Gravitational Waves from Fundamental Axion Dynamics

A totally asymptotically free QCD axion model, where all couplings flow to zero in the infinite energy limit, was recently formulated. A very interesting feature of this fundamental theory is the ability to predict some low-energy observables, like the masses of the extra fermions and scalars. Here we find and investigate a region of the parameter space where the Peccei-Quinn (PQ) symmetry is broken quantum mechanically through the Coleman-Weinberg mechanism. This results in an even more predictive framework: the axion sector features only two independent parameters (the PQ symmetry breaking scale and the QCD gauge coupling). In particular, we show that the PQ phase transition is strongly first order and can produce gravitational waves within the reach of future detectors. The predictivity of the model leads to specific characteristics of the phase transition (like its duration and the nucleation temperature) and the gravitational wave spectrum.


Introduction
The PQ symmetry [1] provides an elegant solution of the strong CP problem: it explains dynamically why the strong interactions preserve CP, while the electroweak ones break it. This solution, regardless of its implementation, manifests itself at low energies through the axion, the pseudo-Goldstone boson associated with the breaking of this (approximate) symmetry [2]. Moreover, the axion is a good dark matter (DM) candidate and can actually account for the whole DM when the PQ symmetry breaking scale f a is around 10 11 GeV. Even if there may be other components of DM, astrophysical observations require anyhow f a to be above ∼ 10 8 GeV (see Ref. [3] for a recent review). This makes testing the PQ idea difficult as current colliders are far from probing those energies.
Cosmology, on the other hand, allows us to have a window on physical processes occurring at such a high mass scale. A classic example is the observation of the cosmic microwave background (CMB) anisotropies. Indeed, these give us information on the inflationary dynamics, which typically occurs at energies much above those within the reach of colliders (see Refs. [4,5] for up-to-date observations). Recently, the experimental discovery of gravitational waves (GWs) [6] has opened another window on very energetic processes predicted by particle physics models. An important example is the possibility to have experimental information on the characteristics of phase transitions, if these are strongly first order (see Ref. [7] for a textbook introduction to this topic). Moreover, the possible observation of GWs due to a first-order phase transition would be a remarkably clear signal of new physics as the finite-temperature symmetry breaking dynamics in the Standard Model (SM) is not of this type.
If the PQ phase transition is strongly first order it could lead to observable GWs. But, as pointed out in Refs. [8,9], the particular features of the phase transition and the corresponding GWs depend on the specific dynamics implementing the PQ symmetry 1 . These features include, for example, the temperature at which the phase transition occurs, its duration and the GW spectrum.
One way one can significantly reduce this ambiguity, and thus have a guide regarding the directions towards which the experimental efforts may be focused, is considering only fundamental, that is truly UV complete, QCD axion models. This is because fundamental field theories, such as asymptotically free [12] or asymptotically safe [13] models, have typically the ability to predict the value of some observables (see e.g. Ref. [14]).
Recently, the first fully calculable and realistic QCD axion model of this sort has been constructed [15]. It features a new non-Abelian interaction and new quarks and scalars. This model is valid up to infinite energy because is totally asymptotically free (TAF): all couplings can flow to zero in the infinite energy limit. The requirement of total asymptotic freedom is a sufficient but non-necessary feature of a fundamental theory, as asymptotic safety can be an alternative (see [16][17][18][19] for phenomenological applications). However, having all couplings approach zero in the UV allows us to trust perturbation theory, at least for sufficiently high energies, and thus obtain a calculable model.
This asymptotically free axion model, as is typically the case in fundamental field theories, predicts some low energy observables: this is because some couplings (specifically the Yukawa and quartic couplings) are compatible with the TAF requirement only if they acquire some specific isolated values at low energy. In other worlds, these couplings are IR attractive.
The requirement of having a classically scale-invariant theory, with no dimensionful parameters in the classical Lagrangian, further reduce the number of independent parameters. Indeed, the observed mass scales are then obtained quantum mechanically rather than with additional mass parameters in the classical Lagrangian. For this reason we investigate here a region of the parameter space of the TAF axion model [15] where the PQ symmetry is broken quantum mechanically through the Coleman-Weinberg (CW) mechanism [20]. This mechanism allows us, among other things, to generate scales via quantum corrections in the regime of validity of perturbation theory and, therefore, have a fully calculable setup.
The main purpose of this work is to study the PQ phase transition and the characteristics of the possible associated GWs in this highly predictive framework. As pointed out in Ref. [21] (see [8,9] for recent applications to QCD axion models) the phase transitions associated with the CW symmetry breaking are typically of first order and can, therefore, lead to observable GWs 2 . This, together with the high predictivity of the TAF axion model mentioned above, can lead to testable implications at GW detectors. Furthermore, a typical feature of CW phase transitions is the presence of a phase of strong supercooling, meaning that the temperature below which the phase transition is effective (the nucleation temperature) turns out to be much smaller than the critical temperature [21] and, in the case of the PQ symmetry breaking, much below f a [8,9].
Since the TAF requirement regards the extrapolation of the theory at arbitrarily high energies, some comments about the behavior of gravity in the UV are now in order. Here we assume that the gravitational interactions are softened (compared to their behavior in Einstein's theory) above and only above a certain energy scale Λ G [14]. While at large lengths all successes of Einstein's theory are reproduced, gravity is assumed to be so weak from the UV down to the PQ scale that its impact on the renormalization group equations (RGEs) can be neglected. This is possible because Λ G can be much below the Planck scaleM P , where quantum gravity effects in Einstein gravity 2 See also Ref. [22] for a previous study of GWs and phase transitions in models with CW symmetry breaking. would become sizeable. Since a phase of strong supercooling occurs in our classically scaleinvariant TAF axion model, the nucleation temperature is much smaller than f a and the relevant values of the fields and their derivatives are also much smaller than f a . This implies that gravity, as far as the production of gravitational waves is concerned, is well-described by Einstein's theory for all values of Λ G satisfying f a Λ G <<M P .
This softened-gravity scenario may be realized, for example, in UV modifications of gravity featuring quadratic curvature terms in the action [23,24] or in non-local extensions of general relativity [25]. Interestingly, the former case also admits a classically scale-invariant formulation, called Agravity [23], in which the Planck and the Fermi scales as well as the cosmological constant are generated quantum mechanically via a gravitational version of the CW mechanism.
The paper is organized as follows. In Sec. 2 we introduce the TAF axion model without a fundamental PQ scale and review results regarding the RG flow in [15] that are most relevant for our purposes. In Sec. 3 it is shown that f a is generated quantum mechanically through the CW mechanism and the mass spectrum is also discussed. In the same section we also show that the axion sector only features one independent mass scale, i.e. f a , and one independent dimensionless parameter (the QCD gauge coupling evaluated at f a ). The PQ phase transition is then studied in Sec. 4, which gives details about the finite-temperature effective potential, the bounce solutions and their actions, the nucleation temperature and the reheating after the supercooling era. In Sec. 5 the predictions of the classically scale-invariant TAF axion model regarding the GW spectrum are worked out. A description of the relevant GW detectors is then provided in Sec. 6, where we also compare the sensitivities of the GW experiments with the predictions of our fundamental model. Finally, in Sec. 7 we provide our conclusions.

A fundamental axion model without a fundamental PQ scale
Here we consider a dimensionless version of the TAF axion model of [15]: we take the limit in which the dimensionful parameters in the microscopic Lagrangian of the TAF axion model of [15] go to zero. The axion sector is invariant under an SU(2) group (henceforth SU(2) a ). Then the full gauge group contains the factor SU(3) c ×SU(2) a , where SU(3) c is the ordinary QCD group. The gauge group should also include extra factors to account for a TAF extension of the SM (see e.g. [14,26,27]). We call such an extension the "SM sector", which of course has to be present, in addition to the axion sector we describe here in order for the complete model to be fully viable. We will not commit ourselves to a specific TAF SM extension here, but we note that in general the SM and axion sectors interact via SU(3) c gauge interactions.
The model features two extra Weyl fermions q andq in the fundamental and antifundamental of SU(3) c ×SU(2) a , with the same PQ charge: {q,q} → e iα/2 {q,q}, where α is a constant. The PQ charges of all particles in the SM sector vanish for simplicity, like in KSVZ-like axion models [28]. We introduce a scalar field A, which spontaneously breaks the PQ symmetry (denoted here U(1) PQ ) and gives mass to the extra quarks (as required by the experiments). Therefore, A is complex and have Yukawa interactions with q andq, Note that V A , just like any other term in the Lagrangian, only contains dimensionless coefficients. The possible mass term m 2 Tr(A † A) in the potential in [15] has been erased or, more generally speaking, it will be assumed that m is much smaller than the effective mass generated by the CW mechanism (we will discuss how this mechanism works in the present model in Sec. 3). The necessary and sufficient conditions for vacuum stability at high-field values (henceforth "high-field stability") are [15] Here we neglect the couplings with the scalars of the SM sector; setting those couplings exactly to zero is consistent at one-loop level. At higher-loop level those couplings can be generated, but they remain small and will have negligible effects on our results. Let us now review the beta-functions of this model [15]. The renormalization group equation (RGE) of the gauge coupling g of a generic gauge group is where t ≡ ln(µ 2 /µ 2 0 )/(4π) 2 , the quantity µ 0 is an arbitrary reference energy and µ is the usual RG scale. The solution to Eq. (2.4) is UV attractive for any g 0 ≡ g(0) and AF requires b > 0. The value g = 0 is a trivial fixed point of Eq. (2.4) and so we take g > 0 without loss of generality. For SU(2) a and SU(3) c the constant b for the corresponding gauge couplings g a and g s reads, respectively, where ∆ is the positive extra contribution due to the fermions and scalars in the SM sector. In the numerical calculation we will use for definiteness the reference value ∆ = 28/3, which is compatible with known TAF SM extensions [15,27]. The RGE of y is instead Like for the gauge couplings, the value y = 0 is a trivial fixed point of Eq. (2.6) and so we take y > 0 without loss of generality. This equation admits a closed-form solution for any b a and b s [15]. Finally, the RGEs of λ 1 and λ 2 are dλ 1 dt = β 1 , and dλ 2 dt = β 2 , where [15] β 1 (g a , y, λ) = 9 2 g 4 a + λ 1 8λ 2 + 6y 2 − 12g 2 a + 14λ 2 1 + 8λ 2 2 − 3y 4 (2.7) and Note that λ 1 = 0 and λ 2 = 0 generically are not zeros of β 1 and β 2 and, therefore, the quartic couplings, or some combinations thereof, could a priori change sign and pass through zero during the RG evolution.
In [15] the system of equations (2.4)-(2.8) was solved and it was found that for any initial condition for the gauge couplings there is one and only one TAF solution satisfying the stability conditions in (2.3) at high-field values. For such solution both y and λ i are IR attractive and are, therefore, predicted at low energies [15]. We will pick up this TAF solution from now on as the high-field stability conditions in (2.3) are necessary to have a viable setup.

Quantum generation of f a
The RGEs dictate that the couplings run with energy. The conditions in (2.3) are necessary at high energy for high-field stability, but they can be violated in the IR. In the present model we find that, while the first condition is always preserved, the second one is violated at small energy. This leads to spontaneous symmetry breaking of U(1) PQ through the CW mechanism [20], because at the energy scale µ PQ where λ ≡ λ 1 + λ 2 = 0 the effective potential develops a flat direction (corresponding to A = A † , such that the two terms in the potential (2.2) become equal: Tr(A † A) = Tr(AA)). We interpret µ PQ as the PQ scale. On the flat direction A = A † the three components A k of A along the Pauli matrices, A = A k σ k /2, can always be transformed through an element of SU(2) a in a way that only one of these components is not vanishing and positive. We call this non-zero component φ. For example, a possible choice for φ is φ = |Re(A 1 )|. In [15] the case in which the explicit mass m is larger than µ PQ was considered. Here we focus instead on the opposite case m µ PQ , so the PQ symmetry breaking is entirely driven by the CW mechanism. A non-vanishing value of φ breaks SU(2) a down to a residual Abelian group U(1) a leading to a massless spin-1 particle (a dark photon), two spin-1 particles with equal mass (which can be described by one complex vector field), two degenerate Dirac fermions with mass two scalars with squared mass and two massless scalars (one is the axion, which as usual acquires a mass through quantum correction, and the other one corresponds to the flat direction). Note that M 2 S ≥ 0 when λ 1 ≥ λ 2 , which turns out to be satisfied at all scales, from the TAF requirement [15].
At one-loop the quantum potential at zero temperature along the flat direction is given by where V 0 is the tree-level potential, where λ is evaluated at the renormalization scale µ = µ 0 exp(8π 2 t), and V 1 is the quantum one-loop correction In this expression the sum over b runs over all bosons (with number of degrees of freedom n b ), that over f runs over all fermions (with number of degrees of freedom n f ), M b,f (φ) are the corresponding background-dependent masses (which, for our model, are given in Eqs. (3.1), (3.2) and (3.3)) and a b and a f are renormalization-scheme dependent quantities. It is understood that the coupling constants in V 1 are evaluated at the same renormalization scale, µ. This part of the potential can be computed explicitly by using the background-dependent masses given above.
whereβ is the beta-function of λ evaluated at µ PQ , namelȳ and the scale f a has been introduced in a way that the CW potential, Eq. (3.7), has its stationary point at φ = f a . We conventionally choose a renormalization scheme such that f a = exp(−1/4)µ PQ . Whenβ > 0 the stationary point at φ = f a corresponds to a minimum. We have numerically verified the positivity ofβ. Then U(1) PQ is spontaneously broken and φ acquires the VEV φ = f a and a mass squared M φ = β f a . The energy scale f a is, therefore, identified with the PQ symmetry breaking scale. The remaining mass spectrum besides M φ corresponds to the axion, the dark photon and the tree-level masses obtained by setting φ = f a in Eqs. (3.1), (3.2) and (3.3). Since f a is at least of order 10 8 GeV and the mass of the extra-quarks M Q (f a ) is around f a , the dark photon satisfies all present bounds, including those of cosmological nature [15]. Note that the arbitrariness of µ 0 tells us that µ PQ (and thus f a ) is a free parameter. The TAF axion sector features only another free parameter, which can be taken to beḡ s ≡ g s (t PQ ) (here and in the rest of the paper a bar indicates that a generic coupling is evaluated at t = t PQ ≡ ln(µ 2 PQ /µ 2 0 )/(4π) 2 ). Indeed, once the gauge couplings are chosen at µ = µ PQ the other couplings y,λ 1 andλ 2 are predicted [15] and one must consider a particular IR value of one of the gauge couplings, sayḡ a , to enforceλ 1 +λ 2 = 0, namely to have CW symmetry breaking. In Fig. 1 we give the couplingsḡ a ,ȳ andλ 1 = −λ 2 as functions ofḡ s to show that all dimensionless quantities in the axion sector are fixed onceḡ s is chosen. As a result, when the PQ symmetry is brokenà la CW one also obtains a prediction for the mass of the extra complex vector field (in addition to the predictions of the extra scalar and fermion masses when m is larger than µ PQ [15]) once f a and g s are chosen.

Peccei-Quinn phase transition
In order to investigate the nature of the Peccei-Quinn phase transition we take into account thermal corrections as well as quantum corrections. We consider the one-loop effective potential Figure 1: The couplings of the model as functions of the QCD gauge coupling at the PQ scale.
where the thermal correction V T to the effective potential is given by [29] (see also [30]) and we have included in V eff (φ, T ) a constant term Λ 0 to account for the observed value of the cosmological constant when φ is at the minimum. It is understood that the coupling constants in V T are evaluated at the same renormalization scale, µ, used in V CW .
Since the background-dependent squared masses are all non-negative V eff has a vanishing imaginary part. This is due to the fact that f a is generated quantum mechanically rather than through an explicit tachyonic scalar mass, which would unavoidably lead to a concave tree-level potential and thus to a complex effective potential for some field values. Therefore, the CW symmetry breaking supports the validity of perturbation theory: indeed, a non-negligible imaginary part (absent in the CW case) generically signals the breaking of the perturbative expansion. Further comments regarding the approximation used will be given below in this section.
In Fig. 2 (left plot) we show V eff as a function of φ for two values of the temperature: the critical temperature T c and T = 0. That figure shows that the transition is of first order. Although we use a fixed value ofḡ s in that figure, other choices of this parameter lead to the same qualitative situation. In the right plot of Fig. 2 we give the dimensionless quantity T c /f a as a function of the only dimensionless parameter of the axion sector,ḡ s .
The absolute minimum of the effective potential is at φ = 0 for T > T c , while, for T < T c , is at a non-vanishing temperature-dependent value. In the latter case the decay rate per unit  Here S d is the action evaluated at the O(d) bounce, which is defined as the solution of the differential problem where a prime denotes a derivative with respect to r. Also, R 4 is the size of the O(4) bounce. Note that S d evaluated at the O(d) bounce can be simplified by using the scaling arguments of [35] to obtain Using the expression above instead of the one in (4.5) makes numerical calculations easier because the derivative of the bounce does not appear in the action. We numerically checked that S 3 /T < S 4 , so Γ is dominated by the O(3) bounce. Therefore, the phase transition is essentially due to thermal effects rather than quantum effects. As an example, in Fig. 3 we give a plot of the O(3) and O(4) bounces for representative values of the QCD gauge coupling and the temperature below T c .
Some words on the approximation used are now in order. We note that the correction to the two-derivative term in the one-loop effective action is small as long as the temperature is small compared to the field values [36] characterising the bounce solution. We have checked that T is small compared to (within 10% of) the relevant field values for all numerical calculations performed in this work. This also implies that we are far from the high-temperature regime for which the perturbative expansion is known to break down [37]. Moreover, note that, generically, the higher-derivative corrections to the one-loop effective action are suppressed when the laplacian applied to the solution of interest (in this case the bounce) is small compared to the background dependent masses times that solution: this follows from the structure of equations of motion without higher derivatives. Since the largest couplings are of order one in our case, those higher-derivative corrections are small when dV eff dφ is small compared to φ 3 in the relevant range of r (where S 3 /T gets its dominant contribution). We numerically checked that this condition is also satisfied (at the 10% level). Therefore, our one-loop approximation for the effective action is reliable.
Also the gravitational corrections to the false vacuum decay are amply negligible in our case. This is because the typical scales of the bounce and the temperature are always below f a (as illustrated in Figs. 2 and 3) and the gravitational corrections are, therefore, suppressed by factors at least as small as f 2 a /M 2 P [38]. The latter quantity is tiny because f a is several orders of magnitude belowM P in order for the axion not to overproduce DM (see Ref. [3] for a review).
Like for other models with CW symmetry breaking [8,9], we find that when T goes below T c the scalar field φ is trapped in the false vacuum φ = 0 until T is much below T c , in other words the universe features a phase of strong supercooling. Then the energy density is dominated by the vacuum energy of φ and the universe grows exponentially like during inflation but with Hubble rate H I = β f 2 a /(4 √ 3M P ), whereM P is the reduced Planck mass that is defined in terms of the Planck mass M P byM P ≡ M P / √ 8π. The bubbles created are diluted by the expansion of the universe and they cannot collide until T reaches the nucleation temperature T n , which corresponds to the temperature when Γ/H 4 I ∼ 1 or, equivalently, using the fact that the decay is dominated by the O(3) bounce, In Fig. 4 (left plot) we give T n as a function ofḡ s for two physically interesting values of the PQ symmetry breaking scale: f a = 10 11 GeV (for which the axion can account for the whole DM) and f a = 10 9 GeV (for which the contribution of the axion to the DM abundance is negligible) 3 . We find that the equation that determines T n in (4.8) does not always admit a solution: there is a minimal value of the couplingḡ s below which there is no solution. As clear from Fig. 4, when the coupling goes to this minimal value T n becomes very small. By comparing the right plot of Fig. 2 with the left plot of Fig. 4 one can see that supercooling takes place, namely T n T c . In Fig. 4 we also give a plot of the bounce actions S 3 /T and S 4 evaluated at T n as functions ofḡ s . That plot shows, as mentioned above, that S 4 > S 3 /T and so the phase transition is dominated by thermal effects. In the figure we consider, as an example, f a = 10 11 GeV, but the other values of f a lead to the same qualitative behavior.
We also note that strong supercooling and the corresponding inflationary period efficiently dilute the density n(T ) of monopoles due to the breaking SU(2) a → U(1) a . In a strong first-order phase transition monopoles may be created by bubble collisions and well-known estimates [39] lead to where p is the probability that the scalar field configuration is topologically non trivial, C = 0.6/ g * (T n ) and g * (T ) is the effective number of relativistic species in thermal equilibrium at temperature T . Even for p ≈ 1, setting g * (T n ) of order 10 2 (a realistic setup given the existing TAF SM sectors) we find that the theoretical bound in (4.9) is amply compatible (unlike in grand unified theories without inflation) with the bound coming from the fact that the mass density of monopoles must not exceed the limit on the total mass density imposed by the observed Hubble constant and deceleration parameter [39]. Indeed, the latter bound is around n(T 0 )/T 3 0 10 eV/M m , where T 0 is today's temperature, M m is the monopole mass and M m ∼ 4πf a /ḡ a , n(T 0 )/T 3 0 n(T n )/T 3 n and the window 10 8 GeV f a 10 12 GeV have been used. The strength of the phase transition is measured, as usual, by the parameter α defined as the ratio between the free-energy density associated with the transition where ∆V eff ( φ , T ) ≡ V eff ( φ , T ) − V eff (0, T ), and the energy density of the thermal plasma. So (4.11) We find a very strong phase transition with α exceeding one by several orders of magnitude. At the end of supercooling the universe should be reheated. This occurs in general thanks to the unavoidable coupling between the axion sector and the SM sector due to gluons. The field φ couples at one loop to gluons through the extra quarks q andq so one gets an effective interaction where G µν is the gluon field strength. This leads to the following rate of the decay of φ into two gluons (4.13) We find Γ φ→gg H I so the reheating is approximately instantaneous. The reheating temperature due to this channel may be computed through But this formula is only valid if the radiation energy density ρ R is not exceeding the vacuum energy density ρ vac due to φ (because ρ vac represents the full energy budget of the system). If this condition is not satisfied we determine T RH as the maximal temperature compatible with ρ R ≤ ρ vac , leading to Our estimate of T RH agrees to very good accuracy with previous determinations [8]. We find a very high T RH . For example, setting g * ∼ 10 2 and {ḡ s , f a } ≈ {0.91, 1.2×10 9 GeV} we obtain T RH ∼ 10 8 GeV, while for {ḡ s , f a } ≈ {0.97, 10 11 GeV} we obtain T RH ∼ 10 10 GeV. Figure 5: The theoretical prediction in the (T n , β/H I ) plane compared with the sensitivity curves of BBO, CE, DECIGO, ET and advanced LIGO (see Sec. 6). The theoretical curves are produced by varyinḡ g s as in Fig. 4. The parameter space enclosed within the experimental curves represents detectable signals.
Finally, another important parameter is β, which measures the inverse of the duration of the phase transition and, for fast reheating, can be computed with the formula

Gravitational waves
When the temperature drops below T n GWs are produced. The dominant source of GWs are bubble collisions that take place in the vacuum. This is because in the era when T reaches T n the energy density is dominated since a long time by the vacuum energy density associated with φ, which leads to an exponential growth of the cosmological scale factor as we have seen. This inflationary behavior as usual dilutes preexisting matter and radiation and, therefore, we neglect the GW production due to turbulence and sound waves in the cosmic fluid [7]. From [40] we find the following GW spectrum due to vacuum bubble collisions where f peak is the red-shifted frequency peak today and is given by [40] f peak ≈ 3.79×10 2 β H(T RH ) T RH 10 10 GeV 100 g * (T RH ) 1/6 Hz. (5.2) Being the reheating almost instantaneous we can approximate the Hubble rate at T RH with its value H I . In Fig. 6 we give h 2 Ω GW as a function of f for some relevant values of the parameters and setting as an example g * = 2×10 2 : this is a reasonable value given that the SM sector has to be extended to satisfy the TAF requirement. In the same plot we also compare these theoretical findings with the experimental sensitivities of GW detectors (see Sec. 6).
Note that any cosmic source of GW background acts as an extra radiation component and, therefore, modifies the expansion rate of the Universe. This means that it is highly constrained by big-bang nucleosynthesis (BBN) measurements of primordial elements [41]. The measurement of the number of effective neutrinos N eff and the observational abundance of Dueterium and Helium, gives us the following bound [42,43]: which is shown explicitly in Fig. 6.

Gravitational wave detectors
GWs serve as a probe for early universe cosmology, and is particularly important for the era prior to BBN. While inaccessible in collider or in other terrestrial experiments due to the high energy scales involved, the strong first-order phase transition associated with the PQ symmetry breaking induces a stochastic GW source and, therefore, gives us a window to experimentally detect a PQ model. Moreover, observing the particular GW spectrum predicted by the TAF axion model may serve as a first indirect verification of the TAF principle. Although there may be other sources of stochastic GWs in the early universe like inflation (Refs. [44][45][46][47]) or unidentified binary black hole mergers (Ref. [48]) the spectrum of GW radiation produced by phase transitions is generically different (see Ref. [7,49] for reviews). For instance, the expected stochastic background from compact binary coalescences, such as binary neutron stars or black holes, has a different dependence on frequency than the one in (5.1), [50]. Various limits on astrophysical sources contributing as stochastic GW background exist [51]. In 2015, when the GW event (named GW150914) from binary back hole mergers was observed [6,52] and its contribution reanalyzed as a source of stochastic background (with 90 % C.L. statistical uncertainty, propagated from the local rate measurement, on the total background) the event was found to be in the frequency range of 0.01 Hz -100 Hz, and of amplitude within the reach of advanced LIGO (see Ref. [48] for details). Such range overlaps with the one of the PQ phase transition in the TAF axion model. This, however, is not a problem: due to its weaker frequency dependence, any multiple network of GW detectors, like LIGO [53,54], VIRGO [55], GEO600 [56], KAGRA [57], and LIGO-India [58], will be able to separate the astrophysical signal from other sources, like that from a cosmological PQ phase transition or other events taking place after cosmic inflation. In the TAF axion scenarios, future GW detectors will be able to probe the model as we will discuss in detail below. Figure 6: The black lines give the GW amplitude corresponding to the values of the TAF model as mentioned in the figure. The areas above the colored lines correspond to the projected sensitivities for various GW observatories as detailed in section 6. The colored region is excluded by the BBN bound discussed at the end of Sec. 5.
All experiments to detect GW has strain noise power spectrum (Ω noise = 6f 2 peak /f 3 S noise ). The quantity S noise (f ) consists of intrinsic noise from the instrument as well as other astrophysical confusion noise S noise (f ) = S ins (f ) + S gcn (f ) (see Ref. [59] for details) that varies from detector to detector. This leads to a signal-to-noise ratio R given by [60] where t obs denotes the experiment's observing time, N = 1 (N = 2) for experiments that perform an auto (cross) correlation measurement of the stochastic gravitational wave background and f min and f max are the minimum and maximum frequency accessible to the detector, respectively. The GW spectrum can be detected by ground-based interferometers; these include advanced LIGO in Hanford and Livingston [53,54], Cosmic Explorer (CE) [61,62], Einstein Telescope (ET) [63][64][65]). Moreover, additional information on the GW spectrum can be obtained through space-based interferometers (BBO [66][67][68], DECIGO [69,70], and LISA [71]).
There are astrophysical sources contributing as confusion noise (see e.g. [72]) but its impact becomes insignificant with time as one is able to subtract this due to information from the individual foreground sources, which increases in number. However, we consider no foreground contamination, therefore, the results can be treated as an upper limit of the experimental reach possible in future.
In Fig. 6, we depict the predicted GW spectra for our benchmark points along with the projected sensitivities of current and future interferometer experiments. Our benchmark points include physically interesting cases with f a ∼ 10 9 GeV and f a ∼ 10 11 GeV for different values ofḡ s . Each GW projected sensitivity is denoted with corresponding legends. Curves are drawn following peak-integrated sensitivity curves 4 from Refs. [59,73]. As clear from Fig. 6, Advanced LIGO, ET, CE, DECIGO and BBO will be potentially able to detect GWs produced by the PQ phase transition in the TAF axion model. Fig. 5 provides sensitivity plots in the (T n , β/H) plane for various experiments and compare them with the theoretical curves from the TAF axion model. We considered the physically interesting cases f a = 10 11 GeV and f a = 10 9 GeV. Fig. 5 shows that ET, CE, BBO and DECIGO will be able to test this scenario through measurements of T n and β/H.

Conclusions
The detection of GWs and many upcoming GWs detectors have recently reinforced the interest in phase transitions predicted by particle physics models. In this paper we have studied the PQ phase transition and the corresponding spectrum of GWs in a QCD axion model where all couplings flow to zero in the infinite energy limit (TAF property) and the PQ symmetry breaking scale f a is generated quantum mechanically through the CW mechanism. This fundamental (i.e. UV complete) model features an extra gauge group SU(2) a , which is spontaneously broken to an Abelian U(1) a subgroup. The low energy spectrum, therefore, includes a dark photon, which has previously been shown to be compatible with current bounds from particle physics experiments and cosmology [15]. This TAF QCD axion model is highly predictive; indeed, the axion sector has only one independent dimensionful quantity, f a , and one independent dimensionless parameter, g s . Therefore, the masses of the particles in the axion sector are all predicted in terms of these two parameters.
We have found that this model features a first-order PQ phase transition, which is very strong. The presence of only few adjustable parameters results in interesting predictions regarding the main quantities associated with the phase transition, T n , β/H I , etc, mainly summarized in the left plot of Fig. 4 and in Fig. 5. We have shown that, like in previous effective axion models with CW PQ symmetry breaking, the phase transition is characterized by a period of strong supercooling, T n T c , when the universe inflated. Thanks to this period the monopole density associated with the breaking SU(2) a → U(1) a is efficiently diluted. Reheating then generically occurs via the unavoidable couplings between the SM and the axion sector due to gluons, which guarantee a rather large reheating temperature. Forḡ s around 1, the model predicts values of T n and β/H I that are within the reach of future GW detectors, such as ET, CE, DECIGO and BBO (see Fig. 5).
The key theoretical tool, which we have used to obtain these results regarding the first-order phase transition, is the calculation of the bounce solutions associated with the tunnelling from the PQ symmetric configuration to the PQ breaking vacuum. Within this formalism we have also checked that the phase transition is mainly due to thermal effects rather than quantum effects: the action of the O(3)-symmetric bounce divided by the temperature is always much smaller than the action of the O(4)-symmetric one.
Finally, the predictivity of the model also interestingly leads to specific features of the GWs spectrum produced by the PQ phase transition. We have compared this theoretical spectrum with the sensitivities of several future detectors such as ET, CE, DECIGO, BBO and advanced LIGO (see Fig. 6), finding conclusively that these experiments will be able to test the fundamental QCD axion model.
We believe that the precision that GW astronomy promises due to the planned worldwide network of GW detectors can make the dream of testing high-scale and fundamental BSM scenarios of UV-completion a reality in near future.