The Parker Instability in Disk Galaxies
Luiz Felippe S. Rodrigues, Graeme R. Sarson, Anvar Shukurov, Paul J. Bushby, Andrew Fletcher
AA CCEPTED FOR PUBLICATION IN THE A STROPHYSICAL J OURNAL
Preprint typeset using L A TEX style emulateapj v. 05/12/14
THE PARKER INSTABILITY IN DISK GALAXIES
L. F. S. R
ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , AND
A. F
LETCHER School of Mathematics and Statistics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UK
Accepted for publication in the Astrophysical Journal
ABSTRACTWe examine the evolution of the Parker instability in galactic disks using 3D numerical simulations. We con-sider a local Cartesian box section of a galactic disk, where gas, magnetic fields and cosmic rays are all initiallyin a magnetohydrostatic equilibrium. This is done for different choices of initial cosmic ray density and mag-netic field. The growth rates and characteristic scales obtained from the models, as well as their dependenceson the density of cosmic rays and magnetic fields, are in broad agreement with previous (linearized, ideal)analytical work. However, this non-ideal instability develops a multi-modal 3D structure, which cannot bequantitatively predicted from the earlier linearized studies. This 3D signature of the instability will be of im-portance in interpreting observations. As a preliminary step towards such interpretations, we calculate syntheticpolarized intensity and Faraday rotation measure maps, and the associated structure functions of the latter, fromour simulations; these suggest that the correlation scales inferred from rotation measure maps are a possibleprobe for the cosmic ray content of a given galaxy. Our calculations highlight the importance of cosmic rays inthese measures, making them an essential ingredient of realistic models of the interstellar medium.
Keywords: instabilities, galaxies: ISM, ISM: cosmic rays, ISM: magnetic fields INTRODUCTIONThe magnetic buoyancy (or magnetic Rayleigh–Taylor) in-stability enhanced by cosmic rays in disk galaxies is knownas the Parker instability (Parker 1966, 1967, 1969). Since cos-mic rays are practically weightless but exert pressure compa-rable to that of thermal gas, turbulence and magnetic fields,their effect on the instability is significant. The growth rateof the most unstable mode is of the order of the Alfv´en cross-ing time, τ A = h / c A (cid:39) yr, based on the disk scale height h =
100 pc, and assuming an Alfv´en speed of c A =
10 km s − ;the corresponding wavelength is of order 1 kpc (Parker 1979).The nature of the instability is quite generic and it is there-fore expected to occur in all magnetized astrophysical disks,including accretion disks and spiral galaxies. The saturationmechanism of the instability and the resulting state of the sys-tem are not well understood. As discussed by Foglizzo &Tagger (1994, 1995), differential rotation weakens the insta-bility. Gas viscosity, random magnetic fields and finite cosmicray diffusion also suppress it (Kuznetsov & Ptuskin 1983), asthey interfere with the sliding of the gas to the bottom of mag-netic loops as well as the streaming of cosmic rays to the topof these loops.There have been several studies of the Parker instability us-ing linear perturbation theory. Giz & Shu (1993) extended theinitial works by Parker by employing a more realistic verticalgravitational acceleration, g z ( z ) ∝ tanh z , which is continuousacross the Galactic midplane, rather than the step-function | g z ( z ) | = constant that was previously used. This allowedthem to consider the vertical symmetry of the instability andidentify a family of unstable modes that cross the midplane(in which case the vertical velocity, U z , and the vertical mag-netic field, B z , are not identically zero at z = [email protected] (LFSR),[email protected] (GRS),[email protected] (AS),[email protected] (PJB),andrew.fl[email protected] (AF). this type, which are often referred to as ‘antisymmetric’, werealso found by Horiuchi et al. (1988), although they used a dif-ferent gravitational acceleration, g ( z ) ∝ z / ( r + z ) / . (Theissue of symmetry about the midplane is discussed in detailin section 3.2.) Kim & Hong (1998) showed that the use ofa more realistic gravity profile also results in a faster growthrate and shorter horizontal wavelength for the instability. Ryuet al. (2003) included the effects of finite cosmic ray diffu-sion; the previous works had taken the diffusivity to be infinitealong the magnetic field and zero perpendicular to it. Theirlinear analysis found that the finite, but large, parallel diffu-sivity of cosmic rays in the interstellar medium (ISM) is wellapproximated by infinite-diffusivity models, and that the per-pendicular diffusivity has no significant effect. Kuwabara &Ko (2006) studied Parker-Jeans instabilities in the ISM usinga model which does not have an external gravitational field butdoes include self-gravity. They confirmed that higher cosmicray diffusivity results in faster growth rates, and that modeswith a wavevector parallel to the magnetic field direction areleast stable.There have also been many numerical studies of the Parkerinstability. To date most numerical models have been ideal:the terms involving the gas viscosity in the Navier-Stokesequation, and the magnetic diffusivity in the magnetic induc-tion equation, are set to zero; notable exceptions are Hanaszet al. (2002) and Tanuma et al. (2003), which did include fi-nite resistivity. Numerical models have generally assumedthat the vertical gravity is uniform and discontinuous at themidplane and/or do not include cosmic rays. Matsumoto et al.(1988) and Basu et al. (1997) used 2D simulations to inves-tigate the nonlinear evolution of the instability, and showedthat the midplane crossing modes dominate. Kim et al. (1998)found that filamentary or sheetlike structures formed when thesystem approaches a new equilibrium after 10 yr in 3D mod-els. Disk rotation was included in the 3D simulations of Chouet al. (1997) and Kim et al. (2001), who showed that the Cori-olis force results in a twisting of the magnetic loops and mayalso help to randomize the initially uniform magnetic field. a r X i v : . [ a s t r o - ph . GA ] N ov L. F. S. R
ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER
3D models including self-gravity and differential rotation, al-lowing for the development of Jeans and swing instabilities aswell as the Parker instability, were used by Kim et al. (2002)to show that a mixture of midplane symmetric and midplanecrossing modes can develop in the ISM. The above simula-tions all excluded a cosmic ray component, which was addedby Hanasz & Lesch (2000, 2003) in 3D, to investigate the trig-gering of the Parker instability by a supernova remnant, andby Kuwabara et al. (2004), who used 2D MHD simulationsto confirm the connection between cosmic ray diffusivity andthe Parker instability growth rate found in their earlier linearanalysis. Mouschovias et al. (2009) allowed for phase transi-tions in a 2D MHD model (without cosmic rays), found thatthe midplane crossing mode was favoured, and argued that anon-linear triggering of the Parker instability by a spiral den-sity shock-front could allow the instability to form gas cloudson a timescale compatible with observations.The instability produces buoyant loops of a large-scalemagnetic field at a kiloparsec scale. These ‘Parker loops’ areexpected to lie largely in the azimuthal direction (the directionof the large-scale field). This corresponds to the ‘undular’modes (with wavevector parallel to magnetic field
BBB ) whichare expected to dominate over the ‘interchange’ modes (withwavevector perpendicular to
BBB ) in most linear analyses (see,e.g. Matsumoto et al. 1993, for a discussion of these modes).(Indeed, some authors identify the Parker instability specifi-cally with these undular modes, rather than simply in terms ofthe fundamental buoyancy mechanism.) There have been nu-merous attempts to detect such Parker loops in observations,but they have met with little success . This may be surpris-ing, since the instability is strong and generic. Therefore, weconsider here the detailed development and possible observa-tional consequences of the instability in galactic disks.We simulate the evolution of a simple model of a sectionof a galactic disk, where gas, magnetic fields and cosmicrays evolve in the gravitational field of stars and dark mat-ter halo. Cosmic rays are described in the advection-diffusionapproximation with a diffusion tensor aligned with the localmagnetic field. The simulation setup is aimed at facilitat-ing comparisons with previous analytical works, as well ascapturing clear signatures of the Parker instability which canlater be compared to more advanced simulations and observa-tions. For these reasons, we use the isothermal approximation(where there is no possibility of confusion with thermally-driven instabilities), and omit rotation and rotational shear(to avoid any confusion with other instabilities, such as themagneto-rotational instability). We also refrain from usingany source terms for the magnetic field and the cosmic rayenergy density, since this could lead to artificial spatial biases.This paper is organized as follows. In Section 2 the modelis described in detail. In Section 3 the results are presented:in § § § § Recently, magnetic loops and filaments have been detected in WMAPand Planck data (Vidal et al. 2015; Planck Collaboration et al. 2015). How-ever, these are likely a consequence of interstellar turbulence driven by super-novae remnants (Planck Collaboration et al. 2015, 2014; Mertsch & Sarkar2013) and not of Parker instability. observing Parker loops from synchrotron emission of edge-on galaxies, or from the Faraday rotation measure signal fromface-on galaxies. Finally, in Section 4 we state our conclu-sions. MODEL DESCRIPTIONAll the numerical calculations in this work were performedusing the open source, high-order finite difference
PENCIL code , designed for fully non-linear, compressible magne-tohydrodynamic simulations. The cosmic ray calculationsemployed the COSMICRAY and
COSMICRAYFLUX modules(Snodin et al. 2006).2.1.
Basic equations and cosmic ray modeling
We solve equations for mass conservation, momentum, andmagnetic induction, for an isothermal gas assuming constantkinematic viscosity ν and magnetic diffusivity η :D ln ρ D t = − ∇ · UUU , (1)D UUU D t = ggg − ∇ ( p g + p cr ) ρ + jjj × BBB ρ + ν (cid:2) ∇ UUU + ∇ ( ∇ · UUU ) + W · ∇ ln ρ (cid:3) , (2) ∂ AAA ∂ t = UUU × BBB + η ∇ AAA , (3)where ρ is the gas density, UUU is the gas velocity,
AAA is the mag-netic vector potential,
BBB = ∇ × AAA is the magnetic field (with ∇ · BBB = jjj = ∇ × BBB ( c / π ) is the electriccurrent density, with c the speed of light in vacuum, W is therate of strain tensor, W i j = ∂ U i ∂ x j + ∂ U j ∂ x s − δ i j ∇ k U k , p g = ρ c is the gas pressure, c s is the adiabatic speed ofsound, p cr is the cosmic ray pressure and ggg is the gravitationalacceleration. We neglect the effects of rotation and rotationalshear.Cosmic rays are modeled using a fluid approximation (e.g.Parker 1969; Schlickeiser & Lerche 1985) where the cosmicray energy density ε cr is governed by ∂ ε cr ∂ t + ∇ · ( ε cr UUU ) + p cr ∇ · UUU = − ∇ · FFF , (4)with FFF the cosmic ray flux defined below. The cosmic raypressure that appears in equation (2) is obtained from theequation of state p cr = ε cr ( γ cr − ) , with γ cr the cosmic ray adiabatic index. We adopt γ cr = / / ≤ γ cr ≤ / FFF is introduced in a non-Fickian formjustified and discussed by Snodin et al. (2006), τ ∂ F i ∂ t = − κ i j ∇ j ε cr − F i , (5) Documentation and download instructions can be found at http://pencil-code.nordita.org/ . ARKER INSTABILITY IN DISK GALAXIES τ can be identified with the decorrelation time of thecosmic ray pitch angles, and κ is the diffusion tensor, κ i j = κ ⊥ δ i j + ( κ (cid:107) − κ ⊥ ) ˆ B i ˆ B j , (6)where a circumflex denotes a unit vector. When τ is negligi-ble, equation (5) reduces to the Fickian diffusion flux, but thelarge diffusivity of cosmic rays makes the non-Fickian effectsboth physically significant and numerically useful. When τ is finite, the propagation speed of cosmic rays remains finite(e.g., Bakunin 2008). This prescription also avoids the pos-sibility of singularities at X-type magnetic null points where κ i j is undefined (Snodin et al. 2006).For numerical solution, the maximum time step δ t forthe cosmic ray physics might be constrained either by anadvective-style Courant condition based on the cosmic rayflux equation (5), or by a diffusive-style Courant conditionbased on the cosmic ray energy density equation (4). (The lat-ter would be appropriate for the Fickian case τ = δ t ≤ C adv δ z / c cr , or δ t ≤ C diff δ z / κ (cid:107) , with δ z the smallest numerical grid spacing, and c cr = ( κ (cid:107) / τ ) / the characteristic speed for the cosmic ray diffusion along thelocal magnetic field; here C adv and C diff are empirical con-stants. For similar advective and diffusive processes in theother equations (based on speeds c s , c A , U , and diffusivi-ties ν , η , respectively), we adopt the constants C adv = . C diff = .
06. For the relatively large value of κ (cid:107) used in thefollowing, this value of C diff results in a significant reduc-tion of δ t , to a value which numerical experiments show tobe unnecessarily small. I.e. an appropriate empirical valueof C diff for the cosmic ray diffusion would be significantlygreater than our standard value of 0.06. (This is perhaps un-surprising, given the anisotropic nature of our cosmic ray dif-fusion; κ (cid:107) only applies along magnetic field lines, where cos-mic ray energy density gradients tend to be small.) On theother hand, our standard value of C adv = . δ t (allowing the cosmic ray physics tocontrol the timestep when necessary, but avoiding unnecessar-ily small timesteps), and achieves this without the introduc-tion of an additional empirical constant. In the calculationsreported here, we use this advective-style constraint.2.2. Numerical domain and initial conditions
We use a local Cartesian box to represent a portion of thegalactic disk, with x , y and z representing the radial, az-imuthal and vertical directions, respectively. The domain isof size 6 kpc ×
12 kpc × . x , y and z directions,with the galactic midplane at the center of the vertical range − .
75 kpc ≤ z ≤ .
75 kpc. We use a grid of 256 × × ( δ x , δ y , δ z ) ≈ ( . , . , . ) . The height of the domain has beenchosen to be much larger than the characteristic vertical scaleof the instability (see below).The initial state is an isothermal disk under magnetohydro-static equilibrium in an external gravity field with the gravita-tional acceleration given by g z ( z ) = − π G Σ tanh (cid:16) zh (cid:17) , (7)where h =
500 pc is the the scale height of warm gas, Σ =
100 M (cid:12) pc − is the total surface mass density of the disk, and G is Newton’s gravitational constant.The relative contributions of the gas, magnetic and cosmicray to the total pressure in the initial state are parameterized using α = p B p g and β = p cr p g , (8)where p B = B / ( π ) is the magnetic pressure. The initialstate is unstable for any values α + β > γ g − = γ g = T = π G Σ m p µ hk B ( + α + β ) , (9)where m p is the proton mass, µ is the mean molecular mass(here µ = k B is Boltzmann’s constant. Under magne-tohydrostatic equilibrium with the gravity field adopted, theinitial gas density is given by ρ ( z ) = Σ g h sech (cid:16) zh (cid:17) , (10)where Σ g is the gas surface density.When the pressure scale heights of magnetic field and cos-mic rays are identical, as imposed by equation (8), the initialenergy density distribution of cosmic rays is ε cr ( z ) = π G ΣΣ g ( γ cr − ) sech (cid:16) zh (cid:17) β + α + β . (11)The initial magnetic field, assumed to be azimuthal, is ob-tained from A x = π h (cid:18) G ΣΣ g α + α + β (cid:19) / arctan (cid:104) tanh (cid:16) z h (cid:17)(cid:105) , A y = A z = , (12)which leads to B y = π G ΣΣ g sech (cid:16) zh (cid:17) α + α + β , B x = B z = . (13)This initial state is not in exact magnetohydrostatic equilib-rium, except for α = β =
0, due to the magnetic diffusion andperpendicular cosmic ray diffusion. For the parameters used,the diffusive time scales are much longer than the timescale ofthe Parker instability, so any diffusive evolution of the initialstate is inconsequential (see Section 2.4).The initial state is perturbed by noise in the velocity field:at each mesh point a random velocity was specified witha magnitude uniformly distributed in the range 0 ≤ U ≤
14 sech ( z / h ) km s − . The velocity perturbation is subsonic,with the rms value of order 1 km s − throughout the domain.This value was chosen to avoid unnecessarily long transientsin the onset of the instability (during which time the initialstate might evolve significantly), while avoiding directly driv-ing the instability.2.3. Boundary conditions
We assume periodicity at the horizontal ( x and y ) bound-aries. The height of the domain, Z = .
75 kpc, is muchlarger than the characteristic vertical scale of the instability,so that the precise boundary conditions at the top and bot-tom of the domain do not significantly affect the results; theseare adopted as follows. Symmetry about the top and bottom L. F. S. R
ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER
Table 1
Various timescales of the problem, with respect to vertical disk crossing times.Time scale Definition MagnitudeSound speed τ s = hc s = (cid:18) h γ g π G Σ (cid:19) / ( + α + β ) / . (cid:18) + α + β (cid:19) / Alfv´en speed a τ A = hc A = (cid:18) h π G Σ (cid:19) / (cid:18) + α + βα (cid:19) / . (cid:18) + α + β α (cid:19) / Viscous τ ν = h ν . τ η = h η . τ cr , (cid:107) = h κ (cid:107) . τ cr , ⊥ = h κ ⊥ . b τ cr , τ = hc cr = h / (cid:113) κ (cid:107) / τ . a Initial value. b Based on parallel diffusion.
Table 2
Values of parameters adoptedQuantity Symbol ValueDisk gas surface density Σ g (cid:12) pc − Disk total surface density Σ (cid:12) pc − Disk scale-height h . τ ν . × cm s − Magnetic diffusivity η . × cm s − CR parallel diffusivity κ (cid:107) . × cm s − CR perpendicular diffusivity k ⊥ . × cm s − boundaries is imposed for U x and U y , and antisymmetry for U z : ∂ U x ∂ z = ∂ U y ∂ z = U z = | z | = Z , (14)making these boundaries impenetrable and stress free. Themagnetic vector potential is symmetric for A z , antisymmetricfor A y and generalized antisymmetric for A x : ∂ A x ∂ z = A y = ∂ A z ∂ z = | z | = Z . (15)These conditions are similar to those for a perfectly conduct-ing boundary (which would require A x = A x has been relaxed to avoid overlydistorting the initial background profile A x ( z ) . Similarly, tominimize distortions to the initial background density profile,generalized antisymmetry is imposed on ln ρ at the top andbottom of the domain: ∂ ln ρ∂ z = | z | = Z . (16)This condition has no obvious physical interpretation but itplaces no strong constraints on the density distribution. Thecosmic ray energy density and flux at the top and bottom boundaries satisfy ∂ ε cr ∂ z = , ∂ F x ∂ z = ∂ F y ∂ z = F z = , (17)which corresponds to vanishing vertical flux in the Fickianlimit.The boundary conditions have been chosen largely for sim-plicity. Because of the impenetrability of the surface of thedomain, there could in principle be reflection of waves and/orthe artificial accumulation of buoyant material, but no suchproblems emerged in the results reported below (which fo-cus primarily upon the early phases of evolution of the Parkerinstability). We have also rerun one of our simulations in adomain with larger vertical size to confirm that the domainis large enough to prevent spurious boundary effects over thetimescales explored (see appendix A).2.4. Parameters and timescales
The parameters used, including our fiducial choices of thevarious diffusivities involved, are shown in table 2; the re-sulting timescales of interest for our disk are shown in ta-ble 1. These values of Σ and Σ g correspond to an ini-tial midplane gas density of ρ ( ) = . × − g cm − ;for α = β =
1, the corresponding initial azimuthal fieldstrength is B y ( ) ≈ µ G and the cosmic ray energy density, ε cr ( ) ≈ . × − erg cm − .Our choice of τ =
10 Myr for the cosmic ray flux correla-tion timescale ensures that it is smaller than the characteris-tic timescale of the Parker instability (see Section 3), as wellas the sound and Alfv´en speed timescales (see table 1). Theimpact of choosing smaller values for τ was examined, andconfirmed to have negligible impact on the results. This valueof τ , together with our other choices of parameters, resultsin the Strouhal number associated with the cosmic ray flux,St = c cr τ / h = √ κ (cid:107) τ / h ≈ .
45, being of order 1. (And notethat this estimate, involving parallel diffusivity with the ver-tical scale height, is arguably an overestimate. Typical hori-zontal lengthscales are longer than h ; see Section 3.1, below.)Thus we do not expect wavelike behaviour in these simula-tions; the diffusion of cosmic rays does not differ significantly ARKER INSTABILITY IN DISK GALAXIES ν and magneticdiffusivity η correspond to a magnetic Prandtl number Pm = ν / η =
2; the effects of varying magnetic Prandtl number areinvestigated in Section 3.3. The values adopted for ν and η are very significantly larger than expected molecular values(and the resulting Pm ∼ O ( ) is very much smaller than thegalactic regime associated with molecular values, Pm (cid:29) τ cr , (cid:107) and τ cr , τ in table 1)are significantly shorter than the timescales of the Parker in-stability, and the sound and Alfv´en speed timescales; we aretherefore in a reasonable regime for modelling galactic disks,yet may expect some deviations from those analytic resultsobtained on the basis of infinitely fast propagation of cosmicrays along fieldlines.The perpendicular diffusivity of cosmic rays, on the otherhand, was somewhat reduced from expected values to avoidthe escape of the cosmic ray energy density before the devel-opment of the instability. As noted in Section 2.2, the initialstate will evolve under the diffusion of magnetic field and theperpendicular diffusion of cosmic rays; but the timescales forthese processes, 4.77 Gyr and 2.39 Gyr, are very long com-pared to the other timescales of interest, making the slow dif-fusive evolution of the background state insignificant. RESULTS AND DISCUSSIONThree-dimensional renderings of the runs of the simulationwith ( α , β ) = ( , ) and ( α , β ) = ( , ) are shown in Fig-ure 1, in panels (a) and (b) respectively. The images show thesnapshots when the kinetic energy in each simulation matchedthe thermal energy (which is fixed due to the isothermal con-dition); this corresponds to simulation time 330 Myr for thecase without cosmic rays (panel a) and 200 Myr for the caseof equal pressure contribution of cosmic rays, magnetic fieldsand thermal pressure (panel b).As expected (and as quantified further below), one can seethat the presence of cosmic rays leads to an increased num-ber of substructures in the density field. The magnetic fieldlines, in green, display the characteristic ‘Parker loops’, whichalso clearly span smaller length scales in the case with cosmicrays.3.1. Dependence on magnetic and cosmic ray pressures ( α and β ) In Figure 2, the time evolution of the root-mean-square ve-locity U rms , and of the kinetic energy E kin , are shown for dif-ferent choices of α and β . A period of exponential growth, as-sociated with the development of the instability, can be clearlyidentified. Increases in the pressure ratio of magnetic fields, α , and of cosmic rays, β , enhance the instability, leading tolarger growth rates.The dependence of the growth rate of U z , Γ , on α and β isexamined in Figure 3, where it is compared with the analyti-cal results of Giz & Shu (1993, GS93) and Ryu et al. (2003,R+03). Under the assumption of equipartition of initial mag-netic pressure ( α = . − ,at β =
0, to 26 . − , at β =
1. To compare our find- ings to previous analytical results, one has to take into ac-count the effect of a finite vertical domain with impenetrableboundaries, which limits the maximum vertical wavelength ofthe perturbation, λ z . In order to illustrate this, we show threesets of analytic growth rate curves for different choices of λ z .For arbitrarily large vertical wavelengths GS93 predict fastestgrowing modes with growth rates 28 . − and 52 . − at β = ,
1, respectively, which differ significantly from ourresults. If, instead, the maximum attainable vertical wave-length in our domain is considered, λ z = . β = Γ / d β , found from the simulation results isclearly smaller. It should be noted that GS93 assumes infinitediffusion of cosmic rays parallel to the field and no perpen-dicular diffusion, and both of these factors would lead us toexpect smaller growth rates (Ryu et al. 2003).The analytical results of Ryu et al. (2003), which take intoaccount finite cosmic ray diffusivity, give a better match forthe slope d Γ / d β , but have smaller growth rates because of theassumption of a constant gravitational field; the importance ofthe gravity profile in this respect is clear from Kim & Hong(1998), who considered uniform, linear and ‘realistic’ fields(the latter being similar to our own).In the lower panel of Figure 3, the weaker dependence ofthe growth rate on α is shown, for a fixed β = .
25. This rela-tively weak dependence is expected for realistic gravity fields(Kim & Hong 1998), but highlights the importance of cosmicrays (and the importance of realistic background states) forthis instability.Of particular interest are the disturbances to the initial stateof the system defined in Section 2.2. These can be expressedin terms of dimensionless quantities, defined as s ( z ) = ρ ( z ) − ρ ( z ) ρ ( z ) , (18) uuu ( z ) = UUU ( z ) c s , (19) bbb ( z ) = BBB ( z ) − BBB ( z ) B , y ( z ) , (20) e cr ( z ) = ε cr ( z ) − ε cr , ( z ) ε cr , ( z ) , (21)where the subscript 0 denotes the quantity at t =
0. Thesehave the same form as the perturbations defined by Giz & Shu(1993) and provide an estimate for the degree of non-linearityof the system.In Figure 4 slices at constant height ( z = . β . The selected snapshots correspond to theinstant when the kinetic energy in the domain reaches 70 percent of thermal energy, to allow comparisons of the spatialstructure on an equal footing (i.e. to account for the differenttimescales of three runs), with the actual time indicated at thetop of each column. The system is still in the linear phase(with amplitudes typically being below 20%) but a significantamount of non-trivial structure is already visible. The spatialdistribution of the substructures does not change with time ,only their amplitudes vary. The observed structures are com- For a movie, see http://doi.org/8j3 .
L. F. S. R
ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER B (a) Density B (b) Density
Figure 1.
Renderings of the simulations at the time where the kinetic energy is close to the thermal energy. Panel (a) shows the outcome of the choice α = . β = . t = α = . β = . t = . × − gcm − . The green lines represent magnetic field lines, color-coded by field magnitude following the color-bar in the plots, which is units of 0 . µ G.The slices show the density in units of 10 − gcm − . Only the top half (i.e. the volume between the midplane and the top boundary) of the simulation domain isshown. ARKER INSTABILITY IN DISK GALAXIES
50 100 150 200 250 300 350 t [Myr] -1 u r m s [ k m s − ] α =1 . , β =0 . , Γ=22 Gyr − α =0 . , β =0 . , Γ=21 Gyr − α =1 . , β =0 . , Γ=26 Gyr − α =1 . , β =0 . , Γ=30 Gyr − α =1 . , β =1 . , Γ=38 Gyr −
50 100 150 200 250 300 350 t [Myr] -3 -2 -1 E k i n / E t h m α =1 . , β =0 . α =0 . , β =0 . α =1 . , β =0 . α =1 . , β =0 . α =1 . , β =1 . Figure 2.
The top panel shows the evolution of the the root mean square ofthe velocity, U rms , for selected choices of α and β as indicated. The bottompanel shows the evolution of the kinetic energy in the whole simulation box(normalised by the fixed thermal energy). The horizontal dashed line corre-sponds to equipartition of thermal and kinetic energy; the analysis in the textis for states before this point, when E kin = . E therm . prised of over-dense sheets of inflowing gas surrounded byoutflowing, magnetized, cosmic ray-rich, under-dense gas.In Figures 5 and 6, we examine the power spectra of per-turbations to the magnetic and velocity field, comparing thecases of ( α , β ) = ( , ) and ( α , β ) = ( , ) , at the snapshotwhen E kin = . E thm . In the absence of cosmic rays, thetwo dominant modes are ( k x , k y ) = ( . − , . − ) and ( . − , . − ) . The k x value of the first mode actuallycorresponds to the domain size in the x -direction (i.e. 2 π / L x ).There is significant spread in k x , which is consistent with theweak dependence of the growth rate with wavenumber notedin the analytic calculations. There is less spread in k y however,with energy being concentrated close to the wavenumber as-sociated with the fastest growing mode of Giz & Shu (1993), k y , GS , which is the dotted curve plotted. Our calculations favorslightly smaller wavenumbers (larger wavelengths), which isunderstandable given the inclusion of non-zero diffusivities.In the presence of cosmic rays, there is a general de-crease of the wavelength along the azimuthal direction. Tofacilitate the comparison, the two most prominent modesin the β = ( k x , k y ) = ( . − , . − ) and ( . − , . − ) . While most of the energy is stillclose to the relevant azimuthal wavenumber from Giz & Shu(1993), k y = k y , GS (here plotted as the dashed line), there ismore spread, showing that other modes are also being sig-nificantly excited. Again, the dominant modes tend to havewavenumbers slightly smaller than the value for the idealcase. β Γ [ G y r − ] (a) α =1 modelGS93, λ z → ∞ GS93, λ z = 6 kpc GS93, λ z = 3 . R+03 α Γ [ G y r − ] (b) β =0 . modelGS93, λ z → ∞ GS93, λ z = 6 kpc GS93, λ z = 3 . R+03
Figure 3.
Dependence of the growth rate of the velocity field with pressureratio. Panel (a): varying cosmic ray pressure ratio, β (with fixed α = α (with fixed β = . λ z , are shown with different linestyles, as indicated. The fastestgrowing modes predicted by Ryu et al. (2003, R+03) are shown by the bottomdash-dotted curve. Closer to the midplane there is less overall power and lessrelative power in the fastest growing modes; this is discussedin Section 3.2, below.3.2.
Dependence with height (z)
The evolution of the rms of some components of the quan-tities defined in equations (18)–(21) is shown in Figure 7, fordifferent z . The vertical components of the velocity and mag-netic fields show clear exponential growth phases. There issmall, but non-negligible, variation in the growth rates calcu-lated for varying z (with growth rates slightly smaller for low z ). The reason for this variation is that, as discussed in the pre-vious Section, these solutions do not correspond to a singleunstable mode; we are not here solving the linearized prob-lem. (And, in fact, even if we had a single mode, our back-ground state is slowly evolving away from the initial state.)Nevertheless, these variations do not prevent a meaningfuloverall growth rate for the instability being calculated. Thegrowth is clearest in the components u z and b z shown, as thesecomponents are most directly driven by the instability. Othervariables (including s , also shown) are more subject to non-linear effects, including the mixing of different modes and theevolution of the background state. The latter is a significanteffect, as the average density in the midplane ( z =
0) at the L. F. S. R
ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER y [kpc] x [ k p c ] β =0 . , t =0 .
33 Gyr s y [kpc] x [ k p c ] u z y [kpc] x [ k p c ] b z y [kpc] x [ k p c ] β =0 . , t =0 .
24 Gyr s y [kpc] x [ k p c ] u z y [kpc] x [ k p c ] b z y [kpc] x [ k p c ] e cr y [kpc] x [ k p c ] β =1 . , t =0 .
19 Gyr s y [kpc] x [ k p c ] u z y [kpc] x [ k p c ] b z y [kpc] x [ k p c ] e cr Figure 4.
Slices at z = . s , to the vertical component of the velocity field, u z and to the magnetic field, b z andto the cosmic rays energy density, e cr . The snapshots at the instant when E kin = . E thm . time E kin = . E thm is on average 3% / /
5% higher thanthe initial state for the cases α = β = / . /
1, due tothe net sinking of gas produced by the Parker instability. Sim-ilarly, there is a net rise of cosmic rays, with the cosmic rayenergy density at that time being 15% /
11% smaller than theinitial state, for α = β = . / s , u x , u y , b x , b y and e cr symmetricand u z and b z antisymmetric under the mapping z → − z ; orantisymmetric (the ‘even’ modes, allowing midplane cross-ings), with s , u x , u y , b x , b y and e cr antisymmetric and u z and b z symmetric under this mapping. Note that these perturba-tions are with respect to the initial state, which in these termsis itself symmetric (‘odd’). For some systems, (e.g. Horiuchiet al. 1988; Basu et al. 1997) it is found that the antisymmet-ric mode of instability is preferred; and it is notable that thisallows an effective factor of two decrease in the horizontalseparation of high density patches (since the antisymmetricdensity perturbation allows adjacent peaks in the ( y , z ) -planeto be at ( y , z ) and ( y + λ y / , − z ) ; rather than at ( y , ± z ) and ( y + λ y , ± z ) , as for the symmetric case). This has implicationsfor the establishment of realistic spacings in galactic den-sity concentrations (i.e. giant molecular clouds); and also for the timescales of evolution of these features, since the lineargrowth rates depend upon the wavelength. Some authors (e.g.Basu et al. 1997; Mouschovias et al. 2009) consider the non-linear evolution of the most unstable pure symmetry mode.However, Giz & Shu (1993), Kim & Hong (1998) notethat for a system with a ‘realistic’ gravity profile, such as westudy, the growth rates of the most unstable modes show nopreference for either type of symmetry. And in any event,as the solution evolves nonlinearly (and with interaction withthe symmetric initial state), the antisymmetric modes excitesymmetric components, and so purely antisymmetric pertur-bations cannot be sustained in the nonlinear regime.A simple visual inspection suggests that our nonlinear so-lutions do not exhibit either pure symmetry or antisymmetrywith respect to reflection in the midplane. Note that our weakrandom initial velocity is asymmetric, so does not seed eitherpure symmetry of perturbation. Nevertheless, if one symme-try was genuinely preferred, we would expect the weak initialcomponent of the other symmetry simply to decay. (And werea pure symmetry solution to be obtained from a pure symme-try initial condition, it could not be considered robust to theinevitable fluctuations of a real system.)To better examine this, we studied how the magnetic energyis distributed between symmetric and antisymmetric modes.For simplicity, we restrict the analysis to the x and z compo-nents, which are zero at the start of the simulation and avoid ARKER INSTABILITY IN DISK GALAXIES k x [kpc − ] α = . , β = . k y [ k p c − ] z = −
500 pc k x [kpc − ] k y [ k p c − ] z = −
150 pc k x [kpc − ] k y [ k p c − ] z =150 pc k x [kpc − ] k y [ k p c − ] z =500 pc k x [kpc − ] α = . , β = . k y [ k p c − ] z = −
500 pc k x [kpc − ] k y [ k p c − ] z = −
150 pc k x [kpc − ] k y [ k p c − ] z =150 pc k x [kpc − ] k y [ k p c − ] z =500 pc Figure 5.
Two-dimensional power spectrum of b z = B z / B y , . The power spectra were computed over slices of fixed z , as indicated in each panel. The top rowshows the power spectra for α = β =
0. The bottom row shows the power spectra for α = β =
1. Different colors correspond to different powers, asshown in the colorbar. Contour lines (drawn without taking the log) highlight major the peaks. The short- and long-dashed orange curves show the k y associatedwith the fastest growing mode obtained from the calculations of Giz & Shu (1993) for α = β = ,
1, respectively. k x [kpc − ] α = . , β = . k y [ k p c − ] z = −
500 pc k x [kpc − ] k y [ k p c − ] z = −
150 pc k x [kpc − ] k y [ k p c − ] z =150 pc k x [kpc − ] k y [ k p c − ] z =500 pc k x [kpc − ] α = . , β = . k y [ k p c − ] z = −
500 pc k x [kpc − ] k y [ k p c − ] z = −
150 pc k x [kpc − ] k y [ k p c − ] z =150 pc k x [kpc − ] k y [ k p c − ] z =500 pc Figure 6.
Same as Figure 5 but for the vertical component of the velocity perturbation, u z = U z / c s . the complication of the symmetry of the background state.The symmetric part is B ( s ) x ( x , y , z ) = ( B x ( x , y , z ) + B x ( x , y , − z ) ) , (22) B ( s ) z ( x , y , z ) = ( B z ( x , y , z ) − B z ( x , y , − z ) ) , (23)and the antisymmetric, B ( a ) x ( x , y , z ) = ( B x ( x , y , z ) − B x ( x , y , − z ) ) , (24) B ( a ) z ( x , y , z ) = ( B z ( x , y , z ) + B z ( x , y , − z ) ) . (25)From these, we compute the total energy associated witheach magnetic field component, E ( s / a ) x / z = (cid:82) ( B ( s / a ) x / z ) / π d V , inthe simulation domain and construct the following diagnostic0 L. F. S. R ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER
50 100 150 200 250 300 350 t [Myr] -3 -2 -1 u z β =0 . z =0 . , Γ =18 . − z =0 . , Γ =18 . − z =1 . , Γ =21 . −
50 100 150 200 250 t [Myr] -3 -2 -1 u z β =1 . z =0 . , Γ =27 . − z =0 . , Γ =28 . − z =1 . , Γ =36 . −
50 100 150 200 250 300 350 t [Myr] -3 -2 -1 b z z =0 . , Γ =18 . − z =0 . , Γ =19 . − z =1 . , Γ =22 . −
50 100 150 200 250 t [Myr] -3 -2 -1 b z z =0 . , Γ =27 . − z =0 . , Γ =29 . − z =1 . , Γ =37 . −
50 100 150 200 250 300 350 t [Myr] -3 -2 -1 s z =0 . z =0 . z =1 .
50 100 150 200 250 t [Myr] -3 -2 -1 s z =0 . z =0 . z =1 . Figure 7.
Root mean square evolution of the perturbations for models with α = z , as indicated in eachpanel. The panels in the left column shows the results in the absence of cosmic rays, while the right-hand column shows the case of β =
240 260 280 300 320 t [Myr] D β =0 . B x B z
120 135 150 165 180 195 t [Myr] D β =1 . B x B z Figure 8.
Time evolution of D (equation 26), in the cases of β = β = α = D = D = − quantity D = E ( s ) x / z − E ( a ) x / z E ( s ) x / z + E ( a ) x / z . (26)Thus, D = D = − D is shown. In the β = > D > − .
15 ). The presence of cosmic rays makes anypreference even less clear.We therefore conclude that pure symmetry with respect toreflections in the midplane should not be expected in the non- linear case. In this respect, the situation is rather like thatwith the horizontal structures, noted above; in the nonlin-ear regime, a spectrum of modes, of differing wavenumbers ( k x , k y ) and of differing midplane-symmetries, must be ex-pected. We note that this was appreciated by Giz & Shu(1993), in their discussion of continuum modes (to which oursolutions should correspond; since, as noted by Kim & Hong(1998), these have growth rates significantly higher than thediscrete modes, so must be expected to dominate). They ob-served that the analysis of individual modes is of limited rel-evance, with an initial value treatment (such as that adoptedhere) being more appropriate; and with the finally realisedstate depending on initial conditions and nonlinear effects.As Giz & Shu (1993) also note, however, the modal analysisremains valid with respect to conclusions about exponentialgrowth (as confirmed by Figure 7), since we are still workingwith a complete set of functions.3.3. Dependence on diffusivities ( ν , η , κ (cid:107) ) The previous results were for our fiducial choices of con-stant kinematic viscosity ν fiducial ≈ . × cm s − andconstant magnetic diffusivity η fiducial ≈ . × cm s − ,corresponding to magnetic Prandtl number Pm = ν / η = . ν fiducial and 0 . η fiducial . There are smallincreases in the growth rate, ∼ − (or (cid:46) / / ν and η by factors 2 and 3, respectively, leads to decreasesin the growth rate (cid:46) ARKER INSTABILITY IN DISK GALAXIES η [10 cm s − ] Γ [ G y r − ] ν [10 cm s − ]0 . . . . . ν [10 cm s − ] Γ [ G y r − ] η [10 cm s − ]0 . . . . . P m Γ [ G y r − ] κ [10 cm s] Γ [ G y r − ] Figure 9.
The dependence of the growth rate of U rms in the model with α = β = η ; the viscosity, ν ; the magneticPrandtl number, P m , and the cosmic ray diffusivity, κ (cid:107) . The fiducial valuesare shown as black circles. rate with increasing diffusivities might be expected on ener-getic grounds; but the choice of diffusivities clearly has onlya small effect on the growth rate, particularly when comparedto its variation with α and (particularly) β .Together, the variations with ν and η show no clear depen-dence on the magnetic Prandtl number in the range Pm = . (cid:29)
1; but along with the agreement with the ideallinearized studies, this nevertheless suggests that the instabil-ity is not strongly influenced by diffusion.In the lowest panel of Figure 9 we show the impact on thegrowth rate of varying the cosmic ray diffusivity along themagnetic field, κ (cid:107) . In agreement with the results of Ryu et al.(2003), the difference between realistic values and even largervalues (i.e. moving towards the approximation of infinite par-allel diffusivity) is small: doubling the value of κ (cid:107) from ourfiducial value leads to a 9% increase in the growth rate. (De-creasing κ (cid:107) from the fiducial value gives a slightly larger ef-fect; but this direction is moving away from galactically rea-sonable values.)3.4. Observational signatures
Having a satisfactory model where Parker instability is thedominant effect, we now look for simple possible observa-tional signatures and their dependence on the cosmic ray con-tent. We start by examining the archetypal signature of theParker instability: the so-called Parker arches/loops. Naively,one would expect the loops to be present in (polarized) syn-chrotron intensity maps of edge-on galaxies. Thus, we showin Figure 10 polarized synchrotron intensity at λ = y [kpc] z [ k p c ] ∆ x =6 kpc P I [ . × µ G e r g c m − ] y [kpc] z [ k p c ] ∆ x =47 pc P I [ . × µ G e r g c m − ] y [kpc] z [ k p c ] ∆ x = 6 kpc (with small scale random field) P I [ . × µ G e r g c m − ] Figure 10.
Polarization view of the model with α = β =
1. The shadedcontours in the three panels show the polarized intensity, PI , in arbitraryunits; dashes are oriented perpendicularly to the angle of polarization andtheir lengths show the degree of polarization. For each panel, the integrationwas over the interval indicated above that plot. The bottom panel show theeffect of the inclusion of a gaussian random field. y [kpc] z [ k p c ] ∆ x =6 kpc ( PI − › PI fi z ) / › PI fi z [%] Figure 11.
Fluctuations in the polarized intensity relative to the average ateach z are shown, for the solution analyzed in Figure 10. appendix B for the details of the calculation). In the top panel,which shows the result of integrating over a thin slice of width ∆ x =
47 pc, the polarization angles are affected by the undu-lations of the magnetic field lines. The polarized and totalsynchrotron intensities, however, are dominated by the diskstratification, and the variations at the same disk height are al-most negligible. Observing an isolated thin slice of an edge-on galactic disk is not possible. Thus, in the second panelthe integration is performed over the whole simulation box( ∆ x = z is shown: a very weak signal ( (cid:46) z -averageare found ( ∼ . ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER we check the effect of the presence of a turbulent random fieldcomponent, by adding to each point a gaussian random fieldwith a root mean square value equal to the magnitude of themodel field at that height z . The presence of the random fieldremoves any remaining (weak) signal of the Parker instabilityin the polarization.The prospect of obtaining direct evidence of Parker loopsfrom this kind of synchrotron observations is further com-plicated by the curvature of the galaxy, which would furtherweaken the signal after the integration. Figure 10 also as-sumes the observations are made at a wavelength of 1 cm,which minimizes Faraday rotation effects.An alternative approach is to consider the rotation measure(RM) to probe the fluctuations in the density and magneticfield produced by the Parker instability (Figure 4). Assumingthe case of a face-on galaxy, the RM can be computed from φ = ( .
812 rad m − ) (cid:90) n e ( z ) − B z µ G d z , (27)where n e is the number density of electrons in the path (forsimplicity, we assumed the plasma was completely ionized).In Figure 12, Faraday maps obtained from the simulationare shown. The integration was performed over the whole do-main (i.e. a distant external radio source is assumed). Thereare sharp peaks in the rotation measure signal, with a maxi-mum | φ | ≈
85 rad m − in the β = | φ | ≈
74 rad m − in the β = . | φ | ≈
91 rad m − in the β = E kin = . E thm ) and that with the current setupno steady state is reached. More informative is the spatial dis-tribution of the RM structures, which changes negligibly withtime. In Figure 13 the structure function (SF) of the RM isplotted, computed parallel and perpendicular to the directionof the initial magnetic field ( y -direction), i.e. D (cid:107) ( l ) = (cid:104) [ φ ( x , y ) − φ ( x , y + l )] (cid:105) , (28) D ⊥ ( l ) = (cid:104) [ φ ( x , y ) − φ ( x + l , y )] (cid:105) . (29)There are systematic differences in the SFs along and per-pendicular to the large scale field, and features which are de-pendent on β . The first peak of the SF is generally larger inthe direction along the field than in the direction perpendicu-lar to it: for β =
0, the first peak occurs at l peak ≈ D ⊥ ,while l peak ≈ . D (cid:107) . When β is increased, the peaks ofthe SF parallel to the field are shifted towards smaller scales,with l peak ≈ .
95 kpc for β =
1. In the direction perpendicularto the initial field, the SF is flattened with the increase in thecosmic ray density.In Figure 14 this is examined for varying α . There is neg-ligible variation on the peak positions if only the magneticfield strength is varied. Although this may be surprising, itis consistent with the relatively weak effect of varying mag-netic buoyancy shown in Figure 3, and again highlights theimportance of incorporating cosmic rays into these studies. CONCLUSIONSWe have performed a set of three-dimensional numericalsimulations exploring the growth rates and length scales asso-ciated with the Parker instability, in a section of a galactic diskwith realistic vertical structures, and including cosmic raysvia a fluid approximation. Aiming for clarity and to avoidany bias (via unintended forcing), our system is allowed to evolve passively from its initial state; i.e. there is no artificialmaintenance of cosmic rays and magnetic fields by imposedsources. This evolution of the background state represents aslight difference from the linearized systems investigated inrelated analytical studies, but the evolution is very slow com-pared with the timescales of the Parker instability; and we dif-fer from such studies in more fundamental ways, in any event,incorporating finite diffusivities and retaining nonlinear inter-actions.Nevertheless, the growth rates and length scales we obtain— and their dependences on the model parameters — qual-itatively agree well with such linearized studies. The quan-titative differences can be clearly understood in terms of thedifferences between the models (the linearized studies beingideal, and in some cases using alternative background pro-files), or simply due to the nonlinear nature of our solutions.Our calculations invariably produce non-trivial spatialstructures, significantly involving multiple modes in the in-stability, and with these modes being fully three-dimensional(rather than pure undular or interchange modes) and with nopreference for either pure symmetry about the midplane. Thisunderlines the importance of three-dimensional calculations,if signatures of the instability structures, for comparison withobservations, are sought. And the clear effect of the cosmicrays on the growth rate and wavelengths of the instability (asseen when varying our parameter β ), highlights the impor-tance of including this ingredient of the ISM in such calcula-tions.Although our calculations use high values for viscous andmagnetic diffusivities, in comparison with galactic disks, ourinvestigation with varying values suggests that the instabilityis not being significantly affected by these diffusivities.For simplicity, but in contrast to real galaxies, our calcu-lations do not involve rotation. Rotation is normally found tostabilize against the Parker instability, but we note that Kim &Hong (1998) argue that this will be less the case for a cosmicray-driven instability, cf. instabilities with purely magneticand gas buoyancy. And preliminary runs with rotation (notreported here) confirm that the effect is minor. Similarly, ourcalculations do not involve rotational shear (differential rota-tion), which is significant in disk galaxies. Foglizzo & Tagger(1994, 1995) find differential rotation to weaken the Parkerinstability; the specific effects on our calculations should beconsidered in future work.As a first step towards investigating the observational impli-cations, we calculate synthetic polarized intensity and Fara-day rotation measures maps from our simulations, and com-pute the structure functions associated with the latter. Wefind that it is very unlikely that radio observations of edge-ongalaxies would be able to clearly detect any structures pro-duced by the instability, because of the averaging along theline-of-sight and the masking of any signal by the presence ofthe random (small-scale) magnetic field. However, our resultssuggest there may be strong signatures in Faraday rotationmeasures of face-on galaxies. And our preliminary analysissuggests that, when combined with independent data aboutthe disk scale height, the correlation scales inferred from suchrotation measure maps may be a useful observational diagnos-tic for the cosmic ray content of galaxies.A comprehensive study based on the present model, solelydedicated to the observational signatures of Parker instabilityin galaxies, is underway. ARKER INSTABILITY IN DISK GALAXIES y [kpc] x [ k p c ] β =0 . , t =0 .
33 Gyr y [kpc] x [ k p c ] β =1 . , t =0 .
19 Gyr φ [ r a d m − ] Figure 12.
Faraday rotation measure maps for different β as indicated (for α = E kin = . E thm ). l [kpc] D ( l ) [ r a d m − ] α =1 . , β =0 . α =1 . , β =0 . α =1 . , β =0 . α =1 . , β =1 . l [kpc] D ( l ) [ r a d m − ] α =1 . , β =0 . α =1 . , β =0 . α =1 . , β =0 . α =1 . , β =1 . Figure 13.
Structure functions computed from Faraday rotation measuremaps for runs with α = β varying as indicated in the legend. The toppanel shows the structure function along the x -axis (i.e. perpendicular to theinitial magnetic field) while the bottom panel shows the structure functionalong the y -direction. ACKNOWLEDGEMENTSWe thank the referee for many constructive comments onthe original version of this paper. We also thank Ann Mao foruseful discussions. LFSR has been supported by STFC (grantST/L005549/1) and acknowledges support from the Euro-pean Commission’s Framework Programme 7, through theMarie Curie International Research Staff Exchange SchemeLACEGAL (PIRSES-GA-2010-269264). This work waspartially supported by the Leverhulme Trust (PRG-2014-427). This work made use of the facilities of N8 HPC pro-vided and funded by the N8 consortium and EPSRC (GrantNo.EP/K000225/1). The Centre is co-ordinated by the Uni-versities of Leeds and Manchester. This research has made l [kpc] D ( l ) [ r a d m − ] α =0 . , β =0 . α =0 . , β =0 . α =0 . , β =0 . α =1 . , β =0 . l [kpc] D ( l ) [ r a d m − ] α =0 . , β =0 . α =0 . , β =0 . α =0 . , β =0 . α =1 . , β =0 . Figure 14.
Structure functions computed from Faraday rotation measuremaps for runs with β = .
25 and α varying as indicated in the legend. Toppanel shows the structure function along the x -axis (i.e. perpendicular to theinitial magnetic field) while the bottom panel shows the structure functionalong the y -direction. use of NASA’s Astrophysics Data System.APPENDIX A. COMPARISON WITH A LARGER DOMAINTo ensure that the vertical size of the domain is largeenough to avoid boundary effects, the simulation was also runin a larger domain, with a vertical range − < z < . U rms is, however, ∼
6% higher in the taller domain for the reasonsdiscussed in section 3.1.4 L. F. S. R
ODRIGUES , G. R. S
ARSON , A. S
HUKUROV , P. J. B
USHBY , A. F
LETCHER
100 120 140 160 180 200 220 t [Myr] -1 u r m s [ k m s − ] L z = 3 . L z = 6 .
100 120 140 160 180 200 220 t [Myr] -3 -2 -1 E k i n / E t h m L z = 3 . L z = 6 . Figure 15.
Evolution of the root mean square velocity (left panel)and kinetic energy (right panel), comparing the standard simulation do-main, ( L x , L y , L z ) = ( , , . ( L x , Ly , L z ) = ( , , ) , for α = β = When comparing figure 16 with figures 5 and 6, it can beseen that the enlargement of the vertical dimension leads onlyto very small changes in the power spectra; and these can beattributed to small variations in the random initial conditions. B. POLARIZATION PROPERTIESWe discuss here the details of the calculation of the polar-ized intensity in section 3.4.First the Stokes parameters are computed, I ( y , z ) = (cid:90) ∆ x s ( rrr (cid:48) ) d x (cid:48) , (B1) Q ( y , z ) = (cid:90) ∆ x p s ( rrr (cid:48) ) cos (cid:2) ψ ( rrr (cid:48) ) (cid:3) d x (cid:48) , (B2) U ( y , z ) = (cid:90) ∆ x p s ( rrr (cid:48) ) sin (cid:2) ψ ( rrr (cid:48) ) (cid:3) d x (cid:48) , (B3)where the intrinsic polarization degree was assumed to be p = .
75, and the local polarization angle is obtained from ψ ( rrr ) = π + arctan (cid:20) B z ( rrr ) B y ( rrr ) (cid:21) (B4) + .
81 rad (cid:18) λ (cid:19) (cid:90) ∆ xx (cid:20) n e ( rrr (cid:48) ) − (cid:21) (cid:20) B x ( rrr (cid:48) ) µ G (cid:21) (cid:18) d x (cid:48) (cid:19) , To compute the emissivity it was assumed that the energy den-sity of cosmic rays is proportional to the number density ofcosmic rays, thus the quantity s ( rrr ) = ε cr ( rrr ) (cid:2) B y ( rrr ) + B z ( rrr ) (cid:3) . (B5)corresponds to the emissivity in arbitrary units.From equations (B1)–(B3) it is then possible to compute thequantities shown in Figures 10 and 11. The polarized intensityis given by PI = (cid:113) ( Q + U ) , (B6)the observed polarization angle is Ψ =
12 arctan (cid:18) UQ (cid:19) , (B7)and the polarization degree p = PI / I . (B8)REFERENCES Bakunin, O. G. 2008, Turbulence and Diffusion: Scaling Versus Equations(Springer-Verlag Berlin Heidelberg)Basu, S., Mouschovias, T. C., & Paleologou, E. V. 1997, ApJ, 480, L55Chou, W., Tajima, T., Matsumoto, R., & Shibata, K. 1997, PASJ, 49, 389Foglizzo, T., & Tagger, M. 1994, A&A, 287, 297, astro-ph/9403019——. 1995, A&A, 301, 293, astro-ph/9502049Giz, A. T., & Shu, F. H. 1993, ApJ, 404, 185Hanasz, M., & Lesch, H. 2000, ApJ, 543, 235——. 2003, A&A, 412, 331, astro-ph/0309660Hanasz, M., Otmianowska-Mazur, K., & Lesch, H. 2002, A&A, 386, 347Horiuchi, T., Matsumoto, R., Hanawa, T., & Shibata, K. 1988, PASJ, 40, 147Kim, J., & Hong, S. S. 1998, ApJ, 507, 254, astro-ph/9807073Kim, J., Hong, S. S., Ryu, D., & Jones, T. W. 1998, ApJ, 506, L139,astro-ph/9808244Kim, J., Ryu, D., & Jones, T. W. 2001, ApJ, 557, 464, astro-ph/0104259Kim, W.-T., Ostriker, E. C., & Stone, J. M. 2002, ApJ, 581, 1080,astro-ph/0208414Kuwabara, T., & Ko, C.-M. 2006, ApJ, 636, 290, astro-ph/0506137Kuwabara, T., Nakamura, K., & Ko, C. M. 2004, ApJ, 607, 828,astro-ph/0402350Kuznetsov, V. D., & Ptuskin, V. S. 1983, Ap&SS, 94, 5Matsumoto, R., Horiuchi, T., Shibata, K., & Hanawa, T. 1988, PASJ, 40, 171Matsumoto, R., Tajima, T., Shibata, K., & Kaisig, M. 1993, ApJ, 414, 357Mertsch, P., & Sarkar, S. 2013, J. Cosmology Astropart. Phys., 6, 41,1304.1078Mouschovias, T. C., Kunz, M. W., & Christie, D. A. 2009, MNRAS, 397,14, 0901.0914Parker, E. N. 1966, ApJ, 145, 811——. 1967, ApJ, 149, 517——. 1969, Space Sci. Rev., 9, 651——. 1979, Cosmical Magnetic Fields: their Origin and their ActivityPlanck Collaboration et al. 2015, ArXiv e-prints, 1502.01582——. 2014, ArXiv e-prints, 1409.6728Ryu, D., Kim, J., Hong, S. S., & Jones, T. W. 2003, ApJ, 589, 338,astro-ph/0301625Schlickeiser, R., & Lerche, I. 1985, A&A, 151, 151Snodin, A. P., Brandenburg, A., Mee, A. J., & Shukurov, A. 2006, MNRAS,373, 643, astro-ph/0507176Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212, astro-ph/9807150Tanuma, S., Yokoyama, T., Kudoh, T., & Shibata, K. 2003, ApJ, 582, 215,astro-ph/0209008Vidal, M., Dickinson, C., Davies, R. D., & Leahy, J. P. 2015, MNRAS, 452,656, 1410.4438
ARKER INSTABILITY IN DISK GALAXIES k x [kpc − ] α = . , β = . k y [ k p c − ] z = −
500 pc k x [kpc − ] k y [ k p c − ] z = −
150 pc k x [kpc − ] k y [ k p c − ] z =150 pc k x [kpc − ] k y [ k p c − ] z =500 pc k x [kpc − ] α = . , β = . k y [ k p c − ] z = −
500 pc k x [kpc − ] k y [ k p c − ] z = −
150 pc k x [kpc − ] k y [ k p c − ] z =150 pc k x [kpc − ] k y [ k p c − ] z =500 pc Figure 16.
Two-dimensional power spectra computed for a simulation with an enlarged domain in the z-direction ( L z = α = β =
1. The top rowshows the power spectrum of u z = U z / c s while the bottom row shows the power spectrum b z = B z / B y ,0