MNRAS
000
, 1–21 (2020)
Preprint 8 May 2020
Compiled using MNRAS L
A
T
E
X style file v3.0
Nonlinear dynamical tides in white dwarf binaries
Hang Yu
1
,
2
?
, Nevin N. Weinberg
1
and Jim Fuller
2
1
Department of Physics, and MIT Kavli Institute, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
2
TAPIR, Walter Burke Institute for Theoretical Physics, Mailcode 350-17 California Institute of Technology, Pasadena, CA 91125
ABSTRACT
Compact white dwarf (WD) binaries are important sources for space-based gravitational-wave
(GW) observatories, and an increasing number of them are being identified by surveys like
ZTF. We study the effects of nonlinear dynamical tides in such binaries. We focus on the
global three-mode parametric instability and show that it has a much lower threshold energy
than the local wave-breaking condition studied previously. By integrating networks of coupled
modes, we calculate the tidal dissipation rate as a function of orbital period. We construct
phenomenological models that match these numerical results and use them to evaluate the spin
and luminosity evolution of a WD binary. While in linear theory the WD’s spin frequency can
lock to the orbital frequency, we find that such a lock cannot be maintained when nonlinear
effects are taken into account. Instead, as the orbit decays, the spin and orbit go in and out
of synchronization. Each time they go out of synchronization, there is a brief but significant
dip in the tidal heating rate. While most WDs in compact binaries should have luminosities
that are similar to previous traveling-wave estimates, a few percent should be about ten times
dimmer because they reside in heating rate dips. This offers a potential explanation for the
low luminosity of the CO WD in J0651. Lastly, we consider the impact of tides on the GW
signal and show that LISA and TianGO can constrain the WD’s moment of inertia to better
than
1%
for deci-Hz systems.
Key words:
instabilities - white dwarfs - stars: oscillations (including pulsations)- binaries
(including multiple): close - gravitational waves
1 INTRODUCTION
As binary white dwarfs (WDs) with short orbital periods inspi-
ral due to the emission of gravitational waves (GWs), they can
evolve into a variety of interesting systems, including AM CVn
stars (Nelemans et al. 2001), R Cor Bor stars (Clayton 2012),
and rapidly rotating magnetic WDs (Ferrario et al. 2015). Merg-
ing WDs may also explode as type Ia supernovae (Webbink 1984;
Iben & Tutukov 1984; Toonen et al. 2012; Polin et al. 2019b) or
in other types of luminous thermonuclear events (Shen et al. 2018;
Polin et al. 2019a). Compact WD binaries emit GWs with frequen-
cies of
≈
1
−
100
mHz, which makes them prominent sources for
proposed space-based GW observatories such as the Laser Inter-
ferometer Space Antenna (LISA, Amaro-Seoane et al. 2017), Tian-
Qin (Luo et al. 2016), and TianGO (Kuns et al. 2019).
The tidal interaction between the binary components spins
them up and heats their interiors. As they inspiral, the tide becomes
progressively stronger and eventually their spin frequency nearly
equals the orbital frequency. However, they never become perfectly
synchronous because of the continual GW-induced orbital decay.
The degree of spin asynchronicity affects the tidal heating rate and
?
E-mail: hangyu@caltech.edu
luminosity of the WDs (Iben et al. 1998; Fuller & Lai 2012a, 2013;
Piro 2019) and the outcome of their potential merger (Raskin et al.
2012; Dan et al. 2014; Fenn et al. 2016).
The dominant mechanism of tidal dissipation is most likely
the excitation of internal gravity waves, either in the form of stand-
ing waves (i.e., g-modes; Fuller & Lai 2011; Burkart et al. 2013),
or traveling waves (Fuller & Lai 2012a,b, 2013, 2014). As we
will show, for orbital periods between approximately
10
min and
150
min, which describes many of the observed WD binaries, the
resonant g-modes excited by the tide have such large amplitudes
that they cannot be considered small, linear perturbations to the
background star. On the other hand, the amplitudes are not so large
that the modes break due to strong nonlinearities. The tidal dynam-
ics and dissipation in this intermediate, weakly nonlinear regime
are complicated and depend on details of the nonlinear coupling
between g-modes driven directly by the tide and the sea of sec-
ondary modes they excite.
In this Paper, we apply the weakly nonlinear tidal formalism
developed in Weinberg et al. (2012) to study tides in WD binaries.
Our study fills the gap between those that assume the excited modes
are linear standing waves (e.g., Fuller & Lai 2011; Burkart et al.
2013) and those that assume they break and form strongly nonlin-
ear traveling waves (Fuller & Lai 2012a,b, 2013, 2014). In Sec-
tion 2, we present the background WD model we use throughout
c
©
2020 The Authors
arXiv:2005.03058v1 [astro-ph.SR] 6 May 2020
2
H. Yu, N. N. Weinberg and J. Fuller
Table 1.
The mass
M
, radius
R
, effective temperature
T
eff
, and moment
of inertia
I
WD
, of our WD model. We will often express results in terms
of the primary’s natural units of energy
E
0
≡
GM
2
/R
= 1
.
10
×
10
50
erg
and frequency
ω
0
≡
√
GM/R
3
= 0
.
346
rad s
−
1
.
M
R
T
eff
I
WD
0
.
6
M
8
.
75
×
10
8
cm
9000 K
0
.
257
MR
2
much of our analysis. In Section 3, we describe the mode coupling
and tidal driving equations that governs the mode dynamics and
in Section 4 we describe our numerical method for solving these
equations. In Section 5, we present our solutions of the mode dy-
namics and show how tidal dissipation and synchronization varies
with orbital period in the weakly nonlinear regime. We also com-
pare our results with the previous studies that assumed the tide was
either linear or strongly nonlinear. In Section 6, we describe the
observable electromagnetic and GW signatures of the tidal interac-
tion, including the tidal heating luminosities, GW phase shifts, and
projected constraints on the WD moment of inertia. In Section 7,
we summarize our key results and conclude.
2 BACKGROUND MODEL
We use
MESA
(version 10398; Paxton et al. 2011, 2013, 2015,
2018) to construct a WD model, whose key parameters are sum-
marized in Table 1. To construct this model, we adopt parameters
similar to those used by Timmes et al. (2018). Specifically, we start
with a pre-main sequence star with an initial mass of
2
.
8
M
and
metallicity
Z
= 0
.
02
and let it evolve to a CO WD with mass
M
= 0
.
6
M
and effective temperature
T
eff
= 9000
K
. We in-
clude element diffusion, semiconvection, and thermohaline mixing
throughout the evolution. We use
GYRE
(Townsend & Teitler 2013;
Townsend et al. 2018) to compute the model’s eigenmodes and con-
struct our mode networks.
In the upper panel of Figure 1, we show the propagation dia-
gram of our WD model. The solid line is the buoyancy frequency
N
, where
N
2
=
g
2
(
1
c
2
e
−
1
c
2
s
)
,
(1)
c
2
e
=
d
P/
d
ρ
is the equilibrium sound speed squared,
c
2
s
= Γ
1
P/ρ
is the adiabatic sound speed squared, and
Γ
1
is the adiabatic index.
All other quantities have their usual meaning. The dashed line is
the Lamb frequency
S
l
for
l
= 2
, where
S
2
l
=
l
(
l
+ 1)
c
2
s
r
2
.
(2)
For the short-wavelength g-modes that comprise the dynamical
tide, the square of the radial wavenumber
k
2
r
=
ω
2
c
2
s
(
S
2
l
ω
2
−
1
)(
N
2
ω
2
−
1
)
,
(3)
where
ω
is the angular eigenfrequency of the mode. A g-mode
propagates where
k
2
r
>
0
, i.e., in regions where
ω <
N
and
ω < S
l
, and is evanescent where
k
2
r
<
0
.
The lower panel of Figure 1 shows the composition profile of
our model. As is typical of stars supported by degeneracy pressure,
the buoyancy is due largely to composition gradients, with peaks in
N
associated with sharp transitions in the internal composition.
0.1
1
10
100
Frequency [min
−
1
]
N
/
2
π
S
l
= 2
/
2
π
0.30
0.80
0.93
0.97
0.99
r/R
10
15
10
17
10
19
10
21
10
23
P
[
dyn cm
−
2
]
0.01
0.1
1.0
Mass Fraction
1
H
4
He
12
C
16
O
Figure 1.
Propagation diagram (top panel) and composition profile (bottom
panel) of our WD model from the center to the surface. Note that we have
oriented the bottom x-axis such that the radius increases to the right.
3 FORMALISM
3.1 Equation of motion
Consider a primary star of mass
M
and a secondary star of mass
M
′
and choose a coordinate system whose origin is at the cen-
ter of the primary and co-rotates with it. We assume that the orbit
is circular and that the spin angular momentum of the primary is
aligned with the orbital angular momentum. For simplicity, we do
not account for the effect of rotation on the mode dynamics except
through the Doppler shift of the tidal driving frequency. The equa-
tion of motion governing the Lagrangian displacement field
ξ
(
r
,t
)
of a perturbed fluid element at location
r
at time
t
is then (see, e.g.,
Weinberg et al. 2012, hereafter WAQB12)
ρ
̈
ξ
=
f
1
[
ξ
] +
f
2
[
ξ
,
ξ
] +
ρ
a
tide
,
(4)
where
f
1
and
f
2
represent the linear and leading-order nonlinear
internal restoring forces, and
a
tide
=
−∇
U
−
(
ξ
·∇
)
∇
U
(5)
is the tidal acceleration. The tidal potential can be expanded as
U
(
r
,t
) =
−
∑
l
≥
2
,m
W
lm
GM
′
D
(
t
)
[
r
D
(
t
)
]
l
Y
lm
(
θ,φ
)e
−
i
m
(Ω
orb
−
Ω
s
)
t
,
(6)
where
Y
lm
is the spherical harmonic function, and
D
,
Ω
orb
, and
Ω
s
are the orbital separation, the orbital angular frequency, and the
spin frequency of the primary, respectively. We focus on the leading
order quadrupolar (
l
= 2
) tide, whose non-vanishing
W
lm
coeffi-
cients are
W
2
±
2
=
√
3
π/
10
and
W
20
=
−
√
π/
5
. It is useful to
define
=
(
M
′
M
)(
R
D
)
3
=
(
y
1 +
y
)(
Ω
orb
ω
0
)
2
,
(7)
where
y
=
M
′
/M
is the mass ratio. The quantity
character-
izes the overall tidal strength and will be useful when we want to
MNRAS
000
, 1–21 (2020)
Nonlinear tides in white dwarf binaries
3
distinguish the system’s dependence on the tidal strength from its
dependence on the driving frequency
2(Ω
orb
−
Ω
s
)
.
In order to solve Equation (4), we expand the six-dimensional
phase space vector as
[
ξ
(
r
,t
)
̇
ξ
(
r
,t
)
]
=
∑
q
a
(
t
)
[
ξ
a
(
r
)
−
i
ω
a
ξ
a
(
r
)
]
,
(8)
where
q
a
(
t
)
,
ω
a
, and
ξ
a
(
r
)
, are the amplitude, frequency, and dis-
placement of an eigenmode labeled by subscript
a
. The frequency
and displacement are found by solving the linear, homogeneous
equation
f
1
[
ξ
a
] =
−
ρω
2
a
ξ
a
,
(9)
which we normalize as
2
ω
2
a
∫
d
3
rρ
ξ
∗
a
·
ξ
b
=
GM
2
R
δ
ab
≡
E
0
δ
ab
.
(10)
Each eigenmode has a unique set of three quantum numbers: its
angular degree
l
a
, azimuthal order
m
a
, and radial order
n
a
. The
summation in Equation (8) runs over all mode quantum numbers
and both signs of eigenfrequency in order to include each mode and
its complex conjugate
1
. Using the orthogonality of the eigenmodes,
Equation (4) can now be expressed as a set of evolution equations
for the mode amplitudes
̇
q
a
+ (
i
ω
a
+
γ
a
)
q
a
=
i
ω
a
[
U
a
+
∑
b
U
∗
ab
q
∗
b
+
∑
bc
κ
abc
q
∗
b
q
∗
c
]
,
(11)
where
U
a
(
t
) =
−
1
E
0
∫
d
3
rρ
ξ
∗
a
·∇
U,
(12)
U
ab
(
t
) =
−
1
E
0
∫
d
3
rρ
ξ
a
·
(
ξ
b
·∇
)
∇
U,
(13)
κ
abc
=
1
E
0
∫
d
3
r
ξ
a
·
f
2
[
ξ
b
,
ξ
c
]
.
(14)
The linear and nonlinear tidal coefficients
U
a
and
U
ab
characterize
the strength of the coupling of modes to the tide, and the three-
mode coupling coefficient
κ
abc
characterizes the strength of the
coupling of modes to each other.
We further simplify Equation (11) by noting that the three-
mode coupling involving the equilibrium tide cancels significantly
with the nonlinear tide (i.e.,
∑
c
∈
eq
κ
abc
q
∗
c
' −
U
ab
; WAQB12).
We therefore ignore
U
ab
and the equilibrium tide and focus on the
dynamical tide.
2
The latter is dominated by the linear driving of the
most-resonant
l
a
=
|
m
a
|
= 2
modes, for which
|
∆
a
/ω
a
|
1
,
where
∆
a
=
ω
−
ω
a
is the linear detuning and
ω
= 2(Ω
orb
−
Ω
s
)
is the linear driving frequency. We refer to such linearly resonant
1
If the amplitudes
q
a
+
and
q
a
−
correspond to eigenfrequencies
ω
a
and
−
ω
a
, respectively, then the reality of
ξ
requires
q
a
+
=
q
∗
a
−
, where the
asterisk denotes complex-conjugation.
2
We note that the nonlinear driving by equilibrium tide might be unsta-
ble depending on the residual coupling. Roughly, the growth rate for the
equilibrium-tide-driven instability is
Γ
(eq)
nl
∼
κ
(eq)
abc
(Ω
orb
−
Ω
s
)
. If the
residual coupling
κ
(eq)
abc
∼
1
after accounting for the cancellation with
U
ab
,
then we have
Γ
(eq)
nl
∼
10
−
8
s
at
P
orb
= 50 min
. This, while smaller than
the nonlinear growth rate of the dynamical tide [Equation (47)], could be
greater than the damping rate of the resonant
l
= 2
modes at the same
period. We defer the study of this effect to future work.
10
30
100
300
n
a
10
-9
10
-7
10
-5
10
-3
γ/ω
0
γ
diff
γ
turb
α
2
.
4
×
10
−
13
n
2
a
0
.
28(
l
a
/n
2
a
)
0.003
0.01
0.03
ω
a
/ω
0
Figure 2.
Linear damping rate of
l
a
= 2
modes. The blue and green circles
represent the respective contributions of radiative diffusion
γ
diff
and con-
vective turbulence
γ
turb
. The orange circles are the inverse group traveling
time
α
of each mode (see Section 3.5). The solid and dashed black lines
show the WKB scaling relations [Equations (21) and (42)].
modes as parent modes. By contrast, the other modes in our net-
works (the daughters, granddaughters, etc.) are primarily excited
through three-mode parametric resonances rather than direct driv-
ing by the tide since they have large
|
∆
a
|
and smaller
U
a
than the
parents [since they have larger
n
a
and
l
a
; see Equations (23) and
(30)]. In our mode network calculations, we therefore solve a re-
duced set of amplitude equations in which the parent modes
{
a
}
satisfy
̇
q
a
+ (
i
ω
a
+
γ
a
)
q
a
=
i
ω
a
U
a
+
i
ω
a
∑
bc
κ
abc
q
∗
b
q
∗
c
,
(15)
and the daughter modes
{
b,c
}
satisfy
̇
q
b
+ (
i
ω
b
+
γ
b
)
q
b
=
i
ω
b
∑
ac
κ
abc
q
∗
a
q
∗
c
,
(16)
and similarly for the granddaughters, great-granddaughters, etc.
The energy of a mode
E
a
(
t
)
is related to its amplitude by
E
a
(
t
) =
q
∗
a
(
t
)
q
a
(
t
)
E
0
.
(17)
This neglects the energy in the three-mode coupling,
1
3
∑
b,c
κ
abc
q
a
q
b
q
c
+ c
.
c
.,
(18)
where c.c. stands for complex conjugate. As we show in Sec-
tion 3.4, this energy is much less than
∑
a
E
a
and therefore, we
will use Equation (17) to represent the mode energy.
3.2 Power-law relations for the coefficients
In Appendix A we describe our calculations of
γ
a
,
U
a
, and
κ
abc
in detail. For the tidal synchronization problem, we are mostly
interested in binaries with orbital periods in the range
P
orb
=
[10
,
100]
min, which corresponds to
l
a
= 2
parent modes with ra-
dial orders in the range
n
a
'
[10
,
100]
[see Equation (20)]. For
such high-order modes, we find that the coefficients follow simple
power-law relations in
n
a
and
l
a
.
MNRAS
000
, 1–21 (2020)
4
H. Yu, N. N. Weinberg and J. Fuller
We find that the eigenfrequencies of our WD model are ap-
proximately given by
ω
a
'
0
.
28
l
a
n
a
ω
0
= 0
.
098
l
a
n
a
rad s
−
1
,
(19)
i.e., the mode periods are given by
P
a
=
2
π
ω
a
'
1
.
1
n
a
l
a
min
,
(20)
where
ω
0
=
√
GM/R
3
is the dynamical frequency of the WD.
In Figure 2, we show the linear dissipation rates
γ
a
. The dissi-
pation is dominated by electron conduction and radiative diffusion
and for
n
a
l
a
(as is true of all modes in our networks) is approx-
imately given by
γ
a
'
2
.
4
×
10
−
13
n
2
a
ω
0
= 8
.
4
×
10
−
14
n
2
a
s
−
1
.
(21)
By comparison, the dissipation due to turbulent convective damp-
ing (green dots) is much smaller for the modes we are interested in
(see Appendix A for details).
From Equations (6) and (12), we can write the linear tide co-
efficient as
U
a
=
W
lm
Q
a
(
M
′
M
)(
R
D
)
l
+1
e
−
i
m
(Ω
orb
−
Ω
s
)
t
,
(22)
where the overlap integral
Q
a
=
1
MR
l
∫
d
3
rρ
ξ
∗
·∇
(
r
l
Y
lm
)
,
'
2
.
1
n
−
3
.
7
a
δ
l
a
l
δ
m
a
m
.
(23)
In Figure 3, we show
Q
a
, calculated using the method described in
Appendix A3, and the numerical fit above. Note that the overlap is
non-zero only if
l
a
=
l
and
m
a
=
m
.
In Figure 4 we show the three-mode coupling coefficients as
a function of the parent mode’s radial order
n
a
. For high-order
modes, we find
κ
abc
'
41
(
T
0
.
18
)(
n
a
l
a
)
2
,
(24)
where
a
is the parent mode. Here
T
is an angular integral that de-
pends only on each mode’s angular quantum numbers and van-
ishes if the modes do not satisfy the angular selection rules: (i)
|
l
b
−
l
c
| ≤
l
a
≤
l
b
+
l
c
, (ii)
l
a
+
l
b
+
l
c
is even, and (iii)
m
a
+
m
b
+
m
c
= 0
. Otherwise, it is of order unity for the typi-
cal triplets that we consider, e.g.,
T
'
0
.
18
for
l
a
=
l
b
=
l
c
= 2
and
(
m
a
,m
b
,m
c
) = (2
,
−
2
,
0)
. In addition to these angular selec-
tion rules, the modes couple significantly only if their radial orders
satisfy
|
n
b
−
n
c
|
.
n
a
(Wu & Goldreich 2001,WAQB12).
3.3 Nonlinear instability threshold
In the absence of nonlinear interactions, a mode driven by the linear
tide has an energy
E
a,
lin
E
0
=
ω
2
a
∆
2
a
+
γ
2
a
U
2
a
,
(25)
where
∆
a
=
ω
−
ω
a
and
ω
=
m
(Ω
orb
−
Ω
s
)
. In linear theory,
the parent’s energy and dissipation rate are smallest when the par-
ent is half-way between resonances, i.e., when the detuning is at a
maximum
|
∆
a
|
=
|
∂ω
a
/∂n
a
|
/
2
'
ω
a
/
2
n
a
(
γ
a
for the pe-
riods of interest). The linear energy of a parent half-way between
resonances is
E
a,
lin
E
0
= 7
.
9
×
10
−
18
(
2
y
1 +
y
)
2
(
P
orb
50 min
)
−
9
.
3
(26)
10
30
100
n
a
10
-8
10
-6
10
-4
|
Q
a
|
Q
a
2
.
1
n
−
3
.
7
a
0.01
0.03
ω
a
/ω
0
Figure 3.
Linear tidal overlap
Q
a
(blue circles). The black line is the fit
given by Equation (23).
10
30
100
n
a
10
2
10
3
10
4
10
5
|
abc
|
abc
abb
Approximate
abc
41(
n
a
/l
a
)
2
0.01
0.03
ω
a
/ω
0
Figure 4.
Three-mode-coupling coefficient
κ
abc
as a function of the par-
ent mode’s radial order. The blue circles are the coupling computed with
daughter pairs that have the smallest frequency detuning with respect to
the parent mode, and the orange circles further restrict it to self-coupled
daughters (
b
=
c
). The green circles connected with solid lines use the ap-
proximate expression for the coupling coefficient integrand given by Equa-
tion (A18). Note that for a given parent mode
a
, its
κ
abc
for different daugh-
ter pairs satisfying the selection rules are all approximately equal as long as
|
n
b
−
n
c
|
.
n
a
.
assuming a non-rotating WD such that
P
a
'
P
orb
/
2
.
Now consider a simple three-mode system consisting of a par-
ent mode driven by the tide coupled to a resonant daughter pair. If
E
a,
lin
> E
th
the parent is unstable even at maximum
∆
a
, where
the energy threshold (see, e.g., WAQB12 and Essick & Weinberg
2016, hereafter EW16)
E
th
E
0
=
1
4
κ
2
abc
(
γ
b
γ
c
ω
b
ω
c
)
[
1 +
(
∆
bc
γ
b
+
γ
c
)
2
]
,
(27)
with
∆
bc
=
ω
+
ω
b
+
ω
c
the nonlinear detuning. Note that if
ω >
0
,
we have
ω
b
,ω
c
<
0
according to our sign convention.
In Figure 5, we show
E
a,
lin
(dotted line) and the minimum
MNRAS
000
, 1–21 (2020)
Nonlinear tides in white dwarf binaries
5
10
20
30
50
100
200
P
orb
[min]
10
-21
10
-18
10
-15
10
-12
10
-9
E/E
0
E
a,
lin
E
th
E
brk
Numerical Search
Figure 5.
Linear energy
E
a,
lin
for a parent half-way between adjacent res-
onance peaks (Equation 26; dotted black line), the nonlinear threshold en-
ergy
E
th
based on the analytic scaling relations (Equation 29; blue line),
and the wave breaking energy
E
brk
(Equation 44; orange line) as a function
of orbital period assuming a non-rotating WD. The black crosses show the
minimum
E
th
obtained from a numerical search for daughter pairs.
E
th
from a numerical search of daughter pairs (black crosses) as-
suming a non-rotating WD. We also show an analytic estimate of
the minimum
E
th
(blue line), whose calculation we describe below.
We see that for
P
orb
.
150 min
, even maximally detuned parent
modes are parametrically unstable. In fact, since
E
a,
lin
E
th
over much of this range, we will see that a single parent excites
many unstable daughter pairs.
The daughter pairs that minimize
E
th
are those that satisfy the
κ
abc
selection rules, have
|
n
b
−
n
c
|
< n
a
, and nonlinear detunings
∆
bc
≈
γ
b
+
γ
c
since they minimize the sum in quadrature in the
brackets of Equation (27); see also EW16. We can obtain an ana-
lytic estimate of the minimum
E
th
by using the scaling relations
given above and an estimate for the minimum
∆
bc
. Following the
argument given by Wu & Goldreich (2001), we obtain an estimate
for the minimum
∆
bc
by noting that for a fixed parent mode
a
,
there are
∼
n
2
a
daughter pairs satisfying
|
n
b
−
n
c
|
.
n
a
at fixed
l
b
and
l
c
. As we allow the angular degree
l
b
and
l
c
to vary, we ob-
tain an extra factor of
l
b
l
a
of mode pairs that satisfy the condition
|
l
b
−
l
c
|
< l
a
< l
b
+
l
c
. The eigenfrequencies of these potential
daughter modes span a range of order
n
a
|
∂ω
b
/∂n
b
| '
n
a
ω
b
/n
b
.
Therefore, the typical minimum three-mode detuning assuming a
non-rotating WD is of order
3
∆
bc
≈
1
l
b
l
a
n
2
a
n
a
ω
b
n
b
'
0
.
07
l
a
ω
0
l
2
b
n
3
a
,
(28)
where in the second approximation we first eliminated
n
b
in terms
of
(
ω
b
,l
b
)
using Equation (19) and then assumed that
|
ω
b
|'
ω
a
/
2
.
The factor of 0.07 is from a fit to our numerical search for daugh-
ter pairs that minimize
E
th
. By using the scaling relations in Sec-
tion 3.2 and Equation (28) and setting
∆
bc
'
γ
b
+
γ
c
, it follows
that the minimum threshold energy
E
th
E
0
'
1
.
0
×
10
−
20
(
0
.
18
T
)
2
(
P
orb
50 min
)
−
3
,
(29)
3
For a rotating WD, the detuning is smaller by yet another factor of
l
b
because the degeneracy between different combinations of
m
b
and
m
c
is
lifted (i.e., rotational splitting).
where this assumes a non-rotating WD and
P
a
=
P
orb
/
2
. As Fig-
ure 5 shows, this
E
th
estimate is in good agreement with that from
the numerical search for daughter pairs. The daughters that mini-
mize
E
th
typically have
l
b
'
6
.
0
(
P
orb
50 min
)
−
5
/
4
and
n
b
'
281
(
P
orb
50 min
)
−
1
/
4
.
(30)
3.4 Energy and angular momentum transfer
In this Section we derive the tidal power
̇
E
tide
and the tidal torque
τ
tide
in the inertial frame. We define the torque to be from the orbit
to the WD and when
τ
tide
>
0
the tide spins up the WD. Given the
interaction Hamiltonian
H
int
=
−
2
E
0
∑
ω
a
>
0
(
q
∗
a
U
a
+
q
a
U
∗
a
)
,
(31)
the tidal torque acting on the WD is
τ
tide
=
∂H
int
∂
Φ
'−
4
∑
ω
a
>
0
Re
[
q
∗
a
∂U
a
∂
Φ
]
E
0
=
−
4
∑
ω
a
>
0
m
a
Im [
q
∗
a
U
a
]
E
0
,
(32)
where a factor of two arises from the sum over modes and their
complex conjugate and another because we restricted the sum
to positive frequencies, and the last equality follows because
∂U
a
/∂
Φ =
−
i
mU
a
. Note that we dropped the term
∝
U
ab
in the
interaction Hamiltonian because only the linearly resonant parents
have a significant direct coupling to the tide (Section 3.1).
The associated tidal power, assuming a circular orbit, is given
by
̇
E
tide
= Ω
orb
τ
tide
.
(33)
In general, this will power a combination of mode energy,
tidal heating, and WD spin energy. However, as we illustrate
in Section 5.1 (see, e.g., Figures 7 and 8), in steady state the
time-averaged total mode energy is approximately constant and
∑
a
̇
E
a
'
∑
a
( ̇
q
∗
a
q
a
+
q
∗
a
̇
q
a
)
E
0
'
0
. Using Equation (15) we
thus have, in a time-averaged sense,
∑
a
ω
a
Im [
q
∗
a
U
a
]
'−
∑
b
γ
b
q
∗
b
q
b
.
(34)
The summation on the left-hand side is only over parent modes
since only they feel a strong, direct driving by the tide (Section 3),
whereas on the right-hand side it is over all modes from all gener-
ations. We also dropped the three-mode dissipation terms as they
contribute little to the total dissipation.
4
For the most-resonant parent modes with
ω
a
≥
0
and
l
a
= 2
,
the azimuthal order
m
a
=
m
= 2
and
ω
a
'
ω
=
m
(Ω
orb
−
Ω
s
)
.
4
There are two dissipation terms that arise directly from three mode
coupling: the first originates from
∑
a
( ̇
q
∗
a
q
a
+
q
∗
a
̇
q
a
)
and contributes
∑
abc
2(
ω
a
+
ω
b
+
ω
c
)
κ
abc
Im
[
q
a
q
b
q
c
]
and the second comes from the
nonlinear piece in the total mode energy [Equation (18)] and contributes
∑
abc
2
γ
a
κ
abc
Re[
q
a
q
b
q
c
]
. Since the detuning
(
ω
a
+
ω
b
+
ω
c
)
∼
γ
a
for the most unstable daughters, the two terms are comparable. Since
|
∑
abc
κ
abc
q
a
q
b
q
c
|
∑
a
E
a
(see last paragraph of this Section), the
nonlinear dissipation is much smaller than the lower-order contribution
∑
a
γ
a
E
a
.
MNRAS
000
, 1–21 (2020)
6
H. Yu, N. N. Weinberg and J. Fuller
We can therefore relate the tidal torque and power to the total dis-
sipation rate inside the star,
̇
E
diss
, as
τ
tide
=
m
ω
̇
E
diss
,
(35)
̇
E
tide
=
Ω
orb
(Ω
orb
−
Ω
s
)
̇
E
diss
,
(36)
where
̇
E
diss
= 4
∑
ω
b
>
0
γ
b
q
∗
b
q
b
E
0
'
4
∑
ω
b
>
0
γ
b
E
b
.
(37)
If we assume that the WD rotates with uniform angular velocity
Ω
s
,
then the tidal torque spins it up at a rate
̇
Ω
s
=
τ
tide
I
WD
,
(38)
and the orbital frequency changes at a rate
̇
Ω
orb
=
̇
Ω
orb
,
gw
+
3
τ
tide
μD
2
,
(39)
where the GW induced orbital decay rate
̇
Ω
orb
,
gw
=
96
5
(
G
M
c
c
3
)
5
/
3
Ω
11
/
3
orb
.
(40)
Here
μ
is the reduced mass and
M
c
=
y
3
/
5
(1 +
y
)
−
1
/
5
M
is
the chirp mass. In our study, we obtain
̇
E
diss
from our mode net-
work simulations according to Equation (37) and thereby determine
τ
tide
,
̇
Ω
s
, and
̇
Ω
orb
.
5
Our large nonlinear networks display complicated dynam-
ics. Nonetheless, some insights can be gained by considering the
nonlinear equilibrium of simple three-mode systems (WAQB12,
EW16). For such a system, the parent mode’s saturation is
E
a,
s
=
E
th
and for
E
a,
lin
E
a,
s
, the daughter mode’s equilibrium is
E
b,
s
E
0
'
√
γ
c
ω
b
γ
b
ω
c
∣
∣
∣
U
a
2
κ
abc
∣
∣
∣
'
6
.
7
×
10
−
16
(
2
y
1 +
y
)(
P
orb
50 min
)
−
7
.
7
.
(41)
Comparing Equations (29) and (41), we see that
E
b,
s
'
E
c,
s
E
a,
s
. As a result, the leading order drive to the granddaugh-
ters will be via daugher-granddaughter three-mode coupling rather
than parent-granddaughter coupling at higher nonlinear orders. We
therefore include multiple generations in our networks but only ac-
count for three-mode coupling between adjacent generations.
6
Energy is stored not only in each individual mode (
∝
E
a
)
but also in the three-mode couplings [
∝
Re [
κ
abc
q
a
q
b
q
c
]
; see
Equation (18)]. However, the latter makes a negligible contribu-
tion to the total energy at saturation. We can easily see this for a
three-mode system, since at saturation the total nonlinear energy is
∼
κ
abc
√
E
a,
s
E
b,
s
E
b,
s
. To see roughly why this also holds
for our large mode networks, note that at saturation the nonlinear
forces approximately balance the linear forces. By Equation (11),
this implies
|
∑
bc
κ
abc
q
∗
b
q
∗
c
|
approximately equals
|
U
a
|
for parent
5
Alternatively, we can compute
τ
tide
by taking the time average of Equa-
tion (32) in steady state. Although we verified that the two methods yield
consistent results, in practice we find that
Im[
q
∗
a
U
a
]
of individual parents
is much more oscillatory than
∑
b
γ
b
E
b
and thus Equation (37) provides a
more numerically accurate estimate of the torque.
6
Four-mode coupling can be important for the
p
-
g
instability (Venumad-
hav et al. 2014; Weinberg 2016). However, that is a nonresonant instability
whereas here we focus on the resonant parametric instability.
modes and
(
γ/ω
a
)
|
q
a
|
for the other modes. Since both are
|
q
a
|
,
it follows that
|
∑
bc
κ
abc
q
a
q
b
q
c
|
E
a
and therefore the nonlin-
ear energy is a small contribution to the total energy.
3.5 Standing waves vs. traveling waves
The relations above and our mode network calculations assume
that the modes are all standing waves. In order to be a stand-
ing wave, a mode’s linear damping time must be longer than its
group travel time through the propagation cavity (which here spans
much of the WD radius; see Figure 1),
T
a
= 2
∫
d
rv
−
1
grp
=
2
∫
d
r
(
d
ω
a
/
d
k
r
)
−
1
, where
v
grp
is the mode’s group velocity. Oth-
erwise, it is a traveling wave. Defining the inverse group travel time
α
a
= 2
π/T
a
, we find
α
a
'
0
.
28
l
a
n
2
a
ω
0
= 0
.
097
l
a
n
2
a
rad s
−
1
.
(42)
In Figure 2, we compare
α
a
to the linear damping rate of modes.
We find that the standing wave condition
γ
a
< α
a
is satisfied for
n
a
.
1200
(
l
a
2
)
1
/
4
,
i.e.,
P
orb
.
1320
(
l
a
2
)
−
3
/
4
min
,
(43)
which is true of all the modes in our networks.
Another necessary condition for standing wave is that the
shear
|
d
ξ
r
/
d
r
| ' |
k
r
ξ
r
|
be everywhere less than unity, where
ξ
r
Y
lm
is the radial component of the physical Lagrangian displace-
ment
ξ
. If a g-mode’s shear exceeds unity, it is strongly nonlinear
and overturns the local stratification and breaks (see, e.g., Goodman
& Dickson 1998; Barker 2011).
Fuller & Lai (2012a; hereafter FL12) and Burkart et al. (2013;
hereafter BQAW13) use this local wave-breaking condition to ad-
dress the onset of nonlinear tidal effects in WD binaries. They show
that at sufficiently short orbital periods, the tide excites internal
gravity waves that are initially linear deep within the WD but be-
come nonlinear and break as they approach the stellar surface.
7
We first evaluate the wave-breaking condition
assuming a
standing wave
, i.e., a g-mode. Using the approach described in Ap-
pendix A1, we find that a g-mode’s shear exceeds unity if its energy
exceeds [see Equation (A4)]
E
brk
E
0
= 3
.
6
×
10
−
12
(
P
orb
50 min
)
−
2
.
(44)
In Figure 5, we show
E
brk
as a function of
P
orb
. We find that
E
a,
lin
first exceeds
E
brk
at
P
orb
'
10
min. Moreover, even highly res-
onant parent modes are unlikely to break before
P
orb
'
10
min.
That is because the parent is parametrically unstable (
E
a,
lin
>
E
th
) out to
P
orb
≈
150
min (Figure 5) and excites secondary
modes which prevent it from reaching the wave-breaking limit (see
Section 5.1).
We now evaluate the wave-breaking condition
assuming a
traveling wave
rather than a standing wave. Specifically, we use
the approach described in FL12 (see Appendix C for a brief synop-
sis) to find the traveling-wave solution of the linear inhomogeneous
tidal equations [Equations (C1) - (C3)]. Just as the standing wave
assumption is valid only if max
|
k
r
ξ
r
|
<
1
, the traveling wave as-
sumption is valid only if max
|
k
r
ξ
r
|
>
1
. In Figure 6 we show
7
It is interesting to note that whereas the the local wave-breaking occurs
at the surface, the global three-mode coupling happens mostly in the core
region. See Appendix A4 and Figure 4. This is different from the case of
solar models (WAQB12).
MNRAS
000
, 1–21 (2020)
Nonlinear tides in white dwarf binaries
7
3
5
10
20
30
50
100
P
orb
[min]
10
−
4
10
−
3
0.01
0.1
0.1
1
10
max
r
(
|
k
r
ξ
r
|
)
Numerical Results
0
.
008
×
(
P
orb
/
50 min)
−
3
Figure 6.
Maximum local shear
|
k
r
ξ
r
|
from the traveling wave solution for
a non-rotating WD with
T
eff
= 9000 K
.
max
|
k
r
ξ
r
|
computed under the traveling-wave assumption for a
non-rotating WD with
T
eff
= 9000 K
. We find that the upper en-
velope of the shear
∝
P
−
3
orb
(which we explain in Appendix C) and
reaches unity at
P
orb
'
10 min
, consistent with the results assum-
ing a standing wave.
The weakly nonlinear regime of this study therefore spans a
large range of orbital periods (
10
.
P
orb
/
min
.
150
). Evaluat-
ing the global, multi-mode dynamics in this regime is essential for
understanding the impact of tidal dissipation on WD binaries.
Our analysis assumes that all modes, not just the parent, are
standing waves and thus below the wave breaking threshold. Since
E
th
∝
ω
3
while
E
brk
∝
ω
2
, higher generation (i.e., lower fre-
quency) modes have even smaller ratios of
E
th
to
E
brk
than the
parent. Thus, the excited daughters, granddaughters, etc. will likely
become parametrically unstable and saturate before breaking (see
also Appendix F in EW16). In practice, because we truncate our
networks at the fifth generation and include only the most resonant
pairs for each generation (since including more modes does not sig-
nificantly increase the calculated
̇
E
diss
; see Section 5.2), some high
generation modes in our network can have shears that momentar-
ily exceed unity. However, at any given time these represent only a
very small fraction of the excited modes and thus they are unlikely
to modify the overall dynamics and dissipation.
It is also worth noting that the shear can be a sensitive function
of the WD temperature. For example, in Appendix A1 we show that
for
T
eff
= 18000 K
the maximum shear is about an order of mag-
nitude larger than for
9000 K
. On the other hand, this is compen-
sated by the tidal synchronization which decreases the driving fre-
quency (see Section 5.4 and Appendix C). The orbital period where
the dynamical tide transitions from weakly to strongly nonlinear is
therefore still
≈
10
min when the effects of both temperature and
synchronization are taken into account.
4 NUMERICAL IMPLEMENTATION
The modes in our networks oscillate near their eigenfrequencies
and have small linear detunings
∆
a
(parents) or nonlinear detun-
ings
∆
bc
(daughters, granddaughters, etc.). We can therefore fac-
tor out the fast-oscillations by transforming coordinates to
c
a
=
q
a
exp(
i
ω
a
t
)
, similar to the approach of previous mode network
studies (Brink et al. 2005, EW16). The parent mode amplitude
Equation (15) is then
̇
c
a
+
γ
a
c
a
=
i
ω
a
|
U
a
|
e
−
i
∆
a
t
+
i
ω
a
∑
bc
κ
abc
c
∗
b
c
∗
c
e
i
∆
′
bc
t
,
(45)
and similarly for the other modes, where
∆
′
bc
=
ω
a
+
ω
b
+
ω
c
.
We initialize our networks by starting each mode at its linear tidal
energy and a random phase. We implemented the calculations in
Python
and used the
NUMBA
package (Lam et al. 2015) to en-
hance the computational performance.
Initially, the amplitudes of the unstable daughters will grow
exponentially at a characteristic rate (see WAQB12)
Γ
nl
'
ω
a
κ
abc
√
E
a
E
0
.
(46)
This allows us to define a characteristic nonlinear growth timescale
T
nl
≡
1
Γ
nl
∼
0
.
12
(
P
orb
50 min
)
3
.
7
yr
,
(47)
where the numerical value is for a parent mode at an initial energy
E
a,
lin
.
We find that the mode networks saturate and reach a nonlinear
equilibrium over a few nonlinear growth times
T
nl
. This is much
shorter than the GW-induced orbital decay timescale
T
gw
=
Ω
orb
̇
Ω
orb
,
gw
= 4
.
8
×
10
7
(
P
orb
50 min
)
8
/
3
yr
,
(48)
where the numerical value assumes a typical WD binary with
M
=
M
′
= 0
.
6
M
. The timescale
T
nl
is also shorter than the time it
takes for the GW orbital decay to change the three-mode detuning
∆
bc
by an amount
(
γ
b
+
γ
c
)
'
2
γ
b
(see Section 3.3),
T
det
=
γ
b
̇
Ω
orb
,
gw
'
150
(
P
orb
50 min
)
19
/
6
yr
.
(49)
Therefore, the particular parametrically unstable pairs that are most
resonant and thus have the lowest
E
th
do not change on a timescale
of a few
T
nl
. We therefore only construct our mode networks once
for each
P
orb
we consider.
In order to construct our mode networks, we search for the
daughters, granddaughters, etc. with the lowest threshold energies.
Numerically, we find that a network’s total energy dissipation rate
̇
E
diss
converges once we include five mode generations constructed
as follows. The first generation (parents) includes the two most lin-
early resonant modes. The second generation (daughters) includes
the three lowest threshold daughter pairs of each parent. Since the
two parent modes both oscillate at the tidal driving frequency, they
usually have the same pair of most-resonant daughter modes and
thus the second generation typically has 6 modes instead of 12. The
third through fifth generations include the single lowest threshold
pair of each mode from the previous generation. A typical network
consists of 92 modes, with
(2
,
6
,
12
,
24
,
48)
modes in each gener-
ation (since modes sometimes appear in more than one pair, some
networks have slightly fewer than 92 modes). We find that increas-
ing the number of modes and generations does not significantly
change the computed
̇
E
diss
(see Section 5.2).
A collective instability can occur if daughters form large sets
of mutually coupled pairs (WAQB12). Collectively unstable daugh-
ters initially grow much more rapidly than the isolated pairs de-
scribed above. However, in our problem the collective instability
threshold
E
th
,
col
is higher than the isolated pair instability thresh-
old
E
th
. EW16 found that the parents, whose linear energy might
be well above
E
th
,
col
, reach a nonlinear equilibrium at an energy
MNRAS
000
, 1–21 (2020)