Mixing and transport of metals by gravitational instability-driven turbulence in galactic discs
Antoine C. Petit, Mark R. Krumholz, Nathan J. Goldbaum, John C. Forbes
MMon. Not. R. Astron. Soc. , 1–11 (2014) Printed 10 November 2018 (MN L A TEX style file v2.2)
Mixing and transport of metals by gravitationalinstability-driven turbulence in galactic discs
Antoine C. Petit, , Mark R. Krumholz, , Nathan J. Goldbaum , & John C. Forbes ICFP undergraduate program, Physics department, École Normale Supérieure, 45 Rue d’Ulm 75005, Paris, France Department of Astronomy and Astrophysics, University of California, 1156 High Street, Santa Cruz, CA 95064, USA
10 November 2018
ABSTRACT
Metal production in galaxies traces star formation, and is highly concentrated towardthe centers of galactic discs. This suggests that galaxies should have inhomogeneousmetal distributions with strong radial gradients, but observations of present-day galax-ies show only shallow gradients with little azimuthal variation, implying the existenceof a redistribution mechanism. We study the role of gravitational instability-driventurbulence as a mixing mechanism by simulating an isolated galactic disc at high res-olution, including metal fields treated as passive scalars. Since any cylindrical field canbe decomposed into a sum of Fourier-Bessel basis functions, we set up initial metalfields characterized by these functions and study how different modes mix. We findboth shear and turbulence contribute to mixing, but the mixing strongly depends onthe symmetries of the mode. Non-axisymmetric modes have decay times smaller thanthe galactic orbital period because shear winds them up to small spatial scales, wherethey are erased by turbulence. The decay timescales for axisymmetric modes are muchgreater, though for all but the largest-scale inhomogeneities the mixing timescale isstill short enough to erase chemical inhomogeneities over cosmological times. Thesedifferent timescales provide an explanation for why galaxies retain metallicity gradi-ents while there is almost no variation at a fixed radius. Moreover, the comparativelylong timescales required for mixing axisymmetric modes may explain the greater di-versity of metallicity gradients observed in high redshift galaxies as compared to localones: these systems have not yet reached equilibrium between metal production anddiffusion.
Understanding the dynamics of the flow of metals throughand around galaxies is a key problem in the study of galaxyformation. Metals trace the history of the gas flows in galax-ies such as the young Milky Way (Ivezić et al. 2012) andrecord the buildup of stellar populations in their progeni-tors at high redshift (Tremonti et al. 2004; Erb et al. 2006).Moreover, metals are not just passive tracers. They changethe chemical and radiative cooling properties of the inter-stellar medium, altering how stars form (Krumholz et al.2009, 2011; Krumholz & Dekel 2012; Krumholz 2012, 2013;Glover & Clark 2012a,b).It is possible to observe the spatial distribution of met-als in the Milky Way (Henry et al. 2010; Balser et al. 2011;Luck & Lambert 2011; Yong et al. 2012), in nearby galaxies(Vila-Costas & Edmunds 1992; Considere et al. 2000; Pi-lyugin et al. 2004; Kennicutt et al. 2011), and even in thehigh redshift universe (Cresci et al. 2010; Jones et al. 2013).The distribution differs from one galaxy to another in thelocal universe, but generally, a radial gradient of the order of − . dex kpc − is observed; the negative sign means thatmetallicity decreases at large radii. High redshift galaxiesare far less regular, and show patterns that vary from flat or slightly positive gradients to some even steeper than in thelocal universe.Metal production in a galaxy is continuously fed bythe stellar feedback (Phillipps & Edmunds 1991). Withoutlarge-scale mixing, the density of metals should be directlyproportional to the stellar density. However it has been ob-served that outer galaxy regions have a metallicity that issignificantly higher than one would expect if the only metalspresent were those produced by stars at the same galactocen-tric radius, while the inner regions of galaxies have metal-licities much smaller than one would expect if all locally-produced metals remained part of the disc (Bresolin et al.2009, 2012; Werk et al. 2011). Moreover, galaxies show gas-phase metallicity inhomogeneities of at most a few tenths ofa dex at fixed galactocentric radius (Rosolowsky & Simon2008; Bresolin 2011; Sanders et al. 2012; Berg et al. 2013),while star formation is far patchier. All of these observationsimply that metal transport must play an important role indetermining the metallicity distribution of galaxies.Unfortunately, metal transport theories remain veryprimitive. The main problem is that it is very difficult torun high resolution simulations of entire galaxies from red-shift z (cid:39) − , when the bulk of the stars in present-daydiscs formed, to the present-day. Simulations that do include c (cid:13) a r X i v : . [ a s t r o - ph . GA ] M a r A. Petit et al. explicit metal tracking over cosmological times are generallyforced by resolution limits to adopt a temperature floor of ∼ K, and are thus unable to provide a realistic treatmentof transport through the multi-phase interstellar medium(ISM) of a galactic disc (e.g. Wiersma et al. 2009; Pilking-ton et al. 2012a,b; Brook et al. 2012; Few et al. 2012a,b;Minchev et al. 2013, 2014). The alternative is parameterized1D models (e.g. Chiappini et al. 2001; Spitoni & Matteucci2011; Forbes et al. 2013), but for these cases the rate ofmetal transport is a free parameter, not a prediction of themodel.In order to study mixing and transport within galacticdiscs, the only solution is to do local studies. A large num-ber of such studies have been published focusing on the ra-dial migration of stars through a galactic disc (e.g., Brunettiet al. 2011; Bird et al. 2012; Di Matteo et al. 2013; Grandet al. 2014, to name only a few), but much less work has beendone focusing on the gas. de Avillez & Mac Low (2002) per-formed early studies of supernova-driven turbulent mixing,and more recently Yang & Krumholz (2012) used shearing-box simulations to show that turbulence driven by ther-mal instability is very efficient at mixing metals. They alsodemonstrated that the multi-phase nature of the ISM hasdramatic consequences for the mixing of metals, implyingthat only simulations with enough resolution to allow sucha multi-phase medium to form are likely to produce reli-able results. However, both de Avillez & Mac Low’s andYang & Krumholz’s studies were limited to small portionsof a galaxy, making it impossible to study transport at thegalactic scale, driven by geometrical effects and galaxy-scaleprocesses like spiral arms. Kubryk et al. (2013) simulatedan entire isolated disc galaxy to measure the effects of abar on radial mixing of gas and stars, and Kubryk et al.(2015) used the results of this simulation to construct asemi-analytic model of element mixing. However, these simu-lations suffered from much the same resolution limitations asthe cosmological ones, in that they used a ∼ K temper-ature floor and thus lacked a multi-phase ISM. The highest-resolution study performed to date is that of Grand et al.(2015), who studied stellar and gas migration in an isolatedgalaxy simulated at ∼ pc resolution, which is still notsufficient to capture the cold phase of the atomic ISM.Our goal in this paper is to study the transport of met-als in a global simulation of a large spiral galaxy, includingthe effects of spiral structure and thermal and gravitationalinstability. To this end, we simulate an isolated galaxy at pc resolution, roughly an order of magnitude higher in res-olution than previous studies, and with a cooling floor thatis low enough to allow development of a full multi-phaseatomic medium. In Section 2 we describe the setup of oursimulations, and how we treat the metal fields. Section 3presents a quantitative analysis of mixing. In Section 4 wediscuss the astrophysical implications of our work. Finally,we summarize in Section 5. To study the turbulent mixing, we simulate an isolatedgalaxy using the Adaptive Mesh Refinement (AMR) code ENZO (Bryan et al. 2014). ENZO includes fluid dynamics,gravity, sink particles, and radiative cooling implementedwith the cooling and chemistry library Grackle . We do notinclude magnetic fields since Yang & Krumholz (2012) haveshown that they have a very a small effect on the turbulentmixing on galactic scales. We model star formation suchthat, when the gas density in any cell reaches a thresholddensity of 50 particles per cubic centimeter, there is a prob-ability equal to (cid:15) SF m g, cell /M min that a part of its gas massis transformed into a star particle of M min that representsa star cluster. Here, (cid:15) SF = 0 . is the star formation effi-ciency, m g, cell is the mass of gas in the considered cell, and M min = 1000 M (cid:12) is the minimum mass of a star particle(Goldbaum et al., 2014, in preparation). This simulation didnot include supernova feedback, both in order to limit thecomputation time and to provide a baseline for the effects ofgravitational instability alone, without the extra turbulenceprovided by supernovae. Simulations including supernovaeare in progress, and will be described in future work.Our initial conditions follow the setup for isolated galax-ies defined by the AGORA project (Kim et al. 2014). Fulldetails on the simulation setup are given in Kim et al. (2014),Springel et al. (2005), and Goldbaum et al. (2015, in prepa-ration), so we simply summarize here the most relevantproperties. We start with a dark matter halo taken fromthe AGORA project, with a circular velocity of km s − .We initialize the baryons as a cylindrical gas cloud with amass M g = 4 . × M (cid:12) in a gaseous halo with the samemass. The gas fraction of the galaxy is 20%. The remainingmass is composed of star particles. The gas density in thegalaxy decreases exponentially both in radius and in height,i.e., the original gas profile is ρ g ( r, z ) = ρ g, exp (cid:18) − rR g (cid:19) exp (cid:18) − | z | h g (cid:19) , (1)where ρ g, = 10 − g cm − is the initial gas density in thecenter of the galaxy and h g = 343 pc and R g = 3 . kpc arethe initial vertical and radial scale length. The gas densityfollows this profile until it reaches the halo. The boundarybetween the halo and the galaxy is determined by the con-dition ρ g T g = ρ h, T h , (2)where ρ h, = 10 − g cm − is the density of the halo, T g =10 K is the initial temperature of the disc, and T h = 10 K is the initial temperature of the halo. The density andtemperature within the halo are uniform. The velocity in thedisc is initially purely circular and fits a flat rotation curveat large radius, giving an orbital period t orb = 175 Myr at . The velocity is equal to zero in the halo initially.We place the galaxy in a computational box formed bya cube root grid with 10 levels of refinement by a factorof 2 per level. The cell size on the finest level is dx = 20 pc, so we are able to marginally resolve the scale heightof the galaxy after the gravitational collapse. The grids arerefined on criteria of gas mass in a cell, particle mass in acell and Jeans length (Truelove et al. 1997). We resolve theJeans length by at least 32 cells on all levels except the finestlevel; on the finest level we add artificial pressure support https://grackle.readthedocs.org/c (cid:13) , 1–11 ixing of metals to ensure that the Jeans length is always resolved by atleast 4 cells. We also refine at the beginning to ensure theaccuracy of the initial conditions (Goldbaum et al., 2014, inpreparation).To create a turbulent thin disc, we first run the sim-ulation for 150 Myr. Over this timescale the gas collapsesinto a thin disc with a height of 200 pc, and a turbulentregion in the center with a roughly 10 kpc radius appears.See Goldbaum et al. (2014, in preparation) for full detailson the evolution of the galaxy during this time. The densitydistribution of the inner region of the galaxy after this 150Myr of evolution is displayed in Figure 1. To study how the turbulence mixes metals in a galaxy, wetake the simulation as it stands after 150 Myr of evolutionand add a set of passive scalar fields, or “colors”. The initialconditions for a realistic metal field depend on the full starformation and merger history of a galaxy, and are obviouslynot available to us for our artificial initial conditions. Evenin a fully cosmological simulation there might be very sig-nificant variations in metal distribution from one galaxy toanother. However, we can overcome this problem by initial-izing the tracer with a very general pattern. If we considera field χ ( r, θ ) defined on a disc that goes to zero (or anyother constant, since the offset does not change the mixing)at the edge, χ can be decomposed in a sum of the productof Bessel and trigonometric functions χ ( r, θ ) = ∞ (cid:88) n =1 ∞ (cid:88) m =0 [ a nm J cnm ( r, θ ) + b nm J snm ( r, θ )] , (3)where J cnm and J snm are the Fourier-Bessel functions definedby J cnm ( r, θ ) = J m ( z nm r/R ) cos( mθ ) (4) J snm ( r, θ ) = J m ( z nm r/R ) sin( mθ ) , (5) a n,m and b n,m are coefficients defined by a nm = (cid:104) χ | J cnm ( r, θ ) (cid:105)|| J cnm ( r, θ ) || (6) b nm = (cid:104) χ | J snm ( r, θ ) (cid:105)|| J snm ( r, θ ) || , (7) R is the radius of the disc we are considering, z nm is the n thpositive zero of J m , and J m is the m th Bessel function ofthe first kind. The inner products appearing in the equationabove are defined by (cid:104) f | g (cid:105) = ( πR ) − (cid:90) fg dA = (cid:90) π (cid:90) R fgπR r dr dθ. (8)Note that the functions are J cnm and J snm are or-thogonal under this inner product, i.e., (cid:104) J cnm | J cn (cid:48) m (cid:48) (cid:105) ∝ δ nn (cid:48) δ mm (cid:48) and similarly for J snm . Furthermore, recallthat J m ( z nm r/R, θ ) cos( mθ ) and J m ( z nm r/R, θ ) sin( mθ ) areeigenvalues for the Laplacian. Thus we would expect that,for a pure diffusion process, the evolution of any initial metalfield χ could be described simply as a decrease in amplitudeof the various Fourier-Bessel modes into which it can be de-composed, with no transfer of power between them. In thisrespect using a Fourier-Bessel decomposition is the naturalextension to cylindrical coordinates of the decomposition of a metal field into Fourier modes for a shearing periodic boxintroduced by Yang & Krumholz (2012).Given this analysis, to study mixing we take the stateof the simulation at 150 Myr and add 25 different color fieldswith the following initial conditions: χ nm ( r, θ, z, t = 0) = J m ( z nm r/R ) cos( mθ ) = J cnm ( r, θ ) , (9)where R = 8 kpc is the radius of the turbulent area, n = 1 . . . and m = 0 . . . , and ( r, θ, z ) ∈ [0 , R ] × [0 , π ] × [ − h g , h g ] ; χ nm is set to zero outside of this cylinder. We onlyconsider the cosine terms because the fields represented bythe sine terms are identical up to a 90 ◦ rotation about the z axis. Since the problem is, in a statistical sense, invariantunder such a transformation, the modes represented by thesine terms should not mix any differently than the ones rep-resented by the cosine terms, and thus there is no additionalinformation gained by including them.Once they are added, the color fields are advected by thegas motion as passive scalars. Numerically, we handle advec-tion of metal tracers in exactly the same manner as advec-tion of mass. Thus if the hydrodynamic solver returns a massflux F M across a particular cell face, the flux of “metal mass"in field n, m is simply χ nm, upwind F M , where χ nm, upwind is theupwind value of the metal field χ nm , determined using thesame PPM interpolation scheme as for the hydrodynamicquantities. This scheme is conservative, in the sense thathydrodynamic evolution leaves the total “metal mass" inthe simulation domain, (cid:82) ρχ nm dV , unchanged, even as thecolor field is advected between computational cells. Whenstars form, we leave the color concentration χ nm in the star-forming cell unchanged, which amounts to assuming thatthe metal concentration in the stars is identical to that ofthe gas cells in which they form. In principle we could thentrack the metal contents of the stars (e.g., Feng & Krumholz2014), but we have not done so here, because our simulationdoes not run long enough for significant stellar migration tooccur, and it is not clear if our resolution of the N -bodydynamics of the stars is sufficient to follow stellar migrationaccurately in any event.We pause to make two final points regarding our numer-ical method. Some previous simulations of metal transport inisolated galaxies have used smoothed particle hydrodynam-ics (SPH) methods (e.g. Kubryk et al. 2013; Grand et al.2015). These codes have an advantage over our method inthat, because they allow the authors to track individual par-ticles, it is possible to distinguish between “radial migration”of gas, i.e., bulk radial transport of gas, and “mixing”, i.e.,homogenization of the metal distribution but without bulkradial redistribution of mass. We cannot make this distinc-tion, and thus we will use the generic term mixing to referto any process that homogenizes the metal distribution.On the other hand, our Eulerian method also offers asignificant advantage over SPH when it comes to simula-tion of metal transport. SPH codes require a prescriptionfor sub-grid scale turbulent mixing in order to follow metaltransport, because without such a model SPH artificiallysuppresses mixing below the resolution scale (e.g., Wadsleyet al. 2008; Greif et al. 2009). These sub-grid models con-tain a number of free parameters, and their accuracy hasnot been extensively calibrated. In contrast, our simulations c (cid:13) , 1–11 A. Petit et al.
10 5 0 5 10 x (kpc) y ( k p c ) -5 -4 -3 -2 -1 D e n s i t y ( g c m ) Figure 1.
Projection of the density of the inner turbulent region of the galaxy. mix naturally at the grid scale, and do not require us to usean explicit subgrid mixing model.
After inserting the color fields as described above, we allowthe simulation to evolve for another 100 Myr. During thistime each of the passive scalar fields χ nm evolves, and inthis section we discuss the nature of this time evolution. Wehave available the value of each scalar field as a function oftime and position, χ nm ( r, θ, z, t ) , but in order to simplifythe analysis we neglect their vertical structure and only an-alyze their column density-weighted averages. Specifically,we define χ nm ( r, θ, t ) = (cid:82) χ nm ( r, θ, z, t ) ρ ( r, θ, z, t ) dz (cid:82) ρ ( r, θ, z, t ) dz , (10)where ρ is the gas volume density, and from now on whenwe refer to χ nm we mean the vertically-integrated quantity.We perform the integration, and all other analysis presentedin this paper, using the analysis tool yt (Turk et al. 2011).Since our simulations do not include supernova feedback,the disk remains thin and there is no gas expulsion on thevertical direction. For this reason, transport of tracers along the vertical direction is negligible and the integration over z does not significantly alter the results.We quantify the evolution of the color fields by com-puting the projection of χ nm ( r, θ, t ) on the different Fourier-Bessel basis functions. The projection is done on a disc ofradius R = 8 kpc with the inner product defined on equation8. We define these projections by P cnm,n (cid:48) m (cid:48) ( t ) = (cid:104) χ nm ( t ) | J cn (cid:48) m (cid:48) (cid:105) (11) P snm,n (cid:48) m (cid:48) ( t ) = (cid:104) χ nm ( t ) | J sn (cid:48) m (cid:48) (cid:105) . (12)The quantity P s,cnm,n (cid:48) m (cid:48) ( t ) indicates how much of the initialpower from the mode n, m has been transferred to the mode n (cid:48) , m (cid:48) (or, for n = n (cid:48) and m = m (cid:48) , how much of the origi-nal power remains) at time t . Note that, even if the metalfield remained in a fixed pattern and the galaxy rotated as asolid body, this rotation would exchange power between the J cnm and J snm . To remove this effect, rather than consider-ing P cnm,n (cid:48) ,m (cid:48) and P snm,n (cid:48) m (cid:48) separately, we instead compute P nm,n (cid:48) m (cid:48) , where P nm,n (cid:48) m (cid:48) = (cid:0) P cnm,n (cid:48) m (cid:48) (cid:1) + (cid:0) P snm,n (cid:48) m (cid:48) (cid:1) . (13)This quantity is invariant under rotation, so if the galaxyrotated as a solid body and the spread of the metal fieldwere described simply by a diffusion equation in the rotat-ing frame, then we would have P nm,n (cid:48) m (cid:48) ( t ) = δ nn (cid:48) δ mm (cid:48) f ( t ) ,with f (0) = 1 and f ( t ) strictly decreasing with t for t > . c (cid:13) , 1–11 ixing of metals P n , n ( t ) n = 1n = 2n = 3n = 4n = 5 Figure 3.
Evolution over time of P n ,n ( t ) , the fraction of theoriginal power remaining in each of the m = 0 modes (equation13) for n varying from 1 to 5. In reality we shall see that this is not how the color fieldsevolve, but this case nonetheless provides a useful zerothorder model against which our findings can be compared.
A realistic metallicity field will contain contributions fromtwo different types of modes: the axisymmetric ones, i.e.,those with m = 0 , carry information on the average radialdependence, and the non-axisymmetric ones, i.e., the m (cid:54) = 0 modes, which describe spiral structure or other deviationsfrom uniformity at fixed radius. Since 2-armed spirals are themost common type of non-axisymmetric pattern in galaxies,we will focus on the case m = 2 . We discuss the m = 0 modesin this Section, and the m = 2 ones in the subsequent one.Figure 2 shows the vertically-integrated color fields χ nm for m = 0 , n = 1 , at t = t (the time when the color fieldsare first added) and t = t + 100 Myr. As a consequence ofthe choice of R , the gas in which the tracers are deposited isfully turbulent. The original pattern is slowly destroyed bythe turbulence, and the destruction looks faster for n = 2 than for n = 1 .Confirming this visual impression, Figure 3 shows thedecrease of P nm,nm ( t ) , the power remaining in the originalmode, with time for the axisymmetric modes n ∈ { . . . } .We observe that the mode m = 0 , n = 1 remains al-most completely stationary for the time scale we consider,while the other modes slightly decrease. The rate at whichthe modes decrease is inversely correlated with n , which isnot surprising: higher n modes correspond to smaller spa-tial scales, and we expect that turbulence should mix outsmaller-scale inhomogeneities faster than larger-scale ones.After 100 Myr, the mode that has decreased the most, n = 5 has lost roughly 70% of its original power, but it is still in alinear phase.We can next investigate how power is transferredbetween modes, which is described by the quantity P nm,n (cid:48) m (cid:48) (t). In Figure 4, we show P ,n (cid:48) m (cid:48) ( t ) at t = t +100 Myr. Physically, this quantity shows how the power that n m l og [ P n m ] Figure 4.
Power transferred from n = 2 , m = 0 to all modes, P ,n (cid:48) m (cid:48) , for n (cid:48) varying from 1 to 15 and m (cid:48) from 0 to 10 after100 Myr of mixing was originally in the n = 2 , m = 0 mode has been transferredto other modes over 100 Myr of evolution. We can see fromthe figure that turbulent diffusion transfers power to bothhigher n and higher m , and that it does so approximatelyisotropically. This figure is also consistent with the resultshown in Figure 3 that the mixing is slow for the n = 2 , m = 0 mode. The majority of the power remains in theoriginal mode, and only ∼ percent of it has spread toother modes. Recall, however, that, at the time illustrated inFigure 4, the simulation has evolved for only 100 Myr, whichis less than a full orbital period of the disc (at 8 kpc). It istherefore not surprising that only a relatively small fractionof the total power has been transferred out of the mode. Star formation that is concentrated in spiral arms is likely toproduce non-axisymmetry in the metal field, and we there-fore next consider non-axisymmetric modes. Which modedominates in a real galaxy will likely depend on whetherthe spiral arms are grand design or flocculent. For this ex-ample we choose to focus on m = 2 , the mode that shoulddominate for a two-armed spiral. We present results for allother m modes below. In Figure 5, we show some snapshotsof the modes m = 2 , n = 1 , at t and t + 100 Myr. Wecan see that the mixing looks more efficient than in the ax-isymmetric case. The pattern has been destroyed both bythe turbulence and the differential rotation. We will studyin this part both of these mechanisms.In Figure 6, we can see that, after a brief transient, P nm,nm ( t ) decreases linearly with time before oscillating be-tween 0 and 0.2 of the original value for the higher valuesof n . The modes vanish much faster than the axisymmetricones. However, most of the reduction in power is a resultof transfer to other modes rather than outright destruction.Indeed, if we compute the sum of the power remaining overall the modes, (cid:80) n (cid:48) ,m (cid:48) P nm,n (cid:48) m (cid:48) ( t ) , which is equivalent tocomputing the norm in the real space || χ nm ( t ) || , we findthat it decreases by at most of order over the 100 Myr c (cid:13)000
Power transferred from n = 2 , m = 0 to all modes, P ,n (cid:48) m (cid:48) , for n (cid:48) varying from 1 to 15 and m (cid:48) from 0 to 10 after100 Myr of mixing was originally in the n = 2 , m = 0 mode has been transferredto other modes over 100 Myr of evolution. We can see fromthe figure that turbulent diffusion transfers power to bothhigher n and higher m , and that it does so approximatelyisotropically. This figure is also consistent with the resultshown in Figure 3 that the mixing is slow for the n = 2 , m = 0 mode. The majority of the power remains in theoriginal mode, and only ∼ percent of it has spread toother modes. Recall, however, that, at the time illustrated inFigure 4, the simulation has evolved for only 100 Myr, whichis less than a full orbital period of the disc (at 8 kpc). It istherefore not surprising that only a relatively small fractionof the total power has been transferred out of the mode. Star formation that is concentrated in spiral arms is likely toproduce non-axisymmetry in the metal field, and we there-fore next consider non-axisymmetric modes. Which modedominates in a real galaxy will likely depend on whetherthe spiral arms are grand design or flocculent. For this ex-ample we choose to focus on m = 2 , the mode that shoulddominate for a two-armed spiral. We present results for allother m modes below. In Figure 5, we show some snapshotsof the modes m = 2 , n = 1 , at t and t + 100 Myr. Wecan see that the mixing looks more efficient than in the ax-isymmetric case. The pattern has been destroyed both bythe turbulence and the differential rotation. We will studyin this part both of these mechanisms.In Figure 6, we can see that, after a brief transient, P nm,nm ( t ) decreases linearly with time before oscillating be-tween 0 and 0.2 of the original value for the higher valuesof n . The modes vanish much faster than the axisymmetricones. However, most of the reduction in power is a resultof transfer to other modes rather than outright destruction.Indeed, if we compute the sum of the power remaining overall the modes, (cid:80) n (cid:48) ,m (cid:48) P nm,n (cid:48) m (cid:48) ( t ) , which is equivalent tocomputing the norm in the real space || χ nm ( t ) || , we findthat it decreases by at most of order over the 100 Myr c (cid:13)000 , 1–11 A. Petit et al. y [ kp c ] y [ kp c ] χ n ( t ) Figure 2.
Snapshots at t (left) and t = t + 100 Myr Myr (right) of the m = 0 , n = 1 (top) and m = 0 , n = 2 (bottom) color fields χ nm . of evolution. Figure 7, which shows P ,n (cid:48) m (cid:48) evaluated at t = t + 100 Myr, reveals where the power that is removedfrom the n = 1 , m = 2 mode is transferred. We can seethat most of the power is still in modes with m = 2 , while asmaller fraction has leaked into m (cid:54) = 2 modes (about 20%). The main cause of the spread of the power in the m = 2 lineshown in Figure 7 is the shear induced by the differentialrotation. We illustrate this point in Figure 8, which showsthe real-space reconstruction of the color field that resultsfrom summing up only the m = 2 modes in Figure 7, i.e., χ sh n ( t ) = (cid:88) n (cid:48) (cid:20) P cn ,n (cid:48) ( t ) || J cn (cid:48) || J cn (cid:48) + P sn ,n (cid:48) ( t ) || J sn (cid:48) || J sn (cid:48) (cid:21) . (14)As is clear from Figure 8, this procedure picks out thesheared pattern created by the differential rotation. As first We can estimate the maximum values of n and m accessible toour simulation from our resolution of 20 pc. Since the disc radiusis 8 kpc, we have roughly 400 cells per disc radius. If we assumethat we need ∼ cells per wavelength to resolve the mode, thenwe should be able to obtain reliable estimates up to n and m values of ∼ . P n , n ( t ) n = 1n = 2n = 3n = 4n = 5 Figure 6.
Same as Figure 3, but for the m = 2 rather than m = 0 modes. pointed out by Yang & Krumholz (2012), this a very im-portant effect for diffusion. As we have already seen whenconsidering axisymmetric modes, turbulence is much moreefficient at mixing away inhomogeneities on small physical c (cid:13) , 1–11 ixing of metals y [ kp c ] y [ kp c ] χ n ( t ) Figure 5.
Same as Figure 2, but now for the m = 2 , n = 1 (top) and m = 2 , n = 2 (bottom) color fields. n m l og [ P n m ] Figure 7.
Same as Figure 4, but for the n = 1 , m = 2 modes. scales (higher n and m ) than on large ones. Since the pri-mary effect of differential rotation is to transport power tohigher m , even if the shear by itself doesn’t mix efficiently,it increases significantly the speed of the diffusion becausethe radial pattern is not preserved by differential rotation. y [ kp c ] χ s h ( t ) Figure 8.
Projection of χ ( t ) after 100 Myr of mixing on themodes m = 2 , n = 1 .. using P c,s ,n (cid:48) m (cid:48) as defined in equation 14.It reproduces the shear induced by the differential rotation of thegalaxy.c (cid:13)000
Projection of χ ( t ) after 100 Myr of mixing on themodes m = 2 , n = 1 .. using P c,s ,n (cid:48) m (cid:48) as defined in equation 14.It reproduces the shear induced by the differential rotation of thegalaxy.c (cid:13)000 , 1–11 A. Petit et al. y [ kp c ] ( χ ( t ) − χ s h ( t )) Figure 9.
Difference between the full color field χ ( t =t + 100 Myr) at time t + 100 Myr , as shown in the top rightpanel of the figure 5, and the partial reconstruction χ sh12 definedby equation 14 and shown in Figure 8. Intuitively, this field showsthe change in the color field due to turbulence rather than shear. While the transport of power along the m = 2 row is causedby differential rotation, the power driven to other modesis a manifestation of the turbulence. To demonstrate this,in Figure 9 we show the χ ( t ) − χ sh12 ( t ) after 100 Myr. Thisquantity is simply the residual that results when we subtractFigure 8 from the upper right panel of Figure 5. We can usethe power associated with this residual field to quantify theimportance of turbulent mixing.Let us consider the quantity T nm = 1 − (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) (cid:80) n (cid:48) P n (cid:48) m,n (cid:48) m (cid:80) ( n (cid:48) ,m (cid:48) ) P nm,n (cid:48) m (cid:48) . (15)Intuitively, the numerator in the fraction is the total powerin modes with the same m as the original one, and thusrepresents the fraction of the original power that is in thesheared field. The denominator is simply the total powersummed over all modes, and the ratio is therefore the frac-tion of the original power that remains in a coherent, shear-distorted pattern. One minus this quantity is the fraction ofthe original power that has been transported or destroyedby turbulence. Thus one may intuitively think of the quan-tity T nm as the fraction of the original power that has beendiffused by the turbulence. We show T nm for m = 2 in Figure 10. In the case of m = 2 , n = 5 , which has the highest dissipation rate for While equation (15) formally involves a sum over all n (cid:48) and m (cid:48) up to infinity, in practice we must of course terminate thecomputation at finite values of n (cid:48) and m (cid:48) . We choose to stop at n (cid:48) = 25 and m (cid:48) = 10 because the sum seems to converge withthese limits, as shown in the Appendix. T n ( t ) n = 1n = 2n = 3n = 4n = 5 Figure 10.
The power transferred from m = 2 to m (cid:54) = 2 modesby turbulence, T nm (equation 15), as a function of the time for m = 2 and n = 1 . . . . the modes we are considering in this part, roughly 10% ofthe power is transferred to m (cid:54) = 2 modes. However, thisis a lower limit on the power dissipated by the turbulence,since turbulence as well as shear can transfer power between m = 2 modes.We can schematically summarize the joint effects of tur-bulence and shear in Figure 11. Turbulence is the only agentcapable of truly mixing the disc and wiping out inhomo-geneities. However, it operates fairly slowly for patters withlarge spatial scales. In comparison, for non-axisymmetricmodes, m (cid:54) = 0 , turbulence is greatly aided by shear. Sheartransports power out of low n modes and into higher n ones; physically, differential rotation transforms any non-axisymmetric pattern into a tightly-wound spiral on atimescale comparable to the orbital period. This works inconjunction with the turbulence to greatly increase the rateof mixing, by moving power out of low n modes where it ishard to dissipate and into high n modes where dissipationis more rapid (Yang & Krumholz 2012). We saw in the previous two sections that some modes canbe very stable whereas others can be completely diffused ontimescales well under an orbital period. We now make thisanalysis quantitative by associating to the destruction ofeach mode a timescale. Figures 3 and 6 show that the powerremaining in the original mode, P nm,nm , is roughly constantduring an initial transient, then undergoes a phase where itdecays roughly linearly with time, and finally stabilizes andbegins oscillating once the majority of the power is gone.This oscillation is caused by the fact that the original modeis no longer dominant and so the power from other modescan feed back into the original one. Examining comparableplots for other values of m , we find that these three phasesare generic. It therefore seems reasonable to estimate a decaytimescale for each mode n, m by computing the slope D mn in the linear phase. D nm increases with n and m (Figures3 and 6), we can also notice that the slope is very small for c (cid:13) , 1–11 ixing of metals nm Effect of shear Effect of the turbulence
Figure 11.
Schematic representation of the actions of the shearand the turbulence on the spectrum of the metal field. n M i x i ng T i m e s ca l e [ y r] m = 0 m = 1 m = 2 m = 3 m = 4 10 -1 M i x i ng T i m e s ca l e [ G a l ac ti c p e r i od ] Figure 12.
The timescale for metallicity inhomogeneities to de-cay, τ nm (equation16), as a function of n for m = 0 . . . . Onthe right side, we normalize the timescale to the galactic orbitalperiod, t orb = 175 Myr at . τ is not plotted because thesimulation time was too small to reach the linear phase for thismode. the m = 0 modes. The timescale is simply the inverse of thisslope, i.e., τ nm = D − nm . (16)Figure 12 shows τ nm versus n and m . We see that forthe m = 0 modes these times are bigger than or similar tothe galactic orbital period t orb = 175 Myr at , whereasfor the non-axisymmetric modes, all are smaller than thegalactic orbital period, sometimes by an order of magnitude.
This analysis of the destruction timescale for metallicity in-homogeneities has strong astrophysical implications. First,we shown that a non-axisymmetric pattern in the metallicityfield is smoothed in less than a galactic orbital period. Thisimplies that non-axisymmetries in the metal field driven byspiral patterns, bar, or similar phenomena in the star for-mation distribution will be suppressed very rapidly. The un-derlying physical mechanism driving this is that illustratedin Figure 11 and discussed by Yang & Krumholz (2012): forany non-axisymmetric mode, differential rotation winds themetal pattern up into a tight spiral on orbital timescales, andthis small-scale pattern is then easily destroyed by turbu-lence. This mechanism likely explains why observed galaxymetallicity variations in the gas at fixed galactocentric ra-dius are so small.However, the axisymmetric patterns seems more stable,particularly for n (cid:54) where τ n is bigger than t orb4 . Themuch larger diffusion timescale for the axisymmetric modeslikely explains why radial metallicity gradients in galaxiespersist even as azimuthal ones are wiped out. Moreover, thediffusion time is also related to the timescale required forthe metal distribution in a galaxy to reach equilibrium be-tween star formation, which drives the metal distributionaway from homogeneity, and turbulent mixing, which drivesit toward homogeneity. Since the low n axisymmetric modesmay dominate the overall metallicity distribution, the re-laxation time required for the metallicity to reach a steadystate would be bigger than several orbits. This might explainwhy metal gradients for high redshift galaxies are far morevaried than those in nearby galaxies. Most z (cid:39) galaxies, in-cluding the Milky Way, experienced their last major mergerbetween z = 1 − , and thus have been in their present con-figuration for ∼ − Gyr. This is a timescale much largerthan their orbital period, suggesting that most present-daygalaxies have had time to reach equilibrium between metalproduction and diffusion. On the other hand, high redshiftgalaxies are often at most a few orbital periods old and aretherefore still in a non-equilibrium state. As a result, theirmetal gradients have not undergone significant smoothing,and instead reflect the patchy distributions of star formationwithin them without much smoothing.Finally, we end this discussion with a caution. We haveargued that, because at present only isolated galaxy sim-ulations can reach the resolutions requires to capture themulti-phase ISM, such simulations are the only way to de-rive realistic rates of metal transport. However, the limi-tation of this approach is that we have run a closed boxsimulation without gas infall or external perturbations. Ina more realistic cosmological environment, the quantitativemixing timescales we have derived might be modified by theprocesses we have been forced to omit. However the differ-ences between axisymmetric and non-axisymmetric shouldbe preserved even in a cosmological context, since the maineffect is purely geometrical. τ is at least bigger than τ and probably bigger than severalgalactic orbital periods.c (cid:13) , 1–11 A. Petit et al.
In this work we simulate the diffusion of inhomogenous metaldistributions in a galactic disc as a result of gravitationalinstability-driven turbulence. To make our study as generalas possible, we note that, in cylindrical geometry, the metalfield can always be decomposed into Fourier-Bessel func-tions. We therefore study how different Fourier-Bessel modesdecay and mix as a result of shear and turbulence. We findthat the efficiency of mixing strongly depends on whetherone is considering an axisymmetric or a non-axisymmetricinhomogeneity. In the former case, the metal field is verystable, destruction of the original pattern requires at leastseveral orbital periods for large-scale modes, and is causedonly by turbulence. In the latter case, the original patternvanishes in less than a galactic orbital period. This differencein timescale is due to the effects of shear. Shear acceleratesdiffusion by winding up the inhomogeneities into tight spi-rals on small spatial scales, effectively transporting powerfrom large to small scales in an orbital period. Once thistransport is complete, turbulence is then able to diffuse thesmall-scale power quite rapidly.This difference between modes in terms of dissipationtime has strong implications for our understanding of theobserved distributions of metals in galactic discs. In partic-ular, the far greater rapidity of mixing for non-axisymmetricmodes than for axisymmetric ones helps explain why galax-ies show consistent radial gradients in metallicity but little tono variation at fixed radius. The long mixing timescales wefind for radial modes also suggest that metal distributions inhigh-redshift galaxies are most likely not yet in equilibriumbetween metal production and mixing. This provides a likelyexplanation for the much greater diversity of metallicity dis-tributions seen at high redshift: these are more reflective ofthe patchy distribution of star formation in these galaxies,while the comparatively uniform behavior of low- z galaxiesarises from the balance between metal production and dif-fusion. Indeed the most important modes are not only thosewith the longest dissipation timescales, but also the mostfed ones by stellar metal production. We leave considerationof that problem to future work. ACKNOWLEDGEMENTS
This work was supported by NSF grants AST-0955300 andAST-1405962, NASA ATP grant NNX13AB84G, and NASATCAN grant NNX14AB52G (MRK, NJG, and JCF). Thesimulations for this research were carried out on the UCSCsupercomputer Hyades, which is supported by the NSF(award number AST-1229745).
APPENDIX A: CONVERGENCE
For computational reasons, we must truncate the Fourier-Bessel expansions we use in our analysis at finite values of n (cid:48) and m (cid:48) . We select as our limits n (cid:48) = 25 and m (cid:48) = 10 , andin this Appendix we show that these limits are high enoughto ensure convergence when we reconstruct fields such as χ sh and quantities like T nm that depend on them. In FigureA1, we plot the quantity T computed using expansionstruncated at various values for n and m . As shown in the Time [Myr] T ( t ) N n = 5, N m = 4 N n = 7, N m = 5 N n = 9, N m = 6 N n = 11, N m = 7 N n = 13, N m = 8 N n = 15, N m = 9 N n = 17, N m = 10 N n = 19, N m = 10 N n = 21, N m = 10 N n = 23, N m = 10 N n = 25, N m = 10 Figure A1. T , computed using a sum truncated at the indi-cated values of n (cid:48) and m (cid:48) . Figure, this quantity is very well converged by the time wereach n (cid:48) = 25 , m (cid:48) = 10 .Another mean to check the convergence is to computethe difference between the norm || χ nm || and the sum overthe first terms of P nm,n (cid:48) m (cid:48) . Indeed, we have the equality: (cid:90) χ d S = ∞ (cid:88) m =0 ,n =1 a nm || J cnm || + b nm || J snm || . (A1)We compute the ratio of the sum over n = 25 and m =10 and || χ || for the 25 color fields studied after 100 Myr. Itappears that the fields lost less than 2% for n = 1 , m = 0 ,and less than 20% for the highest n and m studied. REFERENCES
Balser D. S., Rood R. T., Bania T. M., Anderson L. D.,2011, Astrophys. J., 738, 27Berg D. A., Skillman E. D., Garnett D. R., Croxall K. V.,Marble A. R., Smith J. D., Gordon K., Kennicutt Jr.R. C., 2013, ApJ, 775, 128Bird J. C., Kazantzidis S., Weinberg D. H., 2012, MNRAS,420, 913Bresolin F., 2011, ApJ, 730, 129Bresolin F., Ryan-Weber E., Kennicutt R. C., Goddard Q.,2009, Astrophys. J., 695, 580Bresolin F., Kennicutt R. C., Ryan-Weber E., 2012, Astro-phys. J., 750, 122Brook C. B., et al., 2012, MNRAS, 426, 690Brunetti M., Chiappini C., Pfenniger D., 2011, A&A, 534,A75Bryan G. L., et al., 2014, Astrophys. J. Suppl. Ser., 211, 19Chiappini C., Matteucci F., Romano D., 2001, ApJ, 554,1044Considere S., Coziol R., Contini T., Davoust E., 2000, As-tron. Astrophys.Cresci G., Mannucci F., Maiolino R., Marconi A., GnerucciA., Magrini L., 2010, Nature, 467, 811Di Matteo P., Haywood M., Combes F., Semelin B., SnaithO. N., 2013, A&A, 553, A102 c (cid:13) , 1–11 ixing of metals Erb D. K., Shapley A. E., Pettini M., Steidel C. C., ReddyN. A., Adelberger K. L., 2006, Astrophys. J., 644, 813Feng Y., Krumholz M. R., 2014, Nature, 513, 523Few C. G., Courty S., Gibson B. K., Kawata D., Calura F.,Teyssier R., 2012a, MNRAS, 424, L11Few C. G., Gibson B. K., Courty S., Michel-Dansac L.,Brook C. B., Stinson G. S., 2012b, A&A, 547, A63Forbes J. C., Krumholz M. R., Burkert A., Dekel A., 2013,eprint arXiv:1311.1509, p. 20Glover S. C. O., Clark P. C., 2012a, MNRAS, 421, 9Glover S. C. O., Clark P. C., 2012b, MNRAS, 426, 377Grand R. J. J., Kawata D., Cropper M., 2014, MNRAS,439, 623Grand R. J. J., Kawata D., Cropper M., 2015, MNRAS,447, 4018Greif T. H., Glover S. C. O., Bromm V., Klessen R. S.,2009, MNRAS, 392, 1381Henry R. B. C., Kwitter K. B., Jaskot A. E., Balick B.,Morrison M. A., Milingo J. B., 2010, Astrophys. J., 724,748Ivezić v., Beers T. C., Jurić M., 2012, Annu. Rev. Astron.Astrophys., 50, 251Jones T., Ellis R. S., Richard J., Jullo E., 2013, Astrophys.J., 765, 48Kennicutt R. C., et al., 2011, Publ. Astron. Soc. Pacific,123, 1347Kim J.-h., et al., 2014, Astrophys. J. Suppl. Ser., 210, 14Krumholz M. R., 2012, ApJ, 759, 9Krumholz M. R., 2013, MNRAS, 436, 2747Krumholz M. R., Dekel A., 2012, ApJ, 753, 16Krumholz M. R., McKee C. F., Tumlinson J., 2009, Astro-phys. J., 699, 850Krumholz M. R., Leroy A. K., McKee C. F., 2011, Astro-phys. J., 731, 25Kubryk M., Prantzos N., Athanassoula E., 2013, MNRAS,436, 1479Kubryk M., Prantzos N., Athanassoula E., 2015, A&A,Luck R. E., Lambert D. L., 2011, Astron. J., 142, 136Minchev I., Chiappini C., Martig M., 2013, A&A,Minchev I., Chiappini C., Martig M., 2014, A&A,Phillipps S., Edmunds M. G., 1991, Mon. Not. R. Astron.Soc. (ISSN 0035-8711), 251, 84Pilkington K., et al., 2012a, MNRAS, 425, 969Pilkington K., et al., 2012b, A&A, 540, A56Pilyugin L. S., Vílchez J. M., Contini T., 2004, Astron.Astrophys., 425, 849Rosolowsky E., Simon J. D., 2008, ApJ, 675, 1213Sanders N. E., Caldwell N., McDowell J., Harding P., 2012,ApJ, 758, 133Spitoni E., Matteucci F., 2011, A&A, 531, A72Springel V., Di Matteo T., Hernquist L., 2005, MNRAS,p. 623Tremonti C. A., et al., 2004, Astrophys. J., 613, 898Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H.,Howell L. H., Greenough J. A., 1997, ApJ, 489, L179Turk M. J., Smith B. D., Oishi J. S., Skory S., SkillmanS. W., Abel T., Norman M. L., 2011, ApJS, 192, 9Vila-Costas M. B., Edmunds M. G., 1992, Mon. Not. R.Astron. Soc. (ISSN 0035-8711), 259, 121Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008,MNRAS, 387, 427Werk J. K., Putman M. E., Meurer G. R., Santiago- Figueroa N., 2011, Astrophys. J., 735, 71Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C.,Tornatore L., 2009, MNRAS, 399, 574Yang C.-C., Krumholz M., 2012, Astrophys. J., 758, 48Yong D., Carney B. W., Friel E. D., 2012, Astron. J., 144,95de Avillez M. A., Mac Low M.-M., 2002, ApJ, 581, 1047 c (cid:13)000