Resolving Structure in the Debris Disk around HD 206893 with ALMA

To measure the disk structure, we generated synthetic model images of debris disks with varying geometries that were compared with the interferometric data. For each model image, we compared the data with the model in the visibility domain and calculated a χ² metric to evaluate the goodness of fit between the model and the data.

To generate sky-projected model images of debris disks, we use a ray-tracing code described in Flaherty et al. (2015), ¹⁰ which is a Python adaptation of an earlier IDL code by Rosenfeld et al. (2013). The code translates parametric models of the density and temperature of dust and gas to a grid of density and temperature that is then rotated and integrated along the line of sight to produce a sky-projected map of the millimeter intensity. While there is a fully 3D version of the code available for eccentric disks, we use the 2D version that assumes a circular, axisymmetric flux distribution with a linearly increasing scale height, since there is no a priori evidence for deviations from axisymmetry like a clear offset of the disk center from the star position (and indeed, the assumption is verified a posteriori by the lack of residuals after the best-fit model is subtracted from the data). We assume a dust opacity of 2.3 cm² g⁻¹ (Beckwith et al. 1990), yielding optically thin emission for the dust in all of the models within the parameter space explored by the MCMC chains. Since density and temperature are degenerate for optically thin emission, we assume that these large grains are in blackbody equilibrium with the central star and calculate a dust grain temperature of

$$\begin{eqnarray}&&{T}_{\mathrm{dust}}={\left(\frac{{L}_{* }}{16\pi \sigma {R}^{2}}\right)}^{1/4},\end{eqnarray} \tag{ 1 }$$

where L_* is the bolometric luminosity of the star, σ is the Stefan–Boltzmann constant, and R is the distance of the dust grain from the central star. We adopt a stellar luminosity of L_* = 2.83 L_⊙ for HD 206893 (Delorme et al. 2017). We assumed several different possible functional forms for the surface density of the disk; the functional forms for the surface density are summarized in Table 2.

For example, the surface density in the disk follows a power-law distribution where

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{c}{R}^{p}\end{eqnarray} \tag{ 2 }$$

between radii of R_in and R_out, where Σ_c is a normalization for the total dust mass in the disk, M_disk (see Table 2 for definition), and p is a power-law index that we initially set to a value of 1.0. While the value of p is interesting from the perspective of debris disk evolution, there is a well-known degeneracy between p and the location of the outer radius R_out (see, e.g., Section 4.2.2. of Ricarte et al. 2013), and our data are of sufficiently limited sensitivity and angular resolution that we chose to focus on R_out rather than p. We chose a value of 1 as a middle-of-the-road estimate bounded by theoretical predictions of 0 for collisional evolution models that assume a collisional lifetime of the largest planetesimal longer than the age of the system (Schüppler et al. 2016; Marino et al. 2017), and 7/3 for those that assume the opposite (Kennedy & Wyatt 2010).

As it became clear that the local minimum and maximum evident in the radial surface density profile (Figure 2) would require a more complex radial surface density distribution than the initial flat-disk model, we explored two families of more complicated models: models with broken power laws that switched power-law indices at one or two transition radii within the inner and outer radius of the disk, and models that included a sharp radial gap in the surface density profile at a particular radius. When analyzing surface brightness features with modest S/N, it is often not possible to determine the exact shape of the feature, but broken power laws and sharp gaps are both common families of solutions assumed in the literature (see, e.g., Ricci et al. 2015) and adequately represent the extremes of an abrupt and deep gap versus a broad and shallow gap.

After generating the sky-projected images, we apply a primary-beam correction (multiplying by the primary beam) and then convert the model image into synthetic model visibilities using the MIRIAD task uvmodel (Sault et al. 1995), which we then compare with the data in the uv plane to calculate a χ² value as a goodness-of-fit test. Comparing data in the uv plane is desirable both because uncertainties are well characterized in the uv plane (unlike in the image domain, where the uncertainties are unknown and correlated between pixels, and nonlinearities can be introduced by the CLEAN process) and because it is agnostic to the choice of imaging parameters and allows us to take full advantage of the range of baseline lengths sampled.

We fit the models to the data using an affine-invariant MCMC sampler (Goodman & Weare 2010) implemented in Python in the software package emcee (Foreman-Mackey et al. 2013). The goodness of a fit between the synthetic and observed visibilities are evaluated by a log-likelihood metric $\mathrm{ln}{ \mathcal L }=-{\chi }^{2}/2$. The MCMC code directs an ensemble of walkers in an exploration of parameter space, according to the calculated probability that a given walker position (representing a single set of model parameters) provides a better fit to the data than the previous walker position. After a "burn-in" phase during which the walkers search downhill for the χ² minimum, the process results in a set of model parameters describing the different "steps" that each walker makes, which can then be agglomerated into a marginalized posterior probability distribution for each parameter (see Figures 7, 9, and 8).

We performed several MCMC runs in order to investigate a variety of model formalisms. Initially we varied eight parameters: the inner radius (R_in), the distance between the inner part of the disk and the outer part of the disk (ΔR, which is related to the disk outer radius as R_out = R_in + ΔR), the mass of the disk (M_disk), the flux density of the central star (F_*), the position angle of the disk major axis (PA), the inclination of the disk relative to the observer's line of sight (i), and the position offset in R.A. (Δx) and decl. (Δy) of the star−disk system relative to the pointing center of the interferometer. Subsequently, we introduced a gap inner radius (R_in,Gap) and width (ΔR_Gap), and for the power-law models, one or two transition radii (Rt1 and Rt2) with power-law indices for the disk segments between the transition radii (pp1, pp2, and pp3). All parameters were sampled linearly except for disk mass, which was sampled logarithmically, effectively equivalent to assuming a log-uniform prior, and position angle and inclination, which were sampled as $\cos i$ and $\cos \,\mathrm{PA}$ to avoid undersampling of the extrema.

The initial run revealed two issues that made us suspect that a more complex distribution of material was necessary: First, the stellar flux was approximately double the anticipated value based on the blackbody approximation, and the models showed a much sharper central peak than we saw in the data, suggesting that there must be some diffuse emission around the central star. In addition, the models seemed to explore two regions of parameter space: initially they explored a region of parameter space with a small inner radius, which had a more reasonable stellar flux but left a ring of negative residuals farther from the star, and then eventually they landed in a region of parameter space where the disk inner radius was larger than expected—too large for the inner edge to be carved by the brown dwarf, though it is certainly possible that other unseen companions could be carving the disk edge—and the stellar flux was too high. This behavior suggests that there might actually be two edges to the disk: one close to the star, accounting for the diffuse flux around the star, and another farther out. Therefore, we explored a model with a radial gap in the disk. We assumed uniform priors that required the inner radius of the gap to fall within the radial extent of the disk ([R_in, R_out]) and that required the width of the gap to be smaller than the total width of the disk ([0, ΔR]). The functional form of the gap is a top hat: we assumed that the gap was completely empty, and the edges of the gap are step functions.

In an attempt to remain agnostic about the functional form of the minimum in the surface brightness, we explored both a power law with an empty gap and a double power law, motivated by the set of functional forms assumed by Ricci et al. (2015). However, we found that the double power-law transition radius fell on the peak in the surface brightness around 115 au, which meant that we were not able to evaluate whether a break in the power-law surface density could reproduce the observed surface brightness profile as well as a disk with a gap. We therefore explored two additional profiles: a double power law with a radial gap, and a triple power law with two transitional radii. All of the models were consistent in preferring a dip in the radial surface density profile near 75 au and a peak near 115 au, and the comparison between the latter two functional forms demonstrates that an empty gap with sharp edges yields comparable results to a more shallow power-law inflection point and is preferred with modest statistical significance. Table 2 presents a summary of the functional forms, surface density normalization, free parameters, and best-fit ln prob values for each of the seven classes of models that we fit to the data. For the remainder of the paper, we focus on the comparison between the flat disk (which ignores the local maximum and minimum of the surface density), the double power law with a gap, and the triple power law, since the latter two were the models that best (statistically and by eye) reproduced the features of the observations.

The limits of the priors for all parameters are listed in Table 4. For the flat disk we used 16 walkers and ran the chain for 2000 steps beyond the burn-in period. For the double power law with gap we used 30 walkers and ran the chain for 2000 steps beyond the burn-in period. For the triple power law, we used 30 walkers and ran the chain for 3000 steps beyond the burn-in period. In all cases, the burn-in period was estimated by eye based on where the ln prob values leveled off to a relatively constant maximum. We also followed up with an autocorrelation analysis showing that while the autocorrelation time was still rising by the end of each chain, all parameters had leveled off so that the fractional error in the mean was no more than a few percent, and the standard error of the mean was stable. We show some sample plots from the double power-law model with a gap and the triple power-law model in the Appendix. The best-fit models, as well as the median and uncertainties given by the 16th and 84th percentiles of the posterior distribution, are presented in Table 3. Figure 3 shows the radially averaged visibility profile for the data, compared with the best-fit models for the flat disk, double power law with gap, and triple power law. The visibilities have been deprojected assuming an inclination of 44° and position angle of 60° (the best-fit values for the double power-law model).

Zoom In Zoom Out Reset image size

Figures 4–6 show the tapered data image (left) compared with the best-fit model (middle), sampled at the same baseline separations and orientations and imaged with the same parameters as the data, and the residuals (right). The lack of significant (>3σ) residuals for these models shows that all models adequately fit the data. The main difference visible in the best-fit models is the distribution of flux around the star, toward the center of the disk. For the model without a gap or power-law minimum, the stellar flux is higher (by a factor of two) and the distribution of flux in the center of the system is strongly and centrally peaked. For the models with a gap or power-law break, the stellar flux is lower (more in line with expectations from the blackbody estimate) and the emission is more diffuse around the star, since the model with a gap incorporates disk emission extending throughout the inner regions of the disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the models with a gap or power-law break used five more parameters than the model without a gap, we expect them to be able to provide a better fit simply owing to the larger number of degrees of freedom. In order to appropriately penalize the additional parameters when evaluating changes in the goodness-of-fit parameter, we used the AIC, a form of the Akaike information criterion, and the BIC, the Bayesian information criterion. The BIC penalizes the use of additional parameters more than the AIC, and it also takes into account the sample size (which the AIC does not). A BIC score of ΔBIC > 10 implies "very strong" evidence of a statistically improved fit (Kass & Raftery 1995).

When calculated, the ΔBIC comparing the double power law with a gap to the flat disk is 37.6, implying "very strong" evidence that the double power law with a gap is a better fit to the data than the flat disk. The AIC also returns a probability of 1.1 × 10⁻⁸ that the flat disk is a better fit than the model with a gap. The corresponding BIC and AIC values for the triple power-law gap compared with the flat-disk model are 44.2 and 3.1 × 10⁻⁷, also highly significant. With so many models, there are a large number of potential comparisons we could make, but it is perhaps worth noting that for any given model type (single power law, or double power law) adding a radial gap provides a significantly better fit than the same model without the gap, even taking into account the increase in the number of parameters (the ΔBIC value for the flat power law with a gap compared to the flat power law without a gap is 16, the value for the single power law with a gap compared to the single power law without a gap is 23, and the value for the double power law with a gap compared to the double power law without a gap is 5.9). These values corroborate the conclusion that we have statistically significant evidence for a radial gap in the disk.

As for the shape of the gap, the double power law with a gap and the triple power law both provide an excellent fit to the data using the same number of parameters, though the double power law with a gap does have the larger ln prob value. On the basis of the difference in ln prob, the probability that the double power law with a gap provides a better fit to the data than the triple power law is 0.03, which we take as suggestive but not conclusive evidence that a sharp, empty gap might provide a better fit to the data than a more shallow transition between two power laws.

The results of our MCMC fits (with posterior distributions shown in Figures 7–9) demonstrate that the HD 206893 disk is well described by a broad underlying flux distribution extending from radii of <51 au to ${194}_{-2}^{+13}$ au (this constraint on R_out is measured from the posterior distribution of R_out = R_in + ΔR, which is not shown in Figure 8), with a gap beginning at a radius of ${63}_{-16}^{+8}$ au and exhibiting a width of ${31}_{-7}^{+11}$ au. There are no statistically significant (>3σ) residuals, which indicates that the observations are well described by an azimuthally symmetric and circular distribution of flux. Given the limited S/N of the data (the brightest parts of the dust disk are approximately at the 5σ level on average), we can only rule out local (or point-to-point) spatial variations in flux of ≳50% on spatial scales comparable to the beam FWHM of ∼30 au or larger.

Zoom In Zoom Out Reset image size

When we compare our derived disk size to the sample of spatially resolved disks investigated by Matrà et al. (2018), we note that the value of R for HD 206893 according to the definition laid out in that paper would be approximately 115 au. Of their sample of 10 F stars, only one (HR 8799) has a larger radius—and that star has a particularly unusual λ Boo spectrum that is in some ways more consistent with an A star spectral type. The disk around HD 206893 is therefore among the largest known debris disks around F stars, although comparable to that around HD 170773 (Sepulveda et al. 2019).

Our model makes several simplifying assumptions that may or may not be realistic: the completely empty gap (rather than allowing the depth of the gap to vary), and the sharp edges of both the disk inner and outer radii and the inner and outer edges of the gap. Better characterizing these features—the depth of the gap, the sharpness of the edges—would be interesting from the perspective of understanding debris disk evolution and interactions between the disk and the planet. Unfortunately, this initial reconnaissance of the disk structure is necessarily limited in angular resolution and sensitivity, and the quality of the data does not justify a more complex model (as evidenced by the lack of residuals). However, armed with this new knowledge of the critical spatial and flux scales for the system, future observations at higher angular resolution and sensitivity will be able to elucidate these details.

One of the key parameters of interest is the inner radius of the disk. Previous models of the spatially unresolved SED had estimated that the location of the inner radius should fall around 50 au (Moór et al. 2011), estimated from the derived dust temperature of 49 ± 2 K, although comparisons between spatially resolved observations and SED estimates tend to reveal a large amount of scatter, with spatially resolved measurements biased toward larger radii than SED-based estimates (Booth et al. 2013; Pawellek et al. 2014; Morales et al. 2016). In this case, the measured radius is smaller than the SED-estimated radius and appears to be spatially unresolved with a radius of <51 au. This is good news for obtaining dynamical measurements of the mass of the brown dwarf: an inner radius of 50 au (as suggested by the SED models in Moór et al. 2011) is too large to be plausibly truncated by a brown dwarf with a projected separation of 10 au, at least in the absence of extreme eccentricity, which is ruled out by the combination of astrometric, radial velocity, and direct imaging constraints explored by Grandjean et al. (2019). We should note that there is some tension between a derived inner radius of <51 au and the 50 K temperature constraint from the SED (Moór et al. 2011), since even blackbody-like grains should be hotter than 50 K at distances interior to ∼40 au from the central star. The temperature probed by the SED is likely dominated by the bulk of the opacity, i.e., ∼micron-sized grains in a collision-dominated debris disk. Small micron-sized grains are expected to be located farther from the star than those of blackbody-like grains at the same temperature. The measured disk extent and the disk SED are largely consistent. However, the morphology of the ALMA data seems to suggest that at least some millimeter-size grains must be present within the central 1–2 beams, corresponding to ≲37 au in radius (assuming a beam FWHM of 06 and a diameter of 1.5 beams for the central flux component). Future investigation is needed, particularly the detailed shape of the infrared SED and the exact location of the inner disk edge.

The bad news is that the current data do not have sufficient angular resolution to allow us to provide constraints on the mass of the BD companion that are more restrictive than the constraints from Grandjean et al. (2019). Our best-fit value for the inner radius of the triple power-law model (35 au), when interpreted in the context of chaotic zone theory and assuming a semimajor axis of 13 au (from the MCMC fits to the combined astrometric, RV, and direct imaging constraints in Grandjean et al. 2019), provides an upper limit on the mass of the brown dwarf of <1170 M_Jup, assuming zero eccentricity. The extent of the chaotic zone depends on the semimajor axis of the companion's orbit, as well as the mass ratio μ between the companion and the central star. To derive the upper limit on the mass, we use the functional form derived by Morrison & Malhotra (2015) for the outer extent of the chaotic zone for a high-μ companion: Δa ≃ 1.7μ^0.31 a_p . It is possible to place somewhat better constraints through a more detailed dynamical analysis, although even that is difficult owing to the likely presence of an additional companion causing the RV drift detected by Grandjean et al. (2019). Marino et al. (2020) conduct just such an analysis and show that the location of the inner edge is consistent with being carved by the brown dwarf, though it could be farther out than expected based on chaotic zone theory owing to secular resonances with the companion responsible for the RV drift. Since the location of the inner edge of the disk is so steeply dependent on the companion mass, even a modest improvement in the angular resolution will yield substantially improved constraints on the dynamical mass measurement of the brown dwarf companion.

Another important constraint provided by the circumstellar disk measurement is the inclination of the planetesimal belt. While the nearly face-on configuration precludes meaningful constraints on the vertical structure, and therefore the dispersion of inclinations, the average inclination is constrained by our analysis to be 45° ± 4°. This range of inclinations is somewhat inconsistent with the plausible range of 20°–41° for the brown dwarf companion reported by Grandjean et al. (2019), as well as the 30° ± 5° inclination of the stellar pole derived by Delorme et al. (2017). Interestingly, the version of the Grandjean et al. (2019) dynamical MCMC analysis that combines the constraints from radial velocity, astrometry, stellar proper motion, and direct imaging results in a posterior distribution for the inclination of the companion of 45° ± 3° and a mass of ${10}_{-4}^{+5}$ M_Jup, though they note that the median χ² on the RV data is six times higher than when they do not include the stellar proper-motion constraints, indicating that there may be some inconsistency in the ability of their model to reproduce both the RV and stellar proper-motion data—perhaps providing additional support for the presence of another companion in the system. At this point, while we cannot rule out mutual inclination for the star and brown dwarf companion, coplanarity seems plausible.

The AIC and BIC comparison between the best-fit models provides strong evidence for a gap or local minimum in the surface density near a radius of 80 au. This marks the fourth detection of a radial gap in a broad debris belt, after HD 107146 (Ricci et al. 2015; Marino et al. 2018), HD 15115 (MacGregor et al. 2019), and HD 92945 (Marino et al. 2019). While the sample size is not large, there are only a couple of other debris disks with broad (ΔR/R ≳ 1) debris belts that have been imaged with sufficient resolution and sensitivity at millimeter wavelengths to detect substructure, including β Pic (Dent et al. 2014; Matrà et al. 2019) and AU Mic (MacGregor et al. 2013), the latter of which in fact does exhibit tentative (AIC 3.7σ, ΔBIC = 4.3) evidence for a gap in its edge-on planetesimal belt (Daley et al. 2019). The newly detected gap in the HD 206893 debris disk is therefore part of an emerging trend of gapped structure in broad debris disks, which will be exciting to confirm and explore with future high-resolution observations.

The cause of the gapped structure in debris disks is not well understood. The main categories into which proposed mechanisms fall include (1) gas-dust dynamics, (2) inheritance from the protoplanetary disk phase, and (3) dynamical interactions with a planet.

Lyra & Kuchner (2013) note that a robust clumping instability can organize dust into multiple eccentric rings even in the absence of a planet. The main prerequisite for the effect to operate is a gas-to-dust mass ratio ≳1. Given our upper limit on the CO flux density of 40 mJy km s⁻¹, the approximate upper limit on the CO mass, assuming LTE and an excitation temperature of 30 K, is roughly 8.4 × 10⁻⁵ M_⊕, using molecular data for CO v = 0 drawn from the Cologne Database for Molecular Spectroscopy (Endres et al. 2016). With the best-fit dust mass from Table 3 of 0.021 M_⊕, the upper limit on the gas-to-dust mass ratio is 1.0 × 10⁻⁴ if we assume that the gas is dominated by CO, or 0.060 if we assume a protoplanetary-like composition with a CO/H₂ abundance ratio of 10⁻⁴, corresponding to a mass ratio of ∼0.0014. Given these stringent limits, it is unlikely that gas dynamics are responsible for the double-ringed structure in this system.

One way in which gas dynamics could conceivably play a role despite the lack of gas in the present-day debris disk is if the gapped structure is inherited from the protoplanetary disk phase. The Disk Structures at High Angular Resolution Project (DSHARP) has characterized the locations of rings and gaps in a sample of 20 (18 of which are single) bright, nearby protoplanetary disks (Huang et al. 2018). The outer dust radius of HD 206893 would place it among the larger disks in the sample (3 of 20 DSHARP disks are as large or larger), and the central gap radius of 78 au would place it roughly in the 79th percentile of the 52 protoplanetary disk gaps identified by DSHARP. The structure of the HD 206893 system is therefore generally consistent with the distribution of disk radii and gap locations identified by DSHARP, though perhaps on the larger side especially given its spectral type. It is also worth noting that the DSHARP sample is biased toward systems with high dust luminosities and is almost certainly not representative. Of course, the origin of the gaps and rings in protoplanetary disks is also not well understood, so if the structure is inherited, then at most we have gained some evidence that the dust rings in DSHARP do correspond to planetesimal rings. Proposed mechanisms for generating the gaps in protoplanetary disks include chemical effects like snow lines (e.g., Banzatti et al. 2015; Okuzumi et al. 2016), magnetohydrodynamic instabilities of various flavors along with resulting gas pressure gradients (e.g., Johansen et al. 2009; Bai & Stone 2014; Simon & Armitage 2014; Flock et al. 2015; Lyra et al. 2015), and, of course, one or more (proto)planets (e.g., Papaloizou & Lin 1984; Bryden et al. 1999; Nelson et al. 2000). At this point, there is enough evidence pointing to the presence of planets in the disk at young ages to make it likely that at least some of the gaps and rings in protoplanetary disks are indeed caused by planets (e.g., Isella et al. 2018, 2019; Pinte et al. 2018; Teague et al. 2018).

The extent to which disk structure is likely to be retained between the protoplanetary and debris disk phase is still largely an open question. Most theoretical studies of debris disk structure assume that the protoplanetary disk dust profile is not retained and make predictions for debris disk structure that depend only on the locations of colliding planetesimals that produce the dust, along with the masses and orbital properties of the larger bodies that dynamically sculpt it (plus, for smaller grains, radiation pressure, stellar winds, and sometimes Poynting–Robertson drag, depending on the collision rate).

In the absence of a reservoir of gas whose mass is comparable to that of the dust, it is difficult to break the radial symmetry of the disk without planets. Even within the category of planet-sculpted gaps, there are multiple theoretical considerations that point to different masses and orbital properties of the underlying planetary system. However, the most straightforward explanation is that a single additional unseen planet-mass companion, located within the gap, is responsible for clearing the dust at the location of the gap. In that case, the semimajor axis and width of the gap encode the location and mass of the planet. Ignoring the presence of the brown dwarf at 10 au and using the best-fit values for R_gap and ΔR_gap of 67 and 32 au, respectively, the three-body chaotic zone theory would predict that the gap is carved by a planet with semimajor axis a_p = 80 au and mass ratio μ = 1.0 × 10⁻³, corresponding to a mass of 1.4 M_Jup, assuming an orbital eccentricity of zero. To derive these estimates, we use the general functional form from Morrison & Malhotra (2015) where Δa = C μ^β a_p , where C and β are constants tabulated separately for the inner and outer chaotic zone width for a range of ratios of the radius of the planet R_p to the extent of the Hill sphere R_H; we use the values tabulated for R_p /R_H = 0.001. If the planet has significant eccentricity, then the width of the gap would be expected to increase as 1.8e^1/5 μ^1/5, at least above a critical eccentricity of 0.21μ^3/7 = 0.004, implying that a lower-mass planet on an eccentric orbit could be responsible for carving the observed 31 au gap width (Mustill & Wyatt 2012). Propagating uncertainties through to the mass estimate is nontrivial, not only due to the relative uncertainties on our measured values of the gap radius and depth but also due to uncertainties on the parameters C and β, the asymmetric nature of the chaotic zone, and the degeneracy between model parameters like p and the gap properties, so the value of 1.4 M_Jup should be considered a ballpark estimate subject to many systematic and relative uncertainties.

More exotic possibilities have been invoked to explain the gapped structure in the HD 107146 disk, notably the possibility that both the inner disk edge and the gap could be carved by a single planet with mass comparable to the debris disk and eccentricity between 0.4 and 0.5, located interior to the inner edge of the broad debris belt (Pearce & Wyatt 2015). That configuration was ultimately ruled out for the HD 107146 system by Marino et al. (2018). It is similarly unlikely to apply to the HD 206893 system, because the location of the inner disk edge appears to be adjacent to the chaotic zone of the directly imaged brown dwarf, and the brown dwarf's mass is so much larger than that of the disk that the apsidal anti-alignment described for low-mass companions in Pearce & Wyatt (2015) would not occur. Another possibility, of course, is that multiple planets could be carving the gap, in which case the mass of each planet would be substantially smaller than the 1.4 M_Jup that we derive for the single-planet case.

The chaotic zone theory on which we base our estimate of planet mass is fundamentally a three-body result, but adding the dynamical influence of the brown dwarf (effectively treating the central star as a low-μ binary system) makes the system fundamentally a four-body problem. We therefore conducted N-body simulations with and without the brown dwarf (the latter to verify that the gap depth and width were not significantly affected by the presence of the brown dwarf) to investigate whether the orbital properties of the brown dwarf have a detectable impact on the width of the gap.

In order to simulate gravitational interaction between the star, brown dwarf, putative planet, and dust particles in the disk, we used the N-body software REBOUND (Rein & Liu 2012) with hybrid integrator MERCURIUS (Rein et al. 2019), which switches from a fixed to variable time step when any particle comes within a certain distance of a massive particle, here chosen to be three Hill radii. We chose the fixed time step to be 4% of the brown dwarf's initial period, which is less than 8% of every dust particle's period when the simulation begins. Dust particles are treated in the test particle limit.

Particles are initially randomly distributed with a uniform distribution of semimajor axes between 8 and 158 au, a uniform distribution of eccentricity from 0 to 0.02, and a uniform distribution of inclination from 0° to 10°. We placed 10⁴ disk particles, adequate to sample the 150 au span of the disk and recover smooth images of the final disk density distribution. We assumed a stellar mass of 1.32 M_☉ and a stellar radius of 1.26 R_☉ (Delorme et al. 2017). The brown dwarf is placed on an orbit with semimajor axis 10.4 au and assumed to have a radius of 1.3 R_Jup. We assumed a planet mass of 1.4 M_Jup and a bulk density of 1.3 g cm⁻³. We placed the planet on a circular orbit at 80 au. We integrate the system up to 10 Myr, which is sufficiently long to capture many orbital timescales for all the bodies. We then assign a weight to each particle based on its initial semimajor axis to mock up a surface density that follows the best-fit double power-law model parameters. We ran a set of simulations, varying the brown dwarf mass between 12 and 50 M_Jup (the range of plausible values from Grandjean et al. 2019) with a 1.3 M_Jup planet. We also investigated the effects of placing the brown dwarf on orbits with eccentricities of 0.03, 0.1, and 0.3, as well as introducing an orbital inclination to the plane of the planetesimal disk of 5° or 10°.

Zoom In Zoom Out Reset image size

A summary of the properties of the brown dwarf and the measured gap width and depth based on the final particle distribution is presented in Table 5. The gap width and depth were estimated by fitting a power law to the initial surface brightness distribution, subtracting the final profile from the fit to the initial profile, excluding the "Trojan" particles at the center of the gap, and fitting a top-hat function to the result. Trojans were defined as any particle with a semimajor axis that falls within 1 Hill radius of the semimajor axis of the planet. Two sample radial surface brightness profiles are presented in Figure 10. The left panel shows the result for a low-mass (12 M_Jup) brown dwarf with zero eccentricity, and the right panel shows the result for a high-mass (50 M_Jup) brown dwarf with an eccentricity of 0.3. In both panels, the solid vertical lines represent the location of the brown dwarf (near 10.4 au) and planet (near 79 au), and the adjacent dashed lines represent the extent of each body's chaotic zone, as calculated from the formulae presented in Morrison & Malhotra (2015) for the zero-eccentricity case and Mustill & Wyatt (2012) for the high-eccentricity case. While the top-hat function fit to the surface brightness profile yields a slightly broader width of 34 au than expected from the chaotic zone extent of 32 au, the figures demonstrate that the gap edge is not well defined and the chaotic zone approximation is a good estimate for where the surface brightness starts to increase away from the influence of the companion.

Zoom In Zoom Out Reset image size

We therefore find that the N-body simulations confirm that the chaotic zone approximation is appropriate for most of the parameter space covered by current estimates of the properties of the brown dwarf companion. Even a high-mass brown dwarf (50 M_Jup) with substantial eccentricity does not significantly perturb the planet's orbit (and therefore the width and location of the gap estimated by chaotic zone theory). This result is consistent with a back-of-the-envelope estimate that places the secular timescale at <10 Myr and the maximum eccentricity of the outer planet at 0.06 for a brown-dwarf-to-planet mass ratio ≫1. The depth of the gap for the case of 15 M_Jup is roughly 94%, which is slightly shallower than the assumed 100% depth of the gap in our ray-tracing code. Assuming a shallower gap would likely lead to a broader estimated width of the gap in the MCMC code, indicating that our estimate of the planet mass may be slightly too low. However, since the N-body simulation includes Trojans the depth may or may not be realistic (for example, the model does not include collisional evolution or the formation of the planet's core, either of which might deplete particles within the gap).

Finally, we then fixed the mass of the brown dwarf at 15 M_Jup, on a circular orbit at 10.4 au, and ran a set of simulations varying the planet mass between 0.5 and 50 M_⊕. We also investigated the effects of placing the planet with eccentricities of 0.03, 0.1, and 0.3, as well as introducing an orbital inclination to the plane of the disk of 5°, 10°, and 15°. A summary of the properties of the planet and the measured gap width and depth based on the final particle distribution is presented in Table 6. The gap width calculated by the top-hat function fig continues to yield results higher than the chaotic zone estimate for all planet masses and eccentricities. However, the changes in gap width generally conform to chaotic zone theory, as gap width increases with planet mass to the 0.27 power, close to the formula from Morrison & Malhotra (2015). The gap width increases proportionally to eccentricity to the 0.21 power, close to the formula from Mustill & Wyatt (2012). As expected, as mass increased, the depth of the gap increased, as well as the width. As eccentricity increased, the depth of the gap decreased. Inclining the planet's orbit had little effect on the width of the gap, although the depth of the gap slightly decreased.