© ESO 2021

1 Introduction

β Pictoris (also known as β Pic, HD 39060, HIP 27321, or HR 2020) is an A6V-type star hosting both a disk and a planetary system. It is one of the most frequently studied such systems in the past several decades. The debris disk around the star was originally detected in the form of infrared excess by the IRAS satellite (Aumann 1985), and shortly thereafter became the first circumstellar disk to be spatially resolved (Smith & Terrile 1984). In addition to being nearby at 19.44–19.75 pc (van Leeuwen 2007; Brown et al. 2018) and bright with a magnitude of V = 3.86 mag of the primary star, the β Pic disk is seen largely edge-on, with high resulting column densities along the line of sight, all of which benefits detectability relative to other disk host systems. On the celestial sphere, the disk extends between the north-east (NE) and south-west (SW) directions, and we correspondingly refer to the two arms of the disk as the NE and the SW arms.

The disk has been studied intensively both in terms of its dust (e.g. Weinberger et al. 2003; Liseau et al. 2003; Okamoto et al. 2004) and in terms of its gas (e.g. Slettebak 1975; Dent et al. 2014; Cavallius et al. 2019). The dust reveals fundamental morphological features and population substructures within the disk. For example, the NE arm of the disk is known to be substantially brighter in dust continuum than the SW arm (Kalas & Jewitt 1995), which might imply eccentric components within the disk. Moreover, a warp has been observed in the disk (Heap et al. 2000) that could be characterised as a secondary disk component superimposed on the primary one (Golimowski et al. 2006). This could be interpreted as revealing the existence of a giant planet on a slightly inclined orbit with respect to the main disk plane, dynamically inducing the secondary disk (Mouillet et al. 1997). In 2009 and later, this was confirmed when planet β Pic b was discovered (Lagrange et al. 2009, 2010), which was found to be on an orbit compatible with creating the warp (Chauvin et al. 2012). The planet has since been redetected and characterised through several different methods (e.g. Snellen et al. 2014; Snellen & Brown 2018; Hoeijmakers et al. 2018; GRAVITY Collaboration 2020), and a second planet (β Pic c) has recently been announced (Lagrange et al. 2019), further underlining the richness of the β Pic disk and planetary system.

Gas in the disk has been observed both through transmission spectroscopy of material (in some cases, comets) passing in front of the star (e.g. Vidal-Madjar et al. 1994; Roberge et al. 2006; Kiefer et al. 2014), and through emission spectroscopy farther out in the disk (e.g. Dent et al. 2014; Brandeker et al. 2016; Cataldi et al. 2018). In the SW arm of the disk, CO and CI gas have been observed to be concentrated in a clump between 50 and 100 AU projected separation. Because of the dissociation timescales in the high-energy radiative circumstellar environment of β Pic, much of the localised gas is thought to have been released very recently. This supports a scenario in which a large planetesimal was disrupted in the relatively recent past and had its fragments scattered over a range of eccentric orbits (Jackson et al. 2014). Once per orbit, all of the fragments converge back on the collisional point, causing frequent collisions, and ensuring a fresh supply of gas released in the process, thus potentially explaining the observed properties of the gas. While the hypothesised disruption appears to have taken place at a separation in the 50–100 AU range, the eccentric distribution of the aftermath will largely reside at much larger separations. Further clues to explain this event might therefore be found in the outer parts of the disk.

As indicated above, the β Pic planet and disk system has been studied in considerable detail, but as a result of limitations in the field of view (and/or contrast and sensitivity), such studies have typically been limited to separations within a few hundred AU from the parent star. A wider-field coronagraphic image presented in Larwood & Kalas (2001) has revealed that the disk extends much farther than this, although the relatively modest telescope size (2.2 m) limited the sensitivity for characterizing the outermost parts of the disk. In late 2003 and early 2004, a deep-imaging sequence was launched with the FORS1 camera at the 8.2 m European Southern Observatories (ESO) Very Large Telescope (VLT), in program 072.C-0299 with R. Liseau as principal investigator. The data were of very good average quality, but difficulties associated with the outer regions of the bright stellar point spread function (PSF) prevented an accurate analysis of the disk at the time, which led to the data set being shelved. Since that time, considerable progress has been made in the field of high-contrast processing and PSF modeling and subtraction. Meanwhile, the FORS data set from 2003–2004 remains the deepest imaging sequence for β Pic ever acquired (as far as we can identify), and the best chance for characterizing the outer (~500–2000 AU) regions of the disk. Thus, we were motivated to revisit this archival data set in the context of modern PSF subtraction techniques and analyse the PSF-free outer disk. This paper is the outcome of that effort.

The paper is organised as follows: we outline the properties of the archival data and the initial steps of data reduction in Sect. 2. In Sect. 3, we describe the PSF subtraction schemes tested and the final procedure for generating an image that is (to first order) free from the impact of the stellar PSF, with a particular emphasis on mitigating self-subtraction of the disk flux in Sect. 3.2. The results of this procedure and corresponding analysis are presented in Sect. 4. We then discuss some implications of these results in Sect. 5, and finally summarise our conclusions in Sect. 6.

2 Observations and data reduction

All observations for this programme were acquired in the V band with the FORS1 camera at the ESO VLT (the ‘Antu’ unit telescope) between 28 November 2003 and 23 February 2004. The programme was executed in nine observing sessions, but because one had poor ambient conditions, eight data sets are included in the analysis. Each data set contains five image frames of the target β Pic and five frames of a reference star, α Pic. α Pic was chosen as a reference star because it happens to be similar to β Pic in many respects: the brightness and spectral type of β Pic are V = 3.86 mag and A6V, while the corresponding quantities for α Pic are V = 3.30 mag and A8VnkA6. Being separated by about 1h in right ascension, the two stars are also placed well in the sky to be observed back-to-back within an observational set. Each β Pic frame consisted of a single 250 s exposure, and each α Pic frame consisted of a 154 s exposure (matching the integration times after the relative brightnesses). As a result, across the 40 frames per object, the total effective on-source exposure time was 2.8 h for β Pic and 1.7 h for α Pic. With overheads and including both targets, the total VLT time spent on the programme was approximately 7.5 h. This was a partial completion because the initially approved allocation was 13 h.

The high-resolution imaging mode of FORS1 was used, which corresponds to a pixel scale of 0.09975 arcsec pixel⁻¹, and thus a square field of view ~200′′ across. For all astrometric values quoted in this paper, we have taken the mild distortion of the FORS1 field into account using the equation for this purpose given in the FORS manual. During the observation, an opaque bar of width ~47′′ was placed on top of the bright central star (for both α Pic and β Pic) in order to block out as much of the stellar PSF as possible and to simplify detection of faint material in the circumstellar environment. The bar was placed horizontally in the field and was fixed with respect to the sky (not with respect to the telescope pupil), which creates challenges in the data-processing step, as we discuss in Sect. 3.

The fundamental data reduction steps were provided by the ESO FORS pipeline, which includes steps such as bias- and dark subtraction and flat-field corrections. The remaining procedure was performed with a custom reduction procedure. While the star was always placed at approximately the same position during the observations, there were small drifts between frames and small offsets between observation sets, so that in order to register and recentre the frames to a common frame of reference, we used the bright PSF spider arms sticking out from below the mask. By determining centroids along each arm and making linear fits to the results, we determined the stellar photocentre of each frame as the location in which the linear fits to the two sets of arms overlapped. The relative position angle of the x-pattern formed by the spider arms between different frames also gave a high-quality estimate of the relative parallactic angles from frame to frame. The parallactic angle variation within each individual session was 5°–9°, while the variation between sessions was in excess of 120°. These angles were later used in the PSF subtraction scheme described in the next section. To prepare the data for this PSF subtraction, the frames were shifted to a common centre using spline interpolation.

3 PSF subtraction

3.1 Subtraction strategy

The α Pic observations acquired adjacently to each β Pic observation could in principle have been almost ideal for PSF referencing because the two stars are similar in brightness and spectral type and because the observations were obtained relatively closely on the sky and in time. The PSF subtraction scheme would then simply have consisted of finding ideal matching between representations of the α Pic PSF and the β Pic PSF, for example by pairing representations that match closely in time, or by constructing an idealised PSF model for each β Pic frame based on principal component analysis (PCA) using the α Pic frames as a training set (e.g. Soummer et al. 2012; Amara & Quanz 2012). We tried several such strategies, but unfortunately, it turns out that α Pic simply does not function as a suitable PSF reference for β Pic to a sufficient level of precision so that a high final image depth would be reached. As Fig. 1 shows, a pair of α Pic and β Pic frames can display very different PSF substructure patterns, even if they are taken relatively closely together in time. An alternative display of these differences is shown in Fig. 2. This is probably not due to changing seeing or similar, but rather due to scattered or reflected light in the optical system. Throughout this paper, we use terms such as “PSF substructure” loosely to refer to any light registered on the detector that originates from the bright central star, regardless of whether it stems from the seeing halo, from the diffraction pattern, or from reflections or scattering within the optical system. The mismatch between the target frames and the reference frames leads to large residuals and inevitably to a poor image depth in the final collapsed image.

With α Pic no longer considered as PSF reference, the only viable option was to use β Pic itself for PSF referencing. This can be done using various implementations of angular differential imaging (ADI; see Marois et al. 2006), in which several frames of the same target at different rotation states can be used as PSF references for each other. Normally, any observation for which ADI is foreseen would ideally be acquired in so-called pupil-stabilised mode. In pupil-stabilised mode, the telescope pupil is kept fixed on the detector throughout the observation, while the astronomical sky is allowed to rotate relative to the detector as the Earth itself rotates. Because the FORS observations of β Pic were performed before any implementations of modern ground-based ADI¹ existed and did not foresee any such operations, they were executed in standard field-stabilised mode with the sky fixed on the detector. However, a pupil-stabilised framework can be readily produced during post-processing by simply determining the parallactic angle of each frame, and digitally rotating all frames to a common angle through interpolation. As long as the exposures are not long enough to substantially smear the pupil in the field stabilised observations, this framework is a strong analogy for a set of pupil-stabilised observations.

The smearing affects all localised static PSF features (e.g. static or quasi-static speckles), but in the context of these observations, the most notable features subject to smearing are the diffraction spider arms. Given the 250 s exposures, some smearing is inevitable; fortunately, the sky placement of β Pic means that the rotation rate of the field relative to the telescope pupil is quite uniform so that the smearing is largely the same in all frames. This facilitates reproducing and subtracting in an ADI framework. The issue might have been more serious if the target had resided at a declination closer to −24°, where the minimum zenithal angle from the telescope is small. This would have resulted in a highly non-uniform rotation rate, which in turn would produce non-uniform smearing and complicate the PSF subtraction. Hence, while smearing does affect the PSF subtraction quality, it is a manageable effect in the context of our observations. On this topic, we also note that there is in fact a small hidden benefit from using field-stabilised observations for ADI purposes rather than pupil-stabilised observations: in pupil-stabilised imaging, field objects such as planets or disks exhibit smearing to the same extent that PSF features exhibit smearing in field stabilised imaging. Smearing of the field objects is a known complication in conventional ADI (e.g. Lafrenière et al. 2007)and can be a strong effect for near-zenith targets, potentially reducing the signal-to-noise ratio (S/N) of the planet or disk fluxes. In field-stabilised mode, this effect is thoroughly mitigated. Despite this minor benefit, we assess (in retrospect) that pupil-stabilised mode would have been preferable for these observations, and that it is generally preferable for most ADI applications.

Following the reasoning above, the general procedure for ADI processing of the FORS data is (1) to register the centre and relative parallactic angle of each frame using the spider arms of the PSF. (2) For each target frame, to choose a subset of the other β Pic frames to act as PSF references. (3) To shift and rotate all included frames to a common centre and a common pupil angle. (4) To construct an optimal PSF reference frame from the available PSF representations and subtract it from the target frame. (5) To rotate each PSF-subtracted target frame back to a common field angle and collapse them into a final frame, using a median for the collapse. (6) To filter out any remaining smooth radial residuals from the PSF subtraction by taking the azimuthal median at any given radius from the central star and subtracting it from all pixels at that radius (for all radii). (7) To evaluate any self-subtraction of the β Pic disk that may have arisen in the process; substantial self-subtraction of extended flux sources is a common issue in ADI processing. The key steps in this processing are steps 2, 4, and 7, that is, how to best choose the PSF reference set, how to best construct an idealised master PSF from that set, and how to best characterise the self-subtraction that results as a consequence of this procedure.

Under favourable circumstances, these choices may have been nearly trivial. In total, the β Pic observations span a wide range of parallactic angles, and because the β Pic disk is seen nearly edge-on, its morphology is primarily radial, with a relatively limited azimuthal span. Hence, a natural option for any given target frame would be to select almost all other frames for the reference library, removing only those frames with a small (few degrees) relative parallactic angle to the target frame. This is equivalent to imposing a protection angle in conventional ADI-based techniques (e.g. Lafrenière et al. 2007), and mitigates or removes self-subtraction by excluding reference frames where the β Pic disk would overlap with itself in the target frame. The large number of remaining reference frames can then in principle be used to construct a high-quality optimised reference using some reasonably sophisticated scheme such as Karhunen-Loève image projection (KLIP, see Soummer et al. 2012).

Unfortunately, this is not a feasible strategy for this set of observations, for several reasons. Firstly, just as the observations are not stable between α Pic and β Pic as discussed above, they are not stable either between different observing sessions of β Pic itself. The PSF substructure changes from night to night in similar ways as it changes between α Pic and β Pic, so that while the contrast comes out better when using only β Pic than when using α Pic as reference,it is still very far from ideal. This problem is strongly amplified by the specifics of the opaque bar that was used as a form of coronagraphic mask during all the observations. While modern coronagraphic techniques would typically fix the mask to the pupil of the telescope, the FORS bar is fixed to the field. This has the consequence that when rotating any pair of exposures to a common parallactic orientation, the bar maps over two separate angles in the common image space. Because the bar already takes up a very large fraction of the image space for a single orientation, the combined space it takes up for two different orientations, or for an even larger number of orientations across the full sequence, quickly becomes unmanageable. This is particularly important for an algorithm such as KLIP because if a certain pixel position is behind the bar in any individual reference or target frame, the corresponding pixel position must be masked not only in the frame itself, but in every other frame that is used for the PSF modelling and subtraction. In the end, the usable image space only consists of each pixel that avoids the bar in every single target and reference frame, which for all practical intents and purposes is an empty set if all available reference frames are used.

The only practical solution to these case-specific issues is to instead be very selective in terms of reference frames, and in particular to take advantage of the fact that the stability within observing blocks is much higher than between observing blocks. There are five frames per observing block, so that for each target frame, it is possible to use four reference frames without compromising PSF stability by choosing references from other blocks. This is too few frames for algorithms such as KLIP to perform properly, but a classical median-based PSF construction method (Marois et al. 2006) works well for the purpose. All PSF-subtracted target frames can then be combined into a final image, regardless of which observing block they originated from. Using this strategy, we reach a much higher image depth and quality in the final image than with any of the other strategies tested and described above. A downside of this strategy is that the field rotation within an observing block is only ~10°, so that the effective protection angle is not very large in general. As a result, a rather substantial amount of self-subtraction is imposed in the process (see Fig. 3). However, this self-subtraction can be accounted for, as discussed in the next section.

Fig. 1 Images of β Pic (left) and the PSF reference star α Pic (right) at an image stretch focusing on the faint outer parts of the PSF. The two frames are as close together in time as possible with only the slew from one target to the next in between. Significant differences between the two PSF representations still persist, however, making it challenging to use α Pic as reference for β Pic.

Fig. 2 Same as Fig. 1, but here the two images have been rotated to a mutual pupil angle and one has been subtracted from the other. This allows us to distinguish the systematicPSF differences more clearly, which include asymmetrically distributed PSF halos, different sharpness for the spider arms (probably due to different degrees of pupil smearing between α Pic and β Pic in the field-stabilised observations), and most critically for our purposes, a range of artefacts in the outer parts of the image domain, which are not in common in the β Pic and α Pic images.

Fig. 3 ADI-subtracted image of β Pic and its circumstellar debris disk. Self-subtraction has inevitably taken place in the procedure, as shown by the dark traces on each side of the disk. Still, the disk flux is detectable out to very large separations.

3.2 Accounting for self-subtraction

A central method for correcting for self-subtraction and derive subtraction-corrected parameters for disks in ADI processing is through forward-modelling with negative disk injection (see e.g. Thalmann et al. 2014; Janson et al. 2016; Lagrange et al. 2016). In this framework, a parametric model of the disk is constructed and subtracted from each individual frame (at its original orientation) prior to the ADI reduction. The ADI reduction is then performed as normal, and at the end of the procedure, the residuals in the disk region are calculated. If the model matches the actual properties of the disk, the residuals are small, while if the model is a poor match, the residuals are larger. Thus, it is possible to construct a grid of parameter values and run a new ADI reduction for each combination of parameters in the grid. The minimum residual case can be considered as the best parametric fit to the actual data. A perfect fit is fully unaffected by self-subtraction because after subtraction of such a model from the science data, there is no more disk flux that can overlap with itself during the PSF subtraction. We therefore implemented such a scheme on the β Pic data in the manner described below.

The parametric disk model can be constructed in several different ways, depending on the implementation. For example, if the disk consists of a single dust belt with an eccentric ring-like shape such as in HR 4796 (Milli et al. 2017), a code such as GRATER (Augereau et al. 1999) can be used to model a simple dust distribution and derive a flux distribution from radiative transfer calculations. In the case of β Pic, where the inclination is close to edge-on and we suspect that multiple disk components are overlapping, we chose an approach more like previous work that has been done for β Pic at smaller separations by Ahmic et al. (2009), for instance. First, we chose an approximate position angle to represent the disk plane. The two arms of the disk were treated separately with individual angles in this regard because there is a clear angular offset between them, which we return to later. Along the axis defined by the approximate disk midplane, we then characterise the disk as slices of vertical flux distribution parametrised by their width, shape, and amplitude. The analysis was made in a context wherein the image of the disk was rotated such that it was parallel to the horizontal axis of the pixel array (again, separately for the two sides of the disk). Slices of the disk are therefore represented by vertical columns of the pixel array. As a first step, we simply fit Gaussian functions to pixel columns in the self-subtracted ADI image. The free parameters in the fit are the amplitude, width, and centre of the Gaussian, and the minimum χ² solution gives an initial estimate for these quantities. These estimates are inevitably skewed due to the self-subtraction, which affects the apparent width and amplitude of the disk slices, and because the real vertical profile differs somewhat from a Gaussian (as we verify in subsequent steps), but it provides a good starting point for defining the grid values of the negative disk-injection procedure that follows as the next step.

Because the negative injection procedure attempts to minimise residuals in the disk region, we need to define a specific region in the final image space in which the residuals are evaluated. For this purpose, we chose a rectangular box centred on the disk midplane that was 6′′ wide in the vertical direction and stretched from separations of 26′′ to 98′′ in the horizontal direction on the brighter NE side, and from 31′′ to 59′′ on the fainter SW side. In addition to the faintness of the SW arm, there is also a group of faint extended objects (probably a clustering of background galaxies) just outside of 59′′, preventing an accurate analysis beyond that point. The inner separations are slightly larger (by ~2′′) than the separation range covered by the central bar, which arises from the fact that the bar edge is non-orthogonal to the disk plane, leaving a thin wedge where a fraction of the disk is visible, but cannot be sampled symmetrically with respect to the disk plane. We simply disregarded this minor part of the image space in our analysis. The vertical profile is characterised by the function , where p_a is the amplitude parameter, p_w the width parameter, and p_s is the shape parameter, where for example p_s = 1.0 implies an exponential drop-off. In the inner ranges of the β Pic disk, the shape has been found in previous studies to be super-exponential, that is, with values of p_s < 1 (Ahmic et al. 2009). Before we subtracted the corresponding vertical profile from a given location in the disk, we convolved it with a PSF as derived directly from one of the brighter (but non-saturated) background stars in the field. Because an exponential or near-exponential function is very sharply peaked near the midplane, the PSF convolution is an important step even for vertical profiles with a quite large p_w.

Because different parts of the disk interact azimuthally but not radially in a classical ADI reduction, the optimisation problem of finding a best-fit model becomes separable in the radial direction. In principle, the disk could be divided into radial sections where each section is assigned its own set of parameters prior to the run, but this would be inefficient (or even computationally unmanageable) because every set of parameters needs to be run through its own separate ADI reduction, which is a time-consuming procedure. A much more efficient strategy is to apply parameter grid steps uniformly across the whole disk, except for the output image after each ADI reduction, and then evaluate the residuals locally among the output images. In other words, allowing for different parameter sets to yield the minimum residuals in different radial sections of the disk. In the case of the amplitude, each grid step is a uniform factor that is multiplied with the initial estimation of the radial brightness profile. Because the self-subtraction that we evaluated is expected to vary smoothly with separation, it is not necessary toestimate it in every individual column, and indeed, such a procedure would be strongly affected by random noise in the outer ranges of the disk where the flux is low. Hence, while the innermost 20 pixels of the NE side of the disk are very bright and thus could be easily evaluated individually, for the remaining disk we performed the evaluation over ranges of 100 pixels (approximately 12.5 resolution elements for a seeing of ~0.8′′) at a time on the NE side. The total radial extent of the NE side is approximately 1000 pixels, or ~125 resolution elements, of which 740 pixels are used outside of the central bar. Individual column corrections were then assigned through linear interpolation. The procedure was similar on the SW side, but because the disk is less extended and has no very bright part there, evaluations were performed in steps of 25 pixels across that arm. The procedure was performed in several steps in order to ensure that a global minimum could be reached and to gradually fine-tune the parameters, and initially, this was done only for variations in amplitude and width, keeping the shape fixed at p_s = 0.83, which is what was found for this parameter closer in Ahmic et al. (2009). However, because the intermediate best-fit solution still was not perfect and the residuals visually suggested that the vertical shape was off, we ran an additional procedure in which the shape parameter was allowed to vary along with the width and amplitude. In this iteration, the residuals were evaluated globally across the whole disk region at once. This yielded a best-fit p_s = 1.12 ± 0.32, which we use for the remainder of the analysis. Here, the estimated error has been set as the threshold for which the fit residuals in the evaluation region are within 10% of the minimum residuals. In this context, the slightly wider profile than in the inner disk (0.83, Ahmic et al. 2009) is only marginally significant. The S/N per 1′′ circular footprint at the disk midplane at the NE side is 30.0 at 650 AU. At 900 AU, the disk is approximately two times fainter but the PSF residuals are also smaller, leading to an S/N of 24.9. At the SW side, the S/N is 15.1 at 650 AU and 6.9 at 900 AU.

While some local residuals can still be seen in the resulting images for the NE (see Fig. 4) and SW (see Fig. 5) arms of the disk, they are considerably reduced by the fitting procedure. Most likely, the vertical profile is not perfectly represented by the parametric equation that we used, or is a combination of several overlapping components with different parameters. Once provided with the best-fit parameters from the above procedure, we can produce an output image which, to first order, is free from self-subtraction effects. This was done in an iterative process, where the starting point was a tailored ADI procedure in which the best-fit model disk was subtracted from the PSF reference frames, but not from the target frames. The result was an ADI-reduced image in which the disk was visible and free from ADI-induced self-subtraction features. The exact same ADI procedure was then repeated, but in this step, the output image from the previous ADI reduction was used as the input model for the new one. Any negative value in this semi-empirical input model was set to zero because it is an unphysical value. This procedure was then repeated iteratively. With every step the model is expected to be a more accurate representation of the actual disk morphology, and thus subtracting it from the PSF reference frames should yield a yet more accurate output image to be used as a model in the next step. We can assess how well this works by monitoring what happens to the background stars and galaxies in the FORS field of view for which no explicit optimisation has been performed during any step of the process. Successively with increasing iterations, the fidelity of the background sources increased, with decreasing background shadows around them, confirming that the procedure worked as intended.We find that the output results effectively converge after five iterations, where further iterations did not significantly change the output result. Thus, we took the fifth iteration image as our final output, free (to first order) from self-subtraction effects. The latter iterative part of the procedure is similar to the scheme devised by Pairet et al. (2018) for the same purpose of achieving accurate disk morphologies in ADI-subtracted data.

In general, this iterative procedure can be expected to work increasingly well at large separations because self-subtraction effects are smaller in this regime to begin with. This means that the template image that is subtracted from the PSF frames in each iteration is more accurate, such that a stronger convergence can be reached faster. At smaller separations, the input image is less accurate, and the convergence is consequently slower or worse. Because physically unrealistic (negative) regions are masked out in the procedure, corrections to the PSF using a non-perfect template should always be as good or better than using no correction at all, however, so that while convergence may be limited in the innermost regions, there should be little or no risk of any divergent solutions. As mentioned above, we are fortunate to have background stars spread around the field, allowing us to check for such effects. In contrast to extended sources such as disks, point sources have a precisely known morphology, so that any effects imposed by the algorithm can be monitored through how they affect the point sources. We find that point sources far out in the field (even very bright ones) lack any discernible subtraction shadows around them in the final output. The closest-in point sources do have weak shadows around them, but are nonetheless well represented by point-source morphologies. This confirms that the procedure does work best at large separations and becomes increasingly challenging at smaller separations, but that even at the smallest separations probed in this study, the algorithm works and improves on the pre-iteration results.

Fig. 4 Residuals after negative injection of a disk model prior to the ADI processing for the NE arm of the disk. Some residuals remain in the region closest to the star, which may be related to the detailed vertical profile of the disk. However, they are greatly reduced by the fitting procedure.

Fig. 5 Same as Fig. 4, but for the SW arm of the disk.

4 Results and analysis

Despite all the challenges related to the PSF subtraction described in the previous section, the final output image shown in Fig. 6 is of good quality, and constitutes the deepest imaging of the β Pic system acquired to date. We also show a version of the image smoothed with a Gaussian kernel with a full width at half maximum (FWHM) of 0.6′′ and rotated to the conventional sky orientation in Fig. 7.

We extracted a surface brightness profile along the disk midplane through aperture photometry using a circular aperture with a diameter of 0.6′′, which is close to the FWHM of the seeing-limited PSF. Uncertainties were calculated based on the residuals in a disk-subtracted image, where the iterative disk image was subtracted from both the target and PSF frames, instead of from the PSF frames alone. The result was converted into units of mag arcsec⁻² by normalising for the size of the extraction aperture. This conversion also requires a photometric zero-point, which we determined using background stars in the FORS field around β Pic. For this purpose, we selected five stars that are all recorded in the Gaia DR2 catalogue (Brown et al. 2018), bright enough to have a very good significance but not bright enough to be non-linear or saturated in the FORS images, and at a substantial distance from β Pic itself in order to minimise the effect of its PSF halo on both the FORS and Gaia data. The calibration stars are summarised in Table 1. Because the FORS data are in V band, a photometric conversion from Gaia magnitudes is required, which we acquired using the G_BP − G_RP colour according to transformation equations given in the ESA online Gaia documentation.

In Fig. 8 we show the resulting surface brightness profile in calibrated photometric units. The uncertainty in the photometric calibration is dominated by the scatter among the calibration stars and is 0.14 mag. This scattermight arise from the effect of the bright stellar halo on the Gaia data. PSF effects around bright targets are known to have substantial effect on stars around them in DR2 (e.g. Brandeker & Cataldi 2019); this will most likely improve in subsequent Gaia releases. Because the zero-point calibration uncertainty affects all brightness profile data points systematically in the same way, we plot it separately in Fig. 8 rather than addingit in quadrature to the intrinsic photometric uncertainties there. The calibration clearly dominates the uncertainty in the absolute level in the bright inner parts of the disk (particularly in the NE arm), while intrinsic uncertainties dominate in the faint outer parts. In the outer parts of the FORS field of view, the median 3σ sensitivity (per PSF footprint) is 27.7 mag in V band. The previous deepest image of the β Pic disk at any comparable wavelength was presented in Larwood & Kalas (2001), who achieved an R-band 3σ sensitivity of 25.8 mag at these separations. Thus, the high depth of the FORS image makes it a valuable resource for characterising the outer ranges of the β Pic disk.

In addition to the general smooth drop-off in brightness with increasing distance from the star, there is apparent substructure in the form of bumps at specific radial locations. This substructure in the β Pic disk has been noted before in the inner parts of the disk (e.g. Kalas et al. 2000; Wahhaj et al. 2003). One possible interpretation of these bumps is that they represent an ansa from a discrete ring in the edge-on disk, but another possibility that we need to account for, in the context of the outer ranges of the β Pic disk, is an overlap between the disk and a background galaxy. In these ranges, we are probing very faint flux at considerable image depth, which means that we are at the border of being confusion limited with respect to galaxies in the background. Many background galaxies can indeed be seen across the FORS field, therefore it would not be surprising if some of them overlapped or partially overlapped with the area covered by the disk.

One bump that stands out in the surface brightness profiles in Fig. 8 is located in the NE arm of the disk between 750 and 800 AU. This feature has in fact been noted before as feature ‘A’ in Kalas et al. (2000). Because almost a decade of baselines has passed between the original image of feature A from 1995 and the FORS image from 2004, we can distinguish between a disk feature and a background feature on the basis of whether it shares a common proper motion with the star and disk. We find that the photocentre of the bump resides at a separation of ~39.7′′ in 2004, while in 1995, it resided at a separation of ~40.7′′. While astrometric centring on the bump is very difficult due to its overlap with the disk, the precision is certainly much better than the FWHM of the PSF of 0.6′′. This means that the bump has certainly moved with respect to the star. Moreover, the NE proper motion of β Pic of μ_RA = 4.65 mas yr⁻¹ and μ_Dec = 83.10 mas yr⁻¹ (van Leeuwen 2007) means that a static background object would move in a relatively radial direction toward the star, and by a similar amount over a decade as observed for bump A. We thus conclude that the most probable interpretation of bump A is a background galaxy. If this interpretation is correct, future imaging of the bump should place it at an even closer separation to the star, and it should gradually start separating from the NE arm on the southern side of it. Another conspicuous bump in the FORS data occurs at the outer edge of where we measure the SW arm. As we noted previously, the outer edge of our SW arm measurements is set by the fact that a group of faint objects overlap with the extended direction of the arm. These objects appear to be extended and may represent a clustering of galaxies in the background. The measured bump mightmark the first such galaxy, and we thus consider it unlikely that this is a feature within the disk itself. Future imaging for relative proper motion measurements is required to fully establish the nature of this bump because it is smoothly extended and quite faint and therefore difficult to visually distinguish from disk flux in the images. The outermost range (beyond ~1400 AU) in the NE arm also seems to indicate some broad bump, which is subject to the same uncertainties as the outer bump on the SW side.

Because of the suspected (and in one case verified) background nature of the blobs discussed above, we rejected the corresponding areas when making power-law fits to the general flux slopes of the arms. This was made through line fitting to the flux distribution versus separation in logarithmic space, defining α as the power-law index on the projected radial dependence r^α. In the NE arm within 1200 AU we find that α_NE = −3.7 ± 0.1, and when we included the whole range out to nearly 2000 AU, the value remained very similar at α_NE = −3.9 ± 0.2. In the SW arm out to 900 AU, the slope is similar or very slightly steeper at α_SW = −4.0 ± 0.3, but then becomes significantly steeper such that the average slope of the full arm out to 1050 AU (where the outer bump starts) is α_SW = −6.0 ± 0.8. If interpreted as a power-law break around 900 AU, the outer slope out to 1050 AU is extremely steep at α_SW = −18.0 and would imply a sharp outer cut-off for the SW arm, but we note that this latter value is based only on a small section in one of the faintest parts of the disk. The general trends are consistent with the fact that the brightness asymmetry between the NE and SW side, which is clearly present also in the inner parts of the disk (e.g. Kalas & Jewitt 1995), is even more extreme in the outer parts that we probe here. The slopes on either side out to ~1000 AU are very similar to the corresponding slopes at ~100–300 AU calculated by Kalas & Jewitt (1995) as − 3.76 ± 0.05 on the NE side and − 4.07 ± 0.05 on the SW side.

The vertical height and trace of the disk can also be derived from the ADI-corrected PSF-subtracted image, but it is a more challenging prospect than for the surface brightness profile because the brightness profile is measured in the disk midplane where the flux is highest, but the height and trace depend strongly on the behaviour in the wings extending vertically from the midplane, where the flux is lower and more susceptible to the confusion limit mentioned above. We therefore restricted this analysis to within separations of 50′′. We performed an additional round of fitting for the width and trace, where we fit vertical profiles with the p_s = 1.12 exponent determined in Sect. 3.2, and with the surface brightness profile in the midplane fixed as determined above. In this round of fitting, we downsampled the data by averaging groups of 10-pixel columns and individually fitted the height and centroid of each averaged column. The best-fit heights are shown in Fig. 9. The SW arm appears to have a generally broader profile at these separations, with a mean height of 22.3 ± 1.4 AU compared to the mean height of the NE arm of 13.4 ± 0.9 AU.

The disk traces of each arm are shown in Fig. 10 with an alternative display in Fig. 11. They confirm the tilt angle between the arms that we noted previously. We made linear fits to the traces and related them to the average position angle of the inner disk of ~30.8 deg (e.g. Kalas & Jewitt 1995). The NE side has a differential angle of 2.0 deg in a clockwise (CW) direction in the image plane with respectto the reference angle, while the SW side has a differential angle of 5.2 deg in a counter-clockwise (CCW) direction. The total tilt angle between the arms is therefore 7.2 deg. This is higher than in the inner parts of the disk (e.g. 1.3 deg in Kalas & Jewitt 1995). A tilt angle between two arms of a nearly edge-on disk arises naturally if there is strong forward-scattering in the disk and if the disk has a non-zero angular offset from being perfectly edge-on. In this context, because the tilt angles are asymmetric relative to the inner disk, it would require an angular offset between the inner and outer disk planes not only in inclination, but also in the disk position angles. The higher the degree of forward-scattering and the larger the offset angle, the higher the expected tilt angle.

We can thus use the tilt angle to constrain the inclination and scattering anisotropy, but with some degeneracy between them. To do this, we generated simple models with different inclinations and degrees of forward-scattering in the GRATER code (Augereau et al. 1999). In GRATER, the forward-scattering can be quantified by the Henyey-Greenstein index g, which has a value between −1 (maximally back-scattering case) and +1 (maximally forward-scattering case). We input values of 0.2 to 0.9 in steps of 0.1 for g and inclinations i of 85 to 89 degin steps of 1 deg. In addition, because the tilt angle can in principle depend on what the location is relative to parent planetesimal belt that the dust is expected to originate from, we tried three different values of the belt radius r₀ : 400, 500, and 600 AU. The dust density in the model drops off slowly from the ring centre, consistent with the trend observed in the surface brightness profile; in this case represented by a power-law index in the density of −1.8. For each simulated disk, we evaluated the tilt angle between the two arms of the disk by fitting Gaussian centroids to sequential vertical slices of the disk between 30′′ and 50′′ (~600–1000 AU), and then weperformed a linear fit to the resulting centroids to find the average slope of the arm relative to the main disk plane. The model disk is circular and single-component, and is not meant to represent the actual morphology of the disk (which is much more complex), but only to provide a handle on what is required to produce tilt angles of the size observed.

The dependence of the tilt angle on r₀ is found to be generally weak in this context, but the dependence on both i and g is strong. Interestingly, we find that both a higher degree of forward-scattering and a larger offset from edge-on (i.e. smaller inclination) is required for the outer disk than has previously been determined for the inner primary disk. A g < 0.8 or an i> 86 would both have great difficulty to reproduce the tilt angle. A combination of i = 85–86 deg and g = 0.8–0.9 is necessary to acquire a tilt angle of 7.2 deg. While i = 85–86 deg is consistent with the inclination range determined in Kalas & Jewitt (1995), it is lower than the i = 87.7 deg determined for the inner primary component of the disk in Ahmic et al. (2009). The g in the inner range is determined as 0.3–0.5 in Kalas & Jewitt (1995) and as 0.64 and 0.85 for the primary and secondary components of the disk, respectively, in Ahmic et al. (2009). A value of g = 0.8–0.9 for the outer disk is similar to the secondary disk in the inner range, but higher than the primary disk. g = 0.9 would be consistent with the degree of forward-scattering for dust, for example, in comet McNaught (C/2006 P1; see Marcus 2007).

The difference in forward-scattering could therefore be explained by a gradient in the dust properties in the disk. The inclination might, for instance, imply a warp in the disk, where the outer disk has a systematically different orientation than the inner disk, much like how there are two inner disk components with systematically different orientations already known in the system. The asymmetry in the tilt angle of the two arms that we measured with respect to the inner disk plane would further support this scenario. Because a considerable brightness asymmetry between the arms exists and probably implies at least some strongly eccentric subcomponent of the disk, this might also contribute to the asymmetry in the tilt angle: the two arms might then contain relative populations that are asymmetrically distributed outside of the image plane, with one population of dust primarily in front of the star at one side and another population primarily behind the star at the other, affecting their relative scattering properties even if the underlying population of dust is the same.

Fig. 6 Final image after subtracting a best-fit model from the PSF reference frames. The self-subtraction normally imposed by ADI is removed (to first order) by this procedure.

Fig. 7 Smoothed version of the final image, rotated to the conventional sky orientation with north up and east to the left.

Fig. 8 Surface brightness profile for the brighter NE arm (blue) and the fainter SW arm (red) of the disk. The bumps that deviate from the average trends are probably caused by background galaxies (see text). The SW arm measurements start at a larger projected separation than for the NE arm because of a small asymmetry in the placement of the star behind the central bar.

Fig. 9 Vertical height of the brighter NE arm (blue) and the fainter SW arm (red). The SW arm appears to be slightly broader on average at these separations in terms of its scale height above the midplane.

Fig. 10 Disk trace of the brighter NE arm (blue) and the fainter SW arm (red). The disk has a clear tilt, with both arms deviating from the average disk plane at smaller separations, which isthe reference plane in this figure. This implies a high degree of forward-scattering and/or asignificant deviation from an edge-on inclination (see text). The spatial scalings of the x- and y-axes are very different, which causes a strong visual exaggeration of the tilt angles.

Fig. 11 Same as Fig. 10, but with the absolute value of the horizontal offset (i.e. mirroring the NE arm) in order to enhance visibility.

5 Discussion

As we have noted above, there are multiple asymmetries in the β Pic disk that become increasingly extreme farther out in the disk. Asymmetries are present not only in the dust, but also in the second-generation gas that has long been known to exist in the disk (Slettebak 1975; Olofsson et al. 2001; Brandeker et al. 2004). The dissociation timescales of the detected molecules imply that the gas may have been released in a series of events following a disruption of a large planetesimal in the recent past (Cataldi et al. 2018). The disruption might have occurred as a result of a collision with another planetesimal (Jackson et al. 2014; Lawler et al. 2015) or through tidal disruption in a close encounter with a massive planet nested within the disk (Cataldi et al. 2018; Janson et al. 2020). A large number of fragments of different sizes would have been formed in the disruption, and acquired a range of eccentric orbits with different orbital elements, but which all converge once per orbit to the disruption location where they originally formed. Over time, this creates a series of additional collisions at this location, which can maintain a localised release of short-lived gases over several thousand years (several orbital timescales). The observed gas in the disk is predominantly found in the SW arm at a projected separation of ~50–100 AU, which marks the origin of the event and the mutual periastron in this scenario. Meanwhile, the fragments and associated dust would be spending most of their time at apastron, at substantially larger projected separations, at the other side (i.e. the NE arm) of the disk. It might therefore be the case that the excess flux associated with dust that we observe in the NE armcould in part be related to the recent disruption event. In this framework, it might be viewed as one highly eccentric component of the disk superimposed on a less eccentric, more ordered outer disk around β Pic, much like how there are (at least) two distinct components in the inner disk.

Despite about 17 yr of technical developments, it is a far from trivial prospect to improve the image depth and wide-separation contrast at visible wavelengths acquired in the observations presented here. As mentioned in Sect. 2, the observational programme spans about 7.5 h of VLT time, of which 2.8 h are effective integration time on β Pic, under good seeing conditions and in a dark lunar phase. The disk is probed down to a surface brightness of 27–28 mag arcsec⁻². The sky brightness (~22 mag arcsec⁻² at Paranal during the time of observation, see Patat 2003) is obviously an important factor in this regard, but still, the large size of the VLT and high sensitivity of FORS gives it a significant advantage even with respect to space-based facilities. Acquiring a 27 mag arcsec⁻² sensitivity in a 5 × 5 pixel box at 5σ significance with the WFC3 camera at the Hubble space telescope, for example, would require 39 h of effective on-source integration time, and the field of view of 162′′ would be insufficient to cover the whole disk in one shot. From the ground, PSF-limited visible light imaging has benefited from adaptive optics developments (e.g. Schmid et al. 2018; Close et al. 2018), but because the isoplanatic angle is small, the good performance reached within a few arcseconds of separation cannot be extended to the image scales required to cover the β Pic disk. Nonetheless, in possible future imaging efforts for probing the outer disk, some clear improvements might be made based on our experiences with the FORS data set. Most importantly, the bar used to mask the star in the images was fixed to the sky rather than to the pupil. This made any highly sophisticated PSF subtraction schemes essentially impossible. If an observational setting could be used in which the mask would be fixed to the pupil, and ideally, in which the detector would also be pupil stabilised, this would greatly enhance the achievable contrast performance. Observations should be acquired in long uninterrupted sequences because the largest PSF variations were experienced between nights (or between target and reference star); this procedure would also allow us to accumulate a substantial amount of parallactic angle variation across the observation, which further benefits contrast in ADI-related applications. Observations in different photometric bands than V may be particularly useful because the colour of the dust might help constrain dust sizes and similar characteristics. A first indication of dust colour might be provided by the fact that the β Pic disk appears to be at least marginally visible out to ~1900 AU in the Larwood & Kalas (2001) R-band image, even though that image is ~2 mag less deep than the FORS image. This could imply a redder colour in the scattered disk flux than for the central star, which has a typical A-type colour of V − R = 0.12 mag. Polarimetric imaging could be an alternative or complementary characterisation pathway, although its success heavily relies on a reasonable high degree of polarisation for the dust involved.

Having been probed from single-AU distances in the form of transiting comets (e.g. Vidal-Madjar et al. 1994; Zieba et al. 2019) out to ~2000 AU in the FORS images presented here, the observational data for the β Pic debris disk cover a truly vast range, and provide information on its characteristics at distances from the star that are comparable to the inner Oort cloud in the solar system; a unique feature in the known debris disks. We also note that radar measurements tracking micrometeorites that enter the Earth’s atmosphere have led to the detection of a population of particles whose kinematics imply an interstellar origin (Baggaley 2000). Trace-back studies have pointed out β Pic as the most probable origin of this population (Krivov et al. 2004). It might therefore be argued in this sense that the observable range of the β Pic disk even reaches out to scales of tens of parsecs.

6 Conclusions

We have examined a deep FORS archival imaging data set of β Pic in order to determine the fundamental properties of the outer ranges of its debris disk. One of the main challenges in this context is the complex and varying PSF structure of the bright primary, particularly because these observations were acquired before ADI became widely applicable. This meant that special considerations had to be taken in the subtraction procedure, and that the subtraction had to be performed over data sets spanning a rather small range of parallactic angles. As a consequence, substantial self-subtraction was inevitable as a consequence of the procedure, and in response to this, a dedicated negative injection procedure had to be implemented to mitigate these effects.

Aside from the impressive scale spanned by the far-reaching arms of the disk, its most prominent characteristic is the extreme asymmetry between the NE and SW arms of the disk. The asymmetry exists already in the inner disk region, but it becomes increasingly enhanced with increasing separation, and manifests itself both in the surface brightness, in the vertical height, and in the trace of the disk with respect to the midplane. One apparent component of this asymmetry includes discrete bumps in the surface brightness profiles, but we note that at the faint brightnesses that are being probed, the observations are in a partially confusion-limited regime with respect to background galaxies, and we thus conclude that the most likely underlying cause of the bumps is blending with background galaxies. In the case of the most prominent bump of the NE arm, we confirmed this hypothesis by relating its position to previous epoch imaging and noting that its motion resembles that of a static background object. The asymmetries in the total brightness, brightness slope, height, and trace are however all real physical effects within the disk, and reflect asymmetries in the morphology of the underlying dust. Eccentric subpopulations of dust are a likely cause for many of the observed feature, some of which might be related to the localised gas emission previously reported in the SW arm of the disk, which in turn has been theorised to be caused by a recent large disruption event within the disk.

Acknowledgements

M.J. gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation. This study made use of the CDS services SIMBAD and VizieR, and the SAO/NASA ADS service. We thank Matthias Samland for useful discussion, and the anonymous referee for the careful review and useful comments that helped improve the manuscript.

References

All Tables

All Figures

	Fig. 1 Images of β Pic (left) and the PSF reference star α Pic (right) at an image stretch focusing on the faint outer parts of the PSF. The two frames are as close together in time as possible with only the slew from one target to the next in between. Significant differences between the two PSF representations still persist, however, making it challenging to use α Pic as reference for β Pic.
In the text

Fig. 2 Same as Fig. 1, but here the two images have been rotated to a mutual pupil angle and one has been subtracted from the other. This allows us to distinguish the systematicPSF differences more clearly, which include asymmetrically distributed PSF halos, different sharpness for the spider arms (probably due to different degrees of pupil smearing between α Pic and β Pic in the field-stabilised observations), and most critically for our purposes, a range of artefacts in the outer parts of the image domain, which are not in common in the β Pic and α Pic images.

In the text

	Fig. 3 ADI-subtracted image of β Pic and its circumstellar debris disk. Self-subtraction has inevitably taken place in the procedure, as shown by the dark traces on each side of the disk. Still, the disk flux is detectable out to very large separations.
In the text

	Fig. 4 Residuals after negative injection of a disk model prior to the ADI processing for the NE arm of the disk. Some residuals remain in the region closest to the star, which may be related to the detailed vertical profile of the disk. However, they are greatly reduced by the fitting procedure.
In the text

	Fig. 5 Same as Fig. 4, but for the SW arm of the disk.
In the text

	Fig. 6 Final image after subtracting a best-fit model from the PSF reference frames. The self-subtraction normally imposed by ADI is removed (to first order) by this procedure.
In the text

	Fig. 7 Smoothed version of the final image, rotated to the conventional sky orientation with north up and east to the left.
In the text

	Fig. 8 Surface brightness profile for the brighter NE arm (blue) and the fainter SW arm (red) of the disk. The bumps that deviate from the average trends are probably caused by background galaxies (see text). The SW arm measurements start at a larger projected separation than for the NE arm because of a small asymmetry in the placement of the star behind the central bar.
In the text

	Fig. 9 Vertical height of the brighter NE arm (blue) and the fainter SW arm (red). The SW arm appears to be slightly broader on average at these separations in terms of its scale height above the midplane.
In the text

Fig. 10 Disk trace of the brighter NE arm (blue) and the fainter SW arm (red). The disk has a clear tilt, with both arms deviating from the average disk plane at smaller separations, which isthe reference plane in this figure. This implies a high degree of forward-scattering and/or asignificant deviation from an edge-on inclination (see text). The spatial scalings of the x- and y-axes are very different, which causes a strong visual exaggeration of the tilt angles.

In the text

	Fig. 11 Same as Fig. 10, but with the absolute value of the horizontal offset (i.e. mirroring the NE arm) in order to enhance visibility.
In the text