© The Authors 2025

1. Introduction

Protostellar outflows are essential consequences when young stars are forming (e.g. Frank et al. 2014; Bally 2016; Pascucci et al. 2023). They transport away a substantial amount of excess angular momentum from the protostellar system and thus prevent the spinning up of the protostar. In order to fully understand the nature of their gas excitation and the launching mechanism it is important to accurately measure the physical conditions of these outflows, that is the electron density n_e, electron temperature T_e, and ionisation fraction x_e. These key parameters can be constrained observationally by analysing the emitted forbidden emission lines (FELs), which are attributed to the presence of shocks. Protostellar outflows are typically traced by FELs, such as [S II]λλ4068, 4076, [O I]λλ6300, 6363, [S II]λλ6716, 6731, and [N II]λλ6548, 6583 indicating low-excitation gas conditions (x_e ≲ 0.3, n_e ∼ 10² − 10⁶ cm⁻³, T_e ∼ 5000 − 20 000 K).

In this context different methods are used to derive the gas parameters. First, shock models can be used to compare observed FELs with predictions (e.g. Hollenbach & McKee 1989; Hartigan et al. 1994). Although this approach is the natural choice as shocks cause the observed emission, this approach introduces a strong model dependency and degeneracy. The geometry of the shock, the shock type (C-shock, J-shock, radiative shock) and other parameters (shock velocities, preshock densities, magnetic fields) substantially influence the resulting FEL spectra. Second, individual line ratios that are sensitive to one or more gas parameters can be used as diagnostic tools (e.g. Osterbrock & Ferland 2006; Natta et al. 2014; Fang et al. 2018; Nisini et al. 2024). In this framework the extinction and the gas parameters are typically constrained separately and/or consecutively, depending on the available line ratios. For example, electron densities and temperatures are typically inferred from the [S II]λ6716/[S II]λ6731 or [O I]λλ (6300+6363)/[O I]λ5577 line ratios, respectively. The ionisation fraction is then derived from further line ratios (e.g. [N II]λ6583/[O I]λ6300) under the assumption of a certain ionisation model with n_e and T_e as previously determined parameters. This method is in essence a multi-step approach. Often only a range for n_e, T_e, and x_e can be constrained and some available lines remain unused. A third approach is to find the best excitation model, whereby again no detailed shock physics is involved, that describes simultaneously all observed line ratios that are attributed to the same emitting gas (e.g. Hartigan & Morse 2007; Giannini et al. 2015, 2019). Mathematically, this represents a non-linear minimisation problem. Model grid calculations are usually preferred over model fitting, which in turn can be computationally very expensive already for a handful of parameters and a grid with moderate resolution. As a caveat, not all appearing lines can be blindly combined as they can come from different gas phases. Fourth, a compromise between computationally expensive grid calculations and an in-depth analysis of a few diagnostic diagrams is the BE99 method (Bacciotti & Eislöffel 1999), where only one diagnostic diagram is utilised to derive n_e, T_e, and x_e (e.g. Nisini et al. 2005; Podio et al. 2006, 2009; Ray et al. 2007; Coffey et al. 2008; Maurri et al. 2014). The BE99 method utilises six forbidden emission lines ([O I]λλ6300, 6363, [S II]λλ6716, 6731, and [N II]λλ6548, 6583) in the visual part of the electromagnetic spectrum and works in the framework of collisional ionisation equilibrium, where the atoms and ions are assumed to exist mostly in their ground state. In a nutshell, the BE99 diagnostic diagram is constructed as follows. From the observed six FELs, three line ratios with their errors are calculated: O I/N II, O I/S II, and N II/S II (O I := [O I]λλ6300+6363, N II := [N II]λλ6548+6583, and S II := [S II]λλ6716+6731). These three ratios are then compared with theoretical predictions in the (x_e, T_e) parameter space. If the observed line ratio with error margins is consistent with a given set of values in the parameter space (e.g. {x_eⁱ, T_e^j}) this point is accepted as a potential representation of the data. For all three line ratios these points are collected and effectively appear as three extended curves, that is stripes, in the (x_e, T_e) parameter space. The location, where all stripes overlap, gives the determined gas parameters.

All of these methods are based on time-independent models, that is they assume that the emitting gas is in equilibrium. The hydrogen recombination time τ_rec and dynamical timescale τ_dyn are usually given as justification for their usage. Since protostellar outflows are not stationary, as proper motions studies show, these and potentially other timescales play a crucial role.

Furthermore, with new spectroscopic instruments such as VLT¹/XShooter, VLT/UVES², or VLT/MUSE³ the visual part of the spectrum can be studied, and valuable emission lines that are connected to the BE99 lines in the blue or near-infrared are accessible as well. In turn, the BE99 method can be extended in a natural way (i.e. by including more emission line ratios as stripes in the (x_e, T_e) parameter space) potentially providing more information about the present gas conditions.

In this study we explore the idea of extending the BE99 method by utilising additional line ratios in the spectral range of 3500 − 11 000 Å and considering higher ionisation fractions (x_e > 0.3). Furthermore, we numerically solve a network of ionisation reactions, on which the BE99 ionisation balance is based, in the time domain to investigate how the non-equilibrium situation affects the BE99 method. We study the timescales of reaching the equilibrium to inspect how well the equilibrium assumption is justified when applied to protostellar outflows. We test our propositions on the low-excitation outflow of Par Lup 3−4 and the high-excitation 244−440 Proplyd in Orion.

2. The BE99 method in the time domain

One of the main assumption of the BE99 method is that the reaction equilibrium in the involved atoms and ions is reached. In this section we investigate the situation, when the BE99 ionisation balance is not yet reached. We use a simple network of ionisation reactions that lead to the equilibrium ionisation balance of the BE99 method.

2.1. Deriving the time-dependent reaction rates

Four elements are involved in the BE99 technique. These are nitrogen (N), oxygen (O), sulphur (S), and hydrogen (H). Their relevant ionisation states are denoted as Y^k ∈ {N⁰, N⁺, O⁰, O⁺, S⁺, H⁰, H⁺}. Their number densities are denoted as n(Y^k). The electron density is n_e. The total number densities of an element Y, here introduced as n_tot(Y), are given by the sum of all relevant ionisation states. In a low-excitation medium we expect sulphur to be mostly singly ionised. Nitrogen, oxygen, and hydrogen coexist in their neutral and singly ionised form simultaneously. Thus, their total number densities and ionisation fractions are given by

$$\begin{array}{c}\hfill n_{\mathrm{tot}}(Y)=n(Y^0)+n(Y^{+}),\end{array}$$(1)

and

$$\begin{array}{c}\hfill x_Y=n(Y^{+})/[n(Y^0)+n(Y^{+})].\end{array}$$(2)

Three atomic processes are considered in the BE99 technique: (i) direct and inverse charge exchange, (ii) collisional ionisation, and (iii) radiative and dielectronic recombination. Photoionisation is neglected. The BE99 method assumes ionisation equilibrium in the involved reactions. Without this assumption a time-dependent version of the BE99 method can be formulated from the underlying set of reactions, that is

$$\begin{array}{cc}& {\mathrm{N}}^0+{\mathrm{H}}^{+}{\displaystyle \underset{{\delta }_{\mathrm{N}}^{\prime }}{\stackrel{{\delta }_{\mathrm{N}}}{\rightleftharpoons }}}{\mathrm{N}}^{+}+{\mathrm{H}}^0(\mathrm{charge}\mathrm{exchange}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\mathrm{N}}^0+{\mathrm{e}}^{-}\stackrel{C_{\mathrm{N}}}{⟶}{\mathrm{N}}^{+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(4)

$$\begin{array}{cc}& {\mathrm{N}}^{+}+{\mathrm{e}}^{-}\stackrel{{\alpha }_{\mathrm{N}}}{⟶}{\mathrm{N}}^0(\mathrm{recombination}),\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\mathrm{O}}^0+{\mathrm{H}}^{+}{\displaystyle \underset{{\delta }_{\mathrm{O}}^{\prime }}{\stackrel{{\delta }_{\mathrm{O}}}{\rightleftharpoons }}}{\mathrm{O}}^{+}+{\mathrm{H}}^0(\mathrm{charge}\mathrm{exchange}),\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\mathrm{O}}^0+{\mathrm{e}}^{-}\stackrel{C_{\mathrm{O}}}{⟶}{\mathrm{O}}^{+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\mathrm{O}}^{+}+{\mathrm{e}}^{-}\stackrel{{\alpha }_{\mathrm{O}}}{⟶}{\mathrm{O}}^0(\mathrm{recombination}),\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\mathrm{H}}^0+{\mathrm{e}}^{-}\stackrel{C_{\mathrm{H}}}{⟶}{\mathrm{H}}^{+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\mathrm{H}}^{+}+{\mathrm{e}}^{-}\stackrel{{\alpha }_{\mathrm{H}}}{⟶}{\mathrm{H}}^0(\mathrm{recombination}).\hfill \end{array}$$(10)

The respective reaction rates are indicated above or below the arrow of each reaction. They are functions of the T_e (Appendix A).

From the above reactions (Eqs. (3)–(10)) the full time-dependent system consisting of seven coupled ordinary differential equations – one for each involved particle N⁰, N⁺, O⁰, O⁺, H⁰, H⁺, e⁻ – can be derived straightforwardly (Appendix A). In the BE99 method, no assumptions for the hydrogen reactions are involved. In turn, without invoking any knowledge about additional reactions with hydrogen, the hydrogen ionisation fraction x_e remains a free parameter in the BE99 method. In the excitation model, however, the ionisation fraction in the equilibrium will be a function of the electron temperature T_e and not a free parameter. This model then provides the opportunity to study the BE99 method in the time domain.

2.2. The initial conditions

The stated system of ordinary differential equations can be integrated numerically by specifying eight initial number densities n_init(N⁰), n_init(N⁺), n_init(O⁰), n_init(O⁺), n_init(H⁰), n_init(H⁺), n_init(S⁺), and $n_{\mathrm{e}}^{\mathrm{init}}$. To keep our model reasonably simple, we assume the gas in isothermal condition. Since we assume constant element abundances and global charge neutrality over the whole integration, only four initial values instead of eight fully determine the reactions. In order to invoke physically relevant initial conditions we chose $n_{\mathrm{e}}^{\mathrm{init}}$, $x_{\mathrm{e}}^{\mathrm{init}}$, $x_{\mathrm{N}}^{\mathrm{init}}$, and $x_{\mathrm{O}}^{\mathrm{init}}$ as a set of independent initial parameters. Element abundances are taken from Asplund et al. (2009). Their stated values are N/H = 6.76 × 10⁻⁵, O/H = 4.90 × 10⁻⁴, and S/H = 1.32 × 10⁻⁵. We note that these values differ from the values used in the BE99 paper (N/H = 1.1 × 10⁻⁴, O/H = 6.0 × 10⁻⁴, and S/H = 1.6 × 10⁻⁵). Different assumed element abundances and updated atomic data affect the results of any diagnostic method.

The transformations between the two sets of initial conditions can be derived straightforwardly. We find

$$\begin{array}{c}\hfill n({\mathrm{H}}^{+})=n({\mathrm{H}}^0)·x_{\mathrm{e}}/(1-x_{\mathrm{e}}).\end{array}$$(11)

The element abundances are invoked by

$$\begin{array}{cc}& n({\mathrm{N}}^0)+n({\mathrm{N}}^{+})=\mathrm{N}/\mathrm{H}·(n({\mathrm{H}}^0)+n({\mathrm{H}}^{+})),\hfill \end{array}$$(12)

$$\begin{array}{cc}& n({\mathrm{O}}^0)+n({\mathrm{O}}^{+})=\mathrm{O}/\mathrm{H}·(n({\mathrm{H}}^0)+n({\mathrm{H}}^{+})),\hfill \end{array}$$(13)

$$\begin{array}{cc}& n({\mathrm{S}}^{+})=\mathrm{S}/\mathrm{H}·(n({\mathrm{H}}^0)+n({\mathrm{H}}^{+})).\hfill \end{array}$$(14)

Furthermore, the electron density is fixed by global charge neutrality

$$\begin{array}{c}\hfill n({\mathrm{e}}^{-})=n({\mathrm{H}}^{+})+n({\mathrm{O}}^{+})+n({\mathrm{N}}^{+})+n({\mathrm{S}}^{+}).\end{array}$$(15)

With these equations we can derive

$$\begin{array}{c}\hfill n({\mathrm{H}}^0)=\frac{(1-x_{\mathrm{e}})n_{\mathrm{e}}}{x_{\mathrm{e}}+\frac{\mathrm{S}}{\mathrm{H}}+\left(\frac{\mathrm{N}}{\mathrm{H}}\right)x_{\mathrm{N}}+\left(\frac{\mathrm{O}}{\mathrm{H}}\right)x_{\mathrm{O}}}·\end{array}$$(16)

Thus, from Eqs. (11)–(16) the initial conditions can be transformed one-to-one by specifying $n_{\mathrm{e}}^{\mathrm{init}}$, $x_{\mathrm{e}}^{\mathrm{init}}$, $x_{\mathrm{N}}^{\mathrm{init}}$, and $x_{\mathrm{O}}^{\mathrm{init}}$. The ionisation parameters are constrained to $x_{\mathrm{e}}^{\mathrm{init}}$, $x_{\mathrm{N}}^{\mathrm{init}}$, $x_{\mathrm{O}}^{\mathrm{init}}\in [0,1]$, with 0 representing a neutral medium and 1 being fully ionised (in that element). We integrate the full system of the coupled ordinary differential equations numerically until the reaction equilibrium is reached using the Runge-Kutta-solver implemented in Python (Virtanen et al. 2020). An example of a numerical solution with T_e = const = 10 500 K is depicted in Fig. 1. Figure 1 shows how the three ratios H⁺/H⁰, O⁺/O⁰, and N⁺/N⁰ change with time until the reaction equilibrium is reached after about 2 × 10¹⁰ s ≈ 634 years. The ratios H⁺/H⁰ and O⁺/O⁰ show a very similar evolution over the shown integration time, whereas the N⁺/N⁰ ratio behaves differently (an extra maximum at about 10⁸ s). This behaviour is a consequence of the dominant charge exchange reactions between oxygen and hydrogen in the reaction network (see following Section 2.3). Both ratios H⁺/H⁰ and O⁺/O⁰ are locked in – with O⁺/O⁰ ≲ H⁺/H⁰ – shortly after the start of the reactions.

Fig. 1. Time resolved numerical solution of the reaction rates for the ionisation fractions of hydrogen, oxygen, and nitrogen (Eqs. (3)–(10)). The initial conditions at t0 = 0 were $n_{\mathrm{e}}^{\mathrm{init}}=2000{\mathrm{cm}}^{-3}$, $x_{\mathrm{e}}^{\mathrm{init}}=0.1$, $x_{\mathrm{O}}^{\mathrm{init}}=0.1$, and $x_{\mathrm{N}}^{\mathrm{init}}=0.1$. They have already evolved to the values at the starting point of this plot. The reaction equilibrium is reached after τeq ≈ 2 × 1010 s as indicated as pink shaded area. We fix Te = const = 10 500 K during the simulation.

2.3. The reaction equilibrium

An interesting situation occurs when a stable equilibrium is reached after some time t. For the equilibrium we require

$$\begin{array}{c}\hfill \mathrm{d}[Y^k]/\mathrm{d}t=0,\mathrm{for}\mathrm{all}{Y}^k\mathrm{and}\mathrm{for}t>{\tau }_{\mathrm{eq}}.\end{array}$$(17)

By invoking these constraints to the full system of ODEs, we can derive the ionisation balance in the equilibrium. They are

$$\begin{array}{cc}\hfill \frac{{\mathrm{N}}^{+}}{{\mathrm{N}}^0}{|}_{\mathrm{eq}}& =\frac{{\delta }_{\mathrm{N}}·{[{\mathrm{H}}^{+}]}_{\mathrm{eq}}+C_{\mathrm{N}}·{[{\mathrm{e}}^{-}]}_{\mathrm{eq}}}{{\delta }_{\mathrm{N}}^{\prime }·{[{\mathrm{H}}^0]}_{\mathrm{eq}}+{\alpha }_{\mathrm{N}}·{[{\mathrm{e}}^{-}]}_{\mathrm{eq}}},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill \frac{{\mathrm{O}}^{+}}{{\mathrm{O}}^0}{|}_{\mathrm{eq}}& =\frac{{\delta }_{\mathrm{O}}·{[{\mathrm{H}}^{+}]}_{\mathrm{eq}}+C_{\mathrm{O}}·{[{\mathrm{e}}^{-}]}_{\mathrm{eq}}}{{\delta }_{\mathrm{O}}^{\prime }·{[{\mathrm{H}}^0]}_{\mathrm{eq}}+{\alpha }_{\mathrm{O}}·{[{\mathrm{e}}^{-}]}_{\mathrm{eq}}},\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill \frac{{\mathrm{H}}^{+}}{{\mathrm{H}}^0}{|}_{\mathrm{eq}}& =\frac{{\delta }_{\mathrm{O}}^{\prime }·{[{\mathrm{O}}^{+}]}_{\mathrm{eq}}+{\delta }_{\mathrm{N}}^{\prime }·{[{\mathrm{N}}^{+}]}_{\mathrm{eq}}+C_{\mathrm{H}}·{[{\mathrm{e}}^{-}]}_{\mathrm{eq}}}{{\delta }_{\mathrm{O}}·{[{\mathrm{O}}^0]}_{\mathrm{eq}}+{\delta }_{\mathrm{N}}·{[{\mathrm{N}}^0]}_{\mathrm{eq}}+{\alpha }_{\mathrm{H}}·{[{\mathrm{e}}^{-}]}_{\mathrm{eq}}}·\hfill \end{array}$$(20)

We note that only Eqs. (18) and (19) appear in the BE99 method, since the reactions with hydrogen are not constrained by any assumption. Charge exchange with oxygen dominates for T_e ∼ 5000 − 15 000 K, since ${\delta }_{\mathrm{O}},\delta {\prime }_{\mathrm{O}}\sim {10}^{-9}{\mathrm{cm}}^3{\mathrm{s}}^{-1}$ and C_O, α_O ∼ 10⁻¹² − 10⁻¹³ cm³ s⁻¹. Thus, the ionisation states of oxygen and hydrogen are tightly locked via

$$\begin{array}{c}\hfill \frac{{\mathrm{O}}^{+}}{{\mathrm{O}}^0}{|}_{\mathrm{eq}}\approx \left(\frac{{\delta }_{\mathrm{O}}}{{\delta }_{\mathrm{O}}^{\prime }}\right)·\left(\frac{{\mathrm{H}}^{+}}{{\mathrm{H}}^0}\right){|}_{\mathrm{eq}}\approx \left(\frac{8}{9}\right)·\left(\frac{{\mathrm{H}}^{+}}{{\mathrm{H}}^0}\right){|}_{\mathrm{eq}}\mathrm{for}{T}_{\mathrm{e}}\sim {10}^4\mathrm{K}\end{array}$$(21)

as shown in Osterbrock & Ferland (2006). Furthermore, we can assume n_e ≈ n(H⁺) (neutral plasma condition). Thus we have

$$\begin{array}{c}\hfill n({\mathrm{H}}^0)\approx (1-x_{\mathrm{e}})·n_{\mathrm{e}}/x_{\mathrm{e}}\end{array}$$(22)

and the ionisation equilibrium for nitrogen (Eq. (18)) and oxygen (Eq. (19)) can be approximated by

$$\begin{array}{cc}\hfill \frac{{\mathrm{N}}^{+}}{{\mathrm{N}}^0}{|}_{\mathrm{eq}}& \approx \frac{({\delta }_{\mathrm{N}}+C_{\mathrm{N}})·x_{\mathrm{e}}^{\mathrm{eq}}}{({\alpha }_{\mathrm{N}}-{\delta }_{\mathrm{N}}^{\prime })x_{\mathrm{e}}^{\mathrm{eq}}+{\delta }_{\mathrm{N}}^{\prime }},\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill \frac{{\mathrm{O}}^{+}}{{\mathrm{O}}^0}{|}_{\mathrm{eq}}& \approx \frac{({\delta }_{\mathrm{O}}+C_{\mathrm{O}})·x_{\mathrm{e}}^{\mathrm{eq}}}{({\alpha }_{\mathrm{O}}-{\delta }_{\mathrm{O}}^{\prime })x_{\mathrm{e}}^{\mathrm{eq}}+{\delta }_{\mathrm{O}}^{\prime }}\approx \left(\frac{{\delta }_{\mathrm{O}}}{{\delta }_{\mathrm{O}}^{\prime }}\right)·\frac{x_{\mathrm{e}}^{\mathrm{eq}}}{1-x_{\mathrm{e}}^{\mathrm{eq}}}\hfill \end{array}$$(24)

as demonstrated by Bacciotti & Eislöffel (1999). The above Eqs. (23) and (24) are in full agreement with our time-dependent simulations displayed in Fig. 1.

2.4. The ionisation fraction in the equilibrium

Extensive numerical experiments confirm that the ionisation fraction in the equilibrium ($x_{\mathrm{e}}^{\mathrm{eq}}$) is a function of only T_e. In detail, we integrated the time-dependent reaction model shown in Fig. 1 until the equilibrium state for a fine temperature range (T_e = 8000 − 16 500 K) and grid of initial conditions. We plot $x_{\mathrm{e}}^{\mathrm{eq}}$ in Fig. 2 and fit a four parametric function to the curve

$$\begin{array}{c}\hfill x_{\mathrm{e}}^{\mathrm{eq}}(T_{\mathrm{e}})=A·{(1+\exp [-B·(T_{\mathrm{e}}-T_0)])}^C\end{array}$$(25)

Fig. 2. The hydrogen ionisation fraction of the BE99 reactions in the equilibrium $x_{\mathrm{e}}^{\mathrm{eq}}$ as a function of the electron temperature Te.

with A = 0.9785, B = 8.778 × 10⁻⁴, T₀ = 12 000 K, C = −2.1762. The function $x_{\mathrm{e}}^{\mathrm{eq}}(T_{\mathrm{e}})$ looks qualitatively similar to the hydrogen ionisation fraction derived from the Saha-Equation. However, both functions are quantitatively very different (see Appendix B). In turn, the Saha-Equation does not properly describe the BE99 gas conditions in the equilibrium.

In principle, Eq. (25) offers the opportunity to measure the ionisation fraction directly from the electron temperature, i.e. without the need to construct a BE99 diagram. The line ratio [O I]λλ(6300+6363)/[O I]λ5577 would be sufficient for that purpose as it is sensitive to T_e. However, this way of determining the ionisation fraction requires that the proposed ionisation reactions of our model are the dominant ones and the temperature stays constant with time. Using Eq. (25) for determining the ionisation fraction would therefore introduce a model dependency, which the BE99 technique deliberately avoids.

2.5. Timescales: τ_rec, τ_eq, τ_dyn

The BE99 method is based on the assumption that the gas under investigation is in equilibrium. The hydrogen recombination time

$$\begin{array}{c}\hfill {\tau }_{\mathrm{rec}}\sim 1/(n_{\mathrm{e}}{\alpha }_{\mathrm{H}}(T_{\mathrm{e}}))\approx {10}^5\mathrm{yr}/n_{\mathrm{e}}[{\mathrm{cm}}^{-3}]\end{array}$$(26)

with α_H as hydrogen recombination rate coefficient is usually used to justify this assumption (Bacciotti et al. 1995; Osterbrock & Ferland 2006; Hartigan & Morse 2007). The line of reasoning is that the ionisation fractions of O⁺/O⁰ and H⁺/H⁰ are coupled through charge exchange and therefore ultimately with $x_{\mathrm{e}}^{\mathrm{eq}}$. In this regard Osterbrock & Ferland (2006) noted that ‘deviations from ionisation equilibrium are ordinarily damped out in times of this [τ_rec] order of magnitude’. However, oxygen and nitrogen react on different timescales as compared to hydrogen and in addition the hydrogen recombination time is derived for a reaction system involving only hydrogen, i.e. no coupling with other reactants. Thus, the hydrogen recombination time may not provide an appropriate justification, if the equilibrium is reached or not. In the following, we therefore solve the ODEs numerically to study the how fast the equilibrium is reached.

We define that the equilibrium is reached (τ_eq), when the ionisation fractions of hydrogen, nitrogen, and oxygen stay constant with time and deviate only by 5% from the ionisation equilibrium (Eqs. (23) and (24)). With this definition we numerically found a crude estimate:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{eq}}\approx 10\mathrm{yr}\frac{x_{\mathrm{e}}^{\mathrm{init}}}{n_{\mathrm{e}}^{\mathrm{init}}[{\mathrm{cm}}^{-3}]}\exp \left(\frac{145528}{T_{\mathrm{e}}[\mathrm{K}]}\right)·\end{array}$$(27)

We can compare the equilibrium times scale (Eq. (27)) with the hydrogen recombination time (Eq. (26)) by plugging in realistic numbers for T_e, n_e, and x_e. We base our illustration on the Bacciotti & Eislöffel (1999) measurements of the HH 34 jet. The brightest knot J of the HH 34 jet beam is located at about θ = 20″ away from the source (Figs. 6a–f therein). An estimate of τ_eq can be done by taking the values of T_e, n_e, and x_e at the apex of knot J, for which we read T_e ≈ 8000 K, n_e ≈ 500 cm⁻³, and x_e ≈ 0.05. With Eq. (27) we therefore derive τ_eq ≈ 80 000 yr. The hydrogen recombination time τ_rec derived from Eq. (26) gives a value of about 200 yr, which is more than two orders of magnitude smaller than the equilibrium timescale. From that comparison the medium before and behind knot J of the HH 34 jet has not enough time to reach the reaction equilibrium and in turn the BE99 technique would not be applicable here. In comparison, the dynamical timescale (dynamical shock age) of the jet can be estimated via

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dyn}}=\theta ·D/v_{\mathrm{t}}\approx 4.74\times \frac{\theta [″]·D[\mathrm{pc}]}{v_{\mathrm{t}}[{\mathrm{km}\mathrm{s}}^{-1}]}\mathrm{yr}\end{array}$$(28)

with D as distance to the target, v_t as transverse velocity, and θ as projected angle on the plane of sky. For the HH34 knot J the dynamical timescale is of the order of the hydrogen recombination time, i.e. τ_dyn ≈ 200 yr for D = 480 pc and v_t = 220 km s⁻¹. In order to resolve the stated timescale-problem (τ_eq ≫ τ_rec, τ_dyn), we introduce a new timescale in the following Section 2.6 that is more suitable for describing the state of the system.

2.6. Resolving the BE99 method in the time domain

The critical question for the applicability of the BE99 method may actually not be if the gas is in reaction equilibrium, which may never be the case in any outflow region, but instead how far off the gas is from the equilibrium.

In order to answer this question we followed the numerical solution of our model in the time domain. At each point in time we had full knowledge of the true gas conditions and could therefore calculate synthetic spectra in the non-equilibrium and equilibrium state. These (model-dependent) synthetic spectra can be used as input for the (model-independent) BE99 method. This situation corresponds to the observers perspective who does not know if the equilibrium is reached. Applying the BE99 technique out of equilibrium will give measurements of the electron density, temperature, and ionisation fraction ($n_{\mathrm{e}}^{\mathrm{BE}99}$, $T_{\mathrm{e}}^{\mathrm{BE}99}$, $x_{\mathrm{e}}^{\mathrm{BE}99}$) that do not represent the true gas conditions ($n_{\mathrm{e}}^{\mathrm{true}}$, $T_{\mathrm{e}}^{\mathrm{true}}$, $x_{\mathrm{e}}^{\mathrm{true}}$). As the system evolves in time from a non-equilibrium to the equilibrium state, the BE99 method will provide higher accuracy in these three gas parameters (i.e. it converges). In order to quantify this transition phase we introduced the timescale τ_BE. We define τ_BE as the time at which the BE99 method provides a measurement of the ionisation fraction within 5% of the true ionisation fraction, which is comparable to a typical measurement error of excellent observational data:

$$\begin{array}{c}\hfill \parallel 1-x_{\mathrm{e}}^{\mathrm{BE}99}/x_{\mathrm{e}}^{\mathrm{true}}\parallel <0.05\mathrm{for}t>{\tau }_{\mathrm{BE}}.\end{array}$$(29)

We visualise the convergence of the BE99 method in the time domain in Fig. 3: already within τ_BE ≈ 10⁸ s (3 yr) the BE99 method measures the true ionisation within the 5% error margin. In comparison, the hydrogen recombination time derived from Eq. (26) is about τ_rec ≈ 50 yr, implying τ_BE < τ_rec ≪ τ_eq. In conclusion, neither τ_eq nor τ_rec are indicating the relevant timescale for the applicability of the BE99 method – it is rather τ_BE. We tested this hypothesis by integrating 192 models with different input parameters (Fig. 4). All the tested models have a BE99 equilibrium time smaller than 10⁹ s ≈ 32 yr. This result supports the notion that the BE99 equilibrium is reached fast enough – even faster that the hydrogen recombination time.

Fig. 3. Time resolved BE99 method for the model of Fig. 1. The yellow horizontal stripe indicates the region, where the BE99 method measures an ionisation fraction to 5% accuracy with respect to the ionisation fraction. Within τBE ≈ 108 s the BE99 method converges – long before the reaction equilibrium ($x_{\mathrm{e}}^{\mathrm{eq}}=0.0337$) is reached (τBE ≪ τeq).

Fig. 4. Convergence of the BE99 method in the time domain. In total 192 models with different initial values have been integrated. Parameters: Te [K] = 8000, 9000, 10 000; $n_{\mathrm{e}}^{\mathrm{init}}[{\mathrm{cm}}^{-3}]=1000,2000,5000,10000$; $x_{\mathrm{e}}^{\mathrm{init}}=0.05,0.1,0.2,0.3$; $x_{\mathrm{N}}^{\mathrm{init}}=0.01,0.1$; $x_{\mathrm{O}}^{\mathrm{init}}=0.01,0.1$. Each black curve represents one model, where the BE99 method (assuming equilibrium) has been applied following the gas through the non-equilibrium state. For all these models the BE99 method measures the true gas conditions up to 5% accuracy (indicated as yellow horizontal stripe) after ∼109 s.

3. The BE99 method with more line ratios (BE99e)

In the classical BE99 method three line ratios are involved: OI/NII, SII/NII, and OI/SII. All six associated emission lines are located close to each other in the visible part of the spectrum and thus extinction correction is not crucial. In addition, the three stripes of line ratios most probably always meet in one common region in the (x_e, T_e) space when applied to observations of protostellar outflows. From the BE99 diagram it is not possible to see if the equilibrium is reached or if the extinction determination is correct or the BE99 assumptions are sufficiently fulfilled. However, other line ratios of the involved species in the BE99 ionisation equilibrium in the UVB and NIR part of the spectrum can in principle be used to extend the BE99 diagram in a natural way. Figure 5 shows a synthetic emission line spectrum of the model in the equilibrium with the strongest emission lines of [O I], [N II], [S II], [N I], and [O II] in the spectral range of 3500 − 11 000 Å. These lines are also comparably bright in protostellar outflows and thus can be included in the diagnostics. Other transitions such as the [N II]λ5754 line are typically too faint. From the inspection of this spectrum many extensions of the BE99 diagram, hereby labelled as BE99e, are possible. We note that no additional assumption or equation has been invoked by transitioning from BE99 to BE99e. A selection of line ratios for constructing BE99e diagrams are depicted in Fig. 6. Including any combination of the depicted line ratios in the BE99e diagram is straight forward and can give more information about the state of the gas. Observers can in principle choose according to the dataset available which stripes to use for determining the ionisation fraction. They do not necessarily have to rely on the six classic BE99 lines.

Fig. 5. Predicted BE99e emission lines. Black vertical lines are the synthetic line emissivities normalised to the [O I]λ6300 line. Red vertical lines are the reddened line emissivities assuming the extinction law of Cardelli et al. (1989) with AV = 1.0 and RV = 3.1. The equilibrium conditions are Te = 10 500 K, $n_{\mathrm{e}}^{\mathrm{eq}}=672.65{\mathrm{cm}}^{-3}$, $x_{\mathrm{e}}^{\mathrm{eq}}=0.0336$.

Fig. 6. Possible combinations of lines in the BE99e method. Note: In the (Te, xe)-plane the choice of the numerator/denominator for a given line ratio is irrelevant, e.g. OI/NII and NII/OI represent the same stripe.

In order to demonstrate how the new line ratios can be exploited we present predictions of how the various new stripes in the (x_e, T_e) diagram look like in the time domain together with the influence of extinction for the model (Appendix C). As examples we present four BE99e diagrams at five time steps (10⁴ s, 10⁶ s, 10⁸ s, 10¹⁰ s, 10¹² s) and three extinction values (A_V = 0.0, 0.3, 0.6) in Figs. C.1–C.5. We note that emission lines towards shorter wavelengths are more affected by extinction than emission lines in the near-infrared. In order to illustrate this effect we plot in Fig. 5 the extinction correction factor for A_V = 1.0 and R_V = 3.1. The intrinsic and observed line fluxes at wavelength λ are denoted by F_int(λ) and F_obs(λ). We assume that

$$\begin{array}{c}\hfill F_{\mathrm{obs}}(\lambda )=F_{\mathrm{int}}(\lambda )·{10}^{-0.4A_{V}c(\lambda ,R_V)},\end{array}$$(30)

with c(λ, R_V) as normalised reddening curve and A_V as visual extinction (e.g. Osterbrock & Ferland 2006). From Cardelli et al. (1989) we adopt the expression for c(λ, R_V).

3.1. Determining the extinction using the BE99e diagram

Applying the extended BE99e method on observational data requires that the line emissivities are first corrected for extinction via e.g. the Balmer decrement, Paschen decrement, Millers method, or line ratios of iron (e.g. [Fe II]1.32 μm/1.63 μm). All the mentioned methods come with flaws and uncertainties, which the observer tries to minimise. If the gas is in equilibrium and the extinction is corrected for properly all stripes overlap in one region in the (x_e, T_e) space. However, the critical case occurs, when not all stripes are crossing in one region. In such a situation one should try to find a visual extinction $A_V^{\mathrm{BE}}$ in Eq. (30), for which the stripes overlap. If such an extinction value exists the gas is close enough to the equilibrium. This procedure should be repeated until a satisfactory $A_V^{\mathrm{true}}=A_V+A_V^{\mathrm{BE}}$ is found. Also the electron density and electron temperature can be treated in the same way. By changing them one can in principle check if the stripes will cross in one location.

3.2. Insights from the extra stripes in the (x_e, T_e) diagram

We briefly discuss the global behaviour of the new stripes in the (x_e, T_e) diagram depending on the non-equilibrium/equilibrium state and the extinction. Our main findings are:

• In the classic BE99 diagram all three stripes always meet in one location in the (x_e, T_e) space. From the diagram alone it cannot be seen if the equilibrium condition is fulfilled or if extinction is correctly determined. In our models, the classic BE99 method is robust in measuring the electron density and electron temperature (i.e. the influence of extinction and non-equilibrium conditions are marginal).

• A plethora of BE99e diagrams can be constructed with the same underlying assumptions of the original BE99 paper. In contrast, these new diagrams are sensitive to the equilibrium state and the extinction. In the equilibrium state all stripes cross in one location in the (x_e, T_e) space measuring the true gas conditions, if the spectra are correctly dereddened. If the equilibrium is not reached or the spectra have not been properly dereddened the stripes spread over the (x_e, T_e) space not crossing in one location. Instead they potentially cross in many locations which are not representing the true gas conditions. Some stripes of line ratios are more sensitive to non-equilibrium or extinction influence than others. For example the ratios of OI/[S II]λλ4068+4076, NII/[OII]λλ3726+3729, OI/[N I]λλ5198+5200 are very sensitive to both influences. Less sensitive line ratios are e.g. SII/[S II]λλ4068+4076, SII/[N I]λλ5198+5200, SII/[N I]λλ10397+10398.

• With the additional BE99e diagrams a determination of x_e and T_e is possible even if not all classical BE99 lines are available. Line ratios can be exchanged according to the lines available.

• Some emission line ratios represent almost the same stripes in the (x_e, T_e) space. For example OI/[N I]λλ10397+10398 overlaps almost perfectly with OI/[N I]λλ10407+10408.

4. Extension to higher ionisation fractions (BE99e+)

The BE99 and BE99e methods work for a low-ionisation medium (x_e < 0.3). We can, however, extend it to higher ionisation degrees. We indicate this extension by BE99e+. In a high-excitation medium we expect to have several ionisation states of oxygen, nitrogen, sulphur, and hydrogen. We therefore include the following species

$$\begin{array}{c}\hfill Y^k=({\mathrm{O}}^0,{\mathrm{O}}^{+},{\mathrm{O}}^{2+},{\mathrm{N}}^0,{\mathrm{N}}^{+},{\mathrm{N}}^{2+},{\mathrm{S}}^{+},{\mathrm{S}}^{2+},{\mathrm{S}}^{3+},{\mathrm{H}}^0,{\mathrm{H}}^{+}).\end{array}$$(31)

Additional to the reactions rates Eqs. (3)–(10) there are

$$\begin{array}{cc}\hfill {\mathrm{O}}^{+}+{\mathrm{e}}^{-}& \stackrel{C_{{\mathrm{O}}^{+}}}{⟶}{\mathrm{O}}^{2+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill {\mathrm{O}}^{2+}+{\mathrm{e}}^{-}& \stackrel{{\alpha }_{{\mathrm{O}}^{+}}}{⟶}{\mathrm{O}}^{+}(\mathrm{recombination}),\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill {\mathrm{O}}^{+}+{\mathrm{H}}^{+}& {\displaystyle \underset{{\delta }_{{\mathrm{O}}^{+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{O}}^{+}}}{\rightleftharpoons }}}{\mathrm{O}}^{2+}+{\mathrm{H}}^0(\mathrm{charge}\mathrm{exchange}),\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill {\mathrm{N}}^{+}+{\mathrm{e}}^{-}& \stackrel{C_{{\mathrm{N}}^{+}}}{⟶}{\mathrm{N}}^{2+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill {\mathrm{N}}^{2+}+{\mathrm{e}}^{-}& \stackrel{{\alpha }_{{\mathrm{N}}^{+}}}{⟶}{\mathrm{N}}^{+}(\mathrm{recombination}),\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill {\mathrm{N}}^{+}+{\mathrm{H}}^{+}& {\displaystyle \underset{{\delta }_{{\mathrm{N}}^{+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{N}}^{+}}}{\rightleftharpoons }}}{\mathrm{N}}^{2+}+{\mathrm{H}}^0(\mathrm{charge}\mathrm{exchange}),\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{+}+{\mathrm{e}}^{-}& \stackrel{C_{{\mathrm{S}}^{+}}}{⟶}{\mathrm{S}}^{2+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{2+}+{\mathrm{e}}^{-}& \stackrel{{\alpha }_{{\mathrm{S}}^{+}}}{⟶}{\mathrm{S}}^{+}(\mathrm{recombination}),\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{2+}+{\mathrm{e}}^{-}& \stackrel{C_{{\mathrm{S}}^{2+}}}{⟶}{\mathrm{S}}^{3+}+2{\mathrm{e}}^{-}(\mathrm{coll}.\mathrm{ionisation}),\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{3+}+{\mathrm{e}}^{-}& \stackrel{{\alpha }_{{\mathrm{S}}^{2+}}}{⟶}{\mathrm{S}}^{2+}(\mathrm{recombination}),\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{+}+{\mathrm{H}}^{+}& {\displaystyle \underset{{\delta }_{{\mathrm{S}}^{+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{S}}^{+}}}{\rightleftharpoons }}}{\mathrm{S}}^{2+}+{\mathrm{H}}^0(\mathrm{charge}\mathrm{exchange}),\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{2+}+{\mathrm{H}}^{+}& {\displaystyle \underset{{\delta }_{{\mathrm{S}}^{2+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{S}}^{2+}}}{\rightleftharpoons }}}{\mathrm{S}}^{3+}+{\mathrm{H}}^0(\mathrm{charge}\mathrm{exchange}).\hfill \end{array}$$(43)

From these equations it is straight forward to derive a set of ordinary differential equations describing the time-dependent reaction network (Appendix A).

By invoking the equilibrium condition (Eq. (17)) to the full system of ODEs, we can derive the new ionisation balance in the equilibrium for BE99e+. They are

$$\begin{array}{cc}\hfill \frac{{\mathrm{O}}^{+}}{{\mathrm{O}}^0}& =\frac{({\delta }_{\mathrm{O}}+C_{\mathrm{O}})·x_{\mathrm{e}}}{({\alpha }_{\mathrm{O}}-{\delta }_{\mathrm{O}}^{\prime })x_{\mathrm{e}}+{\delta }_{\mathrm{O}}^{\prime }},\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill \frac{{\mathrm{O}}^{2+}}{{\mathrm{O}}^{+}}& =\frac{({\delta }_{{\mathrm{O}}^{+}}+C_{{\mathrm{O}}^{+}})·x_{\mathrm{e}}}{({\alpha }_{{\mathrm{O}}^{+}}-{\delta }_{{\mathrm{O}}^{+}}^{\prime })x_{\mathrm{e}}+{\delta }_{{\mathrm{O}}^{+}}^{\prime }},\hfill \end{array}$$(45)

$$\begin{array}{cc}\hfill \frac{{\mathrm{N}}^{+}}{{\mathrm{N}}^0}& =\frac{({\delta }_{\mathrm{N}}+C_{\mathrm{N}})·x_{\mathrm{e}}}{({\alpha }_{\mathrm{N}}-{\delta }_{\mathrm{N}}^{\prime })x_{\mathrm{e}}+{\delta }_{\mathrm{N}}^{\prime }},\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill \frac{{\mathrm{N}}^{2+}}{{\mathrm{N}}^{+}}& =\frac{({\delta }_{{\mathrm{N}}^{+}}+C_{{\mathrm{N}}^{+}})·x_{\mathrm{e}}}{({\alpha }_{{\mathrm{N}}^{+}}-{\delta }_{{\mathrm{N}}^{+}}^{\prime })x_{\mathrm{e}}+{\delta }_{{\mathrm{N}}^{+}}^{\prime }},\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill \frac{{\mathrm{S}}^{2+}}{{\mathrm{S}}^{+}}& =\frac{({\delta }_{{\mathrm{S}}^{+}}+C_{{\mathrm{S}}^{+}})·x_{\mathrm{e}}}{({\alpha }_{{\mathrm{S}}^{+}}-{\delta }_{{\mathrm{S}}^{+}}^{\prime })x_{\mathrm{e}}+{\delta }_{{\mathrm{S}}^{+}}^{\prime }},\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill \frac{{\mathrm{S}}^{3+}}{{\mathrm{S}}^{+}}& =\frac{({\delta }_{{\mathrm{S}}^{2+}}+C_{{\mathrm{S}}^{2+}})·x_{\mathrm{e}}}{({\alpha }_{{\mathrm{S}}^{2+}}-{\delta }_{{\mathrm{S}}^{2+}}^{\prime })x_{\mathrm{e}}+{\delta }_{{\mathrm{S}}^{2+}}^{\prime }}·\hfill \end{array}$$(49)

With the additional species of O²⁺ and S²⁺ we can include several new line ratios such as SIII = [S III]λ9068 + [S III]λ9530 or OIII = [O III]λ4958 + [O III]λ5007. With that plenty new line ratios in the BE99 diagram can be added straightforwardly including line ratios such as OI/SIII or NII/OIII. The presented equations for the reaction network in the equilibrium represent an extended excitation model such as the one described by Giannini et al. (e.g. 2015).

5. Extension including other species (BE99e++)

In principle, the described method may be extended to include species other than oxygen, nitrogen, and sulphur. In order to determine the ionisation fraction at least two ionisation states of the same species have to be included, that is for example Fe⁺/Fe²⁺ or Ne⁺/Ne²⁺. The ionisation balance Y⁺ⁱ/Y^+i + 1 can be found by writing down all relevant physical reactions as described in Section 4 and solve the reaction network for the equilibrium. The main challenge here may be to identify the relevant reactions. The magnitude of the reaction rates, however, indicate the dominant reactions. The ionisation balance Y⁺ⁱ/Y^+i + 1 is ideally described as a function of x_e and T_e:

$$\begin{array}{c}\hfill Y^{+i}/Y^{+i+1}{|}_{\mathrm{eq}}=f(x_{\mathrm{e}},T_{\mathrm{e}}).\end{array}$$(50)

Line ratios of Y⁺ⁱ and Y^+i + 1 (BE99e++ method) can be included as stripes in the (x_e, T_e) diagram straightforwardly.

5.1. Non-thermal contributions

The outlined theory of extending the BE99 method in various directions is limited by the assumption that all lines are thermally excited. Other effects such as (a) photodissociation of OH molecules in the disk surface layers (e.g. Gorti et al. 2011; Rigliaco et al. 2013), (b) photoionisation (e.g. Bally & Reipurth 2001; Riaz & Whelan 2015; Méndez-Delgado et al. 2021), (c) fluorescent excitation (Bautista 1999), or (d) pumping (Ferland et al. 2012; Nemer et al. 2020) can contribute to the observed line emission and preclude the applicability of the BE99 method. Line ratios such as [S II]λ4068/[O I]λ6300 or [O I]λ6300/[O I]λ5577 can give further clues on the excitation mechanism. For example for the low-velocity component of protostellar outflows of T Tauri stars, Natta et al. (2014) and Fang et al. (2018) found that the observed line ratios are consistent with thermal excitation. For the high-velocity component typically shock models are used to explain the observed line ratios. However, each source is different and thus the assumption of thermal excitation must be justified as well.

5.2. Other caveats and limitations

Element abundances can substantially differ from solar abundances in protostellar outflows and jets due to, for example, grain sputtering or different evolutionary pathways of the observed material. Deviant element abundances change the position of the stripes in any BE99 diagram and therefore the determined gas parameters. Furthermore, collisions with other species such as hydrogen atoms are not yet included in the excitation model. The lack of atomic data prevents this so far. Line emissivities from the excitation model are expected to be affected to an unknown amount.

Another limiting aspect is the critical density (n_crit) of the utilised FELs. Although its definition is not consistent in the literature, its meaning can be described as follows: For a given energy level and temperature the critical density refers to the density of the collision partners, at which the de-excitation rate is equal to the radiative decay rate from that level. In the high-density limit (here: n_e ≫ n_crit) that implies that the line luminosity per atom or ion coming from that level does not increase with the density of the collision partners anymore, that is it is not sensitive to density variations. The frequent collisional excitation and de-excitation essentially set up a Boltzmann population and the emerging emission lines are quenched. A common misconception is that these quenched lines disappear, which they in fact do not. However, line ratios, where quenched lines are involved, may not be used as diagnostics. The utilised emission lines have distinct critical densities n_crit in the range of n_crit = 10³ − 10⁸ cm⁻³ (see Table 3 in Giannini et al. 2019). In turn, line ratios can be affected by quenching or are just not sensitive to certain gas conditions. For instance, the optical sulphur lines, which are used for determining the electron densities, are only sensitive around n_e ∼ 10⁴ cm⁻³. In the high-density limit these lines are quenched and not sensitive to the density anymore. In the BE99 method, this line ratio can in principle be substituted by the [O II]λλ3726/3729 line ratio. However, also the blue [O II] lines have similar critical densities as compared to the optical sulphur lines. In turn, regions with high electron densities, for example close to the driving source, may not be accessible able to any BE99 method.

6. Applications of the extended BE99 methods

In the following section we apply the excitation model described in Section 4 to observational data of two targets with very distinct excitation conditions.

6.1. Application to the (high-excitation) Proplyd 244–440

In order to demonstrate the applicability of additional stripes in the (x_e, T_e) diagram we take a look at recently published MUSE data of the outflow associated with the 244−440 Orion proplyd in Kirwan et al. (2023). The determination of the ionisation fraction via the BE99 method was not possible as stated therein. No [N II] line fluxes are specified in the paper (see Table 1 therein) indicating that they are either not detected or other unspecified problems prevent their usage. However, other line fluxes of oxygen and sulphur have been measured making a determination of the ionisation fraction possible with the outlined extension of the BE99e+ diagrams. For example, the extracted dereddened line fluxes of [O I]λλ6300+6363, [S II]λλ6716+6731, [O II]λ7320, and [O II]λ7330 at knot E3 of Proplyd 244−440 are stated in Table 1 therein. Since no line errors are specified we assumed that the corresponding errors are 2%. Higher/lower errors would thicken/thin down the stripes in the BE99e+ diagram. Even without the [N II] lines at hand it is now possible to generate a BE99e+ diagram based on the line ratios of OI/SII, OI/OII, SII/OII. We present the BE99e+ diagram based on the stated line ratios in Fig. 7. From Fig. 7 we see that: The line ratios do cross in a common region in the (x_e, T_e) space as predicted by our model. We measure a electron temperature of T_e = 18 800 K ± 1000 K and an ionisation fraction of about x_e = 0.58 ± 0.05. This indicates that this region is in a high-excitation state, which is consistent with the expected gas conditions.

Fig. 7. An application of the additional line ratios (BE99e+) to determine the gas parameters of the outflow region associated with knot E3 of 244−440 Orion proplyd (Kirwan et al. 2023).

6.2. Application to the (low-excitation) Par Lup 3−4 outflow

Par Lup 3−4 (16^h 08^m 51.44^s, −39^d 05^m 30.5^s) is a well-studied protostellar outflow of the Lupus III cloud. A detailed source description can be found in Bacciotti et al. (2011). Its molecular outflow was investigated in Santamaría-Miranda et al. (2020). Observations with the medium resolution slit spectrograph X-Shooter revealed that the protostellar outflow associated with Par Lup 3−4 (HH 600) displays low-excitation (Whelan et al. 2014).

6.2.1. Observations and data reduction

In our analysis we make use of the archival X-Shooter observations of Par Lup 3−4 (programme ID: 085.C-0238(A), date of observations: April 7, 2010). X-Shooter has three arms covering a spectral range of about 3000 − 5590 Å (UVB), 5590 − 10 240 Å (VIS), and 10 240 − 24 800 Å (NIR). In the case of Par Lup 3−4 the slit was aligned with the outflows PA of about 130°. The slit widths were $1\stackrel{″}{.}0$, $0\stackrel{″}{.}9$, and $0\stackrel{″}{.}9$ implying a resolving power 5100 (58 km s⁻¹), 8800 (34 km s⁻¹), and 5600 (53 km s⁻¹) for the UVB, VIS, and NIR arm, respectively. The seeing during the observations was about $0\stackrel{″}{.}9-1\stackrel{″}{.}2$. Exposure times were 4 × 900 s in each arm. The X-Shooter spectra of Par Lup 3−4 were analysed in Bacciotti et al. (2011), Giannini et al. (2013), and Whelan et al. (2014) – hereafter B11, G13, and W14, respectively.

We re-reduced the archival X-Shooter data for Par Lup 3−4 using the up-to-date Esoreflex environment (v. 3.6.1, Freudling et al. 2013) to obtain the 1D and 2D spectra. The data reduction includes bias subtraction, flat fielding, wavelength calibration, and optimal spectrum extraction. We checked the 1D and 2D spectra for consistent relative flux calibration in the UVB, VIS, and NIR arms – no relative flux adjustment was necessary as in the data reduction of W14. The pixel scales in the 2D data are $0\stackrel{″}{.}16{\mathrm{pixel}}^{-1}$ in the UVB and VIS arms, and $0\stackrel{″}{.}21{\mathrm{pixel}}^{-1}$ in the NIR arm. The spectral samplings are 0.2 Å in the UVB and VIS arms, and 0.6 Å in the NIR arm. Cosmics were removed with the Lacosmic routine implemented in Python (van Dokkum 2001). No photospheric correction was needed as the photospheric features are marginally present. We adopt the systemic velocity of v_LSR = 5.5 km s⁻¹ stated in B11. We extracted the position-velocity (PV) diagrams for the relevant emission lines by fitting and removing a 2D Gaussian around the emission line region. We present the PV diagrams in Fig. D.1. In total 26 line regions were inspected for our analysis. These include the six optical BE99 lines plus 20 further line regions lying in the spectral range XShooter.

6.2.2. Extracting line fluxes

We are interested in the physical conditions of the outflow close to the driving source, that is the innermost one arcsec region. This region displays a clear detection in most of the BE99 lines. Following W14 and G13 we restrict our analysis to the region of $-0\stackrel{″}{.}5$ to $+0\stackrel{″}{.}5$ with respect to the position of the driving source. We decided to extract flux values in this region confined spectrally by −80 km s⁻¹ and +30 km s⁻¹. Initially, we extracted fluxes in five consecutive pixel rows in spatial direction at the driving source. However, we noticed no substantial change in the relative line fluxes.

Figure 8 displays the flux box for the [O I]λ6300 line of Par Lup 3−4 as an example for the flux extraction. Since we only need relative line fluxes for our analysis we normalise all extracted fluxes to the line flux of the bright [O I]λ6300 line (i.e. we set its flux to F([O I]λ6300) = 100). The extracted normalised flux values with their 1σ-errors or upper limits are stated in Table E.1.

Fig. 8. The position-velocity diagrams for the [O I]λ6300 line of Par Lup 3−4. Line fluxes have been extracted within the sketched box.

For our line diagnostics we adopt A_V = 0.0, which was derived from the analysis of the near-infrared iron lines in the Par Lup 3−4 outflow (G13). We did not utilise line ratios of e.g. [S II]λ4068/[S II]λ10320, [S II]λ4076/[S II]λ10336 to derive the extinction (Miller 1968), since no extinction value is consistent with the observed Miller line ratios. This is puzzling and implies that either some of the [S II] emission in the UVB arm is overluminous or the [S II] emission in the NIR arm is underluminous.

6.2.3. NIR line fluxes in W14

For a consistency check we also list the W14 flux values in Table E.1. Comparing our extracted relative flux values with the ones presented in W14 we see a very good agreement within the errorbars for the VIS and UVB arms. However, the NIR flux values in W14 are systematically too high by a factor of four as compared to our values. W14 investigated two outflows, namely Par Lup 3−4 and ESO Hα 574. In their data reduction for ESO Hα 574 W14 found that the NIR arm fluxes are too low at about 25% when comparing the continuum in the 1D spectra of the VIS and UVB arms. For Par Lup 3−4 this was not the case, which is in full agreement with our analysis. We suspect that W14 nonetheless applied the correction factor found for ESO Hα 574 also on the NIR fluxes of Par Lup 3−4, which fully explains the stated inconsistency. The usage of different pipeline versions of Esoreflex can explain small differences in the relative line fluxes as displayed in Table E.1, but not the systematically too high NIR fluxes in W14.

6.2.4. Timescale of the outflow

We can estimate the dynamical scale and hydrogen recombination time of one of the two outflow lobes of Par Lup 3−4. (We do not know the equilibrium timescale for the outflow, but as our model calculations show this timescale is not decisive). The distance to Par Lup 3−4 is D ∼ 155 pc (Santamaría-Miranda et al. 2020). We assume a jet velocity of v_t = 170 km s⁻¹ and a jet length of $0\stackrel{″}{.}5$ (W14). With that we get a dynamical timescale of about 2 yr. The hydrogen recombination time (Eq. (26)) is about τ_rec ∼ 1 − 10 yr assuming n_e ∼ 10⁴ − 10⁵ cm⁻³. Thus the gas of the Par Lup 3−4 outflows (both lobes assuming symmetry) should have had enough time to stay close enough to the equilibrium.

6.2.5. Physical conditions of Par Lup 3−4

In order to investigate the physical conditions in the Par Lup 3−4 outflow we first check if thermal excitation can explain the observed line ratios (e.g. Osterbrock & Ferland 2006; Giannini et al. 2015, 2019; Nisini et al. 2024). Figures 9a–c show three exemplary diagnostic diagrams as a function of T_e, n_e, and x_e for the observed line ratios of [S II]λ6716/[S II]λ6731, [S II]λ4068/[O I]λ6300, and [N II]λ6583/[O I]λ6300. The observed line ratio of [S II]λ6716/[S II]λ6731 ∼0.5 is close to the high-density limit indicating high electron densities, that is n_e ≳ 5 × 10³ cm⁻³. Such high electron densities are indeed expected to be measured close to the driving source of the outflow (Fang et al. 2018). The ionisation fraction can be constrained from the [N II]λ6583/[O I]λ6300 line ratio. For high electron densities the observed [N II]λ6583/[O I]λ6300 line ratio indicates x_e ≲ 0.03 for T_e = 5000 − 25 000 K. The observed [S II]λ4068/[O I]λ6300 line ratio is consistent with these gas parameters. In conclusion, all observed line ratios are consistent with the thermal excitation.

Fig. 9. Diagnostic diagrams for the observed line fluxes of Par Lup 3−4. The red horizontal stripes indicate the measured line flux ratios including their 1σ error margins. The black and grey curves are theoretical line ratios derived from our underlying excitation model. They represent reproductions of Fig. A.2. in Giannini et al. (2019).

W14 determined the gas parameters n_e, T_e, and x_e via the classical BE99 method. The stated values in Table 1 therein are: n_e = 4300 ± 200 cm⁻³, T_e = 31 000 ± 500 K, and x_e = 0.003 ± 0.0002 for the blue lobe and n_e = 4600 ± 300 cm⁻³, T_e = 24 000 ± 600 K, and x_e = 0.0063 ± 0.0002 for the red lobe. These values are consistent with the values that we derive for the full flux box when applying the classical BE99 method. When adding additional line ratios (BE99e or BE99e+), such as OI/[N I]λλ5198+5200, SII/[N I]λλ5198+5200, and NII/[N I]λλ5198+5200, however, not all stripes cross in one location in the (x_e, T_e) diagram. We extensively tried to bring the stripes together by applying different extinction values, but they never cross in one location. Thus, the gas parameters derived from the classical BE99 method are not in good agreement with the observations and any BE99 method fails.

6.2.6. Fitting line ratios to excitation model

Apart from the diagnostic diagrams and the BE99e and BE99e+ methods we could take the whole spectrum and fit an excitation model to it, that is we utilise all observed lines at once and not separately. For such a multi-line approach with n observed line fluxes (F_i^obs with i = 1, …, n) following Hartigan & Morse (2007) we calculated an extensive grid of excitation models with the three parameters n_e, T_e, and x_e. Our model grid spans: n_e = 2000 cm⁻³ − 200000 cm⁻³ spaced linearly in 2000 grid points (Δn_e = 100 cm⁻³), T_e = 5000 K − 30 000 K spaced linearly in 250 grid points (ΔT_e = 100 K), and x_e = 0.001 − 0.200 spaced linearly in 200 grid points (Δx_e = 0.001). In total 2000 × 250 × 200 = 10⁸ models with all relevant synthetic line fluxes (F_i^model) are calculated to find the best fit to the observed line fluxes. The best model in our model grid is then found by minimising the residuals defined by

$$\begin{array}{c}\hfill C={\displaystyle \sum _{k=1}^{n(n-1)/2}}{\left(\frac{r_k-m_k}{{\sigma }_k}\right)}^2,\end{array}$$(51)

with r_k = log(F_i^obs/F_j^obs), m_k = log(F_i^model/F_j^model), and σ_k as propagated line ratio errors (i, j = 1, …, n with i < j). We also tested non-zero extinction values to check for better fits, but the residuals were always better for A_V = 0.0.

We compare the observed line fluxes to the modelled line fluxes in Fig. 10. It is remarkable how good the excitation model fits to the Par Lup 3−4 data. In fact, all observed lines are consistent with the best excitation model within 10%, and all the non-detections are consistent with the model as well. The multi-line approach indicates the following excitation condition in the Par Lup 3−4 outflow: n_e ∼ 45000 cm⁻³ − 53000 cm⁻³, T_e = 7600 K − 8000 K, and x_e ∼ 0.027 − 0.036. These values substantially differ from the ones derived from the classical BE99 method. The electron density is about 6 times higher than the one derived from BE99. The electron temperature is about three times lower than from BE99 and the ionisation fraction is about four times higher than from BE99. We think that the multi-line approach by Hartigan & Morse (2007) best describes the observations of Par Lup 3−4 due to the following reasons: 1. The [O I]λ5577 line is not detected towards Par Lup 3−4 but should be, if the classical BE99 conditions prevail, 2. The [S II]λ4068 are substantially underpredicted by the BE99 model. The BE99 method may not give proper results for the Par Lup 3−4 outflows because of the high electron densities at the driving source (quenching effects).

Fig. 10. A comparison of the observed emission line ratios (black bars) with the best model (yellow bars) found by the multi-line analysis and the W14 model found via the classical BE99 method (blue bars). [O I]λ6300 is reproduced by both models, the plot therefore zooms in on the weaker lines. We see an overall agreement of the observed data with the two models – especially with the optical BE99 lines. The [O I]λ5577, [N I]λ5198, and [N I]λ5200 emission lines, however, are substantially overpredicted and the [S II]λ4068 and[S II]λ4076 lines are heavily underpredicted in W14. The near-infrared lines of sulphur and nitrogen are well reproduced.

7. Conclusions

We have numerically integrated a set of ionisation reactions in the time domain which represent the BE99 reaction network. We used this model to generate synthetic emission line spectra of the atomic/ionic species that are involved in the BE99 equilibrium (NI, NII, SII, OI, OII). From these spectra we were able to study how the system evolves from a non-equilibrium to the equilibrium state. Our analysis demonstrates that the equilibrium of the BE99 reactions is reached on timescales shorter than the hydrogen recombination time. This, however, is not critical for the application of the original BE99 method, since the BE99 method converges faster than the hydrogen recombination time.

We propose a natural extension of the BE99 method, that is based on the same assumptions as in the original BE99 paper (BE99e), by including many more line ratios as stripes in the (x_e, T_e) diagram in the spectral range of about 3700 − 10 400 Å. If the BE99 conditions are fulfilled all the new stripes in the (x_e, T_e) space meet in one region indicating the true gas parameters. If, however, the system is out of the equilibrium, or extinction has not been corrected properly, or the determined gas parameters are not correct the stripes do not overlap in one but in many regions. Furthermore we put forward an extension of the BE99 method towards higher ionisation fractions (x_e > 0.3, BE99e+). Additional line ratios including lines of for instance [S III] and [O III] can be exploited. Other species such as Fe⁺/Fe²⁺ or Ne⁺/Ne²⁺ can in principle be included straightforwardly.

We tested the outlined ideas on two examples of distinct targets, that is the high-excitation Proplyd 244−440 and the low-excitation outflow Par Lup 3−4. For Proplyd 244−440 we were able to derive the ionisation fraction even without the [N II]λλ6548+6583 lines. Line ratios of [S II]λλ6716+6731 and [O I]λλ6300+6363 with [O II]λ7320 or [O II]λ7330 can be used as substitutes.

We extensively analysed the emission lines of the protostellar outflow Par Lup 3−4. When including the NI = [N I]λ5198+[N I]λ5200 line ratio in the BE99e+ diagram we see that not all stripes overlap in one region in the (x_e, T_e). Our analysis shows that neither extinction nor non-equilibrium effects are likely responsible for that notion. A comparison with a multi-line approach revealed that the most probable explanation is that any BE99 method underestimates the electron density due to quenching of the optical sulphur lines in this case. The application of the BE99 method close to the driving source, where densities are likely very high, that is of the order of the critical density of the sulphur lines, has to be checked carefully by e.g. inspecting other line ratios.

Data availability

The figures published via Zenodo can be downloaded under https://zenodo.org/records/14228515

Acknowledgments

The authors thank the anonymous referee for helping to improve the quality of the paper. This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the YTTHACA Project (Young stars at Tübingen, Tautenburg, Hamburg & ESO – A Coordinated Analysis) under the programme ID EI 409/20-1.

References

Appendix A: Reaction rates and atomic data

A.1. Charge transfer rates

We use the charge transfer rates ${\delta }_{\mathrm{N}}\prime$, δ_N⁺′, δ_S⁺′, and ${\delta }_{{\mathrm{S}}^{2+}}\prime$ stated in Kingdon & Ferland (1996). Since the stated rate for oxygen in therein is only valid for temperatures of T_e = 10 − 10 000 K we use a simple power law fit to the data of Stancil et al. (1999). For T_e = 500 − 50 000 K we find the best fit to be

$$\begin{array}{c}\hfill {\delta }_{\mathrm{O}}^{\prime }(T_{\mathrm{e}})=5.262\times {10}^{-11}{T}_{\mathrm{e}}{[\mathrm{K}]}^{0.3986}[{\mathrm{cm}}^3{\mathrm{s}}^{-1}].\end{array}$$(A.1)

The charge transfer rates δ_Y^k′ and δ_Y^k are connected by the Milne relation (e.g. Osterbrock & Ferland 2006). For example for oxygen the charge transfer reaction is

$$\begin{array}{c}\hfill {\mathrm{O}}^0{(}^3{\mathrm{P}}_{2,1,0})+{\mathrm{H}}^{+}\underset{{\delta }_{\mathrm{O}}^{\prime }}{\stackrel{{\delta }_{\mathrm{O}}}{\rightleftharpoons }}{\mathrm{O}}^{+}{(}^4{\mathrm{S}}_{3/2}^0)+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}).\end{array}$$(A.2)

Here, we assume that only ground state transitions as indicated in brackets in the above Eq. A.2 are relevant. The function δ_O(T_e) can then be calculated via

$$\begin{array}{c}\hfill \frac{{\delta }_{\mathrm{O}}(T_{\mathrm{e}})}{{\delta }_{\mathrm{O}}^{\prime }(T_{\mathrm{e}})}=\frac{g{(}^4{\mathrm{S}}_{3/2}^{\mathrm{o}})\times g{(}^4{\mathrm{S}}_{1/2}^{\mathrm{o}})}{g{(}^3{\mathrm{P}}_0)+g{(}^3{\mathrm{P}}_1)+g{(}^3{\mathrm{P}}_2)}\exp (-\frac{\mathrm{\Delta }E}{k_{\mathrm{B}}{T}_{\mathrm{e}}}),\end{array}$$(A.3)

where ΔE = χ(O⁰)−χ(H⁰) is the energy difference of the two involved ionisation potentials. In the same way we calculate δ_N, δ_O⁺, δ_N⁺, δ_S⁺, and ${\delta }_{{\mathrm{S}}^{2+}}$ from the below reactions

$$\begin{array}{cc}\hfill {\mathrm{O}}^0{(}^3{\mathrm{P}}_{2,1,0})+{\mathrm{H}}^{+}& \underset{{\delta }_{\mathrm{O}}^{\prime }}{\stackrel{{\delta }_{\mathrm{O}}}{\rightleftharpoons }}{\mathrm{O}}^{+}{(}^4{\mathrm{S}}_{3/2}^0)+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}),\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill {\mathrm{N}}^0{(}^4{\mathrm{S}}_{3/2}^0)+{\mathrm{H}}^{+}& \underset{{\delta }_{\mathrm{N}}^{\prime }}{\stackrel{{\delta }_{\mathrm{N}}}{\rightleftharpoons }}{\mathrm{N}}^{+}{(}^3{\mathrm{P}}_{0,1,2})+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\mathrm{O}}^{+}{(}^4{\mathrm{S}}_{3/2}^0)+{\mathrm{H}}^{+}& \underset{{\delta }_{{\mathrm{O}}^{+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{O}}^{+}}}{\rightleftharpoons }}{\mathrm{O}}^{2+}{(}^3{\mathrm{P}}_{2,1,0})+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}),\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\mathrm{N}}^{+}{(}^3{\mathrm{P}}_{0,1,2})+{\mathrm{H}}^{+}& \underset{{\delta }_{{\mathrm{N}}^{+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{N}}^{+}}}{\rightleftharpoons }}{\mathrm{N}}^{2+}{(}^2{\mathrm{P}}_{1/2,3/2}^0)+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}),\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{+}{(}^4{\mathrm{S}}_{3/2}^0)+{\mathrm{H}}^{+}& \underset{{\delta }_{{\mathrm{S}}^{+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{S}}^{+}}}{\rightleftharpoons }}{\mathrm{S}}^{2+}{(}^3{\mathrm{P}}_{2,1,0})+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}),\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill {\mathrm{S}}^{2+}{(}^3{\mathrm{P}}_{0,1,2})+{\mathrm{H}}^{+}& \underset{{\delta }_{{\mathrm{S}}^{2+}}^{\prime }}{\stackrel{{\delta }_{{\mathrm{S}}^{2+}}}{\rightleftharpoons }}{\mathrm{S}}^{3+}{(}^2{\mathrm{P}}_{1/2,3/2}^0)+{\mathrm{H}}^0{(}^2{\mathrm{S}}_{1/2}).\hfill \end{array}$$(A.9)

We take the ionisation potentials from Appendix D of Draine (2011)

A.2. Collisional ionisation rates

The reaction rates for collisional ionisation are taken from Landini & Monsignori Fossi (1990). In the range of T_e = 10⁴ − 10⁸ K (hot plasma) the analytic expression is used

$$\begin{array}{cc}\hfill C(T_{\mathrm{e}})& =A_{\mathrm{col}}{T}_{\mathrm{e}}{[\mathrm{K}]}^{1/2}{(1+0.1\frac{T_{\mathrm{e}}[\mathrm{K}]}{T_{\mathrm{col}}})}^{-1}\hfill \\ \hfill & \exp (-\frac{T_{\mathrm{col}}}{T_{\mathrm{e}}[\mathrm{K}]})[{\mathrm{cm}}^3{\mathrm{s}}^{-1}].\hfill \end{array}$$(A.10)

The exponential factor exp(−T_col/T_e[K]) in the above equation may be read as activation energy barrier for the collisional ionisation with electrons, i.e. roughly the ionisation potential.

For the hydrogen collisional ionisation rates (units: [cm³ s⁻¹]) we take the formula from Dopita & Sutherland (2003)

$$\begin{array}{c}\hfill C_{\mathrm{H}}(T_{\mathrm{e}})=2.5\times {10}^{-10}{T}_{\mathrm{e}}{[\mathrm{K}]}^{1/2}(1+\frac{T_{\mathrm{e}}[\mathrm{K}]}{78945})\exp (-\frac{157890}{T_{\mathrm{e}}[\mathrm{K}]}).\end{array}$$(A.11)

A.3. Recombination rates

Two processes contribute to the total recombination rate α_tot (Landini & Monsignori Fossi 1990). These are the radiative α_R and the dielectronic recombination rates α_D

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{tot}}={\alpha }_{\mathrm{R}}+{\alpha }_{\mathrm{D}}.\end{array}$$(A.12)

In the range of T_e = 10⁴ − 10⁸ K the analytic expressions for α_R are given by Landini & Monsignori Fossi (1990)

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{R}}(T_{\mathrm{e}})=A_{\mathrm{rad}}{T}_4^{-{\mathrm{X}}_{\mathrm{rad}}}[{\mathrm{cm}}^3{\mathrm{s}}^{-1}].\end{array}$$(A.13)

We note that in Landini & Monsignori Fossi (1990) the sign of X_rad is not correct. The radiative recombination rates stated therein are in agreement with Pequignot et al. (1991). We take α_D from Nussbaumer & Storey (1983), which are valid for T_e ∼ 10³ − 6 × 10⁴ K. Explicitly they are given by (units: cm³ s⁻¹)

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{D}}(T_{\mathrm{e}})={10}^{-12}(\frac{a}{T_4}+b+c{T}_4+d{T}_4^2)T_4^{-3/2}\exp (-\frac{f}{T_4}).\end{array}$$(A.14)

For the hydrogen recombination rate (case B) we take the formula presented in Osterbrock & Ferland (2006)

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{H}}(T_{\mathrm{e}})=3.69\times {10}^{-10}{T}_{\mathrm{e}}{[\mathrm{K}]}^{-0.79}[{\mathrm{cm}}^3{\mathrm{s}}^{-1}].\end{array}$$(A.15)

For sulphur we take the dielectronic recombination rates from Landini & Monsignori Fossi (1990). The analytic expression is

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{D}}(T_{\mathrm{e}})& =A_{\mathrm{di}}{T}_{\mathrm{e}}{[\mathrm{K}]}^{-3/2}\exp (-\frac{T_0}{T_{\mathrm{e}}[\mathrm{K}]})\hfill \\ \hfill & (1+B_{\mathrm{di}}\exp (-\frac{T_i}{T_{\mathrm{e}}[\mathrm{K}]}))[{\mathrm{cm}}^3{\mathrm{s}}^{-1}].\hfill \end{array}$$(A.16)

A.4. Five level system

In this paragraph we describe the numerical solution of a five level atomic/ionic system of species Y^k (see Fig. A.1). The energy levels, (internal) number densities, and statistical weights are labelled as E_i(Y^k), n_i(Y^k), and g_i(Y^k) with i = 1, …5, respectively. In a five level system the internal number densities n_i(Y^k) are connected to the total number density n(Y^k) via

$$\begin{array}{c}\hfill n(Y^k)=\sum _{i=1}^5{n}_i(Y^k).\end{array}$$(A.17)

We can assume that n(Y^k) is known from the ionisation balance and the internal level population are the quantities we want to determine. We consider three processes (rate coefficients in brackets): (a) spontaneous decay (A_ul), (b) collisional de-excitation (C_ul), and (c) collisional excitation (C_lu). Here, the indices u and l indicate the upper level and the lower level, respectively. In an equilibrium (${\dot{n}}_i(Y^k)=0$ for all i = 1…5) the rate equation for a fixed energy level i is given by (e.g. Osterbrock & Ferland 2006)

$$\begin{array}{cc}\hfill 0=& \sum _{j>i}(A_{\mathit{ji}}+C_{\mathit{ji}})n_j+\sum _{j<i}{C}_{\mathit{ji}}{n}_j\hfill \\ \hfill & -\sum _{j<i}(A_{\mathit{ij}}+C_{\mathit{ij}})n_i-\sum _{j>i}{C}_{\mathit{ij}}{n}_i.\hfill \end{array}$$(A.18)

We note that the internal level populations reach the equilibrium on very small timescales so that the equilibrium assumption here is justified. In a five level system only four rate equations are linearly independent. Therefore, we can replace the last rate equation (i = 5) with Eq. A.17 to formulate the matrix equation

$$\begin{array}{c}\hfill \left(\begin{array}{ccccc}{M}_{11}\hfill & M_{12}\hfill & M_{13}\hfill & M_{14}\hfill & M_{15}\hfill \\ M_{21}\hfill & M_{22}\hfill & M_{23}\hfill & M_{24}\hfill & M_{25}\hfill \\ M_{31}\hfill & M_{32}\hfill & M_{33}\hfill & M_{34}\hfill & M_{35}\hfill \\ M_{41}\hfill & M_{42}\hfill & M_{43}\hfill & M_{44}\hfill & M_{45}\hfill \\ 1\hfill & 1\hfill & 1\hfill & 1\hfill & 1\hfill \end{array}\right)·\left(\begin{array}{c}{n}_1(Y^{\mathrm{k}})\hfill \\ n_2(Y^{\mathrm{k}})\hfill \\ n_3(Y^{\mathrm{k}})\hfill \\ n_4(Y^{\mathrm{k}})\hfill \\ n_5(Y^{\mathrm{k}})\hfill \end{array}\right)=\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ 0\hfill \\ n(Y^{\mathrm{k}})\hfill \end{array}\right).\end{array}$$(A.19)

The coefficients of the matrix $\widehat{M}$ are functions of A_ul, C_lu, C_ul. Explicitly, they are

$$\begin{array}{cc}\hfill M_{11}& =-[C_{15}+C_{14}+C_{13}+C_{12}],\hfill \end{array}$$(A.20)

$$\begin{array}{cc}\hfill M_{12}& =C_{21}+A_{21},\hfill \end{array}$$(A.21)

$$\begin{array}{cc}\hfill M_{13}& =C_{31}+A_{31},\hfill \end{array}$$(A.22)

$$\begin{array}{cc}\hfill M_{14}& =C_{41}+A_{41},\hfill \end{array}$$(A.23)

$$\begin{array}{cc}\hfill M_{15}& =C_{51}+A_{51},\hfill \end{array}$$(A.24)

$$\begin{array}{cc}\hfill M_{21}& =C_{12}\hfill \end{array}$$(A.25)

$$\begin{array}{cc}\hfill M_{22}& =-[C_{25}+C_{24}+C_{23}+C_{21}+A_{21}],\hfill \end{array}$$(A.26)

$$\begin{array}{cc}\hfill M_{23}& =C_{32}+A_{32},\hfill \end{array}$$(A.27)

$$\begin{array}{cc}\hfill M_{24}& =C_{42}+A_{42},\hfill \end{array}$$(A.28)

$$\begin{array}{cc}\hfill M_{25}& =C_{52}+A_{52},\hfill \end{array}$$(A.29)

$$\begin{array}{cc}\hfill M_{31}& =C_{13},\hfill \end{array}$$(A.30)

$$\begin{array}{cc}\hfill M_{32}& =C_{23},\hfill \end{array}$$(A.31)

$$\begin{array}{cc}\hfill M_{33}& =-[A_{31}+A_{32}+C_{31}+C_{32}+C_{35}+C_{34}],\hfill \end{array}$$(A.32)

$$\begin{array}{cc}\hfill M_{34}& =A_{43}+C_{43},\hfill \end{array}$$(A.33)

$$\begin{array}{cc}\hfill M_{35}& =A_{53}+C_{53},\hfill \end{array}$$(A.34)

$$\begin{array}{cc}\hfill M_{41}& =C_{14},\hfill \end{array}$$(A.35)

$$\begin{array}{cc}\hfill M_{42}& =C_{24},\hfill \end{array}$$(A.36)

$$\begin{array}{cc}\hfill M_{43}& =C_{34},\hfill \end{array}$$(A.37)

$$\begin{array}{cc}\hfill M_{44}& =-[A_{41}+A_{42}+A_{43}+C_{41}+C_{42}+C_{43}+C_{45}],\hfill \end{array}$$(A.38)

$$\begin{array}{cc}\hfill M_{45}& =A_{54}+C_{54}.\hfill \end{array}$$(A.39)

Equation A.19 can be solved numerically for the internal number densities n_i(Y^k). In the case of dominant electron collisions, the collision rates C_ul(n_e, T_e, Y^k) and C_lu(n_e, T_e, Y^k) are given by

$$\begin{array}{cc}\hfill C_{\mathit{ul}}& =\frac{8.629\times {10}^{-6}{n}_{\mathrm{e}}}{\sqrt{T_{\mathrm{e}}[\mathrm{K}]}}·\frac{{\mathrm{{\rm Y}}}_{\mathit{lu}}(T_{\mathrm{e}},Y^k)}{g_u},\hfill \end{array}$$(A.40)

$$\begin{array}{cc}\hfill C_{\mathit{lu}}& =\frac{8.629\times {10}^{-6}{n}_{\mathrm{e}}}{\sqrt{T_{\mathrm{e}}[\mathrm{K}]}}·\frac{{\mathrm{{\rm Y}}}_{\mathit{lu}}(T_{\mathrm{e}},Y^k)}{g_l}·\exp (-\frac{E_u-E_l}{k_{\mathrm{B}}{T}_{\mathrm{e}}}).\hfill \end{array}$$(A.41)