An implicit advection scheme for modeling relativistic shocks with high Lorentz factors

0.66 MB·English

by  A. A. Hujeirat

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),
[e-mail] 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