β PICTORIS' INNER DISK IN POLARIZED LIGHT AND NEW ORBITAL PARAMETERS FOR β PICTORIS b

Polarimetric observations of β Pic were carried out on 2013 December 12 UT, while performing a series of AO performance and optimization tests (Table 1). β Pic was observed for a total of forty-nine 60 s frames, during which the field rotated in parallactic angle by 91°. Between each image the half waveplate (HWP) modulator was rotated by 225. For 25 frames, GPI's two Sterling cycle cryocoolers (Chilcote et al. 2012) were set to minimal power to reduce vibration in the telescope and instrument to improve the AO performance. The external Gemini seeing monitors were not operational during these observations and as a result the seeing throughout the sequence remains unknown. However, β Pic b can easily be seen in the majority of raw detector images, even before data reduction.

Each raw data frame was dark-subtracted, corrected for bad-pixels and then "destriped" to remove any remaining correlated noise in the raw image caused the by the cryocooler vibration (Ingraham et al. 2014). Since the time of these observations, the level of vibration has been significantly mitigated through the use of a new controller, which drives the two coolers 180° out of phase (Hartung et al. 2014). The vibration caused by the two coolers now interferes destructively and the overall effect is significantly damped. With a reduced level of vibration, the destriping algorithm is only needed for very short exposure times, and incorrect use may result in the injection of unwanted noise. We direct the reader to the GPI IFS Data Handbook³¹ for further details on the destriping algorithm and its appropriate use.

In GPI's polarimetry mode (pol. mode), a Wollaston prism splits the light from each lenslet into two spots of orthogonal polarization states on the detector. Flexure effects within the instrument cause these lenslet PSF spots to move from their predetermined locations on the detector, typically by a fraction of a pixel. For each frame, the PSF offset was determined using a cross-correlation between the raw frame and a set of lenslet PSF models measured using a Gemini Facility Calibration Unit (GCAL) calibration frame. The overall method is decribed in Draper et al. (2014) using high-resolution microlenslet PSFs. Here we use a Gaussian PSF model, which is less computationally intensive, but provides similar results. The raw frames were then reduced to polarization datacubes (where the third dimension carries the two orthogonal polarizations) using a weighted PSF extraction centered on the flexure-corrected location of each of the lenslets' two spots (see Perrin et al. 2015).

Each cube was divided by a reduced GCAL flat field image, smoothed using a low pass filter. The flat field corrects simultaneously for throughput across the field and a spatially varying polarization signal. In theory, this polarization signal should be removed during the double differencing procedure later in the pipeline; however, we have found empirically that this polarization signal is best divided out of each cube individually.³²

To determine the position of the occulter-obscured star, a Radon-transform-based algorithm (Pueyo et al. 2015) was used to measure the position of the elongated satellite spots (Wang et al. 2014). Knowledge of the obscured star's location is critical when combining multiple datacubes that must be both registered and rotated. Each datacube was then corrected for distortion across the field of view (Konopacky et al. 2014). The datacubes were then corrected for any non-common path biases between the two polarization spots using the double differencing correction described by Perrin et al. (2015), before being smoothed with a 2-pixel FWHM Gaussian profile.

Instrumental polarization, due to optics upstream of the waveplate, converts unpolarized light from the stellar PSF into measurable polarization that, if left uncalibrated, can mimic an astrophysical signal. This signal was removed from each difference cube individually, first by measuring the average fractional polarization (i.e., the difference of the two orthogonal polarization slices divided by their sum) inside of the occulting mask, where the flux is due solely to star light diffracting around the mask. We assume that this fractional polarization signal is due to polarization of unpolarized stellar flux by the instrument and telescope. For each lenslet, the fractional instrumental polarization was then multiplied by the total intensity at that location, and then subtracted off in a similar manner to the double differencing correction (see Footnote ³²). Using this method we find the instrumental polarization to be ∼0.5%, a similar level to that reported by Wiktorowicz et al. (2014) using the same data set.

The difference cubes were then shifted to place the obscured star at the center of the frame and then rotated to place north along the y-axis and east along the x-axis. All of the polarization datacubes were then combined using singular value decomposition matrix inversion to obtain a three dimensional Stokes cube, as described in Perrin et al. (2015). Non-ideal retardance in GPI's HWP makes GPI weakly sensitive to circular polarization, Stokes V. Measurements of the circular polarization of an astrophysical source would require knowledge of the HWP's retardance beyond the current level of calibration. Therefore, in almost all cases the Stokes V cube slice should be completely disregarded.

The Stokes datacube was then transformed to "radial" Stokes parameters: $(I,Q,U,V)\to (I,{Q}_{r},{U}_{r},V)$ (Schmid et al. 2006; see Footnote ³²). Under this convention, each pixel in the Q_r image contains all the linear polarized flux that is aligned perpendicular or parallel to the vector connecting that pixel to the central star. A positive Q_r value indicates a perpendicular alignment and a negative Q_r indicates a parallel alignment. Note that this sign convention is opposite that used in Schmid et al. (2006), where positive values of Q_r correspond to a parallel alignment. The U_r image holds the flux that is aligned ±45° to the same vector. For an optically thin circumstellar disk, the polarization is expected to be perpendicular and all the flux is expected to be positive in the Q_r image. The U_r image should contain no polarized flux from the disk and can be treated as a noise map for the Q_r image. The final reduced disk image can be seen in Figure 1.

Zoom In Zoom Out Reset image size

Observations of β Pic in spectroscopic mode (spec. mode) were carried out during four separate GPI commissioning runs, as well as during an ongoing astrometric monitoring program scheduled during regular general observing time. In total, we present nine individual sets of observations over seven epochs (Table 1). Two of the observation sets have been previously published: the H-band data set from 2013 November 18 (Macintosh et al. 2014) and the J-band data set from 2013 December 10 (Bonnefoy et al. 2014). Here we have re-reduced all the data in a consistent manner in an effort to reduce systematic biases and maximize the homogeneity of the data set. As with the polarization mode observations, those observations that were taken during the instrument's commissioning were carried out during AO optimization tests and therefore have a range of exposure times and filter combinations.

All data sets were reduced with standard recipes provided by the GPI DRP. Raw data frames were dark subtracted and destriped for microphonics in the same manner as the polarimetry observations. A short-exposure arc lamp image was taken contemporaneously with each science observation to measure the offsets of the lenslet spectra due to flexure within the IFS. The mean shift was calculated for a subset of lenslets across the field of view relative to a high signal-to-noise ratio (S/N) arc lamp image taken at zenith via a Levenberg–Marquardt least-squares minimization algorithm (Wolff et al. 2014).

The raw detector image was then transformed into a spectral datacube, using a box extraction method. For observations obtained with the K1 and K2 filters, thermal sky observations were taken immediately before or after the observation sequence. Sky background cubes were created in the same manner described above and subtracted from science datacubes. Finally, all cubes were corrected for distortion (Konopacky et al. 2014).

Each data set was PSF subtracted using the methods outlined in Pueyo et al. (2015). To minimize systematic biases, the ensemble of data sets was treated as uniformly as possible. The main steps of this data reduction process include: high-pass filtering, to remove the remaining PSF halo; wavelength-to-wavelength and cube-to-cube image registration, to correct for atmospheric differential refraction and sub-pixel stellar motion across the observing sequence; subtracting the speckles using the KLIP principal component analysis algorithm (Soummer et al. 2012) on each wavelength slice in each cube; rotation to align the north angle of each image; and co-adding the resulting cubes in time.

For the epochs in which β Pic b was observed on consecutive nights, relative alignment was tested using both the cross-correlation method described in Pueyo et al. (2015) and the absolute stellar locations based on the satellite spot centroids derived using the GPI DRP. For these epochs we found better consistency in the location of β Pic b using the DRP centroids, which we then chose to adopt for all data sets.

The KLIP algorithm was implemented using both spectral differential imaging (SDI; Marois et al. 2000) and ADI, building for each slice a PSF library that takes advantage of the radial and azimuthal speckle diversity (in wavelength and in PA, respectively). Due to the relative brightness of β Pic b with respect to the neighboring speckles we limited the exploration of KLIP parameter space to two zone geometries and two exclusion criteria (1 and 1.5 PSF FWHM) for each data set. For each slice, the 30 PSFs that were the most correlated in the region where β Pic b is located were used for PSF subtraction, except for the J-band data which required 50 PSFs for satisfactory subtraction.

To determine the optimal number of principal components to use for each data set, we examined both the evolution of the extracted spectrum and the astrometric stability as a function of wavelength as we increased the number of components. This latter test helps us to rule out the pathological cases for which either a residual speckle (i.e., insufficiently aggressive PSF subtraction) or self-subtraction (i.e., over-aggressive PSF subtraction) bias planet centroid estimates. We checked for potential biases by comparing astrometric positions measured when using only a high-pass filter with those measured when applying KLIP. Finally, we checked for self-consistency by injecting six synthetic point sources at the same separation as β Pic b, but at different PAs. Based on these tests we concluded that the astrometric measurements do not feature systematics either introduced by residual speckles or biases associated with KLIP above the uncertainty levels reported in Section 4.1.

The principal objective of our disk modeling is to retrieve basic geometric properties of the disk. The modeling approach adopted here is to combine a simple recipe for the 3D dust density distribution with a parametric model of the polarized scattering phase function and then fit to the data using the affine-invariant sampler in a parallel-tempering scheme from the emcee Markov chain Monte Carlo (MCMC) package (Foreman-Mackey et al. 2013). Parallel tempering uses walkers at different "temperatures" to broadly sample the posterior distributions and is therefore a useful strategy when the likelihood surface is complex.

For a disk seen in edge-on projection, the radial dust density distribution becomes degenerate with the dust scattering properties. This degeneracy is typically broken with the use of physical grain models, which describe scattering properties (including polarization) as a function of wavelength. In practice, observations are fit to grain models either using simultaneous polarization and total intensity information (e.g., Graham et al. 2007), or multicolor images (Golimowski et al. 2006). With only single wavelength polarized intensity images available, we instead use the HG scattering function (Henyey & Greenstein 1941) to describe the scattering efficiency as a function of scattering angle. The shape of the HG scattering function is a function of only one parameter, the expectation value of the cosine of the scattering angle, $g=\langle \mathrm{cos}\theta \rangle ,$ and thus provides a useful tool to approximate grain scattering when using physical models is impractical. The applicability of the HG scattering function to the modeling of our polarized intensity images is discussed in Section 5.1.

Our dust density model, expressed in stellocentric coordinates, $\eta (r,z),$ follows a power law between an inner radius, R₁, and an outer radius, R₂, and has a Gaussian vertical profile with rms height ${h}_{0}r$ and constant aspect ratio, h₀:

$$\begin{eqnarray*}&&\eta (r,z)\propto {\left(\displaystyle \frac{r}{{R}_{1}}\right)}^{-\beta }\mathrm{exp}\left[-\frac{1}{2}{\left(\frac{z}{{h}_{0}\ r}\right)}^{2}\right]\end{eqnarray*}$$

where r is the radial distance from the star, z is the height above the disk midplane and β is the power law index of the dust density. Inside R₁ and outside R₂ the dust density is zero. The dust density distribution is combined with the HG function, $H(\theta ,g),$ to generate a scattered light image of the disk as seen in 2D projection from the observer's frame, where the intensity for a given pixel $(x^{\prime} ,y^{\prime} )$ is calculated as the integral along the line of sight direction $\hat{z^{\prime} }:$

$$\begin{eqnarray*}&&I(x^{\prime} ,y^{\prime} )={I}_{0}+{\displaystyle \int }_{z^{\prime} =-{R}_{2}}^{{R}_{2}}\frac{{N}_{0}}{{r}^{2}}{\eta }_{\phi }(r,z)H(\theta ,g){dz}^{\prime} .\end{eqnarray*}$$

Here, ${\eta }_{\phi }(r,z),$ represents the dust density distribution, but tilted with a disk inclination, ϕ, relative to the observer's line of sight. The scattering angle $\theta =\theta (x^{\prime} ,y^{\prime} ,z^{\prime} )$ is a function of position.

The $1/{r}^{2}$ term accounts for the diminishing stellar flux as a function of distance from the star. The disk's PA, ${\theta }_{\mathrm{PA}},$ is implemented as a coordinate transformation between the stellocentric coordinates and the projected observer's coordinates. The constant, I₀, and the flux normalization, N₀, have been included to account for any possible biases and the conversion between model flux and detector counts, respectively. In summary, our model has a total of nine free parameters: ${R}_{1},{R}_{2},\beta ,{h}_{0},g,\phi ,{\theta }_{\mathrm{PA}},{N}_{0},{I}_{0}.$

Within the current model there exists a degeneracy between forward scattering ($g\gt 0$) with an inclination of $\phi \lt 90^\circ$ and backwards scattering ($g\lt 0$) with an inclination of $\phi \gt 90^\circ .$ In an effort to conserve computation time we chose to assume forward scattering and place a prior constraint on the scattering parameter, $g\gt 0,$ which is consistent with the model of Milli et al. (2014). A summary of the model parameters and their prior distributions can be found in Table 2.

Table 2.  Disk Model Parameters

Parameter	Symbol	Range	Prior Distribution
Inner Radius	R1	1–100 AU	Uniform in $\mathrm{log}{R}_{1}$
Outer Radius	R2	R1–500 AU	Uniform in $\mathrm{log}{R}_{2}$
Density Power Law Index	β	0.5–4	Uniform in β
Scale height aspect ratio	h0	0.01–2	Uniform in $\mathrm{log}{h}_{0}$
HG asymmetry parameter	g	0–1	Uniform in g
Line of Sight Inclination	ϕ	80°–90°	Uniform in $\mathrm{cos}\phi$
Position Angle	${\theta }_{\mathrm{PA}}$	25°–35°	Uniform in ${\theta }_{\mathrm{PA}}$
Flux Normalization	N0	1–1000	Uniform in $\mathrm{log}{N}_{0}$
Flux Offset	I0	−5–10	Uniform in I0

Download table as:  ASCIITypeset image

We fit the model to the GPI disk image using the parallel-tempering sampler from the emcee package. The diffraction limit in the H-band for Gemini south is 0043, equal to about three GPI pixels. We therefore apply a 3 × 3 pixel binning to both the Q_r and U_r images before fitting. This improves the noise statistics and speeds up the execution time of the MCMC fit, without significantly sacrificing spatial information. At each angular separation in the Q_r image, the errors were estimated as the standard deviation of a 3 pixel-wide annulus centered at that separation in the U_r image. The error estimates therefore contained photon noise, read noise and the unsubtracted instrumental polarization.

The MCMC sampler was run for 2500 steps with 2 temperatures, 128 walkers and burn-in of 500 steps. Additional temperature chains were not employed because of the additional computational cost incurred and the lack of evidence that the Markov chain sampler was only selecting local islands of high likelihood. One strength of using ensemble sampling over other types of sampling for MCMC fitting is that large speed-ups are possible via parallel-processing. On a 32-core (2.3 GHz) computer the entire MCMC run took nearly five days to complete.

After the run, the maximum auto-correlation across all parameters was found to be 85 steps, indicating that the chains should have reached equilibrium (Foreman-Mackey et al. 2013 recommend $\sim {\mathcal{O}}(10)$ autocorrelation times for convergence). In addition, the chains were examined by eye and appeared to have reached steady-state by the end of the burn-in phase. The posterior distributions (Figure 4) were estimated from the zero temperature walkers, using only one of out every 85 steps to ensure statistical independence of the surviving samples. The expected covariance between the inclination ϕ and g (Kalas & Jewitt 1995) is reproduced. Degeneracies are also found between ϕ and g, h₀ and g, h₀ and β, h₀ and ϕ, β and g, and β and ϕ, but in all cases, the parameters appear to be well constrained.

Zoom In Zoom Out Reset image size

The 16%, 50%, and 84% percentiles for each parameter are displayed in a table in the upper right corner of Figure 4. Marginalized across all parameters we find a disk inclined relative to the line of sight by $\phi =85\buildrel{\circ}\over{.} {27}_{-0.19}^{+0.26},$ with an inner radius of ${R}_{1}={23.6}_{-0.6}^{+0.9}\;\mathrm{AU}$, an outer radius of ${R}_{2}={139}_{-13}^{+19}\;\mathrm{AU}$ and an aspect ratio of ${h}_{0}={0.137}_{-0.006}^{+0.005}.$ The PA of the disk is ${\theta }_{\mathrm{PA}}=30\buildrel{\circ}\over{.} {35}_{-0.28}^{+0.29},$ where the errors include GPI's systematic error in PA (∼02). Note that the systematic uncertainty in the PA is not reflected in Figure 4. The dust is well fit by forward scattering grains, with a scattering asymmetry parameter of $g={0.736}_{-0.007}^{+0.008}.$ These results are further discussed in Section 5.1.

Figure 5 displays the best fit model and the residuals of the Q_r image minus the model. The best fit model was generated using the median value of each parameter in the marginalized posterior distribution. We find that the highest likelihood disk model successfully reproduces the GPI data. When examined at a different color scale, the residuals image displays similar low-level structure as that of the U_r image (Figure 2). The structure in the NW–SE direction is likely the Q_r counterpart of the residual instrumental polarization that's seen in the U_r image. A possible alternative explanation is that the structure could be due to a mismatch between the true scattering properties of the dust and the HG scattering function at small angular separations. A second structure can be seen along the disk midplane to the NE of the star. This asymmetric brightness feature is possibly due to a local overdensity of dust, that would increase the scattering at that location. Indeed, the β Pic disk is known to have multiple brightness asymmetries (Apai et al. 2015 provide a good summary). However, the feature is detected at similar brightness levels as the residual instrumental polarization and may yet be an uncharacterized artifact of the polarimetry reduction. Deeper observations of the disk will be required to distinguish between a true brightness asymmetry and instrumental effects.

Zoom In Zoom Out Reset image size

We describe here in broad terms our astrometric measurements and estimation of uncertainties, without delving into the details of each individual data set. For each epoch, the entirety of the data set is combined to estimate the planet's position relative to β Pic. The errors on this relative position are a combination of the error on the star's position, the planet's position, GPI's pixel scale and the accuracy to which we know GPI's orientation relative to true North.

For each data set, the stellar position was calculated using two methods. The wavelength slices of each datacube were first registered using the relative alignment procedure described in Section 2.2 and then collapsed into a broadband image. A Radon transform was then performed on the radially elongated satellite spots to find the stellar position (as in Pueyo et al. 2015). The stellar position was also estimated using the geometric mean of the satellite spot locations provided by the GPI DRP. Most H-band data sets show agreement between two methods at the 0.05 pixel level, with the exception of the 2013 December 11 commissioning sequence, during which extensive AO performance tests where being carried out. For K-band data sets the difference between the two methods is no more than 0.05 pixels and for the J-band data set it is 0.2 pixels. We found greater consistency in the relative location of β Pic b between observations obtained on consecutive nights when using the Radon method, and therefore chose to adopt the centroids measured with the Radon method for all measurements. For each data set, we considered the difference between the two methods as our estimate for the uncertainty on stellar position.

The location of β Pic b (in detector coordinates) was estimated at each wavelength channel, at each epoch and in each filter using the modified matched filter described in Pueyo et al. (2015). The uncertainty in β Pic b's location was estimated as the scatter in the position of the planet as a function of wavelength and number of principal components. We found the uncertainly to range from 0.05 pixels, for the data sets with significant field rotation and where the planet was at larger separations, up to 0.15 pixels, for the later epochs where the planet is significantly closer to the stellar host and SDI is less effective.

We estimated GPI's pixel scale using the methods described in Konopacky et al. (2014) by combining all the data presented therein with four new observations of ${\theta }^{1}$ Ori B, taken between 2014 September and 2015 January. We find an updated pixel scale value of 14.13 ± 0.01 mas/lenslet. Konopacky et al. (2014) measured a PA offset of −100 ± 003 during GPI commissioning. Subsequently, version 1.2 of the GPI DRP was updated to incorporate that 1° offset and correct for it automatically. Using the new measurements of ${\theta }^{1}$ Ori B, we find a residual PA offset of −011 ± 025.

Based on the measured location of β Pic b and its parent star in detector coordinates we calculated the relative separation and PA at each epoch. The separation was converted to milliarcseconds using the new platescale estimate and the PA was adjusted by $-0\buildrel{\circ}\over{.} 11.$ The separation and PA from each measurement can be found in Table 1. Uncertainties on these quantities were combined with the errors on the star position and planet position to yield the errors presented in the table.

β Pic b is detected in the Stokes I image as a point source superimposed on the extended PSF halo (Figure 3). After applying PSF subtraction using a python implementation of KLIP/ADI (Wang et al. 2015) to the image, the planet is recovered at extremely high S/N. The planet's position in the Stokes I image was estimated using the StarFinder IDL package (Diolaiti et al. 2000), which requires the user to input a PSF model for precision astrometry. In GPI's polarimetry mode the entire bandpass is seen by each frame and therefore the satellite spots are elongated and cannot be used as a PSF reference, as they are in spectroscopy mode. Instead, we used a GPI PSF generated with AO simulation software (Poyneer & Macintosh 2006).

To estimate astrometric errors we used StarFinder to measure β Pic b's location in the total intensity image from each of the 49 polarization data cubes. The rms difference between the planet location in the individual cubes and the Stokes I image was taken to be the error in the planet location. The error on the location of the star is estimated from the rms scatter of the measured star's position across the set of cubes. This error tracks the motion of the star behind the coronagraph between frames, which we expect to be larger than the errors on the star's position determined by the Radon transform, and therefore likely overestimates the errors.

The position of β Pic b in the polarimetry mode observations can be found in Table 1. As with the spectroscopy mode data, the errors represent a combination of the errors on the star's and planet's positions, GPI's pixel scale and GPI's PA offset on the sky.

Using the ten newly obtained astrometric points (nine from spec. mode and one from pol. mode), combined with the data sets presented by Chauvin et al. (2012) and Nielsen et al. (2014), we fit for the six Keplerian orbital elements of β Pic b plus the total mass of the system using the parallel-tempered sampler from emcee (Foreman-Mackey et al. 2014). While astrometric datapoints have been published in other papers, in an effort to minimize systematics between data sets, we limited ourselves to only these two large data sets where considerable effort has been made to calibrate astrometric errors. The fitting code was previously used in Kalas et al. (2013), Macintosh et al. (2014), and Pueyo et al. (2015). We also fit the radial velocity measurement of the planet from Snellen et al. (2014), which allows us to constrain the line of sight orbital direction and break the degeneracy between the locations of the ascending and descending node.

The model fits seven parameters: the semimajor axis, a; the epoch of periastron, τ; the argument of periastron, ω; the PA of the ascending node, Ω; the inclination, i; the eccentricity, e; and the total mass of the system, M_T. Our orbital frame of reference followed the binary star sign convention used in Green (1985). Under this convention the ascending node is defined as the location in the orbit where the planet crosses the plane of the sky (centered on the star), moving southward in projection. Note that this is 180° different from the convention used in Chauvin et al. (2012). The projected PA of the ascending node on the sky is defined from north, increasing to the east. The argument of the periastron is defined as the angle between the ascending node and the location of the periastron in the orbit, with ω increasing from the ascending node. The epoch of periastron, τ, is defined in units of orbital period, from 1995 October 10 (Julian date 2450000.5). A summary of the orbital parameters and their prior distributions can be found in Table 3.

Table 3.  Orbit Model Parameters

Parameter	Symbol	Range	Prior Distribution
Semimajor axis	a	4–40 AU	Uniform in $\mathrm{log}a$
Epoch of Periastron	τ	−1.0–1.0	Uniform in τ
Argument of Periastron	ω	$-2\pi -2\pi$ rad	Uniform in ω
Position Angle of the Ascending Node	Ω	25°–85°	Uniform in Ω
Inclination	i	81°–99°	Uniform in $\mathrm{cos}i$
Eccentricity	e	0.00001–0.99	Uniform in e
Total Mass	MT	0–3 M⊙	Uniform in ${M}_{\odot }$

Download table as:  ASCIITypeset image

The MCMC sampler was run for 10,000 steps with 10 temperatures and 256 walkers after a "burn-in" of 2000 steps. After the run, the maximum auto-correlation across all parameters was found to be 25 steps, indicating that the chains should have reached equilibrium. The posterior distributions (Figure 6) were estimated using the zero temperature walkers, using only one of out every 25 steps. In Figure 6, the epoch of periastron was wrapped around to be only positive between 0 and 1 and the arguments of the periastron was wrapped around to range from 0° to 360°. A random selection of 500 accepted orbits are plotted on top of the astrometric and radial velocity datapoints in Figures 7 and 8, respectively. While the orbital fits are generally consistent with most of the astrometric datapoints, the majority of the orbital solutions fall more than 1σ from the measured radial velocity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the planet has a semimajor axis of ${9.25}_{-0.45}^{+1.46}\;\mathrm{AU}$, an inclination of 8901 ± 030, and an ascending node at a PA of 3175 ± 015. We take the median of the marginalized posterior distribution to be the best estimator of each parameter's value, and the 68% confidence values as the errors. Following this convention, the eccentricity of the orbit is found to be $e={0.07}_{-0.05}^{+0.11}.$ However, the eccentricity is a positive definite quantity and typical estimators (e.g., the mean and median) will overestimate the true eccentricity when it is small ($e\lt 0.1$). When considering eccentricities of radial velocity planets, Zakamska et al. (2011) consider several different estimators and find that for small eccentricities the mode of the distribution is the least biased. The mode of our distribution falls in the smallest eccentricity bin indicating an eccentricity very close to zero. Therefore, it is perhaps more appropriate to quote the upper limit on the eccentricity $e\lt 0.26$ (95% confidence).

For orbits with higher eccentricities ($e\gt 0.1$), the epoch and argument of periastron have strong peaks at ∼0.5 periods and $\sim 170^\circ ,$ respectively. At lower eccentricities these two parameters remain degenerate, with a large range of acceptable values. Overall, the marginalized distributions reveal that these parameters are still relatively unconstrained. The ensemble of accepted orbits at the end of the run have a reduced ${\chi }^{2}$ of ${1.55}_{-0.05}^{+0.09}.$

Marginalized across all parameters, the total mass of the system is 1.61 ± 0.05 M_⊙. At $\sim 11\;{M}_{\mathrm{Jup}},$ β Pic b contributes less than 1% to the total mass, giving β Pic itself a mass of ${M}_{\star }=1.60\pm 0.05\;{M}_{\odot }.$ This falls slightly below the range estimated by Crifo et al. (1997) (1.7–1.8M_⊙) and just within the range of Blondel & Djie (2006) (1.65–1.87 M_⊙), who both use evolutionary models and the HR diagram to date β Pic. This estimate provides a slightly smaller value than that presented in Nielsen et al. (2014), though still consistent within their errors (${1.76}_{-0.17}^{+0.18}\;{M}_{\odot }$)

By combining the semimajor axis and stellar mass values of each walker at each accepted iteration, we are able to create a posterior distribution for the orbital period, from which we derive that $P={22.4}_{-1.5}^{+5.3}\;\mathrm{years}$. The large upper limit is due to the extended tail in the semimajor axis distribution.

Giant exoplanets may have polarized emission in the NIR either due to rotationally induced oblateness (Marley & Sengupta 2011) or asymmetries in cloud cover (de Kok et al. 2011). For β Pic b, the recently measured rotational period of ∼8 hr would induce a polarization signature due to oblateness of less than 0.1% (below GPI's current sensitivity limit, Wiktorowicz et al. 2014). Therefore, any detected polarization signal would be indicative of cloudy structure.

To estimate β Pic b's polarization, we first created a disk-free linear polarized intensity image by combining the model-subtracted Q_r image with the U_r image ($P=\sqrt{{Q}_{r}^{2}+{U}_{r}^{2}}$). The total polarized flux at the location of β Pic b, within an aperture of radius $1.22\lambda /D,$ was then compared to the flux of 38 independent apertures at the same angular separation. We find that β Pic b's polarized flux is $0.5\sigma$ from the mean flux of the independent apertures, consistent with zero linear polarization signal from the planet (see Figure 3). While this measurement does not provide any evidence for cloud structure, it does not exclude the possibility either; the magnitude of cloud-induced polarization depends on many factors, including the atmospheric temperature and pressure profile, the composition, the nature of the inhomogeneities, rotation, and viewing angle. The PSF variability during the observations makes accurate recovery of the total intensity of the planet difficult, and thus we leave the characterization of an upper limit on the planet's polarization fraction for future work.

With GPI we probe the projected disk between 03 and 17 at high spatial resolution. The work presented here has two advantages over previous attempts to model the disk at similar angular separations. First, the polarized intensity images provide a unique view of the disk, in particular the vertical extent is free of any biases associated with ADI PSF subtraction. Second, the MCMC fitting allows us to fully explore the multi-dimensional parameter-space and place quantitative confidence intervals on the model parameters.

MCMC fitting requires evaluation of the likelihood function for each set of parameters that is examined. The cost of fitting depends on the computational expense of evaluating the model and the dimensionality of the model parameter space. For that reason we have limited our exploration to optically thin scattering, an analytic recipe for the phase function, and a simple model of the dust distribution. We do not consider multi-component disks (as modeled for the outer disk, e.g., Ahmic et al. 2009) and we assume that the disk aspect ratio is constant. Regardless of these simplifications, we find that this model provides an excellent fit to our polarized image.

The HG scattering function is often used to model the total intensity scattering efficiency of dust grains, but has not been used extensively for polarized intensity. This is at least partially due to the fact that in most circumstances where polarized intensity is measured, total intensity is obtained as well, allowing for more sophisticated modeling of the dust scattering. In addition, the scattering efficiency of polarized intensity of small spherical particles approaches zero at very small scattering angles, a feature that is not captured by the HG function. While the exact shape of the HG function cannot fully reproduce the polarized scattering efficiency function for physical models, a quick informal survey of possible grain models indicates that our best fit $g=\langle \mathrm{cos}\theta \rangle \approx 0.7$ can be reproduced by Mie scattering particles with a radius of ∼1 μm and an index of refraction of m = 1.033–0.01i, similar to the porous, icy grains inferred by Graham et al. (2007) for AU Mic. However, as previously mentioned, a true characterization of the physical grain scattering properties will require either an unbiased total intensity image, or polarized intensity images at other wavelengths. We leave the characterization of the dust properties of the inner disk for future work.

The observations of Milli et al. (2014) have a field of view (04–38) that overlaps with our disk detection and therefore provide a good point of comparison. They model the L' emission with a single component disk model similar to ours, albeit with different radial and vertical dust density profiles. Even so, their best fit inclination (i = 86°) and PA (${\theta }_{\mathrm{PA}}=30\buildrel{\circ}\over{.} 8$) agree fairly well with our own. Their data set constrains the sky-plane inclination less precisely and inclinations of 85°–88° provide good fits to their data. The consistency between their measurements and ours builds confidence that the measured angles are not highly sensitive to the assumed scattering properties and radial dust distribution. The PA of the disk seen in our images (${\theta }_{\mathrm{PA}}=30\buildrel{\circ}\over{.} {35}_{-0.28}^{+0.29}$) and those of Milli et al. (2014) appears to be misaligned from both the outer main disk (${\theta }_{\mathrm{PA}}=29\buildrel{\circ}\over{.} 1\pm 1^\circ ;$ Apai et al. 2015) and the warp (${\theta }_{\mathrm{PA}}\approx 32^\circ -33^\circ$). This offset, and how our disk images fit into the context of the whole system, will be further discussed in Section 5.3.

The results of our model fitting reveal an inclined disk with an inner radius of ${R}_{1}\approx 23.5\;\mathrm{AU}$ (12), populated by grains that preferentially forward scatter polarized light. The majority of the detected polarized flux is therefore inside the projected inner radius and the result of forward scattering by the constituent dust grains. Without direct detections of either the inner or outer radius, the constraints on both are governed by the overall shape and spacing of isophotal contours (see Figure 9).

Zoom In Zoom Out Reset image size

The location of the inner edge of the disk seen in our model is a unique feature of this work and has not been found in previous scattered light imaging at similar angular separations. This could be attributed to both the scattering properties of the dust, which make the inner edge difficult to see, and the modeling strategies used in those studies. Milli et al. (2014) also use a HG function to model their dust. However, their model considers a population of parent bodies between 50 AU and 120 AU, with the density falling as separate power laws inside and outside of these radii and they do not define an inner radius in the same manner as in our model. Apai et al. (2015) make surface brightness measurements of the disk between 05 and 150, but find no noticeable change in the brightness profile at 11. In our model, we find that the forward scattering nature of the dust grains means that the inner edge itself contributes minimally to the observed surface brightness at its projected separation. This serves to emphasize the critical role of dust scattering when interpreting the surface brightness as a function of radius; a smooth surface density by itself does not necessarily exclude features in the radial dust profile.

Note that our model has been defined with a sharp cut-off inside the inner radius, and caution should be used when interpreting the exact value. There may be dust inside of the inner radius with a lower surface density. For example, the true dust density inside the inner radius may have a slowly decreasing inner power-law, such as those considered in Milli et al. (2014).

Imaging and spectroscopic studies in the mid-IR have probed similar regions of the debris disk at wavelengths where contrast between the stellar flux and the dust (thermal) emission is more favorable than in the optical and NIR. Okamoto et al. (2004) found spectroscopic evidence for dust belts at 6, 16, and 30 AU. Wahhaj et al. (2003) fit a series of four dust belts to deconvolved 18 μm images and found their best fit radii to be 14, 28, 52, and 82 AU. With the exception of the belts close to ∼15 AU, all of these belts are either well outside or at the very edge of our field of view. The 6 AU belt seen by Okamoto et al. (2004) is below our inner working angle. We note that the Okamoto et al. (2004) belt at ∼16 AU only occurs on the NE side, at roughly the same location as the tentative brightness asymmetry seen in our disk model residuals. We see no evidence of the other belts in projection, but we model the disk with a continuous dust distribution and therefore may not be sensitive to dust at their locations. Mid-IR imaging by Weinberger et al. (2003) indicates emission within 20 AU that is significantly offset in PA from the main outer disk. In our disk image we see no indication of this offset.

Previous studies of β Pic's debris disk in polarized scattered light have been carried out both in the optical (Gledhill et al. 1991; Wolstencroft et al. 1995) and the NIR (Tamura et al. 2006). These observations image the disk at separations of 15''–30'' and 26–64, in the optical and NIR, respectively. At these angular separations the total intensity observations are not limited by the PSF halo and when combined with the polarized images, polarization fraction can be used to model the dust grains (Voshchinnikov & Krügel 1999; Krivova et al. 2000). Tamura et al. (2006) combine the optical measurements with their K-band data and find that the observations could be explained by scattering from fluffy aggregates made up of sub-micron dust grains. Unfortunately, the lack of total intensity images and a non-overlapping field of view make a direct comparison between our observations and this past work difficult.

In general, our orbit fit is consistent with those previously published (e.g., Chauvin et al. 2012; Macintosh et al. 2014; Nielsen et al. 2014), but the longer temporal baseline and increased astrometric precision significantly tighten the constraints on the orbital parameters. In particular, we find that the PA of the ascending node of the planet lies in between the main outer disk and warp feature, consistent with Nielsen et al. (2014).

At first glance, the errors on our orbital elements appear comparable to those in Macintosh et al. (2014). However, our fit includes the total mass of the system as an additional free parameter. Nielsen et al. (2014) modeled the system's total mass as a free parameter in their orbital fit and found that accounting for the uncertainty in the system's total mass resulted in larger uncertainties in the planet's orbital elements. In particular, they find that with a floating system mass the eccentricity distribution has a long tail that peaks at high eccentricities. Due to a degeneracy between semimajor axis and eccentricity, this stretched the semimajor axis distribution to higher values as well. In Figure 6, we find that the eccentricity is now significantly better constrained ($e\lt 0.26$), and while the degeneracy remains, the semimajor axis is constrained to be $\lt 10.7\;\mathrm{AU}$ with 84% confidence.

For each orbit defining our posterior distribution, we calculate the epoch of closest approach and find that it will fall between 2017 November 20, and 2018 April 4 with 68% confidence. With our derived inclination of i = 8901 ± 036, the updated transit probability is ∼0.06%, assuming that the planet will transit if the inclination is within 005 from 90°. This is a reduction by a factor of ∼50 from the estimate in Macintosh et al. (2014), who found i = 90.7 ± 0.7.

Even though the likelihood of a planet transit is small, it is still possible that dust particles orbiting within the planet's Hill sphere (${R}_{\mathrm{Hill}}\approx 1$ AU) will transit. Indeed the transit of a ring system surrounding an exoplanet was recently detected around J1407 (Mamajek et al. 2012). In the outer solar system, satellites around the giant planets have stable orbits within a Hill sphere about the planet out to $\sim 0.5{R}_{\mathrm{Hill}}$ when in prograde orbits and $\sim 0.7{R}_{\mathrm{Hill}}$ in retrograde orbits (Shen & Tremaine 2008). For β Pic b, assuming a planetary mass of 11 M_J, a semimajor axis of 9.25 AU, a circular orbit and a stellar mass of 1.60 M_⊙, we calculate a Hill radius of ∼1.2 AU. Thus, stable orbits within the Hill sphere will transit if the planet's inclination is within 38 and 53 of edge-on, for prograde and retrograde orbits, respectively. Our new constraints on the inclination indicate that these orbits will almost certainly transit. However, the true transit probability will depend not only on the exact location of the dust, but also its orientation relative to the observer. For example, dust that fills the stable orbits and is orbiting face-on relative to the observer will transit, but if it is orbiting edge-on it will not.

The presence of infalling comets (a.k.a. falling evaporating bodies, or FEBs) has been previously inferred by redshifted absorption features in β Pic's spectrum (Beust & Morbidelli 1996). Thébault & Beust (2001) suggested that a massive ($M\geqslant {M}_{\mathrm{Jup}}$) planet within ∼20 AU on a slightly eccentric orbit ($e\gtrsim 0.05-1$), could be responsible for imparting highly elliptical orbits on bodies within a 3:1 or 4:1 resonance, that then plunge toward the star. In this scenario the argument of the periastron of the planet is restricted to a value of ${\omega }^{\prime }=-70^\circ \pm 20^\circ$ from the line of sight. Using our definitions, the equivalent requirement is $\omega =200^\circ \pm 20^\circ .$ Our results neither confirm nor rule out the infalling comet scenario. While the marginalized distribution of the argument of periastron allows for a broad range of values, if the orbit is indeed eccentric, then ω peaks strongly around 170°, just outside of the acceptable values for this scenario. Thébault & Beust (2001) find that if the eccentricity of the massive perturber (here assumed to be β Pic b) is as large as $e\approx 0.1$ then the infalling comets most likely originate in the 3:1 resonance, which occurs between 18 and 22 AU based on our 68% confidence range for β Pic b. The inner edge of the dust in our scattered light images falls at ${R}_{1}={23.61}_{-0.57}^{+0.86}\;\mathrm{AU}$, outside of range of values for the 3:1 resonance. However, as noted above, our inner radius is sharply defined, and there may still be material inside. For smaller perturber eccentricities the infalling comets originate in the 4:1 resonance, which occurs between 22.2 and 26.5 AU. Our disk model does not constrain whether there is an excess of bodies librating in the 4:1 resonance.

A planet on an inclined orbit is thought to be responsible for the warp feature in the region of the disk outside our field of view at ∼80 AU (Mouillet et al. 1997; Augereau et al. 2001). The directly imaged planet β Pic b (Lagrange et al. 2009) has a mass, semimajor axis, and inclination consistent with producing the warp (e.g., Dawson et al. 2011). The updated PA of β Pic b's ascending node, ${\rm{\Omega }}=31\buildrel{\circ}\over{.} 75\pm 0\buildrel{\circ}\over{.} 15,$ is offset by 265 with respect to the flat outer disk (291 ± 01; Apai et al. 2015), consistent with producing a warp tilted by 5° counter-clockwise with respect to the flat outer disk. As illustrated by Apai et al. (2015, Figure 21), our azimuthal viewing angle of the warped disk affects the degree to which the inner disk and the planet's orbit appear aligned with the flat outer disk in projection and also affects the projected height of the warp. Although the planet's updated orbit remains consistent with sculpting the outer regions of disk, several features of the inner regions of the disk that we measured here are unexpected solely from sculpting by β Pic b (Figure 10).

Zoom In Zoom Out Reset image size

First, the inner edge of the disk is at 23.6 AU, about twelve Hill radii from the planet. We performed a simulation using mercury6 (Chambers 1999) of a planet with orbital parameters set to the median values in Figure 3 embedded in a disk of test particles initially spanning 10–40 AU. On the timescale of hundreds to thousands of orbits, the planet clears out the disk to ∼15 AU, with the inner edge persisting at that location over the 20 Myr stellar lifetime (Binks & Jeffries 2014; Mamajek & Bell 2014). An inner edge at ∼15 AU is in agreement with simulations by Rodigas et al. (2014; see their Table 2).

We have not explored disk models with gradual inner edges (e.g., Milli et al. 2014), so there may be material between 15 and 23.6 AU with a lower surface density, or planet bodies that are less collisionally active. If the disk inside 23.6 AU is truly cleared out, an undetected low-mass planet in between β Pic b and the disk's inner edge could be responsible; we find that a planet could exist on a stable orbit in that region.

Second, we expect the inner disk to be centered on the planet's orbital plane. Given a warp located at ∼80 AU, the width of a secular cycle (i.e., the difference in semimajor axis for which the planetesimals are $2\pi$ out of phase in their oscillations about the planet's orbital plane) is only about 1 AU at a radius of 40 AU and the timescale of a secular cycle is about forty times shorter than at the location of the warp. Therefore, close to the planet, a sufficient number of secular cycles should have passed that the parent bodies' free inclination vectors are randomized about the forced inclination from the planet. Under certain conditions, we found that our simulation could produce a parent bodies sky plane inclinations distribution with peaks at $\sim {i}_{\mathrm{planet},\ \mathrm{sky}}\pm {i}_{\mathrm{planet},\ \mathrm{outer}\ \mathrm{disk}}$ (one of which could correspond to $\sim 85^\circ$), where ${i}_{\mathrm{planet},\ \mathrm{sky}}$ is the line of sight inclination of the planet and ${i}_{\mathrm{planet},\ \mathrm{outer}\ \mathrm{disk}}$ is the mutual inclination between the planet and the outer disk. However, we expect that even in these circumstances the measured disk midplane would be aligned with the planet's orbital plane. Moreover, damping by collisions, small bodies, or residual gas—provided that it occurs on a timescale shorter than half a secular timescale—reduces the free inclination, decreasing the disk scale height but keeping it centered about the planet's orbital plane.

Instead, the average plane of the inner disk appears mutually inclined with respect to the planet's orbit. If the polarized intensity images were dominated by scattered light from the outer disk, a mutual inclination with respect to the planet could be consistent, (depending on the semimajor axis of the dominant dust; see Figure 1 from Dawson et al. 2011), but in the current disk model the observed light is dominated by a close-in disk. Contribution from another planet to the forced plane of the disk is a possibility but the available parameter space for an additional planet that tilts the disk toward us, yet is too low mass to escape detection, is quite limited. In the future, we plan to explore a wider range of dust-scattering models to ensure that this result (a disk mutually inclined to the planet's orbit at ∼25 AU) is not dependent on the assumed dust properties.

Finally, the scale height of the disk appears larger than expected from stirring by β Pic b or self-stirring of the parent bodies. In the absence of damping, the total thickness of the disk would be $\sim 2{i}_{p},$ corresponding to a scale height aspect ratio of about 0.06 for a planet inclined by 36 with respect to the primordial plane. Self-stirring to the escape velocity of 10 km planetesimals would contribute only about 0.001 to the aspect ratio; self-stirring to the escape velocity of Pluto sized bodies would be required. In practice, we do not expect most parent bodies participating in the collisional cascade to be stirred to random velocities of the largest bodies (e.g., Pan & Schlichting 2012); their steady-state random velocities depend on the balance between stirring, damping by smaller bodies and each other, and radiation forces. The scale height is also significantly larger than observed further out in the disk (Ahmic et al. 2009)—even at 50 AU (Milli et al. 2014). The robustness of the scale height to the dust scattering model should be explored further; for example, a significant contribution from polarized back scattering (not modeled here) could result in a smaller inferred scale height. If the current inferred large scale height in the very inner disk is robust, a sub-detection planet located between β Pic b and the inner edge of the disk and mutually inclined with respect to β Pic b is a possible explanation.

Nesvold & Kuchner (2015) recently simulated the dynamical and collisional behavior of β Pic's planetesimals and dust grains using SMACK (Nesvold et al. 2013), which models planetesimals across a range of sizes using super particles. They find that collisional damping is not important in shaping the morphology of the disk. Their detailed model also does not predict the surprising observational features discovered here: they find the planet clears a gap only out to 14.5 AU and that the disk is centered about the planet's orbital plane (see their Figure 3). They find that some planetesimals in the inner disk are scattered by each other or the planet to inclinations larger than $2{i}_{p},$ increasing the thickness of the inner disk by about 50%, not enough to account for the ($\sim 200\%$ larger) observed scale height.