Mixing and diffusion in protoplanetary disc chemistry

We develop a simple iterative scheme to include vertical turbulent mixing and diffusion in ProDiMo thermo-chemical models for protoplanetary discs. The models are carefully checked for convergence toward the time-independent solution of the reaction-diffusion equations, as e.g. used in exoplanet atmosphere models. A series of five T Tauri disc models is presented where we vary the mixing parameter {\alpha} mix from 0 to 0.01 and take into account (a) the radiative transfer feedback of the opacities of icy grains that are mixed upward and (b) the feedback of the changing molecular abundances on the gas temperature structure caused by exothermic reactions and increased line heating/cooling. We see considerable changes of the molecular and ice concentrations in the disc. The most abundant species (H2, CH4, CO, the neutral atoms in higher layers, and the ices in the midplane) are transported both up and down, and at the locations where these abundant chemicals finally decompose, for example by photo processes, the release of reaction products has important consequences for all other molecules. 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 unmixed case. We discuss the impact on three spectral observations caused by mixing and find that (i) icy grains can reach the observable disc surface where they cause ice absorption and emission features at IR to far-IR wavelengths, (ii) mixing increases the concentrations of certain neutral molecules observable by mid-IR spectroscopy, in particular OH, HCN and C2H2, and (iii) mixing can change the optical appearance of CO in ALMA line images and channel maps, where strong mixing would cause the CO molecules to populate the distant midplane.


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 gaskinetic 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 2 , 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 −2 .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 DM Tau with given densities A164, page 1 of 19 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe-to-Open model.Subscribe to A&A to support open access publication.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 2 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 T Tauri 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 −2 .
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.

Basic equations
The basic rate equation to solve for each chemical species is where n i [cm −3 ] is the particle density of species i, P i , and L i [cm −3 s −1 ] are the chemical production and destruction rates, and Φ i [cm −2 s −1 ] is the diffusive particle flux according to Hu et al. (2012), Zahnle et al. (2016), andRimmer &Helling (2016): where n = i n i is the total particle density.m i is the mass of gas species i and µ = n i m i /n the mean molecular weight.g z = (GM ⋆ z)/(r 2 + z 2 ) 3/2 = z Ω 2 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 2 s −1 ] 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: 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 α-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 : where c s = √ kT/µ is the local isothermal sound speed and T the temperature.Typical values used in the literature are α ≈ 10 −5 ... 10 −2 using different constraints (see reviews by Lesur et al. 2022 andPinte 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 2 s −1 ], due to Brownian motions, 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) where the effective radii of the molecules r i can be measured from viscosity experiments (Jeans 1967), or calculated from

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 timeindependent functions of z.The diffusive fluxes vanish and we have with molecular and atmospheric scale heights H i = kT/(m i g 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 2 , 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.

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 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: Chemical reactions conserve the element abundances, and thus the first term in Eq. ( 8) vanishes.Concerning the diffusive flux Φ i [cm −2 s −1 ], we first consider the case K ≫ D i where turbulent mixing dominates.From Eq. (2) we then have The diffusive element fluxes are hence given by where K does not depend on molecule.The last term in this equations can be expressed by In general, the total particle density n 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 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 (Φ k → 0 ⇒ ∂ϵ k ∂z → 0), 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), andYoudin &Lithwick (2007).The result of Schräpler & Henning (2004, see their Eq.( 27)), is where St is the Stokes number given by where τ stop is the stopping distance, τ eddy the turnover timescale of the largest turbulent eddy, ρ d the dust material density, v th the thermal velocity, c s the sound velocity, H the pressure scale height, and a the grain radius.The factor √ π/8 is because c 2 s = kT/µ but v 2 th = 8kT/(πµ).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 a 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): Defining the element abundances present in the gas and ice, respectively, as we find the following result that is independent of α mix :

Small Stokes numbers
In the limiting case of small dust radii a 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 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 2 , 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.

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 Equation ( 20) states a remarkable result.If an element can freeze out anywhere in the column, a rather low level of ϵ gas k will be enforced there, and this low level of ϵ gas k must be identical everywhere in the column to make the diffusive fluxes vanish, irrespective of the local ϵ ice k .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 ϵ gas k .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 −2 , such that the mixing timescale according to Eq. ( 32) approaches the disc lifetime.For an assumed disc lifetime of 5 × 10 6 yr, this would suggest α mix < 3 × 10 −5 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 ϵ gas k (z) and ϵ ice k (z).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 µ.

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 andRimmer &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 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 with the existing time-independent chemistry solver in the PRODIMO model, where 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 : 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 for each chemical species i on all grid points k.

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 A164, page 5 of 19 A&A 668, A164 (2022) 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 α settle = 10 −2 .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 2 O-ice, CO 2 -ice, NH 3 -ice, CH 4 -ice, CO-ice, and CH 3 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.

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: 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 −3 .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 H 2 O 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. 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 2 O cause reinforced line cooling at column densities between about 10 19 and 10 22 cm −2 , 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 2 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.

Beating UV photo-desorption by mixing
At large radial distances, ∼100 au, the model with the highest mixing strength, α mix = 10 −2 , exhibits some rather peculiar results (see Fig. 5).The concentrations of the thermally most stable ice species such as H 2 O#, NH 3 #, CO 2 #, OH#, and C 2 H 4 # (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 2 # (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 −2 , 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 2 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 10 photons cm −2 s −1 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 8 photons cm −2 s −1 (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 where n d is the number density of dust particles, π⟨a2 ⟩n d the total dust cross section per cm 3 , Y i the UV photo-desorption yield, and n act = 4π⟨a 2 ⟩ n d 2 n surf the total number of active adsorption sites on all grains per cm 3 .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 15 sites per cm 2 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 ice tot expresses the probability of finding ice species i in the ice mixture present, where the total ice density is n ice tot = 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 layers = n ice tot /(4π⟨a 2 ⟩n d n surf ) ≈ 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 where we have assumed a total initial ice concentration of n ice tot /n ⟨H⟩ ≈ 10 −3.5 , a mean photo-desorption yield of Ȳ ≈ 10 −3 , and a total dust surface per H-atom of 4π⟨a 2 ⟩n d /n ⟨H⟩ ≈ 2 × 10 −23 cm 2 (which can be compared with Fig. 5).We note that χ ≈ 300 is what we find in our model at r = 100 au.χ ∝ 1/r 2 (i.e.τ UVdesorb ∝ r 2 ) is valid until the interstellar UV photons become more relevant than the diluted stellar UV, which happens at about 1700 au in this case.
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 −2 , 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 −2 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.
Figure 7 shows how the up-mixing of H 2 fuels the ion-neutral chemistry in the upper disc regions, resulting in an enhanced production of OH and H 2 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 + hν, water is mainly produced by the followup radiative association reaction OH + H → H 2 O + hν, and H 2 is produced partly by the formation on grains and partly via On the right side of Fig. 7, which shows the results at the same point for the case with mixing α mix = 10 −3 , the H 2 -concentration is enhanced by a factor of about 1000.Here, the starting point for the production of OH and H 2 O is the reaction where  33) to be more efficient in producing OH + , because the mixing leads to a coexistence of H 2 and O + .Since H 2 is required several times along these formation pathways, even a small increase in the H 2 concentration, caused by up-mixing, can lead to a much higher production of OH and H 2 O molecules.The resulting OH and H 2 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 −3 .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 2 is the only molecule that is strongly affected by mixing directly in Fig. 7.All other molecules shown in Fig. 7  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 −3 disc model at r = 0.15 au, at a height where the C 2 H 2 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 −2 on the left plot in Fig. 4. Points such as these are flagged by PRODIMO as the main emission region for the C 2 H 2 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 −3 , 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 rad V = 194, but semioptically thin in the vertical direction, A ver V = 3.1, still resulting in a total UV strength of χ = 20000.The most abundant species are H 2 [0.5],H [5.6(−3)], H 2 O [1.7(−4)], CO [1.3(−4)], and N 2 [3.6(−5)], with no ices.
Figure 9 shows how the up-mixing of CH 4 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 2 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 (CH 4 ) 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 with branching ratios from the UMIST 2012 database (McElroy et al. 2013).This introduces rates of the order of 10 −1 cm −3 s −1 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 3 H.The concentration of acetylene (C 2 H 2 ), 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 2 H 2 and CO 2 ) in similarly high concentrations.
Again, only the abundant molecules shown with big red arrows in Fig. 9 (CO, CO 2 and CH 4 ) are significantly affected by mixing directly.All other molecules just adjust to those circumstances quickly.For example, C 2 H 2 has a chemical relaxation timescale of about n C2H2 /rate ≈ 8(−7) × 1(+13) cm −3 /(3(+4) cm −3 s −1 ) = 270 s at this point, which is much shorter than the mixing timescale ≈ 2 yr (Eq.31).
In summary, the CH 4 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 2 that is mixed out.3.8(-2) Fig. 9.Chemical pathways to create C 2 H 2 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 On the left side, the reaction paths for C 2 O, C ++ , CH 5 + and H 3 CS + are not followed up.

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.

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 2 O #, NH 3 #, CO 2 #, CO #, and CH 4 #) 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 −5 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 −4 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 −3 , 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 −2 shows a population of H 2 O # and NH 3 # 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).

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 2 O, CO 2 , HCN, and C 2 H 2 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 2 O, HCN and C 2 H 2 are enhanced by mixing, CO 2 is reduced.This is because of the upmixing of CH 4 , which leads to a rich hydrocarbon chemistry as discussed in Sect.4.2.3, which suppresses the formation of CO 2 .The figure also shows that the gas in the upper layers cools down via increased line emission from OH and H 2 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 −2 (see Sect. 3).These spectra A164, page 13 of 19 A&A 668, A164 (2022) 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.
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 2 , HCN, C 2 H 2 , and NH 3 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 −1 .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 −2 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 2 O, HCN, and C 2 H 2 ).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 −4 .Only at α settle = 10 −3 do the OH lines start to become significantly stronger, and the C 2 H 2 and HCN emission features at 13.7 µm and 14 µm, respectively, start to be enhanced and become about comparable to the CO 2 feature at 15 µm, as is observed for many T Tauri stars (Pontoppidan et al. 2010).For α settle = 10 −2 , the OH lines become too strong with respect to the H 2 O lines, similar to the findings of Heinzeller et al. (2011).C 2 H 2 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.

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 a maximum 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 13 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   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 −2 , 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.

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 H 2 , CH 4 , 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.

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 2 O, HCN, and C 2 H 2 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 −3 , to produce notable spectral changes.The spectral features of OH, HCN, and C 2 H 2 are enhanced by mixing, whereas the CO 2 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.

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 −3 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 −5 and 10 −4 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 −3 predict more equal strengths of the CO 2 , HCN, and C 2 H 2 features between 13.7 µm and 15 µm, which is difficult to achieve in our models otherwise.However, for α mix ≈ 10 −2 , the OH and C 2 H 2 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 inbetween midplane, allows us to put some constraint on the α viscosity parameter, because for a large α mix , we would expect A164, page 16 of 19 CO to populate the midplane.However, that constraint is weak, about α mix < 10 −2 .
We conclude that models with a modest choice of α mix ≈ 10 −5 to 10 −3 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).
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 = 1 B i,k A i,k n i,k+1 + C i,k n i,k−1 .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.-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.

Fig. 2 .
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 .Fig. 3 .
Figure 4 shows our resulting temperature and chemical structures.Each plot shows a sequence of well-iterated results for α mix ∈ {0, 10 −5 , 10 −4 , 10 −3 , 10 −2 } 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 −2 ].The figure shows a number

Fig. 4 .
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 .
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, T g and T d .The middle row of plots shows the bare grain surface area per H nucleus = 4π⟨a 2 ⟩n d /n ⟨H⟩ after dust settling in blue, and the number of ice layers on each grain N layers in red.The lower plots show the concentrations of selected ice species.

Fig. 6 .
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 n H2O# /n ⟨H⟩ for water ice, is larger than 10 −5.5 .The background colours show the hydrogen nuclei density n ⟨H⟩ from 10 2 cm −3 (black) to 10 16 cm −3 (red).The dashed black contour lines show the radial visual extinction A V,rad = 1 and vertical visual extinction A V,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.A164, page 9 of 19

H2Fig. 7 .
Fig. 7.Chemical pathways to create OH and H 2 O 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 H 2 .Wiggled arrows indicate photo-dissociation.Light green arrows on the left mark the H 2 formation on grains, and the radiative association reactions O + H → OH + hν and OH + H → H 2 O + hν.H ⋆ 2 denotes electronically excited molecular hydrogen.The figure does not follow up the production of CH 2 + , 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-H 2 O (blue) and p-H 2 O (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. 11.Influence of vertical mixing on millimetre CO observations.The top row shows the CO concentration ϵ(CO) = n CO /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 T dust ≈ 23 K, CO is mostly frozen.The second row shows the 13 CO 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 I ν c 2 /(2ν 2 )/k, where I ν [erg cm −2 s −1 Hz −1 sr −1 ] is the spectral intensity, c the speed of light, ν the frequency, and k the Boltzmann constant.

Fig. 11
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 −2 , the two lobes visually merge beyond about 150 au.Cleeves et al. (2016) andPinte 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 −2 , 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 ofSemenov 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.
Figure A.1 demonstrates the numerical behaviour during the course of the iterations:

Table 1 .
Effective molecular radii for diffusion coefficients calculated with the tools provided at Wolfram's webpage 1 at T = 300 K.
Hu et al. (2012)19 P.Woitke et al.:Mixing and diffusion in protoplanetary disc chemistry Particle i Radius r i −1 s −1 ] for collisions between species i and a 'background atmosphere' denoted by a.Some binary diffusion coefficients for a few relevant combinations of H, H 2 , and CO 2 are given inHu et al. (2012).Most atoms and molecules have effective collisional radii r i between about 1 Å and 3 Å, as can be seen in Table1.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).
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 2 (after outmixing, photo-dissociation, and reaction with CH + ), it triggers the production of OH and H 2 O via the formation of OH + .Quick H-addition reactions with H 2 subsequently produce H 2 O + and H 3 O + , before recombinations with free electrons create the neutral molecules OH and H 2 O.This H 2 O formation path is nomally initiated by reactions of O with H 3 + , for example O + H + 3 → OH + + H 2 , but here we find reaction (