Self-gravity in thin protoplanetary discs: 1. The smoothing-length approximation versus the exact self-gravity kernel

This paper introduces an exact self-gravity kernel based on modified Bessel functions for thin, hydrostatically supported protoplanetary discs, which overcomes the limitations of traditional smoothing-length approximations by preserving Newtonian features, ensuring computational efficiency, and revealing a previously unnoticed source of gravitational runaway at infinitesimal distances.

The Gist

Imagine a protoplanetary disc as a giant, spinning cosmic pizza dough. This dough is made of gas and dust, and it's slowly spinning around a baby star in the center. Over time, this dough doesn't just sit there; it tries to clump together to form planets.

For a long time, astronomers trying to simulate this process on computers faced a tricky problem: Gravity.

In the real 3D universe, gravity is strong when things are close and weak when they are far away. But when scientists try to simulate these discs on a 2D computer screen (flattening the pizza dough to save computing power), they have to guess how gravity works in this flat world.

The Old Way: The "Softening" Trick

For years, scientists used a shortcut called a Plummer potential (or "smoothing length"). Think of this like putting a fuzzy, blurry lens over your camera.

• The Problem: If you make the lens too blurry (too much "softening"), you can't see the fine details. You miss the strong gravitational pull between clumps of dust that are very close together. It's like trying to bake a cake but telling the oven to ignore the heat when the batter is right next to the heating element.
• The Fix: If you remove the blur entirely to see the details, the math breaks down. The computer thinks the gravity becomes infinitely strong, causing the simulation to explode or behave wildly.

The authors of this paper say, "We've been using a blurry lens for too long. We need a way to see the whole picture clearly, from the very center of the dough to the very edge, without the math breaking."

The New Way: The "Exact Kernel"

The team (led by S. Rendon Restrepo) has derived a new mathematical formula, which they call an exact self-gravity kernel.

Here is how to understand their breakthrough using a few analogies:

1. The "Smart Lens" vs. The "Blind Spot"
Imagine the old method (Plummer potential) as a pair of glasses that has a blind spot right in the center of your vision. If two dust particles get too close, the glasses just say, "I can't see you, so I'll pretend you aren't there." This leads to errors.
The new method is like a pair of smart glasses that adjust automatically.

• When particles are far apart: The glasses show you the standard 3D gravity (like looking at a mountain from a distance).
• When particles get very close: The glasses smoothly switch to a different mode that understands the disc is flat. It doesn't break; it just changes how it calculates the pull, ensuring the physics stays realistic even when things are squished together.

2. The "Double-Decker" Cake
The paper also deals with dust (the heavy sprinkles) and gas (the fluffy cake).

• In the old models, scientists often treated the dust and gas as if they were the same height or ignored how the dust settles into a thin layer at the bottom.
• The new formula accounts for the fact that dust can settle into a very thin, dense layer (like a heavy chocolate ganache on top of a light sponge). It calculates how the heavy dust layer pulls on itself and on the gas, even if the dust is much thinner than the gas. This is crucial because, in some cases, the dust's own gravity might be just as important as the gas's gravity.

3. The "Runaway" Effect
One of the most exciting findings is that their new formula reveals a hidden danger: Gravitational Runaway.

• The Analogy: Imagine a snowball rolling down a hill. As it rolls, it picks up more snow, gets bigger, and picks up snow even faster.
• The Science: The old "blurry lens" method often stopped this snowball from growing too big because it smoothed out the gravity. The new "exact" method shows that if a tiny clump of dust gets slightly denser, its gravity pulls in even more dust, which makes it denser, which pulls in more dust. This creates a runaway effect that can lead to the rapid formation of planets or even the collapse of the disc. The old methods were likely missing this because they were too "soft."

Why Does This Matter?

This paper is like upgrading the map for a road trip.

• Before: We were driving with a map that had a "Do Not Enter" zone in the middle of the road (the smoothing length). We had to guess what was happening there.
• Now: We have a high-definition, GPS-enabled map that shows every pothole and hill, even the tiny ones.

The Bottom Line:
The authors have created a new mathematical tool that allows computers to simulate planet formation much more accurately. It respects the laws of physics (Newton's laws) at all distances, handles the difference between gas and dust correctly, and reveals that planets might form faster and more violently than we previously thought.

They also note that while this new math is more complex, it can still be calculated quickly enough for modern supercomputers, meaning future simulations of how our solar system (and others) formed will be far more realistic.

Technical Summary

1. Problem Statement

Protoplanetary discs (PPDs) are often modeled using 2D approximations to achieve high resolution in the $(r, \phi)$ plane, neglecting the vertical dimension. A critical challenge in these simulations is accurately modeling Self-Gravity (SG).

• Current Limitations: The standard approach uses a Plummer potential with an ad-hoc smoothing length (SL), $\epsilon$, to avoid numerical singularities.
  • Arbitrariness: $\epsilon$ is often treated as a free parameter (typically $\epsilon \approx 0.6-1.2 H$, where $H$ is the scale height).
  • Physical Inaccuracies: This approach suppresses the Newtonian nature of gravity at short distances, violates Newton's third law (due to symmetry breaking in previous space-varying SL prescriptions), and fails to account for the vertical structure of the disc (gas and dust stratification).
  • Regime Failure: Existing corrections (e.g., Rendon Restrepo & Barge 2023) are valid only for light discs ($Q \gtrsim 5$) and break down in massive discs where gravitational instability (GI) occurs ($Q \lesssim 1$).
• The Goal: To derive an exact, analytical 2D self-gravity kernel valid for Gaussian-stratified thin discs that incorporates both gas and dust components, respects Newton's laws, and remains valid from light to massive discs.

2. Methodology

The authors developed an analytical framework based on the vertical hydrostatic equilibrium of a bi-fluid (gas + dust) disc.

• Vertical Stratification: Instead of assuming a razor-thin (Dirac $\delta$) disc, the authors assume a Gaussian vertical density profile for both gas and dust. They utilize a revised scale height ($H^{sg}$) that incorporates the gravitational influence of both the star and the disc's own mass (gas and dust), derived from Rendon Restrepo et al. (2025).
• Derivation of the Kernel:
  • They start with the 3D SG force per unit volume and perform a vertical double integration over the Gaussian profiles of the source and target fluids.
  • This integration yields a closed-form expression for the 2D force kernel ($K_{ab}$) acting between fluid $a$ and fluid $b$.
  • The resulting kernel is expressed using modified Bessel functions of the second kind ($K_0$ and $K_1$).
• Numerical Implementation:
  • The kernel was implemented in two codes: FargoCPT (using full FFT with logarithmic grids) and Nirvana-iii (using a semi-spectral method allowing linear grids).
  • The authors benchmarked the kernel against analytical solutions for power-law and exponential discs, as well as 3D numerical simulations of Gaussian spheres.

3. Key Contributions

• Exact Bessel Kernel: The derivation of a closed-form kernel (Eq. 8) that is valid for any disc mass (Toomre $Q$) and accounts for gas-dust interactions.
  $$K_{ab} = \frac{1}{\sqrt{\pi}} (H^{sg}_{ab})^{-2} \frac{d_{ab}}{8} \exp\left(\frac{d_{ab}^2}{8}\right) \left[ K_1\left(\frac{d_{ab}^2}{8}\right) - K_0\left(\frac{d_{ab}^2}{8}\right) \right]$$
  Where $d_{ab}$ is the normalized distance and $H^{sg}_{ab}$ is the root-mean-square scale height of the interacting fluids.
• Restoration of Symmetry: The kernel inherently satisfies Newton's Third Law (action-reaction symmetry) because the scale height definition ($H^{sg}_{ab}$) is symmetric with respect to the interacting positions $r$ and $r'$. This eliminates spurious radial accelerations found in previous space-varying SL methods.
• Bi-fluid Capability: The framework naturally handles four interaction types: gas-gas, dust-dust, gas-dust, and dust-gas, without requiring arbitrary additional parameters.
• Computational Efficiency: The kernel is compatible with Fast Fourier Transform (FFT) methods, allowing for $O(N \log N)$ complexity, which is crucial for global simulations.

4. Key Results

• Behavior at Short vs. Long Distances:
  • Long Range: The kernel scales as $\propto 1/s^2$, recovering standard 3D Newtonian gravity.
  • Short Range: The kernel scales as $\propto 1/s$, exhibiting a purely 2D gravitational nature. This is physically consistent with the vertical integration of a thin disc, contrasting with the Plummer potential which either vanishes (if $\epsilon > 0$) or diverges as $1/s^2$ (if $\epsilon \to 0$).
• Massive Discs and Runaway Instability:
  • Unlike previous prescriptions, the Bessel kernel scales inversely with the Toomre parameter ($Q$) in massive discs.
  • Crucially, it allows for a gravitational runaway at infinitesimal distances. In massive discs, a local density perturbation decreases the local scale height, which increases the kernel value, further amplifying the gravity. This feedback loop is physically necessary for modeling gravitational collapse but is suppressed by standard smoothing lengths.
• Validation:
  • Analytical Tests: The kernel successfully reproduces radial accelerations for power-law and exponential discs with relative errors of $10^{-3}$ to $10^{-4}$ as the aspect ratio ($H/r$) decreases.
  • 3D vs. 2D Comparison: When compared to 3D simulations of Gaussian density bumps, the Bessel kernel maintained errors below 3% across various scale heights. In contrast, Plummer potentials with $\epsilon/H = 1.2$ underestimated SG by up to 28%, and $\epsilon/H = 0$ overestimated it by up to 129%.
• Dust Effects: The study shows that for highly settled dust ($\eta = H_d/H_g \ll 1$), the dust-dust self-gravity can be enhanced by a factor of $1/\eta$, making dust SG comparable to gas SG even at low dust-to-gas ratios ($Z=0.01$).

5. Significance and Implications

• Obsolescence of Plummer Potentials: The paper argues that the widely used Plummer potential with a smoothing length is fundamentally flawed for 2D SG simulations, as it either artificially suppresses or amplifies gravity depending on the choice of $\epsilon$.
• Impact on Planet Formation: The new kernel suggests that previous 2D simulations of Gravitational Instability (GI) may have been inaccurate. The ability to capture the short-range runaway mechanism implies that the threshold for disc fragmentation and the formation of bound objects (planets/brown dwarfs) could be significantly different from current models.
• Future Work: The authors note that this paper provides the theoretical and numerical foundation. An accompanying paper will investigate how this exact kernel alters the outcomes of GI-driven planet formation scenarios.
• General Applicability: The method is applicable to any Gaussian-stratified thin disc, bridging the gap between "razor-thin" approximations and full 3D simulations, offering a computationally efficient yet physically rigorous alternative.

In summary, this work provides the first exact, symmetric, and computationally efficient 2D self-gravity kernel for protoplanetary discs, resolving long-standing inconsistencies in the treatment of vertical structure and short-range gravitational interactions.