This issue was drafted by Claude (Anthropic) on behalf of Markus Röllig, based on a joint analysis session; the numbers below were produced and verified in that session.
Summary
SimLine currently accepts at most two collision partners, which must be pH₂ (+ optionally oH₂), blended internally with the thermal ortho/para ratio into a single rate set applied to the one density column of the shell input. For PDR applications this is a real accuracy limit: the surface layers that emit [CII], [OI], [CI] and CH⁺ are predominantly atomic, so excitation there is driven by H (and partly e⁻), not H₂. We have been post-processing KOSMA-τ models with SimLine and can now quantify exactly where the single-partner model breaks — and where it provably cannot be worked around from outside the code.
Current behavior
collrate.f (readspecies):
ncspecies = 1 → one partner block, must contain H₂;ncspecies = 2 → first block must be pH₂, second oH₂ ("First collision partner must be para-H2 !"); combined via the thermal o/p ratio;- anything else → "Cannot yet handle multiple collision partners !" (the yet suggests this was anticipated).
The molecule files in circulation already carry the additional partner blocks as unreachable reference data — e.g. c+_and_h2.lamda contains C⁺+H (Barinovs et al. 2005) and C⁺+e⁻ (Wilson & Bell 2002) after the two H₂ blocks, and cI.lamda / oI.lamda carry H, e⁻, He, H⁺ blocks (LAMDA partner ids 5, 4, 6, 7). So the data side is ready; only the reader and the per-shell rate assembly are missing.
Why an external workaround is not sufficient
Our converter feeds SimLine an "H₂-equivalent" density (n_H₂ + w·n_H). For a two-level ion this is exact if w = k_H/k_H₂(T) (+ electron term): with that correction our SimLine [CII] agrees with an independent escape-probability code (ONION, using H+pH₂+oH₂ throughout) to 0.4% on a test model — where the naive kinematic weight (0.707) had produced a spurious 15% deficit.
But for multi-level species a single scalar density cannot represent partner mixing when the per-transition rate ratios differ. O I is the decisive case: for its two emitting transitions k_H/k_H₂ ≈ 0.7–1.3, but for the ΔJ=2 channel (³P₀↔³P₂) the ratio is 25–90 (H₂ barely drives ΔJ=2; H does). We bracketed the consequence on a KOSMA-τ clump (n_s = 10⁴ cm⁻³, χ = 10⁴, Joblin-type setup) by running the same model in the two pure-collider limits:
| line |
pure-H₂ rates (current) |
pure-H rates |
ONION (mixed partners) |
| [OI] 63 µm |
2.49e-3 |
3.63e-3 |
3.08e-3 |
| [OI] 145 µm |
1.58e-4 |
2.34e-4 |
2.01e-4 |
(erg s⁻¹ cm⁻² sr⁻¹, clump-averaged.) The limits differ by ~46–50%, and the mixed-partner reference sits at their geometric midpoint — exactly as it should. I.e. SimLine's [OI] currently carries a ≈ −20% systematic in PDR conditions, and no input-side trick can remove it.
Proposed extension
- Reader (
readspecies): accept ncspecies > 2; keep the pH₂/oH₂ pair as the first partner (preserving the existing o/p blending), then read additional blocks keyed by the LAMDA partner id (4 = e⁻, 5 = H, 6 = He, 7 = H⁺). Fully backward compatible: existing files with count ≤ 2 behave identically, and the reference blocks already present in circulating files become live.
- Shell input: optional additional density columns after the existing n column (e.g.
n_H2, T, X, vturb, lcorr, v, [n_H, n_e, …]), with absent columns defaulting to zero — old input files remain valid.
- Rate assembly: wherever the collision rate enters (the
crates path), replace k(T)·n by Σₚ kₚ(T)·nₚ per shell. The level-population/ALI machinery is untouched — partners only change the collision matrix.
The o/p pair support shows all the structural patterns (multiple tables, blending) already exist; this generalizes them from a hardcoded pair to a partner list.
Offer
We are happy to contribute this as a PR (we already build and run current master, including mclarge), together with a validation case: the KOSMA-τ model above with per-shell n_H/n_H₂/n_e profiles, where the multi-partner result must land between the two pure limits and can be cross-checked against the ONION reference values quoted above.
🤖 Drafted with Claude Code
Summary
SimLine currently accepts at most two collision partners, which must be pH₂ (+ optionally oH₂), blended internally with the thermal ortho/para ratio into a single rate set applied to the one density column of the shell input. For PDR applications this is a real accuracy limit: the surface layers that emit [CII], [OI], [CI] and CH⁺ are predominantly atomic, so excitation there is driven by H (and partly e⁻), not H₂. We have been post-processing KOSMA-τ models with SimLine and can now quantify exactly where the single-partner model breaks — and where it provably cannot be worked around from outside the code.
Current behavior
collrate.f(readspecies):ncspecies = 1→ one partner block, must contain H₂;ncspecies = 2→ first block must be pH₂, second oH₂ ("First collision partner must be para-H2 !"); combined via the thermal o/p ratio;The molecule files in circulation already carry the additional partner blocks as unreachable reference data — e.g.
c+_and_h2.lamdacontains C⁺+H (Barinovs et al. 2005) and C⁺+e⁻ (Wilson & Bell 2002) after the two H₂ blocks, andcI.lamda/oI.lamdacarry H, e⁻, He, H⁺ blocks (LAMDA partner ids 5, 4, 6, 7). So the data side is ready; only the reader and the per-shell rate assembly are missing.Why an external workaround is not sufficient
Our converter feeds SimLine an "H₂-equivalent" density (n_H₂ + w·n_H). For a two-level ion this is exact if w = k_H/k_H₂(T) (+ electron term): with that correction our SimLine [CII] agrees with an independent escape-probability code (ONION, using H+pH₂+oH₂ throughout) to 0.4% on a test model — where the naive kinematic weight (0.707) had produced a spurious 15% deficit.
But for multi-level species a single scalar density cannot represent partner mixing when the per-transition rate ratios differ. O I is the decisive case: for its two emitting transitions k_H/k_H₂ ≈ 0.7–1.3, but for the ΔJ=2 channel (³P₀↔³P₂) the ratio is 25–90 (H₂ barely drives ΔJ=2; H does). We bracketed the consequence on a KOSMA-τ clump (n_s = 10⁴ cm⁻³, χ = 10⁴, Joblin-type setup) by running the same model in the two pure-collider limits:
(erg s⁻¹ cm⁻² sr⁻¹, clump-averaged.) The limits differ by ~46–50%, and the mixed-partner reference sits at their geometric midpoint — exactly as it should. I.e. SimLine's [OI] currently carries a ≈ −20% systematic in PDR conditions, and no input-side trick can remove it.
Proposed extension
readspecies): acceptncspecies > 2; keep the pH₂/oH₂ pair as the first partner (preserving the existing o/p blending), then read additional blocks keyed by the LAMDA partner id (4 = e⁻, 5 = H, 6 = He, 7 = H⁺). Fully backward compatible: existing files with count ≤ 2 behave identically, and the reference blocks already present in circulating files become live.n_H2, T, X, vturb, lcorr, v, [n_H, n_e, …]), with absent columns defaulting to zero — old input files remain valid.cratespath), replace k(T)·n by Σₚ kₚ(T)·nₚ per shell. The level-population/ALI machinery is untouched — partners only change the collision matrix.The o/p pair support shows all the structural patterns (multiple tables, blending) already exist; this generalizes them from a hardcoded pair to a partner list.
Offer
We are happy to contribute this as a PR (we already build and run current master, including
mclarge), together with a validation case: the KOSMA-τ model above with per-shell n_H/n_H₂/n_e profiles, where the multi-partner result must land between the two pure limits and can be cross-checked against the ONION reference values quoted above.🤖 Drafted with Claude Code