The First Scattered Light Image of the Debris Disk around the Sco-Cen target HD 129590

We present the first scattered light image of the debris disk around HD 129590, a ~1.3 M$_\odot$ G1V member of the Scorpius Centaurus association with age ~10-16 Myr. The debris disk is imaged with the high contrast imaging instrument SPHERE at the Very Large Telescope, and is revealed by both the IRDIS and IFS subsytems, operating in the H and YJ bands respectively. The disk has a high infrared luminosity of $L_{\textrm{IR}}/L_{\textrm{star}}$~5$\times$10$^{-3}$, and has been resolved in other studies using ALMA. We detect a nearly edge on ring, with evidence of an inner clearing. We fit the debris disk using a model characterized by a single bright ring, with radius ~60-70 AU, in broad agreement with previous analysis of the target SED. The disk is vertically thin, and has an inclination angle of ~75$^\circ$. Along with other previously imaged edge-on disks in the Sco-Cen association such as HD 110058, HD 115600, and HD 111520, this disk image will allow of the structure and morphology of very young debris disks, shortly after the epoch of planet formation has ceased.


INTRODUCTION
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 disc, are found within 100 pc (Sacco et al. 2014). The Scorpius Centaurus Association (de Zeeuw et al. 1999, hereafter Sco-Cen) is the nearest OB2 association, with a mean distance of 140 pc, making it perhaps the most promising collection of young stars which 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 Janson et al. 2013;Lafrenière et al. 2011;Hinkley et al. 2015).
Furthermore, such young systems often possess bright circumstellar debris disks, belts of planetesimals anal-ogous 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 spiralling towards 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 et al. 2006;Chiang et al. 2008;Mustill & Wyatt 2012).
One Sco-Cen system with a known circumstellar debris disk is HD 129590 (HIP 72070, T eff = 5945 K, 1.3 M ⊙ 2.8 L ⊙ , Chen et al. 2011). While not originally catalogued 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 et al. 2000;Chen et al. 2011).
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 §4.1 for further details).
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 +17 −12 • . The ALMA data finds a best-fit grain size of 3.2 +0.6 −0.5 µm, and at this size find the outer edge to be 110 +50 −30 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 postprocessing are presented in §2 and §3 respectively. In section §4 we describe our modelling 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 §6.

OBSERVATIONS
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 2560s on-sky integration was collected by Spectrum of HD 129590.
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. 2006). The entire sequence of observations, including acquisition and calibration, lasted 59 minutes spanning an airmass range of 1.037 to 1.042. The primary science frames covered a total field rotation of 36 degrees, and included the meridian crossing of the target.
Our data post-processing was carried out following the process described in Vigan et al. (2015) (here on: 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.

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 normalise 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 neighbours 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. Soummer et al. 2012;Amara & Quanz 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 3. We tested reductions with between 2 and 100 principal components subtracted, and the feature is robust to the number of principal components removed.

IRDIS
Initial pre-processing of the IRDIS data was performed using the ESO SPHERE pipeline. Master dark and flat frames were created, and a waffle frame was used to calibrate the position of the star centers. Each frame was then independently reduced by applying the master dark and flat frames, and realigned using the star center calibrations and the dither positions so that the central star position was consistent between images.
The 160 individual images (80 timesteps, 2 wavelengths) were then input into a PCA algorithm (e.g. Soummer et al. 2012;Amara & Quanz 2012) to remove stellar speckle noise. For this process we used our own modification of the V15 IFS code, which was reconfigured to accept IRDIS data. A clear disk feature was observed, as demonstrated in Figure 3. This disk fea-ture closely matches that observed in the IFS data. As with the IFS data, we tested a range of reductions with between 2 and 100 principal components removed and found the disk feature to be robust.

SED fitting
The infrared excess of HD 129590 is well-studied. Ballering et al. (2013) fit a single cold dust component at a temperature of 89 K, while Jang-Condell et al.
(2015) find a best fit with grain temperature 93.7±0.1K. Chen et al. (2014), however, fit two separate components at 94 K and 72 K. The evidence for the second, cold belt comes from a single photometric point at 70 µm (see SED in Figure 1), which might be equally well explained by a dust model that is more complex than a simple blackbody. We conclude that the SED is best and most simply described by a single blackbody, with a best-fit temperature of 91 K (and radius 16 AU assuming blackbody absorption/emission). Small dust grains emit poorly at wavelengths significantly longer than their physical size, which we parametrise by modifying our blackbody fit by a factor (210µm/λ) β at wavelengths longer than 210µm (Wyatt 2008). This extra parameter is required by the ALMA observations, and β = 0.5 yields the best fit.

Spatial Constraints
The debris disk is highly symmetric in both the IFS and the IRDIS data. In the IRDIS data, there is a clear dark hole within 0.32 ′′ of the star. The brightest emission extends to 0.67 ′′ , but there is evidence of extended emission as far as 1.03 ′′ from the star, in line with the debris disk. Both lobes of the disk can be clearly seen. The disk signal in the IFS data is fainter, with only the front lobe visible. This has a spatial extent of 0.57 ′′ . As in the IRDIS data, there is extended, faint emission to the edge of the the field of view (0.87 ′′ ).

Disk Modelling
We then use injection modelling to characterise 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: where r is the radial distance from the disk center, z is the vertical distance from the disk midplane, and In each case, the image is a co-add of the entire wavelength range of the subsystem. The left hand images have 6 principal components subtracted, while on the right hand side a more aggressive reduction is presented, where 20 principal components have been removed. In both cases, the reduction is full-frame treatment, where the entire field of view is considered simultaneously. The ringed structure observed in the IFS images is an artifact of the fast Fourier transform process used to rescale the various wavelength observations. ρ • , α in , α out , ξ • , β and γ are free parameters. The exponential term defines the disk vertical profile, while the denominator expresses a radial profile which rises as α in , peaks at r • and then falls as α out . At each point in the disc, the scattered light contribution is calculated as: 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 ρ • parameter and denote this updated parameter as ρ ′ • . The phase function, p(θ), is the standard Henyey-Greenstein scattering function: Since we assume an optically thin disc, 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 modelling procedure. In the interests of computational efficiency, we use only a subset of the data for our modelling: images are trimmed to the central 240×240 pixels (∼3 ′′ ), and only every fourth individual exposure is in-cluded.
The goodness of each model is assessed using the normal χ 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 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 disc 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 parameters space, with probabilities assigned as exp −χ 2 2 . 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 of our starting position. The best fit and error values are calculated from the remaining 500,000 samples. For each parameter, the best fit given by the median value of the marginalized distribution and the 1-σ uncertainties 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 ρ ′ • , the forward scattering parameter g, the radius r • , the inclination angle and the position angle of the disk. We fix the radial profile to α in = 3.5 and α out = −3.5 and fix the stellar distance to 141 pc. The parameters defining the vertical profile, namely ξ • , β 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 α in and α out . The best fit in this case is a smaller ring, with softer power law edges -most notably on the outside, where α out = −1.313 +0.011 −0.012 . This is a surprisingly low value: Thébault & Wu (2008) find typical cases of either "smooth edges" with α out -3.5, or "sharp edges" with α 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.
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. Note-Values for ρ ′ • are relative. PA is measured anti-clockwise of North. *ALMA fit parameters are taken from (Lieman-Sifry et al. 2016) and converted to our reference system.

DISCUSSION
The latest generation of dedicated high resolution exoplanet imaging platforms such as GPI and SPHERE has already revealed a number of scattered light debris disc images. Sco-Cen has proved to be a particularly fortuitous region for these searches. Currie et al. (2015b) 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 The initial debris disk fit, where the Henyey-Greenstein parameter g, the disk radius, and the position and inclination angles of the disk have all been allowed to vary. A residual halo of light is observed. Bottom: Here the parameter αout has also been allowed to vary. As such, the residual halo is better modelled. Fit parameters for both models are given in Table 1. In both cases, IRDIS data is shown with 6 PCA modes subtracted. multiple, separate, rings of debris. Perhaps the only debris discs without complex morphology are HD 114082 (Wahhaj et al. 2016) and HD 129590. All of these debris discs 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 disc evading detection, or two temperature discs correspond to two belt discs less often than expected. Additionally, all of the targets with debris disc 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 §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 forwardscattering disk could produce only 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 angu-500mas=71au Figure 4. The data shown in the top left panel of Figure  3, but presented as a contour plot. Contours have relative brightness values of 1, 2, 4, 8 and 16, while negative contours are plotted in red, with relative flux -1 and -2. The self-subtraction wings above and below the bright ring edge are clearly visible. No attempt has been made to take the throughput of the PCA routine into account in this plot. lar separations to the bright disk edge, and so hide faint wing structures nearby. In Figure 5 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.
Our modelling work highlighted the presence of extended emission, at large semi-major 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 α 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 disc.

CONCLUSION
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 edgeon 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 • =59.3 AU or r • =73.3 AU depending on the underlying model, and a forward scattering parameter 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. EM thanks the University of Exeter for support through a Ph.D. studentship. MB acknowledges support from DFG project Kr 2164/15-1. GMK is supported by the Royal Society as a Royal Society University Research Fellow.