Theoretical Framework

The central goal of DebrisPy is to predict the azimuthally averaged surface density profile, \(\bar{\Sigma}(r)\), from a physically motivated distribution of orbital elements.

Note

For a complete derivation and foundational formalism, see Rafikov (2023).

The surface density of a debris disc, \(\Sigma(r, \phi)\), defined in polar coordinates, is directly linked to the disc’s observed luminosity. In this work, we focus on the azimuthally averaged surface density, \(\bar\Sigma(r)\), which captures the radial structure of the disc after integrating out angular dependence,

\[\bar\Sigma(r) = \frac{1}{2\pi} \int_0^{2\pi} \Sigma(r, \phi)\,\mathrm{d}\phi.\]

The theoretical foundation for this calculation is based on the framework introduced in Rafikov (2023), which we implement in a modular and extensible form in DebrisPy. In this section, we summarise the derivation of \(\bar\Sigma(r)\) under this approach.

The radial position of a particle in orbit with semi-major axis \(a\) and eccentricity \(e\) is determined by the Keplerian relation,

\[r = \frac{a(1 - e^2)}{1 + e \cos(\phi - \varpi)},\]

where \(\varpi\) is the apsidal angle relative to a fixed direction. In such a case, the orbit ranges between the periastron, \(r_p = a(1 - e)\), and apoastron, \(r_a = a(1 + e)\), distances. This range is particularly important when determining which particles contribute to the matter density at a specific orbital distance.

To proceed, we first define the mass distribution over semi-major axis. Specifically, we let \(\mu(a)\) be the total mass per unit \(a\), such that \(\mathrm{d}m = \mu(a) \mathrm{d}a\) is the mass of the particles which belong in the \((a, \mathrm{d}a)\) interval. This allows us to introduce the associated surface density profile in \(a\)-space,

\[\Sigma_a(a) = \frac{1}{2\pi a} \frac{\mathrm{d}m}{\mathrm{d}a} = \frac{\mu(a)}{2\pi a}.\]

Although this is not a physical quantity, it plays a central role in defining and calculating \(\bar\Sigma(r)\). The second quantity of central importance is the eccentricity profile of the particles.

The remainder of the derivation depends on whether eccentricity is specified uniquely, \(e = e(a)\), or instead follows a distribution \(\psi_e(e, a)\). In both cases, we focus solely on the azimuthally averaged density, so the results are independent of the apsidal angle, \(\varpi\). This is because the time-averaged spatial distribution of particles at an orbital distance \(r\) (i.e. the time a particle spends at a radial distance) is unaffected by orbital orientation which is defined via \(\varpi\).

Eccentricity as a Unique Function of Semi-Major Axis

We now consider the case where \(e = e(a)\), which is equivalent to assuming the eccentricity distribution is a delta function in \(a\)-space, \(\psi_e(e,a) = \delta(e - e(a))\).

It is helpful to define the boundary condition that determines whether or not particles on orbits with given \((a, e)\) can contribute to the surface density at a particular radial location \(r\), arising from \(r_p \leq r \leq r_a\).

This requirement can be rewritten as a lower bound on the eccentricity \(e(a)\) for a given semi-major axis \(a\). Specifically,

• For \(r < a\), the constraint becomes \(e(a) > 1 - \frac{r}{a}\),

• For \(r > a\), it becomes \(e(a) > \frac{r}{a} - 1\).

Both of these cases can be unified using the relation,

\[e(a) > \kappa(r,a), \quad \text{where} \quad \kappa(r,a) = \left|1 - \frac{r}{a} \right|.\]

In other words, only orbits with \(e(a)\) above this threshold can physically reach \(r\). Following this, we introduce an auxiliary mass function \(\eta(r|a)\), such that \(\eta(r|a)\,\mathrm{d}r\,\mathrm{d}a\) is the mass contributed at \((r, \mathrm{d}r)\) from particles with semi-major axis in the \((a, a + \mathrm{d}a)\) range. Naturally, this auxiliary mass function can be related back to \(\mathrm{d}m(a)\) using the condition

\[\mathrm{d}a \int_{r_p(a)}^{r_a(a)} \eta(r|a)\,\mathrm{d}r = \mathrm{d}m(a) = 2\pi a \Sigma_a(a)\,\mathrm{d}a.\]

Time-averaging over particle orbits implies the contribution to the matter density at \(r\) is proportional to \(\mathrm{d}t/\mathrm{d}r \propto 1/v_r\), where \(v_r\) is the radial velocity. Using this, it can be shown that

\[\eta(r|a) = \frac{C(a) r}{\sqrt{(r - r_p)(r_a - r)}} = \frac{2 \Sigma_a(a) r}{\sqrt{(r - r_p)(r_a - r)}}.\]

The azimuthally averaged surface density, which is our quantity of interest, is then defined as

\[\bar\Sigma(r) = \frac{\eta(r)}{2\pi r}, \quad \text{with} \quad \eta(r) = \int \eta(r|a)\,\mathrm{d}a.\]

Putting this together gives

\[\bar\Sigma(r) = \frac{1}{2\pi r} \int_{0}^{\infty} \eta(r|a)\,\mathrm{d}a = \pi^{-1} \int_{0}^{\infty} \frac{\Sigma_a(a) \theta(r,a)}{\sqrt{(r - r_p)(r_a - r)}}\,\mathrm{d}a,\]

where \(\theta(r,a)\) is a Heaviside step function ensuring the integrand vanishes for \(r\) outside the orbital domain \(r_p(a) \leq r \leq r_a(a)\).

Due to this constraint, the lower limit of the integral must be \(a > r/2\), rather than zero. To further simplify the integrand, we define the dimensionless kernel function

\[\Phi_e(r,a) = \frac{\theta(r,a)}{\sqrt{e^2(a) - \kappa^2(r,a)}}.\]

This kernel function fully maps the \(a - r\) relationship of particles. As a result, we arrive at the final expression used in this case:

\[\bar\Sigma(r) = \pi^{-1} \int_{r/2}^{\infty} \frac{\Sigma_a(a)}{a} \Phi_e(r,a)\,\mathrm{d}a.\]

This is the central result computed by the DebrisPy package in the case of unique eccentricities.

Distribution of Eccentricities per Semi-Major Axis

We now consider the general case where particles at fixed \(a\) follow a distribution \(\psi_e(e,a)\). The main constraint we have on this function is that it must be normalised:

\[\int_0^1 \psi_e(e,a)\,\mathrm{d}e = 1.\]

In this case, the auxiliary function becomes \(\eta(r|a,e)\), since we now need to consider many eccentricities for a particle with a given semi-major axis. The mass conservation condition becomes

\[\mathrm{d}a\,\mathrm{d}e \int \eta(r|a,e)\,\mathrm{d}r = \psi_e(e,a) \,\mathrm{d}m(a)\,\mathrm{d}e.\]

Following similar steps as before, one can show that

\[\bar\Sigma(r) = \pi^{-1} \int_0^{\infty} \mathrm{d}a \int_0^1 \mathrm{d}e\, \frac{\Sigma_a(a)\,\psi_e(e,a)\,\theta(r,e,a)}{\sqrt{(r - r_p)(r_a - r)}}.\]

Conveniently, this can be written in the same form as in the unique eccentricity case, where the kernel is now defined by

\[\Phi_e(r,a) = \int_{\kappa}^{1} \frac{\psi_e(e,a)}{\sqrt{e^2 - \kappa^2(r,a)}}\,\mathrm{d}e.\]

Substituting this kernel into the expression for \(\bar\Sigma(r)\) again gives

\[\bar\Sigma(r) = \pi^{-1} \int_{r/2}^{\infty} \frac{\Sigma_a(a)}{a} \Phi_e(r,a)\,\mathrm{d}a.\]

This is the central result computed by the DebrisPy package in the case of a distribution of eccentricities.