© P. Woitke et al. 2022

1 Introduction

Mixing and diffusion are important processes in chemical models for planetary atmospheres (see e.g. Moses et al. 2011; Hu et al. 2012; Parmentier et al. 2013; Kreidberg et al. 2014; Charnay et al. 2015; Rimmer & Helling 2016; Zhang & Showman 2018). In the deep layers of gas giant atmospheres, for instance, where densities and temperatures are high, and the chemical reactions are fast, thermo-chemical equilibrium holds. With the falling temperatures and decreasing densities in higher layers, the chemical rates decrease and eventually become smaller than the turbulent eddy-mixing rates. Consequently, the molecular abundances become constant above a certain pressure level – an effect known as ‘transport-induced quenching’ in planetary chemistry. Only at even higher altitudes, once the atmosphere becomes transparent to the stellar UV irradiation, do the molecular abundances start to change again quickly due to photo-dissociation and photo-ionisation, along with other high-energy processes induced by cosmic rays, X-rays, and stellar energetic particles (SEPs, see e.g. Barth et al. 2021).

Another well-known effect in planetary chemistry is that above a certain layer called the ‘homopause’, where the gas-kinetic micro-physical diffusion becomes more important than eddy-diffusion due to turbulent mixing (Hunten 1973; Zahnle et al. 2016). This leads to molecule-dependent scale heights and changes in element abundances above the homopause and the transition into the exosphere.

In protoplanetary disc chemistry, despite the striking physical and chemical similarities to planetary atmospheres, mixing and diffusion have been rarely taken into account. The introduction of diffusive transport terms to the chemical rate equations requires a transition from first-order ordinary differential equations to second-order partial differential equations. Since disc chemistry is at least 2D, this means a considerable numerical and conceptual challenge.

Ilgner et al. (2004) pioneered the research on the impact of eddy mixing on disc chemistry in the frame of (1+1) D viscous disc evolution models for the innermost 10 au. They considered the sulphur molecules CS, SO, SO₂, NS, and HCS⁺ in particular, and concluded that vertical mixing is the most important process among all transport processes considered.

Semenov et al. (2006) used a similar approach but included gas-grain processes. They focused on the CO molecule in the outer disc regions. They concluded that CO can be considerably enhanced by vertical and radial mixing in the cold disc regions, where the dust temperature is lower than about 25 K and CO is expected to freeze out. The mixing rates were derived from an α disc model using α = 10⁻². They discussed the effects of radial, vertical, and ‘full 2D’ mixing in a time-dependent disc model. Unfortunately, since their paper only has a single equation, it remains somewhat unclear how they actually solved the time-dependent partial differential equations.

Semenov & Wiebe (2011) studied the outer parts (>10 au) of the protoplanetary disc around DMTau with given densities and temperatures. They distinguished between ‘responsive’ and ‘unresponsive’ molecules and ice species, based on the calculated column densities. Surprisingly, with regard to the results of this paper, CO, OH, and H₂O, as well as water and ammonia ice were classified as ‘unresponsive’ whereas some sulphur-bearing molecules and polyatomic (organic) molecules were found to be ‘hypersensitive’, meaning that the column densities of the latter molecules were affected by more than two orders of magnitude. We note that these criteria are solely based on changes in the vertical column densities, which are not necessarily proportional to the observable line fluxes. In the TTauri disc models presented in this paper, many of the prominent (sub-) millimetre lines are optically thick in the inner (<100 au) disc regions. In the mid-IR, the lines can also be covered by dust opacity. Hence, the column densities, which are mostly set by the midplane concentrations, are not good predictors of line fluxes in many cases.

Heinzeller et al. (2011) used a simple numerical scheme to include advection caused by inwards accretion and winds, and radial and vertical eddy-mixing in the thermo-chemical disc models developed by Nomura & Millar (2005). In their appendix, it becomes clear that an iterative approach was used, in which, for each disc time integration step, the adjacent concentrations were simply taken from the neighbouring cells at the previous time step and then assumed to be constant, while the central concentrations were advanced in time, point by point. We follow this basic idea in this paper. We note, however, that their Fig. 1 seems to show some artefacts in the form of radial oscillations, which could indicate that their numerical scheme is not fully converging. The results of Heinzeller et al. (2011) are complex, as many effects (vertical and radial mixing, winds, accretion) are discussed at the same time. They concluded that vertical mixing has a considerable impact on the disc chemistry and increases the abundances of many species in the warm molecular layer of the disc observable with the Spitzer Space Telescope and the James Webb Space Telescope (JWST). The viscous parameter for the turbulent mixing was set to α = 10⁻².

Common to the disc papers cited here is that the chemical mechanisms and pathways, which cause the observed changes in the models due to turbulent mixing, do not become fully clear. In this paper, we concentrate on vertical mixing only, and thoroughly explain and test our numerical approach. After an introduction to mixing and diffusion in Sect. 2, we explain our thermo-chemical disc model and test its convergence in Sect. 3. We systematically vary the α-parameter that sets the strength of the turbulent mixing and take the feedback of the changing concentrations on the gas temperature structure into account. We then show and discuss our results in Sect. 4, and demonstrate the main chemical mechanism that changes the concentrations. We continue in Sect. 5 to study the impact caused by turbulent mixing on a few observable quantities, such as ice emission and absorption features in the spectral energy distribution (SED), mid-IR molecular line emission features observable with JWST, and CO J = 2 → 1 channel maps observable with the Atacama Large Millimetre Array (ALMA). We conclude by summarising the main chemical mechanisms found and the effects caused by turbulent mixing in protoplanetary discs in Sect. 6.

2 Chemistry with vertical mixing

2.1 Basic equations

The basic rate equation to solve for each chemical species is

$$\frac{\partial n_i}{\partial t}=P_i-L_i-\frac{\partial {\text{Φ}}_i}{\partial z},$$(1)

where n_i [cm⁻³] is the particle density of species i, P_i, and L_i [cm⁻³ s⁻¹] are the chemical production and destruction rates, and Φ_i [cm⁻² s⁻¹] is the diffusive particle flux according to Hu et al. (2012), Zahnle et al. (2016), and Rimmer & Helling (2016):

$${\text{Φ}}_i=-\left(K+D_i\right)n\frac{\partial }{\partial z}\left(\frac{n_i}{n}\right)+D_i{n}_i\frac{\left(\text{μ}-m_i\right)g_z}{kT},$$(2)

where n = ∑_i n_i is the total particle density. m_i is the mass of gas species i and µ = ∑ n_im_i/n the mean molecular weight. g_z = (GM_★ z)/(r² + z²)^3/2 = z Ω² is the gravitational acceleration in the vertical direction, G the gravitational constant, M_★ the mass of the central star, r the distance from the rotation axis, z the height over the midplane, and Ω the Keplerian circular frequency. We have neglected here the ‘thermal diffusion factor’ α_T, which describes a second-order effect that the lightest molecules tend to diffuse towards the warmest places. The original derivation of Eq. (2) can be found in Hunten (1973) with further reading in Kuramoto et al. (2013).

K [cm² s⁻¹] is the eddy-diffusion coefficient caused by turbulent mixing motions in the z-direction, sometimes also denoted by k_zz. K is generally given by a characteristic mixing velocity times a characteristic length:

$$K=\langle v_z\rangle L,$$(3)

where L is the size of the largest eddy (‘mixing length’) and 〈v_z〉 is the root-mean square of the fluctuating part of the vertical gas velocities. According to the a-disc model of Shakura & Syunyaev (1973), the typical vertical length scale is given by the pressure scale height L ≈ H and the turbulent mixing velocity is approximated as 〈v_z〉 ≈ α_mix c_s:

$$K={\alpha }_{\text{mix}}{c}_{s}H,$$(4)

where $c_{\text{s}}=\sqrt{kT/\text{μ}}$ is the local isothermal sound speed and T the temperature. Typical values used in the literature are α ≈ 10⁻⁵… 10⁻² using different constraints (see reviews by Lesur et al. 2022 and Pinte et al. 2022). While theoretical accretion models often favour higher values to explain the mass accretion rates (see e.g. Manara et al. 2016; Mulders et al. 2017), observational contraints about turbulent line broadening (Flaherty et al. 2015) and dust settling (Pinte et al. 2018) suggest lower values. Eddy mixing transports all molecules in the same way. In contrast, the gas-kinetic micro-physical diffusion coefficient D_i [cm²s⁻¹], due to Brownian motions,

$$D_i=\frac{1}{3}\frac{v_{\text{th,}i}}{{\sigma }_{ia}n}$$(5)

depends on the species i because the collisional cross section σ_ia depends on the ‘size’ of the molecule i and the relative velocity depends on its reduced mass, 1/m_red = 1/m_i + 1/m_j. A simple hard-sphere approximation would be (Jeans 1967)

$$\begin{array}{lll}{D}_i\hfill & =\hfill & \frac{1}{3}{\left({\displaystyle \sum _j\frac{n_j\pi {\left(r_i+r_j\right)}^2}{\sqrt{\frac{8kT}{\pi m_{\text{red}}}}}}\right)}^{-1}\hfill \\ \hfill & =\hfill & \frac{1}{3}\sqrt{\frac{8kT}{\pi m_i}}{\left({\displaystyle \sum _j\frac{n_j\pi {\left(r_i+r_j\right)}^2}{\sqrt{1+\frac{m_i}{m_j}}}}\right)}^{-1}=\frac{b_{ia}}{n},\hfill \end{array}$$(6)

where the effective radii of the molecules r_i can be measured from viscosity experiments (Jeans 1967), or calculated from Lennard-Jones potentials¹. b_ia is called the binary diffusion coefficient [cm⁻¹ s⁻¹] for collisions between species i and a ‘background atmosphere’ denoted by a. Some binary diffusion coefficients for a few relevant combinations of H, H₂, and CO₂ are given in Hu et al. (2012). Most atoms and molecules have effective collisional radii r_i between about 1 Å and 3 Å, as can be seen in Table 1. However, when the diffusion of charged particles and/or diffusion in plasma is considered, the Coulomb cross sections (Langevin cross sections) are larger, by orders of magnitude, than the geometrical cross sections (see e.g. Sect. 2.5.2 in Wiesemann 2014).

2.2 Location of the homopause

Equations (4) and (5) show that while the eddy-diffusion coefficient K is assumed to be constant as a function of height z, the molecular diffusion coefficient D_i scales with 1/n (i.e. it increases with height). In fact, in most parts of the disc, we have K≫D_i. However, in the region above the K = D_i layer, called the homopause (see Fig. 1), the molecules, atoms and ions will tend to arrange themselves with individual scale heights, which could be denoted by ‘gravitational de-mixing’. Above the homopause, the concentrations of the heavy molecules drop quickly, whereas the lighter molecules keep on extending further up. This behaviour can be verified by considering the case of a non-reactive (P_i = L_i = 0) mixture of ideal gases at a constant temperature without eddy mixing (K → 0) in equilibrium, where both n(z) and n_i(z) become time-independent functions of z. The diffusive fluxes vanish and we have

$$\begin{array}{l}{\text{Φ}}_i=0=-D_{i}n\frac{\text{d}}{\text{d}z}\left(\frac{n_i}{n}\right)+D_i{n}_i\frac{\left(\text{μ}-m_i\right)g_z}{kT}\hfill \\ \Rightarrow \frac{\text{d}}{\text{d}z}\ln \frac{n_i}{n}=\frac{\left(\text{μ}-m_i\right)g_z}{kT}=\frac{1}{H}-\frac{1}{H_i}\hfill \\ \Rightarrow n_i\left(z\right)=n_i\left(0\right)\exp \left(-\frac{z}{H_i}\right)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}}n\left(z\right)=n\left(0\right)\exp \left(-\frac{z}{H}\right),\hfill \end{array}$$

with molecular and atmospheric scale heights H_i = kT/(m_ig_z) and H = kT/(µg_z), respectively.

Figure 1 shows the location of the homopause based on the PRODIMO model described in Sect. 3 for a couple of α_mix values. A demonstration molecule with mass m_i = 20 amu and radius r_i = 2 Å is considered, in a mixture of H⁺, H, H₂, He⁺, and He, with radii 0.8 Å, 0.8 Å, 1.366 Å, 1.084 Å, and 1.084 Å, respectively. The surprising conclusion from Fig. 1 is that gravitational de-mixing seems likely to play a significant role in the upper disc layers from which, for example, disc winds are launched, which can therefore be expected to contain only little amounts of heavy elements. We have neglected plasma effects here, meaning we have increased Coulomb collisional cross sections, where the homopause would be located higher in the disc.

Fig. 1 Location of the homopause, where the diffusion coefficients caused by eddy mixing and gas-kinetic motions are equal (K=Di). The colours show the hydrogen nuclei density, n〈H〉, as a function of radius r and relative height z/r. The coloured contour lines show where the homopause is situated when K is calculated with different values of log10 αmix. The dashed lines show where the vertical optical extinction AV,ver equals 100 (black), ten (dark grey) and one (light grey). The dust opacities correspond to our un-mixed model described in Sect. 3. The uneven AV contours are caused by the uneven spatial distribution of the ices, which are taken into account when calculating the opacities.

2.3 Impact of mixing on element abundances

An important question is whether the element abundances can be affected by diffusive transport and become spatially dependent quantities. The element abundances are given by

$$n_{\langle \text{H}\rangle }{ϵ}_k={\displaystyle \sum _i{s}_{i,k}{n}_i,}$$(7)

where s_i,k is the stoichiometric factor of element k in molecule i, ϵ_k are the element abundances with respect to hydrogen, and n_〈H〉 is the total hydrogen nuclei particle density. ϵ_H = 1 by definition. The changes in the element abundances follow from Eq. (1) after multiplication with s_i,k for a selected element k and summation over all chemical species i:

$$\frac{\partial }{\partial t}\left(n_{\langle \text{H}\rangle }{ϵ}_k\right)={\displaystyle \sum _i{s}_{i,k}}\left(P_i-L_i\right)-{\displaystyle \sum _i\frac{\partial }{\partial z}\left(s_{i,k}{\text{Φ}}_i\right).}$$(8)

Chemical reactions conserve the element abundances, and thus the first term in Eq. (8) vanishes. Concerning the diffusive flux Φ_i [cm⁻²s⁻¹], we first consider the case K ≫ D_i where turbulent mixing dominates. From Eq. (2) we then have

$${\text{Φ}}_i=-Kn\frac{\partial }{\partial z}\left(\frac{n_i}{n}\right).$$(9)

The diffusive element fluxes are hence given by

$${\text{Φ}}_k={\displaystyle \sum _i{s}_{i,k}{\text{Φ}}_i}=-{\displaystyle \sum _i{s}_{i,k}Kn\frac{\partial }{\partial z}\left(\frac{n_i}{n}\right)}$$(10)

$$=-Kn\frac{\partial }{\partial z}\left(\frac{{\text{Σ}}_i{s}_{i,k}{n}_i}{n}\right)=-Kn\frac{\partial }{\partial z}\left(\frac{n_{\langle \text{H}\rangle }{ϵ}_k}{n}\right),$$(11)

where K does not depend on molecule. The last term in this equations can be expressed by

$$\frac{\partial }{\partial z}\left(\frac{n_{\langle \text{H}\rangle }{ϵ}_k}{n}\right)=\frac{\partial }{\partial z}\left(\frac{n_{\langle \text{H}\rangle }}{n}\right){ϵ}_k+\frac{n_{\langle \text{H}\rangle }}{n}\frac{\partial {ϵ}_k}{\partial z}.$$(12)

In general, the total particle density depends on the chemical composition whereas n_〈H〉 does not, see further discussion of this point in Sect. 2.3.1. However, if we assume n_〈H〉 ∝ n anyway, the first term in Eq. (12) vanishes and we find

$$\frac{\partial }{\partial z}\left(n_{\langle \text{H}\rangle }{ϵ}_k\right)=\frac{\partial }{\partial z}{\text{Φ}}_k=\frac{\partial }{\partial z}\left(K{n}_{\langle \text{H}\rangle }\frac{\partial {ϵ}_k}{\partial z}\right).$$(13)

Equation (13) states that eddy diffusion works against spatial gradients of element abundances. No matter how complicated the chemical processes are, eddy diffusion will always diminish such gradients $\left({\text{Φ}}_k\to 0\Rightarrow \frac{\partial {ϵ}_k}{\partial z}\to 0\right)$, making sure that the element abundances will eventually become constant throughout the disc.

This result seems a bit counter-intuitive at first; when we consider a particular molecule only, for example, that is abundant in deep layers and mixed up where it dissociates, one would expect this process to enrich the higher layers with the elements of which the molecule was composed. However, this is not how eddy mixing works. Mixing rather exchanges molecules for other molecules in the local neighbourhood, and if we properly sum up the total effect on the element abundances, via Φ_k = ∑_i s_i,k Φ_i, the net effect is exactly zero, unless there is a spatial gradient in element abundances in the first place, in which case eddy diffusion will reduce that gradient.

However, this conclusion is only valid as long as K is independent of molecule to be mixed by diffusion. Once we add D_i with D_i < K (above the homopause), we can no longer pull out the net diffusion coefficient (K + D_i) from the sum in Eq. (10). In that case, diffusion will change the element abundances. An example of such long-term effects by gravitational de-mixing was already described in Sect. 2.2.

Another relevant case occurs when we consider the difference between gaseous and ice species. The ice species are transported with the dust grains, and the grains are only moved by the larger turbulent eddies. In order to arrive at an effective dust particle diffusion coefficient, the advective effect of all individual turbulent eddies has to be averaged. This procedure has been carried out, with different methods and approximations, by Dubrulle et al. (1995), Schräpler & Henning (2004), and Youdin & Lithwick (2007). The result of Schräpler & Henning (2004, see their Eq. (27)), is

$$K_{\text{ice}}=\frac{K}{1+St},$$(14)

where St is the Stokes number given by

$$St=\frac{{\tau }_{\text{stop}}}{{\tau }_{\text{eddy}}}=\frac{\frac{{\rho }_{\text{d}}a}{\rho v_{\text{th}}}}{\frac{H}{{\alpha }_{\text{mix}}{c}_s}}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_{\text{d}}a{\alpha }_{\text{mix}}}{\rho H},$$(15)

where τ_stop is the stopping distance, τ_eddy the turnover timescale of the largest turbulent eddy, ρ_d the dust material density, the thermal velocity, c_s the sound velocity, H the pressure scale height, and α the grain radius. The factor $\sqrt{\pi /8}$ is because $c_{\text{s}}^2=kT/\text{μ}$ but $v_{\text{th}}^2=8kT/\left(\pi \text{μ}\right)$.

It is questionable, however, how far Eq. (14) describes the actual dynamical motions of dust grains with large Stokes numbers in discs. The main reason is that such particles will start to migrate inwards (i.e. they develop a systematic velocity with respect to the gas), which goes beyond the mixing approach. Such motions have been shown to be able to change the elemental abundances across ice lines of key volatiles such as CO (Krijt et al. 2018). Furthermore, it is not obvious how to choose α as we actually need the cumulative effect of all grains of different sizes on the diffusive flux, where it also must be debated how the ice is distributed onto the bare grains as a function of size, and how the size dust distribution is altered by ice formation (see e.g. Arabhavi et al. 2022).

Assuming that we can determine K_ice anyway, and assuming that all vertical element fluxes Φ_k must vanish everywhere in the disc in the time-independent limiting case, we can derive a local condition for each element k from Eq. (10):

$${\text{Φ}}_k=-Kn\frac{\partial }{\partial z}\left(\frac{{\text{Σ}}_{\text{gas}}{s}_{ik}{n}_i}{n}\right)-K_{\text{ice}}n\frac{\partial }{\partial z}\left(\frac{{\text{Σ}}_{\text{gas}}{s}_{ik}{n}_i}{n}\right)=0.$$(16)

Defining the element abundances present in the gas and ice, respectively, as

$${ϵ}_k={ϵ}_k^{\text{gas}}+{ϵ}_k^{\text{ice}},n_{\langle \text{H}\rangle }{ϵ}_k^{\text{gas}}={\displaystyle \sum _{i,\text{gas}}{s}_{ik}{n}_i,n_{\langle \text{H}\rangle }{ϵ}_k^{\text{ice}}={\displaystyle \sum _{i,\text{ice}}{s}_{ik}{n}_i,}}$$(17)

we find the following result that is independent of α_mix:

$$\frac{\partial }{\partial z}\left(\frac{n_{\langle \text{H}\rangle }{ϵ}_k^{\text{gas}}}{n}\right)+\frac{K_{\text{ice}}}{K}\frac{\partial }{\partial z}\left(\frac{n_{\langle \text{H}\rangle }{ϵ}_k^{\text{ice}}}{n}\right)=0.$$(18)

2.3.1 Small Stokes numbers

In the limiting case of small dust radii α and/or large gas densities ρ, we have St ≪ 1 (i.e.well-coupled grains), and hence K_ice ≈ K. Thus, from Eq. (18), we obtain

$$\frac{n_{\langle \text{H}\rangle }{ϵ}_k}{n}=\text{const},$$(19)

previously referred to as the case of ‘constant element abundances’. Clearly, when all species are transported in equal ways, and hence mixed with the same efficiency, the element abundances cannot change. However, Eq. (19) actually states a somewhat different result. Diffusion will make sure that, in the long run, the concentrations of elements as measured with respect to the total gas particle density n become constant. At the locations of the major phase transitions H⁺ ↔ H and H ↔ H₂, where the mean molecular weight changes by roughly a factor of two in both cases, the element abundances must change as well, in the opposite direction, approximately ϵ_k ∝ 1/µ. Otherwise, when µ is constant, we have n ∝ n_〈H〉 and indeed ϵ_k = const. Considering a complete vertical column, from ionised down to molecular conditions, Eq. (19) states that the helium abundance, for instance, should be about four times larger in the ionised regions than in the molecular regions.

2.3.2 Large Stokes numbers

In the limiting case of large Stokes numbers (a → ∞ and/or ρ → 0), we have K_ice ≪ K, and thus from Eq. (18), we have

$$\frac{n_{\langle \text{H}\rangle }{ϵ}_k^{\text{gas}}}{n}=\text{const}.$$(20)

Equation (20) states a remarkable result. If an element can freeze out anywhere in the column, a rather low level of ${ϵ}_k^{\text{gas}}$ will be enforced there, and this low level of ${ϵ}_k^{\text{gas}}$ must be identical everywhere in the column to make the diffusive fluxes vanish, irrespective of the local ${ϵ}_k^{\text{ice}}$. This effect is commonly know as ‘cold trap’ (see e.g. Krijt et al. 2016 and Öberg & Bergin 2021). Since a cold surface cannot move, but air does, the gaseous water molecules in a room will continue to diffuse towards that cold surface and condense onto it, until eventually the supersaturation of the air over the cold surface approaches unity. If the surface temperature is kept lower than the air temperature, this will eventually cause an under-saturation of water vapour in the warm gas in the room.

In a protoplanetary disc, we have the same principle situation. Certain molecules will freeze out quickly in the cold midplane, establishing a very low local ${ϵ}_k^{\text{gas}}$. After just a few vertical mixing timescales (see Eq. (32)), we would therefore expect most molecules in the upper layers to simply be mixed down into the midplane and freeze out there, leaving behind very low concentrations also in the upper layers. However, this conclusion would be in obvious conflict with ALMA observations, for example (Pinte et al. 2018; Rab et al. 2020), where we clearly see CO molecules only above the icy midplane at 100 au distances.

We conclude that the cold trap effect cannot work in discs at full effect, possibly because: (i) the actual αmix in discs is much lower than 10⁻², such that the mixing timescale according to Eq. (32) approaches the disc lifetime. For an assumed disc lifetime of 5 × 10⁶ yr, this would suggest α_mix < 3 × 10⁻⁵ at 100 au, and even smaller α_mix at smaller radii; (ii) the effective Stokes numbers (after averaging over dust size distribution) result to be small in Eq. (14), such that the limiting case described by Eq. (20) is not reached; or (iii) a combination of both.

But even if this means that Eq. (20) is not valid, the cold trap effect could still cause significant deviations from the case of constant element abundances (Eq. (19)) over long timescales, although these will be difficult to prove observationally. In principle, one can solve Eq. (18) in an iterative way, to make it consistent with the ice chemistry, in order to obtain ${ϵ}_k^{\text{gas}}\left(z\right)$ and ${ϵ}_k^{\text{ice}}\left(z\right)$.

Further investigations are clearly required here, both on the influence of molecular diffusion on the upper disc layers and on the cold trap effect, where Eq. (18) could serve as a useful starting point. However, in the remaining sections of this paper, we do not discuss these complications any further, but concentrate on the most basic case, where we have vanishing molecular mixing D_i → 0 and no difference between the diffusive transport of gas and dust K_ice → K, in which case the element abundances can be assumed to be constant throughout the disc. We also neglect here the vertical gradients of the mean molecular weight µ.

2.4 Implementation of mixing rates in PRODIMO

A direct numerical solution of the reaction-diffusion equations (Eq. (1)), which is a system of coupled stiff partial differential equations of second order, is difficult and can be computationally very demanding (see Semenov et al. 2006, Semenov & Wiebe 2011, and Hu et al. 2012). Here, we follow an iterative approach similar to Heinzeller et al. 2011 and Rimmer & Helling 2016. We consider height grid points {z_k | k = 1, …, K}. As shown in Appendix A, the divergence of the diffusive particle flux i at point k can be numerically discretised as

$$-\frac{\partial {\text{Φ}}_{i,k}}{\partial z}=A{n}_{i,k+1}-B{n}_{i,k}+C{n}_{i,k-1},$$(21)

where n_i,k is the particle density of species i at height grid point k, and where the coefficients A, B, and C are given by the geometry of the grid points, the eddy and molecular diffusion coefficients, K and D_i, and the total particle densities n = Σn_i at the neighbouring height grid points k −1, k, k + 1. Our approach is then very simple. We consider the neighbouring concentrations n_i,k+1 and n_i,k−1 as fixed during the solution of the chemistry at height point k. Hence, we can solve

$$\frac{\partial n_i}{\partial t}={\overline{P}}_i-{\overline{L}}_i,$$(22)

with the existing time-independent chemistry solver in the PRODIMO model, where

$${\tilde{P}}_i=P_i+A{n}_{i,k+1}+C{n}_{i,k-1}$$(23)

$${\tilde{L}}_i=L_i+B_{n_{i,k}}.$$(24)

Since PRODIMO solves the chemistry and the heating and cooling balance on all grid points in a particular order (first downwards – then outwards, following the two main irradiation directions), the constants n_i,k+1 are taken from the same chemistry iteration (from the last point above just solved), whereas the constants n_i,k−1 are taken from the previous chemistry iteration. In Appendix A.4 we show that we can actually solve simple diffusion problems with this basic iterative scheme. It is not very efficient or elegant, but it allows us to quickly implement vertical mixing in PRODIMO without changing the code structure at all.

The contribution of the diffusive mixing to the production and destruction of species i on grid point k is formulated in terms of two additional pseudo reactions. The ‘mixing in’ (MI) reaction is a spontaneous process with constant rate, with no entries in the Jacobian. The ‘mixing out’ (MO) reaction has a rate that is proportional to n_i,k:

$$\text{MI}:\to mo{l}_i\left(k_{\text{MI}}=A{n}_{i,k+1}+C{n}_{i,k-1}\right)$$(25)

$$\text{MO:\hspace{0.17em}}{\text{mol}}_i\to \left(k_{\text{MO}}=B\right),$$(26)

where k_MI and k_MO are the rate constants applied in the chemical network. The iteration of the chemistry is then continued until the solution relaxes towards a steady state where all particle densities at all grid-points n_i,k do not change any further. In that case, the neighbouring concentrations n_i,k+1 and n_i,k−1 correctly describe the spatial gradients in the system and we have reached the time-independent solution of

$$0=P_i-L_i-\frac{\partial {\text{Φ}}_i}{\partial z}$$(27)

for each chemical species i on all grid points k.

3 The disc model

The physical and chemical setup of the disc model considered in this paper was described in detail in Woitke et al. (2016) and we refer the reader to their Table 3. The central star is a 2 Myr old T Tauri star with a mass of 0.7 M_⊙ and a luminosity of 1 L_⊙, the disc has a mass of 0.01 M_⊙, with a gas/dust ratio of 100, and the dust grains have sizes between 0.05 µm and 3 mm, with a power law size-distribution of index −3.5. The dust component, approximated by 100 size bins, is settled according to the prescription of Dubrulle et al. (1995) with α_setüe = 10⁻². In a consistent α-model, α_mix = α_settle should be used, but here we chose to keep α_settle constant to study the changes that are solely caused by the introduction of mixing to the chemistry, using different choices of α_mix. All other parameters, including disc shape, mass, radial extension, dust material and opacity, as well as the cosmic-ray, UV, and X-ray irradiation parameters, are listed in Table 3 in Woitke et al. (2016).

We use the large DIANA standard chemical network (Kamp et al. 2017) in this work, which has 235 chemical species, with reaction rates mostly based on the UMIST 2012 database (McElroy et al. 2013), but including a simple freeze-out and desorption ice chemistry, X-ray processes including doubly ionised species (Aresu et al. 2011), excited molecular hydrogen, and polycyclic aromatic hydrocarbons in five different charging states; altogether there are 2843 reactions. Surface chemistry (Thi et al. 2020a,b) is not included.

Since the up-mixing of ices formed in the disc midplane is a strong effect in the models presented in this paper, we used the consistent calculation of ice opacities as implemented by Arabhavi et al. (2022), including H₂O-ice, CO₂-ice, NH₃-ice, CH₄-ice, CO-ice, and CH₃OH-ice. We assume that these ices form a layer of equal thickness on top of each refractory grain (q = 0 in Eq. (3) of Arabhavi et al. 2022). The thickness and composition of this ice layer depends on the position in the disc and is adjusted according to the local results of the chemistry. Thus, the dust opacity used in the continuum radiative transfer depends on the results of the chemistry, which requires global iterations to solve. However, since we do up to 1000 global iterations anyway, to solve the reaction-diffusion equations, and each iteration includes continuum radiative transfer, and chemistry and heating and cooling balances, these additional feedbacks are automatically resolved.

We chose to use 80 radial grid points and 120 vertical grid points for our models, a compromise between having acceptable computation times and the need to resolve the vertical structures as good as possible, to study the vertical mixing. The computation time for one global iteration of one model is about 40 min user time on 32 CPUs. For 1000 global iterations, we needed about 1 month of user time. All computations were performed with PRODIMO git version 2f79b977 (11 February 2022).

Although it would be relatively straightforward to include radial mixing as well, we have only included vertical mixing in this paper. We also chose to concentrate entirely on turbulent eddy mixing and drop all terms for gas-kinetic molecular diffusion (D_i = 0), both of which might be investigated in forthcoming papers.

4 Results

4.1 Convergence of the numerical method

Before we discuss the physical and chemical results of our disc models with vertical mixing, we first need to make sure that our simple iteration scheme actually converges. We measure the convergence by means of the following two criteria, which evaluate the differences in the various chemical concentrations between the current and the previous global iteration in different ways:

$$\text{crit\hspace{0.17em}1=}\underset{i}{\max }\left\{\left|{\log }_{10}{m}_{i,j}\text{/}{m}_{i,j-1}\right|\right\}$$(28)

$$\text{crit\hspace{0.17em}}2={\left(\frac{1}{N_{\text{p}}{N}_{\text{sp}}}{\displaystyle \sum _{k=1}^{N_{\text{p}}}{\displaystyle \sum _{i=1}^{N_{\text{sp}}}{\left({\log }_{10}{n}_{k,i,j}/n_{k,i,j-1}\right)}^2}}\right)}^{1/2};$$(29)

where m_i,∙j is the disc volume integrated mass of species i at iteration j, and n_k∙i,∙j is the particle density of species i on a 2D grid point k at iteration j. N_p is the number of 2D grid points, and N_sp is the number of chemical species in the network.

Figure 2 plots these two convergence criteria as a function of iteration number for the model with medium strong mixing, α_mix = 10⁻³. The figure shows that, due to our very simple iteration method, the convergence is rather slow. We need a couple of hundred iterations to bring the relative changes down to about 0.003 dex, but even at 1000 iterations, there are still fluctuations of a similar order of magnitude.

Figure 3 shows some details about how the temperature and a few selected molecular concentrations change during the course of the iterations, considering a vertical cut through the disc at 10 au. We observe that the chemical structure changes considerably within the first few iterations, but then the changes per iteration settle down, and eventually the differences between iteration 400 (dash-dotted lines) and iteration 1000 (full lines) are barely recognisable on the right side of Fig. 3. The temperature structure also changes. The gas and dust temperatures cool down because: (i) more OH and H2O form, which lead to increased line cooling; and (ii) the optical depths increase due to the presence of icy grains in upper regions. We conclude that our method does converge, albeit very slowly, which agrees with the expected convergence behaviour as discussed in Appendix A.4.

Fig. 2 Measuring the convergence of the iterative solution technique to solve the chemistry with vertical mixing, here for the model with medium strong mixing, αmix = 10−3. See text for explanations.

4.2 Impact on the chemical and temperature structure

4.2.1 Overview of principle effects

Figure 4 shows our resulting temperature and chemical structures. Each plot shows a sequence of well-iterated results for α_mix ∈ {0,10⁻⁵, 10⁻⁴, 10⁻³, 10⁻²} using different line styles. The results are shown by three vertical cuts at constant distance to the star, at 1 au, 10 au, and 100 au, each as a function of hydrogen nuclei column density N_〈H〉 [cm⁻²]. The figure shows a number of remarkable changes in the chemical structure that are caused by the inclusion of vertical mixing, which will be detailed in the forthcoming sections.

First, the ices that predominantly form in the midplane are mixed upwards by a distance that corresponds to several orders of magnitude in column density.

Second, molecules that are mainly formed in the warm intermittent disc layers, where the temperatures are a few hundred K, are mixed both downwards and upwards.

Third, the increased molecular abundances of OH and H₂O cause reinforced line cooling at column densities between about 10¹⁹ and 10²² cm⁻², which lowers the gas temperature and improves the conditions for the formation and survival of more molecules. These layers are the regions directly observable with infrared line spectroscopy. We see that the maxima of the OH and H₂O concentrations shift leftwards and upwards by about one to two orders of magnitude in these plots. These changes are larger for more distant regions and/or a larger mixing parameter α_mix.

Forth, neutral atoms are mixed up into the highest layers, where they absorb X-rays and FUV radiation, which leads to additional absorption of UV and X-rays, causing an increase in temperature, and an increase in the shielding for the upper disc layers at larger radii.

Fifth, ions such as O⁺, which are predominantly formed in high layers, are increasingly able to meet molecules and even icy grains, which are mixed upwards. This leads to a very rich ion-molecule chemistry, with new chemical pathways that are normally unimportant in models without mixing.

Fig. 3 Monitoring the convergence of the iterative solution method along a vertical cut at r = 10 au as a function of the vertical hydrogen nuclei column density $N_{\langle \text{H}\rangle }={\displaystyle {\int }_z^{\infty }{n}_{\langle \text{H}\rangle }\left(z\right)\text{\hspace{0.17em}d}z}$ for the model with medium strong mixing, αmix = 10−3. The upper plots show the changing gas and dust temperature as a function of iteration number. The lower plots shows the concentrations of a few selected oxygen-bearing gas and ice species. The left row shows the first few iterations and the right row shows the convergence towards iteration 1000, which is considered as the final solution. The different full, dashed and dotted lines correspond to the different numbers of iterations carried out as indicated.

4.2.2 Beating UV photo-desorption by mixing

At large radial distances, ~100 au, the model with the highest mixing strength, α_mix = 10⁻², exhibits some rather peculiar results (see Fig. 5). The concentrations of the thermally most stable ice species such as H₂O#, NH₃#, CO₂#, OH#, and C₂H₄# (and some others not shown) become apprximately constant and remain high even at the very top of the model, where they face the unattenuated stellar radiation field. This is in contrast to the less thermally stable ices such as CO# and N₂# (and others).

In models without mixing, the ices are generally not stable in the optically thin disc regions because of UV photo-desorption. However, in the model with strong mixing, α_mix = 10⁻², the upwards mixing rates seem to be large enough to replenish the ice species faster than the UV desorption can destroy them. This is in contrast to thermal desorption. Once the ices become thermally unstable, due to the increasing dust temperature in higher layers, they will sublimate very quickly and mixing cannot change that. We see this behaviour very clearly in the upper part of Fig. 6, where the upper boundary of CO# does not surpass the dust temperature T_d = 25 K line, and the H₂O# upper boundary does not surpass the T_d = 150 K line, which corresponds to the desorption energies of these ices.

With an assumed stellar luminosity L_★ = 1 L_⊙ and f_UV = 0.01, we obtain a UV flux of about 6 × 10¹⁰ photons cm⁻² s⁻¹ between 912 Å and 2050 Å at r = 100 au in the optically thin part high above the disc, which is about a factor of χ = 300 larger than the standard Draine flux in the ambient interstellar medium, F_Draine = 1.921 × 10⁸ photons cm⁻² s⁻¹ (Draine & Bertoldi 1996). With the symbols introduced in Sect. 5 of Woitke et al. (2009), the UV photo-desorption rate of an ice species is

$$\frac{\text{d}{n}_{i\#}}{\text{d}t}=-\pi \langle a^2\rangle n_{\text{d}}{Y}_i\chi F_{\text{Draine}}\times \{\begin{array}{cc}\frac{n_{i\#}}{n_{\text{act}}},& n_{\text{tot}}^{\text{ice}}<n_{\text{act}}\\ \frac{n_{i\#}}{n_{\text{tot}}^{\text{ice}}},& n_{\text{tot}}^{\text{ice}}\ge n_{\text{act}}\end{array},$$(30)

where n_d is the number density of dust particles, π〈a²〉n_d the total dust cross section per cm³, Y_i the UV photo-desorption yield, and n_act = 4π〈a²〉n_d 2n_surf the total number of active adsorption sites on all grains per cm³. Following Aikawa et al. (1996), we assume that the uppermost two ice layers are active, meaning that frozen molecules can be hit directly by UV and can escape freely from these layers, with a surface density of n_surf = 1.5 × 10¹⁵ sites per cm² of solid surface.

The first case in Eq. (30) corresponds to a situation where the grain surface is almost unoccupied and the ratio n_i#/n_act expresses the probability that a UV photon absorbed by a grain hits a site occupied by ice species i. In the second case, the surface is more than fully occupied and the ratio $n_{i\#}/n_{\text{tot}}^{\text{ice}}$ expresses the probability of finding ice species i in the ice mixture present, where the total ice density is $n_{\text{tot}}^{\text{ice}}={\text{Σ}}_i{n}_{i\#}$. Hence, thick ice mantles can only be reduced layer by layer via photo-desorption, meaning at a constant rate (‘zero-order photo-desorption’). Figure 5 shows that, according to the assumed dust size distribution and dust/gas ratio, the grains in the midplane are covered by about $N_{\text{layers}}=n_{\text{tot}}^{\text{ice}}/\left(4\pi \langle a^2\rangle n_{\text{d}}{n}_{\text{surf}}\right)\approx {10}^4$ ice layers.

The timescale to completely UV photo-desorb such grains in the optically thin upper disc layers can be obtained by summing up Eq. (30) over all ice species i, applying the second case

$${\tau }_{\text{UV,desorb}}=\frac{n_{\text{tot}}^{\text{ice}}\text{/}{n}_{\langle \text{H}\rangle }}{\frac{\pi \langle a^2\rangle n_{\text{d}}}{n_{\langle \text{H}\rangle }}\overline{Y}\chi F_{\text{Draine}}}\approx \left(35000\text{\hspace{0.17em}yr}\right)\left(\frac{300}{\chi }\right),$$(31)

where we have assumed a total initial ice concentration of $n_{\text{tot}}^{\text{ice}}/n_{\langle \text{H}\rangle }\approx {10}^{-3.5}$, amean photo-desorption yield of Ȳ ≈ 10⁻³, and a total dust surface per H-atom of 4π〈α²〉n_d/n_〈H〉 ≈ 2 × 10⁻²³ cm² (which can be compared with Fig. 5). We note that χ ≈ 300 is what we find in our model at r = 100 au. χ ∝ 1/r² (i.e. τ_UVdesorb ∝ r²) is valid until the interstellar UV photons become more relevant than the diluted stellar UV, which happens at about 1700 au in this case.

In comparison, the turbulent eddy-diffusion timescale is

$${\tau }_{\text{mix}}=\frac{H_p}{{\alpha }_{\text{mix}}{c}_{\text{s}}}\approx \left(18000\text{\hspace{0.17em}yr}\right){\left(\frac{r}{100\text{\hspace{0.17em}au}}\right)}^{1.4}\left(\frac{{10}^{-2}}{{\alpha }_{\text{mix}}}\right),$$(32)

where we assumed a scale height of H_p = (10 au)(r/100 au)^1.15, a temperature structure as T = (20 K)(r/100 au)^−0.5, a constant mean molecular weight² of µ = 2.3 amu, and $c_{\text{s}}=\sqrt{kT/\text{μ}}$.

The conclusion from this timescale comparison is that only for a very large α_mix and large radii, the up-mixing of icy grains can indeed be faster than UV desorption. The two timescales meet at about 30 au in our model with α_mix = 10⁻², which is in agreement with where we start to see icy grains in the uppermost disc layers in our models. In all other cases, the icy grains stay confined to the midplane due to UV photo-desorption. However, the vertical spread of the icy grains is distinctively wider in comparison with un-mixed models, and strongly depends on parameter α_mix. The entirely flat concentrations of the most stable ices in the α_mix = 10⁻² model are still somewhat surprising. It is possible that the strong increase in the gas temperature, from a few tens of K to a few hundred K, causing larger c_s and thus shorter mixing timescales, is responsible for this behaviour.

Fig. 4 Resulting temperatures and chemical concentrations as a function of vertical hydrogen column density N〈H〉, for radial distances 1 au, 10 au, and 100 au, from left to right. The top row shows the gas and dust temperatures, the second row the main hydrogen species, the third row the oxygen chemistry, and the last row the carbon chemistry. Only a few selected species are plotted. In each part of the figure, the various dashed and dotted lines correspond to the mixing parameter αmix, as explained in the legends.

Fig. 5 Vertical spreading of icy grains at r = 100 au for different mixing strengths αmix. Results are shown as a function of the vertical hydrogen nuclei column density N〈H〉. The top row of plots shows the computed gas and dust temperatures, Tg and Td. The middle row of plots shows the bare grain surface area per H nucleus = 4π〈a2〉nd/n〈H〉 after dust settling in blue, and the number of ice layers on each grain Nlayers in red. The lower plots show the concentrations of selected ice species.

Fig. 6 Upper row of plots shows the location of the ices, in a series of models with increasing mixing parameter αmix. The icy regions are hashed when the respective ice concentration, for example nH2O#/n〈H〉 for water ice, is larger than 10−55. The background colours show the hydrogen nuclei density n〈H〉 from 102 cm−3 (black) to 1016 cm−3 (red). The dashed black contour lines show the radial visual extinction AV,rad = 1 and vertical visual extinction AV,ver = 1, and the red dashed lines show two dust temperature contours, 25 K and 150 K. The lower row of plots shows the respective SEDs for selected inclination angles at a distance of 140 pc. The yellow lines indicate the adopted photospheric spectrum of the naked star. We note that the absorption feature around 2.5 µm is due to molecular absorption in the stellar photosphere, whereas all other narrow absorption and emission features are due to ices, except for the broad 10 µm/20 µm silicate emission features.

4.2.3 Chemical pathways

In order to better understand the chemical processes induced by mixing, we examine one point in our disc model with α_mix = 10⁻³ at 10 au, located at the maximum of the concentration of OH, that is, the maximum of the red curve in Fig. 3 (right part) at N_〈H〉 ≈ 1(+20)cm⁻². The notation a(+b) means a × 10^+b. Here we have a hydrogen nuclei density of n_〈H〉 = 2.4(+7) cm⁻³. The gas and dust temperatures are T_g = 347 K and T_d = 126 K, and the point is optically thin with radial and vertical visual extinctions $A_{\text{V}}^{\text{rad}}=0.024$ and $A_{\text{V}}^{\text{ver}}=6.3\left(-4\right)$, respectively. The most abundant (‘available’) species are H [9.5(−1)], H₂ [2.3(−2)], O [3.0(−4)], e⁻ [1.5(−4)], C⁺ [1.3(−4)], N [7.9(−5)], and H⁺ [2.1(−5)], where the numbers in squared brackets are the concentrations n_i/n_〈H〉.

Figure 7 shows how the up-mixing of H₂ fuels the ion-neutral chemistry in the upper disc regions, resulting in an enhanced production of OH and H₂O. The left side of the figure shows the rates and concentrations in the case without mixing, where OH is mainly produced via the radiative association reaction O + H → OH + hv, water is mainly produced by the follow-up radiative association reaction OH + H → H₂O + hv, and H₂ is produced partly by the formation on grains and partly via H + H⁻ → H₂ + e⁻.

On the right side of Fig. 7, which shows the results at the same point for the case with mixing α_mix = 10⁻³, the H2-concentration is enhanced by a factor of about 1000. Here, the starting point for the production of OH and H₂O is the reaction

$${\text{H}}_2+{\text{O}}^{+}\to {\text{OH}}^{+}+\text{H,}$$(33)

where the O⁺ is produced by charge exchange with H⁺, which results from secondary X-ray reactions with super-thermal electrons H + e⁻ — H⁺ + e⁻ + e⁻. Although reaction (33) is only the fourth most important destruction reaction of H₂ (after out-mixing, photo-dissociation, and reaction with CH⁺), it triggers the production of OH and H₂O via the formation of OH⁺. Quick H-addition reactions with H2 subsequently produce H₂O⁺ and H₃O⁺, before recombinations with free electrons create the neutral molecules OH and H₂O. This H₂O formation path is nomally initiated by reactions of O with H₃⁺, for example O + H₃⁺ → OH⁺ + H₂, but here we find reaction (33) to be more efficient in producing OH⁺, because the mixing leads to a coexistence of H₂ and O⁺. Since H₂ is required several times along these formation pathways, even a small increase in the H₂ concentration, caused by up-mixing, can lead to a much higher production of OH and H₂O molecules. The resulting OH and H₂O concentrations are higher by about 2.5 and six orders of magnitude, respectively, at this point, when mixing is taken into account with α_mix = 10⁻³. The reaction cycles are completed by photo-dissociation and dissociative recombination reactions, which again destroy the molecules and replenish the atoms.

It is worth noting that H₂ is the only molecule that is strongly affected by mixing directly in Fig. 7. All other molecules shown in Fig. 7 have the same reaction coefficients for the MI and MO quasi reactions, yet their low concentrations render these rates insignificant. The only other species directly affected by mixing are those atoms and ions listed as ‘abundant’ in the beginning of this paragraph. This finding can be explained by considering the chemical formation and destruction timescales (particle densities over rates). For low-abundance molecules, particle densities are low but rates are similar (i.e. the chemical timescales are short). Only for the most abundant molecules can the chemical timescale exceed the mixing timescale. The other species with low abundances and short chemical timescales rather tend to adjust their concentrations quickly to these changed circumstances.

Next, we look at the hydrocarbon chemistry. We consider a point in our α_mix = 10⁻³ disc model at r = 0.15 au, at a height where the C₂H₂ concentration is steeply rising towards a plateau in the midplane, that is, at a vertical hydrogen column density of about N_〈H〉 ≈ 1.0(+23) cm⁻² on the left plot in Fig. 4. Points such as these are flagged by PRODIMO as the main emission region for the C₂H₂ 13.7 µm emission feature observable with Spitzer/IRS and JWST/MIRI (see Fig. 8). Here, we have a hydrogen nuclei density of n_〈H〉 = 1.0(+13)cm⁻³, the gas and dust temperatures are T_g = 353 K and T_d = 354 K, and the point is optically thick in the radial direction, $A_{\text{V}}^{\text{rad}}=194$, but semioptically thin in the vertical direction, $A_{\text{V}}^{\text{ver}}=3.1$, still resulting in a total UV strength of χ = 20000. The most abundant species are H₂ [0.5], H [5.6(−3)], H₂O [1.7(−4)], CO [1.3(−4)], and N₂ [3.6(−5)], with no ices.

Figure 9 shows how the up-mixing of CH₄ followed by photo-dissociation leads to much higher concentrations of all hydrocarbon molecules. The hydrocarbon chemistry forms an almost self-contained reaction subsystem that is only weakly connected to other chemicals. Comparing the two depicted cases (with and without mixing), we see that the same neutralneutral, ion-neutral, and photo-dissociation reactions connect the different hydrocarbon species in similar ways, forming sustained reaction cycles. The various reactions with O or O₂ in these cycles form C-O bonds that, once formed, can hardly be destroyed again. We denote these oxidisation reactions as ‘burning’ in the following paragraphs. On the left side of Fig. 9 (without mixing), these reactions consume almost all the hydrocarbon molecules, until an equilibrium with neutral carbon atoms is established, on a very low level n_C/n_〈H〉 = 3(−14), where further burning of the H–C molecules is balanced by X-ray destruction of CO and X-ray-induced charge exchange (XA) with C⁺⁺.

On the right side of Fig. 9 (with mixing), methane molecules (CH4) that formed further down in the UV-protected disc layers are mixed up into the semi-optically thin region considered here, where they are quickly photo-dissociated as

$${\text{CH}}_4+hv\to \{\begin{array}{ll}{\text{CH}}_3+\text{H}\hfill & \left(15.5\%\right)\hfill \\ {\text{CH}}_2+{\text{H}}_2\hfill & \left(69.0\%\right)\hfill \\ \text{CH}+{\text{H}}_2+\text{H}\hfill & \left(15.5\%\right)\hfill \end{array}$$(34)

with branching ratios from the UMIST 2012 database (McElroy et al. 2013). This introduces rates of the order of 10⁻¹ cm⁻³ s⁻¹ to the hydrocarbon system. The system reacts to this input by increasing all concentrations of hydrocarbon molecules by many orders of magnitude, from about four decades for C to 15 decades for C₃H. The concentration of acetylene (C₂ H₂), which turns out to be the most abundant hydrocarbon molecule in this case, increases from 3.4(−19) to 8.3(−7), creating a peculiar chemical situation where carbon can be found in opposite redox states (C₂H₂ and CO₂) in similarly high concentrations.

Again, only the abundant molecules shown with big red arrows in Fig. 9 (CO, CO₂ and CH₄) are significantly affected by mixing directly. All other molecules just adjust to those circumstances quickly. For example, C₂H₂ has a chemical relaxation timescale of about n_C2H2/rate ≈ 8(−7) × 1(+13) cm⁻³/(3(+4) cm⁻³ s⁻¹) = 270 s at this point, which is much shorter than the mixing timescale ≈ 2 yr (Eq. 31).

In summary, the CH₄ molecules that are mixed up from deeper layers by eddy-diffusion are quickly photo-dissociated once they reach the borderline optically thin layers. The products of these photo-processes react further to create a rich hydrocarbon chemistry in an oxygen-rich situation, where the various hydrocarbon molecules are ultimately burnt on timescales of a few years at r = 0.15 au, forming C-O bonds. We actually see in Fig. 9 that there is an over-production of CO and CO₂ that is mixed out.

Fig. 7 Chemical pathways to create OH and H2O molecules without mixing (left) and with αmix = 10−3 (right) at one selected point in the disc at r = 10 au located in the upper disc layers (see text). The green numbers represent the particle concentrations and the blue numbers represent the rates [cm−3 s−1] in the time-independent case. The big red ‘IN’ and ‘OUT’ arrows on the right denote the in-mixing and out-mixing rates for H2. Wiggled arrows indicate photo-dissociation. Light green arrows on the left mark the H2 formation on grains, and the radiative association reactions O + H → OH + hv and OH + H → H2O + hv. ${\text{H}}_2^{★}$ denotes electronically excited molecular hydrogen. The figure does not follow up the production of CH2+, CO+, and HCO+ marked in blue. The notation a(−b) means a × 10−b.

Fig. 8 Fig. 8 ~ Continuum-subtracted infrared line emission spectra for a disc inclination of 45° and spectral resolution R = 3000 with varying mixing parameter αmix. It is important to note the different scaling of the y-axis. Most visible spectral features are overlapping o-H2O (blue) and p-H2O (black) emission lines. A few other features emitted from other molecules and ions are annotated in one of the plots. Each molecular spectrum is first convolved with a Gaussian R = 3000 and then overplotted with a certain colour, but not co-added.

Fig. 9 Chemical pathways to create C2H2 and other hydrocarbon molecules without mixing (left) and with αmix = 10−3 (right) at one selected point in the disc at r = 0.15 au (see text). It is important to note the vastly different concentrations of all the hydrocarbon molecules. As in Fig. 7, the green numbers represent the particle concentrations and the blue numbers represent the rates [cm−3 s−1] in the time-independent case. The big red ‘IN’ and ‘OUT’ arrows denote the in-mixing and out-mixing rates. Wiggled arrows indicate photo-dissociation. The light green wiggled arrows are X-ray induced processes. Three reactions appear twice with an added dashed line, when a reaction is important for two reactants or products (C2H2 + C → C3 + H2, C2H2 + C → C3H + H and C3 + hv → C2 + C). On the left side, the reaction paths for C2O, C++, CH5+ and H3CS+ are not followed up.

5 Impact on continuum and line observations

In this section, we demonstrate how turbulent mixing can affect continuum and line observations. The impact is generally quite substantial, and we have selected to demonstrate three of the major effects in this paper. We do not intend to make any predictions for particular targets in this paper.

5.1 Impact on spectral ice features

In our models, the up-mixing of icy grains leads to a substantial increase in the ice concentrations in the observable disc surface (see details in Sect. 4.2.2). Figure 6 shows the spatial distribution of five selected ice species (H₂O#, NH₃#, CO₂#, CO#, and CH₄ #) in the entire disc. In the un-mixed case, most ices only populate the midplane regions below A_V,ver ≈ 10, due to active UV photo-desorption in the upper disc regions. The only counter-example is CO# in the outer regions ≳100 au, but this region is too cold to cause any significant far-IR ice emission features. Since the SED is mostly produced by thermal emission and scattering of the bare (i.e. ice-less) dust grains above A_V,ver = 10, the resulting SED is almost entirely free of ice features, also at grazing inclination angle (≈ 72° in this model), which tests the A_V,rad ≈ 1 regions in absorption. The ice just sits too deep to cause any notable changes in the integrated spectral flux. One would need to observationally search for the ice features by analysing image data, using the intensities coming from the dark lane (see discussion in Arabhavi et al. 2022).

However, even a weak turbulent mixing with α_mix = 10⁻⁵ lifts some of the icy grains upwards, almost reaching the A_V,ver = 1 surface, and some weak ice emission features between 20 µm and 200 µm start to appear in the spectrum, which become more clearly observable in the α_mix = 10⁻⁴ model. We also see that the A_V,ver = 1 contour moves upwards with increasing α_mix due to the opacity of the icy grains. So, scattered starlight can now directly interact with the icy grains.

At α_mix = 10⁻³, the ice emission features become even stronger, and the icy grains reach the A_V,rad = 1 line (i.e. direct starlight can now interact with icy grains). The first three ice absorption features appear when the disc is observed under an inclination angle of ≈ 74°.

Eventually, the extreme model with α_mix = 10⁻² shows a population of H₂O # and NH₃ # up to the very top of the model, accompanied by another considerable lift of the A_V,ver = 1 and A_V,rad = 1 contours. We note that the vertical extension of the snowy grains coincides with the dust temperature T_dust = 150 K contour, which means that the snow extends right up to the water ice sublimation temperature in this model, whereas the UV photo-desorption becomes an irrelevant process (see explanations in Sect. 4.2.2). The same is true for the vertical extension of CO # at larger radial distances, which is limited by the T_dust = 25 K contour. The SED is full of strong ice emission features, and clear absorption features for inclination angles ≳ 72°. We note that the strengths of the ice features depends on how we assume the ices to be distributed onto the pre-existing bare grains, and our assumption of a layer of equal thickness on top of each refractory grain creates the strongest ice opacities (see discussion in Arabhavi et al. 2022).

5.2 Impact on mid-IR gas emission line spectra

Figure 10 shows the dust and gas temperatures, and the concentrations of the IR-active molecules OH, H₂O, CO₂, HCN, and C₂H₂ along a vertical cut at r = 1 au. We see that all molecules shift their location of maximum concentration to the left (i.e. upwards towards smaller column densities), but whereas OH, H₂O, HCN and C₂H₂ are enhanced by mixing, CO₂ is reduced. This is because of the upmixing of CH4, which leads to a rich hydrocarbon chemistry as discussed in Sect. 4.2.3, which suppresses the formation of CO₂. The figure also shows that the gas in the upper layers cools down via increased line emission from OH and H₂O, the concentrations of which are enhanced by mixing.

Figure 8 shows a series of mid-IR gas emission line spectra between 10 µm and 20 µm computed from our series of T Tauri disc models with α_mix = 0 to 10⁻² (see Sect. 3). These spectra were calculated with the Fast Line Tracer (FLITS) developed by Michiel Min (see Woitke et al. 2018) based on the gas and dust temperature structure, the continuum dust opacities, the molecular abundances, and the non-local thermodynamic equilibrium excitation populations calculated by PRODIMO as described in the previous sections. Due to our lack of sufficient collisional data, the populational numbers of the molecules OH, CO₂, HCN, C₂H₂, and NH₃ are assumed to be in local thermodynamic equilibrium. FLITS calculates the entire spectrum on a wavelength grid that corresponds to a velocity resolution of 1 km s⁻¹. We calculated the spectra separately for each molecule listed in the legend of Fig. 8. in contrast, the IR lines of the atoms and ions Ne⁺, Ne⁺⁺, Ar⁺, Ar⁺⁺, Mg⁺, Fe⁺, Si⁺, S, S⁺ and S⁺⁺ were calculated in one go. For more details and explanations, we direct the reader to the work by Woitke et al. (2018). The resulting continuum-subtracted spectra are then convolved with a R = 3000 Gaussian before plotting, to simulate the spectra that we can expect from JWST/MIRI.

As discussed in Woitke et al. (2018), these predicted spectra are broadly consistent with Spitzer/IRS observational data of T Tauri stars concerning line colours and ratios. However, since we assumed a dust/gas ratio of 100 here, instead of 1000 in this paper, and since we only applied a very mild dust settling with α_settle = 10⁻² according to Dubrulle et al. (1995), the dust is still covering most of the line emitting regions, and the predicted IR lines are generally too weak, by a factor of about 10–100. It also means that the lines that are emitted from higher disc layers (e.g. by Ne⁺ and OH) are stronger than those emitted from lower layers (e.g. H₂O, HCN, and C₂H₂). Therefore, we can only discuss the trends here that are caused by the introduction of mixing.

We see little change between α_settle = 0 and α_settle = 10⁻⁴. Only at α_settle = 10⁻³ do the OH lines start to become significantly stronger, and the C₂H₂ and HCN emission features at 13.7 µm and 14 µm, respectively, start to be enhanced and become about comparable to the CO₂ feature at 15 µm, as is observed for many T Tauri stars (Pontoppidan et al. 2010). For α_settle = 10⁻², the OH lines become too strong with respect to the H₂O lines, similar to the findings of Heinzeller et al. (2011). C₂H₂ 13.7 µm is now the strongest spectral feature between 10 µm and 20 µm, after the various OH ro-vibrational features, which is also in disagreement with Spitzer observations.

Fig. 10 Gas and dust temperatures (upper panel), and concentrations of selected IR active molecules (lower panels) along a vertical cut at r = 1 au as a function of vertical hydrogen nuclei column density N〈H〉. The results of five models are overplotted using different line styles, from αmix = 0 to 10−2, see top panel.

5.3 Impact on CO lines and channel maps

The upper row of Fig. 11 shows the resulting CO concentration on linear radial and vertical scales in the same model series. In our disc models, carbon is shown to be mostly present in the form of CO, with amaximum concentration of n_CO/n_〈H〉 ≈ є_C = 10^−3.85, in a vertical layer that is limited by photo-dissociation on top (ionisation parameter χ/n_〈H〉 ≈ 10^−4.5) and by ice formation at the bottom (T_dust ≈ 23 K). The ¹³CO J = 2 → 1 line at 1360.2 µm is mostly emitted from this vertical layer, and 70% of the line flux is emitted radially from about 50 au to 200 au in the disc model with α_mix = 0. Beyond about 250 au, the line becomes optically thin. Figure 11 shows the resulting line fluxes and double-peaked velocity profiles, the velocity-integrated line maps, and selected channel maps in the lowest panel.

We see two major effects of the mixing on the spectral appearance of the CO in millimetre lines. First, the CO line emission becomes weaker with increasing α_mix at small radii. Figure 4 shows that at 100 au (right column), the disc cools down slightly with increasing α_mix. This is because of the additional presence of icy grains in higher disc layers due to mixing (see Sect. 4.2.2.) Because of their enhanced opacity, the icy grains increase the vertical optical depths and also scatter away the starlight more efficiently. As a consequence, the height below which CO freezes out increases, and hence the vertical CO column density decreases. Figure 11 shows a clear trend of the CO emissions from the inner disc regions becoming weaker with increasing α_mix (see line maps), while the vertical layer of maximum CO concentration becomes thinner (see concentrations).

Second, the CO line emissions intensify with increasing α_mix at large radii. This is a direct consequence of downwards CO mixing. For a large α_mix, the CO mixing timescale becomes shorter than the ice adsorption timescale in the outer regions, which leads to a more gradual decline of the CO concentration towards the midplane, and hence larger CO column densities in the outermost disc regions. The lower boundary of the vertical layer of maximum CO concentration becomes less well defined for a large α_mix. This leads to a more ‘diffuse’ spectral appearance in millimetre CO channel maps (see lower row of plots in Fig. 11). Without mixing, we obtain clearly separated line emissions from two lobes: one from the upper disc half, and one from the lower disc half. With strong mixing, α_mix = 10⁻², the two lobes visually merge beyond about 150 au. Cleeves et al. (2016) and Pinte et al. (2018) noticed that for the case of the giant disc of IM Lupi, where the CO lines are detectable out to ~1000 au, it becomes impossible to identify two separate CO layers beyond ~450 au. Pinte et al. (2018) achieved a distinctively better fit when adding diffuse CO emission from the outer midplane. We note, however, that even without mixing, the CO concentration in the midplane increases with radius as the freeze-out timescale of CO increases with falling density, whereas the photo-dissociation timescale is density independent.

In the model with the strongest mixing, α_mix = 10⁻², we furthermore observe that the CO is also vertically more extended towards the top. This is partly caused by an upwards shift of the χ/n = 10^−4.5 contour line due to the increased dust opacities of the icy grains, but at the very top, the CO is actually mixed up quicker than it can be photo-dissociated. In summary, we confirm the finding of Semenov et al. (2006) that the CO density can be enhanced by mixing above and below the layer of maximum concentration, but we disagree with Semenov & Wiebe (2011) that the CO is not responsive to mixing.

Fig. 11 Influence of vertical mixing on millimetre CO observations. The top row shows the CO concentration є(CO) = nCO/n〈H〉 in the disc on linear scales for radius r and vertical height z. The two dashed contour lines enclose the region of maximum CO concentration, where χ is the UV strength and n = n〈H〉 is the total hydrogen nuclei particle density. Below dust temperature Tdust ≈ 23 K, CO is mostly frozen. The second row shows the 13CO J= 2 → 1 line velocity profile observed from a distance of 140 pc at a disc inclination of 45°. The third row shows the integrated line map after continuum subtraction and the bottom row shows a channel map at radial velocity v + 0.94 km s−1. To plot the intensities in Kelvin, we consider Iv c2/(2v2)/k, where Iv [erg cm−2 s−1 Hz−1 sr−1] is the spectral intensity, c the speed of light, v the frequency, and k the Boltzmann constant.

Fig. 12 Schematic sketch of the disc regions responsible for the formation of atoms, molecules, and ices, and their vertical diffusive transport.

6 Summary and conclusions

6.1 Summary of chemical mechanisms

Figure 12 summarises the principle mechanisms that we have identified in our disc models by means of a physical sketch. The mixing is immediately relevant for the most abundant species, such as H2, CH4, CO, the neutral atoms in higher layers, and the ices in the midplane. By ‘immediately relevant’, we mean that the transport rates can locally dominate over the chemical formation and destruction rates. The abundant species, whose concentrations are limited by element abundance constraints, are transported away from their regions of maximum concentration, both up and down, and in the locations where these abundant chemicals are finally destroyed, for example by photo-reactions, their reaction products fuel the local chemistry. The release of these fragments has important consequences for all other, chemically interlinked but less abundant molecules, which quickly adjust their concentrations to these changed circumstances. This generally creates a more active chemistry, with a richer mixture of ionised, atomic, molecular, and ice species, and new chemical pathways that are not relevant in the un-mixed case.

This situation differs somewhat from exoplanet atmospheres as discussed, for example, by Moses et al. (2011), Zahnle et al. (2016), and Barth et al. (2021). In exoplanet atmospheres, the chemical factories are located in the deepest layers, where temperatures and densities are high. From there, molecules are mixed upwards, resulting in constant concentrations at intermediate heights, until photo-dissociation destroys them at high altitudes. In discs, the main chemical factories are located in the warm intermediate layers. From there, molecules are mixed upwards and downwards. The dense midplane is too cold to have an active chemistry. Instead, molecules freeze out there, and the icy grains are mixed up.

6.2 Summary of impact on spectral observations

Concerning the mid-IR molecular spectra, we see increased molecular abundances in the upper layers observable with Spitzer and JWST, in particular for OH, H₂O, HCN, and C₂H₂ in our models when we include mixing. The chemical response to mixing in these layers happens on short timescales, of the order of hours or days, and liberates more heat via exothermic reactions. However, at the same time the increased abundances of the IR active molecules lead to increased line cooling and a decrease in gas temperatures in those high layers. Despite the strong effects on both molecular abundances and the gas temperature structure, the calculated IR spectra are surprisingly robust. We need to assume strong mixing, α_mix ≳ 10⁻³, to produce notable spectral changes. The spectral features of OH, HCN, and C₂H₂ are enhanced by mixing, whereas the CO₂ feature at 15 µm is reduced.

Concerning the CO millimetre emission lines, mixing reduces the line intensities at small radii but enhances them at larger radii, leading to an overall more flat line intensity profile in ALMA line observations. For intense mixing, CO starts to populate the outer midplane regions, which leads to a distinct, more diffuse appearance of the CO in velocity channel maps. However, it is likely that such differences in CO channel maps can also be produced in alternative ways in our models, for example by changing the column density profile N_〈H〉(r) of the disc.

6.3 Conclusions about the α viscosity parameter

The parameter α_mix describes the strength of the vertical turbulence according to the Shakura & Syunyaev (1973) α-model. Our mixed models with α_mix ≳ 10⁻³ show strong ice emission features between 20 µm and 200 µm, which neither Spitzer nor Herschel have seen for T Tauri stars. However, the α_mix ≈ 10⁻⁵ and 10⁻⁴ models predict weak ice features that might require the signal-to-noise ratio of JWST to be detectable. Without mixing, it would be difficult to explain such detections with our models.

Concerning the mid-IR molecular emission lines, our models with α_mix ≈ 10⁻³ predict more equal strengths of the CO₂, HCN, and C₂H₂ features between 13.7 µm and 15 µm, which is difficult to achieve in our models otherwise. However, for α_mix ≈ 10⁻², the OH and C₂H₂ emission features become unrealistically strong.

The detection of layered CO emission in ALMA channel maps, in the form of two distinct lobes emitted by the lower and the upper disc halves, with no direct emission from the in-between midplane, allows us to put some constraint on the α viscosity parameter, because for a large α_mix, we would expect CO to populate the midplane. However, that constraint is weak, about α_mix < 10⁻².

We conclude that models with a modest choice of α_mix ≈ 10⁻⁵ to 10⁻³ might work best to produce realistic continuum and line observations. However, a constant value of α_mix is certainly unrealistic in the first place, as the turbulence is expected to vary with distance and height (see e.g. Fig. 9 of Thi et al. 2019).

Acknowledgements

P.W. acknowledges funding from the European Union H2020-MSCA-ITN-2019 under Grant Agreement no. 860470 (CHAMELEON).

Appendix A Numerical implementation of mixing

Appendix A.1 Numerical scheme

We consider height grid points {z_k | k = 1,…, K}. The basic reaction-diffusion equation (Eq. 1) is numerically discretised as

$$\frac{\partial n_{i,k}}{\partial t}=P_{i,k}-L_{i,k}-\frac{{\text{Φ}}_{i,k+\mathrm{½}}-{\text{Φ}}_{i,k-\mathrm{½}}}{z_{k+\mathrm{½}}-{\text{z}}_{k-\mathrm{½}}},$$(A.1)

where the notation means $x_{k\pm \mathrm{½}}=\frac{1}{2}\left(x_k+x_{k\pm 1}\right)$, and

$$\begin{array}{lll}{\text{Φ}}_{i,k+\mathrm{½}}\hfill & =\hfill & -\left(K_{k+\mathrm{½}}+K_{i,k+\mathrm{½}}\right)n_{k+\mathrm{½}}\frac{\frac{n_{i,k+1}}{n_{k+1}}-\frac{n_{i,k}}{n_k}}{z_{k+1}-z_k}\hfill \\ \hfill & \hfill & +D_{i,k+\mathrm{½}}{n}_{i,k+\mathrm{½}}\frac{\left({\text{μ}}_{k+\mathrm{½}}-m_i\right)g_{k+\mathrm{½}}}{k{T}_{k+\mathrm{½}}},\hfill \end{array}$$(A.2)

$$\begin{array}{lll}{\text{Φ}}_{i,k-\mathrm{½}}\hfill & =\hfill & -\left(K_{k-\mathrm{½}}+K_{i,k-\mathrm{½}}\right)n_{k-\mathrm{½}}\frac{\frac{n_{i,k}}{n_k}-\frac{n_{i,k-1}}{n_{k-1}}}{z_k-z_{k-1}}\hfill \\ \hfill & \hfill & +D_{i,k-\mathrm{½}}{n}_{i,k-\mathrm{½}}\frac{\left({\text{μ}}_{k-\mathrm{½}}-m_i\right)g_{k-\mathrm{½}}}{k{T}_{k-\mathrm{½}}}.\hfill \end{array}$$(A.3)

Assuming that all quantities are given and constant except for n_i,k−1, n_i,k and n_i,k+1 we write

$$-\frac{{\text{Φ}}_{i,k+\mathrm{½}}-{\text{Φ}}_{i,k-\mathrm{½}}}{z_{k+\mathrm{½}}-z_{k-\mathrm{½}}}=A_{i,k}{n}_{i,k+1}-B_{i,k}{n}_{i,k}+C_{i,k}{n}_{i,k-1},$$(A.4)

with

$$\begin{array}{lll}{A}_{i,k}\hfill & =\hfill & \frac{1}{z_{k+\mathrm{½}}-z_{k-\mathrm{½}}}(+\left(K_{k+\mathrm{½}}+D_{i,k+\mathrm{½}}\right)\frac{n_{k+\mathrm{½}}}{n_{k+1}\left(z_{k+1}-z_k\right)}\hfill \\ \hfill & \hfill & -D_{i,k+\mathrm{½}}\frac{1}{2}\frac{\left({\text{μ}}_{k+\mathrm{½}}-m_i\right)g_{k+\mathrm{½}}}{k{T}_{k+\mathrm{½}}}),\hfill \end{array}$$(A.5)

$$\begin{array}{lll}{B}_{i,k}\hfill & =\hfill & \frac{1}{z_{k+\mathrm{½}}-z_{k-\mathrm{½}}}(-\left(K_{k+\mathrm{½}}+D_{i,k+\mathrm{½}}\right)\frac{n_{k+\mathrm{½}}}{n_k\left(z_{k+1}-z_k\right)}\hfill \\ \hfill & \hfill & +D_{i,k+\mathrm{½}}\frac{1}{2}\frac{\left({\text{μ}}_{k\text{+½}}-m_i\right)g_{k+\mathrm{½}}}{k{T}_{k+\mathrm{½}}}\hfill \\ \hfill & \hfill & +\left(K_{k-\mathrm{½}}-D_{i,k+\mathrm{½}}\right)\frac{n_{k-\mathrm{½}}}{n_k\left(z_k-z_{k-1}\right)}\hfill \\ \hfill & \hfill & -D_{i,k-\mathrm{½}}\frac{1}{2}\frac{\left({\text{μ}}_{k-\mathrm{½}}-m_i\right)g_{k-\mathrm{½}}}{k{T}_{k-\mathrm{½}}}),\hfill \end{array}$$(A.6)

$$\begin{array}{lll}{C}_{i,k}\hfill & =\hfill & \frac{1}{z_{k+\mathrm{½}}-z_{k-\mathrm{½}}}(+\left(K_{k-\mathrm{½}}+D_{i,k-\mathrm{½}}\right)\frac{n_{k-\mathrm{½}}}{n_{k-1}\left(z_k-z_{k-1}\right)}\hfill \\ \hfill & \hfill & +D_{i,k-\mathrm{½}}\frac{1}{2}\frac{\left({\text{μ}}_{k-\mathrm{½}}-m_i\right)g_{k-\mathrm{½}}}{k{T}_{k-\mathrm{½}}}).\hfill \end{array}$$(A.7)

The apparent ± asymmetries in the (K + D) and D terms are correct; they reflect the different nature of mixing, which is modelled as divergence of a gradient in concentration (second order), and gravitational de-mixing, which is modelled as divergence of a mass-dependent flux (first order).

Fig. A.1 Test of the iterative numerical method to solve the diffusion equation without chemical sources and sinks. We consider a δ-peak for the initial concentration xmol = ni,k/nk at the start of the iterations in a Gaussian density structure (see text). The upper figure shows the concentrations as a function of the number of iterations carried out. The lower figure shows the relative deviations between iterations (see Eqs. (28) and (29)).

Appendix A.2 Boundary condition

We apply zero-flux boundary conditions at k = 1 and k = K:

$${\text{Φ}}_{i,\mathrm{½}}={\text{Φ}}_{i,K+\mathrm{½}}=0.$$(A.8)

In practice, this means that all terms in Eqs. (A.5) to (A.7) with K_k−½ or D_i,k−½ are dropped at k = 1, and all terms with K_k+½ or D_i,k+½ are dropped at k = K. For the pre-factor, we approximate z_k+½ − z_k-½ ≈ z₂ − z₁ at k = 1 and z_k+½ - z_k−½ ≈ z_K − z_K−1 at k = K.

Fig. A.1 – Test of the iterative numerical method to solve the diffusion equation without chemical sources and sinks. We consider a δ-peak for the initial concentration x_mol = n_i,k/n_k at the start of the iterations in a Gaussian density structure (see text). The upper figure shows the concentrations as a function of the number of iterations carried out. The lower figure shows the relative deviations between iterations (see Eqs. (28) and (29)).

Appendix A.3 Initial condition

The very first disc iteration is computed with zero mixing rates. In all forthcoming iterations, the mixing rates are activated.

Appendix A.4 Test of the numerical solution method

To verify our numerical implementation, we made a quick test to check the capabilities of our simple iterative scheme. We considered a Gaussian density profile on an equidistant spatial grid with 101 points k = 1…101

$$z_k=k-1,$$(A.9)

$$n_k=\exp \left(-{\left(\frac{z_k}{20}\right)}^2\right),$$(A.10)

$$K_k=1,$$(A.11)

$$D_{i,k}=0,$$(A.12)

and chose, for a single test species i = 1, the following initial condition:

$$n_{i,k}=\{\begin{array}{ll}1\hfill & ,k=50\hfill \\ 0\hfill & ,\text{otherwise}\hfill \end{array}$$(A.13)

at the start of the iteration. Then we copied the exact implementation of the vertical mixing in PRODIMO to an external tool

where we omitted the chemical source and sink terms P_i,k and L_i,k. In a loop from k = K down to 1, we then

• 1)

  calculated A_i,k, B_i,k and C_i,k, and

• 2)

  solved the rate equation without chemical sources and sinks $n_{i,k}=\frac{1}{B_{i,k}}\left(A_{i,k}{n}_{i,k+1}+C_{i,k}{n}_{i,k-1}\right)$.

This numerical scheme mimics the downwards sweep in PRODIMO while solving the chemistry in the time-independent case from top to bottom. Figure A.1 shows the resulting concentrations n_i,k/n_k as a function of z, where z = 0 represents the midplane. The results after selected numbers of iterations are plotted with different colours as indicated. This plot can be compared to our results with chemical sources and sinks plotted in Fig. 2 in the main part of the paper.

The figure indicates that up to 5000 iterations are required to reach a constant concentration, as is physically expected in this setup. The final value of that constant concentration is 0.61%, which is close to the predicted value of 0.63% that we get from a numerical integral of n_ik(z) over z in the initial state, divided by the numerical integral of n(z). We note that in real applications that include the chemical source terms, the chemicals will rarely travel many scale heights. Instead, there are normally some dominant but spatially restricted formation and destruction regions, and the model must properly describe the diffusive transport between those regions. Hence, we expect a couple of hundred iterations to be sufficient.

Figure A.1 demonstrates the numerical behaviour during the course of the iterations:

• The first few tens of iterations are featured by an asymmetry, because the resultant n_i,k+1 is already used for the determination of n_i,k, whereas n_i,k−1 is not. Therefore, the information of a certain molecule being mixed downwards is travelling faster downwards than it travels upwards.

• After some hundreds of iterations, however, the upper regions fill in quicker than the lower regions. This is because of the density gradient in the system; a large reservoir in the denser lower disc regions is able to more quickly fill in the low-density regions above than vice versa.

References

All Tables

All Figures

Fig. 1 Location of the homopause, where the diffusion coefficients caused by eddy mixing and gas-kinetic motions are equal (K=Di). The colours show the hydrogen nuclei density, n〈H〉, as a function of radius r and relative height z/r. The coloured contour lines show where the homopause is situated when K is calculated with different values of log10 αmix. The dashed lines show where the vertical optical extinction AV,ver equals 100 (black), ten (dark grey) and one (light grey). The dust opacities correspond to our un-mixed model described in Sect. 3. The uneven AV contours are caused by the uneven spatial distribution of the ices, which are taken into account when calculating the opacities.

In the text

	Fig. 2 Measuring the convergence of the iterative solution technique to solve the chemistry with vertical mixing, here for the model with medium strong mixing, αmix = 10−3. See text for explanations.
In the text

Fig. 3 Monitoring the convergence of the iterative solution method along a vertical cut at r = 10 au as a function of the vertical hydrogen nuclei column density $N_{\langle \text{H}\rangle }={\displaystyle {\int }_z^{\infty }{n}_{\langle \text{H}\rangle }\left(z\right)\text{\hspace{0.17em}d}z}$ for the model with medium strong mixing, αmix = 10−3. The upper plots show the changing gas and dust temperature as a function of iteration number. The lower plots shows the concentrations of a few selected oxygen-bearing gas and ice species. The left row shows the first few iterations and the right row shows the convergence towards iteration 1000, which is considered as the final solution. The different full, dashed and dotted lines correspond to the different numbers of iterations carried out as indicated.

In the text

Fig. 4 Resulting temperatures and chemical concentrations as a function of vertical hydrogen column density N〈H〉, for radial distances 1 au, 10 au, and 100 au, from left to right. The top row shows the gas and dust temperatures, the second row the main hydrogen species, the third row the oxygen chemistry, and the last row the carbon chemistry. Only a few selected species are plotted. In each part of the figure, the various dashed and dotted lines correspond to the mixing parameter αmix, as explained in the legends.

In the text

Fig. 5 Vertical spreading of icy grains at r = 100 au for different mixing strengths αmix. Results are shown as a function of the vertical hydrogen nuclei column density N〈H〉. The top row of plots shows the computed gas and dust temperatures, Tg and Td. The middle row of plots shows the bare grain surface area per H nucleus = 4π〈a2〉nd/n〈H〉 after dust settling in blue, and the number of ice layers on each grain Nlayers in red. The lower plots show the concentrations of selected ice species.

In the text

Fig. 6 Upper row of plots shows the location of the ices, in a series of models with increasing mixing parameter αmix. The icy regions are hashed when the respective ice concentration, for example nH2O#/n〈H〉 for water ice, is larger than 10−55. The background colours show the hydrogen nuclei density n〈H〉 from 102 cm−3 (black) to 1016 cm−3 (red). The dashed black contour lines show the radial visual extinction AV,rad = 1 and vertical visual extinction AV,ver = 1, and the red dashed lines show two dust temperature contours, 25 K and 150 K. The lower row of plots shows the respective SEDs for selected inclination angles at a distance of 140 pc. The yellow lines indicate the adopted photospheric spectrum of the naked star. We note that the absorption feature around 2.5 µm is due to molecular absorption in the stellar photosphere, whereas all other narrow absorption and emission features are due to ices, except for the broad 10 µm/20 µm silicate emission features.

In the text

Fig. 7 Chemical pathways to create OH and H2O molecules without mixing (left) and with αmix = 10−3 (right) at one selected point in the disc at r = 10 au located in the upper disc layers (see text). The green numbers represent the particle concentrations and the blue numbers represent the rates [cm−3 s−1] in the time-independent case. The big red ‘IN’ and ‘OUT’ arrows on the right denote the in-mixing and out-mixing rates for H2. Wiggled arrows indicate photo-dissociation. Light green arrows on the left mark the H2 formation on grains, and the radiative association reactions O + H → OH + hv and OH + H → H2O + hv. ${\text{H}}_2^{★}$ denotes electronically excited molecular hydrogen. The figure does not follow up the production of CH2+, CO+, and HCO+ marked in blue. The notation a(−b) means a × 10−b.

In the text

Fig. 8 Fig. 8 ~ Continuum-subtracted infrared line emission spectra for a disc inclination of 45° and spectral resolution R = 3000 with varying mixing parameter αmix. It is important to note the different scaling of the y-axis. Most visible spectral features are overlapping o-H2O (blue) and p-H2O (black) emission lines. A few other features emitted from other molecules and ions are annotated in one of the plots. Each molecular spectrum is first convolved with a Gaussian R = 3000 and then overplotted with a certain colour, but not co-added.

In the text

Fig. 9 Chemical pathways to create C2H2 and other hydrocarbon molecules without mixing (left) and with αmix = 10−3 (right) at one selected point in the disc at r = 0.15 au (see text). It is important to note the vastly different concentrations of all the hydrocarbon molecules. As in Fig. 7, the green numbers represent the particle concentrations and the blue numbers represent the rates [cm−3 s−1] in the time-independent case. The big red ‘IN’ and ‘OUT’ arrows denote the in-mixing and out-mixing rates. Wiggled arrows indicate photo-dissociation. The light green wiggled arrows are X-ray induced processes. Three reactions appear twice with an added dashed line, when a reaction is important for two reactants or products (C2H2 + C → C3 + H2, C2H2 + C → C3H + H and C3 + hv → C2 + C). On the left side, the reaction paths for C2O, C++, CH5+ and H3CS+ are not followed up.

In the text

	Fig. 10 Gas and dust temperatures (upper panel), and concentrations of selected IR active molecules (lower panels) along a vertical cut at r = 1 au as a function of vertical hydrogen nuclei column density N〈H〉. The results of five models are overplotted using different line styles, from αmix = 0 to 10−2, see top panel.
In the text

Fig. 11 Influence of vertical mixing on millimetre CO observations. The top row shows the CO concentration є(CO) = nCO/n〈H〉 in the disc on linear scales for radius r and vertical height z. The two dashed contour lines enclose the region of maximum CO concentration, where χ is the UV strength and n = n〈H〉 is the total hydrogen nuclei particle density. Below dust temperature Tdust ≈ 23 K, CO is mostly frozen. The second row shows the 13CO J= 2 → 1 line velocity profile observed from a distance of 140 pc at a disc inclination of 45°. The third row shows the integrated line map after continuum subtraction and the bottom row shows a channel map at radial velocity v + 0.94 km s−1. To plot the intensities in Kelvin, we consider Iv c2/(2v2)/k, where Iv [erg cm−2 s−1 Hz−1 sr−1] is the spectral intensity, c the speed of light, v the frequency, and k the Boltzmann constant.

In the text

	Fig. 12 Schematic sketch of the disc regions responsible for the formation of atoms, molecules, and ices, and their vertical diffusive transport.
In the text

Fig. A.1 Test of the iterative numerical method to solve the diffusion equation without chemical sources and sinks. We consider a δ-peak for the initial concentration xmol = ni,k/nk at the start of the iterations in a Gaussian density structure (see text). The upper figure shows the concentrations as a function of the number of iterations carried out. The lower figure shows the relative deviations between iterations (see Eqs. (28) and (29)).

In the text