© The Authors 2023

1 Introduction

Photoionisation is a ubiquitous process in astrophysics. Any gaseous astrophysical environment embedded in a field of ionising photons (E > 13.6 eV for ground state hydrogen) gets photoionised to a level that depends critically on the ratio of the density of ionising photons at the illuminated face of the cloud to the electron density in the gas (i.e. the ionisation parameter U = (Q_ion/4πr²c)/n_H, where Q_ion is the rate of ionising photons emitted by the source and n_H, r, and c are the hydrogen electron density, the distance of the illuminated face of the gas cloud from the ionising source and the speed of light, respectively¹). A few important examples are the gaseous environment of active galactic nuclei (AGNs) and accreting compact sources, the interstellar medium (ISM) surrounding explosive events like gamma-ray bursts (GRBs), the intergalactic medium embedded in the diffuse radiation field, and the HII regions around O and B stars.

Diagnostics on the ionising source and on the physics, kinematics, and chemistry of the ionised gas can be gathered through UV and X-ray spectroscopy, combined with gas ionisation modelling. Several public photoionisation codes do exist; they describe (given a set of input parameters) the stationary state of a gas cloud at equilibrium (e.g. Cloudy, Ferland et al. 2017; XSTAR, Kallman & Bautista 2001; Kallman et al. 2021; SPEX, Kaastra et al. 1996). Although each has its own working structure and atomic libraries, they generally model equilibrium configurations (i.e. they assume gas in local photoionisation equilibrium). In such circumstances, for a given incident spectral energy distribution and metallicity, the ionic population of the gas, and thus the emerging absorption spectrum, is uniquely determined² by the value of U. From an observational measurement of U, Q_ion it is possible to derive an estimate of the product r² · n_H, but not independently of n_H and r.

However, photoionisation equilibrium does not necessarily apply to gaseous environments illuminated by highly variable sources, such as those surrounding GRBs and AGNs, which are known to exhibit luminosity variations of up to several orders of magnitude. When the ionising source varies on timescales shorter than the equilibration timescale of the gas t_eq, equilibrium models can only provide an average description of the gas properties and may yield the wrong diagnostic on the physics and ionisation state of the gas when compared to observations. The equilibration time t_eq critically depends on the free electron density n_e and (in a three-ion approximation) can be expressed as (Nicastro et al. 1999) $$t_{\text{eq}}\approx {\left\{\left({\alpha }_{\text{rec}}\left(n_{{\text{X}}^i},T\right)\cdot n_{\text{e}}\right)\cdot \left(1+\frac{{\alpha }_{\text{rec}}\left(n_{{\text{X}}^{i-1}},T\right)}{{\alpha }_{\text{rec}}\left(n_{{\text{X}}^i},T\right)}+\frac{n_{{\text{X}}^{i-1}}}{n_{{\text{X}}^i}}\right)\right\}}^{-1}$$(1)

where α_rec and $n_{{\text{X}}^i}$ are the recombination coefficient and the fractional abundance of the i-th ion of the element X, respectively, and T is the electron temperature. We note that t_eq accounts for the time required to reach a new equilibrium and, depending on whether the luminosity increases or decreases, is similar to the photoionisation and recombination timescales, respectively (see García et al. 2013; Sadaula et al. 2023). In these cases, only time-evolving photoionisation models (combined with time-resolved spectroscopy) can provide instantaneous diagnostics on the physics and kinematics of the gas, and also offer an efficient way to remove the gas density-distance degeneracy intrinsic in the definition of U, both through changes in relative ion abundance (e.g. Nicastro et al. 1999; Krongold et al. 2007) and electron population of metastable levels (e.g. Kaastra et al. 2004) in response to flux variations. This is the focus of the present paper.

GRB studies of the temporal variation of absorption lines in the optical/UV band have provided a wealth of information on the medium surrounding the burst and the ISM of the host galaxy along our line of sight, including its density, metallic-ity and distance from the burst (see e.g. D’Elia et al. 2009; Heintz et al. 2018). These results were based on ad hoc time-evolving level-population models built under the UV photo-pumping hypothesis for a number of selected FeII transitions rather than a general-purpose multi-wavelength time-evolving photoionisation code. In the X-rays, the limited resolution and signal-to-noise ratio of current observations, together with the lack of self-consistent time-evolving photoionisation models, have so far hampered detailed spectroscopic studies and investigations of the ionisation structure of the intervening gas. X-ray absorption in the ~0.5–10 keV portion of X-ray afterglow spectra was instead modelled phenomenologically, adopting single-zone cold (i.e. neutral) photoelectric absorption models, which typically provided reasonable fits of the data (see e.g. Campana et al. 2012; Asquini et al. 2019). Alternatively, single-zone, equilibrium ionised absorber modelling was attempted on the X-ray spectra of X-ray bright and optically dark GRB afterglows (e.g. Piro et al. 2002), which exhibit substantial absorption in X-rays and the optical. These studies, however, only provided upper limits on the ionisation parameter, consistent with absorption by either dense absorbing clouds of the star-forming region embedding the GRBs or the more distant ISM of the host galaxy.

Time variability represents one of the best tools for exploring the gaseous environment of AGNs as well. The response of nuclear absorbers to central source luminosity variations has provided fundamental insights into the gas distance, density, geometry and launching mechanisms in a number of sources, revealing a close interaction between the accretion-powered luminosity and the outflowing gas at accretion disc scales (e.g. Nicastro et al. 1999; Peterson et al. 2004; Krongold et al. 2007; Bianchi et al. 2009; Parker et al. 2017; Kara et al. 2021). These results, however, were reached either through simple recombination and variability timescale comparisons or through simple time-evolving photoionisation models (e.g. Nicastro et al. 1999).

Flexible, fast and self-consistent time-evolving photoionisation codes, able to produce time-resolved spectra for any kind of ionising spectral energy distribution (SED) and light curve over the whole rest-frame UV and X-ray bands, have only very recently started to become available (Sadaula et al. 2023; Rogantini et al. 2022) and are urgently needed to fill the gap between current models and the new generation of upcoming spectrometers. In this paper we present the novel code Time Evolving PhotoIonisation Device (TEPID). TEPID follows the energy and ion abundance evolution of a gas undergoing time variable photoionisation and, in its current version, computes the corresponding absorption spectra in the UV and X-ray bands. The main outputs of a TEPID run are the gas temperature, cooling and heating rates, ion abundances and column densities. These quantities are fed into a customised version of the PHotoionised Absorber Spectral Engine (PHASE; Krongold et al. 2003), which builds tables of time-resolved UV and X-ray absorption spectra for a number of input parameters (e.g. time, electron density, ion column density, etc.) that can be used in fitting packages (such as XSPEC, Arnaud 1996) to fit the data and constrain the properties of the absorbing gas. TEPID is highly flexible and can be applied for diagnostics to any of the astrophysical scenarios described above. It has been tested over a broad range of electron densities, equivalent hydrogen column densities, and ionising fluxes (i.e. ionisation parameters) and compared, at equilibrium, with existing photoionisation codes. There are two main limitations of the current version of TEPID. First, the time-evolving level population balance is not computed and all ions are considered in their ground states; second, and as a consequence, gas emissivity is not produced and only a first-order radial-opacity radiative-transfer is performed throughout the cloud (see Sect. 5 for a discussion on these limitations).

In this first paper, we present the code and its application to two different astrophysical cases: the highly ionised gaseous environment of AGNs and the post-explosion surroundings of GRBs. The two environments undergo dramatically different time-evolving photoionisations. In the first case, luminosity changes are stochastically distributed, both in time and magnitude, with intensity variations of up to few orders of magnitude on timescales from a kilosecond to hundreds of kiloseconds and longer. In the second case, instead, the ionising light curve is typically described by an initial short high-luminosity phase called the prompt (which for the purpose of this paper can be approximated as a constant), shortly followed by a steep and virtually monotonic decaying phase lasting several hundreds of kiloseconds, during which the luminosity drops by several orders of magnitude (the afterglow).

The structure of the paper is as follows. In Sect. 2, we describe the general structure of the code and its main outputs. Then in Sects. 3 and 4 we present the results for the AGN and GRB cases, respectively. Finally, the discussion and our conclusions are in Sects. 5 and 6.

2 The TEPID code

2.1 Code overview and input parameters

We developed TEPID starting from the prototypical photoionisation code of Nicastro et al. (1999, hereafter n99), further adapted by Krongold & Prochaska (2013, hereafter K13) to the case of GRB X-ray afterglows. The code solves a set of first-order differential equations following the gas energy (cooling-heating) and ion-fraction (ionisation-recombination) balances as a function of time for the ten astrophysically most abundant elements: H, He, C, N, O, Ne, Mg, Si, S, Fe. The code does this self-consistently and with an in-flight optimised time-resolution, returning at each instant and at each position in the cloud (see below) the gas ionisation balance, temperature T and transmitted absorption spectrum.

We follow the ionisation of the cloud in space, from its illuminated face outwards, by implementing a first-order radiative transfer calculation. We divide the gas cloud into a series of optically thin adjacent shells. The time-evolving calculation is performed shell by shell, starting from the illuminated face of the cloud, from t = 0 to t = t_end and saving all relevant quantities at each time resolution element. The calculation is then repeated at all times for the adjacent shell by using as input, at each time, both the quantities saved at the previous spatial iteration and the geometrically diluted and opacity-attenuated ionising flux emerging from the previous shell and impinging the shell under consideration (see Eq. (9)). The calculation is then iterated up to the last shell. The gas temperature and ionisation balance are saved as output of the run for each time step and for each shell. Finally, the ionic population of the gas cloud is fed to a customised version of the PHASE spectral interface (Krongold et al. 2003), which produces the time-resolved absorption spectra, including both continuum and line absorption, to be compared with observations within the xspec spectral fitting package.

For each simulation run the input parameters of the code are the following:

• $n_{\text{e}}^{\text{tot}}$, the total (bound + free) electron density;

• the spectral energy distribution (SED). TEPID works in the spectral range 10⁻⁵−100 keV, from the infrared, which is important for Compton cooling (see Eq. (5)), to the hard X-rays to properly account for ionisation and heating;

• the initial ion abundances $n_{{\text{X}}^i}$ at t = 0. They can be specified either by providing a value of the ionisation parameter U or through a set of equilibrium ion abundances;

• F_ion(t), the radiative flux impinging on the gas cloud as a function of time (i.e. the light curve).

The code assumes spherical symmetry and can work in two geometrical settings. In the plane-parallel scenario the geometrical thickness of the shell is negligible with respect to the distance from the luminosity source (i.e. Δr/r ≪ 1), and thus the initial T, $n_{{\text{X}}^i}$ are constant throughout the shell. When this condition is not fulfilled (i.e. Δr/r ≳ 1), the input T, $n_{{\text{X}}^i}$ are referred to the illuminated face of the cloud and are then propagated radially outwards throughout the gas column. Elemental abundances are also free parameters: we assume solar metallicity (Lodders 2003) throughout the paper.

We note that at each spatial location (i.e. at each ith shell) all our times t are actually relative times $\text{d}t=t-t_{\text{prop}}^i$, where $t_{\text{prop}}^i$ is the light travel time from the illuminated face of the cloud to the ith shell. These are actually also the ‘observer times’; in other words, at the relative time dt the observer sees photons that have been emitted by the source at time $t-t_{\text{prop}}^{\text{tot}}$ (where $t_{\text{prop}}^{\text{tot}}$ is the light travel time of the entire gas cloud) and that have already crossed the entire gas column along our line of sight on their way to us.

From the point of view of the transmitted spectrum along the line of sight, all that matters is that at each i-th shell the calculation is performed at relative times $\text{d}t=t-t_{\text{prop}}^i$. The time coordinate in each shell is set by imposing dt = 0 at t = 0 (i.e. correcting for $t_{\text{prop}}^i$). For dt < 0 (i.e. t < $t_{\text{prop}}^i$), the ith shell has not yet been reached by the first photons of the time-evolving computation and its gas lies in the status described by the initial conditions (i.e. initial temperature and ion abundances; see Sect. 2.2). We refer to Krolik & Kriss (1995) for further discussion on this point.

The new features of TEPID with respect to N99 and K13 are i) the self-consistent computation of the cooling-heating balance; ii) the in-flight optimised time resolution; iii) the time resolved absorption spectra; and iv) the inclusion of collisional ionisation, dielectronic recombination, and updated data for photoionisation and radiative recombination (see Sect. 2.3).

2.2 Initial conditions

For a given set of initial conditions the code starts the calculation by guessing a temperature of the gas and computing corresponding equilibrium abundances and cooling and heating rates. These are then used to iterate the calculation by varying T until heating and cooling balance each other, and thus thermal equilibrium is reached. As an example, Fig. 1 shows the photoionisation equilibrium temperature T as a function of U (top panel) and U/T (i.e. radiation over gas pressure ratio, bottom panel) for an optically thin cloud (log(N_H/cm⁻²) = 18) of gas illuminated by the standard AGN SED tabulated in Cloudy³. The values computed with TEPID (red line) agree reasonably well with those obtained with Cloudy (green) and XSTAR (blue) over more than three orders of magnitude in U (fully encompassing the range of values typically observed in AGN ionised absorbers; see e.g. Laha et al. 2021), and almost two orders of magnitude in U/T (see Sect. 5 for a discussion on the observed differences)⁴.

2.3 Equations

The temporal evolution of the relative abundance of the ith ion of a given element X, $n_{{\text{X}}^i}$, can be written as $$\begin{array}{l}\frac{\text{d}{n}_{{\text{X}}^i}}{\text{d}t}=-\left[F_{{\text{X}}^i}+C_{{\text{X}}^i}{n}_{\text{e}}+{\alpha }_{\text{rec}}\left({\text{X}}^{i-1},T\right)n_{\text{e}}\right]n_{{\text{X}}^i}\hfill \\ +\left(F_{{\text{X}}^{i-1}}+C_{{\text{X}}^{i-1}}{n}_{\text{e}}\right)n_{{\text{X}}^{i-1}}+{\alpha }_{\text{rec}}\left({\text{X}}^i,T\right)n_{\text{e}}{n}_{{\text{X}}^{i+1}}+{\displaystyle \sum _{k=1}^{i-2}{F}_{{\text{X}}^k}^A{n}_{{\text{X}}^k}.}\hfill \end{array}$$(2)

The first term on the right-hand side is the destruction term due to ionisation from the stage i to i + 1, both radiative and collisional, with temporal rates $F_{{\text{X}}^i}$ and $C_{{\text{X}}^i}$, respectively, and recombination to the stage i − 1 with rate α_rec. The second and third source terms are due to ionisation from the stage i − 1 to i and recombination from the i + 1 stage. Finally, the last term accounts for Auger ionisation, which originates from the ionisation of an inner shell electron, and leads to multiple ionisation events. Photoionisation rates are computed by integrating fluxes over the photoelectric cross-sections tabulated in Cloudy (Ferland et al. 2017), while collisional rates are from Dere (2007). The rate α_rec is the sum of radiative and dielectronic recombination rates, which are both functions of the gas temperature T and the free electron density n_e. For hydrogen we use the values from Ferland (2002) and we sum over levels from 1 to 20 to obtain the total recombination rate. For heavier elements the total rates to the ground states are taken from a) Badnell (2006a,b) for ions with 0–14 electrons (i.e. fully stripped to Si-like); b) Shull & van Steenberg (1982) and Mazzotta et al. (1998) for ions with more than 14 electrons (i.e. low charge S and Fe) for radiative and dielectronic recombination, respectively. Auger probabilities are from Kaastra & Mewe (1993).

The gas temperature T is computed by integrating heating and cooling rates of each ion, weighted for the ion abundance $n_{{\text{X}}^i}$ and the element abundance Z_X, plus a Compton energy exchange term Θ: $$\frac{\text{d}T}{\text{d}t}={\displaystyle \sum _{\text{X,}i}\left(H\left({\text{X}}^i\right)-\text{Λ}\left({\text{X}}^i\right)\cdot n_{\text{e}}\right)\cdot {\text{Z}}_{\text{X}}{n}_{{\text{X}}^i}+\text{Θ}\left(F_{\text{ion}},T\right)}$$(3)

The heating rate H(Xⁱ) consists of two terms: $$H\left({\text{X}}^i\right)={\displaystyle {\int }_{v_{{\text{X}}^i}}^{\infty }\frac{F\left(v\right){\sigma }^{{\text{X}}^i}\left(v\right)}{hv}}\left(hv-h{v}_{{\text{X}}^i}\right)\text{d}\left(hv\right)+E_{\text{A}}\left({\text{X}}^i\right)\cdot {\displaystyle \sum _{k=1}^{i-2}{F}_{{\text{X}}^k}^{\text{A}}}\cdot$$(4)

The first term is the heating from photoionised electrons (see e.g. Netzer 2013). The quantity F(v)/hv is the incoming photon flux and it is linked to the ionising flux as $F_{\text{ion}}={\displaystyle \int {}_{{10}^{-5}}^{100\text{\hspace{0.17em}keV}}F\left(v\right)\text{d}v}$, while ${\sigma }^{{\text{X}}^i}\left(v\right)$ is the ion photoelectric cross-section as a function of energy and ($hv-h{v}_{{\text{X}}^i}$) is the energy that goes into heating, which is the energy of the incident (absorbed) ionising photon minus the binding energy of the electron. The integral is performed from the photon energy threshold $h{v}_{{\text{X}}^i}$ up to the high-energy boundary of the SED (i.e. 100 keV). The second term is due to the additional energy E_A carried by the Auger electrons.

The cooling rates for each X_i are taken from Gnat & Ferland (2012) and depend on T and (linearly) on n_e. They include line emission, recombination, collisional ionisation, and thermal bremsstrahlung.

The net Compton rate Θ accounts for the energy exchange between the gas and the radiation field via Compton scattering. It can be written as (Levich & Sunyaev 1970; Ferland 2002) $$\text{Θ\hspace{0.17em}=\hspace{0.17em}}\frac{2\pi }{m_{\text{e}}{c}^2}\cdot \frac{n_{\text{e}}}{n_{\text{e}}^{\text{tot}}}\left({\displaystyle \int {\sigma }_{\text{h}}F\left(v\right)hv\text{\hspace{0.17em}d}\stackrel{⌣}{v}4\text{\hspace{0.17em}KT}{\displaystyle \int {\sigma }_{\text{c}}F\left(v\right)dv}}\right)$$(5)

and includes both heating and cooling terms. In Eq. (5) m_e, σ_h and σ_c are the electron mass and the Compton heating and cooling cross-sections, respectively.

The free electron number density n_e is given by the weighted sum over all the elements X and all the ionic levels i, each multiplied by the number of ionised electrons $$n_{\text{e}}=n_{\text{e}}^{\text{tot}}\cdot \frac{{\displaystyle {\sum }_{\text{X}}\left({\displaystyle {\sum }_{i=2}^{K_{\text{X}}+1}{n}_{{\text{X}}^i}\cdot \left(i-1\right)}\right){\text{Z}}_{\text{X}}}}{{\displaystyle {\sum }_{\text{X}}\left(K_{\text{X}}\cdot {\text{Z}}_{\text{X}}\right)}}$$(6)

where K_X is the atomic number and Z_X the element abundance. Finally, the frequency-dependent optical depth within the mth shell, τ_m(v) is given by $${\tau }_m\left(v\right)=\delta N_{\text{H}}\cdot {\displaystyle \sum _{\text{X}}{\text{Z}}_{\text{X}}\left({\displaystyle \sum _{i=1}^{K_{\text{X}}}{n}_{{\text{X}}^i}{\sigma }^{{\text{X}}^i}}\left(v\right)\right)}$$(7)

where δN_H is the hydrogen-equivalent column density of the shell, which is used as a proxy for the total gas column and is obtained by integrating n_H along the line of sight within the shell. The relation between hydrogen and electron number densities depends on the gas metallicity; for solar values (as we use in this paper) $n_{\text{H}}=n_{\text{e}}^{\text{tot}}/1.2$.

Ionisation and heating rates in the mth shell are computed using the spectrum $$F{\prime }_m=F_m\cdot \frac{1-{\text{e}}^{-{\tau }_m}}{{\tau }_m}$$(8)

where F_m is the incident spectrum at the innermost face of the shell. The spectrum $F{\prime }_m$ accounts for the attenuation of the continuum across the mth shell (see e.g. Ferland et al. 2017). As discussed in Appendix A, using $F{\prime }_m$ instead of F_m allows the propagation of the radiation through the gas column to be followed to much higher accuracy.

Finally, the incident spectrum on the (m + 1)th slab is equal to the transmitted spectrum from the m-th slab ($F_m\cdot {\text{e}}^{-{\tau }_m}$) corrected for geometric dilution (i.e. the squared ratio between the radial locations r_m, r_m+1 of the two slabs) $$F_{m+1}=F_m\cdot {\text{e}}^{-{\tau }_m}\cdot {\left(\frac{r_m}{r_{m+1}}\right)}^2.$$(9)

Equations (2)–(9) all depend explicitly on time, and form a coupled system. Therefore, in our code we integrate these equations simultaneously as a function of time. The following subsections show how these equations are solved to obtain the temporal evolution of the gas.

Fig. 1 Photoionisation equilibrium temperature in an optically thin cloud of gas ionised by the standard AGN SED of Cloudy as a function of the ionisation parameter U (top) and the U/T ratio (bottom). The green, blue and red curves correspond to the values computed with Cloudy, XSTAR and TEPID, respectively.

Fig. 2 Adaptive temporal binning by TEPID. Top panel: step-function light curve. Middle: Sampling frequency ω. Bottom: zoomed-in images of ω around the flux transition phases (see text for details).

2.4 Adaptive time binning

In order to optimise memory allocation and running time, the temporal sampling frequency ω is self-evaluated by the code through an adaptive approach. In the following example we use an ideal two-phase step-function light curve (Fig. 2, top panel), where at t = 0 the gas is in equilibrium with the ionising flux. The flux then increases instantaneously (dF_ion/dt = ∞) by a factor of 10 and stays at this level for 10 ks (shaded light blue region in top panel of Fig. 2). At t = 10 ks the flux abruptly goes back to its initial value and stays at this level for an additional 10 ks (shaded orange region) up to t_end = 20 ks. Given the dependence of the equilibration time t_eq on n_e (Eq. (1)), the gas reacts to the luminosity changes at t = 0 and 10 ks on a timescale inversely proportional to n_e. At high densities the gas adjusts quickly to the abrupt changes in flux and so the time resolution of the integration has to be high, while at low density the gas takes longer to adjust to the new flux values and lower time resolutions are sufficient, even at the flux transition phases. This is done automatically through an algorithm, as shown in the middle and bottom (zooms around the flux transition phases) panels of Fig. 2, that show the frequency oo of the time-integration sampling for two dramatically different volume electron densities, log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6 (red curves) and 10 (blue curves).

The algorithm works as follows:

• The code scans the entire ionising light curve, from t_start to t_end, and divides the time domain into a number of N_gross intervals defined by looking for positive or negative flux changes |ΔF_ion ≥ |10%: every time this condition is met, a new time-interval is defined.

• For each time interval $n_{\text{gross}}^i$ (where $1\le n_{\text{gross}}^i\le N_{\text{gross}}$), with width $\text{Δ}{t}_{\text{gross}}^i$, the equilibration time of each ion is evaluated, using the initial set of equilibrium ion abundances. The longest of these timescales $t_{\text{eq}}^{\text{max,}i}\left(n_{\text{gross}}^i\right)$ (typically << $\text{Δ}{t}_{\text{gross}}^i$) is used within each $n_{\text{gross}}^i$ interval to define finer sub-grids, each with N_fine,i intervals with sizes $\text{Δ}{t}_{\text{fine,}i}^j=Int\left(\text{Δ}{t}_{\text{gross}}^i/t_{\text{eq}}^{\text{max,}i}\right)\times 2^j$, where j = 0, N_fine,i.

• The total number of finer intervals from t_start to t_end is given by $N_{\text{grid}}={\displaystyle {\sum }_{i=1}^{N_{\text{gross}}}{N}_{\text{fine,}i}}$. This series of intervals defines the time-dependent temporal resolution $\text{Δ}{t}_{\text{fine,}i}^j$. of the simulation, and thus the grid over which the relevant outputs of the computation (temperature, cooling and heating rates, fractional ion abundances, optical depth and transmitted spectrum) are stored in memory for the next spatial iteration. The value of $\text{Δ}{t}_{\text{fine,}i}^j$ is computed according to the formula in Gallavotti (1986) to ensure the numerical error of each time step is kept below a certain threshold.

• The inverse of the temporal resolution defines the time-dependent sampling frequency ${\omega }_j^i=1/\text{Δ}{t}_{\text{fine.}i}^j$, i.e. the quantity plotted in Fig. 2 (middle and bottom panels).

We verified that the above algorithm provides adequate resolution by re-running the simulations presented in this paper with ${\omega }_j^i$ increased by a factor of 10, and checking that all the results (i.e. temperature, ionic populations, transmitted spectra) are unaffected.

2.5 General working scheme

The code structure can be summarised as follows:

• After the time optimisation algorithm has defined the sampling frequency ω(t) (see Sect. 2.4), the time integration starts at the inner illuminated face of the gas cloud, r = r_in, and the temperature, ionisation and recombination rates within the first slab are updated at each time step and stored in memory.

• Next, the ion abundances and column densities are printed to an output file.

• Once t = t_end, the simulation moves to the second slab. The time integration is repeated by using as incident spectrum that emerging from the previous adjacent slab, corrected at each temporal resolution step $\text{Δ}{t}_{\text{fine,}i}^j$. for absorption and for geometric dilution according to Eq. (9). The output files for the ion abundances and column densities are updated with the same temporal resolution as above.

• The spatial iteration continues until the input hydrogen-equivalent column density N_H (or outward radius r_out) is reached.

The final outputs of the simulation are T, $n_{{\text{X}}^i}$, $N_{{\text{X}}^i}$, and N_H as a function of time t. We use these quantities to build PHASE table models (Krongold et al. 2003) containing time-resolved UV and X-ray absorption spectra emerging from the cloud. The spectra can then be used in xspec to perform time-resolved spectroscopic analysis of the observed data.

Fig. 3 AGN light curve.

2.6 Temperature balance and ion abundances

To illustrate how temperature and ionisation balance evolve under varying ionising flux, here we adopt a plane-parallel setting, where a thin layer of gas (log(N_H/cm⁻²) = 18) is in photoionisation equilibrium with the radiation field at t = 0. We assume log(U) = 0.5 and the standard Cloudy AGN SED described in Sect. 2.2. We employ the light curve shown in Fig. 3, a trapezoidal function with the flux F_ion increasing linearly from F_min = 1 (in arbitrary units) at t = 0 up to F_max = 10 at t₁ = 4 ks, staying at F_max from t₁ through t₂ = 6 ks and then decreasing again linearly (and with the same slope) to F_min = 1 from t₂ through t₃ = 10 ks, and remaining at F_min for the next 10 ks (up to t_end = 20 ks).

Figure 4 shows the temporal evolution of H, Λ and Θ for two values of log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 8 (top panel) and 12 (middle and bottom panels), while Figs. 5 and 6 show the temporal evolution of the temperature T and of the fractional abundances of the most abundant ions of carbon, oxygen, and iron for values of electron densities in the range log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 4–12. Generally, for all densities, lower ionisation stages (left and middle panels of Fig. 6) are depleted during the high-flux phase and are increased during the low-flux phase, while high-ionisation stages (right panels of Fig. 6) show the opposite trend.

In the highest density case, log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 12, the gas is in photoionisation equilibrium and readjusts practically instantaneously to the incident flux. All quantities follow the light curve closely. In particular, the time derivative of T (i.e. the algebraic sum H − Λ + Θ; see bottom panel of Fig. 4 for a zoomed-in image) is > 0 (net heating) during flux increases (from t₀ to t₁), then is = 0 during the maximum constant flux phase (from t₁ to t₂ where the temperature is constant and equal to 10⁶ K, the equilibrium temperature corresponding to F_max; see Figs. 1 and 5) and, finally, is negative (net cooling) during the flux decreasing phase (between t₂ and t₃). Accordingly, the temporal evolution of the ion abundances (Fig. 6) closely follows the light curve of the ionising source by quickly adjusting to the instantaneous equilibrium values.

For decreasing $n_{\text{e}}^{\text{tot}}$ the evolution of the gas physical properties departs from equilibrium. Ion-by-ion photoionisation rates do not explicitly depend on $n_{\text{e}}^{\text{tot}}$; however, for a given initial log(U) the incident flux F_ion, which is proportional to photoionisation rates, scales linearly with $r^2{n}_{\text{H}}\propto r^2{n}_{\text{e}}^{\text{tot}}$. Thus, at a given distance, lower values of $n_{\text{e}}^{\text{tot}}$ correspond to lower photoionisation rates. Instead, radiative and dielectric recombination rates explicitly depend on n_e and also on T via a polynomial fitting formula (see references in Sect. 2.2).

This leads to two time-dependent effects for decreasing $n_{\text{e}}^{\text{tot}}$. First, both recombination and photoionisation rates are slower, and thus the variations in T and $n_{{\text{X}}^i}$ in response to the flux changes are compressed and can also be delayed with respect to the changes. This, in turn, affects further recombination rates, which are slower not only because of the lower density, but also due to the compressed dynamics of temperature variations. The delayed response of the gas is clearly visible in Figs. 5 and 6, where the maximum T and fractional abundance of fully stripped ions are not simultaneous with the peak of the light curve (i.e. between t₁ and t₂), but are increasingly shifted towards the end of the high-flux phase (t = 10 ks) as density decreases, and the gas stays overheated and overionised at t > 10 ks. This is also evident in the top panel of Fig. 4 (log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 8), where dT/dt assumes negative values at t > 10 ks (i.e. net cooling, since the gas temperature is higher than the equilibrium one). As we will show below, these effects lead to unique spectral signatures in time-resolved spectra of variable ionising sources that, when properly modelled, allow the characterisation of the gas electron density $n_{\text{e}}^{\text{tot}}$.

Fig. 4 Heating H, cooling Λ, and Compton net heating Θ (red, blue, and green lines, respectively) in units of Ks−1 as a function of t. The black line represents the algebraic sum of these terms, corresponding to the time derivative of the temperature. The top and middle panels correspond to log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 8,12, while the bottom panel is a zoomed-in image of the temperature derivative for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 12.

Fig. 5 T as a function of time for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) from 4 to 12 (lightest to darkest line).

Fig. 6 Selected ion abundances for carbon, oxygen, and iron (from top to bottom) as a function of time for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) from 4 to 12 (thinnest to thickest line). The different panels correspond to different ions, as labelled. The y-scale is different for each panel.

Fig. 7 Selected ion abundances (from top to bottom and left to right, carbon, oxygen, iron; see legends) and temperature (bottom right) as a function of NH for a gas in photoionisation equilibrium, computed with TEPID (solid lines), Cloudy and XSTAR (dashed and dotted lines, respectively).

3 Time-resolved spectra of ionised AGN absorbers

In this section, we show an application of TEPID to the case of highly ionised absorbers in AGNs. We set the AGN SED of Sect. 2.6 for the ionising source and an initial equilibrium ionisation parameter of the gas to log(U) = 1.5 (a factor of 10 higher than in Sect. 2.6), typical of both the high-ionisation components typically observed in warm absorbers (see e.g. Krongold et al. 2003, 2009) and the low-ionisation counterparts (e.g. Krongold et al. 2021) of ultra-fast outflows (UFOs; e.g. Tombesi et al. 2010; Luminari et al. 2021). Simulations are run for a range of densities log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 4–12 and up to an equivalent hydrogen column density log(N_H/cm⁻²) = 24. As a consistency check of our radiative transfer (see Eqs. (8) and (9)), we first run TEPID at equilibrium (by inputting a constant light curve) and compare the resulting ion abundances and temperature with Cloudy. Then, we use the same light curve in Fig. 3 to illustrate the time-dependent behaviour.

3.1 Photoionisation equilibrium

Figure 7 shows the ion abundances for carbon, oxygen, iron, and T as a function of N_H for a gas in photoionisation equilibrium with log(U) = 1.5. The solid, dashed and dotted lines correspond to TEPID, Cloudy, and XSTAR, respectively. At low column densities the gas is quite transparent: its ionisation and temperature are that of an optically thin gas up to log(N_H/cm⁻²) ≈ 22.5. The agreement between the ionic populations of TEPID and Cloudy is within a factor of 2 for the main ions (i.e. those with abundance > 10⁻²). At higher N_H, the TEPID radiative transfer correctly accounts for the increasing optical thickness of the gas, and the (decreasing) gas ionisation and temperature are in good agreement with the Cloudy and XSTAR values.

3.2 Time-dependent optical depth

Figure 8, second and third rows, shows the ionic fractions of O and Fe as a function of N_H for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 10 and at two different times during the simulations: t₁ = 5 ks and t₂ = 16 ks (see caption). With an initial log U = 1.5, H, He, and light metals are fully ionised at virtually all times during the simulations, and the gas becomes optically thick only at log(N_H/cm⁻²) ≳ 23. Consequently, the gas column is virtually isothermal at all times during the simulation (i.e. T(t₁) = 3 × 10⁶ and T(t₂) = 10⁶) and has constant light metal abundances (i.e. f_CVII ≃ f_OIX ≃ 1) up to log(N_H/cm⁻²) ≃ 23 (middle panel of Fig. 8).

At log(N_H/cm⁻²) ≳ 22.5–23 (depending on the considered element), the gas becomes increasingly optically thick, and its temperature and ionisation both degree decrease steadily. However, the smoothness of this decrease depends on $n_{\text{e}}^{\text{tot}}$ and on the ionising flux state. At the chosen log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 10 and during the minimum flux state (t = 16 ks), the abundance of the fully stripped oxygen ion starts decreasing slowly at log(N_H/cm⁻²) > 23.8 (rightmiddle panel of Fig. 8), while the He-, Be-, and C-like fractional abundances of Fe start decreasing quickly already at log(N_H/cm⁻²) ≳ 22.5 to the advantage of the Ne- and O-like ions (right bottom panel of Fig. 8).

Fig. 8 Time-dependent optical depth. Top: AGN light curve. The thick lines show the times for which ion abundances are shown (see below). The thin lines show the times corresponding to the absorption spectra in Fig. 9. Bottom, first and second row: Selected ion abundances (see legend) respectively for oxygen and iron as a function of NH. The left and right panels correspond to t = 5, 16 ks. log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 10 in all cases.

3.3 Time-resolved spectra of AGN outflows

Figure 9 shows the absorption spectra extracted from the above simulation for log(N_H/cm⁻²) = 22 and log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6, 10 (blue and orange curves, respectively), at different times t = 0, 2, 8 and 16 ks. For reference, we report the Fe ion abundances as a function of time in the top panels. The left and right panels show, respectively, the spectral regions containing the Fe L (0.8–1.1 keV) and Fe K (6.2–7.2 keV) transitions. At t = 0 (top panels) the gas is in photoionisation equilibrium and therefore the spectra are independent of $n_{\text{e}}^{\text{tot}}$ (orange superimposed on blue).

At high density (log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 10, orange curves) the gas quickly responds to flux variations, though it is never in precise photoionisation equilibrium with the ionising flux (except at t = 0); therefore, the t = 2 ks and t = 8 ks spectra, corresponding to equal values of F_ion symmetrically preceding and following the maximum, respectively, are similar, but not identical. Analogously, the t = 16 ks spectrum (bottom panels) is similar (though not identical) to the initial equilibrium spectrum because both are taken at the minimum ionising flux level.

Conversely, when log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6 (blue curves) the gas is out of equilibrium at all times t > 0 and the temporal evolution of its ionisation balance and temperature is both compressed and delayed, as discussed in Sect. 2.6. This is clearly reflected in the spectra at t = 2, 8 and 16 ks. Little or no spectral variation is observed in the two spectral bands between t = 0 and t = 2 ks, despite the increase in the ionising flux. Between t = 2 and t = 10 ks the gas ionisation degree increases steadily, despite a first constant (t = 2–4 ks) and then steadily decreasing (t = 4–10 ks) ionising flux. Correspondingly the t = 8 spectrum is overionised with respect to its ionising flux and shows lower opacity in its low-ionisation Fe transitions than the spectrum at t = 2 ks. The gas ionisation keeps increasing even at t > 10 ks, while the ionising flux stays constant and equal to the minimum flux in the simulation. Consequently, the low-density spectrum at t = 16 ks is the least opaque of the four spectra of Fig. 9 in its Fe-L energy band. The gas eventually returns to equilibrium only after s¡300 ks has the flux returned to its minimum value.

Following the evolution in time of the ionisation degree of AGN outflows with X-ray telescopes offers a unique opportunity to constrain, via time-resolved spectroscopy, the electron density of the gas, and thus its distance from the ionising source. This, in turn, allows the mass load of the outflow and its energy to be estimated.

Fig. 9 Spectra of AGN outflows. Top: selected iron ion abundances (see legend) as a function of t for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6, 10 (left and right panels). The dashed lines show the times corresponding to the spectra shown below. Bottom: absorption spectra for log($N_{\text{H}}/{\text{cm}}^{-2}$) = 22 and t = 0, 2, 8, 16 ks (top to bottom, corresponding to the times indicated above and in Fig. 8) and log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6, 10 (blue and orange lines). The left and right columns show respectively the Fe L and Fe K bands (i.e. 0.8–1.1 and 6.2–7.2 keV).

Fig. 10 Pre-burst ionisation structure of the surroundings of a Wolf–Rayet star at the centre of a star-forming region with radius 100 pc and constant density $n_{\text{e}}^{\text{tot}}={10}^2{\text{cm}}^{-3}$. From left to right, the panels show C, O and Fe ion abundances and temperature as a function of the distance (in pc) from the star.

4 The ionisation of the surroundings of GRB afterglows

Here we use TEPID to reproduce the highly structured gaseous ionised environment of long GRBs, thought to originate from the death of massive stars, embedded in a dense star-forming region (Perna 2007; Heintz et al. 2018; Fryer et al. 2022). Our aim is to show how proper time-evolving photoionisation modelling of GRB X-ray afterglow spectra can efficiently constrain the physical and chemical properties of the interstellar medium of galaxies hosting a GRB explosion.

4.1 Physical setting of the GRB surroundings

For simplicity we assume that the ISM surrounding the GRB has constant number density $n_{\text{e}}^{\text{tot}}$ (i.e. independent of the radial distance r from the GRB). The inner boundary is set to r_in = 1 pc to encompass the forward shock, from where much of the X-ray radiation is irradiated. We set an outer radius r_out = 100 pc, which represents an upper envelope for the size of a typical OB stellar association (Stahler & Palla 2004), where the massive progenitor star is expected to assemble. Since the density is constant, the total gas hydrogen-equivalent column density is simply given by N_H = (r_out − r_in) · n_H. Unless otherwise specified, we focus on the density range log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1–4. We also neglect contributions from the outflows associated with the prompt phase or the progenitor star. This is justified since at a distance r_in the density of a typical Wolf–Rayet wind is already approximately one order of magnitude lower than the typical intragalactic densities (n ≃ 1 cm⁻³), and thus much lower than the range of densities we explore in our simulations (Chrimes et al. 2022; Eldridge et al. 2006).

Finally, for the pre-burst physical conditions of the gas permeating the star-forming region we assume those of a medium photoionised by a bright and hot Wolf–Rayet star, modelled as a black body with a temperature of 10⁵ K and a bolometric luminosity of 10^5.8 L_⊙ (Perna et al. 2018, but see Sect. 5 for relaxations of this hypothesis). Figure 10 shows the pre-burst most abundant ions of C, O, Fe and the temperature as a function of the radial distance.

4.2 The ionising source

We parametrise the source ionising flux F_ion(t), i.e. the GRB spectrum and its temporal evolution, as in K13, based on the analysis of the broad sample of Margutti et al. (2013) consisting of more than 650 GRBs observed with Swift. This has a different X-ray spectral shape and light curve in the prompt (t ≤ 100 s) and afterglow (t > 100 s) phases. The X-ray prompt phase is described as a power-law spectrum with energy spectral index α_xp = 0.0 (i.e. a flat spectrum), irradiating with constant luminosity up to t₁ =85 s and decaying as t⁻² between t₁ and t₂ = 100 s. At t₂ the GRB transits to the afterglow phase: the spectrum becomes steeper, α_xag = 1.0, and the light curve flattens as t⁻¹. Figure 11, top panel, shows the light curve from t = 0 to t_end = 10⁶, a typical end time of X-ray afterglow observations. The total irradiated isotropic energies in the range 0.3–10 keV during the prompt and afterglow (up to 10⁶ s) phases are 1.4 × 10⁵¹ erg and 5.4 × 10⁵¹ erg, respectively. Given the shape of the light curve and the extremely high initial flux, the internal time-step optimisation algorithm chooses a finer resolution at earlier times, which then gradually decreases with time. We plot the time sampling frequency ω in the bottom panel of Fig. 11. It can be seen that ω decreases from 10³ s⁻¹ at t = 0 to 10⁻³ s⁻¹ at t = t_end = 10⁶ s. The only exception is between 85 and 100 s, where oo increases to better follow the temporal decay transition.

Fig. 11 Luminosity and temporal resolution as a function of time in the GRB case. Top: GRB light curve. t1, t2, and tend correspond respectively to the end of the initial constant luminosity phase, the start of the afterglow and the end time of the simulation. Bottom: resolution ω as a function of time.

4.3 Temperature and ion abundances as a function of time and distance

Figure 12 shows the post-explosion temperature of the surrounding ISM as a function of the distance r from the GRB for a density log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1 and three increasing times (dotted, dashed and solid lines, respectively). When the burst explodes (t = 0) it quickly and fully ionises its immediate surroundings (photoionisation largely dominates over recombination, due to the low density and the 18 order of magnitude flux increase) and the temperature of the ISM reaches very high values (~2 × 10⁵ K at t = 100 s, in this example). The extent of the ionised region expands with time (dotted to solid lines) driven by the ionising flux, which however becomes geometrically diluted and increasingly absorbed for increasing distance. As an example, the flux of ionising photons at r = 20 pc is lower than that at r = 10 pc by a factor (20/10)² = 4, and therefore it takes four times longer to collect the same number of ionising photons (and thus to reach the same degree of ionisation within the medium), provided that F_ion stays constant (i.e. during the prompt phase). During the afterglow the expansion of the ionised region is even slower, since F_ion decrease with time as ∝ 1/t (see Fig. 11).

The temperature smoothly decreases with distance up to critical radii that depend on the time after the explosion and the density of the ISM (and so its cumulative opacity to the incoming radiation). At these radii the temperature profile displays sudden jumps or discontinuities and then continues decreasing smoothly, eventually approaching the temperature of the pre-burst ISM at large distances. These time-dependent critical distances are set by the recombination fronts of abundant elements in high-ionisation states (particularly H and He, but also, to a lesser extent, lower ionisation species of heavier elements). These recombination fronts can be clearly seen in Fig. 13, where the contributions of H and He ions to the heating and cooling rates (top left and right panels) are displayed together with their fractional abundance (bottom left and right panels, respectively) as a function of radial distance from the burst for t = 100 s and log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1. The first recombination front is found at a distance of about 11 pc and is due to recombination of HeIII in HeII (bottom right panel), while the second front is located at ~20 pc and is due to both the recombination of HeII in HeI and, most importantly for cooling, the recombination of H, which becomes virtually all neutral at r ≳ 20 pc (bottom left panel). Correspondingly, heating and cooling rates show peaks and edges at the recombination fronts. The temporal evolution of such fronts, which move outwards as the medium becomes progressively ionised over time, leads to the peaks in the heating and cooling, and thus in the temperature shown in Fig. 12 between ≈ 10 and 30 parsec.

Finally, Fig. 14 shows the abundances of the most abundant ions of oxygen (top panels) and iron (bottom panels) as a function of the distance r from the GRB (bottom X-axis) and the line of sight hydrogen-equivalent column density N_H (top X-axis) for three increasing times (dotted, dashed and solid curves). The panels from left to right correspond to different values of the ISM density, log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = −1,1,4. We include log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = −1 to explore the case of the recently observed hybrid bursts, displaying a long-GRB-like radiative emission but happening outside star-forming regions, and thus embedded in lower density gas (see e.g. Troja et al. 2022).

At all sampled times the ionisation degree of the ISM surrounding the GRB is highly structured and stratified, and it remains asymptotically stratified for years after the afterglow emission has completely faded (the time dictated by the slow gas cooling and recombination rates, which both depend on the gas density). Close to the GRB all the elements are highly ionised, up to hydrogen-like and fully stripped. The degree of ionisation and the temperature decrease with increasing r and approach the pre-burst values at a sufficient distance from the GRB.

However, for a given N_H the relative abundances of high-versus low-ionisation ions strongly increases with $n_{\text{e}}^{\text{tot}}$ since the denser the medium, the smaller (and closer to the GRB) the radial interval enclosing that N_H. This can be seen by expressing the geometric dilution of the flux as a function of N_H, i.e. $F_{\text{ion}}\left(N_{\text{H}}\right)=F_{\text{ion}}\left(r_{\text{in}}\right)\cdot {\left(r_{\text{in}}/r_{N_{\text{H}}}\right)}^2\propto F_{\text{ion}}\left(r_{\text{in}}\right)\cdot {\left(n_{\text{e}}^{\text{tot}}/N_{\text{H}}\right)}^2$; in other words, at a given N_H (enclosed in a radius $r_{N_{\text{H}}}$) the ionising flux is ∝ ${\left(n_{\text{e}}^{\text{tot}}\right)}^2$, and thus less diluted with increasing density. As a result, the transmitted spectra will be different, and thus high-resolution X-ray absorption spectra of GRB afterglows taken at different times after the explosions can yield a direct tomography of the GRB environment, as we show in Sect. 4.4 below.

Fig. 12 Temperature for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1 as a function of the radius r (i.e. the distance from the GRB) for t = 102,104,106 s (dotted, dashed and solid lines, see legend).

Fig. 13 Contribution of hydrogen and helium to the total heating and cooling rates. Top panels: heating (left) and cooling (right) rates as a function of the radius r. The red line in the left panel and the blue line in the right panel correspond to the total rates, while the green and yellow shaded areas correspond to the total hydrogen and helium rates, respectively. The green(yellow) solid and dashed lines correspond to the rates of H(He) I, II, respectively, while the yellow dot-dashed line to He III. Bottom: hydrogen and helium (left and right, respectively) ion abundances (see legend for colour-coding) as a function of r. In all cases log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1, t = 102 s.

4.4 Absorption spectra

Figure 15, left panel, shows the broad-band 0.1–10 keV absorption spectra imprinted on an incident GRB afterglow power-law continuum with photon index Γ = 2 emerging from the first 100 pc of circumburst medium, at t = 10⁴ s and for increasing log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$), from -1 to 4 (and thus increasing log(N_H/cm⁻²), from 19.5 to 24.5). The right panel shows a zoomed-in image of the 0.5–0.7 keV range, where oxygen lines from neutral (OI) to hydrogen-like (O VIII) transitions lie. Higher $n_{\text{e}}^{\text{tot}}$ probes larger and, on average, less ionised columns of circumburst medium (see Fig. 14). This implies emerging spectra attenuated by relatively high fractions of neutral (and so more opaque) gas for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 4 and by increasingly larger fractions of ionised gas as the density decreases.

Figure 16 shows the same spectra for the same range of log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) and a fixed column density of log(N_H/cm⁻²) = 22. This translates into a broad range of radii according to the electron density. For log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = 4(= 3), the radius is only 3.2(= 0.32) pc, and therefore the medium is totally ionised and the spectra are unabsorbed (except for some lines for log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = 3). On the other hand, for log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = −1 the radius is 32 kpc and the absorbing gas is mostly neutral.

The dependence on $n_{\text{e}}^{\text{tot}}$ is a peculiar feature of non-equilibrium photoionisation that cannot be reproduced with equilibrium, steady-state photoionisation models. As shown above, the simultaneous presence of all these ionisation stages is due to the stratification of the ionisation structure of the intervening medium, whose closer-in ionised regions have no time to recombine before the afterglow fades, while farther-out zones receive a progressively geometrically diluted number of ionised photons and never reach high degrees of ionisation. The out-of-equilibrium ratio between different ionisation stages of a given element in the observed spectra is thus unique to the specific physical properties of the medium surrounding the GRB (namely its density profile and extent of the star-forming region in the line of sight direction) and is very different from what would be expected if the gas were in photoionisation equilibrium (see also Lazzati & Perna 2003). This can be clearly seen in Fig. 17, where we compare the oxygen ion abundances predicted by TEPID for log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = 1 and t = 10⁴ s with those expected if the medium were in photoionisation equilibrium with four different values of the ionisation parameter. The red line shows the OI-IX abundances averaged over a gas column log(N_H/cm⁻²) = 21 (i.e. 32 pc), while the blue to green lines correspond to the equilibrium abundances (computed with Cloudy) for log(U) from −3 to 0 and averaged over log(N_H/cm⁻²) = 21. When the gas is in photoionisation equilibrium with the ionising source, ion abundances peak around a single ionisation stage (see also Kallman et al. 2004). Our time-evolving abundances, instead, clearly show the stratification of the medium. At equilibrium the abundance of each ionic stage is simply proportional to the ratio of its recombination to photoionisation rates. Here, instead, recombination is negligible with respect to ionisation (see Fig. C.1) because of the high-ionising flux and of the relatively low values of n_e, and so the (rising) gas ionisation for increasing time is simply given by the continuous photoionisation of the gas due to the GRB radiative output.

The strong dependence of the afterglow transmitted spectrum on $n_{\text{e}}^{\text{tot}}$ offers a valuable opportunity to directly constrain it in observed spectra by fitting afterglow spectra with time-evolving photoionisation models such as TEPID.

Fig. 14 Abundances for some oxygen (top) and iron (bottom) ions (see legend for colour-coding) as a function of the radius r. The dotted, dashed and solid lines identify increasing times from the onset of the GRB, t = 102,104,106 s, respectively, while the left, middle and right panels correspond to increasing ISM density, log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = −1, 1, 4.

Fig. 15 Absorption spectra for a distance r = 100 pc from the GRB and increasing $n_{\text{e}}^{\text{tot}}$ (and thus increasing NH; see legend for colour-coding). Left: broad-band 0.1–10 keV spectra. Right: zoomed-in image of the 0.5–0.7 keV range, showing several oxygen absorption lines from neutral to hydrogen-like states.

Fig. 16 Absorption spectra for column density NH = 1022 cm−2 and increasing $n_{\text{e}}^{\text{tot}}$ (and thus decreasing distance from the GRB; see legend for colour-coding). Left: broad-band 0.1–10 keV spectra. Right: zoomed-in image of the 0.5–0.7 keV range.

Fig. 17 Relative abundances of the oxygen ion stages. Red line: TEPID abundances for log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = 1, t = 104 s averaged up to log(NH/cm−2) = 21. Blue to green lines: Photoionisation equilibrium abundances for log(U) from −3 to 0, averaged up to log(NH/cm−2) = 21.

4.5 Impact of the pre-burst conditions

To quantify the impact of the initial ionisation of the surrounding medium we compare the results obtained under our assumption of an initial Wolf–Rayet-like pre-ionisation with those obtained by assuming an initially neutral medium. As shown in Fig. 10, stellar luminosity produces an ionisation of the surrounding medium up to far smaller distances than those covered by the GRB ionisation, especially for metals (i.e. elements heavier than He), which determine the gas spectroscopic appearance in the X-ray band. This is due both to the soft SED of stars, which can be described as black-body radiation, and thus is rather steeper than the typical GRB power-law SED, and to the much lower luminosity (around 18 orders of magnitude for the physical settings of this paper). As a result, the impact of the initial temperature and abundances is negligible (see Appendix B for details) and the absorption spectra are totally indistinguishable.

5 Discussion

5.1 Current limitations of the code and ongoing upgrades

The main limitation of the current version of TEPID is the assumption that all ions lie in their ground states at all times. In photoionised gas, levels above the ground with ΔE ~ kT become significantly and stably populated at densities higher than those for which collisional excitation rates are comparable to spontaneous de-excitation rates (e.g. Netzer 2013; Kallman et al. 2021). For permitted transitions this implies densities n_e ≫ 10¹⁰ cm⁻³, at which proper time-evolving photoionisation of the intervening gas becomes less and less critical. However, we note that time-evolving effects at such high densities could manifest when the timescale of the luminosity variability is very short, as in the case, for example, of X-ray lags and quasi-periodic oscillations in AGNs and compact sources (see e.g. Kara et al. 2016; Wang et al. 2022). We refer to García et al. (2013) for further discussion on this point.

The lack of a time-evolving electronic level population treatment prevents the calculation of diffuse and reflected line emission (and so an accurate radiative transfer treatment and the prediction of the gas emission spectrum), and also affects the correct ion balance computation, and thus, in turn, the heating–cooling balance. This explains most of the differences seen in Fig. 1 between the equilibrium temperature computed with TEPID and those derived from Cloudy and XSTAR at log(U) ≃ 1–3, together with the further improvement required to handle near-neutral gas states.

Our approximated radiative transfer considers outward photoelectric absorption of the incident continuum throughout the gas cloud. Neither the gas diffuse continua nor the reflected continua are computed. They are only important when the diffuse and reflected radiation fields contribute significantly to the total emitted radiation. In the GRB scenario the incident ionising continuum largely dominates, as the ISM is quickly photoionised by the incoming photons during the prompt and the early afterglow phases (i.e. from hundreds to thousands of seconds; see Fig. 14), while recombination continuum within the circumburst medium is produced over much longer timescales even in the densest ISM and star-forming regions.

To assess the impact of the lack of diffuse spectrum on the ionisation balance in the AGN scenario, we ran TEPID at equilibrium (as discussed in Sect. 3.1) and compared the resulting ion abundances as a function of N_H with those predicted by Cloudy in the range 0 < log(U) < 1.5. Figure 18, top panel, shows, for each metal included in TEPID and as a function of log(U), the maximum N_H up to which the agreement between the ion abundances computed with TEPID and Cloudy is within 50%. To perform a meaningful comparison, for each value of log(U) we limited the analysis to the most abundant ions (i.e. those with abundance ≥0.1). As expected, the agreement improves for increasing ionisation since the gas opacity, and thus its emission, decreases. The mismatch between the Cloudy- and TEPID-computed abundances for Fe and S, which is significantly higher than all the other metals, is due to the lack of updated recombination rates for low ionisation levels (from Fe 1 to 11 and S 1), for which Cloudy uses custom-computed means, while we employed the latest available rates from Mazzotta et al. (1998). To give an idea of the scatter between different photoionisation equilibrium codes, the bottom panel of Fig. 18 shows the limit N_H obtained by comparing Cloudy and XSTAR ion abundances. As can be seen, the agreement is roughly of the same order as that between TEPID and Cloudy.

It should also be noted that we neglect photon scattering, which leads to a sizeable underestimate of the gas opacity for hydrogen-equivalent column densities above 1 /σ_T = 1.5 × 10²⁴ cm⁻², where σ_T is the Thompson cross-section (see Appendix D for further discussion on this point).

Finally, in the current version of TEPID we do not consider free-free absorption within the cloud. This is a significant source of heating at low ionisation degrees and temperatures, and explains the cooler gas equilibrium temperatures predicted by TEPID at log(U) ≲ − 1 compared to Cloudy and XSTAR (Fig. 1).

The next version of the code will self-consistently solve the level population equations and this, in turn, will allow us to implement a physically self-consistent time-evolving radiative transfer by following the approach outlined in van Adelsberg & Perna (2013), and to compute the emerging absorption and emission spectra as a function of time, geometry and dynamics of the gas.

We also note that, for time-dependent effects to be detected in observed data, the ionised gas will have to spend a significant fraction of the observing time of each (time-resolved) spectral slice out of photoionisation equilibrium. The fulfilment of this condition strongly depends both on the observed source properties (variability amplitude and timescale of the ionising continuum, density and ionisation of the absorbing gas) and on the signal-to-noise ratio of the observation, and thus it requires an accurate evaluation on a case-by-case basis.

However, we note that for GRB circumburst environments the gas recombination tends to be negligible with respect to ionisation (as discussed in Sect. 4.4 and Appendix C.1), and thus the transmitted spectrum is mainly shaped by the progressive reduction of the gas opacity over time (see Fig. 14). As we show in Appendix E, this leads to absorption spectra that are dramatically different from those due to equilibrium absorbers, even when integrated over the typical XMM-Newton observing times of GRB afterglows.

Fig. 18 Comparison between TEPID and equilibrium photoionisation codes. Top: maximum value of NH up to which the agreement between TEPID and Cloudy is within 50% as a function of log(U). The colour-coding corresponds to the different metals (see legend). Bottom: same comparison as above, but between Cloudy- and XSTAR-predicted abundances.

5.2 Comparison with other time-evolving photoionisation codes

The upcoming advent of a new generation of X-ray spectrometers (see below) has increased the need for accurate modelling. Over the past few months two time-evolving ionisation codes similar to TEPID have been presented by Sadaula et al. (2023) and Rogantini et al. (2022), optimised on the physics of ionised absorbers in AGNs. As TEPID does, the two codes solve ion and energy balances as a function of time, neglecting the calculation of diffuse and reflected continuum and line emission, and compute time-resolved absorbed X-ray spectra. The code of Sadaula et al. (2023) performs an approximated radiative transfer treatment similar to ours, while that of Rogantini et al. (2022) considers only optically thin gas clouds.

With respect to the above codes, TEPID is probably faster as it uses the in-flight adaptive temporal resolution algorithm described in Sect. 2.4. This reduces the computing time to an affordable level, both when the ionising source SED varies significantly over short timescales (e.g. variable AGNs and X-Ray Binaries; see Figs. 6 and 8) and when the gas is geometrically thick (i.e. Δr/r ≳ 1; see Sect. 2.1), such as for the circumburst medium around GRBs (see Fig. 14).

We compared TEPID and TPHO, the code of Rogantini et al. (2022), by reproducing their ‘flaring light curve’ case (see their Sect. 3.4), in which an optically thin gas, initially in photoionisation equilibrium, is exposed to a factor of 10 increase in the incident ionising radiation within 20 ks. We obtain a remarkable agreement between the ionic abundances, temperature, and energy balance as a function of time, as shown in detail in Appendix F.

5.3 Time-resolved spectral analysis with future X-ray spectrometers

The unprecedented energy resolution and sensitivity of the microcalorimeters on board the next X-ray missions XRISM and Athena (XRISM Science Team 2020; Barret et al. 2023; Piro et al. 2022) will revolutionise our understanding of the gaseous environment of AGNs, compact sources, and X-ray transients. Observations of gas undergoing photoionisation by variable sources will require appropriate time-evolving codes to be properly analysed and to meaningfully constrain the physical parameters of the gas. As shown in Fig. 17, time-evolving ionisation leads to unique ion abundance patterns that cannot be mimicked by any equilibrium code. Fitting single-epoch absorption spectra of a time-variable warm absorber with equilibrium models may lead to erroneous diagnostics (e.g. Nicastro et al. 1999, Krongold et al. 2007) and even suggest the presence of multiple gas layers with different ionisation parameters where only one non-equilibrium gas component is actually present (Rogantini et al. 2022).

6 Conclusions

In this paper, we present the new Time Evolving PhotoIonisation Device (TEPID). TEPID is a highly modular and flexible time-evolving photoionisation code that produces time-resolved absorption spectra of gas illuminated by variable ionising astro-physical sources. TEPID’s time-resolved spectra can be stored in table models to be used within spectral fitting packages (such as xspec) to fit observed data. We applied TEPID to the case of AGN ionised absorbers and the circumburst medium around a GRB. Our main findings are the following:

• The equilibration time t_eq (Eq. (1)) is inversely proportional to the free electron density n_e, and dictates the timescale over which the gas readjust following a variation of the incident ionising luminosity.

• For ionising flux variations typical of AGNs and total (free + bound) electron densities $n_{\text{e}}^{\text{tot}}\gtrsim {10}^{10}c{m}^{-3}$ (see Figs. 4, 5 and 6), the absorbing gas reacts quickly and is close to photoionisation equilibrium at all times. In the highest density case, $n_{\text{e}}^{\text{tot}}={10}^{12}$ cm⁻³, gas temperature and ion abundances are practically in instantaneous equilibrium with the ionising flux. For decreasing $n_{\text{e}}^{\text{tot}}$ (i.e. increasing t_eq), the gas is increasingly far from photoionisation equilibrium and its reaction to flux variations is both delayed and compressed. For the lowest explored density, $n_{\text{e}}^{\text{tot}}={10}^4$ cm⁻³, the gas is insensitive to flux variations.

• The radial temperature and ionisation profiles of a star-forming region illuminated by a powerful GRB explosion strongly depend on the density of the circumburst medium (see Figs. 12 and 14). For all sampled densities, ${10}^{-1}\le n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}\le {10}^4$, the resulting radial ionisation profile is highly structured and always far from what would be expected if the gas were instantaneously in photoionisation equilibrium with the ionising radiation field (see Fig. C.1).

• In both the AGN and the GRB scenarios the emerging X-ray absorption spectra are uniquely associated with the particular value of the gas electron density $n_{\text{e}}^{\text{tot}}$ (see Figs. 9, 15 and 16), which can therefore be tightly constrained through time-resolved spectroscopy. While preliminary and pioneering studies have been performed on low-resolution X-ray observations of such targets (e.g. N99, Krongold et al. 2007, K13), future high-resolution and high-throughput X-ray spectrometers (e.g. the Resolve spectrometer of XRISM and the Athena X-IFU) will allow time-resolved spectroscopy to be systematically extended to sizeable samples of GRB afterglows and outflows from AGNs and X-Ray binaries, enabling a detailed physical characterisation of the ionised gas under a number of different astrophysical scenarios (XRISM Science Team 2020; Piro et al. 2022; see Juráňová et al. 2022 for a feasibility study with the Athena X-IFU microcalorimeter).

Applications of TEPID to both GRB afterglow and AGN X-ray data will be presented in a series of follow-up papers.

Acknowledgements

We thank the referee for their valuable comments, which significantly improved the paper. A.L., F.N., and L.P. acknowledge support from the HORIZON-2020 grant “Integrated Activities for the High Energy Astrophysics Domain” (AHEAD-2020), G.A. 871158. We made use of the NUMPY and MATPLOTLIB packages (Harris et al. 2020; Hunter 2007).

Appendix A Radiation transfer

Figure A.1 shows the total gas opacity τ of the circumburst environment around a GRB with $\text{log}\left(n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}\right)=4$, t = 10⁴s, and a total hydrogen-equivalent column density N_H = 10²³cm⁻². The red and green lines are computed ignoring the continuum absorption within each shell (i.e. using F_m instead of F′_m; see Eq. 8) and setting a linear spatial resolution with steps of ${\text{Δ}}_{N_{\text{H}}}={10}^{19}$,10²⁰cm⁻². In this setting the shell ionisation is simply proportional to the incident spectrum (F_m), and thus it is overestimated since the intra-shell absorption is not accounted for. This, in turn, leads to lower absorption when computing the transmitted spectrum $F_{m+1}=F_m\cdot e^{-{\tau }_m}$. This effect sums up throughout the total gas column, leading to a more ionised medium for thicker slabs: τ is indeed lower for higher ${\text{Δ}}_{N_{\text{H}}}$.

Fig. A.1 Gas opacity τ as a function of energy for the circumburst environment of a GRB for NH = 1023cm−2. The green and red lines correspond to spatial resolution ${\text{Δ}}_{N_{\text{H}}}={10}^{19}$,1020cm−2, respectively, and do not account for the gas absorption within each shell. The blue and purple lines (superimposed throughout most of the energy range) correspond to resolutions of ${\text{Δ}}_{N_{\text{H}}}={10}^{20}$,1021 cm−2 and account for absorption following Eq. 8.

The blue and purple lines, instead, are computed with F′_m (as in Eq. 8) and linear ${\text{Δ}}_{N_{\text{H}}}={10}^{20}$,10²¹cm⁻², respectively. It can be seen that absorption is more accurately accounted for, and thus the resulting τ is higher than in the previous cases. This demonstrates that Eq. 8 follows the propagation of the radiation with greater accuracy, and thus allows ${\text{Δ}}_{N_{\text{H}}}$ to increase up to two orders of magnitude (thus reducing considerably the computing time). Moreover, τ is practically identical for both values of ${\text{Δ}}_{N_{\text{H}}}$, showing that a resolution of 10²¹ cm⁻² is enough to properly sample the gas column.

Appendix B Impact of the pre-ionisation on the GRB case

Here we compare the gas temperature and ion abundances around the GRB obtained under two different assumptions concerning the initial ionisation of the gas: i) a Wolf–Rayet-like ionisation, as described in Sect. 4.1 and adopted throughout the paper, and ii) a totally neutral medium, with a typical ISM temperature of 100 K (Draine 2011). Figure B.1 reports T, H, Λ, Θ, and the ion abundances of O and Fe. The solid and dashed lines correspond to Wolf–Rayet pre-ionisation and initially neutral medium, respectively. In all cases t = 10⁴s, $log\left(n_{\text{e}}^{\text{tot}}\right)=2$. The differences between the two cases are practically negligible since the huge GRB luminosity is able to ionise up to far higher distances than the area affected by the progenitor star, and therefore the resulting ionisation is weakly dependent on the initial conditions.

Fig. B.1 Temperature (top left), energy balance (top right), and oxygen and iron ion abundances (bottom left and right) for t = 104s, $log\left(n_{\text{e}}^{\text{tot}}\right)=2$ assuming both a Wolf–Rayet pre-ionisation and an initially neutral medium (solid and dashed lines, respectively).

Appendix C Oxygen ionic stages

Figure C.1 shows the ratio of the ion abundances of O IX to O VIII (blue solid lines) and of O II to O I (green solid lines).

For comparison, we plot with dashed lines the values expected for a gas in photoionisation equilibrium, which are given by the relation $$\frac{n_{{\text{X}}^{i=1}}}{n_{{\text{X}}^i}}=\frac{F_{{\text{X}}^i}}{{\alpha }_{\text{rec}}\left({\text{X}}^i\right)}$$(C.1)

The ratios are in photoionisation equilibrium only within the innermost hot ionised region, where most of the oxygen is in the O IX and O VIII stage (Fig. 14), while at greater distances the gas is underionised, and thus its ionisation increases with time.

Fig. C.1 Ratio between selected ionic stages of oxygen (O IX/O VIII, blue solid line; O II/O I, green solid line). For comparison, the dashed lines show the ratios expected for a gas in photoionisation equilibrium. The top(bottom) row corresponds to $\text{log}\left(n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}\right)=1\left(=4\right)$, while from left to right t = 102,104,106s.

Appendix D Thompson scattering

The main effect of the electron scattering on our computations is to deviate photons away from the line of sight and deplete the transmitted spectrum. To model this, we ran a simulation in which we increase the gas opacity by using a modified opacity τ′(v) = τ(v) + N_H • σ_T, where σ_T is the Thompson cross-section and τ(v) is the opacity defined in Eq. 7. This framework, which roughly mimics the ‘open geometry’ setting of equilibrium photoionisation codes such as Cloudy and XSTAR, maximises the opacity due to electron scattering since it neglects photons scattered from different lines of sight towards ours. We find very small variations in both the ion abundances and the overall gas opacity. Figure D.1, top panel, compares the ion abundances in the circumburst environment of a GRB obtained using τ (i.e. ignoring Thompson scattering, solid lines) and τ′ (dotted lines). It can be seen that the differences are very small and are mostly concentrated at N_H > 5 • 10²³cm⁻². In all cases, n_e = 10⁴cm⁻³, t = 10⁴s after the burst.

Figure D.1, bottom panel, shows the gas opacity for N_H = 10²⁴ cm⁻² computed using τ and τ′ (red and green lines, respectively). The two opacities overlap in most of the energy band. At low energy (E = 10–20 eV) there is a small difference due to the slightly lower gas ionisation in the second case (see top panel). At E ≳ 10 keV, instead, the difference is due to the additional Thompson opacity, whose value is N_H • σ_T = 10²⁴ • 6.65/10²⁵ = 0.6.

Fig. D.1 Impact of Thompson scattering. Top: From left to right, selected carbon, oxygen, and iron ion abundances (see legends) as a function of NH, computed using the standard opacity τ and the Thompson-enhanced τ′ (solid and dotted lines, respectively). Bottom: Gas opacity for NH = 1024cm−2 as a function of the energy using τ and τ′ (red and green line, respectively).

Appendix E GRB afterglow spectra

Figure E.1 shows the time-integrated TEPID spectrum (in red) of a GRB environment with log(N_H/cm⁻²) = 21, log(n/cm⁻³) = 1 between t=20 and 80 ks, corresponding to a typical XMM-Newton pointing of a GRB afterglow (A. Thakur et al., in prep.). For comparison, the green lines correspond to photoionisation equilibrium solutions computed with the standard version of PHASE (Krongold et al. 2003) for log(U) between -2 and 0 (see legend for colour-coding). In all cases, the underlying continuum is a power law with photon index Γ = 2 and the same normalisation.

It can be seen that also when accounting for the (finite) duration of a typical X-ray observation, the time-evolving ionisation leads to a unique absorption spectrum that cannot be mimicked by the equilibrium solutions, in quite good agreement with what is shown in Fig. 17 for a given instant of time (which would correspond to an ‘instantaneous’ observation with a duration of zero seconds).

Fig. E.1 TEPID time-averaged spectrum for log(NH/cm−2) = 21, log(n/cm−3) = 1 between t=20 and t=80 ks (red lines). The green spectra correspond to equilibrium absorbers, with log(NH/cm−2) = 21 and increasing ionisation parameter log(U) (see legend). The left panel shows the entire 0.1–10 kev X-ray band, while the right panel is a zoom-in on the 0.5–0.7 keV range.

Appendix F Comparison with the TPHO code

Figure F.1 shows the ionic concentrations of oxygen VIII, IX (left) and iron XIX, XX (right), to be compared with those from TPHO reported in Fig. 7 in Rogantini et al. (2022). It should be noted that ionic concentrations are defined in their paper as the product between the ionic abundance and the abundance of that element relative to hydrogen. The thinnest to thickest lines correspond to increasing number density, from $\text{log}\left(n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}\right)=4$ to 10 in steps of 0.5. Figure F.2 shows the temperature (top) and heating and cooling rates (bottom) for the same density range (darkest to lightest line), to be compared with their Fig. 8.

Notwithstanding the different atomic libraries (TEPID is based on the Cloudy database, while TPHO is based on SPEX) and the different structure of the codes, we obtain a remarkable agreement between the ionic abundances, temperature and heating and cooling rates as a function of time.

Fig. F.1 Oxygen VIII, IX (left) and iron XIX, XX (right) ionic concentrations as a function of time.

Fig. F.2 Temperature (left) and heating and cooling rates (right, solid and dashed lines respectively).

References

All Figures

	Fig. 1 Photoionisation equilibrium temperature in an optically thin cloud of gas ionised by the standard AGN SED of Cloudy as a function of the ionisation parameter U (top) and the U/T ratio (bottom). The green, blue and red curves correspond to the values computed with Cloudy, XSTAR and TEPID, respectively.
In the text

	Fig. 2 Adaptive temporal binning by TEPID. Top panel: step-function light curve. Middle: Sampling frequency ω. Bottom: zoomed-in images of ω around the flux transition phases (see text for details).
In the text

	Fig. 3 AGN light curve.
In the text

Fig. 4 Heating H, cooling Λ, and Compton net heating Θ (red, blue, and green lines, respectively) in units of Ks−1 as a function of t. The black line represents the algebraic sum of these terms, corresponding to the time derivative of the temperature. The top and middle panels correspond to log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 8,12, while the bottom panel is a zoomed-in image of the temperature derivative for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 12.

In the text

	Fig. 5 T as a function of time for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) from 4 to 12 (lightest to darkest line).
In the text

	Fig. 6 Selected ion abundances for carbon, oxygen, and iron (from top to bottom) as a function of time for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) from 4 to 12 (thinnest to thickest line). The different panels correspond to different ions, as labelled. The y-scale is different for each panel.
In the text

	Fig. 7 Selected ion abundances (from top to bottom and left to right, carbon, oxygen, iron; see legends) and temperature (bottom right) as a function of NH for a gas in photoionisation equilibrium, computed with TEPID (solid lines), Cloudy and XSTAR (dashed and dotted lines, respectively).
In the text

Fig. 8 Time-dependent optical depth. Top: AGN light curve. The thick lines show the times for which ion abundances are shown (see below). The thin lines show the times corresponding to the absorption spectra in Fig. 9. Bottom, first and second row: Selected ion abundances (see legend) respectively for oxygen and iron as a function of NH. The left and right panels correspond to t = 5, 16 ks. log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 10 in all cases.

In the text

Fig. 9 Spectra of AGN outflows. Top: selected iron ion abundances (see legend) as a function of t for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6, 10 (left and right panels). The dashed lines show the times corresponding to the spectra shown below. Bottom: absorption spectra for log($N_{\text{H}}/{\text{cm}}^{-2}$) = 22 and t = 0, 2, 8, 16 ks (top to bottom, corresponding to the times indicated above and in Fig. 8) and log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 6, 10 (blue and orange lines). The left and right columns show respectively the Fe L and Fe K bands (i.e. 0.8–1.1 and 6.2–7.2 keV).

In the text

	Fig. 10 Pre-burst ionisation structure of the surroundings of a Wolf–Rayet star at the centre of a star-forming region with radius 100 pc and constant density $n_{\text{e}}^{\text{tot}}={10}^2{\text{cm}}^{-3}$. From left to right, the panels show C, O and Fe ion abundances and temperature as a function of the distance (in pc) from the star.
In the text

	Fig. 11 Luminosity and temporal resolution as a function of time in the GRB case. Top: GRB light curve. t1, t2, and tend correspond respectively to the end of the initial constant luminosity phase, the start of the afterglow and the end time of the simulation. Bottom: resolution ω as a function of time.
In the text

	Fig. 12 Temperature for log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1 as a function of the radius r (i.e. the distance from the GRB) for t = 102,104,106 s (dotted, dashed and solid lines, see legend).
In the text

Fig. 13 Contribution of hydrogen and helium to the total heating and cooling rates. Top panels: heating (left) and cooling (right) rates as a function of the radius r. The red line in the left panel and the blue line in the right panel correspond to the total rates, while the green and yellow shaded areas correspond to the total hydrogen and helium rates, respectively. The green(yellow) solid and dashed lines correspond to the rates of H(He) I, II, respectively, while the yellow dot-dashed line to He III. Bottom: hydrogen and helium (left and right, respectively) ion abundances (see legend for colour-coding) as a function of r. In all cases log($n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}$) = 1, t = 102 s.

In the text

	Fig. 14 Abundances for some oxygen (top) and iron (bottom) ions (see legend for colour-coding) as a function of the radius r. The dotted, dashed and solid lines identify increasing times from the onset of the GRB, t = 102,104,106 s, respectively, while the left, middle and right panels correspond to increasing ISM density, log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = −1, 1, 4.
In the text

	Fig. 15 Absorption spectra for a distance r = 100 pc from the GRB and increasing $n_{\text{e}}^{\text{tot}}$ (and thus increasing NH; see legend for colour-coding). Left: broad-band 0.1–10 keV spectra. Right: zoomed-in image of the 0.5–0.7 keV range, showing several oxygen absorption lines from neutral to hydrogen-like states.
In the text

	Fig. 16 Absorption spectra for column density NH = 1022 cm−2 and increasing $n_{\text{e}}^{\text{tot}}$ (and thus decreasing distance from the GRB; see legend for colour-coding). Left: broad-band 0.1–10 keV spectra. Right: zoomed-in image of the 0.5–0.7 keV range.
In the text

	Fig. 17 Relative abundances of the oxygen ion stages. Red line: TEPID abundances for log(${n_{\text{e}}^{\text{tot}}/\text{cm}}^{-3}$) = 1, t = 104 s averaged up to log(NH/cm−2) = 21. Blue to green lines: Photoionisation equilibrium abundances for log(U) from −3 to 0, averaged up to log(NH/cm−2) = 21.
In the text

	Fig. 18 Comparison between TEPID and equilibrium photoionisation codes. Top: maximum value of NH up to which the agreement between TEPID and Cloudy is within 50% as a function of log(U). The colour-coding corresponds to the different metals (see legend). Bottom: same comparison as above, but between Cloudy- and XSTAR-predicted abundances.
In the text

Fig. A.1 Gas opacity τ as a function of energy for the circumburst environment of a GRB for NH = 1023cm−2. The green and red lines correspond to spatial resolution ${\text{Δ}}_{N_{\text{H}}}={10}^{19}$,1020cm−2, respectively, and do not account for the gas absorption within each shell. The blue and purple lines (superimposed throughout most of the energy range) correspond to resolutions of ${\text{Δ}}_{N_{\text{H}}}={10}^{20}$,1021 cm−2 and account for absorption following Eq. 8.

In the text

	Fig. B.1 Temperature (top left), energy balance (top right), and oxygen and iron ion abundances (bottom left and right) for t = 104s, $log\left(n_{\text{e}}^{\text{tot}}\right)=2$ assuming both a Wolf–Rayet pre-ionisation and an initially neutral medium (solid and dashed lines, respectively).
In the text

	Fig. C.1 Ratio between selected ionic stages of oxygen (O IX/O VIII, blue solid line; O II/O I, green solid line). For comparison, the dashed lines show the ratios expected for a gas in photoionisation equilibrium. The top(bottom) row corresponds to $\text{log}\left(n_{\text{e}}^{\text{tot}}/{\text{cm}}^{-3}\right)=1\left(=4\right)$, while from left to right t = 102,104,106s.
In the text

	Fig. D.1 Impact of Thompson scattering. Top: From left to right, selected carbon, oxygen, and iron ion abundances (see legends) as a function of NH, computed using the standard opacity τ and the Thompson-enhanced τ′ (solid and dotted lines, respectively). Bottom: Gas opacity for NH = 1024cm−2 as a function of the energy using τ and τ′ (red and green line, respectively).
In the text

	Fig. E.1 TEPID time-averaged spectrum for log(NH/cm−2) = 21, log(n/cm−3) = 1 between t=20 and t=80 ks (red lines). The green spectra correspond to equilibrium absorbers, with log(NH/cm−2) = 21 and increasing ionisation parameter log(U) (see legend). The left panel shows the entire 0.1–10 kev X-ray band, while the right panel is a zoom-in on the 0.5–0.7 keV range.
In the text

	Fig. F.1 Oxygen VIII, IX (left) and iron XIX, XX (right) ionic concentrations as a function of time.
In the text

	Fig. F.2 Temperature (left) and heating and cooling rates (right, solid and dashed lines respectively).
In the text