Hydrodynamic ablation of protoplanetary disks via supernovae

TL;DR

This study uses 3D simulations to explore how nearby supernovae can strip material from protoplanetary disks, revealing the effects of inclination, mass, and interaction duration on disk mass loss and asymmetry.

Contribution

It provides the first detailed 3D simulation analysis of supernova-induced ablation of protoplanetary disks considering various disk orientations and masses.

Findings

Mass-loss rates range from 1e-7 to 1e-6 Msol/yr, peaking at 1e-5 during rapid stripping.

Inclination angle minimally affects mass loss unless the disk is nearly edge-on.

Surviving disks contain up to 5e-6 of SN material, insufficient to explain early solar system radionuclides.

Abstract

We present three-dimensional simulations of a protoplanetary disk subject to the effect of a nearby (0.3pc distant) supernova, using a time-dependent flow from a one dimensional numerical model of the supernova remnant (SNR), in addition to constant peak ram pressure simulations. Simulations are performed for a variety of disk masses and inclination angles. We find disk mass-loss rates that are typically 1e-7 to 1e-6 Msol/yr (but peak near 1e-5 Msol/yr during the "instantaneous" stripping phase) and are sustained for around 200 yr. Inclination angle has little effect on the mass loss unless the disk is close to edge-on. Inclined disks also strip asymmetrically with the trailing edge ablating more easily. Since the interaction lasts less than one outer rotation period, there is not enough time for the disk to restore its symmetry, leaving the disk asymmetrical after the flow has passed. Of the low-mass disks considered, only the edge-on disk is able to survive interaction with the SNR (with 50% of its initial mass remaining). At the end of the simulations, disks that survive contain fractional masses of SN material up to 5e-6. This is too low to explain the abundance of short-lived radionuclides in the early solar system, but a larger disk and the inclusion of radiative cooling might allow the disk to capture a higher fraction of SN material.

Tables

Table 1. Table 1: SNR parameters

Variable	Value
$M_{ej}$	$20$M☉
$E_{ej}$	${10}^{44}$J
${\rho }_{amb}$	$2.34\times {10}^{-24}$g cm-3
$T_{amb}$	${10}^4$K

Table 2. Table 2: Disk model parameters

Variable	Value(s)
$R_0$	$1$au
$T_0$	$400$K
${\rho }_0$	$\{\begin{array}{c}3.5\times 10{}^{-13}\text{g }\text{cm}{}^{-3}\\ 3.5\times 10{}^{-12}\text{g }\text{cm}{}^{-3}\\ 3.5\times 10{}^{-11}\text{g }\text{cm}{}^{-3}\end{array}$
$M_{\ast }$	$1$M☉

Table 3. Table 3: Summary of the simulations. The disk angle ( i 𝑖 i ) is defined as the difference in angle between the angular momentum vector of the disk and the vector of the flow direction, making 0 ∘ superscript 0 0^{\circ} face-on and 90 ∘ superscript 90 90^{\circ} edge-on. M J subscript 𝑀 𝐽 M_{J} is the mass of Jupiter.

Simulation	${\rho }_0$/ g cm-3	$M_{disk}$/$M_J$	$i$	Dynamic flow?
const00high	$3.5\times {10}^{-11}$	8.22	$0^{\circ }$	n
const00med	$3.5\times {10}^{-12}$	0.822	$0^{\circ }$	n
const00low	$3.5\times {10}^{-13}$	0.0822	$0^{\circ }$	n
const45low	$3.5\times {10}^{-13}$	0.0822	${45}^{\circ }$	n
const90low	$3.5\times {10}^{-13}$	0.0822	${90}^{\circ }$	n
dyn00low	$3.5\times {10}^{-13}$	0.0822	$0^{\circ }$	y
dyn45low	$3.5\times {10}^{-13}$	0.0822	${45}^{\circ }$	y
dyn90low	$3.5\times {10}^{-13}$	0.0822	${90}^{\circ }$	y

Equations

$$X_{i,after}=\left(1-f_{enr}\right)X_{i,before}+f_{enr}{X}_{i,supernova}.$$

$$f_{enr}=\frac{X_{i,after}}{X_{i,supernova}}.$$

$$T(R)=T_0(R/R_0{)}^{-1}.$$

$${\rho }_{mid}(R)={\rho }_0(R/R_0{)}^{-3/2},$$

$$\rho (\text{mn@boldsymbol}r)={\rho }_{mid}\exp \left(\frac{G{M}_{\ast }}{c_s^2}\left[\frac{1}{r}-\frac{1}{R}\right]\right),$$

$$\Omega (\text{mn@boldsymbol}r)={\Omega }_k(R)\sqrt{\frac{R}{r}-\frac{5}{2}{\left(\frac{c_s}{{\Omega }_k(R)R}\right)}^2},$$

$$P_{grav}=\frac{G{M}_{\ast }\sigma }{R^2},$$

$$\dot{M}\approx 10^{-7}{M}_{\text{sun}}{\mathrm{yr}}^{-1}{\left(\frac{r_{disk}}{100\mathrm{au}}\right)}^{1-1.5},$$

$$P_{ram}=\frac{\dot{M}}{4\pi r^2\nu },$$

Full text

1 Introduction

It is widely believed that massive stars can trigger the formation of lower mass stars in their surrounding regions, leading to frequent associations of massive stars, and young low mass stars. For instance, there are $\sim 2000$ low-mass stars within 2pc of the centre of the Trapezium cluster (Hillenbrand & Hartmann, 1998), and stars with circumstellar disks (O’dell et al., 1993; McCaughrean & O’dell, 1996; Bally et al., 1998) within a few tenths of a parsec of the central star, $\theta^{1}$ Ori C. As massive stars are often the source of strong winds, high ionising radiation flux, and eventual supernovae, these regions are potentially hostile to protoplanetary disks and the subsequent formation of planets. The interactions massive stars have on other stars and their disks can be broadly categorised into gravitational, radiation or ablation.

When two stars pass one another, the associated disks become perturbed by the gravitational interaction. If this interaction is strong enough, material from the disk can become unbound (Clarke & Pringle, 1993). Scally & Clarke (2001) use N-body simulations to model the dynamics of the stars in the Orion Nebula. They find that only a small fraction of the stars interact closely on the expected disk dispersal timescales, and so close encounters are not a major contributor to protoplanetary disk dispersal.

Direct observational evidence of the effects of radiation on disks is seen in the form of the proplyds in the Orion Nebula (e.g. O’dell & Wen (1994); McCaughrean (2001)). The radiation from $\theta^{1}$ Ori C causes a photoevaporation flow from the disks which collides with the oncoming stellar wind. This results in a number of cometary-shaped ionization fronts associated with stars around $\theta^{1}$ Ori C. Estimates put the mass-loss rate of the proplyds at $\sim 4\times 10^{-7}M_{\sun}\,$yr*-1* (Henney & O’Dell, 1999). If sustained over the lifetime of the cluster this would mean that the stars would have had unrealistically massive disks (greater than the star’s mass) at some point in the past. However, $\theta^{1}$ Ori C may have switched on only $\sim 10^{5}\,$years ago (Scally & Clarke, 2001), in which case it is unlikely that the combination of photoevaporation and viscous accretion will have significantly depleted the disks (Garcia, 2011). This suggests that the Orion Nebula is in a relativity short-lived stage of its evolution, which is consistent with the lack of proplyd-like objects in other clusters (e.g. Stapelfeldt et al. (1997); Balog et al. (2006)).

There is now a lot of observational evidence that massive stars influence the disks of neighbouring stars. Cometary tails are seen to extend from protoplanetary disks in Hubble Space Telescope images of the Orion Nebula (e.g. Smith et al., 2005; Ricci et al., 2008), while sub-millimeter observations of the Trapezium Cluster reveal that disks within a few tenths of a parsec of $\theta^{1}$ Ori C are truncated (Mann & Williams, 2009, 2010; Mann et al., 2014). However, the degree to which massive stars affect nearby disks remains uncertain. For instance, Mann et al. (2015) claim that there is no evidence for disk truncation in the Flame Nebula (NGC 2024), though it is somewhat younger and hosts a less massive star compared to the Trapezium cluster. In addition, Richert et al. (2015) argue that there is no evidence for the depletion of the inner disks around pre-main-sequence stars in the vicinity of O-type stars, even very luminous O2-O5 stars. In contrast, Povich et al. (2016) claim that rapid disk evolution on $<1\,$Myr timescales has produced a significant population of intermediate-mass pre-main-sequence stars lacking inner dust disks.

Likewise, significant differences in the mass-loss rate and lifetime of disks exist in theoretical studies of their photoevaporation. Some studies claim lifetimes of only $10^{5}\,$yrs (e.g. Johnstone et al., 1998; Störzer & Hollenbach, 1999), substantially shorter than the $10^{6}\,$-$\,10^{7}\,$yrs lifetime of disks subject only to viscous accretion (Haisch et al., 2001; Hernández et al., 2007, and references therein). In contrast, a recent model combining viscous accretion and external photoevaporation predicts that disk lifetimes are shortened by only a factor of a few (Anderson et al., 2013). In addition, if the outer parts of the disk can be totally removed, the inner disk may rapidly disappear if it is starved of material needed to offset accretion onto the star (Anderson et al., 2013). Thus the significance of external photoevaporation on disk lifetimes remains highly uncertain.

Supernova remnants (SNRs) provide another source of agitation for protoplanetary disks. There is evidence in $\pi$ Sco of circumstellar disks similar to those in the Orion nebula (Bertoldi & Jenkins, 1992) and that a supernova occurred nearby $\sim 10^{6}\,$yrs ago (de Geus, 1992). Lada & Lada (2003) compiled a catalog of embedded clusters and find that the vast majority ($70$ - $90\%$) of stars form in clusters of $>100$ members, and the majority of those ($\sim 75\%$) are in clusters containing stars massive enough to give rise to supernovae (Hester & Desch, 2005).

Analytic estimates for the ablation of protoplanetary disks due to a SNR were presented by Chevalier (2000). He found that for typical disk and supernova parameters, partial stripping of the disk can occur, but typically not its complete disruption. This suggests that although the disks are affected, they can survive such events.

The interaction of a SNR with a protoplanetary disk is also of interest from the point of view of injecting short-lived radionuclides (SLRs) into the early solar system (e.g. Cameron & Truran (1977); Jacobsen (2005); Looney et al. (2006)). Analyses of meteorites suggest that the early solar system had a higher abundance of SLRs than the average interstellar medium (ISM) value (see e.g. Diehl et al., 2006, for ${}^{26}Al$). This means the SLRs must have been injected into the solar system at some point during its history, and due to the short half-life of 60Fe (Rugel et al., 2009, $2.6\,$Myr,)) this constrains the injection event to a time close to the solar system’s formation. Given that supernovae are abundant sources of nucleosynthesis and they are often in close proximity to young, low mass stars, they are seen as a good candidate for the source of SLRs in the early solar system. While some SLRs can be explained by production internal to the solar system, 60Fe is practically impossible to produce locally (Lee et al., 1998; Gounelle et al., 2006).

If the formation of the solar system was triggered by a SNR this could provide a source of SLRs (see e.g. Foster & Boss (1996); Boss & Keiser (2010, 2012, 2014, 2015)). While this is a plausible model, it puts extra constraints on the timing of the formation of the disk. This is estimated to take $\sim 0.1-10\,$Myr (see e.g. Shu et al., 1993; Vanhala & Cameron, 1998; Tassis & Mouschovias, 2004). This is on top of the time needed to form small bodies in the solar system, which is $\sim 1-10\,$Myr (Hartmann, 2005). If SLRs are injected during the molecular cloud stage, disk formation must occur quickly such that enough SLRs exist in the early solar system.

In the current paper we investigate injection into an already formed disk as an alternative to this model. Injection would have to occur in the relativity early stages of the disk’s lifetime, before small bodies begin to form. It is currently unknown at what point material can be injected and homogenised in order to be captured by forming meteorites. As all these processes are likely to occur simultaneously, simulations that track both injection and the subsequent formation into meteorites would go a long way to answering this question. However, this is beyond the scope of the current paper.

The abundance of SLRs is typically parametrised by the ratios of 60Fe and 26Al to their stable isotopes. Quitté et al. (2005) and Tachibana et al. (2006) find a 60Fe/56Fe ratio of $\sim 3-7\times 10^{-7}$. However, the technique used (secondary ion mass spectrometry (SIMS)) has been shown to give a significant positive bias (Ogliore et al., 2011). Subsequent studies using a different technique (Tang & Dauphas, 2012) show a lower ratio of $5\times 10^{-8}$, although Mishra & Goswami (2014) use an improved SIMS method and find a ratio of $7\times 10^{-7}$. Therefore, the true ratio is still unclear. However, if the lower estimates for the amount of 60Fe are correct, then the 60Fe/26Al ratio becomes much lower than that generated by supernovae (Gounelle & Meynet, 2012). A possible resolution to this problem is suggested by Goodson et al. (2016), who note that different sized dust grains are injected at different rates and if they preferentially carry different SLRs this could have a significant effect on the 60Fe/26Al ratio.

The amount of 60Fe produced by a supernova varies depending on the mass and metallicity of the exploding star. Calculations by Woosley & Weaver (1995) put the 60Fe mass fraction at $\sim 2\times 10^{-6}$ for a $20M_{\sun}$ supernova. Using initial solar system abundances from Lodders (2003), Looney et al. (2006) calculate a meteoritic 60Fe mass fraction of $1.10\pm 0.38\times 10^{-9}$.

We can define an enrichment fraction, $f_{enr}$, as the mass of supernova material mixed into the disk as a fraction of the disk mass. Then the change in mass fraction, $X$, of a given species, $i$, after enrichment can be calculated as:

[TABLE]

If the initial solar system has no significant 60Fe, then $f_{enr}$ becomes simply:

[TABLE]

Using $X_{i,after}=1.10\times 10^{-9}$ and $X_{i,supernova}=2\times 10^{-6}$ from above this gives $f_{enr}=5.5\times 10^{-4}$. Doing the same calculation for 26Al ($X_{i,after}=3.77\times 10^{-9}$ and $X_{i,supernova}=3\times 10^{-6}$) gives $f_{enr}=1.3\times 10^{-3}$.

The assumption of zero $X_{i,before}$ is a contested one. It has been suggested (e.g. Gaidos et al., 2009; Young, 2016) that successive supernovae and massive stellar winds can allow a molecular cloud to self-enrich, such that stellar systems form with enriched material already embedded. This scenario is explored in simulations by Vasileiadis et al. (2013) and Kuffmeier et al. (2016). Using a turbulent periodic box, the evolution of a massive ($\sim 10^{5}M_{\sun}$) star forming region is followed over $20\,$Myr. It is found that supernova feedback enriches the star forming gas to a level that is largely consistent with that of the early solar system, although the authors rely on numerical diffusion to mix the SN ejecta and cold gas. This model also remains somewhat controversial as it requires an age spread of star formation within the same cloud, which is at odds with observations (Elmegreen, 2000; Jeffries et al., 2011; Soderblom et al., 2014) and simulations (Bonnell et al., 2003; Dale et al., 2012). While a non zero value for $X_{i,before}$ doesn’t necessarily change the scenario investigated in this work, it would increase the mass fractions we see.

Variations in the supernova mass can also produce quite large differences in the mass fraction, not to mention the inherent uncertainties in supernova modelling. For example, Limongi & Chieffi (2006) calculate a 60Fe mass fraction for a $20M_{\sun}$ supernova at $7.8\times 10^{-7}$, $2.5\times$ less than that calculated by Woosley & Weaver (1995). Their 60Fe mass fractions for the next nearest supernova masses ($17M_{\sun}$ and $25M_{\sun}$) vary from the $20M_{\sun}$ mass fraction by a factor of $\sim 1.5$. With this in mind $f_{enr}$ is only constrained to a factor of two or three at best.

Hydrodynamical simulations of the interaction of a SNR with a circumstellar disk were performed by Ouellette et al. (2007). They found that only a small fraction ($\sim 1\%$) of the disk mass is removed, in contrast to the $\sim 13\%$ expected from the analytical prediction of Chevalier (2000). They attribute this to the cushioning and deflecting effect of the bow shock in addition to the compression of the disk further into the gravitational well of the central star. However their simulations are limited to two dimensions and face-on impacts. While an inclined disk might be expected to be disadvantageous for ablation, the geometry of the bowshock that shields the disk is likely to be substantially different and the impacting flow acts with the rotation of the disk on one side, reducing the threshold for ablation. This scenario therefore needs to be investigated. Ouellette et al. also found that only $\sim 1\%$ of the passing ejecta is captured by the disk. This is not enough to explain the observed abundances, although they suggest that dust grains might be a more efficient mechanism for injecting SLRs into the disk.

To investigate this possibility, Ouellette et al. (2010) incorporated dust grains into their simulations. They found that about $70\%$ of dust grains larger than $0.4\,\mu\,$m are injected into the disk. This could potentially be enough to explain the observed abundances, although it relies on higher than observed dust condensation in the ejecta and the ejecta to be clumpy with the disk being hit by a high density clump of ejecta. They estimate the probability of these conditions being met at $0.1-1\%$ and that the solar system may be atypical in this respect.

The mixing effects of clumpy supernova ejecta interacting with a molecular cloud were investigated by Pan et al. (2012). When the mixing is efficient the enrichment fraction, defined as the ratio of ejecta mass to cloud mass, is $\sim 10^{-4}$. This is approximately the correct level, although as the SLRs are injected into the molecular cloud, and it takes of order a Myr to form stars and disks, the abundance of 60Fe could be significantly reduced in this time.

The formation and subsequent evolution of a star-disk system under the effects of a supersonic wind were investigated by Li et al. (2014). They found that a disk can indeed form and survive in these conditions: a $10^{-3}M_{\sun}$ disk remains after being exposed to the wind for $0.7\,$Myr. However, their disk radius after formation is $\sim 1000\,$au, whereas disks near massive stars are typically photoevaporated down to tens of au (Johnstone et al., 1998). The flow speed is also much slower than that adopted by Ouellette et al. (2007), which in turn reduces the ram pressure significantly. Therefore, while their work was useful for investigating triggered star formation, the region of parameter space that they explored is unrealistic for disk ablation via nearby supernova.

Recently Goodson et al. (2016) performed three dimensional simulations of the interaction of a SNR with a large ($8.8\,$pc), clumpy molecular cloud, including dust grains. They found that the majority of large dust grains are injected into the molecular cloud, within $0.1\,$Myr of the supernova explosion. They note that if 60Fe and 26Al preferentially condense onto different sized grains this could explain the discrepancy in the ratio between the two SLRs.

This paper aims to improve and expand on the 2D hydrodynamical calculations of Ouellette et al. (2007) in a number of ways. Firstly, by doing fully 3d simulations, more complex dynamics can be studied, and the angle between the disk and the flow can be varied. Secondly, we consider a range of disk masses ($1\times$, $0.1\times$ and $0.01\times$ their canonical disk mass) as disk density is the main factor in determining the extent to which a disk is ablated. By considering a number of disk masses spanning a range of midplane densities, we can investigate how different disks react and at what point the disk becomes significantly affected. Finally, we adopt time-dependent flow properties from a 1D simulation of the SNR. This provides more accurate flow conditions past the disk, compared to the analytical approximation used by Ouellette et al. (2007), particularly in the early stages of the SNR’s expansion. We aim to better determine the nature of the interaction, the mass-stripping rate of the disk, and the injection rate of SLRs into the disk.

2 Models

2.1 Overview

A stellar disk is simulated in three dimensions, with mass injected onto the grid to simulate the impact of a SNR. The strength of the flow is dependent on the size and distance of the supernova and is time dependent. To calculate this precisely, a simple, one dimensional simulation of a spherical SNR is performed. The evolution of the density, velocity and temperature are recorded at the desired distance from the supernova, which can then be used by the main simulation to control the properties of the wind. This is in many ways simpler than performing an analytical approximation (e.g. Matzner & McKee (1999)), and has the advantage of reproducing all the features of a realistic SNR, particularly in the early expansion.

2.2 The supernova remnant

The supernova has an ejecta mass of $20\,$M*☉* and an ejecta energy of $10^{51}\,$erg (following Ouellette et al. (2007) to represent a typical type II supernova), and is initially confined to a radius of $2000\,$au. The ambient medium is set to a density of $2.34\times 10^{-24}\,$g cm*-3* and a temperature of $10^{4}\,$K. The simulation parameters are noted in Table 1. For all simulations, the mean molecular weight, $\mu=2.4$, and $\gamma=5/3$. The supernova calculations were performed in spherical symmetry on a uniform one-dimensional grid, extending to $r=10^{5}\,$au with 3200 cells. The same calculations were done at half and quarter resolution with identical results.

Figure 1 shows how the fluid variables change over time at a point $0.3\,$pc from the origin of the SNR. $0.3\,$pc is chosen to aid comparison with Ouellette et al. (2007) and is consistent with observations of the Orion Trapezium Cluster. These values will determine the properties of the flow past the circumstellar disk in the 3d simulation. The supernova ejecta takes $55.4\,$yr to reach $0.3\,$pc, which is insignificant compared to the half-lives of the SLRs so there is no significant decay during the travel time. The SNR is still in its free expansion stage at this point. The strength of the flow is parametrised by the ram pressure, $P_{ram}=\rho\nu^{2}$, and is shown in Figure 2. Also shown for comparison is the analytical ram pressure curve used by Ouellette et al. (2007), which is only a close approximation to the numerical results at late times. The interaction is strongest in the early stages, where the approximation is poor.

2.3 The circumstellar disk

The disk model is adapted from Gressel et al. (2013). The temperature is defined to be constant across z, and to be inversely proportional to cylindrical radius, R,

[TABLE]

This gives the disk a constant opening angle. A constant opening angle is typically desirable at it allows the disk edge to align with cell boundaries in a spherical-polar domain. While this is not a concern for these simulations, it aids comparison with works in the literature. Defining the midplane density as

[TABLE]

and enforcing hydrostatic equilibrium vertically defines the three dimensional density structure, which can be derived as

[TABLE]

where $c_{s}$ is the isothermal sound speed, and $r$ is the distance to the centre of the disk. Finally, the angular velocity, $\omega$, is set to achieve radial equilibrium,

[TABLE]

where $\Omega_{k}(R)=\sqrt{GM_{*}}R^{-3/2}$ is the Keplerian angular velocity. The free disk parameters are chosen such that the disk resembles that described in Ouellette et al. (2005). As they observe no significant ablation of their disk, lower disk mass simulations are also performed, which are likely to be more susceptible to ablation. The values are summarised in Table 2. The disk is truncated at an inner radius of $4\,$au and an outer radius of $40\,$au. The temperature of the disk ranges from $100\,$K to $10\,$K from the inner to outer disk boundaries. The outer rotation period is $\sim 250\,$yr and the inner rotation period is $\sim 10\,$yr. The initial ambient medium ($\rho_{amb}$ and $T_{amb}$) is set to the same as that of the supernova simulation for consistency (see Table 1). The disk was evolved in isolation and found to be stable for several outer rotation periods.

It should be noted that much bigger disks (with radii up to $\sim 1000\,$au see e.g. Bally et al. (2015)) have been observed. Such disks have a much larger surface area and the outer regions are less strongly held by the gravitational field, meaning they will lose mass to ablation more readily. However, a $40\,$au disk is used here to aid comparison with previous work. Without good statistics of the disk radius at formation it is also hard to say what is more typical and in any case it likely depends on unknown quantities such as the conditions under which the cluster formed. At larger radii, both the gravitational field strength and the disk’s surface density decrease, making the material more tenuously held and more easily ablated. For a given disk model, only the central density and stellar mass determine the stripping radius: the stripping radius does not depend on the initial disk radius. With this in mind, some insight into the behaviour of larger disks can also be gleaned from our simulations.

3 The simulations

The calculations were performed with the hierarchical adaptive mesh refinement (AMR) code, MG (Falle, 1991). The code uses the Godunov method, solving a Riemann problem at each cell interface, using piece-wise linear cell interpolation and MPI-parallelization. The Eulerian equations of hydrodynamics are solved using the second order upwind scheme described in Falle (1991). Refinement is on a cell-by-cell basis and is controlled by error estimates based on the difference between solutions on different grids, i.e. the difference between the solutions on $G_{n-1}$ and $G_{n}$ determine refinement to $G_{n+1}$. Spatial resolution is doubled on each refinement level. Further details on the AMR method can be found in Falle (2005).

All simulations are performed on a three dimensional Cartesian grid, with the disk situated at the origin, and the plane of rotation aligned with the X-Y plane. For each simulation a “flow injection region” is defined where the values of the grid cells are set to the current SNR flow properties at the beginning of each time step. This region is placed $96\,$au from the origin. The simulation region extends to $\pm 256\,$au in all directions, except in the cases where this would create a large flow injection region. For all simulations the lowest grid level has a resolution of $8\,$au, with 4 additional grid levels giving an effective resolution of $0.5\,$au. This gives $80$ cells per disk radius. Detailed convergence studies for shock-cloud interactions (Pittard & Parkin, 2016) indicate that a resolution of 32-64 cells per cloud radius is needed for signs of convergence, including in the mixing fraction. Although the scenarios are different, this provides some confidence that our simulations are of adequate resolution. They are also at a higher resolution than any previous studies of interactions of this type. A detailed resolution study is left to future work.

There are two reasons the disk is kept in the same orientation and the wind direction is changed (as opposed to vice versa). Firstly, if the plane of rotation is not aligned to the direction of the grid, instabilities can develop and cause the disk to fragment with no outside influence (e.g., Davis (1984); Hahn et al. (2010); Hopkins (2015)). Secondly, the grid is split among processors by dividing the domain along one of the grid directions. Placing the disk such that the divides split the disk across processors helps to distribute the computational load more evenly and allows the simulations to run more efficiently.

A set of eight simulations are performed. Table 3 details the differences between them. Some simulations are done with a constant flow, as this is the most straight forward to compare against analytical approximations. The parameters of the constant flow are defined by the peak ram pressure point of the SNR (see Figures 1 and 2) and occurs $55.8\,$yrs after the disk is first hit by the SNR. The peak is taken as it represents the worst case scenario for the survival of the disk. Specifically this is a density of $5.74\times 10^{-21}\,$g cm*-3*, a velocity of $2.87\times 10^{8}\,$cm s*-1* and a temperature of $2.72\times 10^{5}\,$K. The corresponding ram pressure is $4.74\times 10^{-4}\,$g cm${}^{-1}\,$s*-2*. In contrast the flow in the dynamic simulations follows the temperature, density and velocity of the calculated SNR, shown in Figure 1.

3.1 Defined quantities

In order to quantify the mass of the disk and the amount of SN ejecta that is captured by the disk, we define several quantities. The mass of the disk, $M_{disk}$, is defined as the total mass of gas that is bound to the central star’s gravitational field, i.e. the velocity of gas in a given cell is less than the escape velocity at that point in space. The simulation uses an advected scalar to track the SLR enriched SN ejecta. Disk material is given a scalar value of $0.0$ and the enriched SN ejecta is given a scalar value of $1.0$. Thus, for a given grid cell, the product of the advected scalar, the cell density and the cell volume gives the mass of enriched material in that cell. We define the mass of captured enriched material, $M_{enr}$, as the total mass of enriched material that is gravitationally bound. Thus the enrichment fraction $f_{enr}=M_{enr}/M_{disk}$.

4 Analytical Approximations

Material in the disk will be disrupted if the ram pressure of the wind exceeds the gravitational force per unit area, $P_{grav}$. This can be estimated as (Chevalier, 2000):

[TABLE]

where $\sigma$ is the surface density of the disk. Figure 3 shows the gravitational pressure as a function of radius for the three disk masses considered. By integrating surface density from the centre of the disk to the point at which the gravitational pressure drops below the peak ram pressure of the SNR, an estimate for the extent of the instantaneous stripping can be obtained.

Figure 4 shows the cumulative integrated mass for the disk. If all the material at a gravitational pressure less than the peak ram pressure of the supernova is stripped, then the low mass disk is left with $10.7\%$ of its initial mass, the medium density disk retains $65.1\%$ and the high mass disk retains $99.95\%$.

In their investigation of the effects of inclination on the ablation of disk galaxies, Roediger & Brüggen (2006) provide an argument for why stripping should be independent of inclination angle for small angles. Assuming the gas disk to be infinitely thin, the force due to ram pressure on a surface element $dA$ is $\rho_{wind}\nu_{wind}^{2}\mathrm{cos}(i)$. The ram pressure is effectively reduced by a factor of $\mathrm{cos}(i)$. As the radial gravitational force is balanced with the centrifugal force only the force perpendicular to the plane of the disk contributes. This is also reduced by a factor of $\mathrm{cos}(i)$. This means that the criteria for material to be stripped from the disk is independent of inclination angle. For highly inclined disks, the assumption of an infinitely thin disk will break down. While galactic disks are in a very different area of parameter space to stellar disks, these arguments are equally valid for either case.

5 Results

Figures 5-7 show snapshots of the evolution of the disk with two dimensional slices though the three dimensional grid. Figure 5 shows the evolution of disks of different masses subject to a constant flow. As expected the more massive disks are more resilient to stripping. The low and medium mass disks are significantly deformed by the flow. For the low mass disk this breaks up the disk and forms a turbulent tail. However, in the medium mass case the deformed disk remains bound.

Figures 6 and 7 show the effect of inclination angle for the constant flow and dynamic flow cases respectively. For inclined disks the stripping is asymmetrical. For disk inclinations of $45^{\circ}$ the leading edge fragments and strips from the disk first, as the bow shock partially shields the trailing edge of the disk. As the disk continues to evolve the trailing edge is stripped and the disk becomes more symmetrical again. As the $45^{\circ}$ disk evolves, the disk is heavily stripped and eventually destroyed. At an inclination of $90^{\circ}$ the side of the disk rotating in the direction of the flow is stripped more heavily at first than the side rotating against the flow. Mass is stripped from the $90^{\circ}$ disk much more slowly than for disks with lower inclinations, due mainly to the lower cross-sectional area of the disk to the flow.

The morphology of the constant wind and dynamic wind cases is broadly similar in all cases. The main distinguishing feature is that as the pressure of the dynamic wind starts to decrease, the tail becomes much wider.

The amount of mass bound to the star’s gravitational well is shown as a function of time in Figures 8-10. Simulations with a constant flow exhibit an initial rapid (“instantaneous”) stripping down to the level where the ram pressure is balanced by the gravitational pressure (where the disk mass plateaus), followed by a slower continual stripping caused primarily by the Kelvin-Helmholtz instability. Figure 8 shows this for the three different disk masses simulated, together with the analytical approximations. For the low mass disk, the analytical prediction fits the position of the plateau very well. The medium mass disk begins to plateau but the continual stripping starts to dominate before it does so, pulling the mass below the analytical approximation. The high mass disk loses virtuallly no mass due to instantaneous stripping and effectively starts in the continuous stripping phase. At the end of the simulation both the medium and high mass disks are losing mass at a rate of $\sim 10^{-6}\,$M${}_{\sun}\,$yr*-1*. During continuous stripping (between $t=40\,$yr and $t=70\,$yr), the mass-loss rate of the low mass disk is $\sim 10^{-7}\,$M${}_{\sun}\,$yr*-1*. While stripping is still occurring at the end of simulation, the constant ram pressure simulations are not continued beyond this point as the longer the simulation runs, the more unrealistic the assumption of peak ram pressure is. Running the simulations for $\approx 150-250$yrs provides enough data to establish trends and for comparison to other works.

As the inclination angle of the disk is increased, the stripping is generally slower (Figures 9 and 10). However, the difference between inclination angles of $0^{\circ}$ and $45^{\circ}$ is relatively small. Not only do they both plateau at about the same level, but they do so at about the same time (see Figures 9 and 10). This is because the flow in this instance is strong enough to deform the disk such that the initial inclination angle is no longer relevant. At high inclination angles a different behaviour is observed, and the disk survives significantly longer at an inclination of $90^{\circ}$. The same general trend can be seen in Figure 10 for the dynamic flow. Interestingly, for the face-on flow, the mass plateaus before the peak ram pressure has been reached, indicating that the history of the flow is important in shaping the disk and determining whether it is susceptible to ablation. The low mass disk only survives the dynamic flow when placed edge-on to the flow, retaining nearly $60\%$ of its mass in this case (Figure 10). The edge on disk is losing $\sim 2\times 10^{-7}\,$M${}_{\sun}\,$yr*-1* at the end of the simulation. If this rate were to be sustained the disk would be destroyed in $\sim 300$yr. However, the ram pressure decreases by two orders of magnitude in this time, so there will be some but not complete stripping during this period. The disk is therefore expected to survive significantly longer than we have simulated.

Figure 11 shows the z-direction mass flux of the disk material at the end of simulation const00high. To show the regions where the disk is losing mass the unbound material is shown in the right panel. Mass is lost from two main regions: along the sides of the deformed disk and also though the central hole. While the mass flux in both these areas is similar the central region spans a much smaller area. Figure 12 shows the same mass flux as Figure 11 now multiplied by a factor of $2\pi r$ to account for the differing contributions to the total mass flux. It is clear that the majority of the mass is lost from the edge of the disk. Ablation from the centre of the disk accounts for approximately $30\%$ of the mass being ablated at this point in the disk’s evolution.

Figure 13 shows the amount of SN ejecta that becomes bound to the disk, tracked using an advected scalar. The disks in simulations const00med, const00high and dyn90low are the only disks to survive to the end of simulation and the only ones that could be considered to be enriched by SN ejecta, although all simulations are shown for completeness. The mass of SN ejecta in the disk is $\sim 10^{-9}\,$M*☉* for the two face-on disks and $\sim 10^{-10}\,$M*☉* for the edge-on disk.

Figures 14-16 show the distribution of bound material for the three disks that survive until the end of the simulation. The right panels in these figures show the value of the advected scalar used to trace supernova ejecta: a value of 1.0 indicates pure ejecta, and a value of 0.0 indicates that no supernova ejecta is present. In all 3 figures the bound supernova ejecta is predominantly found on the upstream face (or edge) of the disk. For the edge-on disk the SN ejecta which becomes bound enters an orbit around the edge of the disk. The ejecta concentration on the leading surface is of order 1 part in $10^{4}$ (in some small localised low-density regions the concentration is at the 1 per cent level). Inside the disk the concentration drops to typically 1 part in $10^{5}-10^{6}$. Thus we see some ejecta material being deposited onto the surface of the disk, but little mixing into the disk interior on the timescales of our simulations. While our simulations do not run long enough to capture the subsequent mixing of this material (nor are turbulent motions within the disk properly resolved), Boss (2013) provide good estimates of the timescales over which mixing is likely to happen, finding that homogenization in a $40$AU disk occurs on timescales of $10^{3}$ to $10^{4}$yr. No significant decay of SLRs would be expected over this time frame.

6 Discussion

6.1 Comparison to previous works

The main point of comparison is with Ouellette et al. (2005). Our face-on high density disk can be considered equivalent to their canonical disk. Good agreement is obtained in that virtually no stripping is seen in either case. However they see less than $0.1\%$ mass lost over $2000\,$yr whereas we see $\sim 2.5\%$ lost over $140\,$yr. This difference can be attributed to the fact that our high density disk is only simulated for the constant ram pressure flow where the long term ablation is much stronger, whereas Ouellette et al. (2005) used a time varying flow. Our simulations also differ in that we observe consistency with the analytical predictions of Chevalier (2000), whereas they state that the bow shock shields the disk causing it to be able to survive higher levels of ram pressure.

Li et al. (2014) study the effect of triggered star formation, allowing a disk to form and evolve under the effects of a wind. Their simulations are in a very different region of parameter space to ours: the ram pressure of their wind is $\sim 7$ orders of magnitude lower, so they are observing the interaction on a timescale of Myrs and their disk after formation has a radius of $\sim 1000$au, with a resolution of $23$au. However, they do note that at the end of their simulation the disk radius is much less than that predicted using the method of Chevalier (2000).

As there have been no studies specifically of the ablation of inclined stellar disks, the only comparison that can be made is with galactic disks. Roediger & Brüggen (2006) and Jáchym et al. (2009) both specifically investigated the effect of inclination on the ram pressure stripping of galactic disks. We observe the same general dependence on inclination angle, where stripping is similar for low inclination angles and is only strongly impeded for close to edge-on orientations. We also observe the same asymmetry when disks are not face on. Jáchym et al. (2009) argues that the parallel side is stripped more easily as the wind works in the direction of rotation to push the material off the disk, whereas the material on the anti-parallel side must first be slowed to zero velocity and re-accelerated in the opposite direction if it is to be stripped on that side on the disk (otherwise it will continue to orbit round to the other side of the galaxy where it may then be stripped more easily). Roediger & Brüggen (2006) note that after an outer rotation period the disc becomes symmetrical again, as the entire disk has experienced the side with stronger stripping.

For galactic disks the stripping occurs over a relatively long period of time, with the disk experiencing several rotation periods while being ablated. However, the nature of the stripping depends on how the galaxy is orbiting within the cluster. A galaxy on a circular orbit will experience constant ram pressure throughout its lifetime, whereas a galaxy falling radially will experience a shorter peak of ram pressure. Roediger & Brüggen (2007) calculate the dynamic ram pressure for several different orbital paths a galaxy might take. Defining a significant ram pressure as one that would be expected to strip at least half of the disk mass if maintained constantly (this corresponds to what Roediger & Brüggen (2007) refer to as a "medium" ram pressure and is $\sim 10^{-11}\,$g cm${}^{-1}\,$s*-2* in their case), the galactic disk experiences significant ram pressure for approximately $0.75$ to $2$ times the outer rotation period during each orbit.

For the supernova and stellar disk interaction examined here the ram pressure is only significant for $\sim 1/3$ of the outer rotation period. This means it is possible for the stellar disk to maintain asymmetry after the point where stripping is significant, although it is not clear where the boundary lies, either in terms of ram pressure or rotation periods, in order for this to occur. Also note that this effect requires the flow to impact the disk edge-on or near edge-on for extended periods, whereas a galaxy is likely to experience stripping from a number of different angles at different points in its orbit.

6.2 Internal stripping

Our simulations show that a significant amount of mass is lost though the central hole of the disk during the continuous stripping phase. Theoretical considerations show that large planets can cause gaps to form in the disk (Takeuchi et al., 1996) and observations of protoplanetary disks have since detected the presence of such gaps (Calvet et al., 2002; Andrews et al., 2016). If these gaps are large and deep enough they could provide additional channels for ablation to take place.

Similarly, in disk galaxy simulations cooling causes the disk to fragment and allows the wind to flow though areas of low density within the galaxy (Tonnesen & Bryan, 2009). This causes the instantaneous stripping to proceed faster but to the same radial extent compared to the adiabatic, non-fragmented disk. The differences in long term stripping were not investigated but it is likely to be accelerated in the same way by the additional surface area, and will almost certainly effect the morphology of the flow. However, cooling is not significant in the bowshock surrounding the disk in our simulations. This may change if the SNR shock itself became strongly cooling (due to a greater age at the time of interaction, and/or higher circumestellar densities).

6.3 Continuous stripping

A number of different authors (Nulsen, 1982; Hartquist et al., 1986; Arthur & Lizano, 1997) have provided estimates for the mass-loss rate of a pressure supported globule under the effect of a wind. All their estimates are based on the mass-loss rate of the globule being, due to conservation of momentum, approximately equal to the mass flow rate though the cross-sectional area of the object being ablated, i.e., $\dot{M}=\pi r^{2}\rho_{wind}\nu_{wind}$. Calculating this value for the constant wind used for our simulations gives $\dot{M}=3\times 10^{-8}\,$M${}_{\sun}\,$yr*-1* which is around 100 times lower than the value seen in the simulations ($\sim 10^{-6}\,$M${}_{\sun}\,$yr*-1*). The flaw in this argument when applied to disks is that each parcel of gas already contains a large amount of momentum due to its rotation. For protoplanetary disks the pressure support is relatively weak so the velocity of the gas is close to Keplerian and thus the escape velocity. This means that the oncoming wind needs to transfer a smaller proportion of its momentum to the gas in the disk so that it becomes unbound. In practice this interaction is fairly complex, particularly for non-face-on orientations and there has been no simple analytical formulae proposed, but it is clear that it is not well described by those for simpler, non-rotating objects.

Photoevaporation by a nearby massive star can also play a role in disk dispersal, and is thought to be the cause of the proplyd objects in the Orion nebula (Henney & Arthur, 1997). Störzer & Hollenbach (1999) model the effect of an external source of ultraviolet radiation on circumstellar disks, finding good agreement with their models and observations of the Orion proplyds. They provide a crude fit to their results for the mass-loss rate:

[TABLE]

at a distance of $0.2$pc from a $\theta^{1}$ Ori C-like star. Richling & Yorke (2000) perform a similar analysis and find agreement with the above equation. Calculating this value for a disk of $r_{disk}=40\,$au gives $\dot{M}=3\times 10^{-8}\,$M${}_{\sun}\,$yr*-1*. This is $\sim 100$x less than the continuous stripping rate from the SNR interaction at a distance of $0.3\,$pc for models const00med and const00high. The lifetime of the const00med and const00high disks to photoevaporation, $t_{\rm life}\approx M_{\rm disk}/\dot{M}$, is thus $3\times 10^{4}$ and $3\times 10^{5}$ yrs, respectively. This indicates that photoevaporation can be an important process which may destroy disks before any nearby massive star explodes as a SN. However, not all disks will be exposed to as high a UV flux as occurs in Orion, despite still being subject to the effects of a nearby SN explosion (for instance, a B0V star has an ionizing flux which is $1.5\times$ lower than $\theta^{1}$ Ori C). In addition, more massive disks than we have considered may survive photoevaporation from Orion-like UV fluxes for several millions of years before being subject to a SN explosion.

Another possible source of ablation for stellar disks is the wind from a nearby massive star. This can last much longer than the SNR. For massive stars ($M\gtrsim 25\,$M*☉*) The strongest wind occurs during the Wolf-Rayet (WR) phase which lasts $\sim 0.3\,$Myr. For an isotropic wind, the ram pressure at a distance, $r$ is

[TABLE]

where $\dot{M}$ is the mass-loss rate of the star and $\nu$ is the velocity of the stellar wind. Typical values for the star’s wind during the WR phase are $\dot{M}=10^{-5}\,$M${}_{\sun}\,$yr*-1* and $2\times 10^{-8}\,$cm s*-1*, giving a ram pressure of $1.2\times 10^{-8}\,$g cm${}^{-1}\,$s*-2*. This is several orders of magnitude lower than appears in any of the three disk masses simulated (see Figure 3) so no direct stripping is possible. While there may be some long term, continuous stripping, we expect stripping from photoevaporation to dominate such periods.

In conclusion, photoevaporation and hydrodynamic ablation are processes which can both affect the disk prior to it being subject to a SN explosion. However, we expect some disks to survive these earlier processes. For those that do, supernova induced stripping can far exceed the photoevaporative mass-loss rate and clearly dominates the mass loss rate during the interaction of the disk with the SNR. This is the scenario which we study in this work.

A final issue concerns the fact that the radiation and winds of nearby massive stars will alter the ambient medium into which the supernova explodes. For instance, a steady wind will produce a density profile which declines as $r^{-2}$. Alternatively, the massive star may undergo an outburst and eject a dense shell of gas (perhaps the most famous example of this is the Homunculus Nebula produced by $\eta$ Carinae). However, unless the mass that the supernova remnant encounters is comparable to the ejecta mass, the SNR will remain in its free-expansion phase and thus its dynamics will be largely unaffected by the surrounding material. In our simulations, the swept up material is a tiny fraction of the ejecta mass at $0.3\,$pc ($\sim 0.015\%$ for the SNR values used here).

6.4 Planet formation

As planet formation occurs on timescales much longer than the SNR interaction simulated in this work, the disk is likely to re-establish equilibrium before planets begin to form. As far as planet formation is concerned, the interaction is therefore equivalent to replacing the disk with one of smaller mass.

6.5 Enrichment via supernova

Overall very little mass from the SNR flow becomes bound to the disk: $\sim 10^{-9}\,$M*☉* for face-on disks that survive the stripping (models const00med and const00high) and $\sim 10^{-10}\,$M*☉* for the edge-on disk that survives (model dyn90low). This is about 10 times less than seen by Ouellette et al. (2007). Some of the difference may be due to the differences in the dimensionality of the simulations (3D in our case, 2D in theirs). However, most of it is likely due to differences in the radiative cooling (which is included in their simulations, while ours are adiabatic). Thus the hot supernova ejecta is more able to cool down and mix with the disk in their simulations. The broad conclusion is the same, however: that the enrichment of protoplanetary disks via supernova ejecta by pure hydrodynamic mixing is too inefficient to explain the abundance of SLRs in the early solar system. The highest value of $f_{enr}$ seen in our simulations with surviving disks is $\approx 5\times 10^{-6}$ (in model dyn90low - see Fig. 13), which falls short of the $10^{-4}$ needed using SLR production values from Woosley & Weaver (1995). Interestingly, this is for our low mass, edge-on case. Low mass disks obviously benefit more from the same amount of ejecta material, but are also generally destroyed more easily. However, we believe that the fact that edge-on disks survive significantly longer than face-on disks allows such disks to intercept more of the SN ejecta, and thus have the highest enrichment fraction. So, at least for the particular parameters chosen here, a low mass edge-on disk can be enriched more than a high mass face-on disk. While the disks here do not meet the enrichment requirements deduced from observations, it may be possible that a much larger disk (and hence larger surface area for collection) placed further away from the supernova (such that its outer parts do not get instantaneously stripped) could be more efficient at absorbing SN ejecta. This is left for the subject of future investigation.

We have also only considered gas phase ejecta. It has been suggested (by e.g. Goodson et al., 2016) that dust grains would be able to penetrate and be absorbed more easily into the disk. However, the abundance of dust grains in the SNR is uncertain (see e.g. Bianchi & Schneider, 2007; De Looze et al., 2017). The dust must also survive the reverse shock as the SNR interacts with the surrounding medium (Micelotta et al., 2016) although this may not be significant for the close proximity SNRs required for disk enrichment.

Finally we note that the enrichment which we see in our simulations is limited to the surface of the disk. The timescales in our simulations are too short for mixing of the SLRs throughout the disk (Boss, 2013, note a mixing time of $10^{4}$ yr), and in any case the resolution is such that we do not adequately resolve the motions responsible for mixing material within the disk.

7 Conclusion

We have presented three dimensional simulations of the stripping of stellar disks due to the influence of a nearby SNR, using a physically motivated dynamic flow. We have also investigated the effect of varying the inclination angle and disk mass.While a number of other processes (e.g. photoevaporation and ablation by the radiation field and wind from a nearby massive star) can cause significant stripping, or even complete destruction of the disk, prior to interaction with a SNR, this is not always the case. We therefore assume in our investigation that at least some part of the disk is present when a nearby massive star explodes.

Good agreement is found with the analytical predications of Chevalier (2000). However, this only accounts for part of the stripping as Kelvin-Helmholtz instabilities can cause additional material to be ablated. In the initial, instantaneous stripping phase, a flow at the peak ram pressure can strip $90\%$ of a low mass disk ($M_{d}\sim 0.1M_{J}$) and $30\%$ of a medium mass disk ($M_{d}\sim 1.0M_{J}$) on timescales of $10$-$100\,$yrs (less than one outer rotation period). High mass disks ($M_{d}\sim 10M_{J}$) are largely unaffected by instantaneous stripping.

During continuous, longer-term, ablation disks lose mass at a rate of $\sim 10^{-6}\,$M${}_{\sun}\,$yr*-1*. This value decreases with time as the SNR passes and the flow weakens, but is several orders of magnitude greater than the mass-loss rate due to photoevaporation or stellar wind ablation and will therefore dominate the disk’s mass-loss rate during this time.

We find that the inclination angle only has a large effect on the evolution when the disk is close to edge on (similarly to previous findings from simulations of disk galaxies). When the ram pressure is large compared to the gravitational pressure in the disk, low inclination angle disks are deformed to the point that their evolution converges to that of a face-on (zero degree inclination) disk. In contrast, our edge-on disks show a much steadier rate of mass-loss (instantaneous stripping is much reduced due to the lower cross-section) and can survive significantly longer than their face-on counterparts. Amongst the low-mass disks simulated, only the edge-on disk survived interaction with the SNR (retaining almost 60% of its mass).

The stripping of inclined disks can be quite asymmetrical, and the direction of the stripped tail may not line up with the direction of the flow (also like disk galaxies). However, unlike disk galaxies, the flow may die down before the asymmetries have disappeared. For the SNR parameters chosen the interaction is very short, lasting only a couple of hundred years. This means that stellar disks are unlikely to be observed during this period, but seeing an asymmetric disk may be evidence that it underwent this type of ablation in its past. For disks more distant from the supernova the duration of significant interaction will be longer.

Amongst all of our simulations, the highest ejecta enrichment fraction is $5\times 10^{-6}$, which is too low to explain the presence of SLRs in the early solar system. The highest enrichment is seen in a low mass edge-on disk, suggesting the ideal case for enrichment is a low mass edge-on disk (that would be destroyed if placed face-on) rather than a face-on disk which one might naively assume.

Acknowledgements

The calculations for this paper were performed on the DiRAC 1 Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Leeds.