Shock formation in stellar perturbations and tidal shock waves in binaries

0.23 MB·English

by  Carsten Gundlach

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed18January2013 (MNLATEXstylefilev2.2)
Shock formation in stellar perturbations and tidal shock
waves in binaries
Carsten Gundlach1 ⋆ and Jeremiah W. Murphy2
2 1School of Mathematics, Universityof Southampton, Southampton, SO17 1BJ, UK
1 2Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195-1580
0 (NSF Astronomy and Astrophysics Postdoctoral Fellow)
2

n
a 24March2011,revised23May2011
J

3
ABSTRACT

]
R We investigate whether tidal forcing can result in sound waves steepening into
S shocks at the surface of a star. To model the sound waves and shocks, we consider
. adiabatic non-sphericalperturbationsof a Newtonian perfect fluid star.Becausetidal
h
forcing of sounds waves is naturally treated with linear theory, but the formation of
p
shocksis necessarilynonlinear,we considerthe perturbations in two regimes.In most
-
o of the interior, where tidal forcing dominates, we treat the perturbations as linear,
r while in a thin layer near the surface we treat them in full nonlinearity but in the
t
s approximation of plane symmetry, fixed gravitationalfield and a barotropic equation
a
of state. Using a hodograph transformation, this nonlinear regime is also described
[
by a linear equation. We show that the two regimes can be matched to give rise to a

3 single mode equationwhichis linear but models nonlinearity in the outer layers.This
v can then be used to obtain an estimate for the critical mode amplitude at which a
1 shockformsnear the surface.As anapplication,we considerthe tidal wavesraisedby
8 the companion in an irrotationalbinary system in circular orbit. We find that shocks
9
form at the same orbital separation where Roche lobe overflow occurs, and so shock
4
formation is unlikely to occur.
.
3
0 Key words: hydrodynamics–shockwaves–stars:oscillations–stars:binaries:close
1 – methods: analytical
1
:
v
i
X
1 INTRODUCTION artefacts of the interaction with the artificial atmosphere.
r
On the other hand, these artefacts may also hide genuine
a As far as we are aware it is unknown if the tidal forces in
shocks.
a binary inspiral can create shock waves before the binary
Mass shedding in small amplitude nonlinear pertur-
objects touch, begin mass transfer or plunge. In order to
bations has been demonstrated numerically for neutron
investigate this, we have developed a quantitative criterion
stars rotating near the mass-shedding angular velocity
forthecriticalamplitudeatwhichstellarperturbationsform
(Stergioulas et al. (2006); Dimmelmeyer et al. (2004)), but
shocks that may be interesting in its own right, or in other
the minimum perturbation amplitude for this to occur has
applications.
not been quantified.
This work was originally motivated by the observation
of unsmooth fluid behaviour in the numerical simulation of We have therefore tried to obtain a quantitative cri-
anirrotational,equalmassneutronstar(NS)binarymerger terion for shock formation using a combination of stellar
(Baiotti et al. (2008); Rezzolla et al. (2010)). The simula- perturbation theory and nonlinear planar fluid dynamics.
tionsshowsurface wavesbreakingwhen theinitial dataare We consider a shock formation scenario where essentially
evolved with the cold equation of state (EOS) P = Kρ2 , radial sound waves steepen as they approach the surface
andastrongwindwhentheyareevolvedwiththeequivalent because the density and sound speed approach zero at the
hot EOS P =ρe (with initially constant entropy). As both surface. (Note that while the shape of the tidal bulges ro-
simulationsshouldbeidenticaluntilgenuineshocksform,it tatesaroundthestar,individualfluidelementsmainlymove
seems likely that both the wind and the surface waves are up and down.)
Ifsuchwavesaregeneratedbytidalforcesfromthecom-
panion,theiramplitudeandshapeisdeterminedinthebulk
⋆ Email:[email protected] ofthestar,wherealmost allthemassis.Inthisregime,lin-
(cid:13)c 0000RAS
2 Carsten Gundlach and Jeremiah W. Murphy
ear perturbation theory can beused to obtain the response 2 NONLINEAR EXTENSION OF LINEAR
of the star to the tidal force, treating its proper oscillation PERTURBATION MODES
modes as forced harmonic oscillators. For simplicity, we as-
2.1 Linear adiabatic perturbation equations
sumethatthebackgroundstarisirrotationalandspherical.
We consider linear adiabatic perturbations of a spherically
Ontheotherhand,nearthesurface thefluidgeometry symmetric static perfect fluid star in Newtonian physics,
canbeapproximatedasplane-parallel,andentropyorcom- in the frequency domain. The background is assumed to
position gradients become irrelevant compared to the den- be non-rotating and in hydrostatic equilibrium. d/dr is de-
sity gradient. In this regime we use a hodograph transform noted by a prime. The background quantities are the den-
to cast the nonlinear dynamics into a single linear second- sity ρ(r), pressure P(r), gravitational potential φ(r), grav-
order PDE. A shock forms if and only if the hodograph itational acceleration φ′(r) > 0, sound speed c(r) defined
transform becomes singular: a criterion for this can be ex- by c2 (∂P/∂ρ)s, entropy per rest mass s(r), and Brunt-
≡
aminedwithinthemodelitself(Gundlach & Please(2009)). V¨ais¨al¨a frequency N(r) definedby
Thiscriterion was testedin thenumerical evolution ofnon- N2 ρ′ P′
linear spherical perturbations of an n = 1 polytropic star . (1)
− φ′ ≡ ρ − c2ρ
andfoundtobeaccuratewithin10%(Gabler et al.(2009)).
Once formed, these shocks quickly take a universal, self- The equations of hydrostatic equilibrium for the spherical
similar form (Gundlach & LeVeque(2011)). background star are
P′+φ′ρ=0, (2)
Thetworegimesarelinkedbynotingthatundercertain 2
approximations the fluid variables and their equations of φ′′+ φ′ =4πGρ. (3)
r
motion in thetwo regimes coincide nearthe surface.
The displacement vector of the (polar) nonspherical
perturbation is expanded in spherical harmonics as
In Sec. 2, we derive the combination of linear non-
sphericalandplanarnon-linearfluidmotion.Wesummarise ξ~(~r,t) e−iωt ξr(r)Ylm(θ,ϕ)~er+ξh(r)r~⊥Ylm(θ,ϕ) ,(4)
≡ ∇
the(well-known)linearperturbationequationsforadiabatic
where (cid:2) (cid:3)
non-spherical stellar oscillations in a suitable notation, and
1 ∂ 1 ∂
theirlimitnearthesurfaceifthedensityvanishesthere.We ~ = ~e + ~e . (5)
∇⊥ r θ∂θ rsinθ ϕ∂ϕ
present the hodograph transform and the shock formation
criterion. We cast the hodograph equation in a form where Because theequations are linear, it is customary to makeξ~
it can be identified with the linear perturbation equations complex as above for ease of calcuation. The physical dis-
nearthesurfaceinahigh-frequencyCowlingapproximation. placement is its real part
Wethenevaluatetheshockformationcriteriononsolutions
ξ~ =Reξ~. (6)
ofthelinearperturbationequationsasiftheyweresolutions real
of the hodograph equation. Weneglectaxialdisplacements,whichinanon-rotatingstar
havenorestoringforce,andarezeromodes.The(real)fluid
In Sec. 3, we apply this general formalism to waves velocity is simply the time derivative of the displacement,
raised in a star by the tidal force of its binary companion. or
We can then use standard methods to calculate the reac-
~v (~r,t)=Re[ iωξ~(~r,t)]. (7)
tion to this force by expandingthe perturbations in proper real −
oscillation modes. We obtain the orbital separation dcrit at TheLagrangianperturbationδf ofanybackgroundquantity
which shocksfirstform as functionof themodes ofthestar f(r) is related to theEulerian perturbation f by
1
and the mass ratio q.
δf f +ξ~ ~f =f +ξ f′. (8)
1 1 r
≡ ·∇
In Sec. 4, we carry out the necessary numerical mode Allscalar perturbationsarealso expandedin spherical har-
calculations for stars with polytropic equations of state. monics, for example theEulerian densityperturbation
Sec. 5 reviews our main approximations and states our as-
ρ (~r,t)=eiωtρ (r)Y (θ,ϕ). (9)
trophyical conclusions. 1 1 lm
The assumption of adiabatic perturbations is that
Asimilarcalculation toourSec.3hasbeencarried out δP =c2δρ. (10)
for g-modes in NSs by Lai (1994). Their frequency is lower
thantheorbitalfrequencyatmerger,andsotheorbitalfre- For given spherical harmonic index l, we define the short-
quency moves through resonance as the orbit shrinks, and hand
the full time-dependent driven oscillator problem must be l(l+1)c2
F(r) 1 . (11)
considered. It was assumed that no shock forms, and dissi- ≡ − ω2r2
pative heating was estimated instead. It turns out that the
The radial and (polar) horizontal parts of the Euler
duration of the resonance is too short to give rise to sig-
equation and the mass conservation equation can be com-
nificant heating. By contrast, we focus on p-modes, which
binedtogiveasinglesecond-orderODEforξ .Withalater
havehigher frequencies and are neverin resonance, and es- r
approximation in mind, we write thisas
timate their amplitude adiabatically in the approximation
of a stationary circular orbit. ξ′′+Aξ′ +Bξ =S, (12)
r r r
(cid:13)c 0000RAS,MNRAS000,000–000
Shock formation in stellar perturbations and tidal shock waves in binaries 3
where P(r)=K0ρ(r)1+n10 K(r)=K0ρn10−n1, (22)
⇔
c′ F′ N2 2 φ′
A = 2 + , (13) withn0 >nrequiredforstability.(Notethatncharacterises
c − F − g r − c2 the EOS, n/n characterises the stratification, and n ap-
0 0
c′ F′ N2 2 φ′ pearsintheLane-Emdenequation).Withρ xn0 asx 0,
B = 2 ∼ →
c − F − φ′ r − c2 we have
(cid:18) (cid:19)(cid:18) (cid:19)
2 φ′ ′ F K(r) x1−nn0, (23)
+ (ω2 N2), (14) ∼
−r2 −(cid:18)c2(cid:19) c2 − and so K, K′ and N all diverge as x 0. Clearly this
→
φ′ l(l+1) N2 c′ F′ 2 assumptionisunphysicalonsufficientlysmallscales,andso
S = 1 2 + + φ . (15)
c2 − ω2r2 φ′ − c F r 1 near the surface the stratification must be adjusted away
(cid:18) (cid:19) from a purepower law.
This is complemented by theperturbed Poisson equation To make an alternative quantitative assumption, we
2 l(l+1) merely assume that N2 is finiteat the surface, or
φ′′+ φ′ φ =4πGρ . (16)
1 r 1− r2 1 1 N2 =O(1). (24)
Itssource term is theEulerian density perturbation
In this approximation, the condition of hydrostatic equilib-
ρ 2 φ′ N2 l(l+1) rium gives
ρ = ξ′ + F ξ φ . (17)
1 −F r r − c2 − φ′ r− ω2r2 1 gx
(cid:20) (cid:18) (cid:19) (cid:21) c2 = +O(x2), (25)
− n
Finally,thecompleteperturbationcanbereconstructedus-
where
ing
ξ (r)= 1 c2 ξ′ + 2 φ′ ξ +φ . (18) g≡φ′(R) (26)
h Fω2r − r r − c2 r 1 isthegravitationalaccelerationatthesurface.Itisalsoclear
(cid:26) (cid:20) (cid:18) (cid:19) (cid:21) (cid:27)
that
Our second-order equations for ξ and φ can be de-
r 1
rived from the first-order systems given in, for example, φ =O(1), φ′ =O(1). (27)
1 1
Unnoet al. (1989) and Christensen-Dalsgaard & Mullan
Weshall also need
(1994).
2
φ′′= g+O(x), (28)
−r
2.2 Expansion near the surface and boundary whichfollowsdirectlyfromthebackgroundPoissonequation
conditions assuming n>1. In this limit we have
The boundary conditions for φ are φ rl at r = 0 and n+1
φ′1+(l+1)φ1/r =0 at r=R.1The bo1un∼dary condition for A = x +O(1), (29)
ξr at r = 0 is ξr ∼ rl−1 (Unnoet al. (1989)). To find the B = n ω2 1 2(n+1) + l(l+1) +O(1)
boundary condition for ξr at r=R, we need to expand the −gx − nσ2 nσ4
equationstoleadingorderinx r R60.Inthefollowing, (cid:20) (cid:21)
n
O(x) will be shorthand for O(≡x/R−). ω˜2+O(1), (30)
| | ≡ −gx
We assume that near the surface the EOS is approxi-
n l(l+1)φ (R)
mated bythe Gamma-law EOS S = φ′(R) 1 +O(1)
−gx 1 − nσ2 R
ρe (cid:20) (cid:21)
P(ρ,e)= , (19) n
n S˜+O(1), (31)
≡ −gx
whereP isthepressure,ρisthemassdensity,etheinternal
energy per rest mass, and n a constant. From the first law where
of thermodynamics, thisis equivalentto R R3
σ2 ω2 =ω2 , (32)
P(ρ,s)=K(s)ρ1+n1, (20) ≡ g GM
is a dimensionless mode frequency. Note that ω˜ and S˜ are
wheresistheentropyperrestmass.TheformofK(s)does
definedas constants.
not matter for our purposes (additional input is required Keeping only the leading O(x−1) in A, B and S, (12)
to fix it), and more generally K can be considered a func-
near thesurface becomes
tion of both entropy and composition to take into account
stratification effects in NSs. d2ξr + n+1dξr n ω˜2ξ = n S˜. (33)
For a stratified background star where K(s) is a given dx2 x dx − gx r −gx
function K(r), we find The solution of (33) that is regular at x=0 is
N2 = n K′. (21) 4nω˜2x S˜
− φ′ −n+1 K ξr(x)=Cfn − g + ω˜2, (34)
(cid:18)r (cid:19)
Intheliterature(forexampleShapiro & Teukolsky(2004)),
where we havedefined thefunction
a power-law stratification is often considered which is pa-
rameterised byspecifying theequilibrium pressure as f (z)=2nΓ(n+1) z−nJ (z). (35)
n n
(cid:13)c 0000RAS,MNRAS000,000–000
4 Carsten Gundlach and Jeremiah W. Murphy
(This is a regular even function of z, normalised to obey 2.4 Matching the two approximations
f (0) = 1.) In the full set of equations S˜ is of course not
n We now have two sets of approximation: in the “pertur-
given a priori but is itself proportional to ξ via (15) and
r bation approximation”, everything is linearised around a
(16).Wealso need to fixan overall factor in themode, and
spherical equilibrium solution. In the “hodograph approx-
we choose to makethemode ξ dimensionless and set
imation”, the vertical fluid motion is treated in full nonlin-
ξ (R)=1. (36) earity,butwe neglect horizontal motion, entropygradients,
r
andangulardependence,andapproximatethegravitational
Then C has a definite value (for any given mode and poly-
field as fixed and constant in space and time.
tropicindexn),whichisdeterminedbysolvingthefullequa-
Weexpectthatthereisanoverlapregionjustbelowthe
tions (12) and (16). (In the Cowling approximation, where
surfacewherebothsetsofapproximationsholdat thesame
φ 0, we would haveC =1.)
1 ≡ time.Inthatregion,weshouldthenfindthesameequation
Theregularsolution(34)obeystheboundarycondition
ofmotion.Toseethis,notethattheperturbationequations
nω˜2 S˜ canbeadaptedtoplane-parallel motionbyformally setting
ξr′ − (n+1)g ξr− ω˜2 =0. (37) l = 0 and 1/R = 0 (and hence F = 1 and σ = 0), and
(cid:18) (cid:19) to a constant gravitational field by setting φ′(r) = g and
This boundary condition is equivalent to Eq. (17.69) φ =0. Neglectingentropygradientscorresponds tosetting
1
of Cox (1980), and the boundary conditions derived in N2 = 0 and c2 = ngx (with g constant). With all these
−
Christensen-Dalsgaard & Mullan (1994) but is not equiva- approximations, (12) reducesto
lenttoδP/P =0.Thislatterboundaryconditionisderived
d2ξ n+1dξ n
in Unnoet al. (1989) under the assumption of finite sound r + r ωξ =0. (44)
speed at the surface, see their Eq. (18.31), and is therefore dx2 x dx − gx r
not applicable here. We now work from the other side. Consider a real, λ-
We introduce the dimensionless radius and mode fre- periodic solution of (41) of theform
quency
v(λ,µ)=Re( iω)e−iωt¯ξ¯(x¯), (45)
r R3 −
s≡ R, σα2 ≡ GMωα2. (38) where we havedefined
(Later, when we consider binaries, R and M will refer to µ2 λ
x¯ , t¯ . (46)
R1 andM1.)Wecanthenwritetheapproximationnearthe ≡−4ng ≡ g
surface as
Then ξ¯(x¯) obeys
S˜
ξr(s)=Cfn 4nσ˜2(s−1) + ω˜2. (39) d2ξ¯+ n+1dξ¯ n ω2ξ¯=0, (47)
(cid:16)p (cid:17) x¯2 x¯ dx¯ − gx¯
which is of course formally the same equation as (44), al-
2.3 Nonlinear isentropic perturbations in the
though it represents nonlinear physics. Consider now a so-
constant gravitational field, plane-parallel
lution of (41) that represents a small perturbation about
approximation
thehydrostaticequilibrium solution (25) with v=0, in the
In Gundlach & Please (2009), we considered nonlinear sense that
smooth adiabatic motions in the approximations of planar
v c, δc c. (48)
geometry, thebarotropic EOS | |≪ | |≪
P =Kρ1+n1, (40) Then from thedefinitions (42) and (46) we can infer that
x¯ x, t¯ t. (49)
with K constant,and constant gravitational acceleration g, ≃ ≃
and derived the linear partial differential equation Furthermore, identifying the planar velocity v with the ra-
dialvelocityv ,comparing(45)with(7),andusing(49),we
2n+1 r
vλλ =vµµ+ µ vµ, (41) have
ξ¯ ξ . (50)
in theindependentvariables r
≃
Wehavenow justified thecoincidence of (44) and (47).
µ 2nc, λ v+gt. (42)
≡ ≡ However,theactuallimitoftheperturbationequations
Here suffices denote partial derivatives. The surface is now near the surface is not (44) but (33). They differ in that ω˜
atµ=0,andtheinteriorofthestaratµ>0.Theboundary is not ω,and bythe(constant in x) source term S˜ which is
condition vµ =0 at µ=0 selects theregular solution. not present at all in thehodograph approximation. Tracing
Thecriterionforthenonlinearfluidequationstoforma thedifferencesback totheEuler andPoisson equations, we
shockisthatthetransformationfrom(x,t)to(µ,λ)becomes see that the terms in (30,31) proportional to l(l+1) arise
singular. Weshowed that thisis equivalentto from horizontal motion, and the middle term in (30) arises
(1 v )2 v2 <0. (43) from the spherical (rather than planar) symmetry of the
− λ − µ background. The difference between ω˜ and ω vanishes in
(Obviously,vmustberealinthisformula.)Ifandonlyifthis thehigh-frequency limit σ2 1.
≫
conditionisobeyed,ashockhasformed,andthesolutionof Generally, the source term S in (12) represents the ef-
(41) no longer has physical significance. fect of the perturbed gravitational potential on the fluid
(cid:13)c 0000RAS,MNRAS000,000–000
Shock formation in stellar perturbations and tidal shock waves in binaries 5
displacement.Inthehodographapproximation,suchaterm Introducingthenew shorthands
cannot be accounted for because the mathematics require
σ2
the gravitational field g to be constant. However, the part κ α, (59)
S˜/ω˜2 ofthenear-surfaceapproximation(39)tothemodeξr α ≡ σf2
is constant in space, and so corresponds to the whole near- ψ(n,κ,y) f (y)2+κf′(y)2, (60)
n n
surface region bobbing up and down as one. Clearly, this ≡
Ψ(n,κ) maxψ(n,κ,y), (61)
part of motion has no effect on shock formation. We will ≡ py
therefore identify ξ¯with ξ minus its constant-in-xpart.
r we can write out the shock formation criterion for a single
(We note in passing that in the high-frequency ap-
mode as
proximation S˜/ω˜2 φ′(R)/ω2. If the mode oscillates with
≃ 1
its own proper frequency, the corresponding displacement A C σ2Ψ(n,κ )max Y >1, (62)
φ′(R)/ω2cosωt is precisely what results from the gravita- α α f α θ | lm|
1
tional field φ′(R)cosωt.) Analysis of the function f shows that for κ < 2(n+1),
− 1 n
In summary, to piece together the two approximations ψ(y) has its global maximum at y = 0, so Ψ = 1. For κ >
into a single “nonlinear mode equation”, we solve the stan- 2(n+1),theglobalmaximumisattainedatthepointwhere
dard linear perturbation equations for ξr and φ1 on the fn+κfn′′ =0, with value Ψ>1. In particular, if the mode
whole domain 0 6 r 6 R, and then consider only the Cfn oscillates at its proper frequency, ωf = ωα, we have κ = 1
part of ξr when we evaluatethe shock formation criterion. and henceΨ=1.
2.5 Evaluating the shock formation criterion for a 2.6 Evaluating the shock formation criterion for a
single mode with periodic time dependence single mode with arbitrary time dependence
Consider now a mode ξ~α with proper frequency ωα that Consider now a mode ξ~α driven with an arbitrary time-
is driven with another frequency ωf, resulting in some di- dependentamplitude Ξα(t),that is
mensionless amplitude A . (To consider a mode oscillating
α
freely, we can just set ωf =ωα in what follows.) Hence the Ξr(r,θ,ϕ,t)=RΞα(t) Ylm(θ,ϕ) ξα(r). (63)
actual time-dependent displacement radial displacement is
Hence,near thesurface and in thebobbing frame
given by
λ ω
Ξr(r,θ,ϕ,t)=AαR cos(ωft) Ylm(θ,ϕ) ξα(r). (51) v=RCαΞ˙α g Ylm(θ,ϕ) fn gαµ . (64)
(cid:18) (cid:19) (cid:18) (cid:19)
Near thesurface thisapproximates as
Hencewe can write in dimensionless form
Ξ (r,θ,ϕ,t) A Rcos(ω t) Y (θ,ϕ)
r ≃ α f lm vλ =p′(z)q(y), vµ =p(z)q′(y), (65)
4nω˜2x S˜
C f α + . (52) where
α n − g ω˜2
" r ! α# z
p(z) ω−1Ξ˙ , (66)
From the identification we have discussed, going into the ≡ α α ωα
bobbingframe we obtain q(y) C Y σ(cid:16)2f (cid:17)(y). (67)
≡ α lm α n
ω ω Wethenhavetominimise(1 v )2 v2 overbothz andy,
v(λ,µ)=−AαCαRωfYlmsin gfλ fn gαµ . (53) to see if it reaches a negative−valλue.− µ
(cid:18) (cid:19) (cid:18) (cid:19)
Note that the periodic case of the previous subsection
Herewehave,somewhat arbitrarily, chosen touseω asan
α is recovered with p(z)=sinτ,with τ =√κz.
approximation toω˜ ,andφ′(R)asan anapproximationto
α 1
S˜.Inthisformula,weconsidervasaslowlyvaryingfunction
of the angles.
Substituting (53) into the shock formation criterion 3 PERTURBATIONS RAISED BY TIDAL
(43), and first focussing on the λ-dependence,we can write FORCES
theresult as
3.1 Calculation of the tidal acceleration
(1 Ucosτ)2 (V sinτ)2>0, (54)
− − Consider a binary system with masses M and M in an
1 2
where elliptic orbit. Let theorbital angular velocity and the spins
y ωαµ, τ ωfλ, (55) of the stars with respect to an inertial reference system be
≡ g ≡ g Ω~, S~1 and S~2. In the (non-inertial) reference system that
U A C Y σ2f(y), (56) movesand spinswith star 1 and with origin in its centreof
≡ − α α lm f mass, theEuler equation for star 1 becomes
V A C Y σ σ f′(y). (57)
α α lm f α
≡
∂~v ~P
Wcoseτfin=dUth/a(tUt2h+iscVri2t)e,riaonndisheshnacerptehsetpcreirtieordioicnalolyveirnaτtwlehaesnt ∂t +(~v·∇~)~v+ ∇ρ +∇~φ1=~a, (68)
a full oscillation period is equivalent to where
U2+V2>1. (58) ~a ~φ +~¨r 2S~ ~v S~˙ ~r S~ S~ ~r. (69)
2 0 1 1 1 1
≡−∇ − × − × − × ×
(cid:13)c 0000RAS,MNRAS000,000–000
6 Carsten Gundlach and Jeremiah W. Murphy
where~v,ρandP arethefluidvelocity,densityandpressure 3.2 Response of the star
in star 1, φ and φ are the gravitational potentials gener-
1 2 The response of the star to~a is governed byLai (1994)
ated by star 1 and star 2, respectively,~r (t) is the location
0
ofthecentreofmassofthebinary,and~e0(t)theunitvector ρ∂2 + Ξ~(~r,t)=ρ~a, (76)
in its direction. Using angular momentum conservation in ∂t2 L
theform d2Ω~ =const, we find (cid:18) (cid:19)
1
where isalineardifferentialoperatorcontainingonlyspa-
~¨r0 =d¨1~e0+Ω~ Ω~ ~r0. (70) tial derLivatives.
× ×
It is generally assumed (for example, Unnoet al.
In the following we assume S~ = 0, both because this
1 (1989)) that a star admits a complete set of eigenmodes
is believed to be correct for NS binaries (Bildsten & Cutler
which obey
(1992); Kochanek (1992)) and because it simplifies the cal-
culation, as we can use perturbation theory on a spherical ξ~ (~r)=ρω2ξ~ (~r) (77)
L α α α
backgroundstar1.Thedistancebetweenthecentreofmass
and
ofstar1andthecentreofmass ofthebinaryis ~r (t) d
0 1
| |≡
andthedistancebetweenthecentresofmassofthetwostars ξ~ ,ξ~ δ , (78)
α β αβ
is ~r (t) d d +d . Here M d =M d by virtue of~r h i∝
2 1 2 1 1 2 2 0
be|ing th|e≡cen≡tre of mass. where theinner product is
We approximate φ2 as spherically symmetric, that is 1
ξ~,~η ξ~ ~η∗ ρr2drdΩ (79)
GM h i≡ M ·
φ = 2 , (71) Z
2 − ~r ~r2(t) 1 R
| − | = ξr ηr∗ +l(l+1)ξh ηr∗ ρr2dr. (80)
and expand ~φ2 up to O(r2) in~r. Then M l,m Z0 lm lm lm lm
∇ X (cid:2) (cid:3)
GM In the second line, we have assumed a decomposition into
~a = d¨ + 2 d Ω2 ~e
1 d2 − 1 0 spherical harmonics, and into the radial and horizontal
(cid:16)GM GM (cid:17) parts defined by Eq. (4), and the standard normalisation
− d32~r+3 d32~rk0+O(r2), (72) Y 2dΩ = 1. For simplicity of notation, and because~a is
| |
where Ω denotes the projection into the plane normal to Rpolar, we neglect the axial parts of all vectors. Note that ξ~
Ω~ and⊥0 denotes the projection into the direction ~e . The isbyourconventiondimensionless,Ξ~ hasdimensionlength,
0
O(r0) tkerm in round brackets vanishes by the assumption and~a has dimension acceleration.
that the origin of ~r is the centre of mass of star 1. The Decomposing (76) into modes, we have
remainder can bewritten as
Ξ~(~r,t)=R Ξ (t)ξ~ (~r), (81)
α α
~a=∇~χ+O(r2) (73) Xα
where thedimensionless amplitudes Ξ (t) obey
with α
χ≡−G2Md32r2+ 3G2dM32(xcosφ+ysinφ)2, (74) ddt22 +ωα2 Ξα = R1 hξ~ξ~α,,ξ~~ai ≡fα(t). (82)
(cid:18) (cid:19) h α αi
wherewehavechosenCartesiancoordinatessothattheorbit The general solution is
is in the xy-plane, and where φ(t) Ωdt is the orbital
pd(hta)steh,eΩi(nts)tathnetainnesotaunstoarnbeiotualssoerpbairtaa≡tlioaRnn.guWlaercvaenlowcirtiyt,eatnhde Ξα(t)=Bαe−iωαt+ ω1α Zt0tsinωα(t−s)fα(s)ds. (83)
tidal potential in terms of spherical harmonics as
GM 3GM
χ = 2r2+ 2r2sin2θ[1+cos2(φ ϕ)] 3.3 Circular orbits
− 2d3 4d3 −
3GM 1 For a circular orbit, only the l = m = 2 component of
= 2r2 sin2θ cos2(φ ϕ) cos2θ
4d3 − − − 3 the tidal potential is time-dependent, through φ(t) = 2Ωt,
h (cid:16) (cid:17)i whereΩis theconstant orbitalangular velocity.Neglecting
3GM 2π
= 2 r2Re e−2iφY thetime-independentparts,andsplittingthe3-dimensional
d3 15 22
r gradient intoa horizontal and radial part, we have
GM π (cid:0) (cid:1)
2 r2Y . (75)
− d3 5 20 ~a= GM2 6πRee−2iΩt 2rY ~e +r2~ Y . (84)
q d3 5 22 r ∇⊥ 22
Hence the tidal force is polar, and to leading order in ~r r
is given by quadrupole terms. To linear order in perturba- Hencewe have (cid:0) (cid:1)
tion theory, the deformations caused by the Y , Y and
00 20 GM 6π ξ~ ,(2s,s)
Yin2d2etpeernmdsendtecboeucpalue.seIndciarncdulaΩr oarrbeitcso,ntshteanYt2,0atnedrmthiserteifmoree- f22(t)= d32r 5 e−2iΩt h αξ~α,ξ~α i. (85)
h i
doesnotcauseshocks.Wecanneglectitalsoformoderately
The solution of (76) is then
eccentricorbits.TheY termisalwaystime-dependentbe-
22
cause φ is: physically, the tides rotate around the star, so Ξ~(~r,t)=R B e−iωαt+A e−2iΩt ξ~ (~r), (86)
α α α
that individual fluid elements moveup and down.
Xα (cid:0) (cid:1)
(cid:13)c 0000RAS,MNRAS000,000–000
Shock formation in stellar perturbations and tidal shock waves in binaries 7
wherethedimensionlessamplitudeoftheparticularintegral haveapproximatelyΞ (t) ω−2f (t),asforthecircularor-
α α
≃
is bit case. Hence we expect even the highly eccentric case to
excitethe22,20and00perturbationsatsimilaramplitudes
1 GM 6π ξ~ ,(2s,s)
Aα= ωα2 −4Ω2 d32r 5 h hαξ~α,ξ~αi i. (87) rΞaαd(itu)s,.comparable to Ξ22 in a circular orbit at periastron
Notethatincontrast toLai(1994),weapproximatetheor-
bit as circular, and we do not take into account resonance.
3.5 Roche lobe overflow and resonance
We assume that the tidal force is the dominant source of
oscillations, and that the orbit evolves so slowly that tran- Our shock formation criterion becomes irrelevant once
sients can beneglected, and so we set B =0. Roche lobe overflow occurs (or when two stars of identical
α
Weintroducethebinarymassratioandthedimension- mass and size touch). Defining
less orbital separation
d
λ (97)
q M2, η d . (88) ≡ d11
≡ M1 ≡ R1 where d is the distance from the centre of star 1 to the
11
Wenote Lagrange point L , a simple calculation gives thequintic
1
Ω2 = G(1+q)M1. (89) λ5 2λ4+λ3 (1+3q)λ2+(2+3q)λ (1+q)=0. (98)
circ d3 − − −
circ
The one real solution gives L (the four complex solutions
1
Wecan then write
give the other Lagrange points in the complex plane). Ap-
6π q ξ~ ,(2s,s) proximating the Roche lobe as a sphere centered on star 1,
Aα=r 5 σα2η3−4(1+q) h hαξ~α,ξ~αi i. (90) RochReelsoobneaonvceeriflsoowbtoacicnuerdsaattηan=oλrb.iAtaslqse≫pa1ra,tλio≃n (3q)1/3.
Combiningthiswith(62),andusingmax Y = 15/32π
and ω = 2Ω, we have the dimensionlesθs|sh22o|ck formation ηres =σ−2/3(4+4q)1/3. (99)
f
p
criterion Ourresult(91)hasbeenobtainedundertheassumptionthat
1+q q ξ~ ,(2s,s) ωα >2Ω, and this is always the case until shock formation
3CαΨ(n,κα) η3 σα2η3−4(1+q) h hαξ~α,ξ~αi i >1, (91) or Rochelobe overflow.
where α characterises the l = m = 2 mode that maximises
thecriterion, and where 4 MODE RESULTS FOR POLYTROPIC STARS
κ = σα2η3 . (92) We have used a publicly available code
α 4(1+q) (Christensen-Dalsgaard(1997);Christensen-Dalsgaard et al
(1996)) to calculate mode functions ξ(s) and frequencies
σ and hence determine the overlap integrals and critical
3.4 Elliptic orbits
valuesof η as a function of q.
For elliptic orbits, we have We have carried out the calculation for the isentropic
(n =n) Gamma-law EOS
GM 6π ξ~ ,(2s,s) 0
f22(t) = d3(t2) 5 e−2iφ(t) h αξ~ ,ξ~ i, (93) P =Kρ1+n1 (100)
r α α
h i
f20(t) = − 4Ω32(t) + Gd3M(t2) π5 hξ~αξ~,(,2ξ~s,s)i, (94) awpitphronxim=at1e/2e,qu1a,ti3o/n2saonfdst3a.tenfo=r c1olids oNftSenmuatsteedr.asna=n
(cid:18) (cid:19)q h α αi 1/2 represents possible stiff neutron star equations
4Ω2(t) ξ~ ,(2s,s)
f00(t) = − 3 √π h αξ~α,ξ~α i. (95) onofns-traeltaet.ivnist=ic3a/n2darenldatnivi=stic3daergeengoeroadteaeplpercotrxoimnaptrieosnssurteo,
h i
For slightly elliptic orbits, the time-dependence of Ω and d respectively. While these equations of state are simplistic,
is weak, and l = m = 2 is the dominant time-dependent they have the advantage that the resulting stellar models,
deformation. For highly eccentric orbits, both Ω2 and d−3 andhenceourresults,dependonlyonn,notonM1 andR1.
aresharplypeakedatperiastron.Thishastwoconsequences: Themass,radiusandpolytropicconstantarerelatedbythe
First, depending on how quickly perturbations set up by scaling relation (Shapiro& Teukolsky (2004))
tidalforcesaredamped,itmaybeappropriatetotreateach Rn−3 KnMn−1. (101)
periastron passage as a transient, rather than as part of 1 ∝ 1
periodic excitation. Secondly, f , f and f are now all In these simple stellar models there is no stratification and
22 20 00
comparably time-dependent.In fact, hencethere are nog-modes.
Ournumericalresultsforcircularorbitswithequations
G(1+q)M
Ω2 =(1+e) 1, (96) of state n=1/2, 1, 3/2 and 3 are shown in the figures. For
pa d3
pa the first several modes, we plot the critical value of orbital
where e is the eccentricity and d the periastron orbital separation η against mass ratio. In all cases, shock forma-
pa
separation. The three overlap integrals are also likely to be tion first occurs for the lowest frequency mode, so that is
comparable. Moreover, as the orbital frequency even at pe- the curve that matters for shock formation. On the same
riastron is substantially lower than the mode frequency, we plot, we also show the orbital separation η at which Roche
(cid:13)c 0000RAS,MNRAS000,000–000
8 Carsten Gundlach and Jeremiah W. Murphy
Figure 1. Critical value for shock formation of dimensionless Figure 3. The same plots for the polytropic equation of state
orbitalseparationη=d/R1againstmassratioq=M2/M1.This withn=3/2(Γ=5/3).
appliesforacircularorbit,thepolytropicequation ofstatewith
n = 1/2 (Γ = 3), and l = m = 2 perturbations. The labels
0,1,2,... denote the number of radial nodes in the mode. The
critical value of η forRoche lobe overflow through the Lagrange
point L1 is shown as a dashed line. At large q, η ∝ q1/3 for all
curves.
Figure 4. The same plots for the polytropic equation of state
with n = 3 (Γ = 4/3). (The bottom right corner appears blank
becausewehavenotcalculatedmodesthere.Allourfiguresshow
thesamerangesofη andq.)
Figure 2. The same plots for the polytropic equation of state
anybuoyancy(non-barotropic)effectscanalsobesafelyne-
withn=1(Γ=2).
glectednearthesurfaceaslongastheentropygradientand
any composition gradients are merely bounded.
lobe overflow starts. For all three equations of state, for all Our calculation of the tidal waves in perturbation
values of q, the critical orbital separation for the formation theory is straightforward for circular orbits (Lai (1994)).
of tidal shocks coincides, within our approximations, with For stars with a simple polytropic equation of state P =
thecritical orbital separation for Roche lobe overflow. Kρ1+1/n in irrotational circular binary orbit, we find that
thecriticalorbitalseparationforshockformationessentially
coincides with the one for Roche lobe overflow. In other
words, tidal forces create shocks roughly when the
5 CONCLUSIONS
binarybeginstomerge.Withinourapproximations,
This paper consists of two parts: a quantitative criterion the two curves agree remarkably closely, so that we
for linear perturbation modes to form shocks, and an ap- cannot say which actually occurs first. In any case,
plication of this criterion to modes raised by tidal forces in the p-mode shock formation mechanism we have in-
compact binaries. vestigated here is not the primary mechanism for
Our shock formation criterion relies on the hodograph binary disruption. As discussed above in Sec. 3.4,
transformation to link linear perturbation modes in the in- weexpectthesameresultevenforhighlyellipticorbits.Al-
terior to the fully nonlinear shock formation criterion of though this is a negative result, it should be stressed that
Gundlach & Please (2009) near the surface. This criterion it was not obvious from dimensional analysis: we have also
is exact for plane-symmetricmotion of a polytropic fluidin estimated thedimensionless factors.
a constant gravitational field. The approximation of planar Extending our analysis to more realistic stellar
symmetry near thesurface is natural, and it turnsout that models would require more extensive modelling, in
(cid:13)c 0000RAS,MNRAS000,000–000
Shock formation in stellar perturbations and tidal shock waves in binaries 9
particular a realistic treatment of the surface. (We Rosswog S., Ramirez-Ruiz E., Hix W. R., 2009, AJ, 695,
have shown above in Sec. 2.2 that another simple 404.
stellar model, assuming one polytropic constant for Shapiro S. L., Teukolsky S. A., 2004, Black holes, white
the equation of state and another for the stellar dwarfs and neutron stars, Wiley, Weinheim.
structure, gives rise to a divergent Brunt-V¨ais¨al¨a Stergioulas N., Apostolatos T. A., Font J. A., 2006, MN-
frequency, and so is inconsistent with our assump- RAS,368, 1609.
tions of a perfect fluid surface). However, the fact Unno W. et al., 1989, Nonradial oscillations of stars, 2nd
that our result holds for a polytropic index ranging edition, University of TokyoPress, Tokyo.
fromn=1/2ton=3suggeststhatotherequationsof
statewouldnotshowshockformationbeforemerger
either. Intuitively,the weaknessof the shock forma-
tion mechanism is dominated by a factor of (tidal
force frequency/mode frequency)2 [the factor σ2 in
f
Eq. (62)].
In a related result, Rosswog et al. (2009) give a cri-
terion for the tidal disruption of a WD in orbit around a
much more massive compact object 2 (black hole or NS) as
(in our notation) η q1/3 based on numerical simu-
disrupt
≃
lations. Our results are consistent with this for large q in
that both Roche lobe overflow and tidal shock for-
mation occur at the same η, namely η 1.44q1/3
crit
≃
for q 1.
≫
Finally, this paper is motivated by the obser-
vance of surface shocks in numerical simulations of
binary neutron star mergers just before the stars
touch (Baiotti et al. (2008); Rezzolla et al. (2010)).
Our results indicate that these surface shocks are
probably not physical.
ACKNOWLEDGMENTS
J.W.M. is supported by an NSF Astronomy and Astro-
physicsPostdoctoralFellowshipunderawardAST-0802315.
C.G.wouldliketothankRandyLevequeandtheUniversity
ofWashingtonforhospitalitywhenthisworkwasbegun.We
liketothankEricAgol,WynnHo,TomMaccaroneandJohn
Miller for helpfuldiscussions.
REFERENCES
BaiottiL.,GiacomazziB.,RezzollaL.,2008,Phys.Rev.D,
78, 084033.
Bildsten L., Cutler C., 1992, AJ, 400, 175.
Christensen-Dalsgaard J., 1997, code from
http://astro.phys.au.dk/ jcd/adipack.n/ (accessed
∼
in 2010)
Christensen-Dalsgaard J., Mullan D. J., 1994, MNRAS,
270, 921.
Christensen-Dalsgaard J. et al, 1996, Sci, 272, 1286.
CoxJ.P.,1980,Theoryofstellarperturbations,Princeton
UniversityPress, Princeton.
Dimmelmeyer H., Stergioulas N., Font J. A., 2004, MN-
RAS,352, 1089.
Gabler M., Sperhake U., Andersson N., 2009, Phys. Rev.
D,80, 064012.
GundlachC.,PleaseC.P.,2009,Phys.Rev.D,79,067501.
GundlachC.,LeVequeR.J.,2011,J.Fluid.Mech.,inpress.
Kochanek C. S., 1992, AJ, 398, 234.
Lai D., 1994, MNRAS270, 611.
RezzollaL.etal,2010, Class.QuantumGrav.,27,114105.
(cid:13)c 0000RAS,MNRAS000,000–000