© The Authors 2026

1 Introduction

With the current advancements in numerical capabilities and optimised codes, conducting complex simulations of physical processes within protoplanetary discs (PPDs) has become routine. Despite the increased computational power, high-resolution 3D simulations present practical challenges such as data storage, analysis, or visualisation that remain unresolved. More significantly, they remain computationally challenging, necessitating a resolution trade-off across the three spatial directions. As a result, when studying discs numerically, one of the common approaches involves employing a 2D approximation, which neglects the vertical dimension, and hence offers a significant gain of resolution in the (r, ϕ) plane. In particular, this approximation allows one to conduct sufficiently resolved global simulations including the gravitational force of the disc exerted onto itself, commonly referred as self-gravity (SG).

Deriving the ‘correct’ approximation for the effective gravity in a 2D disc is a subtle and tedious problem that requires one to first define what is meant by 2D approximation. As an extreme case, this could mean living in a ‘2D Universe’, which implies that the laws of nature are fundamentally different from our 3D space; that is, the gravitational field generated by a point mass scales as ∝1/r (Landy 1996). Of course, this does not apply in the case of discs, as they are objects embedded within a 3D space. But based on theoretical expectations, and supported by PPD observations (Pinte et al. 2016), cataclysmic variable accretion discs (Baptista et al. 2016), and galactic discs (Hu & Peng 2014; van Dokkum et al. 2019), we expect discs to have low geometrical aspect ratios. This justifies the flat, or razor-thin, disc assumption for the density vertical profile, which adopts a Dirac δ distribution. In the particular case of axisymmetric discs, this assumption permits one to derive the gravitational potential in the equatorial plane. For instance, this could be achieved by means of a continuous overlap of infinitely flattened homoeoid shells (Binney & Tremaine 2008, Sects. 2.5.1 and 2.6) or through direct integration, resulting in a closed form for the associated potential (Durand 1953; Huré & Hersant 2007; Huré et al. 2008). Although almost flat, PPDs are inherently 3D structures, yet due to computational limitations or for the sake of analytical simplicity, the oversimplified Dirac profile assumption is often made. To eliminate this limitation, the other approach, which we endorse in this paper, consists of averaging vertically the 3D SG force in a thin disc (Li et al. 2009; Müller et al. 2012; Rendon Restrepo & Barge 2023). As a result the obtained 2D force, meant to act in the midplane, encapsulates on average all information from the vertical structure of the disc. We emphasise that the thin disc structure originates from the vertical hydrostatic equilibrium, which, in the simplest case, leads to a Gaussian profile (Armitage 2022).

To account for the disc’s vertical thickness in the 2D SG calculation and to presumably avoid numerical singularities, the gravitational potential of a thin disc was often approximated by a Plummer potential¹; that is, $${\Psi }_{\text{Plumm}}\left(\mathit{r}\right)=-{\int \int }_{\text{disc}}\frac{\Sigma \left({\mathit{r}}^{\text{'}}\right)}{\sqrt{{|\mathit{r}-{\mathit{r}}^{\text{'}}|}^2+{ϵ}^2}}{d}^2{\mathit{r}}^{\text{'}}$$(1) where the smoothing length (SL), denoted as ε, was regarded during a considerable period of time as a free parameter. Over the years, several values have been proposed and theoretical findings have converged towards a softening length proportional to the disc’s scale height (Huré & Pierens 2009; Huré & Hersant 2011; Huré & Trova 2015). The prevailing value frequently used in numerical simulations is ε_g/H_g=0.6−1.2 (Müller et al. 2012). We highlight that Li et al. (2009) identified another prescription that does not utilise a softening prescription. However, it has received little attention, likely due to its apparent complexity and the lack of efficient, compatible numerical methods. In the latest work to date by Rendon Restrepo & Barge (2023), it was shown that while this widely used prescription correctly models the long-range SG interaction, it tends to underestimate the mid- and short-range interaction by up to 100% – a discrepancy that aligns with the removal of the Newtonian behaviour in the presence of softening (Adams et al. 1989; Hockney & Eastwood 2021; Baruteau & Masset 2008; Young & Clarke 2015). As a consequence, these authors introduced a space-varying smoothing length (SVSL), departing from the simple constant SL approach. This SVSL allows the 2D force to better fit the vertically averaged 3D counterpart. Further, they generalised their correction to bi-fluids; that is, the case in which a dust phase is embedded in the gas. They found that two additional SLs are required to account for the gravitational interactions of dust with dust and the crosswise interaction of dust with gas. This could be of great importance, since the low levels of turbulence in the vertical direction – supported both by observations of thin dust layers (Villenave et al. 2020, 2022) and by 3D shearing box simulations (Baehr & Zhu 2021) – may suggest that dust SG is comparable to the SG of gas at small scales.

Even if Rendon Restrepo & Barge (2023) solved the main issue with their new prescription, they nonetheless used a number of mathematical simplifications that are debatable. In particular, they assumed that the relative scale heights between component a and b, expressed via η_ab=H_a(r^′)/H_b(r), are constant over the 2D disc. This, however, breaks the r/r^′ symmetry in the SG force. This symmetry breaking results in a spurious radial acceleration, which needs to be compensated for (Baruteau 2008; Rometsch et al. 2024). Furthermore, the authors disregarded the effect of SG on the vertical direction, which could significantly affect the layering of the disc (Bertin & Lodato 1999; Lodato 2007), and thus also indirectly the 2D SG force. As a consequence, Rendon Restrepo & Barge (2023) prescription is only valid for comparatively light discs, whose Toomre’s parameter satisfies Q_g ≳ 5, and cannot be used self-consistently in simulations of gravitational instability (GI) in which Toomre’s parameter is expected to reach values smaller than unity (Gammie 2001; Takahashi et al. 2016; Baehr et al. 2017; Béthune et al. 2021). In a broader context, this could profoundly impact a broad range of planet formation scenarios, since at some point or another the SG force should come into play in maintaining the gravitational cohesion of the formed objects (Rendon Restrepo & Gressel 2023). Finally, the work of Rendon Restrepo & Barge (2023) should be considered as a preliminary step in view of the dust stratification, since it was presumed to be Gaussian without further details being provided. This can be readily contradicted for strong SG of dust whereby the sedimentation follows a sech²−profile (Klahr & Schreiber 2020, 2021). Furthermore, the mentioned assumption disregards the gravitational contribution of the gas to the vertical layering. We conclude that all issues and ambiguities raised in this paragraph need to be corrected in order to conduct consistent 2D simulations of single or bi-fluids incorporating SG.

In this first paper of a series of two, we aim to provide an exact SG kernel that must be used in 2D global simulations of discs – a formulation that is valid from light to massive discs all at once. In the second accompanying paper, we study how this new kernel affects the GI planet formation scenario. We begin by recalling briefly the vertical stratification of self-gravitating, bi-fluid discs and provide the exact SG kernel in Section 2. Then, in Section 3, we highlight the relevance of this kernel and discuss how it possibly leads to a new runaway situation. In Section 4, we benchmark our kernel against well-established analytical results from thin disc literature and both 2D and 3D numerical tests. Finally, in Section 5 we propose a discussion of numerical methods, the limitations of our approach, the indirect term, and the consistency of the 2D Poisson equation. Important equations and derivations can be found in the appendices.

2 Self-gravity for two-phase discs

Estimating in-plane SG forces in thin discs involves a two-step process. The initial step, often overlooked, involves determining the vertical hydrostatic equilibrium of the system. The resulting stratification is subsequently used, in the final phase, as an input for vertically integrating all forces; including, in our specific case, SG.

Accordingly, in the next sections we recall the results of Rendon Restrepo et al. (2025) regarding the stratification of a bi-fluid disc and the key steps taken by Rendon Restrepo & Barge (2023) for computing SG in 2D, bi-fluid simulations of PPDs. Both results combined enable us to propose an exact SG kernel that works for a broad range of Toomre’s parameters, as well as a bi-fluid. At the same time, it frees us from the problems inherent in a Plummer potential. To enhance readability of this paper, we have compiled the definitions and abbreviations of the main quantities in Table 1.

2.1 The thin self-gravitating, bi-fluid disc

For simplicity, it is often assumed that the vertical structure of the gas and dust phases composing a PPD is set either by the balance between the star’s gravity and by pressure, or by viscous stirring, respectively. However, the vertical component of the disc’s SG should be included in this equilibrium and particularly when modelling class 0/I discs or the outskirts of PPDs where the SG force is predominant over the star’s gravity. This is most likely the case in the Elias 2-27 system (Pérez et al. 2016; Parker et al. 2022) and in IM Lupi and GM Aur (Lodato et al. 2023).

When the gravitational contribution of gas and dust are treated separately in massive discs, their respective layering are no longer Gaussian, but are set by a sech² function (Spitzer 1942; Klahr & Schreiber 2020, 2021). Additionally, when the disc is solely composed of gas, it is possible to connect the stratification of light and massive discs using a ‘biased’ Gaussian stratification. Such an approach incorporates all SG information into a modified scale height dependent on Toomre’s parameter, as shown by Bertin & Lodato (1999). Rendon Restrepo et al. (2025) extended their work and showed that it is not possible to disentangle the contribution of gas and dust to SG, and proposed a very accurate and general solution valid for any SG strength, of gas and/or dust. At first order, their solution is $$\{\begin{array}{c}{\rho }_g(\mathit{r},z)=\frac{{\Sigma }_g}{\sqrt{2\pi }{H}_g^{sg}}exp[-\frac{1}{2}{(z/H_g^{sg})}^2]\hfill \\ {\rho }_d(\mathit{r},z)=\frac{{\Sigma }_d}{\sqrt{2\pi }{H}_d^{sg}}exp[-\frac{1}{2}{(z/H_d^{sg})}^2]\hfill \end{array},$$(2) with $$\{\begin{array}{cc}{H}_g^{sg}\hfill & =\sqrt{\frac{2}{\pi }}{H}_{g}f\left(Q_{\mathrm{bi}-\mathrm{fluid}}\right)\hfill \\ H_d^{sg}\hfill & =\sqrt{\frac{2}{\pi }}{H}_{d}f\left(Q_{\mathrm{bi}-\mathrm{fluid}}\right)\hfill \\ Q_{\mathrm{bi}-\mathrm{fluid}}\hfill & ={(\frac{1}{Q_g}+\frac{1}{Q_d})}^{-1}\hfill \\ f\left(x\right)\hfill & =\frac{\pi }{4x}[\sqrt{1+\frac{8x^2}{\pi }}-1]\hfill \end{array},$$(3) where H_g is the thermal gas scale height and H_d is the dust diffusive scale height. The quantities $Q_g=\frac{c_g{\Omega }_K}{\pi G{\Sigma }_g}$ and $Q_d=\frac{c_d{\Omega }_K}{\pi G{\Sigma }_d}$ are the Toomre’s parameter of gas and dust, respectively. The effective dust sound speed is c_d. Meanwhile, the term f(Q_bi-fluid), present in the definition of both scale heights, indicates that both fluids experience an equal gravitational influence from the star and from their combined mass distributions. Indeed, the crosswise gravity of both fluids is encapsulated by the general Toomre’s parameter, Q_bi-fluid, defined as the harmonic average of its constituents. For a detailed discussion, we invite readers to refer to Rendon Restrepo et al. (2025). Contrary to the razor-thin disc assumption mentioned in the introduction, the approach discussed in this section is able to accommodate any vertical stratification for the gas and dust phases. For instance, in the case of uniquely dust, it could be applicable to systems with limited sedimentation (Lin et al. 2023) to high settling (Villenave et al. 2020, 2022). As we demonstrate in the following section, the vertical profiles exhibited in Eq. (2) can be elegantly exploited for obtaining the exact 2D SG force via a kernel description.

2.2 Exact self-gravity kernel

We consider two fluids (gas and dust) with the Gaussian vertical stratification defined by Eqs. (2)–(3). In 3D, the SG force per unit volume exerted by the PPD on gas and dust parcels are expressed as $$\begin{array}{rl}& {\mathit{f}}_{3D}^{g,\text{tot}}(\mathit{r},z)=-{\rho }_g(\mathit{r},z)[\mathrm{\nabla }{\mathrm{\Phi }}_g+\mathrm{\nabla }{\mathrm{\Phi }}_d],\\ & {\mathit{f}}_{3D}^{d,\text{tot}}(\mathit{r},z)=-{\rho }_d(\mathit{r},z)[\mathrm{\nabla }{\mathrm{\Phi }}_g+\mathrm{\nabla }{\mathrm{\Phi }}_d],\end{array}$$(4) where Φ_g and Φ_d are the distinct gravitational potentials associated with the gas and dust discs, respectively. For the sake of generality and conciseness, we denote the two fluid phases as ‘a’ and ‘b’ for the remainder of the paper. The 2D SG force per unit volume exerted by the disc made of fluid ‘a’ on a volume element of fluid ‘b’ is (see Eq. (4) in (Rendon Restrepo & Barge 2023)): $$\begin{array}{rl}{\mathit{f}}_{2D}^{a\to b}(\mathit{r})=-& {\rho }_b(\mathit{r},z=0){\iint }_{disc}{\rho }_a({\mathit{r}}^{\mathrm{\prime }},z=0)s{\mathit{e}}_s\\ & \left({\int }_{z,z^{\mathrm{\prime }}=-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{e^{-\frac{1}{2}{(z/H_b^{sq}(\mathit{r}))}^2}{e}^{-\frac{1}{2}{\left(z^{\mathrm{\prime }}/H_a^{sq}\left({\mathit{r}}^{\mathrm{\prime }}\right)\right)}^2}}{{(s^2+{(z-z^{\mathrm{\prime }})}^2)}^{3/2}}dzd{z}^{\mathrm{\prime }}\right)d^2{\mathit{r}}^{\mathrm{\prime }}\end{array}$$(5) where s=|r−r^′| stands for the separation between two fluid elements and e_s=(r−r^′)/s is a unit vector in the 2D plane. We highlight that all variables linked to phase ‘a’ and ‘b’ are dependent on r^′ and r positions, respectively. Hence, for the rest of the paper we skip all explicit dependencies on position, except when a distinction is necessary. We highlight that for a disc made of gas and dust, our naming convention leads to four possible combinations (a, b)={(gas, gas), (gas, dust), (dust, dust), (dust, gas)} for computing ${\mathit{f}}_{2D}^{a\to b}$. With these clarifications, Eq. (5) can be rearranged in the following manner: $${\mathit{f}}_{2D}^{a\to b}\left(\mathit{r}\right)=-{\Sigma }_b{\int \int }_{disc}{\Sigma }_a{\mathrm{K}}_{ab}{\mathit{e}}_s{d}^2{\mathit{r}}^{\text{'}},$$(6) where $${\mathrm{K}}_{ab}=\frac{1}{2\pi }\frac{s}{H_a^{sg}{H}_b^{sg}}{\int \int }_{z,z^{\text{'}}=-\infty }^{+\infty }\frac{e^{-\frac{1}{2}{(z/H_b^{sg})}^2}{e}^{-\frac{1}{2}{(z^{\text{'}}/H_a^{sg})}^2}}{{(s^2+{(z-z^{\text{'}})}^2)}^{3/2}}dzd{z}^{\text{'}}$$(7) is the SG force kernel between fluids ‘a’ and ‘b’ ². The column densities are Σ_i. Previously, Rendon Restrepo & Barge (2023) attempted to provide an accurate approximation of the integral defined by Eq. (7) employing a SVSL, which, despite its accuracy, remains burdensome. In particular, it does not respect the r/r^′ symmetry and is not valid in the regime where SG starts to become significant (i.e. Q ≲ 5). Luckily, when the vertical component of SG is disregarded, the double integral defining the SG force kernel can be assessed analytically (Li et al. 2009, Eqs. (20)–(22)) leading us to the derivation of the force kernel in the general case of a self-gravitating bi-fluid as seen below (see Appendices A and B for a detailed derivation): $${\mathrm{K}}_{ab}=\frac{1}{\sqrt{\pi }}{\left(H_{ab}^{sq}\right)}^{-2}\frac{d_{ab}}{8}exp\left(\frac{d_{ab}^2}{8}\right)[K_1\left(\frac{d_{ab}^2}{8}\right)-K_0\left(\frac{d_{ab}^2}{8}\right)],$$(8) where K_α are modified Bessel functions of the second kind and order α. The notation $d_{ab}=s/H_{ab}^{sg}$ stands for the normalised distance between two fluid elements. The root mean square of $H_a^{sg}$ and $H_b^{sg}$ is defined by $$H_{ab}^{sg}=\sqrt{\frac{H_a^{sg2}+H_b^{sg2}}{2}}.$$(9) The above-mentioned distance H_{a b} not only combines the different fluids contribution to the scale height but also takes into account the different disc positions, r and r^′. We highlight that the unique definition of the above root mean square (rms) scale height (Eq. (9)) permits one to lift the r/r^′ symmetry issues associated with a 2D approximation, even when using the Plummer potential prescription. The symmetry issues are avoided because Newton’s third law is fulfilled as long as the force between two density columns is symmetric.

The closed form expressed in Eq. (8) offers several advantages over the one proposed by Rendon Restrepo & Barge (2023). Despite employing special mathematical functions, this formula remains both simple and general. Notably, there is no longer a necessity for the auxiliary functions (λ, δ) introduced when using a SVSL. Furthermore, Eq. (8) encompasses all mutual interactions between different phases embedded in a PPD. Finally, this new kernel formulation achieves the r/r^′ position symmetry, implying Newton’s third law, as well as the interchangeability symmetry between fluid elements ‘a’ and ‘b’. All these aspects reassure us in the accurate representation of the 2D SG when employing our new kernel. We stress that our work differs from Li et al. (2009) due to the incorporation of the bi-fluid analysis (readily adaptable to N-fluids, see Sect. 5.2 for a discussion) and our consideration of how the SG of both components affects their vertical density profile. Indeed, in our investigation, the vertical component of SG is encapsulated in the definition of the revised scale height (Eq. (3)). We stress that these two last aspects are novel. The chosen approach makes our method similar to a 1 D+2 D problem, justified by the fact that the timescales of vertical settling are much smaller than the in-plane dynamic timescales. Further, we want to highlight that due to computational costs, Li et al. (2009) did not encourage to use the Bessel kernel but instead an interpolated formula. Presumably, they did not contemplate that fast Fourier transform (FFT) methods were still valid for such a prescription and that, under certain circumstances, the evaluation of the Bessel kernel is only done once for a whole simulation. This is regrettable because there are a lot of physical lessons to be drawn from their formula, which we explore in the next section.

3 Advantages associated with the Bessel kernel

In this section, we explore the benefits of putting our prescription to use. To achieve this, we examine how it compares to other prescriptions from the literature, assess how SG is rescaled when considering a dust phase, and demonstrate that our kernel is the only one enabling a potential runaway process at infinitesimal distances.

Fig. 1 Normalised SG kernels, $\pi d_g{H}_g^2\mathrm{K}$, for light (top) and massive (bottom) discs. For light discs, the prescriptions of Müller et al. (2012) and Rendon Restrepo & Barge (2023) are insensitive to the Toomre’s parameter. This is in contrast with the Bessel prescription, which scale as 1/Qg in the limit of massive discs. When a constant SL prescription is used, the kernel vanishes at short distances, whereas for a vanishing SL, the kernel diverges, suggesting an overestimation of SG. The Bessel prescription is the unique prescription that is simultaneously inversely proportional to Qg for massive discs, does not vanish at the singularity and permits a smooth transition between massive and light discs.

3.1 Comparison with other kernels from the literature

In order to show the interest of our findings, we depicted in Fig. 1 a comparison between our Bessel kernel and the prescriptions suggested by Müller et al. (2012) and Rendon Restrepo & Barge (2023). We remind that the first approach relies on a constant SL, which we set at ε/H_g ∈[0, 0.3, 1.2], three customary values in several numerical studies (Paardekooper et al. 2011; Young & Clarke 2015; Zhu & Baruteau 2016; Baruteau & Zhu 2016; Vorobyov & Elbakyan 2018). The second approach is intended for weak SG conditions, Q_g ≳ 5, and employs a SVSL. For simplicity, this comparative study assumes a gas-only disc. To ensure finite values at the singularity (d_g=0) we normalised the different kernels by ${\pi }^{-1}{d}_g^{-1}{H}_g^{-2}$. In particular, the last factor allows us to account for modifications in the scale height due to SG whenever it is incorporated in the kernel’s prescription – as is the case only for our Bessel kernel. In top panel of Fig. 1, we compare the kernels designed for application in the weak SG regime, while the bottom panel shows the Bessel kernel applied to massive discs.

As was expected, our plots revealed similar behaviour among all kernels for the long-range interaction: d_g ≳ 3 when the disc is light or d_g ≳ 3 Q_g when the disc is massive. Notably, all kernels scale as $\propto d_g^{-2}$ for large distances, which respects the Newtonian behaviour. Conversely, all these prescriptions reveal distinct behaviours at short range, which we discuss next. First, we highlight that for a vanishing SL the kernel adopts an inverse square law at small distances, which highly overestimates SG. For the constant and non-vanishing SL, two scenarios arise. When the SL approaches zero (e.g. ε/H_g=0.3), the SG vanishes for small distances but is overestimated for d_g ∈ [0.2, 3]. Conversely, larger SL (e.g. ε/H_g=1.2) necessarily underestimate SG for d_g ∈ [0.2, 3]. This last limitation was corrected by Rendon Restrepo & Barge (2023) prescription, as their normalised kernel does not indeed vanish at the singularity, aligning more closely with expectations for gravitational interactions: the closer two bodies are, the stronger their mutual gravitational pull, a fundamental principle that must remain true even in the 2D approximation. However, this latter kernel is insensitive to the Toomre’s parameter of the disc, which is inconsistent. The Bessel kernel proposed in this work addresses this issue, as its value scales inversely with the Toomre’s parameter (bottom panel of Fig. 1), matching the expectation that SG strengthens in massive discs. Furthermore, it does not vanish at the singularity, a feature consistent with the nature of gravity. Finally, the Bessel kernel permits a smooth transition between light (Q_g ≳ 5) and massive discs (Q_g ≲ 1). These are expected features from a coherent 2D SG kernel. We also note that with the Bessel kernel, the effects of SG start to become appreciable when Q_g ≲ 5. It is noticeable that for the short and long range interactions, the Bessel kernel scales as $\propto d_g^{-1}$ and $\propto d_g^{-2}$, respectively, in agreement with an apparent 2D-3D gravity that we discuss next. Finally, we highlight that the point-wise symmetry of our Bessel kernel, enabled by the root mean square scale height definition (Eq. (9)), eliminates the need for correction terms arising from radial self-accelerations (Baruteau 2008).

3.2 The gravitational character of our formalism

It is worth mentioning a property that we find quite intriguing: while the long-range interaction of our kernel mirrors that of a gravitational force in a 3D space, at short range, it exhibits characteristics akin to a purely 2D nature since the force scales as ∝ 1/s (see the introduction). This feature is particularly important for the short-range interaction of SG within the context of PPDs – and it insightfully highlights the inherent discrepancies of the Plummer potential paradigm.

On the one hand, as the SL approaches the disc scale height, ε ∼ H_g, it causes a significant underestimation of the short-range effects of SG, potentially hindering the formation of gravitationally bound objects. However, on the other hand, if the SL vanishes, it is implied that the scaling of the 2D force is 1/s², which in turn results in a substantial overestimation of SG. Specifically, in the context of GI, we anticipate that the use of a vanishing SL could result in an overestimation of the formation of bound objects. This may offer a new perspective on the results reported by Vorobyov & Elbakyan (2018), among others. The realism of GI simulations employing a finite yet small SL, in works like Paardekooper et al. (2011) and Takahashi et al. (2016), appears more challenging to predict. While the overestimation of SG in the range d_g ∈ [0.2, 3] might suggest an amplification of GI, the suppression of SG at small scales could indicate the inability of fragments to further contract gravitationally, thereby promoting their destruction upon encounters with spirals. A more in-depth numerical investigation is required to evaluate the consequences of employing different gravity prescriptions. This will be the focus of our upcoming paper.

3.3 Gas and dust kernel comparison

An important feature of our analysis lies in the fact that our kernel does not apply solely to gas discs but it can also account for an embedded dust fluid. Such a situation results in four different contributions to the SG forces: one is attributed to the gas disc, another to the dust layer and two to the mutual interaction between the gas and dust discs. Even if the crosswise contributions are not commutative, the above four contributions only necessitate the evaluation of three kernels: K_gg, K_dd and K_gd=K_dg. As a consequence, we aim to compare how the kernels involving dust, K_dd and K_dg, compare to the kernel of a disc uniquely made of gas for different dust-to-gas scale heights ratios: $$\eta =H_d^{sg}/H_g^{sg}.$$(10) In the top panel of Fig. 2, we depict the dust-to-dust kernel, while in the bottom panel we study the dust-to-gas kernel, which accounts for the crossed contribution. Both quantities were normalised with respect to the gas-gas kernel, K_gg. Motivated by highly settled discs (Villenave et al. 2020, 2022) and also by younger systems with limited sedimentation (Lin et al. 2023), we constrained our investigation to 0.01 ≤ η ≤ 1. As expected, for both cases we observe that the kernel matches the one of a system only made of gas at long distances. Indeed, at long separations a force SG kernel does not depend on the value of the scale height, a feature that also applies within the SL prescription – see Fig. 13 of Müller et al. (2012). Most interestingly, we notice that at short distances, d_g ≲ η, the dust-dust kernel is proportional to 1/η. For instance for η=0.01, this suggests that there is an enhancement of dust SG by a factor of 100. Acknowledging that we expect column density dust-to-gas ratios of Z=Σ_d/Σ_g=0.01, this suggests that the SG acceleration due to dust is comparable, if not higher, than the one of gas in the limit of short distances. On numerical grounds, this should be however taken with a grain of salt since currently there are little chances that such distances can be resolved by global simulations, except maybe when using adaptive mesh refinement techniques or in the context of shearing box simulations. Regarding the normalised gas-dust kernel, we observe that, at short distances, it rises very quickly (with respect to decreasing η) towards a limiting value of $\sqrt{2}$. Surprisingly, this means that the kernel tends towards the planet-disc kernel (see Eq. (20) of Müller et al. 2012), a feature already pointed out by Rendon Restrepo & Barge (2023). We provide a physical explanation of this result in Sect. 4, below.

3.4 Gravitational runaway

In the following, we explore the ability of the Bessel kernel to potentially lead to a runaway process. Let us consider a disc only made of gas where a perturbation caused a local and infinitesimal mass increase over an infinitesimal distance. First, we focus on a massive disc (Q_g ≲ 1), where the local mass increase will, in turn, decrease the local Toomre’s parameter. But in the specific case of massive discs the kernel at the singularity³ is $K_{gg}\sim {\left(\sqrt{2}{Q}_g{H}_{g}s\right)}^{-1}$. Subsequently, the local SG is amplified, which feeds into the mechanism, resulting in a possible runaway process. Now, consider the same scenario for a light disc, where the Toomre’s parameter is large. In such a case, the kernel at the singularity is independent of the density perturbation: $K_{gg}\sim {\left(\sqrt{\pi }{H}_{g}s\right)}^{-1}$. Thus, this time, any local increase in density leaves the kernel unaffected, which does not amplify the mechanism and possibly shuts it down. We believe that such a runaway process, which mimics a potential collapse in the vertical direction (through the reduction of the vertical scale height), is uniquely captured by our prescription.

We found that out of all prescriptions presented in Sect. 3.1, only ours permits one to maintain a gravitational collapse at infinitesimal distances. Indeed, for Müller et al. (2012) kernel, the rationale of the above paragraph would never lead to a runaway process at infinitely small scales since this kernel vanishes at the singularity, preventing any amplification of gravity caused by the kernel. A runaway is only possible if the length of the perturbation is of the order ∼ H_g for the Müller et al. (2012) kernel but this time the runaway is solely driven by the overdensity term present in the total force calculation (Σ_g term in Eq. (6)). Similarly, for the SVSL formalism proposed by Rendon Restrepo & Barge (2023), even if the kernel does not vanish at the singularity, the runaway would be driven by the overdensity since the kernel is insensitive to the local Toomre’s parameter increase. All the aspects raised above will be addressed in more detail in the accompanying paper. In the following section, we propose to validate our kernel against analytical and numerical tests.

Fig. 2 Normalised kernels associated with dust with respect to distance for different dust-to-gas scale heights, η. The dust-dust, Kdd, and dust-gas, Kdg, kernels are insensitive to η at large distances, an aspect also in agreement with the SL prescription (see Fig. 13 of Müller et al. 2012). The dust-dust kernel scales with 1/η at short range, which suggests that the SG acceleration due to dust can be comparable of the one of gas even for Z=Σd/Σg=0.01. At short distances, the normalised dust-gas kernel converges quickly towards $\sqrt{2}$ (i.e. with respect to decreasing η). As a consequence this kernel adopts a nature akin to the kernel associated with a planet-disc interaction.

4 Analytical and numerical validation

We found it puzzling that the mathematical expression of the Bessel kernel (Eq. (8)) is very similar to the kernel of a planet-disc interaction (Müller et al. 2012, Eq. (20)) by a scaling factor of $\sqrt{2}$ in the normalised distance. Convinced that this similarity is not coincidental, we proceed to demonstrate it within our formalism. To elaborate, let us consider a disc made of two fluids: one gas component and a second ‘fluid’ component that we consider as a planet. From a purely gravitational perspective, a planet is nothing more than a fluid with a Dirac distribution in all space directions. Accordingly, the fluid planet is a point and its density in the 2D plane can be expressed as: $${\Sigma }_p\left({\mathit{r}}^{\text{'}}\right)=\frac{m_p}{\pi R_p^2}\chi \left(\frac{\parallel {\mathit{r}}^{\text{'}}-{\mathit{r}}_{c,p}\parallel }{2R_p}\right),$$(11) where χ is the rectangular function, R_p and r_{c, p} are the radius and location of the planet, respectively. As the planet’s radius is significantly smaller than all other distances in the problem, the rectangular function could be naturally approximated by a Dirac distribution. Furthermore, the planet inherently possess a zero vertical extension, which leads to the next rms scale height for the system {gas disc + planet}: $H_{gp}=H_g/{\sqrt{2}}^4$⁴ (see Eq. (9)). With all these observations, we can calculate the gravitational force density of the planet acting on the gas disc: $$\begin{array}{rl}{\mathit{f}}_{2D}^{p\to g}& =-{\mathrm{\Sigma }}_g{\iint }_{disc}{\mathrm{\Sigma }}_p{K}_{gp}{\mathit{e}}_s{d}^2{\mathit{r}}^{\mathrm{\prime }}\\ & =-{\mathrm{\Sigma }}_g\frac{m_p}{s}\sqrt{\frac{2}{\pi }}\frac{1}{H_g}{c}^2\exp \left(c^2\right)[K_1\left(c^2\right)-K_0\left(c^2\right)]{\mathit{e}}_{c,p}\end{array}$$(12) where $c^2=\frac{s^2}{4H_g^2}$ and e_{c, p}=(r−r_{c, p}) /|r−r_{c, p}| is the unit vector oriented from the planet to the fluid element located at r. The above expression is exactly the force density of a planet acting on the gas disc (Müller et al. 2012, Eq. (25)), which validates the consistency of our approach. In other words, the planet-disc interaction kernel is nothing more than a special case of the Bessel prescription emphasised in this article. This observation, which began as a curiosity, constitutes our first analytical test.

Another basic analytical test for assessing the correctness of our approach lies in the retrieve of the potentials associated with razor-thin discs, i.e. when the scale height of the disc tends to zero, $H_g^{sg}\to 0$. This is achieved simply by using the well-known result $$\delta \left(z\right)={lim}_{H_g^{sg}\to 0}\frac{1}{\sqrt{2\pi }{H}_g^{sg}}{e}^{-\frac{1}{2}{(z/H_g^{sg})}^2},$$(13) which, when used in Eq. (7), permits one to get the force kernel of a razor-thin disc: $$\begin{array}{rl}{\mathrm{K}}_{ab}& =\underset{H_g^{sg}\to 0}{lim}\frac{1}{2\pi }\frac{s}{H_g^{sg}{H}_g^{sg}}{\iint }_{z,z^{\mathrm{\prime }}=-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{e^{-\frac{1}{2}{\left(z/H_g^{sg}\right)}^2}{e}^{-\frac{1}{2}{\left(z^{\mathrm{\prime }}/H_g^{sg}\right)}^2}}{{(s^2+{(z-z^{\mathrm{\prime }})}^2)}^{3/2}}dzd{z}^{\mathrm{\prime }}\\ & =s{\iint }_{z,z^{\mathrm{\prime }}=-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\delta (z)\delta \left(z^{\mathrm{\prime }}\right)}{{(s^2+{(z-z^{\mathrm{\prime }})}^2)}^{3/2}}dzd{z}^{\mathrm{\prime }}\\ & =\frac{1}{{|\mathit{r}-{\mathit{r}}^{\mathrm{\prime }}|}^2}\end{array}$$(14) Such a limit naturally applies to the kernel under its closed form (Eq. (8)) by the use of the dominated convergence theorem, which permits one to commute the vertical integrals and the limits. This second simple analytical validation test shows that we expect to recover analytical and numerical results that apply to flat discs. Consequently, in order to benchmark our Kernel and check our numerical implementation, we decided to compare in the next sections the SG forces provided by our prescription against exact derivations for specific matter distributions and/or 3D numerical tests. We primarily focus in the next sections on the accuracy of our spectral solver and the comparison between 3D simulations and different 2D gravity prescriptions, leaving for a future article the benchmark of the multi-fluid treatment. We start by presenting our numerical set-up and codes.

4.1 Numerical method and codes

The tests of this section were performed employing both the FargoCPT (Rometsch et al. 2024) and Nirvana-III (Ziegler 2004, 2016) codes. FargoCPT is a 2D finite difference code with a staggered mesh, employing an advection scheme akin to finite volume methods. It solves the hydrodynamics equations using operator splitting and a second-order upwind scheme. FargoCPT and its variants, such as FARGOCA (Lega et al. 2014), FARGOADSG (Baruteau & Masset 2008), and FARGO3D (Benítez-Llambay & Masset 2016), are built upon the Fargo code presented in (Masset 2000). It is formally accurate up to second-order in space and first-order in time and relies on artificial viscosity for shock smoothing (Rometsch et al. 2024). In this code the Bessel kernel was computed by means of full FFT methods, which requires one to satisfy H/r= const. (See Sect. 5.2). Thus, we used a logarithmic grid in the radial direction and a linear grid in the azimuthal direction. We also made use of zero-padding for density in order to guarantee artificial periodicity in the radial direction.

The Bessel kernel was also implemented in the updated version of the Nirvana-iii code, as presented in Gressel et al. (2020). Unlike our implementation in FargoCPT, this time we used a semi-spectral method, which removes the constraint on linear scale heights, allowing for linear radial grids and greater generality (See Sect. 5.2). Finally, in Sect. 4.4, we also employed Nirvana-iii to perform 3D simulations, enabling a comparison between the outcomes of 2D and 3D simulations. In particular, Nirvana-III features a new approach (Gressel & Ziegler 2024) based on James’ method for the computation of the gravitational potential at the poloidal domain boundaries (see James 1977; Moon et al. 2019). We start by showing that the Bessel prescription permits the retrieval of the radial accelerations of flat power law discs.

4.2 Power law discs

One of the very few exact analytical solutions of the gravitational potential in flat discs concerns axisymmetric, finite-size, and power law density discs (Huré & Hersant 2007; Huré et al. 2008; Baruteau 2008). For such a disc, the authors decomposed the exact potential in terms of series, which permits an evaluation of the SG radial component⁵: $$\frac{g_r^{\left(n\right)}\left(\varpi \right)}{g_{0,H}}=A(1+{\beta }_{\rho }){\varpi }^{{\beta }_{\rho }}+{\sum }_{k=0}^n(2k{a}_k{\varpi }^{2k-1}-(2k+1)\frac{b_k}{{\varpi }^{2(k+1)}}),$$(15) where ϖ=r/r_out, β_ρ is the surface density power index and g_{0, H}=−2π G Σ(r_out). The constant A and the sequence coefficients can be found in the original papers. Naturally, the above series converges asymptotically to the exact gravitational force. In their study, the authors considered an infinitely thin disc and evaluated the potential in the equatorial plane. Conversely, our results are applicable to a Gaussian vertical profile but they can be safely compared to theirs in the limit of the razor-thin disc approximation, i.e. small aspect ratios: $H_g^{sg}/r\ll 1$.

In Fig. 3, we depicted the force derived from Huré et al. (2008) potential (solid blue line) for n=100 and compared it with the force provided by our Bessel kernel for the power law index β_ρ=−1.5 and decreasing aspect ratios. Our numerical set-up is $$\{\begin{array}{cc}{\Sigma }_0\left(r\right)\hfill & =20000{(r/1\mathrm{AU})}^{{\beta }_{\rho }}\mathrm{g}{\mathrm{cm}}^{-2}\hfill \\ H/r^q\hfill & =[0.1,0.01,0.001]\hfill \\ (r_{\text{in}},r_{\text{out}})\hfill & =(20,250)\mathrm{AU}\hfill \\ (N_r,N_{\phi })\hfill & =(1200,2800)\hfill \end{array},$$(16) where N_r and N_ϕ are the number of cells in the radial and azimuthal direction, respectively. Using FFT in the radial direction with FargoCPT required setting q=1, while for Nirvana-iii, which offers more flexibility, we used q=0. We chose these two distinct values of q to illustrate that our framework accommodates to any disc flaring. For FargoCPT simulations, the relative error decreases with a decreasing disc aspect ratio, aligning with our Bessel kernel’s ability to achieve the flat disc limit. For smaller aspect ratios, the accuracy is limited to approximately ∼ 10⁻³ and ∼ 10⁻⁴ for FargoCPT. The same trend is observed with Nirvana-iii, but the relative error reaches 10⁻⁴ more rapidly as the aspect ratio decreases. Notably, when H=0.1, the error is around 10⁻³, two orders of magnitude smaller than the equivalent H/r=0.1 set-up. This indicates that the common assumption of H/r= const. in thin disc simulations may introduce significant errors. We highlight that the exact radial SG force cancels at ϖ ≃ 0.07, which leads to a statistical bias in the relative error in this region. Finally, we checked that the accuracy trend with respect to the disc aspect ratio remains consistent for other power law indices ranging from [−3, 3], validating our first test.

Fig. 3 Radial SG force of the power law disc for different scale heights. The accuracy of our kernel increases with decreasing aspect ratios. Note that for ϖ ≃ 0.07 the relative error is beyond unity, which is a bias in the statistical measure due to a vanishing radial acceleration. We recall that in FargoCPT and Nirvana-iii we used a logarithmic and uniform radial grid, respectively.

4.3 The exponential disc

Another analytical result that is useful for comparison with our Bessel kernel relates to exponential discs: $$\Sigma \left(r\right)={\Sigma }_{0}exp(-r/r_0).$$(17) Indeed, the gravitational potential associated with these discs in the equatorial plane is known (Freeman 1970; Binney & Tremaine 2008, Sect. 2.6), which permits one to exactly derive the radial acceleration due to the disc’s own gravity: $$g_r\left(r\right)=-g_{0,F}y[I_0\left(y\right)K_0\left(y\right)-I_1\left(y\right)K_1\left(y\right)],$$(18) where y=r /(2 r₀) and g_{0, F}=2 π G Σ₀. The notations I_α and K_α refer to modified Bessel functions (of order α) of the first and second kind, respectively. Once again, our goal here is to compare the exact acceleration of the exponential discs (Eq. (18)) with the acceleration obtained numerically by the use of the Bessel prescription in the limit of small disc aspect ratios. For the numerical side, our set-up is $$\{\begin{array}{cc}{\Sigma }_0\left(r\right)\hfill & =200\mathrm{g}{\mathrm{cm}}^{-2}\hfill \\ H/r^q\hfill & =[0.1,0.01,0.001]\hfill \\ r_0\hfill & =250\mathrm{AU}\hfill \\ (r_{\text{in}},r_{\text{out}})\hfill & =(0.01,2500)\mathrm{AU}\hfill \\ (N_r,N_{\phi })\hfill & =(1200,2800)\hfill \end{array}.$$(19) We maintained the same set-up as in the previous section, with q=1 for FargoCPT and q=0 for Nirvana-iii. The results of this comparison are depicted in Fig. 4. The conclusions are similar to the ones drawn for the power law disc (see Sect. 4.2), with the difference that this time the radial SG force vanishes at y=0. The only difference is that FargoCPT more closely approximates the analytical solution than Nirvana-III for H/r^q=0.001. This is because FargoCPT's logarithmic grid better resolves the inner disc, where most of the mass lies.

Fig. 4 Radial SG force of the exponential disc for different scale heights. The accuracy of our kernel increases with decreasing aspect ratios and piles up at ∼ 10−4. Note that for y=0 the relative error is beyond unity, which is a bias in the statistical measure due to a vanishing radial acceleration.

4.4 Gaussian discs test

To assess the accuracy of 2D SG prescriptions, it is common to resort to the Gaussian spheres test (Chan et al. 2006; Li et al. 2009). This test involves averaging vertically the 3D gravitational field associated with 3D Gaussian spheres positioned at different locations and comparing it with the results obtained using a 2D kernel. A benefit of this benchmark is its capability to assess the azimuthal component of the gravitational force. Given the analytical challenges in vertically averaging the SG force of a Gaussian sphere, we chose to conduct this test solely through numerical methods, as is explained in the following paragraphs.

Our test entails computing the 3D gravitational field of Gaussian discs in z= const. planes and serving as perturbation to a background Gaussian vertically stratified density. The resulting 3D gravitational field is then vertically integrated and used for computing the reference field amplitude in the equatorial plane: $${\mathcal{G}}_{\text{ref}}^{3D}=\sqrt{{\langle \frac{1}{\sqrt{2\pi }H}{e}^{-\frac{1}{2}{(z/H)}^2}{\mathcal{G}}_r\rangle }_z^2+{\langle \frac{1}{\sqrt{2\pi }H}{e}^{-\frac{1}{2}{(z/H)}^2}{\mathcal{G}}_{\phi }\rangle }_z^2}.$$(20) We highlight that ${\mathcal{G}}_{\text{ref}}^{3D}$ denotes the 3D gravitational field obtained with Nirvana-iii, followed by vertical integration denoted by 〈…〉_z, with the appropriate factors accounted for. This reference field is then compared with the results of the 2D SG prescriptions, $${\mathcal{G}}_{2D}=\left|{\int }_{\text{disc}}{\Sigma }_g{\mathrm{K}}_{gg}{\mathit{e}}_s{d}^2{\mathit{r}}^{\text{'}}\right|,$$(21) specifically the Bessel kernel, on the one hand, and a Plummer potential with two different SL values, ε/H=[0, 1.2], respectively, on the other hand. While 𝒢_2D is not strictly the absolute value of a 2D gravitational field, as it is not derived from a potential (as detailed in Sect. 5.5), it remains the quantity that can be faithfully compared with its 3D counterpart, ${\mathcal{G}}_{\text{ref}}^{3D}$. To assess the potential limitations of each prescription based on the disc vertical extent, we conducted this test for three different scale heights: H = [0.01, 0.1, 0.4].

For this test, the density consists of a constant background distribution, $$\{\begin{array}{cc}{\rho }_0(\mathit{r},z)\hfill & =\frac{{\Sigma }_0}{\sqrt{2\pi }H}exp(-\frac{1}{2}{\left(\frac{z}{H}\right)}^2)\hfill \\ {\Sigma }_0\hfill & =200\mathrm{g}{\mathrm{cm}}^{-2}\hfill \\ H\hfill & =[0.1,0.01]\mathrm{AU}\hfill \end{array},$$(22) with perturbations introduced in the form of the aforementioned Gaussian discs, $$\rho (\mathit{r},z)={\rho }_0(\mathit{r},z){\prod }_{i=1}^{N=3}[1+A_{i}exp(-\frac{1}{2}\frac{{|\mathit{r}-{\mathit{r}}_{c,i}|}^2}{R_i^2})].$$(23) The parameters A_i, r_{c, i}, and R_i, which represent the relative amplitude, centre location, and radius of the Gaussian, respectively, are detailed in Table 2. For the evaluation of the 3D gravitational field employing the NIRVANA-iii code, we used a cylindrical grid with next computational domain extents: $$\{\begin{array}{ccc}r\hfill & \in \hfill & [2,10]\mathrm{AU}\hfill \\ \phi \hfill & \in \hfill & [0,2\pi ]\hfill \\ z\hfill & \in \hfill & [-5H,5H]\hfill \end{array}.$$(24) We used two different number of cells for the 3D set-up: $$(N_r,N_{\phi },N_z)=\{\begin{array}{cc}(256,1024,64)\hfill & \text{if}H=0.01\mathrm{AU}\hfill \\ (256,1024,128)\hfill & \text{if}H=0.4,0.1\mathrm{AU}\hfill \end{array}.$$(25) This distinct choice permitted us to reduce the anisotropy in our 3D grid, which facilitated the convergence of the multigrid solver. For the 2D simulations counterpart, we kept the radial and azimuth resolution unchanged compared to the 3D set-up.

In Fig. 5 we depicted the relative difference between the reference simulation and it is 2D analogues for different SG prescriptions, $$\epsilon =|{\mathcal{G}}_{2D}-{\mathcal{G}}_{3D}^{ref}|/max\left({\mathcal{G}}_{3D}^{ref}\right),$$(26) associated with the above set-up and for H=0.4. The amplitude of the 2D gravitational field, computed using the semi-spectral method implemented in Nirvana-iii, is depicted by 𝒢_2D. The choice of the above error stems from the fact that the gravitational field amplitude approaches 0 in two regions, which makes the standard relative norm diverge, which bias the interpretation. We quantified the deviation by presenting the following L_∞ norm: $${\parallel \epsilon \parallel }_{\infty }=max\left(\epsilon \right).$$(27) The above L_∞ norm is presented in Table 3 for various 2D SG prescriptions. Our findings indicate that for the three different scale heights, the error associated with the Bessel kernel remains within a few percent. Although the maximum error is found near the Gaussian discs positions, it is not centred on them. These minimal deviations, even for high scale heights such as H=0.4, are due to our method’s versatility, which is not limited to razor-thin discs but is applicable to any geometrical aspect ratio. Conversely, the error associated with the Plummer potential varies with the disc scale height and decreases as the scale height diminishes. Specifically, for ε/H=1.2 and H=0.4, the maximum error reaches 28% and is located in the disc of the smallest radius. Indeed, this prescription fails to resolve gravitationally bumps with dimensions smaller than the SL. Notably, this set-up generally underestimates SG. On the other hand, for ε/H=0 and H=0.4, the maximum error reaches 129%. In this case the highest errors are located in the regions of highest density, and this prescription generally overestimates SG. When the scale height decreases to H=0.01, all three prescriptions yield similar results, with an overall maximum error of just a few percent.

With this test, we have demonstrated that the Plummer potential prescription fails to provide an accurate estimation of SG in regions of high density, where SG is overestimated for ε/H=0, and in regions smaller than the scale height for ε/H=1.2, where SG is underestimated. Conversely, the Bessel prescription consistently delivers accurate results across the entire range of scale heights.

5 Discussion

In this section, we explore the potential consequences of our approach for tackling GI, as well as the numerical methods of estimating SG, and we consider alternative methods that offer greater generality. Then, we identify the limitations of our approach and present solutions for scenarios where we deviate from the initial assumptions. Finally, we delve into a detailed discussion on the Poisson’s equation for 2D systems and its implicit implications for 3D systems. In this context, we propose an effective 2D gravitational potential consistent with our formalism.

5.1 Consequences for planet formation

Our detailed study of the Bessel kernel, along with the tests conducted in Sects. 3–4, has revealed and quantified that the Plummer potential can both under- and overestimate the importance of SG. These findings could have profound implications for planet formation theories, particularly the GI, as most simulations characterising it have been conducted using a thin disc approximation. This approach typically employs a SL prescription for the gravitational potential or involves solving a 2D Poisson equation. We take this opportunity to clarify that solving a 2D Poisson equation, as is often done in the flat disc limit, is equivalent to using a Plummer potential with ε/H=0, as is explained in Sect. 5.5.

On the one hand, we affirm that a finite SL suppresses the Newtonian character of gravity, potentially quenching gravitational collapse or enhancing clump disruption during spiral encounters (Young & Clarke 2015). On the other hand, we also affirm that a vanishing SL artificially amplifies gravity. In this context, we recall that an ongoing issue is the convergence problem, which stems from the lack of consensus about the threshold cooling parameter that separates the regime of gravito-turbulence from fragmentation in 2D simulations (Young & Clarke 2015; Kratter & Lodato 2016). The issue first emerged in 3D SPH simulations of gravito-turbulence (Meru & Bate 2011), where resolution-dependent artificial viscosity likely played a significant role (Lodato & Clarke 2011). In 2D grid-based simulations, some aspects of the problem were later attributed to inner edge effects (Paardekooper et al. 2011), though the issue remains only partially addressed in these set-ups. Consistent with these earlier findings, Paardekooper (2012) demonstrated that fragmentation is a stochastic process, allowing discs to fragment at broad beta-cooling ranges and blurring the boundary between gravito-turbulence and fragmentation regimes. Given that these simulations were conducted thanks to 2D shearing box simulations with a similar framework where SL equals zero, we think that our Bessel formulation would be key to addressing the convergence problem. More specifically it will likely narrow and better constrain the beta cooling parameter range where stochastic fragmentation occurs.

Fig. 5 Relative difference between the 2D SG prescriptions and the reference gravitational field (see Eq. (26)) for H = 0.4 AU. From left to right: (left panel) the sharp potential with ε/H = 0, (middle panel) the variable Bessel kernel, and (right panel) the soft potential with ε/H=1.2. The colours of the contour levels at 10−5, 10−3, 0.1, 0.5 are indicated in the colour bar. Softening prescriptions exhibit errors exceeding tens of percent at disc positions, whereas the Bessel prescription maintain errors below 3%.

5.2 Numerical assessment with FFT methods

The calculation of SG using the Bessel kernel remains compatible with FFT techniques. This compatibility arises from expressing the radial and azimuthal components of the SG force as convolution products (Baruteau & Masset 2008; Rendon Restrepo & Barge 2023). The explicit forms of the kernel components, ready for use in convolution operations, are provided in Appendix C. The main advantage of this method lies in the acceleration of the calculation since it requires N_r N_ϕ log (N_r N_ϕ) operations, in contrast with direct summation that requires (N_r N_ϕ)² operations. This gain in performance requires one nonetheless to obey some constraints. First, in order to enable a description in the form of a convolution product, it is necessary that the radial grid is logarithmic. This requires that the disc aspect ratio of the self-gravitating disc, $H_{ab}^{sg}/r$, remains constant in the whole numerical domain (see Appendix D for a discussion when deviating from this condition). Following the procedure of Rendon Restrepo & Barge (2023, Sect. 4.2), we found that this indirectly implies that the bi-fluid Toomre’s parameter, Q_bi-fluid, and the standard disc scale height, H_a/r, should remain constant throughout the entire numerical domain – and also in time. This stringent requirement precludes the use of FFT methods in realistic set-ups. Consequently, if FFT methods are necessary, the simplest approach is to use the thermal scale height without the SG correction: H_sg ≃ H. Second, the use of Fourier transforms demands periodicity in all directions, which requires one to double the radial domain with a zero-padding in order to perform an aperiodic convolution (Hockney & Eastwood 2021). This operation does not affect the computational endeavour since the Hermitian property of Fourier transforms halves the total number of transforms and, as a consequence, compensates the doubling of the domain. Finally, this zero-padding does not affect the SG computation as checked in Sect. 4.

Computing the Bessel kernel requires one to use special functions, which is computationally expensive. In FargoCPT, the Bessel kernel calculation is roughly 20 times slower than the original kernel calculation. If the Bessel kernel is recomputed at every iteration, this increases the cost of SG from around 40% to around 130% of the transport step cost, or from 21% to 45% of the total computational costs when using a locally isothermal equation of state. For the latter, however, the SG kernel can be precomputed for the whole simulation because it does not change with time. Then, there is no difference in computational cost after the initialisation. In case of a non-isothermal equation of state, the kernel can be recomputed only every so often to save computational costs. This can be justified by the fact that the scale height of the disc only changes slowly, at least when the cooling timescale is short compared to the dynamical timescale. For a more realistic stratification in presence of SG, there is however also a physical issue. In this case, the aspectratio of the disc is no longer constant and the method is formally no longer applicable and only yields an approximation to the SG accelerations.

In general, we expect the disc flaring to not be constant, which renders the application of our method with FFT impractical. For instance, the $H_{ab}^{sg}/r$= const. constraint (necessary for full spectral methods) becomes prohibitive when investigating massive discs. Indeed this constraint implies that the Toomre’s parameter needs to be constant in the whole domain and frozen in time. Therefore, a potential decrease in the Toomre’s parameter due to a gravitational runaway process will not be captured by these methods. As a consequence, the other computational possibility, which was in fact implemented in Nirvana-iii, consists of resorting to a pseudo-spectral method based on a combination of FFT in the azimuthal direction followed by a direct summation in radial direction, whose complexity is in $N_r{N}_{\phi }log\left(N_{\phi }\right)+N_{\phi }{N}_r^2$ (Li et al. 2009). What is new in this implementation is that every few time steps the scale height is updated according to Eq. (3). This method is a good trade-off between computational efficiency and generality, since it is adapted to any scale height prescription beyond problems governed only by SG. Another possible solution that eliminates the H_{a b}(r)/r= const. limitation while preserving the computational efficiency, with N log N operations, involves using a fast Hankel transform algorithm in the radial direction.

We now address the numerical treatment of the multi-fluid case. In this scenario, it is necessary to pre-compute the self-kernels for each fluid, as well as the cross-gravitational terms between different fluid species. This can be quite burdensome. However, since the Kernels ultimately depend only on the root mean square scale height, we propose that some of the involved fluids are assumed to have the same scale heights. This significantly reduces the number of required Kernels. For instance, consider a situation in which two dust species have the same scale height. In this case, only four Kernels and six convolutions would need to be computed at each time step:

Volume forces acting on the gas: $${\Sigma }_g{\mathcal{F}}^{-1}[{\widehat{\Sigma }}_g*{\widehat{\mathrm{K}}}_{gg}]\text{and}{\Sigma }_g{\mathcal{F}}^{-1}[({\widehat{\Sigma }}_{d1}+{\widehat{\Sigma }}_{d2})*{\widehat{\mathrm{K}}}_{dg}]$$(28)

Forces acting on the dust 1: $${\Sigma }_{d1}{\mathcal{F}}^{-1}[({\widehat{\Sigma }}_{d1}+{\widehat{\Sigma }}_{d2})*{\widehat{\mathrm{K}}}_{dd}]\text{and}{\Sigma }_{d1}{\mathcal{F}}^{-1}[{\widehat{\Sigma }}_g*{\widehat{\mathrm{K}}}_{gd}]$$(29)

Forces acting on the dust 2: $${\Sigma }_{d2}{\mathcal{F}}^{-1}[({\widehat{\Sigma }}_{d1}+{\widehat{\Sigma }}_{d2})*{\widehat{\mathrm{K}}}_{dd}]\text{and}{\Sigma }_{d2}{\mathcal{F}}^{-1}[{\widehat{\Sigma }}_g*{\widehat{\mathrm{K}}}_{gd}].$$(30)

Here, ^∩ denotes the Fourier transform, ℱ⁻¹ the inverse Fourier transform, and * the convolution product. Accuracy of the approximation depends on the chosen scale height for calculating the dust kernels. A wise choice would be to take the harmonic mean of the involved fluids. In the above example, this would be: $H_d={(\frac{1}{H_{d1}}+\frac{1}{H_{d2}})}^{-1}$.

5.3 Limitations

The use of the Bessel kernel proposed in this paper requires that the stratification of the disc be strictly Gaussian. This is of course an assumption that depends on the physical ingredients incorporated in the simulated disc and can be altered; for example, by considering the stellar irradiation (Wu & Lithwick 2021) or in the presence of magnetic fields (e.g. Johansen & Levin 2008; Gaburov et al. 2012). In the former case, the system is constantly evolving in the vertical direction which suppresses a potential vertical hydrostatic equilibrium. However the Bessel kernel does not require hydrostatic equilibrium. It only requires a Gaussian stratification, which is possible even when the equilibrium is not ensured in the vertical direction. Once again all the time evolution information should be incorporated in the scale height. For the latter, a simple rectification would be to replace the sound speed of gas by c_g(1+1/β)^1/2, where β is the plasma beta parameter (Johansen & Levin 2008). This would only modify the standard scale height H_g used in Eq. (3). In the case of more sophisticated stratifications that significantly deviate from the Gaussian one, a reasonable approach to maintain the applicability of our prescription would be to decompose the vertical profile of the disc as a finite linear combination of Gaussian functions. It is important to note that this method will result in an increasing number of numerical evaluations by FFT methods. If such a case occurs, it is surely wiser to conduct 3D simulations directly.

During our numerical tests, we observed the generation of NaN values in the SG forces when dealing with small disc aspect ratios. This occurrence was traced back to the overflow of the exponential function employed in the Bessel kernel definition. This issue is simply addressed by using the Taylor expansion of the kernel at infinity: $${\mathrm{K}}_{ab}\left(X\right)=\frac{1}{2\sqrt{2}}\frac{H_{ab}^{-2}}{d_{ab}}[\frac{1}{X^{1/2}}-\frac{3}{8}\frac{1}{X^{3/2}}+\frac{45}{128}\frac{1}{X^{5/2}}+O\left(\frac{1}{X^{7/2}}\right)],$$(31) where $X=d_{ab}^2/8$. We applied the above approximation when the variable X was larger than 50, successfully resolving our issue without any impact on the accuracy of the method.

5.4 Indirect term

The presence of asymmetries in the disc can lead to an offset between the barycenter of mass and the reference frame, located at the star’s position. This offset is accounted through an indirect term: Φ_ind(r)=−a_* · r, where a_* is the acceleration of the central object (Zhu & Baruteau 2016; Regály & Vorobyov 2017; Rendon Restrepo & Barge 2022). In order to be compatible with this investigation, we found relevant to derive the indirect term associated with the pull exerted by fluid ‘a’ on the star using our approach: $${\Phi }_{a,\text{ind}}\left(\mathit{r}\right)=\sqrt{2}r{\int \int }_{disc}{\Sigma }_a\left({\mathit{r}}^{\text{'}}\right)cos(\phi -{\phi }^{\text{'}}){\mathrm{K}}_a\left(d_a\right)d^2{\mathit{r}}^{\text{'}},$$(32) where $${\mathrm{K}}_a\left(d_a\right)=\frac{1}{\sqrt{\pi }}{\left(H_a^{sg}\right)}^{-2}\frac{d_a}{4}exp\left(\frac{d_a^2}{4}\right)[K_1\left(\frac{d_a^2}{4}\right)-K_0\left(\frac{d_a^2}{4}\right)]$$ and $d_a=\parallel {\mathit{r}}^{\text{'}}\parallel /H_a^{sg}$. In Appendix E we propose a detailed derivation where you can also find the gradient of this potential. This term is particularly important since it was recently pointed out that the indirect term and the SG forces are inseparable ingredients that should be accounted, or disregarded, simultaneously in simulations (Crida et al. 2025).

5.5 Flat disc Poisson equation and gravitational potential

For numerical simulations addressing SG in thin discs, it is common to solve next Poisson’s equation (Gammie 2001), $$\Delta \Phi =4\pi G\Sigma \delta \left(z\right),$$(33) or to assess an integral (cf. Paardekooper 2012), $$\Phi (x,y,z=0)=-G\int \frac{\Sigma \delta \left(z^{\text{'}}\right)}{\sqrt{{(x-x^{\text{'}})}^2+{(y-y^{\text{'}})}^2+z^{\text{'}2}}}d{x}^{\text{'}}d{y}^{\text{'}}d{z}^{\text{'}}$$(34) in order to evaluate the gravitational potential, Φ. It is worth remembering that solving the above equations places us automatically in a paradigm where the disc is considered as infinitely thin, suggesting that the disc is cold and/or the disc SG dominate in the vertical direction. Therefore, when the above equation is used in 2D simulations, there is a possibility of overestimating the impact of SG, potentially resulting in inconsistent set-ups. For instance, in the context of GI, Toomre’s parameters are often in the range of 0.5−1. This may pose a contradiction with the razor-thin disc approximation implied by Eqs. (33)–(34), which rather implicitly suggests Q ⋘ 1. Further, there is several evidence that PPDs are thin but far from being flat, even in the dust phase, as supported by observations (Lin et al. 2023). These two statements encourage us to think that the method presented in this paper is appropriate. Indeed, in contrast to Eqs. (33)–(34), which are limited to razor-thin discs, the Bessel prescription is applicable to any Gaussian stratification of thin discs, ranging from razor-thin to thick discs, in alignment with observations.

In this context, it is legitimate to verify the possibility of existence of a Poisson’s equation with our formalism. Therefore, it is useful to rewrite the equivalent gravitational field in 2D: $$\begin{array}{rl}{\mathit{G}}_{2D}^{a\to b}& =-{\int }_{z=-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-\frac{1}{2}{\left(z/H_b^{sq}\right)}^2}}{\sqrt{2\pi }{H}_b^{sg}}{\mathbf{\nabla }}_{3D}{\mathrm{\Phi }}_a(\mathit{r},z)dz\\ & =-{\int }_{z=-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-\frac{1}{2}{\left(z/H_b^{sq}\right)}^2}}{\sqrt{2\pi }{H}_b^{sg}}{\mathbf{\nabla }}_{2D}{\mathrm{\Phi }}_a(\mathit{r},z)dz\end{array}$$(35) where we used the even symmetry of the volume density with respect to the z=0 plane. In the above expression, Φ_a is the potential generated by the mass distribution of the unique fluid ‘a’. Therefore, the total 2D SG force per unit volume exerted by a disc made of two fluids ‘a’ and ‘b’ on a volume element of fluid ‘b‘ is $${\mathit{f}}_{2D}^{tot\to b}\left(\mathit{r}\right)={\Sigma }_b({\mathit{G}}_{2D}^{a\to b}+{\mathit{G}}_{2D}^{b\to b})={\Sigma }_b{\mathit{G}}_{2D}^{tot,b}.$$(36) This expression is another way of writing Eqs. (4)–(5). In general the 2D Nabla operator and the integral cannot commute since the scale heights $H_a^{sg}$ and $H_b^{sg}$ are r-dependent functions. As a consequence, it is immediate that in our framework the SG force is not conservative since $\mathit{\nabla }\times {\mathcal{G}}_{2\mathit{D}}^{\text{tot}\to b}\ne \mathbf{0}$. This is not immediately dramatic – or in contradiction with the underlying physics – as there is no theoretical basis suggesting that the weighted average of a conservative force must retain its conservative nature.

It would be tempting to consider that a force derived from a Plummer potential would better fit the physics since it is conservative. However, we remind that prior computations of SG with a softening prescription do not respect this condition since the gravitational acceleration is not a gradient of the Plummer potential. Specifically, this is due to the absence of a ∇ ε², where ε is the SL, term in the definition of the acceleration. This issue, already highlighted by Rendon Restrepo & Barge (2023, Sect. 2.2), encourages to think that it is not necessary to use conservative 2D force for fitting the vertical integral of a conservative 3D force.

Nonetheless, there is one scenario in which the vertically averaged SG force is conservative, and that is when $H_a^{sg}$ and $H_b^{sg}$ are spatial constants. In such a case, the gravitational field can be expressed in a conservative form: $${\mathit{G}}_{2D}^{a\to b}=-{\mathit{\nabla }}_{2D}{\Phi }_{2D}^{ab}\left(\mathit{r}\right),$$(37) where the combined potential can be written in a closed form, $${\Phi }_{2D}^{ab}\left(\mathit{r}\right)=-{\int \int }_{disc}\frac{{\Sigma }_a\left({\mathit{r}}^{\text{'}}\right)}{2\sqrt{\pi }{H}_{ab}}exp\left(\frac{d_{ab}^2}{8}\right)K_0\left(\frac{d_{ab}^2}{8}\right)d^2{\mathit{r}}^{\text{'}}.$$(38) The detailed derivation can be found in appendix F. We checked that the gradient of the above potential permits one to retrieve the SG force obtained in Eq. (6), which is consistent. The right hand side integrand of Eq. (38) can be interpreted as a Green’s function, which scales as a logarithm or inverse law at the singularity or infinity, respectively. This outcome is in agreement with the observations pointed out in Sect. 3.2, which suggests that this Bessel potential permits a smooth transition for gravity between 2D and 3D. We tested several combinations of differential linear operators to the potential defined by Eq. (38), in order to find a Poisson’s like equation. This task was unsuccessful but we still think that it is possible provided that more sophisticated mathematics tools are explored, like elliptic operators. This is beyond the scope of our investigation. While the constraint of a constant scale height is unrealistic (Paneque-Carreño et al. 2023) and limits the practical use of Eq. (38) in global simulations, it could still be valuable for 2D shearing box simulations. Especially if its Fourier transform is determined. Finally, we want to highlight that in practice, simulations involving two distinct fluids with different scale heights, H_g ≠ H_d, require one to estimate four potentials: ${\Phi }_{2D}^{gg},{\Phi }_{2D}^{dd},{\Phi }_{2D}^{gd}$, and ${\Phi }_{2D}^{dg}$.

6 Conclusion

In this work, we propose an exact SG kernel ready to be used in 2D numerical simulations of PPDs made of gas and dust. The closed form of this Bessel kernel results from a direct vertical double integration and its expression is $${\mathrm{K}}_{ab}=\frac{1}{\sqrt{\pi }}{\left(H_{ab}^{sg}\right)}^{-2}\frac{d_{ab}}{8}exp\left(\frac{d_{ab}^2}{8}\right)[K_1\left(\frac{d_{ab}^2}{8}\right)-K_0\left(\frac{d_{ab}^2}{8}\right)].$$(39) This prescription applies when the vertical stratification of both fluids is Gaussian and all information regarding the vertical hydrostatic equilibrium due to the star’s gravity and the disentangled gravity interaction of both fluids is encapsulated in the definition of the different scale heights, $H_{ab}^{sg},H_a^{sg}$, and $H_b^{sg}$ (Eqs. (3)–(9)). In contrast with other kernels, the dependence of these scale heights with respect to the generalised Toomre’s parameter, Q_bi-fluid (Eq. (3)), allows for a continuous transition from a regime of weak SG to a regime dominated by SG. Additionally, this kernel captures correctly the short-range interaction of SG, permitting a possible gravitational runaway and faithful 2D simulations incorporating SG. Furthermore, this kernel is intrinsically symmetric with respect to r−r^′, which adheres to Newton’s third law. To our knowledge, this is the only 2D kernel in the literature of PPDs that simultaneously respects the above three conditions.

On numerical grounds, the Bessel kernel is compatible with FFT methods, provided that $H_{ab}^{sg}/r=$ const. holds across the computational domain – a requirement our Nirvana-iii implementation relaxes through a semi-spectral method. We tested the analytical consistency and numerical implementation of this kernel by retrieving radial accelerations from the literature of flat discs; namely, the power law and exponential discs. The relative errors are of the order of 10⁻⁴−10⁻³, which also validates the accuracy of our FFT method despite the required zero-padding. Furthermore, we compared the outcomes of 3D simulations with those of 2D simulations. This comparison highlighted the accuracy of our Bessel prescription and demonstrated that the softening prescription either underestimates or overestimates SG, with errors reaching up to 129%.

In agreement with earlier findings, our results indicate that for short distances the Plummer potential can lead to a significant underestimation of SG. Conversely, a Plummer potential with SLs approaching zero, or equivalently solving a 2D Poisson equation, overestimates SG effects. Both of these behaviours stem from the nature of SG within the Plummer potential framework, where the SL defines the scale below which gravity is screened. In its absence, gravity retains its purely 3D form (scaling as 1/r²), which is not always physically appropriate in a 2D context. In the Bessel prescription, a characteristic length also appears, but it represents the distance at which the gravitational force transitions – while retaining its Newtonian character – from a 3D scaling to a 2D scaling (1/r). As a consequence, we assert that the Bessel kernel proposed in this work should be used in order to conduct fully consistent 2D simulations including SG. The SL approach, or 2D Poisson equation, widely used in PPDs studies may significantly misestimate SG, necessitating a re-evaluation of these studies. We provide evidence for this statement in the accompanying article, in which we investigate the implications of our kernel to the GI paradigm of planet formation thanks to 2D global simulations.

Acknowledgements

We thank the referee, Clément Baruteau, for his insightful and constructive comments. Funded by the European Union (ERC, Epoch-ofTaurus, 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. T.R. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) research group FOR 2634 (grants KL 650/29-2 and 30-2) “Planet Formation Witnesses and Probes: Transition Disks”.

References

Appendix A An useful function

Let us define the function sequence: $$F_n\left(x\right)={\int }_{-\infty }^{+\infty }\frac{e^{-\frac{{\xi }^2}{2}}}{{(x+{\xi }^2)}^{n+\frac{1}{2}}}d\xi$$(A.1) for n ∈ ℕ and $x\in {\mathbb{R}}_{+}^{*}$. This function satisfies the following recursive relation: $$\begin{array}{rl}{F}_n(x)& =-\frac{2}{2n-1}{F}_{n-1}^{(1)}(x)\\ & =(-\frac{2}{2n-1})(-\frac{2}{2n-3})\dots (-\frac{2}{3})(-\frac{2}{1})F_0^{(n)}(x)\\ & =(-1{)}^n\frac{2^{2n}}{(2n)!}n!F_0^{(n)}(x)\end{array}$$(A.2) But: $$\begin{array}{rl}{F}_0(x)& ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{e^{-\frac{{\xi }^2}{2}}}{{(x+{\xi }^2)}^{\frac{1}{2}}}d\xi \\ & ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{e}^{-\frac{x}{2}{\sinh }^{2}t}dt\text{with}\xi =\sqrt{x}\sinh t\\ & =e^{\frac{x}{4}}{K}_0\left(\frac{x}{4}\right)\end{array}$$(A.3) where K₀ is a modified Bessel function of the second kind and order 0. Combining Eqs A.2 and A.3, we get: $$\begin{array}{rl}{F}_1(x)& =(-1{)}^1\frac{2^2}{2!}1!\frac{d{F}_0(x)}{dx}\\ & =\frac{1}{2}{e}^{\frac{x}{4}}[K_1\left(\frac{x}{4}\right)-K_0\left(\frac{x}{4}\right)]\end{array}$$(A.4) where K₁ is a modified Bessel function of the second kind and order 1.

Appendix B Analytic expression of the self-gravity kernel

In order to find a close form for the SG kernel, it is useful to compute the following quantity: $$I_{ab}={\int \int }_{z,z^{\text{'}}=-\infty }^{+\infty }\frac{e^{-\frac{1}{2}{z}^2/H_b^2}{e}^{-\frac{1}{2}{z}^{\text{'}2}/H_a^2}}{{(s^2+{(z-z^{\text{'}})}^2)}^{3/2}}dzd{z}^{\text{'}}$$(B.1) In this whole appendix Sect. we skip the “sg” superscripts for clarity. Using a variable substitution we get $$I_{ab}=H_a{H}_b{\int \int }_{u,v=-\infty }^{+\infty }\frac{e^{-\frac{1}{2}{u}^2}{e}^{-\frac{1}{2}{v}^2}}{{(s^2+{(H_{b}u-H_{a}v)}^2)}^{3/2}}dudv.$$(B.2) The above integral can be evaluated using the following variable substitution (Li et al. 2009), $$\xi =\frac{H_{b}u-H_{a}v}{\sqrt{H_a^2+H_b^2}},{\xi }^{\text{'}}=\frac{H_{a}u+H_{b}v}{\sqrt{H_a^2+H_b^2}}$$(B.3) which permits one to get $$\begin{array}{rl}{I}_{ab}& =\frac{H_a{H}_b}{{(H_a^2+H_b^2)}^{\frac{3}{2}}}{\iint }_{\xi ,{\xi }^{\mathrm{\prime }}=-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{e^{-\frac{{\xi }^2}{2}}{e}^{-\frac{{\xi }^{\mathrm{\prime }2}}{2}}}{{(\frac{d_{ab}^2}{2}+{\xi }^2)}^{3/2}}d\xi d{\xi }^{\mathrm{\prime }}\\ & =\sqrt{2\pi }\frac{H_a{H}_b}{{(H_a^2+H_b^2)}^{\frac{3}{2}}}{\int }_{\xi =-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{e^{-\frac{{\xi }^2}{2}}}{{(\frac{d_{ab}^2}{2}+{\xi }^2)}^{3/2}}d\xi \\ & =\frac{\sqrt{\pi }}{2}\frac{H_a{H}_b}{H_{ab}^3}{F}_1\left(\frac{d_{ab}^2}{2}\right)\end{array}$$(B.4) where the function F₁ is defined by Eq. A.1 and its closed form given by Eq. A.4. Finally, combining the above result with Eq. 7, we get $${\mathrm{K}}_{ab}=\frac{1}{\sqrt{\pi }}{\left(H_{ab}^{sq}\right)}^{-2}\frac{d_{ab}}{8}exp\left(\frac{d_{ab}^2}{8}\right)[K_1\left(\frac{d_{ab}^2}{8}\right)-K_0\left(\frac{d_{ab}^2}{8}\right)],$$(B.5) which is the SG force kernel of a thin Gaussian-stratified disc.

Appendix C Formulation as a convolution product

The formulation of the SG forces as a convolution product relies on rewriting the Bessel Kernel, and ultimately the normalised distance $d_{ab}^2$, in terms of shifted arguments. We start with the expression of $d_{ab}^2$: $$d_{ab}^2=2\frac{r^2+r^{\text{'}2}-2r{r}^{\text{'}}cos(\phi -{\phi }^{\text{'}})}{H_a^2(r^{\text{'}},{\phi }^{\text{'}})+H_b^2(r,\phi )}$$(C.1) We can rewrite it as $$d_{ab}^2=2\frac{r/r^{\text{'}}+r^{\text{'}}/r-2cos(\phi -{\phi }^{\text{'}})}{H_a^2(r^{\text{'}},{\phi }^{\text{'}})/\left(r{r}^{\text{'}}\right)+H_b^2(r,\phi )/\left(r{r}^{\text{'}}\right)}$$(C.2) By substituting r=e^u and r^′=e^u′ in the numerator, we obtain the correct form for the convolution. However, the denominator requires additional constraints: $H_a^2(r^{\text{'}},{\phi }^{\text{'}})/\left(r{r}^{\text{'}}\right)$ and $H_b^2(r,\phi )/\left(r{r}^{\text{'}}\right)$ must be solely functions of r^′/r, which occurs only if the scale heights are linear with the polar radius. Under this constraint, the normalised distance is expressed as $$\begin{array}{rl}{d}_{ab}^2& =2\frac{e^{u-u^{\mathrm{\prime }}}+e^{-(u-u^{\mathrm{\prime }})}-2\cos (\phi -{\phi }^{\mathrm{\prime }})}{h_a^2{e}^{u^{\mathrm{\prime }}-u}+h_b^2{e}^{u-u^{\mathrm{\prime }}}}\\ & =4\frac{\cosh (u-u^{\mathrm{\prime }})-\cos (\phi -{\phi }^{\mathrm{\prime }})}{h_a^2{e}^{u^{\mathrm{\prime }}-u}+h_b^2{e}^{u-u^{\mathrm{\prime }}}}\end{array}$$(C.3) where h_a=H_a(r^′)/r^′= const. and h_b=H_b(r)/r= const.

The radial and azimuthal components of the Bessel kernel, expressed in a form suitable for convolution on a logarithmic radial grid, are given by $$\{\begin{array}{c}{\mathrm{K}}_r(X,\phi )=\frac{L_{ab}\left(d_{ab}\right)}{\pi d_{ab}^2}{\left[\frac{e^{-X}}{h_{ab}^2}\right]}^{3/2}[e^X-cos\left(\phi \right)]\hfill \\ {\mathrm{K}}_{\phi }(X,\phi )=\frac{L_{ab}\left(d_{ab}\right)}{\pi d_{ab}^2}{\left[\frac{e^{-X}}{h_{ab}^2}\right]}^{3/2}sin\left(\phi \right)\hfill \end{array}$$(C.4) where X=u−u^′. The square of the geometric mean, the square of the normalised distance and the auxiliary function L_ab are given by $$\begin{array}{rl}{h}_{ab}^2& =H_{ab}^2/\left(r{r}^{\mathrm{\prime }}\right)=\frac{h_a^2{e}^{-X}+h_b^2{e}^X}{2}\\ d_{ab}^2& =2\frac{\cosh X-\cos \phi }{h_{ab}^2}\\ L_{ab}\left(d_{ab}\right)& =\sqrt{\pi }\frac{d_{ab}^2}{8}\exp \left(\frac{d_{ab}^2}{8}\right)[K_1\left(\frac{d_{ab}^2}{8}\right)-K_0\left(\frac{d_{ab}^2}{8}\right)]\end{array}$$(C.5) Note that the expression for L_ab differs slightly from the standard form of the Bessel kernel presented in Eq. 8.

Appendix D Limitations of the full spectral method for non-constant disc aspect ratios

The high performance of the full spectral method, makes necessary to employ it despite not fully satisfying its application hypothesis, namely the constant h=H/r condition. To quantify the error arising from deviating from this condition, we compared a constant flaring disc set-up to reference set-ups featuring radial power law flaring and local temperature gaps. For the reference set-ups, we used the semi-spectral method implemented in Nirvana-iii, as detailed in Sect. 5.2. Our numerical set-up is $$\{\begin{array}{cc}{\Sigma }_0\left(r\right)\hfill & =20000{(r/1\mathrm{AU})}^{-1.5}\mathrm{g}.{\mathrm{cm}}^{-2}\hfill \\ r_0\hfill & =100\mathrm{AU}\hfill \\ (r_{\text{in}},r_{\text{out}})\hfill & =(20,250)\mathrm{AU}\hfill \\ (N_r,N_{\phi })\hfill & =(1200,2800)\hfill \end{array}.$$(D.1) To investigate the impact of flaring, we chose $$H/r=0.05{(r/1\mathrm{AU})}^p,$$(D.2) with p ∈ [0.125, 0.25, 0.50]. For the second set-up, testing the impact of a temperature gap, we used $$H/r=0.05\sqrt{1-0.95exp(-{\left(\frac{r-100\mathrm{AU}}{w}\right)}^2)},$$(D.3) with w ∈ [10, 40] AU.

Our results, presented in Figure D.1, demonstrate that as the h= const. condition deviates, the error increases across the entire disc. The error ranged from approximately a few percent for p = 0.125, to a few tens of percent for p = [0.25, 0.5]. In contrast, the impact of the gap was localised and minimal, with the worst-case error being of the order of 10%.

We conclude that using a full spectral method remains relatively accurate for disc flarings not exceeding 0.125, and that local temperature variations do not significantly impact the error. Therefore, for production runs of GI, the highest concern should not be transient structures like spirals or clumps but rather the global disc flaring. If the full FFT method is to be used on set-ups where the disc flaring is not zero, it may be wiser to simply choose a higher disc aspect ratio h when used in the numerical method for computing SG.

Fig. D.1 Radial SG force of the power law disc for different disc flarings and when a temperature gap is introduced.

Appendix E Indirect term

The indirect term associated with the pull exerted by fluid ‘a’ on the star is $$\begin{array}{rl}{\mathrm{\Phi }}_{a,\text{ind}}& =-{\mathit{a}}_{\odot }\cdot \mathit{r}\\ & ={\iint }_{\text{disc}}{\int }_{z=-\mathrm{\infty }}^{\mathrm{\infty }}\frac{{\mathrm{\Sigma }}_a\left(r^{\mathrm{\prime }}\right)e^{-\frac{1}{2}{\left(z^{\mathrm{\prime }}/H_a^{sg}\right)}^2}}{\sqrt{2\pi }{H}_a^{sg}\left(r^{\mathrm{\prime }}\right)}\frac{(r^{\mathrm{\prime }}{\mathit{e}}_{r^{\mathrm{\prime }}}+z^{\mathrm{\prime }}{\mathit{e}}_{z^{\mathrm{\prime }}})}{{(r^{\mathrm{\prime }2}+z^{\mathrm{\prime }2})}^{3/2}}\cdot \mathit{r}{d}^2{\mathit{r}}^{\mathrm{\prime }}\\ & ={\iint }_{disc}\frac{{\mathrm{\Sigma }}_a\left(r^{\mathrm{\prime }}\right)}{\sqrt{2\pi }}r{r}^{\mathrm{\prime }}\cos (\phi -{\phi }^{\mathrm{\prime }}){\left(H_a^{sg}\right)}^{-3}{F}_1\left(d_a^2\right)d^2{\mathit{r}}^{\mathrm{\prime }}\end{array}$$(E.1) where $d_a=\parallel {\mathit{r}}^{\text{'}}\parallel /H_a^{sg}$ and F₁ is the function defined by Eq. A.1. Using the results of appendix A, we get $${\Phi }_{a,\text{ind}}\left(\mathit{r}\right)=\sqrt{2}r{\int \int }_{disc}{\Sigma }_a\left({\mathit{r}}^{\text{'}}\right)cos(\phi -{\phi }^{\text{'}}){\mathrm{K}}_a\left(d_a\right)d^2{\mathit{r}}^{\text{'}},$$(E.2) where ${\mathrm{K}}_a\left(d_a\right)=\frac{1}{\sqrt{\pi }}{\left(H_a^{sg}\right)}^{-2}\frac{d_a}{4}exp\left(\frac{d_a^2}{4}\right)[K_1\left(\frac{d_a^2}{4}\right)-K_0\left(\frac{d_a^2}{4}\right)]$. We emphasise that thus far, we have not used the vertical integral component of our reasoning. Indeed, it is worth noting that the vertical integration employed in Eq. E.1 is inherently encompassed within the definition of the 3D indirect term potential. Thus, to align this indirect term with our approach, we must calculate the following quantity: $$\begin{array}{rl}-{\int }_{z=-\mathrm{\infty }}^{\mathrm{\infty }}{\rho }_a(\mathit{r},z){\mathbf{\nabla }}_{3D}{\mathrm{\Phi }}_{a,\text{ind}}dz& =-[{\int }_{z=-\mathrm{\infty }}^{\mathrm{\infty }}{\rho }_a(\mathit{r},z)dz]{\mathbf{\nabla }}_{2D}{\mathrm{\Phi }}_{a,\text{ind}}\\ & =-{\mathrm{\Sigma }}_a(\mathit{r}){\mathrm{\nabla }}_{2D}{\mathrm{\Phi }}_{a,\text{ind}}\end{array}$$(E.3) We finally obtain the gradient of the indirect potential: $$\begin{array}{rl}{\mathbf{\nabla }}_{2D}{\mathrm{\Phi }}_{a,\text{ind}}& =\sqrt{2}\underset{disc}{\iint }{\mathrm{\Sigma }}_a\left({\mathit{r}}^{\mathrm{\prime }}\right){\mathrm{K}}_a\left({\mathit{r}}^{\mathrm{\prime }}\right)\\ & [\cos (\phi -{\phi }^{\mathrm{\prime }}){\mathit{e}}_r-\sin (\phi -{\phi }^{\mathrm{\prime }}){\mathit{e}}_{\phi }]d^2{\mathit{r}}^{\mathrm{\prime }}\\ & =\left(\begin{array}{rrr}{A}_{disc}\cos (\phi )& +& B_{disc}\sin (\phi )\\ -A_{disc}\sin (\phi )& +& B_{disc}\cos (\phi )\end{array}\right)\end{array}$$(E.4) with $A_{\text{disc}}=\sqrt{2}{\int \int }_{\text{disc}}{\Sigma }_a\left({\mathit{r}}^{\text{'}}\right){\mathrm{K}}_a\left({\mathit{r}}^{\text{'}}\right)cos\left({\phi }^{\text{'}}\right)d^2{\mathit{r}}^{\text{'}}$ and $$B_{\text{disc}}=\sqrt{2}{\int \int }_{\text{disc}}{\Sigma }_a\left({\mathit{r}}^{\text{'}}\right){\mathrm{K}}_a\left({\mathit{r}}^{\text{'}}\right)sin\left({\phi }^{\text{'}}\right)d^2{\mathit{r}}^{\text{'}}.$$