The first multi-dimensional view of mass loss from externally FUV irradiated protoplanetary discs

TL;DR

This study compares 1D and 2D models of mass loss from externally FUV irradiated protoplanetary discs, revealing that 2D models show higher mass loss rates and different chemical structures, impacting observational predictions.

Contribution

It introduces the first detailed 2D-axisymmetric models of externally photoevaporating discs, highlighting differences from 1D models in mass loss rates and chemical structures.

Findings

2D models show mass loss rates within a factor of four of 1D models.

Half of the mass loss in 2D comes from the disc surface, which 1D models neglect.

Chemical abundances differ significantly, affecting observable signatures.

Abstract

Computing the flow from externally FUV irradiated protoplanetary discs requires solving complicated and expensive photodissociation physics iteratively in conjunction with hydrodynamics. Previous studies have therefore been limited to 1D models of this process. In this paper we compare 2D-axisymmetric models of externally photoevaporating discs with their 1D analogues, finding that mass loss rates are consistent to within a factor four. The mass loss rates in 2D are higher, in part because half of the mass loss comes from the disc surface (which 1D models neglect). 1D mass loss rates used as the basis for disc viscous evolutionary calculations are hence expected to be conservative. We study the anatomy of externally driven winds including the streamline morphology, kinematic, thermal and chemical structure. A key difference between the 1D and 2D models is in the chemical abundances. For instance in the 2D models CO can be dissociated at smaller radial distances from the disc outer edge than in 1D calculations because gas is photodissociated by radiation along trajectories that are assumed infinitely optically thick in 1D models. Multidimensional models will hence be critical for predicting observable signatures of environmentally photoevaporating protoplanetary discs.

Tables

Table 1. Table 1: A summary of the species included and initial gas abundances (taken from Asplund et al., 2009 ) for the reduced network used in this paper, which consists of 33 species and 330 reactions (Bisbas et al., 2012 ) . The sum of hydrogen atoms in atomic and molecular hydrogen is unity. The other abundances are with respect to atomic hydrogen.

Gas
Species	Initial abundance	Species	Initial abundance
H	$4\times {10}^{-1}$	H2	$3\times {10}^{-1}$
He	$8.5\times {10}^{-2}$	C+	$2.692\times {10}^{-4}$
O	$4.898\times {10}^{-4}$	Mg+	$3.981\times {10}^{-5}$
H+	0	H2+	0
H3+	0	He+	0
O+	0	O2	0
O2+	0	OH+	0
C	0	CO	0
CO+	0	OH	0
HCO+	0	Mg	0
H2O	0	H2O+	0
H3O	0	CH	0
CH+	0	CH2	0
CH2+	0	CH3	0
CH3+	0	CH4	0
CH4+	0	CH5+	0
e-	0
Dust
${\sigma }_{\text{FUV}}$	$2.7\times {10}^{-23}$ cm2	Dust cross section in wind
$\delta$	$3\times {10}^{-4}$	Dust-to-gas mass ratio in wind
$f_{\text{PAH}}$	1.0	PAH abundance relative to ISM
${\delta }_{\text{PAH}}$	$2.6\times {10}^{-2}$	PAH-to-dust mass ratio
Other
$\zeta$	$5\times {10}^{-17}$ s-1	Cosmic ray ionisation rate
		e.g. Bergin et al. (1999)

Table 2. Table 2: A summary of our model parameters. Columns are, from left to right: the model ID, stellar mass, disc mass, disc outer radius, surface density at 1AU, ambient UV field, location of the critical radius in an equivalent 1D model and the grid size.

Model	$M_{\ast }$	$M_d$	$R_d$	${\Sigma }_{1\text{AU}}$	UV	$R_{\text{crit}}$ (1D)	$R_{\text{grid}}$
ID	($M_{\odot }$)	($M_{\text{jup}}$)	(AU)	(g cm-2)	$G_0$	(AU)	(AU)
A	1	10	100	134	1000	207	385
B	1	20	200	134	1000	368	770
C	1	30	300	134	1000	380	1154

Table 3. Table 3: Mass loss rates 2D and equivalent 1D models, as well as their ratio.

Model ID	$\dot{M}$ 2D	$\dot{M}$ 1D	2D /
	(${10}^{-6}{M}_{\odot }$ yr-1)	(${10}^{-6}{M}_{\odot }$ yr-1)	1D
A	2.2	0.6	3.7
B	4.8	3.1	1.55
C	6.8	5.8	1.17

Table 4. Table 4: The variation of the 2D mass loss rates with spatial resolution. The equivalent 1D mass loss rate is included for reference. Note that we evaluate the mass loss rate in the lower resolution models only where the streamlines originate from the outer half of the disc ( R > R d / 2 𝑅 subscript 𝑅 𝑑 2 R>R_{d}/2 ). This is because in the higher resolution models all of the mass loss originates from the outer half of the disc, and the inner region of the lower resolution models is poorly resolved.

Model ID	$\dot{M}$ 2D 64	2D 128	2D 256	$\dot{M}$ 1D
	(${10}^{-6}{M}_{\odot }$ yr-1)	(${10}^{-6}{M}_{\odot }$ yr-1)	(${10}^{-6}{M}_{\odot }$ yr-1)	(${10}^{-6}{M}_{\odot }$ yr-1)
A	1.9	2.2	2.2	0.6
B	4.8	4.8	4.8	3.1
C	6.7	6.9	6.8	5.8

Table 5. Table 5: Mass loss rates in a low resolution ( 64 × 64 64 64 64\times 64 cell) version of model A, varying the level of healpix refinement (number of rays used to sample 4 ​ π 4 𝜋 4\pi steradians per cell).

healpix level	Rays per cell	$\dot{M}$
		(${10}^{-6}{M}_{\odot }$ yr-1)
0	12	1.9
1	48	2.1
2	192	2.0

Equations

$$\left(\frac{\chi }{1G_0}\right)=\frac{1}{N_l}\sum _1^{N_l}{\chi }_0\exp (-{\tau }_{UV}).$$

$$\Sigma (R)={\Sigma }_{1\text{AU}}{\left(\frac{R}{\text{AU}}\right)}^{-1}$$

$${\Sigma }_{1\text{AU}}=\frac{M_d}{2\pi R_d\text{AU}}.$$

$$H=\frac{c_s}{\Omega }$$

$$T_{\ast }=100{\left(\frac{R}{\text{AU}}\right)}^{-1/2}.$$

$$\frac{v^2\mu m_{\text{H}}}{k_{\text{B}}}-T-n\frac{\text{d}T}{\text{d}n}\ge 0.$$

$$\dot{M}=4\pi R^2\rho \dot{R}\mathcal{F},$$

$$\mathcal{F}=\frac{H_d}{\sqrt{H_d^2+R_d^2}}.$$

$$\dot{M_{i}}=\left\{\begin{array}[]{l l}2\times 2\pi{R}\Delta{l}\rho_{i}\sqrt{v_{R,i}^{2}+v_{z,i}^{2}},&v_{R}>0\,\,\,\,\textrm{or}\,\,\,\,v_{z}>0\\ &\\ 0&\textrm{otherwise}\\ \end{array}\right.$$

$$\rho =10^{-20}{\left(\frac{R}{100\text{AU}}\right)}^{-2}\text{g}{\text{cm}}^{-3}$$

$$R_g=\frac{G{M}_{\ast }\mu m_H}{k_{B}T}.$$

$$\frac{v^2}{2}+\Psi +\int \frac{\text{d}P}{\rho }=\text{Constant}$$

Full text

††pagerange: The first multi-dimensional view of mass loss from externally FUV irradiated protoplanetary discs –References††pubyear: 2016

1 introduction

The detection of almost 4000 very diverse exoplanets, most of which were discovered in the last 5 years (e.g. Winn & Fabrycky, 2015; Winn, 2018), coupled with new observing facilities such as the Atacama Large Millimetre Array (ALMA) has motivated and made possible a revolution in our understanding of planet formation. There is now overwhelming evidence that planets form from discs of material around young stars (for reviews see e.g. Williams & Cieza, 2011; Andrews, 2015; Owen, 2016; Sicilia-Aguilar et al., 2016; Haworth et al., 2016a; Ercolano & Pascucci, 2017). The vast majority of studies of disc evolution and planet formation consider the young star-disc in isolation and such approaches are proving very effective at understanding the myriad physical processes that take place in planet-forming discs.

The young disc-hosting stars themselves typically form in clustered groups of hundreds to hundreds of thousands (e.g. Lada, 2010; Krumholz et al., 2018). In order to fully understand the planet formation process, and in particular the populations of resulting planets, we also have to understand how the natal cluster can affect discs and any planet formation within them. Specifically, any given planet-forming disc can be influenced either by gravitational encounters or irradiation by other cluster members (for a recent assessment of their relative impacts, see Winter et al., 2018). For example Ndugu et al. (2018) accounted for stellar neighbours through simple passive heating of the outer disc in planetary population synthesis models, finding that this is necessary to suppress large populations of cold Jupiters (particularly at low metallicity) which are not observed.

In particular, UV irradiation by stellar neighbours can heat a disc and drive material away from it in a photoevaporative wind. If this mass loss is higher than the rate of viscous spreading, the disc will also be truncated (e.g. Clarke, 2007; Anderson et al., 2013; Haworth et al., 2017; Rosotti et al., 2017; Winter et al., 2018; Shadmehri et al., 2018). Although the effect of external photoevaporation is expected to act most directly upon the outer disc, limiting the disc mass reservoir and radial extent affects the ability of material to move into the inner disc which may also have consequences for planet formation at small orbital radii (e.g. as in the case of Trappist-1 Ormel et al., 2017; Haworth et al., 2018a). Furthermore, truncation of the disc lowers the viscous timescale and affects the redistribution of angular momentum and hence surface density evolution of the entire disc, an effect which does not occur through internal photoevaporation by the host star, where instead an inner hole is eventually produced (Clarke et al., 2001; Owen et al., 2010). The effect of external photoevaproation on planet formation should therefore not be discounted simply because it operates on the disc at larger radii than those at which most planets are being discovered.

The range of star forming environemts in the Galaxy implies that stars may be exposed to ambient ultraviolet fluxes that vary over five orders of magnitude (e.g. Fatuzzo & Adams, 2008); it is therefore possible that some of the exoplanet diversity being discovered is influenced by diversity in the properties of the natal cluster (e.g. Nicholson et al., 2019). Note that although gravitational interactions in the cluster do take place, generally, the effect of external UV irradiation has a dominant influence (Scally & Clarke, 2001; Guarcello et al., 2016; Winter et al., 2018). Internal winds are also driven by the host star (see e.g. the review by Ercolano & Pascucci, 2017) though the relative impact and interplay between internally and externally driven winds is currently unknown.

A major complication is that environmental evaporation of discs is only easy to observe and/or model in extreme UV conditions, for example in the vicinity of O stars. In such a regime disc evaporation is unsubtle, giving rise to “proplyds” that can be observed in the optical, silhouetted against the H ii region (e.g. McCaughrean & O’dell, 1996; O’dell, 1998; O’Dell, 2001; Henney et al., 2002; Wright et al., 2012). Trends in disc properties such as mass and/or radius near such sources are also starting to be inferred (e.g. Ansdell et al., 2017; Eisner et al., 2018).

The majority of star/discs are not in such an extreme UV environment as in the vicinity of an O star. For example Fatuzzo & Adams (2008) compute probability distributions for the UV environments of Galactic star forming regions. Proplyds lie in the tail end of this distribution, even in the limit of no extinction (and hence more pervasive high UV fields). In the weak-intermediate UV regime environmental disc photoevaporation is more subtle, making it difficult to observe.

It is also much more difficult to model environmental photoevaporation in more modest UV environments, with the thermal state of the gas being sensitive to the far-ultraviolet (FUV) radiation field and photodissociation region (PDR) chemistry (e.g. Adams et al., 2004). In particular, a key issue is that the main coolant is the escape of line photons, which cannot be assumed to be optically thin. To estimate this the escape probability into three dimensions has to be computed from every point on a computational domain (e.g. Bisbas et al., 2012). Although within the disc itself this can be approximated by assuming some dominant trajectory (e.g. vertically and/or outward radially Wang & Goodman, 2017; Nakatani et al., 2018; Wang et al., 2018) in a photoevaporative wind this is not necessarily the case. A 3D estimate is required, which makes multidimensional models difficult to develop and, historically, computationally prohibitive. Given that there is growing observational evidence for external disc photoevaporation in more modest environments with the detection of proplyds around a B star in Orion (Kim et al., 2016) and a possible photoevaporative halo around IM Lup (Cleeves et al., 2016; Haworth et al., 2017; Pinte et al., 2018) further models and observations in these lower UV regimes are essential.

The difficulty of including 3D line cooling means that modelling external disc photoevaporation in the FUV regime has been restricted to 1D models to date (though at very high UV field strengths the role of the FUV can be more easily approximated in 2D axisymmetric calculations Richling & Yorke, 2000). These assume that mass loss occurs predominantly from the disc outer edge (see e.g. Adams et al., 2004, for a discussion on this), since this part of the disc is a substantial mass reservoir that is least gravitationally bound to the star. The exciting UV radiation is assumed to propagate from the outside inwards and vice versa for the line cooling. Initially these models were computed semi-analytically, locating the critical point of the flow and matching to conditions at the flow base (disc outer edge), with the temperature calculated from a pre-computed table as a function of local density, UV field and extinction (Adams et al., 2004; Holden et al., 2011; Facchini et al., 2016). Although solutions are not always possible with this approach, where they are possible mass loss rates that would be significant for disc evolution are predicted (though all such models only considered Solar type stars).

In recent years direct simulations of evaporating discs have become possible (Haworth et al., 2016b, 2017, 2018a, 2018b). These have still been 1D, but directly compute the PDR chemical/thermal structure and hydrodynamics. For example Haworth et al. (2018b) produced a publicly available grid of steady state mass loss rates for a large variety of stellar, disc and UV parameters. A key general result is that mass loss can be important even in very weak UV environments if the disc is only weakly bound, for example in the case of a very extended disc like IM Lup (Haworth et al., 2017) or if the stellar mass is low as would have been the case for the precursor of a Trappist-1 like system (Haworth et al., 2018a). Although these 1D models are easily calculated with modern facilities, they offer only a weak means of predicting the crucial signatures of externally evaporating discs which we require in order to identify the process in action in a range of UV environments. Multidimensional models are required to effectively predict observables. Furthermore, since the 1D models do predict significant mass loss, even in weak UV environments, they must be validated and improved upon with more realistic multidimensional models.

In this paper we present the first 2D cylindrical models of externally FUV irradiated protoplanetary discs. There is a large parameter space where the nature of external photoevaporation may differ (e.g. as a function of UV field strength, stellar mass, disc mass, disc size, metallicity, and so on) so in this paper we focus on only a handful of first calculations. We aim to gain a first insight into the applicability of 1D mass loss rates, as well as test the validity of assumptions in the 1D models such as the mass predominantly being driven from the disc outer edge. We also aim to study the flow morphology and take a first look at the chemical composition of 2D externally driven winds.

2 Numerical method

We begin with an overview of the methodology that permits us to solve PDR-dynamics in arbitrary geometries with 3D line cooling.

2.1 General Overiew

We use the photochemical-hydrodynamics code torus-3dpdr Bisbas et al. (2015b) to run the models in this paper. In its most general form this is a coupling of the grid-based Monte Carlo radiative transfer and hydrodynamics code torus (Harries, 2000; Haworth & Harries, 2012; Haworth et al., 2015; Harries, 2015; Ali et al., 2018) with the 3D photodissociation region modelling code 3d-pdr (Bisbas et al., 2012)111Though note that the Monte Carlo radiation transport is not actually employed within this paper, as we will discuss in section 2.3.. In this hybridisation all calculations take place on the torus grid. That is, components have been stripped from 3d-pdr and directly incorporated into torus. The original idea behind the coupled code is that the Monte Carlo radiation transport of torus can accurately compute the exciting UV radiation field (which 3d-pdr approximates as, e.g. isotropic, planar or spherical) for the PDR modules, which is important for models of the interstellar medium. However the application of using the PDR thermal structure in hydrodynamical problems such as external disc photoevaporation has been a powerful by-product that has dominated the use of torus-3dpdr to date.

PDR-dynamics calculations are performed by iteratively computing hydrodynamics and PDR equilibirium updates (i.e. via operator splitting), which we now summarise separately.

2.2 Hydrodynamics

The hydrodynamics scheme is a finite volume method which employs a van Leer (1979) flux limiter, Rhie & Chow (1983) interpolation and is total variation diminishing. In these calculations the disc masses are low enough that self-gravity is unimportant, so we include only a point source gravitational potential. The hydrodynamics scheme was first presented in torus and detailed in Haworth & Harries (2012).

2.3 Equilibrium PDR chemistry and thermal balance

The PDR aspect of our calculations require solving a chemical network iteratively with thermal balance (for full details see Bisbas et al., 2012, 2015b). The fact that the PDR calculations in 3d-pdr take place in 3D for arbitrary geometries is crucial for our models of externally photoevaporating protoplanetary discs. To compute the line cooling escape probabilities, $4\,\pi$ steradians is sampled from the centre of every cell on the computational domain using a healpix scheme (Górski et al., 2005) which divides the sky into zones of equal solid angle. The size of these zones is dependent on the level of healpix refinement $l$, with a number of zones $N_{l}=12\times 4^{l}$. The optical depth, escape probability, etc. along each of these rays is then computed to estimate the line cooling. These rays are also used to estimate the local UV, “$\chi$”, (i.e. we are not using the Monte Carlo radiative transfer of torus) as attenuated from an isotropic ambient UV field $\chi_{0}$ according to

[TABLE]

where $\tau_{UV}$ is the UV optical depth, which is related to the visual extinction by $\tau_{UV}=3.02{A_{V}}$. $G_{0}$ is the Habing unit of UV field strength Habing (1968). Although our models are 2D-axisymmetric, healpix sampling of the sky is still done in 3D utilising the symmetry of the problem.

In its current form our healpix scheme uses a long characterstics ray tracing, for which the scaling is relatively poor. Multiple rays have to be traced from every single cell on the grid, which is computationally expensive. Furthermore, since there can be significant memory requirements we implement domain decomposition (the grid memory is distributed amongst cores) which means rays may not have access to the information stored on all parts of the grid along their trajectory. Most domains hence have to act as servers (sending relevant information to the domain doing the ray tracing) and take turns doing ray-tracing, rather than tracing rays simultaneously. This is added to the fact that the PDR calculation takes place over many iterations to converge on the equilibrium chemical and thermal structure, so the ray tracing process has to take place many times for any given single PDR update. Of course this then has to be done multiple times iteratively with hydrodynamics updates. The 3D line cooling/PDR calculation therefore overwhelmingly dominates the computational expense of these models.

Given the poor scaling we generally run these calculations on a small number of cores. For example on a $256^{2}$ cell 2D cylindrical grid (where the ray tracing takes place in 3D, utilising the symmetry of the problem) a single PDR update at the lowest level of healpix refinement can have a wall time of around 12-24 hours on 5 MPI threads (with remaining threads on the node being used for shared memory parallelisation). In future we intend to improve on the wall time per PDR calculation by improving the scaling to larger numbers of cores. In light of the above, we are presently limited to calculating the steady state flow structure from an irradiated disc, and also at modest resolution both spatially and in terms of healpix rays. Nevertheless in this paper we do consider convergence with both of these quantities in section 3.4.

The details of the PDR calculations themselves are given in all of our series of papers on 1D calculations, the latest of which was presenting the fried grid of mass loss rates (Haworth et al., 2018b). Nevertheless, the initial abundances are also summarised in Table 1 of this paper. We use the same reduced PDR network of 33 species and 330 reactions. This has been tailored to give temperatures accurate to within around 10 per cent of the UMIST 2012 chemical network database of 215 species and over 3000 reactions (McElroy et al., 2013). The chemistry and thermal balance are solved iteratively until convergence. The main coolants are lines from C, C+, O and CO, with some additional contribution from the dust. Heating processes include photoelectric heating from atomic layers of polycyclic aromatic hydrocarbons, C ionisation, H2 formation and photodissociation, FUV pumping, cosmic rays, turbulent and chemical heating and gas-grain collisions.

The only difference to the microphysics in these calculations compared to those of the fried grid is that we use an interstellar medium (ISM)-like polycyclic aromatic hydrocarbon (PAH)-to-dust mass ratio of $2.6\times 10^{-2}$, an order of magnitude higher than that used in the fried grid. PAHs are important, possibly being the dominant heating mechanism in the PDR (see Figure 2 of Facchini et al., 2016, which summarises all of the heating and cooling processes of the code) but the PAH abundance in discs is very uncertain (Geers et al., 2006; Oliveira et al., 2010; Perez-Becker & Chiang, 2011). In prior work we chose lower values of the PAH abundance to be conservative as to the mass loss rate, but here the higher (ISM like) value leads to a warmer wind that can help achieve a steady state more quickly and hence at less computational cost. Here we are most interested in comparing 1D and 2D models, so the absolute mass loss rate is less important. We note that any 1D models in this paper are also bespoke models rather than being taken from the fried grid, and hence are also tailored to have the same ISM-like PAH abundance.

Finally we note that the microphysics in these calculations is only accurate in the photodissociation regime. At deeper, colder, denser layers within the disc that are extremely optically thick to the UV the chemistry becomes much more complicated and is not captured by our model (see e.g. Walsh et al., 2010; Woitke et al., 2016; Cleeves et al., 2017; Kamp et al., 2017; Wakelam et al., 2019). Although this is unimportant for computing the dynamical evolution of the disc (since such zones are cold regardless of the chemical complexity) we do not attempt to properly capture the composition deep within the disc. Effects such as freeze out upon, for example, the CO abundance distribution can be retrospectively applied for predicting observables.

2.4 Further details

As discussed above, our calculations involve iteratively solving hydrodynamic and PDR updates via operator splitting. At this stage we are only interested in the eventual steady state wind solutions and the associated mass loss rate, so we only do a PDR update every $N$ timesteps. In practice we start off with $N$ being $\sim 10^{4}$, but as the calculation approaches its final steady state increase $N$ to even larger values ($\sim 10^{6}$). This means that the pathway to achieving the steady state cannot be studied in detail, but drastically reduces our computational expense. This is the approach we have historically used in our 1D models, for which we confirmed that the frequency with which these updates are made does not affect the final solution (Haworth et al., 2016b). Note that our models assume chemical and thermal equilibrium. In section 5.2 of Haworth et al. (2016b) we supported this assumption by showing that the thermal timescale is much shorter than the flow timescale.

2.5 Validation

torus-3dpdr has been benchmarked against 1D semi-analytic solutions in the context of externally irradiated discs in Haworth et al. (2016b). The code was also validated against a series of dynamical and radiation hydrodynamic tests in Haworth & Harries (2012), Harries (2015) and Bisbas et al. (2015a). The PDR components have also been checked against the Röllig et al. (2007) benchmarks in Bisbas et al. (2015b)

The healpix scheme used by torus-3dpdr is known to work in 3D applications, however it has been modified for this work to utilize the cylindrical symmetry of our simulation grid when propagating through 3D space. To check this modified scheme we computed the UV field impinging upon a sphere embedded in a low density medium using both the 3D and 2D cylindrical scheme, ensuring that the same UV field in the sphere resulted.

2.6 Disc construction/boundary conditions

In these calculations we define a disc that is imposed and hence acts as a boundary condition and isn’t allowed to dynamically evolve (in the first instance we are interested in obtaining steady state solutions to compare with 1D models). For this purpose we consider a truncated power law for the imposed disc surface density profile

[TABLE]

where for a disc of a given size $R_{d}$ and mass $M_{d}$ the surface density normalization is

[TABLE]

We impose the disc up to one scale height

[TABLE]

where $c_{s}$ is the sound speed and $\Omega$ the Keplerian angular velocity. This scale height is computed assuming an isothermal disc temperature of 20 K. We do allow the disc itself to be subsequently heated in the PDR calculation, but do not then adjust the height above the mid-plane to which the disc is imposed (i.e. we always impose it up to 1 scale height as though the disc temperature were 20K). We retrospectively checked the disc mid-plane temperature profile and it is in the range 20–23 K in the bulk of the outer disc, excluding the outermost few cells which can be heated to higher temperatures. Where the external radiation field becomes extremely weak the source star may set the temperature through passive irradiation, which we assume is of the form

[TABLE]

We take the maximum of the external radiation and stellar-induced temperatures, with a floor value of 20 K. Note that the assumed temperature of 20 K in setting the imposed scale height is accurate to within a few degrees in the mid-plane except for the outermost $\sim 8$ cells, though of course at higher $z$ the disc gets significantly warmer, as we will illustrate in section 3.

2.7 Model parameters

We consider 3 calculations in this paper, the parameters of which are summarised in table 2. At this stage we are considering only a 1 M⊙ star, and all of our models consider a $10^{3}$ G0 field. We do however vary the disc mass and radius such that the surface density normalisation is constant. Of course there will likely be interesting variations in the nature of external disc photoevaporation in different stellar mass and UV regimes that will need to be explored in future work.

In 1D models the required size of the grid is tightly constrained. A unique solution for a steady transonic flow satisfies a criticality condition (simultaneous vanishing of terms in the combined momentum and continuity equations, e.g. Parker, 1965; Clarke & Alexander, 2016) at some point in the flow, given by that at which

[TABLE]

is first satisfied for velocity $v$, mean particle mass $\mu$, temperature $T$ and number density $n$ (Facchini et al., 2016). 1D models will only converge on this solution if the critical point is contained within the grid. Our 2D grid is chosen to be larger than this critical radius for the 1D models in each case.

The 3 main calculations are predominantly presented using a fixed $256\times 256$ cell cylindrical $R-z$ grid and the lowest level of healpix refinement (12 rays per cell, propagated in 3D utilising the axisymmetry of our model). In section 3.4.2 we also explore convergence, running $128\times 128$ and $64\times 64$ cell versions, and for the lowest resolution case also using the next 2 levels of healpix refinement (48 and 192 rays per cell), though we are unable to practically increase both the spatial and healpix resolution further at this stage.

2.8 Computing the mass loss rate

One of the primary quantities of interest from these calculations is the mass loss rate, which is important for the dynamical evolution and lifetime of the disc. In 1D calculations the mass loss rate is computed following Adams et al. (2004) using

[TABLE]

where $\mathcal{F}$ is the fraction of solid angle subtended by the disc outer edge

[TABLE]

In our new 2D calculations the mass loss rate is computed through a spherical surface of radius $r$ at the end of the calculation in a postprocessing step as follows. We generate a number of evaluation points on a circular contour in the ($R,z$) plane of radius $r$ centred on the origin. From any given evaluation point we draw a tangent to this circle and calculate the length of this tangent that lies within the cell that contains it. We ensure that each $i^{th}$ cell only contributes once to the mass loss rate estimate, which is given by

[TABLE]

where $v_{R,i}$ and $v_{z,i}$ are the $R$ and $z$ components of the velocity and the preceding factor of 2 arises because the models in this paper are reflective about the mid-plane.

We tested this scheme for a spherically symmetric diverging flow of constant velocity and inverse square density profile (hence a constant “mass loss rate” as a function of radius). For a spherically diverging velocity of $10^{-3}$ km s*-1* and density profile of the form

[TABLE]

the analytic mass loss rate is $4.46\times 10^{-13}$ M⊙ yr*-1*. The mass loss rate we obtain by summing equation 9 over a spherical surface, as a function of the radius of that spherical surface, is given in Figure 1. This was for the case of a fixed 385 AU 2D cylindrical grid with $256^{2}$ cells. Our discretised estimates agree with the analytic value to within 1 per cent for all spherical radii through which we compute the mass loss rate.

3 Dynamics and chemistry of multidimensional external disc photoevaporation

The analysis of our models proceeds as follows. In section 3.1 we make an initial overview of the steady state flow, thermal and kinematic structure. In section 3.2 we then quantify the location in the disc from which the mass loss originates. With an understanding of the flow structure, we then compute and compare 1D and 2D mass loss rates in 3.3. We then discuss convergence of our models in 3.4 and finish with a first look at the chemical composition of externally irradiated discs in 2D in 3.5.

3.1 Initial overview of flow structure

We begin with an initial broad overview of the steady state flow structure in our models. Figure 2 shows the density and temperature structure, as well as a selection of flow streamlines. Figure 3 also shows the Mach number, colour coded such that subsonic flow is blue, supersonic red and the approximate sonic surface is white. In all cases the cold outer disc is surrounded by a warm ($\sim 50-150$ K) layer of subsonic material of typical density $\sim 10^{8}$ cm*-3*. The edge of this “halo”222We now refer to the disc region that we impose as “the disc” and any CO bearing region outside of this as “the halo”, motivated by the terminology used in the interpretation of IM Lup (Cleeves et al., 2016; Haworth et al., 2017) corresponds to the attainment of an optical depth $\sim 1$, which coincides with the surface of photodissociation of CO (see section 3.5) and the transition to a hotter ($\sim 200-350$ K) supersonic flow. Above the disc surface, interior to about half of the disc radius (a region that contributes negligibly to the mass loss: see section 3.2) the flow is subsonic up to a large height. At small $R$ a hot ($\sim 1500$ K in the $10^{3}$ G0 cases) low density ($10^{-3}$ cm*-3*) cavity is blown along the $z-$axis which is subject to small scale changes over time due to small scale shearing instabilities. The rest of the flow and hence mass loss rate is steady, which we will illustrate further in section 3.4.

The streamline morphology, as illustrated in the upper panels of Figure 2 as well as in more detail for model B in Figure 4, is similar in each of our three cases. Streamlines emanating from the disc surface leave vertically, before turning over and becoming quasispherical at larger radii. Streamlines emanating from the disc outer edge flare such that the expansion is greater than spherical out to approximately 1.5 disc radii, before again turning over to become more spherical.

3.2 Where is mass lost from the disc?

Before discussing the total mass loss rates in our models we first discuss the radial location in the disc from which the bulk of the mass loss originates, as well as how it propagates in the wind. This is a prerequisite to understanding any differences in the mass loss rate between the 1D and 2D models.

The gravitational radius $R_{g}$ is that at which the sound speed is equal to the escape velocity (or equivalently the thermal energy equals the binding energy) and hence matter is unbound

[TABLE]

(e.g. Hollenbach et al., 1994; Johnstone et al., 1998). Scenarios where $R_{g}$ is greater than the disc radius are referred to as sub-critical, and a key assumption of the 1D models of externally evaporating discs in the sub-critical regime is that mass loss is dominated by that from the disc outer edge. 1D models then compute the mid-plane flow structure and calculate a total mass loss rate assuming a spherically diverging wind (as in equations 7 and 8).

In our models, which are all sub-critical (we verified this by checking maps of $R_{g}/R_{d}$) we can compute exactly where the mass loss originates from in the disc and compare with the 1D assumptions above. Recall from the description of our method in 2.8 that we compute the mass loss rate in the 2D models though a spherical surface. From each grid cell on this surface (where we know the local contribution to the total mass loss rate) we trace a streamline back to its origin at the disc. An illustration of this for an example spherical surface and handful of streamlines in model C is given in Figure 4.

Using this streamline tracing approach, in the upper panel of Figure 5 we plot the cumulative mass loss rate through the spherical surface as a function of the $R-$coordinate at which the streamline originates from the disc. The lower panel shows the same information, normalised to the total mass loss rate and disc outer radius. These plots demonstrate that in all cases around 50 per cent of the mass loss rate comes directly from the disc outer edge. Although some mass loss does come from the surface, the mass loss rate is entirely set in the outer half of the disc and is mostly from within about 10 per cent of the disc outer edge. That is, in the lower panel of Figure 5 the fraction of cumulative mass loss is 1 by $R/R_{d}=0.5$. Unsurprisingly larger (less bound) discs lose mass from a larger fraction of the disc surface. Overall the assumption from 1D models that the mass loss is predominantly from the disc outer edge therefore seems reasonable to within a factor two or so. We will compare the actual 1D and 2D mass loss rates in section 3.3.

We further illustrate the origin of mass loss in our models in Figure 6. This shows the spatial ($R-z$) location of the streamline end points near the disc outer edge. These points are colour-coded by the cumulative mass loss rate. Starting from the dark blue points, one can see what fraction of the mass loss originates from certain locations. For example, at the top of the disc outer edge the points are green, corresponding to roughly half of the total mass loss rate.

3.3 Comparing 1D and 2D mass loss rates

The key quantity of interest from external photoevaporation for disc evolutionary models is the mass loss rate. Mass loss rates from 1D models are being used in viscous evolutionary calculations to study disc evolution in different environments, where the role of environment is being predicted to be important (e.g. Clarke, 2007; Anderson et al., 2013; Haworth et al., 2018a, b; Rosotti et al., 2017; Winter et al., 2018; Shadmehri et al., 2018). It is crucial to confirm the validity (or otherwise) of these mass loss rates in multidimensional models.

We ran 1D counterparts to our 2D models to compare the mass loss rates (e.g. similar to the fried grid models Haworth et al., 2018b, only with the required higher PAH abundance used here). These 1D models would be defined as having the same disc mass and radius, as well as stellar and irradiating UV field properties as the 2D models. The 1D models do differ in that the scale height at the disc outer edge is allowed to increase due to heating at the disc outer edge (as adopted in the fried grid). Note that in the 2D models we calculate the mass loss rate through a number of spherical surfaces over the supersonic flow and take the average (it should be constant in a perfect steady state and in section 3.4.1 we will show that this is the case to within $\pm 3$ per cent). For reference, the mid-plane flow variables for model A, which are typical of each of the models, are shown in Figure 7. The 2D models result in more rarefied, faster flows in the supersonic wind along the mid-plane.

Table 3 compares the total mass loss rate in our 2D models with their 1D analogues. We are probing only a very limited parameter space in this paper so broad generalisations must be treated with caution. However, in all of the models here the mass loss rate is larger in the 2D cases. So the contribution of mass loss from the disc surface and faster flow more than compensates for the 1D models being based entirely on the peak density part of the flow along the mid-plane.

Recall that 1D models have been predicting that external photoevaporation of discs is important for disc evolution (and potentially planet formation) in a wide range of regimes. Our 2D models are suggesting that in fact the 1D mass loss rates are probably underestimates (albeit by factors of order unity) and so the impact is, if anything, likely stronger than previously anticipated.

3.4 Convergence

We now briefly assess the convergence properties of our models.

3.4.1 Sensitivity of mass loss rate to size of spherical surface

We calculated the mass loss rate in our models through spherical surfaces of different sizes and used the mean. In 2.8 we showed that the variation for a perfect spherically diverging flow is only about 1 per cent over 200 AU variation in surface radius. Figure 8 shows the percentage variation of the mass loss rate, relative to the mean as a function of the radius of the spherical surface. Generally the mass loss rate slightly decreases with surface size, but is flattening out at larger radii (as the number of cells contributing to the estimate increases). Overall the variation is never more than $\pm 3$ per cent. That this is higher than the 1 per cent level in the perfect case discussed in 2.8 will in part be due to small deviations from a perfectly steady state, which make a bigger perturbation to the total mass loss rate at smaller radii (where there are a smaller number of cells on the surface) than at large radii.

3.4.2 Conservation along streamlines

We checked the variation of specific angular momentum and Bernoulli parameter

[TABLE]

which should be invariant along streamlines in a steady state flow. Here $v$, $\Psi$ and $\rho$ are the local velocity, gravitational potential and density. Note that we evaluate the pressure term $\textrm{d}P$ in the Bernoulli parameter in an approximate way by just taking the difference over a local length interval of the streamline. We checked a posteriori that the flow is barotropic along a streamline (i.e. that the density is single valued function of pressure). The specific angular momentum and Bernoulli parameter are illustrated in the case of model B in Figure 9. Both quantities are tightly conserved along streamlines, with deviations from perfect conservation arising from a mixture of the approximate way in which we evaluate the $\textrm{d}P$ term and small deviations from steady state. At the transition to a supersonic flow (at an optical depth $\sim 1$ and, in this case, the CO dissociation front) there is a large pressure gradient, which introduces numerical errors into our evaluation of the $\frac{\textrm{d}P}{\rho}$ term and corresponds to a jump in Bernoulli parameter. However once through to the supersonic flow the Bernoulli parameter is again conserved. We have checked that the mass and momentum fluxes are conserved within the region of rapidly changing thermodynamic variables.

3.4.3 Convergence with spatial resolution

There are two types of resolution element in these calculations: grid cells which quantify the spatial resolution and solid angle elements making up the healpix sampling of the sky from each point on the grid. Owing to the computationally expensive nature of these calculations we make a first exploration of resolution convergence with spatial resolution by coarsening our grid from $256\times 256$ cells to $128\times 128$ and $64\times 64$ and hence using relatively cheap models. We also explore the effect of using higher healpix resolution (which makes the most expensive part of the calculation and order of magnitude, or more, more expensive) in the case of the $64\times 64$ cell model.

In our fiducial models all of the mass loss comes from the outer half of the disc ($R>R_{d}/2$, see Figures 5, 6). At lower resolution when the pressure gradient in the inner disc is very poorly resolved material can artificially be driven from the inner disc. Therefore, when computing the mass loss rate in the lower resolution models we only consider that from streamlines that originate from $R>R_{d}/2$.

The mass loss rate at different spatial resolution is given in Table 4. In each case the variation in mass loss rate with resolution is small. Although we cannot and do not attempt to claim to have fully converged with spatial resolution (something which can be extremely challenging to do and is rarely genuinely achieved), it is promising that the mass loss rate appears to be insensitive to the spatial resolution over the range we have considered.

3.4.4 Convergence with HEALPIX resolution

In the case of the $64\times 64$ cell spatial grid model we increased the level of healpix refinement from $l=0$ (the fiducial value here) to $l=1$ and $l=2$, where the number of healpix zones is $N=12\times 4^{l}$. That is we consider 12, 48 and 192 rays per cell (the $l=2$ case is roughly an order of magnitude more expensive than the $l=0$ case, which is why we explore convergence with healpix rays on a coarse spatial grid). The mass loss rate for different levels of healpix refinement in the case of model A are given in Table 5. As with spatial resolution, we cannot claim convergence, but there is no strong variation in the mass loss rate with healpix refinement. The reason for the slight increase in mass loss rate that can happen at higher healpix resolution is that there is slightly higher temperature/pressure and hence driving of a wind from the disc surface layers.

3.5 Chemical composition

To solve for the dynamics our calculations necessarily solve a PDR chemistry network and therefore naturally yield the chemical structure of the disc/wind. Figure 10 illustrates some features of the composition of our models. The right hand panels are maps of the UV field strength, which attenuates sharply near the halo boundary and rapidly drops to a negligibly low value within the disc itself. Ionisation and dissociation contributing to the chemistry of the disc itself is therefore more sensitive to the cosmic ray ionisation rate than the ambient UV field (for more information on disc chemistry including assessments of the cosmic ray contribution see e.g. Cleeves et al., 2014; Cleeves et al., 2015; Wakelam et al., 2019). Though as we will shortly see parts of the wind are molecular, so may be identified as part of the disc in molecular line observations and certainly would be influenced by the ambient UV field.

The left and middle columns of of Figure 10 illustrate in which zones certain species dominate in abundance over one another. In the left panels the blue regions are molecular hydrogen and the red atomic and so the H–H2 transition surface is at their interface. Similarly in the middle hand panels we illustrate which of CO, C I and C II dominate over one another in certain zones. The CO halo region extends beyond our imposed disc structure. This is then enveloped by a layer of C I, beyond which the medium becomes warmer and C II is most abundant of these 3 species. Comparing the chemical distribution with the dynamical properties of Figures 2 and 3, the carbon cycle is more strongly correlated with the dynamics than the H–H2 transition.

One of the predictions of 1D models is the the azimuthal velocity in the flow is sub-Keplerian (Facchini et al., 2016; Haworth et al., 2016b) which might be identifiable observationally. Indeed observational hints of sub-Keplerian rotation were found in the case of the disc IM Lup by Pinte et al. (2018). Figure 11 shows the spatial distribution of $v_{\phi}/v_{\textrm{kep}}$ for each model. Marked on with dashed lines are the surfaces at which C I transitions to being dominant over CO (the CO surface) and where C II transitions to being dominant over C I (the C I surface). Deviations from Keplerian rotation in the wind can be very strong, but in the models considered here this is mainly in the atomic and ionised carbon zones. Atomic carbon may offer a promising alternative for probing environmentally driven winds and detecting sub-Keplerian rotation, though we will explore observational signatures using synthetic observations in a subsequent paper.

Finally, in Figure 12 we compare the chemical profile of some species in the 1D models with the mid-plane chemical profile from our 2D models. The panels show models A-C from top to bottom. Solid lines are the 2D distribution and points the 1D. There is a common theme, at least in these models, which is that the radial extent of molecular gas is significantly reduced in the 2D models compared to 1D. This is because the 1D models assume that exciting radiation only propagates inwards radially and that other trajectories are infinitely optically thick. Heating of the mid-plane from trajectories above the mid-plane is responsible for dissociating gas more effectively than the 1D models predict, decreasing the CO extent. The consequences of this are that (at least) 2D models will be required to properly predict observable signatures of externally driven winds, and also that they will likely be harder to detect than we might have expected from prior 1D models, at least at $10^{3}$ G0 UV field strengths.

4 Summary and conclusions

We have computed the first 2D PDR-dynamical models of externally FUV irradiated protoplanetary discs, including direct computation of the PDR microphysics and line cooling in 3D. Previously, only 1D (usually semi-analytic) calculations of this have been possible. In this paper we aim to make a first study of the multidimensional flow structure, calculate mass loss rates and compare with the 1D calculations that are important for providing mass loss rates for disc viscous evolutionary models. We draw the following main conclusions from this work.

1. Mass loss rates from our 2D models are slightly higher than in equivalent 1D cases, to within a factor $\sim 4$. Differences can arise from a few factors. The 1D model mass loss rates are enhanced relative to the 2D because they consider the mid-plane only and assume that this densest part of the flow represents what would propagate spherically outwards from an entire scale height at the disc outer edge. Nevertheless, this is more than compensated for by the fact that we find that only about half of the mass loss comes from the very disc outer edge (it is all assumed to originate from there in the 1D models) with the rest coming from the disc surface. We conclude that 1D estimates of mass loss rates can reasonably be used in other applications, especially since they are expected to be conservative estimates. This is an important result since 1D evolutionary models find that external photoevaporation is important for disc mass and radius evolution. It thus controls the viscous timescale and hence the surface density evolution of the entire disc and so could even affect inner planet formation.

2. Aside from the mass loss, a key difference between the 1D and 2D models is the chemical structure. This arises because the 1D calculations assume that exciting radiation only propagates radially inwards and that other trajectories (i.e. from above the mid-plane) are infinitely optically thick. This means that CO is dissociated at smaller radii in 2D models compared to 1D. Although 1D mass loss rates are reasonable (and conservative) multidimensional models such as those in this paper are going to be essential for predicting observables of externally evaporating discs using synthetic observations.

3. Sub-Keplerian rotation in the wind was a signature of external photoevaporation predicted by 1D models (Facchini et al., 2016; Haworth et al., 2016b). We do still anticipate that sub-Keplerian rotation could be a viable signature of external photoevaporation; however at the intermediate UV field strengths considered here it is anticipated to be difficult to detect in CO (at lower UV field strengths we anticipate that stronger sub-Keplerian deviations will be in CO-bearing parts of the flow). Rather, at intermediate UV field strengths our models suggest that atomic carbon lines could offer a promising alternative tracer of the wind and sub-Keplerian rotation, which we will explore with subsequent synthetic observations.

Acknowledgements

We thank the referee for their swift and constructive review that improved the manuscript. We thank James Owen for useful discussions and his comments on this manuscript. TJH is funded by an Imperial College London Junior Research Fellowship. CJC has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC- 2013-ADG. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. The work also used the COSMOS Shared Memory system at DAMTP, University of Cambridge operated on behalf of the STFC DiRAC HPC Facility. This equipment is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. DiRAC is part of the National E-Infrastructure. It also used the DiRAC Data Analytics system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Serve on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1. The main DiRAC project ID for these calculations was dp100 (PI: J. E. Owen).