HORSESHOE DRAG IN THREE-DIMENSIONAL GLOBALLY ISOTHERMAL DISKS

We study the tidal interaction between a gaseous globally isothermal (hence barotropic) protoplanetary disk in rotational equilibrium around a star of mass ${M}_{\star }$ and an embedded protoplanet of mass ${M}_{p}={{qM}}_{\star }$, on a fixed circular orbit in the disk's midplane, with orbital radius r_p and orbital frequency ${{\rm{\Omega }}}_{p}$. We give the position of a given fluid element either in spherical coordinates $(r,\theta ,\phi )$ in a frame corotating with the planet or in cylindrical coordinates, $R=r\mathrm{sin}\theta$ being the cylindrical radius and $z=r\mathrm{cos}\theta$ being the altitude. In each case the frame is centered on the star and such that the planet's ϕ-coordinate equals zero. The angular speed of a fluid element in a non-rotating frame is denoted ${\rm{\Omega }}(r,\theta ,\phi )$. We denote a quantity in the unperturbed disk (i.e., prior to the insertion of the planet) by a subscript 0, and a quantity evaluated at $r={r}_{p}$ by a subscript p. The surface density is defined as

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={\displaystyle \int }_{-\infty }^{+\infty }\rho (R,z){dz},\end{eqnarray} \tag{ 1 }$$

where ρ is the volumic density. In the unperturbed disk, its midplane value follows a power law:

$$\begin{eqnarray}&&\rho (r,\theta =\pi /2,\phi )\propto {r}^{-\xi }.\end{eqnarray} \tag{ 2 }$$

The (uniform) sound speed is denoted c_s, the pressure scale height is $H(r)={c}_{s}\sqrt{{r}^{3}/{{GM}}_{\star }}$, and the disk aspect ratio is $h(r)=H(r)/r$. Since the disk is globally isothermal, the aspect ratio increases as the square root of the radius: $h(r)={h}_{p}{(r/{r}_{p})}^{1/2}$. One can work out the exact dependence of the volumic density and rotational velocity on r and θ, for the case in which the temperature is a power law of the spherical radius, as described in Appendix A. From these dependences, one can infer that the surface density also obeys a power law of radius: ${\rm{\Sigma }}\propto {r}^{-\alpha }$, where α is given by Equation (88).

We introduce the vorticity of the flow ${\boldsymbol{\zeta }}$ viewed in an inertial frame:

$$\begin{eqnarray}&&{\boldsymbol{\zeta }}={\boldsymbol{\nabla }}\times {\boldsymbol{v}}+2{{\boldsymbol{\Omega }}}_{p},\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{v}}$ is the linear velocity in the frame corotating with the planet. Throughout this work, we denote the PV by

$$\begin{eqnarray}&&{\boldsymbol{w}}=\displaystyle \frac{{\boldsymbol{\zeta }}}{\rho }.\end{eqnarray} \tag{ 4 }$$

In order to use a vocabulary consistent with prior work in two dimensions, we reserve the name vortensity for an adequate vertical integral of the PV, which we shall introduce later.

We define the corotation sheet as the two-dimensional region where ${{\rm{\Omega }}}_{0}(r,\theta ,\phi )={{\rm{\Omega }}}_{p}$. This sheet intersects the midplane at $r={r}_{c}(\phi )$.

We focus on the horseshoe region, defined by the particles that cross over from the inner to the outer disk, or vice versa, after a close encounter with the planet. We call the boundary between the horseshoe region and the rest of the disk the separatrix sheet. The streamlines belonging to this sheet are called either the separatrix streamlines, critical streamlines, or widest horseshoes (Fung et al. 2015).

At large azimuthal elongation from the planet, the horseshoe region can be divided into four regions. The front (rear) part corresponds to fluid elements with $\phi \gt 0$ ($\phi \lt 0$). Within the front and rear parts, we distinguish the upstream and downstream regions. The upstream region is the set of fluid elements that have not yet experienced a close encounter with the planet, whereas the downstream region is the set of fluid elements that have already experienced a close encounter, and therefore crossed the corotation by definition of the horseshoe region. In the front part, the upstream region corresponds to $r\gt {r}_{c}$, and the downstream region to $r\lt {r}_{c}$. Opposite relations hold in the rear part.

The equations that govern the evolution of the flow are the continuity equation and the Euler equations, which read respectively, in cylindrical coordinates and in the frame corotating with the planet,

$$\begin{eqnarray}&&{\partial }_{t}\rho +\displaystyle \frac{1}{R}{\partial }_{R}(R\rho {v}_{R})+\displaystyle \frac{1}{R}{\partial }_{\phi }(\rho {v}_{\phi })+{\partial }_{z}(\rho {v}_{z})=0,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\partial }_{t}{v}_{R}+{v}_{R}{\partial }_{R}{v}_{R}+\displaystyle \frac{{v}_{\phi }}{R}{\partial }_{\phi }{v}_{R}+{v}_{z}{\partial }_{z}{v}_{R}-R{{\rm{\Omega }}}_{p}^{2}\\ &&\quad -2{{\rm{\Omega }}}_{p}{v}_{\phi }-\displaystyle \frac{{v}_{\phi }^{2}}{R}=-\displaystyle \frac{{\partial }_{R}P}{\rho }-{\partial }_{R}{\rm{\Phi }},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\partial }_{t}j+{v}_{R}{\partial }_{R}j+\displaystyle \frac{{v}_{\phi }}{R}{\partial }_{\phi }j+{v}_{z}{\partial }_{z}j=-\displaystyle \frac{{\partial }_{\phi }P}{\rho }-{\partial }_{\phi }{\rm{\Phi }},\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{\partial }_{t}{v}_{z}+{v}_{R}{\partial }_{R}{v}_{z}+\displaystyle \frac{{v}_{\phi }}{R}{\partial }_{\phi }{v}_{z}+{v}_{z}{\partial }_{z}{v}_{z}=-\displaystyle \frac{{\partial }_{z}P}{\rho }-{\partial }_{z}{\rm{\Phi }},\end{eqnarray} \tag{ 8 }$$

where v_R, v_ϕ, and v_z are the cylindrical components of the velocity in the rotating frame, $P=\rho {c}_{s}^{2}$ is the pressure, $j={R}^{2}{\rm{\Omega }}={R}^{2}{{\rm{\Omega }}}_{p}+{{Rv}}_{\phi }$ is the specific angular momentum, evaluated in the non-rotating frame, and Φ is the gravitational potential given by

$$\begin{eqnarray}{\rm{\Phi }} & = & -\displaystyle \frac{{{GM}}_{\star }}{\sqrt{{R}^{2}+{z}^{2}}}-\displaystyle \frac{{{GM}}_{p}}{\sqrt{{R}^{2}+{r}_{p}^{2}-2{{Rr}}_{p}\mathrm{cos}\phi +{z}^{2}}}\\ & & +\displaystyle \frac{{{GM}}_{p}}{{r}_{p}^{2}}R\mathrm{cos}\phi ,\end{eqnarray} \tag{ 9 }$$

where the different terms are respectively the star's potential, the planet's potential, and the indirect term arising from the acceleration of our non-inertial frame, centered on the star.

The right-hand side of Equations (6)–(8) can be recast as the gradient of the effective potential

$$\begin{eqnarray}&&\tilde{{\rm{\Phi }}}={\rm{\Phi }}+\eta ,\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\eta ={c}_{s}^{2}\mathrm{log}\left(\displaystyle \frac{\rho }{{\rho }_{00}}\right),\end{eqnarray} \tag{ 11 }$$

is the fluid enthalpy, and ${\rho }_{00}$ is an arbitrary constant dimensionally homogeneous to a density. Like the gravitational potential, the enthalpy is defined to within an additive constant, and specifying ${\rho }_{00}$ amounts to choosing for which value of the density the enthalpy vanishes. In what follows we set ${\rho }_{00}={M}_{*}{r}_{p}^{-3}$. Multiplying Equations (6)–(8) respectively by v_R, ${v}_{\phi }/R$, and v_z, and summing the results under the assumption of a steady flow, we are left with

$$\begin{eqnarray}&&{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}({E}_{{\rm{kin}}}+\tilde{{\rm{\Phi }}})-{{Rv}}_{R}{{\rm{\Omega }}}_{p}^{2}=0,\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}\equiv {v}_{R}{\partial }_{R}+\displaystyle \frac{{v}_{\phi }}{R}{\partial }_{\phi }+{v}_{z}{\partial }_{z}$ and

$$\begin{eqnarray}&&{E}_{{\rm{kin}}}=\displaystyle \frac{1}{2}({v}_{R}^{2}+{v}_{\phi }^{2}+{v}_{z}^{2}),\end{eqnarray} \tag{ 13 }$$

hence we have

$$\begin{eqnarray}&&{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{B}_{J}=0,\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{B}_{J}={E}_{{\rm{kin}}}+\tilde{{\rm{\Phi }}}-\displaystyle \frac{1}{2}{R}^{2}{{\rm{\Omega }}}_{p}^{2}\end{eqnarray} \tag{ 15 }$$

is a Bernoulli invariant. The J index is in analogy with the Jacobi constant of a test particle, and is also meant to avoid confusion with Oort's second constant, which we shall introduce later. This invariant is conserved along the streamlines of the domains over which Equations (6)–(8) hold, that is to say over the domains that do not contain shocks. We assume this to be the case, in particular, for the co-orbital region of sufficiently low-mass planets.

In previous two-dimensional analyses (Casoli & Masset 2009; Masset & Casoli 2009), the width of the horseshoe region far from the planet can be determined using the fact that Bernoulli's invariant is the same at the stagnation point and far away on the separatrices. From Equation (15) this value of the Bernoulli invariant depends exclusively on the enthalpy of the fluid at the stagnation point and its location. We discuss below whether a similar property can be generalized to the three-dimensional case.

The first question to address is therefore that of the set of stagnation points. A stagnation point, in three dimensions, is found wherever the three components of the velocity simultaneously vanish. In general, the constraint ${v}_{\phi }=0$ is verified on a two-dimensional manifold (the corotation sheet), within which the additional constraint v_R = 0 yields a one-dimensional manifold. Then, in general, the additional constraint v_z = 0 yields a finite number of points within this one-dimensional manifold. We therefore expect to have a finite number of stagnation points in three dimensions. The streamline analysis performed on data of three-dimensional calculations, presented in Section 4.2, corroborates this statement.

In the following we shall assume that in the three-dimensional case the horseshoe region is still bounded by a well-defined set of streamlines, the separatrix sheets, a three-dimensional generalization of the separatrix streamlines of the two-dimensional cases. We will come back to this assumption in the discussion of Section 6.4. The second question is whether all fluid elements of a separatrix sheet go through a stagnation point. While the answer to this question is trivial in the two-dimensional case, it is not that obvious in the three-dimensional case. In particular, one can imagine that a given fluid element of a separatrix, originating at a given altitude, will complete a horseshoe U-turn without ever reaching the altitude of the stagnation point. Assume that such a fluid element exists. By hypothesis, the norm of its velocity is therefore finite everywhere on its associated streamline.⁴ Then consider another fluid element on a neighboring streamline, initially close to the first one. The separation ${\boldsymbol{\xi }}$ between the two fluid elements then obeys $d{\boldsymbol{\xi }}/{dt}={{\boldsymbol{v}}}_{2}-{{\boldsymbol{v}}}_{1}$, where v₂ and v₁ are the velocities of the fluid elements along their trajectories. One can recast this derivative in terms of the curvilinear abscissa s along the trajectory of the first fluid element: $d{\boldsymbol{\xi }}/{ds}=({{\boldsymbol{v}}}_{2}-{{\boldsymbol{v}}}_{1})/{v}_{1}$. Since v₁ remains finite and since the velocity field is continuous, $d{\boldsymbol{\xi }}/{ds}$ can be made arbitrarily small provided the fluid elements are sufficiently close initially. Upon integration over any finite length s, this means that two neighboring streamlines can always remain arbitrarily close to each other, provided their initial separation is adequately chosen. This is in contradiction with the character of a separatrix: two neighboring streamlines on opposite sides of the separatrix sheet always follow divergent paths after a finite length, no matter how close they are initially. This shows that the modulus of the velocity vector along a separatrix streamline must vanish somewhere, which implies the passage through a stagnation point. There is an immediate consequence to this statement: the streamlines of the separatrix sheets have significant vertical excursions in order to pass through one of the stagnation points. By continuity, the horseshoe trajectories are vertically bent near the U-turns, especially those that are close to the separatrix. This has already been observed in three-dimensional calculations (D'Angelo et al. 2003; Fung et al. 2015).

There is another consequence to the fact that all separatrices go through at least one of the stagnation points: the value of Bernoulli's invariant on any separatrix streamline must be the value that it has at the corresponding stagnation point. By spawning streamlines (downstream and upstream) in all directions in the immediate vicinity of a stagnation point, one generates a separatrix sheet with a unique value of Bernoulli's invariant.

As in the two-dimensional case, there may be several stagnation points, each with its own value of the Bernoulli invariant. Each of these points generates a separatrix sheet where Bernoulli's invariant has everywhere the same value as at the point. We note that two separatrix sheets with different values of Bernoulli's invariant cannot cross, so they are disjoint or nested. We also show in Appendix B that a given separatrix sheet cannot be connected to two stagnation points with different values of Bernoulli's invariant.

In what follows, we are interested in the separatrix sheet that lies furthest from the orbit, corresponding to the widest horseshoe streamlines. This is to be compared to the two-dimensional situation, which can feature two X-stagnation points (Casoli & Masset 2009; Guilet et al. 2013), but the one that determines the overall horseshoe dynamics is the one that has the lowest value of the Bernoulli invariant (and hence is connected to streamlines that lie the furthest from the orbit).

An important hypothesis is that the flow is barotropic: all variables of Equation (15) are continuous in a barotropic fluid, hence the value of Bernoulli's invariant is well defined in the vicinity of a stagnation point. In a baroclinic situation (if, for instance, the flow obeyed an energy equation), the enthalpy, and hence Bernoulli's invariant, would not necessarily be continuous.

The fact that Bernoulli's invariant is uniform on a separatrix sheet provides an idea of the vertical shape of the horseshoe region far from the planet. Considering that the flow at a large azimuthal distance from the planet is unperturbed, we can simplify the expression of Bernoulli's invariant given by Equation (15) as

$$\begin{eqnarray}&&{B}_{J}=\displaystyle \frac{1}{2}{v}_{\phi }^{2}+{{\rm{\Phi }}}_{*}-\displaystyle \frac{1}{2}{R}^{2}{{\rm{\Omega }}}_{p}^{2}+{c}_{s}^{2}\mathrm{log}\left(\displaystyle \frac{{\rho }_{0}}{{\rho }_{00}}\right),\end{eqnarray} \tag{ 16 }$$

where ${{\rm{\Phi }}}_{*}({\boldsymbol{r}})$ is the stellar potential. As shown in Appendix A, the azimuthal velocity in a globally isothermal disk in rotational and hydrostatic equilibrium does not depend on the altitude (Equation (84)), so the vertical derivative of Bernoulli's invariant reduces to

$$\begin{eqnarray}&&{\partial }_{z}{B}_{J}={\partial }_{z}{{\rm{\Phi }}}_{*}+{c}_{s}^{2}{\partial }_{z}\mathrm{log}\left(\displaystyle \frac{{\rho }_{0}}{{\rho }_{00}}\right).\end{eqnarray} \tag{ 17 }$$

The right-hand side of Equation (17) cancels out, as can be seen from Equation (8) for an axisymmetric disk in equilibrium, hence Bernoulli's invariant in the unperturbed disk depends only on the cylindrical radius. The horseshoe's separatrix sheets should therefore be cylinders coaxial with the disk's rotation axis, and so should the corotation sheet, since v_ϕ also depends exclusively on the cylindrical radius. The horseshoe region is therefore expected to have a constant width with altitude. This has been seen by Fung et al. (2015), who describe the horseshoe flow as columnar. This property will be corroborated by numerical simulations in Section 4. This allows us to introduce unambiguously the (unique) half-width of the horseshoe region, which we denote by x_s as in previous two-dimensional work.

Note that the cylindrical shape of the corotation sheet is a consequence of the disk being globally isothermal. Expanding in z the relations given in Appendix A, we can see that for the case in which the temperature increases outward, the corotation has its smallest radius at the midplane, and the opposite holds for the more realistic case corresponding to a temperature decreasing outward: its largest radius is at the midplane.

3.3.1. Domain of Interest

As in Casoli & Masset (2009) and Masset & Casoli (2009), we define a domain of interest ${ \mathcal D }$ that encloses the planet and the horseshoe region, and extends from $z=-\infty$ to $z=+\infty$, from $\phi ={\phi }_{-}$ to $\phi ={\phi }_{+}$, where ${\phi }_{-}\lt 0$ and ${\phi }_{+}\gt 0$, and where $| {\phi }_{-}|$ and ${\phi }_{+}$ are chosen sufficiently large so that the planetary potential is negligible at $\phi ={\phi }_{\pm }$. Similarly, we chose the radial boundaries of the domain, R_min and R_max, to cover an extent much wider than the horseshoe region. Our domain therefore encloses not only horseshoe streamlines but also streamlines that do not perform a U-turn, and which belong to the inner or outer disk. We make the assumption that the horseshoe flow reaching ${\phi }_{\pm }$ is relaxed in the sense that the flow velocity is exclusively azimuthal, and the azimuthal derivatives of the hydrodynamic variables vanish. Naturally, a given fluid element starting its horseshoe trajectory at a given altitude z₀ may change its altitude during its journey to the downstream part of the horseshoe flow. As we shall see in numerical simulations, it may also oscillate vertically a few times after the U-turn. We make no restrictive assumption concerning its final altitude, which may be different from z₀. Our only assumption at this stage is that its altitude has converged to a constant value by the time it reaches the azimuth ${\phi }_{\pm }$, which is why we must consider sufficiently large values of $| {\phi }_{\pm }|$.

The torque exerted by the material enclosed in the domain ${ \mathcal D }$ on the planet is

$$\begin{eqnarray}&&{\rm{\Gamma }}=\int \int {\displaystyle \int }_{{ \mathcal D }}\rho {\partial }_{\phi }{\rm{\Phi }}{RdRdzd}\phi .\end{eqnarray} \tag{ 18 }$$

Using Equations (7), then (5), under our assumption of a steady state, this integral can be recast as integrals that account for the angular momentum budget on the faces of the domain:

$$\begin{eqnarray}{\rm{\Gamma }} & = & -{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{\phi }_{-}}^{{\phi }_{+}}d\phi \rho {v}_{R}{jR}\right|}_{{R}_{\mathrm{max}}}\\ & & +{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{\phi }_{-}}^{{\phi }_{+}}d\phi \rho {v}_{R}{jR}\right|}_{{R}_{\mathrm{min}}}\\ & & -{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{R}_{\mathrm{min}}}^{{R}_{\mathrm{max}}}{dR}(\rho {v}_{\phi }j+{RP})\right|}_{{\phi }_{+}}\\ & & +{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{R}_{\mathrm{min}}}^{{R}_{\mathrm{max}}}{dR}(\rho {v}_{\phi }j+{RP})\right|}_{{\phi }_{-}}.\end{eqnarray} \tag{ 19 }$$

The edge terms in z cancel out since the density vanishes at larger altitude. The integrals of Equation (19) involve angular momentum fluxes on the arc-of-cylinder faces of the domain (in $R={R}_{\mathrm{min}}$ and $R={R}_{\mathrm{max}}$) and on the radial faces of the domain (in $\phi ={\phi }_{\pm }$).

3.3.2. Corotation Torque Integral

We now assume that the angular momentum flux on the arc-of-cylinder faces is the flux carried by the spiral wake, which corresponds to the Lindblad torque on the planet, and that the fluxes on the radial faces of the domain correspond, on the contrary, to the corotation torque. We assume that there is a range of aspect ratios $({R}_{\mathrm{max}}-{R}_{\mathrm{min}})/[{r}_{p}({\phi }_{+}-{\phi }_{-})]$ of the domain for which this separation of the torques is correct, and we focus exclusively on the flux on the radial faces, which we interpret as the full, nonlinear corotation torque:

$$\begin{eqnarray}{{\rm{\Gamma }}}_{\mathrm{CR}} & = & -{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{R}_{\mathrm{min}}}^{{R}_{\mathrm{max}}}{dR}(\rho {v}_{\phi }j+{RP})\right|}_{{\phi }_{+}}\\ & & +{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{R}_{\mathrm{min}}}^{{R}_{\mathrm{max}}}{dR}(\rho {v}_{\phi }j+{RP})\right|}_{{\phi }_{-}}.\end{eqnarray} \tag{ 20 }$$

We make a final assumption: there is no shock in the domain of interest ${ \mathcal D }$. As the latter includes the planetary wake in the vicinity of the planet, this domain should be narrow enough so as not to reach the location at which the wake eventually shocks as a consequence of wave steepening. This location can be dangerously close to the planet owing to the differential rotation of the disk (Goodman & Rafikov 2001), so we focus on planetary masses largely below the thermal mass. This assumption is required, strictly speaking, to satisfy our steady-state assumption. Should shocks be present over the domain, they would transfer the planetary torque to the disk material and they would carve a gap (Rafikov 2002) over a timescale longer than the typical time considered here, intermediate between the horseshoe U-turn and libration timescales, precluding the existence of a steady state. Arguably shocks could still appear beyond the radial limits of the domain ${ \mathcal D }$. We assume that radial boundary conditions are tailored to ensure a steady state over ${ \mathcal D }$.

Using Equation (16) and the assumption that the motion is exclusively azimuthal at large azimuthal distances from the planet, we write

$$\begin{eqnarray}&&{\partial }_{R}{B}_{J}={\zeta }_{z}{v}_{\phi },\end{eqnarray} \tag{ 21 }$$

and we change our integration variable to B_J so as to transform Equation (20) into

$$\begin{eqnarray}{{\rm{\Gamma }}}_{\mathrm{CR}} & = & -{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{B}_{m}}^{{B}_{c}}{{dB}}_{J}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{+},R\lt {R}_{c}}\\ & & -{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{B}_{m}}^{{B}_{c}}{{dB}}_{J}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{+},R\gt {R}_{c}}\\ & & +{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{B}_{m}}^{{B}_{c}}{{dB}}_{J}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{-},R\lt {R}_{c}}\\ & & +{{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{B}_{m}}^{{B}_{c}}{{dB}}_{J}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{-},R\gt {R}_{c}},\end{eqnarray} \tag{ 22 }$$

where B_m represents the minimum value of Bernoulli's invariant, at larger distance from corotation (without loss of generality we can assume that it has the same value at the outer and inner sides), B_c represents the value of Bernoulli's invariant at corotation, and B_s its value on the separatrix sheet.

Note that the relation between R and B_J is not one-to-one. Equation (21) shows that the radial derivative of Bernoulli's invariant changes sign at corotation, where B_J is maximum. We therefore had to split each term of Equation (20) in two, specifying whether we consider contributions from inside or outside corotation. We number one to four, in their order of appearance, the four terms on the right-hand side of Equation (22). Terms one and four correspond to downstream material leaving the domain, while terms two and three correspond to upstream material entering it. Regardless of the complexity of the flow in the vicinity of the planet, a fluid element entering the domain leaves it either on the same side of corotation or on the other side.

We examine first the case of a fluid element leaving the domain on the same side of corotation ("circulating" rather than "librating"). For the sake of definiteness we assume it enters the domain in ${\phi }_{+}$, $R\gt {r}_{c}$. It therefore leaves it in ${\phi }_{-}$, $R\gt {r}_{c}$. An example of the stream tube defined by such a fluid element is depicted in blue in Figure 1. The contribution of this term to the second integral of Equation (22) is

$$\begin{eqnarray}&&d{{\rm{\Gamma }}}_{\mathrm{CR}}^{+}={-d{B}_{J}{{dz}}_{+}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{+}}\end{eqnarray} \tag{ 23 }$$

and its contribution to the fourth one is

$$\begin{eqnarray}&&d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}={{{dB}}_{J}{{dz}}_{-}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{-}}.\end{eqnarray} \tag{ 24 }$$

Since we assumed there are no shocks on the domain, the fluid elements have their Bernoulli invariant materially conserved. We can work out the net flux contribution by summing Equations (23) and (24).

Zoom In Zoom Out Reset image size

The fluid element may exit the domain at a radius and altitude different from those at which it entered the domain, and its velocity and enthalpy may also be different. We assume the variations of radius, velocity, and enthalpy to be sufficiently small to allow first-order expansions. We denote by $\delta X$ the small variation of a quantity X at the location of the fluid element, between its arrival at and exit from the domain (i.e., the Lagrangian variation seen by the fluid element), and by ${\delta }^{\prime }X$ the variation of that quantity between the two radial faces, at given radius and altitude. Expanding Equation (16), making use of Equations (6) and (8) together with the assumptions we laid down at large azimuthal distance from the planet (${v}_{R}={v}_{z}=0$, ${\partial }_{\phi }\equiv 0$), we are led to

$$\begin{eqnarray}&&\delta {B}_{J}={\left(\displaystyle \frac{{v}_{\phi }}{R}\right)}_{{\phi }_{-}}\delta j+{\delta }^{\prime }\eta ,\end{eqnarray} \tag{ 25 }$$

with $\delta j={j}_{-}-{j}_{+}$ and ${\delta }^{\prime }\eta =\eta (R,z,{\phi }_{-})-\eta (R,z,{\phi }_{+})$. Since $\delta {B}_{J}=0$, we have

$$\begin{eqnarray}d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-} & = & \displaystyle \frac{{{dB}}_{J}{{dz}}_{-}}{{\zeta }_{z}^{-}}{\left[{\rho }_{-}\left({j}_{+}-\displaystyle \frac{R{\delta }^{\prime }\eta }{{v}_{\phi }}\right)+\displaystyle \frac{{RP}}{{v}_{\phi }}\right]}_{{\phi }_{-}}\\ & = & \displaystyle \frac{{{dB}}_{J}{{dz}}_{-}}{{w}_{z}^{-}}{j}_{+}+{{dR}}_{-}{{dz}}_{-}{R}_{-}P({R}_{-},{z}_{-},{\phi }_{+}),\end{eqnarray} \tag{ 26 }$$

where we have used the first-order expansion ${\delta }^{\prime }P={\rho }_{-}{\delta }^{\prime }\eta$, and where ${{dR}}_{\pm }$ is the radial width of the stream tube at azimuth ${\phi }_{\pm }$. In order to further transform the first term, we use Ertel's theorem⁵ (Ertel 1942). Ertel's theorem, in a reduced form that is of interest here, states that in a barotropic flow we have the conservation law

$$\begin{eqnarray}&&{D}_{t}\left(\displaystyle \frac{{\boldsymbol{\zeta }}}{\rho }\cdot {\boldsymbol{\nabla }}\psi \right)=0,\end{eqnarray} \tag{ 27 }$$

where ψ is any quantity materially conserved by the flow (${D}_{t}\psi =0$). We can take here for ψ an arbitrary function whose derivative in z does not vanish, such as ${z}_{+}$ itself, regarded as a Lagrangian, passive scalar on the fluid elements. Ertel's theorem in this case reads

$$\begin{eqnarray}&&{\left({\boldsymbol{w}}\cdot {\boldsymbol{\nabla }}{z}_{+}\right)}_{{\phi }_{+}}={\left({\boldsymbol{w}}\cdot {\boldsymbol{\nabla }}{z}_{+}\right)}_{{\phi }_{-}},\end{eqnarray} \tag{ 28 }$$

or

$$\begin{eqnarray}&&{w}_{z}^{+}={w}_{z}^{-}\displaystyle \frac{{{dz}}_{+}}{{{dz}}_{-}}+{w}_{R}^{-}\displaystyle \frac{{{dz}}_{+}}{{dR}}.\end{eqnarray} \tag{ 29 }$$

Note that in Equation (29) there is no azimuthal derivative, as per our assumption that these vanish at large azimuthal distance from the planet. The factor ${w}_{R}^{-}$ in the last term is expected to be much smaller than ${w}_{z}^{\pm }$ (we will see in Section 4 that the vorticity tilt angle, especially in the circulating region, is well below 10⁻²), and we expect the second factor of this term (${{dz}}_{+}/{dR}$, which represents the inclination of an initially horizontal sheet of material after the interaction with the planet) to also be a small number. We therefore write

$$\begin{eqnarray}&&\displaystyle \frac{{{dz}}_{+}}{{w}_{z}^{+}}\approx \displaystyle \frac{{{dz}}_{-}}{{w}_{z}^{-}},\end{eqnarray} \tag{ 30 }$$

with an accuracy largely better than the 10⁻² level. The net contribution to the corotation torque of the stream tube that we considered is therefore

$$\begin{eqnarray}d{{\rm{\Gamma }}}_{\mathrm{CR}} & = & d{{\rm{\Gamma }}}_{\mathrm{CR}}^{+}+d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}={{dR}}_{-}{{dz}}_{-}{R}_{-}P({R}_{-},{z}_{-},{\phi }_{+})\\ & & -{{dR}}_{+}{{dz}}_{+}{R}_{+}P({R}_{+},{z}_{+},{\phi }_{+}).\end{eqnarray} \tag{ 31 }$$

We see that the torque contribution of our circulating fluid element only amounts to the pressure on the upstream side: both terms above are evaluated in ${\phi }_{+}$ thanks to the transformation arising from Equation (25). Anticipating what follows, this pressure contribution will cancel out upon the integration over all stream tubes, so that the net contribution will arise exclusively from the horseshoe region, which we now examine.

We consider a stream tube that executes a horseshoe U-turn (corresponding to "librating" fluid elements). For the sake of definiteness we assume it enters the domain of interest in $R\lt {r}_{c}$ and ${\phi }_{-}$ and exits it in $R\gt {r}_{c}$ and ${\phi }_{-}$, like the one depicted in purple in Figure 1. Such a stream tube yields contributions to the third and fourth integrals of Equation (22), which read respectively

$$\begin{eqnarray}&&{d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}\right|}_{\mathrm{in}}={{{dB}}_{J}{{dz}}_{-}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{-},\mathrm{in}}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}\right|}_{\mathrm{out}}={{{dB}}_{J}{{dz}}_{-}\displaystyle \frac{\rho {v}_{\phi }j+{RP}}{{\zeta }_{z}{v}_{\phi }}\right|}_{{\phi }_{-},\mathrm{out}},\end{eqnarray} \tag{ 33 }$$

where the subscript in or out specifies whether we consider a region inside or outside corotation. In the above analysis of a circulating stream tube, from Equations (23) to (31), fluid elements did not cross corotation, so there was no need to specify explicitly the location of each term as they were all evaluated systematically outside corotation. As previously, we can recast Equation (32) as

$$\begin{eqnarray}&&{d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}\right|}_{\mathrm{in}}=\displaystyle \frac{{{dB}}_{J}{{dz}}_{\mathrm{in}}^{-}}{{w}_{z}^{\mathrm{in},-}}{j}_{\mathrm{in}}^{-}-{{dR}}_{\mathrm{in}}^{-}{{dz}}_{\mathrm{in}}^{-}{R}_{\mathrm{in}}^{-}P({\phi }_{-},{R}_{\mathrm{in}}^{-},{z}_{\mathrm{in}}^{-}).\end{eqnarray} \tag{ 34 }$$

We can write Equation (33) in a similar way, and transform it using Equation (25), for which we set $\delta {B}_{J}=0$:

$$\begin{eqnarray}{d{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}\right|}_{\mathrm{out}} & = & \displaystyle \frac{{{dB}}_{J}{{dz}}_{\mathrm{out}}^{-}}{{w}_{z}^{\mathrm{out},-}}{j}_{\mathrm{out}}^{-}+{{dR}}_{\mathrm{out}}^{-}{{dz}}_{\mathrm{out}}^{-}{R}_{\mathrm{out}}^{-}P({\phi }_{-},{R}_{\mathrm{out}},{z}_{\mathrm{out}})\\ & = & \displaystyle \frac{{{dB}}_{J}{{dz}}_{\mathrm{out}}^{-}}{{w}_{z}^{\mathrm{out},-}}{j}_{\mathrm{out}}^{+}+{{dR}}_{\mathrm{out}}^{-}{{dz}}_{\mathrm{out}}^{-}{R}_{\mathrm{out}}^{-}P({\phi }_{+},{R}_{\mathrm{out}},{z}_{\mathrm{out}}).\end{eqnarray} \tag{ 35 }$$

We see that this transformation allows us to express the torque contribution in terms of the upstream pressure $P({\phi }_{+},{R}_{\mathrm{out}}^{-},{z}_{\mathrm{out}}^{-})$ and specific angular momentum ${j}_{\mathrm{out}}^{+}$ on the outer side. Since the transformation of Equation (35) has been obtained by letting $\delta {B}_{J}=0$, the quantity ${j}_{\mathrm{out}}^{+}$ represents the specific angular momentum of the upstream material in the outer disk that has same Bernoulli invariant as our horseshoe stream tube.

When summing all contributions of Equation (22), it is now clear that the contributions of the pressure torques cancel out, since they all amount to summing the upstream pressure field: in the outer disk, Equation (31) for the circulating fluid elements and Equation (35) for the librating ones show that we only consider the upstream pressure ($P({\phi }_{+})$). For the inner disk, transformations similar to those considered previously would yield torque contributions as a function of $P({\phi }_{-})$, such as the contribution given by Equation (32). Since the advective contributions to the angular momentum cancel out for circulating streamlines, as shown by Equation (31), we are only left with the contributions to the angular momentum flux of the librating streamlines. For the rear horseshoe, the total contribution therefore amounts to

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}={\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{B}_{s}}^{{B}_{c}}\displaystyle \frac{{{dB}}_{J}}{{w}_{z}^{\mathrm{in},-}}{j}_{\mathrm{in}}^{-}-{\displaystyle \int }_{-\infty }^{+\infty }{dz}{\displaystyle \int }_{{B}_{s}}^{{B}_{c}}\displaystyle \frac{{{dB}}_{J}}{{w}_{z}^{\mathrm{out},-}}{j}_{\mathrm{out}}^{+}.\end{eqnarray} \tag{ 36 }$$

We now exchange the integral sums, which yields

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}={\displaystyle \int }_{{B}_{s}}^{{B}_{c}}{{dB}}_{J}{\displaystyle \int }_{-\infty }^{+\infty }{dz}\left(\displaystyle \frac{{j}_{\mathrm{in}}^{-}}{{w}_{z}^{\mathrm{in},-}}-\displaystyle \frac{{j}_{\mathrm{out}}^{+}}{{w}_{z}^{\mathrm{out},-}}\right).\end{eqnarray} \tag{ 37 }$$

This expression exclusively features upstream values of the specific angular momentum. As we consider the unsaturated horseshoe drag, we assume that it has the same value as in the unperturbed disk. It is therefore independent of altitude, and can be taken out of the integral over z:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}={\displaystyle \int }_{{B}_{s}}^{{B}_{c}}{{dB}}_{J}\left({j}_{\mathrm{in}}^{-}{\displaystyle \int }_{-\infty }^{+\infty }\displaystyle \frac{{dz}}{{w}_{z}^{\mathrm{in},-}}-{j}_{\mathrm{out}}^{+}{\displaystyle \int }_{-\infty }^{+\infty }\displaystyle \frac{{dz}}{{w}_{z}^{\mathrm{out},-}}\right).\end{eqnarray} \tag{ 38 }$$

We again apply Ertel's theorem with the ancillary function $\psi ={z}_{\mathrm{in}}^{-}$, the initial altitude of a fluid element, which gives

$$\begin{eqnarray}&&{w}_{z}^{\mathrm{in},-}={w}_{z}^{\mathrm{out},-}\displaystyle \frac{\partial {z}_{\mathrm{in}}^{-}}{\partial {z}_{\mathrm{out}}^{-}}+{w}_{R}^{\mathrm{out},-}\displaystyle \frac{\partial {z}_{\mathrm{in}}^{-}}{\partial {R}_{\mathrm{out}}}\approx {w}_{z}^{\mathrm{out},-}\displaystyle \frac{\partial {z}_{\mathrm{in}}^{-}}{\partial {z}_{\mathrm{out}}^{-}},\end{eqnarray} \tag{ 39 }$$

where we have again assumed that the radial tilt of the vortex tubes is negligible. We can therefore write

$$\begin{eqnarray}&&{\displaystyle \int }_{-\infty }^{+\infty }\displaystyle \frac{{{dz}}_{\mathrm{out}}}{{w}_{z}^{\mathrm{out},-}}={\displaystyle \int }_{-\infty }^{+\infty }\displaystyle \frac{{{dz}}_{\mathrm{out}}}{{w}_{z}^{\mathrm{in},-}}\displaystyle \frac{\partial {z}_{\mathrm{in}}^{-}}{\partial {z}_{\mathrm{out}}^{-}}={\displaystyle \int }_{-\infty }^{+\infty }\displaystyle \frac{{{dz}}_{\mathrm{in}}}{{w}_{z}^{\mathrm{in},-}}.\end{eqnarray} \tag{ 40 }$$

We conclude that the integral $\int {w}_{z}^{-1}{dz}$ has the same value in the upstream and downstream flow, far from the planet. The expression of this integral can be made slightly simpler by making use of the inverse of the vertical component of the PV, which we dubbed the load in an earlier paper (Masset & Casoli 2010):

$$\begin{eqnarray}&&l=\displaystyle \frac{\rho }{{\zeta }_{z}},\end{eqnarray} \tag{ 41 }$$

and we denote by L the vertical integral of the load:

$$\begin{eqnarray}&&L={\displaystyle \int }_{-\infty }^{+\infty }l(z){dz}.\end{eqnarray} \tag{ 42 }$$

We can now transform Equation (37) into

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}^{-}={\displaystyle \int }_{{B}_{s}}^{{B}_{c}}{{dB}}_{J}{L}_{\mathrm{in}}({B}_{J})({j}_{\mathrm{in}}^{-}-{j}_{\mathrm{out}}^{+}).\end{eqnarray} \tag{ 43 }$$

The jump in specific angular momentum during a U-turn, ${j}_{\mathrm{in}}^{-}-{j}_{\mathrm{out}}^{+}$, is expressed in terms of the values in the unperturbed disk. It is a function of Bernoulli's invariant exclusively, which we denote $-{\rm{\Delta }}{j}_{0}$. The horseshoe drag arising from the front part reads similarly:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}^{+}={\displaystyle \int }_{{B}_{s}}^{{B}_{c}}{{dB}}_{J}{L}_{\mathrm{out}}({B}_{J}){\rm{\Delta }}{j}_{0},\end{eqnarray} \tag{ 44 }$$

and the net corotation torque is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}={\displaystyle \int }_{{B}_{s}}^{{B}_{c}}{{dB}}_{J}({L}_{\mathrm{out}}-{L}_{\mathrm{in}}){\rm{\Delta }}{j}_{0}.\end{eqnarray} \tag{ 45 }$$

This integral, which only contains quantities of the unperturbed disk, can be evaluated using Equation (21), which can be recast to lowest order as ${\partial }_{R}{B}_{J}=4{A}_{p}{B}_{p}x$ and gives the relationship

$$\begin{eqnarray}&&{B}_{J}={B}_{c}+2{A}_{p}{B}_{p}{x}^{2},\end{eqnarray} \tag{ 46 }$$

where $A=\frac{1}{2}r{\partial }_{r}{\rm{\Omega }}$ and $B=\frac{1}{2r}{\partial }_{r}{j}_{0}$ are respectively Oort's first and second constants. Changing our variable of integration to x, we are led to

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{CR}}}=8| {A}_{p}| {B}_{p}^{2}{r}_{p}{\partial }_{R}{L}_{0}{x}_{s}^{4}.\end{eqnarray} \tag{ 47 }$$

Defining as in earlier works the dimensionless vortensity gradient at the planet's location as

$$\begin{eqnarray}&&{ \mathcal V }={r}_{p}\displaystyle \frac{{\partial }_{R}{L}_{0}}{{L}_{0}},\end{eqnarray} \tag{ 48 }$$

and assuming that the flow's vorticity barely depends on altitude in the unperturbed disk, so that ${L}_{0}\approx {{\rm{\Sigma }}}_{p}/(2{B}_{p})$, we are eventually led to

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{CR}}}=\displaystyle \frac{3}{4}{{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{p}^{2}{ \mathcal V }{x}_{s}^{4},\end{eqnarray} \tag{ 49 }$$

which is the exact same horseshoe drag expression as in the two-dimensional case.

We conclude that the horseshoe dynamics in a three-dimensional, globally isothermal disk, despite the complexity of the three-dimensional flow, retain some of the simplicity of the two-dimensional case, first because the width of the horseshoe region is independent of altitude, which allows us to use unambiguously a unique half-width x_s, and second because the drag expression is the same as in the two-dimensional case, provided an adequate definition of the vortensity is used. Equations (41), (42), and (47) show that the vortensity V in a three-dimensional disk must be defined as

$$\begin{eqnarray}&&V={\left[{\displaystyle \int }_{-\infty }^{+\infty }{\left(\displaystyle \frac{{\zeta }_{z}}{\rho }\right)}^{-1}{dz}\right]}^{-1}.\end{eqnarray} \tag{ 50 }$$

In order to illustrate and check the results of Section 3, we have undertaken numerical simulations of a globally isothermal disk with a low-mass planet embedded on a fixed circular prograde orbit, coplanar with the disk. We used for that purpose the public code FARGO3D⁶ (Benítez-Llambay & Masset 2016), with the setup p3diso, meant to describe a (locally or globally) isothermal gas in orbit around a point-like mass, on a spherical mesh. The use of a spherical mesh may seem not adapted to the problem at hand, since we expect from Section 3 a cylindrical symmetry for a number of variables, such as the azimuthal velocity or Bernoulli's invariant. In a similar manner, the horseshoe separatrices are expected to be cylinders. Our choice of a spherical mesh instead of a cylindrical one fulfills two purposes:

• 1.  
  owing to the strong flaring of a globally isothermal disk ($H\propto {R}^{3/2}$), we avoid very empty regions at lower radius and high altitude;
• 2.  
  we can discard possible mesh effects if we find features with cylindrical symmetry.

Our mesh extends from $-\pi$ to $+\pi$ in azimuth, $0.65{r}_{p}$ to $1.35{r}_{p}$ in radius, and $\pi /2-3{h}_{p}$ to $\pi /2$ in colatitude, where h_p = 0.05 is the disk's aspect ratio at the orbital radius of the planet, so that 99% of the disk's mass at the planet's location lies within the mesh. As we simulate only one hemisphere of the disk, we adopt reflecting boundary conditions at the equator. At high altitude we extrapolate the density and azimuthal velocity fields using the analytic profiles of Appendix A to fill the ghost zones, whereas we use symmetric boundary conditions on the radial velocity and antisymmetric boundary conditions on the velocity in colatitude. Our mesh size respectively in azimuth, radius, and colatitude is ${N}_{\phi }=3290$, N_r = 367, and ${N}_{\theta }=78$, with a uniform spacing, so that the cells are approximately cubic at $r={r}_{p}$, with an edge length of $1.9\times {10}^{-3}{r}_{p}$. The planetary potential has the form

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}({\boldsymbol{r}})=-\displaystyle \frac{{{GM}}_{p}}{\sqrt{{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{p}| }^{2}+{\epsilon }^{2}}},\end{eqnarray} \tag{ 51 }$$

where , a softening length used to avoid a divergence at the planet location, is set to $4\times {10}^{-3}{r}_{p}$, which is 8% of the pressure scale height, or two cell sizes. No viscosity is used in the calculations. The planet-to-star mass ratio is $q={10}^{-5}$, so the planet has 8% of the thermal mass ${h}_{p}^{3}{M}_{\star }$.

The disk's initial conditions are given by ${v}_{r}\equiv {v}_{\theta }\equiv 0$, while the density and azimuthal velocity are given respectively by Equations (83) and (84). Damping boundary conditions (de Val-Borro et al. 2006) are used in the radial direction only, over 10% of the radial extent of the mesh (i.e., over annuli of width $0.07{r}_{p}$), at both the inner and outer edges. Each of our calculations is carried out over 20 orbital periods of the planet. We do not use a temporal tapering of the planetary mass upon the planet's insertion in the disk. The manner in which the planet is turned on can have an impact on the flow at larger time in two dimensions (Ormel et al. 2015a), since some vortensity can be created by shocks and trapped in a closed region around the planet, but no such region exists in three dimensions (Ormel et al. 2015b), and our setup does not have the resolution to capture processes that happen on a scale of the planetary Bondi radius.

Our fiducial calculation is one for which $\alpha =3/2$, which corresponds to a (nearly) vanishing vortensity gradient. Our background disk is therefore identical to that of Fung et al. (2015). The main differences are, in our case, a smaller planetary mass and a much coarser resolution at the planet's location.

In addition to this fiducial calculation, we have run 10 other calculations in which we vary the vortensity gradient, defined by Equation (48). Overall we have 11 runs with a vortensity gradient ranging from −1.5 to 1.5 in increments of 0.3.

We perform all our analyses at t  =  20 orbital periods of the planet, the date t  =  0 corresponding to the planet's insertion in the disk.

We determine the position of the separatrices by integrating the path of fluid elements starting at a given colatitude and radius r, with azimuth ϕ = ±1. The integration is carried out downstream for $\phi (r-{r}_{c})\gt 0$ and upstream otherwise. The trilinearly interpolated value of each component of the velocity field is used for the integration. For each colatitude, the starting radius of the widest horseshoe streamline is found with a dichotomic search to machine accuracy.

We present in Figure 2 the width of horseshoe region as a function of the distance to the midplane. As expected from Section 3.2, the width is nearly constant, to within a few per cent, over the whole vertical extent of the computational domain, which covers three vertical scale heights. Also shown in this figure is the width of the horseshoe region from 2D calculations with different softening lengths of the potential. We see that they match for a softening length of $0.65H$. Interestingly, this value is comparable to the softening length required to match the two- and three-dimensional Lindblad torques (Masset 2002; Kley et al. 2012). These results are to be compared with those of Fung et al. (2015), who considered a planet of marginally sub-thermal mass, and who find the horseshoe width reproduced by 2D calculations with a softening length $\sim 0.35H$. The average value of the front and rear widths measured in our fiducial run is

$$\begin{eqnarray}&&{x}_{s}=0.0147{r}_{p}.\end{eqnarray} \tag{ 52 }$$

Zoom In Zoom Out Reset image size

In order to determine the location of the stagnation points, we proceed as outlined in Section 3.2, i.e., by determining the manifolds of decreasing dimension obtained by imposing successively that ${v}_{\phi }=0$, then also v_r = 0, and finally also ${v}_{\theta }=0$. Care must be taken in this procedure with the staggering of the velocity components in the FARGO3D code. Namely, for each of the N_θ conical slices of our mesh, we determine the location at which v_ϕ cancels out (the intersection of the corotation sheet with the cone $\theta ={\theta }_{k}$, with ${\theta }_{k}=\pi /2-3{h}_{p}(k+1/2)/{N}_{\theta }$ and $k\in [0,{N}_{\theta }-1]$). This location turns out to be unique, for our setup and resolution, and results in a function of r for any given azimuth ${({\phi }_{i})}_{i\in [0,{N}_{\phi }-1]}$: $r={r}_{c}({\phi }_{i},{\theta }_{k})$. Next we linearly interpolate v_r at the locations $[{\phi }_{i},{r}_{c}({\phi }_{i},{\theta }_{k}),{\theta }_{k}]$. For each value of k (within each slice in colatitude), we thus construct a sequence v_rⁱ. We then seek the root(s) in ϕ of the piecewise linear function that has value v_rⁱ for $\phi ={\phi }_{i}$. We thus obtain, for each slice in colatitude, the location at which the bilinearly interpolated values of the azimuthal and radial velocities simultaneously vanish. By varying k, we identify "filaments of horizontal stagnation," nearly vertical curves in which the azimuthal and radial velocities cancel out. These filaments are not streamlines, since they can have some tilt in azimuth and radius, whereas the velocity field on them is purely along the colatitude direction, by construction.⁷ There are no sets of stagnation points either, because in general v_θ does not vanish on these filaments, except at specific locations that are determined in our last step, which consists in constructing for each filament the sequence of ${v}_{\theta }^{k}$, the bilinearly interpolated values of v_θ where the filaments intersect the zone edges in colatitude. As previously, we then seek the roots of the piecewise linear function that coincides with ${v}_{\theta }^{k}$ for $\theta =\pi /2-3{h}_{p}k/{N}_{\theta }$. This gives us the location of the stagnation points, where the three velocity components simultaneously vanish. As a sanity check, we verify that the trilinearly interpolated value of each velocity component is indeed very small at the location of each stagnation point. We overall identify 96 stagnation points in our fiducial run at 20 orbits, most of them relatively far for the planet, especially in the tadpole regions. Of special interest is the filament that lies near the planet, almost at the intersection of the separatrix sheets. This filament contains three stagnation points: one at the midplane, and two others within the first pressure scale height.

We show in Figure 3 the separatrix sheets of the horseshoe region. Although the streamlines of these sheets exhibit significant vertical motion toward stagnation points in the vicinity of the stagnation filament, as expected from the considerations of Section 3.2, they never quite reach these points: we find that an accuracy higher than double precision would be required to better constrain the starting radius of a separatrix streamline that would reach the vicinity of the stagnation points. Our assertion of Section 3.2 that Bernoulli's invariant is uniform on the separatrix sheets, and equal to its value at the corresponding stagnation point, therefore needs to be examined in more detail.

Zoom In Zoom Out Reset image size

Figure 4 shows the value of Bernoulli's invariant in the corotation sheet, in the planet's vicinity. The horizontal stagnation filament is shown in this figure. Bernoulli's invariant is found to vary by less than ${10}^{-6}{r}_{p}^{2}{{\rm{\Omega }}}^{2}$, at the filament location, over the first pressure scale height vertically, and it decays by 6–$7\times {10}^{-6}{r}_{p}^{2}{{\rm{\Omega }}}_{p}^{2}$ over the next two scale heights. This overall variation is to be contrasted with the drop in Bernoulli's invariant between corotation (at large distance from the planet) and separatrices, given by Equations (46) and (52), which amounts to $-8\times {10}^{-5}{r}_{p}^{2}{{\rm{\Omega }}}_{p}^{2}$. The variation over the first scale height is therefore at the level of one per cent, whereas the variation at higher altitude is at most 10% of the drop in Bernoulli's invariant between corotation and separatrices. The values of Bernoulli's invariant at the three stagnation points on the filament are, in order of increasing distance to the midplane,

$$\begin{eqnarray*}{B}_{s}^{1} & = & -1.51343737{r}_{p}^{2}{{\rm{\Omega }}}_{p}^{2}\\ {B}_{s}^{2} & = & -1.51343732{r}_{p}^{2}{{\rm{\Omega }}}_{p}^{2}\\ {B}_{s}^{3} & = & -1.51343750{r}_{p}^{2}{{\rm{\Omega }}}_{p}^{2}.\end{eqnarray*}$$

In order to better assess the extent to which Bernoulli's invariant can be regarded as uniform over the separatrix sheets, we define ${B}_{s}=-1.5134374{r}_{p}^{2}{{\rm{\Omega }}}_{p}^{2}$, an average of the values quoted above, and we compare the isosurfaces ${B}_{J}\equiv {B}_{s}$ to the separatrix sheets. Figure 5 shows a global view of a meridional slice of the disk. The upper half of the figure shows isocontours of Bernoulli's invariant at $\phi =1\;{\rm{rad}}$. The lower half shows the separatrix sheets at same azimuth. The cylindrical symmetry of the separatrix sheet and iso-Bernoulli surfaces is readily apparent, and we find an excellent coincidence between the radii of the separatrices and the radii of the isosurfaces of B_s. Here we use a subscript s to convey that B_J has to be evaluated at a stagnation point. In Section 3.3.2, we used a similar notation to represent the (unique) value of B_J on a separatrix sheet. Our numerical results confirm our expectation of Section 3.2 that the two values coincide, hence it is legitimate to use the same notation B_s for both.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows a detailed comparison of the upstream separatrix distance to corotation with the distance of the isosurface ${B}_{J}\equiv {B}_{s}$. A nearly perfect agreement is found over the first pressure scale height, with residual errors far below the mesh resolution. It is likely that the trilinearly interpolated values of Bernoulli's invariant on the one hand, and of the velocity on the other hand, fulfill Equation (14) with a high accuracy, but we have not investigated the reasons for this nearly perfect match in detail. Nonetheless, this remarkable agreement lends confidence in the fact that Bernoulli's invariant can be regarded as uniform on the separatrix sheet, and equal to the value B_s that it has on the stagnation filament. The agreement, as can be seen in Figure 6, is not so good at higher altitudes (albeit with a discrepancy between Bernoulli's isosurface and the separatrices still below the mesh resolution). This corresponds to the altitudes at which Bernoulli's invariant shows a small departure from the uniform value that it has over the first pressure scale height (see Figure 4), which suggests that the flow has not reached a fully steady state at higher altitudes 20 orbital periods after the insertion of the planet. Our streamline integration shows indeed that the high-altitude fluid elements near the separatrix spend a considerable amount of time near the horizontal stagnation filament, where they execute significant vertical excursions at very small velocity (thereby contributing to homogenizing Bernoulli's invariant along this filament). We will come back to the timescale of high-altitude U-turns in Section 6.2, and we will see that they can last much longer than in the midplane. We note that we could not have taken a significantly larger time after the planet's insertion to perform our analysis, because the first effects of phase mixing are expected around t = 40 orbital periods for the horseshoe zone width that we measured.

Zoom In Zoom Out Reset image size

We note in Figure 4 that the filament of horizontal stagnation is almost vertical over the first pressure scale height. We could have defined stagnation filaments by the requirement ${v}_{R}={v}_{\phi }=0$ (that is, the first two components of the velocity in a cylindrical coordinate system vanish). The filament thus defined could nearly be a streamline connecting the stagnation points, providing direct evidence of why they share the same value of Bernoulli's invariant.

Our next step is to check whether Equation (49) gives a reliable estimate of the horseshoe drag. This expression requires the prior knowledge of x_s, the half-width of the horseshoe region. We use for this purpose the set of simulations that we performed with different vortensity gradients, and we assume that the width of the horseshoe region, as in the two-dimensional case, does not sensitively depend on the vortensity gradient (Casoli & Masset 2009) and is given by Equation (52).

We show in Figure 7 the value measured for the torque in our 11 calculations, as a function of the surface density slope, and we compare it to two analytical expressions: one of them is the total torque estimate obtained in the linear regime by TTW02, which reads

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{TTW}02}=-(1.364+0.541\alpha ){{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{p}^{2}{r}_{p}^{4}{q}^{2}{h}_{p}^{-2},\end{eqnarray} \tag{ 53 }$$

and the other one is the sum of the linear estimate of the Lindblad torque of TTW02 and the nonlinear estimate of the corotation torque of Equation (49), in which we set ${ \mathcal V }=\frac{3}{2}-\alpha$. The Lindblad torque estimate of TTW02 is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{TTW}02}^{\mathrm{LR}}=-(2.340-0.099\alpha ){{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{p}^{2}{r}_{p}^{4}{q}^{2}{h}_{p}^{-2}.\end{eqnarray} \tag{ 54 }$$

The value obtained in our calculations are much more compatible with the second estimate. This result extends to the three-dimensional case the findings of Paardekooper & Papaloizou (2009a) that the corotation torque eventually becomes nonlinear at all planetary masses, even those that are largely sub-thermal. In addition, it shows that the horseshoe drag formula of Equation (49) gives an acceptable estimate of the nonlinear corotation torque.

Zoom In Zoom Out Reset image size

The half-width that we have found for the horseshoe region, given by Equation (52), is only 3.5 times larger than the softening length. One can therefore wonder whether the horseshoe region could be substantially wider in the limit of a vanishing softening length. In order to answer this question, we have undertaken two additional calculations with the same setup as our fiducial one, except for the softening length, which was set respectively to $0.06H$ and $0.10H$, instead of our fiducial value of $0.08H$. We find that the horseshoe half-width displays a nearly linear relationship with the softening length, with a very small slope, and that the width extrapolated for a vanishing softening length is only 0.8% wider than in our fiducial calculation: ${x}_{s}=0.0148{r}_{p}$.

This width is marginally smaller than the widely used estimate derived for two-dimensional disks with a softening length of $0.3H$ for the potential, which reads ${x}_{s}=1.16{r}_{p}{(q/h)}^{1/2}$ (Masset et al. 2006). Here, we find a slightly different numerical factor, about 10% smaller, for the half-width of the horseshoe region:

$$\begin{eqnarray}&&{x}_{s}=1.05{r}_{p}\sqrt{\displaystyle \frac{q}{h}}.\end{eqnarray} \tag{ 55 }$$

These findings are consistent with those of Lega et al. (2015), who also observe a 10% smaller width for the horseshoe region than the standard two-dimensional estimate for sub-thermal-mass planets.

In the limit of a small planetary mass, the location of all stagnation points tends to the corotation of the unperturbed disk, since ${v}_{\phi }=R({{\rm{\Omega }}}_{0}-{{\rm{\Omega }}}_{p})+{v}_{\phi }^{\prime }$, where ${v}_{\phi }^{\prime }$ is the perturbed velocity due to the planet, which scales with q when the planetary mass is sufficiently small. We restrict our discussion to this limiting case, which is tantamount to assuming that the corotation sheet is not significantly distorted by the introduction of the planet, at the location of the stagnation points. Denoting by ${{\boldsymbol{r}}}_{s}$ the location of a stagnation point, we can write Bernoulli's constant at ${{\boldsymbol{r}}}_{s}$ prior to the planet's insertion, which reads

$$\begin{eqnarray}&&{B}_{J,0}={{\rm{\Phi }}}_{*}({{\boldsymbol{r}}}_{s})-\displaystyle \frac{1}{2}{R}_{c}{{\rm{\Omega }}}_{p}^{2}+{\eta }_{0}({{\boldsymbol{r}}}_{s})={B}_{c},\end{eqnarray} \tag{ 56 }$$

and after the planet's insertion, once a steady state is reached in the planet frame,

$$\begin{eqnarray}&&{B}_{J}={{\rm{\Phi }}}_{*}({{\boldsymbol{r}}}_{s})+{{\rm{\Phi }}}_{p}({{\boldsymbol{r}}}_{s})-\displaystyle \frac{1}{2}{R}_{c}{{\rm{\Omega }}}_{p}^{2}+\eta ({{\boldsymbol{r}}}_{s})={B}_{s}.\end{eqnarray} \tag{ 57 }$$

We therefore have

$$\begin{eqnarray}&&{B}_{s}-{B}_{c}={{\rm{\Phi }}}_{p}({{\boldsymbol{r}}}_{s})+{\eta }^{\prime }({{\boldsymbol{r}}}_{s}),\end{eqnarray} \tag{ 58 }$$

where ${\eta }^{\prime }=\eta -{\eta }_{0}$ is the perturbed enthalpy. Using a notation similar to Masset et al. (2006), we write

$$\begin{eqnarray}&&{\tilde{{\rm{\Phi }}}}^{\prime }={{\rm{\Phi }}}_{p}+{\eta }^{\prime },\end{eqnarray} \tag{ 59 }$$

which corresponds to the perturbed value of the effective potential of Equation (10). Using Equation (46), and specifying to the Keplerian case, for which ${A}_{p}=-3{{\rm{\Omega }}}_{p}/4$ and ${B}_{p}={{\rm{\Omega }}}_{p}/4$, we obtain

$$\begin{eqnarray}&&{x}_{s}=\displaystyle \frac{1}{{{\rm{\Omega }}}_{p}}\sqrt{-\displaystyle \frac{8}{3}{\tilde{{\rm{\Phi }}}}^{\prime }({{\boldsymbol{r}}}_{s})},\end{eqnarray} \tag{ 60 }$$

which is the same as Equation (12) of Masset et al. (2006). The value of the perturbed effective potential at the stagnation point is directly related to the width of the horseshoe region, and may be used to infer the latter from numerical simulations, without resorting to streamline analysis.

In the derivation of Equation (49), we had to assume that the vortex tubes had a small tilt angle with respect to the vertical direction in the downstream flow. We examine here how justified this assumption was. Figure 8 shows the tilt angle in the radial direction (the tilt angle in the azimuthal direction does not feature in our torque derivation owing to the assumption that at large azimuthal distance the flow variables have a vanishing azimuthal derivative). We see in this figure that slightly larger values of tilt angle are systematically found in the downstream flow (inner side for $\phi =1\;{\rm{rad}}$, outer side for $\phi =-1$ rad). Nevertheless, the values observed are very small, largely below 10⁻² over most of the horseshoe flow, especially over the first pressure scale height, whose contribution dominates the horseshoe drag. Our assumption was therefore largely satisfied for our fiducial run.

Zoom In Zoom Out Reset image size

We finally mention that, as expected, the PV is not conserved in the flow. This is illustrated in Figure 9 where we see stripes in the downstream distribution of the vertical component of PV (at the midplane), which cannot be accounted for by advection of the initial distribution by the horseshoe flow. These stripes can be traced back to vertical oscillations of the material upon execution of their U-turn. Compression (expansion) of material along a vortex tube (i.e., essentially in the vertical direction, in our case) decreases (increases) the PV (as can be readily seen, for instance, from Ertel's theorem by taking for ψ any curvilinear abscissa along the vortex tube).

Zoom In Zoom Out Reset image size

Using Equations (49) and (55), we can recast the nonlinear corotation torque expression as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}=0.88{ \mathcal V }{{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{p}^{2}{r}_{p}^{4}{q}^{2}{h}_{p}^{-2}.\end{eqnarray} \tag{ 62 }$$

Assuming that to lowest order ${ \mathcal V }=3/2-\alpha$, we have

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}\approx (1.32-0.88\alpha ){{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{p}^{2}{r}_{p}^{4}{q}^{2}{h}_{p}^{-2}.\end{eqnarray} \tag{ 63 }$$

This expression is to be compared with the linear corotation torque of TTW02, which reads

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{CR}}^{\mathrm{lin}}\approx (0.976-0.64\alpha ){{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{p}^{2}{r}_{p}^{4}{q}^{2}{h}_{p}^{-2},\end{eqnarray} \tag{ 64 }$$

showing that the nonlinear corotation torque is ∼1.36 times larger than the linear corotation torque. This statement is reminiscent of the two-dimensional case where a similar ratio is found between the nonlinear and linear estimates of the corotation torques (Masset & Casoli 2010; Paardekooper et al. 2011).

The linear estimate of TTW02 vanishes for $\alpha =1.525$, slightly above the value 3/2 that one would naively expect. The vortensity gradient in a globally isothermal disk is actually not $3/2-\alpha$. Using Equation (84), we can show that it has the expression

$$\begin{eqnarray}&&{ \mathcal V }=\displaystyle \frac{3}{2}\left(1+\displaystyle \frac{3}{2}\;{h}^{2}\right)-\left(1-\displaystyle \frac{3}{2}\;{h}^{2}\right)\alpha ,\end{eqnarray} \tag{ 65 }$$

where we neglect terms of order h⁴ and above. The vortensity gradient therefore vanishes for $\alpha =3/2+(9/2){h}^{2}$ (which in our fiducial disk is $\alpha =1.511$) and so does the nonlinear corotation torque, which scales exactly with the vortensity gradient in our analysis.

One might have expected the width of the horseshoe region to decrease at higher altitudes, since the gas there is subjected to a reduced gravity from the planet, compared to the midplane. Similarly, one might have expected the width in the midplane to be similar to the width of the two-dimensional case in the limit of a vanishing softening length, since in this plane the motion is two-dimensional. This two-dimensional width, however, is much larger than the width reported here: instead of the value given by Equation (61), it is $\sim 1.7{r}_{p}\sqrt{q/h}$ (Paardekooper & Papaloizou 2009b; Ormel 2013). Therefore two questions arise: why is the horseshoe region so narrow at the midplane, and why is it so wide at large altitudes?

6.2.1. Width at the Midplane

Equation (60) shows that the width of the horseshoe region is determined by the value of the perturbed effective potential, given by Equation (59), at the stagnation point. The well of perturbed effective potential is much shallower than the gravitational potential well of the planet, as the latter is partially filled with enthalpy (which is why the horseshoe region is much narrower than predicted by the restricted three-body problem). In the 3D case, the density at the midplane differs from the density in a purely two-dimensional case with same velocity field because of the vertical motion above the midplane, which can compress the gas or allow it to expand. The difference in perturbed enthalpy between the three-dimensional case and a two-dimensional case with same resolution and parameters is shown in Figure 10. This difference is also the difference in perturbed effective potentials between the two cases, since the planetary potential is the same in the two- and three-dimensional cases. In the vicinity of the orbit, this difference is positive, resulting from the trend of the streamlines to bend toward the midplane. The effective potential well is therefore shallower in the 3D case than in the 2D case, and the horseshoe region is consequently narrower. The actual difference in the horseshoe width between the 2D and 3D cases cannot be deduced from Figure 10, because the stagnation point can be at a different location in each case, but the order of magnitude of the enthalpy excess in the 3D case is compatible with a sizable reduction in the horseshoe width with respect to the two-dimensional case.

Zoom In Zoom Out Reset image size

6.2.2. Width at High Altitude

This paper has focused on the unsaturated horseshoe drag. A detailed study of the saturation properties will be presented elsewhere. We can, however, anticipate here the following two points:

• 1.  
  As pointed out in Section 6.2.2, the widest horseshoe U-turns are performed more slowly at large altitude. In a two-dimensional situation, saturation occurs as a result of phase mixing because streamlines located at different distances of corotation have different libration times. In a three-dimensional situation, streamlines at the same distance to corotation but at different altitudes also have different libration times.⁸ We can therefore anticipate that the torque oscillations should exhibit a pattern different from those of two-dimensional disks (see, e.g., Masset & Casoli 2010), and that the torque should relax faster toward its asymptotic value as a result of enhanced phase mixing.
• 2.  
  In the evaluation of the saturated torque value, the parameter x_s (half-width of the horseshoe region) plays a very important role. The torque asymptotic (or saturated) value depends on a competition between libration (which tends to cancel the torque, as it flattens the gradients of vortensity and entropy across corotation) and diffusive processes (viscous and thermal diffusion, which tend to restore the large-scale gradients and the unsaturated torque value). Of these two processes, the one with the shortest timescale will win the competition. The ratio of the diffusive to the libration timescales is proportional to x_s³ (Masset 2001; Masset & Casoli 2010; Paardekooper et al. 2011), hence the ultimate torque value depends very sensitively on x_s. An accurate estimate of x_s such as the one we provide here is therefore crucially important in correctly evaluating the degree of saturation of the corotation torque. Masset & Casoli (2010) suggested that the width of the horseshoe region was expected to depend on the altitude, and that a value at low altitude, where most of the mass resides, could be used to determine the torque value. In view of the results reported here, this suggestion is clearly not adequate, and it leads, not surprisingly, to corotation torque values too saturated and therefore largely underestimated (Bitsch & Kley 2011). Regardless of the improvements required for three-dimensional torque formulae, we show in Appendix C that using for x_s the value given by Equation (61) improves considerably the torque estimate of Masset & Casoli (2010) and essentially reconciles it with other estimates (Paardekooper et al. 2011).

Recent work has highlighted the characteristics of the flow at the scale of Bondi's radius ${r}_{B}={{GM}}_{p}/{c}_{s}^{2}$ and below (Ormel et al. 2015b), and the potential impact of this flow on the horseshoe drag (Fung et al. 2015). Ormel et al. (2015b) show that the flow can enter the Bondi sphere at high altitude and exit near the midplane (in shear-dominated configurations) or enter the Bondi sphere near the midplane and exit at higher altitude (in headwind-dominated configurations). In the former case, Fung et al. (2015) have found, for a planet of slightly sub-thermal mass, that the material expelled near the midplane can reach and cross the horseshoe's boundary. They call this flow the transient horseshoe flow, as it involves fluid elements that participate in a unique close encounter with the planet. This flow destroys the horseshoe boundary (near the midplane), invalidating our assumption of the existence of a critical surface: streamlines of different origins merge near the edge of the downstream horseshoe flow. Our work completely disregards these potentially very important effects, owing to its moderate resolution. As noted by Fung et al. (2015), the huge resolution needed to resolve the flow within the Bondi sphere for planetary masses as low as those considered here is beyond computational tractability on present-day platforms, at least for single-mesh calculations, even with non-uniform resolution. The problem of the impact of the flow at sub-Bondi scales on the horseshoe drag should be tackled by means of nested-mesh calculations, which is a forthcoming feature of the code that we used in this paper. Fung et al. (2015) find that the total torque in their calculation is still negative, but reduced in magnitude with respect to the expected value of the linear Lindblad torque. They perform a horseshoe drag analysis of their outcome by summing the mass flow rate, weighted by the angular momentum jump (which they determine by a streamline integration), over the whole set of horseshoe streamlines. They find that the contribution of the transient horseshoe flow to the torque is negligible, and that there exists a net positive corotation torque due to the standard horseshoe flow, in spite of the disk being barotropic with a vanishing vortensity gradient. It is noteworthy that they find a non-vanishing contribution at all altitudes (see their Figure 15), even above the first pressure scale height, where the flow most resembles the situation that we described here, with a columnar structure and well-defined separatrix sheets. There are ingredients neglected in our analysis that could play an important role, such as the radial tilt of vortex tubes upon U-turns, or the non-conservation of Bernoulli's constant if shocks are present in the planet's vicinity. The impact of the flow at sub-Bondi scale on the planetary torque clearly requires further work, in which such ingredients should probably be incorporated.

Extensions of the present work could also consider the impact of the radial temperature gradient (in so-called locally isothermal disks) or the role of the entropy gradient in flows with an energy equation (in so-called adiabatic disks). The entropy gradient is known to have a major impact on the corotation torque in three-dimensional disks (Paardekooper & Mellema 2008; Lega et al. 2015), but so far its action on the co-orbital flow and the generation mechanism of the entropy-related torque have been studied only in two-dimensional disks (Masset & Casoli 2009; Paardekooper et al. 2010).