Horseshoes and spiral waves: capturing the 3D flow induced by a low-mass planet analytically

ABSTRACT

The key difficulty faced by 2D models for planet–disc interaction is in appropriately accounting for the impact of the disc’s vertical structure on the dynamics. 3D effects are often mimicked via softening of the planet’s potential; however, the planet-induced flow and torques often depend strongly on the choice of softening length. We show that for a linear adiabatic flow perturbing a vertically isothermal disc, there is a particular vertical average of the 3D equations of motion that exactly reproduces 2D fluid equations for arbitrary adiabatic index. There is a strong connection here with the Lubow–Pringle 2D mode of the disc. Correspondingly, we find a simple, general prescription for the consistent treatment of planetary potentials embedded within ‘2D’ discs. The flow induced by a low-mass planet involves large-scale excited spiral density waves that transport angular momentum radially away from the planet and ‘horseshoe streamlines’ within the coorbital region. We derive simple linear equations governing the flow that locally capture both effects faithfully simultaneously. We present an accurate coorbital flow solution allowing for inexpensive future study of corotation torques, and predict the vertical structure of the coorbital flow and horseshoe region width for different values of adiabatic index, as well as the vertical dependence of the initial shock location. We find strong agreement with the flow computed in 3D numerical simulations, and with 3D one-sided Lindblad torque estimates, which are a factor of 2–3 lower than values from previous 2D simulations.

accretion, accretion discs, hydrodynamics, waves, methods: analytical, planet–disc interactions, protoplanetary discs

1 INTRODUCTION

Discs of gas in orbital motion around a massive central body are found in many astronomical contexts. Objects embedded within these discs, e.g. planets embedded in circumstellar discs or stellar-mass black holes within active galactic nucleus (AGN) discs, are in general subject to strong gravitational interactions with the gas and dust that comprise the discs. Indeed, signatures such as gaps, rings, and spirals observed in discs’ emission are highly suggestive of large embedded planets (Grady et al. 2013; ALMA Partnership 2015; Benisty et al. 2015; Andrews et al. 2018; Andrews 2020). These interactions lead to angular momentum exchange between the object and the disc, causing the object to migrate.

Type I migration concerns low-mass embedded objects that excite predominantly linear perturbations to the flow in the disc. Not only is it a crucial ingredient in the theory of planet formation and population synthesis (Mordasini 2018), but it may also play a key role in understanding the rapid merger rate of stellar-mass black hole (sBH) binaries inferred from detections by the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) (Secunda et al. 2019). While the migration of sBHs within very thin AGN discs is arguably better suited to the type I parameter regime, the development of type I migration theory has been driven by the need to understand and predict the dynamics of young planets.

Linear 2D analyses from early work (Goldreich & Tremaine 1979, 1980; Artymowicz 1993; Korycansky & Pollack 1993) presented a picture in which the planet exchanged angular momentum with the disc via two mechanisms: some angular momentum was transported away from the planet by density waves excited at Lindblad resonances, and some separately ‘absorbed’ into corotation resonances. The torque depends only on the properties of the disc and their gradients near the planet, since it is exerted predominantly locally. Type I torque formulae therefore take very simple forms. Tanaka, Takeuchi & Ward (2002), and more recently Tanaka & Okada (2024), extended the 2D linear calculation to 3D locally isothermal discs, arriving at a torque formula in good agreement with the numerical simulations of D’Angelo & Lubow (2010), which drives rapid inward migration on a time-scale of |$\sim 10^5$| yr for an Earth-mass planet.

This unopposed rapid inward migration is worrying; it acts as a potential obstacle to the feasibility of the core accretion model, in which giant gas planets form via the accretion of gas onto a solid core on a time-scale of |$\gtrsim$|10⁶ yr (Lissauer et al. 2009). Furthermore, since the lifetime of the disc is |$\lesssim$|10⁷ yr, it suggests at surface level an additional paradoxical existential threat for extra-solar Earth-like planets. Indeed, planetary population synthesis models have struggled to reproduce the observed distribution of semimajor axes among detected low-mass exoplanets from the linear type I torque estimates (Mordasini et al. 2009).

However, as pointed out by Paardekooper & Papaloizou (2009a, b), the linear treatment of the corotation resonance is inappropriate, unless the disc is sufficiently viscous. Close to corotation, the radial displacement of fluid elements becomes non-linear, as they execute ‘horseshoe turns’ [akin to the behaviour of test particles seen in Dermott & Murray (1981) and Petit & Henon (1986)], exchanging angular momentum with the planet as they are excited onto inner or outer orbits. An asymmetry between the potential vorticity (PV), or vortensity, on incident inner and outer orbits leads to a net torque comparable to that of the Lindblad torque (Ward 1991), not captured in the linear theory.

In addition, the longer time-scale associated with this horseshoe motion means that additional physical effects become relevant, each with an associated strong torque with potential to resolve the aforementioned ‘paradoxes’. This has motivated a wealth of research, both semi-analytical and computational, to quantify the dependence of the corotation torque on diffusive, radiative, and migration-feedback processes, among others (Balmforth & Korycansky 2001; Masset 2001, 2017; Masset & Papaloizou 2003; Paardekooper & Mellema 2006; Paardekooper & Papaloizou 2008, 2009a, b; Masset & Casoli 2009, 2010; Paardekooper et al. 2010, 2023; Paardekooper, Baruteau & Kley 2011; Lega et al. 2014; Fung, Artymowicz & Wu 2015; Masset & Benítez-Llambay 2016; McNally et al. 2020).

The torques that arise often depend strongly on the precise geometry of the coorbital flow, and correspondingly the value of the softening length, b (often used to modify the planetary potential to mimic 3D effects in 2D disc models). The variation of torque components can be more than a factor of 3 for different reasonable choices of b (e.g. Tanaka & Okada 2024, fig. 5); such varied behaviour for different softening lengths is an issue not limited to the study of planetary torques (Müller, Kley & Meru 2012). In this paper, we aim to accurately describe the flow induced by the planet, on top of which torque-inducing physics may be studied. Though there is much to be achieved and understood here, the coorbital flow has received little analytical attention since Paardekooper & Papaloizou (2009b), despite their analysis suffering some limitations (which we point out in Section 3).

A key objective of this work is to address the problem of capturing 3D disc physics with 2D equations. In particular, for low-mass objects embedded within vertically isothermal discs that permit adiabatic perturbations, we show that there exists a particular vertical average of the equations of motion that reproduces the commonly adopted 2D fluid equations. This averaging procedure is closely related to the operator that projects onto the Lubow–Pringle 2D mode (Lubow & Pringle 1993). This paper is concerned with the non-axisymmetric generalization of this mode. The 2D mode is so called as it has the property |$v_z = 0$|⁠, but in general it is z-dependent. It is a member of the wider family of 3D disc modes and describes the spiral wake as well as the horseshoe motion within the coorbital region. Most importantly, our averaging process yields a simple, general prescription for the consistent treatment of planetary potentials embedded within ‘2D’ discs, as well as an interpretation for 2D models.

A quantitative description and understanding of the 2D mode of the planet-induced flow, which captures the impact of the disc’s vertical extent on the spiral wake, also has important observational applications. Detection of disc-embedded exoplanets via their kinematic signatures is a promising method for finding very young protoplanets during their formation stages (Pinte et al. 2023). As Fasano et al. (2024) point out, 3D effects impact the precise structure of the planet’s spiral wake, with important consequences for the observational analysis of this signature, and consequently planet mass estimates.

In Section 2, we project the 3D equations of motion onto this 2D flow. We derive (without further approximation) independent, second-order, linear equations governing the 2D flow in Section 3. In Section 4, we present the flow solution. We discuss our findings in Section 5, and draw our conclusions in Section 6.

2 GOVERNING EQUATIONS

For the remainder of the paper, we will assume two quantities to be small. First, we assume the disc’s aspect ratio |$h = H/r$| to be much smaller than 1. Secondly, we assume the planet’s mass |$M_{\mathrm p}$| to be much smaller than the ‘thermal mass’ (Goodman & Rafikov 2001). That is,

$$\begin{eqnarray} q \equiv \frac{M_{\mathrm p}}{M_\star } \ll \frac{M_\text{th}}{M_\star } = h^3, \end{eqnarray}$$

(1)

where |$M_\star$| is the mass of the central star. This ensures that the equations governing the flow are almost everywhere well approximated as linear (this may be verified post hoc), and that the disc’s vertical structure is almost everywhere determined by the star’s gravity. We therefore have two small parameters, namely h and |$\frac{q}{h^3}$|⁠; we will exploit the fact that both are small to make analytical progress.

2.1 Unperturbed state

Before continuing, it is important to define precisely the model for the disc with which our planet will interact. We will consider the planet to introduce a perturbation to the background disc structure outlined below.

Our background disc is steady, axisymmetric, and comprised entirely of an ideal gas. It is in orbit about a star of mass |$M_\star$|⁠, fixed at the centre of our frame, and is inviscid and non-self-gravitating to a first approximation. The gas that comprises it solves the steady Euler equations and ideal gas equation,

$$\begin{eqnarray} {{\boldsymbol u}}_0\cdot \nabla {{\boldsymbol u}}_0 = - \frac{1}{\rho _0} \nabla p_0-\nabla \Phi _0, \end{eqnarray}$$

(2a)

$$\begin{eqnarray} \nabla \cdot \left(\rho _0 {{\boldsymbol u}}_0\right) = 0, \end{eqnarray}$$

(2b)

$$\begin{eqnarray} \frac{p_0}{\rho _0} = \frac{k_\mathrm{ B}}{\bar{\mu }}T_0, \end{eqnarray}$$

(2c)

where |${{\boldsymbol u}}_0$|⁠, |$\rho _0$|⁠, |$p_0$|⁠, and |$T_0$| are the velocity, density, pressure, and temperature of the background disc. In addition, |$k_\mathrm{ B}$| is Boltzmann’s constant, and |$\bar{\mu }$| the mean molecular mass. We introduce the cylindrical coordinates |$(r, \theta , z)$|⁠, with the z-axis normal to the plane of the disc, so that

$$\begin{eqnarray} \Phi _0 = -\frac{G M_\star }{\sqrt{r^2+z^2}} \end{eqnarray}$$

(3)

is the gravitational potential of the central star, and the velocity of the background state may be written as |${{\boldsymbol u}}_0 = r\Omega (r,z) {\boldsymbol e}_\theta$|⁠.

We assume the background disc’s temperature, |$T_0$|⁠, to be a prescribed function of r only, the result of a relatively fast thermal relaxation in the vertical direction compared to the long time-scale of the evolution of the disc, and ignoring the hotter irradiated outer layers. We remark that the background disc need only be ‘locally isothermal’ to a first approximation for the analysis performed in this paper to hold. Though this background state is idealized, the calculation we perform is robust: we define perturbations relative to the exact background state (though we are perhaps ignorant of its precise structure), but only need knowledge of its leading order behaviour in order to evaluate these perturbations.

Furthermore, the isothermal prescription applies only to the background state: we will consider the planet to perturb this background state adiabatically. This represents a good approximation provided the time-scale for thermal relaxation is longer than the planet’s orbital period. We discuss the validity of this assumption in more detail in Section 5.3.3. We define the isothermal sound speed, |$c_s(r)$|⁠, via

$$\begin{eqnarray} c_\mathrm{ s}^2 = \frac{k_\mathrm{ B}}{\bar{\mu }}T_0. \end{eqnarray}$$

(4)

We define the scale height, H, and aspect ratio, h via

$$\begin{eqnarray} h \equiv \frac{H}{r} \equiv \frac{c_\mathrm{ s}}{r\Omega _\mathrm{ K}}, \end{eqnarray}$$

(5)

where |$\Omega _\mathrm{ K}$| is the Keplerian angular frequency, satisfying

$$\begin{eqnarray} \frac{G M_\star }{r^3} = \Omega _\mathrm{ K}^2. \end{eqnarray}$$

(6)

We assume |$h \ll 1$|⁠, so that the disc is thin. The angular velocity of the gas in the unperturbed disc satisfies

$$\begin{eqnarray} \Omega = \Omega _\mathrm{ K}\left(1 + \mathcal {O}(h^2)\right). \end{eqnarray}$$

(7)

The density and pressure then satisfy

$$\begin{eqnarray} \rho _0 = \frac{p_0}{c_\mathrm{ s}^2} = \frac{\Sigma _0(r)}{\sqrt{2\pi } H}\exp {\left(-\frac{z^2}{2H^2}\right)}, \end{eqnarray}$$

(8)

where |$\Sigma _0(r)$| is the surface density.

2.2 Perturbation equations

We now introduce a planet of mass |$M_{\mathrm p}$| on a Keplerian circular orbit of radius |$r_\mathrm{ p}$| and angular frequency |$\Omega _\mathrm{ p} = \Omega _K(r_\mathrm{ p})\sqrt{1+q}$| to the background disc. We consider the perturbation problem in the frame corotating with the planet, and assume further that the flow in this frame is steady. The fluid velocity in this frame is |${{\boldsymbol v}} = {{\boldsymbol u}}-r\Omega _\mathrm{ p} {{\boldsymbol e}}_\theta$|⁠, and the relevant Euler equations become

$$\begin{eqnarray} {{\boldsymbol v}}\cdot \nabla {{\boldsymbol v}} + 2{{\boldsymbol \Omega }}_\mathrm{ p} \times {{\boldsymbol v}} = - \frac{1}{\rho } \nabla p-\nabla \Phi _t-\nabla \Psi _\mathrm{ p}, \end{eqnarray}$$

(9a)

$$\begin{eqnarray} \nabla \cdot \left(\rho {{\boldsymbol v}}\right) = 0, \end{eqnarray}$$

(9b)

$$\begin{eqnarray} {{\boldsymbol v}} \cdot \nabla \left(p \rho ^{-\gamma }\right) = 0, \end{eqnarray}$$

(9c)

where |$\Phi _t = \Phi _0(r,z)-\frac{1}{2}(r^2- 3r_\mathrm{ p}^2)\Omega _\mathrm{ p}^2$| is the tidal potential, and the planet’s potential, including the indirect term arising from the acceleration of the frame centred on the star, is given by

$$\begin{eqnarray} \Psi _\mathrm{ p} = - \frac{G M_\mathrm{ p}}{\sqrt{r_\mathrm{ p}^2 + r^2-2 r r_\mathrm{ p} \cos {\theta } + z^2}} + q r_\mathrm{ p} \Omega _\mathrm{ p}^2 r \cos {\theta }. \end{eqnarray}$$

(10)

We now define local quasi-Cartesian coordinates |$x = r-r_\mathrm{ p}$|⁠, |$y = r_\mathrm{ p}\theta$|⁠, and perturbed variables

$$\begin{eqnarray} \rho ^{\prime } = \rho -\rho _0, \end{eqnarray}$$

(11a)

$$\begin{eqnarray} p^{\prime } = p-p_0, \end{eqnarray}$$

(11b)

$$\begin{eqnarray} &v^{\prime }_x = v_r, \end{eqnarray}$$

(11c)

$$\begin{eqnarray} v^{\prime }_y = v_\theta -r(\Omega - \Omega _{\mathrm p}) = u_\theta -r \Omega , \end{eqnarray}$$

(11d)

$$\begin{eqnarray} v^{\prime }_z = v_z. \end{eqnarray}$$

(11e)

It is useful to define further the time-derivative following the orbital shear flow

$$\begin{eqnarray} D = -\frac{3}{2}x\Omega _\mathrm{ p}\partial _y, \end{eqnarray}$$

(12)

and upon subtraction of the background state solution from the Euler equations (9), we find the balance at the next order in h (assuming |$x = \mathcal {O}\left(H\right)$| and |$y = \mathcal {O}\left(H\right)$|⁠) is given by

$$\begin{eqnarray} D v^{\prime }_x-2\Omega _{\mathrm p} v^{\prime }_y + \frac{\partial _x p^{\prime }}{\rho _{\mathrm p}} = -\partial _x \phi _{\mathrm p}, \end{eqnarray}$$

(13a)

$$\begin{eqnarray} D v^{\prime }_y + \frac{1}{2}\Omega _{\mathrm p} v^{\prime }_x + \frac{\partial _y p^{\prime }}{\rho _{\mathrm p}} = -\partial _y \phi _{\mathrm p}, \end{eqnarray}$$

(13b)

$$\begin{eqnarray} D v^{\prime }_z + \Omega _{\mathrm p}^2 z \frac{\rho ^{\prime }}{\rho _{\mathrm p}} + \frac{\partial _z p^{\prime }}{\rho _{\mathrm p}} = -\partial _z \phi _{\mathrm p}, \end{eqnarray}$$

(13c)

$$\begin{eqnarray} D \rho ^{\prime } + \rho _{\mathrm p}\partial _xv^{\prime }_x + \rho _{\mathrm p}\partial _yv^{\prime }_y + \partial _z\left(\rho _{\mathrm p}v^{\prime }_z\right) = 0, \end{eqnarray}$$

(13d)

$$\begin{eqnarray} D\left(p^{\prime }-\gamma c_{\mathrm p}^2\rho ^{\prime }\right) + (\gamma -1)\Omega _{\mathrm p}^2 z \rho _{\mathrm p} v^{\prime }_z = 0, \end{eqnarray}$$

(13e)

where |$c_{\mathrm p} = c_s(r_{\mathrm p})$|⁠,

$$\begin{eqnarray} \phi _{\mathrm p} = -\frac{G M_{\mathrm p}}{\sqrt{x^2+y^2+z^2}}, \end{eqnarray}$$

(14)

and we have taken the leading order approximation to the background disc’s density,

$$\begin{eqnarray} \rho _{p}(z) \equiv \frac{\Sigma _{\mathrm p}}{\sqrt{2\pi } H_{\mathrm p}}\exp {\left(-\frac{z^2}{2H_{\mathrm p}^2}\right)}, \end{eqnarray}$$

(15)

for |$\Sigma _{\mathrm p} = \Sigma _0(r_{\mathrm p})$| and |$H_{\mathrm p} = H(r_{\mathrm p})$|⁠. The system of equations (13) describes locally the flow excited by a planet embedded in a stratified 3D disc. We could in theory relax the assumption |$y \lesssim H$| by reintroducing the azimuthally global expression for the planet’s potential, |$\Psi _{\mathrm p}$|⁠, in place of |$\phi _{\mathrm p}$|⁠, and applying periodic boundary conditions at |$y = \pm \pi r_{\mathrm p}$|⁠. The corresponding correction, however, is only of relative size |$\mathcal {O}\left(h\right)$|⁠, so that local system involving |$\phi _{\mathrm p}$| (which we adopt for the remainder of the paper) becomes exact in the limit |$h \rightarrow 0$|⁠.

The system (equations 13) governs the excitation of the spiral density waves by the planet, permits downstream gravity waves (e.g. Lubow & Zhu 2014), and also describes the horseshoe trajectories followed by fluid elements close to the planet’s orbital radius. As noted by several authors, e.g. by Tanaka et al. (2002), the density wave excitation is confined to the region extending only a few scale heights from the planet, though equations (13) do not capture the asymmetry between the inner and outer spiral wakes. Importantly, as we will demonstrate in the next section, equations (13) also have the elegant property that it may be vertically integrated to derive exact 2D flow equations.

2.3 Projection onto a 2D flow

Remarkably, under a particular choice of vertical averaging, the system of equations (13) may be transformed exactly into the familiar linearized 2D flow equations. Associated with this averaging procedure is a definite choice for the ‘softening’ of the planetary potential (specified in equation 22).

The averaging procedure defined below is closely related to the operator which projects onto the 2D mode of the disc found by Lubow & Pringle (1993). Indeed, the 2D equations we find govern the radial and azimuthal evolution of the amplitude of this mode. The 2D mode is a member of a larger family of 3D modes, it exists for |$\gamma < 2$|⁠, and is often referred to as ‘2D’, since it has the property that |$v^{\prime }_z = 0$|⁠, despite possessing a vertical dependence.

Before we proceed, it is helpful to first introduce the adiabatic sound speed, scale height, and aspect ratio as

$$\begin{eqnarray} c_\gamma = \sqrt{\gamma } c_{\mathrm p}, \quad H_\gamma = c_\gamma /\Omega _{\mathrm p}, \quad h_\gamma = H_\gamma /r_{\mathrm p}. \end{eqnarray}$$

(16)

We combine (13d) and (13e) to eliminate |$\rho ^{\prime }$|⁠:

$$\begin{eqnarray} c_\gamma ^2 &\times &(13d) + (13e) \\ &\Rightarrow & D p^{\prime } + c_\gamma ^2 \rho _{\mathrm p}\left(\partial _x v^{\prime }_x + \partial _y v^{\prime }_y\right) \\ &&+\, (\gamma -1)\Omega _{\mathrm p}^2 z \rho _{\mathrm p} v^{\prime }_z + c_\gamma ^2\partial _z\left(\rho _{\mathrm p}v^{\prime }_z\right) = 0,\\ &\Rightarrow & D \frac{p^{\prime }}{\rho _{\mathrm p}} + c_\gamma ^2 \left(\partial _x v^{\prime }_x + \partial _y v^{\prime }_y\right) \\ &&+\,c_\gamma ^2 \exp {\left(\tfrac{z^2}{2 H_\gamma ^2}\right)}\partial _z\left[v^{\prime }_z \exp {\left(-\tfrac{z^2}{2 H_\gamma ^2}\right)}\right] = 0. \end{eqnarray}$$

(17)

We now multiply equation (17) by the factor |$\exp {(-\tfrac{z^2}{2 H_\gamma ^2})} \propto \rho _{\mathrm p}^{1/\gamma }$| and integrate. (Incidentally, the function |$\rho _{\mathrm p}^{1/\gamma }$| is the product of the background density distribution and the vertical profile of the 2D mode. In this way, the vertical integration may formally be seen to be an inner product with the 2D mode.) The result of the integration is the 2D mass conservation analogue

$$\begin{eqnarray} D P^{\prime } + \Sigma _{\mathrm p} c_\gamma ^2\left(\partial _x \bar{v}_x + \partial _y \bar{v}_y\right) = 0, \end{eqnarray}$$

(18)

where we have defined

$$\begin{eqnarray} \bar{v}_x \equiv \left\langle v^{\prime }_x\right\rangle , \quad \bar{v}_y \equiv \left\langle v^{\prime }_y\right\rangle , \quad P^{\prime } \equiv \Sigma _{\mathrm p} \left\langle p^{\prime }/\rho _{\mathrm p}(z)\right\rangle , \end{eqnarray}$$

(19)

for vertical average |$\left\langle \, \cdots \right\rangle$| defined via

$$\begin{eqnarray} \left\langle X \right\rangle \equiv \frac{1}{\sqrt{2\pi }H_\gamma }\int _{-\infty }^{\infty }X\exp {\left(-\frac{z^2}{2H_\gamma ^2}\right)}\mathrm{d}z. \end{eqnarray}$$

(20)

We may similarly apply the averaging procedure |$\left\langle \, \cdots \right\rangle$| defined in equation (20) to equations (13a) and (13b) to obtain the 2D momentum equations

$$\begin{eqnarray} D\bar{v}_x-2\Omega _{\mathrm p} \bar{v}_y + \frac{1}{\Sigma _{\mathrm p}}\partial _x P^{\prime } = -\partial _x \Phi _{\mathrm p}, \end{eqnarray}$$

(21a)

$$\begin{eqnarray} &D\bar{v}_y + \frac{1}{2}\Omega _{\mathrm p} \bar{v}_x + \frac{1}{\Sigma _{\mathrm p}}\partial _y P^{\prime } = -\partial _y \Phi _{\mathrm p}, \end{eqnarray}$$

(21b)

where now |$\Phi _{\mathrm p}$| is precisely defined as |$\left\langle \phi _{\mathrm p} \right\rangle$|⁠, that is,

$$\begin{eqnarray} \Phi _{\mathrm p} \equiv \left\langle -\frac{G M_{\mathrm p}}{\sqrt{x^2 + y^2 + z^2}} \right\rangle = -\frac{G M_{\mathrm p}}{H_\gamma }\frac{\mathrm{e}^{\frac{1}{4}s^2}}{\sqrt{2\pi }}K_0\left(\tfrac{1}{4}s^2\right), \end{eqnarray}$$

(22)

with |$s = \sqrt{x^2 + y^2}/H_\gamma$|⁠, and |$K_0(z)$| the modified Bessel function. Note that in a global 2D disc model, the above expression remains valid upon redefining |$s = \left|{{\boldsymbol r}}-{{\boldsymbol r}}_{\mathrm p}\right|/H_\gamma$| and reintroducing the indirect term (this matter is discussed in greater detail in Section 5.1). This prescription for the potential is therefore generally applicable to 2D disc models. Reassuringly, at large distances, the above expression becomes

$$\begin{eqnarray} \Phi _{\mathrm p} = -\frac{G M_{\mathrm p}}{\sqrt{x^2 + y^2}}\left(1 + \mathcal {O}\left(\frac{H_\gamma ^2}{x^2+y^2}\right)\right), \end{eqnarray}$$

(23)

and close to the planet, our 2D potential has a logarithmic singularity. We therefore have exact 2D equations forced by a 2D potential, valid for arbitrary adiabatic index |$\gamma$|⁠.

We remark that we did not make use of equation (13c) for the vertical velocity; the 2D mode is orthogonal to and ignorant of the permitted vertical motions of the disc (including downstream gravity waves and inertial waves). Studies of these excited gravity waves indicate they may have an important observational signature and impact on the torque on the planet, which is difficult to predict due to the intricacy of the problem (Zhu, Stone & Rafikov 2012; Lubow & Zhu 2014; McNally et al. 2020; Bae, Teague & Zhu 2021; Ziampras, Paardekooper & Nelson 2023). That being said, much of the important dynamics are captured by the 2D mode, and we shall focus the majority of our attention here for the remainder of the paper.

2.4 Relation to 2D Euler equations

2D disc models enjoy a valuable simplicity compared to 3D models. However, the non-linear 2D Navier–Stokes equations comprise only an approximate model for the flow in the disc. The analysis performed above (as well as that in Section 2.7) provides an interpretation for the 2D flow variables in terms of their counterparts in a 3D disc with vertical extent. This relationship holds when the disc’s response (e.g. to forcing by a perturbing embedded planet) constitutes a linear perturbation to its background state. The 2D variables are given by weighted vertical averages of the (adiabatic) pressure and velocity perturbations to the locally isothermal background state of a 3D disc. Specifically, if we define

$$\begin{eqnarray} \bar{\Sigma } \equiv \Sigma _{\mathrm p} + \frac{P^{\prime }}{c_\gamma ^2}, \end{eqnarray}$$

(24)

$$\begin{eqnarray} \bar{{\boldsymbol v}} \equiv \bar{v}_x{{\boldsymbol e}}_x + \left(-\tfrac{3}{2}\Omega _\mathrm{ p} x + \bar{v}_y\right){{\boldsymbol e}}_y, \end{eqnarray}$$

(25)

then we see that our vertically averaged flow equations (18), (21a), and (21b) are linearizations of the familiar 2D barotropic model,

$$\begin{eqnarray} {\bar{{\boldsymbol v}}}\cdot \nabla \bar{{{\boldsymbol v}}} + 2{{\boldsymbol \Omega }}_\mathrm{ p} \times \bar{{{\boldsymbol v}}} = - \frac{1}{\bar{\Sigma }} \nabla \bar{P}-\nabla \bar{\Phi }_\mathrm{ t}-\nabla \Phi _\mathrm{ p}, \end{eqnarray}$$

(26a)

$$\begin{eqnarray} \nabla \cdot \left(\bar{\Sigma } \bar{{\boldsymbol v}}\right) = 0, \end{eqnarray}$$

(26b)

$$\begin{eqnarray} \bar{P} = K \bar{\Sigma }^\gamma , \end{eqnarray}$$

(26c)

where |$\bar{\Phi }_t = -\tfrac{3}{2}\Omega _{\mathrm p}^2 x^2$|⁠, and |$K = c_{\mathrm p}^2\bar{\Sigma }_{\mathrm p}^{1-\gamma }$| uniformly, so that |$\mathrm{d}\bar{P} = c_\gamma ^2 \mathrm{d}\bar{\Sigma }$|⁠. We remind the reader of the formal definitions of |$P^{\prime }$|⁠, |$\bar{v}_x$|⁠, and |$\bar{v}_y$|⁠, given in equation (19), and our convention |$c_\gamma ^2 = \gamma c_{\mathrm p}^2$|⁠. [Note we may treat (26c) as an exact definition of |$\bar{P}$|⁠.]

We note that the error in the system (26) scales as the size of the quadratic non-linear terms, that is, |$(q/h_\gamma ^3)^2$|⁠. This is notably far smaller than the error introduced by any alternative averaging process, which would be of the order of |$q/h_\gamma ^3$|⁠. Recall we assumed |$q \ll h_\gamma ^3$| so that our 3D flow was linear, and that the errors introduced in the linearization of the 3D system also scaled as |$(q/h_\gamma ^3)^2$|⁠. Furthermore, while we have addressed here only a local model of disc dynamics, this approach may be generalized to a global disc. It is sufficient for the background disc to satisfy a locally isothermal equation of state.

It is worth mentioning that the ‘surface density’ |$\bar{\Sigma }$| in the 2D model must really be thought of in terms of the vertically averaged pressure perturbation. Moreover, equation (26b) does not technically enforce physical mass conservation; mass conservation would be derived from a direct (unweighted) vertical integral of the 3D mass conservation equation (13d). Instead, (26b) describes the evolution of the pressure due to a combination of compressive and buoyant motions; however, it takes precisely the same form as a 2D mass conservation equation (and for the remainder of this paper this is how we shall think of it).

In contrast, the 2D model’s flow velocities (⁠|$\bar{v}_x$| and |$\bar{v}_y$|⁠) can be derived from simple weighted vertical averages of their 3D counterparts. Finally, we remind the reader that the ‘planetary potential’, |$\Phi _{\mathrm p}$|⁠, which forces this 2D system should be taken to be the vertical average of the 3D potential specified in equation (22). This prescription is generally applicable to 2D disc models, and is discussed in more detail in Section 5.1.

2.5 Stream function and Bernoulli invariant

The stream function and streamlines of our 2D vertically averaged flow are of particular importance near to corotation. This 2D behaviour is expected to contribute dominantly to the torque exerted on the planet, and for low-mass planets, the corotation torque is expected to scale with the fourth power of the horseshoe region width. To find and compute these streamlines, it is instructive to consider the exact 2D flow which solves (26) (which is well approximated by our averaged flow), and exactly conserves the Bernoulli invariant (which we demonstrate below). Equation (26b) implies the existence of a stream function |$\psi$|⁠. We take the definition

$$\begin{eqnarray} \bar{\Sigma }\bar{{\boldsymbol v}} = -{{\boldsymbol e}}_z \times \nabla \psi . \end{eqnarray}$$

(27)

We may rewrite (26a) as

$$\begin{eqnarray} \left(2 {{\boldsymbol \Omega }}_\mathrm{ p} + \nabla \times \bar{{\boldsymbol v}}\right) \times \bar{{\boldsymbol v}} + \nabla \left(\frac{1}{2}|\bar{{\boldsymbol v}}|^2 + \bar{W} + \bar{\Phi }_\mathrm{ t} + \Phi _\mathrm{ p}\right) = {{\bf 0}}, \end{eqnarray}$$

(28)

where |$\mathrm{d}\bar{W} = \mathrm{d}\bar{P}/\bar{\Sigma }$|⁠. We define the Bernoulli invariant, B, as

$$\begin{eqnarray} B = \frac{1}{2}|\bar{{\boldsymbol v}}|^2 + \bar{W} + \bar{\Phi }_\mathrm{ t} + \Phi _\mathrm{ p}, \end{eqnarray}$$

(29)

and introduce the PV (or vortensity)

$$\begin{eqnarray} \zeta = \frac{2 \Omega _{\mathrm p} + {{\boldsymbol e}}_z\cdot (\nabla \times \bar{{\boldsymbol v}})}{\bar{\Sigma }}. \end{eqnarray}$$

(30)

We may then rewrite equation (28) as

$$\begin{eqnarray} && \bar{\Sigma }\zeta {{\boldsymbol e}}_z \times \bar{{\boldsymbol v}} + \nabla B = 0 \\&&\quad \Rightarrow -\zeta {{\boldsymbol e}}_z \times \left({{\boldsymbol e}}_z \times \nabla \psi \right) + \nabla B = 0\\&&\quad \Rightarrow \nabla B = -\zeta \nabla \psi . \end{eqnarray}$$

(31)

That is to say, the Bernoulli function is constant on streamlines, |$B = B(\psi)$|⁠, and further the PV,

$$\begin{eqnarray} -\frac{\mathrm{d}B}{\mathrm{d}\psi } = \zeta (\psi) \end{eqnarray}$$

(32)

is also conserved. We now further impose that |$\zeta$| is uniformly constant upstream (which holds to leading order within the local approximation), so that

$$\begin{eqnarray} \zeta \equiv \frac{\Omega _{\mathrm p}}{2 \Sigma _{\mathrm p}}. \end{eqnarray}$$

(33)

Combining this with the definition (30) implies that linearly

$$\begin{eqnarray} \partial _x \bar{v}_y-\partial _y \bar{v}_x-\frac{\Omega _{\mathrm p}}{2 c_\gamma ^2} \frac{P^{\prime }}{\Sigma _{\mathrm p}} = 0. \end{eqnarray}$$

(34)

In other words, the PV is uniform. Integrating equation (32) (and setting the arbitrary constant of integration to zero) then gives an expression for the stream function for the 2D flow

$$\begin{eqnarray} -\frac{\Omega _{\mathrm p}}{2 \Sigma _{\mathrm p}}\psi = \frac{1}{2}|\bar{{\boldsymbol v}}|^2 + \bar{W} + \bar{\Phi }_t + \Phi _{\mathrm p}. \end{eqnarray}$$

(35)

That is to say, our 2D vertically averaged flow has streamlines which are the contours of the stream function

$$\begin{eqnarray} \frac{\psi }{\Sigma _{\mathrm p}} = \frac{3}{4}\Omega _{\mathrm p} x^2-\frac{2}{\Omega _{\mathrm p}}\left(\frac{P^{\prime }}{\Sigma _{\mathrm p}}+\Phi _{\mathrm p}\right) + 3 x \bar{v}_y + \mathcal {O}\left(\frac{q^2}{h_\gamma ^6}\frac{c_\gamma ^2}{\Omega _{\mathrm p}}\right). \end{eqnarray}$$

(36)

2.6 On the corotation singularity

Importantly, equation (36) allows us to write down the equation for the radial displacement of fluid elements, |$\bar{\xi }_x$|⁠, which we take to satisfy

$$\begin{eqnarray} \bar{{\boldsymbol v}}\cdot \nabla \bar{\xi }_x = \bar{v}_x, \end{eqnarray}$$

(37)

with |$\bar{\xi }_x = 0$| far upstream. If a given streamline (or contour of |$\psi$|⁠) has radial location |$x_0$| far upstream, we have that

$$\begin{eqnarray} \bar{\xi }_x = x-x_0. \end{eqnarray}$$

(38)

Comparing values of |$\psi$| far upstream to the point |$(x,y)$|⁠, we see that at leading order

$$\begin{eqnarray} \psi &=& \frac{3}{4}\Sigma _{\mathrm p}\Omega _{\mathrm p} x_0^2 = \frac{3}{4}\Sigma _{\mathrm p}\Omega _{\mathrm p} \left(\bar{\xi }_x-x\right)^2 \\&=& \frac{3}{4}\Sigma _{\mathrm p}\Omega _{\mathrm p} x^2-\frac{2\Sigma _{\mathrm p}}{\Omega _{\mathrm p}}\left(\frac{P^{\prime }}{\Sigma _{\mathrm p}}+\Phi _{\mathrm p}\right) + 3 x \Sigma _{\mathrm p}\bar{v}_y, \end{eqnarray}$$

(39)

$$\begin{eqnarray} \hphantom{\psi }\Rightarrow {\bar{\xi }_x}^2-2 x \bar{\xi }_x + \frac{8}{3\Omega _{\mathrm p}^2}\left(\frac{P^{\prime }}{\Sigma _{\mathrm p}}+\Phi _{\mathrm p}\right)-\frac{4 x \bar{v}_y}{\Omega _{\mathrm p}} = 0, \end{eqnarray}$$

(40)

$$\begin{eqnarray} \hphantom{\psi }\Rightarrow \bar{\xi }_x = x \pm \sqrt{x^2 -\frac{8}{3\Omega _{\mathrm p}^2}\left(\frac{P^{\prime }}{\Sigma _{\mathrm p}}+\Phi _{\mathrm p}\right) + \frac{4 x \bar{v}_y}{\Omega _{\mathrm p}}}, \end{eqnarray}$$

(41)

where the choice of |$+$| or |$-$| is determined by whether the fluid element has just undertaken a horseshoe turn or not. Here, |$P^{\prime }$| and |$\bar{v}_y$| may be taken to be the linear solutions to the 2D equations of motion, which we reduce to independent second order equations in Section 3. Equation (40) captures the essence of the corotation singularity, and what is meant by the ‘non-linear’ corotation torque. The linear theory such as that applied in Goldreich & Tremaine (1979) and Tanaka et al. (2002) may be thought to implicitly discard the quadratic term |${\bar{\xi }_x}^2$|⁠, which is non-negligible for |$x = \mathcal {O}(\sqrt{\frac{q}{h^3}}H)$|⁠. In this way, the linear equation for the particle displacement experiences a singularity at |$x = 0$|⁠, which must be resolved non-linearly (though there exist linear equations which remain valid at corotation for |$\bar{v}_x$|⁠, |$\bar{v}_y$|⁠, and |$P^{\prime }$|⁠). The singularity is introduced in linear analyses when the linearized azimuthal momentum equation, in our case equation (21b), is used to directly re-arrange for |$\bar{v}_y$|⁠, namely by integrating with respect to azimuthal coordinate y or |$\theta$| (which may look like dividing by |$\mathrm{ im}(\Omega -\Omega _\mathrm{ p})$| in Fourier space). More specifically, this ‘first integral’ equation for |$\bar{v}_y$| must become non-linear to remain valid near to corotation, just as the azimuthal integral of |$\bar{v}_x$| becomes large enough that inclusion of the quadratic term |$\bar{\xi }_x^2$| in equation (40) is necessary.

Importantly, this means that any materially conserved quantity |$Q(\psi) = Q(\bar{\xi }_x-x)$|⁠, will experience the same singularity as the radial displacement, which appears when the quadratic term in equation (40) is neglected. A gradient of (the materially conserved) PV over the coorbital region will therefore induce a singularity in the linear equations of motion. In order to resolve this singularity, the steady distribution of PV, which is advected by the flow induced by the planet within the horseshoe region, must be ascertained. In this sense, the corotation torque is ‘non-linear’. The corotation torque that arises depends very strongly on the geometry of the base flow, demanding an accurate model for the flow induced by a planet. We find this flow in the case of a low-mass planet in this paper.

2.7 2D mode orthogonality and torque decomposition

For |$\gamma \lt 2$|⁠, the unforced linearized equations of motion in a disc admit a ‘2D mode’ solution with |$v^{\prime }_z = 0$|⁠. This 2D mode is a member of a wider family of modes, whose axisymmetric members are discussed in Lubow & Pringle (1993). Here we demonstrate that the 2D mode is orthogonal to the remaining set of (non-axisymmetric) 3D disc modes, and show how the torque may be decomposed into 2D and 3D components. For |$\gamma \gt 2$|⁠, a different family of modes exists which does not include such a mode with |$v^{\prime }_z = 0$|⁠. Indeed, imposing |$v^{\prime }_z = 0$| in this case yields an unnormalizable mode of infinite energy. Appropriately normalized, the 2D mode may be written as

$$\begin{eqnarray} v^{\prime }_{x,0} = \bar{v}_x(x,y)\sqrt{2-\gamma }\exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}, \end{eqnarray}$$

(42a)

$$\begin{eqnarray} v^{\prime }_{y,0} = \bar{v}_y(x,y)\sqrt{2-\gamma }\exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}, \end{eqnarray}$$

(42b)

$$\begin{eqnarray} p^{\prime }_0 = P^{\prime }(x,y)\sqrt{2-\gamma }\frac{\rho _{\mathrm p}(z)}{\Sigma _{\mathrm p}}\exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}, \end{eqnarray}$$

(42c)

$$\begin{eqnarray} v^{\prime }_{z,0} = 0, \quad \rho ^{\prime }_0 = p^{\prime }_0/c_\gamma ^2, \end{eqnarray}$$

(42d)

where again

(43)

We may then decompose, for example, the radial velocity into its 2D mode component and an orthogonal complement, which we call |$v^{\prime }_{x,c}$|⁠, describing the rest of the disc’s 3D motions, which include for example inertial and gravity waves.

$$\begin{eqnarray} v^{\prime }_x = v^{\prime }_{x,0} + v^{\prime }_{x,c}. \end{eqnarray}$$

(44)

Now, |$v^{\prime }_{x,0}$| and |$v^{\prime }_{x,c}$| are everywhere orthogonal with respect to the inner product involving the density

$$\begin{eqnarray} \left\langle f,g \right\rangle = \frac{1}{\Sigma _{\mathrm p}}\int _{-\infty }^{\infty } \rho _{\mathrm p}(z) f(z) g(z) \, \mathrm{d}z. \end{eqnarray}$$

(45)

This may be seen from the relation for the vertical average |$\left\langle \, \cdots \right\rangle$| in terms of the density-weighted inner product

$$\begin{eqnarray} \left\langle \, \cdots \right\rangle \equiv \left\langle \, \cdots ,\frac{1}{\sqrt{\gamma }}\exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}\right\rangle . \end{eqnarray}$$

(46)

Specifically, the orthogonality follows via

$$\begin{eqnarray} \left\langle v^{\prime }_{x,0}, v^{\prime }_{x,c} \right\rangle &=& \left\langle v^{\prime }_{x,0}, v^{\prime }_{x}-v^{\prime }_{x,0} \right\rangle \\ &=& \sqrt{2-\gamma }\, \bar{v}_x\left\langle \exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}, v^{\prime }_{x}\right\rangle \\ && -\,(2-\gamma)\bar{v}_x^2\left\langle \exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)},\exp {\left(\tfrac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}\right\rangle \\ &=& \sqrt{\gamma (2-\gamma)}\, \bar{v}_x\left\langle v^{\prime }_x\right\rangle -(2-\gamma)\bar{v}_x^2\sqrt{\frac{\gamma }{2-\gamma }}\\ &=&0. \end{eqnarray}$$

As a result, we obtain the Parseval identity,

$$\begin{eqnarray} \left\langle v^{\prime }_x, v^{\prime }_y \right\rangle = \left\langle v^{\prime }_{x,0}, v^{\prime }_{y,0} \right\rangle + \left\langle v^{\prime }_{x,c}, v^{\prime }_{y,c}\right\rangle , \end{eqnarray}$$

(47)

as well as equivalent identities for any two variables drawn from the set |$\lbrace v^{\prime }_x, v^{\prime }_y, p^{\prime }/\rho _{\mathrm p}(z)\rbrace$|⁠. In particular, this allows us to decompose the radial angular momentum flux, |$F_{\mathrm A}$|⁠, into independent components associated with the 2D mode and 3D remainder, since

$$\begin{eqnarray} F_{\mathrm A} &\approx & r_{\mathrm p}\int\!\!\!\int \rho _{\mathrm p}(z) v^{\prime }_x v^{\prime }_y \, \mathrm{d}z\, \mathrm{d}y \\ &=& r_{\mathrm p}\int\!\!\!\int \rho _{\mathrm p}(z) v^{\prime }_{x,0} v^{\prime }_{y,0} \, \mathrm{d}z\, \mathrm{d}y + r_{\mathrm p}\int\!\!\!\int \rho _{\mathrm p}(z) v^{\prime }_{x,c} v^{\prime }_{y,c} \, \mathrm{d}z \, \mathrm{d}y \\ &=& \sqrt{\gamma (2-\gamma)}r_{\mathrm p} \Sigma _{\mathrm p}\int _{-\infty }^{\infty } \bar{v}_x\bar{v}_y \, \mathrm{d}y + F_{\mathrm A}^{\text{3D}} \equiv F_{\mathrm A}^{\text{2D}} + F_{\mathrm A}^{\text{3D}} \end{eqnarray}$$

(48)

It is in this sense that the torque on the planet for a general adiabatic index |$\gamma$| may be separated into 2D and 3D contributions. The factor of |$\sqrt{\gamma (2-\gamma)}$| (approximately 92 per cent for |$\gamma = 1.4$|⁠) may be shown to be the fraction of torque imparted into the 2D mode by an external potential (assumed to be z-independent) at a Lindblad resonance of order |$m \ll 1/h_\gamma$| [see e.g. equation (45) in Lubow & Ogilvie (1998) and appendix B1 of Bate et al. (2002)]. In this way, naïvely approximating the 3D flow velocities by their averaged values inadvertently recovers the total flux in this simplified scenario, as this approximation removes the factor of |$\sqrt{\gamma (2-\gamma)}$| from the flux expression. This follows from the applicability of resonant torque excitation theory to this regime.

This principle is at play in fig. 2 of Zhu et al. (2012). The top left panel depicts the simulated torque density excited in 3D discs which admit adiabatic perturbations to an isothermal background state. Notably, for |$x\gtrsim 2H_\gamma$|⁠, and for each value of |$\gamma$| considered, Zhu et al. (2012) recover torque densities matching the case |$\gamma = 1$| (in which the 2D mode is z-independent and contributes almost all of the torque). The torque density in this ‘outer’ region scales as |$x^{-4}$| (Goldreich & Tremaine 1980), so that their |$\gamma$|-dependent scaling of the x-axis with |$H_\gamma$| and y-axis with |$c_\gamma ^4$| does not meaningfully affect the outer torque density curves. This torque density agreement is achieved despite a considerable fraction of the flux being carried by gravity waves and inertial waves when |$\gamma \gt 1$|⁠. That is, the torque density for |$x\gtrsim 2H_\gamma$| matches that obtained via a naïve approximation of the 3D flow by its vertical average, as defined in equation (19).

There is, however, in the planet–disc interaction problem an important (indeed, a dominant) contribution from azimuthal modes with |$m \sim 1/h_\gamma$|⁠. The vertical profile of the planetary potential plays a key role here too. Indeed, inertial waves and gravity waves excited near the planet also impact the torque exerted on the planet significantly, contributing to |$F_{\mathrm A}^{\text{3D}}$|⁠. The problem of resolving the spectrum of gravity waves excited by an embedded planet and the flux they carry is very challenging, and has received only limited attention (Zhu et al. 2012; Lubow & Zhu 2014; McNally et al. 2020; Bae et al. 2021; Ziampras et al. 2023). Numerical approaches face difficulties resolving the complex and fine structure of the waves, and quantitative analytical approaches struggle since the gravity waves are not separable in the vertical direction.

3 REDUCTION TO INDEPENDENT SECOND-ORDER LINEAR EQUATIONS

In Section 2, we showed how the 3D equations governing the flow near a low-mass planet may be manipulated into the same form as the familiar 2D local equations for mass, PV, and momentum conservation

$$\begin{eqnarray} D P^{\prime } + \Sigma _{\mathrm p} c_\gamma ^2\left(\partial _x \bar{v}_x + \partial _y \bar{v}_y\right) = 0, \end{eqnarray}$$

(49a)

$$\begin{eqnarray} \partial _x \bar{v}_y-\partial _y \bar{v}_x-\frac{\Omega _{\mathrm p}}{2 c_\gamma ^2} \frac{P^{\prime }}{\Sigma _{\mathrm p}} = 0, \end{eqnarray}$$

(49b)

$$\begin{eqnarray} D\bar{v}_x-2\Omega _{\mathrm p} \bar{v}_y + \frac{1}{\Sigma _{\mathrm p}}\partial _x P^{\prime } = -\partial _x \Phi _{\mathrm p}, \end{eqnarray}$$

(49c)

$$\begin{eqnarray} D\bar{v}_y + \frac{1}{2}\Omega _{\mathrm p} \bar{v}_x + \frac{1}{\Sigma _{\mathrm p}}\partial _y P^{\prime } = -\partial _y \Phi _{\mathrm p}, \end{eqnarray}$$

(49d)

where |$D \equiv -\frac{3}{2}x\Omega _{\mathrm p}\partial _y$| is the leading order advective operator arising from the background shear flow. We now derive with no further approximations simple, second-order linear equations from these that describe both the spiral wave excitation and horseshoe streamlines. It is helpful to non-dimensionalize the equations. We first replace

$$\begin{eqnarray} x \rightarrow H_\gamma x, \quad y \rightarrow H_\gamma y, \end{eqnarray}$$

(50)

so that x and y now measure distances in units of ‘adiabatic scale heights’. We also let

$$\begin{eqnarray} \mathcal {D} = -\frac{3}{2}x\partial _{y}, \end{eqnarray}$$

(51)

and introduce the non-dimensional potential |$\hat{\phi }_{\mathrm p}$|

$$\begin{eqnarray} \Phi _{\mathrm p}(x,y) &=& \frac{q}{h_\gamma ^3}c_\gamma ^2\hat{\phi }_{\mathrm p}\left(s\right), \quad s = \sqrt{x^2+y^2},\\&& \Rightarrow \hat{\phi }_{\mathrm p}(s) = -\frac{\mathrm{e}^{\frac{1}{4}s^2}}{\sqrt{2\pi }}K_0\left(\tfrac{1}{4}s^2\right) \sim -\frac{1}{s} + \mathcal {O}\left(\frac{1}{s^3}\right). \end{eqnarray}$$

(52)

We define non-dimensional, scaled x- and y-velocity and ‘enthalpy’ perturbations |$u(x,y)$|⁠, |$v(x,y)$|⁠, and |$W(x,y)$| via

$$\begin{eqnarray} \bar{v}_x \equiv \frac{q}{h_\gamma ^3}c_\gamma u, \quad \bar{v}_y \equiv \frac{q}{h_\gamma ^3}c_\gamma v, \quad P^{\prime } \equiv \frac{q}{h_\gamma ^3}c_\gamma ^2\Sigma _{\mathrm p} W. \end{eqnarray}$$

(53)

Equations (49b), (49c), and (49d) become

$$\begin{eqnarray} \partial _x v-\partial _y u-\frac{1}{2}\chi = - \frac{1}{2}\hat{\phi }_{\mathrm p}, \end{eqnarray}$$

(54a)

$$\begin{eqnarray} \mathcal {D} u-2v + \partial _x \chi = 0, \end{eqnarray}$$

(54b)

$$\begin{eqnarray} \mathcal {D} v + \frac{1}{2}u + \partial _y \chi = 0, \end{eqnarray}$$

(54c)

where |$\chi \equiv W + \hat{\phi }_{\mathrm p}$|⁠. Additionally, ‘mass conservation’, (49a) becomes

$$\begin{eqnarray} \mathcal {D}\chi + \partial _x u + \partial _y v = \mathcal {D}\hat{\phi }_{\mathrm p}, \end{eqnarray}$$

(55)

which may also be derived from (54a), (54b), and (54c). It is clear at this stage that the flow is indeed linear so long as |$q \ll h_\gamma ^3$|⁠.

Remarkably, these equations are independent of |$\gamma$|! They are analogous to the equations considered by Paardekooper & Papaloizou (2009b); however, it is worth noting that the approximations adopted in their subsequent analysis led to three sources of inaccuracy. Firstly, the use of decaying boundary conditions instead of a radiation condition, as well as the use of a softened planetary potential both introduce discrepancies. Additionally, motivated by its validity in the case of test particle dynamics, they further approximated the above equations near to corotation by setting |$\mathcal {D} = 0$| (before then differentiating with respect to x). The resulting system omits large factors at corotation in the fluid case. Note that |$\partial _x\mathcal {D} \rightarrow -\tfrac{3}{2}\partial _y$| as |$x \rightarrow 0$|⁠, which does not vanish at corotation.

Before continuing the derivation, it’s helpful also to note the equation for PV conservation in differential form,

$$\begin{eqnarray} \mathcal {D}(\partial _x v-\partial _y u) + \frac{1}{2}(\partial _x u + \partial _y v) = 0. \end{eqnarray}$$

(56)

To proceed, we take |$\mathcal {D}$|(54b) and combine with (54c), and then take |$\mathcal {D}$|(54c) and combine with (54b). This yields

$$\begin{eqnarray} (\mathcal {D}^2 + 1)u = - \mathcal {D} \partial _x \chi -2 \partial _y \chi , \end{eqnarray}$$

(57a)

$$\begin{eqnarray} (\mathcal {D}^2 + 1)v = - \mathcal {D} \partial _y \chi + \frac{1}{2} \partial _x \chi . \end{eqnarray}$$

(57b)

Now, it may be shown from equations (54a) and (56) that

$$\begin{eqnarray} \mathcal {D} \partial _x \chi + 2 \partial _y \chi = \mathcal {D} \partial _x \hat{\phi }_{\mathrm p} + 2 \partial _y \hat{\phi }_{\mathrm p}-\nabla ^2 u + 3\partial _y(\chi -\hat{\phi }_{\mathrm p}), \end{eqnarray}$$

(58)

and

$$\begin{eqnarray} \mathcal {D} \partial _y \chi -\frac{1}{2} \partial _x \chi = \mathcal {D} \partial _y \hat{\phi }_{\mathrm p}-\frac{1}{2} \partial _x \hat{\phi }_{\mathrm p}-\nabla ^2 v. \end{eqnarray}$$

(59)

equations (57a) and (57b) become

$$\begin{eqnarray} (\mathcal {D}^2 + 1-\nabla ^2)u + 3 \partial _y \chi = \partial _y \hat{\phi }_{\mathrm p}-\mathcal {D} \partial _x \hat{\phi }_{\mathrm p}, \end{eqnarray}$$

(60a)

$$\begin{eqnarray} (\mathcal {D}^2 + 1-\nabla ^2)v = - \mathcal {D} \partial _y \hat{\phi }_{\mathrm p} + \frac{1}{2} \partial _x \hat{\phi }_{\mathrm p}. \end{eqnarray}$$

(60b)

Now, to find an equivalent equation for |$\chi$|⁠, we take |$\mathcal {D}^2$|(54a). We’re then able to use equation (56), the momentum equations (54b) & (54c) and the PV equation (54a) again to deduce

$$\begin{eqnarray} (\mathcal {D}^2 + 1-\nabla ^2)\chi + 3\partial _y u = (\mathcal {D}^2 + 1)\hat{\phi }_{\mathrm p}. \end{eqnarray}$$

(61)

We now define the linearized Riemann invariants¹|$J_\pm$|⁠:

$$\begin{eqnarray} J_\pm \equiv u \pm \chi . \end{eqnarray}$$

(62)

Combining equations (61) and (60a), it follows that

$$\begin{eqnarray} \left[\mathcal {D}^2 + 1 \pm 3 \partial _y-\nabla ^2\right]J_{\pm } = \left[\partial _y-\mathcal {D} \partial _x \pm (\mathcal {D}^2 + 1)\right]\hat{\phi }_{\mathrm p}, \end{eqnarray}$$

(63a)

$$\begin{eqnarray} \left[\mathcal {D}^2 + 1-\nabla ^2\right]v = - \mathcal {D} \partial _y \hat{\phi }_{\mathrm p} + \frac{1}{2}\partial _x \hat{\phi }_{\mathrm p}. \end{eqnarray}$$

(63b)

Now, equation (63b) is the parabolic cylinder equation² derived by Artymowicz (1993) for the azimuthal velocity perturbation, but equation (63a) is novel. Importantly, it provides a simple and non-singular description of the behaviour of the radial velocity and enthalpy at corotation. It describes the linear excitation of the profiles of the Riemann invariants of Goodman & Rafikov (2001), which are conserved non-linearly further from the planet. These equations also accurately capture the behaviour of the coorbital flow, including the horseshoe dynamics.

It is worth noting and crediting the similar equations derived by Heinemann & Papaloizou (2009), who studied the excitation of density waves via turbulence. In the appendix, we discuss the numerical solution of equations (63a) and (63b) to high accuracy, and we present these solutions in the next section.

The inner limit of equations (63a) and (63b) (that is, the simplified equations valid close to corotation taking |$x \ll H$|⁠) may be written as

$$\begin{eqnarray} \left[1 \pm 3 \partial _y-\nabla ^2\right]J_{\pm } = \left[\partial _y-\mathcal {D} \partial _x \pm (\mathcal {D}^2 + 1)\right]\hat{\phi }_{\mathrm p}, \end{eqnarray}$$

(64a)

$$\begin{eqnarray} \left[1-\nabla ^2\right]v = - \mathcal {D} \partial _y \hat{\phi }_{\mathrm p} + \frac{1}{2}\partial _x \hat{\phi }_{\mathrm p}. \end{eqnarray}$$

(64b)

Note that terms including |$\mathcal {D}$| acting on |$\hat{\phi }_{\mathrm p}$| may not be neglected, as in this case |$\hat{\phi }_{\mathrm p}$| has a (logarithmic) singularity at the origin, with the effect that such terms are non-negligible. For comparison, the equivalent equation from Paardekooper & Papaloizou (2009b) expressed in our notation reads

$$\begin{eqnarray} \left[1-\partial _x^2-4\partial _y^2\right]\chi = \phi _{\mathrm p}. \end{eqnarray}$$

(65)

3.1 Connection with resonant wave excitation theory

One may deduce directly from (63b) the locations of the ‘effective’ Lindblad resonances. In a pressureless disc composed of test particles, the Lindblad resonances are located at orbital radii such that the epicyclic frequency, |$\kappa (r)$|⁠, is an integer multiple of the interaction rate, |$\Omega -\Omega _{\mathrm p}$|⁠, of the test particles with the planet. In this way, successive interactions of a test particle with the planet will administer in phase ‘kicks’ to the test particle. In the 2D gas-dynamic case, inertio-acoustic wave disturbances instead oscillate with squared frequency |$\omega ^2 = \kappa ^2 + c^2 \left|{{\bf k}}\right|^2$| for wavevector |${{\bf k}}$| and sound speed c. In this way, disturbances of large wavenumber have a larger frequency, which has the effect that the Lindblad resonances of large order need not have a very small interaction rate. In this way, the Lindblad resonances all stand off and pile up a finite distance from the planet’s orbital radius, namely at |$x = \pm \tfrac{2}{3}H_\gamma$|⁠.

Furthermore, in the gas-dynamic case the contribution of the radial wavenumber |$k_x$| to the wave frequency |$\omega$| leads to the spreading out of each Lindblad resonance radially.

The resonance locations may equivalently be defined in the gas-dynamic context as where the solution for a given azimuthal mode changes from evanescent to wavelike. It is near this turning point, where the wave has zero frequency (and correspondingly the most time to be excited), that the forcing has most influence on the wave excitation. We note that the operator |$\mathcal {D}^2 + 1-\nabla ^2$| is hyperbolic in the wave-permitting region |$\left|x\right| \gt \tfrac{2}{3}H_\gamma$|⁠, and elliptic in the region where all modes are evanescent, namely |$\left|x\right| \lt \tfrac{2}{3}H_\gamma$|⁠. This signposts that the Lindblad resonances must all stand off a distance |$\tfrac{2}{3}H_\gamma$| from the planet’s orbital radius. More specifically, for the azimuthal disc mode proportional to |$\mathrm{e}^{\mathrm{i}m \theta }$|⁠, recalling the original definition |$y = r_{\mathrm p}\theta$| and reintroducing dimensions, (63b) becomes

$$\begin{eqnarray} \left(- \frac{9}{4} \Omega _{\mathrm p}^2 \frac{x^2}{r_{\mathrm p}^2}m^2 + \Omega _{\mathrm p}^2 + c_\gamma ^2\frac{m^2}{r_{\mathrm p}^2}-c_\gamma ^2\partial _x^2 \right)v = \cdots . \end{eqnarray}$$

(66)

We see that were it not for the acoustic terms in the above operator, namely |$c_\gamma ^2 m^2/r_{\mathrm p}^2$| and |$c_\gamma ^2\partial _x^2$|⁠, we would recover the positions of the close-in Lindblad resonances for test particles by setting this operator to 0, that is, |$r_{L,\pm m} = r_{\mathrm p}\left(1 \pm \tfrac{2}{3 m}\right)$|⁠. Appropriately, at these resonances, solving for v would involve division by 0. Including the acoustic terms yields a parabolic cylinder equation for |$v_y$| whose solution is evanescent for |$\left|x\right| \lt \tfrac{2}{3}H_\gamma \sqrt{1 + (h_\gamma m)^{-2}}$|⁠, and wavelike for |$\left|x\right| \gt \tfrac{2}{3}H_\gamma \sqrt{1 + (h_\gamma m)^{-2}}$|⁠. The effective Lindblad resonances close to the planet (with |$h_\gamma m \gtrsim 1$|⁠) therefore have locations

$$\begin{eqnarray} r_{L,\pm m} = r_{\mathrm p} \pm \frac{2}{3}H_\gamma \sqrt{1 + (h_\gamma m)^{-2}}. \end{eqnarray}$$

(67)

4 RESULTS

In this section, numerical solutions to equations (63a) and (63b) are presented and discussed. We compare our flow solution near to corotation with flows computed in 3D simulations, e.g. by Lega et al. (2015), Fung et al. (2015), and Masset & Benítez-Llambay (2016). We further compare the excited wave profiles with those found in 2D studies. We find good agreement with 3D one-sided Lindblad torque estimates, which are typically a factor of 2–3 lower than 2D values.

4.1 Flow streamlines and the horseshoe width

Real-space plots for the solutions of equations (63a) and (63b) for each flow variable (u, v, |$\chi$|⁠, |$J_+$|⁠, |$J_-$|⁠, as well as the stream function |$\psi$|⁠), are shown in Fig. 1. The planet is located at the origin of each figure. We imposed radiation boundary conditions in the x-direction, specifying that our numerical solution includes no incoming waves, and used a Fourier transform method in y. Our numerical method is accurate and tailored specifically to this problem, and the plotted solutions have an uncertainty of |$1 \times 10^{-5}$|⁠. Details on this numerical procedure are discussed in the appendix.

$Graphs of non-dimensionalized linearized Riemann invariants $J_+$ and $J_-$ as well as pseudo-enthalpy $\chi = W + \hat{\phi }_{\mathrm p}$, radial velocity u, and azimuthal velocity perturbation v near the planet. The scalings for the variables are given in equation (53). The planet is located at the centre of each panel. Data were obtained via solution of equations (63a) and (63b) using the method described in the appendix. The bottom right panel depicts streamlines of the 2D mode near the planet.$

Figure 1.

Graphs of non-dimensionalized linearized Riemann invariants |$J_+$| and |$J_-$| as well as pseudo-enthalpy |$\chi = W + \hat{\phi }_{\mathrm p}$|⁠, radial velocity u, and azimuthal velocity perturbation v near the planet. The scalings for the variables are given in equation (53). The planet is located at the centre of each panel. Data were obtained via solution of equations (63a) and (63b) using the method described in the appendix. The bottom right panel depicts streamlines of the 2D mode near the planet.

Open in new tab Download slide

Qualitatively, the vertically averaged flow field streamlines depicted in the bottom right panel of Fig. 1 are comparable to those obtained with a ‘softened’ potential with smoothing length |$b = 0.4H_\gamma$|⁠. Namely, they agree in their prediction for the horseshoe width, and the streamlines are only weakly affected by the presence of the density waves. This is to be expected, as this particular choice for b is known to match well the horseshoe width measured in 3D simulations (Lega et al. 2015; Masset & Benítez-Llambay2016).

In this case, there is however a noticeable difference in the excited density wave’s amplitude. We find via a second analogous calculation that the ‘softened’ potential excites a wave with peak amplitude 30 per cent greater than that of the 2D mode in both the profiles of |$J_+$| and v. This wave correspondingly transports an angular momentum flux inflated by 55 per cent, depicted in Fig. 5. This numerical discrepancy highlights the well-known result that softening prescriptions are unable to simultaneously capture both the corotation torque as well as the Lindblad torque accurately (Masset 2002). Indeed, in fig. 5 of Tanaka & Okada (2024) one observes that torque components vary by up to a factor of 3 as the smoothing length is varied between 0.3 and 0.7.

The horseshoe region semithickness |$x_s$| may be found from our solution by first noting that |$\chi$|⁠, u, and v decay rapidly to 0 as |$y \rightarrow \infty$| in the region |$|x| \lt \frac{2}{3}H_\gamma$|⁠. The leading order far-field horseshoe streamlines therefore have reflectional symmetry about |$x = 0$|⁠. We showed in Section 2.5 that the stream function in our vertically averaged flow may be expressed as

$$\begin{eqnarray} \frac{\psi }{\Sigma _{\mathrm p}} = \frac{3}{4}\Omega _{\mathrm p} x^2{-}\frac{2}{\Omega _{\mathrm p}}\left(\frac{P^{\prime }}{\Sigma _{\mathrm p}}{+}\Phi _{\mathrm p}\right) + 3 x \bar{v}_y + \mathcal {O}\left(\frac{q^2}{h_\gamma ^6}\frac{c_\gamma ^2}{\Omega _{\mathrm p}}\right). \end{eqnarray}$$

(68)

By evaluating the stream function at a stagnation point (located on the line |$x = 0$| in the limit |$q/h_\gamma ^3 \rightarrow 0$|⁠), and then far up- or downstream on the same streamline, we see that

$$\begin{eqnarray} \frac{3}{4}\Omega _{\mathrm p} x_s^2 = -\frac{2}{\Omega _{\mathrm p}}\frac{q}{h_\gamma ^3}c_\gamma ^2 \chi _s, \end{eqnarray}$$

(69)

where |$\chi _s$| is the value of the (non-dimensional) pseudo-enthalpy |$\chi$| at a stagnation point on the separatrix streamline. From our numerical solution, we have

$$\begin{eqnarray} \chi _s = -0.47115, \end{eqnarray}$$

(70)

and in the same limit |$q/h_\gamma ^3 \rightarrow 0$|⁠, the flow has three stagnation points, with coordinates |$x = 0$|⁠, |$y = \pm 0.439 H_\gamma$| and |$y = 0$|⁠. We therefore predict a horseshoe region half width for a low-mass planet of

$$\begin{eqnarray} x_s = 1.12089\sqrt{\frac{q}{h_\gamma ^3}}H_\gamma = 1.12089\sqrt{\frac{q}{h^3}}H\gamma ^{-1/4}, \end{eqnarray}$$

(71)

in good agreement with 2D simulations using a smoothing length |$b = 0.4H_\gamma$| (Paardekooper et al. 2010, equation 44), which match well the horseshoe width measured in 3D simulations for this choice of b (Lega et al. 2015; Masset & Benítez-Llambay2016).

The exact |$\gamma ^{-1/4}$| dependence is a non-trivial result, arising from the coincidence that the 2D mode of the potential, specified in equation (22), depends only on the length-scale |$H_\gamma$| and not H, even though the vertical structure of the background disc varies vertically on the length-scale H.

The phenomenon of three stagnation points appearing on the axis |$x = 0$| seems to be physical, indeed it is observed in 3D simulations (see e.g. figs 3 and 4 in Masset & Benítez-Llambay 2016), where in their isothermal fiducial run, the stagnation points have approximate positions |$y = - 0.36 H$|⁠, |$y = 0.53 H$|⁠, and |$y = 0$|⁠, having been displaced by the far-field radial pressure gradient).

In a 3D disc, in general the horseshoe width will depend on height. In the isothermal case however (when |$\gamma = 1$|⁠), there is no entropy stratification, and in fact close to the orbital radius of the planet, the flow is columnar (akin to Taylor–Proudman columns) since the advection is dominated by Coriolis and body forces. In this way, the horseshoe width becomes approximately height-independent, as is observed in 3D isothermal simulations, e.g. those performed by Fung et al. (2015) and Masset & Benítez-Llambay (2016).

In the adiabatic case, the Taylor–Proudman theorem no longer applies. Indeed, we expect a buoyancy wake within the downstream coorbital region, as observed by McNally et al. (2020). Instead, to gain traction here we appeal to the vertical structure of the 2D mode. As discussed in Section 2.7, this mode behaves as

$$\begin{eqnarray} && v^{\prime }_x \propto \exp {\left(\frac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}, \quad v^{\prime }_y \propto \exp {\left(\frac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}, \quad v^{\prime }_z = 0,\\ && p^{\prime } \propto \exp {\left(-\frac{z^2}{2 H_\gamma ^2}\right)}, \quad \rho ^{\prime } \propto \exp {\left(-\frac{z^2}{2 H_\gamma ^2}\right)}. \end{eqnarray}$$

(72)

Consequently (if we neglect the contributions from the other vertical modes), we expect the horseshoe width to increase with height above the mid-plane as

$$\begin{eqnarray} x_s(z) \approx 1.12 \sqrt{\frac{q}{h^3}}H \gamma ^{-1/4}\sqrt{2-\gamma } \exp {\left(\frac{(\gamma -1)z^2}{2 H_\gamma ^2}\right)}. \end{eqnarray}$$

(73)

The neglect of further 3D modes in this approximation is well-motivated. For example, inertial waves behave very differently to the 2D mode, with enthalpy perturbations typically modulated by |$\sqrt{x}$| near to corotation (Latter & Balbus 2009). Gravity waves are also confined near to corotation to the region above diagonal buoyancy resonances (Lubow & Zhu 2014).

While the exponential increase in horseshoe width (as well as perturbed flow velocities) with height is striking, we note that that the factor |$(\gamma -1)/2\gamma$| is small for reasonable values of |$\gamma$|⁠. Further, we should not be too concerned with the dynamics beyond 2–3 scale heights above the mid-plane, as over 95 per cent of the disc’s mass is contained within the first two scale heights. Moreover, our model’s assumptions, including linearity, and that the background state is vertically isothermal (discussed in Section 5.3) may begin to break down at such heights.

Fig. 2 depicts this vertical dependence for a few values of |$\gamma$|⁠. Curiously, the expression in equation (73) is poorly defined for |$\gamma \geqslant 2$|⁠; however, we do not expect |$\gamma$| to take such values in astrophysical discs. This oddity is due to the change in the character of the disc’s linear modes as |$\gamma$| increases through 2. For |$\gamma \geqslant 2$|⁠, e.g. there is no mode with |$v^{\prime }_z = 0$| (such a mode would contain infinite energy).

$Predicted vertical dependence of horseshoe region width $x_s(z)$ from 2D mode calculation for a planet with $M_{\mathrm p}/M_\star = q = 5\times 10^{-6}$ embedded within a disc with aspect ratio $h = 0.05$.$

Figure 2.

Predicted vertical dependence of horseshoe region width |$x_s(z)$| from 2D mode calculation for a planet with |$M_{\mathrm p}/M_\star = q = 5\times 10^{-6}$| embedded within a disc with aspect ratio |$h = 0.05$|⁠.

Open in new tab Download slide

In the isothermal case, our numerical value for the horseshoe width |$x_s = 1.12\sqrt{\frac{q}{h^3}}H$| is in reasonable agreement with 3D simulations (though inflated by a few per cent). Lega et al. (2015) and Masset & Benítez-Llambay (2016) both estimate a width via 3D simulation of |$x_s = 1.05\sqrt{\frac{q}{h^3}}H$|⁠. The reasons for the discrepancy between our result and theirs are perhaps two-fold. Most significantly, we use a fully local expression for the potential (equation 14), rather than performing a discrete azimuthal mode decomposition that takes into account the finite circumference of the disc. In this way, our calculation is most applicable to ultra-thin discs (e.g. in the problem of black hole migration within an AGN disc). We expect a relative discrepancy between our value and those obtained in 3D simulations of the order of |$h_\gamma \sim 0.05$|⁠. Furthermore, we exclude any non-linearity in our calculation. From equation (68), we expect a next order correction to |$x_s$| of relative size |$\mathcal {O}\left(q/h^3\right)$|⁠, which may also be of the order of a few per cent.

4.2 Velocity profiles at corotation

Of particular importance, especially when studying the dynamics of PV and entropy within the horseshoe region, is the dominant flow structure of the horseshoe streamlines. The radial extent of the horseshoe region is |$\mathcal {O}\left(\sqrt{\frac{q}{h^3}}H\right)$|⁠, whereas the flow velocity varies on the length-scale H. As a result, the flow perturbations induced by the planet within the horseshoe region (which inform the dominant flow structure) are well approximated simply by their values on the line |$x = 0$|⁠, that is, taking

$$\begin{eqnarray} v_{y,\text{hs}} \approx -\frac{3}{2}\Omega _{\mathrm p}x + v^{\prime }_y(0,y), \quad v_{x,\text{hs}} \approx v^{\prime }_x(0,y) . \end{eqnarray}$$

(74)

This was noted by Balmforth & Korycansky (2001), who studied the saturation of the corotation torque, and used matched asymptotics to compute Fourier coefficients for the flow velocity perturbations near to corotation. Fig. 3 depicts the non-dimensional radial and azimuthal velocity distributions at corotation, which satisfy (63a) and (63b). These correspond to taking cross-sections of the flow depicted in Fig. 1 along the line |$x=0$|⁠.

Figure 3.

Profiles of non-dimensionalized azimuthal and radial velocity perturbations at corotation, extracted from the numerical solution of equations (63a) and (63b).

Open in new tab Download slide

4.3 Wave evolution, profiles, and one-sided torque

In this section, we clarify in what sense |$J_+$| and |$J_-$| are linearized Riemann invariants, how they correspond to their approximately conserved non-linear counterparts, and discuss the excited wave profiles shown in Fig. 4 and corresponding one-sided Lindblad torque.

$Profiles of non-dimensional linearized Riemann invariants $J_+$ and $J_-$, pseudo-enthalpy $\chi = W + \hat{\phi }_{\mathrm p}$, and azimuthal velocity perturbation v for $x/H_\gamma = 1$, 2, 3, 4, and 5. $\eta$ is the characteristic coordinate defined in equation (78). The scalings for the variables are given in equation (53).$

Figure 4.

Profiles of non-dimensional linearized Riemann invariants |$J_+$| and |$J_-$|⁠, pseudo-enthalpy |$\chi = W + \hat{\phi }_{\mathrm p}$|⁠, and azimuthal velocity perturbation v for |$x/H_\gamma = 1$|⁠, 2, 3, 4, and 5. |$\eta$| is the characteristic coordinate defined in equation (78). The scalings for the variables are given in equation (53).

Open in new tab Download slide

The non-linear 2D theory developed by Goodman & Rafikov (2001) neglects the azimuthal velocity perturbation, as its linear solution decays as |$|x|^{-1/2}$| away from the planet. The resulting 2D system conserves on characteristic curves the Riemann invariants

$$\begin{eqnarray} R_\pm = u \pm \frac{2}{\gamma -1}(c-c_\gamma), \quad c = c_\gamma \left(\frac{\Sigma }{\Sigma _{\mathrm p}}\right)^{\frac{\gamma -1}{2}}. \end{eqnarray}$$

(75)

For small departures from the background state, the enthalpy perturbation obeys |$W^{\prime }/c_\gamma \approx c-c_\gamma$|⁠. We may write therefore

$$\begin{eqnarray} R_\pm \approx u \pm \frac{W^{\prime }}{c_\gamma }, \end{eqnarray}$$

(76)

equal to |$J_\pm$| up to an unimportant additive contribution of |$\Phi _{\mathrm p}/c_\gamma$|⁠. It therefore makes sense to interpret |$J_\pm$| as linearized Riemann invariants.

A WKB analysis, appropriately imposing outgoing wave boundary conditions, indicates that for each azimuthal Fourier mode (with wavenumber |$k_y \gt 0$|⁠),

$$\begin{eqnarray} \tilde{J}_{+} \sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}|x|^{+\displaystyle\tfrac{1}{2}}\exp {\left(+\mathrm{i}\displaystyle\frac{3}{4}\frac{\Omega _{\mathrm p} k_y}{c_\gamma }x^2\right)}, \quad \text{if} x \gt 0\\ |x|^{-\displaystyle\frac{3}{2}}\exp {\left(-\mathrm{i}\tfrac{3}{4}\displaystyle\frac{\Omega _{\mathrm p} k_y}{c_\gamma }x^2\right)}, \quad \text{if} x \lt 0. \end{array}\right. \end{eqnarray}$$

(77)

The behaviour of |$\tilde{J}_{-}$| can be determined from the symmetry |$J_+(x,y) = -J_-(-x,-y)$|⁠, or in Fourier space, |$\tilde{J}_-(x,k_y) = -\tilde{J}_+^{*}(-x,k_y)$|⁠. Note the strong |$|x|^{-3/2}$| decay of |$\tilde{J}_{+}$| in |$x \lt 0$|⁠. This may be thought of as a consequence of |$J_+$| being approximately conserved on characteristics emanating from |$x = -\infty$|⁠, where it is 0 (so that we have a simple wave). The consequence of this can be seen in Fig. 1, where we see each Riemann invariant nearly vanish one side of the line |$x = 0$|⁠. The decay is not quite as strong as |$|x|^{-3/2}$| in real space close to the planet because of the influence of wave excitation there.

The spiral arm locations are well approximated by the characteristics of the second order operators in equations (63a) and (63b). That is, the density waves in |$x \gt 0$| follow the curve |$\eta (x,y) = 0$|⁠, for characteristic coordinate

$$\begin{eqnarray} \eta = y-\frac{x}{3}\sqrt{\frac{9}{4}x^2-1}-\frac{1}{3}\cosh ^{-1}\left(\frac{3x}{2}\right). \end{eqnarray}$$

(78)

The accurate computation of these wave profiles is an important step in the determination of the total torque on the planet, contributing to the Lindblad torque. They are also needed to compute the locations of eventual shocking due to wave steepening. While the profiles obtained within a strictly local approximation yield cancelling inner and outer torques, the modification to the profiles from any asymmetry will lead to a net torque on the planet. The profiles excited by the vertically averaged |$\Phi _{\mathrm p}$|⁠, shown in Fig. 4, are qualitatively very similar to the profiles obtained in planet-driven wave evolution studies (Goodman & Rafikov 2001; Dong et al. 2011).

Our one-sided torque estimate is however more than a factor of 2 smaller than that obtained in the 2D studies conducted by Goodman & Rafikov (2001) and Dong et al. (2011). We attribute this to the point-mass and softened potentials used to force the 2D system in previous studies artificially inflating the excited wave amplitudes. This is the case even for ‘higher order’ expressions for the planet potential such as those used by Dong et al. (2011), and especially for smaller softening parameters. Indeed, we find that the angular momentum flux carried by the 2D mode, |$F_{\mathrm A}^{\text{2D}} = \sqrt{\gamma (2-\gamma)}r_{\mathrm p}\Sigma _{\mathrm p}\int _{-\infty }^{\infty }\bar{v}_x\bar{v}_y \mathrm{d}y$| (as defined in Section 2.7), far from the planet is

$$\begin{eqnarray} F_{\mathrm A}^{\text{2D}}(\infty) = 0.37 \sqrt{\gamma (2-\gamma)}\frac{(G M_{\mathrm p})^2 \Sigma _{\mathrm p} r_{\mathrm p} \Omega _{\mathrm p}}{c_\gamma ^3}, \end{eqnarray}$$

(79)

corresponding to a one-sided torque

$$\begin{eqnarray} T^{\text{2D}} \equiv F_{\mathrm A}^{\text{2D}} |_0^{\infty } = 0.34 \sqrt{\gamma (2-\gamma)}\frac{(G M_{\mathrm p})^2 \Sigma _{\mathrm p} r_{\mathrm p} \Omega _{\mathrm p}}{c_\gamma ^3}. \end{eqnarray}$$

(80)

For comparison, in the case of a point-mass potential in a 2D disc, Goldreich & Tremaine (1980) calculate a dimensionless numerical torque prefactor of 0.93. We have used here that |$F_{\mathrm A}^{\text{2D}}(0) = 0.03 \sqrt{\gamma (2-\gamma)}\frac{(G M_{\mathrm p})^2 \Sigma _{\mathrm p} r_{\mathrm p} \Omega _{\mathrm p}}{c_\gamma ^3}$|⁠, which is apparent from Fig. 5. The fact that |$F_{\mathrm A}^{\text{2D}}(0) \ne 0$| is due to a small fraction (4 per cent) of the flux carried by the outward-propagating wave being excited in |$x \lt 0$|⁠.

$The angular momentum flux, $F_{\mathrm A}^{\text{2D}}(x)$, carried by the density wave in a 3D disc (calculated by projection onto the 2D mode), compared with the flux excited in a 2D disc by a softened potential with smoothing length $b = 0.4H_\gamma$. The box shows a zoomed-in view of the flux profile, which exhibits the ‘negative torque phenomenon’ (Rafikov & Petrovich 2012). The factor of $\sqrt{\gamma (2-\gamma)}$ in the y-axis scale originates from equation (48), and may be ignored when interpreting the flux in the softened case.$

Figure 5.

The angular momentum flux, |$F_{\mathrm A}^{\text{2D}}(x)$|⁠, carried by the density wave in a 3D disc (calculated by projection onto the 2D mode), compared with the flux excited in a 2D disc by a softened potential with smoothing length |$b = 0.4H_\gamma$|⁠. The box shows a zoomed-in view of the flux profile, which exhibits the ‘negative torque phenomenon’ (Rafikov & Petrovich 2012). The factor of |$\sqrt{\gamma (2-\gamma)}$| in the y-axis scale originates from equation (48), and may be ignored when interpreting the flux in the softened case.

Open in new tab Download slide

In the case |$\gamma = 1$|⁠, the torque estimate (equation 80) is in agreement with the 3D (isothermal) simulations of D’Angelo & Lubow (2010) to within a few per cent (one can find a one-sided Lindblad torque estimate from the symmetric part of the cumulative torque profiles in their figs 7 and 8); it further matches precisely the torque found by Zhu et al. (2012). This is because in the case |$\gamma = 1$|⁠, all gravity waves are suppressed, and inertial waves carry very little flux (Tanaka et al. 2002).

Zhu et al. (2012) and Arzamasskiy, Zhu & Stone (2018) also note the discrepancy of a factor of 2 or 3 between previous 2D and 3D estimates of the one-sided torque. Judging from our result and fig. 2 of Zhu et al. (2012), it appears that in the case |$\gamma = 1.4$|⁠, the 2D mode contributes only 75 per cent of the one-sided Lindblad torque, with inertial waves and gravity waves carrying the majority of the remaining flux.

The angular momentum flux computed from our numerical solution, |$F_{\mathrm A}^{\text{2D}}(x)$|⁠, is plotted in Fig. 5. Of note is that the torque is almost entirely exerted within a few scale heights of the planet, which justifies a local approach to the problem. Further, there is a very small decrease in the flux beyond |$x \approx 3.5 H_\gamma$|⁠, which may be attributed to the ‘negative torque phenomenon’ (Rafikov & Petrovich 2012).

Finally, we are able to comment on the likely vertical profile of the shock front, and offer insight into the modification to the shocking length incorporating 3D effects. We note that for a freely propagating 2D mode, |$v^{\prime }_z = 0$|⁠, so that wave steepening may be considered as a purely horizontal effect. Further, the density wave wake and the gravity wave wake are spatially separated far from the planet (McNally et al. 2020), which suggests considering the wave steepening of an isolated 2D mode would provide a good model for 3D shock formation.

Using equation (42a) and comparing the maximum steepness of our 2D mode profile with that of Goodman & Rafikov (2001) and following their wave steepening calculation, we estimate the location of the initial shock in a 3D disc relative to their 2D estimate as

$$\begin{eqnarray} l_\text{sh}^{\text{3D}} \approx 1.2\frac{\mathrm{e}^{-\tfrac{\gamma -1}{5 H_\gamma ^2}z^2}}{(2-\gamma)^{1/5}}\cdot l_\text{sh}^{\text{2D}}. \end{eqnarray}$$

(81)

Here, we have simply exploited the |$A^{-2/5}$| dependence of the shocking length on the amplitude of the wave, and noted the reduced maximum steepness of our linear 2D mode surface density profile [which is around 65 per cent of that found by Goodman & Rafikov (2001) at |$x = 1.33 H_\gamma$|]. The vertical dependence is weak for physical values of |$\gamma$|⁠; for |$\gamma = 1.4$| it varies by only 5 per cent within the first scale height, |$H_{\mathrm p}$|⁠, above the disc mid-plane, though higher up a shock may be expected to form slightly closer to the planet. At much greater height however, as discussed in Section 4.1, our model’s basic assumptions may begin to break down.

5 DISCUSSION

We have already discussed many key items as and when they arose earlier in this paper, including

the nature and resolution of the corotation singularity in Section 2.6, including the manner in which it has a ‘non-linear’ resolution, and how and why it appears in previous linear analyses;
the relation between the commonly adopted 2D Euler equations and variables and their 3D counterparts (Section 2.4);
the appearance of and displacements in the locations of the Lindblad resonances within the framework of this paper in Section 3.1.

However, it is worth discussing caveats to this work as well as two further topics that warrant additional attention: the ‘rigorous’ softening of planetary potentials to reproduce 3D effects within a 2D disc model, and the role that this 2D mode plays in determining the torque on a low-mass planet.

5.1 Treating planetary potentials in 2D discs

5.1.1 Context: softening prescriptions

When studying planet–disc interactions in 2D, the planet’s potential should be modified to account for the influence of the vertical extent of the disc on the dynamics. Prescriptions such as

$$\begin{eqnarray} \Phi _{\mathrm p} = -\frac{G M_{\mathrm p}}{\sqrt{|{{\bf r}}-{{\bf r}}_{\mathrm p}|^2 + b^2}}, \quad b \sim H/2, \end{eqnarray}$$

(82)

are commonly adopted, as well as other similar prescriptions including the ‘fourth-order’ softened potential described in equation (16) of Dong et al. (2011). Each of these prescriptions, however, struggles to uniformly capture the impact of the disc’s vertical extent on the dynamics within a few scale heights of the planet.

Müller et al. (2012) compared equation (82) to the density-weighted average of a point-mass potential and concluded that the best choice for b varies depending on the distance from the planet, especially for separations |$|{{\bf r}}-{{\bf r}}_{\mathrm p}| \sim H$|⁠. They note, worryingly, that the choice of b can inform the direction of planetary migration, or in fragmentation simulations whether the disc fragments or not. For our purposes, adopting such a prescription can significantly inflate the excited wave flux [by as much as or more than a factor of 2, for both equation (82) and the ‘fourth order’ potential], and can impact strongly the horseshoe region width. As mentioned previously, the variation of net torque components can be by more than a factor of 3 as smoothing length is varied between 0.3 and 0.7 (Tanaka & Okada 2024, fig. 5).

5.1.2 A rigorous choice

In light of this, in this section we emphasize that equation (83) is a far better choice for the 2D potential. In particular, taking inspiration from Lubow & Pringle (1993), we proved in Section 2.3 that in the case of a low-mass planet, it is possible to directly derive 2D fluid equations governing averaged velocities and an effective ‘surface density’ that are forced precisely by equation (83). In this way, we argue that equation (83) (which for simplicity excludes the height-independent indirect term) is the optimal choice for the planetary potential in a 2D disc, especially when considering sub-thermal mass planets.

$$\begin{eqnarray} \Phi _\mathrm{ p} \equiv -\frac{G M_\mathrm{ p}}{H_\gamma }\frac{\mathrm{e}^{\frac{1}{4}s^2}}{\sqrt{2\pi }}K_0 \left(\tfrac{1}{4}s^2\right), \quad s = |{{\boldsymbol r}}-{{\boldsymbol r}}_\mathrm{ p}|/H_\gamma . \end{eqnarray}$$

(83)

(We discuss below how this may be implemented practically in simulations.)

Indeed, equation (83) is uniformly a very good approximation to the potential that forces the global 2D mode of the disc [which would simply involve |$s = |{{\bf r}}-{{\bf r}}_{\mathrm p}|/H_\gamma (r)$|], since far away from the planet

$$\begin{eqnarray} \Phi _\mathrm{ p} \sim -\frac{G M_\mathrm{ p}}{|{{\bf r}}-{{\bf r}}_\mathrm{ p}|}\left[1-\frac{1}{2s^2} + \frac{9}{8s^4} + \mathcal {O}\left(s^{-6}\right)\right]. \end{eqnarray}$$

(84)

That is, far away the potential returns to the familiar point-mass potential with small correction terms. In the same way, we reassuringly recover the point mass potential in the case of a ‘razor-thin’ disc. The higher order corrections may be thought to arise physically from the quadrupolar and higher order multipolar interactions between the vertically extended disc and the planet monopole. In seeking to converge more rapidly to a point-mass potential, the ‘fourth order’ softened potential of Dong et al. (2011) neglects the quadrupolar interaction term.

Near the planet, if we assume the background temperature |$T_0(r)$| [and therefore the scale height |$H(r)$|] varies slowly with radius, then the departure of the potential which forces the global 2D mode from the above expression for |$\Phi _{\mathrm p}$| is also small.

In Fig. 6, we compare graphically the appearance of the vertically averaged potential (equation 83) with a point-mass potential and ‘softened’ potentials of the form (equation 82), as well as their corresponding gradients. Note that |$\Phi _{\mathrm p}$| has a logarithmic singularity at the origin. More precisely,

$$\begin{eqnarray} \Phi _{\mathrm p} \sim \frac{G M_{\mathrm p}}{H_\gamma } \frac{2 \ln s}{\sqrt{2\pi }} \quad \text{as} s \rightarrow 0. \end{eqnarray}$$

(85)

$Comparison of vertically averaged potential $\Phi _{\mathrm p}$ from equation (22) which forces the 2D equations (21a) and (21b) with point-mass and softened potentials of the form (equation 82). The left panel shows acceleration for each prescription and the right panel shows gravitational potential. $s = |{{\boldsymbol r}}-{{\boldsymbol r}}_\mathrm{ p}|/H_\gamma$ measures the in-plane distance from the planet in adiabatic scale heights.$

Figure 6.

Comparison of vertically averaged potential |$\Phi _{\mathrm p}$| from equation (22) which forces the 2D equations (21a) and (21b) with point-mass and softened potentials of the form (equation 82). The left panel shows acceleration for each prescription and the right panel shows gravitational potential. |$s = |{{\boldsymbol r}}-{{\boldsymbol r}}_\mathrm{ p}|/H_\gamma$| measures the in-plane distance from the planet in adiabatic scale heights.

Open in new tab Download slide

A logarithmic singularity is a general feature of any weighted vertical average of the planetary potential.

Despite the singularity experienced by the averaged potential (and correspondingly the enthalpy) at the location of the planet, the linearity of the system is hardly compromised. To demonstrate this, we note that in reality, the singularity of the 3D potential will be truncated at the planet’s radius, |$R_\mathrm{ p}$|⁠. We define the small parameter

$$\begin{eqnarray} \varepsilon = R_\mathrm{ p}/H_\gamma . \end{eqnarray}$$

(86)

For Earth-like planets located at around 1–10 au, we expect |$10^{-5} \lesssim \varepsilon \lesssim 10^{-3}$|⁠. Suppose we truncate the singularity of the 3D planetary potential (excluding the indirect term) at the planet’s radius prior to the vertical averaging as

$$\begin{eqnarray} \phi _{\mathrm p} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-\displaystyle\frac{G M_\mathrm{ p}}{\left|{{\boldsymbol r}}-{{\boldsymbol r}}_\mathrm{ p}\right|}, \quad \left|{{\boldsymbol r}}-{{\boldsymbol r}}_\mathrm{ p}\right| \gt R_\mathrm{ p},\\ -\displaystyle\frac{G M_\mathrm{ p}}{R_\mathrm{ p}}, \quad \left|{{\boldsymbol r}}-{{\boldsymbol r}}_\mathrm{ p}\right| \leqslant R_\mathrm{ p}. \end{array}\right. \end{eqnarray}$$

(87)

In this case,

$$\begin{eqnarray} \Phi _{\mathrm p}({{\bf 0}}) \sim \frac{G M_{\mathrm p}}{H_\gamma }\frac{2\ln \varepsilon }{\sqrt{2\pi }}. \end{eqnarray}$$

(88)

However, |$|2\ln \varepsilon /\sqrt{2\pi }| \sim 7 \lt h_\gamma ^3/q$| by assumption, so that the 2D system (including the enthalpy W) everywhere constitutes an approximately linear perturbation to the background state. We note further that the solutions for u, v, and |$\chi$| shown in Section 4 are very weakly affected by the (integrable) logarithmic singularity in the potential.

This then suggests that for numerical applications it would be very reasonable (namely to avoid numerical divergences and time-step issues), to ‘soften’ the 2D potential (equation 83) by taking |$s^2 = \varepsilon ^2 + |{{\bf r}}-{{\bf r}}_{\mathrm p}|^2/H_\gamma ^2$|⁠, though here a practical choice for |$\varepsilon$| for grid-based simulations in which the planetary radius is unresolved would be of the order of the grid spacing divided by |$H_\gamma$|⁠.

Regarding numerical implementation, we ought to address the subtlety that when directly evaluating the potential for large arguments s, we multiply both exponentially large and small numbers. However, efficient C functions for the exponentially scaled Bessel functions |$\mathrm{e}^xK_0(x)$| and |$\mathrm{e}^xK_1(x)$| (which both appear in the expression for the gradient of |$\Phi _{\mathrm p}$|⁠) already exist within the Cephes library³ (Moshier 1989). Alternatively, the vertically averaged potential in (equation 83) and its derivative are functions of one parameter only, so that a lookup table would provide an efficient method for evaluation.

5.2 Torques on low-mass planets

The torque exerted on the planet by the disc is of particular importance, as it informs the rate and direction of planet migration. The solution presented in this paper constitutes the leading order flow induced by the planet, which is rotationally symmetric and so exerts no net torque on the planet. The next-order correction to this flow is a factor of |$h_\gamma$| smaller, includes geometric asymmetries as well as radial variations in the properties of the background disc, and in general leads to a net torque on the planet. The majority of the torque exerted on the planet is done so locally. In this way, a radially local model of the planet–disc interaction is sufficient to determine the torque on the planet, which must depend linearly on the local gradients of the disc’s properties.

5.2.1 Corotation torque and critical layer at corotation

One key complicating issue is that of the horseshoe streamlines. They form closed loops when the full azimuthal extent of the flow is considered (and when the planet and disc are not migrating or drifting radially relative to each other). In particular, materially conserved quantities (such as Ertel’s PV and the entropy) are advected from inner to outer disc radii, and vice versa. Therefore, in the absence of strong viscosity or dissipation, the properties of the background disc are not restored far (azimuthally) from the planet.

This effect gives rise to an additional torque known as the corotation torque. To determine this torque, it is necessary first to solve for the distributions of these conserved quantities within the critical layer, which is called the horseshoe region, with appropriate prescriptions for viscosity, heating, etc. A key ingredient in this solution is an accurate description of the leading order flow within the corotation region. As mentioned in the introduction, there has already been a wealth of research devoted to the study of corotation torques. However, there is still more to be explored here, including in particular the impact that the vertical extent of the disc has on the torque. Part of this impact is captured in the behaviour of the 2D mode of the leading order flow, whose velocity profiles at corotation are shown in Fig. 3. These profiles vary only slightly over the corotation region, and may be used to enable cheap numerical studies of the influence of a range of physics on the corotation torque.

5.3 Caveats

We point out three main caveats to this work, which are detailed in the subsections below.

5.3.1 Omitted 3D wave modes

The solution presented in this paper is not fully 3D. We solved for the behaviour of the z-dependent 2D mode of the disc, which ought to be superposed with a spectrum of inertial and internal gravity waves. These contribute a significant fraction (upwards of 20 per cent for |$\gamma = 1.4$|⁠) of the one-sided torque on the planet, as discussed in Section 4.3.

5.3.2 Impact of turbulence and magnetic fields

Many physical phenomena may act to disrupt the solution for the disc’s 2D mode. The vertical shear in the 2D mode might be affected by processes that provide greater coupling between different layers in the disc, e.g. turbulence, viscosity, or a vertical magnetic field. Favourable conditions include a large plasma |$\beta \gg 1$| and low level of ionization, as well as a low Reynolds stress, captured by the condition for the viscosity parameter |$\alpha \ll 1$| (Shakura & Sunyaev 1973). It is widely believed that except for the innermost regions, protoplanetary discs do predominantly exhibit |$\beta \gg 1$| (Lesur 2021), and recent numerical and observational studies point towards a smaller turbulent viscosity than previously thought. Non-ideal magnetohydrodynamic effects such as ambipolar diffusion and Ohmic resistivity act to suppress turbulence even in the disc’s surface layers (Turner et al. 2014), and to decouple the fluid from the magnetic field. Indeed, Pinte et al. (2016) and later Villenave et al. (2022) find that observations of the discs around HL Tau and Oph163131 are consistent with small viscosity parameters |$\alpha \sim 10^{-4}$| and |$\alpha \sim 10^{-5}$|⁠, respectively.

5.3.3 Thermodynamic assumptions

The validity of the thermodynamic model adopted in this paper (which simply constitutes a vertically isothermal background state permitting adiabatic perturbations) depends critically on the thermal relaxation time-scale for the disc’s gas. If the cooling is too rapid, then the adiabatic assumption breaks down too slow and the disc’s background state will not have reached a thermal equilibrium as the disc evolves. The disc’s outer layers are typically hotter than its interior due to the stellar irradiation, which is not able to penetrate the optically thick disc interior. The molecular gas that comprises the disc interior is optically thick to its own radiative emission; its primary means of cooling is via the surrounding dust grains. The gas must first transfer its thermal energy to the dust, which is then able to radiate away the excess energy more efficiently (though the disc is not necessarily optically thin to this emission). As pointed out by Malygin et al. (2017), Barranco, Pei & Marcus (2018), Pfeil & Klahr (2019), and Bae et al. (2021), infrequent gas–dust collisions often act as a bottleneck for the thermal relaxation of the gas, especially in its surface layers. Bae et al. (2021) predict that typical gas thermal relaxation time-scales are tens of orbits even at 10 au, corresponding to a cooling parameter |$\beta \equiv 2\pi t_{\text{relax}}/t_{\text{orb}} \gtrsim 10^{1.5}$|⁠.

This slow cooling offers favourable conditions for the adiabatic model governing planet-induced perturbations to represent a good approximation, though the cooling is still more rapid than the disc’s evolution time-scale. Indeed, the downstream buoyancy wake excited by a giant planet orbiting at 90 au is believed to be the source of the tightly wound spirals observed in TW Hya (Teague et al. 2019; Bae et al. 2021). The existence of such a signature adds to the credibility of the adiabatic model, as it implies |$t_\text{relax} \gg N_z^{-1}$| (where |$N_z$| the Brunt–Väisälä frequency), necessary for the buoyancy wake to develop. It is likely, however, that cooling and thermal diffusion play important roles on the much longer time-scale associated with the librating horseshoe motion. While this has only a small effect on the dominant flow behaviour found in this paper, it has important consequences for the entropy and PV distributions within the horseshoe region, and correspondingly implications for the corotation torque.

Insight into the vertical temperature structure of the outer regions of protoplanetary discs is offered by the emission from different CO isotopes that trace different disc altitudes. Observational studies typically find a flat temperature plateau near the disc mid-plane, with an increase in temperature in the disc’s upper layers, a few pressure scale heights above the mid-plane (Dartois, Dutrey & Guilloteau 2003; Law et al. 2024). This is consistent with the physical picture discussed above. The disc’s surface temperature is typically a factor of 2 larger than its mid-plane temperature. While there is no exact analogue for the 2D projection procedure described in Section 2.3 when the temperature of the background state varies with height, the majority of the disc’s mass is contained within the region of temperature plateau. While an oversimplification, this suggests therefore that the vertically isothermal model for the background disc may be expected to give quantitatively reasonable results.

6 CONCLUSIONS

For linear adiabatic perturbations to a vertically isothermal protoplanetary disc, there is a particular vertical average of the 3D gas dynamic equations that exactly yields the familiar 2D linear system, comprised of averaged flow velocities, but an effective surface density, pressure, and planetary potential. This averaging process is closely related to the projection operator onto the Lubow–Pringle 2D mode of the disc (Lubow & Pringle 1993). Importantly, when compared with more traditional softening prescriptions for planetary potentials, adopting the averaged potential (equation 22) provides a more rigorous and accurate, parameter-free method to modify the potential of an embedded planet to account for 3D effects.

In this paper, we presented the solution for the 2D mode of the flow excited by a low-mass planet on a circular orbit, whose features include a spiral wake as well as horseshoe streamlines within the coorbital region. We derived non-singular, independent second-order equations for each flow variable, equations (63a) and (63b). These include novel parabolic cylinder equations for linear combinations of the radial velocity and enthalpy, and offer a correction to the model of the coorbital flow proposed by Paardekooper & Papaloizou (2009b).

The 2D mode (so called as it has the property |$v_z = 0$|⁠, though in general it is z-dependent) is a member of the wider family of 3D disc modes. It provides an interpretation for 2D disc models, and plays a dominant role in the response of the disc to the planet, particularly at corotation. We find that in the limit of an ultra-thin disc, the vertically averaged horseshoe width is |$x_s = 1.12\sqrt{\frac{q}{h^3}}H\gamma ^{-1/4}$|⁠. Taking only the 2D mode contribution to the flow predicts a horseshoe width that grows with height above the disc mid-plane, as described in equation (73) and Fig. 2.

The flow in the corotation region is well approximated by superposing the background shear flow with the perturbed flow fields evaluated at corotation, |$x = 0$|⁠. This coorbital flow solution, depicted in Fig. 3, may be used to inexpensively simulate the impact of various physics on the corotation torque, including diffusive and migration-feedback effects.

Our approach also allowed us to capture accurately the wave angular momentum flux transported by the spiral density wave in a 3D disc, which is reduced by a factor of 2 or 3 compared with previous 2D estimates. We used the profiles of this wave to estimate the location at which it first shocks in a 3D disc, as a function of height above the disc mid-plane.

We demonstrated that the 2D mode is orthogonal to the wider family of permitted 3D motions, which include, for example, gravity waves and inertial waves. As such, the torque on the planet decomposes into the sum of separate 2D and 3D components. We omit from this paper any discussion of the remaining excited spectrum of 3D wave modes. Resolving and understanding the gravity wave spectrum is a challenging and understudied, yet important problem that requires future attention, not only because of the angular momentum that they transport but also for their observational signature.

ACKNOWLEDGEMENTS

This research was supported by a Science and Technology Facilities Council (STFC) PhD studentship (grant number 2750631). We are very grateful to the referee, Clément Baruteau, for a thorough report and very helpful suggestions that have improved the paper.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

The interpretation of |$J_\pm$| as linearized Riemann invariants is clarified in Section 4.3.

More specifically, it becomes the parabolic cylinder equation following an azimuthal mode decomposition or Fourier transform.

https://www.netlib.org/cephes/

REFERENCES

Abramowitz

Stegun

I. A.

1972

Handbook of Mathematical Functions

National Bureau of Standards

Washington, DC

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

ALMA Partnership

2015

ApJ

808

Month:	Total Views:
September 2024	113
October 2024	118
November 2024	43
December 2024	34
January 2025	49
February 2025	22
March 2025	34
April 2025	15
May 2025	8

Article Contents

Horseshoes and spiral waves: capturing the 3D flow induced by a low-mass planet analytically

ABSTRACT

1 INTRODUCTION

2 GOVERNING EQUATIONS

2.1 Unperturbed state

2.2 Perturbation equations

2.3 Projection onto a 2D flow

2.4 Relation to 2D Euler equations

2.5 Stream function and Bernoulli invariant

2.6 On the corotation singularity

2.7 2D mode orthogonality and torque decomposition

3 REDUCTION TO INDEPENDENT SECOND-ORDER LINEAR EQUATIONS

3.1 Connection with resonant wave excitation theory

4 RESULTS

4.1 Flow streamlines and the horseshoe width

4.2 Velocity profiles at corotation

4.3 Wave evolution, profiles, and one-sided torque

5 DISCUSSION

5.1 Treating planetary potentials in 2D discs

5.1.1 Context: softening prescriptions

5.1.2 A rigorous choice

5.2 Torques on low-mass planets

5.2.1 Corotation torque and critical layer at corotation

5.3 Caveats

5.3.1 Omitted 3D wave modes

5.3.2 Impact of turbulence and magnetic fields

5.3.3 Thermodynamic assumptions

6 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: NUMERICAL PROCEDURE

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only