Issue |
A&A
Volume 675, July 2023
|
|
---|---|---|
Article Number | A96 | |
Number of page(s) | 12 | |
Section | Numerical methods and codes | |
DOI | https://doi.org/10.1051/0004-6361/202346178 | |
Published online | 06 July 2023 |
Self-gravity in thin-disc simulations of protoplanetary discs: The smoothing length rectified and generalised to bi-fluids
1
Aix-Marseille Univ., CNRS, CNES, LAM,
Marseille,
France
e-mail: pierre.barge@lam.fr
2
Aix-Marseille Univ., CNRS, Centrale Marseille, IRPHE,
Marseille,
France
3
Leibniz-Institut für Astrophysik Potsdam (AIP),
An der Sternwarte 16,
14482
Potsdam,
Germany
e-mail: srendon@aip.de
Received:
17
February
2023
Accepted:
5
April
2023
Context. To mimic the evolution of protoplanetary discs (PPDs), 2D simulations that incorporate self-gravity must introduce a softening prescription of the gravitational potential. When the disc is only composed of gas, the smoothing length is proportional to the scale height of the gas. On the other hand, when a dust component is included, the question arises as to whether the smoothing length approach can still be used to not only quantify the dust self-gravity, but its gravitational interaction with gas as well.
Aims. We identified grey areas in the standard smoothing length formalism for computing self-gravity in PPDs entirely made up of gas. Our goal is to revisit the smoothing length approach, which can then be generalised to two phases, whereby the dust component may be considered as a pressureless fluid.
Methods. We used analytical developments to approximate the vertically averaged self-gravity when the smoothing length is not assumed to be constant, but by taking a spatial function instead.
Results. We obtained an analytical expression for the space-varying smoothing length, which strongly improves the accuracy of the self-gravity computation. For the first time, this method is generalised to address bi-fluid interactions in a PPD: two additional smoothing lengths are proposed for featuring an isolated dusty disc and gas-dust self-gravity interactions. On computational grounds, we prescribe the use of tapering functions for the purpose of avoiding numerical divergences. We also checked that our method continues to be compatible with standard fast Fourier transform algorithms and evaluated computational costs.
Conclusions. Our space-varying smoothing length allows us to: (i) solve the contradictions inherent in the constant smoothing length hypothesis; (ii) fit the 3D vertically averaged self-gravity with a high level of accuracy; and (iii) render it applicable to a bi-fluid description of PPDs with the use of two additional smoothing lengths. Such results are crucial to enable realistic 2D numerical simulations that account for self-gravity and are essential to improving our understanding of planetesimal formation and type I migration.
Key words: protoplanetary disks / planet–disk interactions / gravitation / hydrodynamics / methods: analytical / methods: numerical
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Despite huge developments in 3D numerical computations and the advent of progressively more sophisticated computational facilities, thin-disc (2D) simulations remain much less expensive and more widely used in the study of protoplanetary discs (PPDs). The 2D approximation relies on the vertical averaging of the 3D physical quantities and governing equations. When dealing with self-gravity (SG), the issue becomes more sensitive since the equations cannot be vertically integrated. In practise, in PPDs studies an equivalent smooth potential, known as the Plummer potential, is used to mimic (as much as possible) the vertically averaged SG force and, in this approximation, a smoothing length (SL) is introduced to account for the disc vertical stratification.
As expected, this approach applies to 2D studies of planet–disc interactions where different values have been suggested (Li et al. 2009; Dong et al. 2011). In particular, Masset (2002) and Huré & Pierens (2009) showed that in the planet case, the SL should be proportional to the scale height of the gas disc, whereas Müller et al. (2012) proposed that ϵg/Hg = 0.7, where ϵg and Hg are the gas SL and the pressure scale height, respectively. Müller et al. (2012) also explored the case of self-gravitating gas discs and found that for vertically isothermal and stratified structures, we should instead use ϵg/Hg = 1.2 to avoid a systematic overestimation of the SG term. They also showed that when accounting for vertical SG, a deviation in the vertical Gaussian distribution occurs and in this case, the SL is also proportional to the disc’s Toomre parameter.
The Müller et al. (2012) prescription clearly leads to small errors (<2%) at large distances, but it also has an important drawback (that is not, however, mentioned by the authors); namely, SG is strongly underestimated at short distances with 100% errors. Indeed, in agreement with Huré & Pierens (2009), the authors found that the accuracy in the approximation of SG terms is strongly improved if the SL is a space varying function (Müller et al. 2012, Fig. 13). However, despite a discernible curve shape and (for unknown reasons), the authors did not test analytical expressions that could adequately fit the numerical curve. Furthermore, they found that the best approximation is obtained when the SL vanishes at the singularity – this is a statement that is inconsistent with the divergence of the error in the SG computations (Figs. 12 and 13 of their paper). The reasons for this inconsistency are subtle mathematical details at the singularity that were not explained by the authors and which, to our knowledge, have gone unnoticed to date.
In a bi-fluid description of PPDs, a dust layer is embedded in the gas disc and the evolution of gas and dust are coupled by aerodynamic forces. In this case, the momentum equations for the gas and the dust, both contain two SG terms: one is attributed to the gas disc and the other to the dust layer. In other words, each fluid parcel (gas or dust) is submitted to the gravity of the gas disc and to the gravity of the dust sub-layer. In contrast to the gas disc, the dust layer has a thickness which is governed not by pressure but turbulent stirring in the vertical direction. The scale height of the dust layer can be estimated using the form: (Dubrulle et al. 1995; Weber et al. 2019), where
, α is the dimensionless constant featuring turbulent α-viscosity (Shakura & Sunyaev 1973), Scz is the vertical Schmidt number, and St is the Stokes number. A vertical averaging over the dust layer faces the same numerical issues than averaging over the gas disc. These problems can be solved in the same way as for a single gas disc and a simple extrapolation provides the dust disc’s SL, namely, ϵd ∝ 1.2Hd. 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 ϵd ≪ ϵg. Thus, any error between these two parameters may lead to an incorrect and significant underestimation of the dust layer contribution to the SG terms. The question becomes more complex when accounting for the crossed SG terms. At the present time, we do not possess any rigorous theoretical approach that would allow us to evaluate this gravitational bi-fluid interaction. We do, however, think that this problem can also be addressed with a SL and, intuitively, we can expect it to lie between the gas SL and the dust SL. To our knowledge, both issues raised here with respect to bi-fluids are new and have not been previously addressed in numerical and theoretical studies of PPDs.
The improvements we perform in this study are crucial to enabling realistic 2D numerical simulations when SG is taken into account and they are also particularly important in planetary formation theories to better understand the formation of planetesimals. Indeed, it makes it possible to study the gravitational fragmentation of the dust layer (Goldreich & Ward 1973), formation of coherent clumps under the streaming instability (Johansen & Youdin 2007), and gravitational clumping of dust particles trapped in a large-scale vortex (Barge & Sommeria 1995) or in a co-orbital trapping scenario. Indirectly, we found that this work could also have implications for studies of type I planet migration.
In this paper, our goal is to provide a method based on the SL formalism to accurately compute the SG terms in thin discs at small and large separations. We aim to introduce the substantial SG interaction of gas and dust, when the latter is considered as a pressureless fluid. In Sect. 2, we carry out a theoretical development from first principles which justifies the use of the Plummer potential. In Sect. 3, we explore the SG estimation through the prism of the SL paradigm. We correct and complete the Plummer potential formalism thanks to a spatially dependent SL and generalise it to the case of bi-fluids. In Sect. 4, we handle practical aspects of the study, such as the numerical divergences, computational costs, and the calculation of the corrected SG thanks to fast Fourier transforms (FFT) methods. Finally, in Sect. 5, we suggest a discussion followed by a conclusion.
2 Self-gravity terms for bi-fluid simulations
In this section, we set up the theoretical background for computing gas and dust SG contributions when solid material is considered as a pressureless fluid in a thin-disc approximation. In this context, we also recall the role of the Plummer potential in 2D SG calculations and generalise this formalism to bi-fluids.
2.1 2D approximation and formal derivation
In 3D, the SG force per unit volume exerted by the PPD on a gas and dust parcels are expressed as:
(1)
where Φg and Φd are the distinct gravitational potentials of the gas and dust discs, respectively. For the sake of generality and concision, we denote both fluid phases as a and b, which reduces each term on the right-hand-side of Eq. (1) to:
(2)
where is force per unit volume the a-disc exerts on an elementary b fluid element. The density of phase a is noted as ρa and
. Assuming a vertical hydrostatic equilibrium and an isothermal approximation in the vertical direction for gas, the volume density can be written as
, where Hg(r) is the gas pressure scale height. We assumed the same vertical Gaussian profile for dust density than for gas, but this time, the vertical equilibrium is instead governed by turbulent stirring, which sets a different scale height for dust: Hd(r). This permits us to write both surface densities as:
(3)
The 2D analogue of the SG force is simply obtained integrating the 3D SG force in the vertical direction:
(4)
where s = ‖r − r′‖ is the separation (or mutual distance) between two fluid elements and es = (r − r′) /s. In this 2D approximation, it is implicitly assumed that the disc is symmetric with respect to the z = 0 plane, which allows us to naturally cancel the vertical component of the above force during the integration. After substituting the variables, we finally obtain:
(5)
is a normalised quantity that we refer to as the ‘self-gravity force correction’ (SGFC)1, with respect to the 3D case, db(r) = s/Hb(r) is the separation normalised with respect to the scale height of fluid b, and ηab = Ha(r′)/Hb(r) is the a-to-b scale height ratio. In the following, we use the notation db = db(r), except when a distinction is necessary. We want to highlight that in our reasoning, ηab is a spatial varying quantity, but, henceforth, we prefer to use:
(7)
where 〈〉 stands for space averaging over the 2D disc. This simplifies the next theoretical developments. We note that the scale heights used in this paper could be time-varying functions, but for concision, the time dependence is not explicit in the equations.
The numerical determination of SG in thin-disc and bi-fluid simulations mainly relies on the computation of the SGFC, but no analytical expression of this integral has been found in terms of standard mathematical functions. This is why a Plummer potential is commonly used to most optimally approximate this integral.
2.2 Fitting self-gravity terms thanks to smoothing lengths
Both vertical averagings in Eq. (6) cannot be performed analytically, but (in practice) we can introduce a Plummer potential, , into the 2D gravitational potential,
, which permits us to approach the 2D SG force defined in Eq. (5):
(8)
and ϵab is the SL between phases a and b. Usually, ϵab is assumed to be constant, so that:
(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 dg; however, this is, in fact, inconsistent with Eq. (10). Indeed, the additional term should be present in the numerator of the aforementioned equation. So, in order to remain mathematically correct and to retain the possibility of making comparisons with the work of Müller et al. (2012), we decided to consider the potential that satisfies only Eq. (10) and not Eq. (9). This is a slight change in the SL paradigm that does not affect the approximation of the 2D SG terms and leads to meaningful 2D results. Based on this clarification, we can define an analogue of the SFGC that is compatible with the SL approach:
(11)
This normalised quantity, which we call the smoothing length force correction (SLFC) should fit the SGFC to correctly estimate SG in 2D simulations. This is only possible by wisely choosing a spatially dependent SL, namely, ϵab(dg), as depicted in the next section. The definition and abbreviation of the main quantities encountered in this paper are given in Table 1.
Definitions and list of abbreviations.
3 Mutual self-gravity interactions based on the smoothing-length approach
We aim here to provide a suitable SL that most closely approaches each of the possible exact gravitational interactions: gas-gas, dust-dust, and dust-gas. We started by retrieving and rectifying Müller et al. (2012) results for a disc only made of gas.
3.1 Contribution of the gas disc on a gas parcel
This case corresponds to a=b=g. From the general definition provided by Eq. (6), the gas SGFC is equal to2:
(12)
where ηg = 〈Hg(r′)〉/〈Hg(r)〉 = 13. We present in Fig. 1 (top panel) the exact SGFC (blue solid line), evaluated numerically, and the equivalent SLFC (dashed lines) for three constant smoothing lengths (CSLs): ϵg(dg)/Hg(r) = [0.0, 0.6, 1.2]. In the bottom panel, we show the respective errors in percentage. These curves support the results from Müller et al. (2012) for the case of a self-gravitating disc and also validates our approach. For long distances, the three aforementioned SLFC exhibit the same behaviour and the error becomes negligible. For instance, the CSL value, ϵg(dg) = 1.2 Hg, offers the smallest error at long distances, namely, less than 2% for dg ≥ 4. On the contrary, at short distances, the SLFC either vanishes or diverges, which leads to significant errors since the exact SGFC converges towards L0 ≃ 1.772. Such behaviour is indeed intrinsic to the CSL formalism at short distances:
(13)
The above limits demonstrate that the CSL formalism either leads to 100% or infinite errors at short distances, when ϵg = const. ≠ 0 and ϵg = 0, respectively. This divergence between the SGFC and SLFC was not highlighted by Müller et al. (2012). A possible reason for this is that the behaviour is suitable since it allows for numerical divergences to be avoided at dg = 0 when ϵg(dg) = const. ≠ 0. Even if Müller et al. (2012) did not shine a light on these divergences at short distances, they nevertheless carried out a numerical evaluation that demonstrates that a space-varying smoothing length (SVSL) is necessary to fit the exact SGFC (reproduced in the top panel of our Fig. 2). From this figure, we can infer that the SVSL should tend towards 0 so that the SLFC matches the SGFC. Yet this is in contradiction to Fig. 1 and Eq. (13), since the SLFC is supposed to diverge analytically for ϵg(dg) = 0 at short distances. In general, no matter the choice of CSL, the gas SG was underestimated by a factor 100 at short distances. This contradiction went unnoticed and we intended to resolve it, as described in the next section.
The analytic expression of the gas SL that provides a better fit to Fig. 2 (top panel) and thus permits the SLFC to fit the exact SGFC should fulfill the following constraints:
(14)
Thanks to analytical arguments completed by a curve fit (made explicit in Appendix A), we found a model for the SVSL which permits us to accurately approach the SGFC:
(15)
where ϵ0,g = [π/L0]1/3 and the numerical values of (α, n) are collected in Table 2. In particular, the power of 2/3 and of ϵ0,g allow us to reach a finite value of L0 for the SLFC at dg = 0. In Fig. 2 (bottom panel), we show the SVSL analytic model. In Fig. 1, we show the respective SLFC and error with respect to the exact SGFC (black cross markers). Within this correction, the SLFC and SGFC overlap and the error is reduced to less than 0.5% across the whole distance range. As a matter of comparison, this correction allows us to decrease the error up to factors 200 and 40 at short and long distances, respectively, compared to the CSL where ϵg(dg) = 1.2Hg(r).
![]() |
Fig. 1 Force corrections for the gas only. We compare the SGFC (exact value, in blue solid line) and the SLFC model for different smoothing lengths based on constant values for the SL (red, green and orange dashed lines) or on function єg(dg) (black cross markers). For a CSL, єg = const., the SLFC and SGFC do not match at short separations: the error curve either tends either towards 0 or towards infinity. On the contrary, for a SVSL, єg = єg(dg), the SLFC matches SGFC with an accuracy better than 0.5%. We found that |
![]() |
Fig. 2 Space-varying smoothing length for gas only. We compare the SVSL obtained numerically by Müller et al. (2012) and that proposed in this work which was obtained by analytical computations and curve fitting. Top: Müller et al. (2012). Bottom: this work. |
Fitting parameters for space varying smoothing lengths and δ models.
3.2 Contribution of the dust layer on a dust parcel
In this Sect., we take a = b = d. From the general definition provided by Eq. (6), the dust SGFC is4:
(16)
where ηd = 〈Hg(r′)〉/〈Hd(r)〉 = 1 and η = Hg(r)/Hd(r) is the gas-to-dust scale height ratio that, in general, is greater than unity. Again, for simplicity, we prefer to use η = 〈Hg(r)〉/〈Hd(r)〉. The above relation naturally applies to the dust SLFC, but it is convenient to write it one of two different ways:
(17)
Both equations are equivalent, but the former provides a physical insight while the latter allows for a simple implementation in hydrodynamical codes. In Fig. 3, the top and bottom panels show the SLFC and the error with respect to the SGFC for different gas-to-dust scale height ratios, respectively. The dust SLFC and SGFC curves overlap over the whole η and distance ranges, so we did not plot these two quantities in the same figure to avoid duplication. The error is again lower than 0.5% in the whole distance range, as expected from the unique gas disc case. We want to highlight that at long separations the dust SLFC matches the gas SLFC, , while at short distances, the dust SG is η times stronger than gas SG and we get
. This latter equality is important because (1) if the gas SLFC is used instead of the dust SLFC and (2) if a CSL is used instead of a SVSL – then the dust SG is underestimated by a factor of ~100 η at short distances. These aspects are of primary interest in an astrophysical context since it is typical to see η ≳ 10, which could have important implications for planetary formation theories (as discussed in Sect. 5.2).
![]() |
Fig. 3 Self-gravity force corrections for dust only. We show the dependence of the dust force corrections for different values of the gas-to-dust scale-height ratio (η). For the whole range of η values, the accuracy is better than 0.5%. The dust SGFC is proportional to η at short distances (dg ≤ 1.5/η). For large values of η, this could favour gravitational clumping. Top: SGFC. Bottom: relative error between SGFC and SLFC (in %): 100 · |SLFC/SGFC − 1|. |
3.3 Crossed contributions
In this section, we take a = d and b = g. We highlight that the SGFC is commutative with respect to phases a and b. This can be demonstrated thanks to Newton’s third law or based on analytical arguments (see Appendix B). Therefore, in the following, we will adopt the notation to refer to either the dust-gas or gas-dust SGFC. In the context of this clarification and Eq. (6), the dust-gas SGFC is:
(18)
Motivated by the SL formalism, we again aim to approach the above double integral by a dust-gas SLFC, but to our knowledge, this has never been done before. In the next sections, we explore, for first time, what constraints should be satisfied in order to construct a consistent SL, which allows us to accurately approach the exact dust-gas SGFC.
3.3.1 Limiting cases for weak and strong layering
When layering is minimal (η = 1), the gas and dust are fully mixed and we immediately retrieve the case of a pure gas disc:
(19)
On the other hand, when layering is strong (η ≫ 1), the dust layer is infinitely thin and we get:
(20)
where Ip(x) = xex [K1 (x) − K0(x)] and K1 and K0 are the modified Bessel functions of the second kind. Interestingly, 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 Lp the planet–disc force correction (PDFC). The PDFC can be approached in the SL formalism and using the same analytical approach conducted in Sect. 3.1 and the results in Müller et al. (2012, Sect. 4.1), and further described in Appendix A. Thus, we find the planet–disc SVSL:
(21)
where ϵ0,p, = (π/2)1/6 and the power 2/3 were obtained by analytical means, while (β, q) by curve-fitting. The latter parameters are gathered in Table 2. We note that, contrary to our initial guess, for strong dust layering, the dust-gas SLFC tends towards the PDFC – rather than towards the dust SLFC. In summary, the dust-gas SL equals the gas SL or planet–disc SL when the gas-to-dust height, η, tends towards unity or infinity, respectively.
3.3.2 Behaviours at short and long separation
To study dust-gas SGFC at short distances, it is convenient to take the following function:
(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, δ is equal to 1 for while for an infinitely thin layer of dust this quantity tends to
. Accounting for these boundary conditions and a noticeable curve shape, we get:
(23)
noting that the couple (γ, 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:
(24)
This implies that at long distances, the dust-gas SVSL is expected to be negligible compared to the square of the distance, , which is undoubtedly satisfied if the SVSL is constant at long distances.
![]() |
Fig. 4 Model functions for defining the dust-gas SVSL. Top: δ function: model and numerical estimation with respect to the gas-to-dust height ratio (η). Bottom: λ model with respect to η. |
3.3.3 Dust-gas space-varying smoothing-length construction
Now, we have on hand all the necessary information needed to build a consistent dust-gas SVSL. We made the choice to look for the dust-gas SVSL in the form of a linear combination of the previous asymptotic cases:
(25)
where λ(η) ∊ [0, 1]. For matching Sect. 3.3.1 constraints, we chose:
(26)
Additionally, from the results given in Sect. 3.3.2 and the Taylor expansion described in Appendix C, we get:
(27)
This λ function is plotted in the bottom panel of Fig. 4, where we checked that the above-mentioned boundary conditions are met. We also plotted, in the top panel of Fig. 5, the dust-gas SVSL, defined in Eq. (25), for different gas-to-dust scale height ratios. For η = 1, the dust-gas SVSL matches the gas SVSL while for modest dust layering (i.e. η ≥ 5), the SVSL tends rapidly to the one of the planet–disc interaction case. This allows us to use the approximation: ϵdg(dg) ≃ ϵp(dg) for η ≥ 5. Finally, we found that our results are mathematically consistent provided that . Comparing this theoretical prediction with the value of L0 obtained numerically, we find that both results match with an accuracy of up to six decimal points.
3.3.4 Summary for the dust-gas space-varying smoothing length
The middle and bottom panels of Fig. 5 show the exact dust-gas SGFC and the error between both estimations for different η values, respectively. The dust-gas SLFC and SGFC curves overlap for the whole η and distance range, so we did not plot these two quantities in the same figure to avoid duplication. This overlap is also reflected in the error, which this time depends on the dust layering: for η ∊ [1, 5[, the error is lower than 2% and for η ∊ [5, 100], the error is lower than 0.5% for the whole distance range. Compared to both cases studied in Sects. 3.1 and 3.2, the error is slightly higher but it remains acceptable. We want to stress that the dust-gas SLFC rapidly tends (with respect to η) toward the planet–disc SLFC, which makes the following simplification possible:
(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, θ) of hydrodynamical codes such as RoSSBi3D (Rendon Restrepo et al. 2022), FARGO (Masset 2000) or Athena (Stone et al. 2008). We start by treating the singularity for vanishing separations, responsible of numerical divergences, then we explain explicitly under which conditions results of Sect. 3 could be used for estimating SG thanks to FFT methods. Finally, we quantify the computational costs for 2D N-fluid simulations with SG.
4.1 Removing numerical divergences
From Eq. (11), it is obvious that is valid at short separations, which makes it necessary to have a product with a tapering function, so as to cancel SG for vanishing separations. The tapering function, f, must be equivalent to
at short distances, with p > 1. On the other hand, it should be equal to unity at large distances. In order to avoid losing the accuracy achieved in the steps described in Sect. 3, the tapering length should be approximately as long as the finest numerical resolution in the grid: rs ~ 3 min(∆r, r∆θ)5,6, where ∆r and ∆θ are the resolution in the radial and azimutal directions, respectively. To satisfy the above constraints, we arbitrarily chose p = 2 and the tapering function as follows:
(29)
where ds = rs/Hg 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 ds = 3/127. We also plotted the same quantity under a CSL assumption (red dotted line). Such a tapering length choice was motivated by the high resolution achieved by Rendon Restrepo & Barge (2022), namely: 146 cells/Hg and 127 cells/Hg in the radial and azimuthal directions. As expected, both Plummer potential gradients exhibit the same behaviour for dg ≳ ds 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 dg ≤ 1.5, but there is no need to resort to a tapering function since the potential gradient is cancelled out analytically at the singularity. This may be the reason why Müller et al. (2012) did not mention that there were 100% errors with the CSL method at short separations.
![]() |
Fig. 5 Force corrections and space-varying smoothing-length for dust-gas. We show the dependence of all quantities for different values of the gas-to-dust scale-height ratio (η). For thin (η ≫ 1) and thick (η = 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 |
![]() |
Fig. 6 Normalised Plummer potential gradient with and without a tapering function. In the SVSL method, the introduction of a tapering function (Eq. (29)) avoids the divergence of the Plummer potential. At the same time, the low tapering length, ds(= 3/127 ≃ 0.024) ensures that most of the SG contribution is not lost over short distances (dg ≤ 1.5). Under the CSL approach (red dotted), the gradient is cancelled out at the singularity, although most of the short distance contribution to the SG is lost. |
4.2 Numerical calculation with FFT methods
An important question regarding our results is whether FFT methods can still be used when we resort to a SVSL. Indeed, the classic SG computation in 2D was based on a CSL. We find it instructive to recall here how SG is computed when including a SVSL and a tapering function. Guided by Baruteau & Masset (2008) and Surville (2013, in french) we get the radial component of the SG mass force (divided by ∑b)7:
(30)
and hg(r) = Hg(r)/r is the gas disc aspect ratio. We aim to write the integral, defined by Eq. (30), as a convolution product that is only possible provided that the tapering function, f, and the ratio δϵαb could be written as and θ − θ′ functions. This constraint is satisfied if the disc aspect ratio is a spatial constant, ha(r) = h0,a, in the whole simulation box. Within this condition, the radial component of the SG force is:
(33)
is the modified radial Green’s function where we performed the variable substitution: r = eX and dr = eX dX. For the azimutal component, the calculation is similar and we obtain:
(34)
is the modified azimutal Green’s function. The rewriting of both integrals, expressed via Eqs. (33) and (34), permits us to finally write out the equations in terms of a set of Fourier transforms:
(35)
where 𝓕 and 𝓕−1 are the Fourier transform operator and its inverse. The symbol * denotes the convolution operator. In practice, such quantities are computed numerically thanks to fast Fourier modules which rapidly accelerate numerical computations. For instance, in the 2D version of RoSSBi3D (Rendon Restrepo et al. 2022), this computation is made possible thanks to the FFTW3 library (Frigo & Johnson 2005). It is important to highlight that the use of Fourier transforms is possible if: (1) a logarithmic mesh is used in the radial direction; this condition stems from the radial variable substitution that enabled us to obtain the formulation with the modified Green’s functions. Furthermore, (2) the function Σ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 a periodicity is artificially ensured thanks to a zero-padding approach: the radial domain is doubled and the function Σa is set to 0 in half of the domain. See Surville (Fig. III.10 2013, in french) for an example. Finally, it is necessary for (3) the disc aspect ratio of gas and dust to be constant with respect to r (but it could be a time-varying function). If such condition is not satisfied, the formulation as a convolution product is not possible.
We take this opportunity to clarify that the aforementioned assumption of a constant disc aspect ratio used in this section only permits us to resort to FFT methods for accelerating the numerical calculations. This condition is not new since it was already implicit in the classical calculation of the SG by FFT. For general scale heights, a direct summation in the radial direction and Fourier transforms in the azimuth direction, ensured by the periodicity, is a straightforward solution.
4.3 Computational costs
For SG simulations, the computational endeavour is a non-negligible aspect, particularly for high-resolution simulations. For concision and without any loss of generality, we only treat the computation of SG in the radial direction and we do not account the Fourier transforms of the modified Green’s functions, since they are computed only once. For the case of a bi-fluid simulation, the next Fourier transforms should be performed at each time step:
From the recapitulation above, we observe that the Fouriertrans-forms of gas and dust densities are duplicated, which reduces the total number of Fourier operations to six during a numerical treatment. Therefore, the numerical endeavour for computing SG for a bi-fluid is three times larger with respect to the case of a single fluid. In general, for N fluids with different scale heights and interacting through SG, the number of required Fourier transforms is:
(40)
Compared to standard self-gravitating simulations on a unique gaseous phase, the computation time spent in the self-gravity module is crudely multiplied by a factor ~N2/2 for large N. Of course, this amount can be decreased assuming that some of the involved fluids have the same scale heights.
5 Discussion
In this section, we treat the possible impact of our findings regarding planet migration and the early stage of planetary formation. Then we identify the limitations of the SVSL approach due to our initial assumptions and to the specificity of the studied problem. Finally, we propose possible ways to improve and test our model.
5.1 Consequences for planet-disc interaction
It is widely known that planets can migrate due to tidal interactions with the gas disc. This is the case of type I migration, whereby the planet can exchange angular momentum with Lindblad and co-rotation resonances (Baruteau & Masset 2013, for a review). In 2D numerical simulations, the value ϵp/Hg(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 ~1.5 Hg. Although we do not question the agreement between the 2D and 3D simulations of the planet–disc interaction, we do believe that our SVSL may be better suited than an adjustment factor. Therefore, it might be constructive to verify whether (1) the results of 3D simulations can be retrieved using our SVSL and (2) to check whether the widely used ϵp/Hg(r) = 0.3−0.6 factor could be retrieved from an analytical basis stemming from our SVSL. We also think that our SVSL approach could be helpful for addressing 2D simulations of low-mass planets embedded in a self-gravitating disc. This would, however, require some improvements: the vertical layering due to SG must be accounted for, as discussed in Sect. 5.3.1, and the vertical stratification of the gas disc due to the planet gravitation must be assessed, as noted by Müller et al. (2012). Both aspects raised in this paragraph require a more detailed work, which is out of the scope of present paper.
5.2 Consequences for planet formation theories
One attractive planetary formation scenario is based on the persistence of gaseous vortices in PPDs. Its main interest is in the strong capture efficiency of the vortices and their ability to confine large concentrations of dust grains that could collapse to form planetesimals or a planetary core (Barge & Sommeria 1995). Even if observational findings seem encouraging in this regard (Varga et al. 2021), numerical simulations have not yet concluded that vortices could form objects that are bound by gravity. The results of the present paper offer the possibility to carry out new numerical simulations that correctly account for SG in the gas and dust components of PPDs. In particular, we expect that these new simulations could answer the questions raised by the vortices. Indeed, as demonstrated in Sect. 3.2, dust SG could be underestimated by a factor of ~2000 for η = 20 at short separations. At the same time, the estimation of the dust-gas SLFC could allow us to quantify (with an acceptable level of accuracy) a possible gaseous envelope capture by dust clumps. From theoretical analyses, it has been found that SG inhibit vortices formation by a Rossby wave instability (Lovelace & Hohlfeld 2013) and this was later confirmed via numerical simulations (Baruteau & Zhu 2016; Regály & Vorobyov 2017; Tarczay-Nehéz et al. 2022). In the latest numerical work to date, Rendon Restrepo & Barge (2022) also found that vortices cannot survive in massive PPD and they provided a stability criterion that vortices should satisfy in order to resist SG destabilising effects. In light of our findings, the results of previous simulations should be checked anew to better understand the extent to which the SVSL affects theoretical predictions on the survival of vortices in self-gravitating PPDs.
5.3 Limitations, improvements, and tests
The limitations of our model are inherent to the initial assumptions we have made about the vertical structure of the disc in relation to gas and dust stratification.
5.3.1 Stratification and disc evolution
The SGFC studied in this paper is based on the vertical integration of Eq. (6) in the particular case of a vertically isothermal disc. However the vertical structure could be affected by different mechanisms, which implies that for any vertical stratification different from the Gaussian stratification, the entire work performed in this paper should be repeated and adapted. For instance, this is the case when including the disc vertical SG for gas (but neglecting the vertical component of the central object gravity), which modifies the Gaussian distribution into: , where Q is Toomre’s parameter (Lodato 2007). A similar layering should occur for the dusty disc and we also expect that the crossed gravitational interaction between both phases could impact their respective vertical stratification. Indeed, gas SG could decrease dust scale height by a factor of ~2 (Baehr & Zhu 2021) and we expect that a strong dust layering will also modify the gas vertical structure in correlation with the dust-to-gas density ratio.
The global vertical structure of PPDs evolves over time due to accretion heating (Schobert et al. 2019) and stellar irradiation (Wu & Lithwick 2021), amongst other effects. Instabilities could also generate time variable structures which could affect (locally) the vertical stratification of the flow. This was reported, for instance, in 3D vortices simulations (Meheut et al. 2010) for rings and gaps generated by poloidal magnetic winds (Suriano et al. 2017) and for spiral density waves (Riols & Latter 2018). The method described in this paper is not limited to a steady vertical stratification and it is compatible with the global and local time evolution of the vertical structure, provided that the time dependence between gas and dust scale heights is known. This could be done, for instance, with the 2D1D strategy adopted by Crida et al. (2009). However, in the particular case of spatially constant aspect ratios, ha(r, t) = ha(t), the computational benefits from the FFT method could be lost. Indeed, this would require us to perform the Fourier transforms of Green’s functions at each time step instead of a unique computation for a steady vertical stratification. For a simulation with only gas, this would result in three Fourier transforms at each time step (instead of two). We believe that for other SG computational methods, the computational costs would be unaffected.
5.3.2 Layering of the dust particles
In this paper, the dust component of PPDs is considered as a pressure-less fluid that is sufficiently mixed by the turbulent motions of the gas disc to be maintained in an equatorial sublayer. This assumption requires that the dust particles and gas aerodynamic coupling are governed by a Stokes number less than unity and that dust is adequately diluted in the gas to avoid frequent mutual interactions (if the dust-to-gas mass-ratio ≲1) (Garaud et al. 2004). We also assumed that small-scale turbulence is maintained in the gas disc by a mechanism that we ended up disregarding. It is interesting to note that outside the bi-fluid pressure-less assumption, turbulent stirring may be replaced by collisional and/or gravitational stirring with the formation of a sub-layer of solid particles, whose scale height is different from the one deduced from turbulent stirring in the introduction. In such cases, the necessary smoothing lengths will be different – but the SVSL approach should remain unchanged.
5.3.3 Additional test
A relevant test for our results consists of a comparison with the vertically averaged SG obtained thanks to a 3D simulation. In order to satisfy our assumptions, the 3D disc should be vertically isothermal and the vertical SG should be discarded. In addition, for the bi-fluid version, the dust layer should be as smooth as possible. Such a simulation is beyond the scope of current paper, but still remains interesting lead that we want to explore in the near future.
6 Conclusions
In this work, we reveal the contradictions and shortcomings of the CSL method commonly used to compute the contribution of SG in 2D numerical simulations. In particular, we found that from short to intermediate separations (dg ≤ 1.5), the SG force is analytically underestimated with an error that reaches 100% at the singularity. We corrected these inconsistencies replacing the CSL, ϵg = const., by a space dependent function, ϵg = ϵg(dg) (SVSL). We found that for a gas disc, the SVSL dependence that better fits the exact SGFC is:
(41)
where ϵ0,g = π1/6 and (α, n) can be found in Table 2. This SVSL approach can be extended to the dust disc using the following dust SL:
(42)
where η = 〈Hg〉/〈Hd〉 is the gas-to-dust height ratio. As a side result, we also found the planet–disc SVSL:
(43)
where and (β, 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:
(44)
where the analytical expression of λ is given in Eq. (27). All these new expressions for the SVSL are valid for any scale height if the stratification is Gaussian and remain compatible with the common FFT method for evaluating SG in 2D hydrodynamical simulations, provided that the gas and dust disc aspect ratios are constant. Finally, the use of a tapering function is required to avoid numerical divergences.
The proposed SVSL correction decreases the error up to factors 200 with respect to the latest CSL prescription proposed by Müller et al. (2012). In particular, our SVSL allows us to match the SGFC with a high accuracy, even at the singularity (dg → 0). Regarding the dust SLFC, we found that it is proportional to η at short separations. This result combined with the improvement brought by the SVSL method demonstrates that dust SG was generally underestimated by a factor ~100 η at short separations (compared to the CSL only based into the gas SL).
Our planet–disc SVSL could affect the torque exerted by the self-gravitating disc on a planet, which suggests that type I planet migration could be impacted. We also think that the improvements we have made in the computation of the SG terms will be decisive in future 2D simulations of PPDs inhabited by a large-scale vortex. Hence, we forecast a much better description of the evolution of the dust-gas mixture in the core of the vortices with our model than with standard ones; we also expect significant consequences in the simulations of planetesimal construction. Indeed, on the one hand, dust SLFC could favour the gravitational binding of the dust clumps trapped in the vortex and, on the other hand, dust-gas SLFC allows us to follow the coupled evolution of dust and gas in the various clumps. We speculate that gas could be dragged with dust during the collapse, before being trapped in a gaseous envelope around a dusty core.
Acknowledgements
We thank the referee for her/his helpful comments that enriched the discussion. We would like to warmly acknowledge Stéphane Le Dizès for fruitfull discussions during the preparation of the paper and financial support during the PhD thesis. S.R.R. thanks Clément Baruteau for useful discussions and Andrej Hermann for proofreading the article. Co-funded by the European Union (ERC, Epoch-of-Taurus, 101043302). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.
Appendix A Gas and planet space varying smoothing length derivation
From the results of Müller et al. (2012), we chose to look for a SVSL of the form:
(A.1)
with the constraint . From a Taylor expansion we get:
(A.2)
The SLFC converges towards L0 if (and only if) k = 2/3 and ϵ = [π/L0]1/3. This model can be highly improved thanks to an additional term in the exponential under the form . The values of (α, n) were obtained by numerical fitting to the exact SGFC. In Fig. A.1 (top panel), we compare the exact SGFC (blue solid line) with the SVSL obtained only by analytical means (orange dashed line) or the SVSL, on which we also introduced the numerical fit (black cross markers). We observe that the unique analytical correction decreases the error to less than 5% for the whole distance range.
The same procedure was used for estimating the SVSL of a planet interacting with a disc. In Fig. A.1 (bottom panel), we compare the SVSL obtained only by analytical means (orange dashed line) with the SVSL, which also includes the numerical fit (black cross markers). Again, the accuracy of the planet–disc SVSL method is excellent compared to the CSL method.
![]() |
Fig. A.1 Gas and planet–disc force corrections. For each force correction, we plotted the exact SGFC evaluated numerically (blue solid line) and the SLFC approximation accounting only for the analytical correction (orange dashed line) or accounting for an additional numerical fit (black cross markers). Top: Gas force correction. Bottom: Planet-disc force correction. |
Appendix B Commutativity of the self-gravity force correction for different phases
Thanks to the principle of action-reaction applied to the whole dusty disc and the whole gas disc8 we get:
(B.1)
which implies that . This commutativity is also easily recovered from Eq. 6:
(B.2)
As a consequence, throughout the whole paper, we only use the notation .
Appendix C Derivation of λ and L0
The Taylor expansion of Eq. 11, where a=d and b=g, in the vicinity of dg ~ 0, is:
(C.1)
However, the SLFC also satisfies (Eq. 22). The equalisation of both equations leads to:
(C.2)
Since δ(η = 1) = 1, it is immediately evident that λ(η = 1) = 1. The constraints of Sect. 3.3.1 require that the λ function is cancelled out for infinite gas-to-dust scale height ratios. Considering that , the above condition is only possible if:
(C.3)
where and
. The mathematical relation C.3 implies the constraint
.
Reference
- Baehr, H., & Zhu, Z. 2021, ApJ, 909, 136 [NASA ADS] [CrossRef] [Google Scholar]
- Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
- Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
- Baruteau, C., & Masset, F. 2013, in Lecture Notes in Physics, eds. J. Souchay, S. Mathis, & T. Tokieda (Berlin: Springer Verlag), 861, 201 [NASA ADS] [CrossRef] [Google Scholar]
- Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927 [NASA ADS] [CrossRef] [Google Scholar]
- Crida, A., Masset, F., & Morbidelli, A. 2009, ApJ, 705, L148 [NASA ADS] [CrossRef] [Google Scholar]
- Dong, R., Rafikov, R. R., Stone, J. M., & Petrovich, C. 2011, ApJ, 741, 56 [NASA ADS] [CrossRef] [Google Scholar]
- Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [Google Scholar]
- Frigo, M., & Johnson, S. G. 2005, Proc. IEEE, 93, 216 [Google Scholar]
- Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292 [NASA ADS] [CrossRef] [Google Scholar]
- Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [Google Scholar]
- Huré, J. M., & Pierens, A. 2009, A&A, 507, 573 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [Google Scholar]
- Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52 [NASA ADS] [CrossRef] [Google Scholar]
- Lodato, G. 2007, Nuovo Cimento Rivista Serie, 30, 293 [NASA ADS] [Google Scholar]
- Lovelace, R. V. E., & Hohlfeld, R. G. 2013, MNRAS, 429, 529 [NASA ADS] [CrossRef] [Google Scholar]
- Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Masset, F. S. 2002, A&A, 387, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Regály, Z. & Vorobyov, E. 2017, MNRAS, 471, 2204 [CrossRef] [Google Scholar]
- Rendon Restrepo, S., & Barge, P. 2022, A&A, 666, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Rendon Restrepo, S., Barge, P., & Vavrik, R. 2022, ArXiv e-prints [arXiv:2207.04252] [Google Scholar]
- Riols, A., & Latter, H. 2018, MNRAS, 476, 5115 [Google Scholar]
- Schobert, B. N., Peeters, A. G., & Rath, F. 2019, ApJ, 881, 56 [NASA ADS] [CrossRef] [Google Scholar]
- Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
- Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
- Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2017, MNRAS, 468, 3850 [NASA ADS] [CrossRef] [Google Scholar]
- Surville, C. 2013, Ph.D. Thesis, Aix-Marseille Université, France [Google Scholar]
- Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [Google Scholar]
- Tarczay-Nehéz, D., Rozgonyi, K., & Regály, Z. 2022, MNRAS, 511, 6055 [CrossRef] [Google Scholar]
- Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178 [NASA ADS] [CrossRef] [Google Scholar]
- Wu, Y., & Lithwick, Y. 2021, ApJ, 923, 123 [NASA ADS] [CrossRef] [Google Scholar]
Müller et al. (2012) defined as force correction a quantity that they defined as Isg. We noticed that in some of their graphics, they mistook sIsg(s) for Isg(s). Therefore, for clarity and consistency with their work, we did not adopted this naming convention.
We think that Müller et al. (2012) implicitly assumed this equality in their work.
All Tables
All Figures
![]() |
Fig. 1 Force corrections for the gas only. We compare the SGFC (exact value, in blue solid line) and the SLFC model for different smoothing lengths based on constant values for the SL (red, green and orange dashed lines) or on function єg(dg) (black cross markers). For a CSL, єg = const., the SLFC and SGFC do not match at short separations: the error curve either tends either towards 0 or towards infinity. On the contrary, for a SVSL, єg = єg(dg), the SLFC matches SGFC with an accuracy better than 0.5%. We found that |
In the text |
![]() |
Fig. 2 Space-varying smoothing length for gas only. We compare the SVSL obtained numerically by Müller et al. (2012) and that proposed in this work which was obtained by analytical computations and curve fitting. Top: Müller et al. (2012). Bottom: this work. |
In the text |
![]() |
Fig. 3 Self-gravity force corrections for dust only. We show the dependence of the dust force corrections for different values of the gas-to-dust scale-height ratio (η). For the whole range of η values, the accuracy is better than 0.5%. The dust SGFC is proportional to η at short distances (dg ≤ 1.5/η). For large values of η, this could favour gravitational clumping. Top: SGFC. Bottom: relative error between SGFC and SLFC (in %): 100 · |SLFC/SGFC − 1|. |
In the text |
![]() |
Fig. 4 Model functions for defining the dust-gas SVSL. Top: δ function: model and numerical estimation with respect to the gas-to-dust height ratio (η). Bottom: λ model with respect to η. |
In the text |
![]() |
Fig. 5 Force corrections and space-varying smoothing-length for dust-gas. We show the dependence of all quantities for different values of the gas-to-dust scale-height ratio (η). For thin (η ≫ 1) and thick (η = 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 |
In the text |
![]() |
Fig. 6 Normalised Plummer potential gradient with and without a tapering function. In the SVSL method, the introduction of a tapering function (Eq. (29)) avoids the divergence of the Plummer potential. At the same time, the low tapering length, ds(= 3/127 ≃ 0.024) ensures that most of the SG contribution is not lost over short distances (dg ≤ 1.5). Under the CSL approach (red dotted), the gradient is cancelled out at the singularity, although most of the short distance contribution to the SG is lost. |
In the text |
![]() |
Fig. A.1 Gas and planet–disc force corrections. For each force correction, we plotted the exact SGFC evaluated numerically (blue solid line) and the SLFC approximation accounting only for the analytical correction (orange dashed line) or accounting for an additional numerical fit (black cross markers). Top: Gas force correction. Bottom: Planet-disc force correction. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.