TIGRESS Shearing-Box Simulations

• TIGRESS shearing-box simulations model a local, high-resolution patch of star-forming ISM by incorporating differential rotation, vertical stratification, magnetic fields, and SN feedback.
• The models employ advanced numerical methods including the Athena MHD code with CTU integrator, HLLD Riemann solver, and adaptive ray-tracing for radiative transfer and radiative cooling/heating.
• Synthetic observables such as dust polarization maps and outflow loading factors are derived to compare simulation outputs with observations from Planck, JWST, and other instruments.

TIGRESS (Three-phase ISM in Galaxies Resolving Evolution with Star formation and Supernova feedback) shearing-box simulations are a class of local, high-resolution, three-dimensional magnetohydrodynamical (MHD) models designed to capture the evolution and feedback-driven dynamics of the multiphase interstellar medium (ISM) in star-forming galactic disks. These simulations employ the shearing-box approximation to represent a corotating patch of a galactic disk that includes differential rotation, vertical stratification, supernova (SN) feedback, radiative heating and cooling, and, in recent extensions, full radiative transfer of ultraviolet (UV) fields. TIGRESS shearing-box simulations have been used to investigate turbulence, magnetic field amplification, galactic outflows, dust polarization, and radiative transfer in the ISM, providing synthetic observables to compare directly with, for example, Planck dust polarization and JWST mid-infrared data (Kim et al., 2019, Kim et al., 2020, Linzer et al., 2024).

1. Governing Equations and Shearing-Box Formalism

TIGRESS simulations solve the equations of ideal MHD in a local Cartesian frame corotating at a galactocentric radius $R_0$, with shearing-periodic boundary conditions in the radial direction and vertical gravitational stratification. The full system includes:

• Continuity:

$$\frac{\partial \rho}{\partial t} + \nabla\cdot(\rho\,\mathbf{v}) = 0$$

• Momentum: Includes pressure, Lorentz force, Coriolis force, tidal/centrifugal terms, external vertical gravity, and (optionally) self-gravity and SN momentum deposition.
• Energy:

$$\frac{\partial E}{\partial t} + \nabla\cdot[(E+P^*)\,\mathbf{v} - (\mathbf{v}\cdot\mathbf{B})\mathbf{B}/4\pi] = -\rho\mathbf{v}\cdot\nabla\Phi_{\rm ext} + {\cal H} - {\cal C} + S_{E, \rm SN}$$

where $E$ includes kinetic, thermal, and magnetic contributions.

• Induction:

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B})$$

with $\nabla\cdot\mathbf{B}=0$ enforced via constrained transport.

• (Optional) Self-Gravity:

$\nabla^2\Phi_{\rm self}=4\pi G\rho$

Shearing-box boundaries in $x$ re-map data across the radial edges with a $y$-offset to account for galactic shear; $y$ and $z$ boundaries are usually periodic (or outflow for winds and radiative transfer applications). The vertical gravity combines a stellar disk and dark matter halo component:

$$g_{\rm ext}(z) = 4\pi G\rho_* z / \sqrt{1+(z/z_*)^2} + \Omega^2 z$$

Shear parameter $q\equiv -d\ln\Omega/d\ln R \approx 1$ for flat rotation curves.

2. Simulation Domain, Resolution, and Numerical Methods

Typical TIGRESS shearing-box domains span $L_x=L_y=512$–$1024\,{\rm pc}$ with vertical extents $L_z=1.5$–$7\,{\rm kpc}$, uniformly gridded with cell sizes $\Delta x = 2$–$8\,{\rm pc}$. The Athena MHD code is used with the CTU (Corner Transport Upwind) integrator, HLLD Riemann solver, and constrained transport for divergence-free $B$. Time integration is directionally unsplit with operator-split radiative cooling/heating and, in some runs, subcycling to resolve short cooling times.

Sink particles, representing star clusters, form in converging, self-gravitating flows that exceed a threshold density. Each sink evolves according to a population synthesis model (usually Starburst99), setting the time-dependent rates for SN feedback and UV luminosity.

Table: Domain and Resolution Examples

Model	$L_x=L_y$ (pc)	$L_z$ (pc)	$\Delta x$ (pc)
R8-4pc	1024	6144	4
LGR4-2pc	512	3072	2

3. Physical Processes: SN Feedback, Radiative Transfer, and ISM Chemistry

Supernova feedback is modeled by injecting either $10^{51}\,$erg as thermal energy (if resolved) or the correct terminal momentum $p_{\rm ST}\approx(2$–$3)\times 10^5\,M_\odot\,{\rm km\,s}^{-1}$ (if unresolved) around the sink cluster's position. OB runaway SNe can be launched with ballistic drift.

Photoelectric heating from FUV, radiative cooling (atomic and ionic lines), and cosmic-ray and gas-dust energy exchanges are included. Recent TIGRESS-NCR runs additionally solve time-dependent photochemistry (H, $\rm H_2$, C, O species) and non-equilibrium ionization.

Radiative transfer (in TIGRESS-NCR (Linzer et al., 2024)) employs adaptive ray-tracing (ART), with rays launched from clusters and split adaptively to ensure angular and spatial resolution. FUV and LyC photons are transported, with frequency-dependent dust and gas opacities and photo-rate and energy density feedback to MHD, chemistry, and cooling modules.

4. Magnetic Field Initialization, Evolution, and Dynamo

TIGRESS boxes are initialized with a uniform azimuthal magnetic field, typically $B_{y,0}\sim2$–$3\,\mu$G, so midplane plasma $\beta\sim10$. SN-driven turbulence rapidly tangles the field, and, with galactic shear, drives a large-scale dynamo. Saturated turbulent and mean field strengths reach $\sim3$–$4\,\mu$G (Kim et al., 2019). The field is decomposed into mean (plane-averaged) and turbulent components for statistical analysis:

$$\overline{B}(z)\equiv \langle B(x,y,z)\rangle_{x,y},\quad \delta B = B - \overline{B}$$

Turbulent velocities and magnetic, sonic, and Alfvénic Mach numbers (${\cal M}_s$, ${\cal M}_A$) are monitored.

5. Synthetic Observables: Dust Polarization and Outflow Diagnostics

a. Dust Polarization Maps

LOS-integrated synthetic dust polarization maps are constructed using local 3D dust, density, and magnetic field, integrating Stokes parameters ($I_\nu$, $Q_\nu$, $U_\nu$) at 353 GHz. The dust opacity, temperature, and polarization fraction are fixed; LOS geometry from multiple midplane observer positions samples the physical variability. The full-sky maps are Healpix-tessellated ($N_{\rm side}=128$).

The angular power spectra are decomposed into $E$- and $B$-modes; the $E/B$ power ratio and TE correlation are computed and compared to Planck observations.

• Median $E/B$ ratio $R=1.4$–1.7 (Planck: $R\simeq1.91\pm0.10$)
• Median TE correlation $r_{TE}=0.20$–0.30 (Planck: $r_{TE}\simeq0.357\pm0.020$) The simulated slopes in radial projection, $\alpha_{EE}\approx-3.1$, $\alpha_{BB}\approx-3.3$, are steeper than Planck ($-2.4$ to $-2.5$) due to inertial-range dominance (Kim et al., 2019).

b. Outflow Loading Factors and PDFs

Horizontally-averaged vertical fluxes of mass, momentum, energy, and metals are extracted at various heights. TIGRESS provides distributions as joint PDFs $f_q(u,w)$ of outflow velocity and sound speed, fit by two-component analytic models (for cool and hot phases) parameterized by $\Sigma_{\rm SFR}$ and $Z_{\rm ISM}$. Analytic expressions determine mass, momentum, energy, and metal loading factors $\eta_q$ as functions of star formation and ISM metallicity (Kim et al., 2020).

6. Validation, Convergence, and Robustness

Resolution studies ($\Delta x=4$ vs. $8$ pc) confirm convergence of slopes and statistical quantities; >90% of SN events resolve the Sedov–Taylor phase at $\Delta x\leq8$ pc. Variants with or without hot/warm gas or projection effects yield consistent $E/B$ and $r_{TE}$ metrics. Loading factors and PDFs recovered from analytic fits match simulation results to within 10–20%. Time snapshots and multiple observer locations quantify physical and statistical variability; $E/B$ and $r_{TE}$ exhibit $\sim40$–$50$ Myr quasi-periodic oscillations and $18$–$45\%$ snapshot-to-snapshot rms variability (Kim et al., 2019, Kim et al., 2020).

7. Applications and Extensions: UV Radiation and ISM Diagnostics

TIGRESS-NCR simulations incorporate frequency-resolved, time-dependent ART for FUV and LyC transfer alongside MHD, non-equilibrium chemistry, and feedback (Linzer et al., 2024). They produce synthetic maps of dust emission (using proxies for heating and column), LyC photon energy density, ionization parameter, and derived emission measures (EM). Key findings include:

• Exponential decline of FUV/LyC fields with height; scale heights set by cluster distribution and attenuation.
• FUV intensity at the midplane $$J_{\rm FUV}(z=0)=6.6\times10^{-14}{\rm\,erg\,cm^{-3}}$$ (solar neighborhood).
• PDFs of normalized FUV intensity show lognormal diffuse peaks plus power-law cluster tails.
• Scaling relations for FUV attenuation versus $n_{\rm H}$ and for LyC mean free path in ionized gas ($\ell_{LyC}\sim100$ pc).
• Validation of analytic/slab models for FUV and LyC distribution against full ART results. These results constrain the reliability of subgrid recipes for dust heating, IR emission, and LyC escape in galaxy-scale and cosmological simulations.

---

TIGRESS shearing-box models have become a reference framework for simulating multiphase, dynamically evolving, feedback-regulated galactic ISM patches, uniquely combining SN-driven turbulence, star formation, large-scale dynamo, radiative transfer, and direct calculation of synthetic observables for comparison with galactic and extragalactic datasets (Kim et al., 2019, Kim et al., 2020, Linzer et al., 2024).