Astronomy&Astrophysicsmanuscriptno.Lorentz c ESO2009 (cid:13) March17,2009 An implicit advection scheme for modeling relativistic shocks with high Lorentz factors Hujeirat, A.A.,Keil,B.W.andHilscher, P.P. 9 0 ZAH,LandessternwarteHeidelberg-Ko¨nigstuhl, 0 Universita¨tHeidelberg,69120Heidelberg,Germany 2 Received.../Accepted... r a M ABSTRACT 7 Aims. Anumericallystableandaccurateadvectionschemeformodelingshockfrontsmovingatultra-relativisticspeedsisfunda- mentallyimportantforunderstandingthethermodynamicsofjetsemanatingfromthevicinityofrelativisticobjectsandtheoriginof 1 thehighgamma-rays. Methods. Inthispaperwepresentaspatiallythird-orderaccurateadvectionschemethatiscapableofcapturingshock-frontswith ] E diverseLorentzfactors.Theconsistencyandaccuracyoftheschemeareinvestigatedusingtheinternalandtotalenergyformulation H ingeneralrelativity. Results. Usingthetotalenergyformulation,theschemeisfoundtobeviableformodelingmovingshocksatmoderateLorentzfac- . h tors,thoughwithrelativelysmallCourantnumbers.InthelimitofhighLorentzfactors,theinternalenergyformulationincombination p withafine-tunedartificialviscosityismuchmorerobustandefficient. - Weconfirmourconclusionsbyperformingtestcalculationsandcomparetheresultswiththeanalyticalsolutionsoftherelativistic o shocktubeproblem. r t Keywords.Methods:numerical–hydrodynamics–relativistic–ultrarelativisticshocks s a [ Received ... / Accepted ... General relativity: neutron stars – Anninos,Fragile, 2003; Komissarov, 2004; DelZannaetal., 1 v blackholes–X-raybursts,Methods:numerical–hydrodynam- 2007; Mignoneetal., 2007; Hujeiratetal., 2008, see also the 5 ics–relativistic referencestherein). 2 Received.../Accepted... These studies have enriched the field of computational astro- 0 physicswith numericaltechniquesandusefulstrategiesfor ac- 3 curatecapturingrelativisticshockfrontsusingboththeinternal 3. 1. Introduction and total energy formulation (e.g., DeVilliers,Hawley, 2003; 0 Mignoneetal.,2007).Nevertheless,mostofthemethodsappear Plasmas moving at relativistic speeds have been observed in 9 toencounterseverenumericaldifficultieswhenevertheLorentz 0 diverse astrophysical phenomena, such as in supernova explo- factorbecomeslarge,i.e., Γ 5.Thereasonisthatsmallveloc- v: sions,injetsemanatingfromaroundneutronstars,microquasars ityerrorsaremagnifiedby Γ≥2 astheyentertheotherequations, and in active galactic nuclei (see Livio, 2004; Marscher, i henceinducingsuper-proportionalerrorsintheestimationofthe X 2006;Fenderetal.,2007,andthereferencestherein).Thebulk othervariables.Recipessuchasfurtherreducingthe gridspac- Lorentzfactor,Γ,oftheobservedjet-plasmas,inmostcases,has r ingandaccordinglythetime stepsize mylead tostagnationof a been found to increase with the mass of the central relativistic the solution procedure, particularly if the methods used are of object. Ejected plasmas from aroundneutronstars have Γ = NS time-explicittype. (1 3),inmicroquasars Γ = (1 4) andinAGNs/QSOs µQSO OΓ −= (5 10). In gamma-rayOburs−ts however,plasmas are The robustness of time-implicit methods at capturing shocks AGN considereOdto−movewithaLorentzfactoroftheorderofseveral propagating with high Lorentz factors has been barely inves- hundreds(e.g.,Piran,2006),i.e., Γ = (100 1000). tigated, hence the purpose of the present paper. Moreover, we GRB O − showthatourtime-implicitmethodincombinationwithamod- Jet-plasmas considered to decelerate via self-interaction or in- ifiedthirdorderMUSCLadvectionschemeiscapableofmodel- teractionwiththesurroundingmediaintheformofshocksand eventuallybecomeefficientsourcesfortheproductionoftheob- ingshockspropagatingwith Lorentzfactors Γ 10with time ≥ stepsizescorrespondingtoCourantnumberslargerthan1. served highenergygamma-raysand probablyfor the energetic cosmicrays(Sikora,1994). In Sec. 2 we describe the hydrodynamical equations solved Modeling the formation of relativistic shock by means of in the present study. The solution method relies on using the HD and MHD simulation has been the of focus of numer- 3D axi-symmetric general relativistic implicit radiative MHD ous studies during the the last two decades (Hawleyetal., solver (henceforth GR-I-RMHD), which has been described 1984b; Alayetal., 1999; Fontetal., 2000; MartiandMu¨ller, in details in a series of articles (e.g., Hujeiratetal., 2008; 2003; DeVilliers,Hawley, 2003; Gammieetal., 2003; Hujeirat&Thielemann.,2009).InSec.(3)wejustfocusonsev- eralnewaspectsoftheadvectionscheme.Theresultsofthever- Send offprint requests to: A. Hujeirat, e-mail: ificationtests ispresentedin Sec.(4)andendupwithSec. (5), [email protected] wheretheresultsofthepresentstudyissummarized. 2 A.Hujeiratetal.:Anumericalschemeforcapturingultrarelativisticshocks 2. The1Dgeneralrelativistichydrodynamical where ˜tot =(ut)2ρh. equations Thus,oEnce ˜tot isfound,thepressurecanbecomputedfromthe E knownconservativevariablesasfollows: In the present study we consider the set of Euler equations in one-dimensionandinflatspacetime: γ 1 ˜tot utD P= − E − . (7) γ (ut)2 ∂D+ (DVr)=0 t r ∂ M +∇ · (M Vr)= ∂ P (1) On the other hand, using the internal energy formulation for t r r r r ∂ tot+∇ · ([ tot+P]−Vr)=0, modeling moving shocks requires the inclusion of an artificial t r E ∇ · E viscosity for calculating their fronts accurately. This modifica- where = 1 ∂ √ g, D((cid:17) ρut) is the relativistic density, tion is necessary in order to recover the loss of heat generally ∇r· √−g r − producedthroughtheconversionofkineticenergyintointernal M istheradialmomentum:M (cid:17) Du,whereD (cid:17) Dh,ut isthe r r t energyattheshockfronts. time-likevelocity1 ,Vr = uµ/ut isthe transportvelocity,“h”is Inthiscase,theRHSofEquation(2)shouldbemodifiedtoin- therelativisticenthalpyh = c2+ε+P/ρ. tot = (ut)2ρh Pis E − cludetheartificialheatingterm: thetotalenergywhichisthesumofkineticandthermalenergies ofthegas.“P”isthepressureofidealgas:P= (γ 1)ρε,εand Q+ =η (∂ utVµ)2, (8) − art art µ γdenotetheinternalenergyandadiabaticindex,respectively. whereη istheartificialviscositycoefficient. art The inclusion of Q+ implies an enhancement of the effective The reader is referred to Sec. (2) of Hujeirat et al. (2008), art pressure.Therefore,thethermodynamicalpressure,andalsothe where the general relativistic hydrodynamical equations and enthalpy,shouldbemodifiedtoincludeanartificialpressureof their derivations in the Boyer-Lindquist coordinates are de- theform: scribed.However,wecontinuetousethesecoordinates,though inthelimitofflatspace. P =P+P =P+η ∂ (utVµ). (9) tot art art µ In the case that the mechanical energy is conserved, then the totalenergyreducestoanevolutionaryequationfortheinternal Note that the artificial pressure enters the momentumequation energy: intheformof P ,whichscalesas ∆(η V).Thisisequiva- art art ∇ ∼ lent to activating a second order viscous operator at the shock d fronts, whose effect is then to transport information from the ∂ d+ ( dVr)= (γ 1)E [∂ut+ (utVr)] (2) tE ∇r· E − − ut t ∇r· downstream to the upstream regions, so to enhance the stabil- ityofthetransportschemeinsuchcriticaltransitions. where d =utP/(γ 1). E − Furthermore, there is a remarkable similarity of how Eq. (2) 2.1.DeterminingtheLorentzfactor deals with d and ut. Taking this similarity into account, the E equationcanbere-formulatedtohavethefollowingform: OurtestcalculationsshowedthattheLorentzfactorcanbebest determinedfromthenormalizationcondition,inwhichboththe ∂t¯d+ r (¯dVr)= γD( r Vr), (3) conservativevariablesand transportvelocities are involved.To E ∇ · E − ∇ · clarifythispoint,thenormalizationconditionreads: where ¯d =Dlog[ d(ut)γ 1]. − E E InEquation(1),theenergyequationdescribesthetime-evolution 1 =uµu µ of the total energy tot. Our numerical solution method relies − =(utVµ)(Mµ) firerssptoonndsienlgecetnintrgietsheiEnprtihmeaJryacvoabriiaanb.leWsaentdhecnalcituelraattiengtothreeicrocvoer-r = uD¯t[Mt+VD¯rMr+VθMθ+VϕMϕ] (10) thecorrectcontributionsofthedependingvariables,suchasthe = uD¯t[{Mt−ggtttϕMϕ}+VθMθ+VϕMϕ] transportvelocityandpressure.Therefore,theadvectionopera- = ut[(D¯ut)+VθM +(Vϕ gtϕ)M ], torcanbedecomposedintotwopartsasfollows: D¯ gtt θ − gtt ϕ whereD¯ = ρhut andwhichcanbedirectlydeterminedfromthe ∂t tot+ r ( totVi)= r (PVi). (4) continuityandtheenergyequations,Vµisthetransportvelocity E ∇ · E ∇ · andgµνaretheelementsofthemetric. Thisequationissolvedasfollows:Solvefor totusingthepres- The final form of Eq. (10) is the quadratic equation: a (ut)2 + surefromthelastiterationlevel.Once tot anEd“D”areknown, but+c = 0, fromwhich the Lorentz factorcan be determined E wemayproceedtoupdatethepressureasfollows: completely. We note that since the state of the conservative variables P= Etot−utD . (5) D, Mµ, d depends strongly on the transport velocity Vµ, it is γ (ut)2 1 thereforeEnecessarytoconsiderVµwhencomputingtheLorentz γ 1 − − factorutasdescribedinEq.(10). Thisvalueof“P”istheninsertedagainintotheRHSofEquation (4)tocomputeacorrectedvaluefor tot. TheeffectofiterationcanbereducedEonthecostofusingpres- 3. Amodifiedthirdorderadvectionschemeforhigh surevaluesfromthelasttimestepusingthefollowingalternative Lorentzfactors formulation: Mosttime-implicitadvectionschemesgenerallyrelyonupwind- ∂ ˜tot+ (˜totVi)=∂P, (6) ing to stabilize the transport near critical fronts. Transport of t r t E ∇ · E quantitieshere can be mediated with CFL numberslargerthan 1 ut isalsothegeneralrelativisticLorentzfactor,whichreducestoΓ unity.Therefore,theadvectionschemeemployedshouldbein- inflatspace. dependentofthetimestepsize.Indeed,themonotoneupstream A.Hujeiratetal.:Anumericalschemeforcapturingultrarelativisticshocks 3 centeredschemesforconservationlaws,orsimplytheMUSCL where ξ is an additional weighting function and LqSR,SL are j schemes,obeytheseconditions(seeHirsch,1988,andtherefer- LUIS corrections of “q” at the interface r (see Fig. 2, which j encestherein).However,inordertomodelmovingshockswith arecomputedinthefollowingmanner: high Lorentz factors accurately, it is necessary to modify the suthcpehsetnrmoearemm, asroliezgtahitoaiontnsthcsehooninudflidtoirobmne,aift.iueor.nt,h-cflearouwlsiamlfirittoyemdcoitnnhdeaitcdicooonwr.dnasntrceeamwittho LLqqSjSjLR==PPllll====jjjj+−22qqllQQikki,,==kkjj,+,kk2==jj−rr2mj−rrrmirjmkm−−r,rmkmk , ::iiff VVjj >≤00. (15) AclassicaladvectionschemeoftheMUSCL-typegivesthefol- i −k lowinginterfacevalues: However, we note that high order advection schemes gener- ally degenerate into lower ones at the shock fronts or discon- qSR =q 1[(1+κ)∆q +(1 κ)∆q ] :if V 0 tinuities, in order to maintain monotonicity of the advection qjS−L1/2 =qj−1+−14[(1 κ)∆qj+−1(1+κ−)∆q ]j−2 :if Vj ≤>0,(11) scheme,whichinturnenhancestheproductionofnumericalen- j 1/2 j 4 − j j+1 j tropy/diffusion at these critical transitions. Therefore, in most − of the test calculationspresented here it was necessary to have where “q” is the transported physical quantity, κ is a switch ξ >1/2,orevenveryclosetounitytoenableconvergence. off/on parameter that specify the accuracy needed and ∆q = j q q .Notethattheschemeisofsecondorderforκ= 1,0,1 j j+1 − − andofthirdorderforκ=1/3. 4. Results In orderto enablean accuratecapturingof shockfrontspropa- Inthissectionwepresentresultsofvariousverificationtestsus- gatingwithhighLorentzfactorsonanon-uniformgriddistribu- ingdifferentenergyformulationsincombinationwiththeadvec- tion,thefollowingmodificationshavebeenperformed: tionschemedescribedintheprevioussection.Thewidelyused hydrodynamical test is the one-dimensional relativistic shock if V 0: qSR =q 1[(1 σ)∆q +(1+σ)∆q ], j ≤ j j−1− 4 − j j−1 tube problem (henceforth RSTP). The advantages of this test where σ=2κ(∆qRR ∆q0R)/[(∆qRR)2+(∆q0R)2+ǫ]. problem is that the obtained numerical solutions can be com- j · j j j (12) if V >0: qSL =q + 1[(1+σ)∆q +(1 σ)∆q ], pared with the correspondinganalytical solution with arbitrary wjhere σ=j2κ(∆jq0L4 ∆qLL)/[(∆qj0L)2+−(∆qLLj)+21+ǫ], largeLorentzfactors. j · j j j For further details on the background of this test problem and thewaytheanalyticalsolutionsareobtainedwereferthereader where to(MartiandMu¨ller,2003). ∆qRR =(xm xm)(q q )/(xm xm ) TheRSTPisaninitialvalueproblem,inwhichthesolutionde- ∆qj0R =(xmj−1−xmj )(qj−2− qj−)1/(xmj−2−xmj)−1 pends strongly on the initial conditions. So the domain of cal- ∆q0jL =(xmj−2−xmj−)1(q j−1−q)j/(xmj−1−xmj) (13) culationsisdividedintoa rightandleftregions(see Fig.2).In ∆qLjjL =(xjmj−1−−jx+mj1)(qjj−−1−qj+1j)/(xjmj−1−−xmj+j1), ittnhhieetirlaeilgfvthetrelrogecgioiitonyn0isx≤sce<txtxo≤v:xawcnei:sshweetevPseReryt=wPLh(2e=/re3.4)0×/31,0−ρ6L,ρ=R1=01an.dThine andwhereǫ isasmallnumbersettoavoiddivisionbyzero. q q L j-2 q q qSL qSR R j+2 j j q j+1 q dx j-1 q j Vj+1 Vj Vj-1 x rm rm rm rm rm j+2 r j+1 r j r j-1 r j-2 r xin xa xc xb xout j+2 j+1 j j-1 j-2 Fig.2. Grid spacing versus position. The actual changes of the vari- btFioeiugs.n.d1a.riAessacnhdemthaetilcocdaetsiocnripotfiothneosfcasleavreqraulanfitnitiiteesv”oql”uamnedctheellsv,eltohceii-r aaunbndliefosformorcx.cTu>rhxeicnin:tiPhtieRal=incto(e2nr/vd3ait)lio×[nx1sa0r−≤e6adaxn:df≤oρrRxxb=,≤1w.xhcer:ePtLhe=g4r0id/3s,pρaLcin=g1i0s The accuracy of above modified scheme may be further in- Using the total energy formulation (TEF) we show in Figures creased by incorporatingthe Lagrangian Upwind Interpolation (3-13)theprofilesofthevelocityandthecorrespondingLorentz Scheme(LUIS).TheLUISstrategyisbasedonbyconstructing factorfordifferentdistributionsoftheinitialconditions. a Lagrangian polynomial of third or fourth order whose main In Fig.3 we show the profiles of the velocity and Lorentz fac- weight is shifted to the right or left, depending on the upwind tor after 0.2sec using the total energyformulation(see Eq. 1). direction. The final combined LUIS-MUSCL scheme reads as Theinterval[0.35 x 0.75]hasbedividedinto400equally ≤ ≤ follows: spaced cells, whereas 220 cells where used to cover the exter- nalintervals,wherethevariablescontinuetoacquiretheirinitial if Vj ≤0: qSjR =ξqSjR+(1−ξ)LqSjR, (14) values. The advection scheme described in (Eq. 14) has been if V >0: qSL =ξqSL+(1 ξ)LqSL, employedtoachieve3rd orderspatialaccuracyandthedamped j j j − j 4 A.Hujeiratetal.:Anumericalschemeforcapturingultrarelativisticshocks Crank-Nicolsonmethodisusedtoachievesecondordertempo- 0.8 Analytic solution 1.45 ral accuracy (Hujeirat,Rannacher, 2001). The time step size is 0.7 620 grid points 1.40 settohaveamaximumthatcorrespondstoCFL=0.4. 0.6 1.35 1.30 TcreroemataseiesndtedtthheuenrPcohLbaunbsgtyneedas.sfTaochftiostrhienocfsroe1la0us0tei,oonwfhptihrloeecetphdreeusroseut,hrweereyviheaalrdviasebtilhnees- Velocity000...435 Lorentz Factor11..2205 1.15 Lorentz factor: ut 3.6 (Fig. 4). However, to enforce conver- 0.2 1.10 ≈ gence, it was necessary to increase the number of iteration per 0.1 1.05 timestepconsiderably. 0.00.0 0.2 0.4 0.6 0.8 1.0 1.000.0 0.2 0.4 0.6 0.8 1.0 x x Comparing the results with those of Fig. (5), we see that the solution obtained is sufficiently accurate, so that doubling the Fig.3. The profiles of the velocity V (left) and the Lorentz factor ut (right, dashed line connecting squares) overplotted on the corre- numberofgridpointsanddecreasingthetimestepsizebyhalf sponding analytical solutions (continuous blue lines) after τ = 0.2. didnotimprovetheaccuracyoftheschemesignificantly. Thenumericalsolutionshavebeenobtainedusingthetotalenergyfor- From Figs. (6) and (7), we see that incorporatingartificial vis- mulation (TEF), the initial pressure (P = 40/3, ρ = 10), (P = L0 L R cosity and artificial heating while using the TEF has a similar (2/3) 10 6, ρ = 1),700non-uniformly distributedgridpointsand effect as that of the normal numerical diffusion. In both cases CFL=×0.4.T−heaRdvectionschemeemployedhereisofthirdorderspatial thevelocityattheshockfrontattainslowervaluesthanexpected andsecondordertemporalaccuracy. fromtheanalyticalsolution. We notethatwithinthecontextoftime-implicitsolutionmeth- 1.0 4.0 Analytic solution ods,thetotalenergymethodappeartohavealimitedcapability 620 grid points 3.5 tocaptureshocksmovingwithhighLorentzfactors(Fig.8).The 0.8 rEceokaninsvoiennrcgirseentahcseaetsr,awutenitlhdikeLecrotehraeesneitnsztwefarinctthaolrEeainnse/ErΓgk2yi.n.AEIisnn,atthhcieosnkcsiaensqeeutieictneicsnee,nretghcye- Velocity00..46 Lorentz Factor32..05 2.0 essary to decrease both the time step size and the grid spacing in addition to increasing the number of iteration per time step. 0.2 1.5 This procedure may easily lead to a stagnation of the solution 0.00.0 0.2 0.4 0.6 0.8 1.0 1.00.0 0.2 0.4 0.6 0.8 1.0 method. x x On the other hand, our time-implicit solution method appears Fig.4. AsinFig(3),theprofilesofVand ut havebeenobtained us- to be morerobustif the internalenergyformulationis used. In ingtheTEF,aninitialpressurePL = 102PL0 andCFL=0.4.620non- Figures(9-13),theresultsofseveraltestcalculationsareshown, uniformlydistributedgridpointshavebeenused which agree well with the analytical solutions. Moreover, the IEF is sufficiently stable, so that running the calculations with 1.0 Analytic solution 4.0 1300 grid points CFL 1 doesnotappearto hinderthe convergenceof the im- 3.5 ≥ 0.8 plicitsolutionprocedure(seeFigs.12and13). 3.0 cTthoherercemocteaffiLjoocrrieednnrttaszwfoafbcatthocekrtahoraftitfitfihcetisaIlwEveFilslciwsosittihhteythnineecoaenrsdaseliytryttioctaoolbfisotnaleiun-ttiuothnnee. Velocity00..46 Lorentz Factor2.5 2.0 Moreover,these fine-tuningproceduresof the parametersmust 0.2 berepeatedeachtimetheratioofthepressureP /P ischanged. 1.5 L R This is a consequence of the conversion of kinetic energy into 0.00.0 0.2 0.4 0.6 0.8 1.0 1.00.0 0.2 0.4 0.6 0.8 1.0 x x heatattheshockfrontswhichcannotbedeterminedaprioriand uniquelybysolvingtheinternalenergyalone. Fig.5.AsinFig(3),theprofilesofVanduthavebeenobtainedusing the TEF, an initial pressure P = 102P , CFL= 0.4 and 1300 non- Nevertheless,usingtheIEF,non-linearphysicalprocessescanbe L L0 uniformlydistributedgridpoints. easily incorporatedinto theimplicitsolutionprocedure.Unlike intheTEFcase,suchprocessesaredirectfunctionsoftheinter- nalenergyandthereforeitismucheasiertocomputetheirentries andincorporatethemintothecorrespondingJacobian.Thishas 5. Summaryandconclusions the consequence that calculations can be performed also with CFL 1, while still maintainingconsistency with the original In this paper we have presented the results of different numer- ≥ ical studies for modeling the propagation of ultra-relativistic physicalproblem. Finally,inFig. (14) we showthe profilesofV andut usingthe shocks using the internal energy (IEF), the total energy (TEF) andpseudo-totalenergy(PTEF)formulations. compactinternalenergyformulation(CIEF,Eq.(3)withhighac- curacies. The methodis foundto display overshootingand un- OurtestcalculationsshowthatTEFismostsuitedformodeling dershootingbothatthecontactdiscontinuityandtheshockfront, thepropagationofshockswithlowtomoderateLorentzfactors, whose appearingwas difficultto preventthroughenlargingthe though with low CFL numbers, hence more suitable for time- numberof grid points, decreasing the time step, enhancingthe explicitsolutionmethods. numberofiterationpertimesteporevenbyactivatingtheartifi- On the other hand, althoughthe PTEF is similar to the TEF, it cialviscosity. was found that a larger number of iterations per time step was It should be noted here that since in the limit of large Lorentz requiredtomaintainconvergence. factors, the CFL number depends weakly on the sound speed UsingtheIEFhowever,itwasverifiedthattheimplicitsolution (see Fig. 15), increasing the global iteration per time step has procedure converges relatively fast even for a Lorentz factor more significant effect on convergencethan merely decreasing ut 1andaCFLnumberofCFL 1. ≫ ≥ thetimestepsize. A.Hujeiratetal.:Anumericalschemeforcapturingultrarelativisticshocks 5 1.0 4.0 1.0 2.4 Analytic solution Analytic solution 620 grid points 700 grid points 3.5 2.2 0.8 0.8 2.0 3.0 Velocity00..46 Lorentz Factor2.5 Velocity00..46 Lorentz Factor11..68 2.0 1.4 0.2 0.2 1.5 1.2 0.00.0 0.2 0.4 0.6 0.8 1.0 1.00.0 0.2 0.4 0.6 0.8 1.0 0.00.0 0.2 0.4 0.6 0.8 1.0 1.00.0 0.2 0.4 0.6 0.8 1.0 x x x x Fig.6. AsinFig(3),theprofilesofVand ut havebeenobtained us- Fig.10.SimilartoFig(3),theprofilesofVanduthavebeenobtained ingtheTEF,aninitialpressureP = 102P ,aCFL=0.4.Anupwind using the IEF with the initial left pressure P = 10P . The rest of L L0 L L0 advectionschemeoffirstorderspatialaccuracyhasbeenused. variablesremainedasinFig.(1)unchanged. 1.0 4.0 1.0 4.0 Analytic solution Analytic solution 620 grid points 700 grid points 0.8 3.5 0.8 3.5 3.0 3.0 Velocity00..46 Lorentz Factor2.5 Velocity00..46 Lorentz Factor2.5 2.0 2.0 0.2 1.5 0.2 1.5 0.00.0 0.2 0.4 x 0.6 0.8 1.0 1.00.0 0.2 0.4 x 0.6 0.8 1.0 0.00.0 0.2 0.4 x 0.6 0.8 1.0 1.00.0 0.2 0.4 x 0.6 0.8 1.0 Fig.7. AsinFig(3),theprofilesofVand ut havebeenobtained us- Fig.11.AsinFig(3),theprofilesofVanduthavebeenobtainedusing ingtheTEFwithaheatingsourceduetoartificialviscosity,aninitial theIEF,aninitialpressureP =102P andaCFL=1.5afterτ=0.215. L L0 pressure P = 102P , CFL= 0.4 and a spatially third order accurate L L0 advectionscheme. 1.0 7 Analytic solution 700 grid points 1.0 6 0.8 6 Analytic solution 700 grid points 5 Velocity000...468 Lorentz Factor435 Velocity000...462 Lorentz Factor432 0.2 2 0.00.0 0.2 0.4 0.6 0.8 1.0 10.0 0.2 0.4 0.6 0.8 1.0 x x 0.00.0 0.2 0.4 0.6 0.8 1.0 10.0 0.2 0.4 0.6 0.8 1.0 Fig.12.AsinFig(3),theprofilesofVanduthavebeenobtainedusing x x theIEF,aninitialpressureP =103P andaCFL=2. L L0 Fig.8. AsinFig(3),theprofilesofVand ut havebeenobtained us- ing the pseudo-total energy formulation (PTEF, see Eq. 6), an initial 1.0 18 pressureP =103P andCFL=0.4. Analytic solution L L0 700 grid points 16 0.8 14 0.8 1.45 Analytic solution 12 000...675 700 grid points 111...343005 Velocity00..46 Lorentz Factor180 Velocity 00..43 Lorentz Factor11..2205 0.2 462 0.2 1.15 0.00.0 0.2 0.4 0.6 0.8 1.0 00.0 0.2 0.4 0.6 0.8 1.0 0.1 1.10 x x 0.0 1.05 Fig.13.AsinFig(3),theprofilesofVanduthavebeenobtainedusing 0.10.0 0.2 0.4 x 0.6 0.8 1.0 1.000.0 0.2 0.4 x 0.6 0.8 1.0 theIEF,aninitialpressurePL=105PL0andaCFL=4. Fig.9.AsinFig(3):theprofilesofVanduthavebeenobtainedusing theinternalenergyformulation(IEF,seeEq.2),aninitialpressureP = L PL0andCFL=0.4. intotheJacobian. Furthermore, turbulence is an important property of most It should be stressed here that the IEF has a major draw back, astrophysical fluid flows. Thus the artificial viscosity in such as it requires a fine-tuned parameters of the artificial viscosity flows naturally will be overwhelmed by turbulent diffusion, to obtain the correct Lorentz factor at the shock fronts. On the hencerelaxingthefine-tuningproblemoftheartificialviscosity. other hand, when modeling astrophysical plasmas, the IEF is favorableovertheTEF,ascoolingandheatingprocessesarein Irrespectiveofthemethodusedforsolvingtheenergyequation, generaldirectfunctionsoftheinternalenergy,thusitiseasierto the Lorentz factor is found to be best computed from the construct the corresponding coefficients and incorporate them normalization condition in which the transport velocities and 6 A.Hujeiratetal.:Anumericalschemeforcapturingultrarelativisticshocks 0.9 2.2 Analytic solution Gammie,C.F.,McKinney,J.C.,To´th,G.,2003,ApJ,589,444-457(HARM) 0.8 620 grid points 2.0 Hirsch,C.,1988,”NumericalComputationofInternalandExternalFlows”,Vol. 0.7 I,II,JohnWiley&Sons,NewYork 0.6 1.8 Hujeirat,A.,Rannacher,R.,2001,NewAst.Reviews,45,425 Velocity000...435 Lorentz Factor11..46 HHHuuujjjeeeiiirrraaattt,,,AAA...,,,2CT0ha0mie5le,enCmzoianPndhn,C,MF,.1.-,K6K8.,,e2i1l0,0B9.,WA.,&2A00,8to,NapepweaArstronomy,13,436 Hawley,J.F.,Smarr,L.L.,Wilson,J.R.,1984a,ApJ,277,296-311 0.2 1.2 Hawley,J.F.,Smarr,L.L.,Wilson,J.R.,1984b,ApJS,55,211-246 0.1 Koide,S.,ShibataK.,Kudoh,T.,1999,ApJ,522,727 0.00.0 0.2 0.4 x 0.6 0.8 1.0 1.00.0 0.2 0.4 x 0.6 0.8 1.0 Komissarov,S.S.,2004,MNRAS,350,1431 Livio,2004,BaltA,13,273L Fig.14.AsinFig(3),theprofilesofVandutafterτ=0.175havebeen Marti,J.M.,Mu¨ller,E.,2003,LRR,6,7 obtainedusingthecompactinternalenergyformulation(CIEF,seeEq. Mignone,A.,Bodo,G.,Massaglia,etal.,2007,ApJS,170,228-242(PLUTO) 3)andtheinitialpressurePL=10PL0. Mizuno,Y.,Nishikawa,J.-I.,etal.,2006,astro-ph/0609004(RAISHIN) Piran,T.,2005,RvMP,76,1143 102 vs Sikora,M.,Begelman,M.C.,Rees,M.J.,1994,ApJ,421,153 vs (Newtonian) vs,limit=(cid:2)(cid:3)(cid:4)1c 101 vcs 100 10-110-2 10-1 100 101 102 103 (cid:0)=(cid:1)p Fig.15.β=V /cfortheNewtoniancase(V =γP/ρ;dashedline)and S S fortherelativisticcase(V =γP/ρh;solidline)versustemperature. S momentaareused. We have also presented a modified MUSCL advection scheme of third order spatial accuracy, which has been constructed to enableaccuratecapturingofrelativisticallymovingshocks.The accuracyoftheschemecan befurtherenhancedbyincorporat- ingtheLagrangianUpwindInterpolationScheme(LUIS). On the other hand, using this advection scheme for modeling astrophysicalfluid flows mayunnecessaryincrease the compu- tationalcosts,inparticularwhen3Daxi-symmetriccalculations withhighresolutionarerequired. Notingthattheappropriateparametersoftheartificialviscosity cannotbedeterminedapriori,obtainingthecorrectLorentzfac- tor may turn out to be a difficult task. In the case of standing shocks, the hierarchical solution scenario (Hujeirat, 2005) can beemployedtogradualenhancetheaccuraciesoftheschemeto obtainthecorrectLorentzfactorsattheshockfronts. In a forthcoming paper, we intend to study the formation of strong shocks in curved spacetime near the surface of ultra- compactneutronstars. AcknowledgmentThisworkissupportedbytheKlaus-Tschira Stiftungundertheprojectnumber00.099.2006. References Biretta,J.A.,1999,AAS,195,5401 Marscher,A.P.,2006,AIPC,856,1 Aloy, M.-A., Ibanez, J.M., Mart, J.M., Mu¨ller, E., 1999, ApJS, 122, 151 (GENESIS) Anninos,P.,Fragile,P.C.,2003,ApJ.Suppl.Ser.,144,Iss.2,243-257 DelZanna,L.,Bucciantini,N.,2002,A&A,390,1177-1186 DelZanna,L.,Zanotti,O.,etal.,2007,A&A,473,11 DeVilliers,J.-P.,Hawley,J.F.,2003,ApJ,589,458 Fender,R.,Koerding,E.,Belloni,T.,etal.2007,astro-ph,0706.3838 Font,J.A.,Miller,M.,Suen,W.,Tobias,M.,2000,Phys.Rev.D,61,044011 Font,J.A.,2003,LRR,6,4