First evidence of external disc photoevaporation in a low mass star forming region: the case of IM Lup

TL;DR

This study provides the first evidence that external photoevaporation significantly affects the evolution of large protoplanetary discs around low-mass stars, demonstrated through modeling of IM Lup's disc and its observed CO halo.

Contribution

It models radiatively driven disc flow in a low-radiation environment, showing external photoevaporation's role in disc truncation even in low-mass star forming regions.

Findings

Material at >400AU is weakly gravitationally bound, enabling mass loss.

Photoevaporation causes disc truncation primarily in the first Myr.

Synthetic observations match the observed CO halo extending beyond 1000AU.

Abstract

We model the radiatively driven flow from IM Lup -- a large protoplanetary disc expected to be irradiated by only a weak external radiation field (at least 10$^4$ times lower than the UV field irradiating the Orion Nebula Cluster proplyds). We find that material at large radii ($>400$AU) in this disc is sufficiently weakly gravitationally bound that significant mass loss can be induced. Given the estimated values of the disc mass and accretion rate, the viscous timescale is long ($\sim10$Myr) so the main evolutionary behaviour for the first Myr of the disc's lifetime is truncation of the disc by photoevaporation, with only modest changes effected by viscosity. We also produce approximate synthetic observations of our models, finding substantial emission from the flow which can explain the CO halo observed about IM Lup out to $\geq1000$AU. Solutions that are consistent with the extent of the observed CO emission generally imply that IM Lup is still in the process of having its disc outer radius truncated. We conclude that IM Lup is subject to substantial external photoevaporation, which raises the more general possibility that external irradiation of the largest discs can be of significant importance even in low mass star forming regions.

Equations

$${\Sigma }_{\text{g}}(R)={\Sigma }_c{\left(\frac{R}{R_c}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{R_c}\right)}^{2-\gamma }\right],$$

$$H(R)=H_{100}{\left(\frac{R}{100\text{AU}}\right)}^{\psi },$$

$$n(R)=\frac{1}{\mu m_{\text{H}}}\frac{{\Sigma }_{\text{g}}(R)}{\sqrt{2\pi }H(R)}.$$

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

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

Full text

††pagerange: First evidence of external disc photoevaporation in a low mass star forming region: the case of IM Lup–References††pubyear: 2016

1 introduction

Protoplanetary discs - the birthplaces of planets - are found around young stars which are themselves formed in clusters. The discs are thus externally irradiated by other cluster members, in particular by the most massive stars. Strong irradiation of discs close to O stars is well established, for example from observations of proplyds in Orion (McCaughrean & O’dell, 1996; Johnstone et al., 1998; O’Dell, 1998; Bally et al., 2000; O’Dell, 2001; Henney et al., 2002). For some time there has also been the theoretical expectation that protoplanetary discs might be significantly affected by more canonical radiation field strengths found in a cluster environment (e.g. Scally & Clarke, 2001; Adams et al., 2004; Holden et al., 2011; Facchini et al., 2016; Haworth et al., 2016). This is now being directly supported by recent observations, such as those by Kim et al. (2016) who identify proplyds irradiatied by a 3000 G0111G0 - the Habing unit - is a measure of the UV field local to our solar system and has the value $1.6\times 10^{-3}$ erg cm*-2* s*-1*. radiation field – approximately 100 times weaker than the field strengths irradiating classical proplyds (Störzer & Hollenbach, 1999). To date, however, there is no direct observational inference of externally driven mass loss from discs at lower, but more typical, radiation field strengths in the range $1<$ G0 $<1000$.

IM Lupi is a roughly Solar mass (Panić et al., 2009) young ($\sim 0.5-1$ Myr; Mawet et al., 2012) M0 star situated at a distance of $\sim 161$ pc (Gaia Collaboration et al., 2016) in the vicinity of the Lupus 2 cloud. Although CO emission likely associated with the disc is detected out to $\sim 1000$ AU, it is only detected in millimetre continuum out to about 313 AU so the gas disc is more extended than the dust (Lommen et al., 2007; Pinte et al., 2008; Panić et al., 2009; Cleeves et al., 2016). The disc is also very massive, with estimates of 0.1 and 0.17 M⊙ from Pinte et al. (2008) and Cleeves et al. (2016) respectively. The mass accretion rate is currently about $10^{-8}$ M⊙ yr*-1* (Alcala’ et al., 2016). The most recent analysis of this system by Cleeves et al. (2016) combined new 12CO, 13CO and C18O ALMA observations with a broad array of modelling resources to provide a very comprehensive chemical and radiative transfer model of IM Lup, which could describe many features of the disc very successfully. They also included the effect of external irradiation on the composition and thermal structure of the disc. Based on their modelling efforts and from geometrical arguments based on HIPPARCOS data they estimate a low UV field incident upon the disc of only about 4 G0. There was, however, a diffuse halo of low velocity CO emission about the disc that their model failed to explain. They suggested that this halo might be a remnant structure rather than being material driven out of the disc by photoevaporation. The photoevaporation interpretation was disfavoured based on the inferred low UV field and outer disk temperatures, which were well below those which had been previously considered by external photoevaporation models (Adams et al., 2004; Facchini et al., 2016). However, since this regime is previously unexplored it is difficult to conclude this with any certainty.

In this letter we use photochemical-dynamical models to investigate the external irradiation of IM Lup by the weak UV radiation field expected. We aim to determine the expected mass loss rate and flow properties and to determine whether the CO halo could be explained by such a flow. Ultimately we aim to determine whether low radiation field strengths can drive efficient mass loss and whether IM Lup offers an opportunity to observationally probe externally driven mass loss in the modest radiation regime.

2 Modelling the external photoevaporation of IM Lup

2.1 Numerical method and disc construction

We directly model the photoevaporative outflow, driven by external irradiation, using a radiation hydrodynamics and photodissociation region chemistry code torus-3dpdr, for which key relevant papers are Haworth & Harries (2012); Harries (2015); Haworth et al. (2015); Bisbas et al. (2015). This code was used to run models of externally irradiated discs in benchmark scenarios where there are semi–analytic solutions in Haworth et al. (2016) - validating the approach. The details of the method are also discussed in the latter paper.

In summary we perform calculations of the photodissociation region (PDR) chemistry in sequence with hydrodynamics using operator splitting. The PDR chemistry network is a reduced version of the umist network (McElroy et al., 2013) including 33 species and 330 reactions, and was derived such that it gives temperatures that do not differ appreciably ($\sim$10 per cent) from the much more substantial (and computationally expensive) full network. We do not include polycyclic aromatic hydrocarbons (PAHs), since although they are a key heating mechanism in PDRs, they are observed to be depleted towards discs (Geers et al., 2006; Oliveira et al., 2010). Our models will therefore yield mass loss rates lower than models that would include PAHs. Because we compute steady state flow profiles we are permitted to perform the PDR calculations relatively infrequently, as the same steady state profile will always eventually result.

Following Adams et al. (2004); Facchini et al. (2016), our models are 1D spherical (see Figure 1 of the latter paper). This is believed to be justified because the mass loss is expected to be dominated from the disc outer edge since: i) the material there is least gravitationally bound and ii) the density falls off vertically in a disc more rapidly than radially. This method also assumes that the incident (exciting) UV field approaches inwards radially and cooling line photons escape outwards radially - so every other direction is infinitely optically thick.

We employ a fixed structure for the disc itself, which acts as an inner boundary condition to the radiatively driven flow. Interior to some outer disc radius $R_{\textrm{d}}$ we do not allow the conditions to evolve over time. For these fixed disc conditions we use the parameters derived by Cleeves et al. (2016). The disc’s gas surface density profile follows that of Lynden-Bell & Pringle (1974), i.e.

[TABLE]

where Cleeves et al. (2016) find $\Sigma_{c}=25$ g cm*-2*, $R_{c}=100$ AU, and $\gamma=1$. The scale height is set by

[TABLE]

where Cleeves et al. (2016) find $\psi=1.15$ and $H_{100}=12$ AU. For the dust we assume a cross section of $\sigma_{\textrm{FUV}}=5.04\times 10^{-23}$ cm*-2*, dust to gas ratio of $d/g=10^{-4}$ and the maximum grain size $s_{\textrm{max}}=1\,\mu$m, which are all representative of the kind of dust parameters in the flow found by Facchini et al. (2016). We assume that the disc outer edge is sufficiently far from the parent star that the temperature there is only 10 K. The outer dynamical boundary condition in our models is free-outflow, no-inflow and the inner condition set by the disc properties at $R_{\textrm{d}}$ as described above. The mid-plane number density in the discs of these 1D spherical models is

[TABLE]

The radial extent of our simulation grid – 1017 cm – was chosen such that the critical radius in the flow (Facchini et al., 2016) is captured, which we check using the approach detailed in section 5.3.2 of Haworth et al. (2016). We use an adaptive grid with a maximum number of cells of 2048 and therefore a maximum resolution of 3.25 AU. We run each model for 1 Myr, though steady state flows are established long before this.

3 Results and discussion

3.1 Disc photoevaporation and evolution

We ran a grid of photoevaporation models for different disc outer radii and incident radiation field strengths. We chose disc radii in 50 AU intervals from 350 to 800 AU and radiation field strengths of 0.5, 1, 2, 4, 8 and 16 G0. We compute the mass loss rate from our models following Adams et al. (2004):

[TABLE]

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

[TABLE]

and $H_{d}$ is the scale height at the disc outer edge $R_{d}$. We compute the average of this quantity over the entire flow (note that $\mathcal{F}$ is constant for a given disc outer radius). A summary of the mass loss rates from our grid of models is shown in Figure 1. For large discs (like IM Lup) where material at the outer edge is not so gravitationally bound, substantial mass loss rates ($\sim 10^{-8}$ M⊙ yr*-1*) can be driven even when the incident radiation field strength is very modest. Note that the 4 G0 field expected to be irradiating IM Lup and driving this mass loss is $\sim 10^{3}$ times weaker than that irradiating the proplyds observed by Kim et al. (2016) and $\sim 10^{4}$ times weaker than the proplyds in the ONC (e.g. Bally et al., 2000; O’Dell, 2001; Henney et al., 2002)

The current mass accretion rate in this system was recently computed using new X–shooter data to be $10^{-8}$ M*⊙* yr*-1* with an uncertainty of 0.35 dex by Alcala’ et al. (2016). The external photoevaporative mass loss rate for UV fields $\geq 4$ G0 is hence expected to be of order or greater than the mass accretion rate.

We fit the mass loss rate as a function of radius, which we feed into the Clarke (2007) secular evolutionary code to examine the disc evolution. The viscous timescale of this disc is of order 10 Myr, so the main evolutionary behaviour is truncation of the disc by external photoevaporation. Figure 2 shows the evolution of the disc outer edge as a function of time for different incident UV field strengths (note that for models that drop below an outer radius of 350 AU we compute additional photoevaporation models to estimate the mass loss rate at these smaller radii). In all cases the disc outer edge rapidly retreats to some stagnation radius in less than 1 Myr, after which it varies in size only slowly. The mean radius over 10 Myr as a function of incident UV field is given in Figure 3, showing strong variation for fields $<8$ G0. A key point is that because the observed CO emission is currently extended out to beyond 1000 AU, even an extremely weak UV radiation field would be expected to truncate this very rapidly. The observed CO emission therefore either has to be part of a photoevaporative flow, or part of some much denser envelope that is resilient against the effects of the incident radiation field. Because IM Lup is very young ($\leq 1$ Myr), its outer edge may still be in the process of retreating towards the stagnation radius. Another interesting point is that due to the disc’s long viscous timescale IM Lup is likely to remain unusually large at the stagnation radius (perhaps > 300 AU) for many Myr, unless some other mechanism further truncates the disc.

3.2 Could external photoevaporation explain the CO halo?

Our models imply that even in the presence of a weak UV field substantial mass loss is being induced from IM Lup by external photoevaporation; however, being 1D, they are difficult to directly compare with the real observed data of Cleeves et al. (2016). Generating synthetic observations from 1D models has the limitation that some assumption about the vertical density, temperature and compositional structure is required. Nevertheless we make an optimistic attempt at comparison. We assume that the disc (the boundary condition of the dynamical models) is hydrostatic. In the flow region we use our simulation results and assume that at a given spherical radius there is a constant density, isothermal, isochemical flow, with scale height set by assuming that $H/R$ beyond the disc outer edge is constant. We produce synthetic data cubes using the comoving frame molecular line radiative transfer components of torus, detailed in Rundle et al. (2010). These cubes are then azimuthally averaged in the same manner used to produce the results in Figure 12 of Cleeves et al. (2016).

Because our synthetic observations are based on 1D models, and there is a large array of possible parameters, we do not aim to fit the CO observations. Furthermore, given that we are comparing with 12CO our synthetic observations will be particularly limited in components of the flow that are optically thick (which vary for each model but we generally find are interior to about 800–900 AU). Rather then, we aim to demonstrate that even weak external photoevaporation is capable of producing substantial emission at large radii, such as that observed in the CO halo of IM Lup.

Figure 4 shows a collection of approximate synthetic emission profiles from our photoevaporation models, as well as emission profiles from a selection of the models from Cleeves et al. (2016). The latter models modify the incident UV field but do not permit radial dynamical evolution and thus impose the surface density profile given by equation 1. As a result the extent of the CO emission is significantly less than that observed. Conversely our external photoevaporation models do show emission comparable in extent and magnitude to the observations.

In Figure 5 we plot the radial extent of the CO emission in our models as a function of the disc outer radius, with different lines representing different incident UV fluxes. The fitted gas extent from Cleeves et al. (2016) is 1200 AU which we take to be “the extent” of IM Lup for our comparison here, though in practice the detection is marginal beyond 1000 AU. From our models the extent is the point at which the flux drops below $2$ mJy beam*-1* km s*-1*, which is the background as calculated using the average of the first and last velocity channels in the synthetic data cube. Most of our models have an extent $1000-1300$ AU. Generally the models that have extent consistent with the observations have disc outer radii which imply that the disc outer edge is still retreating.

If the observed extent were known with higher certainty, we could use it in conjunction with Figures 2 and 5 to constrain the minimum disc outer edge and hence maximum age. For example if we knew that the observed extent was 1200 AU then linearly interpolating Figure 5 would yield minimum disc outer radii of 430, 450 and 530 AU for incident UV fields of 4, 2 and 1 G0 respectively. Using our evolutionary models from Figure 2, these minimum disc outer radii would correspond to approximate maximum IM Lup ages of 0.3, 0.8 and 0.8  Myr respectively – so all would be conceivable given the uncertain 0.5–1 Myr estimate for the age of IM Lup. Future higher sensitivity observations might offer such a constraint.

4 Summary and conclusions

We model the external photoevaporation of the large protoplanetary disc IM Lup. This disc has a large CO “halo” that was identified in recent ALMA observations by Cleeves et al. (2016) which could not be explained by hydrostatic chemical and radiative transfer models which assumed that the surface density at large radius was an extrapolation of the form given in equation 1. We find that although the radiation field irradiating IM Lup is very weak ($<10^{4}$ times the UV field irradiating the proplyds near O stars in Orion), the disc is sufficiently large that the weakly gravitationally bound material at the disc outer edge can be efficiently photoevaporated. Specifically a 4 G0 radiation field induces mass loss of $\sim 10^{-8}$ M*⊙* yr*-1* which is comparable to the current accretion rate onto the star. Having a $\sim$10 Myr viscous timescale, the effect of this mass loss is to rapidly (<1 Myr) truncate the disc outer edge down to some stagnation radius. The stagnation radius ranges from about 600 AU for an irradiating UV field of 0.5 G0 down to about 300-350 AU for fields $8-16$ G0. In the absence of other external influences the disc only evolves slowly away from the stagnation radius over 10 Myr. Once gas from the disc cannot be delivered to the outer edge at a rate sufficient to supply the photoevaporative wind the disc is expected to shrink rapidly.

We also generated approximate synthetic observations from our models, which are able to explain the radial extent of CO emission about IM Lup. Our scenarios that are consistent with the observed extent of CO emission of IM Lup generally imply that its disc outer radius is still in the process of being truncated. More generally we demonstrate that even weak external fields can lead to significant extended emission from large discs, which hydrostatic models are unable to achieve.

Acknowledgements

TJH is funded by an Imperial College London Junior Research Fellowship. This work has also been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. LIC acknowledges the support of NASA through Hubble Fellowship grant HST-HF2-51356.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. This work was undertaken on 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. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.00694. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.