Table Of Content

Draftversion January4,2013 PreprinttypesetusingLATEXstyleemulateapjv.5/2/11 SPECTRAL AND INTERMITTENCY PROPERTIESOF RELATIVISTIC TURBULENCE Jonathan Zrake and Andrew I. MacFadyen Center forCosmologyandParticlePhysics,PhysicsDepartment,NewYorkUniversity,NewYork,NY10003,USA Draft version January 4, 2013 ABSTRACT High-resolutionnumericalsimulationsareutilizedtoexamineisotropicturbulenceinacompressible fluid when long wavelength velocity fluctuations approach light speed. Spectral analysis reveals an 3 inertial sub-range of relativistic motions with a broadly 5/3 index. The use of generalized Lorentz- 1 covariantstructurefunctionsbasedonthefour-velocityisproposed. Thesestructurefunctionsextend 0 the She-Leveque model for intermittency into the relativistic regime. 2 Subject headings: hydrodynamics — methods: numerical — gamma-ray burst: general — turbulence n a J 1. INTRODUCTION relating to the non-linear behavior of the SRHD system 3 remain to be answered. For example, it has not pre- A rapidly accumulating body of astronomical obser- viously been known whether relativistic velocity excita- vations indicates the presence of relativistic turbulence ] tions survive to small scales or are suppressed near the E in a variety of astrophysical systems. Cosmological outer scales of a turbulent flow. If relativistic motions H gamma-ray bursts (GRBs) accelerate relativistic jets persist deeply into the cascade, then universal behavior withenergiesexceedingthekineticenergyofasupernova h. (1051−52ergs)to Lorentz factors γ &100 and containin- may be expected. However, relativistic statistical mea- sures are required to characterize its scaling and inter- p ternal fluctuations with δγ 2 (Piran 2004), while jets ∼ mittency properties. For example, the power spectrum - from active galactic nuclei and micro-quasars are found o ofthree-velocitycannotobeyapowerlawabovethescale to have γ & 10 and γ . 2 respectively. Recent dis- r ℓ wherethevelocitydifferencesasymptoticallyapproach t coveries indicate that a broad class of stellar explosions, γ s light speed. Traditional two-point structure functions includingthoseassociatedwithX-raytransientsandlow- a also break down at scales ℓ > ℓ . Statistical diagnostics luminosity GRBs (Kulkarni et al. 1998; Starling et al. γ [ built on the fluid four-velocity, u = (γc,γv), are thus 2011), and even some supernovae lacking high-energy 2 emission(Soderberg et al.2010;Paragiet al.2010), also needed. v accelerate significant amounts of relativistic material. In this Letter, we present high-resolution simulations 6 These relativistic flows are highly susceptible to turbu- ofisotropicrelativisticturbulenceandintroduceLorentz 6 lence due to the extremely high inferred Reynolds num- covariantscaling diagnostics to analyze its basic proper- 0 ties. Questions to be addressed include (1) Does rela- bersandtheabundantpresenceofshearduetoobserved 4 tivistic turbulence display universal behavior? (2) Can flow asymmetries. . its intermittency be described in terms of known phe- 0 In recent models for fluctuating emission from GRBs, nomenologicalmodels? (3) How deeply into the cascade 1 the relativistic turbulence is proposed to be directly ob- do velocity fluctuations remain relativistic? 2 servable, as relativistic eddies beam their radiation into Inorderto addressthese questions,numericalintegra- 1 narrowcones which sweep by the observers’line of sight : (Narayan& Kumar 2009; Lazar et al. 2009). In addi- tions of the special relativistic hydrodynamics (SRHD) v equationswereconductedonthePleiadessupercomputer tion, the central engine for short-duration GRBs likely Xi involvestheturbulentmergingofneutronstarswhichare at the Ames researchcenter on lattice sizes up to 20483, utilizing more than 2,000,000 CPU-hours. The SRHD r knownsourcesofgravitationalradiation(Hulse & Taylor a 1975). Knowledgeofthepropertiesofturbulencepresent equations are cast in conservation form for the parti- cle number Nν = 0 and energy-momentum tensor in these sources (Melatos & Peralta 2010) may aid with ν Tµν = S∇µ and closed by the adiabatic equation of detection and interpretation of gravitationalwaves from ν ∇ state p = ρe(Γ 1) where Γ = 4/3, e is the specific these sources when interferometers such as advanced − internal energy, ρ is the co-moving mass density and LIGO become operable in the near future. Relativistic Tµν =ρhuµuν+pgµν, and h=c2+e+p/ρ. The source turbulence is also likely generated by phase transitions in the early universe and may be produced in heavy ion termSµ =ρaµ ρ(e/e0)4uµ includes injection of energy − colliders and laser experiments. and momentum at the large scales and the subtraction The fundamental properties of kinematically rela- of internal energy (with parameter e0 = 0.25) to per- tivistic hydrodynamic turbulence have not, until now, mit stationary evolution. Vortical modes at k/2π 3 been investigated. Turbulence in fully relativistic are forced by the four-acceleration field aµ = duµ wh≤ich dτ MHD has been studied for decaying (Zhang et al. smoothly decorrelates over a large-eddy turnover time, 2009; Inoue et al. 2011; Beckwith & Stone 2011) and asdescribedin Zrake & MacFadyen(2012). The numer- driven (Zrake & MacFadyen 2011, 2012) flows. Also, ical scheme employed is the 5th-order accurate weighted simulations ignoring the uid rest mass have been essentially non-oscillatory (WENO) conservative finite- done in the magnetized (Cho 2005) and unmagnetized difference scheme which has been extensively tested for (Radice & Rezzolla 2012) cases. Fundamental questions SRHD (Zhang & MacFadyen 2006) and is implemented 2 J. Zrake and A. I. MacFadyen Figure 1. Slicesoftheco-movingdensity(left)andfour-velocitycomponent outofthepage(right),indicatingthestructureofisotropic relativistic turbulence at a resolution of 20483. Density is scaled to the 1/2 power to elucidate structures in the evacuated regions. The grayscaleisfromminimum(darkest) tomaximum(brightest). in the Mara code (Zrake & MacFadyen 2012) with im- proved smoothness indicators (Shen & Zha 2010). The Table 1 Simulationparameters simulations take place on the periodic cube of size L, andwereevolvedatstationarityforatleast5largeeddy Parameter Symbol Value turnovertimes,withtheexceptionofthe20483runwhich was evolved for 1/2 of a turnover time. Each model at Boxsize L - resolution N3 was initialized from stationary turbulence Drivingscalescale - L/2 at resolution (N/2)3, after which a full turnover time is MReelaantivpiasitriwcisscea-lreelativeLorentzfactor γ¯ℓrγel 03.1.160L requiredtore-establishstationarity. Table1summarizes Meanlab-frameLorentzfactor γ¯ 1.67 the simulation parameters for the 20483 run. Maxlab-frameLorentzfactor γmax 5.73 MeanrelativisticMachnumber M¯ 2.67 2. RESULTS MMeaaxnreplraetsisvuisrteictoMmacehannudmenbseitryratio Mp¯/mρ¯ca2x 02.24.19 Imagesofthefullydevelopedturbulentfieldsareshown Lower1%densityquantile - 0.027ρ¯ in Fig. 1. The co-moving density (left panel) con- Upper1%densityquantile - 4.89ρ¯ tains fluctuations above and below the mean by factors of 16. The highest density structures are sheet-like, Note. — Simulation parameters are taken from the 20483 but not fully two-dimensional and trace the front be- run. The relativistic Mach number is defined to be M = tweencloudsoffluidcolliding withrelativisticvelocities. (γβ)fluid/(γβ)sound. Quantities with over-bars are volume aver- agedateachiteration,andallquantitiesaretime-averagedduring Shocksandcontactdiscontinuitiesareevidentwherethe stationaryevolution. density field changes abruptly. In compressible turbulence, these shocks indicate the geometry of the most behavior is universal, i.e. statistically independent of aggressively dissipating structures (Kritsuk et al. 2007). energy injection and dissipation effects at the largest These shocks are more intense in a relativistic gas be- and smallest scales respectively. This range is referred cause their density contrasts may be arbitrarily large, to as the inertial sub-interval, and is often identified by σ = Γ+1 + Γ e/c2 whereas in the non-relativistic limit Γ−1 Γ−1 self-similar behavior in the power spectrum P(k) of the thedensityratioisboundedby Γ+1. Thetransversefour- velocity field. The 1941 theory of Kolmogorov (K41) Γ−1 velocity component (right panel) reveals a highly struc- (Kolmogorov 1941) for incompressible, non-relativistic turedhierarchyofshearingmotionsdowntothesmallest turbulence predicts that P(k) k−5/3 in the inertial in- ∝ scales. It also exhibits remarkable coherence over some terval. Inrelativisticflowstheamplitudeofvelocityfluc- regions,wherelargecloudsoffluidareinrelativisticmo- tuationscannotexceedc,sothepowerspectrumofveloc- tion into (or out of) the page. Fluid elements separated ity cannot obey a power law to arbitrarily large scales. by 1/100, 1/10 and 1/5 of the image are, on average, However, the physics of a relativistic cascade may still in relative motion with four-velocities of 1/2, 1 and 2 be universal, and it could be revealed by an appropri- respectively. ate choice of generalized diagnostics, such as the power Turbulent flows at sufficiently large Reynolds number spectrum of fluid four-velocity. are known to contain a range of scales over which the Figs. 2aand2bshowthepowerspectrumP(k)offour- Spectral and Intermittency Properties of Relativistic Turbulence 3 100 100 (a) (c) 10-2 k−5/3 10-2 ∝ P(k)10-4 P(k)10-4 Pc(k)∝k−1.8 10-6 10-6 10-8 10-8 (b) (d) ∝k−1.92 ∝k−1.56 100 ) ) k k ( ( P P 5/3k 5/3k10-1 Pc(k)∝k−1.8 10-2 101 102 103 101 102 103 k/2π k/2π Figure 2. (a) Power spectrum of the four-velocity field at resolutions 2563, 5123, 10243, and 20483 (increasing in weight from gray to black). (b)Thesamedata asin(a)butcompensated byk5/3 andoffset, stretched toexaggerate deviations froma5/3law. Anarbitrary verticaloffset isgiventoeach curveinordertoclarifythespectral shapebetween k/2π(=1/L)=10and100. (c)Powerspectrum ofthe Helmholtz decomposed four-velocity field. The compressive component Pc(k) (black) follows a power law with index 1.80. (d) The same dataasin(c)butcompensated byk5/3. velocity at lattice sizes of 2563, 5123, 10243, and 20483, is because a relativistic cascade contains a new dimen- where sionless number, the Lorentz factor, at each scale. As P(k)dk = X u˜k u˜∗k (1) a simple illustration, a cursory application of the K41 · k∈dk dimensional argument to relativistic eddies reveals that the Lorentz factor at a scale ℓ satisfies is computed from the discrete Fourier transform (γ 1)3(γ +1) u˜µk = N13 X uµ(xl,m,n)e−ik·xl,m,n (2) ℓ− γℓ2 ℓ =Cǫ2ℓ2/c6. (3) l,m,n where ǫ is the energy injection rate per unit mass and C of the four-velocity field. Throughout the text, boldface is analogous to the Kolmogorov constant. This relation isusedtodenotethespatialcomponentsofafour-vector. does not admit a single power law solution. Instead, ThenormalizationbyN3, whereN isthe numberoflat- solutions are γ ℓ for large γ, and γ ℓ2/3 for small ℓ ℓ tice points per direction, guarantees that P(k) satisfies γ. ∝ ∝ R P(k)dk = u u . Wefindevidenceforaninertialinter- Unlike the total fluid four-velocity, the compressive h · i val of relativistic velocity fluctuations between 1/10 and four-velocity component uµ is strictly self-similar over c 1/100 of the largest scale. As the resolution increases the inertial interval. This is evidence that a local cas- from 2563 to 10243, a shortinterval obeying the 5/3 law cade of acoustic modes is operating across those scales. emerges. But the subsequent resolution of 20483 reveals The Helmholtz decomposition of spatial components of amorefeaturefulspectrumoffour-velocity,consistingof four-velocity is done in the Fourier domain, with the abrokenpower-lawwhichissteeper (1.92)atthe largest compressive part u˜ck = u˜k k/k and the solenoidal part scalesandshallower(1.56)atmoderatescales. Thismay u˜sk = u˜k u˜ck. As show·n in Figs. 2c and 2d, the − be due to the bottleneck effect which is characterizedby power-law range of dilatational power spectrum P (k) c anaccumulationofpoweratthesmall-scaleendofthein- extends through the energy injection scale. This is be- ertialrange(Falkovich1994). Therehasbeensubstantial cause only the vorticalmodes are being directly excited. efforttodistinguishtrueinertialscalingfromthecontam- Its slope is 1.80, which is in contrast with the super- inationofbottlenecks(Beresnyak & Lazarian2009). We sonic limit of highly compressible turbulence, where all reportthatthepowerspectrumoffour-velocityisbroadly the compressive power is in shocks, and the dilatational 5/3, owing to the fact that each higher resolution adds power spectrum P (k) k−2. The power spectrum of c ∝ additional scales which, on the average, obey 5/3 scal- solenoidal P (k) four-velocity is not at all self-similar, s ing. In relativistic turbulence, asymptotically converged and accounts for the majority of power (82%) through- scaling behavior may not yield a single power law. This out the cascade. We remark that thermally relativis- 4 J. Zrake and A. I. MacFadyen 10.000 Kolmogorov (1941) 2.0 S˜ (ℓ)=2 γ 1 1.000 2 rel− S∥(ℓ)= δuµδxµ p ℓζp∥ (cid:1) (cid:0) ∥ζ31.5 p | δxµδxµ | ∝ ∥/p (cid:2)(cid:0) (cid:1) ζ 1.0 ) ℓ ( 0.100 ˜S2 0.5 She & Leveque (1994) 0.0 S˜∥(ℓ)= δ(cid:9)u ˆℓ2 0.010 2 | · | (cid:1) (cid:0) 2.0 γβ<1 γβ>1 S (ℓ)= δuµ δu p/2 ℓζp ←→ 31.5 p | µ| ∝ ζ 0.001 / (cid:2) (cid:1) 101 102 103 ζp ℓ (grid zones) 1.0 Figure 3. The generalized structure functions of four-velocity S˜2k(ℓ) (longitudinal projection, squares) and S˜2(ℓ) (circles) as a 0.5 functionoftheseparationℓforresolution20483. Theverticalgray bar marks the relativistic scale ℓγ = 0.06L, below which velocity 0.0 fluctuations becomesub-relativistic. 1 2 3 4 5 6 7 exponentp tic turbulence admits a novel mechanism for energy ex- Figure 4. Dependenceofζp ontheexponentpforthepower-law cthheanflgueidb.etwTeheenmtheechvaonritciaclalwaonrdkcdoomnepriensscivoemmproetsisoinnsgoaf scalingoflongitudinalS˜pk(top)andtotalS˜p(ℓ)(bottom)structure functions. ζp is obtained by a least-squares fit of the structure fluid volume goes into internal energy, which for a rel- functionsbetweenℓ=1/10and1/100thedomain. Allthescaling ativistic gas enhances its inertia through the coupling exponents arenormalizedbyζ3. Errorbarsindicatethe standard deviationover6equally-spacedstationarysnapshotsat10243. The of total enthalpy to fluid four-velocity expressed by the intermittency models of K41 and SL94 are shown in dashed and term ρhuµuν in the energy-momentum tensor. This ad- solidlinesrespectively. ditional mechanismof energy transfer adds qualitatively whereδuµ =uµ uµ andℓ=(δxµδx )1/2. Four-velocity new dynamics to the turbulent cascade. In particular,it 2− 1 µ pairs are chosento be simultaneous in the globalcenter- may break the statistical decoupling between the compressive and vortical cascades that was recently discov- of-momentum frame so that δx0 =0. Note that S˜p(ℓ)= ered for non-relativistic compressible turbulence (Aluie 2(γ 1)p/2 where γ is the Lorentz factor of one rel rel 2011). The signature of this effect is proposed to be an hfl|uid ele−ment| asimeasured in the rest frame of the other. enhanced coherency among the phases of similar-scale S˜ (ℓ) becomes v v p in the non-relativistic limit. p 2 1 shearing and dilating velocity modes, that can be deter- We also generalhi|ze t−he lo|ngiitudinal structure function mined from the bi-spectrum of fluid four-velocity. δuµδx Modern descriptions of turbulence have recognized S˜k(ℓ)= µ p (5) that two-point correlations of velocity are not normally p h|(δxµδxµ)1/2| i distributed around the mean fluctuating fluid velocity (Kolmogorov 1962; Kraichnan 1990; Cao et al. 1996). which can also be written as δu ℓˆp , where ℓˆ= ℓ/ℓ, h| · | i Thiseffectis referredtoasintermittency,andrepresents because of pair simultaneity. an important departure from the K41 theory. Intermit- We utilize S˜2(ℓ) to determine the scale ℓγ to which tency may be quantified by examining higher moments relativistic velocity fluctuations persist. The kinematic of the probability distribution function (PDF) of the ve- effectsofrelativityareimportantwhentherelative four- locity increments at various scales ℓ. In particular, the velocityγβbetweenfluidelementsexceeds1,i.e. γrel > h i 1994modelofShe-Leveque(SL94)(She & Leveque1994) √2 (or S˜ (ℓ) = 2√2 2 0.83). As shown in Fig. 3, 2 hashadgreatsuccessinpredictingthescalingexponents this occurs at ℓ =0.−06, o≈r 1/16 of the outermost scale. γ ζp for the velocity structure functions Sp(ℓ) of order S˜ (ℓ)showstworegimesofpower-lawscaling,havingex- p. We have found that relativistic turbulence also dis- 2 ponents of ζ = 0.98 below ℓ and ζ = 1.09 above ℓ . playsintermittency,andthattheSL94modelprovidesan 2 γ 2 γ excellent description, provided that velocity increments The longitudinal structure function S˜2k(ℓ) obeys a power areproperlydefinedtoaccommodaterelativisticfluctua- law with index ζk = 0.89 between 1/10 and 1/100 of 2 tions. Wegeneralizethetotalvelocitystructurefunction the domain. This is steep relative to the K41 theory for as incompressible non-relativistic turbulence that predicts S˜ (ℓ)= δuµδu p/2 (4) ζ =2/3,but is in close agreementwith what was found p µ 2 h| | i Spectral and Intermittency Properties of Relativistic Turbulence 5 insimulationsofhighlycompressiblenon-relativistictur- ical turbulent cascades operating far into the relativistic bulence (Kritsuk et al. 2007). regime. We find that the SL94 intermittency model extends to relativistic turbulence when the structure function is This research was supported in part by the NSF given by Eqn. 4. Fig. 4 shows that the normalized throughgrantAST-1009863andbyNASAthroughgrant exponents ζk/ζk through p = 6.5 are well described by p 3 NNX10AF62G issued through the Astrophysics Theory the SL94 relation ζp = p/9+2 2(2/3)p/3. Error bars Program. Resourcessupportingthis workwereprovided − are obtained from the standard deviation over 6 snap- by the NASA High-End Computing (HEC) Program shots at 10243. The total structure function given in through the NASA Advanced Supercomputing (NAS) Eqn. 5,whichmeasuresthestatisticsofpairwiserelative Division at Ames Research Center. We thank Andrei Lorentz factor, exhibits a far greater degree of intermit- Gruzinov and Paul Duffell for insightful discussions. tency. We propose that a new intermittency model for the statistics of relative four-velocities is required, not- REFERENCES ingthatthesestatisticsareinherentlynon-Gaussianfrom consideration of relativistic kinematics alone. That is to Aluie,H.2011, Phys.Rev.Lett.,106,174502 say, a Gaussian random velocity field does not yield a Beckwith,K.,&Stone,J.M.2011,ApJS,193,6 Gaussian PDF of γ , owing to the non-additive nature Beresnyak,A.,&Lazarian,A.2009, ApJ,702,1190 rel ofLorentzboosts. Asubsequentinvestigationwillexam- Cao,N.,Chen,S.,&Sreenivasan,K.R.1996, Phys.Rev.Lett., ine the shape of the S˜ (ℓ) PDF at fixed scales in order 77,3799 p Cho,J.2005,ApJ,621,324 to more specifically characterize the source of the inter- Falkovich,G.1994, PhysFluids,6,1411 mittency seen in Fig. 4. Hulse,R.A.,&Taylor,J.H.1975,ApJ,195,L51 Inoue,T.,Asano,K.,&Ioka, K.2011,ApJ,734,77 Kolmogorov,A.N.1941,Proceedings:MathematicalandPhysical 3. DISCUSSION Sciences, 434,9 —.1962, JFluidMech,13,82 InthisLetterwehavedescribedanumericalinvestiga- Kraichnan,R.H.1990, PhysicalReviewLetters (ISSN tion of isotropic turbulence in a kinematically and ther- 0031-9007), 65,575 Kritsuk,A.G.,Norman,M.L.,Padoan,P.,&Wagner,R.2007, mally relativistic fluid. It is found that relativistic fluc- TheAstrophysicalJournal,665,416 tuations persist to 1/16 of the outermost scale. Spectral Kulkarni,S.R.,Frail,D.A.,Wieringa,M.H.,etal.1998,Nature, analysisrevealsaninertialsub-rangeofrelativisticveloc- 395,663 ity fluctuations with a broadly 5/3 index. The compres- Lazar,A.,Nakar,E.,&Piran,T.2009,ApJ,695,L10 sive component of the four-velocity field obeys power- Melatos,A.,&Peralta,C.2010,ApJ,709,77 Narayan,R.,&Kumar,P.2009, MNRAS,394,L117 law scaling with index 1.80 over nearly two decades. Paragi,Z.,Taylor,G.B.,Kouveliotou,C.,etal.2010, Nature, The SL94 phenomenology successfully describes inter- 463,516 mittency in relativistic turbulence, provided the correct Piran,T.2004, RevModPhys,76,1143 definition for the velocity increments is used. Radice,D.,&Rezzolla,L.2012,arXiv,astro-ph.HE In closing, we remark that for sufficiently large She,Z.-S.,&Leveque, E.1994,Phys.Rev.Lett.,72,336 Shen,Y.,&Zha,G.2010,IntJNumerMethFl,64,653 Reynolds number, ℓ is larger than the dissipation scale γ Soderberg,A.M.,Chakraborti,S.,Pignata, G.,etal.2010, η so that kinetic energy must be transmitted through a Nature,463,513 sub-relativisticcascadebeforebeingthermalized. Inthis Starling,R.L.C.,Wiersema,K.,Levan,A.J.,etal.2011, casetheratioℓ /Lmaybeindependentofthedissipation MNRAS,411,2792 γ mechanism, and should thus be dictated solely by the Zhang,W.,MacFadyen, A.,&Wang,P.2009,ApJ,692,L40 Zhang,W.,&MacFadyen, A.I.2006,ApJSSeries,164,255 Lorentz factor of the energy bearing eddies. The exact Zrake,J.,&MacFadyen, A.2011,GammaRayBursts2010.AIP functional form of this dependence could be constrained ConferenceProceedings,1358,102 byaparameterstudyovertheouterscaleLorentzfactor. Zrake,J.,&MacFadyen, A.I.2012,ApJ,744,32 Doingsowouldyieldinsightintothenatureofastrophys-