Physical Model

Overview

AsteroidThermoPhysicalModels.jl is a comprehensive toolkit for thermophysical modeling of asteroids. This package allows you to simulate the temperature distribution on an asteroid, and predict the non-gravitational force (i.e., Yarkovsky and YORP effects).

The thermophysical model (TPM) considers the following physical processes:

1. Heat Conduction: Solves a one-dimensional heat conduction equation to model heat transfer from the surface into the interior of the asteroid.
2. Self-Shadowing: Accounts for local shadows cast by topography.
3. Self-Heating: Considers re-absorption of scattered light and thermal radiation from surrounding surfaces.
4. Mutual Shadowing: For a binary asteroid, accounts for eclipses between the primary and secondary bodies.
5. Mutual Heating: For a binary asteroid, considers thermal exchange between the primary and secondary bodies.

Symbols

The following symbols are used throughout the package to represent various physical quantities:

Heat Conduction Equation

Heat conduction within the asteroid is modeled by the following one-dimensional heat diffusion equation:

\[\rho C_p \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial z^2}\]

where:

• $T(z,t)$ is the temperature at depth $z$ and time $t$
• $\rho$ is the density
• $C_p$ is the specific heat capacity at constant pressure
• $k$ is the thermal conductivity

Numerical Solvers

AsteroidThermoPhysicalModels.jl provides three numerical methods to solve the heat conduction equation:

1. Explicit Euler Method (ExplicitEulerSolver)

   • Forward difference in time
   • Conditionally stable: requires $\lambda = \alpha \Delta t / \Delta z^2 < 0.5$
   • First-order accurate in time
   • Fast for small time steps

2. Implicit Euler Method (ImplicitEulerSolver)

   • Backward difference in time
   • Unconditionally stable for any time step
   • First-order accurate in time
   • Requires solving a tridiagonal system

3. Crank-Nicolson Method (CrankNicolsonSolver)

   • Average of forward and backward differences
   • Unconditionally stable
   • Second-order accurate in both time and space
   • Best balance of accuracy and stability

The solver can be specified when creating the thermophysical model:

# Example: Using Crank-Nicolson solver
stpm = SingleAsteroidTPM(shape, thermo_params;
   SELF_SHADOWING = true,
   SELF_HEATING   = true,
   SOLVER         = AsteroidThermoPhysicalModels.CrankNicolsonSolver(thermo_params),
   BC_UPPER       = AsteroidThermoPhysicalModels.RadiationBoundaryCondition(),
   BC_LOWER       = AsteroidThermoPhysicalModels.InsulationBoundaryCondition(),
)

Boundary Conditions

Upper Boundary Condition (Surface)

At the surface, a radiative equilibrium boundary condition is applied:

\[-k \frac{\partial T}{\partial z}\bigg|_{z=0} = (1-R_\text{vis})(F_\text{sun} + F_\text{scat}) + (1-R_\text{ir})F_\text{rad} - \varepsilon \sigma T^4\]

where:

• The left side represents the heat flux from the surface into the interior
• The right side represents the energy balance at the surface (absorbed solar radiation, scattered light, and thermal radiation minus emitted thermal radiation)

Lower Boundary Condition

At the lower boundary, an insulation condition is typically applied:

\[\frac{\partial T}{\partial z}\bigg|_{z=z_\text{max}} = 0\]

Thermal Inertia

Thermal inertia is a physical quantity that represents the ability of a material to resist temperature changes, defined by:

\[\Gamma = \sqrt{k \rho C_p}\]

The unit is tiu (thermal inertia unit) or J·m⁻²·K⁻¹·s⁻¹/².

Thermal Skin Depth

The thermal skin depth is a characteristic length that represents how far a periodic thermal wave penetrates into a material:

\[l = \sqrt{\frac{4\pi P k}{\rho C_p}}\]

where $P$ is the period of the thermal cycle (typically the rotation period of the asteroid).

Non-Gravitational Effects

Yarkovsky Effect

The Yarkovsky effect is an orbital perturbation caused by the asymmetric thermal emission resulting from the day-night temperature difference due to the asteroid's rotation. This effect primarily affects the semi-major axis of the asteroid's orbit.

YORP Effect

The YORP effect (Yarkovsky-O'Keefe-Radzievskii-Paddack effect) is a rotational perturbation resulting from thermal emission due to the asymmetric shape of the asteroid. This effect influences the rotation rate and the orientation of the asteroid's spin axis.

Net force and torque

The photon recoil on each facet $i$ is integrated over the surface as

\[\mathbf{F} = \sum_i \mathbf{F}_i, \qquad \boldsymbol{\tau} = \sum_i \mathbf{r}_i \times \mathbf{F}_i\]

where $\mathbf{F}_i$ is the recoil force on facet $i$ (direct emission and reflection, plus the momentum of photons intercepted by other facets when self-heating is enabled) and $\mathbf{r}_i$ is the position of the facet centre. The net force is the plain vector sum: where a force acts does not enter the motion of the centre of mass, only the torque. The torque is taken about the body-fixed origin, which is therefore assumed to coincide with the centre of mass.

Facets with a roughness model

When a facet $i$ of a ShapeModel carries a roughness model, the recoil is computed on every sub-facet $j$ of that model in its local frame — direct emission and reflection, plus the momentum of photons intercepted by the other sub-facets, which is always included — and the facet force is replaced by their sum. The roughness model is a patch that represents the surface of the facet statistically, so the sum is counted for the area of the facet rather than for the area of the patch:

\[\mathbf{F}_i = \frac{A_i}{A_\mathrm{proj}} \, \mathbf{R}_i^{\mathsf T} \sum_j \mathbf{f}_j, \qquad A_\mathrm{proj} = \sum_j a_j \, (\hat{\mathbf{n}}_j \cdot \hat{\mathbf{z}})\]

where $\mathbf{f}_j$ and $a_j$ are the force and area of sub-facet $j$ in the units of the roughness model, $A_\mathrm{proj}$ is the area of the model projected onto its reference plane ($\hat{\mathbf{z}}$ is the local normal), and $\mathbf{R}_i$ rotates from the body frame to the local frame of facet $i$. The scale given to add_roughness_models! does not enter: the force on one patch grows as $\mathrm{scale}^2$ and the number of patches covering the facet falls as $\mathrm{scale}^{-2}$. The torque uses the facet centre $\mathbf{r}_i$ as the point of action; the torque of the patch about its own centre is smaller by the ratio of the patch size to the body size and is neglected.

When self-heating is enabled, the photons that leave the roughness model towards the sky and are intercepted by other facets $k$ of the global shape are accounted for as for a smooth facet, taking the emission of the patch as isotropic: the power $P_{\mathrm{sky},i} = (A_i / A_\mathrm{proj}) \sum_j E_j a_j f_{\mathrm{sky},j}$ that escapes the model, with $f_{\mathrm{sky},j}$ the sky view factor of sub-facet $j$, contributes $(P_{\mathrm{sky},i} / c) \sum_k f_{ik} \hat{\mathbf{d}}_{ik}$.

Irradiation of a roughness model by the other facets

A facet $i$ that carries a roughness model receives, on its sub-facets $m$, the scattered sunlight and the thermal radiation of every facet $j$ visible from $i$. The same far-field approximation as the view factor is used — both facets are small compared to their distance, so every sub-facet of $j$ sees every sub-facet of $i$ along the single direction $\hat{\boldsymbol{d}}_{ij}$ — and the exchange separates into an emitter side and a receiver side.

The emission of $j$ towards $i$ is that of its own roughness model in the direction $\hat{\boldsymbol{d}}_{ji}$, per unit area of the model's reference plane projected onto that direction:

\[E_j(\hat{\boldsymbol{d}}_{ji}) = \frac{1}{A_\mathrm{proj}\cos\theta_{ji}} \sum_n V_n(\hat{\boldsymbol{d}}_{ji})\,(\hat{\boldsymbol{n}}_n \cdot \hat{\boldsymbol{d}}_{ji})^+ \, E_n \, a_n\]

where $V_n$ tells whether sub-facet $n$ of $j$ is seen from $\hat{\boldsymbol{d}}_{ji}$ (it may be hidden by the model's own walls), $E_n$ is $\varepsilon\sigma T_n^4$ for thermal radiation and $R_\mathrm{vis} F_{\mathrm{sun},n}$ for reflected sunlight (single scattering), and $\cos\theta_{ji}$ is the inclination of $\hat{\boldsymbol{d}}_{ji}$ to the reference plane. A hot sunlit wall that faces $i$ makes $E_j$ larger than the Lambertian $\varepsilon\sigma T_j^4$ of the smooth facet — thermal-infrared beaming — and a shaded wall makes it smaller. A facet without a roughness model radiates as a Lambertian face.

The irradiance reaching $i$ is $F_{ij} = E_j(\hat{\boldsymbol{d}}_{ji}) f_{ij}$, the same form as the smooth-surface term $\varepsilon\sigma T_j^4 f_{ij}$. It is distributed over the sub-facets of $i$ that see $j$, each by its own inclination:

\[F_m \mathrel{+}= \frac{F_{ij}}{\cos\theta_{ij}} \, V_m(\hat{\boldsymbol{d}}_{ij}) \, (\hat{\boldsymbol{n}}_m \cdot \hat{\boldsymbol{d}}_{ij})^+\]

so that the wall facing $j$ is heated and the wall behind it is not — the same rule as for the sunlight, with $F_{ij}$ in place of the solar flux. The power received, $\sum_m F_m a_m$, equals $F_{ij} A_\mathrm{proj}$ up to the shadowing discretisation of the sub-facets. The sub-facet visibilities $V$ of every visible pair are geometric and are computed once when the state is built; the temperatures they weight are those of each time step. For a flat roughness model every sub-facet sees every neighbour with the inclination of the facet, and the sum reduces to the smooth-surface flux of the facet. The global-level fluxes of $i$ are not affected: they remain the smooth-surface baseline.

Direction-Dependent Radiance of a Rough Facet

A facet with a roughness model does not radiate as a Lambertian surface: the sunlit wall of a crater is hotter than the shadowed one, and which wall an observer sees depends on the viewing direction. This is thermal-infrared beaming. The radiance of such a facet towards an observer direction $\hat{\mathbf d}$ is evaluated as a post-processing step from the recorded sub-facet temperatures (roughness_radiance, directional_radiance):

\[L_i(\hat{\mathbf d}) = \frac{1}{A_\mathrm{proj}\,(\hat{\mathbf z}\cdot\hat{\mathbf d}_\mathrm{local})} \sum_j V_j(\hat{\mathbf d}_\mathrm{local})\,(\hat{\mathbf n}_j\cdot\hat{\mathbf d}_\mathrm{local})^+\,a_j\,\frac{\varepsilon\,B(T_j)}{\pi}\]

where $\hat{\mathbf d}_\mathrm{local}$ is the observer direction in the local frame of the facet, $V_j$ is 1 when sub-facet $j$ is visible from that direction (not hidden by the crater walls), $A_\mathrm{proj}$ is the projected area of the roughness model, and $B(T)$ is $\sigma T^4$ for the total radiance or the Planck function for the spectral radiance at a given wavelength. The denominator is the projected area of the patch as seen by the observer, so that an isothermal, unshadowed patch reduces to the Lambertian $\varepsilon B(T)/\pi$.

The brightness temperature $T_b(\hat{\mathbf d})$ is that of a blackbody with the same Lambertian radiance, $B(T_b)/\pi = L_i(\hat{\mathbf d})$ (brightness_temperature); emissivity is not divided out, so a smooth grey facet at temperature $T$ has $T_b = \varepsilon^{1/4} T$. Facets without a roughness model radiate as smooth Lambertian surfaces from their recorded surface temperature. The result is one value per facet, which a ray-caster can place on an image.

Binary Asteroid Systems

For binary asteroid systems, AsteroidThermoPhysicalModels.jl provides comprehensive modeling of thermal interactions between the primary and secondary bodies.

Coordinate Systems

The package uses the following coordinate conventions:

• Primary-fixed frame: The reference frame fixed to the primary body
• Secondary-fixed frame: The reference frame fixed to the secondary body
• r₁₂: Position vector from primary to secondary in the primary-fixed frame
• R₁₂: Rotation matrix from primary to secondary frame
• R₂₁ = R₁₂ᵀ: Rotation matrix from secondary to primary frame

Coordinate Transformations

For binary systems, coordinate transformations are handled automatically by the unified API:

# The unified API handles all transformations internally
update_flux_all!(btpm, r☉₁, r₁₂, R₁₂)

The package internally computes:

• Sun position in secondary frame: r☉₂ = R₁₂ * (r☉₁ - r₁₂)
• Primary position in secondary frame: r₂₁ = -R₁₂ * r₁₂

These transformations ensure accurate calculation of mutual shadowing and heating effects.