Modeling Photoionized Turbulent Material in the Circumgalactic Medium. III. Effects of Corotation and Magnetic Fields

To carry out nonequilibrium MHD simulations of the magnetized CGM, we used Models of Agitated and Illuminated Hindering and Emitting Media (MAIHEM). ³ MAIHEM is a cooling and chemistry package built using FLASH (ver. 4.5), an open-source hydrodynamics code (Fryxell et al. 2000). This package models the hydrodynamics in three dimensions (3D) and evolves a nonequilibrium chemistry network of 65 ions, including hydrogen (H i and H ii), helium (He i–He iii), carbon (C i–C vi), nitrogen (N i–N vii), oxygen (O i–O viii), neon (Ne i–Ne x), sodium (Na i–Na iii), magnesium (Mg i–Mg iv), silicon (Si i–Si vi), sulfur (S i–S v), calcium (Ca i–Ca v), iron (Fe i–Fe v), and electrons. This includes solving for dielectric and radiative recombinations, collisional ionizations with electrons, charge transfer reactions, and photoionizations by a UV background.

MAIHEM was originally developed in Gray et al. (2015), later improved on with the inclusion of an ionizing background in Gray & Scannapieco (2016), and further updated with several charge transfer reactions, radiative recombination rates, and dielectronic recombination rates from Aldrovandi & Pequignot (1973), Shull & Van Steenberg (1982), and Arnaud & Rothenflug (1985) in Gray & Scannapieco (2017). Physical cooling processes are included down to 5000 K.

Most recently, Buie et al. (2020b, hereafter B20) modified the MAIHEM code to simulate this gas in a dark matter (NFW; Navarro et al. 1996) gravitational potential, given by

where is the dark matter density normalization, M_halo is the mass of the halo, and R_s = R_vir/c is the scale radius, with R_vir the virial radius of the halo and c the concentration parameter of the halo. The turbulence, which is used as our driving feedback mechanism, is artificially driven according to the following equation:

where is the original acceleration term resulting from the direct Fourier transform of the solenoidal modes (∇ · F = 0, where F is the artificial forcing term used in the momentum equation; Pan & Scannapieco 2010), r is the radius to a cell, and α is a dimensionless parameter that controls how the acceleration behaves with radius. These solenoidal modes depend on the driving scale of turbulence, total injected energy, and autocorrelation time for the Ornstein–Uhlenbeck (OU) process that these are modeled as (Uhlenbeck & Ornstein 1930).

Gray & Scannapieco (2016) detail the equations solved in MAIHEM that are invariant under the transformation x → λ x, t → λ t, ρ → ρ/λ, meaning that the final steady-state abundances depend only on the product nL of the mean density and the driving scale of turbulence, the one-dimensional (1D) velocity dispersion of the gas, σ_1D, and the ionizing extragalactic UV background (EUVB).

We utilize the unsplit staggered mesh algorithm, which solves multidimensional ideal and nonideal MHD problems on a Cartesian grid, and the equations for these calculations are available in Lee (2013). This study specifically uses the ideal MHD solvers, and we also make use of a hybrid Riemann solver that uses the Harten–Lax–van Leer (HLL) solver (Einfeldt et al. 1991) in places with strong shocks or rarefactions and the Harten–Lax–van Leer–Contact solver (Toro et al. 1994; Toro 1999) in smoother flows to stabilize the code as turbulence ensues. We refer the reader to Gray et al. (2015) and Gray & Scannapieco (2016) for further details.

We conducted a suite of numerical simulations that model an 800 kpc box with periodic boundaries using Static Mesh Refinement to accurately capture important structures that develop in the halo as turbulence ensues. The simulations are run in a Milky Way–mass dark matter halo of 10¹² M_⊙, with a virial radius of 220 kpc and a concentration parameter of 10. Note that the virial radius remains constant with time, and the simulations are not run in a comoving domain.

Cosmological parameters are solely used to simulate the gravitational acceleration of the baryonic matter and are from the Planck 2018 Collaboration (Aghanim et al. 2018). They are h = 0.674, Ω_m = 0.315, Ω_b = 0.049, and Ω_Λ = 0.685, where h is the Hubble constant in units of 100 km s⁻¹ Mpc⁻¹ and Ω_m , Ω_b , and Ω_Λ, are the total matter, baryonic, and vacuum densities, respectively, in units of the critical density.

The resolution was set to that in B20, meaning that for R ≳ 300 kpc the domain was at a resolution of 64³, which translates to 12.5 kpc; for 300 kpc ≳ R ≳ 250 kpc the resolution was 128³, which translates to 6.2 kpc; for 250 kpc ≳ R ≳225 kpc the resolution was 256³, which translates to 3.1 kpc; and for 225 kpc ≳ R the resolution was 512³, which translates to 1.6 kpc.

Note that Buie et al. (2018, hereafter B18) used a smaller, subparsec resolution and a smaller turbulent driving scale as compared to the runs presented here. Both sets of simulations inject turbulence with stirring modes between L/3 ≤ 2π/k ≤ L, but in B18 L was set to be on the scale of the box, 100 pc, while here it is chosen to be 30 kpc. This is about the size of the Milky Way disk, as we assume that such turbulence would be primarily driven by outflow/inflow processes on the scale of the disk. However, it is misleading to match these physical scales directly, because for a fixed ionization parameter the results of the isotropic turbulence simulations, such as presented in B18, depend on the product of the density and the driving scale nL, rather than directly on L. As the densities that were used in this previous work were orders of magnitude larger than those found in the study described here, nL is comparable between the two sets of simulations, even though L is not.

Likewise, it is the resolution in units of L rather than in physical units that is the key measure of resolution. In B18 the resolution was Δx = L/128, while in the current case it is slightly lower, Δx ≈ L/20, where 20 results from the resolution divided by the driving scale of turbulence. However, B20 also conducted a resolution study by doubling the resolution to 1024³ or 0.8 kpc and found that while the cool structures that developed were smaller in scale, the overall properties and conditions in the halo, such as temperatures, densities, and resulting column densities, remained in agreement with the original resolution counterpart. MHD simulations as large as 1024³ that are run for 9 Gyr are prohibitively expensive, and we predict that in such simulations the resulting cool structures would be smaller but the overall properties would be largely similar to the runs presented here.

We initialize the medium with a Keplerian circular velocity modified by the observational findings from the Mg ii studies of Ho et al. (2017). Specifically, we use their Equation (A2) to inform how the velocities should fall off along the minor axis, which is chosen to be the y-axis in our runs. This gives a rotational velocity of

where G is the gravitational constant, is the dependence of the circular velocity on scale height, and f_rot is an overall scaling parameter we set to 60%. In this case the spin parameter of the gas, λ, is defined as

where J is the total angular momentum within a sphere of virial radius R along with its gas mass M_gas. V is the circular velocity at the virial radius, . The chosen f_rot results in λ ≈ 0.03 in our halos, in line with numerical studies of galactic dark matter halos (Bullock et al. 2001; Macciò et al. 2008; Bryan et al. 2013).

In all our runs, the gas was initialized with a fractional ion abundance of 0.3 Z_⊙ in collisional ionization equilibrium (CIE) at the virial temperature of the halo, T = 1.2 × 10⁶ K, and average sound speed of 166 km s⁻¹. Following initialization, the medium is allowed to evolve according to the nonequilibrium chemistry. We note that our simulations assume the metallicity to be constant throughout the halo, while observations support a CGM that is significantly chemically inhomogeneous (Zahedy et al. 2019; Buie et al. 2020a). We also choose a redshift-zero (Haardt & Madau 2012, hereafter HM2012) EUVB whose specific intensity was normalized to 8.23 × 10⁻²⁴ erg cm⁻² s⁻¹ Hz⁻¹ sr⁻¹ at the Lyman limit to irradiate the gas in the runs.

We conduct runs that build off of the previous High run from B20, which was able to match many ion observations of nearby star-forming galaxies’ CGMs. That run was conducted by setting α in Equation (2) to −1, and it used the same driving scale for turbulence, 30 kpc. This ensured that turbulent stirring was strongest toward the center and fell off with radius. Observationally, nonthermal motions that are derived from Voigt profiles are typically smaller for low ionization state ions, which are more often found near the host than in higher ionization state ions, which may reside at a range of radii from the central galaxy (Tumlinson et al. 2013; Churchill et al. 2015; Werk et al. 2016; Faerman et al. 2017). However, this may be largely due to the fact that low ionization state ions occur preferentially in colder material with higher physical densities than higher ionization state ions, meaning that at the same column depth they will probe smaller physical scales, wherein turbulence is much smaller. Additionally, Buie et al. (2020a) conducted a Markov Chain Monte Carlo investigation of absorbers in the CGM of nearby galaxies and found that the turbulent velocities had no definitive radial trend, but instead covered a range in σ_1D between 11 and 60 km s⁻¹, consistent with our turbulent stirring.

The turbulent stirring employed here and in the previous B20 study also allows us to more accurately capture the stronger central feedback of actively star-forming galaxies, as well as the impact of gas accreted from the intergalactic medium, which is likely to deposit the most energy per unit volume near the center of the halo. This choice in turbulent stirring led to an average 1D velocity dispersion σ_1D = 41 km s⁻¹. In addition to this turbulent stirring treatment, the aforementioned corotation prescription was also considered in the run that we call the Hydro Rot run.

The second run we conducted further included a seed magnetic field in the z-direction of 0.1 μG run in ideal MHD. This choice of initial magnetic field strength translated to a plasma β, defined as the ratio of thermal pressure to magnetic pressure, of ≈170. It also used the same turbulent stirring and corotation prescriptions and is called the MHD Rot run.

Previously in B20, the simulations were run for 3 Gyr until they reached a global equilibrium such that the change in total energy was less than 0.25% per time step. With the addition of corotation, we ran them for 9 Gyr as interesting features developed after 3 Gyr, and we wanted to capture the full saturation period of the magnetic field in the MHD run.

Note that these simulations are run in a static domain for 9 Gyr and do not emulate the mass growth that would likely be experienced by a real galaxy in this evolutionary time. Such mass growth would likely lead to increased turbulent motions as new gas enters the CGM. Furthermore, these interactions would likely increase the magnetic field in these areas owing to the twisting of magnetic field lines.

To give the reader visual context for the changes in structure as our halos evolve, we show projections of the mass-weighted number density from face-on and edge-on perspectives for the Hydro Rot and MHD Rot halos in Figure 1, at 3, 6, and 9 Gyr. Furthermore, Figure 2 shows the evolution of the gas mass within R_vir in our new simulations, compared to the High simulation from B20, the mass between 12 kpc and R_vir, and the mass in gas with T < 10⁵ K. As in our previous simulations, gas within r ≈ 30 kpc cools in the initial ≈500 Myr, forming a low-pressure region in the center of the halo, which promotes an accretion flow of gas to the inner 12 kpc, which may be seen from the middle panel of Figure 2. Such cooling may also be inferred from the rapid rise in cool gas observed in the right panel of this figure. Gas near the virial radius is pushed beyond it by the initial turbulence, a feature indicated by a decrease in mass shown in the left panel of this figure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The accretion flow results in an accretion shock that forms near the core and gradually moves outward, dispersing the gas near the core throughout the halo. One observes this feature by looking at the middle panel of Figure 2, which shows a rise in gas mass between 1 and 2 Gyr. We also see the initialization of the turbulence, followed by the accretion shock, both of which push about 14% of the baryonic mass beyond R_vir in the first 2 Gyr.

At approximately 2 Gyr, the accretion shock leaves the halo. Following this, between 2 and 4 Gyr, both of the new runs follow similar evolutionary trends in their mass confined to the intermediate region between 12 kpc and R_vir, with the MHD case showing lower mass during this time. This results primarily from the mass loss that occurs as it leaves the halo, which may be seen from the downward slope in the gold line of the left panel. This may in part be influenced by additional pressure support from the amplified magnetic field, described in more detail in Section 3.3. The amount of gas cooler than 10⁵ also steadily rises in both runs for the remainder of the evolution.

As ambient gas in the MHD halo is at lower densities, it is unable to cool as efficiently as compared to the Hydro Rot run. This is shown by the total cooling in Figure 3, in addition to the lower amount of cold mass found in this run, which may be seen in the right panel of Figure 2. Mass loss is prolonged in the MHD Rot run, continuing on for nearly 7 Gyr. The total gas mass fluctuates at ≲1.1 × 10¹¹ M_⊙ in the MHD Rot run, compared to the Hydro Rot halo, which instead fluctuates between 1.2 × 10¹¹ M_⊙ and 1.3 × 10¹¹ M_⊙ for the remainder of its evolution.

Zoom In Zoom Out Reset image size

We find that the initial corotation, along with the continuously driven turbulence, adds heat to the hot component, thereby raising its pressure and entropy. We see these as subsonic motions on the order of ≈10 km s⁻¹, and they persist throughout the evolution of both runs. This increased pressure and entropy cause hot gas to buoyantly rise beyond the virial radius. Some of this gas travels far enough to reenter the 800 kpc box from opposite sides, which starts at ≈4 Gyr and is a consequence of our periodic boundary conditions. We consider this to be a numerical artifact, and it does not influence the results of this study. The MHD Rot run possesses a hotter ambient medium with higher levels of entropy and pressure; thus, more mass leaves its halo. This rising pressure and entropy are discussed in more detail in the following section.

When looking at the mass-weighted average temperature, shown in Figure 4, we see a positive slope that approaches T_vir in the MHD case and < T_vir in the Hydro Rot case. Gas within the MHD Rot’s halo maintains these high temperatures over a longer duration, approximately 1 Gyr longer than in the case of no MHD.

Zoom In Zoom Out Reset image size

Both runs appear to come out of their cooling minima after nearly 4 Gyr of evolution. Soon after, we also see the mass between 12 kpc and R_vir rise in the Hydro Rot halo. This is largely due to the ejection of material from the central region and somewhat from the reaccretion of material outside R_vir. The ejection of material from the center that causes the rise in mass seen in this region is observed 2 Gyr later in the MHD Rot halo, as the magnetic field helps to keep cold gas contained to the center for longer periods of time. Furthermore, these runs maintain their increased masses in this interim region for the remainder of their evolution. The mass-weighted average temperature trends approach a minimum, with log T being ≈5.6 in the MHD Rot halo and 5.5 in the Hydro Rot halo prior to rising for the remainder of the simulation.

In Figure 5, we show the radial profile of the number density, n, at 3, 6, and 9 Gyr. Gas is found at a variety of densities during these times, with most of the mass tracing the diffuse ambient medium and dense structures, such as clouds and filaments. The Hydro Rot run tends to have a wider range of densities, and at 6 Gyr it shows a large amount of mass with n > 10⁻² cm⁻³ spread across the inner 50 kpc of the halo. In the MHD case, however, mass consolidates to smaller radii, reaching densities above 10⁻¹ cm⁻³, along with possessing mass beyond 125 kpc with n ≈ 10⁻³ cm⁻³.

Zoom In Zoom Out Reset image size

After 6 Gyr, the halos become very dynamic, with angular momentum transfer driving dense, cold outflowing gas throughout, some of which travels beyond R_vir by 9 Gyr, while also having infalling, cooling gas. Both raise the density in the ambient medium by 9 Gyr. Also at this time, we can see that both runs have moved their mass outward, with less of it having n ≳ 10⁻¹ cm⁻³ in the inner 50 kpc and more with 10⁻² cm⁻³ < n < 10⁻⁴ cm⁻³ spread throughout the halos.

In Figure 6, we examine the radial profile of the ratio of cooling times to freefall times, t_cool/t_ff, for these runs at 3, 6, and 9 Gyr, where

E_in is the internal energy, and i denotes the ith ion.

Zoom In Zoom Out Reset image size

This ratio gives insight into the self-regulated balance maintained between gas condensing out of the ambient medium and feedback from the host galaxy in the form of an AGN, star formation, and/or supernovae. Once triggered, such feedback may prevent further condensation from the ambient halo (e.g., Silk 1977; Gaspari et al. 2012; McCourt et al. 2012; Voit et al. 2015, 2017).

The origin of this ratio arises from early theoretical work that attempted to understand thermal instabilities in the galactic context (e.g., Hoyle 1953; Field 1965; Rees & Ostriker 1977). When this ratio becomes too low (i.e., t_cool/t_ff ≲ 1), gas may condense out of the halo and rain onto the central galaxy. As gas precipitates, any one of the aforementioned energetic feedback processes may come into play, which raises this ratio in the ambient halo to ≈10. This cycle facilitates the formation of a multiphase medium around galaxies, especially near the central regions of the CGM.

Others have conducted simulations of the CGM, testing different values of t_cool/t_ff ≲ 10 and their impacts on the development of a multiphase medium (e.g., McCourt et al. 2012; Sharma et al. 2012; Gaspari et al. 2017; Choudhury et al. 2019; Stern et al. 2019); however, Esmerian et al. (2021) produced a bimodal distribution for these ratios in their Feedback in Realistic Environments simulations of Milky Way–like galaxies. Looking at our halos, we find a similar bimodal distribution for t_cool/t_ff, resembling what was found in the Esmerian et al. (2021) simulations. The hot and diffuse ambient medium lies between 1 and 100, while the other traces T ≲ 10⁵ gas and lies between 10⁻⁴ and 10⁻¹. We also find a slightly wider distribution of t_cool/t_ff ratios in the MHD Rot case at 9 Gyr. This run possesses gas with ratios beyond 100 that traces low densities around 10⁻⁵ cm⁻³, lying above and below the central disk-like structure.

We find that earlier in the simulations this extraplanar gas (not shown) is hot, with temperatures ≳ T_vir and t_cool/t_ff ratios ≳10. This is likely a consequence of the turbulent stirring, which we set to be strongest near the center, along with this region not having the additional cooling that results from the further mixing brought on by the corotation of gas. After ≈6 Gyr, densities in these regions increase in the Hydro Rot run. During this stage, matter falls toward the center from above and below the central disk-like structure. This differs from the MHD case, in which infalling matter has more angular momentum and thus funnels through the y = 0 plane, within the extended disk. At 9 Gyr, these regions of cooler gas above and below the disk in the Hydro Rot run extend much farther as compared to the MHD case.

Lastly, we examine the entropy and thermal gas pressure profiles, again at 3, 6, and 9 Gyr, in Figure 7. At 3 Gyr, the ambient medium possesses a maximum K ≈ 100 cm² keV in both runs that grows to about 350 and 500 cm² keV by 9 Gyr in the Hydro and MHD Rot halos, respectively, meaning that unless it is mixed into colder material, its cooling time is much longer than a Hubble time (Voit et al. 2002; Oh & Benson 2003). We also find that the additional heating that results from the continuous stirring and initialized corotation previously mentioned in Section 3.1 causes the ambient medium in both runs to be more constrained to ≈10⁻¹⁴ dyn cm⁻² by 9 Gyr. We further find that the MHD Rot run possesses a pressure profile that spans a smaller range as compared to the non-MHD case at all times. This may also be seen in gas transitioning between the cold and hot phases with entropy values of 0.1–10 cm² keV.

Zoom In Zoom Out Reset image size

Cooling and heating gas also lies across a nearly constant pressure range in both runs. We find that gas with T ≲ 10⁵ K possesses entropies <1 cm² keV and has a wider range of pressures between 10⁻¹⁶ and 10⁻¹³ dyn cm⁻². In both runs, some of the mass is able to cool beneath our temperature floor at 5000 K, more of which is seen in our Hydro Rot run. This is where molecular cooling may become dynamically relevant; however, this is not included in our simulations. By 9 Gyr, both runs show lower-entropy gas with K ≲ 0.1 cm² keV that may be found at lower pressures below 10⁻¹⁵ dyn cm⁻². Further analysis shows us that this low thermal pressure gas has temperatures ranging from ≈10⁴ to 10⁵ K, as well as moderate densities around 10⁻³ cm⁻³. Lastly, we see that the ambient medium possesses a nearly isentropic entropy profile at all times in both of our runs, shown in the Appendix.

Finally, the continuous turbulent stirring, along with the initialized corotation, also raises the pressure in the medium beyond R_vir (not shown). We find that pressures in the medium beyond R_vir fluctuate between 6 × 10⁻¹⁵ dyn cm⁻² and 9 × 10⁻¹⁵ dyn cm⁻², slightly lower than those found in the hot component within the halo. As gas reenters the box, it interacts with buoyantly rising gas between 300 and 400 kpc. The periodic fluctuation in pressure beyond R_vir is not seen in the gas inside R_vir. It is unclear whether these pressures outside of the virial radius have any influence on the observable features within it; however, we do observe a slight increase of the temperature in this distant gas.

Next, we examine the magnetic field in the MHD Rot run and show the magnetic, internal, and kinetic energy densities as a function of time in Figure 8. We remind the reader that this run began with a seed B_z = 0.1 μG, which corresponds to u _B = B²/2μ₀ = 4 × 10⁻¹⁶ erg cm⁻³ and a plasma β of ≈170.

Zoom In Zoom Out Reset image size

The magnetic energy density grows very rapidly as the turbulent gas cools toward the center, amplifying u _B to nearly 5 times its initial value within 500 Myr. For comparison, u _B is able to reach ≈12% of the internal energy density and ≈35% of the kinetic energy density in the CGM at this time. Following this, u _B saturates until ≈7 Gyr, where it settles to about 5% of the kinetic energy density for the remainder of its evolution. ∣B∣ also fluctuates about 0.12 μG at this time.

As before, we take a closer look at the simulation at 3, 6, and 9 Gyr, opting to show the density-weighted projections of ∣B∣ from edge-on and face-on perspectives at these times in Figure 9, overlaid with black quivers that indicate the normalized magnetic field vectors.

Zoom In Zoom Out Reset image size

Higher magnetic field strengths between 1 and 10 μG are found in denser structures such as the cool, extended disk-like gas structure centered at y = 0 and cooling filaments found throughout the halo. We generally find that ∣B∣ decreases more steeply above and below the central disk-like gas structure rather than along the plane of it. ∣B∣ strengths nearing a μG may be found out to r ≈ 100 kpc at 3 Gyr, with these strengths extending to 150 kpc by 9 Gyr as increased amounts of dense filaments and clouds are dispersed throughout the halo during this time. A time of 9 Gyr corresponds to when the halo shows outflowing material along with inflowing cooling gas from beyond R_vir. The outflowing gas is able to carry magnetic energy from the central region out to distant radii, further magnetizing the halo. Although higher B-field strengths may be found at more distant radii during this later time in the MHD Rot run’s evolution, we observe minimal difference in the ∣B∣ radial profile throughout its evolution.

In Figure 10 we show the distribution of average β as a function of height and distance. As the halo evolves, the thermal pressure associated with the ambient medium increases as it heats and becomes more diffuse. This is shown by the increased volume with β ≈ 100, while the gas with β values ≲1 is largely found within the extended disk-like structure at all times. Interestingly, gas above and below this structure and beyond 40 kpc is mostly hot and diffuse; however, it also possesses pockets where it is equally supported by thermal and magnetic pressures. These may be colder voids, as higher temperatures would lead to higher β values.

Zoom In Zoom Out Reset image size

The gas with β values ≲1 may be susceptible to magnetic draping, an effect in which a gas cloud moving through a magnetized plasma is able to rapidly build up a magnetic layer, which may shield it from developing instabilities (e.g., Semenov & Bernikov 1980; Dursi & Pfrommer 2008; Cottle et al. 2020). When looking at the quivers in Figure 9, we indeed see magnetic fields that have oriented themselves to be predominately parallel to the surface of the cool structures, such as the extended central disk-like structure. We speculate that this allows for the persistence of the extended disk-like structure at later times, as mixing directly above and below it is inhibited.

Furthermore, the gas that we observed to have low thermal pressures in Figure 7 tends to be supported by a comparable amount of magnetic pressure as compared to thermal pressure. In fact, we find that the gas with K ≲ 0.1 cm⁻² keV typically corresponds to gas with β values ≲1.

Finally, we compare the kinematics and distributions of ions in the two runs, the properties that are most often constrained by observations. In Figure 11, we show the mass-weighted average temperature as a function of radial velocity, V_radial, and radius at 3, 6, and 9 Gyr, overlaid with the amplitude of turbulence and initial v_rot at y = 0 for reference.

Zoom In Zoom Out Reset image size

In general, we find a smaller spread in radial velocities for the MHD case, as the magnetic field seems to restrict the most extreme motions of the gas. We see that velocities are, on average, ≈10 km s⁻¹ lower than those found in the Hydro Rot run. We can also see how the MHD Rot halo progresses toward hotter temperatures, leading to a difference in the cooling that develops at late times. Both runs develop cool, typically dense outflowing gas structures that originate from the central region at 9 Gyr, albeit more pronounced in the MHD Rot case. This is due to the combination of turbulent driving, which is the only source of feedback in these simulations, and angular momentum transfer, which helps to expel cool gas from the central region.

Figure 12 shows ion-weighted projections of the line-of-sight velocities for these runs, illustrating how the MHD influences the motions of commonly observed ions. From the face-on perspective shown in this plot, we do not have to consider the overall corotation velocity component. Here we find the hot component tracing the ambient medium to be largely static, with low line-of-sight velocities near zero, while the gas within cooling structures traces the fastest-moving material. In particular, the low and intermediate ions, Mg ii and Si iv, show the largest line-of-sight velocities, approaching ±200 km s⁻¹ in the dense cores of cooling structures, while the high ion phase, as traced by O vi, shows smaller velocities overall. The kinematics of the gas is a combination of turbulent motions from the consistent driving, thermal motions from its inherent temperature, and motions due to the gas residing in a gravitational potential.

Zoom In Zoom Out Reset image size

In B20, inhomogeneous turbulent stirring was made to be strongest near the center to mimic the star-forming galaxy environment, with the majority of its feedback resulting from the center, as well as the impact of gas accreted from the intergalactic medium, which is likely to deposit the most energy per unit volume near the center of the halo. This stirring profile promoted a convective flow as hot material was driven to distant radii and replaced by gas condensing out of the ambient halo. When lower degrees of stirring were tested, gas simply collapsed toward the center with very limited mixing between the cold and hot components.

In this study, however, corotation provides a natural mechanism for distributing gas and promoting interactions between cold and hot material. Therefore, a change in the amount of turbulent stirring, or in its radial profile by adjusting α in Equation (2), may lead to different features arising. In the case of stronger stirring, one may observe more outflows and possibly more gas that has been collisionally ionized in the central region, while less stirring may result in less of these. In the case of inhomogeneous stirring that instead increased with radius, one might observe slightly hotter distant ambient gas, as the turbulence would eventually dissipate its energy as heat, along with an increased level of variance in densities found in the outer regions (which may produce changes in the resulting ion column densities in these regions). These considerations may be explored in a future study.

Although not shown here, we find gas moving at faster velocities in the MHD case as compared to the hydrodynamic case when the halos are projected through from an edge-on perspective. It may be the case that magnetic fields help to move gas through the central (y = 0) plane, which results in faster motions of this gas.

Moreover, we find that the fastest-moving clouds also have t_cool/t_ff ratios ≲10⁻¹, meaning that they are cooling at least 10× faster than they are falling, which could potentially lead to star formation in their cores (e.g., Turner et al. 2015). These simulations do not consider star formation (and have a maximum resolution of 1.6 kpc); however, this may have implications for higher-resolution simulations that do consider the formation of star clusters.

To further examine the similarities and differences between the cooling structures found in these two runs, we include the projected column density maps for H i at 9 Gyr from both the face-on and edge-on perspectives in Figure 13. Again, 9 Gyr is a point in the evolution that is characterized by large-scale, fairly dense outflowing and inflowing gas clouds and filaments. We see that there are less cool clouds distributed throughout the MHD Rot run’s CGM as compared to the hydrodynamic case.

Zoom In Zoom Out Reset image size

The smaller amount of cool clouds dispersed throughout the MHD halo may be a consequence of evaporation, as the surrounding medium is hotter in this run as compared to the Hydro case. Furthermore, it has been shown that magnetic fields do not significantly increase the lifetime of cool clouds owing to compression (in the case of transverse fields) or the development of cloud tails (in the case of aligned fields), both of which promote mixing with the ambient medium (Cottle et al. 2020). As both simulations evolve, we find an increasing lower limit of H i columns, as well as a covering fraction of 100% for H i, both seen within ≈50 kpc.

This high covering fraction is most likely a consequence of the idealized nature of our simulations as the amount of cool gas near the central region continues to build without additional ionizing sources of feedback, such as a star-forming UV background, AGN feedback, or cosmic rays. Also, our simulations do not consider radiative transfer, which would provide additional means for heating and possibly ionization. These relevant radiative backgrounds (i.e., star-forming UV background, AGN feedback, or cosmic rays) would provide additional heating to the central region of both halos, likely reducing some of the cooling that takes place in the central region, and even ionizing some of the dense H i that accumulates there.

In Figure 14, we show the Si iv/O vi and N v/O vi ratios for the MHD Rot run at 9 Gyr and a single-zone CLOUDY model for reference. Although we do not show these ion ratios for the Hydro case, we find that both of the run’s ion ratios look very similar to what is shown. Moreover, these ion ratios do not change shape throughout the 9 Gyr evolution, except for being slightly less dispersed at earlier times. In the MHD case, gas with Si iv/O vi ratios ≳1 has temperatures ≲40,000 K and β values ≲1, whereas this gas has no magnetic pressure support in the hydrodynamic case.

Zoom In Zoom Out Reset image size

We also performed a cut on the data to see whether the ratios change considerably if we consider “realistic” outer CGM observations, i.e., those with O vi column densities above 10¹³ cm⁻². We assumed a path length of 100 kpc to obtain a conservative lower limit on n_OVI of 3 × 10⁻¹¹ cm⁻³. This cut essentially removes the most diffuse ambient halo gas that would be unlikely to be observed in the spectra of background quasars. The removal of this gas, however, does not alter the overall shape of these distributions.

The majority of COS-Halos detections and limits presented in Figure 12 of Werk et al. (2016) have Si iv/O vi ≲ 10⁻¹ and N v/O vi ≲ 10⁻¹. While these are difficult to explain from equilibrium models, they naturally correspond to gas with T between 10⁵ and 10⁶ K in our simulations. In the MHD Rot run, this region of parameter space also houses gas with β values ≲1, meaning that some of these systems could trace gas in which the dominant pressure is magnetic.

In Figure 15, we show the distribution of CGM mass as a function of radius and angular momentum per unit mass (otherwise known as the specific angular momentum). We further show a Keplerian specific angular momentum profile for comparison, which we will denote as L_Kep. At 3 Gyr, we find that the majority of the mass in the inner 50 kpc of the Hydro Rot halo lies along L_Kep and the majority of this mass in the MHD case is above the Keplerian profile. Such gas is under the influence of turbulent stirring, along with corotation, and this difference indicates larger velocities in the gas residing in the inner halo of the MHD run. This higher angular momentum gas is found near the plane of the central disk structure, centered at y = 0. At 6 Gyr, we observe that the majority of the mass in the inner 50 kpc has decreased in specific angular momentum in both runs. Given that the mass does not seem to have moved to larger distances, this may indicate a decrease in associated velocities. The mass in this inner region of the MHD halo lies along L_Kep, while the mass in this same region of the Hydro Rot halo lies below it.

Zoom In Zoom Out Reset image size

At 9 Gyr, however, the velocities of the majority of the mass within 100 kpc increase again. At this time, the majority of mass in this inner region lies closer to L_Kep in the Hydro case and at or above the Keplerian profiles in the MHD case, indicating that it has been spun up at this later time in its evolution. Thus, it may be inferred that most of the mass within 100 kpc in the hydrodynamic case is traveling at smaller velocities as compared to the gas in the MHD central disk-like structure. Again, gas with the highest specific angular momentum is found along the plane of rotation at y = 0. We can also see this from the last two rows in Figure 1, which show that in both of our runs the densest (n ≳ 10⁻² cm⁻³) material containing most of the mass in the inner ≈100 kpc region lies within ±50 kpc of y = 0.

Given that this dense gas also houses the largest magnetic fields, it is likely that the magnetic fields help to facilitate this transfer of angular momentum throughout the cold gas. This transfer of angular momentum in the MHD case occurs between 6 and 8 Gyr and is likely caused by the combination of large outflows of v _r ≳ 100 km ⁻¹ material along with infalling condensate from the ambient medium that occurs predominantly along the central disk plane centered at y = 0.

Finally, we note that Hodges-Kluck et al. (2016) found that a spherically rotating model, rather than the thick-disk model adopted here, provided a good match to the observations of O vii absorption lines in the Milky Way halo. Adopting such a model would result in a slightly higher total angular momentum of the system, requiring a smaller f_rot in Equation (3) to obtain a λ near 0.03 for our halo. We speculate that spherical rotation may result in slight changes to the mixing and eventual cooling of hot gas, more so in the gas at the outer regions of the halo. However, it is unlikely that these would lead to substantial changes in the structures, column densities, and other properties, especially since the scale heights of disk models that were compared to a spherical rotational model in Hodges-Kluck et al. (2016) were on the order of a few kiloparsecs, rather than the 50 kpc adopted here.