11institutetext: Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
22institutetext: Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE,
Marseille, France 33institutetext: Leibniz-Institut für Astrophysik Potsdam
(AIP), Potsdam, Germany

# Self-gravity in thin-disc simulations of protoplanetary discs:

smoothing length rectified and generalised to bi-fluids

S. Rendon Restrepo 112233<EMAIL_ADDRESS>P. Barge 11<EMAIL_ADDRESS>
###### Abstract

Context. To mimic protoplanetary discs (PPDs) evolution, 2D simulations with
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 gas scale height. On the other hand when a dust component
is included, the question arises as whether the smoothing length approach can
still be used to quantify not only the dust self-gravity but also its
gravitational interaction with gas.

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

Methods. Analytical developments are used to approximate the vertically
averaged self-gravity when the smoothing length is not assumed to be constant
but rather a spatial function.

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
the computational ground, we prescribe the use of tapering functions for
avoiding numerical divergences, checked that our method remains compatible
with standard fast Fourier transform algorithms and evaluated computational
costs.

Conclusions. Our space varying smoothing length permits (i) to solve the
contradictions inherent to a constant smoothing length hypothesis, (ii) to fit
accurately the 3D vertically averaged self-gravity and (iii) is 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
accounting for self-gravity and are important to deepen our understanding of
planetesimals formation and type I migration.

###### Key Words.:

self-gravity – 2D simulations – smoothing length – bi-fluid – hydrodynamics –
Plummer potential

##  1 Introduction

Despite huge developments in 3D numerical computations and the advent of
progressively sophisticated computational facilities, thin-disc (2D)
simulations remain much less expensive and widely used in the study of
protoplanetary discs (PPDs). The 2D approximation lies on the vertical
averaging of the 3D physical quantities and governing equations. When dealing
with self-gravity (SG) the question becomes more sensitive since the equations
cannot be vertically integrated. In practise in PPDs studies an equivalent
smooth potential, called a 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 for accounting 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). Masset (2002); Huré & Pierens (2009) showed that, in the planet case,
the SL should be proportional to the scale height of the gas disc and 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, one should instead
use $\epsilon_{g}/H_{g}=1.2$ to avoid a systematic overestimate of the SG
term. They also showed that when accounting for vertical SG a deviation to the
vertical Gaussian distribution occurs and in this case the SL is also
proportional to the disc Toomre’s parameter.

Müller et al. (2012) prescription clearly leads to small errors ($<2\%$) at
large distances but it also has an important draw-back (never mentioned by the
authors): 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, they didn’t test analytical
expressions that could fit at best the numerical curve. Further, they found
that the best approximation is obtained when the SL vanishes at the
singularity, a statement which 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 unexplained by the authors and, at our knowledge, gone unnoticed until
now.

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 due 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 in the form
$H_{d}=\sqrt{\frac{\tilde{\alpha}}{\tilde{\alpha}+St}}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 $St$ is the Stokes
number. Vertical averaging over the dust layer faces the same numerical issues
than averaging over the gas disc. The problems can be solved in the same way
than for a single gas disc and a simple extrapolation provides the dust disc
SL: $\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}<<\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 present time we do not possess
any rigorous theoretical approach permitting to evaluate this gravitational
bi-fluid interaction. But we think that this problem can also be addressed
with a SL and, intuitively, we can expect that it should lie between the gas
SL and the dust SL. To our knowledge, both issues raised in this paragraph for
bi-fluids are new and were never addressed in numerical and theoretical
studies of PPDs.

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

In this paper, our goal is to provide a method, based on the SL formalism, to
accurately compute SG terms, in thin disks, at short and large separations. We
aim to introduce the substantial SG interaction of gas and dust, when latter
is considered as a pressureless fluid. We begin in Sect. 2 conducting a
theoretical development from first principles which justifies the use of the
Plummer potential. Then, in Sect. 3 we explore SG estimation through the prism
of the SL paradigm. We corrected and completed the Plummer potential formalism
thanks to a spatial dependent SL and generalise it to the case of bi-fluids.
In Sect. 4 we handle practical aspects such as 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 Sect. 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 remind the interest of
the Plummer potential for 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:

$\begin{array}[]{lll}\@vec{f}^{g,tot}_{3D}(\@vec{r},z)&=&\rho_{g}(\@vec{r},z)\left[\<EMAIL_ADDRESS>\@vec{f}^{d,tot}_{3D}(\@vec{r},z)&=&\rho_{d}(\@vec{r},z)\left[\@vec{\nabla}\Phi_{g}+\@vec{\nabla}\Phi_{d}\right]\end{array}$
(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
conciseness, we denote both fluid phases as a and b which reduces each term of
the r.h.s of Eqs. 1 to:

$\begin{array}[]{llll}\@vec{f}^{a\rightarrow
b}_{3D}(\@vec{r},z)&=&\rho_{b}(\<EMAIL_ADDRESS>&=&\displaystyle-G\,\rho_{b}(\@vec{r},z)\\\
&&\displaystyle\iint\limits_{disc}\int\limits_{z^{\prime}=-\infty}^{+\infty}\frac{\rho_{a}(\@vec{r}^{\prime},z^{\prime})}{||\@vec{r}-\@vec{r}^{\prime}||^{2}+(z-z^{\prime})^{2}}\,\@vec{e}_{l}\,d^{2}\@vec{r}^{\prime}\,dz^{\prime}\end{array}$
(2)

where $\@vec{f}^{a\rightarrow b}_{3D}$ is force per unit volume the a-disc
exerts on an elementary b fluid element. The density of phase a is noted
$\rho_{a}$ and
$\@vec{e}_{l}=\left[\,\@vec{r}-\@vec{r}^{\prime}+(z-z^{\prime})\,\@vec{e}_{z}\right]\big{/}\sqrt{||\@vec{r}-\@vec{r}^{\prime}||^{2}+(z-z^{\prime})^{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}(\@vec{r},z)=\rho_{0,g}(\@vec{r})\,e^{-\frac{1}{2}z^{2}/H_{g}(\@vec{r})}$,
where $H_{g}(\@vec{r})$ is the gas pressure scale height. We assume same
vertical Gaussian profile for dust density than for gas but this time the
vertical equilibrium is rather governed by turbulent stirring which sets a
different scale height for dust: $H_{d}(\@vec{r})$. This permits to write both
surface densities as:

$\Sigma_{a}(\@vec{r})=\int\limits_{z=-\infty}^{\infty}\rho_{a}(\@vec{r},z)\,dz=\sqrt{2\pi}\,H_{a}(\@vec{r})\,\rho_{0,a}(\@vec{r})$
(3)

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

$\begin{array}[]{lll}\@vec{f}^{a\rightarrow
b}_{2D}(\@vec{r})&=&\int\limits_{z=-\infty}^{\infty}\@vec{f}^{a\rightarrow
b}_{3D}(\@vec{r},z)\,dz\\\
&=&-G\,\rho_{0,b}(\@vec{r})\iint\limits_{disc}\,\rho_{0,a}(\@vec{r}^{\prime})\,s\,\@vec{e}_{s}\\\
&&\displaystyle\left(\,\,\int\limits_{z,z^{\prime}=-\infty}^{+\infty}\displaystyle\frac{e^{-\frac{1}{2}z^{2}/H_{b}^{2}(\@vec{r})}e^{-\frac{1}{2}z^{\prime
2}/H_{a}^{2}(\@vec{r}^{\prime})}}{\left(s^{2}+(z-z^{\prime})^{2}\right)^{{3}/{2}}}\,dz\,dz^{\prime}\right)d^{2}\@vec{r}^{\prime}\end{array}$
(4)

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

$\@vec{f}^{a\rightarrow
b}_{2D}(\@vec{r})=-\frac{G}{\pi}\,\Sigma_{b}(\@vec{r})\iint\limits_{disc}\,\frac{\Sigma_{a}(\@vec{r}^{\prime})}{H_{g}(\@vec{r})\,s}\,L^{ab}_{sg}(d_{g},d_{b},\eta)\,\@vec{e}_{s}\,d^{2}\@vec{r}^{\prime}$
(5)

where:

$L^{ab}_{sg}(d_{g},d_{b},\eta_{ab})=\displaystyle\frac{1}{2}\frac{d_{b}^{3}(\@vec{r})}{d_{g}(\@vec{r})}\iint\limits_{u,v=-\infty}^{+\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[d_{b}(\@vec{r})^{2}+(u-\eta_{ab}v)^{2}\right]^{3/2}}\,du\,dv$
(6)

is a normalised quantity that we called self-gravity force correction
(SGFC)111Müller et al. (2012) defined as force correction a quantity that they
defined as $I_{sg}$. We noticed that in some of their graphics they mistook
$sI_{sg}(s)$ with $I_{sg}(s)$. Therefore for clarity and consistency with
their work, we did not adopted their naming convention. (with respect to the
3D case), $d_{b}(\@vec{r})=s/H_{b}(\@vec{r})$ is the separation normalised
with respect to b scale height and
$\eta_{ab}=H_{a}(\@vec{r}^{\prime})/H_{b}(\@vec{r})$ is the a-to-b scale
height ratio. In the following we will use the notation
$d_{b}=d_{b}(\@vec{r})$, except when a distinction is necessary. We want to
highlight that in our reasoning $\eta_{ab}$ is a spatial varying quantity, but
for the rest of this article we will rather use:

$\eta_{ab}=\langle H_{a}(\@vec{r}^{\prime})\rangle/\langle
H_{b}(\@vec{r})\rangle$ (7)

where $\langle\rangle$ stands for space averaging over the 2D disc. This
simplifies next theoretical developments. We note that the scale heights used
in this paper could be time-varying functions, but for conciseness, 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 approximate at best this
integral.

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

Both vertical averaging in Eq. 6 cannot be performed analytically but in
practice we can introduce a Plummer potential, $\Psi^{ab}_{\epsilon}$, in the
2D gravitational potential, $\Phi^{ab}_{\epsilon}$, which permits to approach
the 2D SG force defined in Eq. 5:

$\begin{array}[]{lll}\@vec{f}^{a\rightarrow
b}_{2D,\epsilon}(\@vec{r})&=&\displaystyle-\frac{\Sigma_{b}(\@vec{r})}{\pi}\,\@vec{\nabla}_{s}\Phi^{ab}_{\epsilon}(\@vec{r})\\\\[6.0pt]
&\mbox{with}&\displaystyle\Phi^{ab}_{\epsilon}(\@vec{r})=-G\iint\limits_{disc}\Sigma_{a}(\@vec{r}^{\prime})\,\Psi_{\epsilon}^{ab}(s)\,\differential^{2}\@vec{r}^{\prime}\end{array}$
(8)

where:

$\Psi^{ab}_{\epsilon}(s)={\pi}\big{/}\left({s^{2}+\epsilon_{ab}^{2}}\right)^{1/2}$
(9)

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

$||\nabla\Psi^{ab}_{\epsilon}||=\pi
s\big{/}\left({s^{2}+\epsilon_{ab}^{2}}\right)^{3/2}$ (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}$ that is, in
fact, inconsistent with Eq. 10. Indeed, the additional term
$\partial_{s}\epsilon_{ab}^{2}/2$ should be present in the numerator of
aforementioned equation. So, in order to remain mathematically correct and to
keep the possibility to make comparisons with the work of Müller et al. (2012)
we decided to consider the potential satisfying 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:

$\begin{array}[]{ccc}L_{\epsilon}^{ab}(d_{g})&=&H_{g}(\@vec{r})\,s\,||\nabla\Psi_{\epsilon}^{ab}||\\\\[4.0pt]
&=&\displaystyle\frac{\pi
d_{g}^{2}}{\left[d_{g}^{2}+(\epsilon_{ab}(d_{g})/H_{g}(\@vec{r}))^{2}\right]^{3/2}}\end{array}$
(11)

This normalised quantity, which we have called the smoothing length force
correction (SLFC), should fit the SGFC to correctly estimate SG in 2D
simulations. This is only possible choosing wisely a spatially depending SL,
$\epsilon_{ab}(d_{g})$, as it will be depicted in next Sect. The definition
and abbreviation of the main quantities encountered in this paper are given in
Table 1.

Table 1: Definitions and list of abbreviations Abbrev. | Definition/Name | Symbol
---|---|---
SG | Self-gravity |
| Mutual distance or separation | $s=||\@vec{r}-\@vec{r}^{\prime}||$
| Normalised separation for | $d_{b}=s/H_{b}$
| phase b |
SL | Smoothing length |
CSL | Constant smoothing length | $\epsilon_{ab}=const.$
SVSL | Space varying smoothing length | $\epsilon_{ab}(d_{g})$
SGFC | Self-gravity force correction |  $L^{ab}_{sg}$ (Eq. 6)
SLFC | Smoothing length force correction |  $L_{\epsilon}^{ab}$ (Eq. 11)
PDFC | Planet-disc force correction |  $L^{p}$ (Eq. 20)
|  a-to-b scale height ratio | $\eta_{ab}$
| Gas-to-dust scale height | $\eta=\langle H_{g}\rangle/\langle H_{d}\rangle$

##  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
start 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

Figure 1: SG and SL force corrections for gas only.
_Top:_ SGFC (exact value, in blue solid line) and SLFC 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).
_Bottom:_ Relative error between SGFC and SLFC (in %):
$100\cdot|\mbox{SLFC}/\mbox{SGFC}-1|$
For a CSL, $\epsilon_{g}=const.$, the SLFC and SGFC don’t 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).

This case corresponds to a=b=g. From the general definition provided by Eq. 6
the gas SGFC is equal to 222For sake of simplicity in the whole article we
replaced all superscripts $gg$ by $g$.:

$L_{sg}^{g}(d_{g})=\displaystyle\frac{1}{2}d_{g}^{2}\iint\limits_{u,v=-\infty}^{\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[d_{g}^{2}+(u-\eta_{g}v)^{2}\right]^{3/2}}\,du\,dv$
(12)

where $\eta_{g}=\langle H_{g}(\@vec{r}^{\prime})\rangle/\langle
H_{g}(\@vec{r})\rangle=1$ 333We think that Müller et al. (2012) implicitly
assumed this equality in their work.. We present in Fig. 1 top panel above
exact SGFC (blue solid line), evaluated numerically, and the equivalent SLFC
(dashed lines) for three constant smoothing lengths (CSLs):
$\epsilon_{g}(d_{g})/H_{g}(\@vec{r})=[0.0,0.6,1.2]$. In the bottom panel are
showed the respective errors in percentage. These curves testify the retrieval
of Müller et al. (2012) results in the case of a self-gravitating disc which
validates our approach. For long distances the three aforementioned SLFC have
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, less than $2\%$ for $d_{g}\geq 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}\simeq 1.772$. Such
behaviour is indeed intrinsic to the CSL formalism at short distances:

$\begin{array}[]{cccll}\mbox{if}&\epsilon_{g}(d_{g})=0&:&\displaystyle
L_{\epsilon}(d_{g})=\frac{\pi}{d_{g}}&\rightarrow+\infty\\\
\mbox{if}&\epsilon_{g}(d_{g})=const.\neq 0&:&\displaystyle
L_{\epsilon}(d_{g})\underset{d_{g}\rightarrow
0}{\sim}\pi\left(\frac{H_{g}}{\epsilon_{g}}\right)^{3}d_{g}&\rightarrow
0\end{array}$ (13)

Above limits demonstrate that the CSL formalism either leads to 100 % or
infinite errors at short distances when $\epsilon_{g}=const.\neq 0$ and
$\epsilon_{g}=0$, respectively. This divergence between the SGFC and SLFC was
not highlighted by Müller et al. (2012) and a possible reason is that this
behaviour is suitable since it allowed to avoid numerical divergences at
$d_{g}=0$ when $\epsilon_{g}(d_{g})=const.\neq 0$. Even if Müller et al.
(2012) didn’t shine a light on these divergences at short distances they
nevertheless evaluated numerically that a space varying smoothing length
(SVSL) was necessary to fit the exact SGFC, that we reproduced in top panel of
Fig. 2. From this Fig. we can infer that the SVSL should tend to 0 in order
that the SLFC matches the SGFC. But this is in contradiction with Fig. 1 and
Eqs. 13, since the SLFC is supposed to diverge analytically for
$\epsilon_{g}(d_{g})=0$ at short distances. In general, no matter the chosen
CSL, the gas SG was underestimated by a factor 100 at short distances. This
contradiction went unnoticed and we intend to solve it in next paragraph.

Figure 2: Space varying smoothing length for gas only
_Top:_ SVSL obtained numerically by Müller et al. (2012).
_Bottom:_ SVSL obtained in this work by analytical computations and curve
fitting.

The analytic expression of the gas SL which better fits Fig. 2 top panel and
thus permits the SLFC to fit the exact SGFC should fulfil next constraints:

$\begin{array}[]{ccc}\lim\limits_{d_{g}\rightarrow
0}\epsilon_{g}(d_{g})&=&0\\\
\lim\limits_{d_{g}\rightarrow+\infty}\epsilon_{g}(d_{g})&=&\sqrt{2}\,H_{g}(\@vec{r})\\\
L_{\epsilon}^{g}(d_{g}=0)&=&L_{0}\end{array}$ (14)

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

$\displaystyle\epsilon_{g}(d_{g})=\sqrt{2}\,H_{g}(\@vec{r})\,\left[1-\exp\left(-\frac{\epsilon_{g,0}}{\sqrt{2}}d_{g}^{2/3}-\alpha
d_{g}^{n}\right)\right]$ (15)

where $\epsilon_{0,g}=\left[{\pi}/{L_{0}}\right]^{1/3}$ and the numerical
values of $(\alpha,n)$ are gathered in Table 2. In particular, the power $2/3$
and $\epsilon_{0,g}$ allow to reach the finite value $L_{0}$ for the SLFC at
$d_{g}=0$. In Fig. 2 bottom panel we show the SVSL analytic model and in Fig.
1 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 % in the whole distance range. As a matter of
comparison, this correction allows to decrease the error up to factors 200 and
40 at short and long distances, respectively, compared to the CSL where
$\epsilon(d_{g})=1.2H_{g}(\@vec{r})$.

Table 2: Fitting parameters for space varying smoothing lengths and $\delta$ models $(\alpha,n)$ | $(0.04319874,1.14791757)$
---|---
$(\beta,q)$ | $(0.06427627,1.14735482)$
$(\gamma,m)$ | $(0.98362092,0.75552227)$

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

Figure 3: SG and SL force corrections for dust only
_Top:_ SGFC for different values of the gas-to-dust scale-height ratio
($\eta$)
_Bottom:_ Relative error between SGFC and SLFC (in %):
$100\cdot|\mbox{SLFC}/\mbox{SGFC}-1|$.
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}\leq 1.5/\eta$).
For large $\eta$, this could favour gravitational clumping.

In this Sect. we take a=b=d. From the general definition provided by Eq. 6 the
dust SGFC is 444For sake of simplicity in the whole article we replaced all
superscripts $dd$ by $d$.:

$\begin{array}[]{ll}L_{sg}^{d}(d_{g},d_{d})&=\displaystyle\frac{1}{2}\frac{d_{d}^{3}}{d_{g}}\iint\limits_{u,v=-\infty}^{\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[d_{d}^{2}+(u-\eta_{d}v)^{2}\right]^{3/2}}\,du\,dv\\\
&=\displaystyle\frac{1}{2}\eta\left(\eta
d_{g}\right)^{2}\iint\limits_{u,v=-\infty}^{\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[(\eta
d_{g})^{2}+(u-v)^{2}\right]^{3/2}}\,du\,dv\\\ &=\eta L_{sg}^{g}(\eta
d_{g})\end{array}$ (16)

where $\eta_{d}=\langle H_{d}(\@vec{r}^{\prime})\rangle/\langle
H_{d}(\@vec{r})\rangle=1$ and $\eta=H_{g}(\@vec{r})/H_{d}(\@vec{r})$ is the
gas-to-dust scale height ratio that, in general, is greater than unity. Again
for simplicity we rather use in the rest of this article $\eta=\langle
H_{g}(\@vec{r})\rangle/\langle H_{d}(\@vec{r})\rangle$. Above relation
naturally applies to the dust SLFC but it is convenient to write it in two
different ways:

$L_{\epsilon}^{d}(d_{g})=\eta L_{\epsilon}^{g}(\eta d_{g})\quad\mbox{ or
}\quad\epsilon_{d}(d_{g},\eta)=\frac{\epsilon_{g}(\eta d_{g})}{\eta}$ (17)

Both equations are equivalent but the former provides a physical insight while
the latter allows a simple implementation in hydrodynamical codes. In Fig. 3
top and bottom panels are shown 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 for the whole $\eta$ and distances range, 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 gas SLFC, $L_{sg}^{d}(d_{g})=L_{sg}^{g}(d_{g})$, while
at short distances the dust SG is $\eta$ times stronger than gas SG and we get
$L_{sg}^{d}(d_{g}=0)=\eta L_{0}$. This latter equality is important because if
(1) the gas SLFC is used instead of the dust SLFC and (2) if a CSL is used
instead of a SVSL, the dust SG is underestimated by a factor $\sim 100\,\eta$
at short distances. These aspects are of primary interest in an astrophysical
context since usually $\eta\gtrsim 10$ which could have important implications
for planet formation theories as will be discussed in Sect. 5.2.

###  3.3 Crossed contributions

In this Sect. we take a=d and b=g. We highlight that the SGFC is commutative
with respect to phases a and b. This could be demonstrated thanks to Newton’s
third law or by analytical arguments (Appendix B). Therefore, in the
following, we will adopt the notation $L_{sg}^{dg}$ to refer to the dust-gas,
or gas-dust, SGFC. Within this clarification and Eq. 6 the dust-gas SGFC is:

$L_{sg}^{dg}(d_{g},\eta)=\displaystyle\frac{1}{2}d_{g}^{2}\iint\limits_{u,v=-\infty}^{\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[d_{g}^{2}+(u-v/\eta)^{2}\right]^{3/2}}\,du\,dv$
(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 next Sects. we explore for first time what constraints should be
satisfied in order to construct a consistent SL which allows to accurately
approach the exact dust-gas SGFC.

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

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

$\lim\limits_{\eta\rightarrow 1}L_{sg}^{dg}(d_{g},\eta)=L_{sg}^{g}(d_{g})$
(19)

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

$\lim\limits_{\eta\rightarrow+\infty}L_{sg}^{dg}(d_{g},\eta)=\sqrt{2\pi}\,I_{p}(d_{g}^{2}/4)=L^{p}(d_{g})$
(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 second kind. Interestingly, $I_{p}(d_{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)
(see Appendix A), we find the planet-disc SVSL:

$\epsilon_{p}(d_{g})=H_{g}(\@vec{r})\left[1-\exp\left(-\epsilon_{0,p}\,d_{g}^{2/3}-\beta
d_{g}^{q}\right)\right]$ (21)

where $\epsilon_{0,p}=\left(\pi/2\right)^{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. It is interesting to 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

Figure 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$.

Figure 5: SG and SL force corrections for dust-gas and space varying
smoothing-length
_Top:_ Dust-gas SVSL for different values of $\eta$ (gas-to-dust heights
ratio).
_Middle:_ Dust-gas SGFC for the same values of $\eta$
_Bottom:_ Relative error between SGFC and SLFC (in %):
$100\cdot|\mbox{SLFC}/\mbox{SGFC}-1|$.
For thin ($\eta>>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.

To study dust-gas SGFC at short distances it is convenient to define:

$\delta(\eta)=\lim\limits_{d_{g}\rightarrow 0}L_{sg}^{dg}(d_{g},\eta)/L_{0}$
(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_{sg}^{dg}$ matches
$L_{sg}^{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\left[-\gamma(\eta-1)^{m}\right]$
(23)

where 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\limits_{d_{g}\rightarrow+\infty}L_{sg}^{dg}(d_{g},\eta)=\frac{\pi}{d_{g}}$
(24)

This implies that at long distances the dust-gas SVSL should 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

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

$\epsilon_{dg}(d_{g},\eta)=\lambda(\eta)\,\epsilon_{g}(d_{g})+(1-\lambda(\eta))\,\epsilon_{p}(d_{g})$
(25)

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

$\begin{array}[]{ccccc}\lim\limits_{\eta\rightarrow
1}\lambda(\eta)=1&\Longrightarrow&\lim\limits_{\eta\rightarrow
1}\epsilon_{dg}(d_{g},\eta)&=&\epsilon_{g}(d_{g})\\\
\lim\limits_{\eta\rightarrow+\infty}\lambda(\eta)=0&\Longrightarrow&\lim\limits_{\eta\rightarrow+\infty}\epsilon_{dg}(d_{g},\eta)&=&\epsilon_{p}(d_{g})\end{array}$
(26)

Additionally, from Sect. 3.3.2 results and a Taylor expansion (Appendix C), we
get:

$\lambda(\eta)=\displaystyle\frac{\epsilon_{0,g}\left({1}/{\delta(\eta)}\right)^{1/3}-\epsilon_{0,p}}{\epsilon_{0,g}-\epsilon_{0,p}}$
(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 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, $\eta\geq 5$, the SVSL tends rapidly
to the one of the planet-disc interaction case. This permits to use the
approximation: $\epsilon_{dg}(d_{g})\simeq\epsilon_{p}(d_{g})$ for $\eta\geq
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 decimals.

####  3.3.4 Summary for the dust-gas SVSL

In middle and bottom panels of Fig. 5 are shown 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 distances
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 distances
range. Compared to both cases studied in Sects. 3.1 and 3.2 the error is
slightly higher but it remains very acceptable. We want to stress that the
dust-gas SLFC tends rapidly, with respect to $\eta$, to the planet-disc SLFC
which makes possible the simplification:

$L_{sg}^{dg}(d_{g},\eta)=L^{p}(d_{g})\quad\mbox{ or
}\quad\epsilon_{dg}(d_{g},\eta)=\epsilon_{p}(d_{g})\quad\mbox{ if
}\quad\eta\geq 5$ (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 computational costs for 2D, N-fluid simulations with SG.

###  4.1 Removing numerical divergences

Figure 6: Normalised Plummer potential gradient with, and without, a tapering
function.
Under 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\simeq 0.024)$, ensures that most of the SG contribution
is not lost over short distances ($d_{g}\leq 1.5$). Under the CSL approach
(red dotted), the gradient cancels at the singularity but most of the short
distance contribution to the SG is lost.

From Eq. 11 it is obvious that $||\nabla\Psi_{\epsilon}^{ab}||\propto
d_{g}^{-1}$ at short separations which makes necessary 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 not loose the accuracy reached in Sect. 3 the tapering length should
be approximately as large as the finest numerical resolution in the grid:
$r_{s}\sim 3\min(\Delta r,r\Delta\theta)$ 555In cylindrical coordinates.
666The factor 3 was chosen for ensuring a safety margin. where $\Delta r$ and
$\Delta\theta$ are the resolution in the radial and azimutal directions,
respectively. For satisfying above constraints we chose arbitrarily $p=2$ and
the tapering function as following:

$f(d_{g})=1-\exp\left[-\frac{1}{2}\left(\frac{d_{g}}{d_{s}}\right)^{2}\right]\underset{d_{g}\rightarrow
0}{\sim}\frac{1}{2}\left(\frac{d_{g}}{d_{s}}\right)^{2}$ (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 same quantity under a CSL assumption (red
dotted line). Such tapering length choice was motivated by the high resolution
reached by Rendon Restrepo & Barge (2022), i.e 146 cells/$H_{g}$ and 127
cells/$H_{g}$ in the radial and azimuthal directions. As expected both Plummer
potential gradients have 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}\leq 1.5$ but there is no need to resort to a
tapering function since the potential gradient cancels analytically at the
singularity. This last statement 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 to verify if 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 instructive to remind here from ground how SG
is computed, but this time, 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}$) 777We took $G=1$.:

$\begin{array}[]{ll}-\displaystyle\frac{\<EMAIL_ADDRESS>\quad=\displaystyle\iint\limits_{disc}f(d_{g})\frac{\Sigma_{a}(\@vec{r}^{\prime})(\@vec{r}-\@vec{r}^{\prime})}{\left[s^{2}+\epsilon_{ab}(d_{g},\@vec{r})^{2}\right]^{\frac{3}{2}}}\<EMAIL_ADDRESS>\quad=\displaystyle\iint\limits_{disc}f(d_{g})\frac{\Sigma_{a}(\@vec{r}^{\prime})\left(\left(\frac{r}{r^{\prime}}\right)-\cos{(\theta-\theta^{\prime})}\right)\frac{dr^{\prime}}{r^{\prime}}d\theta^{\prime}}{\left[1+\left(\frac{r}{r^{\prime}}\right)^{2}-2\left(\frac{r}{r^{\prime}}\right)\cos{(\theta-\theta^{\prime})}+\delta\epsilon_{ab}^{2}\right]^{\frac{3}{2}}}&\end{array}$
(30)

where:

$\begin{array}[]{lllll}\delta\epsilon_{ab}&=&\displaystyle\frac{\epsilon_{ab}(d_{g},\@vec{r})}{r^{\prime}}&=&\displaystyle\frac{r}{r^{\prime}}h_{g}(r)\frac{\epsilon_{ab}(d_{g},\@vec{r})}{H_{g}(r)}\end{array}$
(31)
$d_{g}=\displaystyle\frac{r^{\prime}}{r}\frac{1}{h_{g}(\@vec{r})}\sqrt{1+\left(\frac{r}{r^{\prime}}\right)^{2}-2\left(\frac{r}{r^{\prime}}\right)\cos{(\theta-\theta^{\prime})}}$
(32)

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

$\begin{array}[]{ll}\displaystyle-\frac{\@vec{\nabla}\Phi_{\epsilon}^{ab}\cdot\@vec{e}_{r}}{\pi}&\\\
\quad\quad=\displaystyle\iint\limits_{disc}\Sigma_{a}(X^{\prime},\theta^{\prime})\,\mathcal{G}_{r}^{ab}(X-X^{\prime},\theta-\theta^{\prime})\,dX^{\prime}d\theta^{\prime}&\end{array}$
(33)

where:

$\begin{array}[]{ll}\mathcal{G}_{r}^{ab}(X-X^{\prime},\theta-\theta^{\prime})&\\\
\quad=\displaystyle\frac{f(X-X^{\prime},\theta-\theta^{\prime})\left[e^{(X-X^{\prime})}-\cos{(\theta-\theta^{\prime})}\right]}{\left[1+e^{2(X-X^{\prime})}-2e^{(X-X^{\prime})}\cos{(\theta-\theta^{\prime})}+\delta\epsilon_{ab}(X-X^{\prime},\theta-\theta^{\prime})^{2}\right]^{\frac{3}{2}}}&\\\
\end{array}$

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

$\begin{array}[]{ll}\displaystyle-\frac{\@vec{\nabla}\Phi_{\epsilon}^{ab}\cdot\@vec{e}_{\theta}}{\pi}&=\displaystyle\iint\limits_{disc}\Sigma_{a}(X^{\prime},\theta^{\prime})\,\mathcal{G}^{ab}_{\theta}(X-X^{\prime},\theta-\theta^{\prime})\,dX^{\prime}d\theta^{\prime}\end{array}$
(34)

where:

$\begin{array}[]{ll}\mathcal{G}^{ab}_{\theta}(X-X^{\prime},\theta-\theta^{\prime})&\\\
\quad=\displaystyle\frac{f(X-X^{\prime},\theta-\theta^{\prime})\,\sin{(\theta-\theta^{\prime})}}{\left[1+e^{2(X-X^{\prime})}-2e^{(X-X^{\prime})}\cos{(\theta-\theta^{\prime})}+\delta\epsilon_{ab}(X-X^{\prime},\theta-\theta^{\prime})^{2}\right]^{\frac{3}{2}}}&\\\
\end{array}$

is the modified azimutal Green function. The rewriting of both integrals,
defined by Eqs. 33 and 34, permits the last writing in terms of Fourier
transforms:

$\begin{array}[]{ccc}\@vec{f}_{2D,\epsilon}^{a\rightarrow
b}(\@vec{r})\cdot\@vec{e}_{r}&=&\Sigma_{b}(\@vec{r})\,\mathcal{F}^{-1}\left[\mathcal{F}\left(\Sigma_{a}\right)\ast\mathcal{F}\left(\mathcal{G}^{ab}_{r}\right)\right]\\\\[6.0pt]
\@vec{f}_{2D,\epsilon}^{a\rightarrow
b}(\@vec{r})\cdot\@vec{e}_{\theta}&=&\Sigma_{b}(\@vec{r})\,\mathcal{F}^{-1}\left[\mathcal{F}\left(\Sigma_{a}\right)\ast\mathcal{F}\left(\mathcal{G}^{ab}_{\theta}\right)\right]\end{array}$
(35)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier transform operator
and its inverse. The symbol $\ast$ denotes the convolution operator. In
practice such quantities are computed numerically thanks to fast Fourier
modules which highly accelerate numerical computation. 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 made possible if (1) a
logarithmic mesh is used in the radial direction. This condition stems from
the radial variable substitution which enabled us to obtain the formulation
with the modified Green’s functions. Furthermore, (2) the function
$\Sigma_{a}$ must be periodic in radial and azimutal directions in order to be
able to use Fourier transforms. Periodicity is obviously satisfied in the
azimutal direction but in the radial direction this is not necessarily the
case. Such periodicity is artificially ensured thanks to a zero-padding: 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 that (3) the disc aspect ratio of gas and dust must
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 the opportunity to clarify that the aforementioned assumption of a
constant disc aspect ratio used in this Sect. only permits to resort to FFT
methods for accelerating 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, direct summation in the radial direction and
Fourier transforms in the azimuth direction, ensured by the periodicity, are a
straightforward solution.

###  4.3 Computational costs

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

Gas into gas

$\mathcal{F}(\Sigma_{g})\quad\mbox{and}\quad\mathcal{F}^{-1}\left[\mathcal{F}\left(\Sigma_{g}\right)\ast\mathcal{F}\left(\mathcal{G}^{g}_{r}\right)\right]$
(36)

Dust into dust

$\mathcal{F}(\Sigma_{d})\quad\mbox{and}\quad\mathcal{F}^{-1}\left[\mathcal{F}\left(\Sigma_{d}\right)\ast\mathcal{F}\left(\mathcal{G}^{d}_{r}\right)\right]$
(37)

Dust into gas

$\mathcal{F}(\Sigma_{d})\quad\mbox{and}\quad\mathcal{F}^{-1}\left[\mathcal{F}\left(\Sigma_{d}\right)\ast\mathcal{F}\left(\mathcal{G}^{dg}_{r}\right)\right]$
(38)

Gas into dust

$\mathcal{F}(\Sigma_{g})\quad\mbox{and}\quad\mathcal{F}^{-1}\left[\mathcal{F}\left(\Sigma_{g}\right)\ast\mathcal{F}\left(\mathcal{G}^{dg}_{r}\right)\right]$
(39)

From above recapitulation we observe that the Fourier transforms of gas and
dust densities are duplicated which reduces the total number of Fourier
operations to 6 during a numerical treatment. Therefore, the numerical
endeavour for computing SG for a bi-fluid is 3 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 amount of Fourier transforms is:

$N^{2}+N=\displaystyle\underbrace{\left(2N^{2}\right)}_{\mbox{All $N$-tuples
combinations}}-\underbrace{\left(N^{2}-N\right)}_{\mbox{Duplicates}}$ (40)

Compared to standard self-gravitating simulations on a unique gaseous phase,
the computation time 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 same scale heights.

##  5 Discussion

In this Sect. 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 well known that planets can migrate due to tidal interactions with the
gas disc. This is the case of type I migration in which 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}(\@vec{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.5H_{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 check
whether (2) the widely used $\epsilon_{p}/H_{g}(\@vec{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 require,
though, to make some improvements: the vertical layering due to SG must be
accounted, 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 planet 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 way (Varga et al., 2021), numerical simulations have not yet concluded
that vortices could form objects bond 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. Particularly,
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 $\sim 2000$ for $\eta=20$ at short separations. At
the same time the estimation of the dust-gas SLFC could permit to quantify
with an acceptable accuracy a possible gaseous envelope capture by dust
clumps. From theoretical analysis it was found that SG inhibit vortices
formation by Rossby Wave Instability (Lovelace & Hohlfeld, 2013) which was
later confirmed by numerical simulations (Baruteau & Zhu, 2016; Regály &
Vorobyov, 2017; Tarczay-Nehéz et al., 2021). 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. At the light of our
findings, previous simulations results should be checked anew for
understanding at which extent the SVSL affects theoretical predictions on
vortices survival in self-gravitating PPDs.

###  5.3 Limitations, improvements and tests

The limitations of our model are inherent to the initial assumptions we 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,
that’s 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 the 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
$\sim$ 2 (Baehr & Zhu, 2021) and we expect that a strong dust layering will
also modify gas vertical structure in correlation with the dust-to-gas density
ratio.

The global vertical structure of PPDs evolves in time due to accretion heating
(Schobert et al., 2019) and stellar irradiation (Wu & Lithwick, 2021), amongst
others. 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, but is compatible
with the global and local time evolution of the vertical structure, provided
that the time dependence of 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}(\@vec{r},t)=h_{a}(t)$, the computational benefits from the FFT method
could be lost. Indeed, this would require 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 3 Fourier transforms at each time step instead of 2. We think 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 which 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 is 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 $\lesssim 1$)
(Garaud et al., 2004). We also assumed that small-scale turbulence is
maintained in the gas disc by a mechanism that we disregarded. 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 in 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. This simulation is out of the scope of
current paper but still an interesting lead that the authors want to explore
in a near future.

##  6 Conclusions

We revealed 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}\leq 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}=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}(\@vec{r})\,\left[1-\exp\left(-\frac{\epsilon_{0,g}}{\sqrt{2}}d_{g}^{2/3}-\alpha
d_{g}^{n}\right)\right]$ (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}(\@vec{r})}{\eta}\,\left[1-\exp\left(-\frac{\epsilon_{0,g}}{\sqrt{2}}(\eta
d_{g})^{2/3}-\alpha(\eta d_{g})^{n}\right)\right]$ (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}(\@vec{r})\left[1-\exp\left(-\epsilon_{0,p}\,d_{g}^{2/3}-\beta
d_{g}^{q}\right)\right]$ (43)

where $\epsilon_{0,p}=\left(\frac{\pi}{2}\right)^{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})$
(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 to match the SGFC with a high accuracy, even at the
singularity ($d_{g}\rightarrow 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 the 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 gravitational binding of the dust clumps
trapped in the vortex and, on the other hand, dust-gas SLFC enables 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 & Zhu (2021) Baehr, H. & Zhu, Z. 2021, ApJ, 909, 136
  * Barge & Rendon Restrepo (2023) Barge, P. & Rendon Restrepo, S. 2023, submitted
  * Barge & Sommeria (1995) Barge, P. & Sommeria, J. 1995, A&A, 295, L1
  * Baruteau & Masset (2008) Baruteau, C. & Masset, F. 2008, ApJ, 678, 483
  * Baruteau & Masset (2013) Baruteau, C. & Masset, F. 2013, in Lecture Notes in Physics, Berlin Springer Verlag, ed. J. Souchay, S. Mathis, & T. Tokieda, Vol. 861, 201
  * Baruteau & Zhu (2016) Baruteau, C. & Zhu, Z. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3927
  * Crida et al. (2009) Crida, A., Masset, F., & Morbidelli, A. 2009, ApJ, 705, L148
  * Dong et al. (2011) Dong, R., Rafikov, R. R., Stone, J. M., & Petrovich, C. 2011, ApJ, 741, 56
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Frigo & Johnson (2005) Frigo, M. & Johnson, S. G. 2005, Proceedings of the IEEE, 93, 216, special issue on “Program Generation, Optimization, and Platform Adaptation”
  * Garaud et al. (2004) Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292
  * Goldreich & Ward (1973) Goldreich, P. & Ward, W. R. 1973, ApJ, 183, 1051
  * Huré & Pierens (2009) Huré, J. M. & Pierens, A. 2009, A&A, 507, 573
  * Johansen & Youdin (2007) Johansen, A. & Youdin, A. 2007, ApJ, 662, 627
  * Li et al. (2009) Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52
  * Lodato (2007) Lodato, G. 2007, Nuovo Cimento Rivista Serie, 30, 293
  * Lovelace & Hohlfeld (2013) Lovelace, R. V. E. & Hohlfeld, R. G. 2013, MNRAS, 429, 529
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165
  * Masset (2002) Masset, F. S. 2002, A&A, 387, 605
  * Meheut et al. (2010) Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31
  * Müller et al. (2012) Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123
  * Regály & Vorobyov (2017) Regály, Z. & Vorobyov, E. 2017, MNRAS, 471, 2204
  * Rendon Restrepo & Barge (2022) Rendon Restrepo, S. & Barge, P. 2022, A&A, 666, A92
  * Rendon Restrepo et al. (2022) Rendon Restrepo, S., Barge, P., & Vavrik, R. 2022, arXiv e-prints, arXiv:2207.04252
  * Riols & Latter (2018) Riols, A. & Latter, H. 2018, MNRAS, 476, 5115
  * Schobert et al. (2019) Schobert, B. N., Peeters, A. G., & Rath, F. 2019, ApJ, 881, 56
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  * Suriano et al. (2017) Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2017, MNRAS, 468, 3850
  * Surville (2013) Surville, C. 2013, PhD thesis, Aix Marseille Université, thèse de doctorat dirigée par Barge, Pierre Astrophysique et Cosmologie Aix-Marseille 2013
  * Tanaka et al. (2002) Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257
  * Tarczay-Nehéz et al. (2021) Tarczay-Nehéz, D., Rozgonyi, K., & Regály, Z. 2021, MNRAS[arXiv:2111.09073]
  * Varga et al. (2021) Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56
  * Weber et al. (2019) Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178
  * Wu & Lithwick (2021) Wu, Y. & Lithwick, Y. 2021, ApJ, 923, 123

##  Appendix A Gas and planet SVSL derivation

Figure 7: Gas and planet-disc force corrections.
_Top:_ Gas force correction.
_Bottom:_ Planet-disc force correction.
For each force correction we plotted the exact SGFC evaluated numerically
(blue solid line), the SLFC approximation accounting only the analytical
correction (orange dashed line) or accounting an additional numerical fit
(black cross markers).

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

$\epsilon_{g}(d_{g})/H_{g}(r)=\sqrt{2}\left[1-\exp\left(-\epsilon_{0,g}d_{g}^{k}\right)\right]$
(45)

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

$L_{\epsilon}^{g}(d_{g})\underset{d_{g}\rightarrow 0}{\sim}\frac{\pi
d_{g}^{2}}{\left[d_{g}^{2}+\epsilon_{0,g}^{2}d_{g}^{2k}\right]^{3/2}}$ (46)

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}^{n}$. The
values of $(\alpha,n)$ were obtained by numerical fitting to the exact SGFC.
In Fig. 7 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.

Same procedure was used for estimating the SVSL of a planet interacting with a
disc. In Fig. 7 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 SGFC for different phases

Thanks to the principle of action-reaction applied to the whole dusty disc and
the whole gas disc 888This rationale also works for an elementary dust volume
and an elementary gas volume. But not between an elementary volume of dust and
a disc of gas (and vice-versa). we get:

$\iint\limits_{disc}\@vec{F}_{sg}^{d\rightarrow
g}(\@vec{r})d^{2}\@vec{r}=-\iint\limits_{disc}\@vec{F}_{sg}^{g\rightarrow
d}(\@vec{r}^{\prime})d^{2}\@vec{r}^{\prime}$ (47)

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

$\begin{array}[]{lll}L_{sg}^{g\rightarrow
d}(d_{g},\eta)&=&\displaystyle\frac{1}{2}\frac{d_{d}^{3}}{d_{g}}\iint\limits_{u,v=-\infty}^{\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[d_{d}^{2}+(u-\eta
v)^{2}\right]^{3/2}}\,du\,dv\\\\[20.0pt] &&\quad\mbox{but}\quad d_{d}=\eta
d_{g}\\\\[5.0pt]
&=&\displaystyle\frac{1}{2}d_{g}^{2}\iint\limits_{u,v=-\infty}^{\infty}\frac{e^{-\frac{u^{2}}{2}}e^{-\frac{v^{2}}{2}}}{\left[d_{g}^{2}+(u/\eta-v)^{2}\right]^{3/2}}\,du\,dv\\\
&=&L_{sg}^{d\rightarrow g}(d_{g},\eta)\end{array}$ (48)

As a consequence, in the whole paper we only use the notation $L_{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:

$\begin{array}[]{lcl}L_{sg}^{dg}(d_{g},\eta)&\underset{d_{g}\rightarrow
0}{\sim}&\displaystyle\frac{\pi}{\left[\lambda(\eta)(\epsilon_{0,g}-\epsilon_{0,p})+\epsilon_{0,p}\right]^{3}}\end{array}$
(49)

But the SLFC also satisfies $\lim\limits_{d_{g}\rightarrow
0}L_{sg}^{dg}(d_{g},\eta)=\delta(\eta)L_{0}$ (Eq. 22). The equalisation of
both equations leads to:

$\lambda(\eta)=\displaystyle\frac{\epsilon_{0,g}\left({1}/{\delta(\eta)}\right)^{1/3}-\epsilon_{0,p}}{\epsilon_{0,g}-\epsilon_{0,p}}$
(50)

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

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

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