# Nonlinear tidal excitation of super-harmonic gravity waves in main-sequence
stars in binary and exoplanetary systems

P. B. Ivanov1, S. V. Chernov 1, A. J. Barker 2
1Astro Space Centre, P.N. Lebedev Physical Institute, 84/32 Profsoyuznaya
Street, Moscow, 117997, Russia
2 Department of Applied Mathematics, School of Mathematics, University of
Leeds, Leeds, LS2 9JT, UK E-mail<EMAIL_ADDRESS>(PBI)E-mail:
<EMAIL_ADDRESS>(SVCh)E-mail<EMAIL_ADDRESS>(AJB)

(Accepted. Received; in original form)

###### Abstract

We study the role of nonlinear effects on tidally-excited internal gravity
waves in stellar radiation zones in exoplanetary or binary systems. We are
partly motivated to study tides due to massive short-period hot Jupiters,
which preferentially orbit stars with convective cores, for which wave
breaking near the stellar centre cannot operate. We develop a theory (and test
it with numerical calculations) for the nonlinear excitation of super-harmonic
“secondary” waves (with frequencies $2\omega_{p}$) by a “primary” tidal wave
(with frequency $\omega_{p}$) near the interface between the radiation zone
and convective envelope. These waves have the same horizontal phase speeds to
leading order, and this nonlinear effect could contribute importantly to tidal
dissipation if the secondary waves can efficiently damp the primary. We derive
criteria involving the orbital and stellar parameters required to excite these
secondary waves to large amplitudes using a local model of the
radiative/convective interface, which we convert to apply to tides in a
spherical star. We numerically evaluate the critical amplitudes required for
this new nonlinear effect to become important using stellar models, comparing
them to the “conventional” criteria for wave breaking in radiative cores and
the application of WKBJ theory near convective cores. The criteria for this
new effect are easier to satisfy than the conventional measures of
nonlinearity in $1.4$ and $2M_{\odot}$ stars on the main-sequence. We predict
nonlinear effects to be important even for planetary-mass companions around
the latter, but this effect is probably less important in stars with radiative
cores.

###### keywords:

planet -star interactions, stars: binaries: close,

††pagerange: Nonlinear tidal excitation of super-harmonic gravity waves in
main-sequence stars in binary and exoplanetary systems–References††pubyear:
2022

##  1 Introduction

Tidal interactions play an important role in exoplanetary and close binary
systems. Although development of the modern theory of tides in celestial
bodies began over a century ago (starting with e.g. Darwin, 1880) there are
still many open problems in this area. Perhaps the most important questions
concern the efficiency of tidal energy dissipation in stars and giant planets,
which are wholly or partly fluid bodies. This is because tidal dissipation can
drive evolution of the spins and orbits of stars and planets in close binary
and in exoplanetary systems (see e.g. the reviews by Ogilvie, 2014; Mathis,
2019). It has long been argued that when the standard equilibrium (also
referred to as quasi-stationary) tides (e.g. Hut, 1981) are considered in
convection zones the dissipative efficiency could be negligibly small in many
applications involving planets and main-sequence stars (e.g. Goldreich and
Nicholson, 1977, who built upon earlier ideas by Zahn 1966). Modern numerical
calculations have confirmed this qualitative statement111Unless the mechanism
considered by Terquem (2021) can work as efficiently as they estimate, though
this is highly uncertain (Barker and Astoul, 2021)., although the detailed
physical picture of the processes involved may be different from what was
originally envisaged (see e.g. Ogilvie and Lesur, 2012; Duguid, Barker and
Jones, 2020a, b; Vidal and Barker, 2020a, b). Thus, the equilibrium (quasi-
stationary) tides are unlikely to explain the observed parameters of main-
sequence close binary and exoplanetary systems, including the orbital
evolution inferred in some hot Jupiter systems. However, equilibrium tides are
still likely to be the dominant mechanism in giant stars (e.g. Verbunt and
Phinney, 1995; Mustill and Villaver, 2012). Another potentially important
tidal dissipation mechanism is dynamical tides (e.g. Cowling, 1941; Zahn,
1977), which involve the resonant excitation of low frequency normal modes (or
waves) of a planet or a star by the time-dependent tidal potential. Usually,
the dynamical tide is thought to consist of internal gravity waves in stably-
stratified (radiative) regions, and inertial waves in convective regions of
rotating bodies.

In this paper, we consider non-rotating main-sequence stars, or sufficiently
slowly rotating ones for which all relevant tidal frequencies are larger than
the Coriolis frequency (like in most hot Jupiter hosts), therefore tidal
forcing cannot directly excite inertial waves. On the other hand, internal
gravity waves existing in stably-stratified radiative zones may be resonantly
excited by tidal forcing, and, in principle, their dissipation can explain
certain observational phenomena (including WASP-12 b’s inferred orbital decay
e.g. Maciejewski et al., 2016; Chernov et al., 2017; Weinberg et al., 2017;
Patra et al., 2020; Yee et al., 2020; Barker, 2020) provided that these waves
are sufficiently strongly dissipated inside the star. The regime in which this
is the case has been referred to as “moderately large dissipation” (MLD, by
e.g. Ivanov et al., 2013; Chernov et al., 2017) or as the “fully damped” or
”travelling wave” regime (by e.g. Barker and Ogilvie, 2010; Barker, 2020).
Qualitatively, this regime is expected when the propagation time of a tidally-
excited gravity wave packet is larger than its damping time222Since the group
speed of a high radial order gravity wave is much smaller than its phase
speed, its propagation time can be much larger than the dynamical time, so the
damping rate (and hence the viscosity or radiative diffusivity damping the
wave) need only be ‘moderately’ large.. However, the ‘standard’ linear
mechanisms to damp gravity waves (radiative diffusion and convective damping
of the evanescent tails in the envelope) are usually not efficient enough to
justify the validity of this regime for close orbits (e.g. Terquem et al.,
1998; Goodman and Dickson, 1998). In this situation various non-linear effects
have been proposed for which this fully damped (MLD) regime may be possible,
such as weakly non-linear mode-mode interactions (e.g. Barker and Ogilvie,
2011; Weinberg et al., 2012; Essick and Weinberg, 2016), or wave breaking near
the centres of solar-like stars possessing radiative cores in which the waves
can be geometrically focused and attain large amplitudes (Goodman and Dickson,
1998; Ogilvie and Lin, 2007; Barker and Ogilvie, 2010; Barker, 2011, 2020).
Note that, although different theoretical approaches (for example, the normal
mode formalism of Ivanov et al. 2013 or the low-frequency asymptotic approach
taken by e.g. Goodman and Dickson 1998) should, in principle, give the same
results for, say, the conditions under which non-linear effects are predicted
to be important, in practice they use various different simplifying
assumptions, and could therefore differ in their quantitative predictions.
Hence, it is worthwhile reproducing results obtained in one formalism with
those of another when possible.

Wave breaking near the centres of solar-type stars with radiative cores can
lead to efficient wave absorption through the formation a critical layer (e.g.
Barker and Ogilvie, 2010; Barker, 2011), which can naturally explain the
occurrence of the fully damped/MLD regime. For example, it has been proposed
that WASP-12 has a radiative core due to being a subgiant, for which this
mechanism may explain the inferred orbital evolution (e.g. Weinberg et al.,
2017). In the current Sun, the criterion for the onset of wave breaking
requires that a planet in a one-day orbit must exceed about 3 Jupiter masses,
though this mass threshold is a strong function of the stellar mass and age
(Barker and Ogilvie, 2010; Barker, 2020). At the end of the main-sequence,
much lower mass planets can cause wave breaking in the stellar core and
potentially be destroyed. However, stars with masses even slightly larger than
1.1 $M_{\odot}$ typically have convective cores on the main-sequence. Although
a similar criterion for wave breaking can be formulated in the radiation zones
of these stars using WKBJ theory (e.g. Sun et al., 2018; Barker, 2020), the
threshold companion masses to cause wave breaking are much larger than the
ones obtained in solar-like stars, and this effect appears unlikely to ever be
important for planetary-mass companions until the star evolves off the main
sequence. It is important to verify the validity of this result, and to
explore whether tides in F-type stars can be in this fully damped regime or
not, because many of the most massive ultra-short period hot Jupiters have
been observed to orbit such stars (such as WASP-18b, Wilkins et al., 2017).

The criterion for wave breaking is equivalent to the condition that the radial
gradient of the total specific entropy, including both the radiative
background and the wave, becomes negative. This condition is always formally
satisfied in the radiation zone at a point sufficiently close to the interface
between the inner radiative region and an outer convective envelope (as we
will show in Section 3). This suggests that non-linear effects acting on tidal
perturbations could be important in this region of a star even though this is
not predicted from WKBJ theory. In this paper, we analyse the nonlinear
dynamics of these perturbations in some detail by considering a region of
small radial extent near such an interface. We use weakly nonlinear theory
(second order perturbation theory) to find a condition for which the
amplitudes of the generated secondary waves by the primary tidal waves become
comparable with the amplitude of the latter. We apply the formalism of Ivanov
et al. (2013) assuming the fully damped/MLD regime to describe the primary
tidal wave (i.e. the first order perturbations). Additionally, for simplicity,
we assume that, after certain modifications, results obtained in planar
geometry can be applied to a spherical star. We also adopt the Boussinesq
approximation for the equations of motion (e.g. Spiegel and Veronis, 1960),
and assume that square of the Brunt-Väisälä frequency in the radiation zone
depends linearly on the distance to the interface provided that this distance
is small (see e.g. Barker, 2011; Ivanov et al., 2013, for a discussion of this
point).

We find, in agreement with results recently obtained in fluid dynamics and
oceanography (see e.g. Wunsch, 2017; Baker and Sutherland, 2020), that first
order tidal perturbations generate (through their nonlinear self-interaction)
super-harmonic second order perturbations with approximately double the
frequency and wavenumber (and hence the same horizontal phase speed). We
derive the conditions required for the secondary super-harmonic waves to
attain approximately the same amplitudes as their primary waves in our
Cartesian model, and then apply this to realistic stellar models assuming that
the amplitude of the primary tidal wave (at first order) is in the fully
damped/MLD regime. Our criterion is written in terms of the quantity
$q/(1+q)$, where $q$ is the mass ratio (secondary perturber mass/primary star
mass), which must be larger than a certain critical value, $C_{crit}$, which
is a function of the orbital period and stellar parameters. We apply this
criterion to a set of main-sequence stellar models with masses $M=1$, $1.4$
and $2M_{\odot}$ with different ages.

We find that the criterion for nonlinear self-interaction of the primary tidal
waves to be important in generating super-harmonic secondary waves near
radiative-convective interfaces can be satisfied much more easily than the
(WKBJ) criterion for wave breaking near convective cores, in most of our
models with $M>M_{\odot}$. It would be interesting to explore with future
numerical simulations whether tidally-excited primary waves that satisfy our
amplitude criterion could be damped efficiently enough by this mechanism to
validate the occurrence of the fully damped/MLD regime in these stars (as
assumed by e.g. Chernov et al., 2017; Barker, 2020). If so, then our results
may have important implications both for hot Jupiter and close binary systems.
In particular, the occurrence of this fully damped/MLD regime allows
straightforward prediction of orbital decay rates in hot Jupiter systems, and
orbital evolution in some close binary systems, as long as the properties of
the stars and orbital properties are known.

The structure of this paper is as follows. In § 2, we recap the main results
required to apply the formalism of Ivanov et al. (2013) to obtain the linear
tidal response. We then provide simple estimates to predict when nonlinearity
might be expected to be important for tidally-excited gravity waves in § 3 at
three particular locations: near the centres of radiative cores, at the
radiative interface with a convective core, and at the interface between a
radiation zone and a convective envelope. This motivates the more detailed
calculations of the generation of super-harmonic gravity waves by weakly
nonlinear interactions in a local Cartesian Boussinesq model of the transition
region between a radiation zone and a convective envelope in § 4 and 5. We
confirm our analytical results by comparing them with numerical calculations
in § 7. We then derive simple criteria to predict when this new nonlinear
effect is likely to become important in § 6 and then apply these to stellar
models in 8. Finally, we present our conclusions and a discussion in § 9. In
the main text we assume that the star is non-rotating. However, we briefly
discuss the most important correction to our results if the star rotates
slowly (which is the appropriate regime for many stars hosting short-period
planets) in Appendix A.

##  2 Linear tidal response in radiation zones

We adopt the following notation and conventions throughout this paper, unless
specified otherwise: $M_{*}$ and $R_{*}$ are the stellar mass and radius,
respectively, $G$ is the gravitational constant, $M_{p}$ is the perturber’s
mass, and the mass ratio is $q=M_{p}/M_{*}$. We define
$\Omega_{*}=\sqrt{GM_{*}/R_{*}^{3}}$ (i.e. the dynamical frequency) and the
mean density $\bar{\rho}=3M_{*}/(4\pi R_{*}^{3})$. We express all quantities
of interest in these natural units. Namely, the density $\rho$ is represented
as $\rho=\bar{\rho}\tilde{\rho}$, the Lagrangian displacement vector
${\boldsymbol{\xi}}=R_{*}\tilde{\boldsymbol{\xi}}$, and all quantities having
the dimension of frequency are expressed in units of $\Omega_{*}$. To
distinguish dimensional and dimensionless quantities we assign tildes to the
latter.

We follow closely Ivanov et al. (2013, hereafter IPCh) and assume that all
quantities of interest may be represented as a discrete Fourier series in time
and over the azimuthal angle $\phi$. In particular, the Lagrangian
displacement vector induced at a point in the star by tidal interactions,
${\boldsymbol{\xi}}$, can be represented as

${\mbox{\boldmath${\xi}$}}=\sum_{m,k}{\mbox{\boldmath${\xi}$}}_{m,k}e^{-i\omega_{m,k}t+im\phi}+c.c.,$
(1)

where $c.c.$ denotes the complex conjugate, it is implied that summation over
the azimuthal number $m$ contains only terms with $m=0$ and $2$, and the
forcing frequency

$\displaystyle\omega_{m,k}=k\Omega_{orb}-m\Omega_{r},$ (2)

where $\Omega_{orb}$ and $\Omega_{r}$ are the orbital frequency and spin
angular velocity of the star, respectively, with $k$ being an integer 333 Note
that equation (2) implies that $\omega_{m,k}$ are defined in the rotating
frame. In the inertial frame we have $\omega_{m,k}=k\Omega_{orb}$.. Hereafter
we assume that the star is effectively non-rotating by setting $\Omega_{r}=0$
in the main text, and briefly discuss the main effects caused by a slow
stellar rotation in Appendix A..

IPCh derived expressions for ${\mbox{\boldmath${\xi}$}}_{m,k}$ in terms of the
eigenvectors of free stellar pulsations, ${\mbox{\boldmath${\xi}$}}_{j}$, as

${\mbox{\boldmath${\xi}$}}_{m,k}={A_{m,k}\over
2}\sum_{j}{Q_{j}\over\omega_{m,k}n_{j}(\omega_{m,k}-\omega_{j}+i\omega_{\nu})}{\mbox{\boldmath${\xi}$}}_{j},$
(3)

where $j$ is an integer, which represents a sum over all of the free modes of
the star. Here $A_{m,k}$ are quantities characterising the amplitudes of the
Fourier components of the tidal potential, which are given explicitly in
Appendix A of IPCh, for example. For our purposes, we only need to know
$A_{2,2}$, which is given by equation (10) below. The quantities $Q_{j}$,
$n_{j}$ and $\omega_{j}$ are the so-called tidal overlap integral (determining
how efficiently a given mode is excited by the tidal potential), norm and
eigenfrequency of a particular free mode, while $\omega_{\nu}$ is its damping
rate.

For a non-rotating star perturbed by the quadrupolar component of the tidal
potential, ${\mbox{\boldmath${\xi}$}_{j}}$ can be expressed in terms of
spherical harmonics $Y^{m}_{2}$ as

${\mbox{\boldmath${\xi}$}}_{j}=e^{-im\phi}\left\\{\xi_{j}(r)Y^{m}_{2}(\theta,\phi){\bf
e}_{r}+\xi_{j,S}(r)r\nabla Y^{m}_{2}(\theta,\phi))\right\\},$ (4)

where the standard spherical polar coordinate system centred on the star
$(r,\theta,\phi)$ is used, ${\bf e}_{r}$ is the unit vector in the radial
direction, and the plane $\theta=\pi/2$ coincides with the orbital plane. The
presence of $e^{-im\phi}$ in front of the brackets in (4) stems from our
initial definition (1), in that the dependence of ${\xi}$ on $\phi$ is already
taken into account there and it is implied that the eigenvectors do not depend
on this angle.

For a non-rotating star we have

$\displaystyle Q_{j}$ $\displaystyle=$ $\displaystyle
2\int_{0}^{R_{*}}\rho(r)r^{3}(\xi_{j}+3\xi_{j,S})dr,$ (5) $\displaystyle
n_{j}$ $\displaystyle=$
$\displaystyle\int_{0}^{R_{*}}\rho(r)r^{2}({\xi_{j}^{2}}+6{\xi_{j,S}^{2}})dr,$
(6)

where $\rho(r)$ is the stellar density. Note that, for clarity, we explicitly
use the dimensional definition for $Q_{j}$ and $n_{j}$ here, contrary to e.g.
IPCh and references therein.

The summation in the expression for ${\mbox{\boldmath${\xi}$}}_{m,k}$ can be
approximately performed under the assumption that the spectrum of eigenmodes
is dense and regular. This is appropriate for the case of high-order
$g$-modes, for example. In this case only modes with frequencies approximately
equal to $\omega_{m,k}$ contribute to the sum over $k$ in (3) and the
difference is determined by the values of the factor
$\omega_{m,k}-\omega_{j}+i\omega_{\nu}$ in the denominator for different
modes, while other quantities can be taken out of the sum over $k$. Let’s
assume that $\omega_{m,k}$ is close to a particular eigenfrequency, with an
index $j=j_{0}$, and then write down
$\omega_{j_{0}}=\omega_{m,k}+\Delta\omega_{j_{0}}$, where
$|\Delta\omega_{j_{0}}|\ll|\omega_{j_{0}}|$ is a frequency offset. Also, we
assume that eigenfrequencies with indices $j\sim j_{0}$ depend approximately
linearly on the difference $l=j-j_{0}$:
$\omega_{j}=\omega_{j_{0}}+(d\omega/dj)l$, where $d\omega/dj$ stands for the
difference $\omega_{j_{0}}-\omega_{j_{0}-1}$. Under these assumptions, the
factor $\omega_{m,k}-\omega_{j}+i\omega_{\nu}\approx
i\omega_{\nu}-(d\omega/dj)l-\Delta\omega_{j_{0}}$, and we can write down

${\mbox{\boldmath${\xi}$}}_{m,k}=-{A_{m,k}Q_{j_{0}}\over
2\omega_{m,k}(d\omega/dj)_{j_{0}}n_{j_{0}}}S{\mbox{\boldmath${\xi}$}}_{j_{0}},\quad
S=\sum^{\infty}_{l=-\infty}{1\over\delta+l-i\kappa},$ (7)

where $\delta={\Delta\omega_{j_{0}}\over d\omega/dj}$ is a dimensionless
frequency offset, $\kappa={\omega_{\nu}\over d\omega/dj}$ is a dimensionless
damping rate, we formally extend summation to infinite limits and set
$\omega_{j_{0}}=\omega_{m,k}$ in the terms outside the sum. The sum $S$ can be
evaluated using complex variable theory (residue calculus) with the result

$S=-\pi\cot(\pi(i\kappa-\delta)).$ (8)

Of particular importance is the case of so-called ‘moderately large
dissipation’ (MLD), which occurs when the damping rate is larger than the
distance between two neighbouring eigenfrequencies i.e. $\kappa>1$. For our
calculations below related to this regime we use $\kappa=1$ for all values of
$\delta$. Note that in this case $S(\kappa=1,\delta)$ is very close to its
limiting value $S(\kappa\rightarrow\infty,\delta)=i\pi$.

##  3 Simple estimates for tidally-excited gravity waves to become nonlinear

It is often assumed that internal gravity waves break (or become essentially
nonlinear) when the maximum magnitude of the (negative part of the) radial
gradient of the Eulerian perturbation of the specific entropy,
$\partial_{r}s^{\prime}$ becomes larger than the background value, $ds/dr$
(see e.g. Ogilvie and Lin, 2007; Barker and Ogilvie, 2010, 2011), since then
there is a portion of the wave in which the entropy profile is overturned
(i.e. with a decreasing radial gradient which would be “convectively
unstable”). The Lagrangian perturbation of specific entropy, $\Delta s$, is
clearly zero for adiabatic perturbations. We note that $\Delta
s=s^{\prime}+{ds\over dr}\xi^{r}$, where $\xi^{r}$ is the radial component of
${\xi}$. Hence, from the condition $|\partial_{r}s^{\prime}|>|{ds\over dr}|$,
we readily find

$\left|{{\partial\over\partial r}}\xi^{r}+{{d^{2}s\over dr^{2}}\over{ds\over
dr}}\xi^{r}\right|>1,$ (9)

where this expression is meant to be evaluated in stellar radiation zones, and
we consider its maximal value over the coordinates $r,\theta,\phi$, and a
tidal forcing period.

In what follows, we consider the simplest case of an approximately circular
orbit, in which only the term with $m=k=2$ is important in the summation in
(1), and the relevant tidal potential component is $A_{2,2}$, as defined in
equation (A2) of IPCh, which has the following explicit form

$A_{2,2}\approx-{1\over 2}\sqrt{{6\pi\over 5}}{GM_{p}\over a^{3}},$ (10)

where $a$ is the orbital semi-major axis. We also have
$\omega_{2,2}=2\Omega_{orb}$, and from (1) and (7) we obtain

${\mbox{\boldmath${\xi}$}}=-{\pi\over 2}A_{2,2}{Q\over{d\omega\over
dj}\Omega_{orb}n}{\mbox{\boldmath${\xi}$}}_{j_{0}}\sin(\Psi),$ (11)

where $\Psi=2(\Omega_{orb}t-\phi)$. It is implied below that all quantities
associated with an eigenmode correspond to one having an approximately
resonant frequency $2\Omega_{orb}$. We do not show the corresponding mode
index $j_{0}$ in expressions for the overlap integral $Q$, the norm $n$, the
eigenfrequencies $\omega$, and their differences $d\omega/dj$ hereafter. Since
we require that $\xi^{r}$ in (9) should take its largest value over a wave
period, we set $\sin\Psi=1$ in (11).

From (4) it follows that the radial component of
${\mbox{\boldmath${\xi}$}}_{j_{0}}$, $\xi^{r}_{j_{0}}$, entering (9) is
expressed in terms of $\xi(r)$ through the factor
$e^{-2i\phi}Y^{2}_{2}(\theta,\phi)$, which is a function of the angle $\theta$
only. This factor takes its largest value at $\theta=\pi/2$, where it is equal
to ${1\over 4}\sqrt{{15\over 2\pi}}$. Note that $\xi(r)$ is assumed to be
known from the solution of the standard problem of free stellar pulsations. We
express $\xi^{r}_{j_{0}}$ in terms of $\xi(r)$, then substitute (10) in (11),
to obtain the radial component

$\xi^{r}={3\pi\over 16}{q\over 1+q}{\Omega_{orb}\over d\omega/dj}{Q\over
n}\xi.$ (12)

We then substitute (12) into (9) to obtain our criterion for nonlinearity

${q\over{1+q}}>C_{crit},$ (13)

where

$C_{crit}={16\over 3\pi}{\left|{d\omega\over
dj}\right|\over\Omega_{orb}}{n\over\left|Q\left({d\xi\over dr}+{{d^{2}s\over
dr^{2}}\over{ds\over dr}}\xi\right)\right|}.$ (14)

Note that (14) should be evaluated in radiative zones, and then its minimal
value should be used in (13) to determine whether or not wave breaking (or
other important nonlinear effects) is predicted.

###  3.1 Wave breaking in the radiative cores of solar-like stars

Firstly, we wish to apply the criterion (13) to the case of a solar-type star
with a radiative core, under the assumption that the minimum value of
$C_{crit}$ as a function of $r$ is reached at the centre. This is the case
that has been studied by e.g. Goodman and Dickson (1998); Ogilvie and Lin
(2007); Barker and Ogilvie (2010); Barker (2011); Weinberg et al. (2012);
Barker (2020). Following Ogilvie and Lin (2007) we assume that the Brunt-
Väisälä frequency, $N(r)$, depends linearly on the distance to the centre of
the star, i.e. $N=A_{centre}r$. We also assume that the entropy gradient is
proportional to $N$ (and $g\sim r$ there also). It then follows from Ogilvie
and Lin (2007) or Barker (2011), that near the centre $\xi$ can be expressed
in terms of Bessel functions as

$\xi(r)=C_{centre}r^{-3/2}J_{5/2}(R_{centre}r),$ (15)

where $C_{centre}$ and $R_{centre}$ can be found by matching (15) to a WKBJ
solution $\xi_{WKBJ}$ valid in the radiative region of a solar-like star for a
g-mode with a large radial order. Such a solution is given e.g. by equation
(88) of IPCh. Close to the centre it has the form

$\xi_{WKBJ}=-C_{WKBJ}{1\over\sqrt{\rho_{centre}A_{centre}}r^{2}}\sin(\sqrt{6}{A_{centre}r\over\omega}),$
(16)

where $\rho_{centre}$ is the central stellar density and $\omega$ is the mode
eigenfrequency. The expression (15) should match (16) at large values of the
argument of the Bessel function, which allows us to obtain

$R_{centre}=\sqrt{6}{A_{centre}\over\omega},\;\;\&\;\;C_{centre}=6^{1/4}\sqrt{{\pi\over
2}}{C_{WKBJ}\over\sqrt{\rho_{centre}\omega}}.$ (17)

We now consider the opposite limit of small values of the argument of (15) to
obtain our nonlinearity criterion. Taking into account (17) we obtain

$\xi\approx{6^{3/2}\over
15}{C_{WKBJ}A_{centre}^{5/2}\over\rho_{centre}^{1/2}\omega^{3}}r.$ (18)

Now, substituting (18) into (14) and noting that ${d^{2}s\over dr^{2}}$
vanishes under our assumptions, we obtain

$C_{crit,centre}={80\over 6^{3/2}\pi}{n\over
Q}\left|\frac{d\omega}{dj}\right|{\rho_{centre}^{1/2}\omega^{3}\over
C_{WKBJ}A_{centre}^{5/2}}\Omega_{orb}^{-1}$ (19)

As discussed in IPCh, for the modes described by WKBJ theory, the constant
$C_{WKBJ}$ can be expressed through the norm $n$ and eigenfrequency $\omega$
as (see their equation 109)

$C_{WKBJ}=\sqrt{{2n\over I}}\omega,\;\;\text{where}\quad
I=\int_{0}^{r_{c}}{dr\over r}N,$ (20)

and where $r_{c}$ is the radius at the base of the convective envelope, while
the frequency difference $d\omega/dj$ is given by their equation (125) with
$k=2$ and $\Lambda=6$, i.e.

$\frac{d\omega}{dj}={4\pi\over\sqrt{6}}{\Omega^{2}_{orb}\over I}.$ (21)

We substitute (20) in (18), then the result of this and (21) in (19), and set
$\omega=2\Omega_{orb}$, to obtain

$C_{crit,centre}={320\over 9\sqrt{2}}{\rho_{centre}^{1/2}\over
I^{1/2}A_{centre}^{5/2}}{1\over\hat{Q}}\Omega_{orb}^{3},$ (22)

where

$\hat{Q}=Q/\sqrt{n},$ (23)

is a ‘normalised’ overlap integral whose value does not depend on the
eigenfunction amplitudes.

In the low frequency asymptotic limit $\Omega_{orb}\rightarrow 0$, one can
show (see IPCh) that under the usual assumption that $N^{2}$ close to the base
of the convective envelope scales linearly with the distance from it,
$\hat{Q}\propto\Omega_{orb}^{17/6}$. Therefore,
$C_{N}\propto\Omega_{orb}^{1/6}$. Comparing with our criterion (13), we can
see that this scaling agrees with equations (A7) and (A9) of Ogilvie and Lin
(2007) and with equations (49) and (50) of Barker (2020). We have also
confirmed that our criterion gives quantitatively similar results to those in
the literature in the low frequency limit.

###  3.2 Wave breaking near convective cores

Stars with masses exceeding approximately $1.1M_{\odot}$ (in particular we
consider F-type or A-type stars) typically possess both convective envelopes
and convective cores with a radiative zone in between. The transition of
$N^{2}$ from its value in the intermediate radiative zone to the convective
core is often very abrupt relative to the wavelength of the tidal waves. In
this case it is reasonable to assume that the WKBJ approximation remains valid
all the way down to the transition radius, $r_{core}$. In such a situation,
the criterion for wave breaking is that the maximum value of
$|d\xi^{r}/dr|>1$, where $\xi^{r}$ is determined by (12). The maximal value
over one oscillation period is considered and $\xi$ entering (12) is given by
eq. (88) of IPCh at $r\rightarrow r_{core}$. We find

$\left|\frac{d\xi^{r}}{dr}\right|_{max}={9\sqrt{2}\over 32}{q\over
1+q}\rho^{-1/2}r^{-5/2}{Q\over n^{1/2}}{\sqrt{N_{core}I}\over\Omega_{orb}},$
(24)

where444Note that, typically, the dependence of $N^{2}$ on $r$ has very sharp
features in the vicinity of $r_{core}$ in stellar models. We believe that
these features should be discarded when evaluating $N_{core}$, since they are
probably unphysical and would be smoothed out by various hydrodynamical mixing
processes, including convective overshoot. $N_{core}=N(r\rightarrow
r_{core})$. Accordingly, from the condition $|d\xi^{r}/dr|_{max}>1$ we obtain

${q\over 1+q}>C_{crit,core},\quad C_{crit,core}={32\over
9\sqrt{2}}\rho^{1/2}r^{5/2}{\Omega_{orb}\over\sqrt{N_{core}I}}{1\over\hat{Q}},$
(25)

where we note that all quantities in (25) should be evaluated at $r\rightarrow
r_{core}$. The approach here is essentially the same as the WKBJ estimate in
equations (51) and (53) of Barker (2020) except that the amplitude of the wave
is determined by directly computing overlap integrals rather than applying the
energy flux obtained by an asymptotic low-frequency analysis (though except
for the shortest orbital periods these approaches should give similar
results).

###  3.3 Wave breaking near the transition from a radiation zone to a
convective envelope?

Formally, we can see that $C_{crit}$ can be arbitrary small just below the
base of a convective zone. In this region, the WKBJ approximation is strictly
invalid, and the corresponding wave solution predicts $\xi$ to be
approximately constant when ${r_{c}-r\over r_{c}}\ll 1$, where $r_{c}$ is
radius at the base of the convective zone. On the other hand, the term
${d^{2}s\over dr^{2}}/{ds\over dr}$ is expected to diverge when ${r_{c}-r\over
r_{c}}\rightarrow 0$, which leads to $C_{crit}$ formally tending to zero in
the same limit. This observation suggests that the situation close to the base
of a convective envelope deserves a special treatment. This is our aim in the
next few sections.

##  4 Local weakly nonlinear analysis of super-harmonic wave generation near
the base of a convective envelope

In this section, we consider the evolution of stellar perturbations near the
base of a convective envelope and assume that the square of the Brunt-Väisälä
frequency has a linear dependence on the distance from the interface at
$r_{c}$, such that

$N^{2}\approx A_{c}(r_{c}-r).$ (26)

We adopt units of length here in terms of $r_{c}$, and units of time in terms
of $\omega_{c}^{-1}$, where $\omega_{c}=\sqrt{A_{c}r_{c}}$. We assume that the
radial wavelength of gravity waves is small enough that we can approximately
describe them using Cartesian geometry with local coordinates $(x,y,z)$, such
that $z=r_{c}-r$, and $x,y$ are local horizontal coordinates. Additionally, we
employ the well-known Boussinesq approximation for the equations of motion
(Spiegel and Veronis, 1960), which is valid for low frequency gravity waves in
a region of small spatial extent. In this approximation, the corresponding
adiabatic non-linear equations of motion in the radiation zone can be written

$\dot{\bf U}+({\bf U}\cdot\nabla){\bf U}=-\nabla P+b{\bf
e}_{z},\quad\nabla\cdot{\bf U}=0,\quad\dot{b}+({\bf U}\cdot\nabla)b=-zU^{z}.$
(27)

Here a dot stands for an Eulerian time derivative, ${\bf U}$, $P$ and $b$ are
perturbations of the velocity and pressure, and the buoyancy
variable555Defined by $b=-g\rho^{\prime}/\rho_{c}$, where $\rho^{\prime}$ is
an Eulerian density perturbation, $\rho_{c}$ is a constant reference density,
and $g$ is the local acceleration due to gravity., respectively. We have
assumed that $N^{2}$ is given by equation (26) and that all dynamical
variables are dimensionless by being expressed in the units indicated above.
Note that from (27) it follows that

$U^{i}_{,k}U^{k}_{,i}=-\Delta P+{\partial b\over\partial z},$ (28)

where summation over repeated indices is assumed from now on, and $\Delta$ is
the Laplacian operator.

We now consider the $z$ component of the equations of motion, differentiate it
with respect to time, and apply the Laplacian operator to the resulting
expression. Then, using (28) to eliminate $\dot{b}$, we obtain

$\Delta\ddot{U}^{z}+\Delta_{\perp}(zU^{z})=-\Delta_{\perp}({\bf
U}\cdot\nabla)b+\dot{T},$ (29)

where $\Delta_{\perp}={\partial^{2}\over\partial
x^{2}}+{\partial^{2}\over\partial y^{2}}$, and

$T={\partial\over\partial z}(U^{i}_{,k}U^{k}_{,i})-\Delta({\bf
U}\cdot\nabla)U^{z}.$ (30)

In a linear analysis, the right hand side of (29) is ignored, which is
formally valid if the solution is of infinitesimally small amplitude. If we
substitute the linear solution to compute the terms on the right hand side of
(29) we can formulate a second order problem, which is referred to as a weakly
nonlinear analysis. Technically, our solution is then reduced to finding
solutions of the same linear equation, but with a forcing term determined by
the solution to the first order problem. We can then compare the first and
second order solutions for $U^{z}$. The condition that they are of the same
order can be considered as a condition for the breakdown of our perturbation
theory.

When looking for the linear solution, without loss of generality we can assume
that $U^{z}$ depends only on $z$, $x$ and $t$, and consider
$U^{z}=v(z)e^{i(\omega t+kx)}+C.C.$, where $C.C.$ stands for the complex
conjugate. Substituting this ansatz into (29) and neglecting the non-linear
terms, we see that it reduces to

$v_{,zz}=(k^{2}-{k^{2}\over\omega^{2}}z)v.$ (31)

Solutions of (31) can be represented in terms of Airy functions $Ai(x)$ and
$Bi(x)$:

$v(z)=C_{1}Ai(-y)+C_{2}Bi(-y),\quad
y={\left({k\over\omega}\right)}^{2/3}(z-\omega^{2}).$ (32)

In order to express the nonlinear source term in terms of the linear solution,
it follows from the continuity equation that $U^{x}_{,x}=-U^{z}_{,z}$ and,
accordingly, in the linear approximation we have $U^{x}={i\over
k}(v_{,z}e^{i\phi}-C.C.)$, where $\phi=\omega t+kx$. Using these relations we
obtain

$U^{i}_{,j}U^{j}_{,i}=2((v^{2}_{,z}-vv_{,zz})e^{2i\phi}+v_{,z}v^{*}_{,z}+vv^{*}_{,zz}+C.C).$
(33)

Since only the term proportional to $e^{2i\phi}$, $(v^{2}_{,z}-vv_{,zz})$, as
well as its complex conjugate, depends on time, and accordingly, contributes
to the source of second order perturbations, we only consider this term below.
Using equation (31), we obtain ${d\over
dz}(v^{2}_{,z}-vv_{,zz})={k^{2}\over\omega^{2}}v^{2}$. A similar calculations
shows that the term $({\bf U}\cdot\nabla)U^{z}$ entering (30) does not depend
on time at second order, and therefore we obtain

$\dot{T}=4i{k^{2}\over\omega}v^{2}e^{2i\phi},$ (34)

where we consider only the term proportional to $e^{2i\phi}$.

Now let us calculate the term $\Delta_{\perp}({\bf U}\cdot\nabla)b$ entering
(29). From the last equation in (27) it follows that to first order we have

$b=i{z\over\omega}(ve^{i\phi}-v^{*}e^{-i\phi}).$ (35)

Now we take into account that ${\bf U}\cdot\nabla={i\over
k}(v_{,z}e^{i\phi}-v^{*}_{,z}e^{-i\phi}){\partial\over\partial
x}+(ve^{i\phi}+v^{*}e^{-i\phi}){\partial\over\partial z}$ to obtain

$\displaystyle({\bf U}\cdot\nabla)b$ $\displaystyle=$
$\displaystyle{i\over\omega}(v^{2}e^{2i\phi}-{v^{*}}^{2}e^{-2i\phi}),$
$\displaystyle\Delta_{\perp}({\bf U}\cdot\nabla)b$ $\displaystyle=$
$\displaystyle-4i{k^{2}\over\omega}v^{2}e^{2i\phi}+C.C..$ (36)

Substituting (34) and (36) in (29) and representing the second order solution
in the form

$U^{z}_{(2)}=we^{2i\phi}+C.C.,$ (37)

we arrive at the equation describing second order perturbations:

${d^{2}\over
dz^{2}}w-4k^{2}w+{k^{2}\over\omega^{2}}zw=-{2ik^{2}\over\omega^{3}}v^{2}.$
(38)

This equation describes the nonlinear generation of super-harmonics by the
self-interaction of the primary wave, which we assume to have been (linearly)
tidally forced. In the next section we will determine approximate solutions to
this equation.

##  5 A weakly non-linear model eigenproblem

We now consider a model problem focussing on the interface between convective
and radiative regions, which may be shown to be relevant for the global
perturbations in a star, as we show later in section 6. We assume that
$U^{z}=0$ at finite distances $z=z_{r}$ into the stable/radiative zone
($z>0$), and $z=-z_{c}$ into the neutral/convective zone ($z<0$),
respectively. We calculate the spectrum of eigenmodes corresponding to
solutions of eq. (31) for this problem and consider a particular ‘primary’
mode with frequency $\omega=\omega_{0}$. The corresponding eigenfunction is
assumed to determine the source term on the right hand side of (38). The
response of super-harmonic second order waves is again solved through
decomposition of the solution over eigenmodes.

###  5.1 Linear solution in the neutral zone

In the neutral zone for $z<0$, the solution can be expressed in terms of
growing and decaying exponentials,

$v=C^{1}_{n}e^{-kz}+C_{n}^{2}e^{kz}=C^{2}_{n}(e^{kz}-e^{-k(z+2z_{c})}),$ (39)

where the last equality follows from our condition $U^{z}(z=-z_{c})=0$. From
equation (39) it follows that

$v_{,z}(z=0)=\kappa
v(z=0),\quad\kappa=k{(1+e^{-2kz_{c}})\over(1-e^{-2kz_{c}})}.$ (40)

This condition should be matched to the solution in the radiative zone. Note
that when $kz_{c}\ll 1$, $\kappa\approx 1/z_{c}$ does not depend on $k$.

###  5.2 Linear solution in the radiative zone and the eigenspectrum

When $z>0$ we represent the solution (32) in the form

$v=C\bar{v},\quad\bar{v}=(\cos(\phi)Ai(-y)+\sin(\phi)Bi(-y)),$ (41)

where $y=(z-\omega^{2})/\lambda_{*}$, and $\lambda_{*}={({\omega\over
k})}^{2/3}$. In the limit $z\rightarrow 0$ it follows from (41) that

$\displaystyle v$ $\displaystyle\approx$ $\displaystyle
C\left({(\cos\phi+\sqrt{3}\sin\phi)\over
3^{2/3}\Gamma(2/3)}-{(\cos\phi-\sqrt{3}\sin\phi)\over
3^{1/3}\Gamma(1/3)}{\omega^{2}\over\lambda_{*}}\right.$ (42)
$\displaystyle\left.\hskip 56.9055pt+{(\cos\phi-\sqrt{3}\sin\phi)\over
3^{1/3}\Gamma(1/3)}{z\over\lambda_{*}}\right).$

This implies that
$v_{,z}\propto{(\cos\phi-\sqrt{3}\sin\phi)\over\lambda_{*}}$. This ratio
should be finite in the considered limit $\lambda_{*}\rightarrow 0$, and,
therefore, $\phi$ should be close to $\pi/6$. We assume that
$\phi=\pi/6+\Delta$, where $\Delta$ is small, and substitute (42) in (40) to
obtain

$\Delta=-{\Gamma(1/3)\over\Gamma(2/3)}{3^{1/6}\over
2}\kappa\lambda_{*}\left(1-{\Gamma(1/3)\over\Gamma(2/3)}{\kappa\lambda_{*}\over
2\cdot 3^{1/3}}-\kappa\omega^{2}\right).$ (43)

Note that the last two terms in the brackets are unimportant in the limit
$\omega\rightarrow 0$ and are neglected hereafter.

We assume that $z_{r}\gg\lambda_{*}$. In the limit $z\gg\lambda_{*}$,

$v\propto\sin\left({2\over
3}\left({(z-\omega^{2})\over\lambda_{*}}\right)^{3/2}+{5\pi\over
12}+\Delta\right),$ (44)

and, from the condition $v(z=z_{r})=0$, we obtain

${2\over
3}{k\over\omega_{n}}{\left(z_{r}-\omega_{n}^{2}\right)}^{3/2}+\Delta(\omega_{n})=\pi\left(n-\frac{5}{12}\right),$
(45)

which is an equation for our eigenfrequencies $\omega_{n}$.

When $\omega_{n}$ is small it can be represented as
$\omega_{n}=\omega^{0}_{n}+\omega^{1}_{n}$, where

$\displaystyle\omega^{0}_{n}$ $\displaystyle=$ $\displaystyle{2\over
3}{kz_{r}^{3/2}\over\pi(n-5/12)},$ $\displaystyle\omega^{1}_{n}$
$\displaystyle=$ $\displaystyle-{3\omega^{0}_{n}\over
2z_{r}}\left({\omega^{0}_{n}}^{2}+{\Gamma(1/3)\over\Gamma(2/3)}{3^{1/6}\kappa\over
2\sqrt{z_{r}}}{\left({\omega^{0}_{n}\over k}\right)}^{5/3}\right).$ (46)

Note that from (46) it follows that when $\omega^{0}_{n}$ and $k$ are given by
the first expression, then $2\omega^{0}_{n}$ is determined by the same
expression with a doubled wavenumber $2k$. That means that there are
eigenfrequencies corresponding to $k$ and $2k$, which differ only by small
corrections, and accordingly there is a secondary mode which is in near-
resonance with the nonlinear source term due to the primary.

###  5.3 Weakly nonlinear generation of secondary waves

In a weakly nonlinear regime we assume that there is no back reaction of the
second order perturbations on the primary wave. In this case, nonlinear terms
entering eq. (29) are assumed to be given by the right hand side of (38). We
represent (29) in the form

$\Delta\ddot{U}^{z}+\Delta_{\perp}(zU^{z})=S,\quad
S={8ik_{p}^{2}\over\omega_{p}}C^{2}_{p}{\bar{v}_{p}}^{2}e^{2i(\omega_{p}t+k_{p}x)},$
(47)

where we assign the index $p$ to all quantities corresponding to the primary
mode, $\omega_{p}$ and $k_{p}$ are related by eq. (46), $v$ is given by (41),
we note that $\Delta_{\perp}={\partial^{2}\over\partial x^{2}}$, and we omit
the index $n$ from here on in this section. We seek series solutions to (47)
of the form

$U^{z}=\sum_{i}a_{i}(t)\phi_{i}(x,z),$ (48)

where $\phi_{i}$ are solutions of the eigenproblem

$\lambda_{i}\phi_{i}=\hat{A}\phi_{i},\quad\hat{A}=\Delta^{-1}\Delta_{\perp}(z\phi_{i})=-k_{s}^{2}\Delta^{-1}(z\phi_{i}),$
(49)

where $\Delta^{-1}$ is the inverse of the Laplacian operator, $k_{s}=2k_{p}$
and $\lambda_{i}=\omega_{i}^{2}$. It is easy to show that the set of
eigenfunctions $\phi_{i}$ are orthogonal with respect to the inner product

$N_{i,j}=\int z\phi_{i}\phi_{j}^{*}dxdz,$ (50)

where integration over $z$ is performed from $0$ to $z_{r}$ and we assume, for
simplicity, periodic boundary conditions in the $x$ direction with a period
$x_{b}\gg k^{-1}$. In order to prove this, we multiply (50) by
$\lambda^{*}_{j}$ to obtain

$\displaystyle\lambda^{*}_{j}N_{i,j}$ $\displaystyle=$ $\displaystyle\int
z\phi_{i}\Delta^{-1}\Delta_{\perp}(z\phi_{j}^{*})dxdz$ (51) $\displaystyle=$
$\displaystyle\int\Delta_{\perp}(z\phi_{i})\Delta^{-1}(z\phi_{j}^{*})dxdz$
$\displaystyle=$
$\displaystyle\int\Delta\Delta^{-1}\Delta_{\perp}(z\phi_{i})\Delta^{-1}(z\phi_{j}^{*})dxdz$
$\displaystyle=$ $\displaystyle\lambda_{i}\int z\phi_{i}\phi_{j}^{*}dxdz$
$\displaystyle+\lambda_{i}\int\nabla(\phi_{i}\nabla(\Delta^{-1}(z\phi_{j}^{*}))-\Delta^{-1}(z\phi_{j}^{*})\nabla\phi_{i})dxdz$
$\displaystyle=$
$\displaystyle\lambda_{i}N_{i,j}-{\lambda_{i}\lambda_{j}^{*}\over
k_{s}^{2}}\int{\partial\over\partial
z}(\phi_{i}\nabla\phi_{j}^{*}-\phi_{j}^{*}\nabla\phi_{i})dx,$

where we use (49), integration by parts, the facts that our boundary
conditions are periodic in the $x$ direction and that
$\Delta_{\perp}\phi_{i}=-k_{s}^{2}\phi_{i}$, and known properties of the
Laplacian. The last term in (51) is zero due to our condition (40), and we
arrive at $(\lambda^{*}_{j}-\lambda_{i})N_{i,j}=0$, which proves the
statement. In a similar way we can prove that the eigenvalues $\lambda_{i}$
are real and positive. Since $\lambda_{i}=\omega_{i}^{2}$, where $\omega_{i}$
are given by eq. (46) with $n=i$, this has already been shown and we omit the
proof.

We substitute (48) in (47), use (49), multiply the result by $\phi_{j}^{*}$,
and integrate over $x$ and $z$ to obtain

$\ddot{a}_{j}+\lambda_{j}a_{j}=-{\lambda_{j}\over k^{2}}{S_{j}\over N_{j}},$
(52)

where

$S_{j}=\int\phi_{j}^{*}Sdxdz,\quad N_{j}=N_{j,j}=\int
z\phi_{j}^{*}\phi_{j}dxdz.$ (53)

Remembering that $S\propto e^{2i\omega_{p}t}$, solutions to (52) are

$a_{j}={\omega_{j}^{2}\over
k_{s}^{2}(4\omega_{p}^{2}-\omega_{j}^{2})}{S_{j}\over
N_{j}}e^{2i\omega_{p}t},$ (54)

where we remember that $\lambda_{j}=\omega^{2}_{j}$.

###  5.4 Excitation of the near resonant secondary mode

As we have mentioned above, there is a secondary mode with $k_{s}=2k_{p}$,
whose frequency is close to double the frequency of the primary mode,
$2\omega_{p}$, and, accordingly, to the frequency of the source $S$. Both the
secondary mode and the source have the same zeroth order frequency
$2\omega_{n}^{0}$, so their difference is due to the corrections given by
$\omega_{n}^{1}$ in equation (46). Equation (54) tells us that the amplitude
of this secondary mode is expected to be much larger than other secondary
modes due to the factor $4\omega_{p}^{2}-\omega_{j}^{2}$ in the denominator.
Let us calculate $a_{j}$ corresponding to this mode. For that, we take into
account the difference between these frequencies only in the denominator, in
all other expressions we neglect the corrections $\omega_{n}^{1}$, and assume
that the eigenfunction of the secondary mode is given by (41) with
$\phi=\pi/6$ multiplied by $e^{2ik_{p}x}$. We also set the amplitude $C$ of
the secondary eigenfunction to unity, since our final expressions do not
depend on its value, and we omit the indices enumerating different secondary
modes hereafter, assigning index $s$ to quantities belonging to the secondary
mode where this matters.

First, we calculate the difference $4\omega_{p}^{2}-\omega_{j}^{2}\equiv
4\omega_{p}^{2}-\omega_{s}^{2}$. Using equation (46), we obtain

$\displaystyle 4\omega_{p}^{2}-\omega_{s}^{2}={12{\omega^{0}_{p}}^{2}\over
z_{r}}\nu,$ (55)

where

$\displaystyle\nu$ $\displaystyle=$
$\displaystyle{3^{1/6}\over{2\sqrt{z_{r}}}}{\Gamma(1/3)\over\Gamma(2/3)}\left(\kappa(2k_{p})-\kappa(k_{p})\right){\left({\omega_{p}^{0}\over
k_{p}}\right)}^{5/3}+3{\omega^{0}_{p}}^{2},$

and where $\kappa(k)$ is given by (40) and all quantities on the right hand
side depend on $\omega^{0}_{p}$ and $k_{p}$. From (40) we have
$\kappa(2k_{p})-\kappa(k_{p})=k_{p}(1-e^{-2k_{p}z_{c}})/(1+e^{-2k_{p}z_{c}})$.
This relation tells us that the correction proportional to $\lambda_{*}^{2}$
tends to zero when $z_{c}k_{p}\rightarrow 0$. In all other expressions below
we can use $\omega^{0}_{p}$ and $k_{p}$, therefore, for simplicity, we set
$\omega\equiv\omega^{0}_{p}$ and $k\equiv k_{p}$ from now on. We substitute
the expression for the source term (47) and (55) in (54) and remember that
$\lambda_{j}$ there should be equal to $\omega_{p}^{2}$. We obtain

$a_{s}={2iz_{r}\over
3\omega\nu}{\int^{z_{r}}_{0}{\bar{v}_{p}}^{3}dz\over\int^{z_{r}}_{0}z{\bar{v}_{p}}^{2}dz}C_{p}^{2}e^{2i\omega
t}.$ (56)

Note that since $\phi^{*}_{s}\propto e^{-2ik_{p}x}$ and $S_{p}\propto
e^{2ik_{p}x}$ the integrands do not depend on $x$, so integration over $x$ is
trivial. In order to evaluate the integrals entering (56) we neglect
$\omega^{2}$ in the expression for $y$ in terms of $z$ (see (41)), and change
the integration variable from $z$ to $y$. Since the integral in the numerator
converges when $z_{r}/\lambda_{*}\rightarrow\infty$, its upper limit of
integration can be extended to infinity. Numerical evaluation of
$\int^{\infty}_{0}{\bar{v}_{p}}^{3}dy$ shows that its value is close to 0.2.
On the other hand the integral in the denominator diverges as
$(z_{r}/\lambda_{*})^{3/2}$ in the same limit. Therefore, we use the limiting
value of the ratio
$\int^{z_{r}/\lambda_{*}}_{0}y{\bar{v}_{p}}^{2}/(z_{r}/\lambda_{*})^{3/2}dy$,
which is close to $0.1$ to obtain
${\int^{z_{r}}_{0}{\bar{v}_{p}}^{3}dz\over\int^{z_{r}}_{0}z{\bar{v}_{p}}^{2}dz}\approx
2{\lambda_{*}^{1/2}\over z_{r}^{3/2}}$. Substituting this in (56) we get

$a_{s}\approx{4i\over 3}{1\over
z_{r}^{1/2}k^{1/3}\omega^{2/3}\nu}C_{p}^{2}e^{2i\omega t}.$ (57)

We expect nonlinear behaviour to be important when $|a_{s}|\gtrsim C_{p}$, and
potentially even for amplitudes quite a bit smaller than this. Equation (57)
tells that the wave amplitude $C_{p}$ in this case should exceed its critical
value

$C_{p}^{crit}={3\over 4}z_{r}^{1/2}k^{1/3}\omega_{p}^{2/3}\nu,$ (58)

for these nonlinearities to be important. Above this amplitude, the secondary
waves have amplitudes that are comparable with or exceed the primary wave
amplitude.

##  6 Relation to the tidal problem and criterion for the transition to
nonlinearity

###  6.1 Relation of our local model results to a spherical star

In the previous section we used dimensionless units in which spatial scales
were expressed in terms of the radius of the base of the convection zone,
$r_{c}$, and temporal ones were expressed in terms of the characteristic
frequency $\omega_{c}=\sqrt{A_{c}r_{c}}$, where $A_{c}$ was defined in eq.
(26). Bearing this in mind, we can relate results obtained from the local
model in the two previous sections to the global normal modes of a sun-like
star. For that we compare (41) with equation (91), and (46) with equations
(103-106) of IPCh, respectively. From the comparison of the first pair of
equations, we see that they coincide when the term proportional to
$\omega^{2}$ is neglected in the argument of the Airy functions in (41) and we
set $k=\sqrt{\Lambda}$, where $\Lambda$ is defined in IPCh as an eigenvalue of
Laplace’s tidal equation, and it can be expressed as $\Lambda=l(l+1)$, where
$l$ is the spherical harmonic degree for a non-rotating star. Note that the
term $\propto\omega^{2}$ was neglected in IPCh because of its smallness, but
it is, however, important for our purposes, since it gives a potentially
important correction in the expression (58) for the quantity $\nu$. Comparing
the expression for $\omega_{n}^{0}$ with eq. (103) of IPCh, we see that they
are equivalent to each other provided we make the following redefinitions in
(46): $n\rightarrow n+2$ and $z_{r}\rightarrow{\left({3\over
2}{I\over\omega_{c}}\right)}^{2/3}$, where $I$ is defined in (104) of IPCh.
Comparing eq. (105) of IPCh and the term proportional to the factor $\kappa$
in the expression for $\nu$ in (58), these equations formally lead to the same
frequency correction when $\kappa=r_{c}B_{c}$, where the quantity $B_{c}$
defined in equation (100) of IPCh. Note, that the numerical values are, of
course, different, since $B_{c}$ was calculated in spherical geometry and for
a particular size of convective envelope.

It is very important to stress that the planar geometry used in the previous
sections can be shown to be fully equivalent to the spherical one only in the
limit $l\rightarrow\infty$. For the tidal problem considered in this paper,
however, $l=2$ for the primary mode. On the other hand, for the secondary mode
$l$ should be equal to $4$, which results in an additional correction in the
expression for the quantity $\nu$ characterising the frequency difference
between primary and secondary modes. In order to calculate it, we start from
the expression for the eigenfrequencies $\omega_{0,n}$ calculated in IPCh in
the WKBJ approximation, and given in their equation (103):

$\omega_{0,n}={\sqrt{l(l+1)}I\over\pi(n+(l+1)/2+1/12)}.$ (59)

Note that for simplicity we neglect their correction
$\delta\phi(\omega_{0,n})$, since it is already taken into account in the
expression for $\nu$ discussed above. We assume $n\gg 1$ for the validity of
the WKBJ approximation, therefore the terms in the denominator can be
considered as small corrections. We can discard them altogether if we note
that when $l=2$ and $4$ for the primary and secondary modes, they differ by
one and this difference can be absorbed into a redefinition of the mode number
$n$. The rest don’t give any contribution at leading order in the small
parameter $\Delta_{l}=\sqrt{{6\over 5}}-1\sim 0.095$ arising from the
difference between the values $2\sqrt{l(l+1)}$ and $\sqrt{l(l+1)}$ in the
numerator in (59) for $l=2$ and $4$, respectively. Taking this into account we
can use

$\omega_{0,n}={\sqrt{l(l+1)}I\over\pi n}.$ (60)

when calculating the frequency difference $2\omega_{p}-\omega_{s}$ and,
accordingly, the additional correction to $\nu$ setting $l=2$ and $4$ for
$\omega_{p}$ and $\omega_{s}$, respectively. For the primary mode we set
$n=n_{p}$ and for the secondary mode $n=n_{s}\equiv n_{p}-k$, where it is
assumed that $k\ll n_{p}$, and $k$ should be chosen in such a way that the
absolute value of $2\omega_{p}-\omega_{s}$ is minimised. Taking into account
only the leading terms in $\Delta_{l}$, we obtain

$2\omega_{p}-\omega_{s}={2\sqrt{6}I\over\pi
n_{p}}f(n_{p})=2\omega_{p}f(n_{p}),$ (61)

where

$f(n)=\textrm{min}_{k}\left|\sqrt{{6\over 5}}-1-\frac{k}{n}\right|.$ (62)

We show the dependency of $f(n)$ on $n$ in Fig. 1.

Figure 1: The function $f(n)$ defined in eq. (62) is shown as a function of
$n$.

Comparing (61) with the expression (55) for $\nu$ and remembering that
$z_{c}={({3\over 2}{I\over\omega_{c}})}^{2/3}$, we see that $\nu$ should
contain the additional term

${\left({2\over
3}\right)}^{1/3}{\left({I\over\omega_{c}}\right)}^{2/3}f(n_{p}).$ (63)

As we mentioned above for the tidal problem $\omega_{p}\approx 2\Omega_{orb}$.
Using (60) with $l=2$ we can express $n_{p}$ in terms of $\Omega_{orb}$ and
consider $f(n)$ to be a function of the orbital frequency. We have

$n_{p}=\left[{\sqrt{6}I\over 2\pi\Omega_{orb}}\right],$ (64)

where $[.]$ implies that only the integer part of the expression is used.

Restoring physical units, using the redefinitions mentioned above, and
substituting (63) in (58), we obtain from (55) and (58):

$\displaystyle\nu$ $\displaystyle=$
$\displaystyle{3^{\frac{1}{6}}\Gamma(\frac{1}{3})\over 2\cdot
6^{\frac{5}{6}}\Gamma\left(\frac{2}{3}\right)}{\left({{2\over
3}{\omega_{c}\over
I}}\right)}^{\frac{1}{3}}B{\left({\omega_{p}\over\omega_{c}}\right)}^{\frac{5}{3}}+3{\left({\omega_{p}\over\omega_{c}}\right)}^{2}$
(65) $\displaystyle\hskip 28.45274pt+{\left({2\over
3}\right)}^{\frac{1}{3}}{\left({I\over\omega_{c}}\right)}^{\frac{2}{3}}f(\Omega_{orb}),$

where $B=r_{c}(B_{c}(2\cdot 6^{\frac{1}{2}})-B_{c}(6^{\frac{1}{2}}))$, and

$C_{p}^{crit}=3^{\frac{3}{2}}2^{-\frac{13}{6}}{\left({\omega_{p}\over\omega_{c}}\right)}^{\frac{2}{3}}{\left({I\over\omega_{c}}\right)}^{\frac{1}{3}}\nu\omega_{c}r_{c}.$
(66)

###  6.2 The criterion for transition to nonlinearity for a single secondary
mode

The above expression (66) should be compared with the amplitude of the radial
component of the linear tidal response. This amplitude,
$U^{z}=i\omega_{p}{\mbox{\boldmath${\xi}$}}^{r}$ can be calculated using (12),
where $\omega_{p}=2\Omega_{orb}$ and we should take the limit $r\rightarrow
r_{c}$ when evaluating $\xi$. The corresponding expression is given by
equation (101) of IPCh:

$\xi\approx{(-1)^{n}\over\Gamma({2\over
3})}\sqrt{\pi}{\left({3\omega_{p}\over\omega_{c}}\right)}^{-\frac{1}{6}}{6^{\frac{1}{12}}\over\omega_{c}^{\frac{1}{2}}\rho_{c}^{\frac{1}{2}}r_{c}^{\frac{3}{2}}}C_{WKBJ}.$
(67)

Using equation (20) for $C_{WKBJ}$ we obtain from (67)

$\xi\approx{(-1)^{n}\over\Gamma({2\over
3})}3^{-\frac{1}{12}}2^{\frac{7}{12}}\sqrt{\pi}{\left({\omega_{p}\over\omega_{c}}\right)}^{\frac{5}{6}}{\left({\omega_{c}\over
I}\right)}^{\frac{1}{2}}{n^{\frac{1}{2}}\over\rho^{\frac{1}{2}}_{c}r_{c}^{\frac{3}{2}}}.$
(68)

We substitute (68) in (12), use $\omega_{p}=2\Omega_{orb}$, and multiply the
result by $2\Omega_{orb}$, to obtain

$\displaystyle U^{z}$ $\displaystyle=$ $\displaystyle
2i\Omega_{orb}{\mbox{\boldmath${\xi}$}}^{r}(r_{c})$ (69) $\displaystyle=$
$\displaystyle i{(-1)^{n}\over\Gamma({2\over
3})}3^{\frac{11}{12}}2^{-\frac{19}{12}}\pi^{\frac{3}{2}}\left({q\over
1+q}\right)$ $\displaystyle\hskip
17.07182pt\times\left({\omega_{c}\over|d\omega/dj|}\right){\left({\omega_{c}\over
I}\right)}^{\frac{1}{2}}{\left({\Omega_{orb}\over\omega_{c}}\right)}^{\frac{17}{6}}\hat{Q}{\omega_{c}\over\rho_{c}^{\frac{1}{2}}r_{c}^{\frac{3}{2}}}.$

Note that since $\hat{Q}$ has dimensions $g^{1/2}\cdot\mathrm{cm}$, equation
(69) has the correct dimensions $\mathrm{cm/s}$. When (69) is larger than than
(66) we assume that a fully non-linear regime sets in. From the condition
$|U^{z}|>C_{p}^{crit}$ we have

$\displaystyle{q\over 1+q}>C_{crit,c}$ $\displaystyle=$ $\displaystyle
3^{\frac{7}{12}}2^{\frac{1}{12}}\pi^{-\frac{3}{2}}\Gamma\left(\frac{2}{3}\right){|d\omega/dj|\over\omega_{c}}$
(70)
$\displaystyle\times{\left({I\over\omega_{c}}\right)}^{\frac{5}{6}}{\left({\Omega_{orb}\over\omega_{c}}\right)}^{-\frac{13}{6}}\nu{\rho_{c}^{\frac{1}{2}}r_{c}^{\frac{5}{2}}\over\hat{Q}}.$

Finally, we use (21) to obtain

$C_{crit,c}=3^{\frac{1}{12}}2^{\frac{19}{12}}\pi^{-\frac{1}{2}}\Gamma\left(\frac{2}{3}\right){\left({I\over\omega_{c}}\right)}^{-\frac{1}{6}}{\left({\Omega_{orb}\over\omega_{c}}\right)}^{-\frac{1}{6}}\nu{\rho_{c}^{\frac{1}{2}}r_{c}^{\frac{5}{2}}\over\hat{Q}}.$
(71)

This is an amplitude criterion for our secondary super-harmonic wave to have a
comparable (or larger) amplitude to the primary tidal wave, above which we
expect consideration of nonlinear effects to be essential.

###  6.3 The criterion for transition to nonlinearity with a dense spectrum of
secondary modes

If the spectrum of eigenmodes is sufficiently dense, the primary mode can
effectively excite a number of neighbouring secondary modes instead of just
the one having a frequency that most closely matches the primary frequency. As
discussed in Section 2, in order for the spectrum to be sufficiently dense,
either the ratio $\delta$ of the value of the frequency offset (proportional
to the quantity $\nu$ defined in (55)) to the distance between neighbouring
eigenfrequencies, $d\omega_{j}/dj$, is large, or the inverse decay time,
$\omega_{\nu}$, is larger than $d\omega_{j}/dj$. While the former condition
can be straightforwardly calculated within the framework of our model, the
latter condition assumes that the wave packet composed of secondary waves with
approximately the same frequency decays during its travel across a star. Since
this has been assumed for the primary wave, it is reasonable to make the same
assumption for these secondary waves as well. In both cases, the modified
values of $C_{crit,c}$ are supposed to be smaller than the values in (71), so
we can consider (71) as a conservative estimate of $C_{crit}$. On the other
hand, our estimate here can be considered a more ‘optimistic’ estimate,
suggesting what can in principle be achieved from the effect discussed in this
paper.

The quantity $\delta$ defined below eq. (7) can be calculated using (21) and
(55). From equation (55), it follows that the distance between the doubled
primary and secondary frequencies, $\Delta\omega_{s}=2\omega_{p}-\omega_{s}$,
can be expressed as

$\Delta\omega_{s}=3^{\frac{1}{3}}2^{\frac{5}{3}}{\left({\omega_{c}\over
I}\right)}^{\frac{2}{3}}\nu\Omega_{orb},$ (72)

where we remember that $z_{r}={({3\over 2}{I\over\omega_{c}})}^{\frac{2}{3}}$
and $\omega_{p}=2\Omega_{orb}$. Now we note that there should be $2\sqrt{6}$
instead of $\sqrt{6}$ in the denominator of (21), see also eq. (125) of IPCh,
and we divide (72) by this expression to obtain

$\delta={2^{\frac{7}{6}}\cdot
3^{\frac{5}{6}}\over\pi}{\left({I\over\omega_{c}}\right)}^{\frac{1}{3}}\nu{\omega_{c}\over\Omega_{orb}}.$
(73)

From our discussion in Section 2, it follows that when $\delta$ and/or
$\kappa$ are large $\Delta\omega_{s}$ in the expressions for the amplitude of
excited modes, as in e.g. (56), should be substituted by
$({d\omega_{s}/dj})(S(\delta,\kappa))^{-1}$, where $S(\delta,\kappa)$ is given
by (8) and we remember that
$d\omega_{s}/dj={2\pi\Omega^{2}_{orb}\over\sqrt{6}I}$ for secondary modes.
Since $\nu$ is proportional to $\Delta\omega_{s}$ it should be changed
accordingly in all expressions, including in (71). Therefore, in order to
account for the excitation of many near-resonant secondary modes we can
express $\nu$ in terms of $\Delta\omega_{s}$ using (72), and make the
substitution

$\nu\rightarrow{\pi\over 3^{\frac{5}{6}}\cdot
2^{\frac{7}{6}}}|S|^{-1}{\Omega_{orb}\over
I^{\frac{1}{3}}\omega_{c}^{\frac{2}{3}}}$ (74)

in (71), where we use the absolute value of $S$, since it can be a complex
quantity. In this way, we obtain

$C^{dense}_{crit,c}=3^{-\frac{3}{4}}2^{\frac{5}{12}}\pi^{\frac{1}{2}}\Gamma\left({2\over
3}\right){\left({I\over\omega_{c}}\right)}^{-\frac{1}{2}}|S|^{-1}{\left({\Omega_{orb}\over\omega_{c}}\right)}^{\frac{5}{6}}{\rho_{c}^{\frac{1}{2}}r_{c}^{\frac{5}{2}}\over\hat{Q}}.$
(75)

It is very important to stress here, that $C^{dense}_{crit}$ given by (75)
should be considered as a qualitative indication of the importance of having a
dense spectrum of secondary waves rather than a quantitative criterion, since
the derivation of this expression is based on a number of simplifying
assumptions. In particular, we assume that the summation of the series in (7)
can be formally extended to infinity.

##  7 Numerical study of weakly nonlinear excitation of super-harmonic
secondary waves

In this section we briefly present some numerical calculations to verify the
analytical results obtained in § 4–5. We first solve numerically the
generalised eigenvalue problem for time-harmonic solutions proportional to
$\mathrm{e}^{i(\omega t+k_{p}x)}$ to the linearised system of equations

$\displaystyle i\omega{\bf U}$ $\displaystyle=$ $\displaystyle-\nabla P+b{\bf
e}_{z}-\gamma{\bf U},$ (76) $\displaystyle i\omega b$ $\displaystyle=$
$\displaystyle-N^{2}(z)U^{z},$ (77) $\displaystyle\quad\nabla\cdot{\bf U}$
$\displaystyle=$ $\displaystyle 0,$ (78)

where $\omega\in\mathbb{C}$ is the eigenvalue,

$\displaystyle N^{2}(z)=\begin{cases}z,\quad\;\;0\leq z\leq z_{r},\\\ 0,\quad-
z_{c}\leq z<0,\end{cases}$ (79)

subject to $U^{z}(z=-z_{c})=U^{z}(z=z_{r})=0$. We have introduced a frictional
damping term $-\gamma{\bf U}$ where $\gamma=10^{-4}$ is a constant (analogous
to $\omega_{\nu}$ in section 2, which is simpler numerically than including
viscosity or thermal diffusion). This problem is discretised in $z$ using a
Chebyshev collocation method with $\mathcal{N}+1$ points (Boyd, 2001), with
$\mathcal{N}=200$, and we choose $k_{p}\equiv k=z_{c}=z_{r}=1$ to define our
choice of space and time units (noting that we then have
$\mathrm{max}[N(z)]=1$). The solution of this eigenvalue problem gives us a
set of eigenvalues $\\{\omega_{i}\\}$ with corresponding eigenvectors
$\\{({\bf U},P,b)_{i}\\}$. We order the eigenvalues
$\omega_{i}\leq\mathrm{max}[N(z)]$ in descending order using their real part
since the modes with the largest frequencies are those that are best resolved
numerically.

We then choose a single mode $i$ as our “primary mode”, set
$\omega_{p}=\mathrm{Re}[\omega_{i}]$, and substitute the corresponding
eigenfunction $\bar{v}_{p}(z)=U^{z}_{i}$ (without loss of generality, this is
normalised so that it is real with a maximum magnitude that matches (41)) into
the right hand side of equation (47), to which we also add the additional
frictional damping term $\gamma\Delta\dot{U}^{z}$ on the left hand side
(equivalent to the term $-\gamma{\bf U}$ in 76). We solve this equation for
the “secondary modes” as an initial value problem with boundary conditions
$U^{z}(z=-z_{c},t)=U^{z}(z=z_{r},t)=0$, initialising our solution with
$U^{z}(z,t=0)=0$ so that the forcing excites all of the modes until they are
damped by friction. To do this, we use Chebyshev collocation in $z$ with
$\mathcal{N}+1$ points (again with $\mathcal{N}=200$) and a second-order
central difference scheme for time integration. This equation is integrated
until a time $t=5000$ and the frequency power spectrum ($|\hat{U}^{z}|^{2}$)
is computed using a discrete Fourier transform of the signal at the location
$z=0.1$.

Figure 2: Mode spectrum showing the real and imaginary parts of the
eigenfrequencies $\omega_{i}$ for $k=1$ (black crosses) and $k=2$ (blue
circles) with $z_{c}=z_{r}=1,\gamma=10^{-4},\mathcal{N}=200$. The cyan dashed
lines are the analytical predictions for $k=1$ using 46, which agree very well
for low frequencies. The red star highlights the mode with
$\omega_{p}=\mathrm{Re}[\omega_{i}]\approx 0.0456$.

(a) Primary $U^{z}$ profile

(b) Secondary $U^{z}$ vs $z$ at $t=5000$.

(c) Secondary $U^{z}$ vs $t$ at $z=0.1$.

(d) Frequency power spectrum of $U^{z}(z=0.1,t)$.

Figure 3: (a): Comparison of (normalised) primary mode eigenfunction $U^{z}$
with the analytical prediction in (41) for the mode with
$\omega_{p}=\mathrm{Re}[\omega_{i}]\approx 0.0456$. (b): Real (solid black
line) and imaginary (dashed blue line) parts of secondary mode $U^{z}$ as a
function of $z$ at $t=5000$, compared with the prediction (red dashed line)
for $a_{s}$ from (57). (c): Real (solid black line) and imaginary (dashed blue
line) parts of secondary mode $U^{z}$ as a function of $t$ at $z=0.1$ in the
stable layer. (d): Frequency power spectrum of $U^{z}$, $|\hat{U}^{z}|^{2}$,
at $z=0.1$ as a function of angular frequency $\tilde{\omega}$. Otherwise same
parameters as Fig. 2.

(a) Primary $U^{z}$ profile

(b) Secondary $U^{z}$ vs $z$ at $t=5000$.

(c) Secondary $U^{z}$ vs $t$ at $z=0.1$.

(d) Frequency power spectrum of $U^{z}(z=0.1,t)$.

Figure 4: Same as Fig. 3 but for a mode with
$\omega_{p}=\mathrm{Re}[\omega_{i}]\approx 0.263$.

(a) Primary $U^{z}$ profile

(b) Secondary $U^{z}$ vs $z$ at $t=5000$.

(c) Secondary $U^{z}$ vs $t$ at $z=0.1$.

(d) Frequency power spectrum of $U^{z}(z=0.1,t)$.

Figure 5: Same as Fig. 3 but for a mode with
$\omega_{p}=\mathrm{Re}[\omega_{i}]\approx 0.0220$.

We first show the spectrum of $\omega_{i}$ from solving the above eigenvalue
problem in the complex plane in Fig. 2, where we highlight a particular low
frequency mode with $\omega_{p}\approx 0.0456$ with a red star. We show the
analytical predictions for $\omega_{i}$ from (46) as the vertical cyan dashed
lines, and our numerical values for the complex mode frequencies for the modes
with $k=1$ as black crosses and $k=2$ as blue circles (the latter shows the
“secondary” modes that can be excited). The agreement is excellent between
(46) and our numerical results for the $k=1$ modes for
$\mathrm{Re}[\omega_{i}]\lesssim 0.1$, with a slight departure, as expected,
when the frequency is larger.

In Fig. 3(a), we show the $U^{z}$ eigenfunction for the particular primary
mode with $\omega_{p}\approx 0.0456$ together with the analytical prediction
in (41) in the stable layer, again indicating excellent agreement. We then
show in Figs. 3(b) –3(d) the spatial structure of the secondary solution at
$t=5000$ (the real and imaginary parts shown as solid and dashed lines,
respectively), the temporal evolution of the secondary $U^{z}$ at $z=0.1$, and
the frequency power spectrum ($|\hat{U}^{z}|^{2}$) of the secondary signal
$U^{z}(z=0.1,t)$. By $t=5000$, modes with multiple frequencies still
contribute to the signal, though the power is maximised in the mode with a
frequency close to $2\omega_{p}$. In Fig. 3(d), we also show the predicted
frequencies of the secondary modes $\omega_{s}$ from solving the above
eigenvalue problem in the case $k=2$, which indicates that many of these modes
are excited and have not yet fully damped, as we expect from our initial
conditions. Figs. 3(b) –3(d) show that the analytical prediction for the
secondary mode amplitude $a_{s}$ in (57), which is shown as the red horizontal
lines in Figs. 3(b) and 3(c), agrees quite well with our numerical results in
this case. The secondary mode amplitudes are dominant near $z\sim 0$, where
the primary mode also takes its maximum value. We observe here that, due to
the near-resonance explained in § 5, the amplitude $U^{z}$ of the super-
harmonic secondary solution can far exceed that of the primary mode, and hence
we predict this nonlinear effect to be important in this system because
$a_{s}\gg C_{p}$.

We also show results for a larger frequency mode with $\omega_{p}\approx
0.263$ in Figs. 4(a)–4(d) and a lower frequency mode with $\omega_{p}\approx
0.0220$ in Figs. 5(a)–5(d), with otherwise the same parameters. The primary
wave solution (and its frequency) for the larger frequency case is clearly
described less well by our analytical prediction, but the agreement is still
good. The amplitude predicted for $a_{s}$ is also comparable with the values
we observe numerically in both of these examples.

Our numerical calculations here have thus verified the analytical results
derived for low frequencies in § 4–5, so we can be more confident in their
application. It would be interesting to explore the fully nonlinear evolution
of the system studied here using numerical simulations to analyse the long-
term dynamics of the generation of super-harmonic secondary waves and their
interaction with the primary (tidal) wave. These simulations will be presented
in a future publication.

##  8 Application of the criteria for transition to nonlinearity to main-
sequence stellar models

In this section, we turn to apply our criteria for the transition to
nonlinearity in several main-sequence stellar models computed using MESA
version 15140666See http://mesa.sourceforge.net. (Paxton et al., 2011, 2013,
2015, 2018, 2019). We calculate the values of $C^{crit}$ for each of the
various scenarios given by eqns. (22), (25), (71) and (75) for a number of
main-sequence stellar models with masses $M=1$, $1.4$ and $2M_{\odot}$ and
different ages. The basic characteristics of our models are given in table 1,
where we show stellar masses in units of $M_{\odot}$, their radii in units of
$R_{\odot}$, ages in years, positions of the base of the convective zone,
$r_{c}$, in units of the stellar radius $R_{*}$, density at the base of
convective zone, $\rho_{c}$, in units of the mean density
$\bar{\rho}=3M_{\ast}/(4\pi R_{\ast}^{3})$, and the mean density in units of
the solar mean density, and finally the values of $\omega_{c}$ and $I$ in
units of $\sqrt{GM_{\ast}\over R_{\ast}^{3}}$, respectively. Additionally, we
show radial profiles of the Brunt-Väisälä frequency and stellar density in
Figs. 6, 7 and 8, respectively.

Table 1: Basic characteristics of our stellar models Model | Mass ($M_{\odot}$) | Radius ($R_{\odot}$) | Age (years) |  $r_{c}$ ($R_{\ast}$) | $\rho_{c}/\bar{\rho}$ | $\bar{\rho}$ | $\omega_{c}$ | $I$
---|---|---|---|---|---|---|---|---
1a | 1 | 0.96 | $2\times 10^{9}$ | 0.73 | $0.12$ | 1.13 | 9.36 | 9.42
1.4a | 1.4 | 1.65 | $1.07\times 10^{7}$ | 0.68 | $0.32$ | 0.31 | 9.01 | 8.52
1.4b | 1.4 | 1.43 | $2.03\times 10^{8}$ | 0.97 | $2.05\times 10^{-5}$ | 0.48 | 49.2 | 10.3
1.4c | 1.4 | 1.72 | $1.73\times 10^{9}$ | 0.94 | $1.89\times 10^{-4}$ | 0.28 | 30.7 | 14.3
1.4d | 1.4 | 1.93 | $2.38\times 10^{9}$ | 0.87 | $2.59\times 10^{-3}$ | 0.19 | 20.5 | 19.0
2a | 2 | 1.63 | $2.67\times 10^{7}$ | 0.99 | $1.78\times 10^{-7}$ | 0.46 | 182.3 | 8.07
2b | 2 | 2.68 | $9.97\times 10^{8}$ | 0.99 | $5.94\times 10^{-7}$ | 0.10 | 113.4 | 16.7
2c | 2 | 5.35 | $1.41\times 10^{9}$ | 0.96 | $3.02\times 10^{-5}$ | 0.013 | 41.6 | 52.5
2d | 2 | 5.75 | $1.44\times 10^{9}$ | 0.89 | $6.09\times 10^{-4}$ | 0.011 | 17.7 | 59.9
Figure 6: Radial profiles of stellar density $\rho$ (black solid line) and
square of the Brunt-Väisälä frequency $N^{2}$ (red dashed line), both in the
natural units defined in § 2, for model 1a. Figure 7: Radial profiles of
$\rho$ and $N^{2}$ in natural units for the models 1.4a, 1.4b, 1.4c and 1.4d
shown as black solid, red dashed, green dotted and red dot-dashed lines,
respectively. Although the lines corresponding to $\rho$ and $N^{2}$
calculated for the same model have the same style, they can be easily
distinguished from each other, since $\rho$ decreases monotonically with $R$.
Figure 8: Same as Fig. 7, but for the models 2a, 2b, 2c and 2d shown as black
solid, red dashed, green dotted and blue dot-dashed lines, respectively.

In Fig. (9) we show the quantities $C_{centre}$, $C_{crit,c}$, and
$C_{crit,c}^{dense}$ for both $\kappa=0$ and $\kappa=1$, as given by equations
(22), (71) and(75), respectively, as functions of orbital period for the solar
model. We observe that according to these criteria, the transition to non-
linearity at the stellar centre always happens for much smaller mass ratios
than any of the other measures i.e. $C_{centre}\ll
C_{crit,c},C_{crit,c}^{dense}$. This criterion is equivalent to the condition
for nonlinearity discussed by (e.g. Goodman and Dickson, 1998; Ogilvie and
Lin, 2007; Barker and Ogilvie, 2010; Barker, 2020), and is the most important
measure of nonlinearity in our solar-mass model.

Note that in order for our criteria to predict the onset of nonlinearity, due
to its definition we require $C_{crit}\lesssim 1$ for any mass ratio. From
Fig. (22), we see that $C_{crit,c}$ given by (71) becomes smaller than one
only for very short orbital periods $P_{orb}\sim 0.2\mathrm{d}$ and the
corresponding mass ratio should be larger than $\sim 0.6$ i.e. nonlinearity is
predicted only for comparable mass ultra-close binaries. Assuming that when
our criteria for the transition to non-linearity are satisfied the primary
mode is then efficiently damped, we conclude here that generation of super-
harmonic secondary modes near the radiative-convective interface in our solar
model is unlikely to provide an important contribution to tidal dissipation.

It is interesting to note that both $C_{crit,c}$ and
$C_{crit,c}^{dense}(\kappa=0)$ are oscillatory functions of orbital period. In
the former case, this is due to the oscillatory character of the function
$f(n)$ shown in Fig. 1. In the latter case, the behaviour is even more
strongly oscillatory, which is due both to the oscillation of $f(n)$ with $n$
and to the sharp periodic changes of $S(\delta,\kappa)$ with $\delta$ when
$\kappa=0$, as indicated by eq. (8). These curves have the same oscillatory
behaviour in the more massive stellar models considered below.

Figure 9: Various criteria for transition to nonlinearity $C_{crit}$ are shown
for a solar mass model as a function of orbital period in days. The black
solid curve corresponds to $C_{centre}$ using the ‘standard’ expression (22)
for wave breaking in radiative cores, the red dashed one corresponds to
$C_{crit,c}$ in (71), the green dotted and blue dot-dashed curves are given by
$C_{crit,c}^{dense}$ in eq. (75), where in the former we set $\kappa=0$ and in
the latter $\kappa=1$.

Figure 10: Same as Fig. 9, but for our models with $M=1.4M_{\odot}$, and the
black solid line shows (25). Top left, top right, bottom left and bottom right
panels correspond to the models a, b, c, d, respectively.

Figure 11: Same as Fig. 10, but for the models with $M=2M_{\odot}$. Figure 12:
Various expressions for $C_{crit}$ as a function of stellar age in years for
stellar models with $M=2M_{\odot}$. Solid curve, dashed, dotted and dot-dashed
curves correspond to (25), (71), (75) with $\kappa=0$ and (75) with
$\kappa=1$, respectively.

Figs. (10) and (11) are similar to Fig. (9), but the curves are calculated for
more massive models with $M=1.4M_{\odot}$ and $M=2M_{\odot}$. Contrary to our
solar model, apart from model 1.4a the more massive stars have convective
cores, as is shown in Figs. 7 and 8. In this case, the criterion for wave
breaking near convective cores predicted by (25) must be used instead of the
criterion for wave breaking at the centre of a radiative core in (22). The
corresponding curve is shown as a solid line in all figures apart from the one
corresponding to 1.4a, where (22) is shown in the same way instead. Unlike the
solar case we see from (10) and (11) that the criteria for transition to
nonlinear behaviour at the outer convective/radiative interface is satisfied
for smaller mass ratios than the criterion for transition at the convective
core (25), for almost all stellar ages and orbital periods777We ignore the
very large amplitude oscillations of $C^{dense}_{crit,c}$ due to the
trigonometric nature of $S(\delta,\kappa=0)$ and assume that only some average
value of this quantity has physical meaning.. The exceptions are models 1.4a,
1.4d and 2d. In the case of 1.4a, we see from Fig. 7 that this model has a
radiative core and an extended convective envelope. This model is rather
similar to the solar model and it is not surprising that it demonstrates
similar behaviour. The models 1.4d and 2d also have quite large convective
envelopes. The larger the transition radius to the convective envelope,
$r_{c}$, the smaller is the density at this radius, $\rho_{c}$, and the larger
is typical frequency, $\omega_{c}$, at this radius. Both the smaller density
and larger $\omega_{c}$ values lead to the smaller values of $C_{crit,c}$ seen
in most of these models. Note, however, that there are two other important
factors determining values of $C_{crit,c}$, namely, values of the overlap
integrals in units of the frequency $\Omega_{*}=\sqrt{GM_{*}\over R_{*}^{3}}$,
and a value for the “average frequency” of a star. The former factor leads to
smaller $C_{crit,c}$ for stars with extended convective envelopes, since the
overlap integrals are typically larger for such stars (e.g. Chernov et al.,
2013, 2017). Smaller values of the average density lead to larger values of
the tidal forcing frequency expressed in units of the natural frequency and,
as a result, smaller values of $C_{crit,c}$. It may be possible that these two
factors counterbalance each other for certain stellar models however.

If we speculate that whenever our criteria for transition to nonlinear
behaviour are satisfied, the primary tidally-excited modes are efficiently
damped, our results concerning more massive models could have important
implications for tidal evolution around such stars. They would imply that
above a critical mass ratio, tidal dissipation rates can be straightforwardly
computed as a known function of orbital period, stellar mass and age
(essentially according to linear tidal theory assuming a fully
damped/MLD/travelling wave regime). In some cases, this ratio can be as small
as $\sim 10^{-3}$ for $P_{orb}\sim 1\mathrm{d}$, as in e.g. models 2a-2c,
indicating that planetary mass companions can cause a strongly nonlinear tidal
response in the star, potentially with efficient tidal dissipation. To
illustrate this further, we show the dependence of $C_{crit,core}$,
$C_{crit,c}$ and $C_{crit,dense}$ on stellar age in Fig. 12, keeping the
orbital period fixed at one day. As shown in this figure, typical values of
$C_{crit,dense}$ are indeed as small as $10^{-3}$ for stellar ages $\sim
10^{9}\mathrm{yr}$, indicating that this mechanism can potentially be
important for tidal evolution of hot Jupiters around such stars.

##  9 Conclusions and Discussion

We have studied the role of nonlinear effects on tidally-excited gravity waves
in the radiation zones of stars, primarily focussing on a new mechanism that
could be important in stars possessing convective cores. Our work was partly
motivated to study tides due to massive short-period hot Jupiters, which are
observed to preferentially orbit stars with convective cores (e.g. WASP-18 b
Wilkins et al., 2017). For these stars, the geometric focussing and consequent
breaking (or nonlinear wave-wave interactions) of gravity waves in the stellar
core (which we have revisited in § 3.1), which can result in efficient tidal
dissipation when this occurs (e.g. Goodman and Dickson, 1998; Ogilvie and Lin,
2007; Barker and Ogilvie, 2010; Barker, 2011; Barker and Ogilvie, 2011;
Weinberg et al., 2012; Essick and Weinberg, 2016; Barker, 2020), cannot take
place for planetary mass companions, unlike in solar-type with radiative
cores.

We have developed a theory for the nonlinear excitation of super-harmonic
‘secondary’ internal gravity waves by a ‘primary’ gravity wave, assuming that
the latter is generated by tidal forcing in a main-sequence star due to an
orbiting companion. This excitation appears to be the most efficient near the
transition radius, $r_{c}$, between a radiative interior and a convective
envelope. At this location, the usual WKBJ approximation for the description
of high radial order gravity waves is strictly invalid. Similar to the
analogous problem recently studied in fluid dynamics and oceanography (see
e.g. Wunsch, 2017; Baker and Sutherland, 2020, and references therein), the
nonlinear self-interaction of an internal gravity wave propagating in a region
with a spatially-varying Brunt-Väisälä frequency can generate super-harmonic
secondary waves in stars, which can be in near-resonance with the primary mode
and potentially be excited to large amplitudes. If this mechanism operates,
this implies that nonlinear effects could be important for tidal waves even in
stars with convective cores with planetary mass companions (as we have
estimated in § 8).

We have adopted a number of simplifying assumptions to model the generation of
super-harmonic secondary waves. Firstly, we have considered a region of small
spatial extent near $r=r_{c}$, and employed a local Cartesian model instead of
global spherical geometry. Secondly, we assumed the square of the Brunt-
Väisälä frequency to be a linear function of the distance from $r_{c}$ into
the radiative layer (see e.g. Barker, 2011; Ivanov et al., 2013, for a
justification and discussion of the limitations of this assumption). Thirdly,
we have adopted the Boussinesq approximation to the equations of motion, which
formally limits us to a region of small spatial extent near the transition
region.

In this framework, we first considered a single ‘secondary mode’ which has the
smallest frequency distance from the ‘primary mode’ and formulated a condition
for nonlinear effects to be important based on when the amplitude of the near-
resonant secondary is comparable to that of the primary (see eq. (58)). Our
analytical calculations were supported by a numerical study in § 7, where we
analysed the generation of secondary waves as a linear problem with a source
term determined by the presence of the primary wave, which confirmed our
analytical results.

We discussed how our model problem can be related to the normal modes of a
spherical star with a radiative interior and convective envelope in § 6. We
assumed the primary gravity waves were excited in a non-rotating star by a
point-like perturber (i.e. planet or close binary star) on a circular orbit to
illustrate this effect. We formulated the criterion (71), which is based on
(58), and specifies $C_{crit,c}$, which is the smallest value of $q/(1+q)$,
where $q$ is the mass ratio (secondary perturber mass/primary star mass),
needed for the secondary mode with the closest frequency to (double) the tidal
forcing frequency, to exceed the amplitude of the primary tidal wave, which
occurs when $q/(1+q)>C_{crit,c}$. We point out later that several modes with
frequencies sufficiently close to the tidal frequency could be effectively
excited under certain assumptions, and provide the corresponding condition
(75), which can be used if both the secondary modes are non-dissipative and
when the condition of “moderately large dissipation” can be applied to the
secondary modes themselves (approximately equivalent to the statement that the
mode damping time is shorter than the group travel time for a gravity wave
packet across the radiation zone). This condition defines the quantity
$C_{crit}^{dense}$, which is analogous to $C_{crit,c}$, but can be applied
when there is a sufficiently dense spectrum of secondary modes, and typically
gives a slightly more optimistic estimate of nonlinearity. In addition, we
formulate a similar criterion for wave breaking near the centre of a star with
a fully radiative interior (see eq. (22) and the same condition, but for a
star with a convective core near the core itself (see eq. (25) in the
formalism of Ivanov et al. (2013). Both (22) and (25) are analogous to similar
criteria reported by (e.g. Goodman and Dickson, 1998; Ogilvie and Lin, 2007;
Barker and Ogilvie, 2010; Sun et al., 2018; Barker, 2020). These conditions
are formulated in terms of $C_{crit,centre}$ and $C_{crit,core}$,
respectively, and for these nonlinear effects to be important we require
$q/(1+q)>C_{crit,centre}$ or $q/(1+q)>C_{crit,core}$, respectively.

We applied our results to stellar models to estimate the importance of the
nonlinear excitation of super-harmonic secondary waves in § 8. We speculate
that the conditions for transition to non-linear behaviour near radiative-
convective interfaces could result in an efficient damping of primary tidal
waves, thus justifying the assumption of “moderately large dissipation”. We
calculated $C_{crit,centre}$, $C_{crit,core}$, $C_{crit,c}$ and
$C_{crit}^{dense}$ for a number of models of main-sequence stars with masses
$M=1$, $1.4$ and $2M_{\odot}$ and different ages. We found the condition for
wave breaking near the centre is always more important in our solar-mass model
(and others with radiative cores) than the condition for nonlinear behaviour
near $r_{c}$, consistent with prior work. The opposite situation occurs
however in more massive stars. Apart from two models with extended convective
envelopes, and, for one with a radiative centre, for all other models, ages,
and orbital periods, $C_{crit,c}$ and the average values of
$C_{crit,c}^{dense}$ are smaller, and, in many cases, much smaller than
$C_{crit,core}$. In Fig. 12 we show, for example, the behaviour of these
quantities for $M=2M_{\odot}$ and a fixed orbital period
$P_{orb}=1\mathrm{d}$, with stellar age. This figure shows that for ages of
order $10^{9}\mathrm{yr}$, $C_{crit,c}$ can be as small as $10^{-3}$ (i.e.
relevant for planetary-mass companions). If our assumption that the transition
to nonlinear behaviour results in efficient dissipation is valid, this implies
that the nonlinear generation of super-harmonics by tidally-excited gravity
waves could be important for hot Jupiters and other massive companions
orbiting these stars.

In this paper we have ignored stellar rotation (with angular velocity
$\Omega_{r}$) except in Appendix A, under the assumption that the star is
rotating slowly relative to the tidal frequency $\omega\approx 2\Omega_{orb}$,
i.e. that $\Omega_{r}^{2}\ll\omega^{2}$, and certainly that
$\Omega_{r}^{2}\ll\Omega_{*}^{2}$. The neglect of rotation is likely to be
valid for predicting the current and future evolution of most, but not all, of
the shortest-period hot Jupiter systems observed. This is because their host
stars typically rotate with periods much longer than their planetary orbits
(e.g. see the table in Appendix C of Barker, 2020), so we can probably neglect
the corresponding frequency shifts due to rotation (as long as their radiation
zones rotate similarly to their surfaces). For example, WASP-12 is inferred to
have a rotation period longer than 23 days and perhaps as long as 38 days
(based on the observed $V\sin i$, see e.g. Patra et al., 2020, and references
therein), whereas the planet WASP-12 b orbits in only 1.09 days. A brief
calculation of the possible rotational correction to our results in Appendix A
confirms that it doesn’t appear to be significant for such orbital and
rotational periods. In addition, inertial waves (restored by Coriolis forces)
will not be excited by planetary tidal forcing in these stars (which would
require $\omega^{2}\leq 4\Omega_{r}^{2}$, in which case we could no longer
describe the modes as a single spherical harmonic in the form of equation 4).
However, the effects of rotation on the mechanism we have analysed should be
explored in more detail in future work to enable us to study more rapidly
rotating planetary hosts, such as young stars or those that are nearly
rotating synchronously with their planetary orbits e.g. $\tau$-Boo or
WASP-128.

Our initial study here should, however, be considered as a preliminary one.
Firstly, we should confirm our speculation that efficient damping of the
primary tidal wave occurs by exciting super-harmonic secondary waves using
direct nonlinear numerical calculations. Secondly, we should extend our
results to spherical geometry, incorporate stellar rotation, as well as more
realistic Brunt-Väisälä frequency profiles near $r_{c}$, and the possibility
of overshooting in the transition region, etc. These issues are left for
future work. Finally, a different mechanism to efficiently damp gravity waves
in stars with convective cores could involve the excitation of a primary
gravity mode to large amplitudes via resonance locking (e.g. Witte and
Savonije, 1999, 2002; Zanazzi and Wu, 2021; Ma and Fuller, 2021). This
mechanism is unlikely to operate effectively in stars with radiative cores due
to the likelihood of wave breaking, as predicted by (22) in this paper, for
example, but it could potentially be important in stars with convective cores
and should be explored further in those stars. The interaction with
nonlinearly generated super-harmonics should also be considered in that
problem.

## Acknowledgements

We would like to thank the referee Michael Efroimsky for a very prompt and
useful report. We are also grateful to Y. Lazovik for important remarks. PBI
and SVCh were supported in part by the grant 075-15-2020-780
(N13.1902.21.0039) ’Theoretical and experimental studies of the formation and
evolution of extrasolar planetary systems and characteristics of exoplanets’
of the Ministry of Science and Higher Education of the Russian Federation. AJB
was supported by STFC grants ST/R00059X/1 and ST/S000275/1.

## Data availability

There are no new data associated with this article.

##  Appendix A The correction to $\nu$ produced by slow rotation

Rotation also causes a difference between the secondary wave frequency,
$\omega_{s}$, and twice the primary wave frequency, $2\omega_{p}$. We estimate
it below assuming that the rotation axis of the star is perpendicular to the
orbital plane888While hot Jupiters may have significant inclinations of their
rotational axes with respect to orbital plane, this is less common for the
shortest-period planets (e.g. Albrecht et al., 2012). and the stellar
rotational frequency, $\Omega_{r}$ (assumed constant), is much smaller than
the frequencies of interest. We neglect effects that are quadratic in
$\Omega_{r}$, and consider $\omega_{s}$ and $\omega_{p}$ as defined in the
inertial frame, which is different from our definitions in the main text,
where it is implied that these frequencies are defined in the rotating frame.
As discussed in Section 2, when the frequencies are considered in the inertial
frame the primary frequency is simply double the orbital frequency,
$\omega_{p}=2\Omega_{orb}$. For a high-order gravity mode with a frequency
$\omega$ the leading order rotational frequency correction is determined by
(see e.g. Christensen-Dalsgaard, 1998, Equations 8.45 and 8.46)

$\omega\approx\omega_{0}\left(1+\left(1-{1\over
L^{2}}\right){m\Omega_{r}\over\omega_{0}}\right),$ (80)

where $\omega_{0}$ is the mode eigenfrequency for a non-rotating star, and $m$
is the azimuthal wavenumber. We remind the reader that $L=\sqrt{l(l+1)}$, and
we use $l=m=2$ for the primary mode, and $l=m=4$ for the secondary mode. This
choice originates from a comparison of our planar problem with the full
spherical problem. Similar to what is done in Section 6.1, we calculate the
difference $2\omega_{p}-\omega_{s}$ using $\omega_{0}=\omega_{p}$ and
$\omega_{0}=2\omega_{p}$ for the primary and secondary mode, respectively, to
obtain

$2\omega_{p}-\omega_{s}\approx-{7\over 15}\Omega_{r}.$ (81)

Comparing (81) with (61) we see that slow stellar rotation can be accounted
for in the expression for $\nu$ by redefining $f(n)$ entering (63) according
to

$f(n)\rightarrow f(n)-{7\over 60}{\Omega_{r}\over\Omega_{orb}}.$ (82)

It is seen from (82) that the additional term is rather small, approximately
$3.8\times 10^{-3}$, for orbital periods of approximately $1$d and rotational
periods of order $30$d, which are appropriate values for many short period hot
Jupiters like WASP-12. Comparing this value with typical values of $f(n)$
shown in Fig. 1, we see that this term doesn’t appear to be significant for
such systems. However, it may be important for faster rotators. Finally, note
that the correction is negative for prograde rotation, and therefore, the
presence of the additional term in $\nu$ could, in principal, make it
negative. In this case its absolute value should be used in our criteria for
predicting non-linear behaviour.

## References

  * Albrecht et al. (2012) Albrecht, S., Winn, J. N., Johnson, J. A., Howard, A. W., Marcy, G. W., Butler, R. P., Arriagada, P., Crane, J. D., Shectman, S. A., Thompson, I. B., Hirano, T., Bakos, G., Hartman, J. D., 2012, ApJ, 757, 18
  * Baker and Sutherland (2020) Baker, L. E., Sutherland, B. R., 2020, J. Fluid Mech., 891, R1
  * Barker (2011) Barker A. J., 2011, MNRAS, 414, 1365
  * Barker (2020) Barker A. J., 2020, MNRAS, 498, 2270
  * Barker and Astoul (2021) Barker A. J., Astoul A. A. V., 2021, MNRAS Letters, 506, L69
  * Barker and Ogilvie (2009) Barker A. J., Ogilvie G. I., 2009, MNRAS, 395, 2268
  * Barker and Ogilvie (2010) Barker A. J., Ogilvie G. I., 2010, MNRAS, 404, 1849
  * Barker and Ogilvie (2011) Barker A. J., Ogilvie G. .I., 2011, MNRAS, 417, 745
  * Boyd (2001) Boyd, J. P., 2001, Chebyshev and Fourier Spectral Methods, 2nd edn., Dover books on mathematics (Mineola, NY: Dover Publications), ISBN 0486411834.
  * Chernov et al. (2013) Chernov, S. V., Papaloizou, J. C. B., Ivanov, P. B., 2013, MNRAS, 434, 1079
  * Chernov et al. (2017) Chernov, S. V., Papaloizou, J. C. B., Ivanov, P. B., 2017, MNRAS, 470, 2054
  * Cowling (1941) Cowling T. G., 1941, MNRAS, Vol. 101, 367
  * Christensen-Dalsgaard (1998) Christensen-Dalsgaard, J., 1998, Lecture Notes on Stellar Oscillations, 4th Edition, users-phys.au.dk/jcd/oscilnotes/
  * Darwin (1880) Darwin G. H., 1880, Philosophical Transactions of the Royal Society of London, 171, 713
  * Duguid, Barker and Jones (2020a) Duguid C. D., Barker A. J., Jones, C. A., 2020, MNRAS, 491, 923
  * Duguid, Barker and Jones (2020b) Duguid C. D., Barker A. J., Jones, C. A., 2020, MNRAS, 497, 3400
  * Essick and Weinberg (2016) Essick R., Weinberg N., 2016, ApJ, 816, 18
  * Goldreich and Nicholson (1977) Goldreich, P., Nicholson, P. D., 1977, Icarus, 30, 301
  * Goodman and Dickson (1998) Goodman J., Dickson E., 1998, ApJ, 507, 938
  * Goodman and Oh (1997) Goodman J., Oh. S, 1997, ApJ, 486, 403
  * Hut (1981) Hut P., 1981, A&A, 99, 126
  * Ivanov et al. (2013) Ivanov, P. B., Papaloizou, J. C. B., Chernov, S. V., 2013, MNRAS, 432, 2339 (IPCh)
  * Ma and Fuller (2021) Ma L., Fuller J., 2021, ApJ, 918, 16
  * Maciejewski et al. (2016) Maciejewski G. et al., 2016, A&A, 588, L6
  * Mathis (2019) Mathis S., 2019, EAS Publications Series, 82, 5
  * Mustill and Villaver (2012) Mustill A. J., Villaver E., 2012, ApJ, 761, 121
  * Ogilvie and Lin (2007) Ogilvie G. I., Lin D. N. C., 2007, ApJ, 661, 1180
  * Ogilvie and Lesur (2012) Ogilvie G. I., Lesur G., 2012, MNRAS, 422, 1975
  * Ogilvie (2014) Ogilvie G. I., 2014, ARA& A, 52, 171
  * Patra et al. (2020) Patra K. C. et al., 2020, AJ, 159, 150
  * Paxton et al. (2011) Paxton B., et al., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Spiegel and Veronis (1960) Spiegel E. A., Veronis G., 1960, ApJ, 131, 442
  * Sun et al. (2018) Sun M., Arras P., Weinberg N. N., Troup N. W., Majewski S. R., MNRAS, 481, 4077
  * Terquem et al. (1998) Terquem C., Papaloizou J., Nelson R., Lin D., 1998, ApJ, 502, 788
  * Terquem (2021) Terquem C., 2021, MNRAS, 503, 5789
  * Verbunt and Phinney (1995) Verbunt F., Phinney E. S., 1995, A&A, 296, 709
  * Vidal and Barker (2020a) Vidal J., Barker A. J., 2020, ApJL, 888(2)
  * Vidal and Barker (2020b) Vidal J., Barker A. J., 2020, MNRAS, 497, 4472
  * Weinberg et al. (2012) Weinberg N., Arras P., Quateart E., Burkart J., 2012, ApJ, 751, 136
  * Weinberg et al. (2017) Weinberg N., Sun M., Arras P., Essick .R, 2017, ApJL 849 L11
  * Wilkins et al. (2017) Wilkins A. .N., Delrez L., Barker A. .J., et al., 2017, ApJL, 836 L24
  * Witte and Savonije (1999) Witte M. G.,Savonije G. J., 1999, A&A, 350, 129
  * Witte and Savonije (2002) Witte M. G.,Savonije G. J., 2002, A&A, 386, 222
  * Wunsch (2017) Wunsch, S., 2017, J. Fluid Mech., 828, 630
  * Yee et al. (2020) Yee S. W., et al., 2020, ApJ, 888, L5
  * Zahn (1966) Zahn J. P., 1966, Annales d’Astrophysique, 29, 313
  * Zahn (1977) Zahn J. P., 1977, A&A, 57, 383
  * Zahn (1989) Zahn J. P., 1989, A&A, 220, 112
  * Zanazzi and Wu (2021) Zanazzi J. J., Wu Y., 2021, AJ, 161, 263