Impact of internal-delensing biases on searches for primordial B-modes of CMB polarisation

Searches for the imprint of primordial gravitational waves in degree-scale CMB $B$-mode polarisation data must account for significant contamination from gravitational lensing. Fortunately, the lensing effects can be partially removed by combining high-resolution $E$-mode measurements with an estimate of the projected matter distribution. In the near future, experimental characteristics will be such that the latter can be reconstructed internally with high fidelity from the observed CMB, with the $EB$ quadratic estimator providing a large fraction of the signal-to-noise. It is a well-known phenomenon in this context that any overlap in modes between the $B$-field to be delensed and the $B$-field from which the reconstruction is derived leads to a suppression of delensed power going beyond that which can be attributed to a mitigation of the lensing effects. More importantly, the variance associated with this spectrum is also reduced, posing the question of whether the additional power suppression could help better constrain the tensor-to-scalar ratio, $r$. In this paper, we show this is not the case, as suggested but not quantified in previous work. We develop an analytic model for the biased delensed $B$-mode angular power spectrum, which suggests a simple renormalisation prescription to avoid bias on the inferred tensor-to-scalar ratio. With this approach, we learn that the bias necessarily leads to a degradation of the signal-to-noise on a primordial component compared to"unbiased delensing". Next, we assess the impact of removing from the lensing reconstruction any overlapping $B$-modes on our ability to constrain $r$, showing that it is in general advantageous to do this rather than modeling or renormalising the bias. Finally, we verify these results within a maximum-likelihood inference framework applied to simulations.


Introduction
The theoretical framework best fitting cosmological observations includes a period of accelerated expansion at very early times called inflation. It is a general feature amongst inflationary models that a background of gravitational waves (tensor perturbations) would have been generated during that era along with fluctuations in the density (scalar perturbations). Although difficult to be detected directly, primordial gravitational waves should have left an imprint on the temperature and polarisation of the cosmic microwave background (CMB) that may be detectable if inflation occurred at a sufficiently high energy scale [1][2][3].
Although scalar perturbations are now known to overwhelm any contribution from tensor perturbations to the temperature anisotropy (T ) or the gradient-like E-mode component of the polarisation, it is possible to form a curl-like component (B-mode) that at recombination is sourced only by tensor perturbations (in linear theory). A detection of a primordial B-mode is therefore widely regarded as a conclusive test of inflation. Furthermore, a measurement of the ratio of primordial tensor-to-scalar power, r, would reveal the energy scale of inflation. Current experimental bounds place r < 0.06 with 95 % confidence [4].
As we search for a signal below this level, it is not sufficient just to improve the sensitivity of experiments, for gravitational lensing of CMB photons by the large-scale distribution of matter in the Universe contaminates the signal by converting part of the primordial E-mode signal into B-mode [5]. The B-modes induced by gravitational lensing were first detected by [6]. These have a power spectrum resembling that of white noise with ∆ P = 5µK arcmin on the large angular scales where the primordial signal is expected to be strongest. Consequently, the uncertainty on measurements of the primordial tensor signal will henceforth be dominated by the large variance associated with the lensing component unless the latter can be removed.
Unfortunately, the removal of lensing B-modes cannot be achieved via multi-frequency cleaning, since lensing merely induces a remapping of photons by small deflection angles (of order a few arcminutes) while preserving the blackbody spectrum of the unlensed CMB. Progress can be made, however, by building an estimate of the lensing B-modes based on observations of E-mode polarisation and some estimate of the lensing deflections on the sky. This can then be subtracted from observed B-mode maps in a procedure known as delensing. Several groups have now successfully applied this technique to real data [7][8][9][10].
The extent to which the lensing signal can be removed depends on the fidelity with which the lensing potential can be estimated. For sufficiently low experimental noise levels, the lensing potential, φ, can be reconstructed internally from the CMB maps themselves by employing either the quadratic estimators of ref. [11] or the more powerful, albeit analytically and computationally complex, maximum-likelihood [12] or Bayesian methods [13,14]. For experimental noise levels available already in the next generation of CMB observatories, a quadratic combination involving Eand B-fields will provide a sizeable fraction of the signal-to-noise, dominating over other estimators in the regime where ∆ P < 5 µK arcmin.
However, ref. [15] showed that whenever there is an overlap in modes between the field we wish to delens and the fields from which a lensing reconstruction is derived, the delensed power is subject to a bias that leads to a suppression in power going beyond that which can be attributed to a mitigation of the lensing effects. More importantly, the variance associated with this delensed spectrum is also reduced, posing the question of whether the bias could help better constrain the tensor-to-scalar ratio, r.
In this paper, we set out to understand the biased delensed spectrum as a function of experimental characteristics and, importantly, primordial B-mode power, assuming no foreground or survey non-idealities. In section 2, we introduce the quadratic estimators for lensing reconstruction, focusing on the EB combination. In section 3, we show how an estimate of the lensing B-modes at the map level can be made by combining observations of E and estimates of φ. Then, in section 4, we introduce analytic models for the angular power spectrum of delensed B-modes in both the cases with and without overlapping modes, and compare them to simulations. (These models are derived carefully in appendix A, and the simulations are described in appendix B.) In section 5, we consider simple models for the covariance of the delensed B-mode power spectrum. We combine the power spectrum models and covariances in a maximum-likelihood framework in section 6, where we simulate inferences of r to compare quantitatively the different ways in which the bias can be dealt with. Finally, in section 7, we briefly study the impact of the bias on an alternative estimator which correlates observed and delensed B-modes.
The results in this paper will often refer to the experimental specifications of the upcoming Simons Observatory (SO) [16], which will feature a large-aperture telescope (LAT) responsible, among other scientific targets, for lensing reconstruction; and separate small-aperture telescopes (SATs) for observations of B-modes on large angular scales. However, the insights developed here apply more widely, and will likely be relevant to any application of internal delensing that uses information from the B-modes for the purpose of lensing reconstruction.

Quadratic estimators for lensing reconstruction
In order to undo the deflections induced by lensing, an estimate of the projected matter distribution on the sky -which determines the lensing potential -is required. For the sensitivities and resolution of current CMB experiments, the best possible estimate is obtained from tracers external to the CMB such as the cosmic infrared background (CIB) or very deep galaxy surveys [17][18][19][20]. The CIB retains a high degree of correlation with the smaller-scale lenses at high redshift that are important for converting E-modes into B-modes, something that is not yet possible with internal lensing reconstruction, which is very noisy on the relevant angular scales. However, this situation will change with planned CMB polarisation surveys, which can provide high signal-to-noise lensing reconstructions on nearly all scales relevant for B-mode delensing [21]. Consequently, reconstructions of φ derived from the CMB fields themselves will ultimately provide the best delensing performance.
Internal reconstruction techniques use the fact that, if we could average over realisations of the CMB while keeping the lensing potential fixed, lensing would break statistical isotropy by inducing correlations between fluctuations on different scales; consequently, the lensing potential can be reconstructed by combining many off-diagonal correlations. In general, the optimal internal reconstruction of φ will be obtained from studying the likelihood function of the lensed CMB temperature and polarisation anisotropies [12][13][14]22]. However, for the noise levels attainable in the near future, the optimal lensing reconstruction has been shown (see, e.g., ref. [23]) to be equivalent to that from the more tractable and computationally efficient "quadratic estimators" of refs. [11,24]. In fact, in the upcoming era of high-resolution CMB experiments such as the Simons Observatory [16] and SPT3G [25], which feature experimental noise levels 1 µK arcmin < ∆ P < 10 µK arcmin, the optimal reconstruction will arise from a combination of external tracers and internal reconstructions involving quadratic estimators [18].
A minimum-variance internal reconstruction of the lensing potential can be obtained from combining different quadratic estimators, as seen in figure 1. For low enough noise levels, the quadratic estimator involving E and B fields is expected to provide the highest signal-to-noise reconstruction -the reason being that, in this case, the dominant source of reconstruction noise involves the Gaussian contraction EE BB , and the lensing and primordial contributions to BB are small [23]. It is because of its relevance in upcoming efforts to delens the CMB that we focus on the EB estimator in this paper. It takes the form where E obs,LAT and B obs,LAT are beam-deconvolved observed fields 1 , A EB L is a normalisation factor, C XX, obs, fid, LAT l is the fiducial, total power observed by the LAT which is used to inversevariance filter 2 the field X, and the geometric coupling Here, ψ l is the angle between l and the x-axis used to define positive Stokes parameter Q, and similarly for ψ l . Notice the use of the fiducial lensedC EE,fid l instead of its unlensed counterpart, as this choice has been shown to mitigate higher-order error terms in estimates of the lensing potential power spectrum [26,27]. If eq. (2.1) were to minimize exactly the variance ofφ EB , it would include an additional term proportional to the primordial B-mode spectrum, C BB,t l , like the one shown but replacing W (L − l, −l)C EE,fid l → −W (−l, L − l)C BB,t,fid |L−l| . We ignore it here, however, as its effect is negligible for values of r compatible with experimental bounds.

Template delensing of the B-mode
Once we have an estimate for the lensing potential, delensing can be performed by using it to build a template approximating the lensing B-modes, B temp , which can then be subtracted from  Figure 1: Lensing reconstruction noise levels associated with the five different quadratic estimators, and their minimum-variance combination, for an experiment with resolution θ FWHM = 1.5 arcmin (full width at half maximum) and isotropic white noise with ∆ T = 6 µK arcmin for temperature and ∆ P = √ 2∆ T for polarisation. The maximum multipole used in the reconstruction is l max = 3000. These specifications are similar to the goals for the SO LAT at 145 GHz. The dashed blue line is the lensing potential power spectrum.
the observed B-modes, 3 B obs, SAT , to yield a delensed field To leading order in ∇φ (what is usually referred to as the "gradient approximation", and is known to be an excellent approximation to the non-perturbative calculation of the lensed Bmode spectrum on large angular scales [28]), the template can be constructed as with W (l, l ) as defined in eq. (2.2). The function f (l, l ) can be chosen to minimise the power spectrum of the delensed field (3.1), which ref. [17] shows to be the case when the fields are filtered with where are Wiener filters for the E-modes and the estimate of the lensing potential. Here, N φφ L is the reconstruction noise onφ. In practice, we use N φφ L ≈ N (0),EB L , the disconnected (Gaussian) contribution to the power spectrum ofφ, and ignore higher-order terms such as N (1),EB L , as the former dominates over the latter by several orders of magnitude on all relevant scales [29]. Notice also that eq. (3.2) uses the observed (lensed and noisy) E obs,LAT , rather than unlensed (or delensed) E-modes. This is to be preferred for template delensing since an approximate cancellation suppresses higher-order contributions in C φφ L to the delensed power spectrum, which would otherwise be as big as 30 % on large scales [30].

Angular power spectrum of delensed B-modes
Most inflationary models predict a spectrum of primordial perturbations that are very nearly Gaussian distributed. This character extends approximately to the fluctuations of the primary CMB, as linear-theory transfer functions are an excellent approximation in the early universe. For such Gaussian fields, the data can be losslessly compressed into a power spectrum, making the latter one of the most powerful tools for studying the CMB. Although the lensed CMB is actually non-Gaussian (as a result of the Gaussian primary CMB fields being displaced by the nearly-Gaussian lensing potential), delensing has been shown to mitigate its non-Gaussian character to the point where the delensed B-modes -if delensed in such a way that the EB delensing bias is avoided -deviate only slightly from Gaussianity [31] and the information lost from working with the power spectrum is relatively small.
More importantly, in practical applications involving masked observations and anisotropic noise, any computation of the exact likelihood of the data becomes intractable and one must work instead with approximate forms, most typically involving estimators of the power spectra and their covariance [32].
For these reasons, most efforts to detect primordial B-modes work with the power spectrum of the delensed field of eq. (3.1): where in the last line we have used the statistical isotropy of the CMB to define C BB,del l , the angular power spectrum of the delensed B-modes. In the remainder of this section we examine this expression in detail. In section 4.1, we evaluate it in the case where the statistical errors in the estimated lensing potential are independent of the lensed CMB fields, exploring the reduction in power associated with removal of part of the lensing signal. Then, in section 4.2, we study eq. (4.1) in the case whereφ is derived from an EB quadratic estimator; new couplings then appear that further suppress the delensed power spectrum beyond a simple removal of lensing power. We conclude the section by proposing an analytic model to capture the behaviour of such "bias" terms.

The unbiased case: reconstruction errors statistically independent of the lensed CMB
If the noise onφ were independent of the B-mode we would like to delens, as would be the case if an external tracer were used forφ, the delensed power spectrum would take the form , there is a contribution from residual deflections -imperfect delensing -given to leading order in lensing by [17] C BB,res where we have defined and assumed that our fiducial model for the lensing power spectrum, C φφ,fid l , is correct. We note that C W l is simply the fiducial power spectrum of the B-mode template. In appendix A.1, we explain how eq. (4.3) can be recovered in the analytic framework developed in section 4.2 to characterise the delensing bias. Notice that, in the limit of no observational noise and a perfect φ, we have W E l → 1 and W φ l → 1 and all of the (leading-order) lensing signal is removed. Higher-order terms mean that template delensing can, in fact, be used to reduce the lensing contribution to the power spectrum to approximately 1 % of its original level. This is in the case where lensed E-modes are employed in the construction of the B-mode template -the alternative of using unlensed/delensed E-modes performs worse with a lensing residual of order 30 % as noted above [30]. Were estimates of the lensing potential accurate enough to reduce the delensed power to the 1 % level, template delensing should be replaced with non-perturbative methods whereby the Wiener-filteredφ is used to remap the observed CMB directly (see, e.g., [7,9,10,19] for demonstrations of this method).
Henceforth, we will refer to eq. (4.2) as the unbiased delensed power spectrum in order to differentiate it from the case where the errors inφ are statistically dependent on the lensed CMB, as is the case when internal delensing with overlapping modes. We now consider the new terms that arise in the delensed power for this latter case.

The biased case:φ obtained from an EB quadratic estimator
We discussed in section 2 that, in the upcoming era of low-noise CMB experiments, the highest signal-to-noise reconstruction using quadratic estimators will arise from combining E and B fields. Crucially, ref. [15] first noticed that, if such a reconstruction is used to delens observations of B-modes, the delensed power spectrum will be biased: while eq. (4.3) continues to quantify the power spectrum due to residual deflections, the total delensed spectrum will show an additional suppression of power beyond that which can be attributed to delensing. In this section, we add to the work of [15,33,34] by providing a detailed calculation of the biased delensed B-mode power spectrum, including terms neglected by eq. (4.3).
In order to isolate the effects of the bias, we assume henceforth that the lensing reconstruction is obtained exclusively from an EB quadratic estimator. Although this is set to be the dominant source of information on lensing for upcoming experiments, in real applications the optimal reconstruction will, in fact, arise as a minimum-variance (MV) combination of several quadratic estimators. To recover the amplitude of the bias on both the delensed B-mode spectrum and the tensor-to-scalar ratio appropriate to those cases, the results quoted in this paper would need to be scaled appropriately by the MV weight pertaining to the EB estimator.
In this biased case, the first term in the delensed B-mode spectrum of eq. (4.1) is still   The last term, which correlates two templates, takes the form The evaluation of the four-and six-point functions that appear in eqs (4.5) and (4.7) is discussed in detail in appendix A. These are expanded in terms of connected n-point functions with n = 2 and 4 (the connected six-point function is higher order in C φφ l ). A subset of these terms combine to give the standard unbiased result (4.3), as shown in appendix A.1. The remaining terms introduce corrections, with the most important of these identified and evaluated in appendix A.2. Combining these results, we show that the biased delensed B-mode power spectrum can be approximated as where the correlation of experimental noise in the SAT and the LAT is denoted as We have also defined whose origin and properties are discussed further below. Equation (4.8) captures the general case where a large-aperture telescope (LAT) focuses on lensing reconstruction while a separate small-aperture telescope (SAT) makes observations of the B-modes on large angular scales. In such a configuration, the experimental noise is uncorrelated between instruments and N X l = 0. On the other hand, the case where a single telescope is used for both purposes can easily be recovered by letting N X l = N BB,LAT l = N BB,SAT l , in which case the biased delensed spectrum reduces to All of the correction terms in eq. (4.8) are proportional to one or two powers of D l ; for D l = 0 it reduces to the unbiased result (4.3). The function D l arises when one contracts the B-mode template over the pair of observed E-modes that enter explicitly. In detail, (4.12) All of the correction terms retained in eq. (4.8) arise from such contractions. For example, since Eand B-modes are uncorrelated, the disconnected contribution to B temp (l 1 )B obs, SAT (l 2 ) gives a term (4.13) Note how this involves the cross-correlation of the experimental noise between the SAT and the LAT, N X l , and the tensor B-mode power. It is the main correction term and, given that it enters the delensed power through −2 B temp (l 1 )B obs, SAT (l 2 ) , suppresses the power further. (The other correction terms we retain increase the delensed power.) We can gain further insight into D l by noting that in the limit of noiseless E-mode observations, and setting the normalisation of the quadratic estimator to A EB L = N (0),EB L (which is correct in the usual limit that the fiducial spectra used to inverse-variance filter the CMB fields in the lensing reconstruction are close to the true total power), D l → C BB,res l /C BB,obs,fid,LAT l [cf. eq. (4.3)]. Consequently, in this limit D l can be interpreted as the ratio of residual lensing power to fiducial power used to inverse-variance filter the B-modes in the lensing reconstruction. More generally, the Wiener filter W E l will reduce D l below this value. In all cases 0 < D l < 1. In figure 2, we see that D l depends strongly on large scales on the fiducial primordial B-mode power contained in the inverse-variance filters and parametrized by r fid . (The figure includes large values of r fid inconsistent with observations but included for illustration.) To understand this behaviour recall that, for r ∼ 0.01, C BB,t l becomes comparable to the lensing power on large angular scales between 10 < l < 100. By raising r fid above this value, we are effectively down-weighting the large-scale lensing B-modes in the fields from which the reconstruction is derived. Although this will not significantly affect the fidelity of the lensing reconstructionwhich is derived from information coming chiefly from smaller angular scales -it will give a sizeable reduction in D l and the bias in the total delensed power on large scales.
The result in eq. (4.8) is a generalisation of that in ref. [33] (eq. A16 there), which assumes r = 0 and separate LAT and SAT observations (N X l = 0). An important insight from our result is that the bias term responsible for the additional suppression of power, eq. (4.13), is proportional to the total observed B-mode power, and not simply its lensing component 4 . Crucially, our result thus predicts a suppression of the primordial signal whenever it is present, as first noted in ref. [15]. This effect is also hinted at in ref. [34], though the fact that their simulations have r = 0 prevents them from making a quantitative statement.
In order to assess the validity of the model in eq. (4.8), we compare it to the output of applying the reconstruction and delensing pipelines to simulations (detailed in appendix B) of  Figure 3: Upper panel : noise-subtracted simulated power spectra of the observed B-mode (green), biased residual after delensing withφ EB (yellow) and residual after delensing with a reconstruction cutoff at l cut = 300 (cyan). The shaded regions represent the variance of each of the delensed spectra, obtained from simulations. Notice that the biased spectrum (yellow) has less variance than its unbiased counterpart (cyan). Also plotted are noise-subtracted theoretical curves modeling delensed spectra with (red, dot-dashed) and without (black, dashed) an l cut .
Middle Panel : ratio of model-to-simulated delensed power spectrum amplitude for the biased (black, dashed) and unbiased (red, dot-dashed) cases. Lower panel : ratio of delensed power spectrum variance in the case where no l cut is imposed to that where l cut = 300, obtained from simulations. All curves are for an experiment with characteristics similar to the Simons Observatory, as described in the text.
an experiment with characteristics similar to those of SO [16]. Specifically, the LAT survey has angular resolution θ FWHM = 1.5 armin and temperature noise level ∆ T = 6 µK arcmin (with the polarisation level ∆ P = √ 2∆ T ), as in figure 1. It is the simulated data from this telescope which we use for the purpose of lensing reconstruction, with l max = 3000. The SAT survey has θ FWMH = 17 arcmin, ∆ T = 2 µK arcmin and ∆ P = √ 2∆ T . Where needed for calculations of the power spectrum variance, we assume the SAT survey covers approximately 2.8 % of the sky (and is fully contained within the wider LAT survey). The simulations have tensor-to-scalar ratio r = 0.01. The results are shown in figure 3. The biased delensed power spectrum (yellow curve in the top panel) is seen to be in excellent agreement with eq. (4.8), shown as the dashed black curve. The differences are smaller than 0.5 % on all scales. For lower noise levels, the level of agreement is expected to worsen as bias terms neglected by our model play an increasingly significant role. Since the most relevant of those terms (the contribution to the six-point function from the first term on the right of eq. A.10) makes a positive contribution to the power, eq. (4.8) is expected to underestimate the amplitude of the biased delensed power spectrum. Indeed, we simulate a single-telescope experiment with θ FWHM = 6 arcmin and ∆ T = 3 µK arcmin and find that the modelled biased spectrum is around 2-3 % lower than its simulated counterpart.
The bias in the delensed spectrum is around 80 % of the signal power, C BB,t l + C BB,res l , in the unbiased delensed spectrum for the configuration in figure 3 (see discussion below). This bias would need to be modelled or otherwise mitigated in a likelihood analysis if inferences on the tensor-to-scalar ratio are not to be biased. One simple way to remove the bias in the delensed spectrum on the scales relevant for searches for primordial B-modes was noted by ref. [15]: exclude from the reconstruction any B-modes that overlap in scales with the B-modes we wish to delens. This follows from the local character of the bias, whose origin is given in eq. (4.12), meaning that it will be avoided at a given multipole as long as that multipole is removed from the B-field used for the lensing reconstruction. For the purpose of primordial B-mode searches, we shall remove reconstruction B-modes on the largest angular scales l ≤ l cut . In this case, eq. (4.12) becomes where Θ(l) is the Heaviside function, and so the bias terms in the delensed power spectrum vanish for l ≤ l cut . Note that the normalisation and noise of the EB quadratic estimator are changed, albeit by a small amount, on all scales by excising the input B-modes at certain scales, and D l (for l > l cut ) needs to be calculated with the modified normalisation and Wiener filter; see eq. (4.10). Equivalently, we can think of effectively setting to infinity the noise on the B-modes used in the reconstruction on the scales that we wish to exclude. In this case, D l goes to zero for l ≤ l cut , and both eq. (4.8) and (4.11) reduce to the unbiased spectrum of eq. (4.2). Averting the bias in this way comes at the cost of a lower signal-to-noise lensing reconstruction as information is being discarded, although the degradation is expected to be minor since most of the lensing information is obtained from smaller scales of the polarisation fields. In order to take this effect into account when computing the delensed power spectrum below the cutoff using eq. (4.3), we must make sure to calculate C W l using a Wiener filter that includes the higher reconstruction noise. When this modification is included, the simple unbiased model agrees very well with the simulated spectrum, as shown in figure 3 for the case of l cut = 300, with differences between the two less than 0.2 % for an SO-like experiment. Imposing a cutoff at l cut = 300 leads to a lower signal-to-noise lensing reconstruction and a small decrease in C W l which, on large angular scales, results in an increase of the residual lensing power spectrum and a degradation of delensing efficiency by approximately 2 %. Several techniques have been suggested to reduce this small degradation, including realisation-dependent methods [34] and an upgrade of the cutoff approach which involves splitting the whole multipole space into "notches" within which the bias can be avoided by excluding modes common to the notch and the reconstruction [15,35].
Although the bias from internal delensing can be easily removed on large scales with the techniques just discussed, an alternative is to model explicitly the biased delensed spectrum as part of any subsequent likelihood analysis. Figure 3 suggests it is worth exploring this approach further since the bias not only suppresses the delensed B-mode spectrum beyond what can be attributed to the (partial) removal of lensing effects, but also reduces its variance. Indeed, in figure 3 the variance of the biased delensed spectrum is roughly 65 % of that associated with the unbiased case (the latter with l cut = 300). If the power spectrum bias were independent of the primordial B-mode power, this reduction in variance would translate directly into improved constraints on the tensor-to-scalar ratio (assuming there is no significant difference in the non-Gaussian correlations between the power spectra at different multipoles; see section 5). However, part of the bias does depend on the primordial spectrum, suppressing the contribution of C BB,t l in the biased delensed spectrum by a factor (D l − 1) 2 (see eq. 4.8).
The suppression of the primordial signal in the biased delensed spectrum acts as a multiplicative bias. We can remove this by renormalising by (D l −1) −2 , in which case the contribution from C BB,t l is the same as for unbiased delensing: (4.15) with C BB,del,unbiased l as defined in equation (4.2). The last three noise terms on the right-hand side vanish for the case of B obs,LAT = B obs,SAT , i.e., when the same maps are used in all parts of the analysis, and are positive when B obs,LAT and B obs,SAT have independent instrument noise (N X l = 0) since 0 < D l < 1. On the other hand, the term proportional to C W l is always positive. Crucially, this means that, in presence of the bias, the renormalised delensed power is generally larger than from the unbiased approach. 5 In other words, the primordial power is always suppressed by a larger fraction than the sources of (lensing and experimental) noise. This is illustrated in figure 5 where we see that, for the experimental configuration adopted in figure 3, the signal-to-noise noise on a primordial component is lower when the bias is allowed to play a part and modelled than in the alternative approaches where the bias is avoided -worse, even, than in the case of no delensing. In section 6, we shall see that this indeed translates to errors on estimates of r from a maximum-likelihood inference pipeline that inevitably increases (in a statistical sense) whenever the bias is modelled instead of avoided. For noise levels lower than those considered here, we expect the approximate expression for the biased delensed spectrum, equation (4.8), to underestimate the true power, as noted above, with the dominant correction coming from the first term on the right of eq. (A.10). In this term, B-modes only enter through a connected four-point function of the form BEEB c , which receives no contribution from tensor modes (ignoring lensing for these). It follows that the additional term adds power but does not couple to C BB,t l , and so further reduces the signal-to-noise on the primordial signal. The renormalised spectrum of eq. (4.15) is lower, all other things being equal, if the same maps are used for all parts of the analysis. The reason is that in this case the experimental  Observed (no delensing) Unbiased delensed Unbiased delensed with lcut = 300 Normalized biased delensed Figure 5: Effective noise power (i.e., the difference between the total delensed power and the primordial contribution) for different estimators from which a primordial B-mode signal (black dashed) is to be extracted. In the biased case, the spectrum is first renormalised so that the amplitude of the primordial power is the same as in the unbiased cases. The experimental set-up is the same as in figure 3.
noise in the reconstruction and in the B-modes to be delensed are correlated, allowing the disconnected couplings (A.15) and (A.16) to suppress the experimental noise contribution to the delensed power spectrum by the same fraction as for the primordial signal. If, on the contrary, the noise is uncorrelated between the reconstruction and the B-modes to delens, it will not contribute to the disconnected term (A.16) and will consequently be suppressed by a smaller amount than the signal -in fact, it will be amplified relative to the unbiased case by coupling (A.15). The degradation in signal-to-noise on the primordial B-mode power in the biased case can be understood by noticing that, in the case N BB,LAT l = N BB,SAT l = N X l , the disconnected couplings (A.15) and (A.16) combine with the terms retained in the unbiased calculation to yield a suppression of the primordial power and experimental noise by an equal fraction, while for the lensing contribution there is an additional (positive) term, eq. (A.11), which reinstates some lensing power while leaving the primordial signal unaffected. This last contribution arises from the connected coupling of the B-mode in one of the templates and the EEB fields in the other template. This is not the only term that restores lensing power, since the partially-connected six-point function (A.6) coupling the lensing part ofφ across templates has a similar effect. However, the latter is already included in the unbiased calculation.

Covariances of delensed power spectra
The precision with which we can isolate a primordial component from the observed (or delensed) B-mode power spectrum is ultimately determined by the covariance of the latter. The non-Gaussian covariance of lensed CMB spectra is well understood by now [36][37][38][39] and models exist that allow for its numerical evaluation. In this work, we use the publicly-available code LensCov [39] to compute the theoretical bandpower covariances for observations corresponding to our choice of cosmology and experimental parameters. We find good agreement between this theoretical computation and results from simulations after applying the appropriate binning and sky fraction corrections, as illustrated by figures 6 and 7. For visualisation purposes, we plot the cross-correlation coefficient of bandpowers, defined as where Cov(C BB l 1 ,C BB l 2 ) is the covariance of the binned angular power spectrum and the angle brackets denote averaging over realisations of either the lensed or delensed CMB and noise. The delensed bandpower covariance is more complicated, but it can begin to be understood by studying the case where the lensing reconstruction is independent of the CMB. Reference [31] presents the following extension of the lensed B-mode power covariance to this delensed case under the assumptions that the E-modes used in the template are cosmic-variance limited on all relevant scales (i.e., for multipoles l 2000) and that φ is Gaussian: The assumption that E-modes are limited by cosmic variance on the relevant scales is not strictly true for the specifications adopted here. Although it is not difficult to generalise eq. (5.3) to account for noise in the E-modes, ref. [31] finds that the covariance without E-mode noise still matches well that obtained from simulations of experiments with the noise levels of current and upcoming experiments. We evaluate equation (5.3) by modifying the LensCov code for the calculation of lensed B-mode covariances. Along with a modified Gaussian variance (the first term on the right of eq. 5.3), we must also replace the derivative terms with those appearing above. We calculate these analytically from the leading-order expression of equation (4.3). The results are shown in figure 6. From this figure, and more acutely from figure 7, it is clear that delensing reduces the non-Gaussian correlation between scales introduced by lensing [36], effectively increasing the number of independent pieces of information available. Figures 6 and 7 also suggest that the simple model of ref. [31] is consistent -up to Monte-Carlo errors from the finite number (25 000) of simulations 6 -with the simulated covariance of delensed B-mode bandpowers below the cutoff, whereφ is effectively independent of B.
Finally, we study the bandpower covariance in the case where no l cut is imposed and thus the delensed B-mode spectrum is biased. An analytic exploration of this case is beyond the scope of this paper, since the framework of appendix A.2 suggests that the covariance receives contributions from 12-point correlators of lensed polarisation fields although some simplifications 6 It can be shown (assuming that the errors on estimates of individual bandpowers are distributed as Gaussian random variables with appropriate correlations) that the fractional error on an element of the correlation matrix, estimated from N simulations, is σ(Rij) where Rij are the true correlation coefficients between bandpowers and the mean bandpowers are known in advance. In the more usual case where the variances are also estimated from the data, the relation is replaced with  [39], as described in the text. Right column: observed (lensed and noisy) bandpowers. Central column: delensed bandpowers in the case where the bias is avoided by either introducing a cutoff at l cut = 300 (top) -thus avoiding the bias for l < l cut -or by working with aφ that is independent of B yet just as correlated with the actual φ as the reconstructedφ EB employed in the simulated analysis (bottom). Left column: biased delensed bandpowers from simulations.
are likely possible. It can, however, be obtained from simulations, as shown in figure 6. Interestingly, at the noise levels considered here, the lower variance induced by the bias does not appear to be accompanied by the larger off-diagonal correlations between bandpowers seen by ref. [40]. The most likely explanation is that the Gaussian experimental noise dominates the variance on the small angular scales where we would expect the bias to induce strong covariance between delensed bandpowers (cf. figure 11 of [40]). This suppresses the correlation and increases the fractional Monte-Carlo error to the point where the signal gets buried in the noise.
In summary, we have seen that, for an experiment with characteristics similar to SO, the EB lensing reconstruction bias brings about a reduction in the variance of the power spectrum of delensed B-modes while the extent of the correlations between bandpowers remains comparable to the unbiased case. If the primordial signal were not suppressed by the bias, analysis of the biased delensed bandpowers (with appropriate modeling of the bias) would lead to improved constraints on the tensor-to-scalar ratio. However, we have already seen in section 4.2, that the bias acts multiplicatively on the primordial power, suppressing the signal more than the Gaussian variance. We thus expect constraints on r to worsen in this case, which we now demonstrate through a simulated maximum-likelihood analysis.

Maximum-likelihood inference of r
Although the likelihood of CMB temperature and polarisation data can in principle be written down exactly, even accounting for the effects of lensing [13,14,22], in practical applications with realistic survey complexities, easily-computable summary statistics (e.g., angular power spectra) and their approximate sampling distributions are typically preferred (see [41] for a review). For power spectra of approximately Gaussian fields, several such likelihood approximations involving the estimated power spectra and their covariance appear in the literature. We choose to work with the approximation developed in ref. [32], and which was used by the BICEP/Keck collaboration to analyse their most recent data [4]. For our application here, it takes the form ) denotes the ith empirical (model) bandpower of the delensed B-mode spectrum. Only multipoles below l max = 280 are used in the inferences that follow. An advantage of using this approximate likelihood is that the non-Gaussian character introduced by lensing can be accounted for through the fiducial delensed bandpower covariance, Cov(C BB,fid i , C BB,fid j ), which only needs to be computed (either analytically or via simulations) and inverted once. As long as it matches the fiducial C BB,fid l , the exact form of this fiducial covariance has a negligible impact on the resulting inferences. Indeed, we verify that there is no appreciable change in the inferences when we vary the fiducial level of primordial B-mode signal. Furthermore, ref. [32] shows that, even if the fiducial model deviates from the truth, the likelihood is still exact in the full-sky, isotropic limit with Gaussian fields.  Hence, all analyses presented henceforth use fiducial bandpower amplitudes and covariances with the same r input . The question remains whether the bandpower covariances ought to be modelled or whether they could simply be obtained from simulations. In section 5 we saw that, while it is possible to write a simple analytic model for the covariance of unbiased delensed bandpowers, doing the same for the biased case is more difficult. To address this question partially, we compare the distribution across simulations of the maximum-likelihood estimates for the tensor-to-scalar ratior ML and their associated uncertainties σ(r ML ) -derived from the second derivative of the log-likelihood function at the maximum-likelihood point -in the cases where the bandpower covariance is obtained analytically (using LensCov) or from simulations. We do this for the cases of no delensing and unbiased delensing, the latter excluding B-modes below l cut = 300 from the simulated fields used for reconstructingφ EB . The simulations used in these comparisons, and for all results presented in this section, assume the experimental set-up described in section 4.2. The likelihood curves for these two cases are shown in figure 8 for a typical realisation, and the distribution ofr ML and σ(r ML ) across 25 000 simulations are given in figure 9. We find only a slight degradation in the errors on r when simulated bandpower covariances are employed, justifying our ensuing use of simulated matrices.
We illustrate the need to account for (or mitigate) the bias in the delensed B-mode power spectrum when delensing withφ EB with the following naive analysis. For the empirical spectrum C BB l used in the likelihood of eq. (6.1), we take the simulated biased delensed spectrum, but all other fiducial and model spectra and covariances in the likelihood are unbiased -that is, calculated assuming thatφ is independent of the CMB but with the same correlation to φ aŝ φ EB . In figure 10, we verify that ignoring the suppression in the biased empirical spectrum leads to inferences on r that are biased low, with a significant number of likelihoods peaking (in the range r ≥ 0) at r = 0, and artificially small typical errors on r.
We have already discussed, in section 4.2, several ways in which the bias can be taken into account: modelling the biased spectrum; renormalizing the biased spectrum to restore unit response to primordial power; or imposing a low-l cutoff on the B fields used for estimatingφ. The first two approaches should yield equivalent results on r, but are expected to be very nonoptimal since the bias reduces the primordial power by a larger fraction than the rest of the power (see figure 5). Mitigating the bias with a cut-off is preferred, although there is a small penalty due to the lower signal-to-noise lensing reconstruction in this case. In figure 10, we quantitatively compare these methods. For our experimental set-up and input value r input = 0.01, we see that modelling the bias inflates the errors on r by typically 15 % 7 compared to when a cutoff at l cut = 300 is imposed, translating to a slightly wider distribution ofr ML and a lower number of detections of non-zero r. The figure also shows that, for our set-up, modelling the bias degrades the sensitivity on r to a level comparable to (or even slightly worse than) the case of no delensing, as suggested by our earlier results for the normalised power in this case (figure 5).

The B obs (B obs − B temp ) estimator
In section 4.2, we saw that the bias from internal delensing with an EB estimator reduces the signal-to-noise on a primordial contribution to the power spectrum of delensed B-modes. We attributed this to the coupling in eq. (A.11), which arises in the cross-correlation of the Bmode lensing template with itself, and that restores some of the lensing power while leaving the primordial signal unaffected. This motivates the idea of considering an alternative power spectrum estimator of the form B obs (B obs − B temp ), i.e., the cross-correlation of the observed B-modes with the delensed B-modes, B del = B obs − B temp . This estimator does not involve the cross-correlation of the template with itself, and so avoids couplings such as that in eq. (A.11).
In the case of no l cut in the lensing reconstruction (such that the B-modes used for reconstruction Orange: biased delensing, based on a reconstructedφ EB that does not impose a low-l B-mode cutoff, with no attempt to model the resulting power spectrum bias in the likelihood. This results inr ML being biased low (and a significant number of simulated likelihoods that peak at r = 0) and with artificially small errors. Green: biased delensing, but with the power spectrum bias modelled as described in section 4.2. This mitigates the bias in r but, for our experimental set-up, has comparable performance to no delensing. Red : a cutoff at l cut = 300 is imposed on the B-modes used in the reconstruction, avoiding the bias to the delensed B-mode power spectrum below l cut . There is no bias in r in this case, and delensing reduces σ(r ML ) as intended. All analyses employ simulated bandpower covariances.
overlap with those to be delensed on all scales) and two separate surveys, LAT and SAT, the expected value of the cross-correlation is In the case of a single survey, this reduces to The last term on the right-hand side is always negative. Renormalising the cross-correlation by (1 − D l ) to have unit response to the primordial power, we see that its renormalised value is reduced below the power of the unbiased, delensed spectrum (equation 4.2) by this negative term. While the non-primordial power in both of these spectra decreases monotonically with the experimental noise, the difference between the two is greatest when C W l D l /(1 − D l ) peaks. This is the case when the experimental noise level in polarisation is approximately 3 √ 2 µK arcmin, corresponding to a delensing efficiency of approximately 40 %. For such an experimental configuration, the term −D l C W l reduces the cross-correlation in eq. (7.2) by an amount equivalent on large scales to white noise with ∆ P ≈ 2 µK arcmin after normalising to have unit response to a primordial signal. It is instructive to compare the B obs (B obs − B temp ) estimator above to another defined along similar lines, but which features no overlap between reconstruction B-modes and Bmodes to be delensed (for example, by introducing an l cut in the reconstruction, such that B obs (B obs − B temp ) is unbiased on scales l < l cut ). This alternative estimator has expectation value which is the same as the unbiased, delensed power spectrum of eq. (4.2). Comparing this to the biased estimator of eq. (7.2), we see that, after appropriate renormalisation, the expectation value of the latter is lower than the estimator in eq. (7.3), or even the standard power spectrum approach of eq. (4.2). Our analytic models for all these different estimators are found to be in excellent agreement with their values in simulations.
In assessing the performance of the cross-correlation estimators for constraining the primordial signal, we remind the reader that, for a for a cross-correlation, the variance is not determined by the expected value (even in the limit of Gaussian fields). For a single survey, treating both B obs and B del as Gaussian fields, the variances of the estimators at hand are . Clearly, the Gaussian variance is smallest for the usual B del B del unbiased estimator. Moreover, we will now see that the (co)-variance of both the cross-estimators is dominated by non-Gaussian effects, which further degrade their constraining power.
In figure 11, we use simulations of an experiment with ∆ P = 3 √ 2 µK arcmin and θ FWHM = 6 arcmin to quantify the variance associated with the different estimators. Although the Gaussian variance of eq. (7.6) appears to capture correctly the simulated behaviour of the unbiased B del B del estimator, the unbiased and renormalised B obs B del estimators show variances well in excess of the Gaussian result. This can be attributed to the predominance of non-Gaussian contributions, which also manifest themselves in the strong correlations between bandpowers shown in figure 12. Intuitively, the reason why the non-Gaussian character is less acute in B del B del than in B obs B del is that the former involves two delensed fields and delensing is known to mitigate the correlations between scales introduced by lensing.
Given these considerations, we believe that the estimator of eq. (7.2) should lead to inferior constraints on r compared to an analysis involving the usual power spectrum estimator of eq. (4.2). A more careful demonstration of this would require working within a maximumlikelihood inference framework, as in section 6. Unfortunately, this is not straightforward since standard likelihood approximations are not applicable for a single cross-spectrum. These likelihood approximations are very convenient in the presence of real-world effects such as masking, and for including non-Gaussian power spectrum covariances. Extending these to a single crossspectrum would require further development.

Discussion
Searches for primordial B-modes in CMB data already necessarily require the removal of the contamination from gravitational lensing. In the near future, experimental characteristics will be such that the EB quadratic estimator is expected to provide a large fraction of the signalto-noise on reconstructions of the lensing potential, which is used to perform the delensing procedure. It is well known that any overlap in modes between the B-field to be delensed and the B-field from which the reconstruction is derived biases the amplitude and variance of the delensed power spectrum [15,34].
In this paper, we have refined the modelling of this bias, based on a cumulant expansion and identification of the dominant terms, and used simulations to verify its efficacy. Our model is consistent with others in the literature [33,34] that are valid in the absence of a primordial (gravitational wave) signal, but extends our understanding to include the primordial component. The power spectrum bias acts multiplicatively on the primordial power, suppressing its amplitude. The bias also suppresses lensing and instrument noise power but, crucially, generally less so than for the primordial signal. This leads to a loss of constraining power for primordial B-modes, and this becomes more acute whenever separate surveys (i.e., with independent instrument noise) are used for lensing reconstruction and the study of B-modes on large angular scales.
Fortunately, the bias is completely local in multipole space in the sense that bias to the delensed B-mode power at multipole l originates from B-modes at that same multipole being used in the lensing reconstruction. The bias can therefore be avoided by eliminating overlapping modes in the B-field employed in the reconstruction, and within the context of searches for primordial gravitational waves this can be done with little loss of performance. We have shown that, for any reasonable l cut , such elimination of overlapping modes is actually preferable over an approach where the bias is modelled, as it yields improved constraining power. These findings are verified by simulating the reconstruction and delensing procedures and performing maximumlikelihood inferences of the primordial power on the resulting delensed spectra.
Our model for the biased, delensed B-mode power spectrum is in excellent agreement with simulations for experimental noise levels comparable or inferior to those of the Simons Observatory, but we note that the agreement worsens for lower noise levels as new terms that we ignore become relevant. In this limit of low noise, we expect our model to over-estimate the signal-to-noise so that our claim that retaining, but modelling, the bias hinders efforts to constrain the primordial signal continues to hold. We therefore recommend that any attempt to delens large-scale B-modes internally building on lensing information gleaned from the B-modes themselves -be it via quadratic estimators or likelihood-based techniques -feature the excision of large-scale B-modes from the lensing inference.
Finally, we have considered alternative cross-correlation estimators for the primordial power involving the observed and delensed B-modes, with a view to removing some couplings in the biased spectrum that degrade the signal-to-noise on r. We provided analytic models for the expectation values of these cross-correlations and verified these against simulations. However, we showed that ultimately these estimators have lower signal-to-noise to primordial power than the auto-power spectrum of the (unbiased) delensed fields, due to their cross-correlation character and their increased non-Gaussian covariance between multipoles. mode template with the observed B-modes, eq. (4.6), and the power spectrum of the template, eq. (4.7), respectively.
It is convenient to expand these four-and six-point functions in terms of connected n-point functions. In the absence of lensing, the CMB would be Gaussian and all connected n-point functions with n > 2 would vanish. As the CMB fields are zero mean, the evaluation of n-point functions would reduce to products of two-point functions. However, in section 5 we saw that lensing distorts the Gaussian primordial statistics, introducing significant non-Gaussianities in the form of non-vanishing, higher-order connected n-point functions, which we shall also refer to as "connected correlators". In particular, the connected four-point function, or trispectrum, induced by lensing lies at the heart of any effort to infer C φφ L from the lensed CMB [26,29,42,43]. Note that Gaussian instrumental noise does not contribute to the connected correlators of the observed CMB fields. Moreover, if we ignore the impact of lensing on the (otherwise Gaussian) gravitational wave contributions to the CMB (which should be a good approximation given their power falls rapidly on intermediate and small scales), these contributions are independent of the lensed, scalar contribution and so do not contribute to the connected correlators either.
We organise the expansion of the n-point functions with the aid of a graphical representation where fields are represented by nodes, drawing those that are observed by the LAT in red, and those observed by the SAT in green. Lines connecting n nodes denote the connected n-point function of the associated fields, which can be evaluated perturbatively to the desired order in lensing. In order to preserve generality, we ensure that the case where a single telescope is used for both construction of the delensing template and the observation of the large-scale B-modes can be recovered by letting N X l = N BB,LAT l = N BB,SAT l , with N X l = 0 otherwise. In this way, the four-point function appearing in eq. (4.6) can be represented as where terms of the form EB have vanished due to parity invariance. The Band E-fields involved in the lens reconstruction are at the top and bottom of the diagrams, respectively, and the E-modes used further in the template are at the midpoint on the left. The two diagrams on the right of eq. (A.1) correspond to a trispectrum E obs E obs B obs B obs c and a product of twopoint functions E obs E obs B obs B obs . Similarly, the six-point function appearing in eq. (4.7) can be decomposed into

A.1 Contributions that appear in the unbiased delensed power spectrum
The standard calculation of eq. (4.3) assumes that the statistical noise in the lensing reconstruction,φ, that appears in B temp is independent of the CMB fields. For this reason, it includes only some of the terms making up eqs. (4.6) and (4.7) or, equivalently, some of the couplings represented by the diagrams of eqs. (A.1) and (A.2). We calculate these here, before assessing the remaining terms in section A.2.
We start with the cross-correlation of the template with the observed B-modes, i.e., −2 B temp (l 1 )B obs (l 2 ) . Only the part involving the connected four-point function is included in the standard calculation, and then only the "primary coupling" in which there is a contraction over the unlensed E-modes across the two legs of the quadratic estimatorφ. This yields where C W l is defined in eq. (4.4). Here, we have introduced the notation, for example,Ẽ[E, φ] to show the functional dependence of the fieldẼ on the unlensed E-modes and φ. Note that the dependence on the unlensed E-modes is linear, and where φ is uncontracted, the unlensed field E is implied. The one other possible trispectrum coupling is subdominant. Ultimately, our justification of this is the good agreement between our model for the delensed B-mode power spectrum and our simulation results. However, a plausibility argument can be made based on the volume of wavevector space that terms can accumulate in the integral in eq. (4.6), similar to the reasoning why the primary coupling should dominate over other "N (1) " couplings in the auto-power spectrum in CMB lensing reconstruction [42]. Generally, terms that couple together two or more of the weights [W (l a , l b )] in the integrand, i.e., the least factorisable terms, will be subdominant relative to other terms where the weights are uncoupled and some of the nested integrals can be separated, since the volume of wavevector space is reduced in the former case. In the specific case of eq. (A.3), the primary coupling produces the same pattern of weights as appears in the rest of the integrand in eq. (4.6), while for the other coupling there is no such factorisation. The primary coupling is responsible for the removal of the lensing signal in B obs that correlates with B temp ; that is, the actual delensing.
We now consider contributions from the six-point function B temp (l 1 )B temp (l 2 ) . The only fully-disconnected term that appears in the standard calculation is Typically, the fiducial observational power spectra used to inverse-variance weight the CMB fields for lensing reconstruction, and in the Wiener filters, will be calibrated from the observed power so to an excellent approximation C EE,obs,LAT l ≈ C EE,obs,fid,LAT l , and similarly for the B-mode spectra. In this case, the integral over l 1 , combined with the (fiducial) normalisation (A EB |l 1 −l 1 | ) 2 , gives the Gaussian reconstruction noise N The only other diagram that is retained in the standard calculation involves a Gaussian correlation of E-modes across templates multiplying a trispectrum made up of two quadratic estimators, φφ c . This trispectrum has been studied in detail by ref. [29]. To O(C φφ L ), it evaluates to a sum of the lensing power spectrum, from the primary coupling, and an N (1) term from the other couplings. Only the primary coupling is included in the standard calculation, and arises from contractions of the form However, we include N (1) here for completeness to find where we have assumed, again, that W E l C EE,obs,LAT l ≈ C EE,fid l . The explicit form of N (1)EB |l 1 −l | is given in eq. (57) of ref. [29]. We find this contribution to be subdominant on the relevant scales to both (A.5) and to the contribution from the primary coupling of the trispectrum, which gives rise to C φφ l in the integrand of (A.7). Combining eqs. (A.5) and (A.6) gives where in the last line we have neglected N ]. This easily follows from noting that all significant bias terms arise from contracting the pair of observed E-modes in at least one lensing template, that is they involve B temp (l) E obs,LAT when attempting to delens B-modes at wavevector l. The result of this contraction is proportional to the observed B-mode at l, B obs,LAT (l), used in the lensing reconstruction (see eq. A.13). If such modes are excluded from the reconstruction, the bias necessarily vanishes.
In the case of a single survey, the correlated noise between the B-modes used in the lensing reconstruction and the B-modes to delens sources a larger bias on the delensed spectrum than in the case of independent surveys (although, as explained in section 4.2, the signal-to-noise on primordial B-mode power is also greater in the former configuration), and acts to emulate an apparent, but spurious, delensing efficiency much greater than expected. This is particularly clear from the cross-spectrum of the template with the observed B-modes. Combining eq. (A.17) -after identifying D l with C BB,res l /C BB,obs,fid,LAT l in the limit where E-mode noise can be neglected in the template (see section 4.2) -with the unbiased result eq. (A.3), we find that this cross-spectrum produces the entire lensing power: irrespective of the actual fidelity of the lensing reconstruction. In figure 14, we quantify further this apparent delensing. We show the difference of the power spectrum of the observed SAT B-modes and the power spectrum of their delensed counterparts as a fraction of the B-mode lensing power as the instrument noise level is varied. This apparent delensing efficiency is shown averaged over degree-scale multipoles, for r input = 0, keeping otherwise the same simulation and reconstruction parameters as in figure 1. It is shown for the single-survey case in blue, while the black and green curves show the unbiased case (C W l /C BB l ) and biased case with N X l = 0, respectively, for comparison. The single-survey apparent efficiency is broken down into the contribution from the four-point function of the observed fields (i.e., twice the crosscorrelation between the template and the observed B-modes; shown in orange) and the six-point function (i.e., the auto-power spectrum of the template; red). The latter contribution does not depend on N X l and so is the same whether the surveys are independent or not (provided N BB,SAT l = N BB,LAT l ), but the contribution of the cross-spectrum is boosted in the single-survey case. In this case, at low noise levels both cross-and auto-spectra are conspicuously close tõ C BB l , resulting in an apparent delensing efficiency close to 100 %. For higher noise levels, the apparent delensing efficiency is even larger, exceeding 100 %.

B Simulations
We simulate observations of the CMB sky using the publicly-available code QuickLens 8 for a fiducial cosmology best fitting the Planck+WP+highL data of ref. [44]. Inspired by the experimental configurations of the upcoming Simons Observatory, we simulate a reconstructionoriented large-aperture-telescope (LAT) survey with noise levels ∆ P = 6 √ 2 µK arcmin = √ 2 ∆ T and beam FWHM of θ FWHM = 1.5 arcmin, together with a small-aperture-telescope (SAT) survey with ∆ P = 2 √ 2 µK arcmin = √ 2 ∆ T and θ FWHM = 17 arcmin. The simulations we generate are on the flat sky with 1024 pixels per side with an inter-pixel separation of 2 arcmin, covering approximately 2.8 % of the sky. Our simulations have periodic boundary conditions and are free of foregrounds. Note that we only simulate the part of the LAT survey that overlaps with the SAT survey.
The procedure for obtaining lensed CMB maps is as follows: first, we generate unlensed T, E, B and φ fields in harmonic space by drawing their Fourier coefficients from zero-mean Gaussian distributions with variance at each (interpolated) angular scale given by the theory  ). The internal delensing bias is larger for a single survey (blue) than for independent SAT and LAT noise (green), in both cases artificially inflating the apparent delensing efficiency compared to the unbiased case (black). The orange line shows the contribution to the single-survey case (i.e., the blue line) from the cross-correlation between the template and the observed B-modes, which is almost exactly equal to the full lensing power (see eq. A. 19). The red line shows the contribution from the auto-power spectrum of the template. power spectra (obtained from CAMB [45]). At this stage, we choose a pixelisation, which in turn sets the maximum frequency we can adequately sample by the Nyquist-Shannon sampling theorem. Then, the unlensed T, E and B fields are rotated into T, Q and U and remapped according to the deflection field d(n) = ∇φ as T = T [n + d(n)], (B.1) and analogously for Q and U , using a bivariate spline interpolation over a rectangular mesh. In order to mimic the effect of observations, we convolve the lensed fields with the transfer function for an experimental beam that is assumed to be symmetric and Gaussian with the required FWHM. As a final step, we add uncorrelated, Gaussian-distributed experimental noise at the map level.