ALMA Resolves C i Emission from the β Pictoris Debris Disk

The CS 10-9 and SiO 11-10 lines remained undetected. We detect and resolve C i ³P₁–³P₀ emission. The left panel of Figure 1 shows the moment 0 map, produced by integrating the data cube along the spectral axis within ±6 km s⁻¹ (with respect to the systemic velocity of the star, assumed to be v_heliocentric = 20.5 km s⁻¹; Brandeker 2011). The figure has been rotated to align the horizontal direction with the main dust disk (position angle of +293; Lagrange et al. 2012). The beam size is 118 × 095 (23 au × 19 au) with a major-axis position angle with respect to the north of 44°. We achieve a 1σ sensitivity of ∼2.3 × 10⁻¹⁷ W m⁻² Hz⁻¹ sr⁻¹ (∼70 mJy beam⁻¹) in the individual channels. We perform photometry by considering a rectangular aperture extending ±115 au in the horizontal direction (measured from the stellar position; see Figure 1) and ±30 au in the vertical direction (measured from the midplane). This yields a total flux of (1.6 ± 0.2) × 10⁻¹⁹ W m⁻² (9.8 ± 1.4 Jy km s⁻¹), where the error is random without any systematic calibration uncertainty taken into account. We did not correct for the primary beam, as it changes the total flux only within the quoted error bars (the same applies for the continuum; Section 3.3). The measured flux is consistent with the value of (1.7 ± 0.4) × 10⁻¹⁹ W m⁻² (10.3 ± 2.3 Jy km s⁻¹) derived from single-dish observations by Higuchi et al. (2017). The same method yields 3σ upper limits of 5.9 × 10⁻²⁰ W m⁻² (3.6 Jy km s⁻¹) for CS 10–9 and 2.9 × 10⁻²⁰ W m⁻² (1.8 Jy km s⁻¹) for SiO 11-10.

Zoom In Zoom Out Reset image size

The right panel of Figure 1 shows the emission profile along the x-axis of the moment 0 map, obtained by integrating within ±30 au in the vertical direction. Both the moment 0 map and the emission profile are suggestive of an asymmetry with a peak on the SW side of the disk. The same asymmetry is seen for the CO emission (Dent et al. 2014; Matrà et al. 2017a) and tentatively also for C ii (Cataldi et al. 2014). However, the SW/NE flux ratio within the rectangular box of Figure 1 is not significant at 0.9 ± 0.2 (we calculate the error on the ratio $\tfrac{a}{b}$ by propagating the error like this: ${\sigma }_{a/b}^{2}=1/{b}^{2}{\sigma }_{a}^{2}+{a}^{2}/{b}^{4}{\sigma }_{b}^{2}$). The SW/NE ratio of the peaks in the emission profile is 1.3 ± 0.5. Thus, there is not necessarily more flux on the SW side of the disk, but the flux seems to be more compact.

Figure 2 shows the position–velocity (pv) diagram (i.e., the data cube integrated along the vertical spatial direction within ±30 au). This figure thus shows the radial velocity of the emission as a function of the projected position along the disk midplane. By using the pv diagram, we can constrain the radial distribution of the emission despite the edge-on orientation of the disk. Indeed, Figure 2 also shows the radial velocity for circular Keplerian orbits at 50 and 220 au around a 1.75 M${}_{\odot }$ star (Crifo et al. 1997), seen edge-on. This is the approximate radial extent of the CO as found by Matrà et al. (2017a). As can be seen, no significant C i emission is detected beyond these lines, suggesting that most C i emission is confined to the same region as the CO.

Zoom In Zoom Out Reset image size

We may use the pv diagram to define a region in p–v space that contains all significant emission. Then, when integrating over the spectral axis, we only take data points within this region into account (rather than everything within ±6 km s⁻¹ as was done to produce Figure 1). The region is illustrated by the dashed lines in Figure 2, and the resulting moment 0 map and emission profile are shown in Figure 3. The SW/NE flux ratio (within the same box as in Figure 1) is now 1.2 ± 0.2. Most importantly, the S/N of the emission profile is significantly improved and the SW/NE asymmetry more clearly visible. We measure a SW/NE peak ratio of 1.6 ± 0.4. Thus, the significance of the peak ratio asymmetry is only marginal. However, in the right panel of Figure 3 we also show the profiles of the CO 3–2 and 2–1 emission (see Figure 2 of Matrà et al. 2017a), for which the SW/NE flux ratios are 1.08 ± 0.08 (CO 2–1) and 1.49 ± 0.13 (CO 3–2) and the peak ratios are 1.42 ± 0.14 (CO 2–1) and 1.51 ± 0.14 (CO 3–2). The C i emission follows the CO 2–1 emission surprisingly well, with the same distinct peak on the SW side of the disk and the same asymmetric extent (out to ∼150 au in the NE and ∼100 au in the SW). This suggests that the asymmetry is indeed real. This is a surprising result. The asymmetry in CO can be readily understood from its short (less than one orbit) lifetime due to photodissociation: if CO is primarily produced in a clump, photodissociation prevents CO from spreading azimuthally, thus preserving the asymmetry. On the other hand, the C produced from CO photodissociation is expected to spread in azimuth within a few orbits. In Section 5.5 we discuss possible solutions to this puzzle.

Zoom In Zoom Out Reset image size

From the moment 0 maps, it is also apparent that the observed emission is slightly tilted, with the NE side below and the SW side above the midplane of the main outer disk (defined by z = 0). A similar tilt is seen for CO (Dent et al. 2014; Matrà et al. 2017a). As is discussed by Matrà et al. (2017a), two reasons might be imagined for this observation: either the gas disk is indeed tilted with respect to the dust disk, or the tilt is a projection effect, resulting from an azimuthally asymmetric gas disk in combination with a slight deviation from a perfect edge-on inclination. Interestingly, the inner dust disk seen in scattered light has a similar tilt (e.g., Milli et al. 2014; Apai et al. 2015; Millar-Blanchaer et al. 2015). This dust component, known as warp or secondary disk, is localized inward of 80 au (Lagrange et al. 2012) and seems thus slightly inward of the gas. Also, it has already been suggested that this inner dust disk is not perfectly edge-on (e.g., Milli et al. 2014; Millar-Blanchaer et al. 2015).

The reader interested in the procedures to estimate the errors quoted in this section is referred to Appendix A.

Assuming that the gas follows circular Keplerian orbits, we can use the pv diagram (Figure 2) to obtain a face-on view of the emission (e.g., Dent et al. 2014; Matrà et al. 2017a). Indeed, each point in the pv diagram corresponds to two points in the xy-plane of the disk (where we define y as the coordinate running along the line of sight), one in front of and one behind the sky plane. There remains the degeneracy of how to distribute the flux of a given pv point among the two points in the xy-plane. As discussed by Matrà et al. (2017a), the degeneracy can be broken if the disk is not perfectly edge-on. However, for simplicity, we follow Dent et al. (2014) and assume that the disk is edge-on. Our primary interest is to illustrate the radial extent of the emission and the position of the clump. We thus choose to place all flux in front of the sky plane, but other physically motivated choices are possible (Dent et al. 2014). Figure 4 shows the deprojection. Points of the pv diagram with $| x| \lt 40\,\mathrm{au}$ were masked because the radial velocity in this region tends toward zero for all radii, i.e., it becomes difficult to assign a radius to a given radial velocity. Emission is seen approximately out to 150 au. The clump appears at a similar position angle as seen in CO (Dent et al. 2014). Details of the deprojection procedure are described in Appendix C. Since the optical depth of the emission is not negligible, the deprojection does not show the true distribution of the emission in the xy-plane, but rather how much emission the observer receives from various locations in the xy-plane.

Zoom In Zoom Out Reset image size

The continuum image at ∼485 GHz is shown in Figure 5. The 1σ noise level is 4.4 × 10⁻¹⁹ W m⁻² Hz⁻¹ sr⁻¹ (1.3 mJy beam⁻¹). The beam size is 119 × 096 (23 au × 19 au) with the position angle of the major axis being 46°. We measure a total flux of (1.12 ± 0.07) × 10⁻²⁷ W m⁻² Hz⁻¹ (112 ± 7 mJy) in the rectangular region (±140 au from the star along x and ±30 au from the midplane along z) indicated in the figure (see Appendix A for details on the error calculation). The measured flux is consistent with a Rayleigh–Jeans extrapolation of the flux measured by ALMA at 870 μm (Dent et al. 2014) and is a factor of ∼2 below what is expected from the infrared to millimeter SED fit by Vandenbussche et al. (2010).

Zoom In Zoom Out Reset image size

As was already seen in the ALMA data by Dent et al. (2014), the continuum is brighter on the SW side of the disk at projected separations between 50 and 100 au. However, the SW/NE integrated flux ratio is only 1.17 ± 0.14. In any case, the asymmetry is analogous to what is seen in thermal mid-IR images between 8.7 and 18.3 μm (Telesco et al. 2005), as well as for C and CO.

In this section, we estimate the total carbon mass in the β Pic disk under some simplifying assumptions. With a 1D model, we want to reproduce the C i flux measured in the present work and the C ii flux measured with Herschel/PACS (Pilbratt et al. 2010; Poglitsch et al. 2010) by Brandeker et al. (2016; ObsID 1342198171) (we sum the flux values of the PACS spaxels listed in their Table 1). Herschel/HIFI also measured the C ii flux (Cataldi et al. 2014). However, the flux measured by PACS is likely a better estimate of the total C ii flux from the disk, because the HIFI half-power beam width is only ∼200 au.

We need to compute the statistical equilibrium (SE) of the level populations, as in general the emission cannot be assumed to be in local thermal equilibrium (LTE). The gas in the β Pic disk is thought to be of secondary origin and thus depleted in hydrogen (e.g., Matrà et al. 2017a). We thus assume that the dominant collider species is electrons and that photoionization of carbon is the main electron source (Fernández et al. 2006; Kral et al. 2016), i.e., the density of ionized carbon equals the electron density. Under these assumptions, the total carbon mass is estimated in two steps. For a given kinetic temperature, we first determine the amount of ionized carbon necessary to reproduce the C ii emission. The second step is to determine the amount of neutral carbon necessary to reproduce the C i flux observed by ALMA, using the electron density derived in the first step.

To solve the SE and radiative transfer, we use pythonradex,¹⁸ a python implementation of the RADEX code (van der Tak et al. 2007) with atomic data from the LAMDA database¹⁹ (Schöier et al. 2005). This code uses an escape probability formalism to solve the radiative transfer and assumes the geometry of a uniform sphere. In general, the radiative transfer depends on the geometry of the emitting region. The gas observed around β Pic is clearly nonspherical and nonsymmetric. The assumption of a uniform sphere is thus a clear limitation of this model. The masses derived here should therefore be considered first-order estimates.

We choose the size of the sphere such that its volume corresponds to the volume of an elliptic torus with semimajor axis of 35 au in the radial direction and semiminor axis of 10 au in the vertical direction (see Figures 1–3). Note that even in the optically thin case, such an assumption on the scale over which the emission is produced is necessary, unless the emission is in LTE. This is because for a given mass the electron density depends on the assumed volume.

We include radiative excitation and de-excitation (hereafter simply (de-)excitation) by the cosmic microwave background (CMB), the stellar radiation (at 85 au), and the dust continuum, where the latter is taken at 85 au as seen in Figure B1 of Kral et al. (2017) and is the dominant component. However, it turns out that including radiative (de-)excitation from these sources does not change our results because (de-)excitation is dominated by collisions and spontaneous emission.

We then consider a wide range of kinetic temperatures T_kin from 40 to 1000 K. For lower temperatures, the C ii line quickly becomes strongly optically thick such that meaningful constraints on the mass are no longer possible (i.e., the model cannot reproduce the observed flux regardless of how much the mass is increased). However, based on our more detailed modeling in Section 4.2, we deem lower temperatures unlikely. Detailed thermodynamical modeling by Kral et al. (2016) also suggests that the temperature is above ∼50 K within 100 au, although in their model the temperature drops to 20 K at 150 au.

Figure 6 shows the derived total carbon mass as a function of T_kin. The figure also shows the individual C⁰ and C⁺ masses. To assess the importance of non-LTE effects, we also show corresponding curves for which LTE has been assumed. Both C i and C ii are, generally speaking, close to LTE. For higher temperatures, the C⁰ mass required to reproduce the observed C i flux is higher in LTE than in non-LTE. This simply happens because in LTE higher levels get populated more quickly with increasing temperature, thus depopulating the level that produces the C i emission. The maximum optical depths encountered for the considered temperature range are 0.5 for C i and 3.9 for C ii.

Zoom In Zoom Out Reset image size

From Figure 6, we determine a total carbon mass between 5 × 10⁻⁴ and 3.1 × 10⁻³ ${M}_{\oplus }$. The lower bound is quite robust to changes of the size of the emitting region, as it is close to the LTE value. For example, increasing or decreasing the radius of the emitting sphere by 50% does not change the lower bound by more than 15%. On the other hand, it is clear that the upper bound is more uncertain, as it can quickly increase if one allows for lower temperatures and/or smaller emitting volumes (i.e., increased optical depth). Another parameter that can strongly affect the optical depth (and thus the upper bound of the mass range) is the assumed line broadening parameter. Here we used b = 0.57 km s⁻¹, derived in Appendix B.

Previous studies estimating the carbon gas mass in the β Pic disk had only the C ii flux and line profile at disposal. These more detailed models derived higher masses by using the spectrally resolved C ii observations from Herschel/HIFI: Cataldi et al. (2014) obtained 1.3 × 10⁻² ${M}_{\oplus }$, while Kral et al. (2016) derived 1.5 × 10⁻² ${M}_{\oplus }$ (the latter study also used an upper limit on the C i flux). In the Kral et al. (2016) model, the total C mass is dominated by ionized carbon that is located beyond ∼100 au. The C i flux is of little importance for the total C mass budget of that model. The discrepancy thus cannot be explained by the fact that we include the C i measurement. On the other hand, most of the carbon in the Cataldi et al. (2014) model is located between 150 and 300 au with an ionization fraction of roughly 50%, i.e., there is a significant contribution of neutral carbon to the total mass. However, our ALMA data show no C i emission beyond ∼120 au. That model indeed overpredicted the C i flux by a factor of ∼20, suggesting that it contains too much neutral carbon, partially explaining the difference in the mass estimates. But Cataldi et al. (2014) also derived a higher mass of ionized carbon compared to our estimate. In addition, their fit to the resolved C ii line profile (see their Figure 2(b)) clearly suggests that C ii emission beyond 150 au is needed to fit the line core (this is also an issue for our 3D models [Section 5.2], which generally do a bad job in fitting the C ii line core). So maybe there is largely ionized carbon gas beyond ∼150 au present and the reason why we derive a lower ionized carbon mass here is because we assumed (based on the C i data) a too small volume (i.e., too high density and thus more excitation, and therefore less mass needed). However, even when increasing the volume of our 1D model, we still derive a lower mass of ionized carbon compared to Kral et al. (2016) and Cataldi et al. (2014) and thus cannot fully resolve the discrepancy, which might be due to the different assumptions of the models. The relatively small HIFI beam (FWHM of ∼200 au) could also play a role. For example, deviations from axisymmetry in the distribution of ionized carbon could introduce additional uncertainty in the mass estimates of Cataldi et al. (2014) and Kral et al. (2016) since C ii emission from large radii is only detected by HIFI if it arises close to the line of sight toward the star.

In this section, we model the 3D spatial distribution of the carbon gas using our new ALMA observations of C i and the previously published spectrally resolved Herschel/HIFI observation of C ii (Cataldi et al. 2014; ObsID 1342238190). For a given, arbitrary distribution of carbon gas, we first solve the ionization balance in each grid cell. We use the ionization rate from Zagorovsky et al. (2010) (scaled with distance to the star) and assume that the gas is optically thin to ionizing photons. The star is the main source of ionizing photons, but ionization from the interstellar radiation field (ISRF) is also included. The ionization balance is solved analytically by assuming that all electrons come from the photoionization of carbon. The ionization fraction thus only depends on the distance to the star, the local gas density, and the kinetic temperature (via the recombination coefficient). Next, we use the derived electron density to solve the SE of the level populations using atomic data from the LAMDA database. The following processes can (de)excite an atom: spontaneous emission, collisions with electrons (we neglect other colliders), and radiative (de)excitation by line photons and the background radiation field, where the latter is assumed to be composed of the CMB, the star, and the dust continuum. For the dust continuum, we employ the field shown in Figure B1 of Kral et al. (2017) (that figure actually shows the field in the midplane of the disk at the position angle of the clump [L. Matrà 2018, private communication], but for simplicity, we only take the radial dependence into account). In principle, a full radiative transfer computation is necessary to solve the SE. However, we take a simplified approach and assume that the line emission is essentially optically thin. Basically, while the emission can become optically thick along the disk midplane (in our best-fit models, the maximum optical depth is typically of the order of 1), it is optically thin in the vertical direction, i.e., the photon escape fraction is high. Thus, we assume that the background field is not attenuated by the gas (all atoms are subject to the full background field) and that (de)excitation by line photons can be neglected. These assumptions make the SE easy to solve. We further discuss this approximation in Appendix D. The background radiation turns out to be unimportant for our models.

Having solved the SE, we compute the emitted spectrum for each grid cell, redshifting or blueshifting it according to its radial velocity. The final step is then to ray-trace the emission along the line of sight to take optical depth into account. The result is a model data cube that can be compared to observations. For simplicity, we consider isothermal models and fix T_kin = 75 K everywhere. We found that for lower temperatures the C i/C ii flux ratio tends to be too high, while the inverse is true for higher temperatures.

To compare to the C ii data, we take the corresponding model data cube, multiply by the HIFI beam, and integrate spatially to get a model HIFI spectrum (Cataldi et al. 2014). For the C i ALMA observations, we convolve the model data cube to the same spatial and spectral resolution as the observations and multiply by the primary beam. The residual moment 0 maps (respectively pv diagrams) shown in Figures 7, 8, and 11 are obtained by subtracting the model moment 0 map (pv diagram) from the data moment 0 map (pv diagram) shown in Figure 1 (Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. Uniform Ring Model

We first consider a simple, symmetric model consisting of a ring with uniform surface density. The number density reads

$$\begin{eqnarray}&&{n}_{\mathrm{ring}}(r,z)\\ =\,\left\{\begin{array}{ll}\tfrac{{\rm{\Sigma }}}{\sqrt{2\pi }H(r)}\cdot \exp \left(-\tfrac{{z}^{2}}{2{(H(r))}^{2}}\right) & \mathrm{if}\ {r}_{\min }\leqslant r\leqslant {r}_{\max }\\ 0 & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 1 }$$

where r and z are cylindrical coordinates (with z perpendicular to the disk midplane), H(r) is the scale height, r_min and r_max define the radial extent of the ring, and Σ is the constant surface density. The scale height in hydrostatic equilibrium for a vertically isothermal disk is given by (e.g., Armitage 2009)

$$\begin{eqnarray}&&H(r)=\sqrt{\displaystyle \frac{{{kTr}}^{3}}{\mu {m}_{p}{{GM}}_{* }}},\end{eqnarray} \tag{ 2 }$$

with k the Boltzmann constant, T the temperature, m_p the proton mass, G the gravitational constant, and M_* the stellar mass. For the mean molecular weight μ, we follow Matrà et al. (2017a) and assume μ = 14. At 85 au and for T = 75 K, the scale height is 4.2 au.

We compute a grid of models over the parameter space listed in the first four rows of Table 1. Note that we also test different values for the (dynamical) stellar mass that affects the orbital velocities (and scale height). To each model, we assign a χ² measure by using expressions analogous to Equation (2) in Booth et al. (2017) (we take the correlation of neighboring data points into account by using an appropriate noise correlation ratio for each data set). We sum the χ² from the Herschel/HIFI data (C ii) and the ALMA data (C i).

Table 1.  Explored Parameter Space and Best-fit Values for the "Ring" and "Ring + Clump" Models Fitting the C i ALMA Data and the C ii Herschel/HIFI Data Simultaneously

Parameter	Unit	min	max	n	Spacing	Best Fit
						"Ring"	"Ring + Clump"
M*	MβPic	0.6	1	3	lin	0.8	0.8
Mring	${M}_{\oplus }$	2.7 × 10−4	2.7 × 10−3	6	log	6.7 × 10−4	6.7 × 10−4
rmin	au	30	70	3	lin	50	50
rmax	au	120	160	3	lin	120	120
Mclump	${M}_{\oplus }$	2.7 × 10−5	2.7 × 10−4	6	log	⋯	1.1 × 10−4
rc	au	70	100	3	lin	⋯	70
σxy	au	20	40	3	lin	⋯	40
Σ	cm−2	⋯	⋯	⋯	⋯	2.4 × 1016	2.4 × 1016
nmid (r = 85 au)	cm−3	⋯	⋯	⋯	⋯	140	140
nclump	cm−3	⋯	⋯	⋯	⋯	⋯	280

Note. The number of values explored for each parameter is denoted by n. We also indicate whether the values are uniformly distributed in linear or logarithmic space. M_ring and M_clump are the masses of the ring and the clump, respectively. The reference mass of β Pic is assumed to be M_{β Pic} = 1.75M_* (Crifo et al. 1997). We also list the constant surface density of the ring Σ, the midplane number density of the ring n_mid at 85 au, and, for the "ring + clump" model, the number density at the center of the clump n_clump resulting from the superposition of the ring and the clump.

Download table as:  ASCIITypeset image

Figure 7 shows the C i residuals in the xz and xv plane, as well as the C ii line emission of the model with the lowest χ². The corresponding model parameters are given in Table 1. The model provides a decent fit to the data, but unsurprisingly it is not able to reproduce the clump observed in the SW. Thus, we consider a more complicated model by adding a clump to the uniform ring in the next section. Such a model has no direct physical justification but is useful to empirically constrain the spatial distribution of the gas and get an estimate of the gas mass.

4.2.2. Uniform Ring with a Clump

The clump is modeled as

$$\begin{eqnarray}{n}_{\mathrm{clump}}(x,y,z) & = & {n}_{0}\cdot \exp \left(-\displaystyle \frac{{(x-{x}_{0})}^{2}+{(y-{y}_{0})}^{2}}{2{\sigma }_{{xy}}^{2}}\right)\\ & & \cdot \exp \left(-\displaystyle \frac{{z}^{2}}{2{(H(r))}^{2}}\right),\end{eqnarray} \tag{ 3 }$$

where x runs along the disk midplane in the sky plane and y along the line of sight (and x² + y² = r²). Furthermore, ${x}_{0}={r}_{{\rm{c}}}\cos ({\phi }_{{\rm{c}}})$ and ${y}_{0}={r}_{{\rm{c}}}\sin ({\phi }_{{\rm{c}}})$ designate the center of the clump, where we have defined r_c as the clump's radial distance to the star and ϕ_c as its azimuthal angle in the x–y-plane. We fix ϕ_c = −32° (Matrà et al. 2017a). The standard deviation of the clump density distribution in the x–y-plane is σ_xy, and the density at the center is n₀. The total carbon number density is then given by ${n}_{{\rm{C}}}={n}_{\mathrm{ring}}+{n}_{\mathrm{clump}}$. We compute models over the parameter space listed in Table 1. Figure 8 is the analog of Figure 7 for the best "ring + clump" model. As can be seen, the model does a slightly better job in modeling the clump.

4.2.3. Eccentric Gas Distribution

The "ring + clump" model of the previous section is purely empirical and does not have a direct physical justification. As mentioned earlier, proposed mechanisms to get such a morphology are either a giant collision or resonant trapping of cometary bodies by a migrating giant planet. However, as we discuss in Section 5.5, based on our new C i data, we deem both of these possibilities unlikely.

Another way to get a morphology with an emission clump on one side of the disk is to relax the assumption of circular orbits and instead consider gas on eccentric orbits. In Section 5.5.6, we discuss how such orbits could arise in the first place. As an example, we here consider a model where the eccentricity of the orbits is uniformly distributed between two values e_min and e_max and where all orbits share the same pericenter and the same argument of periapsis. The pericenter is then a region of higher density and can mimic a clump as seen in the observations.

The derivation of the gas density for this model is given in Appendix E. We again compute a grid of models with the parameters listed in Table 2. Figure 9 shows a face-on view of the midplane carbon gas density, and Figure 10 shows the pv diagram of the modeled C i emission. As is seen in Figure 11, this model is similarly effective in reproducing the observed C i emission, although it has some difficulties in reproducing the C ii line. We emphasize that the purpose of this model is merely to demonstrate that eccentric gas orbits are able to reproduce the general morphology with the clump. A more detailed model will be presented in a forthcoming paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Same as Table 1, but for the Parameter Space Explored by the Models with Eccentric Orbits

Parameter	Unit	min	max	n	Spacing	Best Fit
M*	MβPic	0.6	1	3	lin	0.8
Mtot	${M}_{\oplus }$	2.0 × 10−4	3.3 × 10−3	6	log	6.2 × 10−4
emin	⋯	0	0.4	5	lin	0
emax	⋯	0.1	0.5	5	lin	0.3
rper	au	60	100	3	lin	80
ω	deg	−135	135	7	lin	−45

Note. The total mass of the model is M_tot, the pericenter distance is r_per, and the argument of pericenter is ω.

Download table as:  ASCIITypeset image

An interesting consequence of eccentric orbits is that the gas along the line of sight toward the star has nonzero radial velocity. Thus, emission (or absorption) toward the star is shifted with respect to the systemic velocity. For example, our best-fit model predicts that the emission peak toward the star appears 0.4 km s⁻¹ blueshifted. Another consequence is an additional velocity broadening compared to the circular case, with broadening parameter b_ecc ≈ 0.4 km s⁻¹ (this is smaller than and thus consistent with the broadening measured from the pv diagram; see Appendix B). The velocity shift and broadening depend on the eccentricity and orientation of the disk. To verify this velocity shift and constrain the model, we would need to know the absolute stellar radial velocity to better accuracy than ∼80 m s⁻¹ (for a 5σ detection of the shift).

The column densities of neutral carbon toward the star for the three models discussed in Sections 4.2.1–4.2.3 are (6–7) × 10¹⁶ cm⁻², slightly above the value of (2–4) × 10¹⁶ cm⁻² measured in absorption by Roberge et al. (2000). For ionized carbon, the models range between 5 × 10¹⁶ cm⁻² and 1.2 × 10¹⁷ cm⁻², while Roberge et al. (2006) report 2 × 10¹⁶ cm⁻².

Recently, Kral et al. (2016) presented a model that computes the temperature, ionization, and hydrodynamical evolution of the atomic carbon and oxygen in the β Pic disk. In this model, the atomic gas is produced in a parent belt from photodissociation of CO and then evolves viscously under the influence of the MRI, i.e., it forms an accretion and decretion disk. Kral et al. (2016) predict that the carbon is mostly neutral in the inner parts of the accretion disk.

However, Figures 2 and 4 indicate that no atomic accretion disk has formed (yet), as there seems to be little C i emission inside of ∼50 au (but see also Section 5.4). To confirm this, we compute the C i emission expected from the model and compare it to our data. We thus take the Kral et al. (2016) distribution of neutral carbon and electrons, temperature profile, and scale height (see their Figure 9) and compute the C i emission with the methods described in Section 4.2. Figure 12 shows the moment 0 map and pv diagram of the model and the residuals when subtracting the model from the ALMA C i data. It is clear from this figure that the accretion profile predicted by Kral et al. (2016) produces too much C i emission in the inner regions of the disk. This is also clearly visible in the pv diagram, where too much emission is present at high velocities.

Zoom In Zoom Out Reset image size

For all the 3D models described in Section 4.2, the grid searches indicate that a dynamical mass of 0.8M_* (with ${M}_{* }=1.75\,{M}_{\odot }$ the assumed stellar mass; Crifo et al. 1997) is preferred to reproduce the data. While this result has to be interpreted with care given that we did not derive any error bars, it is interesting to note that Olofsson et al. (2001) derived the same 0.8M_* dynamical mass by modeling spatially and spectrally resolved Na i emission. They ascribed their finding to radiation pressure opposing the gravity of the star. In fact, radiation pressure on Na is so strong that a braking mechanism is required to keep it on the observed Keplerian orbits and prevent it from being blown out (e.g., Liseau 2003; Brandeker et al. 2004). Fernández et al. (2006) suggested that all species affected by radiation pressure are largely ionized and coupled into a single fluid by Coulomb interactions, with carbon acting as a braking agent. In this situation, neutrals are also expected to be coupled (the only exception would be neutrals that do not ionize yet still feel a significant radiation force). Thus, one might actually expect that Na i and C i share the same dynamics. The dynamical mass (respectively the radiation pressure coefficient) is coupled to the composition of the gas (Fernández et al. 2006) and could thus give interesting information on the elemental abundances. However, given the unknown uncertainty of the derived dynamical mass, which could also be influenced by the assumptions of our models, we do not draw any further conclusions at this stage.

All the 3D models presented in Section 4.2 have difficulties in reproducing the core of the C ii line profile (see Figures 7, 8 and 11). The issue is particularly pronounced for the model with eccentric orbits. The fits to the C ii line profile by Cataldi et al. (2014) indicate that the line core requires ionized carbon beyond ∼150 au to be present. As already discussed in Section 4.1, there might thus exist a gas component with a high ionization fraction beyond ∼150 au. This component might also be required to prevent the metals observed there from being blown out by radiation pressure (see Section 5.7). It is possible that our simple models are unable to capture this extended gas component. For example, a more complex density and/or temperature profile might be required to reproduce it. Note that while the Cataldi et al. (2014) model fits the C ii line profile, it would be a bad fit to our new C i data, as it strongly overpredicts the C i flux.

The strong spatial correlation between the C i and CO emission (Figure 3) suggests a scenario where the carbon is mainly produced from photodissociation of CO, i.e., the mass-loss rate of CO equals the production rate of C. Thus, from the CO mass-loss rate and our estimate of the total C mass, we can estimate the time since C (and CO) production started in the β Pic disk, provided that no carbon has yet been removed.

5.3.1. Revised CO Lifetime

5.3.2. C Production Timescale

Kral et al. (2016) predicted the formation of a C accretion disk based on thermo-hydrodynamical modeling. However, we showed in Section 4.3 that no such accretion disk is present. If the C gas production indeed started as recently as estimated in Section 5.3, the absence of an accretion disk is actually not surprising. Indeed, the timescale for viscous accretion from a radius r can be estimated by

$$\begin{eqnarray}&&{t}_{\mathrm{vis}}=\displaystyle \frac{{r}^{2}}{\alpha {c}_{{\rm{s}}}H},\end{eqnarray} \tag{ 4 }$$

with c_s the sound speed and H the scale height. We assume ${c}_{{\rm{s}}}=\sqrt{\tfrac{{kT}}{\mu {m}_{{\rm{p}}}}}$ with T = 75 K and μ = 14, and H as in Equation (2) with r = 85 au. A very high α ≳ 4 (taking the upper bound of the timescale range calculated in Section 5.3) would be required to match the viscous timescale with the time since C production started. In other words, for any reasonable value of α, the gas has not yet had enough time to form an accretion disk.

However, a certain amount of carbon is still needed in the inner regions of the disk, where metals such as Na or Fe are seen in Keplerian rotation. These species are strongly affected by radiation pressure, and C is needed to prevent them from being blown out. As was shown by Fernández et al. (2006), C, which does not feel any radiation pressure around β Pic, is acting as a braking agent in the form of C⁺ via Coulomb interactions. In order to test whether the amount of carbon needed to brake metals is consistent with the new ALMA C i data, we consider an accretion disk model extending from 10 to 50 au with a carbon surface density Σ_C ∝ r⁻¹ and temperature T ∝ r^−0.5 (with T = 75 K at 85 au; Lynden-Bell & Pringle 1974). We compute the emission from this model as described in Section 4.2. We then search for C i emission only in those regions of the data cube where the model predicts emission. However, we want to exclude regions of the data cube that can contain emission from the outer disk. Thus, for those data points with $| v| \lt {v}_{\mathrm{orb}}(50\,\mathrm{au})$ (with v_orb(r) the orbital speed at radius r), we request points to be at least one spatial resolution element inside of the line in the pv diagram defining a thin ring at 50 au (see Figure 2). For points with a larger $| v|$, we can be sure that no contamination from the outer disk is present. We also exclude points with a predicted emission below 10% of the model peak to avoid considering regions of the data cube where only weak emission is expected anyway. We can then measure the flux in this region of the data cube, with the error estimated with a method analogous to what is described in Appendix A. Defining a detection threshold of 3σ, we do not detect significant emission. We derive an upper limit on the C density by adjusting the model such that the probability for a nondetection, i.e., measuring a flux below the detection threshold, would only be 1% if the model was correct. Assuming Gaussian noise, the probability of a model with flux F_mod to remain undetected is given by $\tfrac{1}{2}\left[1+\mathrm{erf}\left(\tfrac{n\sigma -{F}_{\mathrm{mod}}}{\sqrt{3}\sigma }\right)\right]$, with σ the error on the flux and n = 3 in our case. The upper limit model has a midplane C⁺ number density of ∼600 cm⁻³ at 30 au, while ∼100 cm⁻³ is necessary to brake the metals (Fernández et al. 2006). Thus, the data are consistent with enough carbon being present in the inner disk to brake metals. However, to get better constraints, we also look at the C ii observations from Herschel/HIFI by Cataldi et al. (2014). We measure the flux in the wings of the C ii spectrum corresponding to radial velocities between v_orb(50 au) and v_orb(10 au). Error bars are calculated by calculating the flux in spectral regions of the same size without line emission and taking the standard deviation. We first determine the errors on the flux in the left wing and right wing individually and then add them in quadrature to obtain the error on the total measured flux. No significant (larger than 3σ) flux is detected in either the H or V beam of the data. Thus, we adjust the model such that the combined probability to remain undetected in the ALMA C i and the HIFI C ii (H and V beam) data is only 1% (including the ALMA data has a negligible effect, as the HIFI data are more constraining). This model has a midplane C⁺ density at 30 au of ∼380 cm⁻³. We conclude that the currently available data are consistent with carbon being present inside of 50 au at a level that is sufficient to brake metals. We suggest that this gas in the inner region was not produced in the same event as the gas seen at larger distances. If it was, we would expect a full accretion disk, which is not supported by the data (see also Cataldi et al. 2014). Instead, it might be the leftover from a previous gas-producing event.

A key result from our new ALMA data is that C shows the same asymmetry as CO: a clump on the SW side of the disk. This is surprising since one would expect C to spread in azimuth within a few orbits even though C production might happen preferably at the CO clump. This is in contrast to CO, which has a lifetime shorter than an orbital period and thus remains asymmetric. Here we discuss possible solutions to this puzzle.

5.5.1. Recent Event

Perhaps the simplest explanation for the observed C asymmetry is that C production at the clump location started so recently that there was not yet enough time for azimuthal spreading (respectively symmetrization). We expect the timescale for azimuthal symmetrization to be on the order of a few orbits (at 85 au, the orbital period is ∼600 a). We investigate the symmetrization with a 1D toy model, where the only dimension is the azimuthal angle ϕ. To start, we assume that all C is produced in a single point at azimuth ${\phi }_{0}$, at a distance of 85 au, with a rate equal to the CO destruction rate calculated in Section 5.3. In reality, only 30% of the CO emission is found in the clump (Dent et al. 2014). This setup thus maximizes the asymmetry, so the derived symmetrization timescale can be considered an upper limit. We then write the following simple equations describing the temporal evolution of the neutral and ionized carbon densities under the influence of C production, ionization, recombination, and orbital motion:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{0}}{\partial t}=\delta (\phi -{\phi }_{0}){{\rm{\Lambda }}}_{\mathrm{CO}}\displaystyle \frac{{m}_{{\rm{C}}}}{{m}_{\mathrm{CO}}}+{\rho }_{+}{n}_{{\rm{e}}}\gamma -{\rho }_{0}{\rm{\Gamma }}-\omega \displaystyle \frac{\partial {\rho }_{0}}{\partial \phi }\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{+}}{\partial t}=-{\rho }_{+}{n}_{{\rm{e}}}\gamma +{\rho }_{0}{\rm{\Gamma }}-\omega \displaystyle \frac{\partial {\rho }_{+}}{\partial \phi },\end{eqnarray} \tag{ 6 }$$

where ρ₀ and ρ₊ are the azimuthal densities (in kg rad⁻¹) of neutral and ionized C, respectively, Λ_CO is the CO destruction rate (in kg s⁻¹), m_C and m_CO are atomic and molecular masses, respectively, n_e is the electron density, γ is the recombination coefficient, Γ is the ionization rate, and ω is the angular velocity at 85 au. For the electron density, we assume that all electrons come from the ionization of C and that the typical volume extends over Δr = 70 au in the radial direction (best-fit "ring" model; Section 4.2.1) and over Δz = 2H (with H the scale height; see Equation (2)) in the vertical direction (the model remains one-dimensional; we only assume a certain volume to get a value for the electron density). This means that ${n}_{{\rm{e}}}={\rho }_{+}/({m}_{{{\rm{C}}}^{+}}r{\rm{\Delta }}r{\rm{\Delta }}z)$, where ${m}_{{\rm{C}}}^{+}$ is the mass of a C⁺ ion. Then, Equations (5) and (6) are used to compute the change of the azimuthal densities at each time step. Figure 13 shows how the SW/NE mass ratio of neutral carbon evolves with time. The wavy pattern is due to the orbital motion of the gas. The local minima correspond to the times when the first gas produced by the point source leaves the NE side and enters the SW side. After ∼10³ a, the ratio falls below the value of our best-fit "ring + clump" model.

Zoom In Zoom Out Reset image size

We can also use this model to estimate the time it will take to symmetrize the C from the state that we observe today if there is no mechanism that prevents symmetrization. For example, assuming that only 30% of the C production occurs at the clump, with the other 70% uniformly distributed at all azimuths (this corresponds to adding a production term independent of ϕ to Equation (5)), we get the green curve shown in Figure 13. As expected, the symmetrization occurs faster—within ∼10³ a, the SW/NE mass ratio falls below 1.1.

A caveat to this toy model is the possibility that the gas production is not constant over time. Also, for more robust constraints, detailed hydrodynamical calculations would be necessary.

The upper limit on the symmetrization timescale is below the lower bound of the timescale range over which C production occurred, as derived from the C mass (Section 5.3). Thus, it appears unlikely that the asymmetry can be explained by invoking a very recent event.

In summary, the estimated lifetime (∼5000 a) of the carbon-rich gas disk should be long enough to spread out the azimuthal asymmetry, but not long enough to diffuse the disk radially via viscous spreading.

5.5.2. Resonance Trapping by Planet

5.5.3. Initially Eccentric Disk

If the parent planetesimal disk is eccentric because it has been left in that state by some unknown initial condition, the debris disk will initially be asymmetric. Such a disk, if precessed rigidly, can retain its initial configuration over a time much longer than the sound-crossing time. However, over the lifetime of the system (20 Ma), differential precession should have disturbed this initial imprint. Ignoring the presence of other planets and only considering the quadrupole precessional effect of β Pic b (Lagrange et al. 2009), the timescale for order unity change in the relative precession angle is (Murray & Dermott 2000)

$$\begin{eqnarray}&&{t}_{\mathrm{smear}}\sim \displaystyle \frac{{M}_{* }}{{M}_{{\rm{p}}}}{\left(\displaystyle \frac{a}{{a}_{{\rm{p}}}}\right)}^{2}{\left(\displaystyle \frac{2{\rm{\Delta }}a}{a}\right)}^{-1}{P}_{\mathrm{orb}}\sim 14\,\mathrm{Ma},\end{eqnarray} \tag{ 7 }$$

where we evaluate the expression using M_p = 13 M_Jup (Morzinski et al. 2015) for the planet mass and a_p = 9 au (Wang et al. 2016) for its semimajor axis with the ring at a ∼ 85 au and with a fiducial radial width of Δa ∼ 20 au (consistent with that in Table 2). So even without an additional planet, an initially eccentric disk is expected to be largely smeared out.

5.5.4. Eccentric Disk Secularly Forced by a Planet

5.5.5. Giant Impact

Jackson et al. (2014) proposed a scenario where a giant impact between two comparably massed bodies produces a wide spread of debris that have a range of eccentricities but a similar alignment. This is qualitatively plausible to explain the observations (but see Matrà et al. 2017a, for a counterargument regarding the radial width). If such an event occurred in the recent past, it provides an interesting explanation for our deduced carbon production timescale. However, such an event seems exceedingly improbable, as we will show with an order-of-magnitude estimate of the event rate.

Consider N bodies with radius R, in a belt of semimajor axis a and width Δa, performing vertical epi-cycles around the midplane. Viewed by each one of these bodies, the other bodies present a certain optical depth of ${\tau }_{\perp }=(N4\pi {R}^{2})/(2\pi a{\rm{\Delta }}a)$, or a mean collision time of ${P}_{\mathrm{orb}}/{\tau }_{\perp }$. Summed over N bodies, this implies a mean event time of

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}\sim {P}_{\mathrm{orb}}\displaystyle \frac{2\pi a{\rm{\Delta }}a}{{N}^{2}\times 4\pi {R}^{2}}\sim 3\times {10}^{9}\,\mathrm{years}\\ &&\quad \times \,\left(\displaystyle \frac{a}{80\,\mathrm{au}}\right)\left(\displaystyle \frac{{\rm{\Delta }}a}{20\,\mathrm{au}}\right){\left(\displaystyle \frac{N}{1000}\right)}^{-2}{\left(\displaystyle \frac{R}{2000\mathrm{km}}\right)}^{-2},\end{eqnarray} \tag{ 8 }$$

where we assume that there is no gravitational focusing among these low-mass objects, reasonable if their velocity dispersion lies much beyond their surface escape velocities. The scaled value R = 2000 km is appropriate to explain the total gas/dust mass observed, and Δa = 20 au is inspired by the best fit in Table 2. The value N = 1000 is a placeholder. So, to produce an event as recent as 5000 a ago, one would require N ∼ 10⁵, or an absurdly high total mass of ∼500 M_⊕. This argument makes giant impacts exceedingly implausible at this location in the disk. Another difficulty with such a scenario is that giant impacts tend to be completely accretionary at low relative speeds, and even at speeds beyond the surface escape, only a small fraction of mass can be unbound and released into the circumstellar environment (Agnor & Asphaug 2004).

5.5.6. Tidal Disruption

Here we briefly propose an alternative scenario to produce the observed disk. A more detailed calculation will be presented in an upcoming paper. Consider that the outer disk contains a number of Neptune-like planets (N_N ∼ a few). They have a surface escape velocity of v_esc ∼ 23 km s⁻¹. Let us also assume that there are N bodies similar in mass to the Moon (∼0.01 M_⊕) or Mars (∼0.1 M_⊕) and moving in space with a relatively low dispersion velocity, σ ≪ v_esc. There is little gravitational focusing among these bodies, but strong focusing by the Neptunes. The timescale for a physical collision with a Neptune is

$$\begin{eqnarray} & & {t}_{\mathrm{coll}}\sim {P}_{\mathrm{orb}}\displaystyle \frac{2\pi a{\rm{\Delta }}a}{{{NN}}_{N}\times \pi {R}_{N}^{2}}{\left(\displaystyle \frac{\sigma }{{v}_{\mathrm{esc}}}\right)}^{2}\sim 3\times {10}^{7}\,\,\mathrm{years}\\ & & \times \left(\displaystyle \frac{a}{80\,\mathrm{au}}\right)\left(\displaystyle \frac{{\rm{\Delta }}a}{20\,\mathrm{au}}\right){\left(\displaystyle \frac{N}{1000}\right)}^{-1}\\ & & {\left(\displaystyle \frac{\sigma }{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{v}_{\mathrm{esc}}}{23\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-2},\end{eqnarray} \tag{ 9 }$$

where we chose N_N = 5. This is somewhat arbitrary, but the abundance of cold Neptunes is already suggested by microlensing studies, which claim that their abundance rate per star (mostly M dwarfs) is 52% (Cassan et al. 2012). This timescale is becoming astrophysically interesting. But a physical collision with a Neptune will not produce a debris disk around the star. Instead, we focus on encounters that are close enough that the body is tidally disrupted. This means an encounter distance of ∼2R_N and encounter time ${t}_{\mathrm{disruption}}\sim {t}_{\mathrm{coll}}/2\sim 10$ Ma (at 2R_N, while the geometric cross section goes up by a factor of 4, the gravitational focusing factor, (σ/v_esc)², evaluated at R = 2R_N goes down by a factor of 2). So there could have been a few tidal disruption events over the lifetime of β Pic, if there are thousands of moons floating around and if these moons retain low enough dispersion velocities to experience strong gravitational focusing by Neptune. This is still well above the 5000 a event time we infer and still presents a tension.

The end product of the tidal disruption shares many similarities with that from a giant impact. First, the debris will be ejected on a variety of orbits, bringing about a large spread in eccentricity. The semimajor axis will also have a spread. All debris will return to the disruption site, making this location a region of frequent collisions. The size of the disruption site, however, is larger than that in giant impact. As the debris flies away from Neptune, its stellar-centric orbit is altered by Neptune's gravity while it is still within Neptune's Hill sphere. As a result, unlike the narrow nozzle of the size of the impactor radius in the case of giant impacts, here the nozzle has a typical spread of order the planet's Hill radius, which reduces the peak collision rate. However, Matrà et al. (2017a) measured a radial extent of the CO clump of at least 100 au, which is clearly larger than the planet's Hill sphere. More detailed modeling is needed to see whether a tidal disruption event can produce such a radially broad clump. Also, the clump in the deprojection of Matrà et al. (2017a) might appear more extended than it is in reality because the intrinsic velocity dispersion of the gas and the finite velocity resolution of the instrument degrade the resolution of the deprojection in the y-direction. In addition, the deprojection is carried out assuming velocities corresponding to circular orbits. Thus, if the orbits are eccentric in reality, the deprojection will be distorted and not correspond to the actual gas distribution.

The event in β Pic is recent and must also be relatively short-lived, assuming that we are not observing it at a special time. The duration of a tidal disruption flare will depend not only on the above-discussed collisional geometry but also on the distribution of the largest fragments (that remain or reform) after the disruption event. Both of these effects need to be studied in detail. In addition, a short lifetime will also ensure that debris products from previous events do not interfere with the current event. If not, the asymmetry would be washed out by previous debris that likely have a different asymmetry, and we would observe the radially diffused accretion disk from gas produced in previous events.

In summary, our observation disfavors a few proposed scenarios (planet MMR trapping, giant impact, secular forcing). Instead, we propose that the bright, asymmetric debris disk in β Pic could be the result of a recent tidal disruption of a Moon- to Mars-sized object by a Neptune-like planet. Using the C mass derived in Section 4.1 and assuming that the disrupted body has had a CO mass fraction of 10%, its total mass would indeed be ≳3 M_Moon (lower limit because not all carbon might have been released yet).

Since C and CO have similar horizontal emission profiles (see Figure 3), the question arises how similar the spatial and spectral distributions of the emission really are. To answer this question, we first interpolate the CO data cube onto the coordinates of the C i data cube. Then, we use convolution to adjust the spatio-spectral resolution of the data cubes. Finally, we normalize the data cubes and subtract. Figure 14 shows the moment 0 and pv diagram of the residual cube for the CO 2–1 and 3–2 transitions (Matrà et al. 2017a). The CO emission seems similar to the C i emission, with few residuals above 3σ seen. With the data at hand, we are thus not able to detect clear differences in the emission distribution. Note that there is significant difference in the distribution of CO 2–1 versus CO 3–2 emission (Matrà et al. 2017a).

Zoom In Zoom Out Reset image size

Brandeker et al. (2004) and Nilsson et al. (2012) found that Na and Fe have an NE/SW asymmetry reminiscent of what we see for C and CO (compare Figure 3 of this work to Figure 3 in Brandeker et al. 2004). Na and Fe also show a tilt similar to what is seen for C and CO. The radial distribution is quite different, however, with the density being higher closer to the star instead of peaking at 85 au. The asymmetry is seen in both the inner and outer parts of the disk; concerning the inner distribution of Na and Fe, the emission is traced inward to the observational limit of 13 au in the NE, while to the SW the density peak appears to be located much farther out, beyond 50 au. In the outer parts of the disk, Na and Fe can be seen all the way to the limit of the observation at the projected distance of 330 au in the NE, while the disk seems truncated in the SW at 150–200 au.

A possible scenario is that the origin of CO (and thus C and O as its dissociation products) is different from the origin of the metals observed in the optical (Na, Fe, Ca, Ni, Ti, and Cr; Brandeker et al. 2004). While CO likely comes from the disruption of CO-rich cometary bodies at 85 au, the metals observed in the optical could be produced by the so-called falling evaporating bodies (e.g., Vidal-Madjar et al. 2017) close to the star and then diffused outward by the stellar radiation pressure. As the presence of C⁺ would act as a braking agent (Fernández et al. 2006), its spatial distribution would naturally become imprinted on the spatial distribution of the metals. The spatial distribution of Na and Fe is therefore consistent with C⁺ being in an eccentric distribution resulting from a tidal disruption event, as outlined in Sections 4.2.3 and 5.5.6. With the SW clump at 85 au being at the convergence point for a family of eccentric orbits, the C would be more concentrated to the clump location in the SW, while being much more radially distributed in the NE, in agreement with the observed Na and Fe distributions.