The impact of magnetic fields on thermal instability
MMNRAS , 1–17 (2017) Preprint 10 July 2018 Compiled using MNRAS L A TEX style file v3.0
The impact of magnetic fields on thermal instability
Suoqing Ji ( 季 索 清 ), (cid:63) S. Peng Oh and Michael McCourt , Department of Physics, University of California, Santa Barbara, CA 93106, USA. Hubble Fellow.
Accepted 0000. Received 0000; in original form 0000
ABSTRACT
Cold ( T ∼ K) gas is very commonly found in both galactic and cluster halos. There is noclear consensus on its origin. Such gas could be uplifted from the central galaxy by galacticor AGN winds. Alternatively, it could form in situ by thermal instability. Fragmentation intoa multi-phase medium has previously been shown in hydrodynamic simulations to take placeonce t cool / t ff , the ratio of the cooling time to the free-fall time, falls below a threshold value.Here, we use 3D plane-parallel MHD simulations to investigate the influence of magneticfields. We find that because magnetic tension suppresses buoyant oscillations of condensinggas, it destabilizes all scales below l coolA ∼ v A t cool , enhancing thermal instability. This effectis surprisingly independent of magnetic field orientation or cooling curve shape, and sets ineven at very low magnetic field strengths. Magnetic fields critically modify both the amplitudeand morphology of thermal instability, with δρ / ρ ∝ β − / , where β is the ratio of thermal tomagnetic pressure. In galactic halos, magnetic fields can render gas throughout the entire halothermally unstable, and may be an attractive explanation for the ubiquity of cold gas, even inthe halos of passive, quenched galaxies. Key words: galaxies: haloes – galaxies: clusters: general – galaxies: evolution – galaxies:magnetic fields
Modern theories of structure formation predict the halos of galax-ies and clusters to be filled with hot virialized gas in approximatehydrostatic equilibrium (e.g, Birnboim & Dekel 2003). This is eas-ily seen directly in X-ray emission at group and cluster scales, andat galaxy scales has been inferred from stacked SZ observations(Planck Collaboration et al. 2013; Anderson et al. 2015) and moreindirect measures such as confinement of gas clouds and stripping ofsatellite galaxies (Maller & Bullock 2004; Fang et al. 2013; Stockeet al. 2013). In recent years, a much greater surprise from quasarabsorption line observations of the circumgalactic medium (CGM)was the ubiquity and abundance of T ∼ K, photoionized gasthroughout such halos, which is seen in more than half of sight-lines within 100 kpc in projection (Werk et al. 2013; Stocke et al.2013). Photoionization modeling suggests reservoirs of cold gas of ∼ − M (cid:12) , which is ∼ −
10% of the halo mass (Stockeet al. 2013; Werk et al. 2014; Stern et al. 2016). At higher redshifts,when quasars are prevalent, such cold gas is also seen in area-filling, smoothly distributed fluorescent Ly α emission (Cantalupoet al. 2014; Hennawi et al. 2015; Cai et al. 2017). Importantly, deepLy α imaging has found such widespread cold gas to be ubiquitous inall surveyed quasars, out to distances ranging from ∼ −
300 kpc.Understanding the origin of this cold gas is an important chal-lenge for several reasons. The CGM is an important discriminant of (cid:63)
Email: [email protected] galaxy formation models, which have largely been tuned to matchobserved stellar masses and galaxy morphology. Furthermore, only T ∼ − K gas will be observationally accessible in detail forthe foreseeable future. Finally, cold gas both provides fuel for starformation and must be expelled in outflows, and thus is arguablythe most important component to track and understand.The theoretical challenge is similar to that in galaxy clustercores, where the observed cold filaments are thought to arise eitherfrom uplifted cold gas from the central cD galaxy, or are made in situvia thermal instability. In galaxies, models where cold gas is upliftedby galactic winds face problems in understanding how gas canbe entrained without being shredded by hydrodynamic instabilities(e.g., Scannapieco & Brüggen 2015; Zhang et al. 2017), leading tosuggestions that such gas must instead instead spawned by thermalinstability in the wind (Martin et al. 2015; Thompson et al. 2016;Scannapieco 2017). It is also hard to understand the prevalence ofdense cold gas in passive galaxies. These problems may potentiallybe alleviated in recent models which invoke small scale structurein cold gas to enable rapid entrainment and suspension of coldgas droplets (McCourt et al. 2018). Nonetheless, thermal instabilityremains an attractive way to produce cold gas in situ and solve theseproblems.Simulations of local thermal instability in stratified systems inglobal hydrostatic and thermal balance showed that fragmentationinto a multi-phase medium happens if the ratio of the local coolingtime to the free-fall time falls below a threshold value of t cool / t ff ∼ −
10. Subsequent simulations by other groups have corroborated © 2017 The Authors a r X i v : . [ a s t r o - ph . GA ] J u l Ji, Oh, McCourt this finding, and have also extended its relevance to another domain,suggesting that black holes are fed by cold accretion, maintaininga feedback loop that pins t cool / t ff at this threshold value (Gaspariet al. 2012, 2013; Li & Bryan 2014; Li et al. 2015; Prasad et al.2015). Observationally, there is also support for the idea that aninterplay between cooling and feedback leads cluster cores to self-regulate at a threshold value t cool / t ff ∼
10 (Sharma et al. 2012;Voit & Donahue 2015). The exact origin of this numerical value forthe threshold is still debated (Voit et al. 2017), and the results ofthis paper will bear directly upon this. These ideas also lead to thenotion that star formation in galaxies may be precipitation regulated(Voit et al. 2015).Observational constraints on magnetic fields in galaxy halosare relatively sparse and uncertain. In the Milky Way, best-fittingmodels indicate a field of 1 − µ G (see Table 1 and section 17.4of Haverkorn 2015). Assuming equipartition between cosmic-raysand magnetic fields, the scale height from synchrotron emission is5 − ≥ from the radio halos of edge-on galaxies (Dumke& Krause 1998; Heesen et al. 2009). Faraday rotation along radioquasar sightlines yield magnetic fields as large as tens of µ G out to ∼
50 kpc distances from galaxies, out to z ∼ ∼ µ G fields inthe halo of the edge-on spiral NGC 4631 (Mora & Krause 2013). At face value, these field strengths are far greater than thepressure of gas in the cold phase (Werk et al. 2014), and comparableto or possibly higher than the pressure of the ambient hot gas (for which indirect measurements and upper limits suggest P hot (cid:46) − erg / cm ; Fang et al. 2013; Anderson et al. 2015). This isin contrast to the previously studied situation in galaxy clusters,where the magnetic field is highly subdominant: β = P gas / P mag ∼
50. Such high magnetic fields may potentially arise from galacticoutflows, the α − Ω dynamo, or magnetic buoyancy of strong fieldsproduced in the galaxy. Regardless, observations suggest that thegas in the outer parts of galaxy halos is magnetically dominated,and therefore that the hydrodynamic calculations in the literaturemay miss essential aspects of the dynamics by neglecting magneticfields.Indeed, strong magnetic fields can fundamentally change thephysics of thermal instability. Local thermal instability in a stratifiedmedium produces overstable gravity waves, whose non-linear sat-uration depends on t cool / t ff , representing the competition betweenthe driving processes of cooling and the damping effect of buoyantoscillations. Magnetic pressure can alter gas buoyancy, while mag-netic tension can suspend overdense gas. Our goal in this paper isto carefully evaluate such effects in 3D MHD simulations. Previouswork has already suggested that magnetic fields can promote ther-mal instability (Loewenstein 1990; Balbus 1991). However, there aretwo important distinctions between this work and that of Loewen-stein (1990) and Balbus (1991). Firstly, the background state theyassume is a cooling flow. In the purely hydrodynamic case, thermal Assuming equipartition, the total (ordered and irregular) B-field scaleheight is at least 3 − α times larger than the ∼ α ≈ − . For more details on magnetic field observations in galaxy halos, see arecent review paper Han (Han). This could suggest that the cold phase is supported by non-thermal pres-sure such as magnetic fields, turbulence or cosmic rays (see Wiener et al.2017 for further discussion of the last possibility), and further highlight theimportance of considering non-thermal components of the plasma. instability cannot develop (Balbus & Soker 1989; Li & Bryan 2012).By contrast, we assume a background state in hydrodynamic andthermal equilibrium, where thermal instability and a multi-phasemedium can develop, if t cool / t ff is sufficiently low (McCourt et al.2012; Sharma et al. 2012). Secondly, those papers perform a linearstability analysis. Since the damping of thermal instability by gasmotions is fundamentally a non-linear process, it is only by run-ning MHD simulations that we can find the non-linear endstate ofthermal instability.The outline of this paper is as follows. In §2, we describeour methods; in §3, we describe our numerical results; and in §4we provide a physical interpretation. Finally, in §5 we discuss thepotential impact of some missing physics in our simulations, andsummarize our conclusions. We use FLASH v4.3 (Fryxell et al. 2000) developed by the FLASHCenter of the University of Chicago for our simulations. A direc-tionally unsplit staggered mesh (USM) MHD solver, based on afinite-volume, high-order Godunov method combined with a con-strained transport (CT) type of scheme, is adopted to solve thefollowing fundamental governing equations of inviscid ideal mag-netohydrodynamics (Tzeferacos et al. 2012; Lee 2013): ∂ ρ∂ t + ∇ · ( ρ v ) = ∂ ρ v ∂ t + ∇ · (cid:18) ρ vv − BB π (cid:19) + ∇ P tot = ρ g (1b) ρ T (cid:18) ∂ s ∂ t + v · ∇ s (cid:19) = H − L (1c) ∂ B ∂ t + ∇ · ( vB − Bv ) = ∇ · B = p tot = p + B /( π ) is the total pressure including bothgas pressure p and magnetic pressure B /( π ) , g is gravitationalacceleration, s = c p ln ( P ρ − γ ) is entropy per mass, and H and L arethe source terms of heating and cooling. We adopt a plane-parallelgeometry with a symmetric domain about z = z axis: g = g z / a (cid:112) + ( z / a ) ˆ z (2)which describes a nearly constant gravitational acceleration for | z | > a with a smooth transition to zero at z =
0. Here the softeningdistance is set to 1 /
10 of the disk scale height. The free-fall timescaleassociated with the gravitational field is: t ff = (cid:115) z g (3)With this gravitational acceleration we set up the initial conditionin hydrostatic equilibrium, where temperature T is constant, thedensity profile ρ is: ρ ( z ) = ρ exp − aH (cid:169)(cid:173)(cid:171)(cid:115) + z a − (cid:170)(cid:174)(cid:172) , (4)and the scale height H = T / g . For convenience, we rescale ourunit system by setting ρ = T = g = k B = µ m p = MNRAS , 1–17 (2017) mpact of B-fields on thermal instability with a resolution of 256 for standard runs and higher resolutionsfor convergence tests. The domain of the simulation is set to be6 × × H . We adopt hydrostatic bound-ary conditions for the boundaries along z direction to maintainhydrostatic equilibrium in the initial conditions when both cool-ing and heating are turned off, and periodic boundary conditions forother boundaries. We initialize isobaric density perturbations with awhite noise power spectrum and amplitude of 1% and wavenumberranging from 1 to 20 in units of 2 π / L box where L box is the domainsize. The magnetic fields are initialized in the horizontal directionwith β ranging from infinity to order unity. We also performedadditional simulations with initially vertical fields, which are dis-cussed in §4.2. These straight fields are initialized to have uniformstrength throughout the box, so in the initial background state themagnetic fields do not exert any forces (and do not contribute tohydrostatic equilibrium). Note that our setup implies that plasmabeta falls with height in the background state. The value we quoteis β evaluated at z = H .It is well-known that halo gas in galaxies and galaxy clus-ters requires some form of “feedback heating” to prevent a globalcooling catastrophe. Precisely how feedback works is still an activearea of research. Hence, we will follow McCourt et al. (2012) andimplement a simple model which preserves thermal energy at eachradius (or z ) in the domain: to maintain global thermal balance, ineach time step we sum up total amount of cooling over the volumeof a thin slab perpendicular to z axis, and distribute heat uniformlyover this volume. In the region of | z | ≤ a , heating and coolingare turned off due to the uncertainty of feedback prescription atsmall radii. This idealization allows us to focus on the physics ofthermal instability. Later works (e.g. Sharma et al. 2012; Gaspariet al. 2013) have shown that results are not sensitive to the specificimplementation of heating.As an approximate fit to slopes in the regimes of heavy elementexcitation and bremsstrahlung respectively (Fig. 1), we implementcluster-like and galaxy-like power-law cooling functions as the fol-lowing: L = n Λ ( T ) = (cid:40) Λ n T / Θ ( T − T floor ) galaxy cluster Λ n T − Θ ( T − T floor ) galaxy halo , (5)where Θ is Heaviside function. These cooling functions correspondto two cooling timescales: t cool = T / n Λ galaxy cluster56 T n Λ galaxy halo . (6)The different prefactors are a result of the different cooling curveslopes (see equations 8 and 9 in Sharma et al. 2012). The T floor (set to T /
20) limit prevents cold gas from collapsing to very small valueswhich we cannot resolve. Physically, our assumption is that oncegas reaches this floor, its cooling time becomes so short that it willcool all the way to T ∼ K . In addition, the temperature couldreach arbitrarily large values under our simple heating prescription,which is not physical. We therefore limit the temperature to somemaximum T ceiling (set to 5 T ). Physically, these limits of T cool and T ceiling correspond to the turning points in cooling function at whichgas can reach at the two-phase semisteady state.The cooling time in multiphase gas can become very muchshorter than the dynamical time. Computational methods which T (K) -23 -22 -21 Λ ( e r g c m / s ) Figure 1.
Approximate power-law fits for cooling curve (dashed line), withindex of − . limit the time step to some fraction of a cooling time thus becomeprohibitively expensive. We therefore implement the “exact” coolingalgorithm described in Townsend (2009) into FLASH code, wherewe solve the operator-split energy equation t cool dTdt = − T ref Λ ( T ) Λ ( T ref ) (7)by separating variables where T ref is arbitrary reference temperature,which has an analytic solution for (piecewise) power-law coolingfunctions. With this implementation, we can simulate rapid cooling,even when the cooling time dips below the simulation time step. We perform a series of simulations with different initial values of t cool / t ff and magnetic field strengths, and we plot the magnitudeof density fluctuations as a function of t cool / t ff in Fig. 2. In thisfigure, each data point is taken from an individual simulation withthe corresponding value of initial t cool / t ff . The magnitude of densityfluctuations is computed by volumetrically averaging over the region0 . H < | z | < . H . It is defined as: δρρ = ρ − (cid:104) ρ (cid:105)(cid:104) ρ (cid:105) , (8)where the bracket represents the spatial averaging over at the sameheight. This quantity is computed from a snapshot of each simula-tion which has been evolved over ∼ t cool , when thermal instabilitysaturates and the computed quantity enters a stable stage. Note thatin our unit system we change the cooling time by changing Λ ;physically, however, this corresponds to changing the gas density. In Fig. 2, we show δρ / ρ as a function of t cool / t ff for various values ofmagnetic field strength and a galaxy cluster-like cooling curve ( Λ ∝ T − / ). The hydrodynamic simulations have a δρ / ρ ∝ ( t cool / t ff ) − scaling, as found by McCourt et al. (2012). Two striking facts emergefrom this plot. For this cooling curve, the MHD simulations retainthe overall scaling δρ / ρ ∝ ( t cool / t ff ) − , albeit with a normalizationwhich increases as the magnetic field strength increases (i.e. plasma MNRAS , 1–17 (2017)
Ji, Oh, McCourt -1 t cool /t ff -2 -1 › δ ρ / ρ fi ( t c oo l / t ff ) − Figure 2. δρ / ρ as a function of t cool / t ff with galaxy cluster cooling curveand different initial strengths of magnetic field characterized by the valuesof β : blue – β = ∞ , orange – β = β = β = β =
4, brown – β =
3. The dashed lines are power-law fitsto simulation data points. The dot markers are simulations with initiallyhorizontal field, star markers are with vertical field. All the quantities aremeasured at z = H . β decreases). Evidently, magnetic fields promote thermal instability.In Fig. 3, we quantify this : δρ / ρ ∝ β − / for both galaxy andcluster-like cooling curves, with a strong uptick signaling runawaycooling at low β for galaxy-like cooling curves. Secondly, magneticfields have an impact even when they are remarkably weak. Fig. 2shows that the MHD case converges to the hydrodynamic case onlywhen β ∼ δρρ ∼ . (cid:18) t cool t ff (cid:19) − hydro0 . β − / (cid:18) t cool t ff (cid:19) − MHD . (9)While these fits are derived for horizontal fields and a cluster-likecooling curve, we shall see that (as hinted in Fig. 3) they are inde-pendent of cooling curve slope and field geometry down to β (cid:38) afew. The amplitude of density fluctuations is important because itsignals whether spatially extended multi-phase structure can de-velop in the saturated state of thermal instability. In Fig. 4, we showthe cold mass fraction (where T < T /
3) as a function of t cool / t ff .For stronger magnetic fields, the threshold value of t cool / t ff at whichmulti-phase gas can form is higher by up to an order of magnitude.Since t cool / t ff varies as a function of radius, this implies that coldgas can be produced over an increasingly wider radial range formagnetized halo gas. In Fig 4, the cold gas fraction apparently does not scale mono-tonically with β at high β (note the cross-over at t cool / t ff = . β =
278 (green), β = We only plot this scaling for t cool / t ff = .
7, but it continues to hold forother values of t cool / t ff . Strikingly, the cold fraction becomes independent of t cool / t ff for low β and a galaxy cooling curve, an effect we soon discuss. β -1 › δ ρ / ρ fi β − / Figure 3.
The β dependence of δρ / ρ , where both δρ / ρ and β are takenfrom the snapshots at t cool / t ff = . z = H . The blue dots are with galaxy cluster cooling curve, and theorange dots are with the galaxy halo cooling curve. apparent contradiction with the monotonic scaling of δρ / ρ ∝ β − / seen in Fig 3. In fact, at higher t cool / t ff (off the plot scale), there ismore cold gas in the β =
278 simulation than in the β =
772 andhydro simulation, where cold gas is undetectable. This is a conse-quence of the somewhat arbitrary definition of cold gas as T < T / β increases (Fig 5),which is the more physically important measure. The cold gas massbecomes sensitive to the exact definition as the PDF narrows. Also,note that while overall trends should be robust, detailed results aresensitive to both the cooling curve shape (which we have simplified)and also our cooling curve cutoff at T /
20. In reality, of course, thedensity PDF should be bimodal.To reinforce this point that magnetic fields can play a key rolein the development of a multi-phase medium, in Fig. 5 we show thePDF of gas density for different values of β , for the same parameters( z = H , cluster cooling curve, t cool / t ff = . β ) clearly broaden the gas density considerably.This is reminiscent of how decreasing values of t cool / t ff broaden thegas density PDF (McCourt et al. 2012), and buttresses the expec-tation from Fig. 2 and Eq. (9) that higher values of t cool / t ff can becompensated with stronger magnetic fields to produce comparabledensity fluctuations, even if (as we shall see) the morphology ofcold gas is completely different. In Fig. 6a – 6d, we show slice plots of the magnitude of densityfluctuations for simulations utilizing the cluster cooling curve, for t cool / t ff = . z = β = ∞ ( hydro ) , , , (cid:104) δρ / ρ (cid:105) RMS = . × − , × − , . , .
27. Note that since thedynamic range of (cid:104) δρ / ρ (cid:105) RMS is different for each plot, the colorbarscale is also different for each plot. Also superimposed are themagnetic field lines. While the mean β in the entire box does notevolve significantly, it is clear that magnetic fields are amplifiedaround overdense regions. The origin of this amplification can bededuced by evaluating the stretching and compressional terms in MNRAS , 1–17 (2017) mpact of B-fields on thermal instability t cool /t ff -2 -1 c o l d f r a c t i o n Figure 4.
The mass fraction of cold gas (defined as T < T /
3) over entirevolume, with galaxy cluster cooling curve and different initial strengths ofmagnetic field characterized by the values of β : blue – β = ∞ , orange – β = β = β =
27, purple – β =
4, brown – β =
3. The solid lines are with cluster cooling curve, and the dashed onesare with galaxy cooling curve. For β = ∞ and β =
772 cases, there is nodetectable cold gas for simulations with t cool / t ff ≥ . -2 -1 ρ -4 -3 -2 -1 P D F Figure 5.
The probability density distribution of gas density measured at z = H , with galaxy cluster cooling curve and different initial strengths ofmagnetic field characterized by the values of β : blue – β = ∞ , orange – β = β = β =
27, purple – β =
4, brown – β =
3. Simulations used here are with cluster cooling curve and at longcooling time ( t cool / t ff = . β = the flux freezing equation (Ji et al. 2016):12 d | B | dt = stretching (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) B · ( B · ∇) v − | B | ∇ · v compression (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) − | B | ∇ · v . (10)For high β , the effects of stretching and compression are com-parable, while in the low β case, when magnetic tension is strong,amplification is weaker and mostly compressional. In particular, for weak initial fields, β correlates strongly with overdensity and candrop by a factor of ∼
10 or more in regions within high overden-sities; thus, magnetic fields can have a more significant dynamicaleffect than the mean or initial β might indicate. Furthermore, asdiscussed in §4 [see equations (22) and (23), and Fig. 12], whatmatters is not β = P / P B but δ P / P B ∼ O ( ) . The magnetic stressesshould not be compared against pressure forces (which are nearlybalanced by gravity), but with the perturbed pressure forces, withwhich they are comparable.In the hydrodynamic case, Fig. 6a, density fluctuations grow atthe largest available scales, and show clear horizontal stratification.This horizontal stratification is characteristic of low-frequency g-modes, where buoyancy forces inhibit vertical motions. By contrast,once magnetic fields are introduced, the morphology of densityfluctuations changes completely. For very weak initial fields, forexample with β =
278 (Fig. 6b), cold gas is vertically oriented andhas much more small-scale structure. As the magnetic field strengthincreases ( β = ,
3, Figs 6c and 6d), the vertical bias persists, butthe density fluctuations appear on increasingly larger scales. Finally,if we re-run the simulation with strong magnetic fields β =
3, butshorten the cooling time by a factor of 20 ( t cool / t ff = .
28, Fig. 6e),then vertically oriented overdense filaments with much more small-scale structure appear. Later, we shall argue that these variationsarise because the maximally unstable modes scale as ∼ v A t cool ;thus, weaker magnetic fields or shorter cooling times result in morestructure at smaller scales. Eventually, as β ∼ Clearly, magnetic fields change the structure of g-modes seen inthe hydrodynamic case. Since magnetic fields introduce anisotropicstresses, one might expect their impact to change with magneticfield orientation. For example, one might reasonably expect that inthe horizontal field case, magnetic tension supports overdense gasand suppresses buoyancy-driven g-modes. This effect should vanishfor vertically oriented fields. Indeed, when we run simulations withthe same initial plasma β and t cool / t ff = . β =
278 (Fig. 7a, compare with Fig. 6b), the densityfluctuations appear larger scale and more planar. Note that the mag-netic fields become more horizontal in overdense regions; there issome degree of magnetic tension support. By contrast, for β = β and t cool / t ff .In Fig. 2, we indicate with stars the result of a few different simula-tions with vertical fields. We see that in every instance (cid:104) δρ / ρ (cid:105) RMS closely corresponds to that derived from the equivalent horizontalfield simulations with the same β and t cool / t ff . This remarkableresult suggests that the total amount of cold gas is insensitive tothe (in general unknown) field orientation. This is extremely sur-prising, given how strongly the appearance of density fluctuationsvaries with field orientation; we return to this result in §4, below. MNRAS , 1–17 (2017)
Ji, Oh, McCourt x z › δρ/ρ fi = 0 . − . . δ ρ / ρ (a) β = ∞ , t cool / t ff = .
7, cluster cooling curve x z › δρ/ρ fi = 0 . − . . δ ρ / ρ (b) β = t cool / t ff = .
7, cluster cooling curve x z › δρ/ρ fi = 0 . − . . δ ρ / ρ (c) β = t cool / t ff = .
7, cluster cooling curve x z › δρ/ρ fi = 0 . − . . δ ρ / ρ (d) β = t cool / t ff = .
7, cluster cooling curve x z › δρ/ρ fi = 2 . − . . δ ρ / ρ (e) β = t cool / t ff = .
28, cluster cooling curve x z › δρ/ρ fi = 2 . − . . δ ρ / ρ (f) β = t cool / t ff = .
7, galaxy cooling curve
Figure 6.
Slice plots of the magnitude of density fluctuations with different initial β values and cooling curves, superposed with magnetic field lines.MNRAS , 1–17 (2017) mpact of B-fields on thermal instability x z › δρ/ρ fi = 0 . − . . δ ρ / ρ (a) β = x z › δρ/ρ fi = 0 . − . . δ ρ / ρ (b) β =
3, cluster cooling curve and initially vertical field
Figure 7.
Slice plots of the magnitude of density fluctuations with t cool / t ff = . In Fig. 8, we show δρ / ρ as a function of t cool / t ff for the galaxy-likecooling curve Λ ∝ T − . For high β ( β = ∞ , , δρ / ρ ∝ ( t cool / t ff ) − β − / independentof magnetic field orientation. However, at low β ( β = ,
4) andfor horizontal fields, density fluctuations appear to saturate witha non-linear amplitude of δρ / ρ > independent of t cool / t ff .At face value, this result implies when magnetic fields approachequipartition, gravity drops out of the problem and a multi-phasemedium will always form, independent of t cool / t ff . Importantly, thisonly happens for horizontal fields: the vertical field results show thesame δρ / ρ ∝ ( t cool / t ff ) − scaling at both high and low β , similarto the cluster cooling curve case. We show in §4, however, thatcloser examination reveals this to be actually a box size effect. InFig. 8, we also show with a star the result of a simulation with thesame effective resolution but double the vertical extent (6 H × H × H ). The density fluctuation has fallen considerably, and is closeto the ( t cool / t ff ) − scaling still obeyed by the cluster cooling curve.Examination of the slice plots for the smaller and larger boxes showthat both only have one dominant filament, with v A t cool > L box inthe smaller box, suggesting that box size effects (§3.6) might play arole. We will return to this point in §4.3.The slice plots of the density field for β = ∞ , ,
27 forthe galaxy cooling curve are similar to that for the cluster caseand are not shown. More interesting are the slice plots for low β .It is interesting to contrast the galaxy cooling curve case, Fig. 6f( β = t cool / t ff = . δρ / ρ = . β = t cool / t ff = . δρ / ρ = . ∼ The saturation amplitude δρ / ρ ∼ δρ / ρ ∼ √ δ − ∼ . δ ∼ T / T floor ∼
20 approaches unity. -1 t cool /t ff -2 -1 › δ ρ / ρ fi ( t c oo l / t ff ) − Figure 8.
With galaxy cooling curve, δρ / ρ as a function of t cool / t ff withdifferent initial strengths of magnetic field characterized by the values of β :blue – β = ∞ , green – β = β =
27, purple – β =
4, brown – β =
3. The dot markers are simulations with initially horizontal field andstar markers vertical field. β = ,
27 cases (Fig. 6b, 6c), where multiple such structuresare visible. However, the cold gas in the galaxy cooling curve casereaches much higher peak overdensities, and thus has a much morepronounced filamentary structure. Finally, the contrast between Fig.6f and Fig. 6e ( β = t cool / t ff = .
28) is instructive. Althoughboth have similar RMS density fluctuations and cold gas masses,they have very different topologies. Fig. 6e shows much more smallscale fragmentation and appears more similar to Fig. 6b, albeit witha higher density contrast. We will discuss this in §4 in light of thecharacteristic scale v A t cool .Except at low β (discussed further in §4.3), these results areconsistent with a picture where the maximally unstable length scale MNRAS , 1–17 (2017)
Ji, Oh, McCourt (which we shall later identify as ∼ v A t cool ) depends only on β and t cool / t ff , independent of the shape of the cooling curve. We show the results of convergence tests on a smaller box in Fig. 9.Our fiducial simulations are of size 6 H (where H is the scale height)and 256 , and thus have ∼
42 cells per scale height. It is also usefulto phrase resolution in terms of the length scales: l coolhydro ∼ c s t cool ∼ H t cool t ff (11) l coolA ∼ v A t cool ∼ H β / t cool t ff , (12)where we have used H ∼ c s t ff and β ∼ ( c s / v A ) . We discuss therelevance of these length scales further in §4.1. While these lengthscales are generally well resolved, in some regions of parameterspace where t cool / t ff < β (cid:29)
1, they are not. This applies to datapoints in the short cooling time limit at high β in Fig. 2, where oursimulations may be less reliable.In Fig. 9 we present the results of convergence studies for bothhigh ( β = β = β for a fixed value of t cool / t ff = . β =
278 with ∼ , , H (or ∼ , ,
54 cells per l coolA ). Note thatunlike our fiducial simulations, here we only simulate the upper half( z >
0) of the box. We see that while the detailed morphology of theslice plots changes (with more small scale detail apparent at higherresolution), the amplitude of density fluctuations is stable at the ∼ −
20% level (and similar to that in our fiducial box, (cid:104) δρ / ρ (cid:105) = . l coolA , which is well resolved. Similarly, in Fig.9b, we show the results of simulations for β =
3, with ∼ , H (or ∼ ,
42 cells per l coolA ), which also show convergedvalues for (cid:104) δρ / ρ (cid:105) despite differences in detailed morphology. Notethe differing box size for these convergence tests ( H × H × H inFig. 9a, 6 H × H × H in Fig. 9b). As is apparent in Fig. 9b, thischange is necessary because of the much larger unstable modes inthe low β case. We will explore this further in the following section.We have also run several 2D simulations and compared themto our 3D results. In general, we find that 2D simulations producehigher density fluctuations by a factor of 2 – 3. 2D simulationswere already seen to produce somewhat higher density fluctuationsin the hydrodynamic simulations of McCourt et al. (2012); theeffect is more significant in MHD simulations. The reduced degreesof freedom constrains gas motions, particularly in the anisotropicMHD case, resulting in reduced damping of thermal instabilities.For this problem, it is important to use 3D simulations. In Fig. 10, we show the slice plots from simulations of differentbox size (8 × ×
2, 4 × ×
2, and 1 × ×
2) but with the sameeffective resolution ( ∼
64 cells per H , or ∼
32 cells per l coolA ),for t cool / t ff = . β = δρ / ρ for the largest box is within ∼
20% ofthat of our fiducial simulation. However, for smaller boxes at fixedeffective resolutions, δρ / ρ increases monotonically, and is doublethe fiducial value in the smallest box! This is despite the fact that the slice plots appear broadly similar, with ∼ l coolA / L box isan important dimensionless parameter: we require l coolA < L box for δρ / ρ to be independent of box size (for instance, in Fig. 10c, L box < . l coolA , fails this criterion and thus produces larger overdensities). The preceding results show that magnetic fields impact both the am-plitude and morphology of thermal instability in a stratified medium.The amplitude δρ / ρ ∝ β − / appears to be independent of bothfield orientation and cooling curve shape. However, the orientationof condensed filaments depends on both the amplitude and orien-tation of magnetic fields. Furthermore, the characteristic size offilaments changes with both plasma β and t cool / t ff . We now providea physical interpretation of these results. In order to understand the preceding results, it is useful to con-sider the dispersion relation for a magneto-gravity wave (Stein &Leibacher 1974): ω = (cid:18) k ⊥ k (cid:19) ω + k v , (13)where ω ∼ g / H ∼ t − is the Brunt-Väisälä frequency, k ⊥ in-dicates transverse wave numbers perpendicular to the direction ofgravity, and k B ≡ k · ˆ b is the wavenumber parallel to the backgroundfield direction. The two terms on the right-hand side represent g-mode and Alfvén-wave oscillations, where the restoring forces arebuoyancy forces and magnetic tension, respectively. At first blush,it would appear from equation (13) that interchange modes (withvariations perpendicular to the background field), should be unaf-fected by magnetic fields, and could play an important role. Weshall soon see that our results are remarkable insensitive to the ori-entation of magnetic fields with respect to gravity. We can considerthermal instability to be a process which is driven at a frequency of ω drive ∼ t − , and is modified by the two processes represented bythe right hand side of Eq. (13), g-modes driven by buoyancy forcesand Alfvén waves driven by magnetic tension forces. In the absenceof non-linear effects, this simply leads to exponentially growingoverstable modes. Below, we shall argue that non-linear dampinghas effective frequencies which scales as the linear oscillation fre-quencies.We can estimate the magnitude of density fluctuations in thehydrodynamic case as follows (McCourt et al. 2012). The linearizedmomentum equation: ρ d v dt = δρ g (14)gives a characteristic velocity of v ∼ δρρ g t buoy ∼ δρρ Ht ff , (15)where we have used t buoy ∼ t ff for the buoyant oscillation periodand g ∼ H / t , where H is the pressure scale-height. The condition t drive ∼ t damp yields the fastest growing perturbations if we max-imize t damp ∼ L / v , setting L ∼ H . Setting t drive ∼ t cool , we thusobtain: δρρ ∼ A (cid:18) t cool t ff (cid:19) − (16) MNRAS , 1–17 (2017) mpact of B-fields on thermal instability x z › δρ/ρ fi = 0 . › δρ/ρ fi = 0 . › δρ/ρ fi = 0 . -10 -1 -1 δ ρ / ρ (a) β = x z › δρ/ρ fi = 0 . › δρ/ρ fi = 0 . -10 -1 -1 δ ρ / ρ (b) β = Figure 9.
Convergence tests. Panel (9a): domain size of H × H × H , and resolutions of 64 × × × ×
256 and 256 × ×
512 from left to rightcolumns, at t cool / t ff = .
7. Panel (9b): domain size of 6 H × H × H , and resolutions of 256 × ×
128 and 512 × ×
256 from left to right columns, at t cool / t ff = . in good agreement with our numerical results (Eq. (9); Fig. 3),where in planar geometry the dimensionless coefficient A ∼ . l / v A < t cool ,which implies that only length scales: l < l coolA ≡ v A t cool (17)are destabilized by magnetic fields. Note that in the hydrodynamiccase, the largest scale modes are the most unstable, whereas in theMHD case, it is the small scale modes which grow fastest. Thisaccounts for the distinct change in morphology between hydro-dynamic and MHD simulations, even at high β . Contrast Fig. 6a (hydro) and Fig. 6b ( β ∼ l A ∼ H /√ β are destabilized.The length scale l coolA ∼ v A t cool emerges from balancing the first andthird terms of Eq. (13), assuming ω ∼ ω drive ∼ t − , correspond-ing to a critical length scale where magnetic stresses affect thermalinstability. By contrast, the length scale l A ∼ H /√ β emerges frombalancing the second and third terms of Eq. (13), indicating whenmagnetic tension overwhelms buoyancy forces. Physically we re-quire that magnetic tension beats both cooling and buoyancy forces.These two differ by a factor of t cool / t ff : l coolA ∼ v A t cool ∼ c s β − / t cool ∼ H β − / t cool t ff ∼ l A t cool t ff . (18)Our simulations show that the length scale l coolA (Eq. (17)) is therelevant one, not l A (contrast Fig. 6d and Fig. 6e, where l A is thesame in both simulations since β is the same, but l coolA is smaller bya factor of ∼
20 in Fig. 6e, due to the smaller value of t cool / t ff . Thelatter clearly shows much more small scale structure). MNRAS , 1–17 (2017) Ji, Oh, McCourt x z › δρ/ρ fi = 0 . δ ρ / ρ (a) x z › δρ/ρ fi = 0 . δ ρ / ρ (b) x z › δρ/ρ fi = 0 . δ ρ / ρ (c) Figure 10.
Over-cooling effect: simulations of t cool / t ff = . β = H ,but different box widths of 8 × ×
2, 4 × × × × From these considerations, we can understand how magneticstresses change the amplitude of density fluctuations. Tension forcesbalance gravity when δρ g ∼ B λ B , (19)where the magnetic tension force is ( B · ∇) B = ˆ nB / λ B , wherethe unit vector ˆ n points to the center of the radius of curvature, and λ B is the radius of curvature. Note that this is a non-linear criteria,where δ B ∼ B . Setting λ B ∼ l coolA , the maximal scale which canbe influenced by magnetic tension during thermal instability, weobtain: δρρ ∼ A (cid:18) t cool t ff (cid:19) − β − / (20)which is consistent with our numerical results (Eq. (9); Fig. 2 andFig. 3) which show δρ / ρ ∝ β − / .Comparing Eq. (16) and Eq. (20), we see that the MHD densityfluctuations approaches the hydrodynamic case when β ∼ (cid:18) A A (cid:19) ∼
900 (21)where we have taken A , A from Eq. (9). From Fig. 2, we seethat this is indeed the case. Thus, magnetic stresses start to becomeimportant even when the magnetic energy is highly subdominant.The characteristic scale l coolA ∼ v A t cool also explains why largeboxes are needed for converged results in the low β , long coolingtime limit. If the box is smaller than l coolA , then all modes within thebox are destabilized by magnetic fields (for instance, Fig. 10c, where L box < . l coolA ); gas motions are mostly strongly quenched for thesesmall scale modes. Since damping of gas motions is artificiallyboosted, density fluctuations are larger for small boxes. Only whenthe largest unstable mode l coolA is contained within the box do theRMS density fluctuations become independent of box size.It is useful to directly examine evidence from the simulationsthat magnetic tension from horizontal magnetic fields suppressesbuoyant oscillations. In Fig. 11, we present tracer particle plotsfor particles which manage to eventually cool for a range of initialheights above mid-plane. These particles are selected so that themaximum entropy decrease between t = t cool is maximized.We show both the height above mid-plane ( z -coordinate) as a func-tion of time, and spatial trajectories projected on to the x - z plane(shifted to x = t = | z | <
1, where t cool / t ff is low, become overdense andsink. Note that entropy changes take place over many multiples ofthe cooling time, since the departure from thermal equilibrium issmall. As expected from the stable background entropy stratification(and as is also clear from Fig. 6a), particles execute mostly planarmotions. Once a horizontal magnetic field is introduced, we can see that In a stably stratified medium, v / v ∼ ( v / ω BV H ) (cid:28) , 1–17 (2017) mpact of B-fields on thermal instability t/t cool z e n t r o p y x (shifted) z (a) β = ∞ t/t cool z e n t r o p y x (shifted) z (b) β = t/t cool z e n t r o p y x (shifted) z (c) β =
3, horizontal magnetic fields buoyant oscillations are strongly suppressed (Fig. 11b and Fig. 11c).Instead, particles cool and lose entropy at a fixed height until theyreach sufficient over-densities that they cannot be supported by mag-netic tension. At that point, they fall. As the magnetic field strengthis increased, buoyant oscillations are more strongly damped, andparticles cool earlier at higher z , where t cool / t ff is larger (compareleft-hand panels of Fig. 11b and Fig. 11c). When the field is weak,the characteristic length scale l coolA ∼ v A t cool is small. Particles onlytravel a small horizontal distance in a cooling time before they can-not be supported by the magnetic field and fall (right-hand panelof Fig. 11b). This accounts for the small scale, vertically-oriented morphology of overdense filaments (Fig. 6b). By contrast, whenfields are strong ( β = l coolA is large, and particles can travel aconsiderable horizontal distance as they cool (right-hand panel ofFig. 11c), resulting in filaments which are not necessarily verticallyoriented (Fig. 6d).We can also directly verify the importance of magnetic ten-sion in supporting overdense filaments, by examining terms in thelinearized, perturbed Euler equation. In Fig. 12a, we compare thegravitational force on an overdense filament with the perturbedpressure gradient ∂ z δ P tot (where P tot = P g + P B ) and the perturbedmagnetic tension forces ( B · ∇) B z / π , for the β = MNRAS , 1–17 (2017) Ji, Oh, McCourt t/t cool z e n t r o p y x (shifted) z (d) β = t/t cool z e n t r o p y x (shifted) z (e) β =
3, vertical magnetic fields
Figure 11.
Selected particle temporal trajectories colored by entropy (left-hand panel) and spatial trajectories (right-hand panel, projected to x - z plane andhorizontally shifted to x = t = t ∼ t cool (dotted verticalline) at different heights, which helps to trace the pathway of cold gas formation. All the cases presented above are at long cooling time ( t cool / t ff = . tal field case. We note that the magnetic tension term correspondsclosely in morphology to the δρ g term, providing most ∼
70% ofthe support against gravity. On the other hand, pressure gradientscontribute a non-negligible amount ( ∼ l coolA ∼ v A t cool can explain 3 interesting features seen in our simulations:the dependence of filament size with β and t cool / t ff ; the existence ofa minimum box size for converged results; and mostly importantly,the scaling of RMS density fluctuations with β . In the preceding subsection, we made the case that magnetic fieldtension damps buoyant oscillations by allowing for magnetic supportof overdense gas. Naively, one might expect that this effect willvanish for vertical fields, and thermal instability will no longer beenhanced. Instead, we saw that the amplitude of density fluctuationswas roughly independent of field orientation, at fixed β and t cool / t ff .Why is this so?In Fig. 11d and Fig. 11e, we show particle trajectories forthe β = , t cool / t ff = .
7. In both cases, buoyant oscillations are suppressed, with amplitudescomparable to the horizontal field case (compare with Fig. 11band Fig. 11c). However, the spatial anisotropy of motions appearsto be reversed: motions are horizontally (vertically) biased for theweak (strong) field case respectively, opposite to what is seen inthe horizontal field case. This squares with the horizontal (vertical)orientation of overdense filaments in the weak (strong) vertical fieldcase (Fig. 7), compared to the opposite orientation in the horizontalfield case (Fig. 6b and Fig. 6d).Why are buoyant oscillations suppressed, even for verticalfields? Consider an overdense blob falling under the influence ofgravity. Ram pressure leads to a high pressure region upstream ofthe blob. In the hydrodynamic (and horizontal field) case, this highpressure region is dispersed by lateral gas flows, where unlike inthe vertical direction, there are no restoring forces. However, in thevertical field case, horizontal motions are restrained by magnetictension forces. The high pressure region expands, bending the fieldlines until the horizontal pressure gradients balance magnetic ten-sion forces. We have verified this force balance explicitly in oursimulations. However, this now implies (i) vertical pressure gradi-ents (ii) bent field lines offset from vertical which can exert magnetictension forces. Both of these can support overdense blobs againstgravity and thereby suppress buoyant oscillations.
MNRAS , 1–17 (2017) mpact of B-fields on thermal instability x z δρ g z δP tot ( B · ∇ ) B z / (4 π ) -10 -1 -10 -2 -2 -1 (a) β = x z δρ g z δP tot ( B · ∇ ) B z / (4 π ) -10 -1 -10 -2 -2 -1 (b) β = x z δρ g z δP tot ( B · ∇ ) B z / (4 π ) -10 -10 -1 -1 (c) β =
3, horizontal fieldMNRAS , 1–17 (2017) Ji, Oh, McCourt x z δρ g z δP tot ( B · ∇ ) B z / (4 π ) -10 -10 -1 -1 (d) β =
3, vertical field
Figure 12.
Slice plots of the terms of over-density gravity, total pressure gradient and magnetic tension.
We can write the vertical and horizontal force balance equa-tions as: ∂ z δ P tot + ( B · ∇) B z π = δρ g (22) ∂ ⊥ δ P tot = ( B · ∇) B ⊥ π , (23)where δ P tot is the total hydromagnetic pressure perturbation. Toorder of magnitude, B ⊥ / π ∼ δ P tot . Thus, δ z δ P tot ∼ δ z ( B ⊥ / π ) and Eq. (22) can be written as ( B · ∇)( B ⊥ + B z ) π = ( B · ∇) B π = δρ g , (24)which is identical to the force balance equation for horizontal fields.This is why δρ / ρ is similar for both horizontal and vertical fields: forthe vertical field case, the lateral confinement of the over-pressuredregion associated with a falling over-density creates vertical stresseswhich are comparable to those if field were horizontal, and whichsuppress buoyant oscillations. This also implies that a similar lengthscale l coolA ∼ v A t cool is also singled out.We can also understand the change in morphology of the fila-ments as β decreases. When β is high, lateral confinement is weak,and gas motions are largely horizontal due to the strong stratifica-tion imposed by buoyancy forces, similar to the hydrodynamic case.Thus, filaments are mostly horizontal (Fig. 7a). However, for strongmagnetic fields (low β ), an over-pressured region remains laterallyconfined, and instead expands in the vertical direction. Thus, fila-ments are largely vertical (Fig. 7b). One can think of the magneticfields as semi-rigid walls, across which gas is not free to flow. Sincehorizontal pressure balance is no longer required, alternating stripsof high and low pressure gas develop, corresponding to high andlow density strips respectively, each in a separate state of hydrostaticequilibrium.In Fig. 12b, we show the perturbed forces for the vertical fieldcase. We see that for high β , ∂ z δ P tot largely matches δρ g . Quan-titatively, perturbed pressure gradient dominate vertical magnetictension support by a factor of 5. The field lines are bent to roughlyhorizontal around overdense regions, and the influence of magnetictension is indirect (by confining over-pressured regions). By con-trast, for β = As previously discussed in §3.4, the amplitude of density fluctu-ations is independent of the cooling curve, and is completely de-termined by β and t cool / t ff , as in Eq. (9). These two parametersspecify l coolA ∼ v A t cool ∼ H /√ β ( t cool / t ff ) in the hot medium, whichsets the balance between cooling and magnetic stresses in the hotmedium. The interesting exception for volumetric heating (where Λ ∝ T α ) is the apparent independence of δρ / ρ to t cool / t ff at low β . As mentioned in §3.4, this is a box size artifact; RMS densityfluctuation drops close to the ( t cool / t ff ) − scaling in a larger box.Interestingly, this only happens for the galaxy cooling curvewith horizontal fields at low β . Examination of the box shows asingle dominant filament with a highly rarefied and overheated hotmedium, with a significantly long cooling time such that l coolA ∼ v A t cool > L box . The density and entropy PDFs become bimodal,indicating genuine fragmentation into a multi-phase medium. Bycontrast, in the larger box with identical t cool / t ff and β , the hotmedium does not evolve significantly, with l coolA < L box ; therefore,the entropy and density PDFs remain unimodal.We conclude that in a system with sufficient separation ofscales, then the saturated state of thermal instability obeys the scal-ing of δρ / ρ ∝ β − / ( t cool / t ff ) − (Eq. (9)). However, systems where l coolA > ∼ H are prone to “overcooling”.Under isobaric conditions, emissivity fluctuations scales as δ(cid:15) / (cid:15) ∝ ( α − ) δ T / T , implying stronger emissivity fluctuations forthe galaxy ( α = −
1) as opposed to cluster ( α = .
5) cooling curves.In a finite box, the former is prone to a thermal runaway where l coolA grows beyond the box size and all available modes are destabilizedby magnetic fields.We speculate that this apparent box size effect may be of directphysical relevance to galaxy halos. It is common to find β ∼ few, Note that in the initial hot medium, l coolA ∼ . H < L box ; l coolA subse-quently increases, largely because of the increase in cooling time.MNRAS , 1–17 (2017) mpact of B-fields on thermal instability t cool / t ff > ∼
10 in such halos. By the usual hydrodynamic criteria t cool / t ff <
10 (Sharma et al. 2012), such halo gas should be ther-mally stable. With the addition of magnetic fields, not only is such aconfiguration thermally unstable, but because l coolA > ∼ H (and even l coolA > ∼ r ), it can becomes violently thermal unstable, approaching astate where all large scale modes are thermally unstable and gravitydrops out of the problem. Essentially, with strong magnetic fields,thermal instability is no longer quenched by non-linearity inducedby falling overdense blobs, and proceeds with the same vigor as inan unstratified medium.This has obvious implications for manufacturing cold gas inhalo outskirts, where it is seen (e.g. in Ly α florescence; Borisovaet al. 2016), and in circumstances where thermal instability wouldotherwise not be an obvious culprit. In this first study, we have necessarily made some simplification orleft out some ingredients. These include:(i)
Thermal conduction
Thermal conduction quenches thermal insta-bility parallel to the magnetic field on scales below the Field length, λ F ≈ (cid:112) κ ( T ) T / n Λ ( T ) , where κ ( T ) ≈ n e k B v e λ e is the conductioncoefficient, and v e , λ e are the electron thermal velocity and meanfree path respectively. If we compare this to the characteristic scale λ coolA we have found for thermal instability, we get (cid:32) λ F λ coolA (cid:33) = (cid:18) c s v A (cid:19) m p m e λ e c s t cool = . f β Λ − T / , (25)where f is a suppression factor accounting for reduction in con-duction below the Spitzer rate, due perhaps to small scale magneticfield tangling or electron scattering by plasma instabilities, and β ≡ β / Λ − = ( Λ ( T )/ − ) erg s − cm , and T = ( T / ) K.Since λ F < λ coolA , thermal conduction does not obviously play asignificant role, although the separation of scales is not very largeand the true impact of conduction will have to be assessed in sim-ulations with field-aligned conduction. Note that the suppressionfactor f is unknown and potentially large; there is no evidence forthermal conduction anywhere in dilute plasmas apart from the solarwind, which has a large scale, ordered magnetic field (Bale et al.2013).(ii) Geometry
We have only simulated planar geometry, and the impactof spherical geometry could be important, once λ coolA becomes largeand approaches the radius r . Indeed, once λ coolA / r ∼ λ coolA / H ∼( t cool / t ff ) β − / >
1, we might expect magnetic fields to destabilizeall available modes, similar to the “over-cooling effect” seen insmall boxes with planar geometry (see §3.6). This will have to Previous such simulations have found little difference in simulations withand without anisotropic conduction (McCourt et al. 2012). Note that thissimulations were initialized with extremely weak fields ( β ∼ ), so theydid not see the dynamical effects we report here. In addition, Wagh et al.(2014) employed anisotropic thermal conduction and tested different fieldgeometry including tangled fields, finding little difference between hydroand MHD cases as well. This is most likely because the field strength isalso weak ( β ∼ ), and v A t cool is much smaller than the length scalesdiscussed here. be explored in future work. Note that in the hydrodynamic case,spherical geometry was previously thought to boost the thresholdvalue of t cool / t ff required for a multi-phase medium to develop dueto differing compression on a sinking fluid element (Sharma et al.2012), but it has subsequently been shown that once an apples-to-apples comparison of where the cold gas appears is performed, thereis little to no difference between spherical and Cartesian geometry(Choudhury & Sharma 2016). The latter authors found instead thatthe amount of cold gas depends on how t cool / t ff varies with radius.(iii) Magnetic field topology
We have seen that straight horizontal andvertical magnetic fields have comparable effects in enhancing ther-mal instability. This suggests that magnetic field orientation doesnot play a significant role and that a tangled field would give similarresults, although the case λ B < λ coolA (where λ B is the magneticcoherence length) deserves more careful study. It may well be thatthe characteristic length scale should be min ( λ A , λ coolA , λ B ) .(iv) Hydrostatic and thermal equilibrium
Our model is by constructionalways in global hydrostatic and thermal equilibrium. While this isclearly an excellent approximation for galaxy groups and clusters,it is potentially more questionable at galaxy scales, particularly forlower mass “cold mode” halos which may lack a hot hydrostatic halo,or in systems where galactic winds drive a strong outflow. Whilelocal thermal instability is suppressed in a cooling flow, even weakmagnetic fields enable thermal instability to once again operate(Loewenstein 1990; Balbus 1991; Hattori et al. 1995). It would beinteresting to study if magnetic fields can play a similar role inmodifying thermal instability in galactic winds. Thermal instabilityin galactic winds has been the subject of much recent interest asa means of explaining cold gas clouds outflowing at high velocity;since they are born at high velocity, they can evade the destructionby Kelvin-Helmholtz instabilities which would take place if theywere accelerated by the wind (Martin et al. 2015; Thompson et al.2016; Scannapieco 2017).(v)
Large density perturbations
Our simulations start with tiny densityperturbations, from which thermal instability develops and is inde-pendent of a sufficiently small magnitude of initial perturbations.However, the density perturbations in galaxy halos might be po-tentially large due to strong turbulence or galactic winds. Pizzolato& Soker (2005) and Singh & Sharma (2015) suggest that thermalinstability can develop even at large t cool / t ff when initial densityperturbations are large, but this regime has not been explored in ourstudy. Therefore, it is necessary in future studies to relax the assump-tions about hydrostatic equilibrium by including driven turbulenceinto simulations, and investigate how initially large density pertur-bations produced from strong turbulence can modify the evolutionof thermal instability in MHD cases. In this paper, we explore the effects of magnetic fields on thermalinstability in a stratified medium with 3D MHD simulations. Astro-physical plasmas are invariably magnetized and our main messageis that this is a regime where pure hydro is clearly insufficient; MHDeffects simply cannot be swept under the rug. The main applicationwe have in mind is the formation of cold gas in galaxy ( β ∼ few)and cluster ( β ∼
50) halos, but these results should be broadly ap-plicable to other environments. We find that magnetic fields have astriking impact on both the amplitude and morphology of thermalinstability, even at surprisingly low magnetic field strengths. Thisgreatly broadens the radial range in which halo gas can be thermallyunstable. It is a potentially promising explanation for the ubiquityof cold gas in absorption line spectra of the CGM.
MNRAS , 1–17 (2017) Ji, Oh, McCourt
Our main results are as followings:(i) Magnetic fields have a strong influence in enhancing the amplitude of density fluctuations due to thermal instability, and its abilityto fragment into multi-phase medium This enhancement appearseven at very weak magnetic fields, and the MHD results match thehydrodynamic results only at β ∼ δρ / ρ ∝ β − / ,where β ≡ P gas / P B [see equation (9)]. This enhancement arisesbecause MHD stresses weaken the effective gravitational field, andthus reduce the magnitude of buoyant oscillations (which otherwisedisrupt condensing cold gas) by supporting overdense gas againstgravity. This can be seen in both particle tracer plots and slice plotsof linearized perturbed forces.(ii) Magnetic fields have a strong influence on the morphology ofcold gas filaments. The strength and orientation of magnetic fieldsstrongly influence the size and orientation of cold filaments, evenfor very weak fields, which can be understood in terms of the MHDstresses. In particular, the spacing between filaments scales with thecharacteristic scale l coolA = v A t cool (see below).(iii) Thermal instability in the presence of magnetic fields imprintsa characteristic scale l coolA = v A t cool on cooling gas. This scaleis evaluated for conditions in the hot medium (unlike the shatteringlength c s t cool , which is evaluated for conditions in the cool medium;McCourt et al. 2018), since it determines the characteristic scalein the hot medium at which buoyant oscillations are suppressedand gas cooling is enhanced. The appearance of this characteristicscale can explain three features in our simulations: the scaling ofdensity fluctuations δρ / ρ ∝ β − / [by setting the the scale on whichmagnetic tension operates; B / λ coolA ∼ δρ g , Eq. (19)], the scalingof filament sizes with β and t cool / t ff in simulations, and the “over-cooling” effect seen when the box size is smaller than this scale,implying that all modes in the box are destabilized by magneticfields.(iv) To zeroth order, the amplitude of density fluctuations and cold gasfraction are independent of magnetic field orientation and coolingcurve shape. The independence from field orientation arises becausemagnetic tension can be shown to either support overdense gasdirectly (for initially horizontal fields) or indirectly (by confiningpressurized hot gas, which in term supports the overdense gas ortilts the field lines so that magnetic tension support is important).The independence from cooling curve arises because the effect ofmagnetic fields essentially depends on ∼ B / λ coolA , where λ coolA is evaluated in the hot phase, not the cooler phase. An importantexception is the case when mass dropout is sufficiently large that theproperties of the hot phase evolve; thus, density fluctuations in low β , galaxy cooling curve (which is more unstable to fragmentationto multi-phase gas) requires a larger box (to avoid “over-cooling”)for converged results.We close by noting that while magnetic fields clearly enhancethermal instability in a stratified medium, there are also other mech-anisms for doing so. Any process that damps buoyant oscillationcan in principle enhance thermal instability. Thus, for instance, aflattening of the entropy gradient (Voit et al. 2017), which weakensbuoyant restoring forces, and also expel low entropy gas to largeradii (where heating is weak), can potentially have similar effects,though quantifying the impact of these effects requires detailednumerical simulation. In the MHD case, MHD stress oppose thebending of field lines and act as a drag force on an oscillating blob(similar to how the drag force on a moving blob is enhanced bymagnetic fields; McCourt et al. 2015), thus damping buoyant oscil-lations. The key feature we find that make it particularly conclusivefor enhancing thermal instability is that this mechanism becomes important at remarkably low field strengths, and independent offield orientation. This implies that it could potentially be importantthroughout the halo. Furthermore, magnetic fields are ubiquitousand long-lived. By contrast, other processes are less universal andeither spatially or temporally confined. Halos are generally strongstratified by entropy and exhibit flat entropy profiles only in theircores, if at all. Turbulence decays in an eddy turnover time, and thestrong supersonic turbulence needed to overcome buoyancy restor-ing forces is only present in halos during periods of strong AGNor starburst activity. The effects of magnetic fields are much moreuniversal, and we ignore them at our peril. We thank Mark Voit for helpful conversations. We also thank ourreferee, Prateek Sharma, for helpful and insightful comments whichimproved the paper. We acknowledge NASA grants NNX15AK81G,NNX17AK58G and HST-AR-14307.001-A for support. We alsoacknowledge support from the University of California Office of thePresident Multicampus Research Programs and Initiatives throughaward MR-15-328388. This research has used the Extreme Scienceand Engineering Discovery Environment (XSEDE allocations TG-AST140086). We have made use of NASA’s Astrophysics DataSystem and the yt astrophysics analysis software suite (Turk et al.2010). SPO thanks the law offices of May Oh and Wee for hospitalityduring the completion of this paper.
REFERENCES
Anderson M. E., Churazov E., Bregman J. N., 2015, MNRAS, 452, 3905Balbus S. A., 1991, The Astrophysical Journal, 372, 25Balbus S. A., Soker N., 1989, The Astrophysical Journal, 341, 611Bale S., Pulupa M., Salem C., Chen C., Quataert E., 2013, The AstrophysicalJournal Letters, 769, L22Beck R., Wielebinski R., 2013, Magnetic Fields in Galaxies. Springer,Berlin, p. 641, doi:10.1007/978-94-007-5612-0_13Bernet M. L., Miniati F., Lilly S. J., 2013, ApJ, 772, L28Birnboim Y., Dekel A., 2003, MNRAS, 345, 349Borisova E., et al., 2016, preprint, ( arXiv:1605.01422 )Cai Z., et al., 2017, ApJ, 839, 131Cantalupo S., Arrigoni-Battaia F., Prochaska J. X., Hennawi J. F., MadauP., 2014, Nature, 506, 63Choudhury P. P., Sharma P., 2016, Monthly Notices of the Royal Astronom-ical Society, 457, 2554Cox D. P., 2005, ARA&A, 43, 337Dumke M., Krause M., 1998, in Breitschwerdt D., Freyberg M. J., Truem-per J., eds, Lecture Notes in Physics, Berlin Springer Verlag Vol.506, IAU Colloq. 166: The Local Bubble and Beyond. pp 555–558,doi:10.1007/BFb0104782Fang T., Bullock J., Boylan-Kolchin M., 2013, ApJ, 762, 20Fryxell B., et al., 2000, ApJS, 131, 273Gaspari M., Ruszkowski M., Sharma P., 2012, ApJ, 746, 94Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS, 432, 3401Han J., , Annual Review of Astronomy and Astrophysics, pp 111–157Hattori M., Yoshida T., Habe A., 1995, Monthly Notices of the Royal As-tronomical Society, 275, 1195Haverkorn M., 2015, in Lazarian A., de Gouveia Dal Pino E. M., Melioli C.,eds, Astrophysics and Space Science Library Vol. 407, Magnetic Fieldsin Diffuse Media. p. 483 ( arXiv:1406.0283 ), doi:10.1007/978-3-662-44625-6_17Heesen V., Krause M., Beck R., Dettmar R.-J., 2009, A&A, 506, 1123Hennawi J. F., Prochaska J. X., Cantalupo S., Arrigoni-Battaia F., 2015,Science, 348, 779 MNRAS , 1–17 (2017) mpact of B-fields on thermal instability Ji S., Oh S. P., Ruszkowski M., Markevitch M., 2016, Monthly Notices ofthe Royal Astronomical Society, 463, 3989Lee D., 2013, Journal of Computational Physics, 243, 269Li Y., Bryan G. L., 2012, The Astrophysical Journal, 747, 26Li Y., Bryan G. L., 2014, ApJ, 789, 54Li Y., Bryan G. L., Ruszkowski M., Voit G. M., O’Shea B. W., Donahue M.,2015, ApJ, 811, 73Loewenstein M., 1990, The Astrophysical Journal, 349, 471Maller A. H., Bullock J. S., 2004, MNRAS, 355, 694Martin C. L., Dijkstra M., Henry A., Soto K. T., Danforth C. W., Wong J.,2015, ApJ, 803, 6McCourt M., Sharma P., Quataert E., Parrish I. J., 2012, MNRAS, 419, 3319McCourt M., O’Leary R. M., Madigan A.-M., Quataert E., 2015, MonthlyNotices of the Royal Astronomical Society, 449, 2McCourt M., Oh S. P., O’Leary R., Madigan A.-M., 2018, MNRAS, 473,5407Mora S. C., Krause M., 2013, A&A, 560, A42Pizzolato F., Soker N., 2005, ApJ, 632, 821Planck Collaboration et al., 2013, A&A, 557, A52Prasad D., Sharma P., Babul A., 2015, ApJ, 811, 108Ruszkowski M., Oh S. P., 2010, ApJ, 713, 1332Scannapieco E., 2017, ApJ, 837, 28Scannapieco E., Brüggen M., 2015, ApJ, 805, 158Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420, 3174Singh A., Sharma P., 2015, MNRAS, 446, 1895Stein R. F., Leibacher J., 1974, Annual review of astronomy and astrophysics,12, 407Stern J., Hennawi J. F., Prochaska J. X., Werk J. K., 2016, ApJ, 830, 87Stocke J. T., Keeney B. A., Danforth C. W., Shull J. M., Froning C. S., GreenJ. C., Penton S. V., Savage B. D., 2013, ApJ, 763, 148Thompson T. A., Quataert E., Zhang D., Weinberg D. H., 2016, MNRAS,455, 1830Townsend R., 2009, The Astrophysical Journal Supplement Series, 181, 391Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Nor-man M. L., 2010, yt: A Multi-Code Analysis Toolkit for AstrophysicalSimulation Data, Astrophysics Source Code Library (ascl:1011.022)Tzeferacos P., et al., 2012, High Energy Density Physics, 8, 322Voit G. M., Donahue M., 2015, ApJ, 799, L1Voit G. M., Bryan G. L., O’Shea B. W., Donahue M., 2015, ApJ, 808, L30Voit G. M., Meece G., Li Y., O’Shea B. W., Bryan G. L., Donahue M., 2017,ApJ, 845, 80Wagh B., Sharma P., McCourt M., 2014, MNRAS, 439, 2822Werk J. K., Prochaska J. X., Thom C., Tumlinson J., Tripp T. M., O’MearaJ. M., Peeples M. S., 2013, ApJS, 204, 17Werk J. K., et al., 2014, ApJ, 792, 8Wiener J., Oh S. P., Zweibel E. G., 2017, MNRAS, 467, 646Zhang D., Thompson T. A., Quataert E., Murray N., 2017, MNRAS, 468,4801MNRAS000