ART$^2$ is a 3D multi-wavelength Monte Carlo radiative transfer (RT) code that couples continuum and emission lines to track the propagation of photons and their interactions with the interstellar medium (ISM). The previous version of ART$^2$, which included continuum and Ly$\alpha$ line, has been extensively applied to hydrodynamics simulations to study multi-band properties of galaxies and ISM. Here, we describe new implementations of non-local thermodynamic equilibrium RT of molecular and atomic fine structure emission lines, and the parallelization of the code using a number of novel methods. The new ART$^2$ can efficiently and self-consistently produce a full spectrum that includes both continuum and lines such as [CII], [NII], [OIII], Ly$\alpha$, and CO. These essential features, together with the multi-phase ISM model and the adaptive grid, make ART$^2$ a multi-purpose code to study multi-wavelength properties of a wide range of astrophysical systems from planetary disks to large-scale structures. To demonstrate the capability of the new ART$^2$, we applied it to two hydrodynamics simulations: the zoom-in Milky Way Simulation to obtain panchromatic properties of individual galaxies, and the large-scale IllustrisTNG100 Simulation to obtain global properties such as the line intensity mappings. These products are essential for a broad array of studies such as the correlations between physical and panchromatic properties and their evolution. By enabling direct comparison between numerical simulations and multi-band observations, ART$^2$ provides a crucial theoretical framework for the interpretations of existing observations, the plan for future surveys, and the synergy between multi-band galaxy surveys and line intensity mappings. Therefore, ART$^2$ is a powerful and versatile tool to bridge the gap between theories and observations of cosmic structures.

The "Disk of satellites" (DoS) around Milky Way is a highly debated topic with conflicting interpretations of observations and their theoretical models. We perform a comprehensive analysis of all dwarfs detected in the Milky Way and find that the DoS structure depends strongly on the plane identification method and the sample size. In particular, we demonstrate that a smaller sample size produces a higher anisotropy of the spatial distribution and a stronger clustering of the angular momentum of the satellites. Moreover, we calculate the evolution of the 11 classical satellites with proper motion measurements and find that the thin DoS they currently reside in is transient. Furthermore, we analyze two cosmological simulations using the same initial conditions of a Milky Way-sized galaxy, an N-body run with dark matter only and a hydrodynamic one with both baryonic and dark matter, and find that the hydrodynamic simulation produces more anisotropic distribution of satellites than the N-body one. Our results suggest that an anisotropic distribution of satellites in galaxies can originate from baryonic processes in the hierarchical structure formation model, but the claimed highly-flattened, coherently-rotating DoS of the Milky Way may be a small-number selection effect. These findings may help resolve the contradictory claims of DoS in other galaxies and the discrepancy among numerical simulations.

Observations of globular clusters show that they have universal lognormal mass functions with a characteristic peak at $\sim 2\times 10^{5}\, {\rm{M_{\odot}}}$, but the origin of this peaked distribution is highly debated. Here we investigate the formation and evolution of star clusters in interacting galaxies using high-resolution hydrodynamical simulations performed with two different codes in order to mitigate numerical artifacts. We find that massive star clusters in the range of $\sim 10^{5.5} - 10^{7.5}\, {\rm{M_{\odot}}}$ form preferentially in the highly-shocked regions produced by galaxy interactions. The nascent cluster-forming clouds have high gas pressures in the range of $P/k \sim 10^8 - 10^{12}\, \rm{K}\,\rm{cm^{-3}}$, which is $\sim 10^4 - 10^8$ times higher than the typical pressure of the interstellar medium but consistent with recent observations of a pre-super star cluster cloud in the Antennae Galaxies. Furthermore, these massive star clusters have quasi-lognormal initial mass functions with a peak around $\sim 10^{6}\, {\rm{M_{\odot}}}$. The number of clusters declines with time due to destructive processes, but the shape and the peak of the mass functions do not change significantly during the course of galaxy collisions. Our results suggest that gas-rich galaxy mergers may provide a favorable environment for the formation of massive star clusters such as globular clusters, and that the lognormal mass functions and the unique peak may originate from the extreme high-pressure conditions of the birth clouds and may survive the dynamical evolution.

It has been suggested that the satellite galaxies of the Milky Way reside in a highly-flattened, kinematically-coherent plane called Disk of Satellites (DoS). The origin of the DoS, however, has been hotly debated, and a number of conflicting claims have been reported in the literature on whether or not the DoS is consistent with predictions from the standard Lambda Cold Dark Matter ($\Lambda$CDM) cosmological model. Here we investigate this issue by comparing a high-resolution, hydrodynamic $\Lambda$CDM simulation of a Milky Way sized galaxy with its dark matter only counterpart. We find the following results: (1) The abundance and distribution of satellite galaxies around a host galaxy is significantly different in the hydro simulation compared to its N-body counterpart; (2) No clear coherent rotation is found in the satellite system, as the fractions of corotating and counter-corotating satellites remain comparable across cosmic time; (3) The satellite distribution evolves significantly with time, from nearly isotropic at high redshift to anisotropic at the present day; (4) The DoS properties strongly depend on sample selection and plane identification methods. Our results imply that the spatially-thin and coherently-rotating DoS reported in Milky Way and other galaxies may be a selection effect of small sample size.

We study the impact of baryons on the distribution of dark matter in a Milky Way-size halo by comparing a high-resolution, moving-mesh cosmological simulation with its dark matter-only counterpart. We identify three main processes related to baryons -- adiabatic contraction, tidal disruption and reionization -- which jointly shape the dark matter distribution in both the main halo and its subhalos. The relative effect of each baryonic process depends strongly on the subhalo mass. For massive subhalos with maximum circular velocity $v_{\rm max} > 35 km/s$, adiabatic contraction increases the dark matter concentration, making these halos less susceptible to tidal disruption. For low-mass subhalos with $v_{\rm max} < 20 km/s$, reionization effectively reduces their mass on average by $\approx$ 30% and $v_{\rm max}$ by $\approx$ 20%. For intermediate subhalos with $20 km/s < v_{\rm max} < 35 km/s$, which share a similar mass range as the classical dwarf spheroidals, strong tidal truncation induced by the main galaxy reduces their $v_{\rm max}$. Moreover, the stellar disk of the main galaxy effectively depletes subhalos near the central region. As a combined result of reionization and increased tidal disruption, the total number of low-mass subhalos in the hydrodynamic simulation is nearly halved compared to that of the $\textit{N-}$body simulation. We do not find dark matter cores in dwarf galaxies, unlike previous studies that employed bursty feedback-driven outflows. The substantial impact of baryons on the abundance and internal structure of subhalos suggests that galaxy formation and evolution models based on $\textit{N}$-body simulations should include these physical processes as major components.