Multiphysics Coupling and Regime Selection

Compact-object calculations often combine fluid motion, magnetic fields, radiation, neutrinos, elasticity, reactions, and dissipative stresses. The right approximation is determined not by density alone but by ratios of microscopic and macroscopic scales. A consistent hybrid model states those ratios, uses overlapping asymptotic limits to connect regimes, and makes exchange terms cancel exactly. Otherwise the same pressure or energy can be counted twice, or disappear at a model switch.

Required background. Hilbert Stress–Energy from Metric Variation supplies total stress and exchange; Differential Equations supplies asymptotic limits, stiffness, and coupled evolution.

Helpful background. The Einstein Equation supplies spacetime coupling; Relativistic Matter and Compact Objects gives the physical hierarchy; and Computational Relativity supplies multiphysics evolution methods.

Dimensionless regime map

Let $L$ and $t_{\rm dyn}$ be the macroscopic length and time of the phenomenon being modeled. They need not equal the grid spacing and time step. Useful control ratios include

$${\rm Kn}=\frac{\lambda_{\rm mfp}}{L}, \qquad \tau_\epsilon=\int_{\rm path}\kappa_\epsilon\,d\ell, \qquad {\rm Da}=\frac{t_{\rm dyn}}{t_{\rm react}},$$ $${\rm De}=\frac{t_{\rm relax}}{t_{\rm dyn}}, \qquad \sigma_{\rm m}=\frac{b^2}{w}, \qquad \beta=\frac{2p}{b^2}.$$

The optical depth is an integral of opacity along a specified path, with particle energy $\epsilon$ evaluated locally along that path. Only in a homogeneous medium with $\kappa=\lambda_{\rm mfp}^{-1}$ and path length $L$ does $\tau=1/{\rm Kn}$ follow. Knudsen number tests local scale separation; optical depth to escape is generally nonlocal. Magnetization and plasma beta use the comoving magnetic magnitude $b^2$, Heaviside–Lorentz units, and material enthalpy $w=\varepsilon+p$. No one number selects the entire model: one neutrino energy group can be thick and another thin.

When ${\rm Kn}\ll1$, a fluid or controlled moment expansion is plausible; when ${\rm Kn}\gtrsim1$, kinetic angular structure matters. If reactions relax toward a stable local equilibrium, ${\rm Da}\gg1$ permits that equilibrium approximation; ${\rm Da}\ll1$ means that the corresponding composition is approximately frozen. A small Deborah number lets a dissipative flux track its Navier–Stokes value after its initial transient, while ${\rm De}\gtrsim1$ requires retaining its dynamics. Distinguish scattering, absorption, energy exchange, and flavor-changing reaction times: they enforce different equilibria. For neutrino applications, see Foucart 2023, §§ 3.2, 4.1, and 4.2.7.

Selecting transport and material models

For radiation or neutrinos, $\tau\gg1$ supports equilibrium or diffusion if source relaxation is also fast; $\tau\sim1$ requires transport through decoupling; and $\tau\ll1$ supports free streaming. A grey optical depth may be inadequate when opacity is strongly energy dependent. Two-moment closure is efficient when angular detail is secondary, while Monte Carlo or discrete ordinates are better suited to crossing beams and detailed spectra.

Elastic response requires a shear modulus, temperature below melting, and a deformation timescale shorter than plastic or defect relaxation. A fluid model is recovered as shear modulus tends to zero. Viscosity and heat conduction require small Knudsen number and small departures from equilibrium. A transient model retains relaxation; modern general-frame first-order theories provide another option, as explained in Causal Viscosity and Heat Conduction. The two smallness conditions are distinguished in Rocha et al. 2024, preprint § II.D.

For ideal MHD, microscopic kinetic lengths must be small compared with $L$, and magnetic diffusion must be slow compared with the evolution. In the simple scalar-conductivity model, the magnetic diffusivity is $\eta_{\rm m}=1/\sigma_{\rm cond}$, the diffusion time is $t_\eta\sim L^2/\eta_{\rm m}$, and the magnetic Reynolds number is $R_{\rm m}=vL/\eta_{\rm m}$. Thus $R_{\rm m}\gg1$ supports flux freezing on the advection time $L/v$; it does not exclude narrow reconnection layers. The conductivity’s short electric-relaxation time $\sigma_{\rm cond}^{-1}$ must not be confused with the long magnetic-diffusion time.

Force-free electrodynamics further requires negligible material inertia and pressure, a sufficient current supply, $E\cdot B=0$, and $B^2\gt E^2$ in any local orthonormal frame. Large $\sigma_{\rm m}$ helps justify this limit but does not enforce the other conditions. A current sheet can need resistive or kinetic continuation even when its exterior is force-free. See the derivations on resistive MHD and force-free electrodynamics.

Exchange terms and double counting

Write the total stress as

$$\begin{aligned} T_{\rm total}^{ab} &=T_{\rm material}^{ab}+T_{\rm EM}^{ab} +R_\gamma^{ab}+\sum_sR_{\nu_s}^{ab},\\ T_{\rm material}^{ab} &=T_{\rm perfect}^{ab}+\Delta T_{\rm elastic}^{ab} +\Delta T_{\rm diss}^{ab}. \end{aligned}$$

Here the elastic and dissipative terms are corrections to the selected perfect-fluid contribution, not additional complete material stress tensors. The split is valid only when the named pieces do not overlap. Exchange equations must cancel pairwise so $\nabla_aT_{\rm total}^{ab}=0$. Charge and lepton-number sources require the same cancellation.

A common double counting error includes trapped-neutrino energy and pressure in an equation-of-state table while also evolving the full transported neutrino stress. Another adds elastic shear stress to a material stress tensor whose calibration already includes it. A third applies radiative cooling to matter but fails to deposit the emitted energy in a transported radiation field. The cure is a declared variable partition: which degrees of freedom live in equilibrium matter, which are explicit fields, and which exchange source moves energy between them.

For example, $\nabla_aR^{ab}=-G^b$ must be paired with a $+G^b$ matter source. Neutrino number sources also determine the electron-fraction update; see Foucart 2023, §§ 4.2.3–4.2.4. Check the combined update as well as individual modules. In flat spacetime, integrated energy changes by boundary fluxes. In a general evolving spacetime, covariant conservation is not a conserved coordinate-energy integral: the balance includes geometric terms, and a conserved energy associated with time translation requires the appropriate symmetry.

Related methods and asymptotic overlap

A hard switch at $\tau=1$ is rarely physical because optical depth is nonlocal and energy dependent. Diffusion and collisionless streaming are controlled at opposite limits; neither is generally accurate in the decoupling layer, so $\tau\sim1$ is not automatically their overlap region. Prefer transport equations valid through that layer and recovering both limits. Blend reduced models only where both are demonstrably accurate. Blending weights should not create or remove conserved quantities. If variables differ, reconstruct a common stress tensor and conserved currents before projecting into the destination model.

Switch hysteresis can prevent rapid toggling when a state fluctuates around a threshold, but it introduces history dependence that must be specified. Floors and caps likewise define an unresolved model. Their contributions to energy, momentum, charge, and composition must be accumulated and reported.

At a fluid–kinetic interface, match particle and stress moments while recognizing that a few moments do not reconstruct a unique distribution. At an elastic–fluid interface, match normal traction and velocity, with tangential traction determined by the fluid’s lack of shear support. At an MHD–force-free interface, match electromagnetic fields and total flux while preventing material energy from vanishing when inertia is dropped.

Coupled validation and rerun boundary

The benchmark inputs are a one-dimensional stratified slab with prescribed opacity transitioning smoothly from $\tau\gg1$ to $\tau\ll1$, a material heat capacity, equilibrium emissivity satisfying detailed balance, and a radiation pulse or steady luminosity. Compare a full transport reference with the hybrid model in the diffusion interior, transition layer, and streaming exterior. Add a passive reaction variable with known relaxation time to test equilibrium versus frozen composition independently.

This is a proposed comparison, not an executed benchmark. An implementation should record source revision, numerical method and arithmetic, opacity and equation-of-state versions, thresholds, blending width, tolerances, input files, and output hashes. Use a stationary flat-spacetime slab first so that boundary luminosity and the material-plus-radiation energy supply an unambiguous budget; add gravitational transport only after that coupling check passes.

Run $N=64,128,256$ cells. At each resolution, halve the switching-layer width in physical units and separately hold it fixed to distinguish modeling from truncation error. Report $L_1$ errors against the transport reference, total energy and reaction-charge drift, interface flux jumps, number of switch events, and sensitivity to threshold and hysteresis. Check pure asymptotic limits by making the entire domain thick, thin, elastic, fluid, ideal, or resistive as applicable.

One failure mode is a nonconservative state reset when changing variables. Others are switch chatter, opacity-dependent grid artifacts, double-counted trapped pressure, source imbalance, inconsistent interpolation between tables, and a false plateau caused by compensating errors in two modules. A module passing its isolated test does not establish correctness of the coupled exchange.

Exercises

A zone has $\tau\gg1$ but reaction time longer than the dynamical time. Should its radiation be assumed in full chemical equilibrium with matter?

Solution

No. Large optical depth supports angular isotropy and diffusion, but chemical equilibrium also requires ${\rm Da}=t_{\rm dyn}/t_{\rm react}\gg1$. With slow reactions, particles may be trapped yet retain a nonequilibrium number or chemical potential. Transport and reaction times test different assumptions.

References

• Foucart, François. “Neutrino Transport in General Relativistic Neutron Star Merger Simulations.” Living Reviews in Computational Astrophysics 9, 1 (2023). DOI.
• Rocha, Gabriel S., David Wagner, Gabriel S. Denicol, Jorge Noronha, and Dirk H. Rischke. “Theories of Relativistic Dissipative Fluid Dynamics.” Entropy 26, 189 (2024). DOI. Open preprint.