Impact of thermal effects on the evolution of eccentricity and inclination of low-mass planets

TL;DR

This paper investigates how thermal effects influence the evolution of eccentricity and inclination of low-mass planets in protoplanetary discs, revealing conditions for damping or exponential growth based on planetary luminosity.

Contribution

It introduces a linear perturbation theory framework to analyze thermal effects on planetary orbital evolution, highlighting the impact of planetary luminosity on eccentricity and inclination.

Findings

Eccentricity and inclination are damped in non-luminous planets.

High planetary luminosity causes exponential growth of eccentricity and inclination.

Growth rate of eccentricity exceeds that of inclination at large luminosities.

Abstract

Using linear perturbation theory, we evaluate the time-dependent force exerted on an eccentric and inclined low-mass planet embedded in a gaseous protoplanetary disc with finite thermal diffusivity $\chi$. We assume the eccentricity and inclination to be small compared to the size of the thermal lobes $\lambda\sim(\chi/\Omega)^{1/2}$, itself generally much smaller than the scalelength of pressure $H$. When the planet is non-luminous, we find that its eccentricity and inclination are vigorously damped by the disc, over a timescale shorter by a factor $H/\lambda$ than the damping timescale in adiabatic discs. On the contrary, when the luminosity-to-mass ratio of the planet exceeds a threshold that depends on the disc's properties, its eccentricity and inclination undergo an exponential growth. In the limit of a large luminosity, the growth rate of the eccentricity is 2.5~times larger than that of the inclination, in agreement with previous numerical work. Depending on their luminosity, planetary embryos therefore exhibit much more diverse behaviours than the mild damping of eccentricity and inclination considered hitherto.

Equations

$$(x_p,y_p,z_p)=\left[x_p^0-E\cos ({\Omega }_{p}t),2E\sin ({\Omega }_{p}t),I\sin ({\Omega }_p{t}^{\prime })\right],$$

$${\partial }_t\rho +\nabla .\left(\rho \mathbf{V}\right)=0$$

$${\partial }_t\mathbf{V}+\mathbf{V}.\nabla \mathbf{V}+2{\Omega }_p{\mathbf{e}}_z\times \mathbf{V}=-\nabla ({\Phi }_t+{\Phi }_p)-\frac{\nabla p}{\rho },$$

$${\Phi }_t=-\frac{3}{2}{\Omega }_p^2(x-x_p^0{)}^2+\frac{{\Omega }_p^2{z}^2}{2}.$$

$${\partial }_t\mathfrak{e}+\nabla .\left(\mathfrak{e}\mathbf{V}\right)=-p\nabla .\mathbf{V}-\nabla .{\mathbf{F}}_H+S,$$

$${\mathbf{F}}_H=-\rho \chi \nabla \left(\frac{\mathfrak{e}}{\rho }\right),$$

$$\rho$$

$$\mathfrak{e}$$

$$p$$

$$u$$

$$v$$

$$w$$

$$x_p^0=-\frac{{\partial }_x{p}_0}{3{\Omega }_p^2{\rho }_0}.$$

$$E\ll \lambda \text{mbox}andI\ll \lambda .$$

$$p=(\gamma -1)\mathfrak{e},$$

$$c_s=\gamma p_0/{\rho }_0$$

$${\partial }_t{\rho }^{\prime }-\frac{3}{2}{\Omega }_{p}x{\partial }_y{\rho }^{\prime }+{\rho }_0({\partial }_x{u}^{\prime }+{\partial }_y{v}^{\prime }+{\partial }_z{w}^{\prime })=0,$$

$${\partial }_t{u}^{\prime }-\frac{3}{2}{\Omega }_{p}x{\partial }_y{u}^{\prime }-2{\Omega }_p{v}^{\prime }=-{\partial }_x{\Phi }_p-\frac{{\partial }_x{p}^{\prime }}{{\rho }_0}+\frac{({\partial }_x{p}_0){\rho }^{\prime }}{{\rho }_0^2},$$

$${\partial }_t{v}^{\prime }-\frac{3}{2}{\Omega }_{p}x{\partial }_y{v}^{\prime }+\frac{1}{2}{\Omega }_p{u}^{\prime }=-{\partial }_y{\Phi }_p-\frac{{\partial }_y{p}^{\prime }}{{\rho }_0},$$

$${\partial }_t{w}^{\prime }-\frac{3}{2}{\Omega }_{p}x{\partial }_y{w}^{\prime }=-{\partial }_z{\Phi }_p-\frac{{\partial }_z{p}^{\prime }}{{\rho }_0}+\frac{({\partial }_z{p}_0){\rho }^{\prime }}{{\rho }_0^2},$$

$$\begin{array}{r}{\partial }_t{p}^{\prime }-\frac{3}{2}{\Omega }_{p}x{\partial }_y{p}^{\prime }+\gamma p_0({\partial }_x{u}^{\prime }+{\partial }_y{v}^{\prime }+{\partial }_z{w}^{\prime })=\\ \chi \Delta p^{\prime }-\chi \frac{p_0}{{\rho }_0}\Delta {\rho }^{\prime }+(\gamma -1)S_p,\end{array}$$

$$\tilde{\xi }(x,k_y,k_z)=\int \int {\xi }^{\prime }(x,y,z)e^{-j(k_{y}y+k_{z}z)}dydz,$$

$${\xi }^{\prime }(x,y,z)=\frac{1}{4{\pi }^2}\int \int \tilde{\xi }(x,k_y,k_z)e^{j(k_{y}y+k_{z}z)}d{k}_{y}d{k}_z,$$

$${\partial }_x$$

$${\partial }_y$$

$${\partial }_z$$

$$\nabla$$

$$\Delta$$

$${\partial }_t\tilde{\rho }-j\frac{3}{2}{k}_y{\Omega }_{p}x\tilde{\rho }+{\rho }_0\tilde{\nabla }.\tilde{\mathbf{V}}=0$$

$${\partial }_t\tilde{u}-j\frac{3}{2}{k}_y{\Omega }_{p}x\tilde{u}-2{\Omega }_p\tilde{v}+\frac{{\partial }_x\tilde{p}}{{\rho }_0}=-{\partial }_x{\tilde{\Phi }}_p$$

$${\partial }_t\tilde{v}-j\frac{3}{2}{k}_y{\Omega }_{p}x\tilde{v}+\frac{1}{2}{\Omega }_p\tilde{u}+\frac{j{k}_y\tilde{p}}{{\rho }_0}=-j{k}_y{\tilde{\Phi }}_p$$

$${\partial }_t\tilde{w}-j\frac{3}{2}{k}_y{\Omega }_{p}x\tilde{w}+\frac{j{k}_z\tilde{p}}{{\rho }_0}=-j{k}_z{\tilde{\Phi }}_p$$

$${\partial }_t\tilde{p}-j\frac{3}{2}{k}_y{\Omega }_{p}x\tilde{p}+\gamma p_0\tilde{\nabla }.\tilde{\mathbf{V}}-\chi {\Delta }^{\prime }\tilde{p}+\frac{\chi c_s^2}{\gamma }{\Delta }^{\prime }\tilde{\rho }=(\gamma -1){\tilde{S}}_p.$$

$$\left({\partial }_t-j\frac{3}{2}{k}_y{\Omega }_{p}x\right)(\tilde{p}-c_s^2\tilde{\rho })-\chi {\Delta }^{\prime }\left(\tilde{p}-\frac{c_s^2}{\gamma }\tilde{\rho }\right)=(\gamma -1){\tilde{S}}_p.$$

$$\mathcal{K}(\tilde{\rho })=\mathcal{L}\left({\tilde{\Phi }}_p+\frac{\tilde{p}}{{\rho }_0}\right),$$

$$\mathcal{K}(\tilde{\rho })=\frac{\tilde{\rho }}{{\rho }_0}$$

$$\mathcal{L}(Y)=-\frac{{\partial }_x^{2}Y}{{\Omega }_p^2{D}_2}+\frac{3k_y\tilde{\omega }}{{\Omega }_p^3{D}_2^2}{\partial }_{x}Y+\left[\frac{k_y^2}{{\Omega }_p^2{D}_2}\left(1-\frac{6}{D_2}\right)+\frac{k_z^2}{{\tilde{\omega }}^2}\right]Y,$$

Full text

††pubyear: 2019††pagerange: Impact of thermal effects on the evolution of eccentricity and inclination of low-mass planets–References

1 Introduction

Most analytic studies of planet-disc interactions have long been limited to the barotropic approximation. The inclusion of thermodynamics was introduced more than a decade ago, but was mainly restricted to non-linear effects tackled through numerical simulations, with a focus either on giant planets (D’Angelo et al., 2003) or on the non-linear dynamics of the corotation torque (Paardekooper & Mellema, 2006). The inclusion of non-barotropic effects in a linear analysis has been worked out by Tsang (2014), who performed a study of the corotation torque in a non-barotropic disc. Studies of the role of thermal diffusion itself, however, was even longer restricted to non-linear effects, especially in the modelling of the saturation of the corotation torque of intermediate-mass planets (above a few Earth masses), either through numerical simulations or toy models of the horseshoe dynamics (Masset & Casoli, 2010; Paardekooper et al., 2011; Jiménez & Masset, 2017). The impact of thermal diffusion on the interaction between the disc and low-mass planets using linear perturbation theory has been studied more recently. It was firstly noticed in numerical simulations of non-luminous111Throughout this work we use this expression for planets that do not inject heat into the surrounding gas., low-mass planets embedded in radiative222Discs in which thermal diffusion is effected by radiative transfer. discs by Lega et al. (2014), who found that the torque on a low-mass planet on a circular orbit can be significantly more negative when thermal diffusion is included than when it is not. This effect, dubbed cold finger by these authors, was later described by Masset (2017) who provided an analytic expression for the corresponding torque component. Little is known, however, on the impact of thermal diffusion on the gravitational interaction between an eccentric or inclined non-luminous planet and a gaseous disc. Eklund & Masset (2017) performed numerical experiments in isothermal and radiative discs, and found the eccentricity and inclination damping to be much stronger in the latter than in the former, but this effect, which was not the primary focus of that work, was not systematically quantified, and may have been misrepresented by the low resolution with which it was captured.

While the release of heat into the ambient gas by accreting planets has been studied for more than two decades (e.g. Pollack et al., 1996), such studies considered the planet to be at rest in a uniform medium in order to use 1D grids, and the feedback of heat release on the planetary orbit was ignored. It was until recently that the heat release was incorporated, in a highly simplified manner, to numerical simulations of planet-disc interactions. Benítez-Llambay et al. (2015) have found that luminous planetary embryos in the Earth-mass range may undergo outwards migration if their luminosity is above a threshold that should be easily overcome if they are subjected to fast pebble accretion. The effect at the origin of this outward migration shares many similarities with the “cold finger” effect of Lega et al. (2014), to the point that both effects can be unified into a single description (Masset, 2017). Eklund & Masset (2017) performed a follow-up study of the work of Benítez-Llambay et al. (2015), by relaxing the constraint of a circular and coplanar orbit. They found that planetary embryos, if sufficiently luminous, undergo a growth of eccentricity and inclination. A similar result holds in 2D calculations: Chrenko et al. (2017) found that luminous embryos embedded in 2D discs experience an eccentricity growth. They called this effect the hot trail effect. The disc’s response to the heat release indeed adopts a trailing, cometary shape for eccentricities well below the disc’s aspect ratio (Eklund & Masset, 2017). Unlike the response carried by density waves, the effect of heat release can therefore be captured by a simple calculation of dynamical friction even in the subsonic regime (Masset & Velasco Romero, 2017). Note however that, again, the net force arising from thermal effects in the case of dynamical friction is positive (i.e., a thrust) if the planet is sufficiently luminous (Velasco Romero & Masset, 2019).

From the above we see that thermal effects on low-mass planets have been studied both using linear perturbation theory and through numerical simulations for planets in circular orbits and for perturbers in unsheared, homogeneous media, adequate for the description of planets with a sizeable eccentricity or inclination. However, there has not been any analytical study of the regime of small eccentricities and inclinations, when the disc response is not that of a simple “cometary” trail captured by a dynamical friction calculation. The purpose of this work is to provide this missing part, so as to give analytic expressions for the excitation or damping of the eccentricity and inclination of a low-mass planet that can be used in future models of planetary formation, and to shed some light on the different behaviours observed so far in numerical experiments. We present our governing equations in section 2 and work out the density response in section 3. The force arising from thermal effects is then worked out in section 4, and the time evolution of the planet’s orbital parameters is derived in section 5. We finally discuss our results in section 6 and draw our conclusions in section 7.

2 Governing equations

We consider a planet of mass $M$ embedded in a protoplanetary gaseous disc on a prograde, slightly eccentric and slightly inclined orbit. The central star has a mass $M_{\star}$, the disc has a surface density $\Sigma(r)$ and an angular velocity $\Omega(r)$, where $r$ is the distance to the central star. We make the assumption that the disturbances due to thermal effects are small compared to the pressure lengthscale $H$ of the disc. This assumption has been used and discussed by Masset (2017). It allows us to perform our analysis in a 3D shearing box (Narayan et al., 1987). Our frame is therefore essentially a Cartesian box of dimensions much smaller than the planet’s semi-major axis, which contains the planet and co-rotates with its guiding centre. We use the conventional notation for the axes: $x$ is directed along the gradient of unperturbed velocity, $y$ is directed along the unperturbed motion and $z$ is perpendicular to the disc’s midplane. We will refer to the material at $x>0$ ($x<0$) as the outer (inner) disc, implying that the central object lies on the negative side of the $x$-axis. We denote with $E$ the epicyclic excursion of the planet: $E=ea$, where $e$ is the eccentricity and $a$ the semi-major axis. Similarly, we denote with $I$ its vertical excursion, and have $I=ia$, where $i$ is the inclination. We note $\Omega_{p}$ the planet’s orbital frequency and restrict ourselves to the Keplerian case, for which the epicyclic frequency matches the orbital frequency. In our frame, the planet location is therefore:

[TABLE]

where $x_{p}^{0}$ is the offset between the planet’s guiding centre and its corotation, $t$ is the time measured from a passage at periastron and $t^{\prime}=t-t_{\mathrm{AN}}$, where $t_{\mathrm{AN}}$ is a time of passage at the ascending node.

Our governing equations are the continuity equation, Euler’s equations, and the equation on the internal energy density. The continuity equation is:

[TABLE]

where $\rho$ and $\boldsymbol{V}=(u,v,w)^{T}$ are respectively the density and the velocity of the gas. The Euler equation reads:

[TABLE]

where $\Omega_{p}\boldsymbol{e}_{z}$ is rotation rate of the frame, $\Phi_{p}$ is the planetary potential, $p$ the pressure and $\Phi_{t}$ is the tidal potential given by:

[TABLE]

Finally, the equation for the density of internal energy, which we denote with $\mathfrak{e}$ in order to avoid confusion with the eccentricity, reads:

[TABLE]

where $S=S_{d}(\boldsymbol{r})+S_{p}(\boldsymbol{r})$ is a source term arising from the disc local heating $S_{d}(\boldsymbol{r})$ and from the release of energy into the gas by the planet $S_{p}(\boldsymbol{r})$. In what follows we assume the disc to be inviscid and neglect the perturbation of $S_{d}(\boldsymbol{r})$. We discuss in Appendix B how to extend our results to laminar, viscous discs. The heat flux $\boldsymbol{F}_{H}$ is given by:

[TABLE]

where $\chi$ is the thermal diffusivity.

Hereafter we consider the perturbations of density $(\rho^{\prime})$, velocities $(u^{\prime},v^{\prime},w^{\prime})$ and pressure ($p^{\prime}$) and we express the perturbed quantities as the sum of the unperturbed quantities (with a [math] subscript) and the corresponding perturbations (primed) as:

[TABLE]

The unperturbed velocity in Eq. (11) corresponds to the Keplerian shear and comes from Eq. (3) for the unperturbed quantities, which also yields:

[TABLE]

The offset between the guiding centre and corotation depends on the pressure gradient. In what follows we assume $x_{p}^{0}\ll E$ and $x_{p}^{0}\ll I$. We discuss how our results are affected when we relax this hypothesis in appendix A. We also assume that the epicyclic and vertical excursions $E$ and $I$ are much smaller than the characteristic size $\lambda$ of the disturbance:

[TABLE]

We consider the gas to be ideal. It obeys the relationship:

[TABLE]

where $\gamma$ is the adiabatic index. We denote with $c_{s}$ the adiabatic sound speed, given by:

[TABLE]

Upon linearization in $\rho^{\prime}$, $u^{\prime}$, $v^{\prime}$, $w^{\prime}$ and $\mathfrak{e}^{\prime}$, the governing equations can be recast as follows. The continuity equation reads:

[TABLE]

the three components of the Euler equation are:

[TABLE]

[TABLE]

[TABLE]

and the energy equation becomes:

[TABLE]

where $\Delta\equiv\partial^{2}_{x^{2}}+\partial^{2}_{y^{2}}+\partial^{2}_{z^{2}}$ is the Laplacian operator. Following Masset (2017), we neglect the last term of the right-hand side (R.H.S.) of Eq. (18) and the third term of the R.H.S. of Eq. (20), as the size of the perturbation is assumed small compared to the pressure lengthscale.

We perform a Fourier transform in the $y$ and $z$ directions with the following convention in sign and normalization 333Throughout this work, we denote with $j$ the imaginary number with positive imaginary part such that $j^{2}=-1$, so as to avoid confusion with the inclination.:

[TABLE]

where $\xi^{\prime}$ is the perturbation of any variable and $\tilde{\xi}$ its Fourier transform. The derivative operators become:

[TABLE]

and Eqs. (17)-(21) can be rewritten as follows:

[TABLE]

Using Eq. (29) to substitute the velocity divergence ($\tilde{\boldsymbol{\nabla}}.\tilde{\boldsymbol{V}}$) in Eq. (33) we get:

[TABLE]

This is our main equation. We work out a significant simplification in the following section.

3 Density response

3.1 Simplification of the main equation

In order to give an order of magnitude of the different terms that feature in the density response, we take the Fourier transform in time of the perturbation of the different quantities, and substitute the partial derivative with respect to time in Eqs. (29) to (32) by a multiplication by $-j\omega$. Using Eqs. (30)-(32), we can write the three components of the velocity in terms of $\tilde{\Phi}_{p}+\tilde{p}/\rho_{0}$. Upon substitution in Eq. (29), we are led to:

[TABLE]

where the linear operators ${\cal K}$ and ${\cal L}$ are defined respectively by:

[TABLE]

and

[TABLE]

where $D_{2}$ is the non-dimensional quantity:

[TABLE]

and $\tilde{\omega}$ is the frequency in the frame of matter:

[TABLE]

Using Eq. (13), the linear operator $\mathcal{K}(\tilde{\rho})$ can be recast as:

[TABLE]

We now specify to the case of a perturbation not triggered by a gravitational potential, such as the one arising from the release of heat by the planet. In that case:

[TABLE]

We can work out the order of magnitudes of the different terms of this identity using the following approximations: $|x|\sim|k_{y}^{-1}|\sim|k_{z}^{-1}|\sim\lambda\ll H$, $|\partial_{x}\tilde{p}|\sim|\tilde{p}|/H$, $|\partial_{x^{2}}^{2}\tilde{p}|\sim|\tilde{p}|/H^{2}$. We then consider three cases:

Over most of the perturbation, we have $\tilde{\omega}/\Omega_{p}=O(1)$ and $D_{2}=O(1)$. In that case the dominant term of Eq. (37) is the last one that has the order of magnitude $|\tilde{p}|/[\rho_{0}(\Omega_{p}^{2}\lambda^{2})]$, while the dominant term of Eq. (40) is the first one. We then have:

[TABLE]

hence $\tilde{p}$ is negligible compared to $\tilde{\rho}c_{s}^{2}$. 2. 2.

Whenever $\Omega_{p}/\tilde{\omega}\gg 1$, the dominant term in Eq. (37) is the last one and has order of magnitude $\tilde{p}/(\rho_{0}\lambda^{2}\tilde{\omega}^{2})$, whereas the dominant term of Eq. (40) is $\mathrm{max}[1,(x_{p}^{0}/\lambda\times\Omega_{p}/\tilde{\omega}]\tilde{\rho}/\rho_{0}$, hence:

[TABLE] 3. 3.

Whenever $D_{2}\ll 1$, the dominant term in Eq. (37) is $\tilde{p}/(\lambda\Omega_{p}D_{2})^{2}$ and the dominant term in Eq. (40) is $\mathrm{max}[1,x_{p}^{0}/(\lambda D_{2}^{2})]\tilde{\rho}/\rho_{0}$, therefore:

[TABLE]

The relation $|\tilde{p}|\ll c_{s}^{2}|\tilde{\rho}|$ is therefore verified by any disturbance not triggered by a gravitational potential, i.e. that verifies Eq. (41), when the size of the disturbance is much smaller than the pressure lengthscale.

3.2 Net thermal effects

Hereafter we show that the thermal can be decomposed into two contributions: one that arises from the inclusion of thermal diffusion, even if the planet is non-luminous, and another one that arises from the planet’s luminosity itself. For this purpose we perform the following decomposition of the perturbation:

[TABLE]

where the $a$ subscript refers to the adiabatic solution, that of Eqs. (29) to (33) with $\chi=0$ and $\tilde{S}_{p}=0$, while the quantities without subscript refer to the general solution of these equations (i.e. with $\chi\neq 0$ and $\tilde{S}_{p}\neq 0$). In Eqs. (45) and (46), the quantities with a $t$ subscript, which stands as thermal, therefore appear as the difference between the adiabatic solution and the solution with heat release ($\tilde{S}_{p}\neq 0$) and thermal diffusion ($\chi\neq 0$). Since the thermal diffusivity neither appears in the operator $\mathcal{K}$ nor $\mathcal{L}$, we have:

[TABLE]

and we can write, subtracting Eq. (47) from Eq. (35) and using the linearity of these operators:

[TABLE]

The size of the thermal disturbance being assumed small compared to the pressure lengthscale, we then have, as shown above:

[TABLE]

Noting that Eq. (33) reduces to

[TABLE]

we recast Eq. (34) using the decomposition of Eqs. (45) and (46):

[TABLE]

where we have used Eq. (49) to get rid of the instances of $\tilde{p}_{t}$. The second term of the R.H.S. appears as an additional source term. In order to proceed, we make the same approximation as Masset (2017) and write:

[TABLE]

This approximation is valid for low-mass planets ($M/M_{\star}\ll(H/r)^{3}$) which have a low, largely subsonic velocity with respect to the ambient gas. Assuming that the relative variations of $\rho_{0}$ are much smaller than those of $\Phi_{p}$ over the perturbed region, we have, going back to real space:

[TABLE]

Using Poisson’s equation for the planetary potential:

[TABLE]

and specifying from now on to a singular planetary heating term:

[TABLE]

where $L$ is the planet’s luminosity, we can recast Eq. (53) as

[TABLE]

where the modified source term $S^{\mathrm{tot}}_{p}$ is given by:

[TABLE]

where $L_{c}$ has the value:

[TABLE]

In what follows, Eq. (56) or its Fourier transform in $y$ and $z$, which reads:

[TABLE]

constitutes our main equation.

3.3 Evaluation of the source terms

The source term is given by:

[TABLE]

As said in section 2, we neglect from now on $x_{p}^{0}$ with respect to $E$ and $I$. Taking the Fourier transform in $y$ and $z$ of Eq. (60), we arrive at:

[TABLE]

which reads, keeping only the terms up to first order in $E$ and $I$ as per our hypothesis of Eq. (14):

[TABLE]

We recognize the source term of the circular, coplanar case ($E=I=0$) in the first term of the R.H.S. The two subsequent terms arise from the eccentricity and the last one from inclination. If we define:

[TABLE]

and

[TABLE]

then we have:

[TABLE]

where the “circular” term $S_{c}(x,k_{y},k_{z})$ is $(L-L_{c})\delta(x)$. Since Eq. (59) is linear, we have:

[TABLE]

where $\tilde{\rho}_{c}$ ($\tilde{\rho}_{e}$, $\tilde{\rho}_{i}$) is the solution of Eq. (59) where the source term $\tilde{S}^{\mathrm{tot}}_{p}$ of the R.H.S. is replaced by $\tilde{S}_{c}$ ($\tilde{S}_{e}$, $\tilde{S}_{i}$). The term $\tilde{\rho}_{c}$ is the response to $\tilde{S}_{c}$ and has been studied by Masset (2017). We focus in the following on the response to the “eccentric” and “inclination” source terms $\tilde{S}_{e}$ and $\tilde{S}_{i}$. We note that in real space Eq. (66) is simply:

[TABLE]

and we mention that $\rho_{e}^{\prime}$ and $\rho_{i}^{\prime}$ are not dimensionally homogeneous to densities, as the factors $E$ and $I$ have unit of length.

These source terms can be expressed in a slightly different manner. If we define:

[TABLE]

and

[TABLE]

then the source terms can be written respectively as:

[TABLE]

and

[TABLE]

We denote with $\tilde{\rho}_{e}^{\pm}\exp(\pm j\Omega_{p}t)$ the solution to Eq. (59) where the source term of the R.H.S. is replaced by $S_{e}^{\pm}\exp(\pm j\Omega_{p}t)$. Similarly, we denote with $\tilde{\rho}_{i}^{\pm}\exp(\pm j\Omega_{p}t^{\prime})$ the solution of the same equation where the source term of the R.H.S is replaced by $S_{i}^{\pm}\exp(\pm j\Omega_{p}t^{\prime})$. We therefore have:

[TABLE]

and

[TABLE]

Noting that $S_{e}^{+}(x)=-S_{e}^{-}(-x)$, it is straightforward to show that:

[TABLE]

where the $*$ superscript denotes the complex conjugate. We therefore focus only on the solution of Eq. (72) and will use Eq. (74) whenever the expression of $\tilde{\rho}_{e}^{-}$ is needed. In a similar manner, we have:

[TABLE]

and we use the relationship

[TABLE]

to infer $\tilde{\rho}_{i}^{-}$ whenever needed.

3.4 Non-dimensional form of the main equation

We write Eqs. (72) and (75) in non-dimensional form. For that purpose we introduce new variables:

[TABLE]

[TABLE]

[TABLE]

where $k=(k_{y}^{2}+k_{z}^{2})^{1/2}$. This allows us to rewrite Eq. (72) as:

[TABLE]

where $\tilde{\rho}_{e}^{+^{\prime\prime}}$ is the second derivative of $\tilde{\rho}_{e}^{+}$ with respect to $X$ and where:

[TABLE]

Similarly, we have:

[TABLE]

where:

[TABLE]

We call $T_{X_{c},K}$ ($W_{X_{c},K}$) the solution of Eq. (80) [Eq. (82)] with $s_{1}=1$ ($s_{2}=1$) whose real and imaginary parts tend to zero in $\pm\infty$, and we have:

[TABLE]

and

[TABLE]

We describe hereafter how we obtain the solutions $T_{X_{c},K}$ and $W_{X_{c},K}$.

3.5 Jump conditions at origin

From the considerations above, the function $T(X;X_{c},K):\mathbb{R}\rightarrow\mathbb{C}$ verifies:

[TABLE]

where for the sake of clarity we omit the full dependence of $T$. We work out the jump conditions that $T$ must fulfil in $X=0$. We assume $T$ is finite (but possibly discontinuous). We call $F_{\varepsilon}(X)$ the following integral of $j(X_{c}-X)T+KT$, where $\varepsilon>0$:

[TABLE]

Since the integrand is finite, $F$ is a continuous function of $X$ or $\varepsilon$ in zero. We have, using Eq. (86):

[TABLE]

where $H(X)$ is Heaviside’s step function. We now define:

[TABLE]

which, like $F$, is a continuous function of $X$ or $\varepsilon$ in zero. We have:

[TABLE]

Evaluating Eq. (88) in $X=\varepsilon$ and letting $\varepsilon\rightarrow 0$, we get:

[TABLE]

while evaluating Eq. (90) in $X=\varepsilon$ and letting $\varepsilon\rightarrow 0$, we get:

[TABLE]

Eqs. (91) and (92) are the jump conditions in $X=0$. As their R.H.S. are real, the imaginary part of the solution and of its derivative are continuous in $X=0$, whereas their real parts undergo the jumps above. When $X_{c}=0$, Eq. (86) reduces to:

[TABLE]

which is the equation of the circular problem, and has solution $R_{K}+jI_{K}$ (see Masset, 2017, Eqs. (75)–(79)). Similarly, the function $W(X;X_{c},K):\mathbb{R}\rightarrow\mathbb{C}$ verifies the differential equation:

[TABLE]

and can be easily shown to be continuous in $X=0$, while its first derivative has the following jump:

[TABLE]

3.6 Numerical procedure

We use a shooting method with a Runge-Kutta scheme at fifth order and start our integration at a large distance $\pm X_{0}$ from the planet. We integrate towards corotation ($X=0$), so that the integration started in $-X_{0}$ goes forward and that started in $X_{0}$ goes backwards. Our boundary conditions in $\pm X_{0}$ consist of the real and imaginary parts of the solution. Its first-order derivative can be chosen as in Masset (2017) (Eqs. A3 and A4), but we find that our solution is insensitive to this choice (provided $X_{0}$ is large enough), and simpler choices, such as a null derivative, yield virtually the same solution. In the following we present the method employed to obtain the solution $T_{X_{c},K}$. It can be applied straightforwardly to $W_{X_{c},K}$ by amending the jump conditions in $X=0$. We denote with $\boldsymbol{I}=\{\Re[T(X_{0})],\Im[T(X_{0})],\Re[T(-X_{0})],\Im[T(-X_{0})]\}^{T}$ our input vector. Upon integration, we obtain the jumps in $X=0$ corresponding to our choice of $\boldsymbol{I}$. These jumps constitute our output vector:

[TABLE]

where, here only, the $\Delta$ symbol represents the difference in $X=0$ between the values numerically obtained by the backwards integration with $X\geq 0$ and the forward integration with $X\leq 0$. Our aim is that the output vector matches the jump conditions, i.e.

[TABLE]

with:

[TABLE]

as required by Eqs. (91) and (92) and the fact that the imaginary part of the solution and its first derivative is continuous in $X=0$. Since Eq. (86) is linear in $T$ over the domains $X>0$ and $X<0$, the output vector is a linear function of the input vector, i.e. there exists a $4\times 4$ matrix ${\cal M}$ with real coefficients such that:

[TABLE]

We construct the matrix ${\cal M}$ column by column, using successively $\boldsymbol{I}=(1,0,0,0)^{T}$, $\boldsymbol{I}=(0,1,0,0)^{T}$, $\boldsymbol{I}=(0,0,1,0)^{T}$ and $\boldsymbol{I}=(0,0,0,1)^{T}$. The first two cases correspond to $T(-X_{0})=0$ and $T(X_{0})=1$ or $T(X_{0})=j$, respectively, while the last two cases correspond to $T(X_{0})=0$ and $T(-X_{0})=1$ or $T(-X_{0})=j$, respectively. The output vectors give directly the columns of the matrix ${\cal M}$. Once the latter is built, we determine the boundary conditions $\boldsymbol{I}_{T}$ for the sought solution $T$ using:

[TABLE]

The solution $T$ is then constructed by a fifth and last integration using the boundary conditions $\boldsymbol{I}_{T}$. It fulfils the jump conditions in $X=0$ by construction. The solution $W_{X_{c},K}$ must verify the following differential equation:

[TABLE]

It is constructed in almost the same manner, except that we use different jump conditions given by:

[TABLE]

4 Force expression

Having determined the density response, we can now calculate the force exerted on the planet by the perturbed density. Its three components read:

[TABLE]

The gravitational potential of the planet is:

[TABLE]

which reads at first order in $E$ and $I$:

[TABLE]

with

[TABLE]

where $R=(x^{2}+y^{2}+z^{2})^{1/2}$. In what follows we first work out a separation of the terms into those arising from the eccentricity, which give rise to the horizontal components of the force, and those arising from the inclination, which give rise to the vertical component of the force, then we carry out the calculation of the force components.

4.1 Separation of eccentricity and inclination contributions

Substituting Eqs. (67) and (107) in Eq. (103), we are led to:

[TABLE]

Expanding this expression to first order in $E$ and $I$, we are left with:

[TABLE]

Since Eq. (59) does not include partial derivatives in $z$ of odd order, the solution has the same parity in $z$ as the source term. Since $S_{e}$ is even in $z$, so is $\rho^{\prime}_{e}$, and since $S_{i}$ is odd in $z$, so is $\rho^{\prime}_{i}$. Furthermore, Eqs. (109) and (110) show respectively that $\Phi_{e}$ is even in $z$ and $\Phi_{i}$ is odd in $z$. This implies that the part of the integrand in factor of $I$ in Eq. (112) vanishes: the force in the $x$ direction depends only on the eccentricity. The same is true of the force in the $y$-direction. On the other hand, since the derivative in $z$ changes the parity in $z$ of the function to which it is applied, the component $F_{z}$ of the force reduces to:

[TABLE]

The vertical component of the force therefore only depends on the inclination. In this section we have used the expansion of the planetary potential in real space to separate the force expression into contributions arising from the eccentricity or the inclination. In what follows, we evaluate the force using the expression of the potential in Fourier space, using for that purpose Parseval-Plancherel’s theorem, which reads, with the conventions of Eqs. (22) and (23):

[TABLE]

4.2 Force component in the $x$-direction

Using Eqs. (112) and (114) we can write the $x$-component of the force as:

[TABLE]

with

[TABLE]

and

[TABLE]

We have omitted the constant term $\rho_{c}\partial_{x}\Phi_{c}$, which vanishes in the shearing sheet for symmetry reasons. The expression of $\tilde{\Phi}_{c}$ is (Masset, 2017, Eq. 40):

[TABLE]

The expression of the potential in the eccentric case can be inferred from Eq. (118) through a shift in $x$ and $y$:

[TABLE]

Using our hypothesis of Eq. (14), we expand this expression to first order in $E$. This yields:

[TABLE]

where the dependence on $x$, $k_{y}$ and $k_{z}$ has been omitted for improved legibility, and where:

[TABLE]

and

[TABLE]

Using Eqs. (116) and (118), we have:

[TABLE]

Using Eq. (74), and performing the change of variable $x\rightarrow-x$ on one of the terms of the integrand, this expression can be recast as:

[TABLE]

We start the evaluation of $F_{x}^{(2)}$ using an integration by parts, which leads, using Eq. (117):

[TABLE]

We note that $\tilde{\Phi}_{e}^{+}(-x)=-\tilde{\Phi}_{e}^{-}(x)$ and that (Masset, 2017, Eq. 116 and Appendix A):

[TABLE]

so that we can use a change of variable similar to that used in evaluating $F_{x}^{(1)}$. This yields:

[TABLE]

Using Eqs. (115), (121), (124) and (127), we can write:

[TABLE]

with:

[TABLE]

and:

[TABLE]

where for a more concise notation we define $\tilde{\rho}_{e}^{R}\equiv\Re(\tilde{\rho}_{e}^{+})$, $\tilde{\rho}_{c}^{R}\equiv\Re(\tilde{\rho}_{c})$, $\tilde{\rho}_{e}^{I}\equiv\Im(\tilde{\rho}_{e}^{+})$ and $\tilde{\rho}_{c}^{I}\equiv\Im(\tilde{\rho}_{c})$. An integration by parts has been used to write the second term of the integrand of Eq. (130).

4.3 Force component in the $y$-direction

The $y$-component of the force is calculated following an approach very similar to that of section 4.2. It is expressed as:

[TABLE]

with

[TABLE]

and

[TABLE]

Eq. (132) is transformed in a similar way as the expressions of $F_{x}^{(1)}$ and $F_{x}^{(2)}$. We use the fact that $\tilde{\Phi}_{c}$ is even in $x$, and make use again of Eq. (74) to write:

[TABLE]

Using Eqs. (121) and (122), we can rewrite Eq. (133) as:

[TABLE]

Since $\tilde{\rho}_{c}^{R}$ ($\tilde{\rho}_{c}^{I}$) is even (odd) in $x$, the contribution of the imaginary part of $\tilde{\rho}_{c}$ to the first term of the integral vanishes, and so does the contribution of its real part to the second term. From Eqs. (134) and (135) we can write:

[TABLE]

where the two components read respectively:

[TABLE]

and:

[TABLE]

where we have performed an integration by parts to write the second term of the integrand of $F_{y}^{C}$.

4.4 Vertical component of the force

Prior to the evaluation of the force, we need to work out the Fourier transform of the potential to first order in $I$, using our hypothesis of Eq. (14). It can be obtained from Eq. (118) through a shift of magnitude $I\sin\Omega_{p}t^{\prime}$ in the $z$ direction. It reads:

[TABLE]

where

[TABLE]

As done previously with the other components, we split the $z$-component of the force into two contributions:

[TABLE]

with

[TABLE]

and

[TABLE]

The expression of $F_{z}^{(1)}$ is formally similar to the expression of $F_{y}^{(1)}$ and can be transformed in a similar manner, using Eq. (76):

[TABLE]

where we have $\tilde{\rho}_{i}^{R}\equiv\Re(\tilde{\rho}_{i}^{+})$ and $\tilde{\rho}_{i}^{I}\equiv\Im(\tilde{\rho}_{i}^{+})$. Noting that $\tilde{\Phi}_{i}^{-}=-\tilde{\Phi}_{i}^{+}=k_{z}\tilde{\Phi}_{c}/2$, it is straightforward to recast Eq. (133) as:

[TABLE]

Given the parities in $x$ of the functions $\tilde{\Phi}_{c}$, $\tilde{\rho}_{c}^{R}$ and $\tilde{\rho}_{c}^{I}$, the imaginary part of $\tilde{\rho}_{c}$ yields a null contribution in the integral above, so that $\tilde{\rho}_{c}$ in the integrand can be replaced by $\tilde{\rho}_{c}^{R}$. We can now write, using Eqs. (144) and (145):

[TABLE]

with:

[TABLE]

and

[TABLE]

4.5 Reduction of the force components to a non-dimensional

form

In this section we transform the six integrals that define the cosine and sine amplitude of the three force components. Following Masset (2017), we define:

[TABLE]

and the dimensionless form of the wavevectors $k_{y}$ and $k_{z}$:

[TABLE]

We furthermore define:

[TABLE]

Note that $K_{0}$ differs from $K$ defined at Eq. (78): we have the relationship

[TABLE]

In addition to the non-dimensional functions respectively solutions of the Eqs. (93) and (94), we shall also need the solution of the circular problem, which is:

[TABLE]

where $R_{K}$ and $I_{K}$ are solutions of the system of differential equations (75) and (76) of Masset (2017), and where:

[TABLE]

We mention that $R_{K}+jI_{K}$ is the particular solution $T_{X_{c},K}$ with $X_{c}=0$.

Using Eqs. (81), (84), (129), (149), (150) to (154), we can write:

[TABLE]

where

[TABLE]

is a dimensionless constant, and where for conciseness we note $T^{R}\equiv\Re(T_{X_{c},K})$. Eq. (155) can be recast as:

[TABLE]

with

[TABLE]

Similarly, we have:

[TABLE]

where

[TABLE]

are dimensionless constants. We determine these constants as follows. For definiteness, we specify hereafter to the case of $f_{x}^{C}$, the other five constants being obtained in a strictly similar fashion. We define:

[TABLE]

The integrand of this expression is the same as that of Eq. (156), except that we have explicitly written the dependence of $T$ on $X_{c}$ and $K$, and the integration is performed over $X$ only. The different functions that feature in the integrand of this expression are obtained following the method of section 3.6. The integral is a function of $K_{y}$ and $K_{z}$ only. Indeed $K_{0}$ (Eq. 151), $X_{c}=2K_{0}/3K_{y}$ (Eqs. 77 and 79) and $K$ (Eq. 152) are functions of $K_{y}$ and $K_{z}$ only. Once the functions of the integrand are evaluated, we perform the integral of Eq. (169) using Simpson’s method. The dimensionless coefficient $f_{x}^{C}$ is then obtained by evaluating:

[TABLE]

We can recast this expression as:

[TABLE]

This, in practice, corresponds to the way we numerically evaluate this integral. We use a 2D Simpson method over the plane $(K_{y},K_{z})$, with $(K_{y},K_{z})\in[10^{-7},10^{6}]\times[10^{-7},10^{6}]$, and with a constant spacing in $\log K_{y}$ and $\log K_{z}$. We check that our domain of integration is sufficiently large that the integrand $K_{y}K_{z}\phi_{x}^{C}(K_{y},K_{z})$ is negligible near its edges (see Fig. 1). We apply the same reasoning to all the components and we obtain:

[TABLE]

We show in Fig. 1) that we consider at least all the area embedding the contributions higher than $10^{-12}$ times the peak value. Moreover, we have checked that these values are converged to at least four significant digits with respect to the resolution in $K_{y}$ and $K_{z}$, with respect to the resolution with which we solve the differential equations in $X$ for $T$ and $W$, and with respect to the value of $X_{0}$ that we use to find these functions (see section 3.6).

5 Time evolution of eccentricity and inclination

Having worked out the time-dependent force acting on the planet, we can now evaluate the time derivatives of the eccentricity and inclination. Their expressions are respectively (Tanaka & Ward, 2004, and refs. therein):

[TABLE]

and:

[TABLE]

which, upon averaging over one orbital period, read:

[TABLE]

and

[TABLE]

Defining the thermal time $t_{\mathrm{thermal}}$ as:

[TABLE]

where $F_{00}$ is the value of $F_{0}$ when the planet is non-luminous ($L=0$), and where

[TABLE]

we have:

[TABLE]

with

[TABLE]

It is instructive to compare $t_{\mathrm{thermal}}$ to the damping time $t_{\mathrm{wave}}$ defined by Tanaka & Ward (2004), which reads, with our notation:

[TABLE]

We find:

[TABLE]

where the approximation corresponds to the particular case $\gamma=1.4$. The dimensionless parameter $\ell$ is zero when the planet is non-luminous. In that case, we see that the eccentricity and inclination of the planet are damped, on a timescale that is much shorter than that of Tanaka & Ward (2004). Namely, we have in that case:

[TABLE]

and

[TABLE]

Typically, $H/\lambda_{c}\sim 10$ at a few astronomical units in a protoplanetary disc around a solar-type star (Masset, 2017). This implies that the eccentricity and inclination damping on a non-luminous embryo is completely dominated by thermal effects, which are nearly one order of magnitude more important than those arising from wave launching. This strong damping is in agreement with the findings of Eklund & Masset (2017) who found that the eccentricity and inclination of a non-luminous low-mass planet embedded in a radiative disc was indeed much stronger than that expected from Tanaka & Ward’s formulae. When $\ell>1$, i.e. when the planet’s luminosity is larger than the critical luminosity $L_{c}$, Eqs. (184) and (185) show that the eccentricity and inclination grow exponentially. When $\ell$ is not too close to one, thermal effects largely dominate over those arising from wave launching, and the eccentricity grows at a rate $1.46/0.58\approx 2.5$ times larger than the inclination. This is in correct agreement with the findings of Eklund & Masset (2017). We note that $\ell$ can be of order of a few for an Earth-mass embryo with a mass doubling time of the order of $10^{5}$ yr in a typical protoplanetary disc at a few astronomical units. Naturally, the parameter space is so large that significantly more work is required to assess the relative importance of thermal effects over wave effects. Such study is largely beyond the scope of this work, but we anticipate that, in general, thermal effects should be dominant for planetary embryos in planet-forming regions of protoplanetary discs.

6 Discussion

In the following we discuss a few aspects of our findings, namely their relation with the non-linear corotation torque and horseshoe dynamics, and how they compare to the results previously obtained for planets in circular orbits and perturbers in non-sheared, homogeneous media. We also compare the different timescales for migration, eccentricity and inclination evolution under either thermal or wave-supported disturbances. We discuss the different regimes of eccentricity (opposing shear-dominated to headwind-dominated cases) and we finally briefly discuss the mass range over which the present results are expected to apply.

6.1 Relation to horseshoe dynamics

The processes that we have presented here, which are captured by a linear analysis, occur in the vicinity of the orbit, a location that also gives rise to the corotation torque. It is well known that, in an inviscid disc, the corotation torque always becomes non-linear upon a timescale that depends on the planetary mass (Paardekooper & Papaloizou, 2009). However, the horseshoe dynamics should have little impact, if any, on the processes presented here. The response time of the thermal force, which is the time required for the heat released by the planet to diffuse over a lengthscale $\lambda$, is of the order of the dynamical timescale $\Omega_{p}^{-1}$. This timescale is much shorter than the timescale for the horseshoe U-turns of low-mass planets, which can amount to tens of orbits (see e.g. Masset, 2017). Note also that the horseshoe motion corresponds to minute perturbations of the unperturbed sheared flow. We have seen in section 3.1 that the perturbations of this flow can be neglected: the effects that we present here are described by an advection-diffusion equation, where the advection part comes from the unperturbed flow and the diffusion occurs in the vicinity of a point-like source that describes an epicycle. A full description of the dynamics of the coorbital regions that includes simultaneously thermal effects and the non-linear horseshoe dynamics is probably too complex to be tackled analytically and should be undertaken through numerical experiments. It is however reasonable to expect that such experiments would confirm that thermal effects are essentially decoupled from the horseshoe dynamics, for the following reasons:

• •

The numerical experiments of Eklund & Masset (2017), which were performed prior to the existence of this work and were therefore not aimed at checking analytical predictions with a high accuracy, do show that planets are subjected to effects similar to those described here (i.e. excitation of eccentricity and inclination at large luminosity, as opposed to a strong damping of these quantities when they are non-luminous). The planetary masses considered in these experiments range from a fraction of an Earth mass to about ten Earth masses. The horseshoe region is resolved, and the horseshoe dynamics captured, for planets above one Earth-mass in those experiments444In these experiments the horseshoe region of a one Earth-mass planet spans radially 11 cells..

• •

Arguments similar to those exposed above can be applied to the damping of eccentricity that occurs as a consequence of wave-launching. In that case it is primarily due to wave excitation at coorbital Lindblad resonances (Artymowicz, 1993). Those share their location with the horseshoe region. Yet, there is again a large difference between the response time of the waves (which is typically the dynamical timescale $\Omega_{p}^{-1}$) and the horseshoe U-turn time so that one may expect that the eccentricity and inclination damping may be quite insensitive to the horseshoe dynamics. There is a large body of numerical studies of the eccentricity damping in isothermal discs that support this view, as they confirm the analytical predictions of Tanaka & Ward (2004), even for planets subjected to a horseshoe drag (e.g. Cresswell et al., 2007).

We finally mention that planets that have a luminosity $L$ larger than the critical luminosity $L_{c}$ experience a growth of eccentricity and inclination. In these circumstances, the coorbital corotation torque is quenched (Hellary & Nelson, 2012).

6.2 Comparison to earlier results

The analysis presented here is the last one of a series of three analytical studies devoted to the role of thermal diffusion and luminosity feedback in situations of interest for planet-disc interaction.

• •

The first one dealt with the simplest setup, that of a perturber moving across a 3D medium at rest (a setup that gives rise to dynamical friction when the luminosity feedback is not taken into account), with special emphasis on the low Mach number regime. In that case, the two additional forces (the force arising from the perturber’s luminosity and that arising from the mere inclusion of thermal diffusion, when the perturber is non-luminous) have been studied in two separate publications. The former, dubbed heating force, has been studied by Masset & Velasco Romero (2017, [MV17]), while the latter, named cold thermal force, has been evaluated analytically by Velasco Romero & Masset (2019, [VM19]) and corroborated numerically in that same work. The net force arising from thermal effects is a drag when the perturber’s luminosity $L$ is smaller than a critical luminosity $L_{c}$ that has the exact same expression as Eq. (58), and a thrust otherwise.

• •

The second analysis dealt with planets on circular orbits, and analysed the new torque components that arise from the heat release by the planet (or heating torque) and from the inclusion of thermal diffusion (or cold thermal torque). While primarily intended as an analytical follow-up of the work of Benítez-Llambay et al. (2015, [BLl+15]), this analysis allowed to identify the effect found by Lega et al. (2014, [L+14]), and dubbed the cold finger effect by these authors, with the cold thermal torque. Again, this analysis showed that the net thermal torque changes sign for the perturber’s luminosity $L=L_{c}$, which has again the exact same value as that given by Eq. (58).

These two cases and the case presented in this work are summarized in Fig. 2. The critical luminosity $L_{c}$ of Eq. (58) appears as a universal watershed at which the net thermal effects change sign. This has been discussed by Masset (2017, section 5.7 and Fig. 3). In an adiabatic case, an enthalpy (or temperature) peak surrounds the planet (for low-mass planets the peak is nearly the negative of the potential well). The introduction of thermal diffusion tends to reduce large temperature gradients and flattens this peak. The release of energy into the surrounding gas by a luminous perturber rebuilds the peak. When the luminosity of a perturber with low Mach number is equal to the critical luminosity $L_{c}$, the temperature peak that surrounds the perturber has same amplitude and shape as that of the adiabatic case (even though the underlying physical processes responsible for the peak are quite different) and the thermal effects cancel out.

For this specific case with $L=L_{c}$, one recovers in an unsheared homogeneous gas the force studied by Ostriker (1999, [O99]), whereas for the cases of a planet on circular orbit and the more general case with finite but small eccentricity and inclination, one recovers555The studies of TTW02 and TW04 are strictly speaking for isothermal discs, rather than adiabatic discs. None the less, one can infer from these works the behaviour in an adiabatic disc by substituting the isothermal sound speed with the adiabatic sound speed. respectively the torque of Tanaka et al. (2002, [TTW02]) and the damping rates of Tanaka & Ward (2004, [TW04]). Due to the large value of thermal forces on low-mass planets in planet-forming regions of protoplanetary discs, only when $L\approx L_{c}$ do we nearly recover the results for adiabatic discs. The vertical green band at $L_{c}$ (in the electronic version) has been intentionally represented quite narrow to illustrate this effect.

6.3 Comparison of timescales

The effects presented here are considerable, and dominate over those arising from wave-launching for low-mass planets except when $L\approx L_{c}$. Eqs. (189) and (190) show that the eccentricity and inclination damping timescales of a non-luminous planet are a factor of $\lambda/H$ shorter than those given by Tanaka & Ward (2004). We note that the same factor appeared between the migration timescale arising from the cold thermal torque and that arising from wave-launching (Masset, 2017, Eq. [137]). The timescale for eccentricity and inclination damping being, in adiabatic discs, a factor of $(H/r)^{2}$ shorter than the migration timescale Artymowicz (1993); Tanaka & Ward (2004), we therefore expect the same ratio to hold between the timescales for the evolution of the eccentricity and inclination and that for the evolution of the semi-major axis, under thermal forces. From Eqs. (133), (138) and (145) of Masset (2017), we see that the migration timescale of a non-luminous planet subjected to the thermal torque is, in order of magnitude:

[TABLE]

which is indeed a factor $\sim(H/r)^{-2}$ larger than the thermal time of Eq. (182). We summarize these different relationships in Fig. 3.

Note that the estimates of the timescales from thermal forces are for a non-luminous planet. These characteristic timescales for a planet with a luminosity largely in excess of $L_{c}$ would be even shorter666They would then be growth timescales for the eccentricity and inclination, and a timescale of outward migration for the semi-major axis..

6.4 Different regimes of eccentricity

Planets with super-critical luminosity ($L>L_{c}$) experience an exponential growth of eccentricity and inclination with time. At some point the hypothesis of Eq. (14) that the epicyclic and vertical excursions $E$ and $I$ are much smaller than the size of the thermal disturbance ceases to be valid. It is straightforward to realise that in these circumstances the thermal disturbance tends towards that triggered by a perturber in an unsheared medium. The velocity $v_{p}$ of the perturber with respect to the ambient gas is indeed larger than $\lambda_{c}\Omega_{p}$. The response time of the thermal disturbance is then (Masset & Velasco Romero, 2017):

[TABLE]

The response time being shorter than the shear timescale, the shear becomes unimportant and the thermal force tends towards that of unsheared media (Masset & Velasco Romero, 2017; Velasco Romero & Masset, 2019). This regime has been named the headwind-dominated regime by Eklund & Masset (2017, section 4.5), as opposed to the regime of low eccentricities and inclinations that we studied in the present work, which was referred to as the shear-dominated regime. In the same vein, this effect has been named the hot trail effect by Chrenko et al. (2017, [C+17]). The same kind of transition from shear-dominated to headwind-dominated occurs for the pressure-supported wake, except that it occurs at epicyclic or vertical excursions comparable to the pressure length scale $H$. The calculation of the eccentricity and inclination damping rates of Tanaka & Ward (2004) requires that $e\ll H/r$ and $i\ll H/r$, and this, in general, is the case of estimates of damping rates based on a sum of resonances (Ward, 1988; Artymowicz, 1993), as for eccentricities larger than the disc’s aspect ratio the use of series on resonances becomes impractical (Papaloizou & Larwood, 2000). When the eccentricities are larger than the disc’s aspect ratio, a dynamical friction calculation is much more convenient (Papaloizou, 2002; Muto et al., 2011, resp. [P02] and [M11]). The case of inclinations larger than the disc’s aspect ratio is slightly different, as the planet spends a fraction of its orbit outside of the disc, but is also conveniently dealt with using a dynamical friction calculation (Rein, 2012). We depict the different regimes of eccentricity in Fig. 4.

The asymptotic values reached by the eccentricity in the numerical experiments of Eklund & Masset (2017) typically fall within the $[\lambda,H]$ interval. In these conditions the equilibrium eccentricity is given by the balance between the time-varying force of Tanaka & Ward (2004) for the pressure-supported disturbance, and by the heating force in an unsheared medium (Masset & Velasco Romero, 2017; Velasco Romero & Masset, 2019).

6.5 Dependence of thermal forces on the perturber’s mass

As discussed by Masset & Velasco Romero (2017) in the context of unsheared media, the thermal forces are expected to suffer a decay with respect to their linearly predicted value when $M_{p}>M_{c}$, where the critical mass $M_{c}$ is given by:

[TABLE]

When the perturber’s mass is much smaller than $M_{c}$, the heat diffusion time across the planetary Bondi radius is much smaller than the acoustic time across the Bondi radius, which guarantees that the energy released by the planet in its immediate neighbourhood reaches outside the Bondi sphere, where our linear analysis is valid, as an excess of internal energy. This no longer needs to be the case when $M_{p}\gtrsim M_{c}$, however. Using arguments based on the yield of the heating force, Masset & Velasco Romero (2017) argue that the latter has to be cut off for $M_{p}>M_{c}$. The numerical value of $M_{c}$, in planet-forming regions of protoplanetary discs, is of the order of an Earth-mass, although this value can vary greatly as a function of the position in the disc and as the discs evolves and cools (Masset & Velasco Romero, 2017; Masset, 2017). The numerical experiments of Lega et al. (2014), Benítez-Llambay et al. (2015) and Eklund & Masset (2017) are all compatible with a cut-off of thermal effects above masses commensurable with an Earth-mass. An accurate determination of the cut-off law probably requires high-resolution calculations that resolve the Bondi sphere, as this effect cannot be captured by means of a linear calculation. We mention none the less that thermal effects are so vigorous that even in numerical experiments of Eklund & Masset (2017), a sizeable impact of the planet’s radiative feedback on its eccentricity and inclination is found up to approximately 5-10 Earth masses.

7 Conclusions

We have worked out the time-dependent force acting on a planetary embryo embedded in a gaseous protoplanetary disc, using linear perturbation theory, when thermal diffusion in the gas is taken into account, with or without heat release by the planet into the surrounding gas.

We find that this thermal force has a strong impact on the eccentricity and inclination of the embryo, with an outcome that depends on the embryo’s luminosity $L$. When the latter is smaller than the critical luminosity $L_{c}$ defined at Eq. (58), the thermal force leads to a damping of eccentricity and inclination, typically stronger by up to an order of magnitude than the damping due to wave-launching considered so far (Ward, 1988; Artymowicz, 1993, 1994; Tanaka & Ward, 2004). The sign of thermal forces reverses when the embryo’s luminosity is $L_{c}$. Only in the regime where $L\approx L_{c}$ does the damping due to wave-launching play a role, as it is otherwise masked by the strong effect of the thermal force. For luminosities significantly greater than $L_{c}$, as can be expected for Earth-like embryos with mass doubling times shorter than $100$ kyr at a few astronomical units in discs similar to the Minimum Mass Solar Nebula, the eccentricity and inclination grow exponentially over a short time scale of the order of hundreds of orbits only. The outcome of such growth has been studied numerically by Eklund & Masset (2017) and Chrenko et al. (2017) (in 2D discs). The critical luminosity $L_{c}$ to get a reversal from damping to excitation is the same for the eccentricity and inclination. It is also the same as the critical luminosity at which the thermal torque on a planet on a circular orbit (Masset, 2017) reverses from negative (dominated by the cold thermal torque, see Lega et al., 2014) to positive (dominated by the heating torque, see Benítez-Llambay et al., 2015), and the same as well at which the thermal force on a luminous perturber moving across a uniform medium reverses from drag to thrust (Velasco Romero & Masset, 2019). While the thermal force on planets in discs has not been studied numerically in a systematic fashion, the recent numerical simulations of Velasco Romero & Masset (2019) corroborate with a high accuracy the value of the critical luminosity in unsheared, homogeneous media.

The numerical simulations of Lega et al. (2014), Benítez-Llambay et al. (2015) and Eklund & Masset (2017) all seem to indicate that the effects of thermal forces are cut off above a few Earth masses, for the set of parameters considered in these studies. However, this effect has not been studied in detail and should probably be tackled through numerical simulations.

The effects that we present here should have important consequences on various stages of planet formation, such as the phase of oligarchic growth or the giant impact stage, when those occur in the gaseous disc. It could have also consequences on the trapping in mean motion resonances of Earth-sized protoplanets. A detailed study of such effects requires to track the thermal and luminous history of embryos, taking into account the accretion of solid bodies and possibly the effect of mutual collisions.

8 Acknowledgments

For this work, SF was supported by the Programa de Apoyo a Proyectos de Investigacion e Innovacion Tecnologica (PAPIIT) No IA101619 and FM was supported by the Programa de Apoyo a Proyectos de Investigacion e Innovacion Tecnologica (PAPIIT) No IN101616.

Appendix A Impact of corotation offset

Our derivation assumes that the distance $x_{p}^{0}$ between corotation and the epicycle’s guiding centre is small compared to the radial and vertical excursion of the planet. We hereafter discuss what happens when this assumption is relaxed. It is convenient, in this case, to have the origin of the $x$ axis at the planet’s guiding centre, so that the potential terms of Eqs. (118), (121) and (122) keep the same form, whereas the unperturbed azimuthal velocity has now the form:

[TABLE]

With a finite value of $x_{p}^{0}$, the relationships of Eqs. (74) and (76) are no longer verified, as they relied on the symmetry in $x$ of the response. The expressions of the force components worked out in sections 4.2 and 4.3 can be expressed in a slightly more general fashion as follows:

[TABLE]

and

[TABLE]

where we have used the relations

[TABLE]

and

[TABLE]

Using Eqs. (118), (121) and (122) we write

[TABLE]

and

[TABLE]

to get rid of all instances of $\tilde{\rho}_{c}$ and $\tilde{\Phi}_{c}$ in the time derivative of the eccentricity given by Eq. (178), and eventually obtain:

[TABLE]

This, with our notation, is equivalent to the original derivation of Tanaka & Ward (2004) who checked that their force expression was compatible with the time derivative of the eccentricity given by a sum on the first-order Lindblad resonances (Goldreich & Tremaine, 1980; Artymowicz, 1993). Although the expressions of the “circular” terms $\tilde{\rho}_{c}$ and $\tilde{\Phi}_{c}$ are required to obtain the individual force components of Eqs. (172) to (175), they cancel out in the expression of the time derivative of the eccentricity. Assessing how the latter varies with $x_{p}^{0}$ is therefore tantamount to assessing how $\tilde{\rho}_{e}^{\pm}$ varies with $x_{p}^{0}$. Eq. (80) becomes:

[TABLE]

where

[TABLE]

Calling $\tilde{\rho}_{e,0}^{+}$ the solution for $X_{p}^{0}=0$ and writing $\delta\tilde{\rho}_{e}^{+}=\tilde{\rho}_{e}^{+}-\tilde{\rho}_{e,0}^{+}$, we have:

[TABLE]

From Fig. 1, it is evident that the spatial frequencies that most contribute to the response are $|K_{y}|\sim|K_{z}|\sim 1$ (hence $k_{y,z}\sim\lambda_{c}^{-1}$), and for those $|X_{c}|\sim 1$ and $X_{p}^{0}\sim x_{p}^{0}/\lambda_{c}$. From Eq. (204), we deduce that as long as $|X_{p}^{0}|\ll|X_{c}|$ (i.e. $x_{p}^{0}\ll\lambda_{c}$), the correction $\delta\tilde{\rho}_{e}^{+}$ is small compared to the symmetric solution $\tilde{\rho}_{e,0}^{+}$. The time derivative of the eccentricity is therefore correct as long as $|x_{p}^{0}|\ll\lambda_{c}$, regardless of whether it is smaller or larger than the epicyclic excursion $E$, and the growth of eccentricity when $L>L_{c}$ is not a finite amplitude instability. The numerical experiments of Eklund & Masset (2017, e.g. their Fig. 14) agree with this statement. Note that similar arguments apply to the damping of eccentricity arising from wave launching. The expression of Tanaka & Ward (2004) should be valid all the way to eccentricities much smaller than the offset between corotation and guiding centre. On physical grounds, the eccentricity varies because the perturber is subjected to a force that depends on its position on the epicycle. As long as the corotation offset is not too large, shifting corotation amounts to adding a constant force which has no impact on the eccentricity budget.

Strictly similar considerations apply to the excitation or damping of the inclination.

Appendix B Extension to a viscous disc

Our derivation has not considered how viscous heating would affect the energy budget in Eq. (5). If we consider a laminar disc with kinematic viscosity $\nu$, the dominant source term for the viscous heating is:

[TABLE]

and the dominant term arising from the perturbation induced by the planet is, from Eq. (11):

[TABLE]

where, for lengthscales typical of that of the perturbation, $|\partial_{x}v^{\prime}|\sim|v^{\prime}|/\lambda$. Eq. (30)–(32) imply that, in order of magnitude for $|x|\sim\lambda$, $v^{\prime}\sim GM/(R^{2}\Omega_{p})$, where $R$ is the distance to the planet. We can then write the order-of-magnitude relationships $|v^{\prime}|\sim(R_{B}/R)H^{2}\Omega_{p}/R\sim(\rho^{\prime}/\rho_{0})H^{2}\Omega_{p}/\lambda$, where $R_{B}=GM/c_{s}^{2}$ is the planetary Bondi radius. Eq. (206) is then dominated by its second term, which has order of magnitude $\rho^{\prime}\nu\Omega_{p}^{2}(H/\lambda)^{2}$. This source term has to be compared to the divergence of the heat flux, which has the order of magnitude $\chi c_{s}^{2}\rho^{\prime}/\lambda^{2}\sim\rho^{\prime}\chi\Omega_{p}^{2}(H/\lambda)^{2}$. In discs that have a Prandtl number $\mathrm{Pr}\equiv\chi/\nu\gg 1$, the viscous dissipation induced by the planetary perturbation is negligible compared to the heat flux. Our analysis should remain valid in such discs. The numerical exploration of Eklund & Masset (2017) took place in a disc with $\mathrm{Pr}\sim 5$ and yielded results compatible with this expectation.