Numerical study of dynamical friction with thermal effects -- I. Comparison to linear theory

TL;DR

This study uses numerical simulations to analyze dynamical friction with thermal effects, comparing results to linear theory, and finds that thermal diffusion significantly alters the drag force on perturbers, with implications for planetary formation.

Contribution

It provides a detailed numerical analysis of thermal effects on dynamical friction, including the heating force and non-adiabatic drag, validated against linear theory and implemented with GPU-accelerated nested meshes.

Findings

Heating force aligns with analytic estimates at low and high Mach numbers.

Non-luminous drag behaves like dry or solid friction in non-adiabatic gases.

Luminosity threshold for net acceleration is lower than typical planetary embryo luminosities.

Abstract

We investigate by means of numerical simulations with nested Cartesian meshes the force exerted on a massive and luminous perturber moving at constant speed through a homogeneous and opaque gas, taking into account thermal diffusion in the gas and the radiative feedback from the perturber. The force arising from the release of energy into the ambient medium by the perturber, or heating force, is directed along the direction of motion and induces an acceleration of the perturber. Its value is compared to analytic estimates in the low and high Mach number limits, and found to match those accurately. In addition, the drag exerted on a non-luminous perturber significantly departs from the adiabatic expression when thermal diffusion is taken into account. In the limit of a vanishing velocity, this drag tends to a finite value which we determine using linear perturbation theory and corroborate with numerical simulations. The drag on a non-luminous perturber in a non-adiabatic gas therefore behaves like dry or solid friction. We work out the luminosity threshold to get a net acceleration of the perturber and find it to be generally much smaller than the luminosity of accreting low-mass planetary embryos embedded in a gaseous protoplanetary disc at a few astronomical units. We also present in some detail our implementation of nested meshes, which runs in parallel over several \emph{Graphics Processing Units} (GPUs).

Figures

Equations

$${\partial }_t\rho +\nabla \cdot (\rho \mathbf{v})=0$$

$${\partial }_t\mathbf{v}+\mathbf{v}(\nabla \cdot \mathbf{v})=-\frac{\nabla p}{\rho }-\nabla \Phi$$

$${\partial }_{t}e+\nabla \cdot (e\mathbf{v})+p\nabla \cdot \mathbf{v}+\nabla \cdot \left(\chi \rho \nabla \frac{e}{\rho }\right)=L\delta (\mathbf{r})$$

$$a_{M,L,\chi }=\iiint \frac{G{\rho }_{M,L,\chi }(x,y,z)z}{(x^2+y^2+z^2{)}^{3/2}}dxdydz$$

$$F_{\mathrm{heating}}=M{a}_{M,L,\chi }-M{a}_{M,0,\chi }.$$

$${\rho }_{0,\mathrm{L},\chi }(x,y,z)=-\frac{\gamma (\gamma -1)L}{4\pi \chi c_s^{2}r}\exp \left(\frac{z-r}{2\lambda }\right),$$

$$\lambda =\chi /\gamma V$$

$$a_{0,L,\chi }=\frac{\gamma (\gamma -1)GL}{2\chi c_s^2}\text{sign}(V),$$

$$\tau =\frac{{\lambda }^2}{\chi }=\frac{\chi }{{\gamma }^2{V}^2}$$

$$a_{0,L,\chi }=\frac{2(\gamma -1)GL}{\sqrt{\pi }\chi V^2}f\left(\frac{r_{\min }V}{4\chi }\right),$$

$$f(ϵ)={\int }_{ϵ}^{\infty }\mu \left(\sqrt[]{u}\right)du,$$

$$\mu (u)=u{e}^{-u^2}+\frac{\sqrt[]{\pi }}{2}(u^{-2}-2)\text{erfc}(u).$$

$$a_{0,L,\chi }\approx \frac{(\gamma -1)GL}{\chi V^2}\left[-1.96-\log \left(\frac{r_{\min }V}{4\chi }\right)\right].$$

$$S(\Delta t)⟶T(\Delta t),$$

$$S(\Delta t/2)⟶T(\Delta t)⟶S(\Delta t/2),$$

$$\Phi ({\mathbf{r}}_{i,j,k})=\frac{GM}{\sqrt[]{|{\mathbf{r}}_{i,j,k}-{\mathbf{r}}_0{|}^2+{ϵ}^2},}$$

$$R_B=GM/c_s^2,$$

$$F=\sum _l\frac{z_{l}GM{\rho }_l{V}_l}{(r_l^{\prime 2}+{ϵ}^2{)}^{3/2}},$$

$$f_{i-\frac{1}{2}}^x=\chi {\rho }_0\left(\frac{e}{\rho }-\frac{e_{i-1}}{{\rho }_{i-1}}\right)\frac{S_x}{\Delta x}$$

$$f_{j-\frac{1}{2}}^y=\chi {\rho }_0\left(\frac{e}{\rho }-\frac{e_{j-1}}{{\rho }_{j-1}}\right)\frac{S_y}{\Delta y}$$

$$f_{k-\frac{1}{2}}^z=\chi {\rho }_0\left(\frac{e}{\rho }-\frac{e_{k-1}}{{\rho }_{k-1}}\right)\frac{S_z}{\Delta z}.$$

$$e_{i,j,k}-=\Delta t\left(f_{i+\frac{1}{2}}^x-f_{i-\frac{1}{2}}^x+f_{j+\frac{1}{2}}^y-f_{j-\frac{1}{2}}^y+f_{k+\frac{1}{2}}^z-f_{k-\frac{1}{2}}^z\right)V^{-1},$$

$$h(r^{\prime })=\frac{1-2r^{\prime }-e^{-2r^{\prime }}}{2r^{\prime }\text{texttwosuperior}}+1,$$

$$\Delta F_{0,L,\chi }=\frac{(\gamma -1)GML}{\chi V^2}\log 2.$$

$$\xi \sim \frac{4}{\sqrt{\pi }\gamma }f\left(\frac{r_{\text{min}}V}{4\chi }\right).$$

$${\partial }_t{\rho }^{\prime }+{\rho }_0\nabla \cdot {\mathbf{v}}^{\prime }=0$$

$${\partial }_t{\mathbf{v}}^{\prime }+\frac{\nabla p^{\prime }}{{\rho }_0}=\nabla {\Phi }_p$$

$${\partial }_t{p}^{\prime }+\gamma p_0\nabla \cdot {\mathbf{v}}^{\prime }+\frac{p_0}{{\rho }_0}\chi \Delta {\rho }^{\prime }-\chi \Delta p^{\prime }=0,$$

$$-i\omega \tilde{\rho }+i{\rho }_0{k}_j{\tilde{v}}_j=0,j=1,2,3$$

$$-i\omega {\tilde{v}}_j+i{k}_j\frac{\tilde{p}}{{\rho }_0}=-i{k}_j{\tilde{\Phi }}_p$$

$$-i\omega \tilde{p}+\gamma p_{0}i{k}_j{\tilde{v}}_j-\frac{p_0}{{\rho }_0}\chi k^2\tilde{\rho }+\chi k^2\tilde{p}=0.$$

$${\tilde{\Phi }}_p=-\frac{4\pi G{\tilde{\rho }}_{\mathrm{p}}}{k^2},$$

$$\tilde{p}=\frac{4\pi G{\rho }_0{\tilde{\rho }}_p+{\omega }^2\tilde{\rho }}{k^2},$$

$$\tilde{\rho }=\frac{4\pi G{\rho }_0\left(\frac{i\omega }{k^2}-\chi \right)}{i\omega \left(c_s^2-\frac{{\omega }^2}{k^2}\right)+\chi \left({\omega }^2-k^2{c}_s^2/\gamma \right)}{\tilde{\rho }}_p,$$

Full text

††pubyear: 2018††pagerange: Numerical study of dynamical friction with thermal effects – I. Comparison to linear theory–A

1 Introduction

Recent developments in the field of planet-disc interactions show the importance of heat release into the surrounding gas on the dynamical evolution of hot planetary embryos. The radiative feed back from an embryo yields the appearance of a low density region in its vicinity which can help it withstand or even revert inward migration (Benítez-Llambay et al., 2015; Masset, 2017). This low density region can also excite the embryo’s eccentricity (Eklund & Masset, 2017; Chrenko et al., 2017; Masset et al., 2017) and inclination (Eklund & Masset, 2017; Masset et al., 2017) to sizeable values.

The case of dynamical friction is somewhat simpler, as it does not involve a sheared medium such as a nearly Keplerian disc, but a uniform medium at rest. Yet in this case the release of radiative energy by the perturber into its surroundings also has a strong impact on the back reaction exerted on the perturber, to the point that it may yield a net acceleration, the total force being then a thrust rather than a drag (Masset & Velasco Romero, 2017; Park & Bogdanović, 2017). Masset & Velasco Romero (2017, hereafter MV17) give analytical expressions in the linear regime for the extra force resulting from the radiative feedback (dubbed heating force) in the limit of low and high Mach numbers. They find that in the low Mach number limit, the heating force has an asymptotic, finite value independent of the perturber’s velocity, while in the supersonic regime, the heating force decays as ${\cal M}^{-2}$, ${\cal M}$ being the Mach number. In both cases they find the heating force to be positive.

The aim of the present work is to corroborate these asymptotic behaviours by means of numerical simulations and to investigate the regime of intermediate Mach numbers, for which MV17 did not provide an analytical estimate. In addition, the expression of the total force requires a knowledge of the net force acting on a non-luminous perturber when thermal diffusion is taken into account. MV17 made the simple assumption that this force was comparable to that of the adiabatic case, worked out by Ostriker (1999) and numerically confirmed by Sánchez-Salcedo & Brandenburg (1999). We relax here this assumption and study in some detail the force acting on a cold perturber in a medium with thermal diffusion.

A clean numerical “test bed” for this problem can be built in the rest frame of the perturber. The gas is initially set to have a velocity opposite to that of the perturber in the gas frame, and it is injected at any instant in time at the upstream boundary with that velocity, while heat diffusion is modelled with a uniform thermal diffusivity over the grid. Owing to the large range of length scales involved in this problem, it is impractical to simulate the system on a single mesh, and we must resort to nested meshes.

MV17 argue that their linear analysis should be valid when the ratio $GM/\chi c_{s}$ is small, where $\chi$ is the thermal diffusivity of the gas and $c_{s}$ the adiabatic sound speed. All the simulations presented here have $GM/\chi c_{s}\leq 10^{-1}$ (for most of them this ratio is $10^{-2}$). This ratio can be regarded as the ratio between the time scale $R_{B}^{2}/\chi$ for diffusion of heat across the perturber’s Bondi radius $R_{B}=GM/c_{s}^{2}$ and the acoustic time scale across the Bondi radius111Throughout this work we consider the Bondi radius $R_{B}=GM/c_{s}^{2}$ instead of the effective one $R_{B_{\text{eff}}}=GM/(c_{s}^{2}+V^{2})$. Whenever the former is unresolved so is the latter and we therefore only capture the flow in regions where the relative perturbation of density arising from the perturber’s gravity is small.. In the present work, the first of two papers studying the impact of thermal effects and heat release on dynamical friction, we do not attempt to capture the dynamics of the flow within the perturber’s Bondi radius, so as to deal only with weakly non-linear flows. Our mesh resolution will therefore be of the order of, or larger than the Bondi radius, and the thermal and acoustic time scales across the Bondi radius are not resolved. Calculations resolving the flow within the Bondi radius and corresponding to arbitrary values of $GM/\chi c_{s}$ will be presented in the second paper.

In section 2 we present a description of the problem at hand, remind the main results of MV17 and introduce the conventions used throughout this paper. We run our simulations with the hydrocode FARGO3D222http://fargo.in2p3.fr that we have modified to handle nested meshes. In section 3 we detail how this new capability of the code was implemented. In section 4 we present our results, first for the heating force, then for the drag force on a non-luminous perturber in presence of thermal diffusion and finally for the net (total) force on a luminous perturber. In section 5 we discuss our findings and give our conclusions.

2 Problem description

We introduce a point-like perturber of mass $M$ and luminosity $L$ in an initially uniform gaseous medium with adiabatic sound speed $c_{s}$. The perturber is kept fixed whilst the gas has initially the uniform velocity $V={\cal M}c_{s}$. Without loss of generality, we assume the perturber to lie at the origin of a Cartesian frame $(xyz)$ and the gas velocity to be directed along the $z$ axis. The heat released by the perturber in the gas diffuses according to Fourier’s law with the uniform and constant thermal diffusivity $\chi$. Resorting to this simple modelling of thermal effects is an oversimple description of the radiative processes at work in the perturber’s vicinity. While such a simplification is necessary to make the problem analytically tractable, a numerical treatment could consider a more realistic approach (e.g. Benítez-Llambay et al., 2015; Park & Bogdanović, 2017). However, sticking to the same simple description of thermal effects as that of the analytical work of MV17 allows a clean and accurate comparison between analytical expectations and the results of simulations, let aside the moderate computational cost of such prescription.

2.1 Governing equations

We denote respectively with $\rho$, $\boldsymbol{v}$ and $e$ the gas density, velocity and density of internal energy, and add a subscript “[math]” to refer to unperturbed quantities. The governing equations are the continuity equation, Euler’s equation for the velocity and the energy equation, which read respectively:

[TABLE]

where $\Phi=-GM/|\boldsymbol{r}|$ is the perturber’s potential and $L\delta(\boldsymbol{r})$ the energy source arising from the radiative feed back on the gas. In these equations $p=(\gamma-1)e$ represents the pressure and $\gamma$ the adiabatic index. The fourth term on the left hand side of Eq. (3) arises from thermal diffusion.

2.2 Conventions used in this work

To avoid ambiguity when necessary, we use the indices $M,L,\chi$ for the solution of Eqs. (1) to (3). We also use these indices to refer to the net acceleration $a$ imparted to the perturber:

[TABLE]

The acceleration $a_{M,0,0}$ is for instance the acceleration imparted to a non-luminous333Throughout this paper we will indifferently call a perturber with $L=0$ either a non-luminous or a cold perturber. perturber of mass $M$ in an adiabatic gas ($\chi=0$). When multiplied by $M$, this corresponds to the drag force worked out by Ostriker (1999). Similarly, $a_{0,L,\chi}$ is the acceleration imparted to a massless heat source of luminosity $L$ by the perturbation of density triggered in its surroundings. When the relative perturbations of the flow variables are small, they scale linearly with $M$ and $L$, and so does $a_{M,L,\chi}$ by virtue of Eq. (4).

Fig. 1 summarises the different naming conventions that we will be using throughout this work. The total force acting on a perturber $Ma_{M,L,\chi}$ is decomposed into different components: the adiabatic drag $Ma_{M,0,0}$, the cold thermal force that is the difference between the force exerted on the cold perturber of same mass in a disc with thermal diffusion and an adiabatic disc, and the heating force, which is defined by:

[TABLE]

Only when the flow’s perturbation is nearly linear can this force be expressed as $Ma_{0,L,\chi}$.

2.3 Summary of analytical results

We recap hereafter the results obtained by MV17 that we want to confirm using numerical simulations. MV17 obtained the expression of the perturbation of density induced by heat release in the limit of a low Mach number, in the linear regime, and derived the corresponding heating force. They also worked out the expression of the heating force in the case of a large Mach number. All their calculations assume the flow to be in steady state in the perturber’s frame.

The perturbation of density due to the heat release given by a linear analysis in the low Mach number limit and in steady state reads, in the case $V>0$:

[TABLE]

where

[TABLE]

is the cut-off distance arising from the competition between diffusion and advection and represents the characteristic size of the hot, low density “plume” that surrounds the perturber. The acceleration imparted to a test mass in $r=0$ due to this perturbation of density is:

[TABLE]

This acceleration is directed along the motion of the perturber and is independent on the velocity $V$. The response time

[TABLE]

is an estimate of the settling time for the heating force.

In the highly supersonic regime, MV17 give for the acceleration arising from the energy release by the perturber the expression:

[TABLE]

where $r_{\mathrm{min}}$ is the minimum radius for the force integration and where $f(\epsilon)$ is given by:

[TABLE]

the function $\mu(u)$ being itself defined as:

[TABLE]

When the argument of the function $f$ is small (i.e. when the minimal length scale is much smaller than the cut-off distance $\chi/V$), the acceleration can be approximated as:

[TABLE]

This acceleration too is along the direction of motion of the perturber. We have designed numerical setups to test Eqs. (6), (8), (10) and (13).

3 Numerical implementation

We make use of the FARGO3D code (Benítez-Llambay & Masset, 2016) with a refinement capability through the addition of nested meshes similar to that of the NIRVANA code (D’Angelo et al., 2003). Our implementation provides the possibility of creating a hierarchy of nested meshes with an arbitrary number of levels of refinement, each finer level doubling the resolution of its coarser predecessor along each dimension. Within this system of nested meshes, each mesh passes individually through the hydrodynamical solver of FARGO3D, the full set of meshes being stitched together by means of inter-level communications at specific stages of the simulation to ensure communicating concurrent information. This nested meshes version of FARGO3D is also implemented in CUDA and MPI, providing the possibility of running in parallel over multiple GPUs. Results already obtained with the FARGO3D code with nested meshes in a different context (that of low mass planets embedded in protoplanetary gaseous discs) have been presented elsewhere (Fung et al., 2017). However, the description of our nested mesh implementation in that work was limited. We provide here a more complete description of the different stages required by our implementation. The reader not interested in the implementation details may skip this section and go directly to section 3.4.

3.1 Time step scheduling

In many circumstances the Courant condition implies that the maximum timestep allowed on a given level scales with its resolution. In such a case, for each time step $\Delta t$ on a given level $\ell$, the common practice consists in performing two time steps $\Delta t/2$ on the next level $\ell+1$, in a recursive manner (D’Angelo et al., 2002, 2003; Pepliński et al., 2008). This sub-cycling technique is the one we employed in the present work, because it is the most efficient for the problem at hand. However, for different setups such as a planet embedded in a sheared flow, it may not be the best approach. We have developed a general criterion which automatically determines the best sub-cycling scheme for all levels. Since this paper is intended as a reference for our implementation of nested meshes in FARGO3D, we give the detail of our sub-cycling criterion in appendix A.

3.2 Inter-level communications

In order to coordinate individual meshes into a unique system we need to communicate information in both directions: upwards or coarser to finer, and downwards or finer to coarser. In the upwards direction we have one type of communication (CommUp), in which we communicate data from the coarser level to fill the ghost cells of the finer level. For the downwards direction we have two types of communications: CommMean, where a given number of layers of cells in the coarse level are filled by averaging the underlying finer cells and CommFlux, where the fluxes of the finer level at the boundary between levels are communicated to the coarse level and take precedence over the fluxes of the coarse mesh.

In the presentation below we consider communications between level $\ell$ (referred to as the coarse or coarser level) and level $\ell+1$ (referred to as the fine or finer level). A full mesh consists of the active cells and the ghost or buffer cells, that surround the active mesh and allow to enforce boundary conditions or communications between meshes.

3.2.1 Filling the ghost cells of a finer level: CommUp

The ghost cells of a finer level $\ell+1$ are updated before entering the hydrodynamical solver. These ghost cells are filled by a trilinear interpolation of $27$ cells at the coarser level (as shown in fig 2). For a fine cell hosted by a coarse cell $(i,j,k)$, the trilinear interpolation starts in the $z$ direction by interpolating values at $k-1$, $k$ and $k+1$ to the position of the fine cell $k_{fine}$, therefore reducing the initial $27$ values into $9$, which are further reduced in the $y$ direction by interpolating values at $j-1$, $j$ and $j+1$ into position $j_{fine}$ to obtain $3$ values at $i-1$, $i$ and $i+1$, which are eventually interpolated at position $i_{fine}$ to obtain the updated value of the fine cell. This sequence is depicted in Fig. 2. For each of these interpolations, van Leer’s slope limiter is used (van Leer, 1977).

This type of communication is needed at the beginning of the cycle on a finer level, which is when both finer and coarser levels are at the same date. Therefore upwards communications are omitted between the two sub-cycles of the finer level in order to avoid communicating outdated information from the coarser level. This requires to double the width of the ghost buffer (its ultimate width depends on the algorithms used).

3.2.2 Overwriting coarse cells from a fine mesh: CommMean

At the end of the time step on a coarse level, we overwrite several layers of cells covered by the finer level on the contour of the latter, as depicted in Fig. 3. In this implementation it suffices to overwrite four layers given that the number of stages of our solver prevents that information on the coarse level inward of 4 layers reaches the non-refined region. The overwritten cells are filled by averaging the underlying finer cells. Fig 4 shows that for centred quantities we replace the coarse value with the average of the $8$ finer cells that it contains. For staggered quantities we average over the four finer cells aligned at the staggered dimension.

3.2.3 Matching fluxes: CommFlux

In order to enforce the conservation of conserved quantities to platform accuracy across all levels we need to match their fluxes at the frontiers between levels. We use the fluxes from the contour of the fine level to determine what goes into the coarse cells that surround this level. In practice this is achieved by sending the fluxes on the contour of level $\ell+1$ to level $\ell$, and to use these fluxes on the corresponding interfaces of the coarse level in stead of the fluxes given by the hydrodynamical solver on the coarse level. We communicate and overwrite only fluxes at the frontier.

At a frontier between different levels, the flux $f$ of a coarser cell is therefore overwritten by the sum of the fluxes of the $4$ finer cells that share interface with the coarser cell. The fluxes are communicated at the end of a full cycle of the finer level. In the case of sub-cycling the flux to be communicated is the sum of the fluxes of each sub-cycle $f=f_{1}+f_{2}$.

To ensure concurrent fluxes between levels, we have slightly modified the structure of the source step described by Benítez-Llambay & Masset (2016), as we have split the so-called source step SubStep1 (which corresponds to the update of the velocity field under the action of the potential and fictitious forces) into two separate source steps of half the time step at the given level, one before the transport step and the other one after it. In stead of the original scheme:

[TABLE]

we have the modified scheme:

[TABLE]

where $S$ and $T$ stand for the source and transport step, respectively. The evaluation of fluxes is done with the partially updated velocity field $\boldsymbol{v}^{p}$, obtained upon the end of the first substep.

3.3 Sub-cycling and staging of communications

We present a brief description of our implementation of sub-cycling, showing where the communications are performed. Here we use $\Delta t$ to denote the time step of a sub-cycle on the finer level, we use the superscripts $p$ and $q$ to denote partially updated fields and $f$ for the fully updated fields, while we use subscripts to number the sub-cycles. The names of the different substeps are those of Stone & Norman (1992) and Benítez-Llambay & Masset (2016).

• •

Finer level $(\ell+1)$:

• –

First Sub-cycle:

SubStep1$(\Delta t/2)\rightarrow\boldsymbol{v}^{p}_{1}$.

• *

SubStep2,SubStep3$(\Delta t)\rightarrow\{\rho,e\}^{p}_{1}$

• *

FluxIntegration$(\boldsymbol{v}^{p}_{1},\Delta t)\rightarrow f_{1}$

• *

Transport$(f_{1})\rightarrow\{\rho,e\}^{f}_{1},\boldsymbol{v}^{q}_{1}$.

• *

SubStep1($\Delta t/2)\rightarrow\boldsymbol{v}^{f}_{1}$.

• –

Second Sub-cycle:

SubStep1($\Delta t/2)\rightarrow\boldsymbol{v}^{p}_{2}$.

• *

SubStep2,SubStep3$(\Delta t)\rightarrow\{\rho,e\}^{p}_{2}$

• *

FluxIntegration$(\boldsymbol{v}^{p}_{2},\Delta t)\rightarrow f_{2}$

• *

Transport$(f_{2})\rightarrow\{\rho,e\}^{f}_{2},\boldsymbol{v}^{q}_{2}$.

• *

SubStep1($\Delta t/2)\rightarrow\boldsymbol{v}^{f}_{2}$.

• –

CommFlux: $f_{1}+f_{2}\rightarrow$ Coarse level

• •

Coarse level $(\ell)$:

• –

First Cycle:

SubStep1($\Delta t)\rightarrow\boldsymbol{v}^{p}$.

• *

SubStep2,SubStep3$(2\Delta t)$

• *

FluxIntegration$(\boldsymbol{v}^{p},2\Delta t)\rightarrow f$.

• *

OverWriteFluxes$(f_{1}+f_{2})\rightarrow f$.

• *

Transport$(f)\rightarrow\{\rho,e\}^{f},\boldsymbol{v}^{q}$.

• *

SubStep1($\Delta t)\rightarrow\boldsymbol{v}^{f}$.

• *

CommUp$(\ell\rightarrow\ell+1)$: $\{\rho,\boldsymbol{v},e\}^{f}\rightarrow\{\rho,\boldsymbol{v},e\}^{f}_{2}$

• *

CommDown$(\ell\leftarrow\ell+1)$: $\{\rho,\boldsymbol{v},e\}^{f}\leftarrow\{\rho,\boldsymbol{v},e\}^{f}_{2}$

We then have concurrent velocities and concurrent fluxes at the end of a coarse cycle, as shown in figure 5.

3.4 Implementation of the dynamical friction setup

We have built a setup with Cartesian meshes where a perturber of mass $M$ and luminosity $L$ is kept fixed at $\boldsymbol{r_{0}}=(0,0,0)$. The perturber is immersed in a gas with an initially uniform density $\rho_{0}=1$, adiabatic sound speed $c_{s}=1$ and an adiabatic index $\gamma=1.4$. The initial energy is then $e_{0}=c_{s}^{2}\rho_{0}/\gamma(\gamma-1)$. The gas is advected at constant velocity $V$ in the $z$ direction and has a constant thermal diffusivity $\chi$. Taking advantage of the symmetry of the problem in the $x$ and $y$ directions, we reduce our spatial domain to the positive quadrant $x>0$, $y>0$. The ground or coarsest level $\ell=0$ sets the main domain with $x_{min}=0$, $x_{max}=l$, $y_{min}=0$, $y_{max}=l$, $z_{min}=-l$ and $z_{max}=l$, where the extension $l$ is chosen to fulfil a number of conditions that we specify below. Each nested mesh of level $\ell$ covers the sub-domain: $x_{\ell}\in[0,2^{-\ell}x_{max}]$, $y_{\ell}\in[0,2^{-\ell}y_{max}]$ and $z_{\ell}\in[2^{-\ell}z_{min},2^{-\ell}z_{max}]$.

The gravitational potential follows Plummer’s law:

[TABLE]

where $\epsilon$ is the softening length. As said in the introduction, we restrict ourselves in the present work to the study of perturbations with a small relative amplitude, so as to compare our results to those given by the linear analysis of MV17. As a consequence, we renounce to describe the dynamics within the perturber’s Bondi radius

[TABLE]

and adopt for the cells of the finest level a size comparable to or larger than $R_{B}$. A study involving much more resolved calculations showing the flow within the Bondi sphere will be presented in a forthcoming work. Because of this limitation which implies that the perturbation of density is not large even in the vicinity of the perturber, and because the perturber lies on the vertex of the mesh in all our calculations, we generally use $\epsilon=0$, unless explicitly specified.

We monitor the $z$ component of the force that the gas exerts onto the perturber

[TABLE]

where $r^{\prime}_{l}=|\boldsymbol{r}_{l}-\boldsymbol{r}_{0}|$ is the distance from the centre of cell $l$ to the perturber. In practice evaluating this sum over all the cells is done level by level, and care has to be taken to discard the values of coarse zones covered by a finer level.

Thermal diffusion is implemented by evaluating the energy fluxes given by Fourier’s law. Upon multiplication by the surface area of the interface they are associated to (either $S_{x}$, $S_{y}$ or $S_{z}$), they read:

[TABLE]

In the above expressions, we only write indices which differ respectively from $i$, $j$ or $k$, for the sake of legibility. The energy is then updated as:

[TABLE]

where $V$ represents the volume of the cell.

We implement the release of heat by the luminous point-like perturber in a much similar manner to that of Benítez-Llambay et al. (2015): since the perturber falls exactly at a mesh vertex, the energy $L\Delta t$ released by the perturber over a time step is split equally between the eight adjacent zones and used to increment their internal energy.

The boundary conditions of our setup are as follows:

• •

In $x_{min}$ and $x_{max}$:

• –

For $\rho,v_{y},v_{z},e$, we use symmetric boundary conditions.

• –

For $v_{x}$, we use antisymmetric boundary conditions.

• •

In $y_{min}$ and $y_{max}$:

• –

For $\rho,v_{x},v_{z},e$, we use symmetric boundary conditions.

• –

For $v_{y}$, we use antisymmetric boundary conditions.

• •

In $z_{min}$ and $z_{max}$:

• –

For $\rho,v_{z},e$, we use in/out-flow boundary conditions (we set these variables respectively to $\rho_{0},V,e_{0}$ in the ghost zones).

• –

For $v_{x},v_{y}$, we use symmetric boundary conditions.

4 results

All simulations presented here were run up to $t_{\mathrm{end}}=100\tau$, where $\tau$ is given by Eq. (9) and depends on the wind’s speed $V$. We chose this duration to ensure that the flow relaxes to a nearly steady state. For the calculations involving the net force, we take care that the acoustic sphere triggered by the insertion of the perturber in the gas does not reach the edge of the coarser mesh over the duration of the simulation444When this sphere is reflected on the outer edge of the mesh, the net force exerted on the perturber changes and does no longer represent the dynamical friction.. Calculations focusing specifically on one of the thermal forces are obtained by subtracting the outcomes of two different calculations555A calculation with a luminous perturber and a calculation with a non-luminous perturber for the heating force, and a calculation with thermal diffusion and an adiabatic one for the cold thermal force.. Our calculations being in the nearly linear regime, a reflection of the acoustic sphere on the edge of the computational domain does not have a measurable impact on the thermal forces. For reasons of computational cost we occasionally adopt a smaller mesh, on the edges of which the acoustic sphere is reflected in the course of the simulation, when focusing exclusively on thermal forces.

We quote the mass values in units of $c_{s}^{3}/\sqrt{G^{3}\rho_{0}}$, those of thermal diffusivity in units of $c_{s}^{2}/\sqrt{G\rho_{0}}$ and those of luminosity in units of $c_{s}^{5}/G$. For the sake of legibility, we will not write explicitly the units after each value that we quote. They correspond to the actual values when $c_{s}=1$, $G=1$ and $\rho_{0}=1$, which is our choice in the simulations. This scale free unit system allows to capture indistinctly the dynamical friction on a black hole in a gaseous medium or on a proto-planetary embryo in a proto-planetary disc.

4.1 Heating force

We undertook a set of simulations spanning 80 values of the Mach number $\cal{M}$ $\in[0.02,5]$ in a geometric sequence, with thermal diffusivity $\chi=0.1$ and a perturber of mass $M=10^{-3}$. For each Mach number, two calculations were run: one with luminosity $L=2\cdot 10^{-4}$ and another one with a cold perturber. The heating force is obtained by subtracting the forces of the cold runs from the forces of the hot runs. In this set of calculations we used a box with $l=100\lambda$ thus allowing acoustic waves to bounce off the boundaries for low Mach numbers (in which case $100\tau c_{s}>100\lambda$). The heating force being the result of a subtraction, it is impervious to this acoustic bounce for flows with nearly linear perturbations such as the ones we consider here. We use here $6$ mesh levels in total, including the ground one, each mesh with sizes $N_{x}=16$, $N_{y}=16$ and $N_{z}=32$. Fig. 6 shows the heating force as a function of the Mach number. It confirms expectations for the two regimes studied by MV17: a finite, asymptotic value at low Mach number, and a force inversely proportional to the speed squared in the supersonic regime. The heating force is positive in all cases. We study in more details these two regimes below.

4.1.1 Low Mach number

In the limit of a low Mach number, the relative contribution to the net heating force from all the material within a radius $r$ is (MV17):

[TABLE]

where $r^{\prime}=r/2\lambda$. The fraction of the net force that we can obtain in a simulation is therefore $\approx h(r^{\prime}_{\mathrm{max}})-h(r^{\prime}_{\mathrm{min}})$, where $r^{\prime}_{\mathrm{max}}$ corresponds to the maximal length scale (the box size) and $r^{\prime}_{\mathrm{min}}$ corresponds to the minimal length scale (the resolution of the finest level).

For the simulations presented in figure 6, the finest resolution amounts to $\Delta x=\frac{25}{128}\lambda$, therefore $r^{\prime}_{min}=\frac{25}{256}\approx 9.8\cdot 10^{-2}$ and $r^{\prime}_{max}=l/2\lambda=50$. We therefore expect to get $h(r^{\prime}_{max})-h(r^{\prime}_{min})\approx 92$% of the full heating force due to our finite length scale range. This is consistent with the force that we obtain in the low Mach number limit, as can be seen in the inset of Fig. 6.

In the low Mach number limit, the density perturbation is described by Eq. (6). In Fig. 7 we show a map of the density perturbation arising from the heat release and in Fig. 8 we plot three different cuts of the density perturbation: along, against and perpendicularly to the direction of advection ($z$) for the case ${\cal M}=0.023$ with a ground level of size $64\times 64\times 128$ and $11$ levels of refinement of size $32\times 32\times 64$. Here the first level covers a quarter of the ground level in a given direction, in stead of half of the ground level, which allows to cover a larger domain. We observe a good agreement between numerical results and the analytic expectation. Discrepancies are observed only near the boundaries, where the perturbation of density is negligible, and at distances to the perturber comparable to the resolution of the finest level.

4.1.2 High Mach number

As said above, our setups are such that $l/\lambda$ is fixed, and so is therefore $r_{\mathrm{min}}V/4\chi\equiv r_{\mathrm{min}}/(\lambda\times 4\gamma)$. The only dependence on velocity of the heating force in the supersonic regime, given by Eq. (10), is thus through the denominator in $V^{-2}$, which is in agreement with the plot of Fig. 6. We further investigate the dependence of that force on the resolution, or $r_{\mathrm{min}}$, for a fixed Mach number in the supersonic regime. We reproduce a number of times the same calculation, adding each time one level of refinement. The smallest scale is therefore halved from one simulation to the next. In the limit where $r_{\mathrm{min}}V/4\chi$ is sufficiently small, we see from Eq. (13) that the variation of the force between two successive calculations should be:

[TABLE]

This increment is in excellent agreement with the results depicted in Fig. 9.

The density perturbation arising from the release of heat exhibits the conical shape typical of the supersonic case, as shown in Fig. 10. The perturbation, however, is not strictly bound to the interior of the Mach cone in the vicinity of the apex: thermal diffusion has a typical spread in $\sqrt{\chi t}$ and is therefore “infinitely fast” in the vicinity of the perturber, allowing it to beat the acoustic perturbations that define the cone.

4.1.3 Intermediate Mach number values

The analytic results of MV17 have been established both for the low and high Mach number regimes. No general expression has been given for intermediate values of the Mach number. Fig. 1 shows that the ${\cal M}^{-2}$ dependence in the supersonic regime is valid almost until ${\cal M}=1$. In Fig. 7 we compare the value of the heating force at the transonic point to the value given by Eq. (10) with $V=c_{s}$, as a function of the minimal distance $r_{\mathrm{min}}$. The analytic expectation is in reasonable agreement with the measured values. The ratio $\xi$ of the force value at the transonic point to its value at low Mach number is therefore:

[TABLE]

When the plume’s size is much larger than the minimal length scale $r_{\mathrm{min}}$, this ratio is larger than one so that the force has a peak at the transonic point. Only when the minimal length scale is comparable to the plume’s size does the peak vanish. In that case, by virtue of Eq. (23), the asymptotic value at low Mach number is also significantly below its nominal value.

4.2 Dynamical friction in presence of thermal diffusion

In this section we study the force acting on a non-luminous perturber, when heat exchange can occur between neighbouring fluid elements, i.e. when we relax the adiabatic assumption of Ostriker (1999). In a first step we evaluate the correction to the force due to thermal diffusion using of linear perturbation theory and in a second step we check that expression using numerical simulations. In this whole section we therefore consider a cold perturber, the only difference with respect to the setup considered by Ostriker (1999) being the introduction of thermal diffusion.

4.2.1 Linear analysis of the cold thermal force in the limit of low Mach numbers

We start by linearising our governing equations [Eqs. (1) to (3)] which read:

[TABLE]

where prime quantities denote perturbations (e.g. $\rho^{\prime}=\rho-\rho_{0}$). For the sake of legibility we hereafter drop the primes. Using the same notation and same convention for Fourier transforms as MV17 we get, upon Fourier transforming Eqs. (26) to (28) in space and time:

[TABLE]

We can make use of Poisson’s equation, which reads:

[TABLE]

where $\rho_{p}(\boldsymbol{r},t)$ is the density of the perturber. Using Eqs. (29) and (30), it follows that:

[TABLE]

which upon substitution in Eq. (31) yields the Fourier transform of the density perturbation:

[TABLE]

where $c_{s}=\sqrt{\gamma p_{0}/\rho_{0}}$ is the adiabatic sound speed. When $\chi=0$ this expression can be simplified into:

[TABLE]

which is equivalent to Eq. (4) of Ostriker (1999), and requires a time-dependent solution involving retarded Green’s functions to give the drag force. In the general case, we can separate Eq. (34) into two components:

[TABLE]

where

[TABLE]

is the new component that arises in the presence of thermal diffusion. This thermal component appears as the convolution product of the perturber’s density $\rho_{\mathrm{thermal}}$ by a Green’s function $K(\boldsymbol{r},t)$, the latter being the inverse Fourier transform of the ratio by which $\tilde{\rho}_{p}$ is multiplied in Eq. (37). We now consider a perturber with a trajectory at constant speed $V$ along the $z$ axis that passes through $z=0$ at $t=0$. We evaluate the field $\rho_{\mathrm{thermal}}(\boldsymbol{r},t=0)$ in order to evaluate the force that it induces on the perturber in $\boldsymbol{r}=0$. The perturber’s density field has the expression:

[TABLE]

where ${\cal H}$ is Heaviside’s step function, which we introduce to consider only the past ($t<0$) trajectory of the perturber. The convolution of $\rho_{p}$ by the Green’s function $K$, which reads:

[TABLE]

then takes the particularly simple form:

[TABLE]

The force exerted on the perturber at $t=0$ is along the $z$ axis and has the expression:

[TABLE]

Only the norm $k$ of the wave vector features in Eq. (37) so that the Green’s function has spherical symmetry. The inner triple integral in Eq. (41) corresponds to the force exerted in $\boldsymbol{r}=0$ by a density perturbation equal to a Green’s function centred on the location $(0,0,z^{\prime})$. Since this function has spherical symmetry about that point, one may use Gauss’ flux theorem to express this integral simply in terms of the mass enclosed within the sphere centred on $(0,0,z^{\prime})$ with radius $-z^{\prime}$. This mass can be easily worked out in the limit of a small Mach number. In that case the Green’s function is evaluated at time $t=-z^{\prime}/V$ for radii $|\boldsymbol{r}|<|z^{\prime}|$, hence for radii $|\boldsymbol{r}|\ll tc_{s}$. In this limit we have in Fourier space $kc_{s}\gg\omega$ and we can simplify Eq. (37) as:

[TABLE]

with

[TABLE]

In real space we have therefore:

[TABLE]

In these conditions, the integral of $K(\boldsymbol{r},t)$ on a sphere of radius $R$ has the expression:

[TABLE]

Using Gauss’s flux theorem we can therefore write:

[TABLE]

Substituting this result in Eq. (41), we eventually obtain:

[TABLE]

and, generalising to an arbitrary direction of motion:

[TABLE]

This force is a drag force, like the dynamical friction in an adiabatic gas, as it is opposed to the motion. It has two remarkable properties: it is independent of the thermal diffusivity of the gas and of the perturber’s velocity. The heating force has the same last property in the low Mach number regime.

It is also straightforward to realise, if we substitute in the analysis of section 3.1 of MV17 the convective derivative $V\partial_{z}$ with $\partial_{t}$, that the Fourier transform of the density perturbation arising from heat release on one hand, and the Fourier transform of the density perturbation arising from thermal diffusion [our Eqs. (42) and (43)] on the other hand, have same form, except for a change of sign. In the limit of low Mach numbers, the perturbation of density $\rho_{\mathrm{thermal}}$ has therefore the same shape as that of the hot, low-density plume associated to the heating force and corroborated by the numerical simulations presented in section 4.1.1, except for a change of sign: it is here a cold and dense plume. This plume lagging the perturber, it is intuitive that it exerts an additional drag on it.

Following this analogy, we can further say that $\rho_{\mathrm{thermal}}$ has, in the low Mach number limit, the same form as the density perturbation that would be triggered by a massless heat sink with luminosity $-L_{c}$. The luminosity $L_{c}$ can be found by equating Eqs. (8) and (48), which yields:

[TABLE]

This luminosity has the same expression as that found in another context by Masset (2017) for the cold thermal torques exerted on a low-mass planet on a circular orbit in a gaseous protoplanetary disc.

From the above we expect a perturber with a luminosity $L_{c}$ to be subjected to a heating force that nearly cancels the cold thermal force in the low Mach number regime. Such a perturber should therefore experience a drag force comparable to the force in an adiabatic gas.

4.2.2 Numerical corroboration

In order to confirm the magnitude and direction of the cold thermal force, we ran simulations of a non-luminous perturber for different values of the thermal diffusivity ($\chi=0.01,0.05$ and $0.1$), spanning 44 values of the perturber’s velocities in the subsonic regime ($\cal{M}$ $\in[0.02,0.7]$). The dimensions of the box as well as the final time were adjusted for $\chi=0.1$ (the most demanding case). Here we no longer restrict ourselves to the study of the heating force. It is therefore important to avoid the reflection of the acoustic sphere triggered at $t=0$ on the edges of the computational domain. We therefore chose $l=3200\lambda$ which prevents this reflection (except for the first two velocities where a bounce occurs near the end of the simulations). We have a ground level with size $64\times 64\times 128$, and $11$ levels of refinement with the first level covering a quarter of the ground level. Each level $\ell>0$ has size $32\times 32\times 64$. With this system of meshes we expect to retrieve $\approx 99.1\%$, $98.4\%$ and $92.3\%$ of the net force666Assuming the cold thermal force to have the same dependence on the largest and smallest scales as that given by Eq. (23) for the heating force. respectively for the cases of $\chi=0.1,0.05$ and $0.01$. The results of these runs are shown in Fig. 12. They confirm the analytical expectations with a large degree of accuracy. The net force in runs with thermal diffusion is found to tend toward a finite, negative value in agreement to within a few percents with the prediction of Eq. (48). This value is largely independent of the value of the thermal diffusivity. The minute variation from one run to the next with respect to this parameter can be attributed to changes in length scale coverage and to the onset of non-linearities in the flow near the perturber for the smallest value of $\chi$.

We also show the perturbations of density and temperature in the diffusive and adiabatic cases in Fig. 13, and the difference of these fields between these two cases. It confirms that the introduction of thermal diffusion induces a positive perturbation of density and a negative perturbation of temperature with respect to the adiabatic case.

4.3 Net force in the general case

Having studied in separation the heating force and the cold thermal force, we are now in a position to study the net force on a perturber with mass $M$ and luminosity $L$ in a gas with thermal diffusion, when the flow is nearly linear. The net force is then given by the expression of Fig. 1.

4.3.1 A simple expression of the thermal forces at low Mach number

If we define:

[TABLE]

we can recast the expression of the thermal forces in the regime of low Mach number corresponding to Eqs. (48) and (8) respectively as:

[TABLE]

and

[TABLE]

so that

[TABLE]

We also note that the adiabatic force is, in terms of $F_{0}$, expressed as:

[TABLE]

for all Mach numbers, where $I$ is given by Eqs. (14) or (15) of Ostriker (1999).

4.3.2 Numerical analysis

In Fig. 14 we present the results for different sets of simulations using 11 levels of refinement to reach $r^{\prime}_{min}=0.01$ in a box of $r^{\prime}_{max}=1600$, therefore retrieving $99\%$ of the net heating force. We show the total force for the adiabatic case ($\chi=0$), and for cases with thermal diffusivity $\chi=0.1$ for three different values of the luminosity: $L=0$, $L=L_{c}$ and $L=2L_{c}$. We confirm our expectation that the net force nearly coincides with the adiabatic value of Ostriker (1999) when $L=L_{c}$ in the low Mach number regime (see last paragraph of section 4.2.1). Our analysis of the cold thermal force of section 4.2 was limited to the low Mach number regime. We have not worked out an expression of that force in the supersonic regime. However, our simulations seem to indicate that in this regime the cold thermal force has a negligible role, since the net force without heating (blue curve in the electronic version of this manuscript) tends to the adiabatic force (green curve in the electronic version of this manuscript). We also find that the net force is positive, and corresponds to a thrust on the perturber, when the luminosity is in excess of $L_{c}$ and the Mach number is sufficiently low (e.g. below $\approx 0.48$ when $L=2L_{c}$).

Masset & Velasco Romero (2017) argued that, since both the adiabatic force and the heating force follow a $V^{-2}$ law in the supersonic regime and have opposite signs, the net force may be positive at all speeds in the supersonic regime when the perturber’s luminosity is sufficiently large. Their analysis disregarded the cold thermal force and assumed the force on a non-luminous perturber to coincide with the adiabatic estimate of Ostriker (1999). Our present findings, however, suggest that the cold thermal force is negligible in the supersonic regime, so that their expectation should be correct. We have therefore undertaken additional sets of calculations with higher luminosities, presented in Fig. 15. We see that for a luminosity $L=4L_{c}$, the net force is positive essentially over the whole subsonic domain, while it is always positive for $L=8L_{c}$, for the set of parameters adopted. This behaviour is consistent with the density profile depicted by Fig. 16 in which we observe a deficit of material behind the perturber inside the Mach cone.

5 Discussion and conclusions

We have conducted extensive numerical simulations in order to check the impact of thermal effects on the dynamical friction onto a massive perturber, both with and without intrinsic luminosity (i.e. with or without radiative feedback). Owing to the wide interval of length scales required to capture accurately thermal effects, we had to resort to nested meshes, that we implemented in the GPU version of the FARGO3D code. As a side feature, the present manuscript is intended to stand as a concise reference to our nested mesh implementation.

Our simulations confirmed with a high degree of accuracy the analytic predictions of Masset & Velasco Romero (2017). Namely, when the perturber is luminous, the extra force arising from the density perturbation due to the radiative feedback, or heating force, is found to behave as expected in the regimes of low and high Mach numbers. In the regime of low Mach numbers, it tends toward a finite value when the velocity tends to zero. This somehow counter-intuitive behaviour occurs owing to the following scalings. The region that gives rise to the net force has a typical size $\sim\lambda$ and therefore a volume $V\propto\lambda^{3}$. The perturbation of density in this region scales as $\rho\propto{\lambda}^{-1}$, therefore leading to a perturbation of mass $\delta M\propto\lambda^{2}$. It then follows that the net force arising from this volume scales as $GM\delta M/\lambda^{2}$ and is therefore independent of the cut-off distance $\lambda$. The only effect of changing the velocity on the perturbation of density is to change the scale $\lambda$. It therefore follows that the heating force does not depend on the velocity. However, as the perturber’s velocity tends to zero, the cut-off distance increases (as diffusion can go a long way against the wind) and it takes more time to establish the perturbed density field. In the low Mach number regime, the acceleration due to the heat release is a function of $t/\tau$ only, where $\tau$ is defined by Eq. (9). This function tends toward the asymptotic value of Eq. (8) and $\tau$ tends to infinity when the Mach number vanishes.

In the high Mach number regime the heating force decays as $V^{-2}$, and exhibits a logarithmic dependence on the smallest length scale. In all cases the heating force is along the direction of motion and yields a thrust on the perturber.

The present work also shows that the gravitational drag on a non-luminous perturber, when thermal diffusion is present, is significantly at odds with the expression worked out by Ostriker (1999) for an adiabatic gas. The offset between the two expressions, that we dub the cold thermal force, shares a number of similarities with the heating force in the subsonic regime. In particular, it tends towards a finite value in the limit of a vanishing velocity. The interpretation of this new force component is very similar to that of the cold thermal torque discussed by Masset (2017) in a different context (that of low mass protoplanets in a gaseous protoplanetary disc). The introduction of thermal diffusion flattens out the temperature peak that would otherwise surround the perturber if the gas were adiabatic. Allowing thermal diffusion is therefore tantamount to introducing an additional, negative perturbation of temperature in the vicinity of the perturber. Pressure balance in the low Mach number regime then implies that this negative perturbation is accompanied by a positive perturbation of density. Masset (2017) argues that when the perturber has a very subsonic velocity with respect to the ambient gas, the introduction of thermal diffusion yields a perturbation on the gas that is equivalent to the perturbation that would be induced by a massless heat sink with luminosity $-L_{c}$, where $L_{c}$ is given by Eq. (49), regardless of the geometry of the background flow (whether it is sheared, like in a protoplanetary disc, or not, like in the situation considered here). This is exactly what we find in the present analysis and the simulations presented here constitute the first quantitative confirmation of this effect. One potentially important consequence of this effect is the eccentricity damping of cold embryos in protoplanetary discs. It should be much more vigorous than expected from adiabatic calculations. This has been noticed in numerical simulations by Eklund & Masset (2017), but not in a quantitative manner, owing to resolution constraints.

The cold thermal force being finite at vanishing velocity, it behaves like solid or dry friction, and halts the perturber in a finite time. As the response time of the plume tends to infinity in the limit of vanishing velocity, this can only be true, however, if the stopping distance is a fraction of the cut-off distance $\lambda$ for the initial velocity. It is easy to check that a sufficient condition for that is $L_{c}\gg c_{s}^{5}/G$.

The finite, negative force in the limit of a vanishing velocity for a non-luminous perturber that we found here had not been considered by Masset & Velasco Romero (2017). Their conclusions need therefore to be amended. A perturber in slow motion is accelerated only if its luminosity is larger than $L_{c}$. Otherwise, the net force acting on it is negative, and its velocity decays. For the two fiducial disc models considered by these authors, the critical luminosity for a planetary embryo with mass $\chi c_{s}/G$ coincides with the accretion luminosity for a mass doubling time of roughly $100$ kyrs. For the disc parameters considered by Benítez-Llambay et al. (2015) and Eklund & Masset (2017), the accretion luminosity for a mass doubling time of $100$ kyrs is twice the critical luminosity of an embryo with critical mass $\chi c_{s}/G$. Low mass embryos subjected to fast pebble accretion (Lambrechts & Johansen, 2012), for which mass doubling times largely shorter than $100$ kyrs can be expected, should therefore be subjected to a heating force largely in excess of the cold thermal force.

Although we have not worked out an analytical expression of the cold thermal force in the supersonic regime, our numerical experiments suggest it to be negligible. A tentative explanation could be as follows: the cold force primarily arises from the flattening, by thermal diffusion, of the temperature peak associated to a nearly hydrostatic envelope in the slow regime. No such envelope can build up in the supersonic regime, so that the introduction of thermal diffusion in that regime induces much less significant perturbations in the perturber’s vicinity.

We have found that a sufficiently large luminosity of the perturber warrants a net positive force at all Mach numbers. We comment that the gravitationally induced drag diverges logarithmically as time goes on in the supersonic regime while the heating force remains bound, so that ultimately the net force will revert sign in this regime and will become negative, unless the medium is bound, such as in a thin protoplanetary disc (Muto et al., 2011). In that case, the net force may remain positive at all times and the perturber’s velocity increases as $t^{1/3}$. The hot plume’s size being inversely proportional to the perturber’s velocity, the acceleration ceases when the plume size becomes comparable to the perturber’s physical size, or, if thermal diffusion is effected by radiative transfer, when it becomes comparable to the photons’ mean free path, whichever occurs first.

In the present analysis we have restricted ourselves to weakly perturbed flows, as we have not tried to resolve the flow within the Bondi sphere of the perturber. High resolution simulations resolving the flow in the Bondi sphere and evaluating its impact on the thermal forces will be presented in a forthcoming work.

Appendix A Sub-cycling criterion

In order to show how the time stepping schedule can be optimised, we take as an example a low-mass planet embedded in a nearly Keplerian protoplanetary disc, so that the gas is nearly at rest in the planet’s vicinity. The standard sub-cycling technique may not be the most appropriate in this situation: owing to the shear, the flow velocity (which imposes the strongest constraint on the time step, in a geometrically thin disk) is smaller for the higher resolution meshes, which have a smaller radial extent from corotation. For the sake of definiteness consider a mesh with one level of refinement only in the planet vicinity. Assume that the CFL condition, evaluated separately on the coarse mesh, yields a time step limit $\delta t_{0}$, and evaluated separately on the fine mesh, yields a time step limit $\delta t_{1}=K\delta t_{0}$, where we assume $1/2<K<1$. For simplicity we assume that the fine mesh has the same number of zones as the coarse one, so that an integration on either mesh has same computational cost $C_{0}$. We can now consider two possibilities:

• •

We perform two timesteps on the finer mesh for one timestep on the coarser mesh (sub-cycling). The global time step is therefore limited by $\delta t_{0}$, since $K>1/2$, and the total cost of an integration on the whole domain is $3C_{0}$.

• •

We use the same time step for both meshes (flat case). The global time step is therefore limited by $\delta t_{1}$ and the total cost of an integration on the whole domain is $2C_{0}$.

We naturally wish to maximise the time elapsed in the simulation for a given amount of wall clock time, so we search the procedure that has the largest time step to cost ratio. For the first case contemplated here, it is $\delta t_{0}/3C_{0}$, whereas for the second case it is $K\delta t_{0}/2C_{0}$. We see that if $K>2/3$, the flat procedure (no sub-cycling) is more efficient than the other one. As anticipated above, it is in the case where the timestep limit does not differ much on the fine and coarse meshes that a flat procedure is desirable. This simple example can be generalised to a case of $n$ levels of refinement, each corresponding to a computational cost $(C_{\ell})_{0\leq\ell<n}$, and each having an individual time step limit $(\delta t_{\ell})_{0\leq\ell<n}$. A time step at a given level $\ell$ can correspond either to one (flat case) or two (sub-cycling) time steps at level $\ell+1$, independently of how many timesteps are performed at level $\ell+2$ for one time step at level $\ell+1$. We therefore construct a sequence $(\xi_{\ell})_{0\leq\ell<n-1}$, such that $\xi_{\ell}=1$ if the same time step is used for levels $\ell$ and $\ell+1$, and $\xi_{\ell}=2$ otherwise. The timestep on the whole domain that satisfies the CFL condition for all levels is:

[TABLE]

where the product is meant to be $1$ when $\ell=0$, and the cost of one full time step is

[TABLE]

We seek the sequence $(\xi_{m})_{0\leq m<n-1}\in\{1,2\}^{n-1}$ for which the ratio $\Delta t/C$ is maximal. A given sequence can be regarded as a binary word of length $n-1$ composed of 1’s and 2’s. We therefore have to evaluate Eqs. (55) and (56) for the $2^{n-1}$ possible cases. This optimisation is performed at the beginning of each coarse time step, with a negligible computational overhead. For the setup considered in the present paper, this technique falls back to using $\xi_{\ell}=2$ for $0\leq\ell<n$.