Time Evolving Photo Ionisation Device (TEPID): a novel code for out-of-equilibrium gas ionisation

Photoionisation is one of the main mechanisms at work in the gaseous environment of bright astrophysical sources. Many information on the gas physics, chemistry and kinematics, as well as on the ionising source itself, can be gathered through optical to X-ray spectroscopy. While several public time equilibrium photoionisation codes are readily available and can be used to infer average gas properties at equilibrium, time-evolving photoionisation models have only very recently started to become available. They are needed when the ionising source varies faster than the typical gas equilibration timescale. Indeed, using equilibrium models to analyse spectra of non-equilibrium gas may lead to inaccurate results and prevents a solid assessment of the gas density, physics and geometry. We present our novel Time-Evolving PhotoIonisation Device (TEPID), which self-consistently solves time evolving photoionisation equations (thermal and ionisation balance) and follows the response of the gas to changes of the ionising source. The code can be applied to a variety of astrophysical scenarios and produces time-resolved gas absorption spectra to fit the data. To describe the main features of TEPID, we apply it to two dramatically different astrophysical scenarios: a typical ionised absorber observed in the X-ray spectra of Active Galactic Nuclei (e.g. Warm Absorbers and UFOs) and the circumburst environment of a Gamma-Ray Burst. In both cases, the gas energy and ionisation balances vary as a function of time, gas density and distance from the ionising source. Time evolving ionisation leads to unique ionisation patterns which cannot be reproduced by stationary codes when the gas is out of equilibrium. This demonstrates the need for codes such as TEPID in view of the up-coming high-resolution X-ray spectrometers onboard missions like XRISM or Athena.


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 2 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 1 ).A few important examples are the gaseous environment of active galactic nuclei (AGNs) and accreting compact sources, 1 A widely used alternative definition for the ionisation parameter is ξ = L ion /n tot e r 2 , where L ion and n tot e are the ionising luminosity and the total electron density (Tarter et al. 1969).
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 determined2 by the value of U. From an observational measurement of U, Q ion it is possible to derive an estimate of the product r 2 • 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) where α rec and n 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 timeresolved 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, metallicity and distance from the burst (see e.g.D 'Elia et al. 2009;Heintz et al. 2018).These results were based on ad hoc timeevolving level-population models built under the UV photopumping 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 firstorder 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.A141, page 2 of 24 Luminari, A., et al.: A&A, 679, A141 (2023)

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 selfconsistently 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 tot e , the total (bound + free) electron density; -the spectral energy distribution (SED).TEPID works in the spectral range 10 −5 -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 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 X i are constant throughout the shell.When this condition is not fulfilled (i.e.∆r/r ≳ 1), the input T, n 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 dt = t − t i prop , where t i prop 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 tot prop (where t tot prop 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 dt = t − t i prop .The time coordinate in each shell is set by imposing dt = 0 at t = 0 (i.e.correcting for t i prop ).For dt < 0 (i.e.t < t i prop ), the ith shell has not yet been reached by the first photons of the timeevolving 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).

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 −2 ) = 18) of gas illuminated by the standard AGN SED tabulated in Cloudy3 .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)4 .

Equations
The temporal evolution of the relative abundance of the ith ion of a given element X, n X i , can be written as (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 X i and C 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 Silike); 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 X i and the element abundance Z X , plus a Compton energy exchange term Θ: The heating rate H(X i ) consists of two terms: The first term is the heating from photoionised electrons (see e.g.Netzer 2013).The quantity F(ν)/hν is the incoming photon flux and it is linked to the ionising flux as F ion = 100 keV 10 −5 keV F(ν)dν, while σ X i (ν) is the ion photoelectric cross-section as a function of energy and (hν − hν 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ν X i up to the highenergy 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) 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 where K X is the atomic number and Z X the element abundance.
Finally, the frequency-dependent optical depth within the mth shell, τ m (ν) is given by 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 H = n tot e /1.2.Ionisation and heating rates in the mth shell are computed using the spectrum where F m is the incident spectrum at the innermost face of the shell.The spectrum F ′ 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 ′ 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 • e −τ m ) corrected A141, page 4 of 24 for geometric dilution (i.e. the squared ratio between the radial locations r m , r m+1 of the two slabs) 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.

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 ω of the time-integration sampling for two dramatically different volume electron densities, log(n tot e /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 i gross (where 1 ≤ n i gross ≤ N gross ), with width ∆t i gross , the equilibration time of each ion is evaluated, using the initial set of equilibrium ion abundances.The longest of these timescales t max,i eq (n i gross ) (typically << ∆t i gross ) is used within each n i gross interval to define finer sub-grids, each with N fine,i intervals with sizes ∆t j fine,i = Int(∆t i gross /t max,i eq ) × 2 j , where j = 0, N fine,i .-The total number of finer intervals from t start to t end is given by N grid = N gross i=1 N fine,i .This series of intervals defines the time-dependent temporal resolution ∆t j fine,i 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 ∆t j fine,i 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 timedependent sampling frequency ω i j = 1/∆t j fine,i , 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 ω i j increased by a factor of 10, and checking that all the results (i.e.temperature, ionic populations, transmitted spectra) are unaffected.

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 ∆t j fine,i 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 hydrogenequivalent column density N H (or outward radius r out ) is reached.The final outputs of the simulation are T, n X i , N X i , and N H as a function of time t.We use these quantities to build PHASE table models (Krongold et al. 2003)  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.

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 −2 ) = 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 1 = 4 ks, staying at F max from t 1 through t 2 = 6 ks and then decreasing again linearly (and with the same slope) to F min = 1 from t 2 through t 3 = 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 tot e /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 tot e /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 tot e /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 0 to t 1 ), then is = 0 during the maximum constant flux phase (from t 1 to t 2 where the temperature is constant and equal to 10 6 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 2 and t 3 ).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 tot e the evolution of the gas physical properties departs from equilibrium.Ion-by-ion photoionisation rates do not explicitly depend on n tot e ; however, for a given initial  log(U) the incident flux F ion , which is proportional to photoionisation rates, scales linearly with r 2 n H ∝ r 2 n tot e .Thus, at a given distance, lower values of n tot e 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 tot e .First, both recombination and photoionisation rates are slower, and thus the variations in T and n 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. A141, page 7 of 24 between t 1 and t 2 ), 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 tot e /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 tot e .

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. 2003Krongold et al. , 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 tot e /cm −3 ) = 4-12 and up to an equivalent hydrogen column density log(N H /cm −2 ) = 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.

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 −2 ) ≈ 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 −2 ).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.and the gas becomes optically thick only at log(N H /cm −2 ) ≳ 23.Consequently, the gas column is virtually isothermal at all times during the simulation (i.e.T (t 1 ) = 3 × 10 6 and T (t 2 ) = 10 6 ) and has constant light metal abundances (i.e.f CVII ≃ f OIX ≃ 1) up to log(N H /cm −2 ) ≃ 23 (middle panel of Fig. 8).

Time-dependent optical depth
At log(N H /cm −2 ) ≳ 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 tot e and on the ionising flux state.At the chosen log(n tot e /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 −2 ) ≳ 23.8 (right middle panel of Fig. 8), while the He-, Be-, and C-like fractional abundances of Fe start decreasing quickly already at log(N H /cm −2 ) ≳ 22.5 to the advantage of the Ne-and O-like ions (right bottom panel of Fig. 8).

Time-resolved spectra of AGN outflows
Figure 9 shows the absorption spectra extracted from the above simulation for log(N H /cm −2 ) = 22 and log(n tot e /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 A141, page 9 of 24 Fig. 9. Spectra of AGN outflows.Top: selected iron ion abundances (see legend) as a function of t for log(n tot e /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 H /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 tot e /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).therefore the spectra are independent of n tot e (orange superimposed on blue).
At high density (log(n tot e /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 tot e /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 A141, page 10 of 24 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 ≈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.

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.

Physical setting of the GRB surroundings
For simplicity we assume that the ISM surrounding the GRB has constant number density n tot e (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 tot e /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 −3 ), 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 5 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.

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 1 =85 s and decaying as t −2 between t 1 and t 2 = 100 s.At t 2 the GRB transits to the afterglow phase: the spectrum becomes steeper, α xag = 1.0, and the light curve flattens as t −1 .Figure 11, top panel, shows the light curve from t = 0 to t end = 10 6 , 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 6 s) phases are 1.4 × 10 51 erg and 5.4 × 10 51 erg, respectively.Given the shape A141, page 11 of 24 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 3 s −1 at t = 0 to 10 −3 s −1 at t = t end = 10 6 s.The only exception is between 85 and 100 s, where ω increases to better follow the temporal decay transition.

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 tot e /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 5 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) 2 = 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 tot e /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 tot e /cm −3 ) = −1, 1, 4. We include log(n tot e /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 A141, page 12 of 24 Luminari, A., et al.: A&A, 679, A141 (2023) 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 tot e /cm −3 ) = 1, t = 10 2 s. 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 highversus low-ionisation ions strongly increases with n tot e 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.
e /N H ) 2 ; in other words, at a given N H (enclosed in a radius r N H ) the ionising flux is ∝ (n tot e ) 2 , and thus less diluted with increasing density.As a result, the transmitted spectra will be different, and thus highresolution 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.

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 4 s and for increasing log(n tot e /cm −3 ), from -1 to 4 (and thus increasing log(N H /cm −2 ), 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 tot e 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 tot e /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 tot e /cm −3 ) and a fixed column density of log(N H /cm −2 ) = 22.This translates into a broad range of radii according to the electron density.For log(n tot e /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 tot e /cm −3 ) = 3).On the other hand, for log(n tot e /cm −3 ) = −1 the radius is 32 kpc and the absorbing gas is mostly neutral.
The dependence on n tot e is a peculiar feature of nonequilibrium 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 A141, page 13 of 24 Luminari, A., et al.: A&A, 679, A141 (2023) 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 = 10 2 , 10 4 , 10 6 s, respectively, while the left, middle and right panels correspond to increasing ISM density, log(n tot e /cm −3 ) = −1, 1, 4. 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 tot e /cm −3 ) = 1 and t = 10 4 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 −2 ) = 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 −2 ) = 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 tot e offers a valuable opportunity to directly constrain it in observed spectra by fitting afterglow spectra with timeevolving photoionisation models such as TEPID.

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 A141, page 15 of 24 Luminari, A., et al.: A&A, 679, A141 (2023) temperature and abundances is negligible (see Appendix B for details) and the absorption spectra are totally indistinguishable.

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 10 cm −3 , at which proper time-evolving photoionisation of the intervening gas becomes less and less critical.However, we note that timeevolving 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 nearneutral 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 24 cm −2 , 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.A141, page 16 of 24 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.

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.

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).

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 astrophysical 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 tot e ≳ 10 10 cm −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 tot e = 10 12 cm −3 , gas temperature and ion abundances are practically in instantaneous equilibrium with the ionising flux.For decreasing n tot e (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 tot e = 10 4 cm −3 , the gas is insensitive to flux variations.
-The radial temperature and ionisation profiles of a starforming 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 ≤ n tot e /cm −3 ≤ 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. -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 tot e (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.Figure E.1 shows the time-integrated TEPID spectrum (in red) of a GRB environment with log(N H /cm −2 ) = 21, log(n/cm −3 ) = 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.

Appendix E: GRB afterglow spectra
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 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 log(n tot e /cm −3 ) = 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. 1 .
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 .
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).

Fig. 4 .
Fig.4.Heating H, cooling Λ, and Compton net heating Θ (red, blue, and green lines, respectively) in units of K s −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 tot e /cm −3 ) = 8, 12, while the bottom panel is a zoomed-in image of the temperature derivative for log(n tot e /cm −3 ) = 12.

Fig. 6 .
Fig. 6.Selected ion abundances for carbon, oxygen, and iron (from top to bottom) as a function of time for log(n tot e /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 N H for a gas in photoionisation equilibrium, computed with TEPID (solid lines), Cloudy and XSTAR (dashed and dotted lines, respectively).

Figure 8 ,
Figure8, second and third rows, shows the ionic fractions of O and Fe as a function of N H for log(n tot e /cm −3 ) = 10 and at two different times during the simulations: t 1 = 5 ks and t 2 = 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,

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 N H .The left and right panels correspond to t = 5, 16 ks.log(n tot e /cm −3 ) = 10 in all cases.
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 tot e = 10 2 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.

Fig. 11 .
Fig. 11.Luminosity and temporal resolution as a function of time in the GRB case.Top: GRB light curve.t 1 , t 2 , and t end 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.

Fig. 15 .Fig. 17 .
Fig. 15.Absorption spectra for a distance r = 100 pc from the GRB and increasing n tot e (and thus increasing N H ; 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. 18 .
Fig. 18.Comparison between TEPID and equilibrium photoionisation codes.Top: maximum value of N H up to which the agreement between TEPID and Cloudy is within 50% as a function of log(U).The colourcoding corresponds to the different metals (see legend).Bottom: same comparison as above, but between Cloudy-and XSTAR-predicted abundances. C.1).

Fig. D. 1 .
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 N H , computed using the standard opacity τ and the Thompson-enhanced τ ′ (solid and dotted lines, respectively).Bottom: Gas opacity for N H = 10 24 cm −2 as a function of the energy using τ and τ ′ (red and green line, respectively).
Figure E.1 shows the time-integrated TEPID spectrum (in red) of a GRB environment with log(N H /cm −2 ) = 21, log(n/cm −3 ) = 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.17for a given instant of time (which would correspond to an 'instantaneous' observation with a duration of zero seconds).

Figure
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.7inRogantini 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 log(n tot e /cm −3 ) = 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
Fig. F.2. Temperature (left) and heating and cooling rates (right, solid and dashed lines respectively).