The Disk model

The Disk class computes the physical structure of a protoplanetary disk: surface densities, scale heights, temperatures, densities, and dust settling.

Overview

A Disk is built from a set of parameters (params.disk) and a dust model. It provides methods to compute:

• Gas structure: surface density, scale height, midplane/atmosphere temperatures, vertical density profile.

• Dust structure: dust surface density per grain size (with settling), dust scale heights, dust number densities.

• Optical properties: extinction efficiency, visual extinction \(A_V(z)\).

The disk model supports two surface density modes:

• Parametric (sigma_compute='param'): power-law with exponential taper.

• Custom (sigma_compute='custom'): interpolated from a user-provided file (e.g., from an MHD simulation).

Key concepts

Surface density and settling

The gas surface density follows:

\[\Sigma_g(r) = \Sigma_0 \left(\frac{r}{r_0}\right)^{-p} \exp\left[-\left(\frac{r}{r_c}\right)^{2-p}\right]\]

Following the prescription by Cuzzi et al. (1993), Youdin & Lithwick. (2007), Dong et al. (2015), dust settling reduces the scale height of larger grains relative to the gas:

\[H_d(a, r) = H_g(r) \frac{1}{\sqrt{1 + T_{s,mid} \frac{S_c}{\alpha}}}\]

where \(\mathrm{T_{s,mid}}\) is the dimensionless stopping time in the midplane, \(\alpha\) is the turbulence parameter, and \(S_c\) is the Schmidt number.

In params.disk, \(\alpha\) and \(\S_c\) correspond to alpha and schmidtnumber, respectively.

Vertical temperature structure

The temperature transitions from midplane (\(T_\mathrm{mid}\)) to atmosphere (\(T_\mathrm{atm}\)) following a smooth profile controlled by the parameter delta.

API Reference

class astromugs.modeling.Disk.Disk(params: DiskParams, dust: object | None = None)

Bases: object

Static flared disk model for radiative transfer and chemistry simulations.

Provides methods to compute gas and dust physical quantities (density, temperature, scale height, surface density, extinction) on various coordinate grids.

Parameters:

• params (DiskParams) – Dataclass holding all disk model parameters.

• dust (object, optional) – Dust model object providing grain size distributions and opacities.

scaleheight(r)

Compute the gas scale height.

Parameters:

r (float or array_like) – Distance(s) from the star in AU.

Returns:

hgas – Gas scale height in cm.

Return type:

float or ndarray

Notes

When a custom temperature profile is available (from the custom surface density file), the scale height is computed directly from H = sqrt(k_B T r^3 / (mu m_H G M_*)) at each radius. Otherwise, a power-law profile is used via h0 or tmidplan_ref.

surfacedensity(r)

Compute the gas surface density.

Parameters:

r (float or array_like) – Distance(s) from the star in AU.

Returns:

sigma_g – Gas surface density in g/cm^2.

Return type:

float or ndarray

omega2(r)

Compute the square of the Keplerian angular velocity.

Parameters:

r (float or array_like) – Distance(s) from the star in AU.

Returns:

Keplerian angular velocity squared in s^-2.

Return type:

float or ndarray

temp_mid(r)

Compute the midplane temperature radial profile.

Parameters:

r (float or array_like) – Distance(s) from the star in AU.

Returns:

Midplane temperature in K.

Return type:

float or ndarray

Notes

When a custom temperature profile is loaded from the surface density file, the temperature is interpolated from that table. Otherwise, the analytic power law T = T_ref * (r/r_ref)^(-q) is used.

temp_atmos(r)

Compute the atmospheric temperature radial profile.

Parameters:

r (float or array_like) – Distance(s) from the star in AU.

Returns:

Atmospheric temperature in K.

Return type:

float or ndarray

temp_altitude(r, z)

Compute the vertical temperature profile following Williams & Best (2014).

Parameters:

• r (array_like) – Radial distances from the star in AU.

• z (array_like) – Normalized altitude coordinates (in units of scale height).

Returns:

Temperature array of shape (len(r), len(z)) in K.

Return type:

ndarray

verticaldensity_gauss(r)

Compute the gas number density assuming a Gaussian vertical profile.

Parameters:

r (array_like) – Radial distances from the star in AU.

Returns:

Gas number density in cm^-3.

Return type:

ndarray

Notes

Valid only for vertically isothermal disks. The surface density is divided by the mean molecular mass to yield number density rather than mass density.

density(x1, x2, x3=None)

Compute the gas mass density assuming hydrostatic equilibrium.

Parameters:

• x1 (array_like) – Radial coordinate in AU (spherical: radii; nautilus: radii).

• x2 (array_like) – Second coordinate (spherical: colatitude in rad; nautilus: normalized altitudes in units of scale height).

• x3 (array_like, optional) – Third coordinate (spherical: azimuthal angle in rad). Required for spherical grids.

Returns:

rhog – Gas mass density in g/cm^3.

Return type:

ndarray

Notes

The Gaussian vertical profile is an isothermal approximation. For a non-vertically isothermal profile, the density can be computed iteratively and slightly deviates from a Gaussian profile.

numberdensity(x1, x2, x3=None)

Compute the gas number density assuming hydrostatic equilibrium.

Parameters:

• x1 (array_like) – Radial coordinate in AU.

• x2 (array_like) – Second coordinate (colatitude in rad or normalized altitude).

• x3 (array_like, optional) – Azimuthal angle in rad (spherical grids only).

Returns:

ng – Gas number density in cm^-3.

Return type:

ndarray

Notes

If vertically isothermal, the profile is Gaussian. Otherwise the density is computed iteratively and deviates from a Gaussian.

viscous_accretion_heating(x1, x2, x3=None)

Compute the viscous accretion heating rate.

Parameters:

• x1 (array_like) – Radial coordinate in AU.

• x2 (array_like) – Colatitude in rad.

• x3 (array_like, optional) – Azimuthal angle in rad.

Returns:

q_visc – Volumetric viscous heating rate in erg/(cm^3 s).

Return type:

ndarray

Notes

Currently implemented for spherical coordinates only. The heating is distributed vertically as a Gaussian and truncated beyond lim_h scale heights.

surfacedensity_d(r)

Compute the dust surface density per grain-size bin.

Parameters:

r (float or array_like) – Distance(s) from the star in AU.

Returns:

• sigmad (ndarray) – Dust surface density for each grain-size bin in g/cm^2, shape (n_sizes, *r.shape).

• sig_single (float or ndarray) – Total (summed over sizes) dust surface density in g/cm^2.

scaleheight_d(r)

Compute the dust scale height for each grain size.

Parameters:

r (array_like) – Radial distances from the star in AU.

Returns:

• hd (ndarray) – Dust scale height for each grain size in cm, shape (n_sizes, *r.shape).

• hd_single (float or ndarray) – Dust scale height for the single representative grain in cm.

Notes

When settling is enabled, the dust scale height is reduced relative to the gas scale height based on the Stokes number, the Schmidt number, and the turbulent alpha parameter.

density_d(x1, x2, x3=None)

Compute the dust mass density for each grain-size bin.

Parameters:

• x1 (array_like) – Radial coordinate in AU.

• x2 (array_like) – Second coordinate (colatitude in rad or normalized altitude).

• x3 (array_like, optional) – Azimuthal angle in rad (spherical grids only).

Returns:

rhod – Dust mass density in g/cm^3, shape (n_sizes, len(x1), len(x2)[, len(x3)]).

Return type:

ndarray

numberdensity_d(x1, x2, x3=None)

Compute the dust number density for each grain-size bin.

Parameters:

• x1 (array_like) – Radial coordinate in AU.

• x2 (array_like) – Second coordinate (colatitude in rad or normalized altitude).

• x3 (array_like, optional) – Azimuthal angle in rad (spherical grids only).

Returns:

dens – Dust number density in cm^-3, shape (n_sizes, len(x1), len(x2)[, len(x3)]).

Return type:

ndarray

Examples

Access densities at the third radius for the second grain species:

>>> n_d = disk.numberdensity_d(x1, x2, x3)
>>> n_d[1, 2, :]

numberdensity_d_single(x1, x2, x3=None)

Compute the dust number density for a single representative grain.

Useful when RADMC-3D uses multiple grain sizes but the chemistry solver (e.g. Nautilus) requires a single grain population.

Parameters:

• x1 (array_like) – Radial coordinate in AU.

• x2 (array_like) – Second coordinate (colatitude in rad or normalized altitude).

• x3 (array_like, optional) – Azimuthal angle in rad (spherical grids only).

Returns:

dens_single – Dust number density in cm^-3, shape (len(x1), len(x2)).

Return type:

ndarray

Notes

If nb_species == 1, the result equals numberdensity_d[0].

q_ext(lam, A=2, q_c=4)

Compute the extinction efficiency for each grain size and wavelength.

Parameters:

• lam (array_like) – Wavelengths in microns.

• A (float, optional) – Extinction efficiency in the geometric-optics regime (default 2).

• q_c (float, optional) – Extinction efficiency at the critical wavelength (default 4).

Returns:

qext – Extinction efficiency, shape (n_sizes, len(lam)). Dimensionless.

Return type:

ndarray

av_z(lam, nd, r, z)

Compute the vertical visual extinction profile.

Parameters:

• lam (array_like) – Wavelengths in microns.

• nd (ndarray) – Dust number density array of shape (n_sizes, len(r), len(z)) in cm^-3.

• r (array_like) – Radial distances from the star in AU.

• z (array_like) – Normalized altitude coordinates (in units of scale height).

Returns:

avz – Visual extinction (A_V) of shape (len(r), len(z)), in magnitudes.

Return type:

ndarray

avz_ism(lam, A=2, q_c=4)

Compute visual extinction using the ISM H-column-density-to-Av conversion.

Parameters:

• lam (array_like) – Wavelengths in microns.

• A (float, optional) – Extinction efficiency in the geometric-optics regime (default 2).

• q_c (float, optional) – Extinction efficiency at the critical wavelength (default 4).