The First Scattered-light Image of the Debris Disk around the Sco–Cen Target HD 129590

Observing the youngest stellar systems, shortly after planet formation has ceased, provides a glimpse of nascent circumstellar environments. However, only a small handful of stars (≲3–5) with ages ≲10–20 Myr, just after dissipation of the gaseous primordial disk, are found within 100 pc (Sacco et al. 2014). The Scorpius–Centaurus association (hereafter Sco–Cen; de Zeeuw et al. 1999) is the nearest OB2 association, with a mean distance of 140 pc, making it perhaps the most promising collection of young stars that can be observed shortly after the period of active planet formation. This region, containing stars with ages ∼10–16 Myr (Pecaut & Mamajek 2016) allows the best constraints to be placed on the orbital zones of planet formation, and on the early thermal histories of young planets (Ireland et al. 2011; Janson et al. 2013; Lafrenière et al. 2014; Hinkley et al. 2015).

Furthermore, such young systems often possess bright circumstellar debris disks, belts of planetesimals analogous to the Edgeworth–Kuiper Belt in our own solar system. The presence of dust in these systems suggests that planetesimals are responsible for the dust generation (Wyatt 2008). The dust that is generated through planetesimal collisions in these debris disks is inherently transient, being either blown out by stellar winds, or spiraling toward the host star via Poynting–Robertson drag. The persistence of dust in these systems implies that it is constantly being regenerated. Although populations of planetesimals may stir themselves in some cases (e.g., Kennedy & Wyatt 2010), this dust regeneration may be enhanced by perturbations from massive planets, dynamically exciting the planetesimals onto eccentric orbits, thus causing them to collide. As well as revealing the presence of massive planetary perturbers, the location of these large quantities of dust may hint at the location of giant planets in the system: gaps between dust belts may highlight where planets lie (Su & Rieke 2014), and sharp edges to debris rings can constrain the masses of planets shepherding these edges (e.g., Quillen 2006; Chiang et al. 2008; Mustill & Wyatt 2012).

One Sco–Cen system with a known circumstellar debris disk is HD 129590 (HIP 72070, ${T}_{{\rm{eff}}}=5945$ K, 1.3 M_⊙ 2.8 L_⊙; Chen et al. 2011). While not originally cataloged in the de Zeeuw et al. (1999) or Rizzuto et al. (2011) catalogs of Sco–Cen stars, HD 129590 has been listed as a G1V member of the Sco–Cen subgroup Upper Centaurus Lupus (Hoogerwerf 2000; Chen et al. 2011). HD 129590 has an estimated distance of 141 ± 7 pc (Astraatmadja & Bailer-Jones 2016; Gaia Collaboration et al. 2016a, 2016b) and a high infrared luminosity of ${L}_{{\rm{IR}}}/{L}_{{\rm{star}}}$ ∼ 5 × 10⁻³, twice the value observed for β Pictoris (Jang-Condell et al. 2015; Mittal et al. 2015), as can be seen in Figure 1. Although Chen et al. (2014) fit the SED as two distinct belts, Jang-Condell et al. (2015) suggest the system likely contains a single belt of dust. This is in agreement with the predictions of Ballering et al. (2013; see Section 4.1 for further details).

Zoom In Zoom Out Reset image size

The high fractional luminosity of HD 129590 makes it an extremely promising target for scattered-light imaging of circumstellar material (Currie et al. 2014; Draper et al. 2016). ALMA data have previously been obtained for this object (Lieman-Sifry et al. 2016), with the disk being marginally resolved along the major axis. Lieman-Sifry et al. (2016) constrain the inclination angle to >50° with a best-fit value of 70° and the position angle to $-{{59}_{-12}^{+17}}^\circ$. The ALMA data find a best-fit grain size of ${3.2}_{-0.5}^{+0.6}$ μm, and at this size find the outer edge to be ${110}_{-30}^{+50}$ au and the inner edge to be < 40 au. No gas was detected in that work, and so we do not expect gas to have a strong influence on the dust dynamics in this disk.

In this Letter, we present the first scattered-light image of the debris disk around HD 129590, using the SPHERE instrument on the VLT. The observations and data post-processing are presented in Sections 2 and 3, respectively. In Section 4, we describe our modeling of the disk, where we use an optically thin disk model to conclude that the disk has a radius of ∼60–70 au. Finally, we conclude in Section 6.

HD 129590 was observed on 2016 May 4 with the SPHERE instrument at the VLT (Beuzit et al. 2008), as part of a larger planet-finding survey within Sco–Cen. The data were taken in the IRDIFS mode, whereby light is split through a dichroic beamsplitter, and passed simultaneously to both the differential imager and spectrograph instrument (IRDIS; Dohlen et al. 2008) and the integral field spectrometer (IFS; Claudi et al. 2008). A total of 2560 s on-sky integration was collected by each instrument, with the N_ALC_YJH_S coronagraph in place. The sequence consisted of 80 × 32 s individual exposures for the IRDIS data and 40 × 64 s exposures for the IFS data. We used IRDIS in dual-band imaging (DBI; Vigan et al. 2010) mode with the H23 filter pair, at wavelengths λ = 1588.8 nm, ${\rm{\Delta }}\lambda =53.1$ nm and λ = 1667.1 nm, ${\rm{\Delta }}\lambda =55.6$ nm, while the IFS instrument (Zurlo et al. 2014; Mesa et al. 2015) was used in the YJ mode, which spans the range 0.95–1.35 μm. Platescales of the instruments are 12.25 mas/pix for IRDIS and 7.46 mas/pix for the IFS (Maire et al. 2016), and the inner working angle of the coronagraph is 015.

In addition to these science observations, flux calibration images were collected with the target displaced from the coronagraph. Star center calibration frames (waffle frames) were also collected, by imposing a sinusoidal pattern on the deformable mirror. This creates four starspot images, with equal displacements from the central star, in each corner of the frame. Together, these allow the star position to be accurately measured to ∼0.1 pixels (1.2 mas; Vigan et al. 2016) behind the occulting mask. The observations were carried out in pupil-stabilized mode to allow angular differential imaging analysis (ADI; Marois et al. 2015). The entire sequence of observations, including acquisition and calibration, lasted 59 minutes spanning an airmass range of 1.037–1.042. The primary science frames covered a total field rotation of 36°, and included the meridian crossing of the target.

Our data post-processing was carried out following the process described in Vigan et al. (2015, hereafter V15). We used both the ESO data reduction and handling pipeline (DRH; Pavlov et al. 2008) and the publicly available code described in V15, as well as some custom routines.

3.1. IFS

For the IFS data, basic calibrations were first created using the DRH: dark fields, master flat-fields, IFS spectra positions, initial wavelength calibrations, and an IFU flat were all generated. We then used a custom routine to calculate an accurate parallactic angle and time for each image, and to normalize the data based on its exposure time. Bad pixel and cross-talk corrections were applied (for details see V15). The DRH was then used to interpolate these frames both spectrally and spatially. A sigma-clipping routine was applied to remove remaining bad pixels which deviated from their neighbors by more than 3.3σ. Finally, the wavelengths for each image were recalibrated, due to small systematic errors in the DRH pipeline, as described in V15. This process results in a set of calibrated images in the x–y plane, at 39 wavelengths spanning 0.95–1.35 μm and at 40 distinct timesteps. No frame selection was performed, although each individual frame was visually inspected to confirm there were no data issues.

Speckle subtraction was performed using a custom PCA code (e.g., Amara & Quanz 2012; Soummer et al. 2012). This code simultaneously uses the spectral and temporal (parallactic angle) diversity of speckles to remove starlight scattered within the image plane by the telescope optics, but not genuine astrophysical sources. A disk feature was revealed, as shown in Figure 2. We tested reductions with between 2 and 100 principal components subtracted, and the feature is robust to the number of principal components removed.

Zoom In Zoom Out Reset image size

3.2. IRDIS

4.1. SED Fitting

4.2. Spatial Constraints

4.3. Disk Modeling

We then use injection modeling to characterize the disk more rigorously. We choose to take a Bayesian MCMC fitting approach, as performed in Wahhaj et al. (2014).

Synthetic disk images are created using the GRaTeR radiative transfer code (Augereau et al. 1999). We use an optically thin disk model with density defined as

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{\rho }_{\circ }\exp \left({\left[\tfrac{-| z| }{{\xi }_{\circ }}{\left(\tfrac{r}{{r}_{\circ }}\right)}^{-\beta }\right]}^{\gamma }\right)}{\sqrt{{\left(\tfrac{r}{{r}_{\circ }}\right)}^{-2{\alpha }_{\mathrm{in}}}+{\left(\tfrac{r}{{r}_{\circ }}\right)}^{-2{\alpha }_{\mathrm{out}}}}},\end{eqnarray} \tag{ 1 }$$

where r is the radial distance from the disk center, z is the vertical distance from the disk midplane, and ${\rho }_{\circ }$, ${\alpha }_{\mathrm{in}}$, ${\alpha }_{\mathrm{out}}$, ${\xi }_{\circ }$, β, and γ are free parameters. The exponential term defines the disk vertical profile, while the denominator expresses a radial profile that rises as ${\alpha }_{\mathrm{in}}$, peaks at r_◦ and then falls as ${\alpha }_{\mathrm{out}}$.

At each point in the disk, the scattered-light contribution is calculated as

$$\begin{eqnarray}&&F\propto \displaystyle \frac{\rho \times p(\theta )}{{d}^{2}},\end{eqnarray} \tag{ 2 }$$

where θ is the scattering angle, namely, the angle through which light from the star is scattered so that it reaches the Earth. d is the distance from the star to the grid point. The measured flux is influenced by several factors, such as the stellar luminosity, stellar distance, and the telescope gain. We fold these into a single overall normalization by modifying the ${\rho }_{\circ }$ parameter and denote this updated parameter as ${\rho }_{\circ }^{\prime }$.

The phase function, $p(\theta )$, is the standard Henyey–Greenstein scattering function:

$$\begin{eqnarray}&&p(\theta )=\displaystyle \frac{1}{4\pi }\displaystyle \frac{1-{g}^{2}}{{[1-2g\cos \theta +{g}^{2}]}^{\tfrac{3}{2}}}.\end{eqnarray} \tag{ 3 }$$

Since we assume an optically thin disk, scattered-light contributions are added along each line of sight to create a synthetic disk image.

These synthetic images were convolved with a Gaussian to mimic the effects of the point-spread function. Each model in turn was then rotated and subtracted from each frame of the raw IRDIS data, and the PCA sequence was repeated with six modes subtracted to generate a residual image. By injecting negative disk images, we accurately take into account the throughput of the PCA post-processing, which varies with separation from the star. This process is akin to a forward-modeling procedure. In the interest of computational efficiency, we use only a subset of the data for our modeling: images are trimmed to the central 240 × 240 pixels (∼3''), and only every fourth individual exposure is included.

The goodness of each model is assessed using the normal ${\chi }^{2}$ statistic: each pixel in the residual image is divided by its local sigma value and squared. Values for σ are calculated as in Wahhaj et al. (2014): we first convolve the PCA processed image of the disk with a Gaussian with an FWHM of 2 pixels. This is then subtracted from the image to remove extended spatial components, to leave an image containing only noise information. The rectangular test region is divided into annuli, and in each annulus the standard deviation of the noise image is found, so as to capture the variance of noise with distance from the star. As noted in Wahhaj et al. (2014), some disk signal still contributes to the standard deviation. This is inevitable, and will lead to conservative error calculations on our parameters.

We initially use a downhill minimization routine to find a best fit. We then use these best-fit parameters to initiate a Metropolis Hastings MCMC, using the emcee.py package (Foreman-Mackey et al. 2013). The MCMC chain generates a sampling of our parameter space, with probabilities assigned as $\exp \left(\tfrac{-{\chi }^{2}}{2}\right)$. We use uniform priors in this work. 510,000 random samplings are generated for each fit, and the first 10,000 are discarded to ensure that the results are independent from our starting position. The best fit and error values are calculated from the remaining 500,000 samples. For each parameter, the best fit is given by the median value of the marginalized distribution and the 1σ uncertainties are chosen to enclose 34% of the samples on either side of the median. In addition, we calculate the best-fit parameters for each third of the random samplings, and find them to be consistent. This confirms that the initial parameters do not affect the results, and that they are a genuine sample of the probability distribution.

We initially fit a single ring of dust. For this model, we have five free parameters: the overall scaling of the model ${\rho }_{\circ }^{\prime }$, the forward-scattering parameter g, the radius ${r}_{\circ }$, the inclination angle, and the position angle of the disk. We fix the radial profile to ${\alpha }_{\mathrm{in}}=3.5$ and ${\alpha }_{\mathrm{out}}=-3.5$ and fix the stellar distance to 141 pc. The parameters defining the vertical profile, namely, ${\xi }_{\circ }$, β, and γ, are fixed to values of 1, 0, and 2, respectively.

This fit is presented in the top row of Figure 3 with parameters shown in Table 1. The residual image for this model shows a clear halo outside the ring, and as such we fit the data again, but this time we also fit the values of ${\alpha }_{\mathrm{in}}$ and ${\alpha }_{\mathrm{out}}$. The best fit in this case is a smaller ring, with softer power-law edges, most notably on the outside, where ${\alpha }_{\mathrm{out}}\,=-{1.313}_{-0.012}^{+0.011}$. This is a surprisingly low value: Thébault & Wu (2008) find typical cases of either "smooth edges" with ${\alpha }_{\mathrm{out}}-3.5$ or "sharp edges" with ${\alpha }_{\mathrm{out}}$ as steep as −8. This second fit is shown in the lower panel of Figure 3, with parameters given in Table 1. For both of these models, the Henyey–Greenstein parameter takes a relatively high value (0.52 and 0.43, respectively). Although it has been shown for cases with a wide range of viewing angles that a single component Henyey–Greenstein parameter gives a poor fit (Stark et al. 2014; Hedman & Stark 2015), the geometry of HD 129590 means that this parameter is poorly constrained: the faint edge is severely affected by speckle noise, meaning backscattering is hard to constrain.

Zoom In Zoom Out Reset image size

Table 1.  Fit Parameters for the Two Disk Models in Figure 3

Parameter	Fit 1	Fit 2	ALMAa
${\rho }_{\circ }^{\prime }$	1.106 ± 0.006	1.440 ± 0.008	⋯
g	0.522 ± 0.002	0.4272 ± 0.0012	⋯
${\text{}}{r}_{\circ }$[au]	73.3 ± 0.2	59.3 ± 0.2	⋯
${\text{}}{i}_{\mathrm{tilt}}$[°]	76.87 ± 0.05	74.56 ± 0.05	>50
PA[°]	121.58 ± 0.02	121.80 ± 0.02	${121}_{-12}^{+17}$
${\alpha }_{\mathrm{out}}$	−3.5	$-{1.313}_{-0.012}^{+0.011}$	⋯
${\alpha }_{\mathrm{in}}$	3.5	3.15 ± 0.03	⋯
${\xi }_{\circ }$	1	1	⋯
β	0	0	⋯
γ	2	2	⋯
${d}_{\mathrm{star}}[\mathrm{pc}]$	141	141	⋯

Note. Values for ${\rho }_{\circ }^{\prime }$ are relative. PA is measured counterclockwise of north.

^aALMA fit parameters are taken from Lieman-Sifry et al. (2016) and converted to our reference system.

Download table as:  ASCIITypeset image

There is some residual structure that we are unable to model with GRaTeR. In particular, we find no notable improvement in the fit when we vary the disk eccentricity or offset from the star, and the data do not place meaningful constrains on the vertical fit profile. Our position angle is in very good agreement with that found by Lieman-Sifry et al. (2016) with ALMA observations, and our inclination angle is within the ALMA constraints.

The latest generation of dedicated high-resolution exoplanet imaging platforms such as GPI and SPHERE have already revealed a number of scattered-light debris disk images. Sco–Cen has proved to be a particularly fortuitous region for these searches. Currie et al. (2015) detected a dust ring around HD 115600, which was shown to be eccentric. HD 110058 has a wing-tilt asymmetry (Kasper et al. 2015), while HD 106906 appears to be misaligned with the wide planetary companion (Kalas et al. 2015; Lagrange et al. 2016). HD 111520 (Draper et al. 2016) has a dramatic brightness asymmetry, while both HIP 67497 (Bonnefoy et al. 2017) and HIP 73145 (Feldt et al. 2017) have been shown to have multiple, separate rings of debris. Perhaps the only debris disks without complex morphology are HD 114082 (Wahhaj et al. 2016) and HD 129590.

All of these debris disks are presented in Chen et al. (2014) as hosting multiple distinct debris belts. Multiple debris belts have only been seen in scattered light around two of the above targets: either there are several very small, close in belts of dust evading detection, or two temperature disks correspond to two belt disks less often than expected. Additionally, all of the targets with debris disk images have very high excess infrared luminosity values, and all except HIP 67497 have been observed with ALMA (Lieman-Sifry et al. 2016).

As discussed in Section 4.1, the debris belt around HD 129590 is predicted to lie at ∼16 au, based on the infrared excess temperature and simple blackbody constraints. Our observations show the ring to peak at ∼4 times this separation, implying an abundance of small dust at higher temperatures than the blackbody temperature. This is to be expected for a luminous star (2.8 L_⊙) where the blowout size is a few microns.

We attempt to classify HD 129590 under the categories outlined in Lee & Chiang (2016). The pronounced difference in brightness between the front and back of the disk is reminiscent of the "moth" or "double wing" structures simulated in that work, but we do not observe the bright, extended wings seen in, e.g., HD 61005 (Hines et al. 2007). A strongly forward-scattering disk could produce a fainter wing structure, as light will be preferentially scattered away from the viewer. In addition, an ADI based code such as that used here will self-subtract flux in regions at small angular separations to the bright disk edge, and so hide faint wing structures nearby. In Figure 4, contours of equal brightness are plotted, to demonstrate the full dynamic range of the disk. There is no indication of an extended wing, and indeed the self-subtraction lobes (highlighted in red in the image) show a high degree of symmetry. As such, there is no evidence for moth-like wings, and it is not possible to make detailed inferences about the presence, or otherwise, of a planet in this system.

Zoom In Zoom Out Reset image size

Our modeling work highlighted the presence of extended emission at large semimajor axis. This could be a dust halo, caused by radiation pressure blowing out small grains. Alternatively, this may be a scattered disk of planetesimals, and as such more closely resemble the wings mentioned above—albeit much more compressed into the plane than Lee & Chiang (2016) find. Such a ring of planetesimals might explain our surprisingly low value for ${\alpha }_{\mathrm{out}}$. A higher-resolution ALMA image could differentiate these two scenarios, while future imaging with a space facility such as the Hubble Space Telescope would allow the lowest surface brightness material (such as a dust halo) to be imaged, as well as defining the outer disk edge far better than ground-based AO imaging. Polarimetric differential imaging, meanwhile, would place constraints on the dust grain size and scattering properties, and may even reveal the faint edge of the disk.

HD 129590 is a G1V member of the Sco–Cen association, with an infrared excess twice that observed for β Pictoris. This work presents the first scattered-light images of the debris disk responsible for this infrared excess. The debris disk is revealed to be a nearly edge-on disk, with evidence for inner clearing. We use the GRaTeR radiative transfer code to model the disk as an optically thin ring, inclined to the line of sight by ∼75°. Our best-fitting model has a characteristic radius of ${r}_{\circ }=59.3$ au or ${r}_{\circ }=73.3$ au depending on the underlying model, and a forward-scattering parameter of g = 0.52 or g = 0.43. When the power-law edges were freed, these were found to take values of 3.15 inside the ring and −1.313 outside the ring. These values imply a strongly forward-scattering ring, with a soft outer edge. Even with this model, there is an indication in the final panel of Figure 3 of some residual structure, implying that there is some morphology more complex than a simple ring present in this disk.

Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 095.C-0549 and 097.C-1019. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. E.M. thanks the University of Exeter for support through a Ph.D. studentship. M.B. acknowledges support from DFG project Kr 2164/15-1. G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow.