<figure><img src="image_1.png" /></figure>

HAL
open science

# Self-gravity in thin-disc simulations of protoplanetary discs: The smoothing length rectified and generalised to bi-fluids

S. Rendon Restrepo, P. Barge

► To cite this version:

S. Rendon Restrepo, P. Barge. Self-gravity in thin-disc simulations of protoplanetary discs: The smoothing length rectified and generalised to bi-fluids. Astronomy & Astrophysics, 2023, 675, 10.1051/0004-6361/202346178. insu-04479077

HAL Id: insu-04479077

https://insu.hal.science/insu-04479077

Submitted on 5 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Self-gravity in thin-disc simulations of protoplanetary discs:
The smoothing length rectified and generalised to bi-fluids

S. Rendon Restrepo<fnref n="1" />,<fnref n="2" />,<fnref n="3" />* and P. Barge<fnref n="1" />

Received 17 February 2023 / Accepted 5 April 2023

## ABSTRACT

**Context.** To mimic the evolution of protoplanetary discs (PPDs), 2D simulations that incorporate self-gravity must introduce a softening prescription of the gravitational potential. When the disc is only composed of gas, the smoothing length is proportional to the scale height of the gas. On the other hand, when a dust component is included, the question arises as to whether the smoothing length approach can still be used to not only quantify the dust self-gravity, but its gravitational interaction with gas as well.

**Aims.** We identified grey areas in the standard smoothing length formalism for computing self-gravity in PPDs entirely made up of gas. Our goal is to revisit the smoothing length approach, which can then be generalised to two phases, whereby the dust component may be considered as a pressureless fluid.

**Methods.** We used analytical developments to approximate the vertically averaged self-gravity when the smoothing length is not assumed to be constant, but by taking a spatial function instead.

**Results.** We obtained an analytical expression for the space-varying smoothing length, which strongly improves the accuracy of the self-gravity computation. For the first time, this method is generalised to address bi-fluid interactions in a PPD: two additional smoothing lengths are proposed for featuring an isolated dusty disc and gas-dust self-gravity interactions. On computational grounds, we prescribe the use of tapering functions for the purpose of avoiding numerical divergences. We also checked that our method continues to be compatible with standard fast Fourier transform algorithms and evaluated computational costs.

**Conclusions.** Our space-varying smoothing length allows us to: (i) solve the contradictions inherent in the constant smoothing length hypothesis; (ii) fit the 3D vertically averaged self-gravity with a high level of accuracy; and (iii) render it applicable to a bi-fluid description of PPDs with the use of two additional smoothing lengths. Such results are crucial to enable realistic 2D numerical simulations that account for self-gravity and are essential to improving our understanding of planetesimal formation and type I migration.

**Key words.** protoplanetary disks – planet–disk interactions – gravitation – hydrodynamics – methods: analytical – methods: numerical

## 1. Introduction

Despite huge developments in 3D numerical computations and the advent of progressively more sophisticated computational facilities, thin-disc (2D) simulations remain much less expensive and more widely used in the study of protoplanetary discs (PPDs). The 2D approximation relies on the vertical averaging of the 3D physical quantities and governing equations. When dealing with self-gravity (SG), the issue becomes more sensitive since the equations cannot be vertically integrated. In practise, in PPDs studies an equivalent smooth potential, known as the Plummer potential, is used to mimic (as much as possible) the vertically averaged SG force and, in this approximation, a smoothing length (SL) is introduced to account for the disc vertical stratification.

As expected, this approach applies to 2D studies of planet–disc interactions where different values have been suggested (Li et al. 2009; Dong et al. 2011). In particular, Masset (2002) and Huré & Pierens (2009) showed that in the planet case, the SL should be proportional to the scale height of the gas disc, whereas Müller et al. (2012) proposed that $\epsilon_g/H_g = 0.7$, where $\epsilon_g$ and $H_g$ are the gas SL and the pressure scale height, respectively. Müller et al. (2012) also explored the case of self-gravitating

gas discs and found that for vertically isothermal and stratified structures, we should instead use $\epsilon_g/H_g = 1.2$ to avoid a systematic overestimation of the SG term. They also showed that when accounting for vertical SG, a deviation in the vertical Gaussian distribution occurs and in this case, the SL is also proportional to the disc’s Toomre parameter.

The Müller et al. (2012) prescription clearly leads to small errors (<2%) at large distances, but it also has an important drawback (that is not, however, mentioned by the authors); namely, SG is strongly underestimated at short distances with 100% errors. Indeed, in agreement with Huré & Pierens (2009), the authors found that the accuracy in the approximation of SG terms is strongly improved if the SL is a space varying function (Müller et al. 2012, Fig. 13). However, despite a discernible curve shape and (for unknown reasons), the authors did not test analytical expressions that could adequately fit the numerical curve. Furthermore, they found that the best approximation is obtained when the SL vanishes at the singularity – this is a statement that is inconsistent with the divergence of the error in the SG computations (Figs. 12 and 13 of their paper). The reasons for this inconsistency are subtle mathematical details at the singularity that were not explained by the authors and which, to our knowledge, have gone unnoticed to date.

<footnote n="1">¹ Aix-Marseille Univ., CNRS, CNES, LAM, Marseille, France
e-mail: pierre.barge@lam.fr</footnote>

<footnote n="2">² Aix-Marseille Univ., CNRS, Centrale Marseille, IRPHE, Marseille, France</footnote>

<footnote n="3">³ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
e-mail: srendon@aip.de</footnote>
In a bi-fluid description of PPDs, a dust layer is embedded in the gas disc and the evolution of gas and dust are coupled by aerodynamic forces. In this case, the momentum equations for the gas and the dust, both contain two SG terms: one is attributed to the gas disc and the other to the dust layer. In other words, each fluid parcel (gas or dust) is submitted to the gravity of the gas disc and to the gravity of the dust sub-layer. In contrast to the gas disc, the dust layer has a thickness which is governed not by pressure but turbulent stirring in the vertical direction. The scale height of the dust layer can be estimated using the form: $H_d = \sqrt{\frac{\tilde{\alpha}}{\tilde{\alpha}+S_t}} H_g$ (Dubrulle et al. 1995; Weber et al. 2019), where $\tilde{\alpha} = \alpha/Sc_z$, $\alpha$ is the dimensionless constant featuring turbulent $\alpha$-viscosity (Shakura & Sunyaev 1973), $Sc_z$ is the vertical Schmidt number, and $S_t$ is the Stokes number. A vertical averaging over the dust layer faces the same numerical issues than averaging over the gas disc. These problems can be solved in the same way as for a single gas disc and a simple extrapolation provides the dust disc's SL, namely, $\epsilon_d \propto 1.2H_d$. Since the dust layer is much thinner than the gas disc (due to dust settling), the two SLs are very different from one another with $\epsilon_d \ll \epsilon_g$. Thus, any error between these two parameters may lead to an incorrect and significant underestimation of the dust layer contribution to the SG terms. The question becomes more complex when accounting for the crossed SG terms. At the present time, we do not possess any rigorous theoretical approach that would allow us to evaluate this gravitational bi-fluid interaction. We do, however, think that this problem can also be addressed with a SL and, intuitively, we can expect it to lie between the gas SL and the dust SL. To our knowledge, both issues raised here with respect to bi-fluids are new and have not been previously addressed in numerical and theoretical studies of PPDs.

The improvements we perform in this study are crucial to enabling realistic 2D numerical simulations when SG is taken into account and they are also particularly important in planetary formation theories to better understand the formation of planetesimals. Indeed, it makes it possible to study the gravitational fragmentation of the dust layer (Goldreich & Ward 1973), formation of coherent clumps under the streaming instability (Johansen & Youdin 2007), and gravitational clumping of dust particles trapped in a large-scale vortex (Barge & Sommeria 1995) or in a co-orbital trapping scenario. Indirectly, we found that this work could also have implications for studies of type I planet migration.

In this paper, our goal is to provide a method based on the SL formalism to accurately compute the SG terms in thin discs at small and large separations. We aim to introduce the substantial SG interaction of gas and dust, when the latter is considered as a pressureless fluid. In Sect. 2, we carry out a theoretical development from first principles which justifies the use of the Plummer potential. In Sect. 3, we explore the SG estimation through the prism of the SL paradigm. We correct and complete the Plummer potential formalism thanks to a spatially dependent SL and generalise it to the case of bi-fluids. In Sect. 4, we handle practical aspects of the study, such as the numerical divergences, computational costs, and the calculation of the corrected SG thanks to fast Fourier transforms (FFT) methods. Finally, in Sect. 5, we suggest a discussion followed by a conclusion.

## 2. Self-gravity terms for bi-fluid simulations

In this section, we set up the theoretical background for computing gas and dust SG contributions when solid material is considered as a pressureless fluid in a thin-disc approximation.

In this context, we also recall the role of the Plummer potential in 2D SG calculations and generalise this formalism to bi-fluids.

### 2.1. 2D approximation and formal derivation

In 3D, the SG force per unit volume exerted by the PPD on a gas and dust parcels are expressed as:

$$ 
\begin{aligned}
f_{3D}^{\text{g,tot}}(\mathbf{r}, z) &= -\rho_g(\mathbf{r}, z) [\nabla\Phi_g + \nabla\Phi_d], \\
f_{3D}^{\text{d,tot}}(\mathbf{r}, z) &= -\rho_d(\mathbf{r}, z) [\nabla\Phi_g + \nabla\Phi_d],
\end{aligned}
\quad (1) $$

where $\Phi_g$ and $\Phi_d$ are the distinct gravitational potentials of the gas and dust discs, respectively. For the sake of generality and concision, we denote both fluid phases as a and b, which reduces each term on the right-hand-side of Eq. (1) to:

$$ 
\begin{aligned}
f_{3D}^{a \to b}(\mathbf{r}, z) &= -\rho_b(\mathbf{r}, z) \nabla \Phi_a \\
&= -G \rho_b(\mathbf{r}, z) \\
&\qquad \iint_{\text{disc}} \int_{z'=-\infty}^{+\infty} \frac{\rho_a(\mathbf{r}', z')}{\|\mathbf{r} - \mathbf{r}'\|^2 + (z - z')^2} e_l \, d^2\mathbf{r}' \, dz',
\end{aligned}
\quad (2) $$

where $f_{3D}^{a \to b}$ is force per unit volume the a-disc exerts on an elementary b-fluid element. The density of phase a is noted as $\rho_a$ and $e_l = [\mathbf{r} - \mathbf{r}' + (z - z')e_z] / \sqrt{\|\mathbf{r} - \mathbf{r}'\|^2 + (z - z')^2}$. Assuming a vertical hydrostatic equilibrium and an isothermal approximation in the vertical direction for gas, the volume density can be written as $\rho_g(\mathbf{r}, z) = \rho_{0,g}(\mathbf{r}) e^{-\frac{1}{2}z^2/H_g(\mathbf{r})}$, where $H_g(\mathbf{r})$ is the gas pressure scale height. We assumed the same vertical Gaussian profile for dust density than for gas, but this time, the vertical equilibrium is instead governed by turbulent stirring, which sets a different scale height for dust: $H_d(\mathbf{r})$. This permits us to write both surface densities as:

$$ \Sigma_a(\mathbf{r}) = \int_{z=-\infty}^{\infty} \rho_a(\mathbf{r}, z) \, dz = \sqrt{2\pi} H_a(\mathbf{r}) \rho_{0,a}(\mathbf{r}). \quad (3) $$

The 2D analogue of the SG force is simply obtained integrating the 3D SG force in the vertical direction:

$$ 
\begin{aligned}
f_{2D}^{a \to b}(\mathbf{r}) &= \int_{z=-\infty}^{\infty} f_{3D}^{a \to b}(\mathbf{r}, z) \, dz \\
&= -G \rho_{0,b}(\mathbf{r}) \iint_{\text{disc}} \rho_{0,a}(\mathbf{r}') s \, e_s \\
&\qquad \left( \int_{z,z'=-\infty}^{+\infty} \frac{e^{-\frac{1}{2}z'^2/H_g(\mathbf{r})} e^{-\frac{1}{2}z'^2/H_a(\mathbf{r}')}}{(s^2 + (z-z')^2)^{3/2}} \, dz \, dz' \right) d^2\mathbf{r}',
\end{aligned}
\quad (4) $$

where $s = \|\mathbf{r} - \mathbf{r}'\|$ is the separation (or mutual distance) between two fluid elements and $e_s = (\mathbf{r} - \mathbf{r}')/s$. In this 2D approximation, it is implicitly assumed that the disc is symmetric with respect to the $z=0$ plane, which allows us to naturally cancel the vertical component of the above force during the integration. After substituting the variables, we finally obtain:

$$ f_{2D}^{a \to b}(\mathbf{r}) = -\frac{G}{\pi} \Sigma_b(\mathbf{r}) \iint_{\text{disc}} \frac{\Sigma_a(\mathbf{r}')}{H_g(\mathbf{r}) s} L_{sg}^{ab}(d_g, d_b, \eta) e_s \, d^2\mathbf{r}', \quad (5) $$
where

$$L_{\text{sg}}^{ab}(d_g, d_b, \eta_{ab}) = \frac{1}{2} \frac{d_b^3(\mathbf{r})}{d_g(\mathbf{r})} \iint_{u,v=-\infty}^{+\infty} \frac{e^{-\frac{u^2}{2}} e^{-\frac{v^2}{2}}}{[d_b(\mathbf{r})^2 + (u - \eta_{ab}v)^2]^{3/2}} du dv \quad (6)$$

is a normalised quantity that we refer to as the 'self-gravity force correction' (SGFC)<fnref n="1" />, with respect to the 3D case, $d_b(r) = s/H_b(r)$ is the separation normalised with respect to the scale height of fluid b, and $\eta_{ab} = H_a(r')/H_b(r)$ is the a-to-b scale height ratio. In the following, we use the notation $d_b = d_b(r)$, except when a distinction is necessary. We want to highlight that in our reasoning, $\eta_{ab}$ is a spatial varying quantity, but, henceforth, we prefer to use:

$$\eta_{ab} = \langle H_a(r') \rangle / \langle H_b(r) \rangle, \quad (7)$$

where $\langle \rangle$ stands for space averaging over the 2D disc. This simplifies the next theoretical developments. We note that the scale heights used in this paper could be time-varying functions, but for concision, the time dependence is not explicit in the equations.

The numerical determination of SG in thin-disc and bi-fluid simulations mainly relies on the computation of the SGFC, but no analytical expression of this integral has been found in terms of standard mathematical functions. This is why a Plummer potential is commonly used to most optimally approximate this integral.

## 2.2. Fitting self-gravity terms thanks to smoothing lengths

Both vertical averagings in Eq. (6) cannot be performed analytically, but (in practice) we can introduce a Plummer potential, $\Psi_{\epsilon}^{ab}$, into the 2D gravitational potential, $\Phi_{\epsilon}^{ab}$, which permits us to approach the 2D SG force defined in Eq. (5):

$$f_{2D,\epsilon}^{a \to b}(\mathbf{r}) = -\frac{\Sigma_b(\mathbf{r})}{\pi} \nabla_s \Phi_{\epsilon}^{ab}(\mathbf{r}) \\ \text{with } \Phi_{\epsilon}^{ab}(\mathbf{r}) = -G \iint_{\text{disc}} \Sigma_a(\mathbf{r}') \Psi_{\epsilon}^{ab}(s) d^2\mathbf{r}', \quad (8)$$

where:

$$\Psi_{\epsilon}^{ab}(s) = \pi / (s^2 + \epsilon_{ab}^2)^{1/2}, \quad (9)$$

and $\epsilon_{ab}$ is the SL between phases a and b. Usually, $\epsilon_{ab}$ is assumed to be constant, so that:

$$\|\nabla \Psi_{\epsilon}^{ab}\| = \pi s / (s^2 + \epsilon_{ab}^2)^{3/2}. \quad (10)$$

For a gas disc, Müller et al. (2012, Fig. 13) have shown that the SL that gives the best fits to the SGFC is a spatial function of $d_g$; however, this is, in fact, inconsistent with Eq. (10). Indeed, the additional term $\partial_s \epsilon_{ab}^2/2$ should be present in the numerator of the aforementioned equation. So, in order to remain mathematically correct and to retain the possibility of making comparisons with the work of Müller et al. (2012), we decided to consider

<table><caption>**Table 1.** Definitions and list of abbreviations.</caption><thead><tr><th>Abbrev.</th><th>Definition + name</th><th>Symbol</th></tr></thead><tbody><tr><td>SG</td><td>Self-gravity</td><td>s = ||<b>r</b> - <b>r</b>'||</td></tr><tr><td></td><td>Mutual distance or separation</td><td>d<sub>b</sub> = s/H<sub>b</sub></td></tr><tr><td></td><td>Normalised separation for phase b</td><td></td></tr><tr><td>SL</td><td>Smoothing length</td><td></td></tr><tr><td>CSL</td><td>Constant smoothing length</td><td>ε<sub>ab</sub> = const.</td></tr><tr><td>SVSL</td><td>Space varying smoothing length</td><td>ε<sub>ab</sub>(d<sub>g</sub>)</td></tr><tr><td>SGFC</td><td>Self-gravity force correction</td><td>L<sub>sg</sub><sup>ab</sup> (Eq. (6))</td></tr><tr><td>SLFC</td><td>Smoothing length force correction</td><td>L<sub>ε</sub><sup>ab</sup> (Eq. (11))</td></tr><tr><td>PDFC</td><td>Planet-disc force correction</td><td>L<sup>P</sup> (Eq. (20))</td></tr><tr><td></td><td>a-to-b scale height ratio</td><td>η<sub>ab</sub></td></tr><tr><td></td><td>Gas-to-dust scale height</td><td>η = &lang;H<sub>g</sub>&rang;/&lang;H<sub>d</sub>&rang;</td></tr></tbody></table>

the potential that satisfies only Eq. (10) and not Eq. (9). This is a slight change in the SL paradigm that does not affect the approximation of the 2D SG terms and leads to meaningful 2D results. Based on this clarification, we can define an analogue of the SFGC that is compatible with the SL approach:

$$L_{\epsilon}^{ab}(d_g) = H_g(\mathbf{r}) s \|\nabla \Psi_{\epsilon}^{ab}\| \\ = \frac{\pi d_g^2}{[d_g^2 + (\epsilon_{ab}(d_g)/H_g(\mathbf{r}))^2]^{3/2}} \quad (11)$$

This normalised quantity, which we call the smoothing length force correction (SLFC) should fit the SGFC to correctly estimate SG in 2D simulations. This is only possible by wisely choosing a spatially dependent SL, namely, $\epsilon_{ab}(d_g)$, as depicted in the next section. The definition and abbreviation of the main quantities encountered in this paper are given in Table 1.

## 3. Mutual self-gravity interactions based on the smoothing-length approach

We aim here to provide a suitable SL that most closely approaches each of the possible exact gravitational interactions: gas-gas, dust-dust, and dust-gas. We started by retrieving and rectifying Müller et al. (2012) results for a disc only made of gas.

### 3.1. Contribution of the gas disc on a gas parcel

This case corresponds to $a=b=g$. From the general definition provided by Eq. (6), the gas SGFC is equal to<fnref n="2" />:

$$L_{\text{sg}}^{g}(d_{\text{g}}) = \frac{1}{2} d_{\text{g}}^{2} \iint_{u,v=-\infty}^{\infty} \frac{e^{-\frac{v^{2}}{2}} e^{-\frac{u^{2}}{2}}}{[d_{\text{g}}^{2} + (u - \eta_{\text{g}} v)^{2}]^{3/2}} du dv, \quad (12)$$

where $\eta_{\text{g}} = \langle H_{\text{g}}(\mathbf{r}') \rangle / \langle H_{\text{g}}(\mathbf{r}) \rangle = 1^3$. We present in Fig. 1 (top panel) the exact SGFC (blue solid line), evaluated numerically, and the equivalent SLFC (dashed lines) for three constant smoothing lengths (CSLs): $\epsilon_{\text{g}}(d_{\text{g}})/H_{\text{g}}(\mathbf{r}) = [0.0, 0.6, 1.2]$. In the

<footnote n="1">¹ Müller et al. (2012) defined as force correction a quantity that they defined as $I_{\text{sg}}$. We noticed that in some of their graphics, they mistook $sI_{\text{sg}}(s)$ for $I_{\text{sg}}(s)$. Therefore, for clarity and consistency with their work, we did not adopt this naming convention.</footnote>

<footnote n="2">² For sake of simplicity in the whole article we replaced all superscripts gg by g.</footnote>

<footnote n="3">³ We think that Müller et al. (2012) implicitly assumed this equality in their work.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Force corrections for the gas only. We compare the SGFC (exact value, in blue solid line) and the SLFC model for different smoothing lengths based on constant values for the SL (red, green and orange dashed lines) or on function $\epsilon_g(d_g)$ (black cross markers). For a CSL, $\epsilon_g = \text{const.}$, the SLFC and SGFC do not match at short separations: the error curve either tends either towards 0 or towards infinity. On the contrary, for a SVSL, $\epsilon_g = \epsilon_g(d_g)$, the SLFC matches SGFC with an accuracy better than 0.5%. We found that $L_0 = \sqrt{\pi}$ (Appendix C). Top: SGFC and SLFC. Bottom: relative error between SGFC and SLFC (in %): $100 \cdot |\text{SLFC}/\text{SGFC} - 1|$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Space-varying smoothing length for gas only. We compare the SVSL obtained numerically by Müller et al. (2012) and that proposed in this work which was obtained by analytical computations and curve fitting. Top: Müller et al. (2012). Bottom: this work.</figcaption></figure>

bottom panel, we show the respective errors in percentage. These curves support the results from Müller et al. (2012) for the case of a self-gravitating disc and also validates our approach. For long distances, the three aforementioned SLFC exhibit the same behaviour and the error becomes negligible. For instance, the CSL value, $\epsilon_g(d_g) = 1.2H_g$, offers the smallest error at long distances, namely, less than 2% for $d_g \ge 4$. On the contrary, at short distances, the SLFC either vanishes or diverges, which leads to significant errors since the exact SGFC converges towards $L_0 \approx 1.772$. Such behaviour is indeed intrinsic to the CSL formalism at short distances:

$$ \text{if } \epsilon_g(d_g) = 0: \qquad L_\epsilon(d_g) = \frac{\pi}{d_g} \qquad \rightarrow +\infty $$

$$ \text{if } \epsilon_g(d_g) = \text{const.} \neq 0: \qquad L_\epsilon(d_g) \underset{d_g \to 0}{\sim} \pi \left(\frac{H_g}{\epsilon_g}\right)^3 d_g \qquad \rightarrow 0. \qquad (13) $$

The above limits demonstrate that the CSL formalism either leads to 100% or infinite errors at short distances, when $\epsilon_g = \text{const.} \neq 0$ and $\epsilon_g = 0$, respectively. This divergence between the SGFC and SLFC was not highlighted by Müller et al. (2012). A possible reason for this is that the behaviour is suitable since it allows for numerical divergences to be avoided at $d_g = 0$ when $\epsilon_g(d_g) = \text{const.} \neq 0$. Even if Müller et al. (2012) did not shine a light on these divergences at short distances, they nevertheless carried out a numerical evaluation that demonstrates that a space-varying smoothing length (SVSL) is necessary to fit the exact SGFC (reproduced in the top panel of our Fig. 2). From this figure, we can infer that the SVSL should tend towards 0 so that the SLFC matches the SGFC. Yet this is in contradiction to Fig. 1 and Eq. (13), since the SLFC is supposed to diverge analytically for $\epsilon_g(d_g) = 0$ at short distances. In general, no matter the choice of CSL, the gas SG was underestimated by a factor 100 at short distances. This contradiction went unnoticed and we intended to resolve it, as described in the next section.

The analytic expression of the gas SL that provides a better fit to Fig. 2 (top panel) and thus permits the SLFC to fit the exact SGFC should fulfill the following constraints:

$$ \begin{array}{l} \lim_{d_g \to 0} \epsilon_g(d_g) = 0, \\ \lim_{d_g \to +\infty} \epsilon_g(d_g) = \sqrt{2} H_g(\mathbf{r}), \\ \mathcal{L} \epsilon_g(d_g = 0) = L_0. \end{array} \qquad (14) $$
<table><caption>**Table 2.** Fitting parameters for space varying smoothing lengths and $\delta$ models.</caption><tr><td>(α, n)</td><td>(0.04319874, 1.14791757)</td></tr><tr><td>(β, q)</td><td>(0.06427627, 1.14735482)</td></tr><tr><td>(γ, m)</td><td>(0.98362092, 0.75552227)</td></tr></table>

Thanks to analytical arguments completed by a curve fit (made explicit in Appendix A), we found a model for the SVSL which permits us to accurately approach the SGFC:

$$ \epsilon_g(d_g) = \sqrt{2} H_g(\mathbf{r}) \left[ 1 - \exp\left(-\frac{\epsilon_{g,0}}{\sqrt{2}} d_g^{2/3} - \alpha d_g^n \right) \right], \quad (15) $$

where $\epsilon_{0,g} = [\pi/L_0]^{1/3}$ and the numerical values of $(\alpha, n)$ are collected in Table 2. In particular, the power of $2/3$ and of $\epsilon_{0,g}$ allow us to reach a finite value of $L_0$ for the SLFC at $d_g = 0$. In Fig. 2 (bottom panel), we show the SVSL analytic model. In Fig. 1, we show the respective SLFC and error with respect to the exact SGFC (black cross markers). Within this correction, the SLFC and SGFC overlap and the error is reduced to less than 0.5% across the whole distance range. As a matter of comparison, this correction allows us to decrease the error up to factors 200 and 40 at short and long distances, respectively, compared to the CSL where $\epsilon_g(d_g) = 1.2H_g(\mathbf{r})$.

## 3.2. Contribution of the dust layer on a dust parcel

In this Sect., we take $a = b = d$. From the general definition provided by Eq. (6), the dust SGFC is<fnref n="4" />:

$$ \begin{aligned} L_{\text{sg}}^{d}(d_g, d_d) &= \frac{1}{2} \frac{d_d^3}{d_g} \iint_{u,v=-\infty}^{\infty} \frac{e^{-\frac{u^2}{2}} e^{-\frac{v^2}{2}}}{[d_d^2 + (u - \eta_d v)^2]^{3/2}} du dv \\ &= \frac{1}{2} \eta (\eta d_g)^2 \iint_{u,v=-\infty}^{\infty} \frac{e^{-\frac{u^2}{2}} e^{-\frac{v^2}{2}}}{[(\eta d_g)^2 + (u - v)^2]^{3/2}} du dv \\ &= \eta L_{\text{sg}}^{g}(\eta d_g), \end{aligned} \quad (16) $$

where $\eta_d = \langle H_d(\mathbf{r}') \rangle / \langle H_d(\mathbf{r}) \rangle = 1$ and $\eta = H_g(\mathbf{r})/H_d(\mathbf{r})$ is the gas-to-dust scale height ratio that, in general, is greater than unity. Again, for simplicity, we prefer to use $\eta = \langle H_g(\mathbf{r}) \rangle / \langle H_d(\mathbf{r}) \rangle$. The above relation naturally applies to the dust SLFC, but it is convenient to write it one of two different ways:

$$ L_{\epsilon}^{d}(d_g) = \eta L_{\epsilon}^{g}(\eta d_g) \quad \text{or} \quad \epsilon_d(d_g, \eta) = \frac{\epsilon_g(\eta d_g)}{\eta}. \quad (17) $$

Both equations are equivalent, but the former provides a physical insight while the latter allows for a simple implementation in hydrodynamical codes. In Fig. 3, the top and bottom panels show the SLFC and the error with respect to the SGFC for different gas-to-dust scale height ratios, respectively. The dust SLFC and SGFC curves overlap over the whole $\eta$ and distance ranges, so we did not plot these two quantities in the same figure to avoid duplication. The error is again lower than 0.5% in the whole distance range, as expected from the unique gas disc case. We want to highlight that at long separations the dust SLFC matches the

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Self-gravity force corrections for dust only. We show the dependence of the dust force corrections for different values of the gas-to-dust scale-height ratio ($\eta$). For the whole range of $\eta$ values, the accuracy is better than 0.5%. The dust SGFC is proportional to $\eta$ at short distances ($d_g \le 1.5/\eta$). For large values of $\eta$, this could favour gravitational clumping. Top: SGFC. Bottom: relative error between SGFC and SLFC (in %): $100 \cdot |\text{SLFC}/\text{SGFC} - 1|$.</figcaption></figure>

gas SLFC, $L_{\text{sg}}^d(d_g) = L_{\text{sg}}^g(d_g)$, while at short distances, the dust SG is $\eta$ times stronger than gas SG and we get $L_{\text{sg}}^d(d_g = 0) = \eta L_0$. This latter equality is important because (1) if the gas SLFC is used instead of the dust SLFC and (2) if a CSL is used instead of a SVSL – then the dust SG is underestimated by a factor of ~100 $\eta$ at short distances. These aspects are of primary interest in an astrophysical context since it is typical to see $\eta \gtrsim 10$, which could have important implications for planetary formation theories (as discussed in Sect. 5.2).

## 3.3. Crossed contributions

In this section, we take $a = d$ and $b = g$. We highlight that the SGFC is commutative with respect to phases a and b. This can be demonstrated thanks to Newton’s third law or based on analytical arguments (see Appendix B). Therefore, in the following, we will adopt the notation $L_{\text{sg}}^{dg}$ to refer to either the dust-gas or gas-dust

<footnote n="4">⁴ For sake of simplicity in the whole article we replaced all superscripts dd by d.</footnote>
SGFC. In the context of this clarification and Eq. (6), the dust-gas SGFC is:

$$L_{\text{sg}}^{dg}(d_{\text{g}}, \eta) = \frac{1}{2} d_{\text{g}}^2 \iint_{u, v=-\infty}^{\infty} \frac{e^{-\frac{u^2}{2}} e^{-\frac{v^2}{2}}}{[d_{\text{g}}^2 + (u - v/\eta)^2]^{3/2}} du dv. \quad (18)$$

Motivated by the SL formalism, we again aim to approach the above double integral by a dust-gas SLFC, but to our knowledge, this has never been done before. In the next sections, we explore, for first time, what constraints should be satisfied in order to construct a consistent SL, which allows us to accurately approach the exact dust-gas SGFC.

### 3.3.1. Limiting cases for weak and strong layering

When layering is minimal ($\eta = 1$), the gas and dust are fully mixed and we immediately retrieve the case of a pure gas disc:

$$\lim_{\eta \to 1} L_{\text{sg}}^{dg}(d_{\text{g}}, \eta) = L_{\text{sg}}^{g}(d_{\text{g}}). \quad (19)$$

On the other hand, when layering is strong ($\eta \gg 1$), the dust layer is infinitely thin and we get:

$$\lim_{\eta \to +\infty} L_{\text{sg}}^{dg}(d_{\text{g}}, \eta) = \sqrt{2\pi} I_p(d_{\text{g}}^2/4) = L^p(d_{\text{g}}), \quad (20)$$

where $I_p(x) = xe^x[K_1(x) - K_0(x)]$ and $K_1$ and $K_0$ are the modified Bessel functions of the second kind. Interestingly, $I_p(d_{\text{g}}^2/4)$ is a force function in the limiting case of a planet interacting with a gas disc, as exposed in Müller et al. (2012, Sect. 4). Accordingly, we recall the quantity $L^p$ the planet–disc force correction (PDFC). The PDFC can be approached in the SL formalism and using the same analytical approach conducted in Sect. 3.1 and the results in Müller et al. (2012, Sect. 4.1), and further described in Appendix A. Thus, we find the planet–disc SVSL:

$$\epsilon_p(d_g) = H_g(r) \left[ 1 - \exp\left(-\epsilon_{0,p} d_g^{2/3} - \beta d_g^q\right) \right], \quad (21)$$

where $\epsilon_{0,p} = (\pi/2)^{1/6}$ and the power $2/3$ were obtained by analytical means, while $(\beta, q)$ by curve-fitting. The latter parameters are gathered in Table 2. We note that, contrary to our initial guess, for strong dust layering, the dust-gas SLFC tends towards the PDFC – rather than towards the dust SLFC. In summary, the dust-gas SL equals the gas SL or planet–disc SL when the gas-to-dust height, $\eta$, tends towards unity or infinity, respectively.

### 3.3.2. Behaviours at short and long separation

To study dust-gas SGFC at short distances, it is convenient to take the following function:

$$\delta(\eta) = \lim_{d_{\text{g}} \to 0} L_{\text{sg}}^{dg}(d_{\text{g}}, \eta)/L_0. \quad (22)$$

We did not succeed to express this function in terms of standard mathematical functions, but we estimated it numerically in top panel of Fig. 4 (blue solid line). As expected, $\delta$ is equal to 1 for $\eta = 1$ ($L_{\text{sg}}^{dg}$ matches $L_{\text{g}}^{g}$) while for an infinitely thin layer of dust this quantity tends to $\sqrt{2}$. Accounting for these boundary conditions and a noticeable curve shape, we get:

$$\delta(\eta) = \sqrt{2} + (1 - \sqrt{2}) \exp[-\gamma(\eta - 1)^m], \quad (23)$$

<figure><img src="image_6.png" /><figcaption>Fig. 4. Model functions for defining the dust-gas SVSL. Top: $\delta$ function: model and numerical estimation with respect to the gas-to-dust height ratio ($\eta$). Bottom: $\lambda$ model with respect to $\eta$.</figcaption></figure>

noting that the couple $(\gamma, m)$ can be found in Table 2. This analytic model is also shown in top panel of Fig. 4 (orange dashed line). In contrast, for long separations, the dust-gas SGFC should satisfy:

$$\lim_{d_{\text{g}} \to +\infty} L_{\text{sg}}^{dg}(d_{\text{g}}, \eta) = \frac{\pi}{d_{\text{g}}}. \quad (24)$$

This implies that at long distances, the dust-gas SVSL is expected to be negligible compared to the square of the distance, $\epsilon_{dg}(d_g)/H_g(r) = o(d_g^2)$, which is undoubtedly satisfied if the SVSL is constant at long distances.

### 3.3.3. Dust-gas space-varying smoothing-length construction

Now, we have on hand all the necessary information needed to build a consistent dust-gas SVSL. We made the choice to look for the dust-gas SVSL in the form of a linear combination of the previous asymptotic cases:

$$\epsilon_{dg}(d_g, \eta) = \lambda(\eta) \epsilon_g(d_g) + (1 - \lambda(\eta)) \epsilon_p(d_g), \quad (25)$$

where $\lambda(\eta) \in [0, 1]$. For matching Sect. 3.3.1 constraints, we chose:

$$\lim_{\eta \to 1} \lambda(\eta) = 1 \quad \Rightarrow \quad \lim_{\eta \to 1} \epsilon_{dg}(d_g, \eta) = \epsilon_g(d_g), \\ \lim_{\eta \to +\infty} \lambda(\eta) = 0 \quad \Rightarrow \quad \lim_{\eta \to +\infty} \epsilon_{dg}(d_g, \eta) = \epsilon_p(d_g). \quad (26)$$
Additionally, from the results given in Sect. 3.3.2 and the Taylor expansion described in Appendix C, we get:

$$ \lambda(\eta) = \frac{\epsilon_{0,g} (1/\delta(\eta))^{1/3} - \epsilon_{0,p}}{\epsilon_{0,g} - \epsilon_{0,p}}. \qquad (27) $$

This $\lambda$ function is plotted in the bottom panel of Fig. 4, where we checked that the above-mentioned boundary conditions are met. We also plotted, in the top panel of Fig. 5, the dust-gas SVSL, defined in Eq. (25), for different gas-to-dust scale height ratios. For $\eta = 1$, the dust-gas SVSL matches the gas SVSL while for modest dust layering (i.e. $\eta \ge 5$), the SVSL tends rapidly to the one of the planet-disc interaction case. This allows us to use the approximation: $\epsilon_{dg}(d_g) \approx \epsilon_p(d_g)$ for $\eta \ge 5$. Finally, we found that our results are mathematically consistent provided that $L_0 = \sqrt{\pi}$. Comparing this theoretical prediction with the value of $L_0$ obtained numerically, we find that both results match with an accuracy of up to six decimal points.

### 3.3.4. Summary for the dust-gas space-varying smoothing length

The middle and bottom panels of Fig. 5 show the exact dust-gas SGFC and the error between both estimations for different $\eta$ values, respectively. The dust-gas SLFC and SGFC curves overlap for the whole $\eta$ and distance range, so we did not plot these two quantities in the same figure to avoid duplication. This overlap is also reflected in the error, which this time depends on the dust layering: for $\eta \in [1, 5]$, the error is lower than 2% and for $\eta \in [5, 100]$, the error is lower than 0.5% for the whole distance range. Compared to both cases studied in Sects. 3.1 and 3.2, the error is slightly higher but it remains acceptable. We want to stress that the dust-gas SLFC rapidly tends (with respect to $\eta$) toward the planet-disc SLFC, which makes the following simplification possible:

$$ L_{sg}^{dg}(d_g, \eta) = L^p(d_g) \quad \text{or} \quad \epsilon_{dg}(d_g, \eta) = \epsilon_p(d_g) \quad \text{if} \quad \eta \ge 5. \qquad (28) $$

This approximation could simplify the numerical treatment.

## 4. Numerical treatment

The main goal of our study is to implement an accurate SG computation method for multi-fluids in 2D numerical codes. This could be beneficial for the 2D version ($r, \theta$) of hydrodynamical codes such as RoSSBi3D (Rendon Restrepo et al. 2022), FARGO (Masset 2000) or Athena (Stone et al. 2008). We start by treating the singularity for vanishing separations, responsible of numerical divergences, then we explain explicitly under which conditions results of Sect. 3 could be used for estimating SG thanks to FFT methods. Finally, we quantify the computational costs for 2D N-fluid simulations with SG.

### 4.1. Removing numerical divergences

From Eq. (11), it is obvious that $||\nabla\Psi_e^{ab}|| \propto d_g^{-1}$ is valid at short separations, which makes it necessary to have a product with a tapering function, so as to cancel SG for vanishing separations. The tapering function, $f$, must be equivalent to $d_g^p$ at short distances, with $p > 1$. On the other hand, it should be equal to unity at large distances. In order to avoid losing the accuracy achieved

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Force corrections and space-varying smoothing-length for dust-gas. We show the dependence of all quantities for different values of the gas-to-dust scale-height ratio ($\eta$). For thin ($\eta \gg 1$) and thick ($\eta = 1$) dust layers, the SVSL tends towards the one of the gas and of the planet, respectively. At short distances the dust-gas SLFC is restricted to the interval $[L_0, \sqrt{2}L_0]$. The accuracy of the SL method for the gravitational interaction of the gas disc with the embedded dust is better than 2% for the whole separation range. Top: SVSL. Middle: SGFC. Bottom: relative error between SGFC and SLFC (in %): $100 \cdot |\text{SLFC}/\text{SGFC} - 1|$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Normalised Plummer potential gradient with and without a tapering function. In the SVSL method, the introduction of a tapering function (Eq. (29)) avoids the divergence of the Plummer potential. At the same time, the low tapering length, $d_s(= 3/127 \approx 0.024)$ ensures that most of the SG contribution is not lost over short distances ($d_g \le 1.5$). Under the CSL approach (red dotted), the gradient is cancelled out at the singularity, although most of the short distance contribution to the SG is lost.</figcaption></figure>

in the steps described in Sect. 3, the tapering length should be approximately as long as the finest numerical resolution in the grid: $r_s \sim 3 \min(\Delta r, r\Delta\theta)^{5,6}$, where $\Delta r$ and $\Delta\theta$ are the resolution in the radial and azimuthal directions, respectively. To satisfy the above constraints, we arbitrarily chose $p = 2$ and the tapering function as follows:

$$f(d_g) = 1 - \exp\left[-\frac{1}{2}\left(\frac{d_g}{d_s}\right)^2\right] \underset{d_g \to 0}{\sim} \frac{1}{2}\left(\frac{d_g}{d_s}\right)^2, \quad (29)$$

where $d_s = r_s/H_g$ is the normalised tapering length. In Fig. 6, we show the normalised gradient of the gas Plummer potential with (orange dashed line) and without (blue solid line) the tapering function correction for $d_s = 3/127$. We also plotted the same quantity under a CSL assumption (red dotted line). Such a tapering length choice was motivated by the high resolution achieved by Rendon Restrepo & Barge (2022), namely: 146 cells/H_g and 127 cells/H_g in the radial and azimuthal directions. As expected, both Plummer potential gradients exhibit the same behaviour for $d_g \gtrsim d_s$ but only the corrected one converges towards 0 for vanishing distances, which should avoid numerical divergences. Regarding the CSL assumption, we clearly see that SG is underestimated for $d_g \le 1.5$, but there is no need to resort to a tapering function since the potential gradient is cancelled out analytically at the singularity. This may be the reason why Müller et al. (2012) did not mention that there were 100% errors with the CSL method at short separations.

## 4.2. Numerical calculation with FFT methods

An important question regarding our results is whether FFT methods can still be used when we resort to a SVSL. Indeed, the classic SG computation in 2D was based on a CSL. We find it instructive to recall here how SG is computed when including a SVSL and a tapering function. Guided by

Baruteau & Masset (2008) and Surville (2013, in french) we get the radial component of the SG mass force (divided by $\Sigma_b)^7$:

$$\begin{aligned} & -\frac{\nabla \Phi_e^{ab} \cdot \mathbf{e}_r}{\pi} \\ &= \iint_{\text{disc}} f(d_g) \frac{\Sigma_a(\mathbf{r}')(\mathbf{r} - \mathbf{r}')}{[s^2 + \epsilon_{ab}(d_g, \mathbf{r})^2]^{\frac{3}{2}}} \mathbf{e}_s \cdot \mathbf{e}_r \, \mathrm{d}\mathbf{r}' \, \mathrm{d}\theta' \\ &= \iint_{\text{disc}} f(d_g) \frac{\Sigma_a(\mathbf{r}') \left(\left(\frac{r}{r'}\right)^2 - \cos(\theta - \theta')\right) \frac{\mathrm{d}\mathbf{r}'}{r'} \, \mathrm{d}\theta'}{\left[1 + \left(\frac{r}{r'}\right)^2 - 2\left(\frac{r}{r'}\right) \cos(\theta - \theta') + \delta\epsilon_{ab}^2\right]^{\frac{3}{2}}}, \end{aligned} \quad (30)$$

where:

$$\delta\epsilon_{ab} = \frac{\epsilon_{ab}(d_g, \mathbf{r})}{r'} = \frac{r}{r'} h_g(r) \frac{\epsilon_{ab}(d_g, \mathbf{r})}{H_g(r)}, \quad (31)$$

$$d_g = \frac{r'}{r} \frac{1}{h_g(\mathbf{r})} \sqrt{1 + \left(\frac{r}{r'}\right)^2 - 2\left(\frac{r}{r'}\right)\cos(\theta - \theta')}, \quad (32)$$

and $h_g(\mathbf{r}) = H_g(\mathbf{r})/r$ is the gas disc aspect ratio. We aim to write the integral, defined by Eq. (30), as a convolution product that is only possible provided that the tapering function, $f$, and the ratio $\delta\epsilon_{ab}$ could be written as $r/r'$ and $\theta - \theta'$ functions. This constraint is satisfied if the disc aspect ratio is a spatial constant, $h_a(\mathbf{r}) = h_{0,a}$, in the whole simulation box. Within this condition, the radial component of the SG force is:

$$-\frac{\nabla \Phi_e^{ab} \cdot \mathbf{e}_r}{\pi} = \iint_{\text{disc}} \Sigma_a(X', \theta') \mathcal{G}_r^{ab}(X - X', \theta - \theta') \mathrm{d}X' \mathrm{d}\theta', \quad (33)$$

where:

$$\mathcal{G}_r^{ab}(X - X', \theta - \theta') = \frac{f(X - X', \theta - \theta') [e^{(X - X')} - \cos(\theta - \theta')]}{[1 + e^{2(X - X')} - 2e^{(X - X')} \cos(\theta - \theta') + \delta\epsilon_{ab}(X - X', \theta - \theta')^2]^{\frac{3}{2}}},$$

is the modified radial Green's function where we performed the variable substitution: $r = e^X$ and $dr = e^X \, dX$. For the azimuthal component, the calculation is similar and we obtain:

$$-\frac{\nabla \Phi_e^{ab} \cdot \mathbf{e}_\theta}{\pi} = \iint_{\text{disc}} \Sigma_a(X', \theta') \mathcal{G}_\theta^{ab}(X - X', \theta - \theta') \mathrm{d}X' \mathrm{d}\theta', \quad (34)$$

where:

$$\mathcal{G}_\theta^{ab}(X - X', \theta - \theta') = \frac{f(X - X', \theta - \theta') \sin(\theta - \theta')}{[1 + e^{2(X - X')} - 2e^{(X - X')} \cos(\theta - \theta') + \delta\epsilon_{ab}(X - X', \theta - \theta')^2]^{\frac{3}{2}}}$$

is the modified azimuthal Green's function. The rewriting of both integrals, expressed via Eqs. (33) and (34), permits us to finally write out the equations in terms of a set of Fourier transforms:

$$\begin{aligned} f_{2\mathrm{D},\epsilon}^{a \to b}(\mathbf{r}) \cdot \mathbf{e}_r &= \Sigma_b(\mathbf{r}) \mathcal{F}^{-1}[\mathcal{F}(\Sigma_a) * \mathcal{F}(\mathcal{G}_r^{ab})] \\ f_{2\mathrm{D},\epsilon}^{a \to b}(\mathbf{r}) \cdot \mathbf{e}_\theta &= \Sigma_b(\mathbf{r}) \mathcal{F}^{-1}[\mathcal{F}(\Sigma_a) * \mathcal{F}(\mathcal{G}_\theta^{ab})], \end{aligned} \quad (35)$$

<footnote n="7">⁷ We took $G = 1$.</footnote>

<footnote n="5">⁵ In cylindrical coordinates.</footnote>

<footnote n="6">⁶ The factor 3 was chosen for ensuring a safety margin.</footnote>
where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier transform operator and its inverse. The symbol $*$ denotes the convolution operator. In practice, such quantities are computed numerically thanks to fast Fourier modules which rapidly accelerate numerical computations. For instance, in the 2D version of RoSSBi3D (Rendon Restrepo et al. 2022), this computation is made possible thanks to the FFTW3 library (Frigo & Johnson 2005). It is important to highlight that the use of Fourier transforms is possible if: (1) a logarithmic mesh is used in the radial direction; this condition stems from the radial variable substitution that enabled us to obtain the formulation with the modified Green's functions. Furthermore, (2) the function $\Sigma_a$ must be periodic in radial and azimuthal directions in order to be able to use Fourier transforms. Periodicity is obviously satisfied in the azimuthal direction, but in the radial direction this is not necessarily the case. Such a periodicity is artificially ensured thanks to a zero-padding approach: the radial domain is doubled and the function $\Sigma_a$ is set to 0 in half of the domain. See Surville (Fig. III.10 2013, in french) for an example. Finally, it is necessary for (3) the disc aspect ratio of gas and dust to be constant with respect to $r$ (but it could be a time-varying function). If such condition is not satisfied, the formulation as a convolution product is not possible.

We take this opportunity to clarify that the aforementioned assumption of a constant disc aspect ratio used in this section only permits us to resort to FFT methods for accelerating the numerical calculations. This condition is not new since it was already implicit in the classical calculation of the SG by FFT. For general scale heights, a direct summation in the radial direction and Fourier transforms in the azimuth direction, ensured by the periodicity, is a straightforward solution.

### 4.3. Computational costs

For SG simulations, the computational endeavour is a non-negligible aspect, particularly for high-resolution simulations. For concision and without any loss of generality, we only treat the computation of SG in the radial direction and we do not account the Fourier transforms of the modified Green's functions, since they are computed only once. For the case of a bi-fluid simulation, the next Fourier transforms should be performed at each time step:

Gas into gas:

$$ \mathcal{F}(\Sigma_g) \text{ and } \mathcal{F}^{-1}[\mathcal{F}(\Sigma_g) * \mathcal{F}(G_r^d)], \quad (36) $$

Dust into dust:

$$ \mathcal{F}(\Sigma_d) \text{ and } \mathcal{F}^{-1}[\mathcal{F}(\Sigma_d) * \mathcal{F}(G_r^d)], \quad (37) $$

Dust into gas:

$$ \mathcal{F}(\Sigma_d) \text{ and } \mathcal{F}^{-1}[\mathcal{F}(\Sigma_d) * \mathcal{F}(G_r^{dg})], \quad (38) $$

Gas into dust:

$$ \mathcal{F}(\Sigma_g) \text{ and } \mathcal{F}^{-1}[\mathcal{F}(\Sigma_g) * \mathcal{F}(G_r^{dg})]. \quad (39) $$

From the recapitulation above, we observe that the Fourier transforms of gas and dust densities are duplicated, which reduces the total number of Fourier operations to six during a numerical treatment. Therefore, the numerical endeavour for computing SG for a bi-fluid is three times larger with respect to the case of a single fluid. In general, for $N$ fluids with different scale heights

and interacting through SG, the number of required Fourier transforms is:

$$ N^2 + N = \underbrace{(2N^2)}_{\text{All } N\text{-tuples combinations}} - \underbrace{(N^2 - N)}_{\text{Duplicates}}. \quad (40) $$

Compared to standard self-gravitating simulations on a unique gaseous phase, the computation time spent in the self-gravity module is crudely multiplied by a factor $\sim N^2/2$ for large $N$. Of course, this amount can be decreased assuming that some of the involved fluids have the same scale heights.

## 5. Discussion

In this section, we treat the possible impact of our findings regarding planet migration and the early stage of planetary formation. Then we identify the limitations of the SVSL approach due to our initial assumptions and to the specificity of the studied problem. Finally, we propose possible ways to improve and test our model.

### 5.1. Consequences for planet–disc interaction

It is widely known that planets can migrate due to tidal interactions with the gas disc. This is the case of type I migration, whereby the planet can exchange angular momentum with Lindblad and co-rotation resonances (Baruteau & Masset 2013, for a review). In 2D numerical simulations, the value $\epsilon_p/H_g(r) = 0.3-0.6$ is often used to match the total torque exerted on a planet in 3D simulations (Masset 2002; Tanaka et al. 2002).

Similarly to our results of Sect. 3.1, the planet–disc SLFC is not well captured by a CSL for separations inferior to $\sim 1.5 H_g$. Although we do not question the agreement between the 2D and 3D simulations of the planet–disc interaction, we do believe that our SVSL may be better suited than an adjustment factor. Therefore, it might be constructive to verify whether (1) the results of 3D simulations can be retrieved using our SVSL and (2) to check whether the widely used $\epsilon_p/H_g(r) = 0.3-0.6$ factor could be retrieved from an analytical basis stemming from our SVSL. We also think that our SVSL approach could be helpful for addressing 2D simulations of low-mass planets embedded in a self-gravitating disc. This would, however, require some improvements: the vertical layering due to SG must be accounted for, as discussed in Sect. 5.3.1, and the vertical stratification of the gas disc due to the planet gravitation must be assessed, as noted by Müller et al. (2012). Both aspects raised in this paragraph require a more detailed work, which is out of the scope of present paper.

### 5.2. Consequences for planet formation theories

One attractive planetary formation scenario is based on the persistence of gaseous vortices in PPDs. Its main interest is in the strong capture efficiency of the vortices and their ability to confine large concentrations of dust grains that could collapse to form planetesimals or a planetary core (Barge & Sommeria 1995). Even if observational findings seem encouraging in this regard (Varga et al. 2021), numerical simulations have not yet concluded that vortices could form objects that are bound by gravity. The results of the present paper offer the possibility to carry out new numerical simulations that correctly account for SG in the gas and dust components of PPDs. In particular, we expect that these new simulations could answer the questions
raised by the vortices. Indeed, as demonstrated in Sect. 3.2, dust SG could be underestimated by a factor of ~2000 for $\eta = 20$ at short separations. At the same time, the estimation of the dust-gas SLFC could allow us to quantify (with an acceptable level of accuracy) a possible gaseous envelope capture by dust clumps. From theoretical analyses, it has been found that SG inhibit vortices formation by a Rossby wave instability (Lovelace & Hohlfeld 2013) and this was later confirmed via numerical simulations (Baruteau & Zhu 2016; Regály & Vorobyov 2017; Tarczay-Nehéz et al. 2022). In the latest numerical work to date, Rendon Restrepo & Barge (2022) also found that vortices cannot survive in massive PPD and they provided a stability criterion that vortices should satisfy in order to resist SG destabilising effects. In light of our findings, the results of previous simulations should be checked anew to better understand the extent to which the SVSL affects theoretical predictions on the survival of vortices in self-gravitating PPDs.

### 5.3. Limitations, improvements, and tests

The limitations of our model are inherent to the initial assumptions we have made about the vertical structure of the disc in relation to gas and dust stratification.

#### 5.3.1. Stratification and disc evolution

The SGFC studied in this paper is based on the vertical integration of Eq. (6) in the particular case of a vertically isothermal disc. However the vertical structure could be affected by different mechanisms, which implies that for any vertical stratification different from the Gaussian stratification, the entire work performed in this paper should be repeated and adapted. For instance, this is the case when including the disc vertical SG for gas (but neglecting the vertical component of the central object gravity), which modifies the Gaussian distribution into: $\cosh\left(\frac{z}{QH_g}\right)^{-2}$, where Q is Toomre's parameter (Lodato 2007). A similar layering should occur for the dusty disc and we also expect that the crossed gravitational interaction between both phases could impact their respective vertical stratification. Indeed, gas SG could decrease dust scale height by a factor of ~2 (Baehr & Zhu 2021) and we expect that a strong dust layering will also modify the gas vertical structure in correlation with the dust-to-gas density ratio.

The global vertical structure of PPDs evolves over time due to accretion heating (Schobert et al. 2019) and stellar irradiation (Wu & Lithwick 2021), amongst other effects. Instabilities could also generate time variable structures which could affect (locally) the vertical stratification of the flow. This was reported, for instance, in 3D vortices simulations (Meheut et al. 2010) for rings and gaps generated by poloidal magnetic winds (Suriano et al. 2017) and for spiral density waves (Riols & Latter 2018). The method described in this paper is not limited to a steady vertical stratification and it is compatible with the global and local time evolution of the vertical structure, provided that the time dependence between gas and dust scale heights is known. This could be done, for instance, with the 2D1D strategy adopted by Crida et al. (2009). However, in the particular case of spatially constant aspect ratios, $h_a(r, t) = h_a(t)$, the computational benefits from the FFT method could be lost. Indeed, this would require us to perform the Fourier transforms of Green's functions at each time step instead of a unique computation for a steady vertical stratification. For a simulation with only gas, this would result in three Fourier transforms at each time step (instead of

two). We believe that for other SG computational methods, the computational costs would be unaffected.

#### 5.3.2. Layering of the dust particles

In this paper, the dust component of PPDs is considered as a pressure-less fluid that is sufficiently mixed by the turbulent motions of the gas disc to be maintained in an equatorial sub-layer. This assumption requires that the dust particles and gas aerodynamic coupling are governed by a Stokes number less than unity and that dust is adequately diluted in the gas to avoid frequent mutual interactions (if the dust-to-gas mass-ratio $\le 1$) (Garaud et al. 2004). We also assumed that small-scale turbulence is maintained in the gas disc by a mechanism that we ended up disregarding. It is interesting to note that outside the bi-fluid pressure-less assumption, turbulent stirring may be replaced by collisional and/or gravitational stirring with the formation of a sub-layer of solid particles, whose scale height is different from the one deduced from turbulent stirring in the introduction. In such cases, the necessary smoothing lengths will be different – but the SVSL approach should remain unchanged.

#### 5.3.3. Additional test

A relevant test for our results consists of a comparison with the vertically averaged SG obtained thanks to a 3D simulation. In order to satisfy our assumptions, the 3D disc should be vertically isothermal and the vertical SG should be discarded. In addition, for the bi-fluid version, the dust layer should be as smooth as possible. Such a simulation is beyond the scope of current paper, but still remains interesting lead that we want to explore in the near future.

## 6. Conclusions

In this work, we reveal the contradictions and shortcomings of the CSL method commonly used to compute the contribution of SG in 2D numerical simulations. In particular, we found that from short to intermediate separations ($d_g \le 1.5$), the SG force is analytically underestimated with an error that reaches 100% at the singularity. We corrected these inconsistencies replacing the CSL, $\epsilon_g = \text{const.}$, by a space dependent function, $\epsilon_g = \epsilon_g(d_g)$ (SVSL). We found that for a gas disc, the SVSL dependence that better fits the exact SGFC is:

$$ \epsilon_g(d_g) = \sqrt{2} H_g(\mathbf{r}) \left[ 1 - \exp\left(-\frac{\epsilon_{0,g}}{\sqrt{2}} d_g^{2/3} - \alpha d_g^n \right) \right], \quad (41) $$

where $\epsilon_{0,g} = \pi^{1/6}$ and $(\alpha, n)$ can be found in Table 2. This SVSL approach can be extended to the dust disc using the following dust SL:

$$ \epsilon_d(d_g, \eta) = \frac{\sqrt{2} H_g(\mathbf{r})}{\eta} \left[ 1 - \exp\left(-\frac{\epsilon_{0,g}}{\sqrt{2}} (\eta d_g)^{2/3} - \alpha (\eta d_g)^n \right) \right], \quad (42) $$

where $\eta = \langle H_g \rangle / \langle H_d \rangle$ is the gas-to-dust height ratio. As a side result, we also found the planet-disc SVSL:

$$ \epsilon_p(d_g) = H_g(\mathbf{r}) \left[ 1 - \exp\left(-\epsilon_{0,p} d_g^{2/3} - \beta d_g^q \right) \right], \quad (43) $$

where $\epsilon_{0,p} = (\pi/2)^{1/6}$ and $(\beta, q)$ can be found in Table 2. The crossed gravitational interaction of the gas with the embedded dust can also be evaluated through the SL method.
We constructed this dust-gas SL from a linear combination of aforementioned planet-disc and gas SVSL:

$$ \epsilon_{dg}(d_g, \eta) = \lambda(\eta) \epsilon_g(d_g) + (1 - \lambda(\eta)) \epsilon_p(d_g), \quad (44) $$

where the analytical expression of $\lambda$ is given in Eq. (27). All these new expressions for the SVSL are valid for any scale height if the stratification is Gaussian and remain compatible with the common FFT method for evaluating SG in 2D hydrodynamical simulations, provided that the gas and dust disc aspect ratios are constant. Finally, the use of a tapering function is required to avoid numerical divergences.

The proposed SVSL correction decreases the error up to factors 200 with respect to the latest CSL prescription proposed by Müller et al. (2012). In particular, our SVSL allows us to match the SGFC with a high accuracy, even at the singularity ($d_g \to 0$). Regarding the dust SLFC, we found that it is proportional to $\eta$ at short separations. This result combined with the improvement brought by the SVSL method demonstrates that dust SG was generally underestimated by a factor $\sim 100 \eta$ at short separations (compared to the CSL only based into the gas SL).

Our planet-disc SVSL could affect the torque exerted by the self-gravitating disc on a planet, which suggests that type I planet migration could be impacted. We also think that the improvements we have made in the computation of the SG terms will be decisive in future 2D simulations of PPDs inhabited by a large-scale vortex. Hence, we forecast a much better description of the evolution of the dust-gas mixture in the core of the vortices with our model than with standard ones; we also expect significant consequences in the simulations of planetesimal construction. Indeed, on the one hand, dust SLFC could favour the gravitational binding of the dust clumps trapped in the vortex and, on the other hand, dust-gas SLFC allows us to follow the coupled evolution of dust and gas in the various clumps. We speculate that gas could be dragged with dust during the collapse, before being trapped in a gaseous envelope around a dusty core.

Acknowledgements. We thank the referee for her/his helpful comments that enriched the discussion. We would like to warmly acknowledge Stéphane Le Dizès for fruitfull discussions during the preparation of the paper

and financial support during the PhD thesis. S.R.R. thanks Clément Baruteau for useful discussions and Andrej Hermann for proofreading the article. Co-funded by the European Union (ERC, Epoch-of-Taurus, 101043302). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

Baehr, H., & Zhu, Z. 2021, ApJ, 909, 136

Barge, P., & Sommeria, J. 1995, A&A, 295, L1

Baruteau, C., & Masset, F. 2008, ApJ, 678, 483

Baruteau, C., & Masset, F. 2013, in Lecture Notes in Physics, eds. J. Souchay, S. Mathis, & T. Tokieda (Berlin: Springer Verlag), 861, 201

Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927

Crida, A., Masset, F., & Morbidelli, A. 2009, ApJ, 705, L148

Dong, R., Rafikov, R. R., Stone, J. M., & Petrovich, C. 2011, ApJ, 741, 56

Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237

Frigo, M., & Johnson, S. G. 2005, Proc. IEEE, 93, 216

Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292

Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051

Huré, J. M., & Pierens, A. 2009, A&A, 507, 573

Johansen, A., & Youdin, A. 2007, ApJ, 662, 627

Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52

Lodato, G. 2007, Nuovo Cimento Rivista Serie, 30, 293

Lovelace, R. V. E., & Hohlfeld, R. G. 2013, MNRAS, 429, 529

Masset, F. 2000, A&AS, 141, 165

Masset, F. S. 2002, A&A, 387, 605

Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31

Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123

Regály, Z. & Vorobyov, E. 2017, MNRAS, 471, 2204

Rendon Restrepo, S., & Barge, P. 2022, A&A, 666, A92

Rendon Restrepo, S., Barge, P., & Vavrik, R. 2022, ArXiv e-prints [arXiv:2207.04252]

Riols, A., & Latter, H. 2018, MNRAS, 476, 5115

Schobert, B. N., Peeters, A. G., & Rath, F. 2019, ApJ, 881, 56

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137

Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2017, MNRAS, 468, 3850

Surville, C. 2013, Ph.D. Thesis, Aix-Marseille Université, France

Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257

Tarczay-Nehéz, D., Rozgonyi, K., & Regály, Z. 2022, MNRAS, 511, 6055

Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56

Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178

Wu, Y., & Lithwick, Y. 2021, ApJ, 923, 123
## Appendix A: Gas and planet space varying smoothing length derivation

<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Gas and planet-disc force corrections. For each force correction, we plotted the exact SGFC evaluated numerically (blue solid line) and the SLFC approximation accounting only for the analytical correction (orange dashed line) or accounting for an additional numerical fit (black cross markers). Top: Gas force correction. Bottom: Planet-disc force correction.</figcaption></figure>

From the results of Müller et al. (2012), we chose to look for a SVSL of the form:

$$ \epsilon_g(d_g)/H_g(r) = \sqrt{2} [1 - \exp(-\epsilon_{0,g}d_g^k)], \quad (A.1) $$

with the constraint $L_{\epsilon}^g(d_g = 0) = L_0$. From a Taylor expansion we get:

$$ L_{\epsilon}^{g}(d_g) \sim \frac{\pi d_g^2}{[d_g^2 + \epsilon_{0,g}^2 d_g^{2k}]^{3/2}}. \quad (A.2) $$

The SLFC converges towards $L_0$ if (and only if) $k = 2/3$ and $\epsilon = [\pi/L_0]^{1/3}$. This model can be highly improved thanks to an additional term in the exponential under the form $-\alpha d_g^k$. The values of $(\alpha, n)$ were obtained by numerical fitting to the exact SGFC. In Fig. A.1 (top panel), we compare the exact SGFC (blue solid line) with the SVSL obtained only by analytical means (orange dashed line) or the SVSL, on which we also introduced the numerical fit (black cross markers). We observe that the

unique analytical correction decreases the error to less than 5% for the whole distance range.

The same procedure was used for estimating the SVSL of a planet interacting with a disc. In Fig. A.1 (bottom panel), we compare the SVSL obtained only by analytical means (orange dashed line) with the SVSL, which also includes the numerical fit (black cross markers). Again, the accuracy of the planet-disc SVSL method is excellent compared to the CSL method.

## Appendix B: Commutativity of the self-gravity force correction for different phases

Thanks to the principle of action-reaction applied to the whole dusty disc and the whole gas disc<fnref n="8" /> we get:

$$ \iint_{\text{disc}} F_{\text{sg}}^{d \to g}(\mathbf{r}) d^2 \mathbf{r} = - \iint_{\text{disc}} F_{\text{sg}}^{g \to d}(\mathbf{r}') d^2 \mathbf{r}', \quad (B.1) $$

which implies that $L_{\text{sg}}^{d \to g}(d_g, \eta) = L_{\text{sg}}^{g \to d}(d_g, \eta)$. This commutativity is also easily recovered from Eq. 6:

$$ \begin{align} L_{\text{sg}}^{g \to d}(d_g, \eta) &= \frac{1}{2} \frac{d_d^3}{d_g} \iint_{u,v=-\infty}^{\infty} \frac{e^{-\frac{u^2}{2}} e^{-\frac{v^2}{2}}}{[d_d^2 + (u - \eta v)^2]^{3/2}} du dv \nonumber \\ &\text{but } d_d = \eta d_g \tag{B.2} \\ &= \frac{1}{2} d_g^2 \iint_{u,v=-\infty}^{\infty} \frac{e^{-\frac{u^2}{2}} e^{-\frac{v^2}{2}}}{[d_g^2 + (u/\eta - v)^2]^{3/2}} du dv \nonumber \\ &= L_{\text{sg}}^{d \to g}(d_g, \eta). \nonumber \end{align} $$

As a consequence, throughout the whole paper, we only use the notation $L_{\text{sg}}^{dg}$.

## Appendix C: Derivation of $\lambda$ and $L_0$

The Taylor expansion of Eq. 11, where $a=d$ and $b=g$, in the vicinity of $d_g \sim 0$, is:

$$ L_{\text{sg}}^{dg}(d_g, \eta) \sim_{d_g \to 0} \frac{\pi}{[\lambda(\eta)(\epsilon_{0,g} - \epsilon_{0,p}) + \epsilon_{0,p}]^3}. \quad (C.1) $$

However, the SLFC also satisfies $\lim_{d_g \to 0} L_{\text{sg}}^{dg}(d_g, \eta) = \delta(\eta)L_0$ (Eq. 22). The equalisation of both equations leads to:

$$ \lambda(\eta) = \frac{\epsilon_{0,g} (1/\delta(\eta))^{1/3} - \epsilon_{0,p}}{\epsilon_{0,g} - \epsilon_{0,p}}. \quad (C.2) $$

Since $\delta(\eta = 1) = 1$, it is immediately evident that $\lambda(\eta = 1) = 1$. The constraints of Sect. 3.3.1 require that the $\lambda$ function is cancelled out for infinite gas-to-dust scale height ratios. Considering that $\delta(\eta \to \infty) = \sqrt{2}$, the above condition is only possible if:

$$ \frac{\epsilon_{0,p}}{\epsilon_{0,g}} = \left(\frac{1}{2}\right)^{1/6}, \quad (C.3) $$

where $\epsilon_{0,g} = (L_0/\pi)^{1/3}$ and $\epsilon_{0,p} = (\pi/2)^{1/6}$. The mathematical relation C.3 implies the constraint $L_0 = \sqrt{\pi}$.

<footnote n="8">⁸ This rationale also works for an elementary dust volume and an elementary gas volume – but it does not work between an elementary volume of dust and a disc of gas (and vice-versa).</footnote>