Horseshoes and spiral waves: capturing the 3D flow induced by a low-mass planet analytically
Joshua J. Brown Gordon I. Ogilvie


firstpage–lastpage
Pointwise and uniform bounds for functions of the Laplacian on non-compact symmetric spaces
    [
    September 9, 2024
===========================================================================================




§ 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 which 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 which transport angular momentum radially away from the planet, and ‘horseshoe streamlines’ within the co-orbital region. We derive simple linear equations governing the flow which locally capture both effects faithfully simultaneously. We present an accurate co-orbital flow solution allowing for inexpensive future study of corotation torques, and predict the vertical structure of the co-orbital 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 to 3 times lower than values from previous 2D simulations.





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








§ INTRODUCTION



Discs of gas in orbital motion around a massive central body are found in many astronomical contexts. Objects embedded within these discs, for example planets embedded in circumstellar discs, or stellar-mass black holes within AGN discs, are in general subject to strong gravitational interactions with the gas and dust which comprise the discs. Indeed, signatures such as gaps, rings and spirals observed in discs' emission are highly suggestive of large embedded planets <cit.>. 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 which excite predominantly linear perturbations to the flow in the disc. Not only is it a crucial ingredient in theory of planet formation and population synthesis <cit.>, but it also may play a key role in understanding the rapid merger rate of stellar-mass black hole (sBH) binaries inferred from LIGO detections <cit.>. Whilst 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 two-dimensional analyses from early work <cit.> 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. <cit.>, and more recently <cit.>, extended the 2D linear calculation to 3D locally isothermal discs, arriving at a torque formula in good agreement with the numerical simulations of <cit.>, which drives rapid inward migration on a time-scale of ∼ 10^5 years 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 ≳ 10^6 years <cit.>. Furthermore, since the lifetime of the disc is ≲ 10^7 years, 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 semi-major axes among detected low-mass exoplanets from the linear type I torque estimates <cit.>.

However, as pointed out by <cit.>, 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 <cit.> and <cit.>), 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 <cit.>, 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 <cit.>.

The torques which 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 by more than a factor of 3 for different reasonable choices of b (e.g. <cit.>, figure 5); such varied behaviour for different softening lengths is an issue not limited to the study of planetary torques <cit.>. 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 <cit.>, despite their analysis suffering some limitations (which we point out in section <ref>).

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 which permit adiabatic perturbations, we show that there exists a particular vertical average of the equations of motion which reproduces the commonly adopted 2D fluid equations. This averaging procedure is closely related to the operator which projects onto the Lubow-Pringle 2D mode <cit.>. This paper is concerned with the non-axisymmetric generalisation 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 <cit.>. As <cit.> 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 <ref>, 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 <ref>. In section <ref>, we present the flow solution. We discuss our findings in section <ref>, and draw our conclusions in section <ref>.



§ GOVERNING EQUATIONS



For the remainder of the paper, we will assume two quantities to be small. Firstly, we assume the disc's aspect ratio h = H/r to be much smaller than 1. Secondly, we assume the planet's mass M_p to be much smaller than the `thermal mass' <cit.>. That is,

    q ≡M_p/M_⋆≪M_th/M_⋆ = h^3,

where M_⋆ 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 2 small parameters, namely h and q/h^3; we'll exploit the fact that both are small to make analytical progress.



 §.§ Unperturbed state


Before continuing, it's 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_⋆, fixed at the centre of our frame, and is inviscid and non-self-gravitating to a first approximation. The gas which comprises it solves the steady Euler equations and ideal gas equation,


    u_0·∇u_0 = - 1/ρ_0∇ p_0 - ∇Φ_0,
       ∇·(ρ_0 u_0) = 0,
       p_0/ρ_0 = k_B/μ̅T_0,


where u_0, ρ_0, p_0 and T_0 are the velocity, density, pressure and temperature of the background disc. In addition, k_B is Boltzmann's constant, and μ̅ the mean molecular mass. We introduce the cylindrical coordinates (r, θ, z), with the z-axis normal to the plane of the disc, so that 

    Φ_0 = -G M_⋆/√(r^2+z^2)

is the gravitational potential of the central star, and the velocity of the background state may be written as u_0 = rΩ(r,z) e_θ.

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 idealised, 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'll 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 <ref>. We define the isothermal sound speed, c_s(r), via

    c_s^2 = k_B/μ̅T_0.

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

    h ≡H/r≡c_s/rΩ_K,

where Ω_K is the Keplerian angular frequency, satisfying

    G M_⋆/r^3 = Ω_K^2.

We assume h ≪ 1, so that the disc is thin. The angular velocity of the gas in the unperturbed disc satisfies

    Ω = Ω_K(1 + 𝒪(h^2)).

The density and pressure then satisfy

    ρ_0 = p_0/c_s^2 = Σ_0(r)/√(2) Hexp(-z^2/2H^2),

where Σ_0(r) is the surface density.



 §.§ Perturbation equations


We now introduce a planet of mass M_p on a Keplerian circular orbit of radius r_p and angular frequency Ω_p = Ω_K(r_p)√(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 v = u - rΩ_p e_θ, and the relevant Euler equations become


    v·∇v  + 2Ω_p ×v = - 1/ρ∇ p - ∇Φ_t - ∇Ψ_p,
       ∇·(ρv) = 0,
       v·∇(p ρ^-γ) = 0,


where Φ_t = Φ_0(r,z) - 1/2(r^2- 3r_p^2)Ω_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

    Ψ_p = - G M_p/√(r_p^2 + r^2 - 2 r r_p cosθ + z^2) + q r_p Ω_p^2  r cosθ.

We now define local quasi-Cartesian coordinates x = r - r_p, y = r_pθ, and perturbed variables


    ρ' = ρ - ρ_0,
       p' = p - p_0,
       v'_x = v_r,
       v'_y = v_θ - r(Ω - Ω_p) = u_θ - r Ω,
       v'_z = v_z.


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

    D = -3/2xΩ_p_y,

and upon subtraction of the background state solution from the Euler equations (<ref>), we find the balance at the next order in h (assuming x = 𝒪(H) and y = 𝒪(H)) is given by


    D v'_x - 2Ω_p v'_y + _x p'/ρ_p = -_x ϕ_p,
       D v'_y + 1/2Ω_p v'_x + _y p'/ρ_p = -_y ϕ_p,
       D v'_z + Ω_p^2 z ρ'/ρ_p + _z p'/ρ_p = -_z ϕ_p,
       D ρ' + ρ_p_xv'_x + ρ_p_yv'_y + _z(ρ_pv'_z) = 0,
       D(p' - γ c_p^2ρ') + (γ - 1)Ω_p^2 z ρ_p v'_z = 0,


where c_p = c_s(r_p),

    ϕ_p = -G M_p/√(x^2+y^2+z^2),

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

    ρ_p(z) ≡Σ_p/√(2) H_pexp(-z^2/2H_p^2),

for Σ_p = Σ_0(r_p) and H_p = H(r_p). The system of equations (<ref>) describes locally the flow excited by a planet embedded in a stratified 3D disc. We could in theory relax the assumption y ≲ H by reintroducing the azimuthally global expression for the planet's potential, Ψ_p, in place of ϕ_p, and applying periodic boundary conditions at y = ± r_p. The corresponding correction however is only of relative size 𝒪(h), so that local system involving ϕ_p (which we adopt for the remainder of the paper) becomes exact in the limit h → 0.

The system (<ref>) governs the excitation of the spiral density waves by the planet, permits downstream gravity waves (e.g. <cit.>), and also describes the horseshoe trajectories followed by fluid elements close to the planet's orbital radius. As noted by several authors, for example by <cit.>, the density wave excitation is confined to the region extending only a few scale heights from the planet, though (<ref>) does not capture the asymmetry between the inner and outer spiral wakes. Importantly, as we'll demonstrate in the next section, (<ref>) also has the elegant property that it may be vertically integrated to derive exact 2D flow equations.



 §.§ Projection onto a 2D flow



Remarkably, under a particular choice of vertical averaging, the system of equations (<ref>) 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 (<ref>)).

The averaging procedure defined below is closely related to the operator which projects onto the 2D mode of the disc found by <cit.>. 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 γ < 2, and is often referred to as `2D', since it has the property that v'_z = 0, despite possessing a vertical dependence.

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

    c_γ = √(γ) c_p,    H_γ = c_γ/Ω_p,    h_γ = H_γ/r_p.

We combine (<ref>) and (<ref>) to eliminate ρ':

    c_γ^2 × (<ref>) + (<ref>)
     D p' + c_γ^2 ρ_p(_x v'_x + _y v'_y) 
     
        + (γ - 1)Ω_p^2 z ρ_p v'_z + c_γ^2_z(ρ_pv'_z) = 0


    D p'/ρ_p + c_γ^2 (_x v'_x + _y v'_y) +
     c_γ^2 exp(z^22 H_γ^2)_z[ v'_z exp(-z^22 H_γ^2)] = 0.

We now multiply (<ref>) by the factor exp(-z^22 H_γ^2)∝ρ_p^1/γ and integrate. (Incidentally, the function ρ_p^1/γ 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

    D P' + Σ_p c_γ^2(_x v̅_x + _y v̅_y) = 0,

where we have defined

    v̅_x ≡⟨ v'_x⟩,   v̅_y ≡⟨ v'_y⟩,    P' ≡Σ_p ⟨ p'/ρ_p(z)⟩,

for vertical average ⟨ ⋯⟩ defined via:

    ⟨ X ⟩≡1/√(2)H_γ∫_-∞^∞Xexp(-z^2/2H_γ^2) z.

We may similarly apply the averaging procedure ⟨ ⋯⟩ defined in (<ref>) to equations (<ref>) and (<ref>) to obtain the 2D momentum equations


    Dv̅_x - 2Ω_p v̅_y + 1/Σ_p_x P' = -_x Φ_p,
       Dv̅_y + 1/2Ω_p v̅_x + 1/Σ_p_y P' = -_y Φ_p,


where now Φ_p is precisely defined as ⟨ϕ_p ⟩, that is,

    Φ_p ≡⟨ -G M_p/√(x^2 + y^2 + z^2)⟩ = -G M_p/H_γ^1/4s^2/√(2)K_0(14s^2),

with s = √(x^2 + y^2)/H_γ, and K_0(z) the modified Bessel function. Note that in a global 2D disc model, the above expression remains valid upon redefining s = |r - r_p|/H_γ and reintroducing the indirect term (this matter is discussed in greater detail in section <ref>). This prescription for the potential is therefore generally applicable to 2D disc models. Reassuringly, at large distances, the above expression becomes

    Φ_p = -G M_p/√(x^2 + y^2)(1 + 𝒪(H_γ^2/x^2+y^2)),

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 γ.

We remark that we did not make use of equation (<ref>) 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 <cit.>. 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.



 §.§ 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 <ref>) 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 (for example 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

    Σ̅≡Σ_p + P'/c_γ^2,
    v̅≡v̅_xe_x +  (-32Ω_p x + v̅_y)e_y,

then we see that our vertically averaged flow equations (<ref>), (<ref>) and (<ref>) are linearisations of the familiar 2D barotropic model,


    v̅·∇v̅  + 2Ω_p ×v̅ = - 1/Σ̅∇P̅ - ∇Φ̅_t - ∇Φ_p,
       ∇·(Σ̅v̅) = 0,
       P̅ = K Σ̅^γ,


where Φ̅_t = -32Ω_p^2 x^2, and  K = c_p^2Σ̅_p^1-γ uniformly, so that P̅ = c_γ^2 Σ̅. We remind the reader of the formal definitions of P', v̅_x and v̅_y, given in equation (<ref>), and our convention c_γ^2 = γ c_p^2. (Note we may treat (<ref>) as an exact definition of P̅.)

We note that the error in the system (<ref>) scales as the size of the quadratic non-linear terms, that is, (q/h_γ^3)^2. This is notably far smaller than the error introduced by any alternative averaging process, which would be of order q/h_γ^3. Recall we assumed q ≪ h_γ^3 so that our 3D flow was linear, and that the errors introduced in the linearisation of the 3D system also scaled as (q/h_γ^3)^2. Furthermore, whilst we've addressed here only a local model of disc dynamics, this approach may be generalised to a global disc. It's sufficient for the background disc to satisfy a locally isothermal equation of state.

It's worth mentioning that the `surface density' Σ̅ in the 2D model must really be thought of in terms of the vertically averaged pressure perturbation. Moreover, equation (<ref>) doesn't technically enforce physical mass conservation; mass conservation would be derived from a direct (unweighted) vertical integral of the 3D mass conservation equation (<ref>). Instead, (<ref>) 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 (v̅_x and v̅_y) can be derived from simple weighted vertical averages of their 3D counterparts. Finally, we remind the reader that the `planetary potential', Φ_p, which forces this 2D system should be taken to be the vertical average of the 3D potential specified in (<ref>). This prescription is generally applicable to 2D disc models, and is discussed in more detail in section <ref>.



 §.§ 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's instructive to consider the exact 2D flow which solves (<ref>) (which is well approximated by our averaged flow), and exactly conserves the Bernoulli invariant (which we demonstrate below). Equation (<ref>) implies the existence of a stream function ψ. We take the definition

    Σ̅v̅ = -e_z ×∇ψ.

We may rewrite (<ref>) as

    (2 Ω_p + ∇×v̅) ×v̅ + ∇(1/2|v̅|^2 + W̅ + Φ̅_t + Φ_p) = 0,

where W̅ = P̅/Σ̅. We define the Bernoulli invariant, B, as

    B = 1/2|v̅|^2 + W̅ + Φ̅_t + Φ_p,

and introduce the potential vorticity (PV) (or vortensity)

    ζ = 2 Ω_p + e_z· (∇×v̅)/Σ̅.

We may then rewrite equation (<ref>) as

    Σ̅ζe_z ×v̅ + ∇ B = 0
     - ζe_z ×(e_z ×∇ψ) + ∇ B = 0
    ∇ B = -ζ∇ψ.

That is to say, the Bernoulli function is constant on streamlines, B = B(ψ), and further the PV,

    -  B/ψ = ζ(ψ)

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

    ζ≡Ω_p/2 Σ_p.

Combining this with the definition (<ref>) implies that linearly

    _x v̅_y - _y v̅_x - Ω_p/2 c_γ^2P'/Σ_p = 0.

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

    -Ω_p/2 Σ_pψ = 1/2|v̅|^2 + W̅ + Φ̅_t + Φ_p.

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

    ψ/Σ_p = 3/4Ω_p x^2 - 2/Ω_p(P'/Σ_p+Φ_p) + 3 x v̅_y + 𝒪(q^2/h_γ^6c_γ^2/Ω_p).




 §.§ On the corotation singularity


Importantly, (<ref>) allows us to write down the equation for the radial displacement of fluid elements, ξ̅_x, which we take to satisfy

    v̅·∇ξ̅_x = v̅_x,

with ξ̅_x = 0 far upstream. If a given streamline (or contour of ψ) has radial location x_0 far upstream, we have that

    ξ̅_x = x - x_0.

Comparing values of ψ far upstream to the point (x,y), we see that at leading order

    ψ = 3/4Σ_pΩ_p x_0^2 = 3/4Σ_pΩ_p (ξ̅_x-x)^2 
     = 3/4Σ_pΩ_p x^2 - 2Σ_p/Ω_p(P'/Σ_p+Φ_p) + 3 x Σ_pv̅_y,


    ξ̅_x^2 - 2 x ξ̅_x + 8/3Ω_p^2(P'/Σ_p+Φ_p) - 4 x v̅_y/Ω_p = 0,


    ξ̅_x = x ±√(x^2 -8/3Ω_p^2(P'/Σ_p+Φ_p) + 4 x v̅_y/Ω_p),

where the choice of + or - is determined by whether the fluid element has just undertaken a horseshoe turn or not. Here, P' and 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 <ref>.

Equation (<ref>) 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 <cit.> and <cit.> may be thought to implicitly discard the quadratic term ξ̅_x^2, which is non-negligible for x = 𝒪(√(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 v̅_x, v̅_y and P'). The singularity is introduced in linear analyses when the linearized azimuthal momentum equation, in our case equation (<ref>), is used to directly re-arrange for v̅_y, namely by integrating with respect to azimuthal coordinate y or θ (which may look like dividing by im(Ω - Ω_p) in Fourier space). More specifically, this `first integral' equation for v̅_y must become non-linear to remain valid near to corotation, just as the azimuthal integral of v̅_x becomes large enough that inclusion of the quadratic term ξ̅_x^2 in equation (<ref>) is necessary.

Importantly, this means that any materially conserved quantity Q(ψ) = Q(ξ̅_x - x), will experience the same singularity as the radial displacement, which appears when the quadratic term in (<ref>) is neglected. A gradient of (the materially conserved) potential vorticity 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.



 §.§ 2D mode orthogonality and torque decomposition


For γ < 2, the unforced linearized equations of motion in a disc admit a `2D mode' solution with v'_z = 0. This 2D mode is a member of a wider family of modes, whose axisymmetric members are discussed in <cit.>. 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 γ > 2, a different family of modes exists which does not include such a mode with v'_z = 0. Indeed, imposing v'_z = 0 in this case yields an unnormalisable mode of infinite energy. Appropriately normalised, the 2D mode may be written as


    v'_x,0 = v̅_x(x,y)√(2-γ)exp((γ-1)z^22 H_γ^2),
       v'_y,0 = v̅_y(x,y)√(2-γ)exp((γ-1)z^22 H_γ^2),
       p'_0 = P'(x,y)√(2-γ)ρ_p(z)/Σ_pexp((γ-1)z^22 H_γ^2),
       v'_z,0 = 0,   ρ'_0 = p'_0/c_γ^2,


where again

    v̅_x ≡⟨ v'_x⟩,   v̅_y ≡⟨ v'_y⟩,    P' ≡Σ_p ⟨ p'/ρ_p(z)⟩.

We may then decompose, for example, the radial velocity into its 2D mode component and an orthogonal complement, which we call v'_x,c, describing the rest of the disc's 3D motions, which include for example inertial and gravity waves. 

    v'_x = v'_x,0 + v'_x,c.

Now, v'_x,0 and v'_x,c are everywhere orthogonal with respect to the inner product involving the density

    ⟨ f,g ⟩ = 1/Σ_p∫_-∞^∞ρ_p(z) f(z) g(z)   z.

This may be seen from the relation for the vertical average ⟨ ⋯⟩ in terms of the density-weighted inner product

    ⟨ ⋯⟩≡⟨ ⋯,1/√(γ)exp((γ-1)z^22 H_γ^2)⟩.

Specifically, the orthogonality follows via

    ⟨ v'_x,0, v'_x,c⟩   = ⟨ v'_x,0, v'_x - v'_x,0⟩
       = √(2-γ) v̅_x⟨exp((γ-1)z^22 H_γ^2), v'_x⟩
            - (2-γ)v̅_x^2⟨exp((γ-1)z^22 H_γ^2),exp((γ-1)z^22 H_γ^2)⟩
       = √(γ(2-γ)) v̅_x⟨ v'_x⟩ - (2-γ)v̅_x^2√(γ/2-γ)
       =0.

As a result, we obtain the Parseval identity,

    ⟨ v'_x, v'_y ⟩ = ⟨ v'_x,0, v'_y,0⟩ + ⟨ v'_x,c, v'_y,c⟩,

as well as equivalent identities for any two variables drawn from the set {v'_x, v'_y, p'/ρ_p(z)}. In particular, this allows us to decompose the radial angular momentum flux, F_A, into independent components associated with the 2D mode and 3D remainder, since

    F_A ≈ r_p∬ρ_p(z) v'_x v'_y   z  y 
     = r_p∬ρ_p(z) v'_x,0 v'_y,0  z  y + r_p∬ρ_p(z) v'_x,c v'_y,c  z   y 
     = √(γ(2 - γ))r_p Σ_p∫_-∞^∞v̅_xv̅_y   y + F_A^3D≡ F_A^2D + F_A^3D

It's in this sense that the torque on the planet for a general adiabatic index γ may be separated into 2D and 3D contributions. The factor of √(γ(2-γ)) (approximately 92% for γ = 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 ≪ 1/h_γ (see for example equation (45) in <cit.> and appendix B1 of <cit.>). 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 √(γ(2-γ)) from the flux expression. This follows from the applicability of resonant torque excitation theory to this regime.

This principle is at play in figure 2 of <cit.>. 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≳2H_γ, and for each value of γ considered, <cit.> recover torque densities matching the case γ = 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 <cit.>, so that their γ-dependent scaling of the x-axis with H_γ and y-axis with c_γ^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 γ > 1. That is, the torque density for x≳2H_γ matches that obtained via a naïve approximation of the 3D flow by its vertical average, as defined in (<ref>).

There is however in the planet-disc interaction problem an important (indeed, a dominant) contribution from azimuthal modes with m ∼ 1/h_γ. 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_A^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 <cit.>. 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.



§ REDUCTION TO INDEPENDENT SECOND ORDER LINEAR EQUATIONS



In section <ref>, 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


    D P' + Σ_p c_γ^2(_x v̅_x + _y v̅_y) = 0,
       _x v̅_y - _y v̅_x - Ω_p/2 c_γ^2P'/Σ_p = 0,
       Dv̅_x - 2Ω_p v̅_y + 1/Σ_p_x P' = -_x Φ_p,
       Dv̅_y + 1/2Ω_p v̅_x + 1/Σ_p_y P' = -_y Φ_p,


where D ≡ -3/2xΩ_p_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 which describe both the spiral wave excitation and horseshoe streamlines. It's helpful to non-dimensionalise the equations. We first replace

    x → H_γ x,    y → H_γ y,

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

    𝒟 = -3/2x_y,

and introduce the non-dimensional potential ϕ̂_p

    Φ_p(x,y) = q/h_γ^3c_γ^2ϕ̂_p(s),    s = √(x^2+y^2),
    ϕ̂_p(s) = -^1/4s^2/√(2)K_0(14s^2) ∼ -1/s + 𝒪(1/s^3).

We define non-dimensional, scaled x- and y- velocity and `enthalpy' perturbations u(x,y), v(x,y) and W(x,y) via

    v̅_x ≡q/h_γ^3c_γ u,   v̅_y ≡q/h_γ^3c_γ v,    P' ≡q/h_γ^3c_γ^2Σ_p W.

Equations (<ref>), (<ref>) and (<ref>) become


    _x v - _y u - 1/2χ = - 1/2ϕ̂_p,
       𝒟 u - 2v + _x χ = 0,
       𝒟 v + 1/2u + _y χ = 0,


where χ≡ W + ϕ̂_p. Additionally, `mass conservation', (<ref>) becomes

    𝒟χ + _x u + _y v = 𝒟ϕ̂_p,

which may also be derived from (<ref>), (<ref>) and (<ref>). It's clear at this stage that the flow is indeed linear so long as q ≪ h_γ^3.

Remarkably, these equations are independent of γ! They are analogous to the equations considered by <cit.>; however, it's 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 𝒟 = 0 (before then differentiating with respect to x). The resulting system omits large factors at corotation in the fluid case. Note that _x𝒟→ -32_y as x → 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,

    𝒟(_x v - _y u) + 1/2(_x u + _y v) = 0.

To proceed, we take 𝒟(<ref>) and combine with (<ref>), and then take 𝒟(<ref>) and combine with (<ref>). This yields


    (𝒟^2 + 1)u = - 𝒟_x χ - 2 _y χ,
       (𝒟^2 + 1)v = - 𝒟_y χ + 1/2_x χ.


Now, it may be shown from equations (<ref>) and (<ref>) that

    𝒟_x χ + 2 _y χ = 𝒟_x ϕ̂_p + 2 _y ϕ̂_p - ∇^2 u + 3_y(χ - ϕ̂_p),

and

    𝒟_y χ - 1/2_x χ = 𝒟_y ϕ̂_p - 1/2_x ϕ̂_p - ∇^2 v.

Equations (<ref>) and (<ref>) become


    (𝒟^2 + 1 - ∇^2)u + 3 _y χ = _y ϕ̂_p - 𝒟_x ϕ̂_p,
       (𝒟^2 + 1 - ∇^2)v = - 𝒟_y ϕ̂_p + 1/2_x ϕ̂_p.


Now, to find an equivalent equation for χ, we take 𝒟^2(<ref>). We're then able to use equation (<ref>), the momentum equations (<ref>) & (<ref>) and the PV equation (<ref>) again to deduce

    (𝒟^2 + 1 - ∇^2)χ + 3_y u = (𝒟^2 + 1)ϕ̂_p.

We now define the linearized Riemann invariants[The interpretation of J_± as linearized Riemann invariants is clarified in section <ref>.] J_±:

    J_±≡ u ±χ.

Combining equations (<ref>) and (<ref>), it follows that


    [𝒟^2 + 1 ± 3 _y - ∇^2]J_± =  [_y - 𝒟_x ± (𝒟^2 + 1)]ϕ̂_p,
       [𝒟^2 + 1 - ∇^2]v =  - 𝒟_y ϕ̂_p + 1/2_x ϕ̂_p.


Now, equation (<ref>) is the parabolic cylinder equation[More specifically, it becomes the parabolic cylinder equation following an azimuthal mode decomposition or Fourier transform.] derived by <cit.> for the azimuthal velocity perturbation, but equation (<ref>) 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 <cit.>, 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 <cit.>, who studied the excitation of density waves via turbulence. In the [appx]appendix, we discuss the numerical solution of equations (<ref>) and (<ref>) to high accuracy, and we present these solutions in the next section.

The inner limit of equations (<ref>) and (<ref>) (that is, the simplified equations valid close to corotation taking x ≪ H) may be written as


    [1 ± 3 _y - ∇^2]J_± =  [_y - 𝒟_x ± (𝒟^2 + 1)]ϕ̂_p,
       [1 - ∇^2]v =  - 𝒟_y ϕ̂_p + 1/2_x ϕ̂_p.


Note that terms including 𝒟 acting on ϕ̂_p may not be neglected, as in this case ϕ̂_p has a (logarithmic) singularity at the origin, with the effect that such terms are non-negligible. For comparison, the equivalent equation from <cit.> expressed in our notation reads

    [1 - ∂_x^2 - 4∂_y^2]χ = ϕ_p.




 §.§ Connection with resonant wave excitation theory


One may deduce directly from (<ref>) 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, κ(r), is an integer multiple of the interaction rate, Ω - Ω_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 ω^2 = κ^2 + c^2 |k|^2 for wavevector 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 = ±23H_γ.

Furthermore, in the gas-dynamic case the contribution of the radial wavenumber k_x to the wave frequency ω 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's 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 𝒟^2 + 1 - ∇^2 is hyperbolic in the wave-permitting region |x| > 23H_γ, and elliptic in the region where all modes are evanescent, namely |x| < 23H_γ. This signposts that the Lindblad resonances must all stand off a distance 23H_γ from the planet's orbital radius. More specifically, for the azimuthal disc mode proportional to ^ m θ, recalling the original definition y = r_pθ and reintroducing dimensions, (<ref>) becomes

    (- 9/4Ω_p^2 x^2/r_p^2m^2 + Ω_p^2 + c_γ^2m^2/r_p^2 - c_γ^2_x^2 )v = ⋯.

We see that were it not for the acoustic terms in the above operator, namely c_γ^2 m^2/r_p^2 and c_γ^2_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,± m = r_p(1 ±23 m). 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 |x| < 23H_γ√(1 + (h_γ m)^-2), and wavelike for |x| > 23H_γ√(1 + (h_γ m)^-2). The effective Lindblad resonances close to the planet (with h_γ m ≳ 1) therefore have locations

    r_L,± m = r_p ±2/3H_γ√(1 + (h_γ m)^-2).




§ RESULTS



In this section, numerical solutions to equations (<ref>) and (<ref>) are presented and discussed. We compare our flow solution near to corotation with flows computed in 3D simulations, for example by <cit.>, <cit.> and <cit.>. 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.




 §.§ Flow streamlines and the horseshoe width


Real-space plots for the solutions of equations (<ref>) and (<ref>) for each flow variable (u, v, χ, J_+, J_-, as well as the stream function ψ), are shown in figure <ref>. 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 × 10^-5. Details on this numerical procedure are discussed in the [appx]appendix.

Qualitatively, the vertically averaged flow field streamlines depicted in the bottom right panel of figure <ref> are comparable to those obtained with a `softened' potential with smoothing length b = 0.4H_γ. 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 <cit.>.

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% 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%, depicted in figure <ref>. 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 <cit.>. Indeed, in figure 5 of <cit.> 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 semi-thickness x_s may be found from our solution by first noting that χ, u and v decay rapidly to 0 as y →∞ in the region |x| < 2/3H_γ. The leading order far-field horseshoe streamlines therefore have reflectional symmetry about x = 0. We showed in section <ref> that the stream function in our vertically averaged flow may be expressed as

    ψ/Σ_p = 3/4Ω_p x^2 - 2/Ω_p(P'/Σ_p+Φ_p) + 3 x v̅_y + 𝒪(q^2/h_γ^6c_γ^2/Ω_p).

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

    3/4Ω_p x_s^2 = -2/Ω_pq/h_γ^3c_γ^2 χ_s,

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

    χ_s = -0.47115,

and in the same limit q/h_γ^3 → 0, the flow has three stagnation points, with coordinates x = 0, y = ± 0.439 H_γ and y = 0. We therefore predict a horseshoe region half-width for a low-mass planet of

    x_s = 1.12089√(q/h_γ^3)H_γ = 1.12089√(q/h^3)Hγ^-1/4,

in good agreement with 2D simulations using a smoothing length b = 0.4H_γ <cit.>, which match well the horseshoe width measured in 3D simulations for this choice of b <cit.>.

The exact γ^-1/4 dependence is a non-trivial result, arising from the coincidence that the 2D mode of the potential, specified in equation (<ref>), depends only on the length-scale H_γ 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 for example figures 3 and 4 in <cit.>, 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 γ = 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, for example those performed by <cit.> and <cit.>.

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 <cit.>. Instead, to gain traction here we appeal to the vertical structure of the 2D mode. As discussed in section <ref>, this mode behaves as

    v'_x ∝exp((γ-1)z^22 H_γ^2),    v'_y ∝exp((γ-1)z^22 H_γ^2),    v'_z = 0,
    
        p' ∝exp(-z^22 H_γ^2),   ρ' ∝exp(-z^22 H_γ^2).

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

    x_s(z) ≈ 1.12 √(q/h^3)H γ^-1/4√(2-γ)exp((γ-1)z^22 H_γ^2).

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 √(x) near to corotation <cit.>. Gravity waves are also confined near to corotation to the region above diagonal buoyancy resonances <cit.>. 

Whilst the exponential increase in horseshoe width (as well as perturbed flow velocities) with height is striking, we note that that the factor (γ - 1)/2γ is small for reasonable values of γ. Further, we shouldn't be too concerned with the dynamics beyond 2 to 3 scale heights above the mid-plane, as over 95% of the disc's mass is contained within the first 2 scale heights. Moreover, our model's assumptions, including linearity, and that the background state is vertically isothermal (discussed in section <ref>) may begin to break down at such heights.

Figure <ref> depicts this vertical dependence for a few values of γ. Curiously, the expression in (<ref>) is poorly defined for γ⩾ 2; however we don't expect γ to take such values in astrophysical discs. This oddity is due to the change in the character of the the disc's linear modes as γ increases through 2. For γ⩾ 2, for example there is no mode with v'_z = 0 (such a mode would contain infinite energy). 


In the isothermal case, our numerical value for the horseshoe width x_s = 1.12√(q/h^3)H is in reasonable agreement with 3D simulations (though inflated by a few percent). <cit.> and <cit.> both estimate a width via 3D simulation of x_s = 1.05√(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 (<ref>), 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 (for example 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 order h_γ∼ 0.05. Furthermore, we exclude any non-linearity in our calculation. From equation (<ref>) we expect a next order correction to x_s of relative size 𝒪(q/h^3), which may also be on the order of a few percent. 



 §.§ 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 𝒪(√(q/h^3)H), 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

    v_y,hs≈ -3/2Ω_px + v'_y(0,y),    v_x,hs≈ v'_x(0,y)

This was noted by <cit.>, who studied the saturation of the corotation torque, and used matched asymptotics to compute Fourier coefficients for the flow velocity perturbations near to corotation. Figure <ref> depicts the non-dimensional radial and azimuthal velocity distributions at corotation, which satisfy (<ref>) and (<ref>). These correspond to taking cross-sections of the flow depicted in figure <ref> along the line x=0.




 §.§ 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 figure <ref> and corresponding one-sided Lindblad torque.

The non-linear 2D theory developed by <cit.> 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

    R_± = u ±2/γ - 1(c - c_γ),    c = c_γ(Σ/Σ_p)^γ - 1/2.

For small departures from the background state, the enthalpy perturbation obeys W'/c_γ≈ c - c_γ. We may write therefore

    R_±≈ u ±W'/c_γ,

equal to J_± up to an unimportant additive contribution of Φ_p/c_γ. It therefore makes sense to interpret J_± as linearized Riemann invariants.

A WKB analysis, appropriately imposing outgoing wave boundary conditions, indicates that for each azimuthal Fourier mode (with wavenumber k_y > 0),

    J̃_+∼
    |x|^+12exp(+34Ω_p k_yc_γx^2),   if   x > 0
    
    |x|^-32exp(-34Ω_p k_yc_γx^2),   if   x < 0.

The behaviour of J̃_- can be determined from the symmetry J_+(x,y) = -J_-(-x,-y), or in Fourier-space, J̃_-(x,k_y) = -J̃_+^*(-x,k_y). Note the strong |x|^-3/2 decay of J̃_+ in x < 0. This may be thought of as a consequence of J_+ being approximately conserved on characteristics emanating from x = -∞, where it is 0 (so that we have a simple wave). The consequence of this can be seen in figure <ref>, 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 (<ref>) and (<ref>). That is, the density waves in x > 0 follow the curve η(x,y) = 0, for characteristic coordinate

    η = y - x/3√(9/4x^2 - 1) - 1/3cosh^-1(3x/2).


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're also needed to compute the locations of eventual shocking due to wave steepening. Whilst 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 Φ_p, shown in figure <ref>, are qualitatively very similar to the profiles obtained in planet-driven wave evolution studies <cit.>.

Our one-sided torque estimate is however more than a factor of 2 smaller than that obtained in the 2D studies conducted by <cit.> and <cit.>. 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 <cit.>, and especially for smaller softening parameters. Indeed, we find that the angular momentum flux carried by the 2D mode, F_A^2D = √(γ(2-γ))r_pΣ_p∫_-∞^∞v̅_xv̅_y  y (as defined in section <ref>), far from the planet is

    F_A^2D(∞) = 0.37 √(γ(2-γ))(G M_p)^2 Σ_p r_p Ω_p/c_γ^3,

corresponding to a one-sided torque

    T^2D≡ F_A^2D|_0^∞ = 0.34 √(γ(2-γ))(G M_p)^2 Σ_p r_p Ω_p/c_γ^3.

For comparison, in the case of a point-mass potential in a 2D disc, <cit.> calculate a dimensionless numerical torque prefactor of 0.93. We've used here that F_A^2D(0) = 0.03 √(γ(2-γ))(G M_p)^2 Σ_p r_p Ω_p/c_γ^3, which is apparent from figure <ref>. The fact that F_A^2D(0) ≠ 0 is due to a small fraction (4%) of the flux carried by the outward-propagating wave being excited in x < 0.

In the case γ = 1, the torque estimate (<ref>) is in agreement with the 3D (isothermal) simulations of <cit.> to within a few percent (one can find a one-sided Lindblad torque estimate from the symmetric part of the cumulative torque profiles in their figures 7 and 8); it further matches precisely the torque found by <cit.>. This is because in the case γ = 1, all gravity waves are suppressed, and inertial waves carry very little flux <cit.>.

<cit.> and <cit.> 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 figure 2 of <cit.>, it appears that in the case γ = 1.4, the 2D mode contributes only 75% 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_A^2D(x), is plotted in figure <ref>. 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's a very small decrease in the flux beyond x ≈ 3.5 H_γ, which may be attributed to the `negative torque phenomenon' <cit.>.


Finally, we're 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'_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 <cit.>, which suggests considering the wave steepening of an isolated 2D mode would provide a good model for 3D shock formation.

Using equation (<ref>) and comparing the maximum steepness of our 2D mode profile with that of <cit.> and following their wave steepening calculation, we estimate the location of the initial shock in a 3D disc relative to their 2D estimate as

    l_sh^3D≈ 1.2^-γ - 15 H_γ^2z^2/(2-γ)^1/5· l_sh^2D.

Here, we've 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% of that found by <cit.> at x = 1.33 H_γ). The vertical dependence is weak for physical values of γ; for γ = 1.4 it varies by only 5% within the first scale height, H_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 <ref>, our model's basic assumptions may begin to break down.



§ DISCUSSION



We've 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 <ref>, 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 <ref>);
    
  * the appearance of and displacements in the locations of the Lindblad resonances within the framework of this paper in section <ref>.

However, it's worth discussing caveats to this work as well as two further topics which 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. 



 §.§ Treating planetary potentials in 2D discs




  §.§.§ 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

    Φ_p = -G M_p/√(|r - r_p|^2 + b^2),    b ∼ H/2,

are commonly adopted, as well as other similar prescriptions including the `fourth order' softened potential described in equation (16) of <cit.>. Each of these prescriptions however struggle to uniformly capture the impact of the disc's vertical extent on the dynamics within a few scale heights of the planet.

<cit.> compared (<ref>) 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 |r - r_p| ∼ 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 (<ref>) 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 (<cit.>, figure 5).



  §.§.§ A rigorous choice


In light of this, in this section we emphasise that (<ref>) is a far better choice for the 2D potential. In particular, taking inspiration from <cit.>, we proved in section <ref> 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' which are forced precisely by (<ref>). In this way, we argue that (<ref>) (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.

    Φ_p ≡ -G M_p/H_γ^1/4s^2/√(2)K_0(14s^2),    s = |r-r_p|/H_γ.

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

Indeed, (<ref>) is uniformly a very good approximation to the potential which forces the global 2D mode of the disc (which would simply involve s = |r-r_p|/H_γ(r)), since far away from the planet

    Φ_p ∼ -G M_p/|r-r_p|[1 - 1/2s^2 + 9/8s^4 + 𝒪(s^-6)].

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 <cit.> 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 Φ_p is also small.


In figure <ref> we compare graphically the appearance of the vertically averaged potential (<ref>) with a point-mass potential and `softened' potentials of the form (<ref>), as well as their corresponding gradients. Note that Φ_p has a logarithmic singularity at the origin. More precisely,

    Φ_p ∼G M_p/H_γ2 ln s/√(2)  as  s → 0.

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_p. We define the small parameter

    ε = R_p/H_γ.

For Earth-like planets located at around 1 to 10 au, we expect 10^-5≲ε≲ 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

    ϕ_p = 
            -G M_p/|r-r_p|,   |r-r_p| > R_p,
    
            -G M_p/R_p,   |r-r_p| ⩽ R_p.

In this case,

    Φ_p(0) ∼G M_p/H_γ2lnε/√(2).

However, |2lnε/√(2)| ∼ 7 < h_γ^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 χ shown in section <ref> 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 timestep issues), to `soften' the 2D potential (<ref>) by taking s^2 = ε^2 + |r-r_p|^2/H_γ^2, though here a practical choice for ε for grid-based simulations in which the planetary radius is unresolved would be on the order of the grid spacing divided by H_γ.

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 ^xK_0(x) and ^xK_1(x) (which both appear in the expression for the gradient of Φ_p) already exist within the Cephes library[<https://www.netlib.org/cephes/>] <cit.>. Alternatively, the vertically averaged potential in (<ref>) and its derivative are functions of one parameter only, so that a lookup table would provide an efficient method for evaluation.



 §.§ 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_γ 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.



  §.§.§ 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. In the absence of strong viscosity or dissipation therefore, 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's 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 [s:intro]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 figure <ref>. 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.



 §.§ Caveats


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



  §.§.§ 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% for γ = 1.4) of the one-sided torque on the planet, as discussed in section <ref>.



  §.§.§ 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, for example turbulence, viscosity, or a vertical magnetic field. Favourable conditions include a large plasma β≫ 1 and low level of ionisation, as well as a low Reynolds stress, captured by the condition for the viscosity parameter α≪ 1 <cit.>. It's widely believed that except for the innermost regions, protoplanetary discs do predominantly exhibit β≫ 1 <cit.>, and recent numerical and observational studies point towards a smaller turbulent viscosity than previously thought. Non-ideal MHD effects such as ambipolar diffusion and Ohmic resistivity act to suppress turbulence even in the disc's surface layers <cit.>, and to decouple the fluid from the magnetic field. Indeed, <cit.> and later <cit.> find that observations of the discs around HL Tau and Oph163131 are consistent with small viscosity parameters α∼ 10^-4 and α∼ 10^-5 respectively. 



  §.§.§ 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 won't 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 which 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 <cit.>, <cit.>, <cit.> and <cit.>, infrequent gas-dust collisions often act as a bottleneck for the thermal relaxation of the gas, especially in its surface layers. <cit.> predict that typical gas thermal relaxation time-scales are tens of orbits even at 10 au, corresponding to a cooling parameter β≡ 2 t_relax/t_orb≳ 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 <cit.>. The existence of such a signature adds to the credibility of the adiabatic model, as it implies t_relax≫ N_z^-1 (for N_z the Brunt-Väisälä frequency), necessary for the buoyancy wake to develop. It's likely however that cooling and thermal diffusion play important roles on the much longer time-scale associated with the librating horseshoe motion. Whilst this has only a small effect on the dominant flow behaviour found in this paper, it has important consequences for the entropy and potential vorticity 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 which 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 <cit.>. 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. Whilst there's no exact analogue for the 2D projection procedure described in section <ref> 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. Whilst an oversimplification, this suggests therefore that the vertically isothermal model for the background disc may be expected to give quantitatively reasonable results.



§ CONCLUSIONS



For linear adiabatic perturbations to a vertically isothermal protoplanetary disc, there is a particular vertical average of the 3D gas dynamic equations which 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 <cit.>. Importantly, when compared with more traditional softening prescriptions for planetary potentials, adopting the averaged potential (<ref>) 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, (<ref>) and (<ref>). 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 <cit.>.

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√(q/h^3)Hγ^-1/4. Taking only the 2D mode contribution to the flow predicts a horseshoe width which grows with height above the disc mid-plane as described in equation (<ref>) and figure <ref>. 

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 figure <ref>, 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 under-studied, yet important problem which requires future attention, not only because of the angular momentum which they transport, but also for their observational signature.



§ ACKNOWLEDGEMENTS


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



§ DATA AVAILABILITY


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

mnras




§ APPENDIX A: NUMERICAL PROCEDURE





In this section we outline how the numerical solutions plotted in figure <ref> were obtained. We sought solutions with very high resolution and accuracy so that in future work we may use them to predict accurately the intricate dynamics which inform the corotation (and indeed wave) torques. We consider the problem in the (x,k_y) plane, having Fourier transformed with respect to y. For brevity we denote k_y = k. It is helpful to define

    x' = √(3|k|)x,    a = 1 + k^2/3 |k|.

Equations (<ref>) and (<ref>) become the canonical parabolic cylinder equations,


    [_x'^2 + 14x'^2 - (a ± sgn(k))]J̃_±
                = - sgn(k)[13 + 12 x'_x'∓(14x'^2 - 13 k)]ϕ̃_p.
    
       [_x'^2 + 14x'^2 - a]ṽ = 1/2√(3)[x'√(|k|) - 1√(|k|)_x']ϕ̃_p,


forced by the Fourier transform of the potential ϕ̂_p, which we denote ϕ̃_p. It may be shown (namely by evaluating the Fourier transform prior to the vertical averaging) that

    ϕ̃_p(x,k) = -√(2/)∫_-∞^∞K_0(|k|√(x^2+z^2))^-z^2/2 z,

which is far cheaper to evaluate accurately than the direct expression for the transform of ϕ̂_p.

Now, as x' → +∞, the solution will consist of the homogeneous waves U( a',x'^-/4) and U(- a',x'^/4). Here a' is understood to denote the order of the relevant parabolic cylinder equation, that is, either a± or a, and we adopt U(a,x) as defined in chapter 19 of <cit.>. We use the parabolic cylinder function (PCF) U, as E(a',x) and W(a',x) are not defined in general for a' ∈ℂ. From now on we take k ⩾ 0, noting our real real-space variables will have conjugate-symmetric Fourier transforms. With this in mind, and noting that ∀δ > 0, the PCF U(a,z) ∼ z^-a-1/2^-1/4z^2 as z →∞ in (z) ⩽34 - δ, we identify the following solutions as in- or outgoing waves:


    U( a',x'^-/4)    is an outgoing wave  as  x' →∞
       U(- a',x'^/4)    is an ingoing wave as  x' →∞
       U( a',-x'^-/4)    is an ingoing wave as  x' → -∞
       U(- a',-x'^/4)    is an outgoing wave as  x' → -∞


We want no incoming waves present in the solution. That is, our desired particular integral solution for J̃_± or ṽ, (which we denote η̃_p) behaves as


    η̃_p ∼ b_1U( a',x'^-/4)   as  x' →∞
       η̃_p ∼ b_2U(- a',-x'^/4)   as  x' → -∞.


We eliminate incoming waves in our numerical solution with the following algorithm. Suppose we compute a numerical solution η̃, imposing arbitrary initial data at x' = 0. It follows that, without loss of generality


    η̃∼ (b_1 + c_1) U( a',x'^-/4) 
                    + c_2 U(- a',x'^/4)   as  x' →∞,
    
       η̃∼ (b_2 + d_1)U(- a',-x'^/4) 
                    + d_2U( a',-x'^-/4)   as  x' → -∞.


 We may then fit the numerical solutions via linear regression to the homogeneous PCF wave solutions far from x' = 0 (where the forcing has decayed sufficiently) to extract values for c_2 and d_2. We note that c_2 and d_2 both depend linearly on the (arbitrary) initial data imposed at x' = 0. Integrating the system twice more with two different choices for initial data allows us to solve this linear system for the correct initial data (which yields c_2 = d_2 = 0). We then perform one final integration with this choice of initial data, giving a numerical solution for η̃_p. 

 We carried out this algorithm with a 4th order Runge-Kutta ODE solver to compute the solutions to equations (<ref>). This method provides a very high numerical accuracy, and we were able to achieve a relative error of no more than 1× 10^-5 for all 0.01 ⩽ k ⩽ 500.

It's worth noting some explicit analysis concerning the k=0 case, which is not permitted in the previous algorithm, and further how we can treat the singularities in the system in both real and Fourier-space without adopting softening prescriptions or cutoffs. In the case k=0, ϕ̃_p doesn't exist, and indeed J̃_± diverges too. However, since

    [1-_x^2]W̃|_k=0 = _x^2ϕ̃_p,   and  [1-_x^2]ṽ|_k = 0 = 1/2_xϕ̃_p,

for W̃ the transformed enthalpy perturbation, we may write

    W̃(x,0) = ∫_-∞^∞^-|x'-x|[√(/2)|x'|erfcx(|x'|/√(2)) - 1] x',


    ṽ(x,0) = 1/2∫_-∞^∞^-|x'-x|√(/2)erfcx(|x'|/√(2))sgn(x') x',


    ũ(x,0) = 0,

where erfcx(x) is the scaled complementary error function. Now, as well as the divergence of the potential as k → 0, we have that for large k,

    ϕ̃_p(x,k) ∼ - √(2 )^-|k x|/|k|,

that is, as we should expect, the inverse Fourier transform of ϕ̃_p diverges at x = y = 0 (as the 2D real-space potential has a logarithmic singularity there).

We may however eliminate the need to evaluate the solutions at very many large values of k (in order to obtain high accuracy), and overcome the singular behaviour at k = 0, by defining variables which are well-behaved everywhere in Fourier and real space. The inverse Fourier transform of these variables will then quickly give accurate results. First note that

    ∫_-∞^∞ln(x^2+y^2)^- k y y = -2 /|k|^-|k x|.

It therefore makes sense to define regularised J_± variables

    J_±,reg = J_±∓[ϕ̂_p - 1/√(2)ln(x^2+y^2/1+y^2)]

so that

    J̃_±,reg = J̃_±∓[ϕ̃_p + √(2 )(^-|k x| - ^-|k|/|k|)].

These variables are then indeed well-behaved (with no singularities) everywhere in Fourier and real space. In particular, for k = 0, 

    J̃_±,reg(x,0) = ∓[W̃(x,0) + √(2 )(1 - |x|)],

where W̃(x,0) may be found by evaluating (<ref>). With these regularised variables, we are able to quickly achieve the earlier quoted absolute uncertainty in the flow solution of 1 × 10^-5.




