High-resolution ALMA and HST imaging of κCrB: a broad debris disc around a post-main-sequence star with low-mass companions Free

ABSTRACT

1 INTRODUCTION

The debris discs that have been observed around other stars are inferred to be dust, produced in the collisional cascades of exo-asteroid and exo-comet belts. Impacts between large belt objects form smaller planetesimals and observable quantities of debris dust at infrared and sub-millimetre/millimetre wavelengths (Wyatt 2008; Hughes, Duchêne & Matthews 2018). Debris discs have been observed with features such as cavities, gaps, warps, clumps, radial asymmetries, and gas (such as carbon monoxide, CO) as a result of belt collisions and/or interactions between planets and belts (Kral et al. 2017; Marino et al. 2018; Faramaz et al. 2019; Matrà et al. 2019; Lovell et al. 2021). Consequently, such features mean that studies of debris discs can constrain broader aspects of planetary systems such as their architecture and evolution that may not otherwise be observable.

Exo-Jupiter planets, which we define here as a class of giant planets with |$0.1\!-\!5\, M_{\rm {Jup}}$| and orbital radii between 1 and 5 au, are observed around ∼5 per cent of stars (Winn & Fabrycky 2015). Such planets are regularly determined to have eccentric orbits (Chiang & Laughlin 2013), an outcome predicted to emerge from planetary system instabilities (Jurić & Tremaine 2008) and/or the Kozai mechanism (Nagasawa, Ida & Bessho 2008). Debris discs at tens of au can be depleted by eccentric exo-Jupiters and simultaneously influence other closer-in planets dynamically (Gomes et al. 2005; Raymond et al. 2011, 2012), yet the presence of exo-Jupiters appears thus far to have little to no impact on the occurrence of detectable debris discs (Bryden et al. 2009; Moro-Martín et al. 2015; Yelverton, Kennedy & Su 2020). Despite this, exo-Jupiters and other outer massive planets can strongly influence the structure of planetesimal belts over long Myr–Gyr time-scales (see e.g. Mustill & Wyatt 2009) that can be imaged with high-resolution observatories, such as ALMA and Hubble Space Telescope (HST). Recently, Lovell et al. (2021) presented new ALMA images of the debris disc around q¹ Eri, an F9 main-sequence star that hosts an exo-Jupiter, and connected the disc morphology to a plausible planetary system architecture. Studying the morphology of debris discs at high resolution around stars hosting exo-Jupiters can therefore offer avenues to explore the planetary architecture and exoplanetary system evolutionary histories.

Debris discs are brighter and have higher detection fractions around main-sequence A-type stars, spanning their main-sequence lifetimes (with fractional occurrence rates of |${\sim }20{{\ \rm per\ cent}}$| over the range 5–850 Myr; see e.g. Rieke et al. 2005; Su et al. 2006; Wyatt et al. 2007). Therefore, the planetary systems of retired A-type stars provide useful laboratories to investigate long-term (⁠|${\gtrsim }\,$|Gyr) planet–disc interactions, for the old evolved systems that host at least one known planet and at least one observable debris disc. As part of a study of 36 such evolved post-main-sequence stars, Bonsor et al. (2014) presented evidence of four sub-giants hosting a far-infrared excess indicative of a debris disc (observed at 100 and 160 |$\mu$|m with Herschel), of which three host at least one confirmed giant planet (HR 8461, HD 131496, and |$\kappa \,$|CrB, all of which were once A-type stars). At a distance of 30.1 pc (Gaia Collaboration et al. 2018), of the known sub-giants that host at least one planet and debris disc, |$\kappa \,$|CrB is the closest (by a factor of 2–4) and hosts the brightest disc. Although all four known sub-giants with debris discs are at the base of the giant branch (see Fig. 1), |$\kappa \,$|CrB’s closer distance and far-infrared disc brightness make it the most ideal system to investigate the long-term effects of planet–belt interactions around a post-main-sequence star.

Figure 1.

Hertzsprung-Russell (HR) diagram showing three stellar tracks (from the models of Hurley, Pols & Tout 2000) for stars with Z = 0.025 and masses of 1.8, 1.5, and 1.3 M_⊙, respectively, for the solid blue, dash–dotted amber, and dashed green lines and the four known post-main-sequence sub-giants that host debris discs (Bonsor et al. 2014). Error bars represent 1σ and 3σ uncertainties (assuming |$\sigma {=}2{{\ \rm per\ cent}}$|⁠, see Section 3.3). All four are located near the base of the giant branch based on their ‘Priam’ and ‘FLAME’ modelled effective temperatures and luminosities (i.e. from Gaia Collaboration et al. 2018). All four are consistent with A-type stellar evolution tracks.

Open in new tabDownload slide

|$\kappa \,$|CrB (HD 142091, HIP 77655) is a K1 sub-giant star (an evolved A-type main-sequence star), with a mass |${\sim }1.8\, {\rm M}_\odot$| and age >2.5 Gyr (Johnson et al. 2008). |$\kappa \,$|CrB has been known to be an evolved sub-giant for many decades (based on spectral photometry, surface gravity, and absorption features; see e.g. Keenan & Pitts 1980; McWilliam 1990). It hosts an exo-Jupiter, |$\kappa \,$|CrB b, with semimajor axis |$a_{\rm {pl}}=2.8\,$|au, mass |$M_{\rm {pl}}\, \sin {i}=1.6\, M_{\rm {Jup}}$|⁠, and eccentricity e_pl = 0.07 (see Johnson et al. 2008). |$\kappa \,$|CrB’s debris disc was first imaged at 100 and 160 |$\mu$|m at resolutions of 6.7 and 11.4 arcsec, respectively, with Herschel PACS (Photodetector Array Camera and Spectrometer) observations (Bonsor et al. 2013). The data were fitted with a disc with modelled inclinations in the range 59°–60° and position angles (PA) in the range 142°–145°, and either a wide (20–220 au) dust belt or two belts centred at 41 and 165 au. Although broad morphological constraints were placed on the disc, the low-resolution imaging meant that a number of questions were left open; is this a single or multiple component debris disc; and how have possible disc–planet interactions carved the disc?

We present here the first ALMA and HST images of |$\kappa \,$|CrB, providing strict constraints on this planetary system, with observational resolution ≳3 times higher than previous Herschel images in an attempt to unravel the evolution of this planetary system. In Section 2, we detail our methodology for reducing our ALMA and HST data sets, in Section 3 we present our observational analysis of this system, in Section 4 our methodology to model the data, in Section 5 derive upper limits on the CO flux, and in Section 6 we derive limits on the putative planetary architecture. In Section 7, we provide a detailed discussion of our results and their possible interpretations, and summarize our key findings and conclusions in Section 8.

2 DATA REDUCTION

2.1 ALMA data reduction

κ CrB was observed for ∼195 min on source in four scheduling blocks with ALMA in Band 6 (1.1–1.4 mm) during Cycle 7, using 49 antennas with minimum and maximum baselines ranging from 15.0 to 312.7 m (project 2019.1.01443.T, PI: Lovell). Each scheduling block observed for ∼49 min on source, and was conducted on 2019 December 4, 2019 December 20, and two times on 2019 December 21. In all cases, the correlator had three spectral windows centred on frequencies of 228.527, 215.527, and 217.527 GHz, for continuum observations with a bandwidth of 2.000 GHz and channel widths of 15.625 MHz, as well as a spectral window centred on a topocentric frequency of 230.569 GHz with a bandwidth of 1.875 GHz and channel widths of 976.562 kHz for CO J = 2–1 spectral line observations. The average frequency of these observations corresponds to our quoted wavelength of 1.34 mm (which herein we quote as 1.3 mm). The visibility data set was calibrated using the casa software version 5.6.1-8 with the standard pipeline provided by the ALMA Observatory. The plotms task in casa was used to examine the visibility amplitudes as functions of both time and uv-distance, and as a result, it was deemed that no data flagging would be necessary. The measurement set was purged of calibrator data and time averaged into 10 s bins, channel averaged into four channels per spectral window, and the casa statwt task was run to estimate the visibility weights based on the measured variance. Continuum imaging was conducted using the casa tclean algorithm with natural weighting (providing a good trade-off between resolution and signal-to-noise ratio), shown in Fig. 2 (left-hand panel) in which emission from both the central star and the inclined debris disc are clearly detected. The synthesized beam size in this image is 2.1 × 1.4 arcsec (with a beam position angle |$\rm {PA}_{\rm {Beam}}{=}173.1^{\circ }$|⁠, anticlockwise from north), i.e. a physical resolution of 63.2 × 42.1 au, which has clearly resolved the disc over multiple beams.

Figure 2.

High-resolution images of the κ CrB system. In both panels, north points up, east points left, a 7.5 arcsec dashed ring is plotted, dash–dotted lines indicate the PA (146° ± 5°), and stellar locations are marked in the image centres with asterisk symbols (red left, white right). Left-hand panel: Band 6 (1.3 mm) natural weight cleaned image, showing the star, the disc, and outer emission. This image has an rms of |${\sim }12.0\, \mu$|Jy beam⁻¹, and contour lines of ±2σ, ±3σ, ±5σ, and |${\pm }7\sigma$| emission. Right-hand panel: HST scattered light image at optical wavelengths showing dust nebulosity coincident with the ALMA emission, with an image rms of |${\sim }8\, \mu$|Jy arcsec⁻², with ALMA Band 6 contours overlaid.

Open in new tabDownload slide

2.2 HST data reduction

κ CrB was observed with HST/STIS (Space Telescope Imaging Spectrograph) in three orbits on ut 2014-02-03. With the 50CORON aperture, this 1024 × 1024 CCD camera delivers unfiltered visible-light images with 0.05077 arcsec pixel⁻¹. The star was occulted by the WEDGEA0.6 and WEDGEA2.0 coronagraphic positions, which have widths 0.6 and 2.0 arcsec, respectively. The former provides close-in imaging, but saturation limits the total integration time. In each orbit, we acquired 12 images with 2.5 s integration each at WEDGEA0.6. For deeper, wide-field imaging, we acquired nine images with 130 s integrations each at WEDGE2.0 during each orbit. To facilitate subtraction of the optical point spread function (PSF), the telescope PA was changed by 17° between the first orbit and the second orbit, and 11.83° between the second and third orbits. A PSF reference star, HD 139323 (K3V), was observed with STIS at the WEDGE0.6 and WEDGE2.0 positions in a fourth orbit.

We registered the geometrically corrected, cosmic ray rejected files (*sx2.fits) by subtracting all images for a given wedge position from the very first image obtained in each orbit. By iteratively shifting the images by 0.01 pixel, we determined the x, y shift that minimized the PSF residuals after subtraction. To subtract the sky background, we determined the median pixel value in a 100 × 100 pixel box located in the top left-hand corner of each image. Every image was then divided by the integration time to provide units of counts sec⁻¹ pixel⁻¹. The images at each wedge position in each orbit were then median combined.

Each final κ CrB image for each wedge position in each orbit was then subtracted by the final image for each wedge position in the orbit that targeted HD 139323. The x, y position and flux scaling of HD 139323 relative to κ CrB was iteratively adjusted to minimize the subtraction residuals in the annular region containing the PSF halo and telescope diffraction spikes. Additionally, the WEDGEA2.0 data for κ CrB have an electronic or scattered light artefact appearing as a positive-valued horizontal band across the entire detector to the left and right of the occulted star. The band encompasses ∼37 rows, and it is uniformly brighter (median value 0.07 counts s⁻¹ pixel⁻¹) to the left of the star than to the right (0.04 counts s⁻¹ pixel⁻¹) as measured in the PSF-subtracted images. The steps to convert between counts s⁻¹ pixel⁻¹ and |$\mu$|Jy arcsec⁻² are provided by Schneider et al. (2016), for which we note that a gain of 4 and the pixel size stated above yield 1 counts s⁻¹ pixel|$^{-1}{=}177\, \mu$|Jy arcsec⁻². We subtracted the horizontal band using these values for rows 662:702 and 666:700 to the left and to the right of the star, respectively. Before PSF subtraction, the intersections of the diffraction spikes were used to determine the stellar locations behind the occulting wedge positions. The stellar location was then used as a centre of rotation for the PSF-subtracted images to be oriented with north up and east left. These north-rotated images were then median combined to provide two final PSF-subtracted images of κ CrB at the WEDGEA0.6 and WEDGEA2.0 positions representing 90 and 3510 s cumulative integration time, respectively.

Neither STIS image of κ CrB shows evidence for distinct morphological features such as rings or belts surrounding the star. The WEDGEA0.6 data have negative values up to 1.7 arcsec radius, and then positive values out to 4.1 arcsec radius. The negative halo may be a processing artefact: Changes in the thermal environment of the telescope between orbits change the PSF such that the PSF reference star is mismatched with that of the science target. The morphology of the positive halo is elliptical with a major axis in the south-east–north-west (SE–NW) direction. However, because the WEDGEA0.6 position is near the edge of the detector, the κ CrB image is vignetted at 2.7 arcsec radius in the PA range 45°–90°.

Our analysis will therefore focus on the WEDGEA2.0 data that are not vignetted, and are deeper, but have a larger inner working angle (IWA; ≥1.7 arcsec or 51 au) than the WEDGEA0.6 data. This image (Fig. 2, right-hand panel) does not have a negative-valued inner halo. Instead, there is a positive halo with a flux deficit north-east (NE) of the star, producing an overall fan-like morphology. A key concern is whether or not the deficit could be due to the occulting bar of the camera that crosses the field from the NE to the south-west (SW). This is an unlikely explanation because the bar would also produce a light deficit symmetrically across the star between the NE and SW sides. Also, the positive halo cannot be due to a residual of mismatched PSFs because such residuals would be azimuthally symmetric. Debris discs inclined to the line of sight can exhibit fan-like morphologies when surrounded symmetrically by dust, with an asymmetric scattering phase function (Kalas & Jewitt 1996). We thus conclude that the fan-like light morphology after PSF subtraction is due to light scattered by the dust grains surrounding κ CrB.

3 OBSERVATIONAL ANALYSIS

3.1 ALMA analysis

All analysis presented in this section is derived from the image of κ CrB as shown in the left-hand panel of Fig. 2, and we start by producing profiles of the emission. However, in order to produce accurate radial profiles, we first require estimates of the disc |$\rm {PA}$| and inclination (i). We start by measuring the |$\rm {PA}$| measured anticlockwise from north centred on the star by plotting the total flux within 5° wedges within a distance of 250 au (to enclose all disc emission from the stellar centre as a function of angle between 0 and 180° in 2° increments) and found this was maximized for a |$\rm {PA}{=}146^\circ{\pm }5^\circ$|⁠, consistent with previous work (i.e. 142°–145° as modelled by Bonsor et al. 2013). We define the line parallel with the PA through the centre of the star as the disc major axis, plotted in both panels of Fig. 2 along with its associated uncertainty. Given the inclined nature of this disc, the brightness of the emission coincident with the star, and the beam size, emission along the minor axis is fainter, and not well resolved. To measure the inclination, we therefore subtracted the unresolved stellar emission using a 2D Gaussian, defined by the beam major and minor axes, beam PA, and peak pixel intensity, from the centre of the Band 6 image and plotted elliptical rings around the residual emission, for which |$i{=}\arccos {\rm {(major\, extent / minor\, extent}}$|⁠). We estimate the disc inclination as 61° ± 5°, by varying the inclination such that all ellipses remained within the |$3 \sigma$| contour lines (also consistent with previous work, i.e. ∼59° as modelled by Bonsor et al. 2013). Whilst these inclination and PA values come with imaging uncertainty, we deem these sufficiently accurate for the purposes of our observational analysis. In Section 4, with more detailed modelling we show that these are consistent with best fits of the interferometric visibilities.

Fig. 3 shows the azimuthally averaged radial profile centred on the star, using our derived inclination of 61° and PA of 146°. Clear from this profile is that the peak emission is centred on the star, and that disc emission extends from ∼50 to 205 au (1.7–6.8 arcsec) based on the shoulder feature that falls below the 3σ level at ∼205 au, but peaks at around |${\sim }110\,$|au (3.7 arcsec). To place constraints on the stellar and disc emission, overplotted are the profiles of (i) a point source with the same brightness as the central peak located at the origin (as a blue dotted line), and (ii) a narrow ring of emission with the same brightness as the average belt emission at a deprojected radius of 110 au (as an amber dotted line), with both having been convolved with the average beam extent. Since the star is unresolved, the emission peak at this location can be used to constrain the 1.3 mm stellar flux as |$191{\pm }8\, \mu$|Jy, which combined with the 10 per cent flux calibration uncertainty of ALMA results in |$F_{\rm {\star ,\, B6}}{=}191{\pm }20\, \mu$|Jy.

Figure 3.

Azimuthally averaged radial profile for the Band 6 data (red, solid), along with curves showing emission coincident with the star and the disc at a deprojected distance of 110 au convolved with the beam (blue and amber, respectively, dotted), their combination (green, dashed), and a black dash–dotted line representing |$3 \sigma$| emission. The shaded regions represent |$1 \sigma$| and |$2 \sigma$| errors.

Open in new tabDownload slide

The combined disc and star profile in Fig. 3 shows that whilst a simple point source located at the stellar location and a narrow ring at 110 au can broadly account for the emission internal to ∼150 au, beyond this distance the outer edge emission is too bright, i.e. a narrow ring cannot fully interpret the disc towards κ CrB. This could be the case if the disc is broad. We discuss further the potential origins of this emission signature at the end of this section. It can be seen in the left-hand plot of Fig. 2 that there is a faint clump of emission at an approximate relative RA–Dec. location of [−1.0, −5.0] (i.e. at an approximate deprojected distance of |${\sim }205\,$|au), with a few tens of |$\mu$|Jy beam⁻¹ brightness. Given that the radial profile is consistent with noise at all radii beyond ∼205 au, this clump appears too faint to have significantly increased our estimate of the disc’s extent. We discuss the possible origins of this emission in Section 7.1.3.

The same conclusions from Fig. 3 can be drawn from the integrated flux profile shown in Fig. 4. Here, we show that the integrated flux rises to a peak value at ∼276 au, consistent within 1|$\sigma$| with the integrated flux beyond 180 au. Although the average disc intensity falls below 3σ beyond ∼205 au in Fig. 3, given that the integrated flux profile is broadly constant beyond 180 au, we interpret 180 au as the disc outer edge in the millimetre. Fig. 4 further demonstrates that at the depth of our observations there is no detectable asymmetry in the total flux in the two disc ansae, given that the integrated flux profiles in the NW and SE halves of the disc are consistent over all radii. The flux profile can additionally be used to measure the total 1.3 mm flux from κ CrB based on the peak of the ‘Total’ value of |$0.96\pm 0.04\,$|mJy. Combined in quadrature with an assumed |$10{{\ \rm per\ cent}}$| ALMA calibration uncertainty, this results in a value of |$F_{\rm {tot,\, B6}}{=}0.96\pm 0.10\,$|mJy (which we include in Table 1, and later use to model the flux distribution in Section 3.3).

Figure 4.

Band 6 flux profile as a function of deprojected radius of |$\kappa \,$|CrB. Highlighted regions represent |$1 \sigma$| error. Here, the total flux can be measured. Split either side of the disc minor axis, this plot demonstrates that the total flux is symmetric about the star at the depth observed.

Open in new tabDownload slide

3.1.1 Summary of ALMA analysis

In summary, we conclude from this ALMA analysis that having detected the debris disc towards |$\kappa \,$|CrB, we have observationally constrained its PA on the sky as 146° ± 5°, its inclination as i = 61° ± 5°, the radial peak location of its disc emission at ∼110 au, starting at |${\sim }50\,$|au, and extending outwards to at least ∼180 au.

There are a number of ways that |$\kappa \,$|CrB’s disc emission could arise based on the resolution of our imaging, which we will explore by modelling the dust belt in Section 4. Whilst our observational analysis suggests that the broad disc emission arises from a belt of planetesimals, it remains unclear how broad the underlying planetesimal distribution may be. We will test this with simple models that parametrize the surface density distribution of dust as a function of the radial distance from the star. These models will be able to (1) constrain the extent of the belt, and by comparing different parametrizations, (2) search for the presence of any underlying radial sub-structure in the disc. We note that whilst broad debris discs can be modelled in many ways, here we will use just two parametrizations. We either model the disc with (i) a radial power-law distribution (i.e. ‘model SPL’, with free parameters for the inner and outer edges) or (ii) a radial Gaussian distribution (i.e. ‘model SGR’, with free parameters for the radial peak and width/standard deviation of the Gaussian). In conjunction, these will offer sufficient insight to explore the debris disc and underlying planetesimal belt of |$\kappa \,$|CrB, which we fully described in Section 4.

3.2 HST analysis

All analysis in this section is derived from the HST scattered light image of |$\kappa \,$|CrB as shown in the right-hand panel of Fig. 2. We start by noting that all modelling of the scattered light data will be conducted in Section 4.2, and our comparison of the mm and scattered light observations will be discussed in Section 7. To first investigate the morphology of the disc that can be inferred from the image, we produced a radial profile of the scattered light by averaging emission within 10° wedges either side of the star, assuming a PA of 146° and an inclination of 40° (which we discuss further and justify in Section 4.2). This profile is shown in Fig. 5. We note that the derived profiles with 10° wedge angles are consistent with a wide range of assumed inclinations and PA, including those derived from the ALMA image. Such consistency demonstrates that there are large uncertainties present when estimating the geometry of discs from scattered light emission. The radial profile (designated ‘NW–SE Data’) allows us to estimate the peak brightness, and the inner and outer edges of the dust-scattered light. Together these can constrain the dust emission extent. Further, since the surface brightness of scattered light images is proportional to the dust opacity, surface density of micrometre-sized grains, the scattering phase function, and stellar irradiation, by scaling the surface brightness by r², the resulting profile corrects for the stellar irradiation (which follows an inverse square law; see e.g. Stolker et al. 2016), and is thus a better approximation of the surface density distribution of the micrometre-sized grains. We therefore plot this same radial profile scaled by r². Finally, we estimate these two emission profiles with decaying power-law models, which we discuss further below, and return to in Section 4.2.

Figure 5.

Radial profile for the HST scattered light emission (blue, solid), the same but scaled by r²/10 000 (green, dotted), and modelled power-law profiles (amber, dashed), showing the decaying nature of the emission, over a broad radial distance from the star. The errors in scattered light data are |${\pm }3\sigma$|⁠, and the vertical red dashed line is at the radial location within which coronagraphic obscuration dominates the emission.

Open in new tabDownload slide

Analysing Fig. 5 in detail, we start by noting the WEDGEA2.0 IWA at 51 au limits our ability to analyse emission within this radius (shown with a red dashed line). The brightness profile shows an initial ‘shoulder’ of emission (between 51 and 67 au), a sharp brightness peak (at 73 au, with intensity |${\sim }100\, \mu$|Jy arcsec⁻²), and a decaying intensity profile outwards to ∼280 au. In terms of the r²-corrected profile, the dust distribution has the same shoulder, and follows a relatively flat power-law profile between ∼73 and 280 au, beyond which this becomes consistent with noise. Based on the unknown true inner edge of the dust, we conclude that the scattered light from circumstellar dust has a wide extent spanning ∼51–280 au.

The radial profiles allow us to set two additional broad constraints on the dust distribution, and albedo which we discuss in turn here. First, starting at the 73 au peak, by overplotting an intensity profile |$I(r) = I_0 \times (r/R_{\rm {in,\, {\it HST}}})^{-\gamma }$| (where I₀ is set to the inner edge intensity, i.e. |$I_0{=}100\, \mu$|Jy arcsec⁻²) and I(r) × r², we show that the brightness and surface density fall-offs are consistent with simple power laws. We show these for γ ∼ 0.2 for the surface density profile, and equivalently γ ∼ 2.2 for the surface brightness profile, from which we conclude that over the extent of the 73–280 au emission, the fall-off in scattered light emission is pronounced, with a dust surface density consistent with being relatively flat. Secondly, from the peak scattered light brightness, we estimate the required albedo of the dust [using equation (2) of Marshall et al. (2018), i.e. S = fFA/((2|$\pi$| × ϕ × dϕ × cos i)(1 − A)), where S is the surface brightness in mJy arcsec⁻², f is the fractional luminosity, F is the stellar flux in mJy, ϕ is the disc semimajor axis in arcsec, dϕ is the disc width in arcsec, i is the inclination, and A is the albedo). We find |$A{\sim }10{{\ \rm per\ cent}}$|⁠, based on F = 2923 mJy, f = 5 × 10⁻⁵ (see Section 1), dϕ = 2–8 arcsec, and i = 40° (we note A varies by just |${\sim }2{{\ \rm per\ cent}}$| for a higher inclination of 60°). This albedo is broadly consistent with our expectation based on the broader flux distribution (see Section 3.3, i.e. since the dust appears to be bright in its thermal emission) and consistent with the dust albedo of other A-type stars (e.g. Fomalhaut; see Kalas, Graham & Clampin 2005).

Although faint, we note here the tentative possibility of a scattered light asymmetry, with the dust emission more extended to the NW of the disc than the SE, most easily seen in Fig. 2 (right-hand panel). At the same intensity of emission, the NW extends out to ∼7.0 arcsec, whereas in the SE, this same intensity emission extends only as far as ∼5 arcsec. Although we do not extend analysis of this tentative asymmetry further, we suggest future observations in scattered light may better constrain this with more sensitive measurements.

3.3 Flux distribution

Fig. 6 shows the flux density distribution of κ CrB, including the flux density derived from the Band 6 ALMA data (see Section 3.1), for which wavelengths, fluxes, and magnitudes (up to K-band photometry) are given in Table 1. To investigate where the emission in this system originates from, we fitted the complete flux distribution with a stellar photosphere plus a ‘modified’ blackbody, applying the methodology of Kennedy & Wyatt (2014) and Yelverton et al. (2019). We find that the distribution is best fitted with a single, outer, cool temperature dust component, and a star with a stellar effective temperature of |$T_{\rm {eff}}{=}4850{\pm }100\, \rm {K}$|⁠, and a stellar luminosity and radius of |$L{=}12.0\, L_{\odot }$| and |$R{=}4.9\, R_{\odot }$|⁠, respectively (with 2 per cent calibration uncertainties). We note that this modified blackbody model is simply a normal blackbody over all wavelengths, until a specific wavelength, λ₀, whereby the flux density is then modified by a term |$(\lambda /\lambda _0)^{-\beta _0}$|⁠, and that the stellar spectrum model covers the wavelength range from 0.05 |$\mu$|m to 1.0 cm. We further note that Gaia Collaboration et al. (2018) estimated these stellar parameters (using the ‘Priam’ and ‘FLAME’ algorithms) as |$4750\, \rm {K}$|⁠, |$12.7\, L_{\odot }$|⁠, and |$5.1\, R_{\odot }$|⁠, respectively, which are all broadly in agreement with these model-fitted values. Error flags on the WISE1 and WISE2 data points meant that these were excluded from the fitting procedure (though these are still plotted in Fig. 6).

Figure 6.

Flux density distribution for |$\kappa \,$|CrB, where the total flux measurements are dominated by emission due to the star up to |${\sim }30\, \mu$|m, and by the disc emission between |${\sim }30\, \mu$|m and 1.3 mm. The disc is modelled here as a single component modified blackbody as described in Section 3.3, which successfully interprets the mid-infrared and new ALMA Band 6 data.

Open in new tabDownload slide

The single-component model finds the disc is best fitted with a blackbody temperature of |$T_{\rm {bb}}{=}72{\pm }4\, \rm {K}$| (corresponding to an uncorrected blackbody radius of r_bb = 51 ± 6 au), a relatively faint but significant fractional luminosity of 4.9 ± 0.5 × 10⁻⁵, and modified blackbody parameters β₀ ∼ 1.6 and |$\lambda _0{\sim }400\, \mu$|m. The profile turn-over at 400 |$\mu$|m can be seen clearly in Fig. 6, though we note here that with only one data point beyond |${\sim }160\, \mu$|m, λ₀ and β₀ are poorly constrained and degenerate. It is well known that blackbody radii measured from SEDs (Spectral Energy Distributions) are smaller than the radii determined from resolved millimetre wavelength images (see e.g. Pawellek & Krivov 2015). We also find this to be the case here: whilst the uncorrected blackbody radius is ∼51 au, the peak millimetre-emission radius is ∼110 au, and so a factor of ∼2 different. In Pawellek & Krivov (2015), the ratio of these two radii measurements, Γ, is empirically quantified in equation (8) with respect to the stellar luminosity and dust composition. For a range of compositions (see table 4 of Pawellek & Krivov 2015) and the previously determined stellar luminosity of |$L{=}12.0\, L_{\odot }$|⁠, we find a Γ ratio of 2.1–2.3. Since there is good consistency between Γ × r_bb = 107–117 au and the ∼110 au peak millimetre-emission radius, this demonstrates that our measurements align well with those of discs around similarly luminous stars.

We note that unresolved millimetre emission coincident with the location of the star could exist in the warmer inner regions of this system, for example, from a separate inner warm/hot dust component near to the star. Indeed, Absil et al. (2013) detected a near-infrared K-band excess at the 1 per cent level, which was interpreted as due to the evolutionary stage of the star. However, we are unable to verify the presence of a millimetre excess at 1.3 mm co-located with the star. To demonstrate this, we note that our modelled stellar flux distribution has a value of 220 |$\mu$|Jy at λ =1.3 mm, which is consistent within 1.5σ with the measured flux at the location of the star |${\sim }191{\pm }20\, \mu$|Jy as measured in Section 3.1 and shown in Fig. 6. Consequently, any near-infrared excess towards |$\kappa \,$|CrB (such as that measured by Absil et al. 2013) must be sufficiently faint to have avoided detection at millimetre wavelengths. We thus conclude that across the complete flux distribution of κ CrB, our modelling only finds significant evidence for emission due to a photospheric component and a single disc at many tens of au.

4 MODELLING κ CrB’S DUST OBSERVATIONS

Following our observational analysis in Section 3, here we assess the morphology of κ CrB’s debris disc with parametric models. We start by discussing our Band 6 model set-up in Section 4.1.1, their results and implications in Section 4.1.2, explore how this relates to a possible model of the dust emission in scattered light in Section 4.2, and summarize our key modelling findings in Section 4.3.

4.1 Modelling the millimetre emission

4.1.1 Millimetre model set-up

The models described in Section 3.1.1 share a common set of seven free parameters: dust mass (M_dust), main belt radius (r₁), inclination (i), |$\rm {PA}$|⁠, photospheric Band 6 flux (⁠|$F_{\star ,\, \rm {B6}}$|⁠), and the Band 6 measurement set phase centre offsets in RA and Dec. (x_off, y_off, respectively). We modelled the system with two belt distributions, ‘model SPL’, defined by an inner edge (r₁), outer edge (r₂), and power-law profile, with an index (p_SPL) for a total of nine free parameters, and ‘model SGR’, defined by a peak emission radius (r₁) and Gaussian width (σ_SGR, the standard deviation of the Gaussian distribution) for a total of eight free parameters. Although we measured the Band 6 emission coincident with the star as ∼191 |$\mu$|Jy in Section 3.3, since stellar emission in the millimetre can be poorly constrained, we leave the flux located near the star (⁠|$F_{\star ,\, \rm {B6}}$|⁠) as a free parameter. The models have the same vertical Gaussian density distribution, defined by the scale height, h = H/r, where H is the height of dust above the mid-plane at a radius r (Marino et al. 2016). Given our observations only have a low signal-to-noise ratio and the disc is only moderately inclined, we fix h = 0.05 (consistent with the expected range of h = 0.02–0.12 for debris discs; Hughes et al. 2018), and further note that by investigating model fits with h left as a free parameter, we always yielded unconstrained values. Given the lack of any observed asymmetries in our images, we do not model particle eccentricities. Additionally, all these models have identical (and fixed) minimum and maximum grain radii of |$a_{\rm {min}}{=}1.0\, \mu \rm {m}$| and |$a_{\rm {max}}{=}1\, \rm {cm}$|⁠. These are set at the lower limit by the approximate blowout size, which is just a few micrometres for astrosilicate-dominated grains around a star with |$1.8\, {\rm M}_\odot$| and |$12\, L_\odot$| (see e.g. Wyatt 2008, equation 14; D_blowout = 0.8(L_⋆/M_⋆)(2700/ρ)), and at the upper limit by neglecting emission from grains a factor of 10 larger than the Band 6 wavelength of 1.3 mm. We further assume a dust grain density of ρ = 2.7 |$\rm {g}$| |$\rm {cm}^{-3}$|⁠, a size distribution with power-law exponent α = 3.5 (Dohnanyi 1969), and a weighted mean dust opacity of |$\kappa {=}1.88\, \rm {cm}^2\, \rm {g}^{-1}$|⁠, based on a mix of compositions with mass fractions of 15 per cent water-ice, 70 per cent astrosilicates, and 15 per cent amorphous carbon (see both Li & Greenberg 1998; Draine 2003). We note this composition was likewise used to model the >Gyr old debris discs towards HD 107146 and HD 10647 (see Marino et al. 2018; Lovell et al. 2021, respectively). Since we do not investigate the material composition or size distribution, a_min, a_max, ρ, α, and κ remain fixed throughout.

4.1.2 Millimetre modelling results

We used the radmc-3D package (Dullemond et al. 2012) with the MCMC (Markov Chain Monte Carlo) parameter estimator emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013) to compute the dust temperature of our disc models. These dust temperatures were used to create 2D images at |$1.3\,$|mm and subsequently used to compute model visibilities at identical uv-baselines to our ALMA observations. This was done identically to Lovell et al. (2021), using the tools developed in Marino et al. (2018) and Lovell et al. (2021), where a complete discussion of model fitting is presented.

All model best-fitting results are reported in Table 2, with their corresponding model and residual images shown in Fig. 7. We start by noting that the four parameters determined in Section 3.1, i.e. the main belt radius (r₁), |$\rm {PA}$|⁠, inclination i, and Band 6 stellar flux |$F_{\star \, \rm {B6}}$|⁠, are all found to within |$3 \sigma$| of both best-fitting Band 6 model parameter sets. Note that whilst the power-law model SPL has a sharp inner edge found to be at |$93{\pm }8\,$|au, given this model is broad, when this is imaged and convolved the same ALMA beam as our measurements, this has a peak emission radius at |${\sim }110\,$|au. These models also suggest that there is |${\sim }0.016\, {\rm M}_\oplus$| of dust that will be dominated by the ∼centimetre-sized grains in the model. We discuss the implications of this dust mass further in Section 7.1.1. In addition, our modelling suggests that there may be a phase centre offset given all x_off and y_off parameters are |${\gtrsim }1\sigma\!-\!2 \sigma$| different to 0; however, since there are no situations in which any of these exceed |$3 \sigma$| (in any single axis), we do not discuss this any further. Finally, we note that in both residual images, a faint, unresolved clump is present at relative RA, Dec. position of [−1.0, −5.0], with a flux density |$30\!-\!40\, \mu$|Jy, consistent with the emission of the same clump in Fig. 2 if due to an unresolved point source. Since this is present in the residuals of both images, our model fits have likely not been strongly influenced by the presence of such a clump. We discuss the possible origins of this clump emission further in Section 7.1.3.

Figure 7.

Best-fitting parameter model and residual images (best-fitting model subtracted from the data) for both Band 6 models investigated. Upper plot contours (model images) are +3σ, 5σ, and 7|$\sigma$|⁠, and lower plots (residual images) are ±2σ and 3|$\sigma$| (red positive, blue negative). In both residual plots, the beam is shown in black in the lower right-hand of the panels, and a circle with radius 180 au inclined to 60° is shown with a black dashed loop.

Open in new tabDownload slide

4.1.3 A broad debris disc

Deriving best-fitting parameter estimations provides us with useful constraints on the physical morphology of the debris disc; however, our modelling was primarily conducted to test whether or not our data provided evidence that the underlying planetesimal belt is extended or narrow. To quantify this, by fitting models that can measure either the width of the emission (if fitted as a Gaussian) or the extent of the emission from the distance between the inner and outer edges (if fitted as a power law), both our models show evidence that the emission is indeed broad. Discussing these models in turn, we first find that the Gaussian disc attempts to fit an extremely broad disc with a radially symmetric profile with a width (standard deviation) of |$\sigma _{\rm {SGR}}=98{\pm }19\,$|au. Given the width of model SGR is so broad (symmetric about ∼131 au, with a width ∼98 au), this simultaneously attempts to fit the inner regions where emission is likely dominated by the star, and the outer regions where emission may be entirely due to noise, and thus can be seen to oversubtract and undersubtract emission in various regions of the residual image (see Fig. 7). Whilst this model supports the hypothesis that the disc is broad, it may not be a good parametrization of this system. We quantify this further in Section 4.1.4 by comparing models SGR and SPL statistically.

In contrast, the model describing the belt with a decaying power law (SPL) provides a more promising interpretation, whilst also showing the emission to arise from a broad disc. We quantify this breadth based on the modelled disc outer edge (r₂), which tends to unconstrained high values from which we can set a lower limit. For this reason, we quote |$r_2{\gt }151\,$|au in Table 2, based on where |${\gt }99.7\, {{\ \rm per\ cent}}$| of all r₂ values finished in their MCMC chains. Comparing with the well-constrained inner edge of 93 ± 8 au, this confirms that the outer edge is indeed significantly distant from the inner edge, and thus the disc has a broad extent. We constrain this extent by finding the difference between the MCMC chains of the fitted inner and outer belt radii and find that over |$99.7\, {{\ \rm per\ cent}}$| of these are separated by 50.1 au (comparable with the resolution scale of our imaging, thus at least partially resolved), and assess this as a strict, significant lower limit on the belt extent. Based on the peak emission radius (∼110 au), this model would therefore place a fractional width on this disc of |${\gt }46{{\ \rm per\ cent}}$|⁠. We note here that if we instead use the 50th percentiles in the MCMC chains to find the median extent (instead of the 99.7th percentile), we measure a much broader extent of 135 au (since the median outer edge radius is at a wider value of 232 au). As such, the disc could well be much broader than the 50.1 au lower limit, which higher signal-to-noise ratio observations may measure.

To demonstrate further the broad extent of the underlying models, we show in Fig. 8 the surface density profiles of the two models, based on the 50th percentile of the MCMC chains after burn-in, along with their corresponding ±1σ and |${\pm }2 \sigma$| results. Both models show that their 50th percentile surface densities exceed 0 over a broad range of radii, i.e. that broad radial ranges are required for the modelled dust in order to interpret the Band 6 observations.

Figure 8.

Surface density plots for the two models. These are plotted based on a random sample of 5000 MCMC chains (after burn-in), and show the 50th percentile plots, ±1σ and 2|$\sigma$|⁠.

Open in new tabDownload slide

4.1.4 Millimetre model statistical comparison

Although qualitatively model SPL appears to be a better representation of the data, we additionally compare models SPL and SGR statistically. We do so by assessing their respective χ² and |$\rm {BIC}$| (Bayesian Information Criterion) values, with the former assessing the goodness of fit of the models, and the latter being a comparative assessment based on the size of our interferometric data sets and number of parameters in our model, i.e. |${\rm BIC} = \chi ^2 + N_{\rm {par}}\ln {N_{\rm {dat}}}$| (see Schwarz 1978; Kass & Raftery 1995, for which a difference in BIC values of 10 is considered strong evidence to support a given hypotheses, e.g. that one model is more statistically significant than another). For our Band 6 data, there are 2 × N_Vis = N_Dat = 12, 104, 328 independent data points, accounting for the real and imaginary components of the interferometric visibilities. We report in Table 3 these relative values, and show with reference to the SGR model that whilst the additional parameter included in the SPL model improves the fit to the data (i.e. model SPL has a lower χ²), this may not be significant (i.e. model SGR has a lower |$\rm {BIC}$|⁠). Whilst this may suggest that the increase in the number of parameters from 8 to 9 between models SGR and SPL is not favoured, Lovell et al. (2021) demonstrated that statistical comparisons of disc model fits to low signal-to-noise ratio interferometric data sets can be significantly skewed by large N_Dat. Consequently, as |$\Delta \rm {BIC}$| may not be a meaningful comparative diagnostic, this would suggest that model SPL is statistically favourable of the two models.

4.1.5 Alternative models?

Whilst the data can be fitted with a smooth wide disc, we cannot rule out the possibility that there is a gap or narrow radial sub-structure in the disc that gives rise to the observed emission structure. Albeit with less significant fits in additional tests, i.e. poorer χ² and |$\rm {BIC}$| values, and visibly poorer residuals, we found that the data can nevertheless be modelled by two narrow belts. Such a scenario comprising multiple planetesimal belts would not be unprecedented as it has been found that gaps in wide debris discs may be common (see e.g. Marino et al. 2020). However, given the low signal-to-noise ratio, the resolution of our observations, and that such a model was indeed a poorer representation, we do not assess this possibility any further in this work.

4.2 Modelling the scattered light emission

Whilst the emission in the millimetre traces well the large planetesimals within a disc, due to stellar radiation forces affecting much smaller grains, tracing the parent bodies with small particles is much more challenging. This means that the observed morphologies of discs in scattered light and in the millimetre can be substantially different, even if their dust formed in the same parent belts (e.g. see Beta Pic as imaged in scattered light with HST, and in the millimetre with ALMA; see Golimowski et al. 2006; Dent et al. 2014, respectively). Nevertheless, constraints on disc morphologies such as their PA and inclination can be obtained from scattered light imaging, and comparisons with discs as observed in other wavelengths probe different disc processes. As such, we outline here our investigation into the morphology of the scattered light emission for which we model the HST data in order to compare this with the emission observed in the millimetre images.

We start by noting that our scattered light HST model is simply a variant of the broad single power-law millimetre model (model ‘SPL’) described in Section 4.1.1, but instead set up to produce images of scattered light emission at 0.6 |$\mu$|m, with a range of anisotropic scattering factors (g) based on the Henyey–Greenstein (HG) scattering function (Henyey & Greenstein 1941). In this investigation, we fixed the inner and outer edges of the disc to 73 and 280 au, respectively (based on the peak brightness of the HST observations, and the outer edge at which scattered light emission was observed in Section 3.2), the PA to 148°, and produced a suite of models, based on variations of the disc inclination (ranging from 25° to 60° in 5° increments), scattering factors (ranging from 0.1 to 0.9 in 0.1 increments), and power-law indices (ranging from −2 to 2 in 0.5 increments). To find the best fit of the model suite, we introduced a dummy pre-factor that each model was multiplied by before the data were subtracted (used to scale the emission brightness to maximize the goodness of fit, which was also varied in order to minimize the χ² of each model). To avoid our goodness-of-fit estimates being driven by emission far from the star and by noisy inner regions, a null map was produced based on the data and convolved with all images, such that models contained identical imaging artefacts as the raw data. A noise estimate from the raw data was obtained from emission unaffected by real scattered light emission, i.e. beyond 280 au, with which the reduced χ², |$\chi ^2_\nu$|⁠, of each residual image (HST data minus model) could then be measured. We note here that our null model with zero emission everywhere has |$\chi ^2_{\rm {\nu ,\, null}}=3.75$|⁠, whereas the suite of models ranged in their measured |$\chi ^2_{\nu }$| values between 1.69 and 2.70. Further, by simulating grids of noise, we estimate the uncertainty in measured |$\chi ^2_{\nu }$| values as ∼0.10.

We show in the left-hand panel of Fig. 9 the HST data (unchanged from Fig. 2), in the centre panel the best-fitting model residuals, and in the right-hand panel the best-fitting model residuals for a disc inclined to 60° (consistent with that measured in the millimetre). We note that the overall best-fitting model has |$\chi ^2_{\nu } {=} 1.69$| and the best 60° inclined model has |$\chi ^2_{\nu } {=} 2.05$|⁠. This investigation found generally that to minimize χ², models required low inclinations in the range i = 30°–40°, highly anisotropic HG scattering factors in the range g ∼ 0.8–0.9, and relatively flat surface density power-law indices in the range γ = −0.5 to 0.5 (consistent with our imaging analysis earlier, see Section 3.2). We note here that (1) variations in the disc PA (between 142° and 152°) did not yield any significant improvement in the measured |$\chi ^2_{\nu }$| values, and (2) increasing the value of the scale height, h (to either 10 per cent, 15 per cent, or 20 per cent), did not affect the preferred inclinations of best-fitting model discs. The second of these points suggests that there is no measurable dependence between the scale height and the disc inclination. However, we do not rule out that different model parametrizations could yield different constraints on these parameters, and in particular different phase scattering functions, which we briefly investigate later in this section, and discuss further in Section 7.1.2.

Figure 9.

Scattered light images used to model the HST emission; in all three panels, north is up and east is left. Left-hand panel: HST image (as per Fig. 2). Centre and right-hand panels: residual images (data-model) for disc models with HG scattering factors of 0.9, power-law indices of γ = 0.0, and either inclinations of i = 35° or i = 60°, respectively. In the i = 35° model, |$\chi ^2_\nu {=}1.69$| (the lowest of all models), whereas in the i = 60° model, |$\chi ^2_\nu {=}2.05$|⁠, i.e. a less good fit. Also plotted are white dash–dotted ellipses demonstrating 200 au rings inclined at either 35° (centre panel) or 60° (right-hand panel).

Open in new tabDownload slide

Discussing these two residual images in turn, in the case of the best fit (centre panel), there remains emission in the region between 73 and 280 au; however, the broad azimuthal region covering the SE to the NW is reasonably well modelled. In addition, in the NE where the HST image artefacts are significant (but the absolute scattered light emission is lower), the model accounts for this with its high scattering factor, which naturally lowers the emission brightness on the far side of the disc. Therefore, whilst this minimal |$\chi ^2_{\nu }=1.69$| suggests that a better fit to the data can be found, this model does at least broadly account for the morphology of the scattered light emission, and is statistically favoured to the best-fitting model with i = 60° inclination. For the purposes of this work, we deem this i = 35° best-fitting model to be sufficiently accurate for comparative purposes. However, detailed modelling of |$\kappa \,$|CrB’s scattered light emission could be conducted in the future to investigate this in greater detail and understand the uncertainty distribution of best-fitting parameters. We explore some of the details that require further investigation in Section 7.1.2.

In the case of a model identical except for having a more inclined i = 60° disc, the measured |$\chi ^2_{\nu }$| is significantly higher (2.05), and it can be seen in the image residuals (right-hand panel of Fig. 9) that emission in the azimuthal region between the SE and NW is not as well modelled: there is brighter excess emission in the SW, and a deeper subtraction along the disc major axis in the SE and NW. This demonstrates that the model fits in scattered light favour an inclination that is shallower than the inclination measured and modelled in the millimetre (and those in the far-infrared; see Bonsor et al. 2013). By investigating a range of different HG scattering factors, inclinations, and power-law indices, we have shown that this is a trend seen in all disc models: with highly anisotropic scattering factors, relatively flat power-law indices, and low ∼30°–40° inclinations favoured.

To test whether the millimetre-scattered light inclination discrepancy arises from the HG scattering parametrization, we further investigated whether a more inclined scattered light model could fit the observations by considering spherical grains and Mie theory, consistent with the models presented in Pawellek et al. (2019). Mie theory was investigated since strong forward scattering is expected (see the appendix of Pawellek et al. 2019). In doing so, we again found that the scattered light models prefer low inclination values between 30° and 40°, consistent with the modelling approach presented above, and still discrepant with that observed with ALMA and Herschel. We explore the implications for this inclination discrepancy further in Section 7.1.2.

4.3 Summary of modelling

To summarize, by comparing the outputs of parametrized models fitted to the ALMA Band 6 visibilities, we find that |$\kappa \,$|CrB’s debris disc most likely arises from a radially broad distribution of planetesimals, which we constrained as having an extent |${\gt }50\,$|au. In particular, we showed that whilst radial power-law and Gaussian models can fit the data, the power-law model is favourable. However, given that the disc emission could have been broadened due to the presence of unresolved multiple planetesimal belts, future deep millimetre observations of the |$\kappa \,$|CrB system could be undertaken to examine alternative models.

Additionally, we produced simple parametric models of the scattered light emission to investigate the dust emission seen by HST. We showed with that the models that best explain the observed dust emission require the emission to be (i) highly anisotropic with HG scattering factors between g = 0.8 and 0.9, (ii) to have a relatively flat low power-law index between γ = −0.5 and 0.5, and (iii) to have low inclination angles between i = 30° and 40°. These models come with a range of uncertainties and are not able to fully interpret the HST emission, and thus could be greatly improved with detailed modelling and new scattered light data. Nevertheless, we have shown the disc to have a less inclined morphology in scattered light to that in the millimetre. We discuss the implications of the scattered light and millimetre modelling further in Section 7.

5 CO J = 2–1 ANALYSIS

Alongside measurements of the dust continuum, the ALMA Band 6 observations as introduced in Section 2 included higher spectral resolution measurements over the CO J = 2–1 spectral line frequency. Therefore, we produced a subset of the casa measurement set containing only data covering this spectral line with rest frequency f_CO = 230.538 GHz. To remove continuum emission from any |$\rm {CO}$| signal, we used the casa tool uvcontsub with a fit order of 1 external to the region defined by f_CO ± Δf, where Δf was set to a width |$30\,$|MHz, avoiding fitting to channels where |$\rm {CO}$| could be present. We then produced a data cube in the barycentric reference frame using the casa tclean algorithm with zero clean iterations, since there was no significant emission per beam present in any channel.

We note that the relative line-of-sight velocity of any CO gas emission (if present) would however fall across multiple data channels due to the variation of the observed radial velocity (RV) around the disc. This can however be accounted for since the RV of the star, the disc morphology, and stellar mass are all known, for which |$v_{\rm {rad}}{=}-24.98{\pm }0.14\, \rm {kms^{-1}}$| (see Gaia Collaboration et al. 2018). Given the disc is broad, any gas that is co-orbiting with the disc will span a range of velocities and thus fall over a range of channels. To measure the line profiles from our CO data cube, we therefore applied the spectro-spatial filtering method of Matrà et al. (2015, 2017a) by producing a Keplerian mask to shift pixels between cube channels based on a stellar mass of 1.8 M_⊙ (see Johnson et al. 2008) and the same disc PA and inclination as derived in Section 3.1 in both of the possible disc rotation directions (either clockwise and anticlockwise) and integrated all emission (in both rotation directions separately) between projected radii of either 0–60 au or 60–240 au from the star. These two distances cover the two regions dominated in the continuum by either (i) the star from 0 to 60 au, or (ii) by the disc from 60 to 240 au. Both spectra are shown as a function of RV in Fig. 10, which show that no significant CO J = 2–1 line emission is present towards κ CrB, in either of the spatial regions, and in either of the rotation directions.

Figure 10.

CO spectral data integrated between a deprojected distances from the star, where a significant detection would coincide with the stellar RV marked in grey. Top: covering a range between 60 and 240 au, and bottom: covering a range from 0 to 60 au.

Open in new tabDownload slide

We found the rms of the profiles (in the two integrated regions) across all velocities, and measured |$\rm {rms_{inner}{=}0.32\, \rm {mJy}}$| in the inner 0–60 au region (averaged in both clockwise and anticlockwise directions), and |$\rm {rms_{outer}{=}2.5\, \rm {mJy}}$| in the outer 60–240 au region (averaged in both clockwise and anticlockwise directions). We thus define a |$3 \sigma$| upper bound CO line peak flux (based on the 60–240 au region) of |$F_{\rm {CO,3 \sigma }}{=}7.5\, \rm {mJy}$|⁠. Although our channel spacing equates to a velocity width, |$\Delta v {=} 1.269\, \rm {km\, s}^{-1}$|⁠, the effective spectral bandwidth is ∼2.667 times larger than this (since adjacent channels in ALMA data are not fully independent from each other¹). Therefore, we place a |$3 \sigma$| limit on the integrated CO flux as |$\delta F_{\rm {CO},3 \sigma } {=} 2.667 \times \Delta v \times \rm {rms} {=} 25.4\, \rm {mJy\, kms^{-1}}$|⁠. Since |$\kappa \,$|CrB has more recently been interpreted as having a lower stellar mass (i.e. |$1.32\, {\rm M}_\odot$|⁠, see White et al. 2020, and also Fig. 1), we also checked the robustness of this CO flux limit by re-measuring this flux with a lower stellar mass of |$M_\star {=}1.32\, {\rm M}_\odot$|⁠, which also yielded a consistent non-detection. We discuss the CO J = 2–1 observations further in Section 7.3, where we use our CO flux upper limit to constrain the total CO gas mass, and the possible origins of CO.

6 THE κ CrB PLANETARY SYSTEM

Although there is a known RV planet, a low-eccentricity exo-Jupiter with an ∼1250 d orbit (see Johnson et al. 2008, and our introduction), long-term RV measurements of κ CrB over the space of 8 yr found this to have a gradient in its RV trend that may have been the result of acceleration from an outer massive body (i.e. a large planet or brown dwarf; see Bonsor et al. 2013). If this RV trend is indeed due to an outer body, to evade being detected (for a given orbital radius), it would need to be sufficiently faint (i.e. to have been undetected by the adaptive optics (AO) of Keck II observations; see Bonsor et al. 2013, for a complete discussion), and have sufficiently low mass (i.e. to avoid producing a significant proper motion anomaly in its Gaia eDR3–Hipparcos linear motion solution; see Brandt 2021, for a full discussion). Since the publication of Bonsor et al. (2013), the ∼8 yr of RV data (2004–2011) has expanded to ∼17 yr (2004–2021) with new APF observations (Automated Planet Finder telescope, at the Lick Observatory; see Radovan et al. 2014; Vogt et al. 2014) and HIRES (High Resolution Echelle Spectrometer) observations (the high-resolution echelle spectrometer, at the Keck Observatory; see Vogt et al. 1994). Given the doubled baseline of observations, we updated the constraints on the putative outer body using the radvel package (Fulton et al. 2018).

radvel allows users to define the orbital parameters for N bodies, for which we investigated the goodness of fit of models with N = 0, 1, or 2 planets. In Table 4, we report our model comparison results (assuming a stellar mass of |$M_\star =1.8{\pm }0.3\, {\rm M}_\odot$| and letting the jitter parameter vary freely) and start by noting that the N = 0 planet model is strongly ruled out, but that N= 1 and N= 2 planet models are more comparable. We note that the parameters of the inner planet in both N= 2 and N = 1 models remain strongly consistent with the previously confirmed exo-Jupiter planet in the literature (e.g. Johnson et al. 2008), and further, the inner planet of the N= 2 planet model has consistent values with the single derived planet in the N = 1 fit, as shown in Table 5. Whilst this model comparison investigation suggests that a two-planet configuration is statistically favoured, we argue against such a configuration on dynamical grounds. In the modelled scenario, whilst the inner planet remains in the same location, the model derives the outer planet mass as |$M_{\rm {pl\, c}}\sin {i}{=}0.39{\pm }0.10\, M_{\rm {Jup}}$|⁠, the semimajor axis as a = 6.0 ± 0.4 au, and the eccentricity as e = 0.78 ± 0.13. Such an eccentricity and orbital radius result in this body having a periastron at |${\sim }1.3\,$|au, i.e. inside the 2.8 au orbit of the exo-Jupiter, κ CrB b, which we deem an unstable configuration over the 2.5 Gyr age of the system. From this, we conclude that whilst the RV baseline has doubled from 8 to 17 yr, radvel is unable to determine the orbit of a stable N = 2 planetary configuration from the RV data. We show in Fig.  A1 of Appendix A the RV plots of the N = 2 parameter model.

7 DISCUSSION

Here, we will piece together the observational and modelling analyses presented earlier to tie up our understanding of the κ CrB planetary system as a whole, and place this in the context of other planetary systems. Thus far, our analysis of this 2.5 Gyr sub-giant can be summarized as follows:

(a) from the millimetre imaging and modelling, there is a single, 61° inclined, broad planetesimal belt with a PA of ∼146° (anticlockwise relative to north), with a peak emission radius of |${\sim }110\,$|au, which with an ALMA resolution of 63.2 × 42.1 au produces emission that extends over a wide radial region from the star of ∼50–180 au;

(b) from the scattered light imaging, there is significant dust emission that extends from a radius of at most 51 au out to 280 au (peaking at 73 au) with an estimated albedo of |$A{\gtrsim }10{{\ \rm per\ cent}}$|⁠, and although has a PA on the sky consistent with that measured at millimetre wavelengths, is significantly less inclined;

(c) from the RV monitoring, there is an inner planet on a 1248 d orbit, and a significant long-term RV trend, which if due to an outer body must have a mass between |$0.4\hbox{ and }120\, M_{\rm {Jup}}$|⁠, and an orbital radius between 8 and 66 au.

7.1 Understanding κ CrB’s planetesimal belt

7.1.1 The mass of κ CrB’s planetesimal belt

We start by assessing the mass of κ CrB’s planetesimal belt based on our millimetre analysis. We note that our models of the Band 6 emission both found a dust mass of |${\sim }0.016\, {\rm M}_\oplus$| with maximum model dust grain sizes of |${\sim }2\,$|cm, which can be used to constrain the size of planetesimals, and thus the total belt mass. For example, by setting the 2.5 Gyr system age equal to the collisional lifetime of the planetesimals in the belt (see equation 3 of Wyatt 2008), we set a lower limit on the maximum planetesimal size of bodies in the collisional cascade of |${\gtrsim }1\,$|km, and from this we can estimate a lower limit on the total belt mass of |${\gt }1\, {\rm M}_\oplus$| (assuming an α = 3.5 size distribution, as per Dohnanyi 1969). The above assumes a fixed value for |$Q^\star _{\rm {D}}$| of |${\sim }70\,$|J kg⁻¹, a relative velocity of |$660{\pm }130\,$|m s⁻¹ (using equation 10 of Matrà et al. 2019, with a fixed value of h = 0.05), a grain density of |$\rho =2.7\,$|g cm⁻³, and a width of 50 au, i.e. consistent with the lower-limit power-law model extent. This total belt mass and lower limit on planetesimal size is consistent with other debris discs at a similar age, like 1.4 Gyr q¹ Eri (Lovell et al. 2021), suggesting that although this disc is indeed long-lived, it is not an outlier in having retained detectable dust over ≳Gyr ages.

We next compare the modelled Band 6 dust mass of |$\kappa \,$|CrB with the values determined for the SONS (SCUBA-2 Observations of Nearby Stars) sample of main-sequence debris discs (Holland et al. 2017) in Fig. 11 (top panel). Whilst the SONS dust masses were determined at 850 |$\mu$|m with a dust opacity of |$\kappa _\nu {=}1.7\,$|cm² g⁻¹, to compare these with the dust mass derived for |$\kappa \,$|CrB’s debris disc at 1.3 mm, it is important to use a consistent dust opacity value; thus, we scaled the SONS dust masses by assuming an opacity power-law index of β = 1 (see e.g. Hildebrand 1983; Beckwith et al. 1990, which we note only alters the SONS dust masses by |${\sim }50{{\ \rm per\ cent}}$|⁠). We find that the dust mass of |$\kappa \,$|CrB’s debris disc is lower than most of the discs around less evolved A-type stars; however, this is unsurprising given that |$\kappa \,$|CrB is nearly an order of magnitude older than all SONS A-type stars, and has thus had longer to collisionally grind down. In fact, by comparing the dust mass of κ CrB with the mass evolution of Holland et al. (2017), i.e. the ∝ t^−0.5 trend in their fig. 34, it can be seen that the dust mass measured here is entirely in accordance with this prediction based on collisional evolution models (see e.g. Wyatt et al. 2007; Sibthorpe et al. 2018).

Figure 11.

Top panel: M_dust–L_⋆ for the SONS survey of debris discs (showing the complete data set in blue as presented in Holland et al. 2017), including the new data point for |$\kappa \,$|CrB (amber). Bottom panel: R_disc–L_⋆ relationship plot (showing the complete data set in blue as presented in Matrà et al. 2018, for resolved disc radii), including the new data point for |$\kappa \,$|CrB, along with its 50 au |$3\sigma$| extent (in red) and 135 au median extent (in amber).

Open in new tabDownload slide

Further, by comparing the disc peak emission radius, extent, and stellar luminosity of |$\kappa \,$|CrB with other main-sequence debris discs (using the data presented by Matrà et al. 2018), we show in Fig. 11 (bottom panel) that |$\kappa \,$|CrB’s debris belt peak radius remains strongly consistent with debris discs around main-sequence A stars. We present on this plot both the strict and median extents of the disc (i.e. the 50 and 135 au values demonstrated by the broad disc model SPL), demonstrating that both these measures of the disc extent are also consistent with main-sequence A stars.

Neither of the conclusions that |$\kappa \,$|CrB’s debris disc is consistent with the dust masses and disc radii of main-sequence A-type debris discs are particularly surprising given how early |$\kappa \,$|CrB is into its giant branch evolution (see e.g. Fig. 1 in which early stage sub-giants have comparable luminosities to main-sequence A-type stars) and that A-type stars are readily detected with debris discs across the span of their main-sequence ages (see e.g. Rieke et al. 2005; Su et al. 2006; Wyatt et al. 2007). Nevertheless, we have shown that for the total disc mass, total dust mass, and disc radius, that |$\kappa \,$|CrB’s disc remains consistent with the population of less evolved debris discs around A-type main-sequence stars in the context of collisional evolution, despite its older |${\sim }2.5\,$|Gyr age. Future observations of sub-giants/giant branch stars to search for the presence of debris disc dust may uncover whether or not the standard collisional evolution models remain applicable around stars that have evolved similarly to or beyond the stage of |$\kappa \,$|CrB. We further note that the levels of dust present are significantly higher than those in the Solar system’s Kuiper belt (see e.g. Wyatt 2008; Booth et al. 2009). From this, we highlight the same point made in Bonsor et al. (2013) that this strongly implies that there very likely has been no ‘Late Heavy Bombardment’ type event in the planetary system of |$\kappa \,$|CrB.

7.1.2 The inclination of κ CrB’s planetesimal belt

Unusual about the observations of |$\kappa \,$|CrB is the inclination discrepancy between (i) the millimetre Band 6 ALMA and far-infrared Herschel images (which both suggest this to be ∼60°–61°), and (ii) the scattered light image, which is better fitted with models inclined at ∼30°–40°. Assuming that there are no underlying issues with either data set, here we explore three hypotheses that may explain this inclination discrepancy.

Radial scattering anisotropy dependence? One possibility to explain the inclination discrepancy is that the scattering factor (that determines the level of dust emission anisotropy) has a radial dependence. This could have the effect that emission at increasingly greater distances from the star becomes preferentially scattered, which consequently reduces the assessment of the scattered light dust inclination when modelled with a fixed scattering factor (as was conducted in Section 4.2). This effect is likely to be most pronounced in broad discs where there are a wide range of radial distances over which differences in scattering factors must be accounted for. Indeed, based on the best-fitting models using single-valued HG functions presented in Section 4.2, to model |$\kappa \,$|CrB’s scattered light emission with a single phase function, such radially dependent scattering factors are likely required (given the large residuals that remain after model subtraction). The problems associated with modelling scattered light observations of debris discs with single-valued HG parametrizations are well known (see e.g. Stark et al. 2014; Milli et al. 2017). Whilst in addition to the HG models, we inspected models based on Mie theory, it may still be the case that the inclination discrepancy suggests that a better model is required, for example, with a radially varying phase scattering function. We note that confirming whether or not this can fully account for the measured inclination discrepancy requires additional modelling work beyond the scope of this study with more complexity than the simple scattered light models presented here.

A dust size-dependent scale height? Another possibility is that |$\kappa \,$|CrB’s micrometre-sized dust is being vertically ‘puffed up’ more than the dust observed at millimetre/far-infrared wavelengths, which could have the effect of lowering the estimated inclination at sub-micrometre wavelengths. Scale heights in debris discs have a number of origins, such as being from their primordial planetesimal distribution or from planetesimal scattering by large embedded bodies (for a discussion on this in the case of the edge-on disc, AU Mic; see Daley et al. 2019). The combination of radiation pressure and collisional processes can additionally induce larger scale heights in smaller grains than larger ones. For example, since smaller grains are blown on to eccentric orbits by radiation pressure, these have higher impact velocities relative to larger grains, which due to momentum conservation result in these receiving larger impact ‘kicks’ during collisions with larger grains/planetesimals. The result of this was modelled by Thébault (2009), in which small grains close to the blowout size were shown to be more vertically puffed up than the largest grains, i.e. to average inclinations of ∼8° versus just ∼0.5° for the smallest and largest test particles, respectively. Therefore, whilst this effect may not fully account for the ∼20°–30° difference in the modelled scattered light inclination and that observed in the far-infrared/millimetre, this may contribute to the discrepancy. We therefore recommend further collisional modelling of the |$\kappa \,$|CrB system to constrain the extent to which collisions and radiation pressure can induce dust size-dependent scale heights in this debris disc.

We note in addition that whilst radiation pressure may not be responsible for the inclination discrepancy, radiation pressure does naturally explain how the observed scattered light halo extends out to ∼280 au, whereas the millimetre dust is not observed significantly beyond ∼180 au. In other words, whilst radiation pressure may not dominate the vertical distribution of the micrometre-sized dust, it does have a significant effect on the radial structure.

Misaligned precessing belts? The previous two hypotheses share the feature that their dust belts originate from the same distribution. Here, we suggest the possibility that these discs are instead part of different belt distributions, i.e. that the millimetre/far-infrared and sub-micrometre belts have physically different inclinations, with the parent belt (traced in the millimetre) precessing at a rate different from that of the smaller dust in scattered light halo (in a mechanism similar to that discussed in Sende & Löhne 2019). Whilst this might be at odds with the perceived alignment in the PA of the dust observed with ALMA, Herschel, and HST (which we find are all within the range of 142°–152°), for example, since any precession would be likely to misalign these belts in both inclination and PA, this could be a chance encounter where these distributions are aligned in PA, but not in inclination. On the other hand, higher resolution imaging and more detailed modelling may indeed find that the PA are less well aligned than our existing modelling procedure determined. Ultimately, there remain uncertainties with all three of these scenarios, and we briefly outline further work necessary to understand the nature of this inclination discrepancy.

In our analysis of the scattered light emission, the greatest uncertainty that remains is the location of the inner edge of the dust. Since this lies within the HST imaging artefacts, the scattered light dust inclination derived in this work is model dependent: A direct measurement of the belt inner edge with scattered light imaging could however determine this uniquely. We therefore suggest new scattered light observations of |$\kappa \,$|CrB, e.g. with SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch instrument) or GPI (Gemini Planet Imager), to resolve and directly measure this inner edge, derive more accurately the scattered light inclination, and re-evaluate the hypotheses outlined here.

7.1.3 Asymmetries in |$\kappa \,$|CrB’s debris disc?

In our modelling section, we showed that the millimetre-emission residual maps (Fig. 7) all contained the same |$30\!-\!40\, \mu$|Jy clump to the south of the belt (at a relative RA, Dec. location of −1.0, −5.0). Perhaps simply, given that all of the residual maps likewise contain negative residual emission at the same level (e.g. see RA, Dec. −5.0, −2.5), we cannot rule this out as being due to noise. On the other hand, this emission could be real, and we explore three hypotheses, namely that this is emission from a background millimetre galaxy (MMG) or dust, either from a planetesimal collision or in a circumplanetary orbit, and suggest future work that could distinguish between these scenarios.

A background galaxy? We first assess the more likely of the three options, that this is emission from an MMG. For example, based on the galaxy counts of Simpson et al. (2015), we find that the probability that an MMG with ≥0.035 |$\mu$|Jy falls within a 60° inclined ellipse with a semimajor axis of 5.0 arcsec exceeds order unity; in other words, it is not surprising that at the depth of our observations a significant clump of emission is coincident with the debris disc. Such a hypothesis can be tested since |$\kappa \,$|CrB is a high-proper motion star, whereas the location of an MMG is fixed on the sky. However, given the proper motion of |$\kappa \,$|CrB is μ_{RA, Dec} = (−8.79 ± 0.18, −347.77 ± 0.20) |$\rm {mas}\, \rm {yr}^{-1}$| (see Gaia Collaboration et al. 2018) and the resolution of our imaging (∼2.1 × 1.4 arcsec), this hypothesis can only be tested at the 3|$\sigma$| level 18 yr on from these 2019 observations, in 2037 December.

A planetesimal collision? Secondly, we explore the interpretation that this clump has a planetesimal origin and estimate the size of the body required to produce this, by first inferring the dust mass required to produce 35 |$\mu$|Jy emission as |${\sim }4\times 10^{-4}\, {\rm M}_\oplus$|⁠, i.e. using the same ratio of the total disc flux to the total dust mass. To arise from the complete destruction of a single planetesimal, such a body would need a diameter D|${\gtrsim }1200\,$|km, assuming (1) a grain density of |$\rho =2.7\,$|g cm⁻³ and (2) that the planetesimal breaks up only into the millimetre–centimetre-sized grains that show up as a clump in the emission. Although this is consistent with the lower-limit maximum planetesimal size estimated earlier, this is a factor of ∼1000 larger. We note that this planetesimal size is a lower limit since the catastrophic break-up of a planetesimal is likely to produce a cascade of different-sized particles, whereas we have only considered the millimetre–centimetre-sized grains. Without repeating the analyses of Krivov & Wyatt (2021) and Lovell et al. (2021), we simply re-state here that if such D|${\gtrsim }1200\,$|km-sized bodies are part of a collisional cascade with a Dohnanyi (1969) size distribution, this would suggest a total debris disc mass that is substantially larger than the |${\sim }1\, {\rm M}_\oplus$| estimated earlier, for example, with a mass of many tens of M_⊕, based on |$M_{\rm {disc}} \propto D^{1/2}_{\rm {max}}$|⁠. Since a more realistic collision at this location would produce significant levels of micrometre-sized dust as part of a cascade that would be blown out (a possibility discussed by Grigorieva, Artymowicz & Thébault 2007), this could provide a simultaneous explanation to the inclination discrepancy discussed in Section 7.1.2. For example, whilst the millimetre-sized grains would spread around the orbit of the parent planetesimal, the micrometre-sized grains would be preferentially blown out, and thus cause the observed NE–SW asymmetry seen in the scattered light emission.

Massive circumplanetary rings? Finally, we consider the possibility that this is unresolved circumplanetary dust orbiting an outer body with the same dust mass as inferred from the planetesimal origin hypothesis, i.e. |${\sim }4\times 10^{-4}\, {\rm M}_\oplus$|⁠. This level of mass is just ∼1/30 that present in the moon, or 100–200 times that present in Saturn’s rings, and consequently, would not require a particularly massive body to retain this level of dust. Our modelled limits on the size of the millimetre debris disc can provide additional constraints on the mass of such a body based on the outer edge debris disc dynamical limits, as plotted in Fig. 12 and which we discuss in detail in Section 7.2. For now we note that the limits set by the debris disc model show that the mass of any planet external to the planetesimal belt must be |${\lesssim }3\, M_{\rm {Jup}}$| at the 205 au distance of the clump to the star (assuming this orbits in the same plane as the belt). In other words, whilst it may not be necessary to explain the outer edge of this debris disc via planet sculpting, if this clump is due to circumplanetary dust around a massive outer planet (which we emphasize would be a different body from the one causing the RV trend), then it may additionally constrain the disc outer edge location. Better constraints on the disc outer edge could therefore more tightly constrain the outer planetary system region.

Figure 12.

Shown here are the possible places that planets can/cannot reside (unhatched/hatched regions, respectively) within the |$\kappa \,$|CrB system, based on all the constraints set out in this work, in mass–semimajor axis space (for an assumed inclination of 61°) and the location of the known planet, |$\kappa \,$|CrB b. In green is the area that a companion must exist in order to produce the observed RV trend, and in amber is the area that planets can exist without affecting any observations. In grey are the areas that planets cannot exist, as ruled out by Bonsor et al. (2013), i.e. due to both RV and AO limits. In red are the areas that planets cannot exist based on this work, ruled out by dynamical arguments, i.e. the chaotic zones of |$\kappa \,$|CrB b, and putative planets on the inner and outer edges of the debris disc.

Open in new tabDownload slide

Since the planetesimal collisions and circumplanetary dust hypotheses have a common origin in the emission being due to dust, these are difficult to distinguish, even with long-term, high-resolution millimetre monitoring. Although dust formed in a collision would gradually spread around its orbit over time, whereas dust in a circumplanetary ring system would remain coincident with its host planet, we should expect the shape of the millimetre emission of these two scenarios to diverge over time. However, at ∼205 au, orbital time-scales are very long and it may take many hundreds to thousands of complete orbits for dust formed in a planetesimal collision to spread sufficiently to be measurable (see Jackson et al. 2014, for a discussion on collisional spreading time-scales). Therefore, since it could take many thousands of years to accurately distinguish between these origins in the millimetre, detailed analysis of this emission in other wavebands is likely to be necessary.

7.2 Planet–disc interactions in κ CrB’s planetary system

Having discussed the planetesimal belt in isolation in Section 7.1, here we extend our analysis of the wider system by considering the impact that the plausible planet architecture has on the planetesimal belt structure by considering possible planet–disc interactions (as studied extensively over recent years for the case of axisymmetric and asymmetric discs; see e.g. Wyatt 2005; Mustill & Wyatt 2009; Shannon et al. 2016; Lynch & Lovell 2021).

We start by introducing Fig. 12, which shows the parameter space within which any planets/companions in this system can reside, assuming that the RV trend is caused by a massive outer companion. We define three categories, namely where an outer body must reside (in green), can reside (in amber), and cannot reside [in grey or red, depending on these constraints being based on Bonsor et al. (2013) or this work, respectively]. The grey regions are set by (i) the previous Keck II AO limits and (ii) the RV data (using the 2004–2021 baseline), which strictly constrains the types of planets that can reside within 8 au (see Section 6). The amber region is set by equation (1), above (below) which planets would (would not) induce an RV acceleration trend. The red limits are set by three dynamical constraints. First, that other planets cannot reside within the orbital instability region set by |$\kappa \,$|CrB b [i.e. within its chaotic zone, where we note the inner and external chaotic zones of a planet as Δa_cz = ± 1.5a_pl(M_pl/M_⋆)^0.28 (see e.g. Quillen & Faber 2006), where a_pl is the planet semimajor and M_pl is the planet mass]. Secondly, that planets cannot be present in the parameter space within which these would cause an observable proper motion anomaly between the Gaia DR2 and eDR3 epochs (‘PMA lim’). Thirdly, that planets cannot be present where these would have chaotic zones that would intersect with the inner and outer edges of the debris disc. The green region is then bounded by all such conditions, and defines the parameters of a body that would cause the RV trend. This has two important implications. First, that if the RV trend is due to an outer body, this must be massive, with |$0.4\!-\!120\, M_{\rm {Jup}}$|⁠, and lie between |$8\hbox{ and }66\,$|au, e.g. this could be consistent with a range of bodies, such as an eccentric exo-Saturn or a brown dwarf at tens of au. Secondly, that additional planets can be present without inducing an RV acceleration trend in the RV data if these reside below the RV trend limit, i.e. in the amber regions, meaning that inner terrestrial planets, and even giant planets at tens of au, can stably reside in this system without being detected. Nevertheless, we note that if the RV trend is not due to an outer body, the constraints on where planets may reside differ from Fig. 12. For example, although the red and grey AO limit regions would remain, any regions bounded by the RV data would be unconstrained.

7.2.1 Dynamical constraints: more than two planets?

Since scattered light sub-micrometre emission does not trace the parent planetesimal belt as accurately as millimetre emission, the constraints we place on the planetary architecture are derived from the Band 6 ALMA data. The millimetre modelling in Section 4 showed that whilst 1.3 mm emission begins at |${\sim }50\,$|au, this is consistent with a broad belt between 93 and 151 au (though the outer edge is a lower limit, e.g. see discussion of model SPL in Section 4.1.2). The presence and morphology of this disc provide constraints on the planetary architecture in three ways. First, the presence of a large cavity between the star and the disc inner edge suggests that this was caused by planet sculpting from unseen inner planets. Secondly, since the chaotic zones of planets cannot overlap with the location of the planetesimal belt, there is likely a dearth of planets (above a given mass) in this region. Thirdly, for the disc to be producing visible levels of dust, it is possible that this disc is being stirred by a planet. There is however only one known planet, the inner exo-Jupiter |$\kappa \,$|CrB b, and this planet alone could not be responsible for the observed disc architecture. For example, we rule out the possibility that belt stirring can be achieved by this inner exo-Jupiter since the secular stirring time-scale (see equation 15 of Mustill & Wyatt 2009) is ≳ 1000 times the age of this system. The combination of these three conditions (in addition to the known RV acceleration trend, see Section 6) thus argues in favour of the presence of additional outer planets, which we discuss further below.

We first consider a scenario in which planets underwent no significant migration, and note first that the limits set by the debris disc stability in Fig. 12 are based on the requirement that planetary chaotic zones (i.e. Δa_cz ± 1.5a_pl(M_pl/M_⋆)^0.28; see Quillen & Faber 2006) cannot overlap with the planetesimal belt that spans from 93 to 151 au. The disc inner edge thus sets two constraints on the type of companion/s that may be responsible for causing the observed cavity via (i) sculpting the belt inner edge at 93 au and (ii) depleting planetesimals from a few au outwards to 93 au. Assuming that planetesimals originally formed at all distances from the star, in a scenario where no planets migrated, we find that in order to explain the RV trend, the extent of the cavity, and the disc inner edge location, multiple additional planets within the inner cavity are required. For example, the widest orbiting companion within the cavity that could carve the inner edge (i.e. a massive body at 66 au, see Fig. 12) can only clear planetesimals from 46 au out to 93 au, i.e. such a companion cannot clear planetesimals sufficiently far inwards and thus additional planets would be required to do so in the inner tens of au. Indeed, there are no scenarios in which a single companion in addition to |$\kappa \,$|CrB b can carve the morphology of the debris disc (without requiring any planetary migration, which we discuss further in Section 7.2.2).

To constrain the number and mass of any such planets at tens of au responsible for carving the cavity, we consider two limiting scenarios based on the putative RV companion residing at either end of its constrained location. If at 8 or 66 au, the putative RV companion would have chaotic zones either between 5 and 13 au or between 46 and 93 au, respectively. By using equations (4) and (5) of Shannon et al. (2016) and the system age of 2.5 Gyr, we can estimate possible planetary architectures by assuming that these planets would be responsible for the disc cavity either between 5 and 46 au (i.e. in the case that there is a massive companion at 66 au) or between 13 and 93 au (i.e. in the case that there is a massive companion at 8 au). Respectively these were found to either require 13 planets each with a mass >0.7 M_⊕ or 8 planets each with a mass >2 M_⊕. In either of these two scenarios, all of these bodies would be in the amber region of Fig. 12, and the putative RV companion in the green region. We note here, however, that since the cavity may have formed much earlier in the system’s lifetime [and the equations of Shannon et al. (2016) are based on the age of the system], fewer planets with a greater total mass could have been responsible for the observed cavity, i.e. the above planet mass estimates should be seen as lower limits. Since the putative RV companion could however be anywhere between 8 and 66 au (rather than at either limit), these two example calculations serve to demonstrate that in addition to κ CrB b and the putative RV planet, a significant number of additional planets may be required to produce the observed disc morphology if planet sculpting from non-migrating planets is the main cause of the observed cavity.

7.2.2 Planet migration: just two planets?

An alternative hypothesis is that this system formed with two planets: the inner eccentric exo-Jupiter, κ CrB b, and the outer RV companion that underwent significant migration. For example, this outer body could have formed in the inner regions of this system, and migrated outwards to 66 au following e.g. a planet–planet interaction with κ CrB b, or another dynamical instability, clearing the full cavity and carving the 93 au inner belt edge. Indeed, it has been shown that planet migration over many tens of au is possible in systems with N ≥ 2 planets in which one of these migrates out through a disc (see e.g. Martin et al. 2007). As such, a body migrating through the belt would interact significantly with the belt planetesimals, either accreting, ejecting, or trapping these in resonances. Potentially problematic for such a scenario is that planetesimals trapped in resonances could induce significant sub-structure in the debris disc (e.g. such as clumps; see Wyatt 2003, depending on parameters such as the migration rate and planet mass), which we do not observe. Although at the depth of our millimetre observations we see no evidence for any significant asymmetries in the emission co-located with the belt that might indicate early planet migration, future deep observations may uncover sub-structure related to such a dynamical history. On the other hand, if such a migration happened far earlier in |$\kappa \,$|CrB’s lifetime, the population of planetesimals in the resonant zones may have depleted sufficiently that such asymmetric imprints no longer remain observable, i.e. since such observable asymmetries can be smoothed over. For example, in systems where planetesimals start with moderate eccentricities, asymmetries become less pronounced over time (see e.g. Reche et al. 2008). Thus, deeper images of |$\kappa \,$|CrB’s debris discs may still be insufficient to exclude this type of migration history.

7.2.3 What will become of this debris disc?

The fate of |$\kappa \,$|CrB’s debris disc as the star ascends the giant branch is intimately bound with the location of its planets. Bonsor & Wyatt (2010) showed that as a result of post-main-sequence stellar evolution and the adiabatic stellar mass-loss of A-type stars at the end of the giant branch stage, the enlargement of planet and planetesimal orbits (which can expand by the same fraction of stellar mass lost) has implications for the depletion of debris discs. For example, the fractional increases in the orbital distances are of the order of 2–3 towards A-type stars, thus further increasing the size of the chaotic zones of planets and the opportunity for resonance-overlap instabilities, but also the radii of the planetesimal orbits. This means that if a massive companion is responsible for the sculpting of the inner edge of this debris disc, the disc may be significantly disrupted by such a body as the star loses mass towards the end of the giant branch stage. Such disruptions present significant opportunities for the inward scattering of planetesimals and thus their break-up/sublimation near the stellar surface; thus, the |$\kappa \,$|CrB system could very well be one of the precursors to the population of bright ‘polluted’ white dwarfs seen at much later evolutionary epochs (see Zuckerman et al. 2003, 2010; Koester, Gänsicke & Farihi 2014; Veras et al. 2014). On the other hand, the presence of the inner planet |$\kappa \,$|CrB b (which at such an evolved epoch may instead reside at ∼6–8 au) could efficiently eject such inwardly scattered planetesimals, and thus the flux of planetesimals on star-grazing orbits may be consequently very small. New modelling based on the constrained planetary system architecture of this work could therefore explore the level of pollution expected at the white dwarf stage.

7.2.4 How does κ CrB’s planetary architecture compare to others?

Planet formation in the κ CrB system must have been reasonably efficient. Even just considering the total planet mass of the confirmed exo-Jupiter κ CrB b, more than twice the amount of material formed planets in the young protoplanetary disc of κ CrB than did so in the Solar system (though we note that κ CrB may be nearly twice as massive as the Sun). This however could be as high as ∼120 times that of the Solar system, e.g. if the RV acceleration trend is due to a planet at the upper mass limit of Fig. 12. Whilst compared with our G-type Sun, this suggests that planet formation may have been substantially more efficient around κ CrB, this may instead reflect the larger average mass of material available to form planets within protoplanetary discs that is available around intermediate-mass stars (see e.g. Ansdell et al. 2016; Pascucci et al. 2016). Alternatively, planet formation efficiency around κ CrB may also be substantially lower than other intermediate-mass stars. For example, around the 30 Myr A-type star, HR 8799, at least four giant planets formed with a total mass of ∼30 M_Jup (Marois et al. 2010). If the putative RV companion is towards the lower end of the constrained mass bounds of Fig. 12, then planet formation may have been over an order of magnitude more efficient around HR 8799 than κ CrB. Nevertheless, with only loose constraints on the total planet mass in the κ CrB system, there remains a large uncertainty on whether the planet formation efficiency is indeed discrepant with known planetary systems around stars that could have formed with comparable levels of primordial material.

It is further instructive to compare the location of planets around κ CrB with other systems. We first do so with reference to the snowline. For the simple assumption that dust within planetary systems behaves as a blackbody, we can estimate the location of the snowline from |$r_{\rm {snowline}}=(278.3/T_{\rm {snow}})^2L^{0.5}_\star$|⁠. For a sublimation temperature of T_snow = 170 K and κ CrB’s luminosity |$L_\star =12\, L_\odot$|⁠, this equation would place the snowline at 9.3 au, though we note that earlier in the system’s history when planets were forming, this will have been closer to 6 au. This would mean that the known planet κ CrB b is inside of this, with the putative outer RV companion outside of this. This means that irrespective of the mass of planets formed around κ CrB, the broad architecture of this is different from that of HR 8799 in which all four giant planets orbit with radii significantly beyond their snowlines. Although this may indicate that κ CrB and HR 8799 have distinct planet formation histories, the observed architectural differences may also be due to the young age of the HR 8799 system in which further orbital architecture evolution will occur, for example, via planet migration (e.g. through planet–planetesimal belt interactions) or via planet–planet scattering events. In contrast, whilst the mass levels of the bodies in the Solar system and κ CrB are at least a factor of 2 different, the location and type of bodies may be entirely consistent. For example, Fig. 12 demonstrates that planets consistent with an exo-Earth (lower left amber region), an exo-Jupiter (κ CrB b), an exo-Saturn (the RV trend companion at the lower left of the green region), and multiple exo-Uranus/exo-Neptune analogues (in the middle amber region) are all consistent with the RV data, and need only a few more Earth-mass planets tens of au beyond those to carve the disc cavity and inner edge. In other words, the planetary architecture of κ CrB may be a slightly broader and more massive version of the Solar system, albeit with a significantly more massive planetesimal belt. As such, with tighter constraints on the planetary architecture of κ CrB, Solar system formation models may provide significant insight into how this planetary system formed and evolved.

7.3 CO gas: constraints on giant branch sublimation?

Having fully explored the planetary system dust observations and RV planet modelling, we attempt to constrain the planetary system composition. In Section 5, by analysing the emission near the J = 2–1 transition line frequency, we demonstrated that if any CO gas is present in the |$\kappa \,$|CrB planetesimal belt (or closer to the star), it must be below our detection limits. However, the upper-limit flux (⁠|$25.4\, \rm {mJy\, kms^{-1}}$|⁠) can be used to determine an upper limit on the total mass of CO gas present based on the gas excitation conditions (i.e. as set by the kinetic temperature, density of electrons, and the radiation environment; see Matrà et al. 2015, for further details). For the CO J = 2 rotational level and based on the system being in non-local thermodynamic equilibrium, we used the code of Matrà et al. (2015, 2018) including fluorescence to estimate the CO gas mass at four temperatures (i.e. 10, 20, 50, and 100 K, consistent with the expected temperatures of outer, cool debris discs) for the disc with a peak emission radius at |${\sim }110\,$|au, and the star with a stellar flux consistent with the SED (i.e. see Section 3.3). Fig. 13 shows the outcome of these calculations for each of the four temperatures over a range of electron densities. From this plot, we find an upper limit on the CO gas mass in the range |$(4.2\!-\!13)\times 10^{-7}\, {\rm M}_\oplus$| (i.e. based on the 100 K minima and maxima, which correspond to moderate and high electron collision partner densities). We note that electron densities range from low values (i.e. |${\sim }10^{-4}\,$|cm⁻³) where excitation is radiation-dominated through to high values (i.e. |${\sim }10^{10}\,$|cm⁻³) where excitation is collision-dominated.

Figure 13.

Interpreted limits on the CO gas mass (based on the upper-limit CO flux) plotted as a function of the unknown density of electron collision partners, for four assumed temperatures (10, 20, 50, and 100 K, in blue/solid, amber/dash–dotted, green/dashed, and purple/dotted, respectively).

Open in new tabDownload slide

From this we can estimate how much CO may be present in the solid planetesimals via two scenarios: (i) that all CO is released via collisions or (ii) that all CO is released via sublimation. Considering the first of these, we apply equation (2) of Matrà et al. (2017b) to bound the |$\rm {CO}$| and |$\rm {CO_{2}}$| ice fraction of planetesimals. We estimate this limit between |$f_{\rm {CO+CO_{2}}}{\lt }85\!-\!95{{\ \rm per\ cent}}$|⁠, based on the 110 au belt radius, lower limit on the disc extent of 50 au, fractional and stellar luminosities consistent with Section 3.3, a stellar mass of |$1.8\, {\rm M}_\odot$|⁠, a CO photodissociation time-scale of 120 yr (Visser, van Dishoeck & Black 2009), and the range of CO gas masses determined above. We note that although this range of ice fractions sets only loose constraints on the composition of planetesimals, it is strongly dependent on both the belt extent and photodissociation time-scale, for which a belt with a narrower extent and well-shielded CO (and thus longer CO photodissociation time-scales) would reduce |$f_{\rm {CO+CO_{2}}}$| substantially. Nevertheless, the range of ice fractions determined here remains consistent with Solar system comets, readily found in the range of |$f_{\rm {CO+CO_{2}}}{\sim }5\!-\!20{{\ \rm per\ cent}}$| (Le Roy et al. 2015).

An alternative scenario could be that any CO that is present in the |$\kappa \,$|CrB debris disc is released via the sublimation of CO-ice, as the star ascends the giant branch and planetesimal and grain surfaces are warmed up. However, based on equation (15) of Jura (2004), we find at a temperature of 70 K that the total area of dust in the |$\kappa \,$|CrB debris disc is a factor of 4 × 10¹⁰ too low to liberate detectable levels of CO, i.e. consistent with our non-detection. In estimating this value, we assume a number of parameters. First, that the belt has a temperature of 70 K based on the blackbody temperature of dust at the peak disc radius of 110 au around a |$12\, L_{\odot }$| star. Secondly, that the total cross-sectional area of dust present (see equation 4 of Wyatt 2008) is ∼7.6 au², for which we use the millimetre disc radius, |$r{=}110\,$|au, and disc fractional luminosity, f = 5 × 10⁻⁵. Thirdly, that CO survives over a 120 yr photodissociation time-scale such that there is no shielding of CO (see e.g. Visser et al. 2009). Fourthly, that the mass production rate is consistent with Kuiper belt object sublimation model predictions, i.e. with a mass-loss rate per unit surface area as a function of temperature, |$\dot{\sigma }_0{=}3.8{\times } 10^{8}$| g cm⁻² s⁻¹ K^−1/2 (see Jura 2004). Finally, that to be observable, a CO gas mass of |${\sim }1{\times }10^{-6}\, {\rm M}_{\oplus }$| would be required with a CO-ice fraction of 10 per cent, and that all liberated CO begins in an ice phase. We therefore note that given the very low levels of CO that can be liberated via this sublimation mechanism around |$\kappa \,$|CrB, any CO that is present is likely to be dominated by that produced via collisions, as is the case for main-sequence debris discs.

However, the CO-ice sublimation temperature dependence is sharp, e.g. temperatures beyond 100 K can liberate CO at a much faster rate than is possible at 70 K. The evolutionary time-scale for ascending the giant branch can take tens to hundreds of Myr, and so a period of sublimation-dominated CO gas production around |$\kappa \,$|CrB has not yet been reached; however, more luminous stars may have already reached this stage. We note that for the average debris disc radii of AFGK-type stars this condition is only met when L_⋆ is in the hundreds of L_⊙, and so is not possible at the sub-giant stage. Given the relationship between temperature and sublimation rate, and between the reduction in planetesimal size and temperature (see equations 15 and 21 of Jura 2004, respectively), such high disc temperatures increase the contribution to the total CO available via sublimation from larger planetesimals. Whereas at temperatures ≲100 K, the total CO mass liberated via sublimation is dominated by the contribution from approximately micrometre-sized grains. As such, on the giant branch, there may arise a transition in the dominant mechanism by which CO is released, from collisions to sublimation, and towards a CO production rate sufficiently high that CO becomes readily observable. It may thus be the case that whilst CO could be rare to detect at the sub-giant stage, studies of giant stars detect CO gas in abundance. We note that very few observations of debris discs around giant stars have been made, and thus disc evolution during this stage should be probed further. We therefore suggest new observations and further detailed modelling work to better constrain CO production during giant branch ascent, accounting for the CO release contribution from both collisions and sublimation, the effect of CO photodissocation, and for debris discs with varying radii and widths.

8 CONCLUSIONS

New images and analysis of the sub-giant star |$\kappa \,$|CrB have been presented at 1.3 mm with ALMA and at 0.6|$\, \mu$|m with HST, which demonstrate the existence of a debris disc with peak emission radius of |$110\,$|au with an extent between (at least) 50 and 180 au in the millimetre, and between (at least) 51 and 280 au in scattered light. Given |$\kappa \,$|CrB has now evolved off the main sequence, this is therefore the most evolved debris disc resolved with ALMA to date, and demonstrates the feasibility of imaging post-main-sequence debris discs.

By modelling the Band 6 ALMA data, we have shown that the continuum emission is consistent with being due to a single broad planetesimal belt with an extent |${\gt }50\,$|au and inner edge starting at 93 au. We show that the debris disc of |$\kappa \,$|CrB is consistent with the population of discs around main-sequence A-type stars (which this star evolved from) and has a mass consistent with expectations from collisional evolution models, i.e. sub-giant branch evolution may not yet have significantly influenced its debris disc properties.

We presented an updated analysis of the stellar RV over a baseline of 17 yr and re-confirmed the presence of an inner eccentric exo-Jupiter at ∼2.8 au, and an acceleration trend in the residual RV data. If the RV acceleration trend is due to an outer body, we found that it must have a semimajor axis between 8 and 66 au and a mass between |$0.4\hbox{ and }120\, M_{\rm {Jup}}$|⁠, i.e. this would need to be at least a wide-orbiting giant planet. As such, combined with the known exo-Jupiter, this suggests that planet formation towards |$\kappa \,$|CrB was efficient, forming two or more giant planets.

Further, we used the constraints set by the RV data and the debris disc in the millimetre to consider the architecture of the complete planetary system, showing that the cavity between the star and disc can be explained by a single massive outer companion and a string of lower mass bodies. This analysis demonstrates the unique constraints disc observations can place on planets.

We extended our single broad belt model to investigate the HST scattered light emission of |$\kappa \,$|CrB and found this is best fitted with highly anisotropic dust (with HG factor g = 0.8–0.9) but with a significantly lower inclination (i.e. 30°–40°) than the millimetre (i.e. ∼60°), despite the HST and ALMA images showing broadly consistent PA and emission scale sizes. We explored three possibilities: that our scattered light model is insufficiently detailed, that the sub-micrometre dust has a much greater vertical scale height than the millimetre dust, and that given the old age of this system, these are due to two misaligned precessing distributions.

Since there is a clump in the south of the disc in the millimetre images, we explored a number of scenarios as to the origin of this. Although it seems most plausible that this is either noise or a background MMG, we discussed the potential that this is also a clump of dust, and by constraining its mass, discuss the implications of this being either due to a planetesimal collision in the belt or dust in a circumplanetary ring system.

By further analysing emission near the CO J = 2–1 spectral line, we also report no evidence of CO gas. With our uncertainties, we constrain the total CO gas mass as |$M_{\rm {CO}}{\lt }(4.2\!-\!13) \times 10^{-7}\, {\rm M}_\oplus$|⁠. We find this is consistent with levels expected from planetesimal collisions (for planetesimals with Solar system comet-like compositions), or those associated with the sublimation of pure CO-ice as |$\kappa \,$|CrB begins to ascend the giant branch.

Deeper observations in the future will allow tighter constraints to be placed on the morphology of the belt at sub-micrometre and millimetre wavelengths, search deeper for any gas that may be present as this system begins to warm up as the star evolves off the main sequence, and tighten the constraints placed here on the wider planetary system.

SUPPORTING INFORMATION

RVData_HD142091.csv

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their insightful comments that greatly improved the clarity of this manuscript. We also thank Luca Matrà for providing code with which CO gas mass limit estimates were derived. We gratefully acknowledge the support and expertise provided by ALMA and HST staff involved in the data collection and post-processing/quality assurance undertaken in advance of our analysis. JBL is supported by an STFC postgraduate studentship. GMK is supported by the Royal Society as a Royal Society University Research Fellow. SM is supported by a Research Fellowship from Jesus College, Cambridge. AB is grateful to the Royal Society for a Dorothy Hodgkin Fellowship. PK thanks support from GO-13362 provided by NASA through a grant from STScI under NASA contract NAS5-26555.

DATA AVAILABILITY

This work makes use of the following ALMA data: |$\rm ADS/\, JAO.ALMA\, \#2019.1.01443.T$|⁠. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The HST data used for this study are part of the Cycle 21 with proposal ID GO-13362 and are publicly available on the Barbara A. Mikulski Archive for Space Telescopes (https://archive.stsci.edu/hst/search.php). This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France. This research has made use of NASA’s Astrophysics Data System. All RV data used to model the planet fits have been included as a supporting file alongside this paper.

Footnotes

REFERENCES

APPENDIX A: FITTED PLANET PARAMETERS

Here, we provide the plots of best-fitting N = 2 model derived by radvel in Fig. A1.

Figure A1.

radvel predictions of the N = 2 planet model, based on all available RV data over 17 yr.

Open in new tabDownload slide

APPENDIX B: PROPER MOTION ANOMALIES: THE LONG-PERIOD LIMIT

In this appendix, we develop a simple analytic approximation for the expected proper motion anomaly for a long-period binary on a circular orbit. We relate this to observations of the same system in the Gaia DR2 and eDR3 catalogues (same start date, with a baseline of 22 and 34 months, respectively).

B1 Long-period assumption

B2 Proper motion anomaly