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

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⁴ times lower than the ultraviolet 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 time-scale is long (∼10 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 that can explain the CO halo observed about IM Lup out to ≥1000 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.

accretion, accretion discs, protoplanetary discs, circumstellar matter, stars: individual: IM Lup, photodissociation region (PDR)

1 INTRODUCTION

Protoplanetary discs – the birthplaces of planets – are found around young stars that 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, Hollenbach & Bally 1998; O'Dell 1998; Bally, O'Dell & McCaughrean 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, Clarke & Bisbas 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 G₀¹ 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 < G₀ < 1000.

IM Lupi is a roughly solar mass (Panić et al. 2009) young (∼0.5–1 Myr; Mawet et al. 2012) M0 star situated at a distance of ∼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 ∼1000 au, it is detected only 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⁻⁸ M_⊙ yr⁻¹ (Alcala’ et al. 2017). The most recent analysis of this system by Cleeves et al. (2016) combined new ¹²CO, ¹³CO and C¹⁸O ALMA observations with a broad array of modelling resources to provide a very comprehensive chemical and radiative transfer model of IM Lup that 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 ultraviolet (UV) field incident upon the disc of only about 4 G₀. 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 disc temperatures, which were well below those that 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 (PDR) chemistry code torus-3dpdr, for which key relevant papers are Haworth & Harries (2012), Harries (2015), Haworth et al. (2015) and 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 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 (∼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) and Facchini et al. (2016), our models are 1D spherical (see fig. 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 that acts as an inner boundary condition to the radiatively driven flow. Interior to some outer disc radius R_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.

\begin{equation} \Sigma _{\textrm{g}}(R) = \Sigma _{\rm c}\left(\frac{R}{R_{\rm c}}\right)^{-\gamma }\exp \left[-\left(\frac{R}{R_{\rm c}}\right)^{2-\gamma }\right], \end{equation}

(1)

where Cleeves et al. (2016) find Σ_c = 25 g cm⁻², R_c = 100 au and γ = 1. The scaleheight is set by

\begin{equation} H(R) = H_{100}\left(\frac{R}{100\,\textrm{au}}\right)^{\psi }, \end{equation}

(2)

where Cleeves et al. (2016) find ψ = 1.15 and H₁₀₀ = 12 au. For the dust, we assume a cross-section of σ_FUV = 5.04 × 10⁻²³ cm⁻², dust-to-gas ratio of d/g = 10⁻⁴ and the maximum grain size s_max = 1 μ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_d as described above. The mid-plane number density in the discs of these 1D spherical models is

\begin{equation} n(R) = \frac{1}{\mu m_{\mathrm{H}}} \frac{\Sigma _{\textrm{g}}(R)}{\sqrt{2\pi }H(R)}. \end{equation}

(3)

The radial extent of our simulation grid – 10¹⁷ 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 G₀. We compute the mass-loss rate from our models following Adams et al. (2004):

\begin{equation} \skew4\dot{M} = 4\pi R^2 \rho \dot{R} \mathcal {F}, \end{equation}

(4)

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

\begin{equation} \mathcal {F} = \frac{H_{\rm d}}{\sqrt{H_{\rm d}^2+R_{{\rm d}}^2}} \end{equation}

(5)

and H_d is the scaleheight 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 Fig. 1. For large discs (like IM Lup) where material at the outer edge is not so gravitationally bound, substantial mass-loss rates (∼10⁻⁸ M_⊙ yr⁻¹) can be driven even when the incident radiation field strength is very modest. Note that the 4 G₀ field expected to be irradiating IM Lup and driving this mass-loss is ∼10³ times weaker than that irradiating the proplyds observed by Kim et al. (2016) and ∼10⁴ times weaker than the proplyds in the Orion Nebula Cluster (e.g. Bally et al. 2000; O'Dell 2001; Henney et al. 2002).

Figure 1.

Log mass-loss rate as a function of disc outer radius for different incident UV fields.

The current mass-accretion rate in this system was recently computed using new X shooter data to be 10⁻⁸ M_⊙ yr⁻¹ with an uncertainty of 0.35 dex by Alcala’ et al. (2017). The external photoevaporative mass-loss rate for UV fields ≥4 G₀ 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 time-scale of this disc is of order 10 Myr, so the main evolutionary behaviour is truncation of the disc by external photoevaporation. Fig. 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 Fig. 3, showing a strong variation for fields <8 G₀. 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 has to be either 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 (≤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 time-scale, 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.

Figure 2.

The outer disc radius of IM Lup as a function of time according to our evolutionary models that include external photoevaporation.

Figure 3.

The mean radius of IM Lup over 10 Myr as a function of the incident UV field strength.

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 the scaleheight 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 as used to produce the results in fig. 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 ¹²CO, 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 them interior to about 800–900 au). Rather, 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.

Fig. 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.

Figure 4.

Azimuthally averaged emission profiles from our approximate synthetic observations, compared to the observed data points (with 1σ error bars) from Cleeves et al. (2016). The upper panel varies the disc outer radius for a radiation field of 4 G₀. The middle panel varies the incident radiation field strength upon a 450 au disc. The bottom panel shows a collection of models from Cleeves et al. (2016).

In Fig. 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⁻¹ km s⁻¹, 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.

Figure 5.

The approximate radial extent of CO emission in our models as a function of disc outer radius. The horizontal line represents the radial extent from the best-fitting model of Cleeves et al. (2016).

If the observed extent were known with higher certainty, we could use it in conjunction with Figs 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 Fig. 5 would yield minimum disc outer radii of 430, 450 and 530 au for incident UV fields of 4, 2 and 1 G₀, respectively. Using our evolutionary models from Fig. 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 that 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⁴ 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 G₀ radiation field induces mass-loss of ∼10⁻⁸ M_⊙ yr⁻¹, which is comparable to the current accretion rate on to the star. Having an ∼10 Myr viscous time-scale, 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 G₀ down to about 300–350 au for fields 8–16 G₀. In the absence of other external influences, the disc evolves only 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 that 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.

Acknowledgments

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 the Associated Universities, Inc.

1

G₀ - the Habing unit - is a measure of the UV field local to our solar system and has the value 1.6 × 10⁻³ erg cm⁻² s⁻¹.