Cosmological Particle Production and Pairwise Hotspots on the CMB

Heavy particles with masses much bigger than the inflationary Hubble scale $H_*$, can get non-adiabatically pair produced during inflation through their couplings to the inflaton. If such couplings give rise to time-dependent masses for the heavy particles, then following their production, the heavy particles modify the curvature perturbation around their locations in a time-dependent and scale non-invariant manner. This results into a non-trivial spatial profile of the curvature perturbation that is preserved on superhorizon scales and eventually generates localized hot or cold spots on the CMB. We explore this phenomenon by studying the inflationary production of heavy scalars and derive the final temperature profile of the spots on the CMB by taking into account the subhorizon evolution, focusing in particular on the parameter space where pairwise hot spots (PHS) arise. When the heavy scalar has an $\mathcal{O}(1)$ coupling to the inflaton, we show that for an idealized situation where the dominant background to the PHS signal comes from the standard CMB fluctuations themselves, a simple position space search based on applying a temperature cut, can be sensitive to heavy particle masses $M_0/H_*\sim\mathcal{O}(100)$. The corresponding PHS signal also modifies the CMB power spectra and bispectra, although the corrections are below (outside) the sensitivity of current measurements (searches).


Introduction
Cosmic inflation (for a review see [1]) is one of the leading paradigms to explain the generation of primordial fluctuations that seed the cosmic microwave background (CMB) anisotropies and the large-scale structure (LSS) of the Universe. Similar to particle collider experiments, the inflationary era, characterized by the Hubble expansion parameter H * , generically provides an environment to produce new particles with masses ∼ H * which can in turn imprint their signatures on primordial fluctuations. Given H * during inflation can as high as 5 × 10 13 GeV [2], this provides us with a unique chance to probe very high-energy physics operating at scales H * .
During inflation, the expanding spacetime leads to a time-dependent Hamiltonian, and particles with masses m ∼ H * can be produced non-adiabatically even when starting with the Bunch-Davies vacuum at very early times [3]. While such minimal non-adiabatic production amplitude gets exponentially suppressed as e −πm/H * for m H * , the slowly rolling inflaton field itself can assist in producing even heavier particles if a direct coupling exists between the inflaton and the heavy particles. The produced particles can (a) generate a non-trivial spatial correlation between curvature perturbations ζ( x 1 )ζ( x 2 ) · · · at different locations or, (b) if sufficiently heavy, they can directly modify the curvature perturbation ζ( x) around their individual locations, as will be explained below. In the absence of isocurvature fluctuations, curvature perturbations ζ then remain conserved on superhorizon scales by causality [4,5]. As they re-enter the horizon, the sub-horizon physics, such as the baryon acoustic oscillations (BAO) or the structure formation process, further process these signals and the imprints of the heavy particles are recorded on the CMB and the LSS. The phenomena of inflationary particle production, along with the associated signatures in the power spectrum and non-Gaussianity of ζ (option (a) mentioned above), has been studied extensively in the literature. A relatively recent example of this is the "Cosmological Collider Physics" program [6,7] that aims to detect on-shell mass-spin signatures of heavy particles via nonanalytic momentum dependence and angular dependence of cosmological correlators. In particular, it has been shown recently [8][9][10][11][12] that particle masses m φ 0 , where φ 0 ≈ 60H * [2] is the slowly rolling inflaton velocity, can give rise to distinctive oscillatory signatures in primordial non-Gaussianity. For high-scale inflation, such a scale ∼ 60H * is interesting from a particle physics point of view since it can be close to the scale of Grand Unification. Motivated by this, in this work we ask what kind of mechanisms can produce even heavier particles with their typical masses m φ 0 , and if they can lead to novel phenomenology. A shift-symmetric, derivative coupling of the inflaton to the heavy particles involves factors of ∂ µ φ and therefore, φ 0 is the maximum energy budget that can be used to produce heavy particles and correspondingly, production for masses m φ 0 is extremely suppressed. However, if we allow shift-symmetry breaking interactions, this conclusion is avoided. A simple example of this involves a heavy particle whose mass is dependent on the inflaton field value. Then as the inflaton rolls along its potential, the heavy field mass may pass through a minimum, and depending on that efficient particle production can occur at that instant while it shuts off at other times. The above mechanism with shift-symmetry violation has been noted before, see e.g. [13][14][15][16][17] and [18] for a review. Here we incorporate two differences. First, the minimum of the heavy field mass is always H * , similar to [19] and therefore, the slow-roll inflaton trajectory need not be significantly perturbed. The number of produced heavy particles is also significantly smaller compared to scenarios where heavy particle masses pass through values H * . Second and more importantly, following their production, the heavy particles actually get significantly massive (and non-relativistic) such that they can directly modify the curvature perturbation by giving rise to a non-zero one-point function ζ( x) around their almost static locations. These localized perturbations, after reentering the horizon, give rise to spots/structures localized in position space. Given such a localized nature of the signal that is not repeated across various scales in the sky, it is natural to wonder how a direct position space search for such signatures compares to the traditional approach involving power spectrum, bispectrum or trispectrum of density perturbations.
To investigate this, in this work we present a concrete scenario in which such signatures show up as pairwise hotspots (PHS) on the CMB map. We focus on an example involving a heavy scalar field with an inflaton-dependent mass. We compute the number density of the heavy scalars following their production, and this determines the number of hotspots that appear on the CMB surface after an appropriate volume dilution. For each individual heavy scalar, we also compute the associated curvature perturbation ζ( x) during inflation and then propagate it forward through subhorizon evolution upon its horizon reentry. In particular, by taking into account the necessary Sachs-Wolfe and integrated Sachs-Wolfe effects, we obtain the final localized temperature profile corresponding to each heavy scalar. Owing to their non-adiabatic production from time-dependent vacuum, the heavy particles come in pairs and therefore we finally get pairwise signals on the CMB sky. As we will show later, the two hotspots of a given pair generally overlap with each other.
To make the comparison between position and momentum space searches mentioned above, we study both how the PHS can give rise to localized distortions in the CMB TT-spectrum and also localized signals on the position space. For the position space study, we simulate the effects of PHS using the HEALPix 1 package and perform a simplified "cut-and-count" analysis to identify the PHS from the background of standard CMB fluctuations. For the power spectrum based search, we numerically compute the CMB TT-spectrum in the presence of PHS using CLASS [20,21]. Our simplified analysis indeed demonstrates situations where a position space cut-and-count analysis performs better than the power spectrum based search. In particular, we show that the position space search may probe the minimum value of the scalar mass up to O(100)H * -comparable to the kinetic energy of the inflaton. At the same time, the search crucially relies on typical values of the heavy mass away from the minimum and that can be as big as O(10 4 )H * . In this sense, this search is sensitive to very high mass scales and can go beyond Grand Unification scales as well. While we do not consider the presence of instrumental noise and astrophysical foregrounds, we also do not use the detailed morphology of the signal in our search that is expected to mitigate the former two contributions. We will address this important aspect in a future work. We also do not consider in detail the resulting non-Gaussianity signatures for which a dedicated template study would be needed since the signatures show up on narrow windows in momentum space.
Some earlier literature have discussions on the production of heavy particles that relate to our study. The example model we use containing an inflaton-dependent scalar mass has similar behaviors to the examples studied in [19,22]. However, Ref. [19] and subsequently Ref. [23] focus on N -point correlation functions to study the signatures of periodic particle productions. In contrast, we focus on the position space signatures of heavy particles produced at a single instant of time during inflation. Therefore, in our set up the hotspots are characterized by a single scale and distortions in the momentum space ζ-correlators show up on a specific scale without any periodicity. Regarding the calculation of the temperature profile of the hotspot, Ref. [24] has a similar discussion although they assume the existence of pre-inflationary heavy particles that appear in our visible horizon. However, since inflation dilutes any pre-existing relic, one drawback of such an assumption is that the final number density of the hotspots strongly depends on the unknown amount of inflation we have had, strongly affecting the phenomenology. Our set up is free from this issue since particles are produced during inflation itself when the CMB-observable modes exit the horizon. Furthermore, the hotspots/defects they consider appear on much larger scales with different phenomenology. The idea of generating the PHS signals on the CMB map was discussed in Ref. [25] with the motivation of developing a quantum-entangled system for Bell's inequality test of the Universe's primordial fluctuations. Here in the context of a full model we also compute the number of such hotspots and study their properties. We also go beyond the computation of the primordial curvature perturbation and include the important effects of subhorizon transfer functions that are essential for the position space search.
The plan of the paper is as follows. In Sec. 2 we describe our model and review some necessary aspects of particle production. We then compute how many hotspots can be produced while ensuring small backreaction to the slow-roll inflaton dynamics. In Sec. 3 we discuss the detailed phenomenology of such hotspots by taking into account the subhorizon evolution. In particular, we derive both the central temperature of the hotspots as a function of the comoving horizon size at the time of particle production, and the spatial profiles of individual hotspots. These profiles are then injected into mock CMB maps generated via HealPix as described in Sec. 4. There we also outline our position space search strategy. After studying the effects of PHS on the CMB-TT power spectrum, we present our bounds on the mass-coupling parameter space. We conclude in Sec. 5. In Appendix A, we give some details of the "in-in" computation of the curvature perturbation due to a massive particle, while Appendix B contains some benchmark CMB images relevant for our study. trajectory, the heavy field mass is ∼ µ φ 0, the slowly rolling inflaton velocity. However, the mass passes through a minimum M0 µ during a single instant of (conformal) time η * and non-negligible particle production takes place for M0 φ 0.
2 Heavy particle production during inflation

Set-up
We denote the background inflationary spacetime metric by, where a(t) = e H * t is the scale factor given in terms of the inflationary Hubble scale H * which approximately remains constant during inflation. We will consider a scenario in which an otherwise ultra-heavy particle σ momentarily becomes lighter as the inflaton passes through some field value, and is significantly produced only around that time. To model such events of particle production where the mass of the heavy field passes through a minimum, we consider the following Lagrangian, This defines an inflaton-dependent effective mass for σ, that passes through its minimum at φ = µ/g. While in the context of more UV complete scenarios such as [19,26], and in scenarios with random mass functions [27,28], the full effective mass can have a more complicated dependence on φ, around the time of particle production, we can still approximate M 2 (φ) as a quadratic function of φ as in Eq. (2.3) quite generally. Using this parameterization, we will capture the physics of particle production through the parameters g and M 0 . One of the goals of this work is to use simulations of the CMB map to derive constraints on these two parameters. As we will see shortly, for the kind of signatures that we are looking for, we will typically need 1 g 10. While such values of g can potentially give rise to large radiative corrections to the inflationary potential, in scenarios with additional symmetries, such as supersymmetry, the radiative corrections can be in control, see e.g. [19].
Before getting into a detailed calculation, we briefly describe the qualitative feature of the above Lagrangian (2.2) whose relevant scales are shown in Fig. 1. We consider scenarios in which at a generic point on the inflationary trajectory, gφ ∼ µ M 0 H * and hence, the σ-mass is very large and cosmological production σ particles is exponentially suppressed and negligible. However, as the inflaton passes through the field value φ * = µ/g, M (φ) becomes significantly smaller momentarily, M (φ) ≈ M 0 . As a result, for M 0 ∼ φ 0 , where φ 0 ≈ 60H * is the slow-roll inflaton velocity determined by the scalar power spectrum [2], the kinetic energy of φ can be harnessed to produce σ-particles when φ ≈ φ * , corresponding to conformal time η * . Here and in the following, we use ' * ' to denote quantities evaluated at the time particle production. By construction, this mechanism of heavy particle production lasts for a short period of time and the imprints of the heavy particles would be most prominent on some particular k-modes that exit the horizon around η * . Therefore, the associated signature would show up as a violation of scale invariance. Our aim is to devise detection strategies to isolate such violations. While these effects would be observable in both the CMB and the LSS, in this work, we mostly focus on the effects on the CMB, leaving the effects on LSS for a future study.

Review of particle production
For the sake of completeness, we now review how to compute particle production given Eq. (2.2). Our approach will be somewhat pedagogical and readers can skip to the final result in Eq. (2.19) if needed. We start with the equation of motion (EOM) for σ, written in terms of the conformal time η ≡ dt/a(t). In the above and the following, s denote derivatives with respect to η. To further simplify the EOM into a harmonic oscillator form, we define u = σ/η, in terms of which we get, with a time-dependent frequency, (2.6) To estimate how large the effects of particle production and violation of adiabaticity is, we can compute ω (η)/ω 2 up to terms O(H 2 * /M 2 0 ) which we neglect. Let us now focus on our model where the time dependent mass is given in Eq. (2.3). Denoting µ ≡ gφ * and using the fact that during inflation the homogeneous solution under the slow-roll approximation goes as, we can rewrite the mass term as (η ∝ −e −Ht ), This gives a measure of adiabaticity violation, where we have assumed M 0 gφ| ln (η/η * )|/H * to get the last relation. We have also set k = 0 to estimate the maximum amount of adiabaticity violation. This then shows that particle production is most efficient just around η * where the mass is minimized. The same conclusion can be obtained via the steepest descent method along the lines of [29].
Given this we will study the simplified equation of motion around η * to study particle production. Therefore, we approximate (2.11) Using this Eq. (2.5) can be written around η = η * as, where (2.13) To estimate particle production, we will solve Eq. (2.12). We follow the standard procedure (see e.g., [3]) to estimate the Bogoliubov-β coefficient by studying the amplitude of a negative frequency mode that develops after starting with a positive frequency solution at very early times. To this end, let us first obtain the adiabatic solutions for small and large times. For τ → ∞ the ±−frequency solutions are, (2.14) For τ → −∞ on the other hand, ±−frequency solutions read as, The exact solutions to Eq. (2.12) can be written as parabolic cylinder functions W − κ 2 2 , ± √ 2τ [22,30]. To obtain β we will look at the asymptotic behavior of W . It can be checked that the following linear combination gives the positive frequency solution at τ → −∞, where σ = √ 1 + e −πκ 2 − e −πκ 2 /2 . The Bogoliubov-β coefficient can be read off from the negative frequency part of this solution as τ → +∞ which is given by, Thus we have the standard Bogoliubov α and β coefficients as, It can be checked that they satisfy, |α| 2 − |β| 2 = 1 as required. Thus for a given k-mode we get [19,22], where we have approximated, φ ≈ −φ H * η * . This implies that for any comoving momentum k, particle production is Boltzmann suppressed if the minimum value of the time-dependent mass of the particle, here M 0 , is bigger than gφ. This can be understood as coming from the adiabaticity violation factor e − ω 2 ω . Focusing on k = 0 for simplicity, ω 2 = (gφ − µ) 2 + M 2 0 is the time-dependent frequency. The quantityω/ω 2 parameterizes the violation of adiabaticity andω/ω 2 is maximized when gφ − µ = M 0 / √ 2 M 0 . This implies at the time of particle production, ω M 0 andω gφ which gives, |β| 2 ∼ e − ω 2 ω e −M 2 0 /(gφ) . We also mention that Eq. (2.19) is only valid for |β| 2 < 1 and cannot be applied to scenarios where M 2 0 < 2H 2 , which lie outside our parameter space of interest.

Estimating particle production
We can now calculate the physical number density of produced particles at η * after a phase space integral over physical momenta k phys = k|H * η * |, (2.20) As we will show in the following, each of the particles leads to a localized "disturbance" in the spacetime metric that can manifest itself as either a hot or a cold spot on the CMB. The hot/cold nature depends on the size of the disturbance and is primarily controlled by η * . As we will detail in Sec. 3, for the benchmark choices of η * we will focus on, the disturbances will be hotspots on the CMB and therefore we will use this terminology from now on.
Notation. Except for places where |η * | appears, we will now absorb the minus sign in η * and use it to denote the size of the comoving horizon at the time of particle production, instead of the (negative) conformal time.
The number of such hotspots within the CMB surface, about η 0 = 13.8 Gyr [31] (conformal age of the Universe) distance away with a thickness ∆η rec , is given after an appropriate volume dilution factor, where a * is the scale factor at the time of particle production. In the last relation, we have used the fact that a mode with comoving momentum k crosses the horizon at time t when k = a(t)H * (t). Thus, k * = a * H * is the mode that exits the horizon at the time of particle production. The ratio involving the thickness of the surface of recombination is given by, ∆η rec /η 0 ≈ 10 −3 with ∆η rec ≈ 19 Mpc [32]. Therefore, the number of hotspots inside the surface of last scattering depends on our model parameters g, M 0 and η * as, (2.23) Here we have used the horizon crossing relation k * η * = 1. The exponential sensitivity to M 0 comes from the non-adiabaticity factor e − ω 2 ω mentioned above, but we see that even for M 0 H * , there can still be an appreciable number of hotspot production. The scaling ∝ 1/η 3 * denotes the volume dilution of hotspots after their production. Finally, g controls the number of hotspots and, as we will see below, the temperature of individual hotspots as well.

Backreaction constraints
Before exploring PHS phenomenology and detectability in more detail, we review the existing constraints on the setup. First we want to ensure that the production of heavy particles does not affect the background inflating spacetime [19,29]. The EOM obeyed by the inflaton reads as, In the presence of the interaction implied by the φ−dependent mass term in Eq. (2.2) we get where V φ is the background inflationary potential that drives the inflationary expansion and σ 2 is the vacuum expectation value with respect to the early time vacuum. To ensure the usual slow-roll dynamics for the inflaton, 3Hφ ≈ − Using Eq. (2.9) we can rewrite the above condition as, where we have used the fact that n ∼ M σ 2 is the number density of the produced σ particles around η * , given by eq. (2.20).
To ensure that particle production does not deplete the kinetic energy of inflaton, we require the energy density of the produced σ particles to be smaller thanφ 2 , i.e., Using Eq. (2.9) once again, we can rewrite the above constraint as, which for ln(η/η * ) ∼ 1 reduces to constraint (2.28). The constraint in relation (2.28) is shown in Fig. 2 where the left hand side is a percent of the right hand side and above the green line labeled "% Backreaction" the backreaction will be more than a percent. Various contours of N HS are also shown which determine how many hotspots can be produced with negligible backreaction on the inflationary dynamics. The constraints should be thought of as conservative, as the left hand sides of Eqs. (2.28), (2.30) involve the number density at |η * |, the maximal value, and have no dilution effects. As presented, the constraints are primarily on g and M 0 /H * , with only a very mild dependence on η * .

Hotspot phenomenology
Having estimated the strength of particle production, now we focus on the effect of the produced particles on curvature perturbation and, eventually, on the CMB. We do this in two steps. First, we compute the curvature perturbation due to a massive particle. Then, we include the effects of sub-horizon transfer function to see the final form of the temperature anisotropy.

Curvature perturbation due to a massive particle
To compute the curvature perturbations due to σ production, we employ a sudden transition approximation, i.e., we assume that the particles become non-relativistic immediately after their production. This is justified through Eq. (2.19), which shows that particle production is exponentially suppressed unless physical momenta of the particles kη * H * M 0 , for M 2 0 gφ (as we will choose later). Furthermore, after their production the particle masses increase, whereas the physical momenta redshift further. Therefore, the action due to a single heavy particle can be written in a non-relativistic form as, where x µ = (t, 0) is the location of the static massive particle, chosen to be at the origin. To see how this contributes to the curvature perturbation, we parametrize the metric fluctuations as [33], with the condition ∂ i γ ij = 0 and γ ii = 0. In this gauge there is no inflaton fluctuation, δφ = 0. The quantities N and N i have vanishing time derivative and therefore are not dynamical. In particular, N obeys the algebraic relation [34], In the last line we have expressed N in terms of the comoving curvature perturbation ζ [35], as appropriate in the δφ = 0 gauge. We can now rewrite the action integral (3.1) into (3.5) The last integral determines the effect of a heavy particle on the curvature perturbation. In terms of the conformal time η, the ζ-tadpole term is The above interaction term is non-trivial due to the time dependence of the mass M (η), and it gives rise to a non-zero one point function of ζ, ζ HS . We will see below that after including the effects of transfer function, ζ HS above would give rise to hotspots (HS) on the CMB for our parameter choices, as opposed to cold spots. This one-point function ζ HS can be calculated using the "in-in" formalism along the lines of [25] (for a derivation see Appendix A), where = −Ḣ * /H 2 * . We can understand the above result intuitively as follows. Due to causality, the massive particle can not change the curvature perturbation on scales larger than the horizon size at the time of particle production, determined by η * . Hence there is no effect for |η| > |η * |. 2 After production, the particle mass keeps changing gradually due to the slow-roll of the inflaton field and this temporal change is recorded through a spatial variation of curvature. For example, a spherically symmetric shell that exits the horizon at |η 1 | < |η * | is perturbed according to M (η 1 ). Another shell that exits at a later time |η 2 | < |η 1 | is perturbed according to M (η 2 ). Thus a hotspot shows up with a comoving size |η * |, and the radial curvature profile inside that region is determined by the time evolution M (η) during inflation and after η * .
For our model, M (η) varies as in Eq. (2.9). Therefore, for r η * the hotspot profile can be rewritten after dropping the subleading contributions due to M 0 and usingφ 0 = Given the size of the standard inflationary quantum fluctuation at any point, ζ 2 , the visibility of the hotspot over such fluctuations is mostly controlled by the coupling g. The above relation also indicates that for visibility of the hotspot signal, we typically need g 1, as will also be seen below through our numerical simulations. Given the above relation, we see that while M 0 determines the particle production rate, the hotspot perturbation is determined by the time-dependent mass away from the minimum. This is so because already after one e-folding, M (η) gets dominated by the first term in Eq. (2.9). Therefore, the hotspot profile for r < η * /few is determined by g.

Temperature anisotropy from curvature perturbation
Having calculated the curvature perturbation due to the massive particle, we now connect it to the observed temperature fluctuation on the sky. We write the metric fluctuations in the Newtonian gauge as, (3.10) We can write the temperature fluctuation on the sky corresponding to a Fourier mode k as, δT ( k,n, η 0 )/T ≡ θ( k,n, η 0 ), wheren denotes the direction on the sky, η 0 is the present (conformal) age and T is the average CMB temperature. Then we can rewrite this in terms of Legendre polynomials as, where the moments θ l (k, η 0 ) get contribution from Sachs-Wolfe (SW), integrated Sachs-Wolfe (ISW) and Doppler effects. Given the spherically symmetric profile of the hotspot, the Doppler contribution to θ l (k, η 0 ) is very small, as can be checked explicitly. Therefore, we focus on the SW and the ISW contributions which are given by, where j l (x) is the spherical Bessel function and η rec denotes the time of recombination. We have also isolated the two functions f SW (k) and f ISW (k) to characterize the SW and ISW contributions after mode re-entry with ζ( k) denoting the primordial perturbation as before. In practice, we extract f SW (k) by running CLASS for a purely ΛCDM setup. We numerically fold in the spherical Bessel function with the CLASS output of θ 0 (k, η), Ψ(k, η) perturbations from which the primordial fluctuation ζ( k) is already factored out. This leaves us with f SW (k). Similarly, numerically integrating Ψ (k, η), Φ (k, η) with the spherical Bessel function and the optical depth (second expression in Eq. (3.12)) gets us f ISW (k). The f SW (k), f ISW (k) extracted in this fashion can be applied to ζ HS (k) as subhorizon physics is independent of the origin of the initial curvature perturbation.
With f SW (k) and f ISW (k) in hand, we are now ready to compute the temperature perturbation due to the heavy particle. We rewrite the momentum space expression of ζ HS derived in Eq. (A.7) as, where the function f (x) = gH 2 φ0 (Si(x) − sin(x)). Thus we can write the temperature fluctuation at our location x 0 as, (3.14) Using the plane wave expansion, the distance to the hotspots, 16) and the identity between Legendre polynomials and spherical harmonics Y lm , Eq. (3.18) gives the temperature profile around the hotspot as a function of the angle between the direction of observationn and the direction of the hotspotn HS . We will use this profile in the subsequent analysis of hotspots in CMB simulations.
To obtain the central temperature of a hotspot we can setn ≈n HS implying P l (n ·n HS ) ≈ 1. For small enough η * , the ISW contribution can be neglected. In that case, a simple expression of the central temperature can be obtained by noting f SW (k) = T SW (k)j l (kη 0 − kη rec ) with T SW (k) being the SW transfer function (after factoring out ζ( k)). Upon further using the approximation For larger η * , we need to keep the ISW contribution as well. The full result for the central temperature for various values of η * is shown in Fig. 3. We see the net effect of the subhorizon evolution is to give rise to a hotspot, as opposed to a cold spot. In Fig. 4 we will show the angular profile for a given hotspot as a function of the angular distance from the center of the hotspot.

Separation between hotspots
In the above, we have seen how a single heavy particle gives rise to a hotspot in curvature perturbation. However, due to momentum conservation, heavy particles are produced in pairs. Hence as long as the heavy particles are sufficiently separated, what we expect to see is not a single hotspot but rather a pair of them.
To get an idea about the separation between the two members of a hotspot pair, we recall from Eq. (2.19) that the particles are produced either semi-or non-relativistically. Furthermore, within one Hubble time after particle production, the effective mass of a particle increases from ≈ M 0 to gφ H * | ln(η/η * )|, as seen from Eq. (2.9). Hence, most of the particle propagation due to non-zero physical momentum k phys = ka −1 ∼ (gφ) 1/2 happens immediately after production. The resulting comoving distance travelled by a massive particle up to conformal time η f can be estimated as, To get the above we have used M ≈ M 0 at production and k phys ≈ kη * H * M 0 . Thus we get the intuitive result that for semi-relativistic production, the particles can travel a horizon-sized distance in one Hubble time, whereas for non-relativistic production, the particles mostly remain where they were produced. Given this, we will model the separation between the particles of the same pair to be a random uniform distribution between 0 and η * .

Effect of localized hotspots on the CMB power spectrum
Before discussing in detail how to search for the above-mentioned hotspots, we can ask whether the CMB power spectrum is already good enough for that purpose. To get an analytical understanding, we consider a toy model in which a hotspot is represented by a delta function on the CMB sky. Therefore, the CMB temperature fluctuation reads as (3.21) Here, ∆T inf denotes the standard inflaton-sourced CMB fluctuation and (θ * , φ * ) is the location of the hotspot with a temperature T h (∝ g). By doing the standard spherical harmonics decomposition we get, Therefore, the temperature power spectrum estimator C TT reads as, where the interference term vanishes between the two uncorrelated perturbations. Conventionally, one looks at D ≡ ( + 1)C /(2π) rather than C TT for the power spectrum: Including the g dependence in T h and multiplying by the number of hotspots yields a PHS contribution to D TT ∼ g 2 N HS 2 that increases at higher -modes. However, a more realistic model of the hotspot must include their finite profile size. For hotspots with size ∼ η * , -modes with √ 4πη −1 * correspond to epochs after the particles were produced and should make negligible contribution to D TT . Production of modes with k 1/η * is also Boltzmann suppressed as seen from Eq. (2.19). Combining the growth at low from the toy model with vanishing effect from large , we expect the PHS impact on D TT to be most prominent at * . Depending * , and g, such an effect could be hidden within the uncertainties on the power spectrum. So, while the power spectrum will constrain the PHS model to some degree, it is worthwhile to consider other detection possibilities as well.

Detection strategy
Here we discuss the possibility of identifying PHS using a simplified search for the hotspot signals in position space. We begin with some comments about the various backgrounds to the PHS signal and then motivate two different search strategies.
There are three types of backgrounds to consider when searching for PHS: i.) detector noise, ii.) the astrophysical foreground, and iii.) background from the standard primordial fluctuations. Detector noise depends on the details of the instrument and can be reduced further at future CMB experiments. The foreground emissions arising from compact objects such as galaxies, galaxy clusters, gas, and dust that exist between the last scattering surface and us, can also produce localized signals. Since the compact sources usually produce higher frequency signals than the CMB photons, the comparison between different frequency maps has been a powerful tool to subtract the foreground in CMB studies [36]. Therefore, it is possible to remove the foregrounds by comparing the sky maps taken at different frequency bands. Moreover, the shape of the compact objects are usually not spherically symmetric and may be distinguished from the signal hotspots that we are looking for. Therefore, in the following, we only consider the background from the primordial, almost Gaussian fluctuations when studying the PHS signal. We will assume the CMB maps are masked to reduce the astrophysical foregrounds and badly-conditioned pixels and retain only 60% of the sky for the analysis. The number is similar to the sky fraction used in the Planck analysis [31].
We consider two types of signature from the PHS.

1.) Position-space search:
The PHS show up as localized objects in the sky, motivating a search in position space. Even in idealized scenarios, the primordial quantum fluctuation can still generate a background for the PHS and, unlike the astrophysical foregrounds, background signals from the primordial fluctuation follow the same blackbody distribution as the PHS signal. We need to rely on the hotspots' detailed properties, such as their temperature and distribution, to distinguish them from the background. We therefore consider a "cut-andcount" analysis in position space by applying temperature cuts on the entire sky to veto the standard CMB fluctuations. While past position space searches have utilized wavelet analysis based approach [37][38][39] or a more dedicated template fit (e.g., [40]), in this work we explore how much an even simpler "cut-and-count" method can be useful.
2.) Power spectrum analysis: As discussed in Sec. 3.4, the presence of PHS also modifies correlation functions of CMB fluctuations. The existing measurements of the CMB power spectrum, the searches for bispectrum and trispectrum can constrain the number of PHS.
Our primary concern will be the corrections to the power spectrum C TT , though we will also briefly mention the rough (scale-dependent) amplitude of the non-Gaussianity parameter f NL for the same number of PHS under consideration. There will also be corrections to power spectra C TE,EE involving E-mode polarization. Since such corrections are of similar orderof-magnitude compared to C TT , to get a first bound on our parameter space, we will not consider polarization spectra in this work.
Before discussing these two search strategies in detail, we step through the numerical simulation recipe we follow to generate the PHS signal and mock CMB maps.

Simulating the PHS signal
In this section, we describe a detailed methodology for simulating CMB maps and PHS signals. The first step is to generate a temperature power spectrum C TT using the Boltzmann code CLASS v2.9 [20,21] with the best fit ΛCDM parameters The next step is to convert this angular power spectrum into sample CMB maps using HEALPix [41].
A key input for the HEALPix simulation is N side , which controls the resolution of the CMB maps: N pixels = 12 N 2 side . The number of pixels can be compared to the number of independently fluctuating temperature patches max =1 (2 + 1). Thus a simulation using N side = 256 corresponds to an max 880 while N side = 2048 corresponds to max 7000. 4 The price for larger N side is significantly slower simulation speed. Balancing resolution and computational cost, we proceed with a hybrid approach. Specifically, for the cut-and-count analysis, which requires a large number of simulated CMB images, we use lower resolution, typically N side = 256 or 512. However, for the power spectrum analysis we do not need to generate multiple images, but we want to capture the impact of the PHS up to ∼ 2000 (Planck 2018 range) -so we use N side = 2048. A few other useful relations to keep in mind are i.) how to convert between mode and angle on the sky: θ = √ 4π −1 max , and ii.) how to convert between η and -mode : * η 0 /η * , with η 0 ≈ 13.8 Gpc being the current conformal age of the Universe. For a given HEALPix CMB map, we next add a population of PHS signals based on the profiles that we derived above in Eq. (3.18). From Eqs. (3.13), (3.20), we see that the angular size, separation and temperature of the PHS signal are controlled by the heavy scalar-inflaton coupling g and comoving horizon size η * at particle production, while the overall number of PHS also depends on M 0 . For our simulations, we pick some benchmark values of η * and g, then use the signal and background analysis to set limits on M 0 . The benchmark parameter points we choose are η * = 160, 100, 50 Mpc , g = 2, 3, 6, 10. (4.2) These η * values roughly correspond to * ≈ 85, 140, 280. The shape of the temperature profiles (i.e., setting g = 1 in Eq. (3.18)) is controlled by η * alone and is shown in Fig. 4 as a function of the angle cos θ =n ·n HS away from the center of the hotspot. The vertical lines in the figure show the angles √ 4π/ * that correspond to the angular sizes of η * in the sky. We pick these benchmark values to show examples with different corrections due to the transfer function. If η * is much smaller than 50 Mpc, suppression from the transfer function will diminish the signal to be below the primordial fluctuation even for a large coupling g. We start the maximum horizon from 160 Mpc so that we can compare the result of the position space search to the C TT analysis with a reasonably small uncertainty from cosmic variance.
For the η * = 160 Mpc benchmark, we use N side = 256 maps, corresponding to max 880 or angular resolution of ≈ √ 4π/ max = 0.004. The δT profile in this case effectively extends to θ ≈ 0.02, as is shown in Fig. 4, therefore we divide the hotspot radius into 5.5 pixels in the HEALPix image, where the half pixel corresponds to the central pixel of the hotspot, so that in the 2D  pixel-space, an individual hotspot spans a 11 × 11 block. The temperature in each hotspot pixel is set to the average value of δT across that pixel, following profiles in Fig. 4 and scaled by the approrpriate benchmark value of g. After adding one hotspot for the PHS signal, we add its partner -a second, identical spot, but with its center offset from the original spot. The offset is determined via a uniform random distribution up to a maximum angle of ≈ √ 4π/ * (shown as vertical lines in Fig. 4). For η * = 160 Mpc, the maximum separation corresponds to 10 pixels (∆θ ≈ 0.04).
The PHS signal generation procedure for the η * = 100 Mpc and η * = 50 Mpc benchmarks is identical. For η * = 100 Mpc, the δT profile effectively extends to θ ≈ 0.02. We again use N side = 256, so the hotspots have a radius of 5.5 pixels and a maximum separation of 6 pixels. For η * = 50 Mpc, the δT profile effectively extends to θ ≈ 0.015. The spots are smaller so we increase the resolution to N side = 512, leading to a spot radius of 8.5 pixels and a maximum separation of 6 pixels. Note that, with these sizes and separations, the profiles of the two hotspots from a given pair will overlap for all benchmark points. Also, as signals for the η * = 100 and 50 Mpc cases are colder and more difficult to search for, we consider an even larger coupling g = 10 in these two cases when presenting the results. Some example images from the PHS signal and combination of signal and CMB background for the η * = 160 Mpc with g = 6 benchmark are shown below in Fig. 5. The images focus on a small patch in the sky within ≈ [−16 • , 16 • ] for both the longitude and the latitude. The upper left and upper right plots are the PHS+CMB and PHS-only plots made with N side = 256. To illustrate how N side changes the resolution, the lower left plot is the simulated CMB map, and the lower right plot is the same signal model but made with N side = 2048. Example images of the η * = 50, 100 Mpc scenarios are shown in Appendix B (Fig. 10 and 11).

Searching for PHS signals with temperature cuts
From the signal panels in Fig. 5, the PHS signal looks very similar to jets at the Large Hadron Collider -localized regions on a 2D plane with ∼ circular shapes that stand out as hotter than the background (more energy/transverse momentum in the case of jets). That analogy motivates us to try a position-space search for the PHS using cuts on temperature, inspired by how jets are hunted for at colliders.
Rather than searching the entire sky all at once for PHS, we divide the CMB maps into patches, then seed them with PHS signals. As part of this process, the patches are projected from spherical coordinate maps to Cartesian coordinates. This procedure -cutting the sky into 'flattened' patches -is not technically necessary but is often adopted in the literature [42][43][44][45][46]. Historically it was used since CMB analyzers using wavelets or machine learning could not cope with spherical coordinates [37][38][39][47][48][49]. This technicality has been overcome in the neutral network analysis [50,51], but we will stick to the flat patches for simplicity in this analysis. Projecting from spherical to Cartesian coordinates does introduce some degree of distortion. To mitigate this effect, we restrict the simulated CMB patches to a region around the equator made by 140 × 140 pixels. Taking a sample 140 × 140 pixel patch, we implement the following cuts to isolate the PHS signal from the CMB background: • Reduce δT CMB contribution from lower -modes by cutting the image into regions with a certain < ∼ * , calculate the average temperature of each region and subtract it from the image.
• Further veto pixels below a temperature cut δT cut that is chosen based on the PHS temperature.
Even though the temperature of PHS we consider in Eq. (4.2) is typically much larger than the primordial δT ≈ 27 µK of a single mode before entering the horizon, the total standard deviation -formed from summing over a range of -is much larger. Specifically, where min is set by the size of the patch and max is set by the pixel resolution, roughly max ≈ 12N 2 side . For the range of for our benchmark signals, σ CMB ∼ O(110 µK), comparable to the PHS signal temperatures ∼ g 2 δT with g = O(1). To reduce background coming from modes with much smaller than the hotspot signals ( * ), we divide the simulated maps into regions that correspond to modes with < ∼ * , calculate the average temperature of each division and subtract it from the original map. As is shown in the rightmost panel of Fig. 6, this process reduces the standard deviation of the δT CMB distribution and helps to suppress the background. Depending on the benchmark this subtraction process leads to O(10 − 30)% improvement in the final constraints on the allowed number of hotspots. As a specific example, for our N side = 256 simulations we divide the 140 × 140 pixel map into 14 × 14 pixel blocks. Subtracting the average temperature in these sub-patches raises min from ≈ 6 to ≈ 60, and lowers σ CMB to 91 µK.
After subtracting contributions from the lower-modes, we apply a δT cut to image in the 140 × 140 pixel region. The T cut is chosen to optimize the signal significance -to be defined in more detail below -but ranges from 140 µK to 340 µK depending on the benchmark point. For such a simple cut, the number of background spots N CMB in the region can be estimated analytically by a Gaussian distribution of temperature fluctuations. For example, in the η * = 160 Mpc analysis, the standard deviation of the δT fluctuation is ≈ 91 µK after subtracting the low-modes. Once applying δT cut ≈ 160 µK to the simulated image, a Gaussian temperature fluctuation only gives ≈ 4% chance for each single pixel to pass the cut. This corresponds to N CMB N pat ≈ 3 × 10 4 for the entire sky, where N pat = 40 (160) is the number of 140 × 140 pixel regions in the whole sky for the η * = 100, 160 Mpc (50 Mpc) analysis. The number is close to the result derived from the numerical analysis (and shown in Table 1).
To determine the significance of a given PHS scenario, we determine: 4) where N n sig+CMB and N CMB are the average number of hotspots identified in the 140×140 pixel maps with n-signal injections and no injections. f sky = 60% is the sky fraction that we assume for the analysis. We determine N n sig+CMB by generating 10 3 CMB background maps and increasing the number of injected PHS signals until Eq. (4.4) equals 1 or 2 (for 1σ or 2σ, respectively) for a given T cut . We then iterate the entire process, varying T cut until we find the optimal value. The optimized temperature cuts for each benchmark PHS point are shown in Table 1. We also display the number of PHS pairs that must be injected for each benchmark in order to achieve Sig = 1 ( N 1σ ) or Sig = 2 ( N 2σ ). In other words, N 1,2σ is the maximum number of PHS pairs in the entire sky that can hide (at 1 or 2 σ) from our search. To get an idea of the signal efficiency under the temperature cut, one can compare 2 N 1,2σ to N n sig+CMB − N CMB . Depending on the benchmark point, we find the efficiency (per PHS pair) ranges from O(5%) (η * = 50 Mpc, g = 3) to O(70%) (η * = 100 Mpc, g = 10), see Table 1. From Table 1, this simple cut and count search is sensitive to O(10 − 10 4 ) PHS, a number that the scalar model discussed in Sec. 2.4 can easily generate.
One signal feature that is obviously missing from our analysis is the pairwise nature of the hotspots. Viewing the signal on its own, as in Fig. 5, the fact that there are always two spots in close proximity is striking, and it seems like one should be able to exploit this to increase the power of a position-space search. Taking a cue from jet physics once again, one simple way to capture this pairwise feature is impose an additional cut requiring any pixel hotter than δT cut to have a second hot pixel (either at the same temperature or slightly different) within some neighborhood in order to be considered a viable signal candidate. In practice, we find that a search looking for pairs of hot pixels within disc-like regions on the sky fares worse than the simple δT cut search, even when optimized over the distance between hot pixels and the individual pixel temperature cuts. Part of the difficulty with the paired pixel search is that the core temperature of the PHS -found by multiplying the profile in Fig. 4 by g -is not significantly bigger than σ CMB . For example, for g = 3, the core temperature of the η * = 160 Mpc hotspots is O(120 µK), while for the η * = 50, 100 Mpc benchmarks the core temperature is only O(30 µK). With PHS temperatures so close to σ CMB , it is impossible to find δT cut values for two pixels that are both efficient for the signal and not plagued by high background fake rate. Dropping the requirement for two hot pixels helps because we can raise the single pixel δT cut ; while the individual PHS may not be particularly hot compared to σ CMB , the profiles overlap, so the odds of one pixel passing an increased δT cut remain high while the background rate falls precipitously. A more sophisticated analysis -perhaps using template fitting or a machine learning to better utilize the shape of the profiles -could enhance the signal sensitivity significantly, and we will pursue this in a future study.

PHS signatures in the CMB power spectrum and bispectrum
As we discussed in Sec. 3.4, the existence of PHS gives an overall increase of the C TT spectrum for modes < ∼ * ≈ η 0 /η * , and the increase is proportional to the number of hotspots N HS and the heavy scalar-φ coupling g 2 . To verify this behavior, we add N HS hotspots to the HEALPix CMB map using the profile from Eq. (3.18) (translated to pixels), then derive the power spectrum. In order to be sensitive to = 2000, we work with a higher resolution, N side = 2048. Some example higher resolution images are shown in Figs. 5, 10, 11 (bottom right panels).
In Fig. 7, we show the D TT spectra obtained from the HEALPix for the ΛCDM model (blue) with the best-fit values of the parameters from Planck 2018 results [31], and also from the ΛCDM+PHS (red) with the number of PHS that gives a 1σ excess in the cut-and-count analysis (Table 1). In the lower panel of each plot, we show the difference between the two scenarios (green) as well as the 1σ error bars from the Planck 2018 result [31] (gray), following the same binning of . We also give the total ∆χ 2 from the green curves in each plots. Our reasoning for showing the deviation from a 1σ excess (instead of 2σ) of the cut-and-count search is that we would like to compare ∆χ 2 η * = 160 Mpc   Table 2. Minimum numbers of signals N ∆χ 2 =6 ( N ∆χ 2 =20 ) in the sky for a ∆χ 2 = 6 (∆χ 2 = 20) excess based on the goodness-of-fit test using C TT distributions.
value to the numbers from the Planck fit which roughly correspond to 1σ variations of the ΛCDM parameters. Studying Fig. 7, we see the green curves roughly follow the behavior expected from Sec. 3.4, though with more complicated dependence coming from the detailed shape of the profile. For g = 6, we see < 1σ deviations for all the bins with the η * = 160 and 100 Mpc benchmarks. These deviations of the C TT spectra contribute to total ∆χ 2 ≈ 1 − 5. As a comparison, the ∆χ 2 between the fits of Planck 2018 TTTEEE+lowl+lowE(+lensing) data with the best-fit ΛCDM parameters and the average parameters is ≈ 26 (6) and is ≈ 20 if including only TT+lowl+lowE data [52]. Given that the best fit parameters in these data sets are within the 1σ variation around the average value of each parameter, the ∆χ 2 of the deviation in D TT from PHS signals with a 1σ excess in the cut-and-count analysis is smaller compared to the ∆χ 2 between ΛCDM best fits and average values  Figure 7. D TT spectrum of the CMB using best fit ΛCDM input parameters in Eq. (4.1) with (red lines) and without (blue lines) PHS signals injected in the sky using a resolution parameter N side = 2048. The differences between the two distributions are shown in green lines, and the shaded regions correspond to 1σ uncertainty, taken from the Planck 2018 data [31]. Here we show examples with g = 6. Several benchmarks of η * and the number of PHS signals corresponding to N1σ bounds in Table 1 are chosen. The corresponding ∆χ 2 excess based on the goodness-of-fit test is shown in each panel.
of Planck data. Therefore the cut-and-count analysis would be able to place a stronger constraint compared to standard TT power spectrum search.
However for η * = 50 Mpc benchmark, the corresponding ∆χ 2 = 38, and the same number of PHS from the cut-and-count study may have been excluded by the power spectrum measurement. To be conclusive, a more dedicated study is necessary to constrain the PHS signal using the Markov Chain Monte Carlo (MCMC) method that would fit the observed data using both the ΛCDM and PHS parameters. In the present work we do not carry out such an MCMC study, but to get a feeling of the relative sensitivities between the cut-and-count search and the C TT measurement, we calculate the number of PHS signals that corresponds to either ∆χ 2 = 6 or ∆χ 2 = 20 in each model and show the results in Table 2. The required N 1σ for the cut-and-count search with g = 2 ∼ 6 are all lower than N ∆χ 2 =6 for the deviation of the TT spectrum. This means the simple cut-and-count study can already provide better constraints to the power spectrum measurement for η * =160 and 100 Mpc. A search for bump features with a specific profile using Planck 2015 data and forecasts for LSS surveys appear in [53]. It would be very interesting to carry out a similar procedure for the present scenario. The above results from our search for hotspots is summarized in Fig. 8.
In addition to correction to the power spectrum, the presence of PHS also generates primordial  Table 1 (Table 2). This shows that both for η * =160 Mpc and η * =100 Mpc, the simple cut-and-count search can perform better than a power spectrum search.
non-Gaussianity encoded in the bispectrum and higher point correlation functions, for k * ∼ 1/η * . From a naive dimensional analysis, the going rate of the magnitude of the dimensionless bispectrum, for a given number density n of PHS per Hubble volume and coupling g, is of order f NL (k * ) ∼ 10 3 g 3 (n/H 3 * ) [19]. For g = 2 and η * = 160 Mpc, we have (n/H 3 * ) ≈ 0.05, corresponding to the 1σ result in Table. 1. This gives f NL (k * ) ∼ 400. Such a large f NL (k * ), however, only shows up in a narrow range of k-modes around k * ∼ 1/η * . As was found in [19], for similar models of particle production, the overlap between a large three point functions with a narrow k-range and the standard equilateral template but with a wider k-range can be small, and new templates are necessary to pick up the non-Gaussian signal. A careful study involving N -point correlation functions (allowing for N > 4) of the model in [19] was performed in [23] using the WMAP results. To make a rough comparison, we can consider small values of ω (frequency of particle production) in Fig. 7 of [23], as appropriate for our scenario of particle production at a single instant and use the relation, to relate the quantity n/H 3 * used by [23] and N PHS in Fig. 8. Using the upper bounds on N PHS in Fig. 8, we then see the constraints derived using our simplified search is better than the sensitivity quoted in [23]. We need to improve both the positions space search and the apply N -point function   Table 1 (Table 2). We compare the 2σ bounds for different production times η * = 160, 100, 50 Mpc in the lower-right panel. The results for η * = 50 Mpc and η * = 100 Mpc cases are almost identical and the corresponding bounds overlap.
searches to the present context to better understand the complementarity between the two types of studies.

Bounds on the heavy scalar production
Using Eq. (2.23), we can translate the N 1,2σ values in Table 1 to the minimum mass of heavy scalar M 0 at the time of particle production. We show the results as the red (1σ) and orange (2σ) curves in Fig. 9. We join the points of different g values to illustrate the trend of the M 0 /H * bound as a function of g. Parameter points that lie below the curves are ruled out. The simple cut-and-count analysis can be sensitive to M 0 values that are O(100) times larger than the Hubble H * during the inflation. At the same time, the central temperature of the hotspots are determined by a typical value of the heavy mass away from the minimum, and is ∼ 10 4 H * . This shows the possibility of looking for extremely heavy particles produced by the inflationary dynamics using our position space search. In the plots we also show, via gray dashed and dotted lines, the bounds that correspond to C TT deviations as given in Table 2, which give us an idea about how power spectrum searches can constrain the heavy scalar masses.
We see both for η * = 100 and 160 Mpc, the position space search performs better. We also show in light blue various contours labelling how many hotspots are produced for different values of M 0 and g. These curves illustrate the exponential sensitivity in the number of hotspots to M 0 /H * , so that slight differences in the bounds -either comparing the position space search to the power spectrum or jumping 1σ ↔ 2σ -correspond to big differences (∼ order of magnitude in some places) in the number of hotspots. For all the bounds, we stay above the "% Backreaction" contour, implying the backreaction on slow-roll dynamics due to particle production is sub-percent. Finally, in the bottom right panel of Fig. 9, we combine different η * benchmarks to see that the bounds do not vary drastically as we change η * .

Conclusion
In this work we have studied how production of very heavy particles during inflation can lead to localized perturbations in the spacetime curvature, and eventually lead to anomalously hot or cold spots on the CMB sky. The key ingredient that leads to heavy particle production is an inflaton-dependent mass of the heavy field σ which becomes small (but still H * ) at a particular (conformal) time η * , i.e., M 2 (η) = M 2 0 + f (φ) where φ is the inflaton, and M 2 (η * ) = M 2 0 f (φ) for typical values of φ. Quite generally when the heavy mass pass through such a minimum, we can approximate the mass-squared function as quadratic in φ, M 2 (η) ≈ g 2 (φ − φ * ) 2 + M 2 0 . This implies the σ − φ coupling is controlled by a single parameter g. At η * , deviations from adiabaticity are O(H * /M 0 ), and following the steps laid out in Ref. [22], we find particle production can be significant even for M 0 ∼ O(100H * ). We also compute the number of hotspots that appear within the CMB surface. The non-adiabatic nature of particle production from a time-dependent vacuum also implies that such spots come in pairs with the mutual separation dependent on the size of the comoving horizon η * .
We then discuss in detail how the primordially produced heavy particles lead to temperature perturbations in the CMB, being careful to include subhorizon evolution dominated by Sachs-Wolfe and integrated Sachs-Wolfe effects. Depending on η * , these effects can be significant, even turning an initial cold spot into a hotspot or vice-versa, and in particular, our benchmarks choices of η * lead to hotspots. Compared with a previous work involving such primordial particle production [19], the particle production in our model occurs only at a particular time η * , rather than periodically. We also show that the pairwise nature of particle production leads to an overlap of the two hotspots originating from the same pair, and this leads to a non-circular combined hotspot profile.
Varying M 0 /H * , g and the production time η * , we explore how the hotspot signal can be searched for, both using the CMB power spectrum and with a position-space search that selects out particularly hot pixels. We find the sensitivity for both methods is similar, but the position-space search is more powerful when g is large (∼ O(3 − 6)) and for larger η * 100 Mpc -scenarios where the spots are larger and hotter. This sensitivity is better than the one suggested by Ref. [23] which used an optimized N -point analysis rather than a position-space search, though direct comparison is somewhat tricky since they work in the context of a model with periodic particle production. Furthermore, the optimal estimator in Ref. [23] assumes no overlap between different hotspot profiles, a feature that is absent in our model and analysis. It would be very interesting to see how an optimized N -point analysis performs for non-periodic particle production and with overlapping hotspot profiles.
There remains several important and interesting future directions. We have not considered in detail the effects on non-Gaussianity coming from such events of particle production. The resulting bispectrum or trispectrum would exhibit "bump" features on specific scales, and therefore new templates of shape function would be needed to detect/constrain such kinds of non-Gaussianity.
The position space search we employ is a first step and involves several approximations. Most significantly, we assume the relevant astrophysical point sources have been masked and instrumental noise has been subtracted out from the CMB map. We believe that aspects of our signal, such as the perfectly circular nature of individual hotspot profiles in the CMB, and their appearance in a limited signal frequency range would help differentiate them from normal astrophysical point sources, though this deserves further study. Our position-space study is also simplistic and does not take advantage of the profile shape (which we have computed theoretically). Further analysis, better exploiting the signal properties, as well as an analysis using wavelet analysis and machine learning to better pick out the signal over the CMB background, is underway and will be discussed in a future work.