New constraints on axion-gauge field dynamics during inflation from $Planck$ and BICEP/Keck data sets

We present new constraints on spectator axion-${\rm U}(1)$ gauge field interactions during inflation using the latest $Planck$ (PR4) and BICEP/Keck 2018 data releases. This model can source tensor perturbations from amplified gauge field fluctuations, driven by an axion rolling for a few e-folds during inflation. The gravitational waves sourced in this way have a strongly scale-dependent (and chiral) spectrum, with potentially visible contributions to large/intermediate scale $B$-modes of the CMB. We first derive theoretical bounds on the model imposing validity of the perturbative regime and negligible backreaction of the gauge field on the background dynamics. Then, we determine bounds from current CMB observations, adopting a frequentist profile likelihood approach. We study the behaviour of constraints for typical choices of the model's parameters, analyzing the impact of different dataset combinations. We find that observational bounds are competitive with theoretical ones and together they exclude a significant portion of the model's parameter space. We argue that the parameter space still remains large and interesting for future CMB experiments targeting large/intermediate scales $B$-modes.


Introduction
A stochastic gravitational wave background (hereafter SGWB) at all frequencies is a generic prediction of the inflationary paradigm [1,2]. The importance of measuring such primordial signal cannot be overstated, since a detection would provide strong evidence for cosmic inflation [3][4][5][6]. To reach this goal, an extensive experimental effort, targeting the SGWB spectrum at different frequencies, is ongoing and will continue throughout the next decade and beyond (see e.g. [7] for a review). The amplitude of the SGWB is usually parametrized by the ratio of the amplitudes of the tensor and scalar modes power spectra, the tensor-to-scalar ratio r. Currently, only upper bounds exist on r, the tightest one being r < 0.034 at 95% C.L. 1 [8] placed using a combination of Planck [9] and BICEP/Keck [10] CMB experiments data, through the imprint primordial tensor modes are known to leave in the B-mode polarization pattern of the CMB [11,12]. Being the most sensitive and the closest in the future timeline among all planned probes of the primordial SGWB [7], a positive gravitational wave (GW) detection is likely to come first from future CMB B-mode experiments, such as the LiteBIRD satellite [13] and the ground-based CMB-S4 [14], if tensor-to-scalar ratio reaches r ∼ 10 −3 at CMB scales.
Merely detecting r, however, does not immediately enable us to discriminate between different possible origins of the primordial SGWB. In the simplest scenario of single field inflation, realized by a slowly rolling scalar field minimally coupled to gravity, the SGWB is produced by the quantum vacuum fluctuations of the metric [1,2]. In this simple setup, the tensor-to-scalar ratio r can be related directly to the energy scale of inflation [15]. Furthermore, the SGWB spectrum produced within this framework is known to be (i) nearly scale-invariant (with a slight red-tilt), (ii) nearly Gaussian and (iii) non-chiral (i.e. parity-conserving). However, the relation between the energy scale of inflation and r, together with all the properties enunciated above, can be violated if an energetically-excited extra particle content is present during inflation, feeding the stress-energy tensor in the perturbed Einstein equation for the tensor modes of the metric. This intriguing possibility makes testing the scale dependence, Gaussianity and the chirality of the SGWB very compelling, if r is detected by future experiments [16].
Sourcing gravitational waves with additional matter fields during inflation can have, however, undesirable side effects: the sources are always (at least) gravitationally coupled to the sector responsible for the accelerated expansion, and therefore enhance not only tensor but also strongly non-Gaussian scalar modes (see e.g. [17][18][19][20]). Avoiding overproduction of such perturbations to comply with the tight bounds enforced by the CMB data on non-Gaussianity, while simultaneously maintaining a visible amplitude for the sourced SGWB signal, is therefore a necessary ingredient of any successful inflationary model aiming to achieve GWs of secondary/non-vacuum origin.
Two of the most studied mechanisms that are capable of successfully realizing the above scenario, involve production of a SGWB from amplification of gauge fields of the Abelian [17,[21][22][23][24][25][26][27][28][29][30][31][32] and non-Abelian [33][34][35][36][37][38][39][40][41][42][43][44][45] kind. Particle production of the gauge fields from inflation has been considered in the contexts of several cosmological phenomena, such as the generation mechanisms of magnetic field [46][47][48][49][50][51][52][53][54][55][56] and the matter-antimatter asymmetry in our universe [57][58][59][60][61][62][63][64], among many others. Specifically, in this work we focus on the sourcing of GWs through the Abelian U(1) gauge field fluctuations and in this context we consider a spectator sector including a generic pseudo Nambu-Goldstone boson (i.e. an axion field) coupled to an Abelian gauge field [25,28]. In this setting, inflation is realized through a standard inflaton sector minimally coupled to gravity and whose energy density is dominant with respect to that of the spectator axion and gauge field, thus allowing controlled production of scalar perturbations. In order to break the conformal invariance in the gauge field sector and allow for the amplification of the fluctuations in the latter with subsequent enhancement of tensor modes during inflation, the axion and the U(1) field are considered to interact through a Chern-Simons term [65,66]. The resulting gauge field amplification (and also its impact on the scalar curvature perturbations) is controlled by the transient rolling of the axion along its potential. In particular, for the specific realizations of this scenario, we will consider two different choices for the spectator axion potential, both leading to localized gauge field amplification: i.e one with a standard cosine-type potential [25,67]; and the other endowed with a string-inspired non-compact 2 axion potential [28]. Excitingly, due to parity violating nature of Chern-Simons interaction, only one of the two helicities of the gauge field fluctuations is amplified by the rolling spectator axion, resulting in a fully chiral SGWB.
Gauge field sources leave distinctive signatures in the primordial SGWB compared to the standard quantum fluctuations of the metric. More specifically, we can exploit the characteristic scale-dependence of sourced tensor modes to distinguish between the two [24,25,28,70]. Another possibility is to look for the strongly non-Gaussian signature in the bispectra of sourced gravitational waves at CMB [17,26,71,72] and interferometer scales [70,73]. Finally, amplification of U(1) gauge field sources during inflation is a parity-breaking process and therefore non-zero parity-violating correlations in the CMB angular power spectra [21,74] and bispectra [23,26,72,75] are expected, together with circularly polarized gravitational waves at scales relevant for laser interferometers [76].
In this work, we are going to focus our attention to the scale-dependent signatures of the spectator axion-U(1) gauge field dynamics during inflation in the angular power spectra at large/intermediate CMB scales 3 , with the aim of deriving constraints on the parameter space of spectator axion gauge field interactions. For this purpose, we utilize the state-of-the art CMB dataset for temperature and polarization provided by the Planck satellite [77] in its latest incarnation, complemented by the latest BICEP/Keck B-mode polarization data [10]. In particular, Planck yields the current best measurement of temperature and (E and B-mode) polarization at the largest scales, which are in fact accessible only from space, while the ground-based BICEP/Keck contributes with the best measurement of intermediate scales B-modes to date. We explore the likelihoods using a frequentist method: the profile likelihood. Despite being widespread within the particle physics community [78,79], this statistical technique is relatively less common in cosmology 4 , compared to Bayesian methods. Nonetheless, profile likelihoods present several advantages over the latter: the global maximum likelihood solution is guaranteed by construction and, moreover, parameters estimates are independent from prior distributions and model parametrization, thus making profile likelihoods immune to "volume effects", appearing during marginalization [81].
The structure of the paper is as follows. We start in Section 2 by reviewing the spectator axion-U(1) model and assessing the bounds on its parameter space imposed by theoretical self-consistency (namely from perturbativity and backreaction considerations, see subsection 2.2). In Section 3 we provide instead the latest observational bounds on the model parameters from Planck and BICEP/Keck data. We offer a detailed interpretation of the constraints for different choices of the model parameters, highlighting differences among the two axion potentials, and analyze the impact of different dataset combinations. Finally, in Section 4, we show that observational constraints on the model are competitive with theoretical ones and update the available parameter space of the model. We conclude by suggesting a possible future path for discriminating the spectator axion-U(1) model we consider from the conventional single field realizations of inflation at CMB scales.
2 Theory: A spectator axion-U(1) gauge field model As we mentioned in the introduction, the particle content we consider during inflation contains a spectator axion-Abelian U(1) gauge field sector along with a canonical inflaton sector that does not exhibit direct interactions with the former and both sectors minimally coupled to gravity. The action describing this system is given by [22,25,82], where the first term represents the standard Einstein-Hilbert action, φ is the inflaton and the spectator sector includes an axion-like field χ, U(1) gauge field A µ with the antisymmetric field strength tensor F µν = ∂ µ A ν − ∂ ν A µ . In the spectator sector L spectator , we consider an axion-like field χ that enjoys a(n) (approximate) shift symmetry, taking into account the gauge invariance of the U(1) field, χ is then expected to interact with gauge fields through the leading order dimension five operator 5 that takes the following form whereF µν ≡ √ −g µνρσ F ρσ /2 is the Hodge dual of the field strength tensor F µν , f is the axion decay constant, λ is a dimensionless coupling constant and µνρσ is an anti-symmetric tensor satisfying 0123 = g −1 . Background evolution. We consider an inflationary setup where the spectator sector fields provide subleading contribution to the total energy density during inflation. This implies that energy densities of the scalar fields in the model (2.1) obey ρ φ ρ χ where ρ X =Ẋ 2 /2 + V (X) (with an overdot denoting derivative with respect to cosmological time t) for X = {φ, χ} and assuming negligible backreaction from gauge field fluctuations ρ A ρ χ (see Section 2.2), we have such that quasi-dS expansion is completely dictated by the inflaton's potential. Furthermore, we assume that V (φ) is flat enough to support sufficiently long quasi dS expansion but otherwise we let it unspecified as the fine details regarding the inflaton's dynamics is irrelevant for the discussion we present below. With this assumption, we will treat Hubble rate H as constant and denote the scale factor during inflation in conformal time as a(τ ) = −1/(Hτ ) 6 , −∞ < τ ≤ 0. On the other hand, if the spectator axion χ is displaced from its global minimum, it can also roll down in its potential, albeit in the slow-roll regime thanks to sufficiently flat potential U (χ) endowed with shift symmetry. Due to the slow-roll assumption, we require that χ's background dynamics should obey during inflation when scales associated with CMB observations exits the horizon. In what follows, we will briefly review the impact of such a slowly rolling spectator axion on the behavior of gauge field fluctuations. For the clarity of the discussion, initially we will not specify the explicit form of the potential U (χ) before we introduce explicit axion models (see Section 2.1) that we analyze in this work. Amplification of gauge field fluctuations. For a dynamical spectator axion field with a time dependent profile χ(t), the interaction term can no longer be treated as a surface term in the action (2.1). As a result, the dispersion relation of A µ becomes modified, leading to copious production of its fluctuations provided that spectator axion has a non-trivial velocityχ = 0. To see this, we first decompose gauge field fluctuations into Fourier modes in Coulomb gauge (A 0 = 0) as [86], Another possibility is the coupling of a shift symmetric scalar to the fermion current via dimension five operator, see e.g. [83][84][85] for phenomenological implications of this scenario. 6 In this work, we disregard terms that are subleading in slow-roll expansion.
where h.c represent hermitian conjugate of the first term in (2.5) and the helicity vectors obey together with the commutation relations of annihilation/creation operators â λ ( k),â † λ ( k ) = δ λλ δ( k − k ). Inserting the decomposition (2.5) in the spectator part of the action (2.1), the equation of motion (EoM) for the gauge field mode functions in a flat FLRW background satisfy where we defined dimensionless time variable x ≡ −kτ and the effective coupling ξ between spectator axion and gauge field. Without any loss of generality, we work with ξ > 0 andχ < 0 so that spectator axion rolls down on its potential from positive large to small values χ ≥ 0. As we mentioned before, the correction that appear in the dispersion relation (2.6) arise through the coupling (2.2) whose parity violating nature is apparent from its alternating sign ±. In particular, when the modes are deep inside the horizon (x = k/(aH) 1), the correction term is negligible and the gauge field obeys the standard dispersion relation. However, as the modes stretches outside the horizon, it becomes dominant for x = k/(aH) 2ξ, leading to instability for one of the circular polarization state of the gauge fields. In our conventions (ξ > 0 &χ < 0), A − state experiences a tachyonic instability while A + stays in its vacuum. The roll of spectator axion fieldχ = 0 therefore induces production of helical gauge fields and assuming the roll of axion with a constant velocity ξ constant, the late time (−kτ → 0) gauge field mode functions amplified this way typically exhibit an exponential amplitude that is regulated by the axion's velocity [86]: (2.7) In the model (2.1), although the spectator sector does not exhibit direct couplings with the visible sector fluctuations such as the inflaton δφ and metric h ij , the influence of the particle production processes in the gauge fields inevitably mediate to the visible sector perturbations through gravitational interactions. Below, we will briefly review the impact of the gauge field sources on these fluctuations. Tensor perturbations sourced by vector fields. To study the influence of gauge field fluctuations on the tensor perturbations, we focus on the transverse traceless metric perturbation g ij = a 2 (τ )[δ ij +ĥ ij (τ, x)] and decompose it into circularly polarized states λ = ± in Fourier space asĥ λ (τ, k) = Π ij,λ ( k)ĥ ij (τ, k) where Π ij,λ is the polarization tensor obeyinĝ k i Π ij,λ ( k) = 0, Π * ij,λ Π ij,λ = δ λλ and Π * ij,λ ( k) = Π ij,−λ ( k) = Π ij,λ (− k). Expanding the action (2.1) up to third order in fluctuations includingÂ i andĥ ij , it can be shown that the mode equation of graviton polarization states h λ is sourced by the transverse, traceless part of the energy momentum tensor that is composed of gauge field fluctuations [22]: where we defined dark "electric" and "magnetic" ijk ∂ jÂk in a flat FLRW universe, in analogy with standard model electromagnetism. Scalar perturbations sourced by vector fields. The impact of particle production processes on the visible scalar fluctuations is also encoded indirectly by the presence of gravitational interactions [19]. In particular, integrating out the non-dynamical scalar metric fluctuations such lapse δN and the shift N i reveals a mass mixing between inflaton δφ and spectator axion δχ fluctuations and opens up a channel that can influence the curvature perturbation 7 R H δφ/φ through the inverse decay of gauge fields: A i +A i → δχ → δφ ∝ R. Dynamics of this contribution can be understood by first studying the influence of particle production on the spectator axion fluctuations δχ through, Then focusing on the inhomogeneous solution of the δχ fluctuations in (2.9), one can compute the conversion of δχ to δφ via to find the the part of curvature perturbation R that is sourced by the amplified gauge fields.
In this work, we are interested to the extent of which gauge field sources can influence tensor perturbations consistent with other observations at CMB scales. In this context, it has been recently realized that if the spectator axion χ rolls for a large-amount of e-folds ∆N χ 1 during which the scales associated with CMB observations exit the horizon 8 , the sourced contributions to the R becomes sizeable due to the sensitivity of gauge field amplitudes and δφ − δχ mass mixing (see eq. (2.10)) on the spectator axion's velocity ξ ∝ |χ| [19]. In particular, insisting on a large secondary contribution to the tensor fluctuations through eq. (2.8) generically leads to an exceedingly large scalar non-Gaussianity at CMB scales [19,89]. The origin of the difficulty in efficiently enhancing tensor perturbations compared to the scalars can be readily seen from eqs. (2.8), (2.9) and (2.10), by realizing that the sourced part of both perturbations arise via non-linear terms including the same amount of gauge fields. Notice however from eqs. (2.9) and (2.10) that the efficiency of the process A i + A i → δχ → δφ ∝ R is highly sensitive to the behavior of the spectator axion's velocity as both gauge field production (see eq. (2.7)) and mass mixing δφ − δχ have dependence onχ. In what follows, we will discuss two spectator axion models that exhibit a localized velocity profile that can overcome the aforementioned limitations on scalar fluctuations.

Transiently rolling spectator axion models
In order to minimize the influence of particle production on the curvature perturbation and to render secondary GWs sourced by gauge fields viable, we will consider models that can lead to localized gauge field production where the spectator axion transiently rolls on potentials of 7 Regarding other direct contributions to R from fluctuations in the spectator axion δχ and gauge fields Ai, in this work we will consider a spectator axion model that rolls down to its minimum long before the end of inflation (see Section 2.1), so that the contribution of δχ on the late time curvature perturbation R can be neglected [25,82]. On the other hand, the contribution from gauge fields is proportional to the absolute value of Poynting vector, a| S| = a| E × B| which is also negligible at late times as the particle production saturates at super-horizon scales and the resulting electromagnetic fields decay as E, B ∼ a −2 [28]. For the purpose of evaluating phenomenological implications of the model (2.1), we therefore adopt the standard relation R ≡ Hδφ/φ in this work. 8 The regime of validity of the perturbative description of gauge field production is also questioned in [87,88]. the following form [25,28]: where µ and Λ are parameters of mass dimension one. The first model (M1) features a standard shift symmetric potential (see e.g. [67]) with the size of the axion modulations is set by Λ. In this model, the motion of the axion is contained within the maximum (χ = πf ) and the minimum (χ = 0) where the slope U (χ) vanishes. Therefore at large (early times) and small field values (late times), axion rolls with very small velocities whereasχ obtains relatively large value at an intermediate time when χ passes through an inflection point χ * = χ(τ * ) with U (χ * ) = 0 where the slope of the potential U (χ) becomes maximal.
In the second model (M2), the axion field range is extended via a monodromy term [68,69] proportional to a soft symmetry breaking mass parameter µ and χ is assumed to probe step-like feature(s) 9 in the "bumpy" regime, Λ 4 µ 3 f . Similarly to the first model, in the plateau like region and towards the global minimum 10 , the spectator axion acquires very small velocities where U → 0 but obtains a transient peak when the slope of the potential U (χ) becomes maximal at the cliff region, in particular at inflection point denoted by χ * (See Fig. 1). 9 In the bumpy regime, depending on the initial conditions (χ f ) spectator axion can probe multiple step-like features during inflation. In this work, we assume that χ traverse only one such region on its potential during which observable scales associated with CMB exits the horizon. 10 As shown in the right panel of Fig. 1, the roll of χ towards the global minimum (χ = 0) can be captured by modifying the monomial term as µ 3 χ → µ 3 f [ 1 + (χ/f ) 2 − 1] so that the axion potential (2.11) interpolates between µ 3 χ and (µ 3 /f )χ 2 from large to small field (χ/f → 0) values respectively. By construction this modification is designed to affect the motion of χ far away from the inflection point χ * where the axion's velocity peaks and therefore we do not expect it to influence χ's velocity profile (2.12) and the resulting phenomenology we discuss in this work (see below) which are derived by assuming the potential shown in (2.11). In short, it only stands to ensure a smooth rollover to the minimum χ → 0 with a small velocity U ∝χ → 0 within the slow-roll approximation we are undertaking.
Assuming slow-roll motion (2.4), for typical field ranges dictated by the scalar potentials in (2.11), the spectator field velocityχ and the effective coupling strength ξ = −λχ/(2Hf ) therefore obtains a peaked time dependent profile [25,28]: where the subscript * denotes the value of a quantity at the time when the axion passes through the inflection point (See Fig. 1). It is clear from (2.12) that ξ obtains its peak value at τ = τ * with a maximal value that can be parametrized in terms of the dimensionless coupling constant λ as . (2.13) The width of the time dependent peak in ξ is mainly controlled by the dimensionless ratio δ that essentially characterize the mass of the spectator axion in its global minimum δ ≈ m 2 χ /H 2 . In particular, for larger δ (heavier axion), restoring force towards the global minimum is larger so that χ traverse the inflection point faster, resulting with a sharper peak in ξ. In other words, δ is a measure for the acceleration (ξ/(ξH) =χ/(χH) ∼ δ) of the spectator axion as it rolls down on its potential. At this point it is worth mentioning that due to the slow-roll approximation (2.4) we are undertaking, we are restricted to δ 3. In this work, to derive observational and theoretical constraints on the spectator axion-gauge field model (2.1), we will focus on the following cases δ = {0.2, 0.3, 0.4, 0.6} corresponding to increasingly sharp rise in the velocity of χ.
The peaked structure of ξ profile introduces a critical scale τ −1 * = k * in the EoM (2.6), corresponding to the scale that exits the horizon when the axion's velocity is maximal (i.e when ξ = ξ * ). Since the tachyonic mass of the U(1) field in (2.6) is maximal around this point, it results in a scale dependent growth of the gauge field fluctuations where only modes whose size is comparable to the horizon size at τ = τ * , i.e k ∼ O(1)a * H * , are efficiently amplified. The scale dependent amplification of gauge field modes can be accurately studied using the semi-analytic techniques discussed in [25,28], which is what we will utilize in our analysis. Below we review the impact of such scale dependent vector field production on the 2-point correlators of tensor and scalar fluctuations during inflation. Scale dependent perturbations from gauge field sources. In addition to the standard vacuum fluctuations driven by the quasi-dS background, the perturbations in the observable sectorX = {R,ĥ ± } pick up a sourced contribution from the enhanced gauge field fluctuations that can be described by the particular solutions of (2.8) and (2.10) (see also (2.9)):X = X (v) +X (s) where the superscripts denote vacuum and sourced modes, respectively. We define the power spectra of R and h λ as Since the origin of the vacuum and sourced part of scalar and tensor perturbations are different, these contributions are statistically uncorrelated and therefore the resulting total power spectra can be simply described by the sum of vacuum and sourced part auto-correlators as where the vacuum contributions are given by the standard expressions: with r v denoting the vacuum tensor-to-scalar ratio and φ ≡φ 2 /(2H 2 M 2 pl ) is the slow-roll parameter controlled by the inflaton sector. Note that due to the parity violation in the gauge field sector (A − A + ), it is sufficient to take into account P On the other hand, for the transiently rolling spectator axion models we described above, the sourced power spectra in (2.16) inherit the scale dependence of the gauge field that leads to a Gaussian spectral shape [25,28], where j = {R, ±}. The functions f c 2,j , σ 2,j , x c 2,j control, respectively, the amplitude, the width, and the position of the peak of the sourced signal, which depend on the background model of the spectator axion through the parameters ξ * and δ we discussed above and therefore on the underlying scalar potential (2.11) in the spectator axion sector. For representative choices of the background parameter δ, we present accurate formulas for the amplitude f c 2,j , width σ 2,j and the location x c 2,j of the peak in terms of the effective coupling ξ * in Tables 3-6 in the appendix. Notice that due to parity violating nature of gauge field production, sourced tensor perturbations satisfy f 2,− f 2.+ and therefore it is maximally chiral. Note also from the Tables 3-6 that the amplitude f c 2,j of the sourced signals is exponentially sensitive to the effective spectator axion-gauge field coupling ξ * that parametrizes the efficiency of particle production in the gauge field sector.
In (2.16), the sourced power spectra (2.18) introduce a Gaussian bump feature on top of the standard quasi scale-invariant spectra (2.17). If the former feature is dominant, the total power spectra becomes highly scale-dependent. Clearly, such a scale dependence should not overwhelm the scalar power spectrum and should be consistent with the CMB temperature (T) and polarization modes (E,B) data. Our main goal in this work is to derive constraints on the rolling spectator axion-U(1) gauge field models from Planck and BICEP/Keck 2018 data in order to see to what extent we can realize a chiral, synthetic component of tensor modes through the spectator axion-gauge field dynamics.
However, before we continue our discussion in this direction, we would like to identify and check the parameter space of the spectator axion-gauge field model (2.1) consistent with backreaction and perturbativity considerations first discussed in [87,88].

Limits on backreaction and perturbativity
Induced by the gauge field amplification in the spectator sector, the derivation of the scale dependent contributions (2.18) to the total power spectra assumes that backreaction of the spectator fields on the background evolution is negligible and vector/scalar fluctuations in the spectator sector stay in the perturbative regime. In this section, we study the limitations on the size of the sourced signals from these effects. In our analysis, we will closely follow [28,87,88] which we refer the reader for many details presented below.

Backreaction constraints
Since we assume that axion like field is a spectator and does not contribute effectively to the total energy density during inflation, we need to make sure that ρ χ V φ 3H 2 M 2 pl is satisfied during the motion of χ. As shown explicitly in [28,88], for both of the transiently rolling spectator axion models we consider, potential energy of the axion U (χ) always dominates over the kinetic energy E kin,χ =χ 2 /2 which reaches its maximal value at the inflection point when τ = τ * . Therefore, to ensure that spectator χ does not contribute to the background energy density, it is sufficient to enforce where the "max" refers to the maximum value of the potential energy during the rollover of χ.
An upper bound on f /M pl : In the compact spectator axion model, assuming χ starts its motion close to the maximum χ in π/f of the potential, the maximal value of the potential energy density is set by the height of the oscillatory potential in (2.11) and is given by On the other hand, in the non-compact axion model, the maximal value of the potential depends on the initial conditions as U (χ)| max µ 3 f [χ in /f + 1] in the Λ 4 µ 3 f regime where χ in /f 3π/2, assuming χ traverse a single cliff-like region in its potential [28] (See Fig. 1). In terms of the dimensionless ratios we defined in (2.11), the first condition (2.19) gives (2.20) Therefore, the backreaction constraints translates into an upper bound on the ratio between two fundamental parameters in our model, namely f /M pl .
A lower bound on f /M pl : Next, we need to make sure that gauge field production in the spectator sector does not influence the background evolution of χ. Since the gauge field mode functions are amplified at the expense of spectator axion kinetic energy, we therefore require that maximum energy density contained in the gauge fields to be smaller than peak kinetic energy of spectator axion: ρ A,max (χ 2 /2) τ =τ * ρ φ χ, * /3 where we defined the slow-roll parameter χ =χ 2 /(2H 2 M 2 pl ). Using the definition of effective coupling in (2.6) and (2.12), this condition can be cast into the following form  where ρ A,max is the maximum value of the gauge field energy density which is typically obtained when τ /τ * ∼ 10 −2 [28,88]. For both spectator models we consider, we computed the quantity ρ A,max /(ρ φ φ ) using the formulas provided in [28,88]. Combining the resulting expressions with (2.21), for a given choice of δ, we derive a lower bound on the ratio f /M pl in terms of where for both models and all the δ values we consider, the precise values of the coefficients c 0 and c 1 can be found in Table 1. Comparing the two models we consider from the table, we see that the lower bound is in general less restrictive in the non-compact axion model (M2) compared to the rolling axion model (M1) with the standard cosine potential. Notice that, increasing δ relaxes the bounds as in this caseχ is maximal for a shorter amount of time, reducing width of the gauge field modes that are effected by the roll of χ and in general the efficiency of particle production. Finally, it should be clear from (2.22) that, reducing r v relaxes the bound further where the allowed region for f /M pl increases.

Perturbativity constraints
The production of scale dependent, chiral GWs of non-vacuum origin in the spectator sector typically demands an exponentially large amplitude in the gauge field sources during the times/at scales when the observable effects are produced. Therefore, one may wonder if large amplitudes obtained by the gauge field fluctuations can drive the system out of the perturbative regime which was the intrinsic assumption we made in deriving the sourced templates of scalar and tensor perturbations in (2.18). In the following analysis, our aim is therefore to establish the regime for which these results are under perturbative control. For this purpose, we consider two main requirements that the spectator models (2.1) (and (2.11)) we focus should fulfill [88]: 1. Higher order loop corrections induced through the interaction (2.2) do not spoil the leading order estimates for the amplified gauge field mode functions. This criterion can be written in terms of the model parameters as [28,88], where prime denotes the two point function without the corresponding delta function and the expression in the nominator/denominator represent leading order loop correction to the gauge field propagator and the corresponding three level result, respectively.  2. For the second criterion, we will demand that the interaction (2.2) does not induce a variance δχ 2 that is larger than the typical classical field excursion χ cl . We therefore require [28,88], where we described the numerator as an integral of the leading order loop contribution to the axion's power spectrum 2π 2 P (1) . In the first model M1 with standard axion modulations, while it is natural to identify χ cl → f [88], due to non-compact nature of axion in the second model (M2), we will use χ cl → χ in 3πf /2 assuming χ rolls over one bump like region in its potential before reaching its global minimum at χ = 0 [28].
As indicated by the expression (2.23), the first criterion is time and scale dependent. To evaluate this expression, we will focus on the mode that is most amplified by the rolling axion, k = 5k * by evaluating the expression at a late time τ /τ * → 0 at which the gauge field mode functions are maximally enhanced. As shown in [28,88], this strategy is sufficient to derive strongest constraints as the growth in P A (and also for P χ ) saturates at late times for the most amplified mode. The second criterion (2.24) on the other hand arise as an integral over modes and hence scale independent. Following the procedure outlined in [28,88], we evaluated these criteria for both models we consider and for δ = 0.2, 0.3, 0.4, 0.6 corresponding the increasingly faster rolling axion. Similar to the backreaction constraints we derived earlier, at fixed δ, the resulting constraints can be interpreted as a lower bound on f /M pl in terms of r v and the effective coupling ξ * as where we provide explicit values of the coefficients p,p in Table 2. As can be inferred from the tables, for all the δ choices we focus, bounds from the renormalization of the gauge field wave functions dominate over the second criterion 2) above, and hence we will ignore the latter hereafter. Another information that we can obtain form Table 2 is that constraints on f /M pl tend to be weaker in the second model M2 compared to M1.
Summary of backreaction and perturbativity constraints. Focusing on the range 3 < ξ * < 6.5 of effective coupling within which interesting phenomenology from spectator axion-gauge field dynamics can arise, we compared the lower bounds obtained on f /M pl from backreaction (2.22) and perturbativity considerations in (2.25) using Tables 1 and 2. In this way, we found that for δ = 0.2, backreaction constraints dominate over P A for ξ * < 5.97 in M1 and for all ξ * range we quoted above in M2. For δ = 0.3, the range of domination of the backreaction constraints reduces to ξ * < 5.72 in M1 while it still dominates over the perturbativity for all ξ * in M2. Increasing δ further makes the perturbativity bound stronger compared to backreaction. For example, for the choice of δ = 0.4, backreaction dominates for ξ * < 5.45 in M1 whereas it is stronger than perturbativity for ξ * < 6 in M2. Finally for δ = 0.6, backreaction is stronger than perturbativity for ξ * < 4.77 in M1 while this range is reduced further to ξ * < 4.42 in M2. Considering the upper limits derived in (2.19), we can then compile all the backreaction and perturbativity constraints on the rolling spectator axion-models as where the coefficients b, p can be read from Tables 1 and 2 for all δ choices we focus and following the discussion we presented above. Focusing on the representative cases of δ = {0.2, 0.6}, in Fig. 2 we illustrate the parameter space f /M pl − ξ * consistent with perturbativity and backreaction bounds. Note that, using the relation (2.13) at fixed δ, the same parameter space can be described in terms of the dimensionless coupling λ of the spectator axion-gauge field interaction (2.2) (see the upper x-axis). Comparing the top and bottom panel plots, we confirm that within the same range of ξ * (λ) the non-compact axion model M2 has a larger parameter space where backreaction and perturbativity constraints are satisfied. Furthermore, for a faster rolling spectator axion χ (larger δ), a larger portion of the parameter space opens up for both models 11 . Finally, for smaller r v allowed region for f /M pl enlarges at fixed ξ * as can be also inferred from the expressions we derived in (2.26). The dots indicated by red color in the graphs locate the point in the parameter space where the lower limits derived from perturbativity considerations become comparable to backreaction constraints. In accordance with our discussion above, beyond this point, the dominant constraints on the lower bound for f /M pl come from P A in (2.25), resulting with a slight change in the slope of the limiting line (see p 1 vs b 1 from Tables  1 and 2). This change in the slope is barely recognizable from Fig. 2 since we are showing the constraints using a linear-log scaling. To sum up our findings in this section, the energy density contained in the spectator axion sector is approximately given by ρ χ ≈ H 2 f 2 (see eq. (2.20)). Therefore, the requirement (2.19) on the sub-dominance of axion's energy density with respect to the inflaton sector, ρ φ 3H 2 M 2 pl , give rise to an upper bound (2.20) on the scale f with respect to M pl . On the other hand, backreaction and perturbativity criteria we discussed above can be translated into a lower bound on the scale f (see eqs. (2.22) and (2.25)). This is because lowering the scale f with respect to M pl leads to a stronger interaction (2.2) between spectator axion and U(1) gauge fields, increasing the efficiency of the particle production in the gauge field sector. At fixed δ and vacuum tensor-to-scalar ratio r v , we showed that this lower bound is dictated by the effective coupling ξ * (λ) between the spectator Abelian gauge and axion field and parametrized by the first inequality in eq. (2.26).

New constraints on scale-dependent GWs sourced by vector fields from
Planck and BICEP/Keck data While the previous section was dedicated to the evaluation of theoretical bounds given by the backreaction and perturbativity considerations, we will devote this section to provide observational constraints on the parameter space of the spectator axion-U(1) gauge field model (2.1), and specifically on the transiently rolling axion models (see Section 2.1) M1 and M2 parametrized by the scalar potentials in (2.11).
CMB 2-point function analysis. For this purpose, we use the latest Planck NPIPEprocessed PR4 maps release [9] and BICEP/Keck 2018 (hereafter BK18) data [10], which together represent the state-of-the-art dataset for constraining primordial scalar fluctuations and tensor modes at the largest cosmological scales we are interested in this paper 12 . Similarly to analysis carried in [25] using a simplified data analysis setup (i.e. by fixing all model's, cosmological and likelihood nuisance parameters except for the ξ * parameter) with WMAP temperature data, by exploiting the full power of present day CMB temperature and polarization datasets, we aim to show the extent to which the models considered in this work can produce an observable amount of sourced gravitational waves, while remaining consistent with current tight constraints on the scalar sector.
In the following, we perform a likelihood analysis leaving the ξ * parameter free (together with cosmological and likelihood nuisance parameters) and fixing δ and k * . As motivated in Section 2.1, for δ, we consider the following representative values δ = {0.2, 0.3, 0.4, 0.6}, ordered according to decreasing amplitude (at fixed ξ * ) and increasing sharpness of the sourced Gaussian bump (2.18). For k * , we choose instead the following three representative values: For the models under consideration, the first value above typically induces a sourced bump in the scalar and tensor spectra at the very largest CMB scales (pertaining to the reionization bump in E and B-modes spectra), the second one affects the recombination bump's multipole range in B-modes, and finally the third one impacts scales around the first acoustic peak again in B-modes. Note from eq. (2.18) that the actual scale at which the sourced bump appears is different from the critical scale k * : k p = k * x c 2,− [δ, ξ * ] > k * because of the momentum conservation law of one-loop interactions that generates the sourced signals 13 . This deviation becomes larger at fixed ξ * for increasing δ, as can be inferred from Tables 3-6. Furthermore, k p is typically larger for the model M2 compared to M1 (at fixed k * , ξ * and δ): we will see in the following that this has important consequences on the model constraints.
CMB 3-point function and parity-violating correlations. As highlighted in the previous literature [25,72], the strict constraints on scalar non-Gaussianity at smaller scales can be evaded in the models under examination, if we consider a sourced bump at large scales in the 12 Note that in principle, the M2 model can also simultaneously produce sizeable amount of GWs at sub-CMB scales, which can be probed for instance through pulsar timing arrays and laser interferometers [28]. 13 In particular, at the time when axion's velocity peaks (τ * = k −1 * ), maximally amplified gauge field modes (running in the loop) are typically inside the horizon obeying q > k * O(1) [25,28]. Due to momentum conservation at each vertex (2 × (Ai + Ai → {h λ , R})) contributing to the one-loop power spectra, the resulting correlations among the external states ({h λ , R}) is therefore maximal for wave-numbers satisfying k = kp > k * . spectra, generated by an axion rolling for only a few e-folds ∆N χ ∼ δ −1 during inflation. For example, this condition is realized if the axion velocity satisfiesχ → 0 (or ξ → 0 as in (2.12)) when modes with O(10 2 ) leave the horizon so that gauge field production is ineffective at those scales. Tensor non-Gaussianity 14 [26] generated in spectator axion-U(1) model can also provide complementary information to CMB 2-point functions, even though the most stringent constraints are still obtained from the latter. The sourced bump in BB is indeed accompanied by a similar bump in the BBB 3-point function, which however has smaller signal-to-noise ratio compared to BB one [72]. Therefore, the analysis in the following will be based solely on CMB 2-point functions, and we leave a thorough analysis including bispectrum constraints to future work. Finally, as we discussed in Section 2.1, the models we consider produce fully chiral gravitational waves: the possibility of detecting such circular polarization with the CMB has been addressed in previous literature [25,90]. However, parity-violating EB and T B correlations in Planck data can constrain only very weakly the chirality parameter [90], therefore we will not consider them further in our analysis.

Data and likelihoods
As anticipated above, we exploit the latest Planck and BK18 public data releases, with likelihoods publicly available for the Cobaya [91] MCMC framework. Specifically, we combine the low-T T Commander likelihood (covering multipoles = 2 − 30) with the high-T T + T E + EE HiLLiPoP likelihood in the range = 30 − 2500 and the low-( = 2 − 150) EE + BB + EB LoLLiPoP likelihood, as described in [9]. In the analysis, we also include the B-mode intermediate-scale constraints from BK18, neglecting correlations with Planck. This is doable because B-modes are noise-dominated and the two CMB surveys have uncorrelated noises and, moreover, they observe very different fractions of sky [8,9].
In the following, we will also study the impact of separate datasets on the constraints: in particular, we name Planck TT the combination of low-T T Commander and HiLLiPoP T T likelihoods, while Planck TEB is the combination of Planck TT, HiLLiPoP T E + EE and LoLLiPoP EE, BB and EB likelihoods.

Methodology: the profile likelihood
In order to provide constraints on the model parameters, we perform a frequentist profile likelihood analysis. Compared to the Bayesian framework widely used in cosmology, frequentist methods have been applied in fewer occasions, despite having several advantages [80,92]. First, they do not require to choose arbitrary priors on the parameters, a practice which may have an important impact on the final bounds in a Bayesian setting. Second, the maximum likelihood estimate (MLE) is invariant under different choices of the model parameterization. Third, frequentist parameter estimates are not affected by so-called "volume effects" [81], which can instead appear, due to the marginalization process, in Monte Carlo Markov Chain analysis 15 .
The profile likelihood for a given parameter of interest θ is obtained by fixing θ to a certain value within the range of interest and maximizing the likelihood with respect to all remaining parameters. The maximisation is then repeated for several different values of 14 See also [75] for the study of mixed scalar-tensor type non-Gaussianity in the spectator axion-gauge field models we consider in this work. 15 Volume effects arise because marginalization enhances regions of the parameter space that contain more probability density volume in the marginalised directions. Moreover the volume of probability density in a certain parameter direction depends on both the choice of priors and the model parameterization. This can often result into the peak of the marginalized posterior distribution being far from the global MLE. the parameter of interest, scanning a wide range of θ values. The minimum of the profile likelihood built in this way coincides, by construction, with the global MLE given the full parameters set. Specifically, we minimize the χ 2 function using the iMinuit multi-dimensional minimizer package [93], a python implementation of the famous Minuit algorithm [94]. A typical example of profile likelihood for our parameter of interest ξ * is shown in Fig. 3 for the M1 model where we focus on δ = 0.6 and k * = 5 × 10 −3 Mpc −1 for four representative values of the vacuum tensor-to-scalar ratio 16 r v = {0.0001, 0.001, 0.01, 0.044}.
The behaviour of the profile likelihood reflects the exponential nature of the gauge field production [86], with a steep growth starting at increasing ξ * value for decreasing r v . ∆χ 2 is instead zero for smaller values of ξ * , since the amount of sourced modes produced by the model is negligible and does not affect the likelihood.
An upper bound on ξ * ≡ ξ * ,limit is obtained by cutting the profile likelihood ∆χ 2 = χ 2 − χ 2 min = 4 in each of the cases considered. We note that the limits derived in this paper cannot be directly compared to the ones derived in [25], because in the latter all parameters were fixed to WMAP ΛCDM best-fit values except ξ * . Therefore the approach in [25] is not guaranteed to reach the global MLE of the likelihood, while in the profile likelihood approach used in this paper we vary all parameters (model + cosmological + nuisance) and the result matches the global minimum of the likelihood up to numerical accuracy 17 .

Observational bounds from the CMB: results and discussion
The upper bounds on the ξ * parameter, obtained from the latest Planck and BK18 datasets, are summarized in Fig. 4 for both the M1 and M2 models. The impact of separate datasets (i.e. Planck TT, Planck TEB and Planck + BK18) on the constraints is singled out in Fig. 5 for M1 and Fig. 6 for M2, considering the two representative values δ = 0.2, 0.6. Finally, in Figures 7 and 8, respectively, we show the theoretical CMB spectra and the total (vacuum + 16 We also attempted building the profile likelihood by fitting rv in addition to ξ * , the cosmological and the nuisance parameters. However, because of the degeneracy between rv = 16 φ and ξ * (2.18), the latter remains essentially unconstrained when fitted together with rv. This happens because both ξ * and φ control the amplitude of the sourced signals, so it is always possible to decrease rv to accommodate for larger ξ * . Therefore, we fixed rv to phenomenologically reasonable values in order to obtain more informative constraints. 17 We also checked that the minimizer was not trapped in any local minimum, by starting minimization from a wide range of different initial parameters sets. sourced) primordial tensor power spectra P h (k) evaluated at ξ * , limit for some representative cases 18 . We discuss the bounds for each model separately below, starting from the M1 model.
Bounds on M1 model. For every δ considered, the upper limit on ξ * becomes tighter as the sourced bump moves from larger to smaller scales (Fig. 4, left panel): the reason is that most of the constraining power is coming from scalar modes sourced in the T T and EE spectra.
Planck large-scale T T modes are indeed cosmic variance-limited, so the trend can be imputed mainly to decreasing cosmic variance at smaller scales. Also, ξ * , limit is tighter for smaller δ at fixed k * : the width of the sourced bump is indeed proportional to 1/δ, as sourced modes are produced only while the axion is significantly rolling (i.e. for a number of e-folds ∆N χ 1/δ [25]). The sharper the bump is, the fewer multipoles are affected, and so the constraints on ξ * will be generally weaker. In addition, decreasing ∆N χ reduces more the production of sourced scalars than that of tensors: the process δA + δA → δχ → δφ is indeed very sensitive to the axion's velocity [25] and therefore to ∆N χ . We now compare constraints from Planck TT, Planck TEB and Planck + BK18. Figure  5 confirms that ξ * , limit is governed by temperature data: limits do not improve significantly when adding Planck E and B-modes, except when the bump is sourced at the first acoustic peak scales (i.e. for k * = 5 × 10 −3 Mpc −1 ). In this case, since the addition of BK18 B-mode data has no significant effect, we conclude that Planck intermediate/small scale E modes are providing the extra constraining power. Adding BK18 data has no effect for wider bumps, but can slightly tighten the upper bound in the case δ = 0.6, when sourced tensors are produced around recombination bump scales, since BK18 is sensitive only to multipoles 30 − 250 (Fig. 7).  Bounds on M2 model. Let us now discuss the constraints on the M2 model ( Figure 4, right panel). The ξ * , limit allowed for M2 is larger than the one for M1 in all cases considered: the axion's velocity profile in the M2 model (2.12) is indeed sharper than the M1 one (see also Fig. 8), and therefore allows for larger sourced tensors production for the same level of sourced scalars [28]. While in the M1 case ξ * , limit is always tighter at larger k * and larger δ, for M2 this holds only for δ < 0.6: similarly to the M1 model, scalars drive the constraints Here we assume r v = 10 −4 . We also report as reference the 95% C.L. error bars and upper limits from Planck (in gray) and BK18 (in red) data.
for wider bumps, while for δ = 0.6 tensors play a crucial role. For a bump at the largest and smallest scales considered (i.e. k * = 7 × 10 −5 and 5 × 10 −3 Mpc −1 ), indeed, constraints do not improve significantly when adding polarization data (Fig. 6), and are again primarily driven by temperature data. At intermediate scales (k * = 5 × 10 −4 Mpc −1 ), instead, constraints substantially improve when adding Planck E and B data and even more when adding BK18 data, confirming that tensor modes are driving the constraints. Furthermore, the peak of the signal moves to larger k p for M2 compared to M1 at fixed k * , as can be seen in Fig. 7 and 8. This is exacerbated at large δ and contributes to the loose ξ * , limit allowed for δ = 0.6. To summarize our findings, when the bump is sourced at recombination bump scales or smaller, adding polarization enhances the constraining power on both models considered. Moreover, the difference between the Planck TT and Planck TEB upper limits at fixed k * is larger for larger δ, since for a sharper bump the scalar constraints allow for a larger value of ξ * , increasing the SNR in polarization (especially in BB) and therefore making polarization data more relevant. On the other hand, for a sourced signal peaking at the largest scales considered, adding polarization data has minor or no impact on both models constraints, compatibly with the low sensitivity of Planck B-modes at the very largest scales [8,9].  Total tensor-to-scalar ratio. Finally, in Fig. 9, we show the values of the total (i.e. vacuum + sourced) tensor-to-scalar ratio r * (k) [25], evaluated at ξ * , limit for both models. All spectra appearing in (3.1) are evaluated at the peak of the sourced signal k = k p = k * x c 2,− [δ, ξ * ]. Compatibly with our previous discussion, sourced signals peaking at the largest scales generally allow for larger r * values 19 . Interestingly, the δ = 0.4 case allows for the highest r * at the largest scales for both models: this is because it represents a good compromise between a signal not so spiky that it cannot compensate for the correct normalization of the total scalar power spectrum with sourced modes, and one that is spiky enough that it allows for large values of ξ * . The production of a sizeable amount of sourced tensor modes, while still complying with scalar constraints, is thus realized. Moreover, Fig. 9 highlights the fact that it is still possible to get significant contribution to r * from sourced modes (reaching r * ∼ O(10 −2 )) in the k * = 5 × 10 −4 Mpc −1 case, even with a vacuum contribution as small as r v = 10 −4 or 10 −3 . On the other hand for k * = 5 × 10 −3 Mpc −1 , the allowed sourced contribution is smaller but can still be significant, especially for the second model M2.
Summary of observational constraints. To summarize the content of this section, using the latest Planck and BK18 CMB datasets, we derived, for the first time in the literature, constraints on the effective coupling parameter ξ * of the spectator axion-U (1) gauge field model in (2.1) for two possible potentials of the transiently rolling spectator axion in (2.11). We used a fully frequentist profile likelihood approach to derive upper bounds on ξ * which are  independent from prior distributions and model parametrization choices and thus immune to volume effects. We provided a detailed interpretation of the behaviour of the upper bound ξ * , limit for different choices of the δ and k * model parameters and compared these results for the two axion potentials under consideration. In conclusion, as can be seen by comparing Figures 2 and 4, the observational bounds reported in this section are competitive with the theoretical bounds from perturbativity and backreaction (Section 2). We will address in the next section the effect of these combined theoretical and observational bounds on the model's parameter space, in the context of particle production in the spectator sector.

Conclusions
At the end of this decade, new CMB probes, such as the LiteBIRD satellite [13] and the ground-based CMB-S4 [14], will target the imprint in the B-mode polarization pattern left by the primordial gravitational waves. In case of a detection, however, it will still be necessary to perform further tests in order to understand the origin of this signal and distinguish between the SGWB generated by quantum vacuum fluctuations of the metric, within the leading paradigm of single-field, slow-roll inflation, and the one possibly sourced by additional matter fields present during inflation. The SGWB properties predicted in these two scenarios can greatly differ, e.g. an almost scale-invariant spectrum from quantum vacuum fluctuations versus a strongly scale-dependent one when matter fields intervene.
In this paper, we relied specifically on the axion-U(1) gauge field model (2.1) for sourcing gravitational waves: this model involves, in addition to the usual scalar field driving inflation, a spectator sector including a gauge field with U(1) symmetry directly coupled to an axion. We considered two choices for the rolling axion potential (M1 and M2 in (2.11)), both capable of giving localized gauge field amplification at large/intermediate CMB scales. In Section 2.2, we provided bounds on the parameter space of the model, and more specifically on the effective coupling ξ * (and λ) between the axion and the gauge field, as implied by self-consistency of the theory, i.e. validity of the perturbative regime and negligible backreaction from the gauge field quanta. The theoretical bounds are summarized in (2.26) and the resultant available parameter space is shown in Fig. 2.
In Section 3, we completed the analysis of the model by deriving upper bounds on ξ * from state-of-the-art CMB spectra, namely from the latest Planck and BICEP/Keck data. We adopt for this purpose the frequentist profile likelihood approach, fully exploiting in this context its immunity to prior choices, model parametrization and volume effects, which instead are known to affect Bayesian estimates. We summarize in Fig. 4 the upper bounds on ξ * from Planck and BICEP/Keck data, for typical choices of the model parameters δ and k * which control, respectively, the width and the position in wavenumber space of the bump feature sourced by gauge fields in scalar and tensor spectra.
Reduced viable parameter space by Planck and BK18. The observational upper limits we obtained on ξ * can further tighten the parameter space of the model (2.1), and are competitive with the theoretical bounds presented in Section 2.2 for the self-consistency of our approach. In particular, as discussed in Section 3.3, the Planck and BICEP/Keck data limits the size of the effective coupling ξ * and hence the height of the maximally chiral, scale-dependent tensor perturbations sourced by the gauge fields at scales k * relevant for CMB observations. Including these observational bounds leads to a further reduction of the available parameter space consistent with backreaction and perturbativity bounds, in the f /M pl − ξ * (λ) plane for the spectator axion gauge field model. For δ = 0.3, 0.4, 0.6, we superimpose these observational constraints with the theoretical bounds (see Fig. 2) and show the resulting parameter space in Fig. 10. Comparing with the perturbativity + backreaction constraints presented in Fig. 2, we can clearly observe that the available parameter space shrinks from a large triangle at fixed r v to a smaller right trapezoid by the observational constraints on ξ * (i.e. ξ * ,limit ), shown by the vertical lines corresponding to bounds obtained at k * = 7 × 10 −5 (solid), k * = 5 × 10 −4 (dot-dashed) and k * = 5 × 10 −3 (dotted) Mpc −1 . As can be confirmed from these plots and from our discussion in the previous section, tightest limits on the area of available parameter space of the models (2.1) come from the smallest scales 20 k * = 5 × 10 −3 Mpc −1 while the allowed region gradually enlarges towards larger scales at fixed r v . Similar to our discussion in Section 2.2, for choices of r v smaller than what is presented in Fig. 10, a larger parameter space can in principle be made available. It is worth stressing that such cases correspond to inflaton sectors endowed with flatter scalar potential Within the available parameter space presented in Figure 10, the total tensor-to-scalar ratio r * can be inferred from Figure 9. In particular, 20 The only exception is the second Model (M2) with δ = 0.6, for which the physical peak kp of the sourced signals occurs at scales where observational constraints by Planck and BK18 are weak (see section 3.3 and footnote 19). For this reason, in Figure 10 we did not include the observational constraints on M2 with δ = 0.6 at k * = 5 × 10 −3 Mpc −1 (see the bottom right panel). we clearly observe that a sizeable sourced contribution to r * from gauge fields is viable while complying with CMB data at all k * values we consider.  and on small scales with ground-based experiment (e.g. CMB-S4). The green and purpleshaded areas represent the range of multipoles to which each kind of experiment is typically sensitive to. Dotted black lines show the single-field slow-roll prediction for r = 0.0046 (predicted by the Starobinsky model [2,95]) and for r = 0.001 (close to detectability limit for future B-mode probes). Solid lines show theoretical CMB spectra candidates evaluated at ξ * , limit . Specifically, Model A (M2 model: k * = 7 × 10 −5 Mpc −1 , r v = 0.001 and δ = 0.4, solid purple) is indistinguishable at > 30 from the vacuum prediction for r = 0.001, while featuring a very distinctive reionization bump. Model B (M1 model: k * = 5 × 10 −4 Mpc −1 , r v = 0.001, δ = 0.2, light blue) features instead a recombination bump very similar to the one for r = 0.0046, but with a quite different behavior from the vacuum prediction at reionization bump and smaller scales. Finally, Model C (M1 model: k * = 5 × 10 −3 Mpc −1 , r v = 0.001, δ = 0.6, light orange) is indistinguishable from r = 0.001 in the whole multipole range accessible from space, but has a distinctive bump feature at > 200, making a ground-based mission necessary to distinguish it from the standard slow-roll prediction. For reference, we also show the 95% C.L. error bars and upper limits from Planck (in gray) and BK18 (in red) data.
The path ahead: relevance of a B-mode satellite mission. As we discussed in Section 3.3, the current observational constraints on the model (2.1) are mainly driven by temperature spectra (i.e. from the sourced contribution to scalar fluctuations). B-mode polarization data at large/intermediate scales from Planck and BICEP/Keck are indeed weakly constrained, and therefore have a minor effect on the model bounds ( Figures 5 and 6). Large-scale temperature data are already cosmic variance-limited in the Planck dataset, so sensitivity to (large/intermediate scale) polarization must be improved to better constrain the axion-U(1) gauge field model. In particular, we argue that a B-mode satellite mission with access to the large and intermediate CMB scales (i.e. the ones pertaining to the reionization bump), such as LiteBIRD, would have unique benefits in distinguishing a vacuum-generated SGWB from a sourced one in the model under consideration [96,97]. Ground-based experiments, such as BICEP/Keck considered in this paper or the planned high-sensitivity CMB-S4, indeed, cannot access the multipoles 30, as they typically have much smaller sky coverage, compared to the almost full-sky measurements available from space, and are affected by Earth's atmospheric contamination at the largest scales. Nonetheless, future ground-based experiments would still be highly beneficial and complementary to a satellite mission: high-sensitivity measurements of intermediate scales B-modes would help discriminating between vacuum and sourced origins of GWs for k * 5 × 10 −4 Mpc −1 . In Fig. 11, we illustrate the potential benefits of measuring both large and small scale B-modes, with a space mission and a ground-based experiment, respectively. A full-sky space mission would also be necessary to obtain improved measurements of EB and T B parity-violating correlations at the largest scales, which are non-vanishing for gauge-sourced SGWB production during inflation, as discussed in Section 1. Furthermore, LiteBIRD will also greatly improve limits on tensor non-Gaussianity at large scales, making a signal of order O(1) potentially detectable in tensor-tensor-tensor equilateral (i.e. f ttt,eq N L ) and squeezed (f ttt,sq N L ) configurations [72].