Discovery of an Edge-on Circumstellar Debris Disk around BD+45° 598: A Newly Identified Member of the β Pictoris Moving Group

Table 1 provides the detailed parameters of our observations of BD+45° 598, which was first observed on UT 2013 September 25 with the NIRC2 instrument at the W.M. Keck Observatory using the $K^{\prime}$ filter (2.124 μm central wavelength) and the 400 milliarcsecond coronagraphic mask. The two subsequent epochs of observation of BD+45° 598 were obtained using this mask, and the 200 milliarcsecond radius of this mask effectively sets the ultimate inner working angle for our observations. Seeing measurements obtained at Keck indicate that the average 0.5 μm seeing for the 2013 and 2014 epochs ranged from 05 to 06, but was poorer (∼08) for the 2015 epoch. All scattered light observations were obtained using the Keck Adaptive Optics ("AO") system (Wizinowich et al. 2000), operating in pupil tracking mode (Marois et al. 2006), to differentiate quasi-static speckles (e.g., Hinkley et al. 2007) from any bona fide astrophysical sources. To remove the residual, uncorrected starlight, we model this speckle pattern using a customized algorithm based on projections of the data onto Karhunen–Loève eigenimages (the Karhunen–Loève Image Processing or "KLIP" algorithm, Soummer et al. 2012). Subtracting this model from each image provides deep contrast relative to the host star. However, because of the ∼46° decl. of BD+45° 598, Table 1 shows that typically only 18°–36° of parallactic angle rotation could be achieved for each of our three epochs of scattered light imaging, meaning a substantial amount of residual scattered starlight contamination remains within ∼05.

In Figure 2 we show an image of BD+45° 598 with our three epochs of data coadded together, as well as a corresponding signal-to-noise (S/N) map of the same image. The disk was immediately identified in our first epoch (Figure 2, top panels), extending to ∼1'' on each side of the star. The disk structure in the images was robust against the number of eigenmodes chosen in the post-processing algorithm, further indicating that the structure is a true astrophysical signal, and not residual, scattered starlight. However, to verify that the signal was indeed a physical structure associated with the star, BD+45° 598 was again detected using the $K^{\prime}$ filter at two subsequent epochs, 2014 October 2 and 2015 January 10 (see Figure 2). We also attempted another observation of BD+45° 598 on 2014 November 11, but this data set suffered from poor atmospheric conditions corresponding to poor Adaptive Optics correction, resulting in a nondetection of the disk.

Zoom In Zoom Out Reset image size

BD+45° 598 was not observed by either the Spitzer or Herschel observatories, and there are only upper limits available at 60 μm and 100 μm from IRAS. So, as an additional check that this structure seen in the Keck images was not an artifact from the image post-processing, BD+45° 598 was observed at various dates between 2016 February 19 and 2016 August 10 for a total of 4.5 hr at 850 μm (Band 3) using the SCUBA-2 instrument at the James Clerk Maxwell Telescope. The SCUBA-2 observations achieved an rms sensitivity of 1.1 mJy per 14.5 arcsecond beam, but no signal at the location of BD+45° 598 was detected. These observations thus provide an upper limit of 3.3 mJy at 850 μm, and we do not include these limits for our analysis. Using the Submillimeter Array (SMA) on 2015 August 10, we observed BD+45° 598 at 224 GHz (1.3 mm) in a six hour "filler track." The array configuration using six antennas resulted in a 5.0 × 27 beam, and using a point-source fit to the imaged visibilities, our observations revealed a ∼3.3σ detection of a source with a flux density of 2.17 ± 0.76 mJy, coincident with the location of the disk seen in scattered light. The observations followed standard procedures, with 3c84 and 0303+472 used as gain calibrators, and Uranus as the flux calibrator.

With the goal of obtaining better resolution on the unresolved SMA detection, we observed the field around BD+45° 598 at 230.5 GHz during the winter periods in 2015–16 and 2016–17 using seven antennae with the IRAM NOEMA interferometer at Plateau de Bure at an altitude of 2552 meters. The first session used the compact BD array configuration, with projected baseline lengths <200 m on 2015 November 27 (2.1 hr) and 2015 December 02 (1.1 hr), the second session used the extended A array configuration, with projected baseline lengths up to 800 m on 2017 January 3 (6.4 hr). All observations were conducted during fair weather conditions (precipitable water vapor ∼2–3 mm) apart from those on 2015 December 2, which were not used in our data processing. The data were processed with the Widex correlator over a bandwidth of 3.6 GHz centered on 230.5 GHz in both polarizations combined to provide total intensity visibilities. Phase calibrators J0136+478 and J0300+470 were used to monitor the atmospheric phase variations once every 30 minutes to remove its effect on the target phase. Absolute flux density calibration relies on the monitoring program of standard calibrators (MWC 349 and LkHa 101) conducted routinely by IRAM. At 1.3 mm, the flux density uncertainty can be up to 20%.

The NOEMA data were processed with the IRAM software GILDAS. The data acquired at the three epochs of observations were combined to provide satisfactory uv-coverage. Measured complex visibilities were Fourier inverted and deconvolved using the CLEAN algorithm (Högbom 1974) to provide the final intensity map shown in Figure 3. The data processing used natural cleaning with an elliptical support large enough to include the source and keep some sidelobes (5'' by 10'' with a position angle close to the disk orientation in the Keck image). One thousand iterations were carried out, but the convergence of the flux density was satisfactorily reached after a few hundred. The image has a beam size of (068 × 052) and is limited in extent by the primary beam of individual antennas (18'' at 230 GHz).

Zoom In Zoom Out Reset image size

In the CLEAN map shown in Figure 3, we find a millimeter-wavelength structure detected with a statistical significance of 4σ–5σ, elongated and aligned with the edge-on disk seen in the Keck image (see discussion of the NOEMA data modeling in Section 3.2). The NOEMA structure is well centered at the optical position of BD+45° 598 provided by the Gaia DR2 catalog (Gaia Collaboration et al. 2018) updated with its proper motion at the observation midpoint (α = 02^h21^m1315055 ± 000013 and $\delta =+46^\circ 00^{\prime} 06\buildrel{\prime\prime}\over{.} 4333\pm 0\buildrel{\prime\prime}\over{.} 0011$ at 2016.0). The Gaia and NOEMA celestial reference frames are aligned at the submilliarcsecond level (Lindegren 2020).

The elongated NOEMA structure in the map of Figure 3 is almost certainly millimeter-sized dust that traces a belt of planetesimals centered on BD+45° 598 corresponding to the Keck scattered light images. We note that the west side of the structure is more prominent than the east side, but our modeling finds that any brightness asymmetry is only at the ∼2σ level. The integrated flux density of the disk at λ = 1.3 mm is 1.11 ± 0.17 mJy (statistical uncertainty). This measurement was done by fitting an ellipsoid to the data in the uv-plane with a priori values conjectured from the image itself. The photospheric flux density of BD+45° 598 is three orders of magnitude smaller at this wavelength and so does not contribute. Finally, using the beam full width at half maximum combined with the distance to BD+45° 598, we can infer that the belt width is smaller than 35 au.

For the purpose of kinematic analysis, we obtained an improved radial velocity of the star by observing with the IGRINS spectrograph (Park et al. 2014) at the McDonald Observatory on 2015 October 23. A total exposure time of 360 s was obtained in the "ABBA" observing pattern, and 200 s were obtained for the purpose of telluric calibration on the star HD14212. Cross-correlation against other objects with similar spectral types provides an absolute radial velocity of 0.18 ± 0.22 km s⁻¹.

Our $K^{\prime}$-band coronagraphic images are presented in Figure 2 for the three epochs where we model the disk geometry. No disk signal was detected during our 2014 November 11 observations, and so we do not model these data. For the three epochs shown in Figure 2, we detect a peak S/N for the disk of 4.8, 4.8, and 4.7, respectively. The bottom panel of Figure 1 also shows an S/N map for the combined data set, indicating a peak S/N of 5.8. Figure 1 (top panel) also shows a measure of the disk surface brightness, indicating a peak surface brightness of ∼1–2 mJy arcsec⁻². However, it should be noted that in the speckle dominated regime at inner working angles ≲05, some "self-subtraction" due to the KLIP algorithm will diminish the surface brightness of the disk by a factor of ∼2.

To constrain the disk geometry, we use an injection modeling process as described in Wahhaj et al. (2014) and Matthews et al. (2017). We fit a circular disk with five free parameters to the data. These parameters are the position angle, inclination angle, disk radius, forward scattering (parameterized by the Henyey–Greenstein parameter g), and an overall scaling term. Synthetic model images of an optically thin ring are generated using the GRaTeR radiative transfer code (Augereau et al. 1999), and then convolved with a Gaussian function to reproduce the resolution of the telescope. To generate disk models with GRaTeR, we use the following parameterization. The disk density is calculated 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 the remaining terms are parameters defining the disk vertical and radial structure. For this work, the only disk parameters we fit are ρ₀ and r₀; values for the other parameters are fixed at α_in = 10, α_out = − 3, ξ₀ = 1, β = 0, and γ = 2, i.e., we fit a thin, narrow ring and do not attempt to constrain the vertical disk structure since this is unresolved in the Keck/NIRC2 data.

The scattered light contribution from each point in the disk is calculated as:

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

where θ is the scattering angle of the dust, and d is the distance from the star to each point, and p(θ) is the scattering phase function. In each fit, the disk is scaled by an arbitrary value ρ₀, which takes into account the effect of various factors including the stellar luminosity, stellar distance, and telescope gain. For p(θ) we use the Henyey–Greenstein scattering function:

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

This calculation is performed over a grid encompassing the entire disk, and contributions along each line of sight are added to create a final disk image. As well as the ρ₀, r₀, and g disk parameters, we fit two geometrical viewing terms, namely the position angle and the inclination of the disk, denoted by PA and i_tilt.

These models are then rotated and subtracted from each individual data frame in the cleaned data cubes. These data cubes with negative disks injected are processed using the same KLIP algorithm discussed in Section 2.1. By injecting negative disk images into the cleaned data, we mimic a forward-modeling procedure whereby the throughput of the KLIP algorithm is accurately taken into account in the disk modeling. To determine the goodness-of-fit of each disk model, we calculate the χ² of a rectangular region enclosing the disk, with central pixels close to the star excluded.

The uncertainty per pixel is calculated as a function of radius. We first take the standard deviation of all those pixels outside the mask, in small annuli centered on the star. We then inflate these errors to account for the correlation between pixels. The inflation term is the square root of the number of pixels per correlated region, and is calculated from both the FWHM of the instrument and the radial elongation of speckles due to the filter width (see Figure 1). At radius r, the inflated error is $\sigma (r)=\sqrt{\mathrm{FWHM}\times {\rm{r}}{\rm{\Delta }}\lambda /\lambda }\times {\sigma }_{\mathrm{pixel}}({\rm{r}})$, where λ and Δλ are the central wavelength and bandpass of the NIRC2 $K^{\prime}$ filter. Using this method we find an initial fit value for the first epoch of data using a downhill minimization algorithm. We then initiate Metropolis Hastings MCMC chains from this initial value to explore the parameter space, and calculate the best fit and uncertainty values. MCMC chains are generated with emcee (Foreman-Mackey et al. 2013), probabilities are assigned to each sample as $\exp \left(-{\chi }^{2}/2\right)$, and we use uniform priors.

We find good fits to each of the three epochs, with the parameters showing good agreement between epochs. Images of the three best-fit models are shown in Figure 2 and best-fit parameters are listed in Table 2, where we list the median and uncertainties, which include 34% of the samples on each side of this median value. The disk is ∼3° from edge on, and has a radius of ∼70 au, although the radius is poorly constrained due to the viewing geometry. Given the relatively low signal-to-noise ratio of the disk, we are not able to detect additional structure with the current data.

Table 2. Best-fit Model Parameters to Our $K^{\prime}$ Scattered Light Data

	2013 Sep 25	2014 Oct 2	2015 Jan 10
PA [°]	${109.9}_{-0.7}^{+0.6}$	${109.5}_{-1.0}^{+0.9}$	110.2 ± 1.2
itilt [°]	${86.9}_{-1.6}^{+3.0}$	${87.1}_{-1.7}^{+3.2}$	${86.1}_{-2.0}^{+2.2}$
ρ0	${1.0}_{-0.4}^{+0.9}$	${0.7}_{-0.3}^{+1.0}$	${0.9}_{-0.4}^{+0.7}$
r0 [au]	${71}_{-17}^{+19}$	${70}_{-27}^{+16}$	${68}_{-11}^{+20}$
g	${0.44}_{-0.23}^{+0.32}$	${0.28}_{-0.17}^{+0.43}$	${0.33}_{-0.16}^{+0.36}$

Note. ρ₀ values are quoted relative to the first epoch.

Download table as:  ASCIITypeset image

As the Keck scattered light data and the IRAM NOEMA 1.3 mm data are sensitive to dust populations at different radii, we have chosen to model the Keck and NOEMA data separately, and not to perform a simultaneous fit. We have modeled the NOEMA image with the thermal emission of an optically thin ring and again use a Markov Chain Monte Carlo (MCMC) approach using emcee (Foreman-Mackey et al. 2013). The fitting is done by generating disk images, which are then convolved with the NOEMA beam, and then subtracted from the data to create a residual map (from which we compute a χ² value). With this approach, the NOEMA radius is constrained to $0\buildrel{\prime\prime}\over{.} {85}_{-0.25}^{+0.18}$ (${62}_{-18}^{+13}$ au) comparable to the Kuiper Belt radius, the inclination is ${87}_{-0.45}^{\circ +0.46}$ and the PA is ${109.85}_{\,-0.09}^{\circ +0.11}$. The belt radius derived from the NOEMA data is smaller than that derived from the Keck scattered light imaging (r₀ ∼ 70 au in Table 2) as expected since scattered light imaging will be sensitive to micron-sized grains that are likely blown to high eccentricity orbits by radiation pressure, while the NOEMA data set will reveal millimeter-sized grains tracing the belt of larger planetesimals. Higher dynamic range maps would be needed for a detailed comparison of the NOEMA radius with the radius parameter r₀ of the grain surface density in the GRaTeR model used to fit the scattered light image. Lastly, with the current data, we cannot place strong limits on the belt width besides our constraint that it must be less that 35 au based on FWHM argument mentioned in Section 2.3.

To construct the SED shown in Figure 4, we collected optical (Tycho-2, APASS, Gaia) and infrared (2MASS, AKARI, WISE, IRAS) photometry for BD+45° 598. We fit the SED with various models, using PHOENIX spectra for the stellar photosphere (Hauschildt et al. 1997), and amorphous silicate spectra for the dust opacity (Li & Greenberg 1997). The silicate models use three parameters; a size distribution n(D) ∝ D^2−3q with a minimum size ${D}_{\min }$, a temperature that corresponds to the blackbody equilibrium temperature T_bb at the assumed stellocentric distance, and a normalization parameter (e.g., mass, area, or fractional luminosity). These models are fitted using multinest (Feroz et al. 2009) as described in Yelverton & Kennedy (2018), and also using MCMC using emcee.

Zoom In Zoom Out Reset image size

We fit two models, both of which have the same stellar component and assume the same silicate optical properties. The first is a theoretical model, in which the dust size distribution is fixed to q = 11/6, the minimum grain size to ${D}_{\min }=1.2$ μm, and the dust blackbody temperature to 40 K (as expected at 62 au, the best-fit parent belt location from modeling the NOEMA data). The only parameter for the dust component in this model is the normalization, and it should fit the data if these fairly standard assumptions are correct. The second model is the same, but parameters are allowed to vary.

The best-fit stellar model in both cases has an effective temperature of 5280 ± 100 K, consistent with our analysis in Section 4. The luminosity is L_* = 1.4L_☉, for which the blowout grain size is 1.2 μm. The dust models are shown in Figure 4; the red line shows the theoretical model given the parent belt location as imaged by NOEMA (F_ν,pb), while the semitransparent lines show a handful of dust models from MCMC fitting that are consistent with the data (F_ν,mcmc). The posterior distributions of the silicate belt parameters are shown in Figure 4, which also shows the parameters of the theoretically expected model. The theoretical model is inconsistent with the posterior distributions for the fitted model, whose size distribution slope is significantly shallower. This discrepancy is caused by the NOEMA flux density being higher than predicted given the mid-IR excess seen with WISE; this pushes the fitted models toward shallower size distributions, whose behavior can be understood in terms of the models presented by Kennedy & Wyatt (2014), e.g., Figure 11 in that work.

The astrosilicate model used here therefore illustrates that there may be a discrepancy between the mid-IR and millimeter-wave data, in terms of our ability to explain both with a dust population at a single stellocentric radius. This might be resolved with deeper millimeter-wave imaging to better locate the belt; a more distant location (i.e., a cooler T_bb) would allow a steeper size distribution to fit the data (Figure 4). Alternatively, the dust spectrum may not come from a single location, for example, comprising both warm inner dust as well as cool outer dust. Nonetheless, we use our fitted model of the infrared excess to calculate a fractional luminosity of ${L}_{\mathrm{IR}}/{L}_{\mathrm{tot}}={6}_{-1}^{+2}\times$ 10⁻⁴, higher than that of the debris disk around AU Mic (L_IR/L_tot ≃ 3.5 × 10⁻⁴; Matthews et al. 2015).

No significant evidence for point sources that may be self-luminous, massive companions to BD+45° 598 was found within our Keck scattered light images. To assess the sensitivity of our near-infrared coronagraphic observations to substellar companions in the vicinity of BD+45° 598, we start with the angular measure of our achieved 5σ contrast relative to the host star (i.e., a contrast curve). Here we present the contrast performance for only the 2013 September 25 epoch, as this data set delivered the best constraints on planetary mass companions that might be dynamically interacting with the disk. We then use this measure of our 5σ contrast, together with the models of Phillips et al. (2020) assuming an age of 23±3 Myr for the BPMG, to calculate the lower mass limit of substellar objects that would be detectable with our observations described in Section 2.1. Figure 5 shows this sensitivity to substellar companions as a function of projected orbital separation from the host star. Figure 5 also shows a detection probability map (right panel) based on our estimated contrast curve, and generated using the Exo-DMC algorithm (Bonavita 2020), demonstrating that our observations are sensitive to ∼8 M_Jup objects at 50 au, and ∼3 M_Jup objects at ∼100 au. In the left panel of Figure 5, we also indicate the upper and lower estimates of the age of the BPMG (27 and 21 Myr, respectively); however, these different ages do not substantially change our overall sensitivity.

Zoom In Zoom Out Reset image size