© The Authors 2024

1 Introduction

Medium-resolution spectroscopy of protoplanetary discs with the Spitzer Space Telescope has revealed a dense forest of H₂O emission lines (Carr & Najita 2008; Salyk et al. 2008), in addition to overlapping line emission features from OH, C₂H₂, HCN, and CO₂. From the emission temperatures of these molecules, which are of a few hundred Kelvin to about 1000 K, and the derived sizes of the line-emitting areas, it was immediately clear that these molecular spectra probe the chemistry in the terrestrial planet-forming regions; see for example the review by Pontoppidan et al. (2014). Over the past 15 yr, more than 100 Spitzer/IRS disc spectra have been analysed using simple slab-models, which has revealed certain properties of these line emissions in terms of molecular column densities, emission temperatures, and line emission areas; see Pontoppidan et al. (2010); Carr & Najita (2011); Salyk et al. (2011). Some statistics are also available, and certain trends have been identified concerning the influence of the accretion luminosity, the disc mass, the depletion of molecules in discs with large inner dust cavities, and chemical differences as a function of stellar type; see, for example, Pascucci et al. (2013); Najita et al. (2013); Walsh et al. (2015); Banzatti et al. (2017).

However, the significance of these measured molecular properties, where in the disc these molecules are present, both radially and vertically, and how this all fits together are open questions. The slab-models used for the analysis of the spectra do not provide such insights, and also have certain principle weaknesses, such as their inability to predict absorption lines, and an uncertainty as to how to include the effects of Keplerian broadening, continuous dust emission, and dust absorption under a given inclination angle.

Ultimately, we are seeking a higher level of understanding, and predictions of molecular emission features that are in accordance with the disc chemistry and the heating and cooling balance based on a model for the physical structure of the gas and dust in the disc that can predict the emission features of all molecules simultaneously. However, thermo-chemical disc models, such as those of Akimkin et al. (2013); Zhang et al. (2013); Walsh et al. (2015); Bruderer et al. (2015); Woitke et al. (2018a); Anderson et al. (2021), have yet to demonstrate that this is possible.

The puzzle of revealing the physical structure of the disc, which produces these emission lines, is complicated by the fact that we presently have little understanding of the element abundances of the gas in the inner disc. Indeed, there are very good reasons to assume that the transport of icy grains in the disc can profoundly change these element abundances during disc evolution (e.g. Ciesla & Cuzzi 2006). For example, pebble drift through the water snowline has been demonstrated in hydrodynamical disc models to increase the oxygen abundance in the inner disc (Bosman et al. 2017; Booth & Ilee 2019).

However, results from the Very Large Telescope (VLT)/MATISSE interferometer (Kokoulina et al. 2021) suggest that the dust in the inner 10 au of the disc around HD 179218 is carbon-enriched in the form of nano-carbon grains. These conclusions are based on the spectral shape of the solid opacities. Studying new VLT/MATISSE and VLT/GRAVITY data of HD 144432, Varga et al. (accepted by A&A) concluded that most of the carbon in the inner disc is in the gas phase, whereas solid iron and FeS are the main sources of dust opacity at wavelengths in the L-, M-, and N-bands. Additional clues about element transport processes have been obtained by linking the mid-infrared (MIR) line spectra to continuum images obtained with the Atacama Large Millimeter Array (ALMA). Expanding on earlier results from Najita et al. (2013), Banzatti et al. (2020) reported an anti-correlation between the infrared water line luminosity tracing the gas within a few astronomical units (au) and the distribution of solid pebbles at 10–200 au in discs. That is, radially compact discs show stronger water lines than large discs. Based on new James Webb Space Telescope (JWST) data, Banzatti et al. (2023) see more evidence for this anti-correlation. van ‘t Hoff et al. (2020) and Nazari et al. (2023) proposed that the inward transport and subsequent sublimation of carbonaceous dust grains and Polycyclic Aromatic Hydrocarbon molecules (PAHs) might enrich the inner disc with carbon.

New evidence for carbon-rich/oxygen-poor inner discs has recently been found in the first published disc spectra obtained with the James Webb Space Telescope (JWST); see Kóspál et al. (2023); Grant et al. (2023); Tabone et al. (2023); Banzatti et al. (2023); Perotti et al. (2023). A rich hydrocarbon chemistry was identified by Tabone et al. (2023) in the inner disc around the M5-star 2MASS-J16053215-1933159, with exceptionally optically thick C₂H₂ emission features around 7 µm and 14 µm, which are so strong that they form a wide quasi-continuum with tens of thousands of overlapping weak lines, without detecting any water lines. This is very clear evidence for a strong oxygen depletion in the inner disc of this star. Tabone et al. (2023) also mention a carbon-enrichment of the inner disc by the destruction of PAHs or carbonaceous grains as a possible cause for the large column densities of hydro-carbon molecules. However, this appears to be inconsistent with the very weak CO fundamental ro-vibrational lines as measured with JWST. To fit these lines, removal of the oxygen appears to be required, rather than enrichment of the carbon (Kanwar et al., in prep.).

In this paper, we take on the challenge of fitting the first publisheD JWST spectrum of a T Tauri star (Kóspál et al. 2023) using the 2D thermo-chemical disc model PRODIMO (Woitke et al. 2009, 2016). In Sect. 2, we report a number of important PRODIMO code changes that result in considerably improved predictions of the various MIR molecular emission features in general. In Sect. 3, we explain how we set up our disc model and how we fitted that model to the spectral energy distribution (SED) of EXLupi and its JWST/MIRI spectrum. Section 4 contains a thorough analysis of the best-fitting disc model in terms of the physical and chemical disc structure, including the calculated dust and gas temperatures. We also derive molecular column densities and emission temperatures from this model that can be compared to the values derived from slab-models, and analyse the chemical pathways in the line-forming regions of the various molecules in this section. Our conclusions are presented in Sect. 5.

2 Code improvements

Woitke et al. (2018a) demonstrated the problems faced by current standard thermo-chemical models for T Tauri discs in explaining the approximately equally strong emission features of H₂O, CO₂, HCN, and C₂H₂. First, a regular gas/dust mass ratio of 100 generally results in overly weak MIR molecular emissions. Second, while the molecular emission features of H₂O and OH can be pushed to the observed flux levels using a larger gas/dust ratio of 1000 (Bruderer et al. 2015; Woitke et al. 2018a; Greenwood et al. 2019), the CO₂ emission feature at 15 µm becomes too strong. Finally, the emission lines of HCN and C₂H₂ mostly come from a radially extended optically thin photodissociation region high above the disc in these models, but the emission lines created in these layers are generally optically thin, and their emission characteristics are in conflict with observations, in particular the molecular column densities are too low in these models; see for example Bruderer et al. (2015) for HCN. Kanwar et al. (2024) thoroughly analysed these dilute photodissociation regions using an extended hydro-carbon chemical network.

Increasing the C/O ratio can indeed boost the HCN and C₂H₂ line fluxes, but this leads to the loss of the CO₂ emission feature (Woitke et al. 2018a), and a very fine-tuned C/O value would be required to explain the ratios of the emission line features of different molecules. In contrast, what the observations tell us is that the emission features of H₂O, CO₂, HCN, and C₂H₂ are optically thick in most cases, with column densities in excess of 10¹⁵cm⁻², which causes the flux levels of these features to be on a similar level in a more robust way. We therefore carefully analysed our PRODIMO models to see how we can possibly create larger amounts of optically thick HCN and C₂H₂ molecules. These efforts resulted in a number of important code changes that are reported in this section. These changes are essential for fitting the first JWST spectrum, as described in Sect. 3.

2.1 A new escape probability theory for discs

Many published line fluxes from models of protoplanetary discs rely on the concept of escape probability; for example Nomura et al. (2007). Different variants of the escape probability formalism have been developed (Woods & Willacy 2009; Woitke et al. 2009; Du & Bergin 2014). In PRODIMO, while we normally perform formal solutions of the 3D line-transfer problem to calculate the line fluxes, an escape probability method is still required to compute the non-local thermodynamic equilibrium (non-LTE) populations and the heating and cooling rates. The previously used escape probability method in PRODIMO is based on Eq. (73) in Woitke et al. (2009), which is a rather crude approximation as it only considers radiative pumping by continuum photons in the radial direction and line photon escape into the upward vertical direction. In particular, the IR pumping from the disc below is incorrectly attenuated by the radial line optical depth ${\tau }_{ul}^{\text{rad}}$, which can be huge, which means we are likely underestimating the true IR pumping. In addition, in the outer optically thin disc regions, the escape probability is currently limited to $P_{ul}^{\text{esc}}\le \frac{1}{2}$, whereas it should be $P_{ul}^{\text{esc}}\to 1$ in the optically thin limit.

This paper introduces a new idea about how to calculate the line-averaged mean intensity ${\overline{J}}_{ul}$ by taking into account continuum radiative transfer effects in the line resonance region:

$${\overline{J}}_{ul}=\frac{1}{4\pi }{\displaystyle \int {\displaystyle \int {\varphi }_v{I}_v\left(n\right)d\text{Ω}dv}},$$(1)

where ϕ_v is the line profile function associated with the spectral line u → 1 and I_v is the spectral intensity. Finding ${\overline{J}}_{ul}$ is key to formulating the effective non-LTE rates and heating and cooling functions associated with line absorption and emission; see Eqs. (80), (81), (85), and (86) in Woitke et al. (2009). Here, we show that, by making certain assumptions about disc geometry and the locality of the radiative transfer quantities explained below, it is possible to obtain a final expression that again reads $${\overline{J}}_{ul}=P_{ul}^{\text{pump }}{J}_v^{\text{cont }}+\left(1-P_{ul}^{\text{esc }}\right)S_{\text{L}},$$(2)

but now with new pumping and escape probabilities, $P_{ul}^{\text{pump }}$ and $P_{ul}^{\text{esc}}$, that depend not only on line optical depths, but also on continuum optical depths, and take into account the pumping from and the escape in all directions. S_L is the line source function and $J_v^{\text{cont }}$ the continuous mean intensity that would result from the line having zero opacity. We write the continuum and line radiative transfer equation as $$\frac{d{I}_v}{ds}=-\left({\kappa }_{\text{C}}+{\kappa }_{\text{L}}{\varphi }_v\right)I_v+{\kappa }_{\text{C}}{S}_{\text{C}}+{\kappa }_{\text{L}}{\varphi }_v{S}_{\text{L}},$$(3)

where the continuum, absorption and scattering opacities [cm⁻¹], and the continuum source function, assuming isotropic scattering, are given by $${\kappa }_{\text{C}}={\kappa }_y^{\text{abs}}+{\kappa }_y^{\text{sca}},$$(4)

$${\kappa }_{\text{C}}{S}_{\text{C}}={\kappa }_v^{\text{abs}}{B}_v\left(T_{\text{d}}\right)+{\kappa }_v^{\text{sca}}{J}_v^{\text{cont}}.$$(5)

The line opacity [cm⁻¹ Hz], Gaussian line profile function [Hz⁻¹], and line source function [erg cm⁻² s⁻¹ Hz⁻¹ sr⁻¹] are $${\kappa }_{\text{L}}=\frac{h{v}_{ul}}{4\pi }\left(n_l{B}_{lu}-n_u{B}_{ul}\right),$$(6)

$${\varphi }_v=\frac{1}{\sqrt{\pi }\text{Δ}{\nu }_{\text{D}}}\exp \left(-{\left(\frac{v-v_{ul}}{\text{Δ}{v}_{\text{D}}}\right)}^2\right),$$(7)

$$S_{\text{L}}=\frac{2h{v}_{ul}^3}{c^2}{\left(\frac{g_u{n}_l}{g_l{n}_u}-1\right)}^{-1},$$(8)

where v_ul is the line centre frequency, v_D = v_u1 Δv/c is the Doppler width, and $\text{Δ}v={\left(v_{\text{th }}^2+v_{\text{turb }}^2\right)}^{1/2}$ is the line width in velocity space. Introducing the optical depth and maximum optical depth as $$t_v(s)={{{\displaystyle \int }}^{\text{}}}_0^s{\kappa }_{\text{C}}+{\kappa }_{\text{L}}{\varphi }_{v}ds,$$(9)

$${\tau }_v=t_v(s_{\text{max}}),$$(10)

the formal solution of Eq. (3) is $$I_v(n)=I_v^0(n){\text{e}}^{-{\tau }_v}+{{{\displaystyle \int }}^{\text{}}}_0^{s_{\max }}\left({\kappa }_{\text{C}}{S}_{\text{C}}+{\kappa }_{\text{L}}{\varphi }_v{S}_{\text{L}}\right){\text{e}}^{-t_v}ds,$$(11)

where n is the direction of the ray (unit vector) and s_max is the distance backward along the ray up to the point where the ray enters the model volume. $I_v^0(n)$ is the incident intensity at that point in direction n.

Assumption 1 (locality). We assume that all radiative transfer quantities, that is κ_C, κ_L, S_L, S_C, and ϕ_v, are independent of location and direction, and are given by their local values. Assuming S_L = const is standard in escape probability theory, because the line photons typically interact only within a small ‘resonance volume’ around the point of interest. Generalising this concept to the continuum transfer properties is a new idea (however, see Hummer & Rybicki 1985).

Using the above strategy, we can capture the most relevant continuum radiative transfer effects in that resonance region, where there is competition between line and continuum absorption and emission. Constant ϕ_v also means that we neglect line velocity shifts; that is, we assume a static medium. The only quantities that are allowed to depend on direction are $I_v^0(n)$ and s_max (and hence τ_v). The formal solution of the RT Eq. (11) then simplifies to $$I_v(n)=I_v^0(n){\text{e}}^{-{\tau }_{\text{C}}-{\tau }_{\text{L}}{\varphi }_v}+\frac{{\tau }_{\text{C}}{S}_{\text{C}}+{\tau }_{\text{L}}{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}+{\tau }_{\text{L}}{\varphi }_v}\left(1-{\text{e}}^{-{\tau }_{\text{C}}-{\tau }_{\text{L}}{\varphi }_v}\right),$$(12)

where we use the continuum and line centre optical depths τ_C = κ_Cs_max and τ_L= κ_L s_max.

Assumption 2 (disc geometry). We consider three major directions: radial (the rays coming from the stellar surface, which are often important for pumping), vertically upwards, and vertically downwards. At any given location in the disc, the star occupies a small solid angle Σ_*. The other two principle directions each cover half of Σ_d = 4π − Σ_*. All rays coming from the stellar surface are represented by a single radial ray. For the upward and downward directions, we approximate the disc as being plane-parallel with µ = cos(θ), where θ is the angle with the vertical. Thus, we find from Eq. (1) that $$\begin{array}{l}{\overline{J}}_{ul}=\frac{{\text{Ω}}_{\star }}{4\pi }{\displaystyle \int {\varphi }_v{I}_v^{\text{rad}}dv}+\frac{{\text{Ω}}_{\text{d}}}{8\pi }{\displaystyle \int {\varphi }_v}{{{\displaystyle \int }}^{\text{}}}_0^1{I}_v(\mu )d\mu dv\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }{\displaystyle \int {\varphi }_v}{{{\displaystyle \int }}^{\text{}}}_{-1}^0{I}_v\left(\mu \right)d\mu dv.\hfill \end{array}$$(13)

Introducing the vertically upwards and downwards (across the disc) line and continuum optical depths, ${\tau }_{\text{L}}^{\uparrow }$, ${\tau }_{\text{L}}^{\downarrow }$, ${\tau }_{\text{C}}^{\uparrow }$, and ${\tau }_{\text{C}}^{\downarrow }$, respectively, and the radially inward line and continuum optical depths, ${\tau }_{\text{L}}^{\text{rad}}$ and ${\tau }_{\text{C}}^{\text{rad}}$, and using plane-parallel geometry τ(µ) = τ(µ = 0)/µ, we obtain from Eqs. (12) and (13) that $$\begin{array}{l}{\overline{J}}_{ul}=\frac{{\text{Ω}}_{\star }}{4\pi }{{\displaystyle \int }}^{\text{}}{\varphi }_v(I_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}-{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v}\hfill \\ +\frac{{\tau }_{\text{C}}^{\text{rad}}{S}_{\text{C}}+{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}^{\text{rad}}+{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v}\left(1-{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}-{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v}\right))dv\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }{{\displaystyle \int }}^{\text{}}{\varphi }_v{{{\displaystyle \int }}^{\text{}}}_0^1{I}_v^{\text{ISM}}{\text{e}}^{-\frac{{\tau }_{\text{C}}^{\uparrow }}{\mu }-\frac{{\tau }_{\text{L}}^{\uparrow }}{\mu }{\varphi }_v}\hfill \\ +\frac{{\tau }_{\text{C}}^{\uparrow }{S}_{\text{C}}+{\tau }_{\text{L}}^{\uparrow }{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}^{\uparrow }+{\tau }_{\text{L}}^{\uparrow }{\varphi }_v}\left(1-{\text{e}}^{-\frac{{\tau }_{\text{C}}^{\uparrow }}{\mu }-\frac{{\tau }_{\text{L}}^{\uparrow }}{\mu }{\varphi }_v}\right)d\mu dv\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }{{\displaystyle \int }}^{\text{}}{\varphi }_v{{{\displaystyle \int }}^{\text{}}}_{-1}^0{I}_v^{\text{ISM}}{\text{e}}^{-\frac{{\tau }_{\text{C}}^{\downarrow }}{\left|\mu \right|}-\frac{{\tau }_{\text{L}}^{\downarrow }}{\left|\mu \right|}{\varphi }_v}\hfill \\ +\frac{{\tau }_{\text{C}}^{\downarrow }{S}_{\text{C}}+{\tau }_{\text{L}}^{\downarrow }{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}^{\downarrow }+{\tau }_{\text{L}}^{\downarrow }{\varphi }_v}\left(1-{\text{e}}^{-\frac{{\tau }_{\text{C}}^{\downarrow }}{\left|\mu \right|}-\frac{{\tau }_{\text{L}}^{\downarrow }}{\left|\mu \right|}{\varphi }_v}\right)d\mu dv,\hfill \end{array}$$(14)

where $I_v^{\star }$ and $I_v^{\text{ISM}}$ are the stellar and interstellar incident intensities, respectively. Introducing the second exponential integral function $E_2(\tau )={\displaystyle {\int }_0^1\exp }\left(-\frac{\tau }{\mu }\right)d\mu$, we find $$\begin{array}{l}{\overline{J}}_{ul}=\frac{{\text{Ω}}_{\star }}{4\pi }{{\displaystyle \int }}^{\text{}}{\varphi }_v(I_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}-{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v}\hfill \\ +\frac{{\tau }_{\text{C}}^{\text{rad}}{S}_{\text{C}}+{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}^{\text{rad}}+{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v}\left(1-{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}-{\tau }_{\text{L}}^{\text{rad}}{\varphi }_v}\right))dv\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }{{\displaystyle \int }}^{\text{}}{\varphi }_v(I_v^{\text{ISM}}{E}_2\left({\tau }_{\text{C}}^{\uparrow }+{\tau }_{\text{L}}^{\uparrow }{\varphi }_v\right)\hfill \\ +\frac{{\tau }_{\text{C}}^{\uparrow }{S}_{\text{C}}+{\tau }_{\text{L}}^{\uparrow }{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}^{\uparrow }+{\tau }_{\text{L}}^{\uparrow }{\varphi }_v}\left(1-E_2\left({\tau }_{\text{C}}^{\uparrow }+{\tau }_{\text{L}}^{\uparrow }{\varphi }_v\right)\right))dv\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }{{\displaystyle \int }}^{\text{}}{\varphi }_v(I_v^{\text{ISM}}{E}_2\left({\tau }_{\text{C}}^{\downarrow }+{\tau }_{\text{L}}^{\downarrow }{\varphi }_v\right)\hfill \\ +\frac{{\tau }_{\text{C}}^{\downarrow }{S}_{\text{C}}+{\tau }_{\text{L}}^{\downarrow }{\varphi }_v{S}_{\text{L}}}{{\tau }_{\text{C}}^{\downarrow }+{\tau }_{\text{L}}^{\downarrow }{\varphi }_v}\left(1-E_2\left({\tau }_{\text{C}}^{\downarrow }+{\tau }_{\text{L}}^{\downarrow }{\varphi }_v\right)\right))dv\hfill \end{array}$$(15)

We now introduce the dimensionless profile function $\varphi (x)=\exp \left(-x^2\right)/\sqrt{\pi }$ with $x=\left(v-v_{ul}\right)/\text{Δ}{v}_{\text{D}}$, and six basic 2D functions, which all produce values of between 0 and 1, to further simplify the result: $${\alpha }_0\left({\tau }_{\text{L}},{\tau }_{\text{C}}\right)={{\displaystyle \int }}^{\text{}}\varphi (x){\text{e}}^{-{\tau }_{\text{C}}-{\tau }_{\text{L}}\varphi (x)}dx,$$(16)

$${\alpha }_1\left({\tau }_{\text{L}},{\tau }_{\text{C}}\right)={{\displaystyle \int }}^{\text{}}\varphi (x)\frac{{\tau }_{\text{C}}}{{\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)}\left(1-{\text{e}}^{-{\tau }_{\text{C}}-{\tau }_{\text{L}}\varphi (x)}\right)dx,$$(17)

$${\alpha }_2\left({\tau }_{\text{L}},{\tau }_{\text{C}}\right)={{\displaystyle \int }}^{\text{}}\varphi (x)\frac{{\tau }_{\text{L}}\varphi (x)}{{\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)}\left(1-{\text{e}}^{-{\tau }_{\text{C}}-{\tau }_{\text{L}}\varphi (x)}\right)dx,$$(18)

$${\beta }_0\left({\tau }_{\text{L}},{\tau }_{\text{C}}\right)={{\displaystyle \int }}^{\text{}}\varphi (x)E_2\left({\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)\right)dx,$$(19)

$${\beta }_1\left({\tau }_{\text{L}},{\tau }_{\text{C}}\right)={{\displaystyle \int }}^{\text{}}\varphi (x)\frac{{\tau }_{\text{C}}}{{\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)}\left(1-E_2\left({\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)\right)\right)dx,$$(20)

$${\beta }_2\left({\tau }_{\text{L}},{\tau }_{\text{C}}\right)={{\displaystyle \int }}^{\text{}}\varphi (x)\frac{{\tau }_{\text{L}}\varphi (x)}{{\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)}\left(1-E_2\left({\tau }_{\text{C}}+{\tau }_{\text{L}}\varphi (x)\right)\right)dx,$$(21)

where, using the Einstein relations, $$\begin{array}{l}{\tau }_{\text{L}}={{\displaystyle \int }}^{\text{}}\frac{{\kappa }_{\text{L}}}{\text{Δ}{v}_{\text{D}}}ds={{\displaystyle \int }}^{\text{}}\frac{h{v}_{ul}}{4\pi \text{Δ}{v}_{\text{D}}}\left(n_l{B}_{lu}-n_u{B}_{ul}\right)ds\hfill \\ ={{\displaystyle \int }}^{\text{}}\frac{A_{ul}}{8\pi \text{Δ}v}{\left(\frac{c}{v_{ul}}\right)}^3\left(n_l\frac{g_u}{g_l}-n_u\right)ds.\hfill \end{array}$$(22)

The results are $$\begin{array}{l}{\overline{J}}_{ul}=\frac{{\text{Ω}}_{\star }}{4\pi }\left(I_v^{\star }{\alpha }_0\left({\tau }_{\text{L}}^{\text{rad}},{\tau }_{\text{C}}^{\text{rad}}\right)+S_{\text{C}}{\alpha }_1\left({\tau }_{\text{L}}^{\text{rad}},{\tau }_{\text{C}}^{\text{rad}}\right)+S_{\text{L}}{\alpha }_2\left({\tau }_{\text{L}}^{\text{rad}},{\tau }_{\text{C}}^{\text{rad}}\right)\right)\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }\left(I_v^{\text{ISM}}{\beta }_0\left({\tau }_{\text{L}}^{\uparrow },{\tau }_{\text{C}}^{\uparrow }\right)+S_{\text{C}}{\beta }_1\left({\tau }_{\text{L}}^{\uparrow },{\tau }_{\text{C}}^{\uparrow }\right)+S_{\text{L}}{\beta }_2\left({\tau }_{\text{L}}^{\uparrow },{\tau }_{\text{C}}^{\uparrow }\right)\right)\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }\left(I_v^{\text{ISM}}{\beta }_0\left({\tau }_{\text{L}}^{\downarrow },{\tau }_{\text{C}}^{\downarrow }\right)+S_{\text{C}}{\beta }_1\left({\tau }_{\text{L}}^{\downarrow },{\tau }_{\text{C}}^{\downarrow }\right)+S_{\text{L}}{\beta }_2\left({\tau }_{\text{L}}^{\downarrow },{\tau }_{\text{C}}^{\downarrow }\right)\right)\hfill \end{array}$$(23)

From this geometric model, we also get a prediction of the continuum; that is, the mean intensity at line centre frequency, as would be present if the line opacity was zero: $$\begin{array}{l}{J}_v^{\text{cont }}=\frac{{\text{Ω}}_{\star }}{4\pi }\left(I_v^{\star }{\alpha }_0\left(0,{\tau }_{\text{C}}^{\text{rad}}\right)+S_{\text{C}}{\alpha }_1\left(0,{\tau }_{\text{C}}^{\text{rad}}\right)\right)\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }\left(I_v^{\text{ISM}}{\beta }_0\left(0,{\tau }_{\text{C}}^{\uparrow }\right)+S_{\text{C}}{\beta }_1\left(0,{\tau }_{\text{C}}^{\uparrow }\right)\right)\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }\left(I_v^{\text{ISM}}{\beta }_0\left(0,{\tau }_{\text{C}}^{\downarrow }\right)+S_{\text{C}}{\beta }_1\left(0,{\tau }_{\text{C}}^{\downarrow }\right)\right).\hfill \end{array}$$(24)

Equation (24) is used to determine S_C. This means that we can calibrate S_C in such a way that our simple three-way RT model with constant continuum quantities results in the correct $J_v^{\text{cont}}$ as known from the proper solution of the continuum radiative transfer problem. This confers a significant advantage as it allows us to eliminate most of the principle problems arising from the simplified three-way RT model in the continuum.

Comparing Eq. (23) to Eq. (2), we find the definitions of our new pumping and escape probabilities: $$\begin{array}{l}{P}_{ul}^{\text{pump }}=\frac{1}{J_v^{\text{cont }}}(\frac{{\text{Ω}}_{\star }}{4\pi }\left[I_v^{\star }{\alpha }_0\left({\tau }_{\text{L}}^{\text{rad}},{\tau }_{\text{C}}^{\text{rad}}\right)+S_{\text{C}}{\alpha }_1\left({\tau }_{\text{L}}^{\text{rad}},{\tau }_{\text{C}}^{\text{rad}}\right)\right]\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }[I_v^{\text{ISM}}\left({\beta }_0\left({\tau }_{\text{L}}^{\uparrow },{\tau }_{\text{C}}^{\uparrow }\right)+{\beta }_0\left({\tau }_{\text{L}}^{\downarrow },{\tau }_{\text{C}}^{\downarrow }\right)\right)\hfill \\ +S_{\text{C}}\left({\beta }_1\left({\tau }_{\text{L}}^{\uparrow },{\tau }_{\text{C}}^{\uparrow }\right)+{\beta }_1\left({\tau }_{\text{L}}^{\downarrow },{\tau }_{\text{C}}^{\downarrow }\right)\right)])\hfill \end{array},$$(25)

$$P_{ul}^{\text{esc}}=1-\frac{{\text{Ω}}_{\star }}{4\pi }{\alpha }_2\left({\tau }_{\text{L}}^{\text{rad}},{\tau }_{\text{C}}^{\text{rad}}\right)-\frac{{\text{Ω}}_{\text{d}}}{8\pi }\left({\beta }_2\left({\tau }_{\text{L}}^{\uparrow },{\tau }_{\text{C}}^{\uparrow }\right)+{\beta }_2\left({\tau }_{\text{L}}^{\downarrow },{\tau }_{\text{C}}^{\downarrow }\right)\right).$$(26)

The interpretation of these results is as follows. $P_{ul}^{\text{pump }}$ is the probability that continuum photons make it to the considered point, given the presence of line opacity. The average is taken over all continuum photons that contribute to $J_v^{\text{cont}}$, whether they are emitted by the star or by any point in the disc, averaged over the local line profile function. $P_{ul}^{\text{esc }}$ is the probability that line photons emitted anywhere in the disc – in the direction of the considered point – do not make it to the considered point, averaged over the local line profile function. $P_{ul}^{\text{esc}}<1$ means that these line photons are either reabsorbed by line opacity, or are reabsorbed or scattered by continuum opacity along the path to the considered point.

As the order of photon emission and reabsorption and the optical depth effects along the photon paths can be reversed, $P_{ul}^{\text{esc}}$ can also be interpreted as the probability that a line photon emitted from the considered point will not have another line interaction in the disc. It is important to realise the emphasis on further line interaction here. $P_{ul}^{\text{esc}}$ is not the probability that line photons will escape the disc. This latter interpretation is only valid when there are no continuum opacities. Continuum opacity will prevent further line interaction, and so it increases $P_{ul}^{\text{esc}}$. The latter could be interpreted as ‘intrinsic escape’, meaning that, once a line photon is absorbed and re-emitted by the dust, that continuum photon will escape the disc somehow, but this is not what the equations tell us. In particular, in the optically thick limiting case with τ_L ≪ τ_C, the result is $P_{ul}^{\text{esc}}\to 1$ and $P_{ul}^{\text{pump }}\to 1$.

The six basic 2D functions α₀(τ_L, τ_C) …β₂(τ_L, τ_C) are pre-calculated in the form of 2D numerical tables. $I_{\nu }^{\star }$ and $I_v^{\text{ISM}}$ are known functions, and ${\tau }_{\text{C}}^{\text{rad}}$, ${\tau }_{\text{C}}^{\uparrow }$, and ${\tau }_{\text{C}}^{\downarrow }$ are available after initialisation of the disc structure and dust opacities in PRODIMO. As $J_v^{\text{cont}}$ is available from the proper solution of the continuum transfer problem, we can use Eq. (24) to determine and precalculate S_C for every line at any point in the disc. The centre line optical depths ${\tau }_{\text{L}}^{\text{rad}}$ and ${\tau }_{\text{L}}^{\uparrow }$ are integrated during the downward and outward sweep of the chemistry and the heating and cooling balance in PRODIMO, but we also need ${\tau }_{\text{L}}^{\downarrow }$ here, which requires another assumption.

Assumption 3 (downward line optical depths). We calculate ${\tau }_{\text{L}}^{\downarrow }$ by assuming that the local line/continuum opacity ratio remains the same on the way down and across the disc: $${\tau }_{\text{L}}^{\downarrow }={\tau }_{\text{C}}^{\downarrow }\frac{{\kappa }_{\text{L}}}{{\kappa }_{\text{C}}},$$(27)

where κ_L and κ_C are the local line and continuum opacities. The old escape probability method in PRODIMO was based on $$P_{ul}^{\text{pump,old }}={{\displaystyle \int }}^{\text{}}\varphi (x){\text{e}}^{-{\tau }_{\text{L}}^{\text{rad}}\varphi (x)}dx,$$(28)

$$P_{ul}^{\text{esc,old }}=\frac{1}{2}{{\displaystyle \int }}^{\text{}}\varphi (x)E_2\left({\tau }_{\text{L}}^{\uparrow }\varphi (x)\right)dx.$$(29)

As can be easily shown, $P_{ul}^{\text{pump }}=P_{ul}^{\text{pump,old }}$ is only valid when the direct stellar illumination dominates; that is, when $J_v^{\text{cont }}=\frac{{\text{Ω}}_{\star }}{4\pi }{I}_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}}$. In all other cases, the new $P_{ul}^{\text{pump }}$ probabilities are larger. $P_{ul}^{\text{esc}}=P_{ul}^{\text{esc},\text{ old }}$ is only valid when ${\text{Ω}}_{\star }\to 0,{\tau }_{\text{C}}^{\uparrow }=0$ and ${\tau }_{\text{C}}^{\downarrow }\to \infty$. Under any other, more realistic circumstances, the new $P_{ul}^{\text{esc }}$ probabilities are also generally larger. Both effects are caused by the inclusion of line–continuum interactions, namely (i) locally produced continuum photons can be absorbed in the line, and (ii) line photons can be absorbed locally in the continuum. We use both limiting cases discussed above to check the numerical implementation of the new escape probability method in PRODIMO.

Figure 1 shows a comparison between old and new escape and pumping probabilities based on our T Tauri standard disc model. This model was introduced by Woitke et al. (2016) and has stellar, dust, and disc parameters as listed in Table A.2. The results of the new treatment show larger $P_{ul}^{\text{pump }}$ and $P_{ul}^{\text{esc}}$ in general. The deviations are particularly obvious in the disc midplane close to the star, where trie continuum is optically thick. We see large and about equal $P_{ul}^{\text{pump }}$ and $P_{ul}^{\text{esc}}$ here, in contrast to the old treatment. More significant for line observations are the changes in the more transparent upper disc regions, where $P_{ul}^{\text{esc }}\sim P_{ul}^{\text{esc,old }}$ and $P_{ul}^{\text{pump }}>P_{ul}^{\text{pump,old }}$ in the inner disc. However, in the outer disc, this relation is reversed and we find $P_{ul}^{\text{pump }}<P_{ul}^{\text{pump,old }}$. The reasons for these changes are more complex, as $J_v^{\text{cont }}$ is not given by the attenuated stellar irradiation but is created mostly by continuum emission from the disc below.

Figure 2 shows the resulting leading heating and cooling processes in the innermost 10 AU of the disc (T Tauri standard model). Figure A.1 shows additional vertical cuts at selected radii with more details. There are two significant changes:

• (i)

  In the disc surface regions inside of 1 au, which are borderline optically thin and relevant for the IR line emission, we find increased H₂O heating and cooling rates, whereas previously there was a balance between UV-driven heating by H₂ formation and H₂O cooling. This is due to the increased pumping and escape probabilities in these regions.

• (ii)

  Thermal accommodation in the midplane is now quite unimportant for both heating and cooling, unless the midplane becomes so cold that it barely contains any molecules other than H₂ due to freeze-out. In the optically thick midplane, dust and gas now exchange energy very efficiently via line photons, more efficiently than by inelastic collisions. And since the formulation of the line heating and cooling in PRODIMO is robust, including all reverse processes and stimulated emission, $J_v^{\text{cont }}=B_v\left(T_{\text{d}}\right)$ implies that a balance between line heating and cooling can only be achieved when T_g = T_d.

Fig. 1 Escape probability ($P_{ul}^{\text{esc}}$, top) and continuum pumping probability ($P_{ul}^{\text{pump }}$, bottom) of the [OI] 63 µm line in the PRODIMO standard T Tauri model. Left and right plots show the results from the previous and the new escape probability theory. The white dashed contours show the particle density of atomic oxygen log10 nOI[ cm−3] = {8, 6, 4, 2, 0, −2}.

Fig. 2 Change in heating and cooling balance with the new escape probability treatment and the new 2D concept to apply UV molecular shielding factors. The plots show the leading heating process (top) and the leading cooling process (bottom) in a zoomed-in disc region responsible for the disc IR line emissions. Processes are only plotted when they fill at least 1% of the area in these plots. The dashed lines mark the visual disc surface Av,rad = 1 and the vertical optical extinction AV = 10. The model includes 103 heating and 95 cooling processes with 64910 spectral lines altogether.

2.2 New escape probabilities for line-flux estimations

In PRODIMO, the observable line fluxes are estimated based on the escape probabilities, summing up the contributions from all spatial cells. The following expression is used to calculate the luminosity [erg s⁻¹] of photon energy escaping in the form of line photons from a single cell in the old model: $$L^{\text{cell }}=\text{Δ}V{n}_{\text{u}}{A}_{\text{ul}}h{v}_{\text{ul}}{{\displaystyle \int }}^{\text{}}\varphi (x){\text{e}}^{-{\tau }_{\text{L}}\varphi (x)-{\tau }_{\text{C}}}dx,$$(30)

where ΔV [cm⁻³] is the volume of a cell k centred around point grid point (ix, iz). The line flux [erg s⁻¹ cm⁻²] is then calculated as $$F_{\text{line }}=\frac{{\displaystyle {\sum }_k{L}_k^{\text{cell }}}}{4\pi d^2},$$(31)

where d is the distance to the object. We consider the vertical escape here (face-on disc), and so the line and continuum optical depths are ${\tau }_{\text{L}}={\tau }_{\text{L}}^{\uparrow }$ and ${\tau }_{\text{C}}={\tau }_{\text{C}}^{\uparrow }$. Using the definition of the line emission coefficient, we have $n_{\text{u}}{A}_{\text{ul}}h{v}_{\text{ul}}=4\pi {\kappa }_{\text{L}}{S}_{\text{L}}$, and we can express the line luminosity of a cell, according to the old implementation, by $$L^{\text{cell }}=4\pi \text{Δ}V{\kappa }_{\text{L}}{S}_{\text{L}}{\text{e}}^{-{\tau }_{\text{C}}}{{\displaystyle \int }}^{\text{}}\varphi (x){\text{e}}^{-{\tau }_{\text{L}}\varphi (x)}dx.$$(32)

This result corresponds to the case of a cold background S_C = 0, where the result is always positive; see area A in Fig. 3. However, line and continuum interfere in more complicated ways than that; there is not a simple superposition, as sketch (d) in Fig. 3 suggests, because the escaping continuum photons are partly blocked by line opacity. This effect is visualised by the lower graph in case (b).

In a more detailed model, we need to compute the upper graph in case (b) and the pure continuum, and then subtract the two to simulate the way that line fluxes are measured from observations. Without continuum, we have area A in Fig. 3:

$$L^{\text{cell }}=4\pi \text{Δ}V{{\displaystyle \int }}^{\text{}}{\kappa }_{\text{L}}{S}_{\text{L}}{\varphi }_v{\text{e}}^{-{\tau }_{\text{L}}\varphi (x)-{\tau }_{\text{C}}}dv.$$(33)

With continuum, the (infinite) result is

$$L^{\text{cell }}=4\pi \text{Δ}V{{\displaystyle \int }}^{\text{}}\left({\kappa }_{\text{L}}{S}_{\text{L}}{\varphi }_v+{\kappa }_{\text{C}}{S}_{\text{C}}\right){\text{e}}^{-{\tau }_{\text{L}}\varphi (x)-{\tau }_{\text{C}}}dv.$$(34)

Considering only the continuum, the (infinite) result is

$$L^{\text{cell }}=4\pi \text{Δ}V{{\displaystyle \int }}^{\text{}}{\kappa }_{\text{C}}{S}_{\text{C}}{\text{e}}^{-{\tau }_{\text{C}}}dv.$$(35)

The continuum-subtracted line flux (area O in Fig. 3) is

$$\begin{array}{l}{L}^{\text{cell }}=4\pi \text{Δ}V{\text{e}}^{-\pi }\hfill \\ \times \left({\kappa }_{\text{L}}{S}_{\text{L}}\underset{C_1\left({\tau }_{\text{L}}\right)}{\underbrace{{\displaystyle \int \varphi (x){\text{e}}^{-{\tau }_{\text{L}}\varphi (x)}dx}}}-{\kappa }_{\text{C}}{S}_{\text{C}}\text{Δ}{\nu }_{\text{D}}\underset{C_2\left({\tau }_{\text{L}}\right)}{\underbrace{{\displaystyle \int \left(1-{\text{e}}^{-{\tau }_{\text{L}}\varphi (x)}\right)dx}}}\right)\hfill \end{array}.$$(36)

Equation (36) shows that there is an additional correction term with respect to the old way of calculating the escaping line luminosity from a cell (Eq. (32)), which is always negative. This new term describes the loss of continuum flux due to line absorption. The integration boundaries in these equations are −∞and +∞, and so both intermediate results (Eqs. (34) and (35)) are infinite, but the subtraction of the two (Eq. (36)) is finite. The function C₁(τ_L) describes the probability that an outgoing line photon is not reabsorbed in the line again, and C₂(τ_L) describes the reduction of continuum photons escaping the object due to line opacity. The two functions C₁(τ_L) and C₂(τ_L) are plotted in Fig. 4.

According to this new formulation of escape probability, the line flux is always smaller when continuum emission effects are included, and can become negative when the continuum source function S_C is as large as the line source function S_L. However, we note that the continuum and the line are usually emitted from different disc regions.

Figure 5 shows a comparison between the line-flux predictions obtained with the old and the new escape probability treatment. All results are shown with respect to the proper full PRODIMO line RT results for face-on inclination. We generally see a good agreement between the results of all three methods, with differences being of order 10–30%. The differences become more noticeable at shorter wavelengths, where the lines are mostly emitted from the inner disc. The new escape probability treatment is superior and more balanced around the y = 1 line, whereas the old escape probability treatment has the tendency to overpredict the line fluxes. Two NH₃ lines at around 10.5 µm turn out to be absorption lines according to the proper line RT, and are not plotted in Fig. 5. Both escape probability methods fail to predict these lines correctly. The old escape probability method predicts (by design) emission lines. The new escape probability method does produce absorption lines as it should, but the negative fluxes are too large. We note that both escape probability methods can predict the submillimetre (submm) lines very well (e.g. low-J CO lines, CN, HCN, HCO⁺ in Fig. 5). The fluxes of the far-infrared (FIR) lines (e.g. Herschel/HIFI water lines, fine-structure cooling lines of oxygen) are predicted better than 10%.

Fig. 3 Effect of continuum emission on line flux. Case (a) shows an emission line on a cold continuum SC = 0. Case (b) is an emission line on a warm continuum SL > SC, and case (c) is an absorption line where the continuum is hot SL < SC. The area A is the flux of photon energy [erg s−1 cm−2] carried away by line photons escaping from the object (after continuum extinction and line self-absorption) and is the same in cases (a), (b) and (c). However, the objective is to derive the area O, which is the line flux as measured by observers. For a warm continuum, clearly O < A, and for an absorption line O < 0 (although A > 0).

Fig. 4 Functions C1(τL) and C2(τL) as introduced in Eq. (36).

2.3 Photo cross sections

PRODIMO now uses the new high-resolution photoionisation and photodissociation cross sections ${\sigma }_{\text{ph}}^{\text{mol}}$ of Heays et al. (2017) obtained from the Leiden database¹.

Fig. 5 Line fluxes predicted by the old escape probability treatment (blue) and new escape probability treatment (orange). All line fluxes are divided by the results from the full line RT for zero inclination (face-on).

2.4 Ultraviolet molecular shielding and self-shielding

Shielding factors are used in PRODIMO to reduce the molecular photoionisation and photodissociation rates:

$$R_{\text{ph}}=4\pi {{\displaystyle \int }}^{\text{}}{\sigma }_{\text{ph}}(v)\frac{J_v}{hv}dv,$$(37)

where σ_ph(v) is the photodissociation cross section and J_v is the local mean intensity. Without molecular shielding, we have $J_v=J_v^{\text{cont }}$ where $J_v^{\text{cont }}$ is the result from the continuum radiative transfer (RT). However, the UV molecular opacities are not included in the continuum RT; they generally cause $J_v<J_v^{\text{cont}}$, because there is no molecular UV emission in the disc. These opacities therefore reduce the photo rates.

There are various molecular shielding effects included, both self-shielding i→i and shielding of a molecule j by another molecule i; see for example Bethell & Bergin (2009). Water shielding was recently discussed by Bosman et al. (2022a,b). In some cases, the shielding factors s_i→j are taken from published 1D semi-infinite slab-models, which resolve the UV lines and have published efficiency factors as a function of the column density of the molecular absorber ‘towards the source’ and (local) temperature:

$$s_{i\to j}=s_{i\to j}\left(N_i,T_{\text{g}}\right)\text{\hspace{1em} with \hspace{1em}}{N}_i={{\displaystyle \int }}^{\text{}}n(i)ds.$$(38)

In general, the shielding factors can be applied as a product, because $\exp \left(-{\tau }_{\text{all}\to j}\right)=\exp \left(-\underset{i}{{{\displaystyle \sum }}^{\text{}}}{\sigma }_{\text{i}}(v)N_i\right)={\prod }_i\exp \left(-{\sigma }_{\text{i}}(v)N_i\right)$, at least approximately. Table 1 summarises the basic self-shielding functions used. In PRODIMO, all photorates are calculated by numerical integration of Eq. (37) over the frequency-dependent photo cross sections and the local UV radiation field J_v The latter is interpolated between the values known on approximately seven UV spectral gridpoints used in the continuum radiative transfer; see Woitke et al. (2009) for details. In cases where we do not have detailed photo cross sections, but only the a, β, and γ coefficients from a chemical database, we fit a generic cross-section ${\sigma }_{\text{ph}}(\lambda )={\sigma }_0\exp \left(-{\left({\lambda }_0-\lambda \right)}^2/\text{Δ}{\lambda }^2\right)$, with two free parameters σ₀ and λ₀, until the photorates are correct for optical extinctions A_V = 0 and A_V = 0.5 using standard interstellar dust extinction properties.

Therefore, it is relatively straightforward to use any frequency-dependent recipes to account for the effect of molecular shielding. In particular, and this is new in PRODIMO, we use the following idea to account for generic self-shielding:

$$s_{i\to i}={{\displaystyle \int }}^{\text{}}{\sigma }_{\text{ph}}(v)\frac{J_v}{hv}{\text{e}}^{-N_i{\sigma }_{\text{ph}}(v)}dv/{{\displaystyle \int }}^{\text{}}{\sigma }_{\text{ph}}(v)\frac{J_v}{hv}dv.$$(39)

Other recipes are applied to the shielding of any molecule by H₂ and C; see Table 1. For example, Eq. (8) in Kamp & Bertoldi (2000) describes the T-dependent reduction of UV light by H₂ lines when detected with the flat opacity of neutral C, which is assumed to provide a first approximation of the shielding effect on any molecule with uncorrected σ_ph(λ) for wavelengths shorter than 1100 Å.

2.5 New 2D treatment of molecular shielding

The question arises as to how to deal with the problem of molecular shielding in 2D, and specifically as to which column density N_i is to be considered, choosing from the radial one, the vertically upward one, a mixture, or even something else. Equation (24) holds an answer to this question. Using the explicit shapes of α₀, β₀ an β₁ for τ_L → 0, and neglecting the α₁ term, which is the continuum emission of the disc in the solid angle of the star, and is usually very small, this equation reads

$$\begin{array}{l}{J}_v^{\text{cont }}=\frac{{\text{Ω}}_{\star }}{4\pi }{I}_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad }}}+\frac{{\text{Ω}}_{\text{d}}}{8\pi }{I}_v^{\text{ISM}}\left(E_2\left({\tau }_{\text{C}}^{\uparrow }\right)+E_2\left({\tau }_{\text{C}}^{\downarrow }\right)\right)\hfill \\ +\frac{{\text{Ω}}_{\text{d}}}{8\pi }{S}_{\text{C}}\left(2-E_2\left({\tau }_{\text{C}}^{\uparrow }\right)-E_2\left({\tau }_{\text{C}}^{\downarrow }\right)\right).\hfill \end{array}$$(40)

Equation (40) states the influence of the three principle sources (the star, the ISM background, and the disc itself) on the calculated $J_v^{\text{cont}}$. Contrary to the IR pumping of molecules, which often comes from the warm extended dust below the point of interest, the UV radiation mainly arrives from two directions, the direct radial and the vertically downward directions. Clearly, if direct illumination from the star dominates, then $J_v^{\text{cont }}\approx \frac{{\text{Ω}}_{\star }}{4\pi }{I}_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}}$ and one should use the radial column densities towards the star to compute the shielding factors. If interstellar UV irradiation dominates, one should take the vertical column densities. However, in the IR-line-emitting disc regions, neither of these two terms is usually the important one, and it is the downward scattering of star light by the small dust in the disc’s upper layers that dominates, which is represented by S_C and expressed by the third term in Eq. (40).

It is not obvious which molecular column densities should be considered for the molecular shielding factors applied to the third term in Eq. (40). The proper approach would be to consider the molecular column density along the average scattered photon path, but this information is not available. One could argue in favour of taking the radial column density at a much higher disc level up to the scattering point and then adding the downward vertical column density from that point when considering single vertical scattering events, but then the point at which that scattering event takes place needs to be determined.

Instead we follow a much more simple approach here. We assume that the second part of the molecular column density along the scattered photon path is dominant, namely the downward vertical column density from the scattering site. Since the vertical column density increases exponentially on the way down, because (a) the density increases and (b) the molecular concentration increases, it is mostly determined locally, and there is no significant difference when simply taking the full vertical column density.

According to this approach, we use the radial column densities for the first, and the vertical column densities for the second and third terms in Eq. (40), which can be reformulated as

$$J_{v,\text{rad}}^{\text{cont }}=\frac{{\text{Ω}}_{\star }}{4\pi }{I}_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rad}}}$$(41) $$J_{v,\text{ver}}^{\text{cont }}=J_v^{\text{cont }}-J_{v,\text{rad}}^{\text{cont }}$$(42) $$J_v=s_{i\to j}\left(N_i^{\text{rad}},T_{\text{g}}\right)J_{v,\text{ rad }}^{\text{cont }}+s_{i\to j}\left(N_i^{\text{ver}},T_{\text{g}}\right)J_{v,\text{ vert }}^{\text{cont }}.$$(43)

In previous PRODIMO models, we used the stellar and interstellar irradiation terms (parts 1 and 2 in Eq. (40)) to calculate weighting factors for the application of shielding factors according to the following scheme:

$$J_{v,\text{ rad }}^{\text{cont }}=\frac{{\text{Ω}}_{\star }}{4\pi }{I}_v^{\star }{\text{e}}^{-{\tau }_{\text{C}}^{\text{rd}}}$$(44) $$J_{v,\text{ver}}^{\text{cont }}=I_v^{\text{ISM}}{E}_2\left({\tau }_{\text{C}}^{\uparrow }\right)$$(45) $$w_{\text{rad}}=J_{v,\text{rad}}^{\text{cont }}/\left(J_{v,\text{rad}}^{\text{cont }}+J_{v,\text{ver}}^{\text{cont }}\right)$$(46) $$w_{\text{vert }}=J_{v,\text{ ver }}^{\text{cont }}/\left(J_{v,\text{ rad }}^{\text{cont }}+J_{v,\text{ ver }}^{\text{cont }}\right)$$(47) $$s_{i\to j}=w_{\text{rad }}{s}_{i\to j}\left(N_j^{\text{rad}},T_{\text{g}}\right)+w_{\text{vert }}{s}_{i\to j}\left(N_j^{\uparrow },T_{\text{g}}\right)$$(48) $$J_v=s_{i\to j}{J}_v^{\text{cont }},$$(49)

where w_rad and w_vert are the radial and vertical weights. However, this approach involves a clearly far more questionable assumption, because the often most important term (the third term in Eq. (40)) is not even considered, and thus $J_v^{\text{cont }}\gg J_v^{\text{rad}}+J_v^{\text{vert}}$. Instead, based on the mixing ratio between $J_v^{\text{rad}}$ and $J_v^{\text{vert }}$, an effective shielding factor was calculated and applied to all of $J_v^{\text{cont }}$. As Eq. (48) ignores the scattering altogether, the old vertical part $J_{v,\text{ ver }}^{\text{cont }}$ is too small, and therefore we switched from radial to vertical column densities too late. As the radial column densities are generally larger than the vertical ones, the effects of molecular shielding were overestimated in the old models.

Figure 6 shows that, as expected, the formation of CO now happens slightly deeper in the disc. In both models, CO is found to be abundant when (i) the ionisation parameter χ/n is low enough (χ is the strength of the UV field; see Eq. (42) in Woitke et al. (2009), and n is the total hydrogen nuclei density [cm⁻³]) and (ii) the dust temperature is high enough T_d ≳ 20 K. These two criteria are related to CO photodissociation and freeze-out. The critical value for the CO photodissociation log₁₀(χ/n) changes from about -4.5 to -5 in the new models.

The feedback on the more complex molecules, such as HCN and C₂H₂, is remarkable, leading to profoundly different chemical conditions in the IR-line-emitting disc regions. In Fig. 6, we concentrate on the upper, warm line-emitting regions. There are also some regions with high concentrations of HCN and C₂H₂ in the midplane, but these molecules are covered by optically thick dust and play no significant role in the IR-line emission in this T Tauri standard disc setup. The formation of these IR-active molecules in the upper disc happens in slightly deeper layers according to the new models, and as the gas densities are larger there, the formation of more complex molecules is generally favoured in the new models.

In particular, the formation of C₂H₂ needs both H₂ and neutral carbon atoms; see e.g. Agúndez et al. (2008). In the lower part of Fig. 6, we overplot H/H₂ = 10 and CO/C = 100 contour lines. These two lines nicely encircle the region of high C₂H₂ concentrations in the model. Similar results were reported by Walsh et al. (2015) and Anderson et al. (2021). In the inner disc, these lines are crossed in our models, and so when H₂ finally forms, there is no longer any neutral C, and we get no C₂H₂; see also Kanwar et al. (2024). In the new model, these two lines are further apart, there is a wider region to form C₂H₂, and that region now also extends to smaller radii.

A similar effect is found for HCN, where the changes are even more pronounced. The formation of HCN requires H₂ and neutral N atoms. In the new models, both transitions H → H₂ and N → N₂ happen deeper in the disc, and they are vertically more separated, leading to increased HCN column densities in the upper layers of the inner disc.

2.6 Heating by photodissociation

One heating process that has so far been partly neglected in PRODIMO is the ‘direct’ heating by photodissociation:

$${\text{Γ}}_{\text{pd}}=\underset{\text{mol}}{{{\displaystyle \sum }}^{\text{}}}{n}_{\text{mol}}{R}_{\text{mol}}^{\text{ph}}\text{Δ}{E}_{\text{mol}},$$(50)

where n_mol [cm⁻³] is the molecular particle density, $R_{\text{mol}}^{\text{ph}}$ [s⁻¹] the photodissociation rate, and ΔE_mol is the mean energy liberated per photodissociation. Similar to photoionisation, the excess photon energy can create superthermal dissociation products, or can lead to vibrationally or electronically excited molecular fragments, and the key question pertains to how much of this energy is thermalised by the gas.

Gu et al. (2020) treated this process with 100% efficiency, assuming that the energy released as local heat is set as the incident photon energy in excess of the respective dissociation threshold: $$\text{Δ}{E}_{\text{mol}}=\frac{{{{\displaystyle \int }}^{\text{}}}_{v_0}^{\infty }\frac{J_v}{hv}{\sigma }_{\text{ph}}^{\text{mol}}(v)h\left(v-v_0\right)dv}{{{{\displaystyle \int }}^{\text{}}}_{v_0}^{\infty }\frac{J_v}{hv}{\sigma }_{\text{ph}}^{\text{mol}}(v)dv},$$(51)

where v₀ is the photodissociation threshold frequency. However, more detailed studies show that not all of this energy is thermalised; for example, excited dissociation products can decay radiatively, or can trigger follow-up chemical reactions.

Glassgold & Najita (2015) published detailed energies per photodissociation ΔE for H₂ (1.2 eV), CO (1.4 eV), H₂O (1.2 eV), and OH (3.6 eV), and we have now implemented these data into PRODIMO. For all other molecules, we conservatively assume ΔE ≈ 2 eV. Previously in PRODIMO, only H₂, with ΔE(H₂) = 0.4 eV (Black & van Dishoeck 1987), was considered for the heating by photodissociation in addition to the heating by photoionisation of H⁻, C, Fe, Si, and Mg. Figure 2 shows that this additional heating process can be quite relevant in the IR-active molecular layer. Interestingly, the reduced UV molecular shielding factors cause larger photorates, which lead to a more active chemistry, where the molecular fragments created by photodissociation react further exothermally, which increases the area where chemical heating is important.

Fig. 6 Change of molecular concentrations caused by introducing our new 2D concept to apply molecular shielding factors.

2.7 Chemical rate network

We used the large DIANA standard chemical network (Kamp et al. 2017) in this work, which has 235 chemical species, with reaction rates mostly taken from the UMIST 2012 database (McElroy et al. 2013), but including a simple freeze-out and desorption ice chemistry (Woitke et al. 2009), X-ray processes partly based on doubly ionised species (Aresu et al. 2011), excited molecular hydrogen (see Kamp et al. 2017), and polycyclic aromatic hydrocarbons in five different charging states (Thi et al. 2019), which is 3066 reactions altogether. H₂ formation of grain surfaces is calculated according to Cazaux & Tielens (2002, 2004) with updates described in Cazaux & Spaans (2009); Cazaux & Tielens (2010). We use our standard element abundances (Table A.1), which are close to solar (C/O = 0.46) except for the depleted metals Na, Mg, Si, S, and Fe.

PAH molecules are included in the radiative transfer, with opacities calculated from their overall abundance and an assumed ratio between charged and neutral PAHs to avoid global iterations; see details in Woitke et al. (2016, 2019). PAHs are also included in the chemical network with five charging states (Thi et al. 2019), and are included in the calculation of heating and cooling rates, as the photoelectric effect on PAHs is one of the most important heating processes in the upper disc layers. The PAH heating rate is calculated from the calculated abundances of the different charging states, their individual photo-cross sections, and the calculated local radiation field.

In recent updates of our chemical network, we added a few hand-picked gas-phase reactions not included in UMIST 2012. One of these reactions turned out to be critical for the fitting of the JWST spectrum of EX Lupi (Sect. 3): $$\text{C}+{\text{H}}_2\text{O\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\to \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}HCO}+\text{H,}$$(52)

from Hickson et al. (2016), who measured surprisingly large rates for this reaction at low temperatures, enhanced by quantum mechanical tunnelling. Although the Kinetic Database For Astrochemistry (KIDA²) does not include this reaction in their standard releases, they list it with Arrhenius parameters {α = 1.07 × 10⁻⁸cm³s⁻¹, β = −1.59 , γ = 0} with reference to Hickson et al. (2016). In our models, whenever a CO-bond is broken, for example by photodissociation or X-ray-induced reactions, this reaction would be very efficient in immediately reforming that bond, which would prohibit the subsequent formation of hydrocarbon molecules and HCN as discussed in Sect. 4.3. Tests show that if we were to include this reaction with the Arrhenius coefficients from the KIDA database, the C₂H₂ column density in the line emitting regions of our EX Lupi model would decrease by a factor of about 300, and peak line flux at 13.7 µm would decrease by a factor of about 5. However, closer inspection of Fig. 3 in Hickson et al. (2016) shows that the Arrhenius coefficients {α = 4.2 × 10⁻¹³cm³s⁻¹ , β = −2.5 , γ = 0} fit the original data points at 50 K, 75 K, and 100 K much better, and are consistent with the measurements of Husain & Kirsch (1971), which show that the rate coefficient of this reaction is very small at room temperature (< 10⁻¹² cm³s⁻¹). We therefore decided to use our own Arrhenius fit of this reaction in this paper.

2.8 New HITRAN line-selection rules

The expansion from the Spitzer/IRS to the JWST/MIRI wavelength range made it necessary to revise our line-selection strategy, by which we import only a limited number of molecular lines from the HITRAN database. We also updated PRODIMO to use the HITRAN 2020 catalogue (Gordon et al. 2022) instead of the HITRAN 2012 catalogue. The imported data are the upper level energies E_u, line wavelengths λ_ul, Einstein coefficients A_ul, and level degeneracies g_u. The new HITRAN 2020 database has a greater number of lines, which, if we would import all lines, would vastly increase PRODIMO’s memory consumption. As explained in Woitke et al. (2018a), we therefore have to limit our line selection to lines with a minimum line strength, defined as $$\text{ line strength }=A_{ul}{g}_u\exp \left(-\frac{E_u}{k(1500\text{K})}\right)\text{. }$$(53)

In combination with further restrictions on wavelength range and upper level energy E_u, see Table 2, we aim to select only the relevant, strong observable lines. Together with the rovibrational CO lines of the first four vibrationally excited states (Thi et al. 2013), ro-vibrational H₂O lines, rotational lines of various molecules, and lines of atoms and ions (Woitke et al. 2018a), we are now tracing 64890 lines altogether in this model. All lines are automatically considered for heating and cooling.

2.9 Dust settling

Considering a turbulent one-dimensional column of gas and dust, Mulders & Dominik (2012) write down the continuity equation for passive dust particles with a size distribution function f (a, z, t) as $$\frac{\partial f\left(a,z\right)}{\partial }=\frac{\partial }{\partial z}\underset{j\text{diff}}{\underbrace{\left(\rho D_{\text{d}}\left(a\right)\frac{\partial }{\partial z}\left(\frac{f\left(a,z\right)}{\rho }\right)\right)}}-\frac{\partial }{\partial z}\underset{j\text{set}}{\underbrace{\left(f\left(a,z\right){\stackrel{°}{\upsilon }}_{\text{dr}}\left(a,z\right)\right)}},$$(54)

where ρ is the gas mass density, t the time, z the vertical coordinate, and a the dust particle radius. Assuming that there is no bulk gas motion, the grains are transported downwards by gravitational settling with their size-dependent equilibrium settling velocity ${\stackrel{\circ }{v}}_{\text{dr}}(a,z)$, and are transported upwards by eddy diffusion, The upward mixing is driven by vertical dust particle concentration gradients $\frac{\partial }{\partial z}\left(\frac{f(a,z)}{\rho }\right)$ and calculated according to the local dust diffusion coefficient D_d. Assuming that the grains have relaxed towards a steady state, we have for each size ∂/∂t f(a, z) → 0 and zero net flux j_diff – j_set → 0, and therefore $$\rho D_{\text{d}}\frac{\partial }{\partial z}\left(\frac{f(a,z)}{\rho }\right)=f(a,z){\stackrel{°}{\upsilon }}_{\text{dr}}(a,z),$$(55)

from which we find $$\frac{\partial }{\partial z}\left(\ln \frac{f\left(a,z\right)}{\rho }\right)=\frac{{\stackrel{°}{\upsilon }}_{\text{dr}}\left(a,z\right)}{D_{\text{d}}}.$$(56)

For large Knudsen numbers and subsonic particle velocities (Epstein regime), the equilibrium settling velocity, also referred to as the terminal fall speed, is given by Schaaf (1963) $${\stackrel{°}{\upsilon }}_{\text{dr}}(a,z)=-g_z{\tau }_{\text{stop }}(a,z)=-z{\text{Ω}}^2{\tau }_{\text{stop }}(a,z),$$(57)

where g_z is the vertical gravitational acceleration, and Ω is the Keplerian circular frequency. The stopping timescale, or frictional coupling timescale, is given by $${\tau }_{\text{stop }}(a,z)=\frac{{\rho }_{\text{m}}a}{\rho (z)v_{\text{th}}},$$(58)

where ρ_m is the dust material density, and the thermal velocity is $$v_{\text{th }}=\sqrt{\frac{8k{T}_{\text{g}}}{\pi \mu }}=\sqrt{\frac{8}{\pi }}{c}_{\text{S}},$$(59)

where k the Boltzmann constant, T_g the local gas temperature, µ the mean molecular weight, c_S the local isothermal speed of sound, and $p=c_{\text{S}}^2\rho$ is the local gas pressure. Equation (58) appears to be widely incorrect in the literature, where c_S is often erroneously used instead of the correct mean of the magnitude of the thermal velocities of gas particles υ_th, as required in Schaaf’s formula; see Woitke & Helling (2003).

Eddy diffusion (‘turbulent mixing’) is assumed to be dominant over gas-kinetic diffusion. The eddy diffusion coefficient is given by a characteristic velocity multiplied by a characteristic length, $$D_{\text{gas }}=\langle v_z\rangle L,$$(60)

where L is the size of the largest eddies (‘mixing length’) and 〈υ_z〉 is the root-mean-square average of the fluctuating part of the vertical gas velocities. In a Kolmogorov-type power spectrum P(k) ∝ k^−5/3, there are different turbulent modes associated with different wave-numbers k or different spatial eddy sizes l. Any given particles of size a tend to co-move with the sufficiently large and slow turbulent eddies, whereas their inertia prevents them from following the short-term, small turbulence modes. In order to arrive at an effective particle diffusion coefficient, the advective effect of all individual turbulent eddies has to be averaged, and thereby transformed into a collective particle-diffusion coefficient. This procedure has been carried out, with different methods and approximations, by Dubrulle et al. (1995), Schräpler & Henning (2004), Youdin & Lithwick (2007), and Mulders & Dominik (2012). The result of Mulders & Dominik (2012) is $$D_{\text{d}}=\frac{D_{\text{gas }}}{1+S{t}^2},$$(61)

where St is the Stokes number given by $$St=\frac{{\tau }_{\text{stop }}}{{\tau }_{\text{eddy }}}$$(62)

where τ_eddy is the eddy turnover or turbulence correlation timescale in consideration of the largest eddy, which is assumed to be equal to the mixing length L.

$${\tau }_{\text{eddy }}=\frac{L}{\langle v_z\rangle }.$$(63)

According to an α disc model (Shakura & Sunyaev 1973), the typical vertical length scale is $L=\sqrt{\alpha }{H}_p$ and $v_z=\sqrt{\alpha }{c}_s$ (Ormel et al. 2007), and therefore

$$D_{\text{gas}}=\alpha c_{\text{S}}{H}_p,$$(64) $${\tau }_{\text{eddy }}=\frac{H_p}{c_{\text{S}}}=\frac{1}{\text{Ω}}.$$(65)

Our basic Eq. (56) for the balance between upward mixing and downward gravitational settling is therefore

$$\frac{\partial }{\partial z}\left(\ln \frac{f(a,z)}{\rho }\right)=-\frac{z\text{\hspace{0.17em}Ω}{\tau }_{\text{stop }}}{\alpha H_p^2}\left(1+S{t}^2\right),$$(66)

where all quantities on the right hand side of Eq. (66), that is, Ω, α, D_gas, T, c_S, µ, ρ, and ρ_m, can in general depend on height z. In addition, St and τ_stop depend on z and on grain size a.

Dubrulle et al. (1995) ignored all z-dependencies on the ride hand side of Eq. (66) and replaced all variables with their mid-plane values. These authors also dropped the (1 + St²) term. In that case, we can carry out the z-integration to find

$$\ln \frac{f(a,z)/f(a,0)}{\rho (z)/\rho (0)}=-\frac{\text{Ω}{\tau }_{\text{stop }}(0)}{\alpha H_p^2}{{{\displaystyle \int }}^{\text{}}}_0^z{z}^{\prime }d{z}^{\prime }=-\frac{z^2}{2H_p^2}\frac{\text{Ω}{\tau }_{\text{stop }}(0)}{\alpha }.$$

Introducing the size-dependent dust scale height H_a via $f(a,z)=f(a,0)\exp \left(-\frac{z^2}{2H_a^2}\right)$, and using $\rho (z)=\rho (0)\exp \left(-\frac{z^2}{2H_p^2}\right)$, the result is

$$\frac{H_p^2}{H_a^2}=1+\frac{\text{Ω}{\tau }_{\text{stop }}(0)}{\alpha }.$$(67)

Riols & Lesur (2018) also drop the (1 + St²) term in Eq. (66), but take into account that τ_stop depends on density, and is therefore z-dependent. Furthermore, these authors explicitly assume a Gaussian for the vertical gas density distribution. In that case, the stopping timescale can be expressed as

$${\tau }_{\text{stop }}={\tau }_{\text{stop }}(0)\mathrm{\hspace{0.17em}exp\hspace{0.17em}}\left(\frac{z^2}{2H_p^2}\right),$$(68)

and integration of Eq. (66), assuming that α, c_S, H_p, and Ω are height-independent, yields

$$\ln \frac{f(a,z)/f(a,0)}{\rho (z)/\rho (0)}=-\frac{\text{Ω}{\tau }_{\text{stop }}(0)}{\alpha H_p^2}{{{\displaystyle \int }}^{\text{}}}_0^z{z}^{\prime }\exp \left(\frac{z^{\prime 2}}{2H_p^2}\right)d{z}^{\prime },$$(69)

which reproduces Eq. (33) in Riols & Lesur (2018). The analytical solution of Eq. (69) is

$$\frac{H_p^2}{H_a^2(z)}=1+\frac{\text{Ω}{\tau }_{\text{stop }}(0)}{\alpha }\left(\frac{2H_p^2}{z^2}\right)\left(\exp \left(\frac{z^2}{2H_p^2}\right)-1\right).$$(70)

The decrease in H_a with respect to Eq. (67) is remarkable; the function (2/x²)^[ exp(x²/2) – 1] is 3.2 at x = 2, 20 at x = 3, and 370 at x = 4. Infrared emission lines from discs are quite typically emitted from layers located at about 3–4 scale-heights (Woitke et al. 2018b). The Riols & Lesur settling description reduces the local dust/gas ratios by orders of magnitude in these regions (see Fig. 7), and even reduces the abundance of the smallest grains at high altitudes. We suggest that Eq. (70) should become the standard in disc modelling; we have implemented it in PRODIMO, and used it in the present study.

However, in order to remain consistent with previous publications, in PRODIMO, we multiply the second terms on the right side of Eqs. (67) and (70) by a factor $\sqrt{1+{\gamma }_0}$, with γ₀ = 2 for compressible turbulence, and incorrectly use c_S instead of υ_th in Eq. (57). These two changes amount to a factor of $\sqrt{8\times 3/\pi }\approx 2.8$, which can be compared to the uncertainty in the settling parameter α. Equation (70) can be generalised to a numerical integration of Eq. (56) for cases where $\rho (z)=\rho (0)\exp \left(-\frac{z^2}{2H_p^2}\right)$ is not valid, and c_S = H_p Ω is not valid, for example in hydrostatic disc models where both the gas temperature and the mean molecular weight depend on z.

Fig. 7 PRODIMO model for the quiescent state of EX Lupi. The top plots show the assumed gas and dust column densities (left) and scale heights (right). The second row shows the corresponding gas densities (left) and gas-to-dust mass ratio after settling (right). The third row shows the UV field strength χ (left), and the X-ray ionisation rate ζX in comparison to the cosmic-ray ionisation rate (right). The bottom row shows the resulting dust-temperature (left) and gas-temperature (right) structures. Additional contour lines show selected values for the vertical optical extinction AV and the radial optical extinction AV,rad as indicated.

2.10 Rounded inner rims

As previous Spitzer and recent JWST observations show (Banzatti et al. 2017, 2023; Kóspál et al. 2023; Grant et al. 2023; Tabone et al. 2023), the MIR line spectroscopy of discs first and foremost probes the physics and chemistry in the innermost disc regions ≲ 1 au, close to where the dust reaches its sublimation temperature, with line-emitting radii sometimes as small as 0.1 au for T Tauri stars. Therefore, the shape of the inner rim becomes a significant factor for understanding the new JWST spectra. Our knowledge about these inner disc regions is still rather limited, because of both the difficulty in observing these regions with sufficient resolution and the complexity of the physics involved. According to radiation hydrodynamics models, the direct stellar irradiation causes the dust around the inner rim to evaporate and form a puffed-up inner rim (Natta et al. 2001; Dullemond et al. 2001; Dullemond & Monnier 2010). The radial transition in terms of gas surface density across the inner rim is not an infinitely sharp jump, but a relatively diffuse region, which means that the gas pressure stabilises the rim (Isella & Natta 2005; Kama et al. 2009; Woitke et al. 2009; Flock et al. 2017). In addition, evolutionary MHD models suggest that the switch from magnetically active inner disc to the dead zone may cause a local pressure maximum and formation of gas–dust rings (Masset et al. 2006; Dzyurkevich et al. 2010; Kadam et al. 2019, 2022), which causes a more gradual buildup of the column density towards the first pressure bump. Inspired by the inner disc structures suggested by these models, here we introduce a parametric approach, assuming

$$\begin{array}{l}\text{Σ}(r)\propto \exp \left[-{\left(\frac{R_{\text{soft}}}{r}\right)}^s\right]r^{-ϵ}\exp \left[-{\left(\frac{r}{R_{\text{tap}}}\right)}^{2-\gamma }\right],\hfill \\ R_{\text{soft}}=\text{raduc}\times R_{\text{in}},\hfill \\ s=\frac{\ln [\ln (1/\text{reduc})]}{\ln (\text{raduc})},\hfill \end{array}$$(71)

where Σ(r) is the mass column density at radius r, raduc is the radius up to which the column density initially increases, in units of the inner disc radius R_in, and reduc is the factor by which the column density is reduced in Eq. (71) at R_in. As in previous publications, R_tap is the tapering off radius, and the self-similar solution is obtained for γ = ϵ, but we rather keep ϵ and γ as independent parameters.

In hydrostatic models, where the hydrostatic equilibrium is solved in each vertical column (Natta et al. 2001; Kama et al. 2009; Dullemond & Monnier 2010), we see large scale-heights at the inner rim that is directly heated by the star, followed by a dip in Hp(r) behind the inner rim, because the inner rim casts a shadow on the subsequent disc where the temperatures are lower; see for example Fig. 9 in Woitke et al. (2009) and Fig. 6 in Woitke et al. (2016). We introduce two more parameters to describe this ubiquitous feature of hydrostatic models:

$$\begin{array}{l}{H}_p(r)=H_0{\left(\frac{r}{r_0}\right)}^{\beta }f(r),\hfill \\ f(r)=1-\text{deduc}\times \exp \left(-{\left(\frac{\text{Δ}r}{w}\right)}^2\right),\hfill \\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Δ}r=\ln \left(r/R_{\text{in}}\right)-\ln (\text{daduc}),\hfill \\ w=\{\begin{array}{ll}\frac{1}{2}\ln (\text{deduc}),\hfill & \text{ifΔ}r<0\hfill \\ \ln (\text{deduc}),\hfill & \text{otherwise’}\hfill \end{array}\hfill \end{array}$$(72)

where daduc is the dip radius in units of the inner radius, and deduc is the reduction of the scale height at the dip radius. The effect of the new disc shape parameters reduc, raduc, deduc, and daduc on column density structure and scale heights is visualised in the upper panel of Fig. 7.

3 The PRODIMO disc model for EX Lupi

3.1 Stellar setup

Our photometric data collection for the quiescent state of EX Lupi is shown in Table 3, and the assumed stellar properties are listed in Table 4. Stellar luminosity, effective temperature, and optical extinction are taken from Sipos et al. (2009). Stellar mass, distance to EXLupi, and disc inclination are taken from the collection of ALMA observations by White et al. (2020). The stellar excess UV parameters are derived from XMM UV photometric observations (see Table 3). The quantity f_UV = L_UV/L_* is the UV excess between 91.2 nm and 250 nm. p_UV is a power-law index explained in Appendix A of Woitke et al. (2016). The X-ray properties are explained in Appendix A of Woitke et al. (2016) and are taken from measurements of Teets et al. (2012), which indicate a significant increase in the X-ray luminosity during burst phase, in particular for the soft X-rays. The accretion rate Ṁ_acc = 4 × 10⁻¹⁰ M_⊙ yr⁻¹ is taken from Sicilia-Aguilar et al. (2015), which is used for viscous gas heating according to Eq. (2) in D'Alessio et al. (1998); see details in Woitke et al. (2019).

Recently, Cruz-Sáenz de Miera et al. (2023) analysed an XSHOOTER spectrum of EX Lupi taken in July 2022, just one month before the JWST observations were carried out. This 0.3–2.2 µm spectrum suggests a much higher mass-accretion rate of Ṁ_acc ≈ (2–3) × 10⁻⁸ M_⊙ yr⁻¹. Another value was reported by Wang et al. (2023), namely Ṁ_acc ≈ 3 × 10⁻⁹ M_⊙ yr⁻¹. The data analysis is intertwined with the determination of the optical extinction A_V of the star. Cruz-Sáenz de Miera et al. 2023 arrive at a large value of A_V ≈ 1.1 using R_V = 3.1. In contrast, most published works on EX Lupi have so far assumed A_V ≈ 0, including Sipos et al. (2009). Varga et al. (2018) used A_V = 0.2, Alcalá et al. (2017) used A_V = 1.1, and Wang et al. (2023) used A_V = 0.1. If A_V is indeed as large as 1.1, the XMM photometric data would imply a ten times greater UV luminosity ( f_UV ≈ 0.1). Extrapolating the slab-model for the continuum emission due to accretion of Cruz-Sáenz de Miera et al. 2023 to shorter wavelengths, Ábrahám (2023, priv. comm.) derive f_UV ≈ 0.02, but this UV excess would not contain the strong emission lines, such as Ly-α, because of stellar activity. A discussion on how our modelling results change if larger values of f_UV and Ṁ_acc are used can be found in Sect. 4.5.