CMB $\mu T$ cross-correlations as a probe of PBH scenarios

We propose a new method for probing inflationary models of primordial black hole (PBH) production, using only CMB physics at relatively large scales. In PBH scenarios, the primordial power spectrum profile for curvature perturbations is characterized by a pronounced dip, followed by a rapid growth towards small scales, leading to a peak responsible for PBH formation. We focus on scales around the dip that are well separated from the peak to analytically compute expressions for the curvature power spectrum and bispectrum. The size of the squeezed bispectrum is enhanced at the position of the dip, and it acquires a characteristic scale dependence that can be probed by cross-correlating CMB $\mu$-distortions and temperature fluctuations. We quantitatively study the properties of such cross-correlations and how they depend on the underlying model, discussing how they can be tested by the next generation of CMB $\mu$-distortion experiments. This method allows one to experimentally probe inflationary PBH scenarios using well-understood CMB physics, without considering non-linearities associated with PBH formation and evolution.


Introduction
Gravitational waves from merging black holes, in conjunction with other probes, can be used for testing our understanding of cosmology. There is the intriguing possibility that part of black holes from merging events have primordial origin, arising from the collapse of overdense regions in the early stages of our universe evolution [1][2][3]. The existence of such overdense regions might be attributed to the dynamics of inflationary cosmology 1 [6,7]. In order for producing primordial black holes (PBH) with astrophysically relevant masses, the spectrum of primordial curvature perturbations at small scales should increase by several orders of magnitude with respect to its values at large, CMB scales. Within the context of single-field inflation, such an amplification can be realized during a short non-attractor period, when the would-be decaying mode influences the evolution of super-horizon fluctuations. In this set-up, a peak in the curvature perturbation spectrum can be produced, and the non-linear process of PBH production and subsequent evolution depends on the details of the peak, as first explicitly discussed in [8]. We refer the reader to the reviews [9][10][11][12] for details and references on these topics.
In this work, we ask whether we can probe the PBH formation mechanism by testing larger scales well away from the peak of the spectrum. In fact, a close examination of the profile of the spectrum in inflationary scenarios capable of producing PBHs shows universal features, that are interesting to investigate. In particular, the spectrum is typically characterized by a very pronounced dip, followed by a rapid growth with a well defined slope towards its peak (see e.g. [13]). The position of the dip is not random, but depends on other global properties of the spectrum, as well as on details of the non-attractor inflationary evolution [14]. It typically occurs at much larger scales than the peak, k dip /k peak 10 −2 . Since such a dip seems to be a rather universal feature of power spectrum in single-field models of PBH formation, probing further the physics associated with it offers an indirect way of testing these models, independently on the details of non-linear PBH formation and evolution mechanisms, at scales much larger than k peak .
We start our analysis in Section 2, by analytically computing some of the key properties of the power spectrum and bispectrum of curvature fluctuations around the dip. For this purpose, we make use of the gradient expansion formalism introduced in [15], which allows one to compute the evolution of fluctuations on super-horizon scales in models with non-attractor phases [16]. Among our new findings, we show that the size of the bispectrum is generally amplified at the scale of the dip position. It has a rich dependence on momenta, with a broad support spanning different bispectrum shapes. When focussing on isosceles and squeezed configurations, it acquires a characteristic momentum dependence that is controlled by the underlying inflationary mechanism.
In Section 3, we then propose a potential probe of such a scale dependent bispectrum around the dip position by utilizing cross-correlations between CMB µ-type distortions and temperature anisotropies at large scales. In fact, µT cross correlations are known to be sensitive to primordial non-Gaussianity [17], hence they represent an appropriate observable for testing the features associated with it around the location of the dip. We show quantitatively how such cross-correlations are influenced by the characteristic scale-dependent bispectrum in isosceles configurations, and we estimate how our results depend on the underlying non-attractor evolution during inflation. We also discuss the prospects for detectability of this signal with future µ-distortion experiments, showing that µT cross-correlations constitute a promising avenue to test PBH scenarios with clean and well-understood CMB physics only. In Section 4 we conclude with future directions, and we then include four technical appendixes.
2 Curvature perturbation at super-horizon scales Models of primordial black hole (PBH) formation based on single field inflation often include a short phase of non-attractor evolution, characterized by a transient growth of would be the decaying mode which influences the super-horizon evolution of cosmological scalar perturbations 2 . In this section, we briefly describe a convenient formalism to study the statistics of curvature fluctuations in such situations. We aim to show that both the power spectrum and the bispectrum are characterized by pronounced features, whose phenomenological implications are then analyzed in Section 3.
We assume an isotropic, conformally flat FLRW background metric where a(τ ) is the scale factor, connecting conformal to physical time through the relation a(τ ) dτ = dt. Within canonical single-field inflation, the dynamics of the curvature perturbation on comoving hypersurfaces, R k (τ ), is governed by the following equation in Fourier space [20] 1 where z = aφ/H is the so-called pump field, with φ the background inflaton profile and H =ȧ/a the Hubble parameter. Primes indicate derivatives along conformal time, and dots derivatives along physical time.
To study the dynamics of fluctuations at super-horizon scales (k → 0) we focus on the gradient expansion formalism introduced in [15]. We refer the reader to [16] for further developments and applications of this formalism to the analysis of power spectrum of scalar fluctuations that leads to PBH formation. In this framework, iterative solutions of (2.2) can be generated at a desired order in k → 0 expansion in terms of simple analytic functions describing the background evolution which then allow us to relate the late time τ = τ * (i.e. at the reheating surface) curvature perturbation R k (τ ) during inflation, to its value at an initial time around horizon exit τ = τ k in terms of a complex, k dependent coefficient: Once expanded up to second order, k 2 gradient expansion, the coefficient α k reads where we defined k dependent fractional velocity of the curvature perturbation as The full k dependence of the expression (2.4) on super-horizon scales is then encoded in v R (See Appendix A) and the functions D(τ k ), F (τ k ) which are given by the following nested integrals of the pump field z(τ ) = aφ/H (See [15,16] for further details): (2.7) Whenever the pump field increases with time -as in standard slow-roll inflation, where z ∝ a(τ )the functions D, F rapidly decrease to zero after horizon crossing (i.e. α k → 1), and the curvature perturbation in (2.3) settles to a constant shortly after horizon exit (R k (τ * ) R k (τ k )). On the contrary, in inflationary models containing phases of non-attractor evolution, z(τ ) transiently decreases and the functions D, F can grow and amplify the curvature perturbation (i.e. |α k | 1 in (2.3)) at super-horizon scales.
Making use of eq (2.3), the power spectrum for the curvature fluctuation R at late times can be related to the power spectrum evaluated at the horizon crossing via , and we split α k into its real and imaginary parts using v R = v R R + i v I R . Up to order O(k 2 ) in the gradient expansion, the real and the imaginary part of the enhancement factor are therefore given by (2.10) We assume that R k (τ k ) is a Gaussian random variable: nevertheless, the superhorizon evolution typically introduces non-linearities. In fact, we can go beyond the linear theory provided in eq. (2.3) to compute the bispectrum of the late time curvature perturbation R k (τ * ). For the purpose of deriving an analytic expression for the bispectrum, we adopt the following non-linear version for the curvature perturbation, derived in [21] (2.11) where the last term represents the non-linear contribution parametrized by the convolution of the Gaussian variable R k (τ k ). Defining the corresponding bispectrum as 12) and using (2.11), the bispectrum results with [21] 13) where permutations are among the three external wave-numbers. Correspondingly, we define the scale-dependent non-linearity parameter f NL as where P R (τ * , k) ≡ 2π 2 k 3 P R (τ * , k). Using (2.8), we then have Finally, plugging the bispectrum in (2.13) into (2.15), we re-write the scale-dependent f NL as Notice from (2.16) that the size and the scale-dependence of the f NL parameter depends on the function F (τ k ) and the quantity α k , whose behavior depend on the background dynamics during inflation and in particular for the case of interest on the properties of the non-attractor regime. As a consequence, the corresponding bispectrum shape and its scale-dependence can be very rich, as we will learn in what follows. Armed with these preliminary, general results, we now study the scale dependence of the power spectrum and bispectrum on a representative set-up capable of producing a large PBH population during inflation.

2.1
The spectral shape of the spectrum: a dip is followed by a rapid growth In order to study the evolution and enhancement in the power spectrum, we consider a representative scenario that instantly connects an initial slow-roll era, with η sr = 0, to a slow-roll violating, non-attractor phase with constant η c ≤ −6 where η denotes the second slow-roll parameter, H η ≡ d ln /dt and the first slow-roll parameter is given by H = −d ln H/dt. In this setup, the pump field z(τ ) is assumed to have a profile: describing collectively the slow-roll and the constant-roll phases, the latter being controlled by the negative O(1) "slow-roll" parameter η c . We define τ 0 as the transition time to the constant-roll era, τ f as the conformal time when the constant-roll era ends, while we relate the quantity z 0 with a constant slow-roll parameter sr via z 0 = −a(τ 0 ) √ 2 sr M pl . For simplicity we describe the scale factor as a = −1/(Hτ ) with a constant Hubble rate H during inflation. We also indicate with H 0 the size of the comoving horizon at the time of the transition to the non-attractor era, and with ∆N = ln (τ 0 /τ f ) the duration (in e-fold numbers) of the non-attractor phase.
We proceed with determining the corresponding growth rate of the power spectrum. For this purpose, we re-write equation (2.8) as and we evaluate the power spectrum at τ * → τ f , i.e. at the end of the non-attractor era. Using equation (2.18), as well as the expressions (2.9) and (2.10), we make use of the analytic formulas for α R k and α I k from Appendix A and B to characterize the shape of the power spectrum. In this way, we plot our results in Figure 1 for two different sets of parameters characterizing the instant transition in the pump field profile of (2.17). We note that for studying physical implications of our findings, it is convenient to introduce a fixed quantity which determines the size of a mode k with respect to the horizon (aH) −1 at time τ = τ k , corresponding to the horizon crossing epoch. We then distinguish modes whose momenta lie in the following ranges: i) modes leaving horizon during the initial slow-roll era, i.e. modes satisfying τ k /τ 0 > 1 or equivalently k/H 0 < c k ≤ 1, and ii) modes that leave the horizon during the non-attractor η c ≤ −6 phase, c k < k/H 0 .
In these regimes, the behaviour of late time P R shown in Figure 1 reflects the accuracy of the gradient expansion formalism in capturing characteristic features of the power spectrum in inflationary backgrounds that transiently violates slow-roll conditions which is required to generate PBHs [22]. In particular, we notice a pronounced dip feature occurring at relatively large scales, associated with modes that still leave the horizon in the initial slow-roll era (i.e. k dip H 0 ) far away from the peak, which occurs at k peak 3H 0 . It is worth mentioning that such dip feature is due to competing contributions that appear in the power spectrum that are weighted by opposite signs. In particular, P R (τ * ) initially decreases going from large to small scales, because negative terms in eq (2.4) dominate at small k. Positive contributions to late time P R become instead more and more important at smaller scales (i.e. as k increase), where the power spectrum starts to rapidly grow, smoothly connecting with the spectrum in the regime of non-attractor evolution. To sum up, a dip feature in the power spectrum forms around the range of scales where negative and positive contributions to P R have comparable magnitude.
Interestingly, the presence of such a pronounced dip in the spectrum is a universal feature, being virtually present in all single field models based on non-attractor evolution that are aiming Figure 1: Scale dependence of the power spectrum for an example ultra slow-roll model (η c = −6) (Left) and constant-roll model (η c = −6.8)(Right) as function of k/H 0 . The vertical gray-dotted lines separate the interval of wave-numbers leaving the horizon during the initial slow-roll era (left side) and the non-attractor era (right side). We define A s = H 2 /(8π 2 sr M 2 pl ).
to generate a sizeable peak in the power spectrum for producing PBH, say of order ∆P R /P R 10 7 (see e.g. [23][24][25][26][27][28][29][30][31][32]). As can be noticed from Figure 1, following the scales corresponding to the dip feature, the power spectrum grows with a spectral index that can be as large as n s − 1 = 4 [13] towards the peak 3 . Then, for modes that leave the horizon deep in the non-attractor era, it decays with a spectral index n s − 1 = −|6 + η c | 4 . The properties of the dip turn out to be related with global features of the spectrum profile, see e.g. [14]. In fact, [14] found that the location of the dip position k dip in momentum space is related with a characteristic pivot scale k associated with the duration of non-attractor evolution by the expression k dip /k O(1) × (∆P R /P R ) −1/4 . For the scenarios we consider with a pronounced peak in the power spectrum, we can expect k to be related with k peak , and (∆P R /P R ) −1/4 10 7 −1/4 10 −2 . A close examination of Figure 1 confirms these arguments and shows a robust relation between the peak scale and the location of the dip feature in momentum space, obeying k dip 10 −2 k peak . (2.20) In light of these considerations, we assume (2.20) for the rest of this work. Since the slope of the spectrum changes abruptly at the dip position, it is interesting to ask whether we can make use of this characteristic feature for designing observables for probing the shape of the spectrum. This is our aim for what comes next.

The scale dependence of the bispectrum
Let us concentrate on modes that exit during the initial slow-roll era, k/H 0 < c k , to investigate the scale-dependence of the bispectrum. We expect the bispectrum to have features and be amplified around the scale k dip corresponding to the position of the dip in the power spectrum: in fact, non-linearities are usually enhanced at the location of rapid changes in the power spectrum (see e.g. [39] for a review). In what follows, we confirm this expectation for the system under consideration.
For scales satisfying k/H 0 < c k (recall the definition of c k around eq (2.19)), the power spectrum at around horizon crossing P R (τ k ) is scale invariant, see the discussion in Appendix A, in particular eq (A.8). We can then simplify the expression for the bispectrum for modes that exit during the initial slow-roll era as (we identify P R (τ k ) = P (0) R to emphasize its scale independence) The scale dependent non-linearity parameter 5 then reads as Notice that in eq (2.21), the scale dependence of bispectrum is characterized by the non-linearity parameter f NL , and the terms in the square brackets weighted by the enhancement factors α k contain information on the overall amplification of the power spectrum at late times. In this sense, Eq (2.21) can be interpreted to describe the late time bispectrum in terms of the power spectrum at late times which exhibits an amplification parametrized by α k . The equilateral limit of (2.21) makes these arguments clear, since we can re-write the late time bispectrum as (here (2.23) Nevertheless, we found that other shapes of the bispectrum turn out to be physically interesting.
In fact, the bispectrum in eq (2.21) has a broad support covering different shapes, and our analytical formulas allow us to have analytic control on the various possibilities. For definiteness, here we discuss the representative equilateral and squeezed configurations, while in Section 3.4 we will analyze other possible shapes corresponding to isosceles triangles in momentum space. We now focus on a representative choice of parameters that can generate an enhancement of order 10 7 in the power spectrum -as required for PBH production. Focusing on such scenarios, we plot in Figure 2 the resulting f NL in the equilateral k 1 = k 2 = k 3 and squeezed configuration k 3 k 1 = k 2 , for scales that exit the horizon during the initial slow-roll era, namely for scales around the k dip . From the left panel in Figure 2, we notice that the non-linearity parameter reaches its maximal value at the location of the k dip . The right panel in Figure 2 also provides valuable information on the scale dependence of f NL in the squeezed limit. We learn that this quantity initially grows positive, to then rapidly decrease to negative values around k dip , and to finally increase again from negative values, almost in a symmetrical fashion. Importantly, as can be seen from the overlap of orange dotted curve with the black solid one in the right panel of Figure 2, the scale dependence of f NL in the squeezed limit agrees very well with the one inferred from Maldacena's consistency condition: f NL = 5(1 − n s )/12. The typical magnitude for the squeezed limit of |f NL | around k dip results of order O(10), compatible with the large value for the spectral tilt around the dip region (see e.g. [16]).
It is worth to point out that the findings we present in Figure 2 agree well with previous literature [38,43,44] that numerically studied non-Gaussianity in inflationary setups that include a transient non-attractor era. However, the main advantage of our formulas (2.18), (2.21) and (2.22) is the analytic control they provide us to study features in the power and bispectrum in a model independent way without restoring to numerical techniques. In particular, these formulas are flexible enough to ensure an accurate description of the power and bispectrum in various inflationary backgrounds that exhibit a transient, slow-roll violating phase characterized by essentially two numbers: the duration ∆N of the non-attractor era and the value of the slow-roll parameter η c < 0 in this phase.

Summary of the results so far, and their implications
Let us briefly summarize our findings so far, and anticipate their implications to be developed next. We consider inflationary scenarios that incorporate a short non-attractor phase aimed to increase the size of the curvature perturbation spectrum at small scales. The spectrum profile is characterized by features, including a dip that precedes a steady growth -see Figure  1 -a phenomenon well studied in previous works. Moreover, for the first time we analytically demonstrate that the bispectrum associated with R acquires pronounced features and strong scale dependence around k dip (see Figure 2).
Can we use the pronounced features associated with the statistics of curvature perturbation at the scale k dip , occurring well away from the peak (see eq (2.20)), for probing these inflationary models with CMB physics only, independently on the details of PBH formation? An affirmative answer would allow us to probe the PBH generation mechanism within clean and nearly-linear perturbation theory at relatively large scales, avoiding complications involved with the physics of PBH formation.
With this aim in mind, we explore the idea to study cross-correlations between temperature fluctuations at large CMB scales, of order k L 10 −4 −10 −2 Mpc −1 , and CMB µ-spectral distortions at smaller scales 6 , say k S 10 3 − 10 4 Mpc −1 . Correlations between such disparate scales can indeed be induced by the squeezed limit of the curvature bispectrum [17], which in our case is amplified at the dip position to the levels of order |f NL | O (10). Given that for scales k peak < 10 5 Mpc −1 a peak in the scalar power spectrum is expected to be tightly constrained from µ distortions alone (see e.g. [13,46]), we focus on the interval k peak = 10 5 − 10 6 Mpc −1 . This interval corresponds to k dip = 10 3 − 10 4 Mpc −1 (see eq (2.20)), the range most relevant for µ distortions, as we review in the next section. We point out that selecting k peak = 10 5 − 10 6 Mpc −1 leads to the formation of PBHs in the range M PBH 1 − 100 M (see e.g. [11,47]). PBHs within this mass range are interesting in their own, and moreover can be considered as seeds for Supermassive Black Holes (SMBH) with M SMBH = 10 6 − 10 8 M , taking into account accretion and merging effects after production [46,48]. However, as one can anticipate from the discussion above, we will not consider any actual details of the PBH production, and in the next section we will concentrate on how a squeezed bispectrum enhances cross-correlations between temperature fluctuationsand µ-type distortions at scales that are much larger than k peak .
We conclude this subsection with some comments on the physical significance of the squeezed limit of the bispectrum in single field inflation. There has been some debate in the recent literature on whether such limit is physical, or whether can be set to zero with a gauge transformation. Considering single-field inflationary scenarios in a non-attractor phase, [49] finds that the physical f NL vanishes. On the other hand, [50] (see also the general discussion in [51,52]) points out that the dynamics of the decaying mode can leave important imprints in the squeezed limit of the bispectrum, when carefully considering the process of matching between non-attractor and attractor regimes towards the end of inflation 7 .
The effects that we are going to study hold whenever there is a coupling between long and short modes induced by non-Gaussianity, and does not require to focus on the infinitely squeezed limit of the bispectrum where the previous debate applies. In our case, as we have seen with examples in Section 2.2, the bispectrum shape is quite complex and has a broader support than the strictly local shape. Hence it can provide such couplings even outside the purely local Ansatz, and we will make use of this fact for our phenomenological considerations in Section 3.4. Moreover, although here we consider single-field systems, the same phenomenon can also occur in multiple field or curvaton-like scenarios for PBH production with pronounced features in the statistics of curvature fluctations. In those cases, the effects of non-adiabatic modes in the squeezed limit of the bispectrum can not be removed by gauge transformations. Hence, in what follows we do not discuss these issues any further, and focus on developing the phenomenological consequences of our idea.

µT correlations as a probe of the PBH generation mechanism
We now discuss how to use the cross-correlation among CMB spectral distortions and temperature fluctuations for probing the statistics of curvature perturbation in inflationary scenarios capable of generating PBHs. We start with a brief review that discuss µ-type CMB distortions and their properties for testing the statistics of primordial fluctuations following previous works. We then apply these methods to inflationary scenarios leading to PBH production, analyzing the prospect of using µT cross correlation for probing the corresponding scale-dependent bispectrum at the dip position. The method we propose allows one to indirectly probe the PBH formation properties at relatively large scales, away from the scales where the curvature spectrum grow and has a peak.

Brief review of CMB µ-distortions
The energy injection caused by the diffusion damping (Silk damping) of the acoustic waves in the pre-recombination photon-baryon plasma heats photons and leads to spectral distortions in the black-body spectrum of the CMB as they re-enter the horizon [59,60]. At very high redshifts, these distortions are erased by both photon conserving (i.e. Compton scattering) and non-conserving (double Compton scattering) processes. For the range of redshifts z f ≡ 5×10 4 < z < 2×10 6 ≡ z i , photon non-conserving processes cease to be efficient but photons can still maintain thermal equilibrium by elastic Compton scattering e − + γ → e − + γ, which conserves photon number. The resulting photon spectrum is then described by a Bose-Einstein distribution with a non-vanishing, chemical potentialμ ≡ µ/(k B T ), leading to the so-called µ-type distortion of the black-body CMB spectrum: where n(ν) is the number density of photons per frequency. Being associated with the dissipation induced in the primordial plasma by the super-horizon primordial fluctuations as they re-enter the horizon and starts oscillating, µ distortions have primordial origin, and thus can be expressed in terms of the primordial power spectrum [17,[61][62][63]. To see this, we relate the size of µ distortion to the heat generated by diffusion damping Q γ [64]: where the square of the sound speed obeys c 2 s 1/3 since the universe is radiation dominated at those redshifts, and . . . p is the square of the photon energy density contrast δ γ averaged over a period of acoustic oscillations. This energy release is then converted into µ-distortions as Denoting the transfer function of photon density contrast δ γ (τ 0 , k) = T γ (k)R k -with R k is the conserved super-horizon curvature perturbation -µ-distortions are then related to the primordial curvature perturbation by 17] is the small-scale limit of the linear photon transfer function and k D denotes the diffusion damping scale. During radiation domination, denoting R ≡ 3ρ b /ρ γ 1 (ρ b being the baryon density), the damping scale reads as The quantity k D therefore depends on redshift, and will appear in many formulas in what follows. In (3.4), W (k) = 3k −3 [sin(k) − k cos(k)] is top-hat filter function in Fourier space that smears the dissipated energy over a volume of radius k −1 s k D (z f ) −1 and k ± = k 1 ± k 2 . Using (3.4), the average (monopole) µ−distortion in the CMB generated by acoustic damping of perturbations in the photon-baryon plasma is then given by the log integral of the primordial power spectrum from k D (z f ) 46 Mpc −1 to k D (z i ) 1.2 × 10 4 Mpc −1 : Eq (3.6), shows how µ−type distortions 8 , generated by the dissipation of the photon density perturbation at small scales, allow one to probe the primordial power spectrum at small scales. In particular, µ-distortions can be directly used for constraining inflationary models generating PBHs using a pronounced peak in their curvature power spectrum at scales of order k peak < 10 5 Mpc −1 , see e.g. the recent [13,46]. But as anticipated in Section 2.3, we can use spectral distortions to indirectly probe also smaller scales k peak > 10 5 Mpc −1 . Our idea is to use the large non-Gaussianity produced at the dip position k dip k peak /100 (see Figure 2), which induces large cross-correlations of µ-distortions with the CMB temperature spectrum. As we show in what comes next, the scale-dependence of the bispectrum around the dip feature leads to larger effects compared to the inflationary scenarios that support a purely local bispectrum, making such cross-correlations a valuable observable for testing mechanisms of PBH generation independently from the details of the PBH production.

Cross-correlation of µ-distortions with CMB temperature anisotropies
We now briefly review how primordial non-Gaussianity induces cross-correlation between µdistortions and CMB temperature anisotropies, closely following [17,66]. We then discuss and compare the resulting µT correlation focusing on two different cases: i) for a primordial scenario that can generate a purely local type bispectrum and ii) for a slow-roll violating inflationary model that features a scale dependent squeezed bispectrum around the dip feature as we discussed in Sections 2.1 and 2.2.
As first shown in [17], a non-zero bispectrum in the squeezed limit makes the distribution of µ anisotropic on the sky and allows cross-correlations between µ anisotropy and CMB temperature anisotropy Θ(n) = δT (n)/T . To see this, we can expand the CMB temperature anisotropies detected by an observer into spherical harmonics as Θ(n) = lm a T lm Y lm (n), where the harmonic coefficients are given by and ∆ l is the radiation transfer function. Similarly, for an observer at the origin, direction dependent distortion anisotropies in µ(n) can be expanded into spherical harmonics as µ(n) = lm a µ lm Y lm (n). Using eq. (3.4), we can then relate the the spherical harmonic coefficients to the primordial curvature perturbation as [17,66] where χ * = τ 0 − τ * 14 Gpc is the comoving distance between the last scattering surface and today. Defining the angular correlators of two anisotropies labeled by {i, j} as The cross correlation of (3.8) and (3.7) then gives (3.10) In this expression we choose the smoothing scale to be the same as the largest scale µ distortions are relevant for, i.e. k s ≈ k D (z f ) and focused on the squeezed limit k + = k 3 → 0, k 1 → k − /2 together with the definition of the bispectrum in (2.12). It is clear from (3.10) that the angular correlator C µT l is sensitive to the integral of the bispectrum in the squeezed configuration, in particular with a larger ratio k 1 /k 3 k − /2k + compared to the CMB anisotropies alone. Noting that the smallest wave-number we can probe from the CMB anisotropy in the sky corresponds to the quadrupole, k 3 ≈ l/χ * = 2/(14 Gpc) 1.4 × 10 −4 Mpc −1 , we have access to the ratio k 1 /k 3 within the following range, Note that this ratio is much greater one can reach only via temperature anisotropies of the CMB for l = 2 − 3000 where k 1 /k 3 = l 1 /l 3 = 1 − 1500. To set the stage for the computation of µT in inflationary scenarios that can produce large PBH populations in the post-inflationary universe (See Sections 2.1 and 2.2), below we first review the standard calculation for a purely local bispectrum.
3.3 C µT l for a local type bispectrum For local type non-Gaussianity, the bispectrum of the curvature perturbation in the squeezed limit can be expressed as [67] NL is the scale-independent primordial non-linearity parameter 10 . Adopting for simplicity a scale invariant primordial power spectrum P (0) R = 2.1 × 10 −9 , and plugging (3.12) into (3.10), we obtain 11 C µT where the multipole dependent quantity b (loc) (the suffix (loc) means that it holds for a local bispectrum) can be expressed in terms of a double integration over short and long momenta involving the transfer function ∆ l of temperature anisotropies and window function W as b (loc) (l) = 6 l(l + 1) To compute this quantity, one can first focus on the Sachs-Wolfe (SW) limit (small l) where ∆ l → j l (k + χ * )/3. Then realizing that k + k s k D (z f ) and W (k + /k s ) → 1 in the squeezed limit k + → 0, the integrals in (3.14) can be carried analytically to conclude b (loc) → 1 and arrive at the standard result 10 See [68,69] for earlier works on the study of cross correlations between µ and Θ in the presence of scale dependent fNL. 11 It should be noted that there are non-primordial contributions to the final observed C µT l from non-linear and projection effects which are found to be negligible [70]. function ∆ l (k + ) of temperature anisotropies in the integral (3.14) over long momenta k + , which can be implemented numerically [66]. It turns out that the final C µT l obtained in this way can be related to C µT,SW l by an analytic fit first proposed in [69]: The resulting C µT l in comparison with C µT,SW l is shown in Figure 3. As first shown in [66], the SW approximation breaks down around l = 10 and C µT l changes sign at l 40 becomes negative for l > 40. This implies that for a local type bispectrum with f (p) NL > 0, while temperature anisotropies and µ distortions are correlated at the largest scales l < 40, they become anti-correlated at small scales.
3.4 C µT l for inflationary scenarios with an enhanced power spectrum We next focus on inflationary models that contains a transient non-attractor era η c < 0 following the initial slow-roll phase, as discussed in Section 2. Using the gradient expansion method of Section 2, eqs (2.21) and (2.22), for modes that exit the horizon in the initial slow-roll era, we obtain the following analytic expression for the bispectrum associated with isosceles triangle configurations: where P R (k) = 2π 2 P (0) R /k 3 is the dimensionful power spectrum at around horizon crossing τ k during the initial slow-roll era. Notice that -as opposed to the parametrization we undertake in eq (2.14) -in (3.17) we express the bispectrum with respect to two copies of the power spectrum, evaluated at horizon crossing. This way of expressing the bispectrum provides an easier comparison with the local bispectrum case we reviewed earlier, since eq. (3.17) can be recast in a form that resembles the standard local type bispectrum for isoceles triangle configurations: where the effective non-linearity parameter is given by In appendix D, we study the squeezed limit k q of (3.19) and found that it quickly reaches to a small constant value f 2)) in the large scale tail of the squeezed limit, q → 0 and q > k. Since, we send q → 0 to reach this initial value, we do not expect that f (0) NL can reflect an accurate initial value of the squeezed configurations of the expression (3.19). To keep our discussion on the observability of µT (see Section 3.5) as general as possible, we will therefore normalize f eff NL with respect to this initial value (as in (D.2)) and rescale it with a fiducial primordial f (p) NL that parametrizes the effective local non-linearity parameter valid on large scales i.e. at CMB scales. Following appendix D, in the squeezed limit k q, we will therefore adopt the following where H 0 k peak /3 100k dip /3 (see Section 2.1) and the size of the coefficients c Notice that the parametrization (3.20) ensures f eff NL → f (p) NL in the large scale limit q → 0. On the other hand, when we approach to smaller scales with q k dip , the scale dependence of the bispectrum becomes important, and its features relevant around the k dip . To illustrate this, we present the scale dependence of f eff NL (q, q, k) (3.20) in the squeezed limit (k/q → 0) in Figure 4.
NL in the large-scale tail of the squeezed limit (q → 0 where q > k 3 ). More importantly, we observe that f eff NL changes sign at k dip and starts to grow large in absolute value in the negative direction.
We now turn our attention to the µT ((3.10)) cross correlation that arise through the scale dependence of such bispectrum around the scales k dip associated with the dip feature. For this purpose, noting the definition in (3.20), we utilize (3.18) in (3.10). We then obtain the angular µT cross correlation as where we defined the quantity b (pbh) , whose expression in the inflationary scenarios containing a non-attractor phase can be re-written as (3.22) We note that the expression above only assumes a bispectrum in the isosceles triangle configurations without specifying a specific triangle shape such as equilateral or squeezed configurations. This quantity b (pbh) involves a double integral over short and long momenta, and coincides with expression b (loc) of equation (3.14) if the effective non-linearity parameter does not exhibit scale dependence, i.e.f eff NL → 1. However, as we mentioned previously, for an inflationary scenario that can generate PBH populations, the effective f eff NL (q, q, k) (3.20) is scale-dependent especially around the dip feature present in the power spectrum (See e.g. Figure 4). This situation might lead to substantial differences between b (pbh) and b (loc) that can allow us to probe the underlying PBH production scenario.
In particular, we expect that these differences can appear for two reasons: i) The (possible) dependence off eff NL (q, q, k) on the small momentum k, corresponding to CMB temperature anisotropy scales. Since k also appears in the l-dependent functions ∆ l , j l in the first of the nested integrals in (3.22), such a dependence can potentially change the multipole dependence of b (pbh) and hence the resulting C µT l correlator.
ii) The dependence off eff NL on the momentum q, corresponding to scales probed by spectral distortions. Such a dependence would only modify the second of the nested integrals in eq (3.22), changing the overall magnitude of b (pbh) but not its l-dependence.
For what respects point i) above, in this work we focus on the bispectrum in the squeezed limit. Focusing on such configurations, we found that at leading order in the small parameter k/q 1, f eff NL is independent from the soft momenta k as can be verified from (3.20) (See also appendix D). This implies that, in the ultra-squeezed limit, the l-dependence of the resulting C µT l is not modified compared to the local bipsectrum case we studied earlier (see Section 3.4). Instead, for what respects to point ii), we found that the dependence off eff NL on the large momenta q is quite strong (See eq. (3.20)) such that it can enhance the overall magnitude of b (pbh) , improving the chances of testing the effects of non-attractor periods in models with PBH production.
Based on these facts, we expect that the l-dependence of eq. (3.22) will be the same as in the local bispectrum case. Therefore, the full multipole (l) dependence of b (pbh) (or C µT l (3.21)) can be obtained starting from a computation made in the large-scale SW limit which can then be extended to smaller scales by means of the analytic fit of eq (3.16) using b (pbh) (l) = ρ(l) b (SW) (pbh) . For this purpose, we first focus on the expression b (pbh) of (3.22) in the SW approximation, assuming ∆ l (k) → j l (kχ * )/3 in this limit. Using the leading order result (D.2) forf eff NL in the squeezed limit, the integrals in (pbh) with (3.16). We observe that compared to the scale independent local bispectrum case (See Figure 3), the overall amplitude of the angular correlator is enhanced significantly and its behavior from large to small scales is inverted due to the negativity b (SW) (pbh) for the k dip value quoted in Figure 6. In other words, for a PBH forming inflationary scenario, µ distortions become anti-correlated with temperature anisotropies at large scales contrary to the case arise for a scale independent local bispectrum. It is worth stressing that this behavior arise because in PBH forming inflationary scenarios, scale dependent bispectrum changes sign around the dip feature and grows large in the negative direction as we show in the right panel of Figure 4. This result is interesting on its own as one can in principle distinguish between these scenarios by just looking at the sign of the C µT l at large scales, i.e. for small l. In summary, our results imply that phases of non-attractor inflation can qualitatively and quantitatively change the properties of the cross-correlations between CMB temperature anisotropies and µ-type distortions.

Prospects of detectability for µT
In order to asses the prospects of detectability for the µT correlator, we estimate the cumulative signal-to-noise ratio, S/N using [66] S N 2 = lmax l=2 (2l + 1) where C T T l is the CMB temperature anisotropy power spectrum and C µµ,N l is the noise level for µ distortions. For an experiment as PIXIE [71] one has C µµ,N l 4π µ 2 min e l 2 /84 2 [17] where µ min denotes the minimum detectable µ distortion (monopole) signal. Using (3.7), we denote the T T correlator in (3.24) as Following our discussion in the previous section, taking into account full transfer effects we express the µT cross correlation as with ρ(l) defined in (3.16). Then we plug (3.26) into (3.24) to re-write (S/N ) 2 ratio as where for µ min we take the PIXIE fiducial value 10 −8 [66], and the sum in (3.27) is carried up to l max = 200. To estimate the S/N , we would require the full knowledge of C T T l in (3.25) using the full transfer function ∆ l . Taking these full transfer effects into account, [66,69] found that the signal to noise ratio can be estimated as 40% of the result that one would obtain by adopting the SW limit ρ l → 1, C T T,SW l = 2πP (0) R /(25 l(l + 1)) 13 in (3.27). Keeping this in mind, we obtain This result implies that |f (pbh) | 2892 should be observable for a PIXIE like experiment. On the other hand, for a spectrometer comparable to PRISM [72] with µ min = 10 −9 , |f (pbh) | 290 is required for the detectability. At this point, it is worth mentioning that for a local type f (p) NL ∼ 10 −2 at CMB scales as predicted by standard slow-roll backgrounds, such a signal would be challenging to detect. However, PBH forming inflationary models can in principle generate a more pronounced dip feature in the power spectrum (and hence in the bispectrum) (See e.g. [28,30]) compared to dip features we can capture using analytic formulas we study in this work. Based on our explorations in this paper, we anticipate that for such dramatic dip features, the overall magnitude of b (SW) (pbh) can be in principle enhanced to compensate for the f (p) NL ∼ 0.01 at CMB scales. A detailed analysis on the expected b (pbh) in the aforementioned models require numerical work and is outside the scope of this work.
On the other hand, considering the 2σ limits on the local type bispectrum by Planck [73] at CMB scales: −11.1 < f (p) NL < 9.3 (f NL = −0.9 ± 5.1, 68%CL) and the typical values of |b (SW) (pbh) | a few (10 2 − 10 3 ) that can be obtained in inflationary scenarios that can generate PBH populations (See figure 5), we conclude that the influence of the scale dependent non-Gaussianity (see Section 3.4) on C µT l , as due as the PBH formation mechanism, could be observable by PIXIE or PRISM for an interesting range of k dip values.

Discussion
In this work, we presented a method for probing inflationary scenarios of PBH formation, using only CMB physics at relatively large scales. We based our arguments on the characteristic properties of the spectrum of curvature perturbation in single-field inflationary models that can generate a large population of PBHs. In these models, the curvature perturbation spectrum is characterized by a pronounced dip followed by a rapid growth towards a peak responsible for PBH formation. By making use of the gradient expansion formalism of [15], we analytically computed the properties of the power spectrum and (for the first time) of the bispectrum around the dip position, which occurs at scales well larger than the peak. The bispectrum turns out to have a rich dependence on momenta, with a broad support spanning different bispectrum shapes. We found that when focussing on isosceles and squeezed configurations, the bispectrum can be enhanced at the position of the dip; it also acquires a characteristic momentum dependence that is controlled by the underlying inflationary mechanism. We proposed to probe such enhanced squeezed bispectrum through the correlations it induces between CMB µ-distortions and CMB temperature fluctuations. We extended the methods first explored in [17] to include the case of scale-dependent non-Gaussianity from PBH formation mechanisms, finding analytical expression for quantities controlling µT correlations, and discussing how future CMB µ-distortion experiments can test this observable. Interestingly, the method we propose would allow one to experimentally probe inflationary PBH scenarios using well-understood CMB physics at scales much larger than the peak of the spectrum, without considering non-linearities associated with PBH formation and evolution. In particular, owing the relation between relevant scales associated with µ distortions and the features present in the PBH forming inflationary scenarios, our findings are relevant for PBH masses within the M (pbh) ∼ 1 − 100M . This implies that µT correlations can be considered as a useful tool to distinguish between astrophysical vs primordial origin of BHs.
Our work can be extended in several directions. From the phenomenological side, it would be interesting to explore more broadly the effects of non-Gaussianities in the region occurring around the dip, and in particular to investigate whether an enhanced trispectrum would lead to distinctive signals in the self-correlations of µ-distortions [17]. Moreover, the ideas we pursued in this work can be complemented with more direct methods for constraining the slope of the scalar power spectrum with µ-distortions [13,46], that rely on the knowledge of the growth rate of the spectrum towards the peak and hence control different ranges of scales. On the theoretical side, it would be interesting to have a better understanding of possible consistency relations connecting the features of the spectrum such as dips and peaks [14], and their consequences for the bispectrum. Finally, extending the ideas we presented in this work to the multiple scalar field case appear as another interesting venue to be explored, in particular to make a comparison with single-field results. Then a pronounced dip feature in the power spectrum might be missing [74][75][76][77] (but see [78]). Similarly, in models that utilize axion gauge-field dynamics during inflation to enhance the curvature perturbation, such a future is generically not present [79,80] (See however [81]). On the other hand, still much work is needed to clarify the properties of the spectrum in this context.

Acknowledgments
It is a pleasure to thank Enrico Pajer for comments on the manuscript. The work of OÖ is supported by the European Structural and Investment Funds and the Czech Ministry of Education, Youth and Sports (Project CoGraDS-CZ.02.1.01/0.0/0.0/15003/0000437). GT is partially funded by the STFC grant ST/T000813/1.
A The curvature perturbation R k and fractional velocity v R The enhancement factor α k in (2.9), (2.10) and the eqs. (2.8) and (2.13) implies that we need an expression for R(τ k ) and the fractional velocity v R in (2.5) for being able to determine the spectral behavior of power and bispectrum of R k . In this appendix, we therefore aim to derive an expression for R(τ k ) and v R for the two phase model we study in the main text. For this purpose, we resort to Mukhanov-Sasaki equation for the canonically normalized variable Q k (τ ) ≡ z(τ )R k (τ ), which is exact to all orders in slow-roll parameters. For constant values of slow-roll parameters , η, and assuming 1 (such that H is approximately constant as we assume), the exact solution for Q k in terms of the Hankel functions is given by where For the initial slow-roll phase (η sr = 0 → ν = 3/2), requiring that all modes are in their Bunch-Davies vacuum for −kτ → ∞ in (A.3), the curvature perturbation is then given by where we have used z = (−Hτ ) −1 √ 2 sr M pl . The solution above immediately implies This result makes it clear why the curvature perturbation settles to a constant solution shortly after horizon exit in standard slow-roll inflation, which can be understood in the −kτ → 0 limit of eq. (A.6). For our purposes, we are interested in the expression in (A.6) evaluated at the initial time τ = τ k at around horizon crossing. With this in mind, we split the fractional velocity into a real and imaginary part as, where we defined a positive number −kτ k = k/H k ≡ c k ≤ 1 to identify the size of the each mode with respect to the horizon size at the initial time, i.e. at τ = τ k . It is clear from this expression that the imaginary part of v R includes an extra factor of c k compared to the real part. We note that unless c k = 1, this translates into an extra suppression for the imaginary part of the fractional velocity.
On the other hand, using (A.5), the power spectrum evaluated at around horizon crossing is given by Next, we need to develop an expression for R(τ k ) and the fractional velocity continuous through the transition at τ = τ 0 . For this purpose, we use a matching procedure for R k and its derivative between the initial slow-roll era, i.e. (A.5) to a general solution during the non-attractor stage: where ν = (3 + η c )/2. Matching R k and R k at the transition τ = τ 0 using (A.5) and (A.9), we obtain , (A.10) where we defined y ≡ −kτ . Then for τ k /τ 0 < 1, the real and the imaginary part of the fractional velocity is given by where we define the functions f α = f α (y, y 0 , ν), α = 1, 2, 3, 4 in terms of the Bessel function of the first and second kind as noting f 4 = f 1 (y, y 0 , ν + 1), f 3 = −f 2 (y 0 , y, ν). Using these expressions, power spectrum evaluated at τ = τ k for modes that leave the horizon during the non-attractor phase is given by The continuity of the real and the imaginary part of the fractional velocity can be confirmed explicitly from the eqs in (A.12) and (A.13) as they reduce to their slow-roll counterparts provided in (A.7) at the transition point τ k = τ 0 .
In this appendix, we present the details on the calculation of the integrals associated with the functions D(τ k ),F (τ k ), in the background model given in eq. (2.17). For this purpose, we define x ≡ τ /τ 0 to re-parametrize the background pump field as for modes leaving the horizon during the slow-roll era, the functions D(τ k ) and F (τ k ) were calculated in [16] and found to be , where τ f /τ 0 ≡ x f = e −∆N with ∆N denoting the duration of the non-attractor era. Next we focus on D(τ k ) and F (τ k ) for modes that leave the horizon during the non-attractor era, i.e. x k < 1. In this case, integrals defined in (2.6) and (2.7) are much simpler to evaluate which require only the behavior of the pump field during the non-attractor era, i.e. the second line in (B.1). Proceeding in this way, for modes that exit during the non-attractor era, x k < 1, we obtain the following expressions Denoting x k = c k (k/H 0 ) −1 , the wave-number dependence of these functions can be parametrized as where the coefficients C,C are functions of the parameter set {c k , η c , ∆N } as can be understood from (B.2). Note that for modes that exit the horizon during the non-attractor era, the formulas in (B.5) are valid up to a maximum wave-number k/H 0 corresponding to the mode that exits the horizon at the end of the non-attractor era which obeys kτ f = 1.
It is worth emphasizing that the coefficients C andC that multiply the k dependent terms can be organized in a hierarchal way in powers (determined by η c ) of a(τ f )/a(τ 0 ) = e ∆N where ∆N is the duration of non-attractor era in number of e-folds. This result reflects the fact that modes that leave during the slow-roll and in the early stages of non-attractor era are enhanced due to the slow-roll violation η c ≤ −6.

C Consistency condition
In this appendix, our aim is to show that the formulas we derived for the non-linearity parameter in (2.16) reduce to the standard expression implied by the consistency condition for vanilla slow-roll models. For this purpose, we will focus on the mode equation (2.2) of the comoving curvature perturbation in Fourier space. It is a well known fact that in a standard slow-roll background, the growing mode solution of eq. (2.2) is conserved on super horizon scales. This can be readily seen from the formal integral solution of (2.2), which can be written up to order O(k 2 ) for small but finite wave-numbers as where we obtain the last term by solving iteratively the inhomogeneous part of eq. (2.2) using the leading growing mode which we identify as R k (τ k ) = R (0) . The constant behavior of R k shortly after its scale crosses the horizon can be readily seen from the solution (C.1), by realizing that -in a slow-roll background where z ∝ (−τ ) −1 -the second and the third term in (C.1) decay respectively as (−τ ) 3 and (−τ ) 2 in the late time limit −τ → 0. Therefore, in a slow-roll background we can immediately identify the second and third term in (C.1) as the decaying modes. In fact, the standard decaying mode is given by the last term as it decays slowly, i.e. ∝ (−kτ ) 2 , compared to the second. Notice that the second and the third term are proportional to the functions we defined as D(τ ) and F (τ ) within the gradient expansion formalism we are undertaken. In the language of the main text, the discussion above suggests that in an always slow-roll background, we can neglect the terms proportional to D(τ ) compared to the term proportional to F (τ ) for a sufficiently late time τ k shortly after horizon crossing. Therefore, the power spectrum at late times can be well approximated by where P (0) R = k 3 |R (0) | 2 /(2π 2 ) = H 2 /(8π 2 sr M 2 pl ) is the scale invariant power spectrum evaluated at τ k assuming a constant H with sr 1. From (C.2), we can relate spectral index to the k dependent part of the α k as where F (τ k ) = τ 2 k /6 in a slow-roll background with z ∝ (−τ ) −1 and used the fact that there is no enhancement in this scenario, i.e. |α k | → 1. Taking the squeezed limit k 1 k 2 = q k 3 of the expression in (2.16) then gives f NL (q, q, k 3 ) = 5 12 4F (τ q )q 2 α q α * where we used (C.3) by noting |α k 3 | 1 as k 3 → 0 and |α q | 1 for a slow-roll background as before. It is worth to point out that by virtue of the gradient expansion formalism we undertake, 1 − n s ≈ 2c 2 k /3 1 at leading order in the small parameter c 2 k with c k = −kτ k < 1 which is consistent with the identification that all the modes we consider are already outside the horizon at the initial time τ k . D f eff NL in the squeezed limit Here, we would like to derive the squeezed limit of the effective non-linearity parameter defined in (3.19) for modes that leave the horizon before the background transitions to the non-attractor regime. Using the formulas we derived for the functions D(τ k ), F (τ k ) and the definitions of the enhancement factor α k (See e.g. eqs. (2.9) and (2.10)), we can derive and expression for f eff NL (3.19) in terms of an expansion over the small ratio k/q 1 in the squeezed limit. In this way, up to next to leading order in k/q, we obtain f eff (D.1) where v R 's are defined as in (A.7) and the background dependent coefficients C F −2 = C F −2 (c k ) can be extracted from (B.2) and (B.4). Similarly, in order to obtain an explicit expression in terms of the hard momenta q, we further dissect the expression in (D.1) using the formulas for α q and F (τ q ) to re-write the effective non-linearity parameter as NL (c q )/3 is the "initial" value of the f eff NL we define in the large scale tail of the squeezed limit, i.e. q k as q → 0 and we used the fact that we universally utilize a small number for all the modes we study assuming c k = c q = constant < 1. Recall that the condition c k = −kτ k < 1 serve for the purpose to ensure that each mode labelled by a wave-number k are outside the horizon for an appropriately chosen time τ k . The coefficients c (n) NL that appear in (D.2) can be written in terms of the coefficients C of the functions F, D in Appendix B and the fractional velocities v R R , v I R we defined in Appendix A as where {η c , ∆N, c q } dependence of these expressions should be understood.