ALMA OBSERVATIONS OF THE DEBRIS DISK AROUND THE YOUNG SOLAR ANALOG HD 107146

The CASA task clean was used to Fourier invert the complex visibilities to create an image of the dust emission and to perform multi-frequency synthesis deconvolution. Figure 1 presents the ALMA map of the λ1.25 mm continuum emission from HD 107146 obtained with natural weighting using a Briggs robust parameter of 2. This weighting scheme gives a synthesized beam of 115 × 084, corresponding to a spatial resolution of about 32 × 23 AU at the distance of HD 107146. Figure 2 shows the same ALMA data but weighted with a Gaussian uv-taper with an on-sky FWHM of 2 arcsec. This suppresses the weight of the longer baselines, with the aim of highlighting the more diffuse emission from the debris disk. The angular resolution of the map of 252 × 225 is about twice as large as the map obtained with natural weighting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The integrated flux density from the disk is 12.5 ± 1.3 mJy at 1.25 mm. This was obtained by integrating the surface brightness over the area showing emission at ≳ 2σ above the background level in the outer taper weighted map. The uncertainty reflects the 10% of uncertainty on the absolute flux scale, whereas the 1σ noise levels on the natural and outer taper weighted maps are 0.030 and 0.057 mJy beam⁻¹, respectively. This flux density is consistent at ≈1σ with the measurement of 10.4 ± 1.4 mJy by Corder et al. (2009) using CARMA observations at a slightly longer effective wavelength of about 1.32 mm.

Figure 3 shows the radial averaged surface brightness profile obtained from the ALMA image with natural weighting. Consistent with previous results from CARMA (Corder et al. 2009) and the SMA (Hughes et al. 2011), the dust emission from HD 107146 extends over an angular diameter of ∼11'' with a decrease in the amount of emission toward the star. While the surface brightness of the inner and outer disk are similar, the ALMA observations reveal a decrease in the surface brightness at intermediate radii. These characteristics cannot be caused by spatial filtering from the interferometer. The shortest baseline in the ALMA data is 15 m, which corresponds to angular size scales of 17'', a factor of ∼1.5 larger than the disk. Also, the primary beam size of the ALMA 12 m diameter antennas is 24''at a wavelength of 1.3 mm. Although these observations are likely not fully sensitive to emission from the largest angular scales of the debris disk (see Wilner & Welch 1994), the fact that the recovered flux is consistent with that measured by Corder et al. (2009) with baselines shorter by a factor of ≈2.5 indicates that the loss of flux in our data is only marginal. Nonetheless, to avoid possible biases due to an incomplete (u, v) coverage of the observations, our analysis of the disk structure was performed on the measured visibilities (see Section 4).

Zoom In Zoom Out Reset image size

We imaged the CO J = 2–1 line using different weighting schemes, each providing different synthesized beams and sensitivities to surface brightness. A range of velocities of ±10 km s⁻¹ from the known radial velocity of the star (2.6 km s⁻¹; Valenti & Fischer 2005) was examined. No detection was obtained in any of these attempts by integrating over the disk size as seen in dust continuum. The tightest upper limit to surface brightness, and therefore to column density in CO, that could be derived from our data is about 10 mK (at 3σ) with an approximate angular resolution of 6''. A discussion on the upper limit on the mass in CO derived by this non-detection is presented in Section 7.

Our models consider a debris disk as an axisymmetric, geometrically thin and optically thin layer of solids. At any radial location in the disk r, the temperature T(a, r) of a grain with size a is derived from the balance between the energy absorbed by the grain and the thermal energy emitted by the grain itself:

$$\begin{eqnarray} && \pi a^2 \int _{0}^{\infty } \frac{L^{\star }_{\nu }}{4\pi r^2} Q^{\rm {abs}}_{\nu }(a) d\nu \nonumber\\ && \quad = 4 \pi a^2 \int _{0}^{\infty } Q^{\rm {abs}}_{\nu }(a) \pi B_{\nu } [T(a;r)] d\nu. \end{eqnarray} \tag{ 1 }$$

In Equation (1), $L^{\star }_{\nu }$ is the stellar monochromatic luminosity, $Q^{\rm {abs}}_{\nu }(a)$ is the frequency dependent absorption efficiency of a grain with size a, and B_ν[T(a; r)] is the Planck function at the temperature of the grain. After defining the stellar and dust properties, through $L^{\star }_{\nu }$ and $Q^{\rm {abs}}_\nu (a)$, respectively, Equation (1) can be solved numerically to derive the temperature T(a; r). In the case of HD 107146, $L^{\star }_{\nu }$ was derived from a black body with effective temperature T = 5841 K and bolometric luminosity Log(L_bol/L_☉) = 0.04 (Carpenter et al. 2008; Hillenbrand et al. 2008). The absorption coefficients $Q^{\rm {abs}}_\nu (a)$ were extracted from the dust models of Ricci et al. (2010a), which consider spherical porous solids made of a mixture of silicates, carbonaceous materials, water ice (optical constants from, Weingartner & Draine 2001; Zubko et al. 1996; Warren 1984, respectively,) and vacuum. These are based on the Pollack et al. (1994) models for dust in protoplanetary disks, and derived by an analysis of astronomical data and theory, as well as the composition of primitive bodies in the solar system.

An implicit assumption used to derive Equation (1) is that the whole disk is radially optically thin to the stellar radiation, so that stellar photons can reach any region in the disk without significant absorption. This assumption was verified a posteriori on the disk models that can reproduce the ALMA data for HD 107146.

The radial profile of the surface brightness of the disk as viewed as face-on is given by

$$\begin{equation} I_\nu (r) = \frac{\Sigma (r)\int _{a_{\rm {min}}}^{a_{\rm {max}}} B_{\nu }[T(a;r)] n(a)a^3 \kappa ^1_{\nu }(a) da}{\int _{a_{\rm {min}}}^{a_{\rm {max}}} n(a)a^3da}, \end{equation} \tag{ 2 }$$

where Σ(r) is the radial surface density of dust, n(a) is the grain size distribution, and $\kappa ^1_{\nu }(a)$ is the single-grain dust opacity coefficient that is related to the absorption efficiency through $(3/4a\rho) Q^{\rm {abs}}_{\nu }(a)$, with ρ = 1.2 g cm⁻³ being the mean solid density for the Ricci et al. (2010a) dust model considered here (see Miyake & Nakagawa 1993).

Sub-millimeter continuum observations can be used to constrain the size distribution of dust particles in debris disks. The grain size distribution is commonly approximated as a power-law function dN = n(a)da, n(a)∝a^−q for grains between a minimum and maximum grain size a_min and a_max, respectively. For grain size distributions that follow a power law with 3 < q < 4 over a broad enough interval around the observation wavelength, i.e., a_min ≪ λ ≪ a_max, Draine (2006) derived a relation between q, the spectral index β of the dust opacity (κ_ν∝ν^β), and the spectral index β_s of small grains (a ≪ λ), β_s = 1.8 ± 0.2 for interstellar grains. The assumption that the smallest and largest grains in the dust population are much smaller and larger than λ ≈ 1 mm, respectively, is in line with the interpretation of debris dust produced by a collisional cascade of large planetesimals getting ground all the way down to ∼μm sized grains. Also, physical models of collisional cascades of planetesimals predict q-values between about 3 and 4 (Dohnanyi 1969; Pan & Schlichting 2012; Gaspar et al. 2012), so that we can use the Draine (2006) equation to derive an estimate for q after inferring β (see Ricci et al. 2012).

The spectral index β of the dust opacity can be constrained by measuring the spectral index α of the millimeter spectral energy distribution (F_ν∝ν^α). By combining the flux densities measured for the HD 107146 debris disk at 0.88 mm, 1.25 mm, 1.32 mm, and 3.0 mm by Hughes et al. (2011), this work, Corder et al. (2009) and Carpenter et al. (2005), respectively, we measured a millimeter spectral index α = 2.42 ± 0.16.

In the case of HD 107146, the results of our modeling presented in Section 4 show that the dust emission at these wavelengths is optically thin and in the Rayleigh–Jeans regime and therefore the spectral index of the dust opacity is given by β = α − 2 (see, e.g., Ricci et al. 2010a). We then obtained q = 3.25 ± 0.09 from the Draine (2006) relation.

The integrals in Equation (2) are formally computed over all sizes of solids in the disk. A good approximation for the minimum grain size a_min is the blow-out grain size a_blow, as grains smaller than a_blow are blown out of the system by the stellar radiation field as soon as they are created. We estimated a_blow ≈ 1.7 μm using Equation (2) from Roccatagliata et al. (2009) with an albedo of ≈0.5 from our dust model, and a stellar mass of 1.0 M_☉ for HD 107146. With the equation used by Ricarte et al. (2013), which is valid for a silicate dust particle on a circular orbit around the central star, we would get a_blow ≈ 2.7 μm. However, our analysis is insensitive to this variation for a_blow since the emission at the ALMA wavelength of 1.25 mm is dominated by grains much larger than few μm.

The largest solids in a debris disk are kilometer-sized planetesimals or even planet-sized objects if these are present in the system. However, for practical reasons, the integrals in Equation (2) can be computed up to just a few centimeters as the emission from these and larger grains is negligible at the wavelength of our observations. For this reason we adopted a_max = 2 cm.

After adopting a parameterization for the dust surface density Σ(r) (see Section 5), Equations (1) and (2) define the debris disk model. The free parameters of our models are the parameters used to describe Σ(r), the disk inclination i defined as the angle between the disk axis and line-of-sight direction (0^o for face-on geometry), the disk position angle (P.A.) defined as the angle east of north to the disk major axis, and the offsets Δα and Δδ of the disk center relative to the phase center of the ALMA observations.

The models were fitted to the complex visibilities rather than the images to avoid non-linear effects of the image deconvolution, correlated noise between image pixels, as well as possible filtering out of the disk large scale emission. For a given set of disk model parameters, an image of the dust emission was generated at the mean frequency (239.5 GHz) of the observations. The image was then multiplied by the primary beam response of a 12 m ALMA antenna assuming a Gaussian-shaped beam with an FWHM of 243. The model image has a size of 2048 × 2048 pixel with 0015 pixels. The model visibilities were obtained through a Fourier transform of the image and sampling the model using the same (u, v) data points in the ALMA observations. The flux density of the sampled visibilities were varied with the frequency of the spectral bands as ν^{2 + β}, where β = 0.4. The value of χ² for the model parameters can then be computed.

The best-fit model parameters and parameter uncertainties were found using emcee, an MIT licensed pure-Python implementation of the Goodman & Weare (2010) Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble sampler⁵ (Foreman-Mackey et al. 2013). Relative to the traditional methods of χ²-minimization over a fixed multi-dimensional grid, the emcee algorithm has the advantage of focusing on the regions of the parameter space around the χ² absolute minimum. The models which were accepted by the MCMC Ensemble sampler were used to probe the probability function for each model parameter obtained through marginalization, i.e., by integrating the posterior distribution over all the parameters except the one of interest. The best-fit values and (asymmetric) uncertainties presented for each parameter in this work reflect the mode of the marginalized probability distribution and the 68.3% confidence level, respectively.

For each surface density parameterization, emcee was run with several hundreds walkers (see Foreman-Mackey et al. 2013) for several hundreds iterations. We also ran several trials by either varying or not varying the starting positions of the walkers to confirm that the final constraints on the model parameters are not affected by the randomly Gaussian distributed walkers and/or changes in the starting values of the parameters.

The simplest parameterization we considered to describe Σ(r) is a power law truncated at an inner and outer radius, R_in and R_out, respectively,

$$\begin{equation} \Sigma (r) = \left\lbrace \begin{array}{@{}ll}\Sigma _{\rm {0}}\left(\frac{r}{1\,\rm {AU}}\right)^{p} & \quad\; {{\rm where}\; R_{{\rm in}} < r < R_{{\rm out}}}\\ 0 & \quad\; {{\rm elsewhere},} \end{array} \right. \end{equation} \tag{ 3 }$$

where r is the distance in the disk from the central star, Σ₀ is the surface density normalized at 1 AU,⁶ and p is the slope of the radial power law. The model thus contains eight free parameters: {Σ₀, p, R_in, R_out, i, P.A., Δα, Δδ}.

Table 2 reports the constraints on the model parameters derived by our analysis (uncertainties are at 1σ). The top row in Figure 4 shows the probability distribution for R_in, R_out, p, as derived using the emcee algorithm (Section 4.2). With these models, the disk extends from about 25 to 152 AU. These estimates for inner and outer radii are both smaller than the values of 50 and 170 AU, respectively, derived by Hughes et al. (2011) from the combined analysis of the spectral energy distribution and SMA interferometric data at 0.88 mm using single power-law models, but are still within the uncertainties, reported as larger than 10 AU in their analysis.

Zoom In Zoom Out Reset image size

Table 2. Constraints on the Model Parameters from the Analysis of the ALMA Visibilities for the Three Classes of Models Discussed in Section 5

Parameter	Single Power Law	Single Power Law with Gap	Double Power Law
Rin [AU]	$25.2^{+2.7}_{-1.8}$	$30.3^{+2.4}_{-0.9}$	$35.8^{+1.3}_{-1.8}$
Rgap [AU]	...	$80.9^{+1.8}_{-2.6}$	...
ΔRgap [AU]	...	$9.0^{+1.0}_{-1.5}$	...
Rout [AU]	$152.0^{+1.0}_{-1.8}$	$150.4^{+1.7}_{-1.1}$	$148.3^{+2.1}_{-0.7}$
Log[Σ0/(g cm−2)]	$-5.58^{+0.13}_{-0.09}$	$-5.23^{+0.20}_{-0.06}$	$-2.33^{+0.21}_{-1.21}$
p	$0.74^{+0.06}_{-0.05}$	$0.59^{+0.04}_{-0.09}$	...
p1	...	...	$-1.09^{+0.74}_{-0.11}$
p2	...	...	$1.37^{+0.08}_{-0.20}$
Rbreak [AU]	...	...	$68.5^{+6.2}_{-2.6}$
i [deg]	20.6 ± 1.9	$20.7^{+2.2}_{-2.3}$	$21.0^{+2.0}_{-1.9}$
P.A. [deg]	$144.7^{+7.1}_{-4.2}$	$141.1^{+7.9}_{-3.3}$	$143.0^{+7.4}_{-3.6}$
Δα [arcsec]	$0.048^{+0.039}_{-0.027}$	$0.029^{+0.021}_{-0.039}$	$0.046^{+0.035}_{-0.023}$
Δδ [arcsec]	$-0.123^{+0.048}_{-0.031}$	$-0.081^{+0.046}_{-0.020}$	$-0.054^{+0.030}_{-0.037}$

Download table as:  ASCIITypeset image

The p value of $0.74^{+0.06}_{-0.05}$ indicates a surface density that is increasing with radius, i.e., p > 0. This is in line with the estimate of Σ(r)∝r^{0.3 ± 0.3} found by Hughes et al. (2011), although their large uncertainty did not allow them to rule out a flat or decreasing radial profile with their data. This result can be understood by noticing that at our wavelength of 1.25 mm I_ν∝B_ν[T(r)] × Σ(r) ∼ T(r) × Σ(r) ∼ r^−0.5 × Σ(r), and that as seen in Figure 1, the surface brightness toward the inner and outer edges of the disk is roughly the same. The offsets Δα and Δδ of the disk center relative to the phase center are consistent, within ≈2σ, with the offsets of the star (Section 2). This is true also for the other classes of models discussed in the next sections.

The value of the minimum of the χ² function is 2519918.4, which corresponds to a reduced $\widetilde{\chi }^2$ of 1.028. The image on the top left corner in Figure 5 contains the synthetic map of the best-fit model, which shows how that model would have looked like if it was observed under the same conditions (i.e., same (u, v) coverage, sensitivity) as our actual ALMA observations, with imaging performed using natural weighting. The contours in the right top corner map in the same figure show the residuals of the data—model subtraction. Although for the majority of the disk surface the residuals are below the 2σ level, some peaks and dips at ∼2σ–3σ are seen. In particular, nearly all these peaks (shown as blue continuous lines in the map) are located close to either the inner or the outer edge of the disk where the disk surface brightness is high, whereas the dips (blue dashed lines) are more toward the middle of the ring where the disk surface brightness is lower (see color map in the top right panel in Figure 5). This behavior is also evident in the radial profile of the disk surface brightness shown in Figure 3. This suggests that a surface density radial profile that can account for a depletion region at intermediate radii between R_in and R_out will likely better reproduce the ALMA data than a single continuous power law. In order to test this, in the next two subsections we present two simple extensions to the single power-law models.

Zoom In Zoom Out Reset image size

Figure 6 shows the comparison between visibility data measured by ALMA and binned over deprojected baseline lengths and the prediction of the single power-law model (blue line). The model shown here has parameters values corresponding to the peak of the probability distributions marginalized over each model parameter as probed by the emcee algorithm (Section 4.2). The small plot on the top of the figure highlights a region where the single power-law model does not reproduce well the real part of the visibility data. This occurs at deprojected baseline lengths of ≈75–130 m, or ≈60–110 kλ. These baseline lengths correspond to angular scales of ≈2–3 arcsec on the sky that is also the approximate range of angular distances from the central star of the decrease in surface brightness in the disk.

Zoom In Zoom Out Reset image size

A simple extension to the single power-law models that can account for a radial depletion of dust is provided by the following surface density functions:

$$\begin{equation} {\Sigma (r) = \left\lbrace \begin{array}{@{}ll}\Sigma _{\rm {0}}\left(\!\displaystyle\frac{r}{1\,\rm {AU}} \!\right)^{p} \!\!&\!\!\!\!\! \quad \;{{\rm where}\; R_{{\rm in}} < r < R_{{\rm gap}} - \displaystyle\frac{\Delta R_{{\rm gap}}}{2},\; {\rm or}\; R_{{\rm gap}} + \displaystyle\frac{\Delta R_{{\rm gap}}}{2} < r < R_{{\rm out}}}\\ 0 \!\!&\!\!\!\!\! \quad \;{\rm elsewhere.} \end{array} \right. } \end{equation} \tag{ 4 }$$

Equation (4) represents a single power law truncated at R_in and R_out as in Equation (3) with the addition of an axisymmetric gap with Σ(r) = 0 centered at R_gap and width ΔR_gap. This models has the same 8 parameters as the single power-law models (Section 5.1) plus two additional parameters that define the gap, namely R_gap and ΔR_gap, for a total of 10 model parameters.

The constraints to the model parameters are reported in the third column of Table 2. The probability distributions estimated for ΔR_gap, R_gap and p are shown in the middle row in Figure 4. The best-fit model presents a gap at a distance of about $80.9^{+1.8}_{-2.6}$ AU (at the center of the gap) and with a radial width $\Delta R_{\rm {gap}} = 9.0^{+1.0}_{-1.5}$ AU. The slope p of the power law of the surface density, i.e., $0.59^{+0.04}_{-0.09}$, is only slightly lower but still compatible at 2σ with the positive value found in the previous subsections.

The best-fit model has a χ² value of 2519875.4, which is lower by Δχ² = 43.0 than the case of a single power-law model with no gap. These models are an extension of the single power-law models presented in Section 5.1, which can be reproduced if ΔR_gap = 0 and/or R_gap = R_in or R_gap = R_out. In this case we can use the F −test to investigate whether this decrease of the χ² value with models with two extra parameters is statistically significant, i.e., if the best-fit model with a gap can be considered as a statistically better representation of the ALMA data than the single power-law models with no gap. The F-test statistics in this case has a value of 8 × 10⁻¹⁰, indicating that the models with gap provide a significantly better fit to the ALMA data. A similarly low value of about 3 × 10⁻⁹ for the relative likelihood of these models is derived using the Akaike information criterion (AIC). This is consistent with the fact that in our analysis about 99% of the models have ΔR_gap > 5.5 AU (see also Figure 4).

The fact that models with a gap better reproduce the ALMA data than single power-law models can be seen also by comparing the map of the data best-fit model residuals between these two classes of models (top two rows in Figure 5). The presence of a gap reproduces the decrease in surface brightness at radii roughly intermediate between the disk inner and outer radii. Also, the surface brightness of the best-fit model with gap is larger than in the no-gap case toward the disk inner and outer radii (see also Figure 3). As a consequence, lower absolute values for the residuals are seen in the residual map of the best-fit single power-law model with gap.

The better match between the ALMA data and the models with gap relative to the single power-law models is evident also in Figure 6. The model with gap (green line) fits well the data even in the region with deprojected baseline lengths ≈75–135 m, contrary to the single power-law models with no gap.

The models described in the previous subsection treat the case of a region in the disk that is fully depleted of dust, i.e., a gap. However, the dust depletion may be only partial, its density showing a local decrease but never reaching Σ(r) = 0 for R_in < r < R_out. A simple way to account for this possible behavior is by considering a double power-law radial profile for the surface density:

$$\begin{equation} \Sigma (r) = \left\lbrace \begin{array}{@{}l l}\Sigma _{\rm {0}}\left(\frac{r}{1\,\rm {AU}}\right)^{p_1} & \quad \;{{\rm where}\; R_{{\rm in}} < r < R_{{\rm break}}}\\ \Sigma _{{\rm break}}\left(\frac{r}{R_{{\rm break}}}\right)^{p_2} & \quad \;{{\rm where}\; R_{{\rm break}} < r < R_{{\rm out}}}\\ 0 & \quad \;{\rm elsewhere,} \end{array} \right. \end{equation} \tag{ 5 }$$

with the condition $\Sigma _{\rm {break}} = \Sigma _{\rm {0}} ({R_{\rm {break}}}/{1\,\rm {AU}})^{p_1}$ to assure continuity at r = R_break, and R_in ⩽ R_break ⩽ R_out.

This surface density is described by 6 free parameters, {Σ₀, p₁, p₂, R_break, R_in, and R_out}, for a total of 10 parameters defining this class of models. Also in this case these models represent an extension of the single power-law models discussed in Section 5.1: Equation (5) can "collapse" into Equation (3) if p₁ = p₂, and/or R_break = R_in or R_break = R_out.

The resulting constraints for the model parameters are reported in the fourth column of Table 2. The probability distributions estimated for R_break, p₁ and p₂ are shown in the bottom row in Figure 4. The best-fit model shows a negative value of p₁, indicating a surface density that is decreasing with radius out to $R_{\rm {break}} = 68.5^{+6.2}_{-2.6}$ AU, close to the value for the gap radius ($R_{\rm {gap}} = 80.9^{+1.8}_{-2.6}$ AU) for models with a gap. Beyond R_break the surface density increases radially with a power law of p₂ ≈ 1.4.

The best-fit model has a χ² value of 2519879.9. This value is nearly identical (difference of only 4.5) to the best-fit value for the single power law with a gap that has the same number of free parameters. The AIC test returns a relative likelihood of 0.10 between the single power law with a gap and the double power law. We conclude that the difference in the quality of the fit between the two models is not significant. Relative to the single power-law models with no gap, the χ² is lower by 38.5. The F-test statistics has a value of 7 × 10⁻⁹, and the relative likelihood of these models from the AIC test is 3 × 10⁻⁸. Therefore, like in the case of a gap, also the addition of a second power law to the surface density radial profile gives a statistically better representation of the ALMA data of HD 107146.

This is supported by the comparison of the residual maps in Figure 5. The residual map for the best-fit double power-law model is very similar to the case of a single power-law model with gap and significantly better than the single power-law model with no gap. The same occurs for the comparison of the surface brightness (Figure 3) and visibility predictions for these models (Figure 6).

We have also run models where a gap is added to the double power law, and a double power law without the condition of continuity at R_break. The general result is that for each functional form considered to analyze the continuum ALMA data of HD 107146, the analysis favors models that show inner and outer radii around ∼30 and ∼150 AU from the star, respectively, a depletion of dust at ∼70–80 AU, and a dust surface density toward the outer edge of the disk that is comparable or larger than the values found in the inner disk. The ALMA data presented in this paper do not allow us to distinguish between these possible different functional forms, and we decided to show in this paper the models with the lowest numbers of free parameters. This also means that current data do not permit a precise characterization of the radial profile of the dust surface density in the depletion region around ≈70–80 AU from the star. Future ALMA observations with better sensitivity and angular resolution than the ALMA Cycle 0 data presented here are needed for this.

As described in Section 5.3 the current data do not constrain the detailed radial profile of Σ(r) in the region of the depletion. In Section 5 we described the constraints to the model parameters in the case of single power-law models with a gap as well as for double power-law models. Both classes of models provide a better fit to the ALMA data than single power-law models with no gap. A visual representation of the constrained radial profiles of the dust surface density for each of the three classes of models is shown in Figure 7.

Zoom In Zoom Out Reset image size

If a gap is present in the HD 107146 disk then the most natural explanation would be the dynamical clearing by a planet. An alternative explanation involving "photoelectric instability," which has been proposed to reproduce rings in debris disks with detected gas (Lyra & Kuchner 2013) is made unlikely by our very tight upper limit to the CO-to-dust mass ratio (Section 7): this mechanism can be efficient only for disks with gas-to-dust mass ratios larger than ∼1, which would require a high value >10⁵ for the gas-to-CO mass ratio. The relatively old age of the star, i.e., ∼80–200 Myr, suggests that the gas-to-CO mass ratio in the system is probably lower than ∼10⁵–10⁴. This is the range of values estimated for the younger TW Hya primordial disk (Favre et al. 2013), where nearly all the gas is under the form of molecular H₂ that gets efficiently dispersed by the process of photo-evaporation within the first ∼10 Myr after the formation of the star. However, the gas-to-CO ratio may be underestimated if molecular CO has been strongly depleted, either because of condensation onto dust grains and/or UV photodissociation (e.g., Dent et al. 1995; Kamp & Bertoldi 2000). The temperature derived by our models for small dust particles across the disk is always ≳ 30 K, indicating that CO condensation onto grains should not be effective, at least as long as grains are surrounded by a CO substrate (see the discussion in Najita & Williams 2005). As for the UV photodissociation, given the relatively late spectral type of HD 107146, most dissociating UV photons are expected to come from the interstellar radiation field rather than from the central star (Kamp & Bertoldi 2000; Greaves et al. 2000). Najita & Williams (2005) applied photodissociation models to HD 107146 and estimated a gas-to-dust mass ratio <1 even in the case of a significant amount of molecular H₂ still present in the system. However, the geometry and temperature of the putative gas in the debris disk is obviously uncertain and can potentially play an important role in those calculations.

In the scenario of a planet-induced gap, a measurement of the gap width provides an estimate for the mass of the putative planet. At a given distance from the central star, a more massive planet carves a wider gap (Quillen 2006), which is set by the region of overlapping mean motion resonances on either side of the planet (Wisdom 1980). When the resonance widths exceed the distance between them, particle motion becomes chaotic. If the solids in the disk are in nearly circular orbits and the planet has low eccentricity, then the chaotic zone boundary is given by da_z ≈ (1.3–2.0) × μ^2/7, where da_z is the difference between the zone edge semi major axis and that of the planet divided by the semi major axis of the planet, and μ is the ratio between the planet and stellar mass. The range of (1.3–2.0) for the possible normalization coefficients has been obtained by numerical simulations incorporating collisions in a diffusive limit (Quillen & Faber 2006; Chiang et al. 2009). However, Gladman (1993) has demonstrated that solids in the outer regions of a chaotic zone will stay bound in that region even with their chaotic dynamics. Solids which are closer to the planet, within a region called crossing zone, will instead cross the orbit of the planet and get scattered out very efficiently. The size of the crossing region, which should be more directly related to the notion of a gap in solids around the planet, is given by da_z ≈ (2.1–2.4) × μ^1/3.

Under the framework of the models with the surface density described as a single power law with gap, da_z ≈ ΔR_gap/(2R_gap). The mass of the putative planet in HD 107146 can therefore be estimated using the values of ΔR_gap and R_gap derived by our analysis and reported in Table 2. We obtain $6.1^{+2.1}_{-3.1} \,M_{\oplus }$ and $4.1^{+1.4}_{-2.1} \,M_{\oplus }$ for values of 2.1 and 2.4 for the normalization coefficient in the da_z − μ relation for the crossing zone, respectively. The uncertainties on these estimates are derived from propagating the uncertainties on ΔR_gap and R_gap. If we used the relation from the chaotic zone theory, we would obtain $5.3^{+2.2}_{-3.1} \,M_{\oplus }$ and $1.2^{+0.5}_{-0.7} \,M_{\oplus }$ for values of 1.3 and 2.0 for the normalization coefficient, respectively. A value of ≈1.9 M_⊕ is obtained from the results of recent numerical calculations which account for disruptive collisions other than the gravitational interaction between a planet and a ring of planetesimals (Nesvold & Kuchner 2014).

Although only more sensitive and higher angular resolution observations can possibly confirm the presence of a gap in the HD 107146 debris disk, this analysis shows the potential of probing (indirectly) with ALMA the presence of terrestrial planets (and more massive planets) at wide separations and embedded in planetesimal disks. Terrestrial planets are still well beyond the reach of direct imaging techniques. In the case of HD 107146, upper limits of ∼10–13M_Jup at distances ≳ 15 AU from the central star have been obtained by high-contrast imaging in the near IR (Apai et al. 2008; Metchev & Hillenbrand 2009; Janson et al. 2013).

An interesting result of our analysis is that, at least in the outer regions of the HD 107146 debris disk, the dust surface density increases with the distance from the central star. A similar behavior has been recently found for the debris disk around the M-type, ≈10 Myr-old young star AU Mic (MacGregor et al. 2013). In the case of HD 107146, single power-law models with or without a gap predict power-law indices for the dust surface density radial profile of ≈0.57 and 0.75, respectively, with uncertainties lower than 0.1. In the case of double power-law models, the power-law index beyond R_break ≈ 70 AU is about 1.3. This behavior is different from all the younger gas-rich primordial disks in which the dust surface density decreases with radius⁷ (Andrews & Williams 2007; Andrews et al. 2009; Isella et al. 2009; Guilloteau et al. 2011).

The millimeter-sized grains traced by ALMA are insensitive to the stellar radiation pressure and trace the location of the planetesimals that produced them by collision. Since planetesimals are formed in primordial disks where surface densities and collisional rates decrease with radius, it is unnatural to think that this feature is an imprint of initial conditions. Rather, it may reflect collisional evolution and depletion of planetesimals with radius.

In the models of icy planets formation by Kenyon & Bromley (2002, 2008) the dust in debris disks is efficiently produced once a population of Pluto-sized objects with radii of ∼1000 km or larger is formed in the disk. These bodies stir smaller leftover planetesimals to disruption velocities, so that the outcome of their collisions is a very efficient production of debris dust particles. Once ∼1–10 km sized planetesimals have been ground, the production of dust is inhibited and the dust production rate decreases significantly. Therefore, in these models dust traces the recent formation of a population of Pluto-sized objects triggering the collisional cascade.

These Pluto-sized objects are formed at different times across the disk, so that the spatial variation of dust surface density can probe how the formation of icy planets propagates throughout the disk. For collisional processes, the growth time of solids in the disk is t_growth∝P/Σ where P is the orbital period (Lissauer 1987) and Σ the surface density of solids. As P∝r^3/2, if the initial surface density Σ of solids decreases with radius (or increases with a profile shallower than a power law with index 3/2), then solids growth is faster closer to the star. For a debris disk extending between ∼30 and ∼150 AU with the mass in dust particles as estimated for HD 107146 at ∼100 Myr, the Kenyon & Bromley (2008) models predict a timescale for the formation of 1000 km planetesimals of just ∼10 Myr in the inner disk (30 AU), but as large as ∼1.4 Gyr in the outer disk (150 AU). Taken at face value, the HD 107146 system would have had enough time to form Pluto-sized objects in the innermost disk regions, but not in the outermost ones, the furthest 1000 km objects being formed at ∼57, 62, 78 AU for HD 107146 ages of 80, 100, and 200 Myr, respectively. If Pluto-sized objects are not formed in the outer regions, dust production via collisional cascade would not be efficient in those regions, contrary to our result of a dust surface density which increases with radius in the outer disk.

These estimates strongly depend on the initial conditions assumed for the distribution of solids at the beginning of the calculation.⁸ In particular, the timescale for the formation of Pluto-sized objects scales with the density of solids with a power of −1.15.⁹ The Kenyon & Bromley (2008) models adopt an initial surface density of solids with radii of ∼1–1000 m with a power law radial profile with index −3/2 and a surface density at 30 AU of ≈0.1–0.2 g cm⁻² comparable to the distribution of solids in the minimum-mass solar nebula model (MMSN; Weidenschilling 1977; Hayashi 1981). In order to decrease the timescale for the formation of Pluto-sized objects in the outer disk down to values lower than the estimated age of HD 107146, i.e., ≲200 Myr, the initial density has to be increased by a factor of ∼10. Two possible ways to achieve this are by increasing the solids density at all radii by this factor, or by making the radial profile nearly flat, i.e., Σ_i∝r^−0.1, while keeping the same density in the inner disk.

In the former case the total mass of these solids would be increased by a factor of 10, giving a mass in 1–1000 m bodies of ∼4.7 × 10³⁰ g ≈800 M_⊕, while in the latter the required mass is ∼2 × 10³⁰ g ≈340 M_⊕. These values are close to the masses in smaller solids (radii ≲ 1 cm) estimated for the brightest disks in nearby, ≈1 Myr old star forming regions (e.g., Andrews & Williams 2005, 2007; Isella et al. 2009; Ricci et al. 2010a, 2010b; Andrews et al. 2013). If at this early stage of disk evolution the gas-to-solid mass ratio is still similar to the value ∼100 as found in the interstellar medium (ISM), the required initial gas mass is in the range ≈0.10–0.25 M_☉. This roughly corresponds to the mass limit for having gravitational instabilities in a disk around a 1 M_☉ star such as HD 107146 (e.g., Boss 1997; Rice et al. 2005).

Our analysis shows that if the primordial disk around HD 107146 was initially very massive and close to the gravitational instability limit, models of solid growth predict that the formation of Pluto-sized objects could have had the time to propagate outward to the disk outermost regions. The inferred rising surface density of small solids could be the result of this propagation, the outer regions containing more dust grains than the inner ones because having initiated the collisional cascade more recently.

Kennedy & Wyatt (2010) use this idea of a "self-stirred" planetesimal belt to predict the surface density of solids across the disk. They divide their simulations in two cases, depending on whether the collisional timescale t_c is lower or larger than the time at which the stirring is initiated, t_stir. For models in which t_c ≪ t_stir, the evolution of solids is violent and leads to rising, very steep surface densities with p > 2, much larger than the value constrained for HD 107146. If instead t_c > t_stir, consistently with the case of HD 107146 after using the Kennedy & Wyatt fiducial parameter values and a stellar mass of 1 M_☉ (see Equation (10) in their paper), the evolution is slower and this results into a more shallower but still rising surface density radial profile. For the model shown on the left panel of Figure 2 in Kennedy & Wyatt (2010), at an age of 100 Myr the surface density shows a rather moderate increase with radius in the disk outer regions, which is more consistent with the results of our analysis on HD 107146. Another interesting aspect of these models is that right after the radius at which the surface density peaks, the surface density rapidly drops by orders of magnitude. This can explain the absence in Figure 5 of any significant residuals beyond the outer radius R_out of our models. A more thorough theoretical investigation adapted specifically to the characteristics of the HD 107146 system is needed to further test models of debris disks produced by the self-stirring of a planetesimal belt.

Another possible scenario to explain a rising surface density with radius involves dynamical interaction between a planetesimal belt and a fully formed planetary system. According to the Nice model, dynamical gravitational interaction between the giant planets and the young Kuiper belt has played a crucial role in shaping the current architecture of the solar system (Gomes et al. 2005; Morbidelli et al. 2005). After several hundreds of megayears from the dispersal of the gas-rich primordial disks, Neptune migrated into the Kuiper belt and dynamically excited the orbit of several Kuiper belt objects, possibly explaining the late heavy bombardment (LHB) of the Moon (Tera et al. 1974). Because of this interaction, a large number of planetesimals got scattered both inward and outward spreading out the distribution of mass. The rate at which planetesimals were scattered inward is much greater than the rate at which they were scattered outward. As a result, the surface density profile got asymmetric with radius, with the bulk of the mass being represented by a surface density with a relatively shallow and rising slope (Σ∝r^1.5) followed by a very steep falloff (Σ∝r⁻⁵) at larger radii (Booth et al. 2009).

The scenarios described here is only one of the infinite possible outcomes of a planetary system–planetesimal belt interaction. It is referred to the case of the solar system, and therefore does not reproduce in detail the characteristics observed for the HD 107146 disk (see Section 8.3). At the same time this example shows the potential of this mechanism to possibly reproduce a relatively shallow surface density of small particles that is increasing with radius. Whether this scenario can explain the main properties of the HD 107146 remains to be tested. Models that account for both the dynamical interaction between a planetesimals belt and a planetary system as well as the debris dust production through the collision of large bodies are needed (e.g., Nesvold et al. 2013). The presence of planets in the system will be investigated through future high angular resolution observations (see discussion above).

We note that an important underlying assumption of our models is that the dust opacity is the same throughout the disk. The result of a dust surface density that is increasing with radius could be modified if the dust opacity κ_ν at the wavelength of the ALMA observations were a function of the location in the disk. Since I_1.25mm(r)∝T(r) × Σ(r) × κ_1.25mm(r), the terms Σ(r) and κ_1.25mm(r) would be degenerate. This means that our ALMA data could in principle be reproduced also by disk models with a κ_1.25mm(r) function that is increasing with radius and with a dust surface density that is flat or even decreasing with the distance from the central star.

Spatial variations of the dust opacity could reflect variations of the dust chemical composition and/or grain physical morphology, as well as variations in the grain size distribution. For example, if the slope q of the grain size distribution were changing from 3.0 in the inner disk to 3.5 in the outer disk, that would produce a variation of a factor of ∼3 in κ_1.25mm. That is close to the variation of the dust surface density between the inner and outer regions of the disk, as constrained by our models assuming a constant κ_1.25mm (see Figure 7). The only way to break the degeneracy between Σ(r) and κ_ν(r) and better constrain the radial profile of the dust surface density is by observing the disk at high angular resolution at multiple wavelengths in the (sub-)millimeter, as done for primordial disks (Isella et al. 2010; Guilloteau et al. 2011; Perez et al. 2012; Trotta et al. 2013). Future ALMA observations at multiple wavelengths have the potential to do this in the case of HD 107146.

The HD 107146 main sequence star has the same spectral type and luminosity class as the Sun. Its younger age allows us to investigate the properties of a planetesimal belt around a young Sun-like star.

While thermal emission from solids in extra-solar planetesimal belts trace small dust particles with sizes ≲ 1 cm, emission from dust produced in the Kuiper belt has not been detected yet and only large, kilometer-sized bodies can be seen in the Kuiper belt. The only way to attempt a comparison between the solar system's Kuiper belt and extra-solar belts is via an observational characterization of large solids and small dust particles in the former and latter systems respectively, using models for the dynamics and physics of collisions to connect the information on solids with very different sizes.

The vast majority of objects discovered in the solar system with a greater average distance than Neptune ("trans-Neptunian objects," TNOs) have a semi-major axis between ≈30 and 50 AU. These TNOs constitute the so-called "classical Kuiper belt." Whereas the inner radius of the classical Kuiper belt is very similar to that constrained for HD 107146, the outer disk spreads out a factor of ∼3 further.

Several TNOs have also been found beyond the classical Kuiper belt, as far as ∼1000 AU from the Sun, therefore at distances even larger than the outer radius of the HD 107146 debris disk. These objects are characterized by relatively large eccentricities and inclinations to the ecliptic. They have likely been scattered out through gravitational interactions with the giant planets rather than being formed in situ, and for this reason are called "scattered disk objects" (SDOs).

Vitense et al. (2010) analyzed a large sample of TNOs known at the time and used an algorithm to remove the observational biases due to the inclination and distance selection effects and estimated the main parameters of the Kuiper belt. They derived a mass of ≈0.12 M_⊕, half of which located in the classical Kuiper belt, the other half in SDOs. In order to predict the spatial distribution of dust particles in the presumed Kuiper belt debris disk, they applied to their "de-biased" Kuiper belt a numerical code that includes the effects of stellar gravity, radiation pressure, and Poynting–Robertson drag, as well as disruptive and erosive collisions. They derived a surface density of dust with a peak around ≈40 AU with a shallow radial profile in the inner disk down to radii ≲ 20 AU, and a steeper fall-off with a power-law index of ≈ −2.0 out to radii of several hundreds AU.

In addition to having a very different radial profile of the dust surface density, the debris disk around HD 107146 has a very different mass than the estimated Kuiper belt debris disk. The mass of kilometer-sized bodies in the Kuiper belt estimated by Vitense et al. (2010) is within a factor of two from the mass in dust estimated for the HD 107146 debris disk. Since the mass of kilometer-sized bodies that generate the dust is much larger than the mass of the dust itself, this implies that the dust in the HD 107146 is much more massive than in the Kuiper belt, which is in line with the much younger age of the HD 107146 system.

Quantifying the mass of the planetesimal belt surrounding HD 107146 from the mass in dust is prone to large uncertainties because of the huge step in solid size. If we could extrapolate the grain size distribution with a power-law index of 3.25 as derived in Section 4.1 all the way to bodies with sizes of 100 and 1000 km, we would find M_belt ∼ 10⁴ and 5 × 10⁴ M_⊕, respectively. These values are unreasonably high, as they would require a mass in gas and solids for the parental primordial disk larger than the mass of the central star (assuming the canonical value of 100 for the gas-to-dust mass ratio), much larger than the masses estimated for primordial disks around ∼1 Myr old pre-main-sequence stars (Andrews & Williams 2005; Isella et al. 2009; Ricci et al. 2010a).

Lower, more realistic values are derived by considering a power law index of 3.6 for the size distribution of large bodies as derived by Vitense et al. (2010) in the case of the Kuiper belt. With this value the mass in bodies as large as 1000 km would be 100 M_⊕. This is only a factor of a few larger than the mass in planetesimals invoked by the Nice model to describe the scattering event that lead to the LHB and to the dispersal of the vast majority of planetesimals in the early solar system (e.g., Booth et al. 2009).