New inflationary probes of axion dark matter

If a light axion is present during inflation and becomes part of dark matter afterwards, its quantum fluctuations contribute to dark matter isocurvature. In this article, we introduce a whole new suite of cosmological observables for axion isocurvature, which could help test the presence of axions, as well as its coupling to the inflaton and other heavy spectator fields during inflation such as the radial mode of the Peccei-Quinn field. They include correlated clock signals in the curvature and isocurvature spectra, and mixed cosmological-collider non-Gaussianities involving both curvature and isocurvature fluctuations with shapes and running unconstrained by the current data analyses. Taking into account of the existing strong constraints on axion isocurvature fluctuations from the CMB, these novel signals could still be sizable and potentially observable. In some models, the signals, if observed, could even help us significantly narrow down the range of the inflationary Hubble scale, a crucial parameter difficult to be determined in general, independent of the tensor mode.


Introduction
An axion is a pseudo-scalar field that has a discrete shift symmetry, a ∼ = a + 2πf a , with f a known as the decay constant.Despite its apparent simplicity, an axion could have a plethora of phenomenological and cosmological applications, ranging from solving the strong CP problem in the Standard Model (SM) of particle physics [1][2][3][4], to serving as a classic cold dark matter (DM) candidate [5][6][7].It could be realized naturally as a light pseudo-Nambu-Goldstone boson from the spontaneous breaking of a global Abelian symmetry, the Peccei-Quinn (PQ) symmetry, U (1) PQ , at the high energy scale f a .Though the axion scenario was proposed more than 40 years ago [1][2][3][4], it continues to serve as one of the most motivated feebly-coupled particles beyond the SM, generating growing theoretical and experimental studies.For some recent reviews on axions as well as discussions of future plans and open questions, see [8][9][10][11][12][13].
A light axion, either a QCD axion or a general axion-like particle, could have a close interplay with the inflationary paradigm, one leading explanation for the origin of our observable universe [14][15][16][17].In particular, if U (1) PQ is spontaneously broken during inflation, the resulting massless axion would be present and undergoes random quantum fluctuations with an amplitude set by the inflationary Hubble scale, H [18][19][20][21][22]. These fluctuations, usually referred to as axion isocurvature perturbations, are independent of the primordial inflaton fluctuations.After inflation, they could be imprinted on the DM perturbations, once the axion becomes (part of) DM, while the inflaton perturbations are inherited by radiation, baryons and other DM components.Given the current non-observation of isocurvature modes in the cosmic microwave background (CMB) observations, there exists a strong constraint on the combination of the axion DM fraction and the size of the axion isocurvature perturbations [23].In the QCD axion scenario, the constraint is often interpreted as the "axion isocurvature problem": inflationary QCD axion DM, if it is the leading component of DM, is incompatible with high-scale inflation (e.g., f a ∼ 10 11 GeV, which will make the relic abundance of QCD axion DM through the minimal misalignment mechanism match the observed DM abundance, requires the inflationary Hubble scale to be below 10 7 GeV [23]).In this paper, we will not intend to provide a new solution to this problem, remedying the relationship between the QCD axion DM and high-scale inflation.Instead, we take a different perspective and address the following question: given the existing constraint on the size of axion isocurvature perturbation, could there exist other non-trivial cosmological observables testing the existence of light axions during inflation?
To answer this question, we make use of the recent developments in cosmological collider physics and primordial standard clock models, which provide new opportunities for probing high energy physics during inflation.
During inflation, particles with masses around and below the Hubble scale are produced on-shell.It has been pointed out that these heavy spectator particles could generate non-analytical momentum dependence and angular dependence of the primordial non-Gaussianities (three-and higher-point correlators) [24,25], due to their quantum motions.Observations of primordial non-Gaussiantity (NG) could provide rich information about the mass and spin of the heavy particles.This program, by its analogy to the massspin measurements at terrestrial particle colliders, has been dubbed "cosmological collider" (CC) physics.It has attracted rapidly growing attention , given its tremendous potential to study new physics at a much higher energy scale, which could never be produced at a particle collider.
When sharp features are present in the inflaton trajectory (see [98][99][100][101][102][103][104][105][106][107] for motivations from CMB data analyses and see [108][109][110] for reviews on primordial features), the energy scale of the cosmological collider can become even orders-of-magnitude higher.In the primordial standard clock models [105,106,[111][112][113][114][115][116][117][118][119][120][121][122][123][124], the inflaton deviates from the attractor solution temporarily due to some sharp feature in the trajectory.The departure could trigger the oscillation of a massive spectator field with a mass m ≫ H.As a consequence, in the primordial spectra (two-and higher-point correlators) of the density perturbations, this oscillation gives rise to a clock signal with the frequency sensitive to the background scale factor evolution a(t) and the particle mass m in Hubble unit, providing a way to probe these two kinds of information.The energy injected by the features can also excite quantum oscillations of heavy fields with masses m ≫ H, again leaving their signatures in the primordial non-Gaussianities [88].
In this article, we show that due to interactions between the inflaton and the PQ scalar field, whose phase becomes the axion after the breaking of U (1) PQ , there arises a whole suite of novel cosmological observables, that allow us to test the existence of the axion during inflation as well as its properties such as the associated PQ breaking scale during inflation.They include correlated clock signals in the curvature and isocurvature spectra; and large (correlated) CC NG signatures in the bispectra involving isocurvature modes with shapes and running that are unconstrained in data analyses so far.These bispectra include the mixed ones between curvature mode(s) and one or two isocurvature modes and pure isocurvature bispectra, with all three modes being the isocurvature ones.As far as we know, Ref. [78] is the only paper that studied the CC signal in an inflationary axion model in the pure isocurvature channel.Early papers on axion isocurvature NG but not involving the production of heavy particles with m ≳ H as in the CC physics include [125,126].In [78], to achieve a sizable CC signal in that model, the axion needs to have a non-zero classical rolling speed, or in other words, the axion needs to have a mass about the inflationary Hubble scale instead of being effectively massless during inflation.While this is plausible through model building, we will not pursue it in our paper.Instead, we explore a wider range of light inflationary axion models with various couplings and show that with either classical primordial features or chemical-potential type coupling, potentially observable new classical or quantum signals in either power spectra or bispectra could arise, in the parameter space consistent with current observations.The paper is organized as follows.In Sec. 2, we will review the basics of axion DM isocurvature.In Sec. 3, we will provide a summary of the three models based on effective operators and their main signatures.In the following three sections, Sec. 4, 5 and 6, we present details of each model and relevant key computations.We conclude and point out some future directions in Sec. 7.

Review of axion dark matter isocurvature
During inflation, if PQ symmetry is spontaneously broken, a massless axion, denoted by a or equivalently θ = a/f I with f I the PQ field's vacuum expectation value (VEV) during inflation, is generated.The condition for symmetry breaking requires f I ≫ H, the Hubble scale during inflation.The amplitude of its fluctuation with ⟨θ⟩ = θ i being the initial misalignment angle, is set by the Gibbons-Hawing temperature [127]: This fluctuation could be independent of that of the inflaton and leads to the DM isocurvature perturbation once the axion becomes (part of) cold DM after inflation.The isocurvature fluctuation in DM S d and the contribution due to the axion component S a are defined by where δ denotes fluctuation over the average value, while the subscript rad means photon radiation.The isocurvature fluctuation in DM due to the axion component is given by [128,129] where assuming δθ ≪ θ i , we obtain the axion isocurvature amplitude S a ≈ 2δθ/θ i .The dark matter isocurvature, S d is a product of S a and γ, the fraction of axion DM relic abundance Ω a h 2 among the total DM abundance Ω d h 2 .Throughout the article, we work in the framework of single field inflation and assume that the axion isocurvature is the only source of DM isocurvature.We do not consider the more complicated multifield inflation scenario in which DM isocurvature could be generated during reheating of multiple inflatons [130].Consequently, the remaining part of DM are assumed to be generated from radiation and will not contribute to the DM isocurvature S d .The expression for Ω a in the last line, adapted from [129], is for QCD axion through the vanilla misalignment mechanism [5][6][7].Note that Ω a depends on f a , the axion decay constant after inflation. 1 As we will discuss in more detail, f I is not necessarily equal to f a in the presence of the PQ field coupling to the inflaton.The isocurvature spectrum is then given by where the amplitude is given by and n i is the isocurvature spectral index.The current CMB measurements constrain uncorrelated DM isocurvature to be [23] 1 Beyond the minimal misalignment mechanism, mechanisms such as early matter domination [131][132][133][134][135], thermal friction from gauge fields [136][137][138], or dynamical PQ scale [139,140] could significantly alternate Ωa, which we will not explore further in the paper.Alternatively, an axion-like particle could also go through the misalignment mechanism and become DM.The axion-like particle's mass is not fixed by the standard model QCD scale, introducing an extra free parameter when calculating its relic density.In all scenarios above, the corresponding decay constants could be significantly higher than in the vanilla case.Elaborate discussions for these possibilities will be left to future studies.

Models and effective operators
We consider the following PQ field model during inflation with the (−, +, +, +) signature: where ϕ is the inflaton with a potential V ϕ (ϕ), χ is the complex PQ scalar field with a potential V χ (χ), and f a is the axion decay constant or equivalently the VEV of the PQ scalar today.V add includes the leading-order high-dimensional operators coupling the PQ field to the inflaton or other heavy spectator fields (e.g., we will consider some heavy fermions charged under U (1) PQ , ψ) during inflation, which could lead to different types of interesting spectrum or CC signals.We require these operators to respect the (approximate) shift symmetry of the inflaton and the PQ symmetry U (1) PQ .The operators and their corresponding cosmological signals in density perturbations are listed in Table 1.
Model New cosmological signals  In each model in Table 1, we have only one high-dimensional operator coupling the PQ field to the inflaton for simplicity.In principle, one could combine multiple operators in one model, which may lead to more complicated new CC signals.Note that V add in model 1 is a dimension-6 operator, in contrast to the dimension-5 operators in the other two models.The reason for this slightly odd ordering is that model 1 could lead to a striking signal, the correlated clock signals in the curvature and isocurvature power spectra, which contain most information and could probably be most easily tested with future cosmological observations.Also, although more suppressed by the cutoff scale Λ, the amplitude of signal in model 1 is not necessarily more suppressed than that in model 2 due to different dependence on other parameters.
In all the models, we focus on the scenario in which the PQ symmetry is spontaneously broken so that χ acquires a VEV, f I , during inflation.As we will show, f I could be shifted from f a in the presence of V add .Then χ could be parameterized as where σ is the heavy radial mode and the massless phase mode a serves as the axion, which could become (part of) the cold DM after inflation.As will be discussed in the models, we expect that the radial mode mass, about the same size of f I , is much above the inflationary Hubble scale.The amplitude of cosmological production of heavy particles, which is σ in our case, is proportional to e −πm/H with m the mass of the heavy particle (see also [141,142]).Thus the corresponding CC signal gets exponentially suppressed when m ≫ H.This clearly challenges the power of CC probe to super-H scale and the inflationary axion scenario.To overcome this difficulty, we consider two possible enhancement mechanisms which have already been proposed in the CC literature.The first one is to employ the primordial features, some shift-symmetry violating features along the inflaton trajectory in the inflationary landscape.These features could excite a heavy field with m ≫ H classically and quantum-mechanically, resulting in sizable signals in the power spectrum and non-Gaussianities [88,111].The other one is to use a chemical-potential type coupling ∂ µ ϕJ µ , with J µ a current made of the heavy fields [51,61,62,64,69,72,142].This coupling still preserves the inflaton shift symmetry but can introduce a new scale φ0 /Λ ≫ H when the inflaton is set to its homogeneous background value [51,61,62,64,69,72,142].This new higher scale allows excitation and production of heavy particles, compensating the Boltzmann suppression with e π φ0 /(ΛH) and boosting the NG signals.In our model 1 and 2, we will rely on the primordial feature mechanism, while in model 3 we employ the chemical-potential coupling.The Lagrangian of our first model is given by -6 -where c is a positive dimensionless coefficient and Λ is the cutoff energy scale of the model. 3hroughout this work, we assume such model parameters to be constants independent of ϕ for simplicity.In the second line above, we use the parametrization in Eq. (3.2) after the PQ symmetry breaking, and V σ (σ) is given by Note that in this model, the axion field does not couple directly to the inflaton.The VEV of the PQ field during inflation is given by where φ0 is the homogeneous background of the inflaton (the standard attractor solution), satisfying The equation above assumes the slow-roll approximation with φ ≪ ∂V ϕ ∂ϕ , and the slow- 18 GeV the reduced Planck scale.In single field inflation, φ0 ≈ (60H) 2 , fixed by the amplitude of the density perturbations.The mass squared of the radial mode is For the model to be a valid effective field theory (EFT), we need to impose the following constraints following approaches in [49,143]: (a) The added inflaton-PQ field coupling in this model is just the leading one of a series (∂ϕ) 2n |χ| 2 /Λ 4n−2 .Requiring this power expansion doesn't spoil the EFT during inflation, we need to impose [144] Λ > φ0 .
(b) Naturalness constraint on the quantum correction to λ from the inflaton loop with the cutoff ∼ Λ: (c) Naturalness constraint on the quantum correction to m 2 σ,eff from the inflaton loop with the cutoff ∼ Λ: (d) The high-dimensional operator only leads to a small correction to the inflaton kinetic term of (∂ϕ) 2 : Eq. ( 4.3) and Eq.(4.5) allow us to solve f I /H and m σ,eff /H as a function of f a /H, q and λ.For the natural value of λ ∼ 1 and different choices of f a /H ⊂ O (10), we show m σ,eff /H as a function of q in Fig. 1.One can see that when q increases, the contribution to m σ,eff /H from the inflaton-PQ field coupling dominates over that from the PQ-field potential, approaching an asymptotic value (2λ) 1/4  φ0 /H at q = 1.Smaller λ allows for larger f a /H to achieve the same m σ,eff /H.In all the models we consider, m σ,eff = √ λf I ≫ H during inflation, the cosmological production of σ, and the induced CC signals will be exponentially suppressed by ∼ e −πm σ,eff /H , as mentioned in the summary of Sec. 3. To enhance the signals, we need to introduce some boosting mechanisms.In model 1 and 2, we will rely on the primordial feature mechanism in which the inflaton potential possesses a sharp feature and can be decomposed as where V ϕ0 is the smooth potential responsible for the featureless attractor solution in Eq. (4.4), while the perturbation V ϕ1 is a small sudden change of the potential localized in some regions of the field space.Then the classical background trajectories can be split into the featureless (with subscript 0) and featured (with subscript 1) components correspondingly: where ϕ 0 (t) is given in Eq. (4.4) and σ 0 (t) = 0. Note that since a does not couple to the inflaton directly, a bkg (t) is a constant.There are many examples of sharp features.Such features can be easily envisioned when we embed the inflation models in a potential landscape when the universe is unstable.
There are also motivations from CMB data analyses [98][99][100][101][102][103][104][105][106][107].A sharp feature could naturally excite the classical oscillation of a massive field, such as the radial mode of the PQ field, which is otherwise hard to probe because its mass is much larger than the Hubble scale.In this paper, as a toy example, we consider V ϕ1 to be a sharp step function: where b ⊂ (0, 1) is a dimensionless quantity; θ(ϕ − ϕ s ) is the Heaviside θ function: it is one when ϕ > ϕ s and zero otherwise, which means that the perturbation potential is a small but sharp downward-step.Similar or other examples of sharp features with a more complicated form have been used to explain the CMB residual anomalies [98][99][100][101][102][103][104][105][106][107].Here we just use this simplest form to illustrate the new clock and CC observables.
The time evolution of ϕ 1 and σ 1 due to V ϕ1 follow the equations of motion (EOMs) below: where we approximate φ2 bkg − φ2 0 ≈ 2 φ0 φ1 .Given Eq. (4.12), the solutions are ) where t s corresponds to the time when the inflaton passes through ϕ s and we ignore higherorder terms suppressed by more powers of m σ,eff /H ≃ µ σ in σ 1 (t).An illustration of this scenario is presented in Fig. 2. The inflaton velocity has a steep jump at t = t s with ( φ2 bkg − φ2 0 )/ φ2 0 | t=ts ≈ 2bV ϕ0 / φ2 0 .This triggers oscillation in the radial mode σ with an amplitude proportional to c and suppressed by m 2 σ,eff or equivalently µ 2 σ when µ σ ≫ 1.

Corrections to the power spectra
Now we move on to consider the perturbations, defined as δϕ = ϕ − ϕ bkg , δσ = σ − σ bkg , and δa = a − a bkg , and compute their two-point correlators.Expanding Eq. (4.1), we have the following leading quadratic terms where R is the scale factor; the overhead dot indicates the derivative with respect to time while ∂ i is the derivative of the spatial coordinate.We ignore the contribution from the inflaton kinetic term, which does not contribute to the power spectra corrections.We also ignore terms of higher order O((σ bkg /f I ) 2 ) that are highly suppressed by small parameters.Finally, terms proportional to δσ 2 are dropped as they only contribute to the two-point correlators at higher order in c/Λ 2 ; the δ φδσ term also contributes to two-point correlators at higher order, but we include it here since they will contribute to the three-point functions which we will discuss later.
Sharp feature signal in curvature power spectrum We first compute the leading correction to the curvature power spectrum due to the step function potential in Eq. (4.12), through the first term in Eq. (4.17): where τ is the conformal time with dτ = dt/R and R(τ ) = −1/(Hτ ).τ s is the corresponding conformal time at the transition point ϕ s .
An insertion of this interaction leads to the diagram in Fig. 3, resulting in a correction δϕ ∂ 2 ϕ V ϕ δϕ to the curvature power spectrum P ζ as where is the reference wavenumber for the sinusoidal signal; the function g(k/k 0 ) determines the envelope of the signal and ω gives the phase.This sinusoidal correction is entirely due to the step feature in time and is independent of the excited oscillating heavy σ field.Thus such a signal does not contain any information about the PQ sector and its coupling to the inflaton.Both g(k/k 0 ) and ω are model-dependent, sensitive to the input perturbation potential V ϕ1 .Their detailed forms are not important for our following discussions, so we will not list them here.Instead we will focus on the amplitude at k ∼ k 0 , which is of order There is a suggestive dip feature in the power spectrum present near ℓ ∼ 25 in the CMB, which could be interpreted as a sharp feature signal with ∆P ζ /P ζ ∼ 0.3 [98,106].Taking this as a rough guideline, we consider bV ϕ0 / φ2 0 ≲ 0.3 for evaluating the more important and interesting clock signals below.The abrupt starting point of the oscillation in σ bkg at t = t s is also a potential source of sharp features.Due to the couplings in the 2nd and 3rd lines of (4.17), this sharp feature contributes to additional sharp feature signals in both curvature and isocurvature power spectra.We will not present any details of such signals in this paper, only to mention that common to all sharp feature signals, these signals would also have the sinusoidal-running scale dependence, and model-dependent envelops.The amplitudes of these signals would be smaller than the ones we considered above, due to the sizes of the couplings and the milder nature of the sharp feature.
Correlated clock signals in the power spectra Now we consider the corrections to either the curvature or isocurvature power spectra due to the resonances between the induced oscillating heavy radial mode and the axion or inflaton, which result in clock signals with rich information for the PQ field and its couplings during inflation.The leading diagrams, from interactions in the second and third lines of Eq. (4.17) with insertions of σ bkg , are presented in Fig. 4. We take σ bkg to be Eq.(4.16) and compute the two diagrams in the same way: where k r = −µ σ /(2τ s ) = µ σ k 0 /2.We remind that the clock signals are generated through the resonance between the harmonic (standard) oscillation of the massive field and the oscillation of the inflaton or axion quantum fluctuation mode [111,112].Because the time dependence of the latter frequency is determined by the background scale factor evolution a(t), the momentum dependence in the phase of the clock signal takes the form of the inverse function of a(t), which in this case is a logarithmic function.The √ 2πµ σ enhancement arises from the resonance [145][146][147][148].Note that ∆P i the same up to an overall rescaling factor cf 2 I /Λ 2 , which originates from the different coefficients of the quadratic terms for (δa) 2 and (δϕ) 2 in Eq. (4.17).
Let's estimate the amplitudes of the corrections at k = k r .In the limit µ σ ≫ 1, where we use µ σ ≈ m σ;eff /H = √ λf I /H and q ≡ cf 2 I /Λ 2 ≪ 1 to satisfy the constraint on the kinetic term of ϕ as discussed around Eq. (4.9).
-12 -Estimates in the more general parameter space are shown in Fig. 5.We fix the initial misalignment angle θ i to be one and λ to be an order one value, either 1 (top row) or 0.3 (bottom row), and present contours of physical quantities in the plane of the two most important input model parameters f a and q.We also fix f a /H = 40, but the result won't change much when varying f a /H as long as it is O (10).For f a /H ≳ 100, µ σ = √ λf I /H > √ λf a /H also becomes ≳ 100 (for λ ∼ O( 1)) and makes it more difficult for the clock signals to be observable in CMB due to the high frequency.Here we assume that the axion is the QCD axion, and the relic abundance is computed according to the last line of Eq. (2.4).
Note that Ω a ∝ f 7 6 a , a higher f a will give rise to a larger axion fraction γ.Combining the left and right panels, we see that, for example, when bV ϕ0 / φ2 0 = 0.3, q ≳ 0.01 and f a ≲ a few ×10 9 GeV, 4 |∆P ζ /P ζ | ≳ 2% with µ σ of order a few 10's, which could be observable, while the current isocurvature constraint in Eq. (2.7) is still satisfied.In this region, the axion DM fraction γ ∼ O(10 −3 − 10 −4 ).
From the computations above, one can see that our model 1 could have a striking signal in the power spectra, namely the correlated clock signals, with the following characteristics: • Sizable corrections to both the curvature and isocurvature power spectra, both taking the form sin (µ σ ln(k/k r ) + phase) with a common mass parameter µ σ , a common reference scale k r , and exactly the same phase.
• The correction to the isocurvature two-point correlator is parametrically larger by a factor of 1/q.This signal is illustrated in Fig. 6.Such a signal, if measured, could be very informative: • The correlated clock signals would point to an oscillating massive field coupled to both the inflaton and the dark matter field, and therefore serve as a more distinctive signature of the axion-inflaton Lagrangian in Eq. (4.1) than other more generic predictions of CC physics associated with a massive field.
• It could help us determine µ σ (or equivalently, m σ;eff in Hubble unit) and q, and learn about the energy scale related to the PQ field and its coupling to the inflaton.
• We could also narrow down the range of the inflation energy scale from these measurements in this particular model.There are degeneracies between f I and H in two observables, namely, the amplitude of the axion isocurvature perturbation (combining Eq. (2.4), (2.7) and (4.3)), and the mass of σ in Hubble units µ σ = √ λ(f I /H) from the frequency of the clock signals, assuming QCD axion dark matter with the misalignment mechanism.The measurement from either the axion isocurvature perturbations or the primordial standard clock/cosmological collider physics alone would not be able to determine H over many orders of magnitude.However, these two degeneracy directions in the f I -H plane are quite orthogonal to each other.Combining these two observables (together with the determination of q from the relative sizes of the curvature/isocurvature clock signals 5 ), we would have a very rare opportunity of being able to significantly narrow down the energy scale of the inflation model, better than either type of the experiments alone in terms of order of magnitude, especially for the low-scale (small-field) inflation models of which the tensor mode would be too small to observe. 6 Interestingly, if this clock signal would be observable in the curvature power spectrum in the near future, which means ∆P ζ /P ζ should be at least a few percent [105,106], the correlated oscillatory signal would be present at the leading order of the isocurvature power spectrum because q ≪ 1.In certain parameter space, the oscillatory signal could even become the leading component of the isocurvature power spectrum over the scale-invariant one. 7When a scale-dependent oscillatory signal becomes important, the constraint on the isocurvature perturbation from the Planck data may need to be re-analyzed with oscillatory templates, because such a signal would not be picked up by a scale-invariant template.
• Both clock signals are examples of classical primordial standard clock signals.The phases of their oscillations directly record the time dependence in the background inflationary a(t), and therefore provide direct evidence for the inflationary scenario.
We refer the readers to the following literature for detailed studies of this aspect [105,106,[111][112][113][114][115][116][117][118][119][120][121][122][123][124].The most stringent constraints on primordial feature models so far are from the CMB data.Many feature models have been compared to Planck data and constrained using the curvature power spectrum and non-Gaussianities (see [110,150] for reviews).With the addition of many other experiments and different types of cosmological data, we expect rapid progress in this area in the near future.Future improvements from observations of CMB E-mode polarization are forecasted in [107].Those from galaxy surveys are studied in [151][152][153][154][155][156][157][158][159][160][161][162][163].More futuristically, one can also look forward to 21cm hydrogen line surveys [164,165] and stochastic gravitational wave background mapping [124,166,167].So, the type of features in the curvature power spectrum, such as the benchmark model we presented in this work, will be critically tested very soon.On the other hand, experimental prospects of primordial features in the isocurvature power spectrum are now an open question.

Three-point correlators
It is well-known that primordial features boost NGs of inflation models [145][146][147][148][168][169][170][171][172][173].With the presence of the axion isocurvature during inflation, there could also be large NGs in three-point correlators involving isocurvature modes.We will compute both in this section.The relevant trilinear perturbation terms for leading-order bispectra include Curvature bispectrum.A sharp feature in time could not only generate a sinusoidal signal in the curvature power spectrum, but also in a curvature bispectrum (which we will denote as ccc) through the first trilinear coupling in Eq. (4.24):The corresponding diagram with an insertion of this interaction in Fig. 7 could be -16 -computed as where ⟨δϕ 3 ⟩ ′ is the normalized three-point inflaton correlator without the delta function: 2 are the model-dependent envelope functions.The g 1 , g 2 functions from the simple step function potential we assume should be taken with a grain of salt since they may not be representative of more realistic models.Yet the amplitude still allows us to estimate the strength of NG measured by the dimensionless parameter of the ccc correlator, f ccc NL , where we use that the curvature perturbation ζ is related to δϕ via ζ = −(H/ φ0 )δϕ and the ccc correlator is normalized by . Though it could be a large signal, this sharp-feature-induced ccc NG does not tell us anything about the isocurvature mode, and we will not study it further.We also ignore high-order diagrams involving an exchange of the radial mode.Mixed isocurvature-curvature bispectrum.The quadratic couplings in Eq. (4.17) and the trilinear ones in Eq. (4.24) combined could contribute to curvature-isocurvature-isocurvature three-point correlators, which we will indicate as the "cii" bispectrum.The leading types of diagrams are shown in Fig. 8 and Fig. 9. Fig. 8 has an insertion of φ1 , while the two diagrams in Fig. 9 have an insertion of the oscillating σ bkg at either two-point or three-point vertex respectively.Note that in model 1, δa's always come in pairs, so only bispectra with an even number of isocurvature modes exist.In other words, only ccc or cii exist.In this section, we will compute the cii three-point correlator with momenta k 1 , k 2 , k 3 , where k 1,2 are taken to be associated with the axion modes and k 3 is assigned to the inflaton ϕ.As we will see, the leading behavior of this correlator is a scale-dependent oscillatory function of the momentum k 1 + k 2 + k 3 .Such a signal would be most easily detected in the equilateral limit where there are more triangle configurations.So we will focus on the equilateral limit with k 1 = k 2 = k 3 ≡ k in this section.Define the normalized three-point function (with a prime superscript) as: The contributions from Fig. 8 and the left panel of Fig. 9 include the following integral in common: where j = 3 or 3 2 ± iµ σ depending on whether φ1 or σ is injected.The massive mode function is given by v k (τ ) = −ie −πµσ/2+iπ/4 √ πH(−τ iµσ the Hankel function of the first kind.This integral corresponds to having the feature injection at τ 2 to excite the heavy σ quantum before it decays to the axion quanta at τ 1 .The integration of the two-point vertex at τ 2 reads: where z 1,2,s = −k 3 τ 1,2,s = −kτ 1,2,s , respectively.We also use the Hankel function of the second kind, H iµ (z) * e −πµ .In the regime with k/k 0 ≳ µ σ ≫ 1 we are interested in, the early-time expansion of the Hankel function could be applied.The twopoint integral could then be approximated as where Γ(a, b) is the incomplete Gamma function.
The three-point vertex at τ 1 has the following kinematic structure: where k 12 ≡ k 1 + k 2 and the operator D is defined as with the partial derivative ∂ k 12 acting on k 12 .The integration involving the three-point vertex in the equilateral limit can be simplified as: iµσ (z 1 ), = (−1) iµσ (z 1 ), ( where we set τ end = 0.In the last step, we apply the equilateral limit with For convenience, we denote the integration of Eq. (4.29) as: where X 3 and X 2 are abbreviations of the three-and two-point vertexes, respectively.Combining Eq. (4.31), (4.34) and applying the early-time expansion in the integration of X 3 , we have the following leading-order analytic approximation for the j = 3 case ( φ1 injection in Fig. 8): One could see that this contribution is a sinusoidal function with frequency set by 1/τ s .The computations are similar for the left panel of Fig. 9, aside from the different forms and amplitudes of the injected features.However, numerically we find that the NG from the σ bkg injection is subdominant compared to that of the φ1 injection.First, the amplitude of |σ bkg |/f I ∼ q(λ φ2 0 /m 4 σ,eff )(bV ϕ0 / φ2 0 ) is usually smaller than φ1 / φ0 = bV ϕ0 / φ2 0 except for very light σ.Moreover, since the frequency of the injected feature, µ σ , is the same as the propagator's mass, the oscillating σ bkg field cannot excite a propagating δσ perturbation with the same mass and light fields simultaneously.In contrast, the φ1 insertion is approximately infinitely sharp by construction and is capable of exciting this -19 -massive field. 8Finally, from the right panel of Fig. 2, one can see that σ bkg is continuous, in contrast to φ1 .In this case, the feature of σ bkg is not as sharp as that of φ1 .The final integrand is a highly oscillatory function with its both ends at z = 0 and z s being zero, making the overall integration small.Another contribution to the diagrams comes from the integral where the three-point vertex happens at τ 2 prior to the two-point one at τ 1 , namely: or in the shorthand notation , where X2,3 stands for the twoand three-point vertexes.Numerically, we find the size of the integral above is subdominant or, at most, comparable to that of Eq. (4.29).Similar to what is explained in [88] (for the diagrams Fig. 2b and 2c in [88]), the suppression is because this integral corresponds to the situation in which the massive field is created before the time of sharp feature.
We now estimate the size of contributions from Fig. 8 to cii bispectrum, up to a numerical factor, as The three-point function above can be translated into the dimensionless form of NG parameter, f cii NL , via the following relation: where we use S d = γS a = 2γδa/(f I θ i ), ζ = −(H/ φ0 )δϕ and the definition of β in Eq. (2.7).
Normalizing to the standard curvature NG parameter to take into account that the isocurvature power is suppressed, the contribution from Fig. 8 to the cii NG in the equilateral limit is estimated to be: The right panel of Fig. 9 also contributes to cii through the integration, in the equilateral limit: 2 )H (2) 2 )H (2) where σ bkg = σ 1 with j = 3 2 ± iµ σ is injected through the three-point vertex at time τ 2 prior to τ 1 .With the early-time expansion, the integration of the two-point vertex over z 1 reads: However, when combined with the integration of the three-point vertex, the logarithmic behavior above makes it difficult to obtain an analytical approximation.Moreover, the positive-and negative-frequency parts of the σ bkg oscillation largely cancels with each other, which is also challenging to be approximated analytically.To have a rough estimation in the early-time regime, it is convenient to replace log z 2 by a constant such as 1.With such a simplification, the integration over the three-point vertex becomes a product of polynomials and sinusoidal oscillations.By comparing with the numerical integration of Eq. 4.41, we find the closed-form approximation of the cii bispectrum to be

.43)
Notice that the approximation above tends to overestimate the three-point function when µ σ is large.When z s ≫ µ σ ≫ 1, this contribution to the size of f cii NL could be estimated as In Fig. 10, we show the numerical results for a benchmark model with a relatively small µ σ = 5 in the equilateral limit.For both types of diagrams in Fig. 8 and Fig. 9 (right panel), the cii signals are sizable.The analytical approximations in Eq. (4.39) and (4.44) using the early-time expansion are also plotted as dashed curves, agreeing decently well with the numerical results in both the amplitude and the phase.For more realistic benchmarks with higher µ σ (which are difficult to be evaluated numerically due to the slow convergence), we expect that the diagrams with φ1 injection are not significantly affected as their amplitudes -21 -are independent of µ σ in the early-time region.In contrast, the importance of the diagrams with σ bkg injection will decrease with a larger µ σ as their amplitudes ∝ |σ 1 | are suppressed by µ −2 σ .In all cases, the phase of the cii signal follows e −3ik/k 0 = e i(k 1 +k 2 +k 3 )τs up to a constant shift.
In the other interesting limit, the squeezed limit with k 1 ∼ k 2 ≫ k 3 , the NG signals due to the sharp feature could also be sizable, which we check by numerical computations.The properties of the contributions from each type of diagram, such as the relative importance, are similar to those in the equilateral limit.Numerical integration also shows that the k 3 /k 0 dependences of signals are similar to their equilateral-limit counterparts.However, it is more difficult to find analytical approximations for both the amplitude and the phase of the signal, which we will not detail here.
The sharp feature signals in the curvature bispectrum and mixed bispectrum studied in this subsection, as well as the sharp feature signal (4.18) in the power spectrum, are all due to the same sharp feature in the model, and so they all run sinusoidally (with 2k versus k 1 + k 2 + k 3 between power spectrum and bispectra) with the same starting scale, parameterized by k 0 .Besides these leading behaviors, the more detailed envelope behavior and subleading phase shifts are more complicated and depend on the nature of the sharp feature, but correlations between different spectra may still deserve further studies.In contrast to the clock signals, as mentioned, only from these sharp feature signals we will not be able to learn any specific information about the massive radial mode; thus these signals are less distinctive as signatures of the PQ field.These comments also apply to the sharp feature signals in the next model.
Regarding the observational prospect of these scale-dependent oscillatory bispectra, we note that, although we have shown that the sizes of these bispectra are greatly enhanced relative to their attractor values due to the features, it remains an open question of how large these bispectra have to be to become observable in future experiments.There are some constraints on such curvature bispectra from the Planck data [23,172,173]. 9In addition, the Planck data also sets some weak constraints on mixed curvature/isocurvature bispectra assuming local shapes, which do not apply to the shapes in our model [175].The future observational prospect is less studied.On the other hand, both the constraint and the future prospect of the oscillatory mixed bispectra are open questions.We leave these questions to future works.In model 2, we consider the following Lagrangian, which introduces a kinetic mixing between the inflaton and the radial mode, after the PQ symmetry breaking: where in the second line, we use the parametrization in Eq. (3.2) and Eq.(4.2).In the third and fourth lines, we employ the following field redefinition to remove the kinetic mixing where the dimensionless parameter x is defined as For the redefined canonically normalized fields, we could use the standard Bunch-Davies solutions for the mode functions.
In principle, the dimension-five operator added here could appear together with the dimension-six operator considered in the previous section.Since they lead to different signals and may come from different UV origins, we study them one by one for clarity and simplicity.As we will see, the signal magnitudes are not determined purely by the operator dimensions.
In the absence of a primordial feature, the homogeneous background solution of the inflaton is the same as in Eq. (4.4) with ϕ 0 replaced by φ0 .Around the PQ field's VEV f I during inflation, the background radial mode satisfies σ0 = 0 , ( which requires Using Eq. (4.2) and Eq.(4.4), one obtains a relation determining f I as well as the effective mass squared of σ: and where we ignore the suppressed contribution to The axion field adopts a trivial background solution a 0 (t) = 0. Now as in model 1, we add a step feature, V ϕ1 as in Eq. (4.12), to the inflationary potential that sources the attractor solution, V ϕ0 .The homogeneous background fields become φbkg (t) = φ0 (t) + φ1 (t) , σbkg (t) = σ0 (t) + σ1 (t) = σ1 (t) , where φ0 (t) and σ0 (t) are the featureless solutions as discussed above while φ1 (t), σ1 (t) are the induced features in the inflaton and radial mode trajectories.The axion field -24 -is not directly coupled to the inflaton and its background evolution remains trivial with a bkg (t) = a 0 (t) as a constant.The linearized EOMs of φ1 (t) and σ1 (t) become: and Plugging in Eq. (4.12), the solutions of φ1 , σ1 read: Compared with the solutions in model 1 shown in Eqs.(4.15) and (4.16), the induced inflaton feature remains the same while the induced oscillating radial mode is suppressed by one less power of its mass.

Corrections to the power spectra
Let's first consider the correction to the power spectra.At the quadratic level, the interactions between the perturbations, δ φ = φ − φbkg , δσ = σ − σbkg and δã = ã − ãbkg = ã, read The first term in Eq. (5.13), − 1 2 , due to the step potential V ϕ1 is also present in model 1.It will lead to the same sharp feature signal in the curvature spectrum, computed in Eq. (4.19).We will not repeat the computation here.On the other hand, in contrast to model 1, there are no quadratic terms of (δ φ) 2 from the dimension-5 operator we add in model 2 since it is only linear in φ.Thus diagram as in the left panel of Fig. 4, leading to a clock signal in the curvature spectrum, is absent here.One may wonder that a higher-order diagram, such as the one in Fig. 11, could lead to a clock signal due to the resonance between σ bkg and the intermediate σ mode.However, it turns out not to be the case: the injection's frequency is the same as that of the intermediate δσ and there is no An example of a higher-order diagram contributing to the curvature power spectrum with the exchange of the heavy radial mode.The external points of inflaton are marked with blue squares.The solid blue and dashed black lines represent the propagator for φ and σ, respectively.The empty dot denotes the interaction proportional to x∂ 2 ϕ V ϕ in the first line of Eq. (5.13), while the shaded blob represents interactions with an insertion of σbkg in the third line of the same equation.
resonant production, which is the same reason for the lack of clock signal in cii in model 1 as discussed in the previous section.Clock signal in the isocurvature spectrum.Similar to model 1, there is still a clock signal in the isocurvature spectrum, due to the diagram in the right panel of Fig. 4. The computation is similar and the only change is σbkg , which is Eq. ( 5.12) for model 2. We find that where we take µ σ ≫ 1 and x ≪ 1 limits to simplify the final result.The amplitude (defined at k = k r ) is then: where we use µ σ ≈ m σ,eff /H ≈ √ λf I /H.This benchmark requires that the axion is only a subdominant fraction of dark matter with γ/θ i ≈ 10 −3 β 0.038 to satisfy the current isocurvature bound in Eq. (2.7).In certain parameter space (e.g.x ≳ 0.1 while fixing the other benchmark parameter values), the perturbative theory for the isocurvature spectrum could break down, and one needs to adopt a different method to compute its amplitude.
Comparing the isocurvature clock signal amplitudes in Eq. (4.23) for model 1 and Eq. ( 5.15) for model 2, we note that they could be similar even though x from the dimensionfive operator (suppressed by one power of the cutoff scale) could be much larger than q from the dimension-six operator (suppressed by two powers of the cutoff).This difference is compensated by different power dependence on other parameters, such as φ0 /H 2 ≃ 60 2 .This serves as an example that the signal sizes are not entirely determined by the dimensions of the added operators.

Three-point correlators
Analogous to model 1, the sharp feature also leads to a potentially large NG in ccc via the diagram in Fig. 7 and the result applies here as well.Similarly, the more interesting bispectrum in model 2 is cii, which we will discuss more in this section.The trilinear interactions of the perturbations read: (5.17) The first integral at τ 1 stands for the three-point vertex happening at late times, while the second integral is the mixing two-point vertex excited by the feature which happens at the earlier time τ 2 .The second integration is straightforward: e πµ σ /2 e ik 3 τs H (2) (5.18) The three-point vertex at τ 1 could be estimated in the same way as Eq.(4.34).
In the equilateral limit, where we apply the early-time expansion with k/k 0 ≳ µ σ , the analytical approximation for ⟨δaδaδ φ⟩ ′ to the leading order in x reads where we omit terms with lower powers of k/k 0 .The corresponding dimensionless amplitude of the NG signal for cii in the equilateral limit is then (5.20) Similar to the previous model in Section 4, the observational prospect is still an open question as it depends on the oscillatory templates and multiple parameters.
6 Model 3: The first two models rely on a primordial inflationary feature to enhance the signals either in the two-point or the three-point correlators.In model 3, we will explore a different enhancement mechanism based on a chemical-potential type coupling between the inflaton and the PQ field, as well as the PQ field coupling to some other heavy spectator fermion fields.We find that in this model, there could also be novel correlated signals in several bispectra involving different numbers of the isocurvature modes.We also consider the cross power spectrum between the curvature and axion mode.
In model 3, we add a different dimension-5 operator, ∂ µ ϕJ µ PQ , with J µ PQ ≡ i(χ † ∂ µ χ − χ∂ µ χ † ), which again simultaneously respects the U (1) PQ symmetry and the inflaton's shift symmetry, to the Lagrangian: where in the third line, we use the field redefinition as to remove the kinetic mixing between a and the inflaton.After the spontaneous breaking of U (1) PQ , the radial and phase modes of χ = (f I + σ)/ √ 2 exp(iã/f I ) are related to the basis we start with, χ = (f with the dimensionless parameter z defined as Both fields σ and ã have canonically normalized kinetic terms.We identify ã as the massless axion mode, with the same Bunch-Davies initial condition as the inflaton.In this new basis, one could see that with the chemical-potential term alone, the axion ã is decoupled from the inflaton.Thus there is no mixed correlation, i.e., curvature-isocurvature spectrum and mixed bi-spectra.This is consistent with the discussions in the literature that ∂ µ ϕJ µ coupling alone could not lead to any chemical potential for a scalar field (the axion in our case) [64,72,142].We extend the model in Eq. (6.1) by including a Dirac fermion ψ that couples to the PQ scalar χ: where D is the covariant derivative, and ψ L , ψ R are the two chiral Weyl components of ψ: The above Lagrangian will be invariant under the global PQ symmetry if the component Weyl fermions are also charged under U (1) PQ with opposite charges.Under a PQ rotation by an angle θ, the charged matter transforms as After the PQ symmetry breaking, ψ acquires an effective mass Performing a chiral rotation of the fermions to remove the phase in the Yukawa term, the fermion kinetic term yields: where we omit irrelevant terms for the CC signals.The heavy fermions ψ above could be charged under the SM color group to induce a coupling between the axion and the SM gluon as needed in the KSVZ axion model [176,177].In the presence of ∂ µ ϕJ µ , we identify ã as the axion with its homogeneous background satisfying ȧ0 = 0 during inflation.This implies a rolls during inflation with a speed ȧ0 = z φ0 .Then its coupling to the non-conserved axial current J µ 5 = ψγ µ γ 5 ψ gives rise to a chemical potential for the fermion When the chemical potential µ c ≳ m ψ ≃ yf I / √ 2, ψ with a particular chirality would be produced efficiently from the inflationary background, leading to sizable NG signals [51,142].Using the redefined basis in Eq. (6.3), the couplings between the fermion and the scalars read This could be achieved if the previous condition on y is satisfied and Λ ≲ 4πf a .In addition, z ≪ 1 so that κ 2 (∂ µ ϕ) 2 | χ| 2 /Λ 2 will not lead to a large negative contribution to the inflaton kinetic term.Finally, to avoid the last term in Eq. (6.1) restoring the PQ symmetry during inflation, we also need In practice, f I ≲ f a is possible with a certain level of fine-tuning.In general we expect f I to be of O( φ0 ) or higher.Since the fermion couples to both ã and ϕ as shown in Eq. (6.10), its loop introduces all kinds of CC bispectra signal in terms of curvature and isocurvature, as shown in Fig. 13.Since the coupling of the inflaton and ã to the fermions is proportional, we expect the bispectra signals to be related to each other in a simple way.For the pure isocurvature bispectra (iii), its NG parameter f iii NL can be estimated using the late-time expansion of -30 -fermion modes [51,61]: where we set N c and N ψ are ψ's number of color and flavor, respectively.For the QCD axion that couples to vector-like quarks, N c = 3 while N ψ ≥ 1.
For the other (mixed) bispectra, we have correlated signals with the same mass parameter and related amplitudes: Here the superscripts of f NL are ordered, with the first index corresponding to the mode with the lowest momentum in the squeezed limit.Depending on the sizes of z and √ β, the CC signals in each mixture vary.Fig. 14 shows quantitative results of f iii NL and f ccc NL .One could see that the amplitude of the curvature bispectrum is generally small, but in some special parameter space, f NL becomes much greater than O(0.01).Note that, compared to the first two models, the non-Gaussianities here are scale-invariant.Also the magnitudes of f NL are generally smaller due to loop effects.Therefore it would be very difficult to observe such non-Gaussian signals, although in general the constraints and detection prospects of these types of non-Gaussianities involving isocurvature external lines remain as open questions.On the other hand, the fermion loop also introduces a mixing between ã and ϕ, and thus the correlation between curvature and isocurvature modes at the one-loop level.However, since the scale-invariance is unbroken in the two-point functions, the fermion loop only appears as the local part, and the mixing is suppressed by m −2 ψ .A quick estimation of the mixed spectrum from the fermion loop reads: and the corresponding mixing angle is when m ψ ≳ H.Such a small correlation is thus very challenging for observations and is unlikely to cause more stringent bounds on β.
One could consider a similar model with the axion coupling to heavy spectator gauge bosons instead of fermions.But that type of model requires more ingredients to generate a sufficiently large CC signal.We leave the explanation and estimates of the signals to Appendix A.

Conclusions and outlook
In this article, we point out that if a light axion is present during inflation and becomes even only a small fraction of dark matter later, its associated isocurvature fluctuations could result in a profusion of novel cosmological correlators, depending on the interaction between the PQ field and the inflaton.These correlators are either two-point or threepoint functions involving different combinations of curvature and isocurvature modes.The corresponding signals include correlated classical clock signals, taking the form of sin(log k) with k the wavenumber, in both the curvature and isocurvature power spectra; and the cosmological-collider type NG in mixed curvature-isocurvature or pure isocurvature threepoint functions, which could be correlated as well.With the help of a classical inflationary feature or chemical-potential enhancement, these signals could be sizable and observable in the future, even in the presence of current strong constraints on the overall amplitude of the isocurvature power.They encode a wealth of information about the interplay between the PQ breaking dynamics and inflation at a very high energy scale, which could be challenging to be probed otherwise.In some models, the signals could be combined to infer the range of the inflationary Hubble scale, an important parameter which suffers from a lack of observables, independent of the tensor mode.
Our paper serves as an early step to explore the application of the booming cosmological collider and primordial standard clock physics to explore the axion dynamics during the inflationary epoch in the very early universe.There exists more to study along the direction: • The paper is based on three simple EFTs, each containing only one high-dimensional operator coupling the PQ field to the inflaton.A general EFT could contain several operators simultaneously.New diagrams could be present and lead to different patterns of the observables.
• In our study, we do not intend to solve the QCD axion isocurvature problem and ameliorate the tension between the high-scale inflation scenarios and QCD axion dark matter.There exist quite a few models to address this problem by including specific types of interactions between the inflaton and the PQ field [22,[178][179][180][181][182][183][184].It would be interesting to see whether these models could still give rise to cosmological collider or primordial classical clock signals in the parameter space where the QCD axion isocurvature problem is solved.
• On the observational side, the signals we point out in the paper require new templates in the data analyses.The most outstanding example is that in our first two models, the leading component in the isocurvature power spectrum could be an oscillatory clock, which may not be captured by the current Planck template.
• It would be informative to study the constraints and forecast the future experimental prospects of the oscillatory isocurvature power spectrum and mixed curvatureisocurvature power spectrum, as well as those of the scale-or shape-dependent oscillatory mixed curvature-isocurvature bispectra.
and both factors could be large compared to H.The oscillating part of ccc reads, when μ = µ/H, ν ≡ (m A /H) 2 − 1/4 ≈ m A /H > 1: In the equations above, we have u ≡ 1 + ∼ 1 when the bare mass m Σ is small.In the first line, all the factors are from simple power counting rules except for an extra power (m A /H) 3/2 , which comes from an explicit full loop integration with late-time expansion and is not covered by the power counting rules [69].We require μ ∼ ν and the exponential factor could not be very large.Otherwise, the computation may not hold as it enters the non-perturbative regime.
On the other hand, replacing the external leg from ϕ to ã will produce bispectra signals involving isocurvature modes, such as the icc type.In this case, both the constraint on β in Eq. (2.7), and the small aF F coupling will suppress the signal.In particular, one needs to replace vertexes in Eq. (A.6) for each ã external leg attached: In addition, when converting δã to the DM isocurvature perturbation S d , one needs to replace the curvature normalization by Note that since one Higgs-gauge vertex σA 2 is replaced by aF F , the extra power dependence (m A /H) 3/2 from the loop integration may change as well.We will not explore it further here.Assuming that the power dependence doesn't change, we estimate Since β ≲ 0.04 and m A /H ≲ 10, we expect the icc type non-Gaussianity to be suppressed by O(10 −2 ), compared to ccc.
The same gauge boson loop will also introduce two-point mixing between adiabatic and isocurvature modes.Equivalently, the correlation between A s ad A i , parameterized by

Figure 2 .
Figure 2. The scenario of a sharp feature given Eq.(4.12).Left: φ1 /H 2 as a function of t; right: induced radial mode oscillations |σ 1 /f I | as a funtion of t.

Figure 3 .
Figure 3. Diagram of the sharp feature signal in the curvature power spectrum.The external points of inflaton are marked with blue squares.The solid blue line represents the propagator for ϕ while the empty dot denotes the interaction in Eq. (4.18).

Figure 4 .
Figure 4. Diagrams of the clock signals in the curvature (left) and isocurvature (right) power spectra.The external points of the axion (inflaton) are marked with red (blue) squares.The solid red (blue) line represents the propagator for a (ϕ) while the shaded blob denotes the interaction in either the second or the third line of Eq. (4.17), with an insertion of the oscillating σ bkg .
P i clock is equal to the oscillating part of |σ 1 /f I | times √ 2πµ σ .In other words, |σ 1 /f I | is suppressed by a factor of √ 2πµ σ ∼ O(10) compared to ∆P i P i clock .The two corrections, (4.20) and (4.21), are

Figure 6 .
Figure 6.Illustration of the correlated clock signals in the power spectra in model 1: ∆P i /P i (red dashed) and ∆P ζ /P ζ (blue solid).Both axes are in linear scales.The decaying envelope is due to (k/k r ) −3/2 as shown in Eq. (4.22) and Eq.(4.23).

Figure 7 .
Figure 7.The curvature bispectrum induced by the sharp feature.The external points of inflaton are marked with blue squares.The solid blue line represents the propagator for ϕ.The empty dot denotes the interaction in Eq. (4.25).

Figure 8 .
Figure 8. Diagram with an insertion of φ1 (represented by an empty dot) to the cii bispectrum in model 1.The filled red dot represents a constant vertex.

Figure 9 .
Figure 9. Diagrams with an insertion of σ bkg (represented by a shaded blob) to the cii bispectrum in model 1.The filled dots (red or black) indicate that the vertices are constants.

Figure 10 .
Figure 10.Sizes of the cii NG for two different types of diagrams in the equilateral limit.We set µ σ = 5, c = 0.1, Λ = 3 φ0 , and β = 5 × 10 −3 .In both panels, blue solid curves stand for diagrams where the feature insertions happen first, while red solid curves are the contributions of alternative diagrams where the features are injected at a later time.The blue dashed curves show the analytical approximation of the corresponding solid curves.Top: the diagram in Fig. 8 with φ1 injection.Bottom: the diagram with the σ bkg injection at the three-point vertex, as shown in the right panel of Fig. 9.

Figure 12 .
Figure 12.The leading diagram contributing to cii.The external point of inflaton is marked with blue squares while those of axions are marked with red ones.The solid blue, red wavy and dashed black lines represent the propagator for φ, a, and σ respectively.The empty dot denotes the interaction proportional to x∂ 2 ϕ V ϕ in the first line of Eq. (5.13).The filled red dot denotes the constant vertex between δa and δσ.Mixed isocurvature-curvature bispectrum.The leading diagram contributing to cii is shown in Fig.12.It contains the δσ(δa) 2 /f I vertex in the second line of Eq. (5.16) and one insertion of the −x∂ 2 ϕ V ϕ δ φδσ vertex in the first line of Eq.(5.13).This leading contribution from this diagram is

2 Figure 13 .
Figure13.The fermion loop diagrams for iii, iic, icc and ccc correlators.The red (blue) squares indicate the external axion (inflaton) mode.The red (blue) dots represent the coupling of fermion to axion (inflaton) in Eq. (6.10).The red wavy, blue, and black lines represent the axion, inflaton, and fermion propagators, respectively.

Figure 14 .
Figure 14.The f NL of pure isocurvature (left) and curvature (right) CC signals induced by the fermion loop.In both plot, we have the benchmark total fermion number N c N ψ = 3 and β ≃ A i /A s = 0.02 to satisfy the constraint in Eq. (2.7).

Table 1 .
Three simple PQ scalar models that lead to novel inflationary signals which involves isocurvature fluctuations.The signals in purple font indicate the signals enhanced due to sharp features in the inflationary potential, while those in orange indicate the signals enhanced due to chemical-potential type couplings.