Gas Temperature Update

Overview

The changetemperature subroutine updates the gas temperature in each PDR grid point once the level populations have converged. It enforces thermal balance by iteratively adjusting the gas temperature until the total heating and cooling rates are equal within a prescribed tolerance.

The routine operates point-by-point over the full PDR domain and is parallelised using OpenMP.

Its main responsibilities are:

  • Compute reaction rates at the current gas temperature.

  • Compute all relevant heating processes.

  • Compare total heating and total cooling rates.

  • Adjust the gas temperature using a hybrid strategy:

    • Exponential temperature stepping (±30%)

    • Binary chop (bisection) once a valid temperature bracket is found

  • Detect convergence and enforce physical temperature bounds.

Thermal Balance Condition

Thermal equilibrium is defined as:

\[\Gamma_{\mathrm{tot}} = \Lambda_{\mathrm{tot}}\]

where:

  • \(\Gamma_{\mathrm{tot}}\) is the total heating rate

  • \(\Lambda_{\mathrm{tot}}\) is the total cooling rate

The imbalance is quantified by:

Mean Imbalance

\[F_{\mathrm{mean}} = \Gamma_{\mathrm{tot}} - \Lambda_{\mathrm{tot}}\]

Relative Imbalance

\[F_{\mathrm{ratio}} = \frac{2\,|F_{\mathrm{mean}}|} {|\Gamma_{\mathrm{tot}} + \Lambda_{\mathrm{tot}}|}\]

Thermal balance is considered achieved when:

\[F_{\mathrm{ratio}} \le F_{\mathrm{crit}}\]

Heating and Cooling

Reaction Rates

For each PDR point, reaction rates are computed using the current gas and dust temperatures:

  • Gas temperature \(T_{\mathrm{gas}}\)

  • Dust temperature \(T_{\mathrm{dust}}\)

  • Local radiation field and column densities along all rays

These rates are passed to the heating routines.

Heating Processes

The routine computes up to 12 distinct heating channels, including (but not limited to):

  • Photoelectric heating

  • Cosmic-ray heating

  • Chemical heating

  • H₂ formation heating

  • Gas–dust thermal coupling

The total heating rate is stored as:

\[\Gamma_{\mathrm{tot}} = \texttt{heating(12)}\]

Cooling Processes

Cooling rates are assumed to have already been computed (e.g. via LVG line cooling) and stored as:

\[\Lambda_{\mathrm{tot}} = \sum_i \Lambda_i\]

Temperature Update Strategy

The temperature update proceeds only after level populations have converged (level_conv = .true.).

Initial Bracketing Phase

During the first thermal iteration:

  • If \(F_{\mathrm{mean}} > 0\) (net heating):

    \[T_{\mathrm{new}} = 1.3\,T_{\mathrm{old}}\]

    and the lower bound is updated:

    \[T_{\mathrm{low}} = T_{\mathrm{old}}\]

  • If \(F_{\mathrm{mean}} < 0\) (net cooling):

    \[T_{\mathrm{new}} = 0.7\,T_{\mathrm{old}}\]

    and the upper bound is updated:

    \[T_{\mathrm{high}} = T_{\mathrm{old}}\]

This phase searches for a temperature interval that brackets the thermal equilibrium solution.

Binary Chop Phase

Once heating and cooling change sign between iterations, the algorithm switches to binary chopping (bisection):

  • If heating dominates:

    \[T_{\mathrm{low}} = T_{\mathrm{current}}, \quad T_{\mathrm{new}} = \frac{T_{\mathrm{current}} + T_{\mathrm{high}}}{2}\]

  • If cooling dominates:

    \[T_{\mathrm{high}} = T_{\mathrm{current}}, \quad T_{\mathrm{new}} = \frac{T_{\mathrm{current}} + T_{\mathrm{low}}}{2}\]

From this point onward, only bisection is used.

Convergence Criteria

A PDR point is marked as fully converged if any of the following are met:

  1. Thermal balance criterion:

    \[F_{\mathrm{ratio}} \le F_{\mathrm{crit}}\]

  2. Temperature stagnation without balance:

    \[|T_{\mathrm{new}} - T_{\mathrm{old}}| \le T_{\mathrm{diff}} \quad \text{and} \quad F_{\mathrm{ratio}} > F_{\mathrm{crit}}\]

  3. Physical bounds are reached consistently:

    • \(T \le T_{\min}\) while cooling dominates

    • \(T \ge T_{\max}\) while heating dominates

Temperature Bounds

The gas temperature is strictly constrained to lie within:

\[T_{\min} \le T_{\mathrm{gas}} \le T_{\max}\]

If convergence is reached at a boundary:

  • The temperature is clamped.

  • Further iteration may be forced if level populations must be recomputed.

Parallelisation

The routine is OpenMP-parallelised over PDR points:

  • Each PDR point is independent.

  • Shared arrays are read-only or explicitly indexed by point.

  • Private variables include:

    • Reaction rates

    • Heating arrays

    • Temporary temperature variables

This ensures scalability for large 3D models.

Summary

The changetemperature routine provides a robust and physically controlled method to enforce local thermal balance in PDR models after level population convergence. Its hybrid stepping–bisection strategy ensures:

  • Rapid initial exploration of temperature space

  • Guaranteed convergence once a valid bracket is found

  • Stability in extreme heating or cooling regimes

  • Strict adherence to physical temperature limits

This routine forms the thermal closure of the coupled chemistry–radiative transfer–thermal balance problem in 3D-PDR.