Primordial gravitational waves from excited states

We show that a scalar excited state with large occupation numbers during inflation leads to an enhancement of tensor modes and a characteristic pattern of order-one oscillations in the associated stochastic gravitational wave background (SGWB) sourced during inflation. An effective excited state, i.e. a departure from the Bunch-Davies vacuum, can emerge dynamically as the result of a transient non-adiabatic evolution, e.g. a sharp feature along the inflationary history. We provide an explicit example in a multifield context where the sharp feature triggering the excited state is identified with a strong turn in the inflationary trajectory. En passant, we derive a universal expression for the tensor power spectrum sourced at second order by an arbitrary number of scalar degrees of freedom during inflation, crucially taking into account the nontrivial structure of the Hilbert space in multifield setups. The SGWB sourced during inflation can overcome the standard scalar-induced SGWB sourced at horizon re-entry of the fluctuations after inflation, while being less constrained by perturbativity and backreaction bounds. In addition, one may entertain the possibility of detecting both since they peak at different frequencies exhibiting oscillations with distinct periods.


Introduction
Since LIGO's "first light" [1], an astounding picture has been rapidly unfolding: our universe is filled with gravitational waves (GWs) sourced by black hole mergers. In fact, black holes are just one among numerous sources expected to contribute to the so-called stochastic gravitational wave background (SGWB); a bath of gravitational waves with frequencies spanning at least 20 orders of magnitude. Another important source is cosmic inflation [2][3][4][5][6], the period of exponential expansion that preceded the hot Big-Bang phase of our universe, responsible for the cosmic microwave background (CMB) anisotropies and the large-scale structure (LSS). Indeed, one of inflation's most remarkable predictions is the existence of a background of primordial GWs [7][8][9][10] with wavelengths ranging from planetary scales up to 10 4 Mpc (the largest observable scale). This corresponds to a range of frequencies between 10 −17 and 10 2 Hz. Current CMB observatories are indirectly searching for primordial GWs with frequencies 10 −17 −10 −13 Hz (in the form of B-mode polarisation), whereas the 10 −9 − 10 −1 Hz band will be accessible to future surveys such as LISA [11], SKA [12] and IPTA [13].
Within the inflationary paradigm, these anisotropies are traced back to quantum fluctuations around a quasi de Sitter background. In single-field, slow-roll inflation, the latter are predicted to evolve adiabatically, leaving the universe filled with inhomogeneities characterised by almost scale-invariant power spectra. This can be understood as the result of the smooth evolution of the background, disfavouring any particular time-slice as wavelengths are continuously stretched from sub-to super-horizon scales due to the exponential expansion of space. Scale invariance (i.e. spectra with constant amplitudes) implies that a detection of B-modes in the CMB polarisation would allow us to deduce the entire GW spectrum from its low-frequency amplitude. To date, CMB observations [14] have constrained the power spectrum of primordial GWs to be far below the signal sensitivity available to the next generation of detectors (one order of magnitude for SKA, three orders of magnitude for LISA) leaving essentially no possibility for a direct detection of the SGWB's primordial component.
However, the simplest models of inflation might not constitute an accurate description of the origin of primordial perturbations throughout the entire band of observable frequencies. CMB surveys cover but a narrow window of scales, rendering the extrapolation of the observed scale invariance of scalar perturbations to higher frequencies a potentially premature assumption. Actually, from a theoretical perspective, given the ultraviolet sensitivity of inflation and the associated difficulties to realise a prolonged stage thereof (see e.g. [15]), it is natural to entertain the possibility that inflation has occurred in successive periods with possibly vastly different properties. As a matter of fact, there exists circumstantial evidence pointing towards this possibility [16][17][18]: the black holes observed by LIGO/Virgo might be primordial black holes [19,20] resulting from the collapse of extreme over-densities on small scales. We thus have phenomenological, theoretical and observational motivations to consider non-standard realisations of inflation in order to describe the origin of cosmological perturbations, including primordial gravitational waves.
The purpose of this article is to uncover the characteristic frequency profile imprinted  Figure 1: In red, the characteristic inflationary gravitational wave energy density from the emergence of an excited state during inflation. This spectrum is the subject of the current work where it is computed explicitly. In green, the scalar induced post-inflationary gravitational wave spectrum resulting from the presence of the same excited state triggered by a sharp feature during inflation. The shape of this contribution has been computed in [21]. The grey dashed line represents the analytical template computed in Sec. 4, and explicitly given in Eq. (6.1), which is independent of the specific model generating the excited state.
in the inflationary-era SGWB spectrum stemming from an excited state that is dynamically generated during inflation at wavelengths well below those relevant for the CMB and LSS. Even if one considers standard Bunch-Davies initial conditions in the remote past, an effective excited state can appear as a general outcome of a momentary departure from the adiabatic evolution of the scalar fluctuations, which isolates a preferred set of scales during inflation, leading to (possibly) sizable, scale-dependent enhancements of the power spectrum of scalar perturbations. As we will show, the latter can cause a significant production of primordial GWs at localised frequencies, both during and after inflation. With the appropriate conditions, excited states during inflation can produce observable signals accessible for instance to LISA or SKA, turning them, along with the next-generation of CMB observatories, into main contributors to the nascent field of multi-messenger primordial cosmology [22][23][24][25][26][27]. Excited states can be generated in many scenarios of the primordial universe, a nonexhaustive (and non mutually exclusive) list of phenomena and references including: brief departures from slow-roll in single-field inflation, caused e.g by a step in the inflaton po-tential or a time-dependent effective parameter [28][29][30][31][32][33][34][35][36][37][38]; particle production due to timedependent masses and resonance phenomena [39][40][41][42]; sharp turns in the multifield landscape [21,[43][44][45][46][47]; multiple-stage inflation [48][49][50][51][52][53], etc. Note that in all of these examples, the excited state emerges dynamically at some moment during the evolution, which is crucially different from initialising the system in a non Bunch-Davies state (as in e.g. [54]). Independently of any particular microphysics, which plays an interesting but secondary role in what follows, the formalism and the results presented in this work have a broad range of applicability; they are derived under the single assumption of an excited state with large occupation numbers. We have kept the formalism as general as possible, taking into account the possible sourcing of GWs by multiple sources, which, in general, can be non-trivially quantum-mechanically correlated. Upon doing so, we will see that there is an interesting effect special to multifield dynamics, which amounts to a significant enhancement of the spectrum, even in two-field inflation, as a consequence of the quantum mixing among the fields.
As we will show, the presence of an excited state of scalar perturbations leads to a significant enhancement of the tensor modes generated during inflation compared to a sourcing by fields in their vacuum state. The intuitive physical reason is clear: an excited state means the presence of particles during inflation; these particles carry some energy beyond the minimal one of vacuum quantum fluctuations, and this energy -more accurately the transverse traceless part of the corresponding energy-momentum tensorsources GWs. Something worth highlighting is that this component of the SGWB is governed by the particle content and dynamics of inflation. It thus offers a probe of the physical processes during inflation that is complementary to the scalar-induced SGWB generated after inflation, which is sourced as the wavelengths of the primordial fluctuations re-enter the horizon after inflation and hence is sensitive only to the statistics of the fluctuations at the end of inflation (see [55] for a recent review).
Some of us have recently studied the scalar-induced GW background generated after inflation. The corresponding primordial curvature power spectrum displays large oscillations on small scales, characteristic of sharp features during inflation leading to substantial particle production [21]. The present paper thus complements this study by characterising the unavoidable inflationary-era SGWB that is also generated in this context. As we shall see, these two complementary signals of GWs appear in the SGWB in the form of two characteristic bumps in Ω GW (the density of primordial GWs), located at different frequencies -see Fig. 1. 1 This can be understood as the result of the emergence of two distinct comoving scales marking the positions of these maxima: the horizon size at the time of particle production, and the wavelength of the maximally enhanced mode of the primordial scalar power spectrum. Moreover, these GW spectra are modulated by oscillations with (distinct) frequencies, determined by the aforementioned scales, and whose relative amplitudes are different in the two cases: oscillations in the inflationary-era SGWB are of order one, while they are of order 10% for the ones generated after inflation. The resulting frequency profile of the total SGWB thus displays a very rich structure, offering a smoking gun signature of nontrivial inflationary dynamics far away from the CMB window.

Layout
We begin in Sec. 2 with a review of generic aspects of the SGWB. In Sec. 3, we derive a universal expression for the tensor power spectrum sourced at second order by an arbitrary number of scalar degrees of freedom during inflation, highlighting previously overlooked effects due to the quantum mixing thereof. We analyse the detailed structure of the spectrum in Sec. 4 providing analytical estimates of all the spectral characteristics. In Sec. 5, we discuss the interplay between theoretical constraints related to backreaction and perturbative control and prospects for detection in future GWs observatories, while Sec. 6 contains our conclusions. Finally, the (approximate) equivalence of the spectra computed with the retarded Green's function and with the in-in formalism is shown in App. A, further elaborating on points discussed in Sec. 3; appendix B contains details relevant for the analytic results obtained in Sec. 4.

Stochastic gravitational wave background basics
Let us start by defining the density parameter of gravitational waves, Ω GW , which will constitute the main observable quantity to be dealt with throughout this work. We are interested in studying the generation and evolution of tensor perturbations in a Friedmann-Lemaître-Robertson-Walker (FLRW) background. This may be done by perturbing the FLRW metric as where τ denotes conformal time. Here, h ij (x, τ ) is the tensor perturbation which is traceless (δ ij h ij = 0) and transverse (∂ i h ij = 0), while ζ is the curvature perturbation in this (comoving) gauge. It will be convenient to express h ij (x, τ ) in terms of Fourier modes h ij (k, τ ): The transverse and traceless conditions then translate to k ih ij = δ ijh ij = 0, allowing us to further decomposeh ij (k, τ ) in terms of spin-2 polarisation tensors as: with the two polarisation tensors e + ij (k) and e × ij (k) taken to be real. 2 Consistent with the fact that the source of GWs considered in this work comprises scalar degrees of freedom, in what follows, we consider a statistically homogeneous, isotropic 2 We use e + ij (k) = 1 √ 2 (mimj −ninj) and e × ij (k) = 1 √ 2 (minj +nimj), where {mi,ni} are the unit vectors orthogonal to k. In this way, e r ij (k) = e r ij (−k) and e r ij e s ij = δ rs . and unpolarised SGWB, and we write the two-point function of h r k as where P t (k, τ ) is the (total) dimensionless power spectrum, and where the brackets denote the statistical average, equivalent in our context to a quantum expectation value. When tensor modes are well inside the horizon and truly behave as free gravitational waves, with (h r k ) = kh r k , one can define their energy density as [57,58] where we have used units such that M 2 Pl = (8πG) −1 , while ≡ d/dτ . This expression allows one to obtain the density parameter as with Ω GW (k, τ ) the density parameter per comoving logarithmic scale. This expression gives an explicit relation between Ω GW and the power spectrum evaluated at a time where the modes of interest are sub-horizon. 3 In order to determine P t (k, τ ), we will consider scenarios where scalar fluctuations are enhanced during inflation, making them a relevant source of tensor perturbations h ij . To proceed, we must consider Einstein's equations linearised with respect to h ij in an FLRW background. This gives the following equation of motion for the Fourier modeh ij : whereS ij (k, τ ) represents the Fourier mode of the transverse and traceless component of the energy-momentum tensor acting as a source term. Just as we did with the tensor perturbation, we may decompose this source term asS ij (k, τ ) = r=+,× e r ij (k)S r k (τ ), giving an equation of motion for the polarisation modes h r k : where H = a /a is the co-moving Hubble parameter. The source S λ k (τ ) depends on the matter content of the universe at different stages of its evolution. During inflation, S λ k (τ ) receives contributions from the curvature perturbation ζ and (if present) other degrees of freedom. This gives rise to an inflationary component of the tensor field that we label h inf ij . On the other hand, after inflation S λ k (τ ) receives contributions from the super-horizon fluctuations ζ, produced during inflation, as they reenter the horizon [59][60][61][62]. The result of this second source is a post-inflationary component of the tensor field that we label h rad ij .
Thus, h ij can be written as the sum of two terms: 4 Inserting this decomposition into Eqs. (2.4)-(2.6) we can schematically (and with obvious notation) write Ω GW as follows: As we will see in a specific example, for excited states, the size of the two contributions Ω rad GW and Ω inf GW depends on different quantities and for most of the parameter space one contribution can easily overcome the other. For this reason, we do not consider the mixed term here. Ω rad GW induced by sharp features creating excited states at a given time along the inflationary history has been computed in [21] and the main result will be briefly reviewed in Sec. 4.2. The purpose of the current work is then to compute, in all generality, the pattern arising in Ω inf GW , i.e. the contribution from tensor modes sourced during inflation, due to the presence of excited scalar fluctuations. Usually the inflation-generated contribution is sub-dominant compared to the one sourced during the post-inflationary era. The reason is that, as we will see in more detail later, while in the latter case the source term in Eq. (2.8) is schematically given by S k ∝ ζ 2 ∝ P ζ , in the former case, considering the adiabatic perturbation Q ζ ∝ √ ζ, it reads S k ∝ Q 2 ζ ∝ P ζ ; hence, naively Ω inf GW ∼ 2 Ω rad GW . However, as we will show, the situation can drastically change if the temporal behaviour of the scalar modes becomes non-standard before horizon crossing and/or additional entropic degrees of freedom become relevant in the source term (see also [63]).

Stochastic background today
We are interested in the stochastic background of GWs as measured in the present cosmic era (τ 0 ). For a given mode in Fourier space, the frequency of GWs today is given by With the scale crossing the Hubble radius at matter-radiation equality being k eq 1.3 × 10 −2 Mpc −1 , all modes with frequencies f 10 −17 Hz have re-entered the horizon during radiation domination, unless a non-standard thermal history is considered between the end of inflation and the radiation-dominated era.
To compute the post-inflationary induced GWs today, it is sufficient to note that the source term is active when the corresponding mode re-enters the horizon and soon decays (as ∝ τ −2 during radiation) leaving a free propagating GW with an energy density ρ GW ∝ 1/a 4 . Thus, if we consider Ω rad GW at τ p , a time after horizon crossing for a given mode such that the source has become negligible, we have Deep inside the radiation era ρ cr (τ p ) ρ r (τ p ) and where g S and g * are respectively the effective number of entropic and relativistic degrees of freedom as a function of the temperature T . Thus, Eq. (2.12) becomes Ω rad GW (k, τ 0 ) = c g Ω r,0 Ω rad GW (k, τ p ), (2.14) with Ω r,0 the energy density fraction in radiation today.
The contribution to the stochastic background sourced during inflation is more conveniently expressed in terms of the tensor power spectrum at the end of inflation. In fact, P inf t at a given time after inflation can be written by means of a transfer function that takes into account the evolution of the tensor modes throughout the cosmic history, i.e., P inf t (k, τ ) = T 2 (k, τ )P t (k, τ end ). By inserting this expression into Eq. (2.6) (with the decomposition (2.10) in mind), we obtain [58,64] 5 For consistency, we have multiplied the transfer function one can find in [58] (obtained under the assumption that ρ r has always red-shifted as 1/a 4 ) by c g ; the factor that takes into account the different number of relativistic degrees of freedom when the modes of interest re-enter the horizon.
The first Standard Model degree of freedom to become non-relativistic is the the top quark that annihilates at about T m t /6 30 GeV. By recalling that the frequency of a GW produced at horizon crossing (during radiation) can be directly related to the temperature of the universe at that time 6 , one deduces that f (T 30 GeV) 8 · 10 −7 Hz.
In the review [58], an analytical result for T 2 (k, τ0) has been computed by considering modes that enter the horizon during radiation and matching their behaviour at a given time τ * (the time when the pure radiation and matter solutions for the scale factor cross) related in a precise way to the time of matterradiation equality. This analytical result is in good agreement with the full numerical one and has no substantial difference with respect to the interpolation in Eq. (2.15). We thus use the analytical result obtained in this way to estimate the prefactor in Eq. (2.19). 6 By using the conservation of entropy gS(T )T 3 a 3 = const and ρ rad = π 2 30 g * (T )T 4 we can rewrite Eq. (2.11) as f (Tp) = 2.5 · 10 −8 Hz g * (Tp) 100 where we used that for T Mev, gS g * .
Thus, if one is interested in frequencies relevant for GWs observatories like, for instance, LISA (10 −5 Hz f 10 −1 Hz) and LVK (1 Hz f 10 4 Hz) 7 , all Standard Model degrees of freedom can be safely treated as relativistic at the time of production. Therefore, g * (T p ) g S (T p ) 106.75, which, together with the present-era values g S (T 0 ) 3.91 and g * (T 0 ) = 2, leads to c g 0.4. Furthermore, as per common practice, in order to avoid propagation of uncertainties on the measurements of the Hubble parameter, we will consider the quantity h 2 Ω GW with H 0 = h · 100 Km/(s · Mpc).
Summarising, it is convenient to write the two contributions to h 2 Ω GW as where the explicit expression for Ω rad GW (k, τ p ) will be specified later in Eq. (4.17) (we will re-consider the "red-shifting" factors r i , r r when discussing observability in Sec. 5). The following sections will focus on computing P t (k, τ end ) ≡ P t (k) in presence of an excited state. For notational convenience we will also omit the τ 0 argument and write Ω GW (k) ≡ Ω GW (k, τ 0 ).

Multifield quantisation
We will be interested in GWs sourced by N scalar fluctuations to second order in perturbation theory. To this end, let us denote the corresponding gauge-invariant quantum operators (at linear order) byQ X (k, τ ) (more on their normalisation below), where X runs from 1 to N , and expand them in a basis of canonical creation/annihilation operators aŝ A crucial aspect of a multi-species system is that the ladder operator basis consists of N vectors labeled by the indices i, j (see e.g. [66][67][68][69]). This is the result of properly taking into account the interaction of N scalar degrees of freedom during inflation: the Hilbert space of the system is the tensor product of the individual Hilbert spaces, which is itself isomorphic to the N complex-dimensional vector-space spanned by ladder operators associated with one degree of freedom each. As usual, creation/annihilation operators and their corresponding vacuum state are defined by specifying a set of mode functions. Deep enough on sub-Hubble scales, one can always consider suitably defined fluctuationsQ X (k, τ ) that are decoupled and hence quantum-mechanically independent, and throughout this paper, we impose thereon the Bunch-Davies (BD) vacuum, i.e. the mode functions are chosen so that when |kτ | 1 their corresponding vacuum state is the Minkowski one (the unique one minimising the Hamiltonian at early times). This is equivalent to imposing the following initial conditions: with δ Xi the Kronecker delta, i.e. one can always choose the N elements of the ladder basis to be "aligned" with the N initially independent fluctuating degrees of freedom (see e.g. Sec 3.2 of [70]). In addition, Eq. (3.3) assumes that each field has been canonically normalised, i.e. S = 1 2 dtd 3 x XQ 2 X + . . . . Then, a specific model determines a system of N coupled differential equations satisfied by the operatorsQ X . In particular, the associated N 2 mode functions Q Xi in (3.1) correspond, for each field index X, to N independent solutions (labelled by the index i) of the N -field system at hand. To explicitly find the mode functions Q Xi , one simply solves the corresponding equations of motion with the N different sets of initial conditions; each one given by fixing the index i in Eq. (3.3) and let X vary.
In Sec. 4, we will focus on GW sourced by excited states, which imply a specific form of Q Xi (k, τ ). For the moment though, let us keep the discussion as generic as possible and derive the second-order scalar induced inflationary tensor power spectrum in a general multifield context.
Due to the SVT decomposition, tensors can only be sourced (to lowest order in perturbation theory) by the transverse, traceless component of the energy-momentum tensor furnished by the kinetic terms of whatever scalar fields comprise the model at hand. Therefore, in the notation of Eq. (3.1), the scalar source in Eq. (2.8) will be given bŷ where e +,× (k, p) ≡ p i p j e +,× ij (k) = p 2 √ 2 sin 2 θ cos(2φ), sin(2φ) , (3.5) with θ the angle between the wavevectors k (of the induced GWs) and p (of the source), while φ is the azimuthal angle of p. Note again that Eq. (3.4) is not restricted to models with canonical scalar fields; for instance, given a non-trivial field-space metric one can always diagonalise the kinetic term of the fluctuations by projecting them on a set of vielbeins (see for instance Sec. 4.7). Although our general formalism is independent of this, in such scenarios, it may be convenient to choose these vielbeins such that one of the fields Q X corresponds to the instantaneous adiabatic fluctuation, which we will do in concrete examples. That is, we will distinguish the first of the N fluctuating degrees of freedom by identifying it with Q ζ ≡ M Pl √ 2 ζ, with the other fields corresponding to instantaneous entropic fluctuations ψ 1 , ψ 2 , . . . ψ N −1 . Equivalently, one may use interchangeably the notations X = {1, 2, . . . , N } or X = {ζ, ψ 1 , . . . , ψ N −1 }.
Viewed as an operatorial statement, the formal solution of the tensor equation of motion (2.8) can then be expressed aŝ The Green's function can be written in terms of the vacuum mode functions. In de Sitter (dS) this becomes 8 with ζ(kτ ) = e −ikτ (1 + ikτ ), (3.8) the standard dS mode functions with Bunch-Davies asymptotics (where for convenience, we have not included the H/ √ 2k 3 factor), which are the same for tensor and scalar modes (hence the use of ζ). Furthermore, the vacuum contribution (homogeneous solution)ĥ λ k;0 (τ ) is uncorrelated to the source term (it comes with its own quanta): Before computing the power spectrum, let us comment on the use of the Green's function in this context, a subject that has been discussed in [71][72][73][74][75][76].

Field evolution in the interaction-picture
Equation (3.6) can be viewed as the first term in the expansion of the field operator in a series over the interaction-picture free fields: where T denotes time-andT anti time-ordering, while ∞ ± ≡ ∞(1 ± i ) accounts for the contour deformation in the infinite past. In this language,ĥ λ k;0 (τ ) is the interaction-picture field, whileĤ int is the interaction picture Hamiltonian: H int (τ ) = 1 (2π) 3 dp a 2 (τ )ĥ p;0 (τ )Ŝ λ −p (τ ), (3.11) withŜ given by Eq. (3.4).
To verify this, we may expand the exponentials to first order in H int , to obtain 12) 8 More precisely, the Green's function associated with the canonical (Mukhanov-Sasaki) variable is the solution of the corresponding homogeneous equation, satisfying the Wronskian condition v 0 which, upon inserting (3.11), readŝ Next, let us isolate the infinite past by inserting an arbitrary time 9τ [74]. After expanding the graviton in the canonical basis, with ζ the dS mode function (3.8) andâ,â † satisfying the commutation relations (3.2), we may use (3.7) to finally obtain The first line of this equation coincides with Eq. (3.6) (with the part of the integral between τ and τ ; for i → 0, the matching of the leading terms is exact). We may thus draw two conclusions: i ) the Green's function solution for the tensor field (3.6) is an approximation of the nonlinear result (3.10) [72], and ii ) it also differs from the latter as far as the implementation of the i prescription is concerned [74,75]. We further elaborate on this discussion in App. A. Both characteristics can be thought of as manifestations of the quantum nature of the inflationary tensor modes since both the contour deformation and the nonlinearity arise from a quantum-mechanical treatment which is appropriate at τ = −∞, where the BD initial condition is imposed. However, since we are considering an effective excited state emerging at some later time, the "quantum" character here translates into a statistical property of the random variable h: all the terms in the expansion (3.10) express the non-Gaussian variable 10 as a series over the Gaussian random fields h 0 , Q 0 much like the familiar local ansatz [77] for the scalar fluctuation (see also App. A.1 of [78] for a related discussion).

The tensor power spectrum
Let us for a moment (see the end of the section) neglect the last line of Eq. (3.15). Then the graviton two-point function can be simply written as where we have dropped the vacuum contributionĥ k;0 , since the sourced GWs will be the dominant component in the scenarios under consideration here. Note that the operators involved in Eq. (3.16) are the interaction-picture fields, which, here, have Gaussian eigenvalues, allowing us to proceed via Wick's theorem. Scalar non-Gaussianity boosted by the excited state, see e.g. [79][80][81][82][83][84][85], will enter at two loops and beyond via insertions of (at least) the ever-present cubic vertices [86]Ĥ A back of the envelope estimation implies that the perturbativity/backreaction conditions discussed in Sec. 5 should automatically grant radiative stability.
Performing the Wick contractions and ignoring the irrelevant disconnected contribution, one thus obtains where the scalar power spectra are given by with the mode functions defined in Eq. (3.1). Next, we may substitute the polarisation vectors (3.5) noticing that e +,× (k, p) = e +,× (−k, −p) = e +,× (−k, p − k). The integrals over φ then yield the polarisation Kronecker delta as πδ λµ . Writing the scalar source in Eq. (3.16) explicitly and using the definition (2.4), we arrive at the main formula for the total power spectrum of tensor modes sourced by scalar degrees of freedom during inflation: where we recall that g k is given in (3.7). We thus see that by considering a multisource scenario, not only does there appear a summation over the distinct scalar degrees of freedom sourcing GWs but, due to the mixing, also over the quanta comprising each source.
Notably, when all contributions are of the same order, this leads to an enhancement of the power spectrum proportional to N 4 . Finally, let us remind the reader that, as we show in App. A, this expression can be obtained via the in-in formalism at one loop. 11 Taken at face value, the integrations in Eq. (3.19) lead to divergences in the hardmomentum limit p → ∞ and the infinite past limit τ 1 → −∞. However, one has to keep in mind that our focus is on scenarios where only a short range of modes, starting from a given time, are enhanced. Thus, for arbitrary large values of the internal momenta (p, |k−p|), the mode functions follow a dS-like evolution and can be treated in the same way as in standard inflationary scenarios (without an excited state). The discussion regarding the exact finite part present in the literature (see [88] for instance) would not affect in anyway our results since the former is orders of magnitude suppressed compared to the contribution coming from the enhanced modes computed in the next section. Hence, for practical purposes, one can simply regularise the integral by introducing a finite cutoff in momentum space. 12 Regarding the time integral: in the next section we will argue that for the cases studied in this work, the timeτ that we used in Eq. (3.15) acquires a physical meaning (see also footnote 9). In the τ 1 <τ domain (which we call the "in region" in the next section), the mode functions follow again the standard dS evolution, rendering this contribution to P t subdominant compared to the one coming from τ 1 >τ . This preferred time thus serves here as another cutoff "shielding" the infinite past. As discussed in Sec. 3.2 and App. A, neglecting the τ 1 <τ contribution renders the power spectrum computed via the retarded Green's function, approximately equal to the one-loop, in-in power spectrum.
For later convenience, let us also define the following two sets of dimensionless momenta and Using these, the geometrical factor in (3.19) becomes

Stochastic gravitational wave background from excited states
We are interested in studying the stochastic background of gravitational waves sourced during inflation associated with the appearance of excited states at a given time along the inflationary history. Since interactions among multiple fields provide a natural playground for excited states to arise, we exemplify our claims in Sec. 4.7 with a two-field model in which a short period of strongly non-geodesic motion along the inflationary trajectory lies at the origin of the excited state. Let us, however, emphasise that all our main results regarding the shape of the signal follow solely from the presence of an excited state and are equally valid when the mechanism triggering it occurs in different (e.g. single-field) scenarios. As we will see, the precise origin (multifield/single-field) as well as the particularities of each model are encapsulated in the explicit functional form of the Bogoliubov coefficients, and play an interesting but secondary role. There exist though a characteristic that is specific to multisourced GWs: we have seen that our master formula (3.19) is a nontrival generalisation of the single source result owing not only to the summation over the various types of fluctuations but also over the various types of quanta. In "democratic" situations like the one studied in Sec. 4.7, this will in turn enhance the tensor modes by orders of magnitude (depending on the number of fields) due to a combinatorial factor.

Dynamically generated excited states
The main results of this paper only depend on the existence of a dynamically generated excited state, whose definition is simple: although all modes are initialised in the Bunch-Davies vacuum, a non-trivial inflationary dynamics is such that after some time, some sub-Hubble k-modes are not in their ground state anymore. A paradigmatic and physically motivated large class of models in which this mechanism may be at play is the one of sharp features, and for concreteness, we formulate things in this language in the following, although our results have a broad range of applicability. By a sharp feature, we mean a sudden change in some background parameter f (N ), where N henceforth denotes the number of e-folds. 13 This could be any background quantity controlling the dynamics of the perturbations during inflation. Our special focus is on sharp features that at the same time (exponentially) enhance the power spectrum for a limited range of scales.
In order to see that this naturally leads to an excited state with a large amount of particle production, we consider the evolution of the various degrees of freedom by splitting the time domain in three regions (see Fig. 2): in for N < N out − δ (with δ the duration of the feature), where the mode functions are placed in the BD vacuum in the infinite past and obey the standard dynamics on a slowly changing inflationary background that we approximate with a de Sitter epoch; f eature for the narrow region when the feature is active, i.e. between N out − δ < N < N out ; and out for N > N out , where the dynamics is back to standard (like in the in region) but now with different initial conditions set at N out by matching to the feature-region solution.
The types of sharp features we are interested in are such that the maximally enhanced modes are deep inside the Hubble radius at the time of the feature. This naturally leads one to distinguish two relevant scales: k * marking the maximum of the scalar power spectrum, and k out , the wavenumber of the mode that exits the horizon at N out . In order to quantify the hierarchy between them, let us introduce the parameter γ ≡ k * k out , (4.1) Figure 2: Time evolution of the curvature perturbation power spectrum in presence of a sharp feature for three different k-modes (normalised to P 0 , i.e. the single-field, slow-roll primordial power spectrum). In blue, an illustrative profile for the background parameter defining the feature. At the onset of the so-called out region, the relevant modes are still deep inside the Hubble radius. By that time, the sharp feature has prepared the system in an excited state, so that the corresponding power spectrum oscillates with a k-dependent phase before freezing outside the horizon. Scalar induced GWs generated after inflation are only sensitive to the snapshot of the power spectrum of ζ at late times, i.e. the primordial curvature power spectrum. Instead, GWs generated during inflation are sensitive not only to the various degrees of freedom entering in the energy-momentum tensor but also to their time evolution during inflation. The three curves have been computed explicitly by using the example of Sec. 4.7 with parameters (η ⊥ , δ) = (28, 0.25).
which will be useful when studying the enhancement of the GW energy density produced by excited states. Before discussing the behaviour of the mode functions appearing in the generic solution (3.1) in these three regions, let us re-write them as follows where H/ √ 2k 3 has been factored out for later convenience. In the in and out regions, solutions take respectively the form (up to a global phase factor) and is the standard de Sitter mode function. Note that we are considering scenarios in which the relevant enhanced modes are deep inside the Hubble radius at the onset of the out region, which, as we will see, is the most relevant time for GW production since gradients suppress the source at subsequent times. There, it is thus a good simplifying approximation to use massless mode functions ζ(kτ ). It is straightforward in principle to include mass effects, but technically cumbersome with the appearance of Hankel functions throughout that would obscure the simple physics we want to describe.
The dS mode function and its complex conjugate provide two independent solutions to the corresponding equation of motion, so that Q out Xi (k, τ ) is necessarily a linear combination of them, with coefficients α Xi (k) and β Xi (k) called Bogoliubov coefficients. The latter are not arbitrary though, as they should be such that at any time, theQ X (x) commute with one another, the same for their conjugate momentaπ Q X (x), and that Q hold. This imposes the following set of relations (with an implicit sum over the repeated index i): which are automatically satisfied by any unitary evolution from the in to the out region. In a single-field situation, they reduce to the well known relation |α| 2 − |β| 2 = 1, which generalises here, for any X, to but we stress that the whole set of relations (4.6), (4.7) should hold. 14 We will explicitly verify them for any time in the example of Sec. 4.7.
In the in region, the Bunch-Davies initial conditions (3.3), with α Xi (k) = δ Xi and β Xi (k) = 0 for all k modes, trivially satisfy (4.6), (4.7), with the mode functions in the τ → −∞ limit behaving as the Minkowski positive frequency modes. Then, the vacuum associated with the operatorsâ i ,â † i appearing in Eq. (3.1) corresponds to the lowest-energy state in the remote past, which we label as the in vacuum. The latter also represents the time-independent state of the system (since, as it is customary, the dynamics is described in the Heisenberg picture). In the out region, the dynamics is back to the standard one, with free fields propagating over a (quasi) dS background. Accordingly, the mode functions which behave as the positive frequency Minkowski modes in the remote past (here meaning, for each k mode, |kτ | 1 although still τ > τ out ) are analogous to the ones of the in region given in Eq. (4.3). Thus, one can write the fields operators in the out region aŝ whereb i ,b † i are a new set of creation/annihilation operators defining the vacuum (and the excited states) of the system in the out region. To find the relation between these new set of operators and the ones corresponding to the in vacuum, one simply equates the expansion (4.9) with the one in Eq. (3.1) after substituting the expressions (4.2)-(4.4) therein, giving:b (4.10) One can check that these operators satisfy the canonical commutation relations (3.2) by virtue of Eqs. (4.6), (4.7). The mean number of particles for each species at the onset of the out region is then given by the expectation value of the operatorsN X =b † Xb X in the state of the system i.e. the in vacuum), leading to which is a generalisation of the standard result n(k) = |β(k)| 2 (see also e.g. [89]). The occupation numbers n X (k) are exponentially greater than unity in the set-ups under investigation here. In fact, as we briefly review later, to have an enhancement of the power spectrum for a given range of modes due to the sharp feature, one needs |α(k)| 1 (removing indices for simplicity). This, together with the quantisation conditions (4.8), implies |β(k)| |α(k)| 1 for the relevant range of modes. In addition, the matching of the mode functions at the onset of the out region selects a k-dependent phase difference between the Bogoliubov coefficients (see for instance Sec. 2.3 of [21] for a more detailed explanation): β e 2ik/kout+iϕ α, (4.12) with k out the scale corresponding to the time of the feature N out , i.e. the mode k out exits the Hubble radius at N out , and ϕ is a phase factor whose k-dependence is mild compared to the rapidly oscillating first term; we assume α real without loss of generality. The exact time dependence of the mode functions in the region of the feature is modeldependent. However, since the features we are considering have the property to enhance the power spectrum of the scalar modes by several orders of magnitude during a short period of time, we parameterise this time dependence with an exponential enhancement, i.e.
, (4.13) with g(k) a model-dependent function of momentum.
Let us stress once more the sources of model dependence of the whole set-up under investigation, highlighting in parallel what is completely generic. Firstly, the precise functional form of the Bogoliubov coefficients in the out region depends on the exact mechanism at hand. However, by inserting the parametrisation (4.4) into Eq. (3.19) we derive generic formulae for the enhancement of the tensor modes that hold for arbitrary k-dependent coefficients. It is then reasonable, since one considers a narrowly peaked curvature power spectrum, to study the enhancement of tensor modes in the limit where these coefficients are peaked around a given scale k * . That will allow us to derive explicitly all the main features of the signal. However, as we stress in Sec. 4.5, most of them are indeed independent from the shape of the Bogoliubov coefficients. The second source of model dependence lies in the mode functions Q feat (k, τ ) in the region of the feature. As already stated, we parametrise the latter by means of Eq. (4.13). In any case, the contribution to the tensor power spectrum coming from this region is subdominant (see discussion in Sec. 4.6.2) and the main conclusion about the shape of Ω GW is robust against the exact form of this parametrisation.

Primordial power spectrum from a sharp feature
Since the GWs sourced in the radiation era depend directly on the scalar power spectrum P ζ , to facilitate the comparison with the inflationary-era sourced GWs, we close this section by briefly summarising the effect of a sharp feature on P ζ (see for instance [21] for more details).
For a sharp feature inducing an excited state, from the definition (3.18) and the solution (4.4) in the out region, we have that (4.14) where the relevant quantities are evaluated at horizon crossing and P 0 = H 2 /(8π 2 M 2 Pl ) is the single-field, slow-roll, dimensionless scalar power spectrum. As usual, dimensionless power spectra are defined as P XY = k 3 2π 2 P XY . Suppressing the indices on the Bogoliubov coefficients for simplicity, one can expand the square in Eq. (4.14) and write the power spectrum as where the appearance of the cosine is a direct consequence of the relation (4.12). For |α| ∼ |β| 1, we have that |β|/|α| 1 so that we can simplify this further to write That is, for a sharp feature leading to a significant enhancement of fluctuations the scalar power spectrum takes the form of an enhanced envelope P ζ modulated by sinusoidal oscillations in k with unit amplitude. The frequency of this oscillation is 2/k out , i.e. it is set by the scale of the feature.
For significant but not exponential particle production, i.e. for |α| ∼ |β| ∼ 1, the amplitude of oscillation is less than unity but still O(1). In that case, the power spectrum can still be boosted compared to its value at CMB scales as long as P 0 is larger at the relevant scales, e.g. if has a smaller value than when CMB modes cross the horizon.

Post-inflationary generated GWs from excited states: brief review
Excited states during inflation will also produce scalar-induced GWs in the post-inflationary era when the relevant fluctuations re-enter the horizon [59][60][61][62]. 15 The equation of motion for the tensor modes is again (2.8) but the source in this case depends to leading order on a four-point function of primordial curvature perturbations. If these are Gaussian, this can be written as a product of two instances of the scalar power spectrum, but in general there will also be a contribution proportional to the trispectrum, see e.g. [91][92][93][94][95]. At CMB scales primordial fluctuations are highly Gaussian and one expects the trispectrum contribution to the GW spectrum to be negligible compared to the term involving the power spectrum. However, for the modes affected by the excited state this is not necessarily the case and the trispectrum contribution may become important. We leave this for future work.
The present-era fraction of energy density in GWs that were sourced in the postinflationary era can then be written as [61,62] The factor c g Ω r,0 relates the GW energy density fraction at the time of its sourcing in the radiation-dominated era to that of today -c.f. Sec. 2.1. The integration kernel T rad (d, s) is computed in terms of time-integrals over the transfer functions and Green's function factors that relate the primordial curvature fluctuations to the source term in Eq. (2.8).
It also includes a kinematic factor that arises from rewriting the integral over momenta in terms of the variables (d, s) and an oscillation average has also been performed. The kernel depends on the equation of state of the universe when the relevant fluctuations re-enter the horizon. Assuming a standard thermal history of the universe, here we take this to be an era of radiation domination, in which case the relevant kernel is given in [97,98]. 17 One important property of T rad (d, s) is the existence of a singularity for s = √ 3 signalling resonant amplification. As a result of this, for a finite-width peak in P ζ at k = k * that is sufficiently narrow, the post-inflationary contribution to the GW spectrum will exhibit a narrow principal peak from resonant amplification at k 2k * / √ 3. In addition, there is generically also a lower broad "bump" around k ≈ k * / √ 3. For broader P ζ , the two peaks in Ω rad GW are increasingly "blurred" into one by the double convolution in (4.17) so that one eventually finds a single broad peak in the GW spectrum around k ∼ k * . 15 See also [90] for the impact of primordial dark matter isocurvature fluctuations on the scalar-induced SGWB generated after inflation. 16 Here we have rescaled the integration variables (d, s) compared to their namesakes in the previous works [21,96] by a subset of the authors, to be consistent with the definition of (d, s) in Eq. (3.21). 17 The corresponding expressions for a different equation of state can be found in [99][100][101][102] and the resulting post-inflationary GW spectrum due to an excited state during inflation was analysed in [103].
Let us now briefly review the effect of an excited state due to a sharp feature during inflation on the radiation-era GW spectrum [21]. We consider the case where the scalar fluctuations are enhanced by the feature, i.e. |α ζi | ∼ |β ζi | 1. As described in Sec. 4.1, this produces a peak in P ζ in virtue of |α ζi | 1 that is further modulated by O(1) oscillations due to |β ζi |/|α ζi | ∼ O(1). The width of the peak ∆k is given by the range of scales affected by the feature. One strategy for understanding the corresponding GW spectrum is to treat the oscillation in P ζ as a series of individual peaks. These can in general be taken as narrow peaks in the sense described in the paragraph above, as every such spike only covers a fraction of ∆k, which for a sharp feature already corresponds to a narrow interval in general. Every such peak in P ζ then produces its own resonance peak in Ω rad GW . There are also further resonance peaks from interactions between different spikes coming from the two factors of P ζ in Eq. (4.17) [21,104].
The sum of these contributions then gives rise to a spectral shape of Ω rad GW which for the scales of maximal enhancement is well-reproduced by the template [21]: The sinusoidal modulation arises from the superposition of the various resonance bumps, with the maxima of the cosine coinciding with the maxima thereof. Note that the frequency of oscillation is increased by a factor √ 3 compared to that of the oscillation in P ζ . As the resonance peaks have finite width, their multiple superposition has an effect of averaging out the modulation and suppressing the amplitude A lin . Even for the minimal case with just an O(1) number of oscillations within the interval ∆k, one finds that at best A lin ∼ O(20%), as can be seen from the green curve in Fig. 1. This decreases further as the frequency of the oscillation is increased. The smooth background spectrum Ω rad GW (k) can be shown to be given by the GW spectrum due to the smooth background of the scalar power spectrum P ζ (k) ∼ |α ζi (k)| 2 , which here is given by a peak of width ∆k with maximum at k = k * .
If this is narrow, i.e. ∆k/k * < 1, the GW spectrum Ω rad GW takes the usual form of a broad lower bump at k ∼ k * / √ 3 and a principal resonance peak at k 2k * / √ 3, see again the green curve in Fig. 1.

GW
Let us now turn back our attention to the GWs sourced during inflation. By inserting Eq. (4.2) into Eq. (3.19) with the change of variables (3.21) -or (3.22)-we can rewrite the tensor spectrum P t (k) ≡ lim kτ →0 P t (k, τ ) as with the geometrical factor given by (4.20) Furthermore, in particular where we have introduced the dimensionless variable z, defined as z ≡ kτ. (with δ accounting for the small but finite duration of the feature), the total result can be organised as follows: The constituent power spectra (let us henceforth suppress the integration limits), indexed by x = {out/feat/mix}, are given by with where Q are the mode functions (4.2) corresponding to each region.
Consistently with the conditions of applicability of our formalism, that we spelled out in (3.3), note that we disregard contributions from the in region. It is instructive to consider separately the contributions coming from the excited states alone, i.e. ∆ out t , the one that takes into account the finite time to "excite" these states, ∆ feat t , and consequently the mixed contribution.
In order to see why excited states naturally lead to an enhancement in the sourced gravitational-wave background, let us first focus on the contribution coming from the out region P out t . This is expected to provide the dominant term in the sum (4.25) since it is there that the excited states have support. This statement will be proved later while computing the contribution from the region of the feature in Sec. 4.6.2. As already mentioned, the model dependence of P out t only lies in the explicit functional form of the Bogoliubov coefficients. In this section we are thus able to derive generic formulas valid for any set of excited states. Furthermore, by considering a scalar power spectrum peaked around a given scale (of which a sharp feature is only a particular case) we provide explicit analytical approximations that capture the main characteristics of the signal in Sec. 4.6.

Integrating over internal time
Let us begin by noticing that the time integral in ∆ out t given in Eq. (4.27), can be performed analytically. Plugging Eq. (4.4) and the de Sitter Green's function (4.21) into Eq. (4.30), and relabeling the Bogoliubov coefficient associated with particle production as α + Xi ≡ α Xi and α − Xi ≡ β Xi , allows us to write P out t in a compact way as follows: (4.31) This generic formula for the tensor power spectrum in presence of excited states is one of the main results of this work. The time integral has been factored out as Note that from this definition it follows that G(x, y, z out ) = G(y, x, z out ), The integral (4.32) can be computed explicitly: 18 and For simplicity, let us focus only on one component of the spectrum corresponding to a single quantum degree of freedom. To that end, we will drop the indices from α ± Xi by implicitly fixing them to, say, X = ζ and i = 1, e.g. α + ζ1 ≡ α and α − ζ1 ≡ β. The outregion spectrum (4.31) will contain eight terms and their complex conjugates, which can be written as where we have defined ρ ≡ β α , (4.38) which keeps track of terms related to particle production induced by the excited state. Note that we have also omitted the third argument in G.

Enhancement of the tensor spectrum from excited scalar states
The prefactor outside the integral in Eq. (4.37) can be rewritten as H 4 /(8π 4 M 4 Pl ) = 8 2 P 2 0 , which, together with Eq. (4.16) for |α| 2 , confirms the expectation that the tensor power spectrum is proportional to 2 P 2 ζ . The sum in parenthesis then determines the enhancement of the tensor power spectrum with respect to this naive expectation. Note that there are six different combinations of G's weighing different components of the spectrum: (4.39) To better appreciate the peculiarity of excited states, let us write the time integral explicitly: (4.40) When the momenta of the scalar modes running in the loop are deep inside the horizon, i.e. when xz 1 and yz 1, the highly oscillating exponential damps the value of the integral for non-excited states where only the (+) contribution is present. In contrast, excited states add contributions with a (−) sign. These are terms in which the constructive interference between positive and negative frequency modes enhances the integral in the large argument regime, i.e. |G(x, −y)| |G(x, y)| for large x and y, and close to the resonant point x y where the frequency of the oscillating piece is small. This effect of constructive interferences is well known from studies of primordial non-Gaussianities in the presence of excited states, resulting in an enhancement of the bispectrum near flattened configurations [79][80][81][82][83][84][85]. In fact, the argument of the square in Eq. (4.31) contains the tree-level, tensor-scalar-scalar correlator (this could be made explicit via the cutting rules of [105][106][107] applied to the correlator itself). The similar enhancement of the latter due to excited states leads to the amplification of the tensor power spectrum that we discuss here.
Let us recall that, in general, to have a potentially observable Ω GW , a significant enhancement of P ζ with respect to the single-field, slow-roll result P 0 is required (assuming that the latter corresponds to the CMB pivot value). As reviewed at the end of Sec. 4.1, in our context, this automatically implies a large amount of particle production, i.e. |α| |β| 1, and |ρ| 1. Thus, let us now focus on this arguably most interesting case. Then, for |ρ| 1, and recalling that |G(x, −y)| generally dominates over |G(x, y)| over most of the integration range, it follows that P out t is dominated by the last line in (4.37): where θ ≡ Ph(β/α). Interestingly, the GW spectrum in this case exhibits many universal features independently of the shape in k of the Bogoliubov coefficients, as can be seen in Fig. 3, where we plot Ω inf GW (k) for various realisations of an excited state with |α| 2 |β| 2 1. The GW spectrum exhibits a maximum that can be checked to occur at k ∼ k out . Furthermore, there are oscillations on the UV tail with a relative amplitude that quickly approaches unity along the tail. The frequency of the oscillation is the same in all cases and is given by ω = 2/k out .
The appearance of the oscillations on the tail can be understood to arise from the properties of the kernel G(x, −y). From the explicit expressions in (4.34)-(4.36), note that the contribution F(x, −y, z out ) exhibits an overall prefactor e −izout that does not depend on the integration variables (x, y). It is this prefactor that is unaffected by the integral over (x, y) that is responsible for terms of type ∼ cos(2z out ) = cos(2k/k out ) in Ω inf GW (k), the factor of 2 coming from the fact that the integration kernel contains is quadratic in G(x, −y).
We will revisit this in more detail when specialising to peaked Bogoliubov coefficients below. The main point of this discussion here and of Fig. 3 is that many of the properties of Ω inf GW that we will describe in more detail for peaked Bogoliubov coefficients below, qualitatively also hold for broad profiles for |α(k)| 2 , like the example with σ = 1 in Fig. 3.
We now restrict attention to setups that produce a peak in P ζ (k) at k = k * , or equivalently in |α(k)| 2 , -as induced here by the excited state. In this case one can give an approximate expression for the integral in Eq.  1. Here we consider |α(k)| 2 = α 2 exp − 1 2σ 2 ln 2 (k/k * ) for different values of σ and we took Ph(β/α) = 0 for simplicity, and we also include |α| 2 and Ω inf GW from the example 4.7 with (η ⊥ , δ) = (28, 0.25) (red). The main finding is that the principal properties of Ω inf GW are independent of the precise shape of |α(k)| 2 , i.e. Ω inf GW exhibits a peak at k ∼ k out k * and a UV tail modulated by O(1) oscillations with frequency ω = 2/k out . Although the precise fall-off of the GW spectrum in the UV depends on the shape of the Bogoliubov coefficients, its behaviour is still universal under the (motivated) assumption that the latter are peaked around k * . In fact, for the two examples with the narrowest peaks in |α(k)| 2 one can check that the GW spectra become near-identical and only differ by a numerical factor. upshot is that the peaked Bogoliubov coefficients select the preferred locus kx = ky = k * , and thus we can use for any function F . We also used that from Eq. (4.16), the amplitude of the scalar power spectrum is given byP ζ 2P 0 |α| 2 . Omitting the bar from now on, we thus have where P ζ * ≡ P ζ (k * ) and κ = k/k * . To be complete, note that the approximation in (4.42) is valid when the contribution to the integrand apart from the |α| 2 -terms varies sufficiently slowly.
In the next section, starting from the approximation (4.42), we will be able to derive analytically the main features regarding the shape of the GW spectrum. We stress that we do not need to use the explicit relationship (4.12) between the phases of Bogoliubov coefficients induced by the presence of a sharp feature to draw our main conclusions. In contrast, this phase difference is at the root of the oscillations in Ω rad GW . We can also use the analytical form of G(x, y, z out ) in (4.34) to understand how the (+) and (−) contributions individually affect the tensor power spectrum. Let us focus on the locus x y and distinguish two cases: k k out and k k * meaning momenta of the GW close to the one of the feature and to the maximum of the scalar power spectrum, respectively. In both situations, the momenta of the sourcing scalar modes are xk k * so that xz out ≡ xkτ out k * τ out 1. Thus, for k k out (|z out | 1), x 1. In this case, K(x, x) ∝ 1/x 2 while F(x, x) ∝ x. In contrast, the same functions with a minus sign in one of the argument -signature of an excited state -scale as K(x, −x) ∝ x 2 and F(x, −x) ∝ e −izout x 2 . In the range k k * , x 1 and |z out | 1, the function K, which is independent of z out , is not relevant and we have F(x, x) ∝ e −i(1+O(1))zout z out and F(x, −x) ∝ e −izout z out . Here, the (+) and (−) contributions have similar size and the higher frequency (due to |z out | 1) makes the oscillations visible in the GW signal (see explicit details below). We summarise this discussion in Fig. 4 from which it is clear that, for a given strength of the scalar source, i.e. parameterised by P Q ∼ P ζ A, tensors are orders of magnitude enhanced in presence of an excited state.
Having discussed the case of large of particle production, one may entertain the possibility that the scalar power spectrum is enhanced via other mechanisms (such as a change in the Hubble flow) so that P 0 for the scales of interest is already much larger than the power spectrum at CMB scales. In this framework, a non-adiabatic transition generating excited states does not necessarily need to enhance the power spectrum by several orders of magnitude in order to achieve observationally relevant values of Ω GW . Thus, given the generic framework of this section, let us briefly comment on the case where a transient non-adiabatic evolution leads to a small amount of particle production, i.e. |β| |α| (or equivalently |ρ| 1). When β = 0, i.e. the modes are in a Bunch-Davies state, only the first term in (4.37) contributes and, for a fixed P ζ , the tensor power spectrum is damped with respect to the case |ρ| 1. With a small amount of particle production ρ 1 (due to the enhancement of scalar modes on sub-Hubble scales), if the enhancement proper to excited states is sufficiently large so that |G(−x, y)/G(x, y)| > ρ −1 , then the would-be small corrections linear in ρ in Eq. (4.37) might become dominant, resulting in different shapes of P t . In the limit x y, Eq. (4.37) becomes It is worth stressing that, in the limit of small ρ, the conclusions drawn from the previous formula about the shape of P t should be only taken as qualitative. In fact, in this limit, the concept of a preferred time τ out becomes ill-defined. That was the key assumption behind the use of the classical Green's function method, see Sec. 3 and App. A, and the possibility to neglect the infinite past contributions therein. Thus, the rigorous treatment of this case requires a proper full in-in "quantum" computation for small ρ that is beyond the scope of the current work.

Detailed structure of the spectrum
We now have all the necessary ingredients to extract the overall observable structure of the tensor power spectrum. We will derive analytically the position and the scaling of the peak, which occurs around k k out , as well as the frequency of the oscillations occurring at smaller scales around k k * .

Out-region contribution: maximum and oscillations
Let us first focus on the out region. By means of the approximation (4.42), Eq. (4.31) (in the limit of large particle production) can be easily integrated. The tensor power spectrum becomes where we have assumed that all fields and quanta contribute equally so that N is a place holder for the number of fields. The explicit expression of the function G appearing above simplifies considerably at the point imposed by the approximation (4.42), i.e. x y 1/κ, and we stress once more that the minus sign in the second argument of G is related to the negative frequency mode present only in the case of excited states. After some manipulations, it is useful to rewrite G in (4.45) as follows: where, as before, k * marks the maximum of the primordial scalar power spectrum P ζ , while γ = −k * τ out = k * /k out parameterises how deep inside the horizon was the maximally enhanced mode k * at the onset of the out region. Recall that we consider γ 1 so that the last term in the previous expression can be safely neglected for all practical purposes.
By inserting Eq. (4.46) into (4.45) and (2.19), we obtain an explicit analytical template for Ω inf GW that matches remarkably well the numerical results for the full spectrum (see Figs. 1, 6 and 8): Further, one may estimate the behaviour of the envelope by sending sin(γκ) → 1 and cos(γκ) → 1. This leads to the following simple expression (4.48) As we are going to show in a moment, the signal has its maximum for κ 1. Thus, the falloff of the spectrum right after its peak follows approximately the simple power-law behaviour κ −3 -see Fig. 8. Note that by momentum conservation, the wavenumbers of the tensor modes generated at second order by scalar perturbations peaked at the scale k * cannot exceed 2k * , i.e. κ ≤ 2, where the signal indeed vanishes. Furthermore, for small enough momenta, i.e. for γκ 1 on the left of the peak of the signal, Ω inf GW behaves like κ 3 , resulting in a symmetric envelope of the signal with respect to its peak, which is also well visible in Fig. 8.
As can be checked a posteriori, the maximum of the spectrum is to a good approximation given by the maximum of the function G, with the prefactors κ −1 µ(κ −1 , κ −1 ) in (4.45) only leading to small corrections suppressed by γ −1 . Thus, to find the maximum, we set to zero the derivative of (4.46) with respect to κ and we look for the solution with the smallest κ. In the limit γ 1 this leads to 4 − 3γκ sin γκ + (γ 2 κ 2 − 4) cos γκ = 0, (4.49) which depends only on γκ. Thus, one can find a universal first root γκ max = c, i.e.
k max = c k out , with c 3.505, (4.50) which is valid for any k * and γ 1. 19 It is now easy to estimate how P out t scales with γ at its maximum by inserting (4.50) in Eqs. (4.46) and (4.45), i.e. by considering the regime γ 1 with γκ max = c: where we used A 2 ∝ P 2 ζ /N 2 . By reintroducing the "redshift factors" defined in Eqs. (2.18)- Figure 5: Scaling of the peak in Ω inf GW as a function of γ. Coloured lines correspond to h 2 Ω inf GW evaluated at the maximum (4.50) (computed for the two-field example in Sec. 4.7 without the approximation of peaked Bogoliubov coefficients), for different enhancement of the primordial scalar power spectrum p ≡ 1/2 ln(P ζ /P 0 ) and fixed = 10 −2 . Each dashed line corresponds to 10 −2 N 2 2 γ 5 Ω rad GW evaluated at k = 2k * / √ 3, which is the location of the maximum of the scalar-induced GWs sourced during radiation. A given p translates into a fixed Ω rad GW maximum amplitude. The figure confirms very accurately the scaling obtained with the simple analytical estimate in Eq. (4.52). (2.19), and putting all together we arrive at relating the amplitude of the two peaks, i.e. the one in Ω inf GW and the one in Ω rad GW respectively, to the slow-roll parameter and γ. The factor 10 −2 comes from the ratio of the redshift factors r i /r r 4 · 10 −2 . Remarkably, despite the expected 2 suppression, the inflationary contribution is boosted by a γ 5 factor, which we confirm numerically in Fig. 5. There, the ratio between each solid line (Ω inf GW ) and the corresponding dashed line (N 2 2 10 −2 γ 5 Ω rad GW ) is indeed always an order one number (taking values between 1 and 2). The relation (4.52) underlines that, for a given P ζ , the amplitude of Ω inf GW can easily be comparable or larger than the one of Ω rad GW at their respective maxima. Further, since the different maxima are located at different scales, i.e. k inf max k out while k rad max k * = γk out , the two contributions might be of the same order and still be distinguishable providing a unique pattern for the total Ω GW -see for instance Fig. 1.
Let us now study the shape of the spectrum around k = k * . Here, we can approximate Eq. (4.46) by going to the limit γκ 1. The first term in (4.46) dominates and P out t oscillates with a constant frequency: Although in this region the signal is suppressed compared to its value at the peak k k out , it exhibits order one oscillations with a frequency ω in k-space given by which is the same frequency that also controls the modulations in the primordial scalar power spectrum, see Eq. (4.16), i.e. similar to scalar fluctuations, tensor modes enhanced during inflation select the same preferred scale k out determining the onset of the out region. As reviewed in Sec. 4.2, the same mechanism leads to oscillations in the post-inflationary scalar-induced GWs, but with a larger frequency, i.e. ω rad = √ 3ω, where the numerical factor stems from the fact that scalar perturbations only source tensor modes once they re-enter their sound horizon. 20 Note that the periodic peak-structure in Ω rad GW is due to a resonance mechanism that "averages" order-one oscillations in P ζ to ∼ 10% oscillations in Ω rad GW . In contrast, the oscillations in Ω inf GW from Eqs. (4.47) and (4.53) are genuinely of order one, offering even better prospects of detection in future GWs observatories.
GWs signals exhibiting oscillations with the same frequency have been shown to arise for instantaneous sources active during inflation in [108]. The mechanism presented in the current work is based on a dynamically emergent excited state, and not only provides an explicit realisation for the appearance of such oscillations but it also yields a richer built-in structure for the GWs spectrum; for instance, the enhancement of the peak due to the constructive interference between positive and negative frequency modes as described above, the particular shape of Ω inf GW as given by Eq. (4.45) and the possible oscillatory counterpart in the post-inflationary Ω rad GW . The same line of thought applies to [109,110] where visually similar shapes have been numerically computed for models with resonances. It is likely that our formalism and results are applicable there. Finally, let us stress once more that all the main results of this section hold generally for excited states sourcing GWs during inflation; these are the three equations in the squared boxes regarding the position of the maximum (4.50), the enhancement of Ω inf GW at its peak (4.52) and the frequency of the oscillations around k k * (4.54). In fact, these outcomes rely only on the functional shape of G that comes from the time integral over the three de Sitter mode functions, together with the assumption that the Bogoliubov coefficients are peaked around a given scale. The latter assumption is guaranteed when the event (the feature) creating the excited state is sharp.

Feature-region contribution
Until now we have presented results for the total tensor power spectrum P t as if the only relevant contribution was the one from the out region, i.e. P t P out t . Thus, let us end this section by showing explicitly that the contribution from the feature region, i.e. the one coming from the time dependence of the mode functions in the region of the feature, is indeed small. In particular, as we are going to show, the contribution is suppressed by the small parameter δ describing the finite duration (in e-folds) of the feature. Let us start by considering the term inside the modulus in Eq. (4.28), using Eq. (4.13): with where g ≡ g(xk) + g(yk) comes from the model-dependent exponent g(k) defined in Eq. (4.13), and we recall that G(0, z) is given in Eq. (4.22). The time-integral (4.56) can be performed analytically: which can simply be obtained by neglecting the variation of G(0, z) in the integrand in (4.56) compared to the exponentially growing factor. Note that to first order in δ the modeldependent parameter g drops out. Let us now also consider the prefactor in parenthesis in Eq. (4.55) and, for simplicity, focus again on a single mode (that is, disregarding the sum over fields by fixing the quanta indices to i = j = 1). Taking the modulus square of (4.55), we obtain where 2|α(p)| 2 1 + p 2 τ 2 out + (1 − p 2 τ 2 out ) cos ϕ + 2pτ out sin ϕ ; (4.60) in the second step we have used Eq. (4.12), i.e. β(k) = α(k)e −2ikτout+iϕ . Putting everything together, namely inserting the expression (4.60) into Eqs. (4.59), (4.28) and (4.26), and then using the narrow-peak approximation (4.42), we arrive at In the first expression we have used the definition γ = −k * τ out , together with γ 1. We have assumed that all mode functions in the multifield sum give the same contribution. This brings a factor N 4 multiplying the single-field result. Equation (4.62) is valid under the additional approximation on the integral I in Eq. (4.58). This overestimates the overall amplitude of P feat t -see Fig. 6, unless δ 10 −2 . 21 In any case, from Eq. (4.62) we learn that the contribution from the feature region is suppressed at least by a factor of δ 2 . Moreover, this expression allows us to understand the oscillatory pattern arising from this contribution. In fact, |G(0, −γκ = k/k out )| 2 provides oscillations with the same frequency as given in (4.54). Overall, from Eq. (4.62) we can infer how the contribution from the feature region scales with the various parameters: where again we used A 2 P 2 ζ /N 2 . To summarise, the feature-region contribution is enhanced by a factor γ 4 and suppressed by, at least, a factor δ 2 . In contrast, the out region -c.f. Eq. (4.51)-scales as γ 5 , and thus, in general, dominates the spectrum.

Explicit example: turning in field space
We provide here a concrete example in which an excited state emerges as the result of a sharp feature along the inflationary dynamics. The setup illustrated in this section relies on a model-independent, multifield mechanism first proposed in [44,45] to seed primordial black holes. All numerical plots in the current paper, e.g. Figs. 6-10, are produced using this mechanism as a benchmark, which serves as an explicit illustration for our results. In this example, the sharp feature triggering the excited state corresponds to a sudden and strong turn of the trajectory in field space. Let us start with the Lagrangian for the general class of non-linear sigma models whose target-space geometry is given by the metric G IJ : where φ = φ 1 , · · · , φ N and V (φ) is a generic multifield potential. Let us restrict for simplicity to two fields (N = 2). A convenient way to organise perturbation theory is the adiabatic-entropic basis defined by the unit vectors T I =φ I /(G JKφ JφK ) 1/2 (tangent to the background trajectory) and N I (orthogonal to the trajectory) with the pair T I , N I selecting a definite orientation. Here,˙≡ d/dt. Expressing the field fluctuations as δφ I = Q φ T I +Q ψ N I and fixing the comoving gauge, i.e. Q φ = 0 & g ij = a 2 e 2ζ δ ij , we may proceed to obtain the dynamics of the two degrees of freedom Q ψ , Q ζ ≡ M Pl √ 2 ζ. The effective action at second order in these linear fluctuations around a given homogeneous background reads 22 [67,111,112] where the last term describes the coupling between the two types of perturbations. Here, η ⊥ ≡ 1 H N I D t T I is a dimensionless parameter measuring the turning rate of the trajectory, with D t A I ≡Ȧ I +Γ I JKφ J A K the covariant time derivative along the background trajectory, which deviates from a geodesic in field-space when η ⊥ = 0. The resulting equations of motion readΠ where for simplicity we are considering H constant here and in what follows, and we have used Π ζ ≡Q ζ + 2η ⊥ HQ ψ . The mass parameter m 2 s turns out to be determined by background quantities as with V ;ss = e I s e J s V ;IJ the projection of the covariant Hessian of the potential along the entropic direction, and R fs the field-space scalar curvature. An alternative and also useful notion of mass is that of the entropy mass given by It can be shown that µ corresponds to the rest-energy of the massive particle state of the spectrum on sub-Hubble scales [113,114]. In addition, on super-Hubble scales, one can integrate Eq. (4.66) once. Upon doing so, Eq. (4.67) reduces toQ ψ + 3HQ ψ + µ 2 Q ψ = 0, which shows that µ also plays the role of Q ψ 's mass there. From the expression (4.68), one sees that if η 2 ⊥ 1, so that at a given time the bending parameter is large enough to overcome the other two contributions, then the entropic field experiences a transient tachyonic instability for k 2 /a 2 |m 2 s |, first noticed in [115]. As a series of recent works has shown [116][117][118][119][120][121], this does not lead to a background instability but rather to a transient exponential growth of fluctuations until (effective sound) horizon crossing. Through the derivative coupling between the two perturbations, whose strength is determined by the same parameter η ⊥ , this growth also affects the curvature perturbation. Thus, both adiabatic and entropic perturbations are subject to the same exponential growth (on sub-Hubble scales) compared to standard setups (see e.g. [122][123][124][125][126] for recent discussions about this regime of strongly non-geodesic motion).
Considering a brief period in which the bending parameter is large provides an explicit example of a transient non-adiabatic evolution leading to the dynamical appearance of an excited state. Using the language of Sec. 4.1, a large bending parameter determines a feature region, corresponding to the time interval in which η ⊥ 1, followed by an out region where the bending is negligible and the system is in an excited state.

Mode functions for strong sharp turns
Equations (4.66) and (4.67) can be re-expressed as equations of motion for Bogoliubov coefficients keeping track of the production of excited states during inflation. We may define time-dependent Bogoliubov coefficients aŝ where ζ(kτ ) = (1 + ikτ )e −ikτ is the massless dS mode function given by Eq. (3.8), whereas ν (−kτ ), with ν = 9/4 − µ 2 /H 2 , is the dS mode function for a massive field of mass µ. In fact, µ here coincides with the entropy mass defined in Eq. (4.69). Both mode functions respect Bunch-Davies initial conditions at kτ = −∞.
Then, it can be shown that α ζi (τ ), α ψi (τ ), β ζi (τ ) and β ψi (τ ) obey the following firstorder equations of motion where the coefficient matrices, satisfying det A = det B = 0 and AB = BA = 0, are given by Noteworthily, the first-order, coupled differential equations (4.72) are valid for any time-dependent η ⊥ (τ ) and preserve the relations (4.6) and (4.7) for the time-dependent Bogoliubov coefficients written in Eqs. (4.70) and (4.71). This can be shown by taking the time derivative of Eqs. (4.6), (4.7) and using Eq. (4.72) to show that they stay invariant as long as the Bogoliubov coefficients satisfy them at a given initial time. Given that ζ(kτ ) and ψ(kτ ) coincide with the standard single-field solutions for massless and massive fields respectively, Eq. (4.72) shows explicitly that turns in field space, parametrised by η ⊥ , induce the excitation of modes through the mixing of Bogoliubov coefficients. More to the point, starting from the Bunch-Davies initial conditions (4.3) in the in region where η ⊥ vanishes, these equations make manifest that the bending in the feature region feeds non-vanishing β coefficients, which become time-independent once the bending is over in the out region, to coincide with the ones defined in Eq. (4.4).
Although using a different method, i.e. directly solving for the mode functions Q Xi , analytical solutions were found explicitly in [44], with results generalised in [21], by considering a top-hat profile for the time dependence of the bending parameter 23 η ⊥ (N ), with height and width given by the constants η ⊥ and δ, and in the regime of a sharp turn δ 1. In particular, by setting Bunch-Davies initial conditions, the mode functions in the out and feature region take the simple forms (4.4) and (4.13), respectively. For the sake of completeness, we write here the corresponding functions in a convenient form. In the current setup, the model-dependent function g in the exponent of Eq. (4.13) is given by: while the Bogoliubov coefficients in Eq. (4.4) can be written as with θ k = 2e −δ/2 κη ⊥ + 2 arctan(κ/S) , (4.77) and Here, ξ parametrises the contributions from the potential and field-space curvature to the entropic mass (4.68), i.e. m 2 s = (ξ − 1)η 2 ⊥ (N )H 2 . (4.80) In this notation, ξ < 1 and η ⊥ = 0 correspond to the appearance of the transient tachyonic instability. These simple expressions are valid for scales k ∼ k * , where modes are exponentially amplified and one finds |α X i | |β X i | 1. Furthermore, from Eqs. (4.77) and (4.79) we observe the typical pattern of oscillations proper to a sharp feature, i.e. with a frequency given by ω ∼ 2/k out as described around Eq. 4.12. Here, the slowly varying k-dependent phase factor is given explicitly by ϕ k 2 arctan(κ/S).
For the numerical results shown throughout the paper we have considered the particular case ξ = −3 [44] corresponding to the case in which the entropy mass (4.69) vanishes (µ = 0) [127]. Apart from being a convenient choice, this case can be well-justified from a holographic perspective, wherein one can construct multifield models of inflation with the help of "fake" super-potentials [128]. As shown in [44], in this type of models a shift symmetry of the super-potential (with respect to one of the fields) enforces the value ξ = −3 throughout the full evolution of the multifield system, before, during and after the turn. Thus, we have normalised k in the previous expressions to the maximum of the primordial power spectrum for ξ = −3.
Moreover, in this example, the parameter defined in Eq. (4.1), quantifying how deep inside the horizon the maximally enhanced mode lies at the beginning of the out region, is given by γ = η ⊥ e −δ/2 η ⊥ , (4.81) while the maximal enhancement of the primordial power spectrum at k = k * , reads (4.82) 5 Theoretical bounds and observational prospects

Backreaction and perturbativity
It is well-known that a sharp feature along the inflationary dynamics (a transient period of non-adiabatic evolution) can lead to strong coupling of fluctuations at a given scale, jeopardising perturbative unitary and thus leading to a loss of theoretical control. Furthermore, the quanta copiously produced during this epoch can backreact and compete with the underlying background evolution rendering unsatisfactory the standard treatment of perturbation theory, see e.g. [80,[129][130][131].
In this section, as a proof of principle, we show how bounds coming from the requirements that unitarity and backreaction be under control can constrain a significant part of the parameter space of a given scenario that would otherwise lead to a signal above the threshold of observability in GWs detectors. In [21] we have derived tentative bounds applicable to the example of Sec. 4.7: perturbative control requires while the constraint for backreaction to be under control is where the C on the right-hand side in both inequalities has to be thought of as a fudge factor taking values O(1) − O(100) and which can be estimated in full computations. Note that, interestingly, the perturbativity bound (5.1) is in agreement with the constraint arising from energy conservation discussed in [132]. Let us characterise the enhancement of the scalar power spectrum at its maximum by the parameter where for definiteness, P 0 = 2.4 × 10 −9 is considered in this section. The bounds (5.1) and (5.2) can be rewritten in terms of p and γ as Perturbativity : Assuming that a signal is being potentially detectable if h 2 Ω GW 10 −14 , then the green area corresponds to the region of parameter space where Ω rad GW has observational relevance while the theoretical framework is under control. This allowed and relevant region of parameter space is extended by an additional strip (in yellow) once the Ω inf GW signal is also considered. For δ > 1 (red region) the feature described in Sec. 4.7 is not sharp and the solutions for the Bogoliubov coefficients does not hold.

Constraining the parameter space
Let us now understand how the bounds from the previous section constrain part of the parameter space that would naively lead to a pronounced signal in Ω GW . As discussed in the previous section, the amplitude of the gravitational-wave energy density generated during inflation and in the post-inflationary evolution will depend on γ and on the enhancement p of the primordial scalar power spectrum. In particular, we have shown that Ω inf GW ∝ 2 P 2 ζ γ 5 at its global maximum (4.50), so that we can rewrite it as where quantities with a tilde are meant to be evaluated at a given pivot point (γ,p), r i is the redshift factor defined in (2.19) and is an order-one number, almost independent of the pivot scale chosen. 24 The amplitude of Ω rad GW (at its maximum) depends instead only on the enhancement of the scalar power spectrum and can be written as with r r the redshift factor defined in (2.19) and c r again an order-one number. A certain enhancement p of the primordial scalar power spectrum P ζ fixes the amplitude of Ω rad GW . In contrast, the amplitude of Ω inf GW , for the same value of p, is still allowed to grow as γ 5 . Thus, one can imagine scenarios (part of the parameter space) characterised by a small enhancement of the scalar fluctuations compared to CMB scales, which lead to an unobservable amount of Ω rad GW but for which Ω inf GW is still large in amplitude. Naturally, in such a case, the backreaction and perturbativity constraints (5.1)-(5.2) are more easily satisfied. Furthermore, let us recall that the maxima of the two contributions Ω inf GW and Ω rad GW are located at different frequencies: from Eq. (4.50) we have k inf max 3.5k out and k rad max = 2/ √ 3 γk out (recall that γ 1). Thus, in the region of parameter space where the two peaks have similar amplitude, one can still hope to see the imprint of both in the full Ω GW .
For a given example of a sharp feature, the underlying parameters defining the mechanism are directly related to γ and p. In the case of a strong turn in field-space, these are given by Eqs. (4.81) and (4.82), i.e. p = δη ⊥ − ln √ 2 and γ = η ⊥ e −δ/2 . By fixing the product η ⊥ δ we fix p, and γ is then roughly given by p/δ, where δ is the parameter controlling the duration of the feature.
In this context, Fig. 7 can be thought of as a roadmap in parameter space for sharp features generating excited states and sourcing GWs. On the y-axis we have the enhancement p of the power spectrum. From Eq. (5.8) it follows that lines of constant p coincide with lines of constant Ω rad GW and that the larger the value of p, the larger Ω rad GW . Lines of constant Ω inf GW can be understood by looking at Eq. (5.6): increasing p and γ leads to an enhancement of Ω inf GW at its peak. Thus, if we take a specific value of Ω GW =Ω as a proxy for the threshold of detection in a given GWs observatory, all the parameter space above the corresponding lines of Ω rad GW and Ω inf GW would correspond, a priori, to a "positive detection" for that particular experiment. However, and this is the purpose of this section, perturbativity bounds could forbid a huge part of this parameter space. As an example, let us focus on the valueΩ 10 −14 that in the mHz band can be seen as the optimistic threshold of detectability for LISA. If we consider first the contribution from Ω rad GW , any parameter corresponding to the the points above the dash-dotted red line in Fig. 7 would give a signal in the LISA sensitivity band. Even so, by taking C 10 in Eq. (5.4) the region highlighted in grey would correspond to points where perturbativity is lost and the theory is not under theoretical control, leaving us with the highlighted green triangular region as the de facto available region of parameter space. Considering the contribution from Ω inf GW adds a strip in parameter space, highlighted in yellow, were Ω inf GW >Ω and still perturbativity is roughly under control. Note that in our specific set-up, the pair (p, γ) also determines the region in which δ < 1, cf. Eqs. (4.81), (4.82). The complementary region is marked in red as our analytical solutions to the dynamics of fluctuations is not accurate there.
This discussion applies equally to the bound (5.5) from backreaction. The resulting Fig. 9 would be the analogous roadmap in parameter space for this case. The only small difference is what happens once a given value of the slow-roll parameter is fixed. While in Fig. 7, is controlling only the size of Ω inf GW ∝ 2 , in Fig. 9, also determines the position of the backreaction constraints through Eq. (5.2). Thus, one should keep in mind that a smaller would ameliorate backreaction but at the same time would lift all lines of constant Ω inf GW . Summarising, by looking at Fig. 7 (Fig. 9) we can select specific points in parameter space for which perturbativity (backreaction) bounds are satisfied and Ω inf GW and/or Ω rad GW are potentially observable in a given experiment. For instance, the gravitational wave signal generated during inflation corresponding to the point in parameter space marked with a cross in Fig. 7 is shown in Fig. 8 together with the LISA sensitivity band. The small value of p implies a subdominant Ω rad GW 10 −5 · P 2 0 e 4·5 5 · 10 −15 compared to Ω inf GW . In Fig. 10 we show Ω rad GW and Ω inf GW for the two points in parameter space marked in Fig. 9. To conclude, scenarios in which a sharp feature leads to a significant enhancement of the power spectrum of the scalar fluctuations should be studied carefully, and at the same time, from a theoretical and phenomenological point of view. In particular, the significant constraining power coming from the -often neglected -requirement of remaining in a regime of theoretical control, demands a detailed determination of the corresponding bounds (5.1), (5.2) once a particular set-up is proposed.

Summary of results
We have derived analytical predictions for the enhanced SGWB sourced during inflation as induced by the presence of a scalar excited state. Excited states with large occupation numbers associated with small scales can dynamically emerge as the result of a sudden transition along the inflationary trajectory. We illustrate this concept by means of an explicit example in a multifield setup but analogous dynamics can emerge in other contexts. Allowing non-trivial dynamics during inflation not only provides an unconstrained and exciting phenomenological playground at scales smaller than the CMB; given the the- GW Ω inf GW , (γ, p) = (28, 6.5) Ω inf GW , (γ, p) = (13, 7) Ω rad GW , (γ, p) = (28, 6.5) Ω rad GW , (γ, p) = (13, 7) LISA sensitivity Figure 10: Ω rad GW and Ω inf GW corresponding to the points in parameter space highlighted with crosses in Fig. 9. For illustrative purpose we choose the same time of the feature N out = 28.5 for both scenarios. The position of the peak of Ω inf GW is then fixed -see Eq. (4.50)and larger γ implies a higher and wider signal. Conversely, increasing γ shifts the peak of Ω rad GW to higher frequencies, while its overall amplitude depends only on p, i.e. the enhancement of the power spectrum of the curvature perturbation.
oretical difficulties in building a consistent high-energy embedding for models displaying a prolonged period of standard slow-roll evolution, it is reasonable to assume that the few e-folds probed by current experiments might have been, in reality, rather exceptional.
Let us schematically summarise our main findings. First, and independently of the current application to the case of an excited state: We have derived a generic formula -Eq. (3.19)-for the tensor power spectrum sourced at second order by a generic multi-scalar system during inflation.
The non-trivial sum arising from properly quantising the multifield system leads to a N 4 enhancement of the tensor power spectrum, with N the number of scalar degrees of freedom, when all fields equally contribute. In particular, the former is given in terms of an integral over the scalar-field momenta and the time evolution. At the time when the excited state has been generated, the relevant enhanced modes are, by definition, deep inside the Hubble radius. Assuming, from that time onwards, an almost de Sitter like evolution of fluctuations, the time evolution can be explicitly integrated. Equation (4.41) provides one with a kernel that can be applied to compute the corresponding GWs background for any excited state with sufficiently large occupation numbers, i.e. whenever the Bogoliubov coefficients satisfy |α| |β| 1.
inflationary-era GW spectrum compared to the post-inflationary one, we find that Ω inf GW can easily be comparable to or larger than Ω rad GW . All in all, our results reinforce the importance of GW cosmology in reconstructing the history of the very early universe. Our work leaves several interesting open questions to pursue in the future. On the phenomenological side, for instance, it is possible to obtain more accurate templates to reconstruct the GW signal by taking into account the subdominant component of the mixed contribution in Eq. (2.10), as well as consider more complicated dynamics once excited states have emerged. In addition, although we have estimated the size of backreaction effects, it is necessary to have a more rigorous understanding of how they affect the evolution of the system. On the theoretical side, it would be interesting to study cases where the particle production is not large. There, the classical treatment we have used in this work is not accurate enough and one has to resort to the full in-in formalism. these excited states have support from τ out onward until the end of inflation at τ = 0 (out region), in order to avoid the complications discussed in Sec. 3.2, coming from the i prescription, we may as well focus on the spectrum therein. Let us thus use the notation of Eq. (4.25) and label P out t the part of the spectrum corresponding to the out region. Furthermore, let us ignore the momentum conserving delta function and also write P out t (k) ≡ lim kτ →0 P out t (k, τ ) as in Sec. 4. The in-in formalism yields two one-loop contributions for the tensor power spectrum: This is nothing but rewriting the modulus squared in the second line of Eq. (3.19) (recall that g k is real) as the sum of the squared 25 real and imaginary parts of the quantity X Q Xi (p, τ 1 )Q Xj (q, τ 1 ). The index α in the right-hand side just reorganises the summation over ij; it takes 2N 2 values (for N scalar fields).
Next, let us note that since the scalar source enters in the same way in both methods (the only difference lies in how the external gravitons connect to the internal vertices), we may use this notation to also reorganise the one-loop diagrams as where we have taken into account the combinatorial factor of 2 coming from the distinct h-contractions. Note that in this way, we may compare (A.5), (A.7) and (A.8) term by term and then resum. Finally, by noting that due to the reality condition on f α , the integrand of (A.8) is symmetric under τ 1 ↔ τ 2 , we may disentangle the integrals in P 1;B t and rewrite it as [75]  Upon doing so, we obtain P out t;A (k) + P out t;B (k) = P out t (k), (A.10) which is the desired result.

B Approximate expression for the GW spectrum
For simplicity we suppress all labels on the Bogoliubov coefficients and proceed as if there was only a single scalar mode. From Eq. (4.41), the spectrum of induced GWs takes the form where we find it convenient to perform the integration using the variables d, s defined in Eq. (3.21), and the kernel F (x, y, k) can be read off from Eq. (4.41).
We also restrict our attention to situations where the excited state induces a peak in the scalar power spectrum, as this is the case most relevant for phenomenology. In terms of the Bogoliubov coefficients, this implies that |α(k)| 2 will exhibit a peak, see Eq. (4.15). For a first analysis, we consider the extreme case |α(k)| 2 = C δ[log(k/k * )], i.e., where |α(k)| 2 is given by a single spike. Inserting this profile into Eq. (B.1) one finds P out, δ t (k) = C 2 k 2 * In practice, however, the peak in |α(k)| 2 will have finite extent. If it is sufficiently narrow, the above single-spike result is expected to still give a good approximation to P out t (k), especially for k ∼ k * . However, if the peak is not narrow enough or for sufficiently small k this is not the case anymore.
To illustrate this point, for a peak of finite width ∆k we write where the Heaviside theta functions enforce the finite extension of the peak, with k 1 < k * < k 2 . Inserting Eq. (B.3) into Eq. (B.1) the effect of the Heaviside theta functions is to restrict the integration domain to We can then confirm that for a sufficiently narrow peak with ∆k k * one recovers the single-spike result. In this case the two instances of P in Eq. (B.1) only have significant overlap when the two arguments are close to each other, i.e. in the vicinity of d = 0. To approximate this, we can thus set d = 0 and integrate over the region of overlap as given in Eq. (B.4). For the d-integral this implies integrating up to d max = ∆k/k, which gives a factor of ∆k/k. This implicitly assumes that we are considering values k ∼ k * , so that ∆k/k is small. The s-integration is to be performed over s ∈ [2k 1 /k, 2k 2 /k], i.e. over an interval of width ∆k/k enclosing the value 2k * /k. For sufficiently small ∆k/k the integrand can be taken as constant over this interval. The integral can thus be approximated by evaluating the integrand at s = 2k * /k and multiplying by the width 2∆k/k. Upon including a factor of 1/2 that accounts for the triangular shape of the (d, s)-integration domain, one recovers the single-spike result (B.2).
We now examine how this analysis is changed if the peak in |α(k)| 2 is not narrow, or if we are interested in small values of k. In both these cases the interval ∆k/k is not necessarily small, which will affect how the d-and s-integrals can be performed. In particular, consider ∆k/k > 1. In this case d max = ∆k/k > 1 and the constraint on the d-integral due to the finite size of the peak is void, as the range of integration in Eq. (B.1) is d ∈ [0, 1]. As long as the integrand only varies slowly with d, we can proceed as in the narrow peak case, set d = 0 and perform the integral. The difference however is that this does not give a factor of ∆k/k anymore. The integral over s is again over an interval of width 2∆k/k, which is no longer small. Still, if the integrand only varies slowly over this range, we can proceed as before. Then, for a peak in |α(k)| 2 and for ∆k/k > 1 the spectrum of induced GWs can be written as in the single-spike case, albeit with one less factor of ∆k/k, i.e. P out t (k) ≈ C 2 k * ∆k κ −1 F κ −1 , κ −1 , k , (B.5) with κ = k/k * . Recall that an assumption for arriving at this result is that F (x, y, k) in Eq. (B.1) varies sufficiently slowly. Here we did not include the Heaviside theta function from Eq. (B.2), as this piece is also affected by the finite width of the peak in |α(k)| 2 . The correct factor is not very important here, as this mainly affects the UV part of the spectrum, where this approximate expression is not expected to hold. In case of sharp features, it is reasonable to assume that k * /∆k 1. In the main text we have also redefined the constant C of this appendix such that C = A(2P 0 ) −1 . In this way A gives a direct estimate of the peak of the primordial scalar power spectrum -see Eq. (4.43). Thus, we finally arrive at the approximation in Eq. (4.42) used throughout the paper. This, when appropriately applied to the various contributions of the GW spectrum, will provide a very good match to numerical results for k k * . As we explain in the main text, this includes the regime of maximal amplitude of the GW spectrum.