© The Authors 2023

1 Introduction

Circumstellar debris discs have been detected around a significant fraction (15–30%) of main sequence stars (Hughes et al. 2018). These discs are thought to result from the continuous collisional grinding of leftovers from the planet-formation process, of which only the tail end (dust particles) is detectable as a photometric infrared (IR) excess (Wyatt 2008). An increasing number of these discs have also been imaged, most commonly by scattered-light observations in the visible and near-IR, but also in thermal emission in the mid-to-far IR to millimetre (mm) wavelength domain¹. These images reveal a variety of structures, such as clumps, warps, spirals, and arcs, which have been interpreted as being sculpted by gravitational perturbations of (usually unseen) planets (e.g. Augereau et al. 2001; Wyatt 2006; Thebault et al. 2012) or companion stars (Thebault et al. 2021), or as being the result of violent and transient collisional events (Jackson et al. 2014; Kral et al. 2015; Thebault & Kral 2018).

Apart from the aforementioned spatial structures, most imaged discs share a common feature, which is that the dust surface density is not uniformly decreasing (or increasing) with stellar distance, but peaks at a given radial location, thus creating ring-like features. This belt-like configuration is so ubiquitous that it has been suggested that “debris ring” should be a more appropriate denomination of the debris disc phenomenon (Strubbe & Chiang 2006). In many cases, scattered-light images also show a dimmer but radially extended halo of dust beyond the location of this main ring or belt (Hughes et al. 2018). The analytical and numerical studies of Strubbe & Chiang (2006, hereafter STCH) and Thébault & Wu (2008, hereafter TBWU) have shown that such halos are in fact expected beyond collisionally evolving belts of parent bodies (PBs). These halos should be made of small, typically 1–10 µm sized grains that are collisionally produced within the main ring and then placed on high-eccentricity orbits by radiation pressure or stellar wind. For systems dense enough for the Poynting-Robertson (PR) effect to be negligible with respect to collisional evolution, the vertical optical depth τ in the outer halo naturally tends towards a radial profile in τ ∝ r^−1.5. At wavelengths dominated by scattered light, this translates into a surface brightness (SB) profile that decreases as ∝r^−3.5. This SB ∝ r^−3.5 slope also holds for the projected midplane SB profile of discs seen edge-on, provided that their vertical scale height z is ∝ r (Strubbe & Chiang 2006). However, we note that these results are valid under several simplifying assumptions, notably isotropic scattering and the absence of small unbound grains that are blown out by radiation pressure or stellar wind.

Thanks to cutting-edge instruments such as SPHERE or GPI (Beuzit et al. 2019; Macintosh et al. 2014), recent years have seen an exponential increase in the number of systems for which the radial profile of halos has been constrained (Adam et al. 2021). For these systems, the canonical radial profiles of STCH and TBWU can be used as benchmarks to judge as to whether or not the outer edge of a debris ring is ‘natural’ or, on the contrary, shaped by additional mechanisms (outer planet, interaction with gas, companion star, etc.). However, the quantity the radial profile of which is constrained by observations can differ from one study to another. In some cases, it is the SB that is retrieved, which can either be the (deprojected) stellocentric SB measured along a radial cut (Schneider et al. 2018) or the projected SB in the case of discs seen edge-on (Golimowski et al. 2006). In other cases, it is the underlying grain number density profile that is estimated, usually by fitting the observed resolved images with radiative-transfer models, such as the GRaTer code (Augereau et al. 1999), in which several free parameters (parent belt location and width, radial slope in the halo, scattering phase function, etc.) are explored and adjusted (Choquet et al. 2018; Bhowmik et al. 2019; Perrot et al. 2019). The radial density profile that is estimated can either be the surface density σ or the volumic number density n.

The fact that – depending on the profile-fitting approach – the notion of radial profile can refer to three distinct parameters (SB, σ, and n) sometimes leads to some confusion when comparing different discs or when comparing observations to canonical theoretical profiles. An additional problem is that the theoretical results of STCH and TBWU constrain the system’s optical depth τ (or geometrical cross section), and not the particle number densities σ and n. If the disc is only made of identical particles, then the radial profiles of τ and σ (and also n × r for non-flared discs) should be equivalent², but this is not the case in situation studied here because there is a very strong size segregation as a function of radial distance in the halo (Thebault et al. 2014). Schematically, at a given stellar distance, r, the cross section is indeed dominated by grains produced in the main ring that have their apoastron Q = r. As the value of Q is imposed by the size-dependent mechanism that is radiation pressure (or, for M stars, stellar wind), it follows that, at a distance r, τ is dominated by grains of size s, such that $$\beta \left(s\right)=\frac{1}{2}\left(1-\frac{r_0}{r}\right),$$(1)

where r₀ is the location of the main belt and ß(s) is the ratio between the radiation pressure and stellar gravity (or stellar wind) forces for the size s.

In contrast, the number-density profiles derived in most GRaTeR-type fits assume that the particle size distribution is the same everywhere in the system, which clearly cannot apply to the size-segregated halo. This means that neither σ nor n derived in this way correspond to actual number densities (neither that of the global population nor that of a given size range). Nevertheless, for fits performed on images at short scattered-light wavelengths, the radial dependence of the GRaTeR-derived σ parameter should, to a first approximation, match that of τ, as long as grains in the halo contribute to the flux in proportion to their cross section. As shown in Sect. 3, this assumption does not hold if small unbound grains make a significant contribution and if the scattering phase function is dependent on size; however, before exploring these important issues, let us present in Table 1 the first coherent census of all systems with observationally constrained halo profiles for which we specify which quantity (SB, n or σ) has had its radial profile constrained. It can be noted that the radial index Γ of the surface density profile has been directly estimated for only three systems, and that the most commonly constrained slope is α_out for the volumic density n and, to a lesser extent, γ_out for the stellocentric SB profile. In most studies, the comparison to the theoretical results of STCH and THWU is done by assuming that the profiles of τ and σ are identical (see previous paragraph) and that the indexes Γ, α_out, and γ_out are related through the relations $$\text{Γ}={\gamma }_{\text{out}}+2$$(2)

and $$\text{Γ}={\alpha }_{\text{out}}+\delta ,$$(3)

where δ is the index of the radial profile of the halo’s scale height, which is usually assumed to be equal to 1 (Thalmann et al. 2013; Ren et al. 2019). For the specific case of discs seen edge-on, we have Γ = γ_out + 1 + δ, which reduces to the same Γ = γ_out + 2 relation for δ = (Strubbe & Chiang 2006). However, we note that these relations are in principle only valid in an idealised case with isotropic scattering as well as for grains that are large enough for their contribution to the flux to be proportional to their geometrical cross section, and, in the case of edge-on discs, for an infinite instrumental resolution in the vertical direction.

All halos listed in Table 1 have been imaged in scattered light and, more generally, the halo phenomenon has so far mostly been investigated in the visible or near-IR and not in thermal emission. There are two main reasons for this. The first is theoretical: halos are supposed to be made of small grains not significantly bigger than the radiation-pressure blowout size s_blow, which are not expected to contribute at long wavelengths where they are poor emitters. The second reason is observational: high-end instruments in the visible or near-IR are those that offer the best spatial resolution by far, allowing the structures of halos to be resolved. However, the recent discovery of two extended halos around HD 32227 and HD 61005 with ALMA in the mm (MacGregor et al. 2018) has shed new light on the matter. Based on order-of-magnitude flux-to-mass conversions, MacGregor et al. (2018) estimated that these halos were mostly made of large mm-sized grains, which are not supposed to populate halos according to the canonical scenario of STCH and TBWU. To explain this puzzling presence, these authors explored some additional mechanisms, such as planetary or stellar perturbations, interaction with the ISM, or the aftermath of a large planetesimal breakup, but could not find a fully satisfying scenario. More recently, Olofsson et al. (2022a) revisited the modelling of the ALMA observations of HD 32297 (jointly with SPHERE polarimetric data) and found that one cannot rule out the possibility that these mm-halos are mostly made of smaller micron-sized grains after all. These authors showed that an overdensity of micron-grains, as expected because of their longer collisional lifetimes, could in principle compensate for their lower emissivity. However, using the simplified STCH and TBWU relations, they found that the flux due to these grains could only account for up to ~30% of the measured ALMA levels, but more sophisticated grain-distribution modelling might change these preliminary results and this crucial issue therefore remains an open question. Additional halo detections in the mm are expected in the near future, as ALMA is slowly catching up with the resolution of near-IR instruments such as SPHERE and GPI, and high-angular-resolution projects are currently being executed (e.g. ARKS large program, PI Marino). In addition, new ground- and spaced-based facilities (e.g. VLT/ERIS and JWST) can now observe debris discs in the mid-IR with unprecedented precision, which might potentially allow halos to be imaged at these wavelengths, underlining the need for a more comprehensive study of their detectability in thermal emission.

In light of these pending issues and recent new developments, we undertake a thorough numerical reinvestigation of the halo phenomenon. We focus on two main problems, the first of which is how the scattered light radial profiles derived by STCH and TBWU might be affected by taking into account the effect of small unbound grains, size-dependent scattering phase functions, and, for edge-on discs, instrument resolution in the vertical direction. Secondly, we explore how the belt+halo phenomenon manifests itself at longer wavelengths, in terms of both the system’s radial brightness profile as well as the halo’s imprint on the integrated spectral energy distribution (SED).

In Sect. 2, we briefly present the collisional evolution code that serves as a basis in our study as well as the typical belt+halo setup that we explore. Sections 3.1 and 3.2 present the results of our numerical exploration regarding radial profiles in scattered light. Results regarding profiles at longer wavelengths as well as system-integrated SEDs are presented in Sect. 3.3. We discuss the implications of our results in Sect. 4 and conclude in Sect. 5.

2 Model

2.1 Basic principle

Our numerical investigation uses the tried and tested particle-in-a-box collisional model initially created by Thébault et al. (2003), and constantly improved over the past two decades. As with similar codes, such as the ACE code of the Jena group (Krivov et al. 2005), particles are sorted into logarithmic size bins, whose populations are evolved following estimates of their mutual impact rates and collisional outcomes. This code has a 1D spatial resolution and is divided into radially concentric annuli. We use the latest version of the code, as described by Thebault & Kral (2019), for which collisional rates are estimated by separate deterministic N-body simulations taking into account the effect of stellar gravity and radiation pressure³, and with a more realistic collisional prescription for small particles in the ≲1 mm range taken from the laboratory experiments by Nagaoka et al. (2014). We refer the reader to Thébault & Augereau (2007) and Thebault & Kral (2019) for a detailed description of the code.

To derive SB profiles at different wavelengths, as well as SEDs, we use the (also tried and tested) GRaTeR radiative-transfer package (Augereau et al. 1999).

2.2 Setup

For the sake of readability of our results and in order to avoid numerically expensive parameter explorations, we consider a reference ‘nominal’ setup, chosen as to be the most representative of the observed belt+halo systems. The considered setup is that of a narrow belt of parent bodies, extending from 50 to 66 au, thus centered at r₀ = 58 au and of width ∆r₀ = 16 au, where all the mass is initially confined. We consider 119 log size bins, between s_max = 20 km and s_min = 0.025 µm. Large parent bodies in the belt (large enough for the radiation pressure effect to be negligible) are assumed to be located on orbits of average eccentricity <e₀>=0.05 and inclination <i₀>=<e₀>/2=0.025, which are typical values for debris-producing discs (e.g. Thébault 2009). Regarding grain composition, we consider the generic case of compact astrosilicates (Draine 2003). As for the central star, we consider an A6V stellar type, identical to the archetypal ß Pictoris case. For this stellar type, the blow-out size s_blow due to radiation pressure is of the order of 2 µm, which is much larger than our s_min, meaning that our code takes into account the potentially important effect of unbound grains.

Another important parameter is the level of collisional activity within the disc, which can be, to first order, parameterised by f_d, the disc’s fractional luminosity in IR. We consider a system with f_d = 8 × 10⁻⁴ as a reference case, which is approximately the average value for the resolved-halo systems presented in Table 1. In the spirit of Thebault & Kral (2019), we also consider a ‘very collisionally active’ case with f_d = 4 × 10⁻³, corresponding to the brightest discs in our census (such as HD 32297 or HD 129590). As explained in Thebault & Kral (2019), in practice we always start with discs whose initial masses are expected to correspond to an f_d that is larger than the ones we are aiming for, and then let the systems collisionally evolve until (1) the shape of the particle size distribution (PSD) no longer changes (steady state), and (2) f_d has decreased to the desired value. All the main parameters for this nominal setup are summarised in Table 2.

Fig. 1 Normalised radial profiles of the vertical optical depth (τ) and the SB, in scattered light (λ = 0.8 µm), for the nominal setup presented in Table 2, as well as for a ‘very bright’ disc case (fractional luminosity fd = 4 × 10−3). The blue area marks the radial extent of the parent body belt. Scattering is assumed to be isotropic here.

3 Results

3.1 Radial profile in scattered light

Figure 1 presents the radial profile, at steady state, of both τ and the SB in scattered light, for the nominal case as well as for the ‘bright disc’ case, assuming isotropic scattering when deriving the SB. Strictly speaking, the displayed SB(r) profiles correspond to a disc seen face-on (i = 0), but they should also be valid for any other viewing configuration as long as the disc is not seen edge-on, that is, as long as, at a given position along any radial cut of the observed disc, there are only grains from a given stellocentric distance that contribute to the flux. To a first approximation, this is true if the disc’s opening angle ψ is less than 90 − i. For these non-face-on cases, the displayed SB profiles are the ones that should be obtained for the deprojected stellocentric luminosity profile along any radial cut.

The average radial slope that is asymptotically reached by the optical depth profiles is −1.48 for the nominal case (f_d = 8 × 10⁻⁴), which is very close to the theoretical value of −1.5, and the average slope of the SB(r) profile is −3.42, again very close to the canonical value of −3.5 when applying the relation Γ = γ_out + 2. However, the optical depth profiles are shallower for the highly collisional bright-disc case (f_d = 4× 10⁻³), with a radial index of −1.21 for τ(r) instead of −1.48. This is mainly due to the fact that, for such high-f_d cases, unbound grains (s < s_blow) begin to significantly contribute to the geometrical cross section, the proportion of unbound grains being directly proportional to the level of collisional activity in the system (see Thebault & Kral 2019). The effect of unbound grains is less pronounced on the flux, because their smaller size, especially the ones in the ≲0.1 µm range, makes them less efficient scatterers. The slope of SB(r) in the very bright case is nevertheless still slightly smaller than for the f_d = 8 × 10⁻⁴ case, with a radial index of −3.25 instead of −3.42. In this f_d = 4 × 10⁻³ run, unbound grains even dominate the flux in the halo beyond ~200 au (Fig. 2). These s < s_blow dust particles will tend to flatten the radial profile because their distribution should follow τ ∝ r⁻¹ instead of ∝ r^−1.5 for bound grains (see Sect. 3.3 in STCH). This flattening effect is less visible for the reference f_d = 8 × 10⁻⁴ case, but the contribution of unbound grains to the flux still exceeds 10% everywhere in the halo (see black solid line in Fig. 2).

We note that, as already noted by Thebault et al. (2012), the asymptotic τ and SB slopes are not reached immediately beyond the PB belt, but after a transition region where the density and flux drop more abruptly. This is because, right outside the PB belt, there is a sudden transition from a ring where all particle sizes are present to an outer region only populated by small grains that are significantly affected by radiation pressure. If we define the limiting size s_PR for these latter grains by the criterion $$e\left(s_{\text{PR}}\right)=\frac{\beta \left(s_{\text{PR}}\right)}{1-\beta \left(s_{\text{PR}}\right)}=e_0,$$(4)

then, for e₀ = 0.05, we get s_PR ~ 11 s_blow. 11 s_blow. For a standard dN ∝ s^−3.5 size distribution, the absence of all s > s_PR grains results in a drop of geometrical cross section of ~30%, which is roughly what we observe in Fig. 1. The width of this transition region at the outer edge of the PB belt is ~0.15r₀ for both the nominal and high-f_d cases.

The radial SB slopes of Fig. 1 were derived by implicitly assuming isotropic scattering. Because the viewing angle should not vary along a given radial cut, they are in principle independent of the scattering phase function (SPF), provided that the SPF does not depend on stellocentric distance. This would, for instance, be the case for a standard Henyey–Greenstein SPF prescription (Henyey & Greenstein 1941) with a constant, size-averaged g parameter. Nevertheless, this assumption might not hold for halos, where the strong grain-size segregation as a function of r, coupled to the fact that the SPF is expected to vary with grain size, should result in a radial variation of the phase function. We therefore explore this effect further by including a more realistic SPF in the models. The simplest option would be to use Mie theory (Mie 1908), for which the dust particles are assumed to be compact spheres. However, several studies (e.g. Rodigas et al. 2015; Milli et al. 2019; Chen et al. 2020; Arriaga et al. 2020) have demonstrated that the assumption of spherical grains does not really hold when trying to reproduce observations of debris disks. Here, we instead use the SPF computed using the distribution of hollow spheres model (DHS, Min et al. 2005). For each grain size, the optical properties (e.g. scattering efficiencies, SPF) are obtained by averaging over a distribution of shapes. As discussed in Min et al. (2016), this model is able to reproduce the properties of irregularly shaped samples and the departure from spherical symmetry is controlled by the maximum filling factor 0 ≤ f_max < 1, which we set to 0.8 (Min et al. 2016). The SPFs are computed using the optool (Dominik et al. 2021) using the ‘DIANA’ standard opacities (Woitke et al. 2016, a mixture of pyroxene and carbon, with a mass ratio of 87% and 13% and a porosity of 25%, optical constants from Dorschner et al. 1995 and Zubko et al. 1996, respectively), at a wavelength of 0.8 µm.

As we can see in Fig. 3, the SB profile is significantly affected by the use of a more realistic SPF prescription. The most striking effect is to be found in the innermost part of the halo, just beyond the PB ring, where the luminosity drop is significantly reduced and even almost completely vanishes at low values of θ. In addition, the profile also gets slightly shallower further out in the halo. For the lowest θ = 15 deg run, the slope is ~−3.12 for the nominal f_d = 8 × 10⁻⁴ case (instead of −3.42 for isotropic scattering) and ~−2.95 for the bright f_d = 4 × 10⁻³ disc. This is due to two concurring effects. Firstly, because the average size of bound grains in the halo decreases with stellar distance, asymptotically tending towards s ~ s_blow (see Eq. (1)), and also because, at all scattering angles θ, the considered SPF increases with decreasing grain size in the s > s_blow domain (Fig. 4), the relative contribution of outer halo regions to the SB is increased. Secondly, for θ ≥ 4deg, the SPF actually peaks in the s < s_blow domain. This enhances the contribution of unbound grains to the flux (Fig. 2) which, as already noted, tends to further flatten the SB profile.

Fig. 2 Radial dependence of the fraction of the flux density at 0.8 µm in scattered light due to unbound grains (s < sblοw). Results are shown for the nominal and ‘bright disc’ cases, as well as for anisotropic scattering, at two different angles, using the distribution of hollow spheres model (DHS, Min et al. 2005) for the SPF.

Fig. 3 Radial profile of the normalised SB obtained using the DHS scattering phase function prescription for three different scattering angles.

Fig. 4 Size dependence of the DHS scattering phase function for four different scattering angles.

3.2 Edge-on configuration

The edge-on case introduces additional specific issues that do not affect the general configuration presented in the previous section. The first one concerns the scattering phase function. The SB is indeed now measured in the disc’s midplane along the projected distance ρ, and is therefore now the sum of contributions coming from different physical stellocentric distances r, for which the scattering angle is not the same. As a result, the resulting luminosity should in principle depend on the SPF even for SPF prescriptions that do not depend on grain size, and hence stellar distance. We test the importance of this effect by considering the classical Henyey-Greenstein prescription, for which the anisotropy of the scattering behaviour of the dust grains is characterised by the dimensionless asymmetry parameter −1 < g < 1: $$f\left(\theta \right)=\frac{1-g^2}{4\pi {\left(1+g^2-2g\cos \left(\theta \right)\right)}^{3/2}}.$$(5)

The g = 0 case represents isotropic scattering, in which photons are scattered in all directions with equal probabilities. For positive (negative) values of g, incident photons are scattered in the forward (backward) direction, and as g increases (decreases), the asymmetry is even more pronounced. In practice, the value of the asymmetry parameter g can be used to control the size of the particles, as grains much smaller than the wavelength of observations are expected to scatter isotropically (hence g = 0), while grains larger than the wavelength should display a strong forward-scattering peak (g > 0).

Figure 5 presents the projected midplane profiles SB(ρ) for the isotropic (g = 0) case as well as for three different values of g, assuming that the disc is not flared and has a constant aspect ratio h = H/r. For the isotropic scattering case, we get SB(ρ) ∝ r^−3.3 in the outer regions, which is close to the standard value of −3.5 that directly follows from the aforementioned relation Γ = γ_out + 1 + δ for a constant h (i.e. δ = 1). We see that the projected SB(ρ) profile strongly depends on g in the innermost ρ < r₀ regions. This is expected, because the scattering angle for the most luminous grains, that is, those in the PB ring, can reach very small values at small ρ, which leads, for strong forward scattering (i.e. high g values), to a significant increase in SB(ρ) for decreasing ρ. In the outer regions beyond the projected outer edge of the PB belt, the dependence on g is weaker but still noticeable. For the most extreme g = 0.95 case, the radial index reaches γ_out ~ −3.7 in the outer regions, which is ~−0.4 steeper than for the isotropic (g = 0) case. This weaker dependence is due to the fact that, in these outer regions, the variation of scattering angle with ρ is more limited than in low ρ regions. Interestingly, taking a more realistic and size-dependent DHS prescription for the SPF tends to move the SB profile back close to the reference r^−3.5 slope (blue line in Fig. 5). This is because of the increased role of small bound and unbound grains, which act to flatten the SB profile (see previous section), and thus act in the opposite direction of the geometrical effect of integrating along the line of sight.

Another potentially important issue for the edge-on configuration is whether or not the disc is resolved in the vertical direction. The midplane profiles presented in Fig. 5 indeed implicitly assume that the disc is resolved in z at all projected distances ρ, so that there is a natural geometrical dilution of the flux in the vertical direction. However, vertical resolution is not necessarily achieved with current instrument facilities (see Sect. 4). We explore the effect of not resolving the disk vertically in Fig. 6 by considering four different cases, ranging from fully resolved everywhere (resolution res ≪ h × r₀) to a very poorly resolved case where only the outermost >4r₀ regions are vertically resolved. As can clearly be seen, the SB(ρ) midplane profile flattens with decreasing vertical resolution. For the lowest resolution (res = 4 × h × r₀) case, the slope for the midplane profile tends towards γ_out = −2.3. This value is ~γ_out(res=0) + 1, which is expected because the vertical dilution term is now absent and the ‘midplane’ luminosity corresponds here to the z-integrated flux.

Fig. 5 Edge-on disc. Projected radial profile of the midplane SB for four different values of the g parameter of the Henyey–Greenstein phase function, as well as for the DHS prescription for the SPF (see main body of the text).

Fig. 6 Edge-on disc. Projected radial profile of the midplane SB for four different values of instrument resolution. Here, h is the disc’s aspect ratio, r0 the centre of the parent body belt, and dr/2 its width (isotropic scattering is assumed).

Fig. 7 Radial profile of the normalised SB at four different wavelengths, estimated with the GRaTeR package for the nominal setup.

3.3 Wavelength dependence and SED

We present here, for the first time, an exploration of the halo phenomenon over a wide range of wavelengths, notably in thermal emission up to the mm domain.

Figure 7 presents, for the non-edge-on configuration, how the SB(r) profiles vary with λ. At long wavelengths (160 and 800 µm), there is a sharp drop at the outer edge of the parent-body belt. This is due to the fact that, at these large λ, the absorption coefficient Q_abs of the small grains that dominate the geometrical cross section is much smaller than 1. As a consequence, the flux within the PB belt is dominated by large s ≫ s_blow grains, whose number density abruptly drops beyond the PB belt because their orbits are only very weakly affected by stellar radiation pressure. Beyond r₀ + ∆r₀/2, small grains take over, but their higher number density cannot compensate for their very low Q_abs ≪ 1 values. We note that, in this dimmer far-IR and mm halo, the radial profile index of the SB is ~−2.2, which is significantly shallower than in scattered light. The likely interpretation for this is that, at these long wavelengths, the temperature Τ of grains in the ~70–400 au regions implies that they emit very close to the Rayleigh-Jeans approximation of the Planck function, for which the flux at a given wavelength λ is ∝T. As the blackbody temperature of grains is in turn proportional to r^−0.5, it follows that the relation between the slope Γ of the vertical optical depth τ and the slope γ_out of the SB(r) profile is Γ = γ_out + 0.5 instead of Γ = γ_out + 2 (see Eq. (2)).

However, at λ = 70 µm, the drop at the outer edge of the PB belt is much more limited, and is only ~30% more pronounced than in scattered light. This is because, at this wavelength, the poorly emitting (i.e. with Q_abs < 1) grains come from a narrow size range of s_blow ≤ s ≤ λ/2π ~ 10 µm, which only accounts for ~60% of the total cross section in the birth ring⁴. Moreover, even grains in this s_blow ≤ s ≤ λ/2π domain still have non-negligible Q_abs values at λ = 70 µm (typically ~0.2–0.3, e.g. Morales et al. 2013). The slope of the profile is ~−2.9, which is slightly steeper than at longer wavelengths; this results from the fact that the Rayleigh-Jeans approximation becomes less accurate at these shorter λ.

The SB profile at λ = 15 µm is very different, with an absence of a luminosity drop at the edge of the PB belt, followed by a very steep decrease with radial distance in the halo. This is because, at this wavelength and at r ≥ 60 AU from an A6V star, we are in the Wien side of the Planck function, for which the flux increases exponentially with T. As a consequence, small grains (bound and unbound), whose temperatures exceed the almost black-body temperature of larger particles, totally dominate the flux⁵, and there is no drop at the outer edge of the PB ring due to the sudden absence of large grains. However, because of the exponential dependence of the flux on T in the Wien domain, the decrease in SB with radial distance becomes steeper and steeper as the grains get colder in the outer halo, going from a radial index of ~−4.3 just outside the PB belt to ~−6 in the 300–400 au region.

In Fig. 8 we look at the wavelength dependence from the perspective of the system-integrated SED. We see that, despite corresponding to grains that are further away from the star, the SED of the halo actually peaks at a shorter wavelength than the contribution from the PB belt. In addition, at all wavelengths shorter than ~90 µm, the relative contribution from the halo to the total flux F_halo/F_disc exceeds 50% and thus dominates that of the PB ring (Fig. 9). These results agree well with the radial profiles, showing that there is no sharp luminosity drop at the PB belt–halo interface for wavelengths short of 70 µm. The only exception is a narrow wavelength range around λ ~ 8–15 µm. This corresponds to a ‘sweet spot’ where the thermal flux dominates scattered light but is in the Wien regime of the Planck function, for which there is a very steep decrease in the flux with radial distance (see above) and the contribution from the halo is therefore much lower. However, we note that the λ ~ 10 µm domain corresponds to wavelengths at which the disc is very faint (Fig. 8). At longer wavelengths (≥90 µm), we logically see a decrease of the halo contribution to the total flux, which is the direct consequence of the low emissivity of its small-grain population at increasing λ. Nevertheless, even in the mm-wavelength domain, the contribution from the halo is never totally negligible. As an example, it is still 5% of the total flux at λ = 800 µm.

From Fig. 9, we also see that, when considering the global F_halo/F_disc ratio, we obtain relatively similar results for both the nominal and bright-disc cases. This is because, in terms of the relative contribution of the halo, the only parameter distinguishing these two cases (once normalised) is the fraction of unbound grains in the system. The higher fraction of submicron grains for f_d = 4 × 10⁻³ will only affect the flux at short wavelengths (λ ≲ 1 µm), where we see that F_halo/F_disc is indeed ~10% higher for this bright-disc case, as well as in the aforementioned λ ~ 8–15 µm domain, where the system’s total flux is dominated by unbound grains within the PB belt (see Thebault & Kral 2019).

Fig. 8 Normalised system-integrated SED for the nominal setup (fd = 8 × 10−4), displaying also the respective contributions coming from the parent body belt (between 50 and 66 au) and the halo (beyond 66 au).

Fig. 9 Relative contribution to the total flux as a function of wavelength coming from the whole halo for both the nominal and ‘very bright’ disc cases. The discontinuity at λ ~ 8 µm corresponds to the transition between the scattered-light-dominated domain and the thermal-emission-dominated one.

4 Discussion

4.1 How universal are the SB ∝ r^−3.5 and τ ∝ r^−1.5 profiles?

As mentioned in Sect. 1, halo radial profiles are sometimes used as a proxy to constrain the level of ‘unexpected activity’ in the outer regions of debris discs: perturbations by (unseen) planets, the effect of companion or passing stars, gas drag, and so on. This is usually done by measuring the SB profiles, or, more often, the extrapolated underlying dust density distribution is compared to the expected ‘normal’ profiles for unperturbed systems derived by Strubbe & Chiang (2006) and Thébault & Wu (2008). However, such direct comparisons raise several issues. The first is that the STCH and TBWU reference radial slopes are only valid under several simplifying assumptions, the main ones being the absence of unbound grains and, for the SB profile, isotropic scattering. Conversely, most observation-based fits of the σ and n profiles also make some strong simplifications, notably that the grain-size distribution is the same everywhere in the system, which is clearly not the case for halos that are on the contrary extremely size-segregated. While this simplification is of limited consequence as long as isotropic scattering is assumed⁶, it becomes problematic when considering more realistic and, in particular, size-dependent SPFs.

Here, we have, for the first time, numerically explored how these different approximations and simplifications could bias our understanding of debris-disc halos. One important result is in regards to the proportion of unbound grains, which always account for more than 10% of the scattered-light luminosity and even dominate the flux in the outer-halo regions for a bright, collisionally active disc (Fig. 2). For this f_d = 4 × 10⁻³ case, the presence of unbound grains is able to significantly flatten both the τ and SB profiles, as these small particles have a shallower radial distribution in τ ∝ r⁻¹. The effect of small grains becomes even more pronounced when considering a realistic size-dependent SPF function, which, for all scattering angles θ ≥ 4° always peaks in the s < s_blow domain. To a first approximation, the combined effect of high collisional activity and size-dependent SPF results in SB slopes that tend towards ~−3 instead of the standard −3.5 value. Conversely, we expect GRaTeR-type fits of the underlying dust density – which do not take into account these effects – to underestimate the index of the σ or n slopes by up to a value of ~0.5.

For the specific case of edge-on discs, this flattening effect, which is due to the increased influence of unbound grains and of the smallest bound ones, is less visible for the midplane SB profile. This is mainly because it is in large part compensated for by the purely geometrical effect that comes from the fact that the flux at a projected distance ρ is now the sum of contributions integrated along the line of sight, for which the flux dependence on ρ becomes weaker with increasing projected distance. However, there is a parameter that potentially has a much greater influence on the midplane profile of edge-on discs, which is the non-resolution of the disc in the vertical direction. Our results indeed show that the consequence of not vertically resolving a disc could lead to a flattening of up to +1 in terms of radial index (with an index of −2.5 instead of −3.5) of the SB radial profile of its halo. This result might prove important because, even if second-generation instruments, such as SPHERE or GPI, provide a pixel scale of about 12–15 milli-arcsec (e.g. Maire et al. 2016), and the angular resolution of ALMA observations continues to improve – down to the au scale for the closest systems –, only a handful of systems, and the closest ones, have been resolved in the vertical direction z. From Table 3 and Fig. 11 of Olofsson et al. (2022b), we see that, amongst the list of constrained-halo systems of Table 1, only β-Pic, Au Mic, HR4796, HD 115600, and HD 61005 have been unambiguously resolved in z at near-IR wavelengths.

With these news results in mind, we can take a renewed look at Table 1. Our analysis has shown that the most reliable halo slope estimates are likely to be those made directly on the observed SB profiles, because they do rely on fewer model-dependent assumptions (regarding the presence of unbound grains or the size-dependence of the SPF). As a consequence, we show in Fig. 10 the dispersion of radial slopes for all systems for which it is the SB profile that has been fitted from observations. It can be seen that 6 out of the 13 considered systems have halos with profiles that are fully within the boundaries of acceptable values found by our numerical investigation; namely Fomalhaut, AU Mic, HD 15745, HD 32297, HD 35841, and HD 53143. On the opposite side of the spectrum, three discs, namely HR 4796, HD 107146, and HD 36546⁷, fall fully outside these boundaries and are therefore likely to be systems for which additional mechanisms are sculpting the outer realms of the discs. For the remaining four systems, some parts of the halo have expected radial slopes while some do not, and an assessment of these systems is not possible without undertaking system-specific investigations that go beyond the scope of the present paper. However, we note that for two of these systems, HD 191089 and HD 202628, the steeper slopes have been measured in a narrow region just outside the PB belt (see the difference between the r₀ and r_max values in Table 1), for which our simulations have shown that a steeper SB profile is in fact expected.

As discussed in Sect. 3, fits of the underlying n and σ profiles obtained with GRaTeR or similar codes should be less reliable than fits of the SB, and also more challenging to interpret in terms of the physical meaning of the fitted quantities, especially when considering the fact that there is a strong size segregation in the halo. Nevertheless, we can consider that the radial dependence of the estimated σ and n profiles⁸ should, to a first approximation, give a rough estimate of the radial dependence of the vertical optical depth in the halos of these system, which can be compared to our results regarding the τ profiles. Making this approximation of an equivalence between the σ (or n × r) and τ profiles, we see in Fig. 11 that a little less than half of the systems have radial density profiles that are within ±50% of the reference τ profiles obtained in Fig. 1. The spread of radial indexes is larger than for the fitted SB slopes displayed in Fig. 10, with ~25% of systems having radial slopes that are more than two indexes below that of our reference simulations. However, we note that for the three systems with the steepest estimated slopes, that is, HD 172555, HD 160305, and HD 115600, the halo profiles are only derived for a relatively limited radial region. As already pointed out, for this narrow region just beyond the PB belt, both the SB and τ profiles are expected to be steeper, meaning that the slope indexes in this region should not be representative of the halo profile further out. In addition, for some systems where both the SB and n or σ radial behaviour have been fitted, we see some incoherent results, with for instance the σ profiles being steeper than the profile of the SB for HD 15115 and HD 32297. This points towards an intrinsic potential pitfall in global ring+halo fits, which we discuss in the following subsection.

Fig. 10 Surface-brightness radial slopes taken from Table 1. The dark blue area corresponds to the expected values in a non-perturbed system according to the present numerical investigation. The light-blue area is the same, but for edge-on discs, taking into account the potential nonresolution of the disc in the vertical direction. For some systems, there are different slope estimates depending on the radial position in the halo and on which side of the disc has been considered. In this case, up to four values are displayed: diamonds and cross symbol signify the radial indexes in the inner- and outer halo, respectively, for one disc side, and the squares and plus symbol are the equivalent indexes for the opposite side. For systems where there is only one global fitted radial index, all four symbols do overlap. Edge-on systems are written in italics.

Fig. 11 Values of the radial profile indexes of vertical optical depth (τ) derived from the n and σ fits displayed in Table 1, when making the simplifying assumption that the radial dependence of τ is the same as that of σ, or is equal to that of n plus one (constant opening angle). The blue domain corresponds to the range of τ(r) indexes between the −1.48 value obtained for our nominal case and the −1.21 value for our bright disc (fd = 4 × 10−3) case (see Fig. 1).

4.2 Suggested procedure for fitting discs with halos

Studies deriving σ or n profiles from observations treat the ring+halo system as a whole, usually assuming that the density profile follows $$n\left(r\right)={\left[{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{in}}+{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{\text{out}}}\right]}^{-1/2}{n}_0,$$(6)

and exploring different values of r₀, α_in, and α_out, as well as of the disc inclination i and, often, the g parameter of the HG scattering phase function. The best fit is then found by comparison to the observed luminosity profile through a classical χ² minimisation procedure. A potential problem is that, because the flux in the bright PB belt is usually much higher than in the halo, the χ² fit is dominated by the narrow bright ring, which means that acceptable global fits can actually be a poor match of the luminosity profile in the halo region. This could explain the finding that, in some cases, fitted n profiles seem to disagree with the observed SB profile in the halo (see previous section). More generally, this could lead to large errors when trying to assess the specific structure of the halo.

To alleviate these potential problems, we suggest that, in future observational studies, the fitting of the radial profile of the system should be done in two steps, in which the main ring and the halo are fitted separately. For non-edge-on cases, the PB belt alone could first be fitted by finding the best possible set of i, r₀, and ∆r₀ values, where ∆r₀ is the FWHM of the ring. Once these parameters are constrained, the halo profile beyond r₀ + 0.5∆r₀ can then be investigated by fitting the α_out index. Such a procedure should ideally take into account the varying size distribution within the halo, the potential effect of unbound grains, and the size-dependence of the SPF, all parameters that could potentially change the correspondence between an observed SB profile and the radial profile of the underlying optical depth distribution. Nevertheless, thoroughly exploring these parameters would probably render the fitting procedure much too cumbersome. A possible compromise could be the ‘semi-dynamical’ model that Pawellek et al. (2019) experimented, which takes a PB belt as constrained from observations in the submm as a starting point, from which smaller grains are produced with their abundances scaled up by the corrective factor found by STCH and TBWU. Synthetic images are then produced, potentially taking into account realistic SPFs, which are compared to observations in scattered light. Nevertheless, this is not a fit per se, although it provides important information as to whether or not the observed ring+halo system behaves according to the predictions for unperturbed systems found by STCH and TBWU. Here, we propose to approach the problem the other way around, taking the constrained SB(r) profile in the halo as a starting point, and ‘reverse engineer’ the τ(r) profile from it. Of course, as the SPF in the halo depends on the grain-size distribution, which in turn depends on r and is not known beforehand, there is in principle one unknown too many. However, an approximation of the SPF(r) dependence could be obtained assuming that, at any given position r in the halo, the geometrical cross section should be dominated by grains of a size s, such as their apoastron is located at r when produced at r₀ in the main PB belt (Eq. (1)). The SPF at position r can then be estimated for this specific grain size s and incorporated into the α_out fitting procedure.

We note that, in principle, sophisticated numerical codes such as ACE, LIDT-DD (Kral et al. 2013), or even the collisional model presented in Sect. 2 of the present paper, can provide self-consistent estimates of the particle size distribution as a function of radial location in the system that are more accurate than what would be obtained by the procedure suggested in the present paragraph. However, using such sophisticated models to do parameter best fits of observed discs would require performing a huge set of CPU-consuming simulations exploring an extended parameter space, which would imply a considerable amount of effort. While this is undertaken in some rare cases, such as Müller et al. (2010) for Vega or Löhne et al. (2012) for HD 207129, it can only be done within the framework of long and especially dedicated numerical studies. What we are aiming for here is different, that is, a relatively easy and approximate but efficient way of constraining main disc parameters without the inherent flaws of only using Eq. (6) coupled to a radiationtransfer code. The procedure we are presenting is typically to be used at the end of observational studies presenting new resolved data.

4.3 Halos in thermal emission

Our numerical exploration has shown that halos, despite being made of small micron-sized grains, should remain relatively bright deep into the mid-IR and even far-IR domain. As an example, for a typical debris ring located at 60 au, the halo to PB belt brightness ratio at 70 µm is still close to what it is in scattered light (Fig. 7). It is only at wavelengths longer than ~100 µm that the halo luminosity drops. We also note that, at these long wavelengths, the radial profile (beyond the sharp drop at the outer edge of the PB belt) of the SB in the halo is significantly shallower than in scattered light (∝r^−2.2 instead of ∝r^−3.5). Another important result of the present study is in regards to halo total fluxes. We show that, except for a narrow domain around λ ~ 10 µm, halos contribute between 50% and 70% of the system’s total flux at all wavelengths short of ~90 µm, and even in the mm-domain the halo still has an integrated luminosity that amounts to a few percent of that of the main parent belt. Perhaps more tellingly, we find that the total, wavelength-integrated thermal emission of the halo is approximately half that of the whole system.

Our results show that, in thermal emission, the optimal wavelength window for observing halos beyond belts in the 50–70 au region is λ ~ 20–70 µm. This domain covers two bands (24 and 70 µm) of the MIPS instrument on the Spitzer telescope as well the 70 µm band of the PACS instrument on the Herschel telescope. While the resolution and sensitivity of Spitzer were probably too limited to explore halos, the situation was more favourable for Herschel-PACS at 70 µm. However, most PACS-resolved discs had estimated sizes smaller than the FWHM of the instrument (5.6″ at 70 µm) and were inferred from image deconvolution (e.g. Booth et al. 2013), preventing from deriving reliable radial profiles. Such radial profiles were obtained for only very few systems, such as Vega (Sibthorpe et al. 2010), ϵ Eridani (Greaves et al. 2014), Fomalhaut (Acke et al. 2012), and HD 207129 (Löhne et al. 2012), but without constraining the slopes of the SB profile in the outer regions. A reinvestigation of the Herschel data for these systems, and in particular an estimate of their SB(r) radial slopes, would definitely improve our understanding of the halo phenomenon by comparison with our present results. However, in most cases, the absence of radial profiles means that, if halos were present, it would not be possible to identify them as a fading extension of a bright belt, but they would nevertheless contribute to the estimated disc size. It is therefore likely that, at 70 µm, a significant fraction of PACS-derived disc radii correspond to a blend between the main collisional belt and the halo, and cannot be reliably used to trace the location of the dust mass reservoirs of these systems. This blending effect would be less important for the 100 µm PACS band, but the resolution here is poorer (6.8″) and our results show that, even at this wavelength, the halo still contributes ~40% of the system’s flux (see Fig. 9). New generations of far-IR instruments would be crucially needed here to untangle the belt and halo contributions by providing reliable radial SB profiles in the 20–70 µm domain. We note that the longest wavelength of observation of the JWST, of namely 28.3 µm, overlaps with the low end of our optimal window for halo observations. Given the unparalleled sensitivity of that instrument, we therefore expect it to provide us with the first resolved images of debris-disc halos in the mid-IR.

Our results also challenge the way disc SEDs are sometimes used to constrain the global particle size distribution (PSD) in resolved systems. The procedure to constrain the PSD is indeed usually to consider the geometrical profile constrained from image fitting and then find the power law index q of the PSD that best fits the SED, under the assumption that the dN ∝ s^qds PSD holds everywhere in the system (e.g. Pawellek et al. 2019). However, this procedure is not adapted to systems for which a large fraction of the SED is due to the halo, which has a PSD that is very different from that of the PB belt. Moreover, the very notion of a single dN ∝ s^qds law for the halo makes little physical sense, because this region is strongly size-segregated. The PSD for the halo should therefore in principle be derived as a function of radial location, but this would imply having reliable estimates of the SED at different locations in the outer regions, which is generally impossible to obtain. A possible intermediate solution would be, as for estimating the SPF in Sect. 4.2, to assume a simplified size-distribution in the halo, where a given radial location r is only populated by mono-sized grains produced in the PB belt and having their apoastron at r (Eq. (1)). With this assumption, and for systems for which the halo’s geometry and SB profile have been constrained from image-fitting by the aforementioned procedure, its contribution to the SED can be unequivocally estimated. This can then be subtracted from the total SED to allow estimation of the PSD index q in the main belt through the usual procedure. Of course, our simplifying assumption for estimating the SED of the halo would imply that its SB is close to ∝r^−3.5 and should in principle not be used for systems where the SB’s profile strongly departs from this radial dependence. However, we believe that, in any case, this procedure is preferable to a global fit that assumes a uniform PSD everywhere in the system. As for the procedure outlined at the end of Sect. 4.2, sophisticated codes such as ACE or LIDT-DD could be in principle used to obtain more accurate estimates of the SED. However, here again, what we are aiming for is a quick and straightforward procedure that can be used in observation-based studies without the flaws that come with the assumption of a constant PSD everywhere in the system.

The detailed study of individual systems goes beyond the scope of the present paper, but we conclude this section by briefly discussing whether or not the two halo detections obtained with ALMA for HD 32297 and HD 61005 (MacGregor et al. 2018) can be explained by the ‘natural’ behaviour of halos at long wavelengths without invoking additional mechanisms. We first note that, contrary to our numerical results, the radial profiles in the mm of these two halos appear relatively steep, with −6.2 and −5.5 for the extrapolated index Γ of the surface density instead of −1.5. In addition, the radially integrated deprojected luminosities F_halo of these halos amount to between 20% and 30% that of the parent-body belt F_belt, which is one order of magnitude more than for our synthetic halos at λ = 1.3 µm (see Fig. 9). This points towards an ‘abnormal’ halo and therefore the need for additional mechanisms at play in the outer regions of these systems⁹. However, we note that the fitting procedure adopted by MacGregor et al. (2018) imposes a density continuity at the belt–halo interface and it is not clear to what extent this assumption, coupled to the almost edge-on orientation of the system, affects the obtained results in terms of density slopes and F_halo/F_belt ratios.

4.4 Unresolved systems

Our results also have consequences for systems that have not been resolved at any wavelength and whose size and radial location r_d are only constrained by their SED. One of the most sophisticated methods to retrieve r_d from the SED is the one proposed by Pawellek et al. (2014) and Pawellek & Krivov (2015). The first step of this procedure is to fit the SED with a modified black body (MBB) model (Backman & Paresce 1993) in order to derive the typical dust temperature T_dust, which is in turn used to derive a blackbody radius r_BB. The ‘true’ disc radius r_d, which implicitly corresponds to that of the PB belt, is then obtained by multiplying r_BB by a factor ${\text{Γ}}_{\left(L_{*}\right)}$ that depends on stellar type. The ${\text{Γ}}_{\left(L_{*}\right)}$ ratio is obtained separately by empirically comparing r_BB to real disc radii for a sample of resolved systems in the far-IR with Herschel. To their credit, Pawellek et al. (2014) were aware of the risk of incorrectly estimating r_d if considering wavelengths at which small, radiation-pressure grains can contribute, and thus chose disc sizes retrieved from Herschel PACS images at λ = 100 µm instead of 70 µm. However, as already mentioned, our results show that, even at this relatively long wavelength, the halo of small grains still makes up ~40% of the system’s total flux¹⁰ and could therefore lead to overestimation of r_d. A more reliable estimate would be obtained by considering ALMA images in the mm domain, but such data were not readily available at the time the ${\text{Γ}}_{\left(L_{*}\right)}$-based procedures were developed. We therefore strongly recommend updating the ${\text{Γ}}_{\left(L_{*}\right)}$ empirical law, taking r_d values determined from ALMA images as a reference, in the spirit of the study by Matrà et al. (2018) or Pawellek et al. (2021).

We note that our results challenge the notion that a single MBB law can accurately fit a system that is actually made of two distinct components, PB belt and halo, which have very different spatial structures and particle size distributions, and whose SEDs peak at different wavelengths. As a consequence, the ‘disc temperature’ T_d derived by the MBB procedure probably has a limited physical meaning and cannot be a reliable estimate of the temperature in the PB belt. As the halo’s SED peaks at a shorter wavelength than that of the PB belt, T_d probably overestimates the actual temperature of the collisionally active region of the disc. However, this does not invalidate the Γ-based procedure, as it is in principle independent of whether or not T_d has a physical meaning. Of importance here is the reliability of the empirical ${\text{Γ}}_{\left(L_{*}\right)}$ ratio estimates, for which T_d (or rather r_BB) can be considered an abstract proxy.

5 Summary and conclusion

We carried out the most thorough investigation of the halo phenomenon to date. We focus in particular on two issues: (1) the robustness of the theoretical τ ∝ r^−1.5 and SB ∝ r^−3.5 radial profiles when taking into account the role of unbound grains, realistic SPF prescriptions, and instrument resolution, and (2) the behaviour of halos in thermal emission out to the mm domain, both on resolved images and on system-integrated SEDs.

For a typical halo produced beyond a collisional belt located at ~60 au, our main results can be summarised as follows:

• –

  The contribution of small unbound grains amounts to at least ~10% of halo luminosities in scattered light, and can even dominate in the outer halo regions for bright discs. For these brightest discs, halo radial profiles can become significantly shallower;

• –

  Size-dependent scattering phase functions (SPFs) also result in flatter radial profiles, which directly follows from the fact that halos are strongly size-segregated regions;

• –

  For systems viewed edge-on, not resolving the disc in the vertical direction can lead to a flattening of the SB(ρ) radial profile index by up to one;

• –

  Comparing these new results to a complete sample of obser-vationally constrained halo SB profiles, we find that roughly half of them have radial profiles that are fully compatible with our predictions, while ~25% have profiles that cannot be explained by our models (usually being too steep). For these systems, additional mechanisms should be at play to shape the outer regions. For a large fraction of the remaining ~25%, halo profiles have been derived in a region of insufficient width to allow definitive conclusions to be reached;

• –

  We obtain comparable results for systems for which it is the underlying dust density distribution whose radial profile has been observationally constrained. However, these density-distribution fits should be less reliable than SB ones. In particular, they are likely to be biased by the fact that they do not discriminate between belt and halo, and that they do not take into account the effect of unbound grains and size-dependent SPFs;

• –

  We suggest that future observational fits of the underlying density distribution in systems with halos should be made in two steps, starting with a geometrical fit of the PB belt, whose parameters are then injected into a separate fit of the halo’s radial slope accounting for size-dependent SPF effects;

• –

  Radially extended halos should also be visible in thermal emission in the λ ~ 20–100 µm range, where the halo-to-main-belt contrast is comparable to what it is in scattered light;

• –

  With the exception of a narrow 8 ≲ λ ≲ 15 µm domain, halos always account for more than 50% of the disc’s total flux up to λ ~ 90 µm;

• –

  Despite being located further out than the PB belt, the halo’s SED peaks at a shorter wavelength than that of the belt, which makes the global system appear hotter;

• –

  Beyond λ ~ 90 µm, halo brightness strongly decreases with wavelength, but halos still contribute to a few percent of the flux in the mm domain. However, this seems to be insufficient to explain the bright halos detected with ALMA around HD 32297 and HD 61005;

• –

  For unresolved discs, the presence of a halo can also bias the procedure inferring their radial location from an analysis of their SED.

The study of individual systems goes beyond the scope of the present paper and is deferred to future studies. These studies will explore additional parameters, such as stellar type, PB belt radial location, and dynamical context (such as known stellar or planetary companions).

Acknowledgements

J.O. acknowledges support by ANID, – Millennium Science Initiative Program – NCN19_171. The authors thank Anthony Boccaletti for fruitful discussions.

References

All Tables

All Figures

	Fig. 1 Normalised radial profiles of the vertical optical depth (τ) and the SB, in scattered light (λ = 0.8 µm), for the nominal setup presented in Table 2, as well as for a ‘very bright’ disc case (fractional luminosity fd = 4 × 10−3). The blue area marks the radial extent of the parent body belt. Scattering is assumed to be isotropic here.
In the text

	Fig. 2 Radial dependence of the fraction of the flux density at 0.8 µm in scattered light due to unbound grains (s < sblοw). Results are shown for the nominal and ‘bright disc’ cases, as well as for anisotropic scattering, at two different angles, using the distribution of hollow spheres model (DHS, Min et al. 2005) for the SPF.
In the text

	Fig. 3 Radial profile of the normalised SB obtained using the DHS scattering phase function prescription for three different scattering angles.
In the text

	Fig. 4 Size dependence of the DHS scattering phase function for four different scattering angles.
In the text

	Fig. 5 Edge-on disc. Projected radial profile of the midplane SB for four different values of the g parameter of the Henyey–Greenstein phase function, as well as for the DHS prescription for the SPF (see main body of the text).
In the text

	Fig. 6 Edge-on disc. Projected radial profile of the midplane SB for four different values of instrument resolution. Here, h is the disc’s aspect ratio, r0 the centre of the parent body belt, and dr/2 its width (isotropic scattering is assumed).
In the text

	Fig. 7 Radial profile of the normalised SB at four different wavelengths, estimated with the GRaTeR package for the nominal setup.
In the text

	Fig. 8 Normalised system-integrated SED for the nominal setup (fd = 8 × 10−4), displaying also the respective contributions coming from the parent body belt (between 50 and 66 au) and the halo (beyond 66 au).
In the text

	Fig. 9 Relative contribution to the total flux as a function of wavelength coming from the whole halo for both the nominal and ‘very bright’ disc cases. The discontinuity at λ ~ 8 µm corresponds to the transition between the scattered-light-dominated domain and the thermal-emission-dominated one.
In the text

Fig. 10 Surface-brightness radial slopes taken from Table 1. The dark blue area corresponds to the expected values in a non-perturbed system according to the present numerical investigation. The light-blue area is the same, but for edge-on discs, taking into account the potential nonresolution of the disc in the vertical direction. For some systems, there are different slope estimates depending on the radial position in the halo and on which side of the disc has been considered. In this case, up to four values are displayed: diamonds and cross symbol signify the radial indexes in the inner- and outer halo, respectively, for one disc side, and the squares and plus symbol are the equivalent indexes for the opposite side. For systems where there is only one global fitted radial index, all four symbols do overlap. Edge-on systems are written in italics.

In the text

Fig. 11 Values of the radial profile indexes of vertical optical depth (τ) derived from the n and σ fits displayed in Table 1, when making the simplifying assumption that the radial dependence of τ is the same as that of σ, or is equal to that of n plus one (constant opening angle). The blue domain corresponds to the range of τ(r) indexes between the −1.48 value obtained for our nominal case and the −1.21 value for our bright disc (fd = 4 × 10−3) case (see Fig. 1).

In the text