Dynamical Axions and Gravitational Waves

In this paper we explore the possibility of observable gravitational waves as a manifestation of the QCD axion dynamics. In particular, we focus on dynamical axion models which solve the strong CP problem, and include the confinement of a QCD-like gauge group at the TeV scale. We study the resulting chiral symmetry breaking phase transition for models with $N_F=3$ and $N_F=4$ light flavors using the linear sigma model. This model describes the scalar meson spectrum and its interactions, with the diagonal field $\varphi$ as the order parameter. We find that the amplitude of the gravitational wave spectrum depends on the mass of the dynamical axion $\eta'$ via the ratio $m_{\eta'}/m_\varphi$. The resulting spectra may be observed at future mid-range gravitational wave experiments such as AION/MAGIS, DECIGO, and BBO. Moreover, the TeV states can be searched for at colliders and their quantum numbers characterized, providing a unique connection between axion physics, gravitational waves and collider searches.


Introduction
First order phase transitions may result in gravitational wave spectra observable at the next generation of interferometer experiments, and have therefore recently received much attention.
Here we will focus instead on chiral phase transitions, which result from the breaking of a global symmetry after a gauge group confines. The order of such phase transitions has been the topic of decades of scientific inquiry. A famous analytic argument made by Pisarski and Wilczek [11] (PW in the following) is based on a linear sigma model, 1 where Σ is a quark condensate. The PW argument relies on an expansion in = 4 − d, where d is the number of space-time dimensions. To leading order in , it is found that there is no infrared stable fixed point for N F > √ 3. such that a first order phase transitions is expected for a large number of light flavors [11]. 2 Of course, thermal phase transitions should be described by an effective theory in d = 3 dimensions, as time-like fluctuations are cut off at finite temperature. The PW argument can therefore be taken as a guide only. Nevertheless, the result that chiral phase transitions with N F ≥ 3 are first order is commonly accepted, and demonstrated on the lattice for N F = 6 [12].
The gravitational wave phenomenology of first order confining phase transitions has implications for pressing questions beyond the Standard Model (SM). One such open question is whether or not QCD violates CP, known as the strong CP problem. The experimentally small CP violating coupling,θ = θ QCD + arg det M < 10 −10 [13,14] is unexplained in the SM. The most common dynamical explanation for the lack of CP violation in the strong sector is the introduction of a global U (1) Peccei-Quinn (PQ) symmetry that allows the CP violating parameter to be rotated to zero via field redefinitions [15]. Upon spontaneous breaking of the U (1) P Q , the resulting associated pseudo-Goldstone boson that couples anomalously to the field strength of QCD is the axion [16,17].
The strong CP problem can also be solved using massless quarks [18]. In such models, a U (1) chiral rotation of the massless quark fields rotatesθ to zero [19]. This U (1) is also spontaneously broken by the chiral condensation when the quarks' color group confines. The resulting pseudo-Goldstone boson that couples anomalously to the field strength of the confining group is the η , composed of massless quarks. In this class of models, the η plays the role of the composite, or dynamical, axion.
The behavior of the η can be studied in the framework of new exotic confining color groups, which are prevalent in model building to address the strong CP problem [20][21][22][23][24][25][26][27][28][29][30][31]. Heavy or visible axion models predominately utilize exotic color groups to alter the typical (m a , f a ) relationship [22][23][24][25][26][28][29][30][31][32][33]. There are generic features of chiral phase transitions in dynamical axion models that set them apart from general confining hidden sector models. In particular, axion models that utilize an exotic color group typically have matter charged under both QCD and the exotic group. If any of that matter is fermionic, the three colors of QCD guarantee at least an approximate N F = 3 flavor symmetry from the point of view of the exotic group, and this is exactly the minimum flavor symmetry required for a first order phase transition. Also, the pions associated with the broken SU (N F ) symmetry will get quadratically divergent mass terms due to interactions with QCD, and so such models have a generic form for the explicit symmetry breaking from QCD in the linear sigma model. Although one could imagine a number of ways to achieve the dynamical axion solution, in this paper we often use the example discussed in [30] as a benchmark.
Gravitational wave signals from confining phase transitions have recently been studied in dark QCD-like models with N F = 3 [34,35]. In this paper, we will focus on the relation between the dynamical axion and the gravitational wave signature of the chiral phase transition, for N F = 3 and N F = 4. Using the linear sigma model as the low-energy effective theory, we find that the gravitational wave predictions depend sensitively on the mass of the dynamical axion. This is an interesting result, which invites further study of the phenomenology of these models.

The dynamical axion in the linear sigma model
To discuss the origin of the dynamical heavy axion we consider a typical framework, an exotic SU (Ñ ) color sector that confines at a scale Λ Λ QCD , connected to the SM through a quark field ψ charged under QCD and SU (Ñ ), shown in Table 1. This quark field guarantees an approximate N F ≥ 3 flavor symmetry upon SU (Ñ ) confinement as ψ is a triplet of QCD. Table 1: The messenger field between the SM QCD and the SU (Ñ ) confining exotic group.
Here we study the chiral symmetry breaking phase transition associated with the SU (Ñ ) confinement. Upon confinement, the chiral SU (N F ) L × SU (N F ) R symmetry is broken to SU (N F ) V , where N F is the number of flavors in the exotic sector. Below confinement, the dynamics of the hidden sector are modeled by the linear sigma model. The resulting pion fields we discuss are bound states made out exotic quarks.
If ψ is massless, this field is important to the solution of the strong CP problem. Chiral rotations on the ψ field can be used to rotate away the CP violating angle of either QCD or SU (Ñ ), though not both independently. 3 Moreover, a vanishing tree-level mass for ψ ensures that ψ is active when SU (Ñ ) confines, protecting the first order nature of the chiral phase transition. Below confinement, ψ will be hidden in the resulting bound states. One of those bound states, the η , will be the pseudo-Goldstone boson of U (1) A , which is broken explicitly by the instantons of SU (Ñ ). This particle will then couple anomalously to GG, and plays the role of the dynamical axion in the theory below the confinement scale.

Low energy effective theory
At low energies, the dynamics of the quark condensate Σ ij ∼ ψ Rj ψ Li can be described by a linear sigma model. The effective potential for the dynamical field Σ which transforms as a ( ,¯ ) under SU (N F ) L × SU (N F ) R is given by 3 Additional model building is needed to solve the strong CP problem completely in these types of models.
For example, either another massless quark field is introduced or a symmetry relates θQCD andθ.
The chiral condensate Σ ij ∼ f Σ δ ij spontaneously breaks the global chiral symmetry SU (N F ) L × SU (N F ) R → SU (N F ) V . This effect is captured by the linear sigma model when: Σ can be decomposed as: where T a ij are the generators of the SU (N F ) L,R symmetry. The π a are the pseudo-Goldstone bosons associated with the broken combination SU (N F ) A . The ϕ and X a fields are massive bound states associated with the preserved SU (N F ) V × U (1) V symmetry.
In the µ Σ → 0 limit, the chiral symmetry of V (Σ) is enhanced to U (N F ) L ×U (N F ) R , which contains an extra U (1) A restored in (2.1), but spontaneously broken by the chiral condensate. The η in (2.3) is the Goldstone boson associated with the spontaneous U (1) A breaking. The µ Σ → 0 limit does not properly describe properties of a confining gauge group, as we know the axial anomaly explicitly breaks U (1) A to Z N F by quantum effects. These quantum effects originate from instantons of the confining SU (Ñ ). Indeed, it was recognized in [11] that the sum of instanton and anti-instantons generate 2N F -point interactions of the form detΣ + h.c.. Thus the µ Σ -term in the linear sigma model captures the U (1) A -breaking instanton effects in the low energy theory. However, identifying µ Σ with the temperature-dependent instanton density [36], one can draw the conclusion that µ Σ (T ) vanishes for T → ∞, as Debye screening shields electric field fluctuations at high temperature. We will discuss this issue in more detail in section 3.1.2.
Crucially, the η is anomalously coupled to GG, making it a dynamical axion. Above the confinement scale, the strong CP problem is solved by the presence of the massless ψ quarks. The η is a bound state composed of ψ quarks, and therefore its connection to the global U (1) that rotates the CP violating phase is explicit. Note that in the presence of a mass term for the ψ quarks, the U (1) A is classically explicitly broken and no longer meets the criteria for a dynamical axion. 4 The η axion obtains a mass due precisely to the explicit U (1) A -breaking effects from the instantons of the confining group. This explicit breaking should result in an axion potential, which typically takes the form: Using (2.4) we can connect the µ Σ parameter in the linear sigma model to the physical axion mass predicted by the axion's couplings to GG, motivating our choices for linear sigma model parameters in our gravitational wave signal analysis in Section 3.2.
The QCD-colored pions resulting from the spontaneous symmetry breaking should also receive masses since QCD explicitly breaks SU (N F ). The QCD-induced mass of the pions can be included in the linear sigma model by adding the explicit flavor symmetry breaking term This potential term is motivated by the fact that the QCD flavor breaking comes from the The form of Q a L,R will depend on the number of flavors and how QCD interacts with the fields inside Σ. This is worked out explicitly for N F = 3 and N F = 4 below.
Evidence of QCD breaking the flavor symmetry comes in the form of quadratically divergent gluon loops driving the masses of the pions up towards Λ. Their zero-temperature loop-induced masses are: In general, if the dynamical axion solves the QCD strong CP problem, the massless quarks that compose the axion should be coupled to QCD. Typically, the massless messenger quark present at SU (Ñ ) confinement will leave behind some pion states charged under QCD. Thus, pion masses of the form (2.6) are a generic feature of dynamical axion models with hidden chiral phase transitions. Octet pions are present in both the N F = 3 and N F = 4 models, while the triplets are only present in the N F = 4 model, as will be discussed in Section 2.2.2. Colored states can be searched for at colliders via their gluon couplings. Collider searches provide a lower limit on the masses of the octet and triplet pions [38][39][40], Using (2.16) and (2.17), the lower limits on pion masses provide a lower limit for the on the confinement scale: The other relevant bound comes from the lightest dynamical axion. Here the dynamical scale for the axion is the chiral symmetry breaking scale

Meson masses
In this section we compute the spectrum of meson states, required to compute the one-loop thermal corrections to the effective potential discussed in the next section.

N F = 3
The N F = 3 case is representative of models with a massless quark field charged as bifundamental under SU (3) QCD and a confining SU (Ñ ), as shown in Table 1. Models in Ref. [25] and [30] have a confining exotic color (gauge) group with an SU (3) L × SU (3) R (global) chiral symmetry. The SU (3) flavor symmetry is explicitly broken by QCD. The effects of this are captured by adding the ξ-term to the linear sigma model: where T a are the SU (3) generators. We find upon minimizing the potential that The masses of the η and the pion are These fields are pseudo-Goldstone bosons associated with the SU (N F ) A × U (1) A symmetries spontaneously broken by the chiral condensate. The η gets its mass from the anomalous U (1) A breaking. In the linear sigma model, the µ Σ -term breaks U (1) A and so m η should be proportional to µ Σ . The octet pions get their mass from interactions with QCD, and so their mass should be proportional to ξ. Given that we know ξ is generated by a quadratically divergent effect, we can use (2.6) to estimate ξ ∼ 3 α 4π Λ 2 . The masses of the heavy bound states are The heavy states have masses proportional to m 2 Σ . They are held together by the binding energy from the confining gauge group and should have masses near the confinement scale The N F = 4 case we examine is given in Table 2. This is representative of models in which the exotic confining group requires a second independent massless field to independently rotate away the θ-angle of the confining SU (Ñ ) group. Examples of axion models that have this flavor structure near exotic confinement are Ref. [21] and Model I of Ref. [30]. Table 2: The massless quark content that gives an approximate SU (4) L × SU (4) R chiral symmetry broken by the chiral condensate. Interactions with QCD explicitly break the SU (4) flavor symmetry.
QCD explicitly breaks the SU (4) flavor symmetry. The effects of this are captured in the addition of the ξ-term to the linear sigma model: and the fourth row and column filled with zeroes. Σ ij ∼ Ψ Rj Ψ Li , where the flavor multiplet Ψ i = (ψ i , χ) contains three ψ i quarks and one χ. Since only ψ is charged under QCD, the explicit symmetry breaking effects of the ξ-term should only effect the bound states containing ψ quarks. We know that when SU (4) → SU (3) × U (1), fields in the Adjoint representation of SU (4) break as where the 15 is representation of the π a of SU(4), the 1 field is η χ , and the 3,3, 8 fields are the QCD-charged pions. Even when the SU (4) → SU (3) breaking does not happen, this breaking pattern gives us a hint as to which representations of SU (3) are living inside a representation of SU (4). Given (2.6), we expect that π a contains mass eigenstates where π 8 are the color octet pions and π 3 are the color triplet pions. We find upon minimizing the potential that The η and the pion masses are (2.20) The pions π 3 and π 8 have masses proportional to ξ. The mass of the QCD-colored pions are given by the QCD's quadratically divergent contribution to the pion self energy, (2.16) and (2.17). Using this, we can estimate ξ ∼ 3αc 4π Λ 2 . The η χ and η ψ are the singlet pseudoscalars that couple to the confining group's GG. The confinement dynamics that break U (1) A , namely, the instanton effects, are described by the µ Σ -term of the linear sigma model. Since both η χ and η ψ couple identically to GG, these instantons can only give mass to one eigenstate. This explains the spectrum. The η ψ mass is proportional to µ Σ since its mass corresponds to explicit U (1) A breaking. Then the linear sigma model predicts the light eigenstate m η χ = 0. The η χ is not exactly massless, but obtains its mass via mixing with the SM pion below QCD confinement.
The mass of the heavy bound states associated with the unbroken generators, X a and ϕ are: Because the 15 X a R are no longer degenerate, the mass eigenstates have been renamed X 1 , X 3 , X 8 for the particles charged as singlets, triplets, and octets under QCD, respectively. All the heavy states have masses proportional to m 2 Σ . They are held together by the binding energy from the confining gauge group and should have masses near the confinement scale m 2 Σ ∼ Λ 2 . The light state η χ will be present at QCD confinement and will couple anomalously to QCD's field strength tensor through triangle diagram involving ψ quarks. Thus, the η χ is also a dynamical axion, and without any additional mass sources will mix with the SM pions and yield an invisible axion with the typical relationship m a f a ∼ m π f π characteristic of invisible axion models [44][45][46][47][48].
Ref. [30] provides an example model where additional mass contributions to the axion potential raise the η χ mass. This model is thus a heavy axion model with an N F = 4 flavor symmetry at SU (Ñ ) confinement.
The new mass contributions come from small-sized instantons in the UV theory that couple only to χ and not to ψ. In this analysis, the effect is captured by adding a µ SSI -term to the linear sigma model: where picks out only the χ component of Ψ i = (ψ i , χ). This µ SSI -term should account for mass contributions to any state that contains a χ quark, and so will affect both η ψ and η χ , 5 raising the lightest mass eigenstate. Including this new µ SSI -term, we find upon minimizing the potential that The η and the pion masses are Note that now both mass eigenstates contain µ SSI contributions, and so depending on the strength of the SSI instantons interacting with the χ field, the lightest dynamical axion can be made heavy. The masses of the heavy bound states associated with the unbroken generators, X and ϕ are now where the µ SSI provides an additional mass source, though all heavy bound states still contain m 2 Σ ∼ Λ 2 contributions as well.

Gravitational waves from a chiral phase transition
In this section we study the finite-temperature one-loop effective potential generated by the scalar mesons, and explore the parameter space leading to a first-order phase transition. The gravitational spectrum of the phase transition is then derived for two distinct cases, which differ in the number of light fermions: N F = 3 and N F = 4.

One loop effective potential at finite temperature
At one loop, we consider the following daisy-resummed thermal corrections to the potentials (2.9) and (2.13), 6 Here the m i are thermal, field-dependent masses. The sum is over all scalar states (with multiplicity n i ) which couple to ϕ. We consider the zero temperature (Coleman-Weinberg) contribution to the potential to be a redefinition of the parameters in our zero temperature potential (2.1), and do not add it explicitly.

Thermal masses
For N F = 3, the field dependent thermal masses are given by, (as also found in [34]), and for N F = 4, they are given by As a result of the spurion analysis described in the previous section, the masses of the X and the π mesons are no longer degenerate in N F = 4 case, although the thermal contributions to the spectrum remain degenerate.

The µ Σ term during the phase transition
The determinant interaction (µ Σ detΣ + h.c.) in the potential (2.1) is known to be generated by instanton interactions [19]. It is therefore reasonable to assume that the term is proportional to the strength of instanton effects [11]. The temperature scaling of this quantity at large temperatures (and small gauge couplings) is well approximated by the dilute instanton gas approximation (DGA), as large-scale instantons are suppressed [51]. However, as T → T c and α S → 1 (near the confinement scale), large-scale instantons are no longer suppressed. At low temperatures and strong coupling the instantons can therefore not be considered wellseparated, and the DGA is no longer a good approximation. How the instanton density scales during and after the confinement phase transition is currently an open question. A related, well studied quantity is the topological susceptibility, where F (θ, T ) is the θ-dependent free energy and where d(ρ, T ) is the instanton density as a function of instanton size ρ and temperature. The finite temperature behavior of the topological susceptibility χ(T ) is of interest to the lattice community as well as the axion and dark matter communities, because the misalignment production of light fields depends sensitively on its behavior at finite temperature (see e.g. [52,53]). The DGA can be used to compute χ(T ) at high temperatures, which results in χ(T ) ∼ T −8 (we include a brief review of this calculation in appendix A). Lattice simulations agree with this prediction above T 1.5 T c , and indicate a flattening off of the temperature evolution around T ∼ T c [52]. Motivated by these results, we will use dµ Σ /dT = 0 for T N < T c in the following, but alert the reader that this issue should be revisited if further lattice results become available. A strong temperature dependence of the parameter µ Σ would mean a larger nucleation rate than calculated below, and an overall decrease of the gravitational wave amplitude.

Gravitational wave spectra
In the linear sigma model detailed above, the dynamics of the phase transition are captured by the diagonal field ϕ, which has vaccuum expectation value ϕ = 0 at high temperatures, and ϕ = f Σ after the transition. As is well known, the tunneling between the two vacua of this field is described by the scalar bounce ϕ c (r, T ), a spherically symmetric classical solution to the Euclidean equations of motion [54].
We use a combination of a shooting algorithm and a finite difference technique to solve this equation at different temperatures, using the 1-loop thermal potential described in the previous section. From this solution, the thermal parameters of the phase transition are derived. Firstly, the nucleation temperature T N is conventionally defined as the temperature for which a particular volume fraction is in the new phase. We will use, where p(t) is the nucleation probability per unit time per unit volume, t N is the nucleation time, and S E is the Euclidean action evaluated at the bounce solution ϕ c (r, T N ). We have assumed radiation domination to relate t N and T N . Secondly, the nucleation rate is captured by the parameter β (conventionally normalized to the Hubble rate), and can also be related to the bounce action, Then, importantly, the latent heat can be defined by, where the symbol ∆ indicates that the quantity is to be evaluated on either side of the wall (with a relative sign), and where ρ N = π 2 g * T 4 N /30 is the equilibrium energy density at T N , assuming radiation domination.
As was also argued in [35], it is likely that the chiral phase transition does not exhibit runaway behavior: that is, the bubble walls do not keep accelerating until the bubbles collide. The field ϕ couples to several other bosonic degrees of freedom which will result in friction on the bubble wall. For a non-runaway transition, then, the gravitational wave spectrum resulting from colliding acoustic waves in the plasma is expected to dominate. The thermal parameters can then be used to find predictions for this gravitational wave spectrum [1,55], 7 Here z p is a simulation derived factor, which we take to be z p = 5 in the following [55]. For v w → 1, the efficiency parameter is well approximated by [58],

N F = 3
Here we discuss the thermal parameters of the N F = 3 phase transition and the resulting gravitational wave spectrum. We will express our results as a function of the ratio of physical masses m η /m ϕ . For N F = 3, this ratio is given in terms of parameters of the linear sigma model by, . (3.14) Here the inequality is derived from the mass spectrum (2.11); realness of the physical masses implies κ + 3λ > 0 and µ 2 Σ ≥ 0 in the linear sigma model. 8 As we will see below, the upper bound on the ratio m η /m ϕ has a nontrivial implication for the gravitational wave spectrum from this class of chiral phase transitions.
We perform the bounce calculation described in the previous subsection, with the oneloop thermal potential described in section 3.1, for 50 parameter points. We use parameters for which the ratio (3.13) lies in the range m η /m ϕ = [1, We use (m X − m π )/GeV = [2 × 10 2 , 2 × 10 4 ], recognizing that m X > m π for any choice of parameters, cf. (2.11).
We plot our results in the Fig. 1, along with projected constraints from various future gravitational wave experiments. The limits plotted here are for power-law spectra, for which the signals in different frequency domains are correlated [67]. The most sensitive experiment in this observational window is the Big Bang observer (BBO), a proposed fourth-generation space based interferometer experiment, in a hexagram configuration [59][60][61]. The Deci-Hertz Interferometer Gravitational Wave Observatory (DECIGO) is a space-based Fabry-Perot interferometer which was proposed as early as 2001 [62,68]. B-DECIGO [63] is the scaled-down version of DECIGO (the B stands for "basic"). Finally, AION and MAGIS are proposed experiments which apply the recently proposed atomic interferometery technique [64][65][66].
It is clear from the right panel of Fig. 1 that the gravitational wave predictions from this model will not be probed by any currently proposed gravitational wave experiments. This result can be primarily explained by the small latent heat released in this transition, compared to the large radiation energy density at the time of bubble nucleation. "Dark" QCD models, which do not couple to the Standard Model (and as such evade the experimental constraints described in section 2.2) may feature phase transitions at lower scales. Such models may therefore in principle predict larger gravitational wave amplitudes.
Although the signal is not observable, one can still make a few interesting observations. From the left panel of Fig. 1, it is seen that the latent heat parameter α and the nucleation rate It is seen that the latent heat parameter α is typically larger, and the rate parameter β/H is typically smaller for larger ratio m η /m ϕ , associated with a stronger gravitational wave signature (the panels use the same color scaling).
parameter β correlate positively and negatively with the ratio of masses m η /m ϕ respectively. These correlations are also true for the ratio f Σ /T c (the ratio of the value of the VEV to the temperature at which both phases are degenerate -sometimes referred to as the strength of the transition). Since the mass of the dynamical axion m η is determined in part by instanton effects as described in the previous subsection, this correlation motivates further study of this parameter at low temperatures.

N F = 4
For N F = 4, we may again express the ratio of physical (zero temperature) masses in terms of parameters of the linear sigma model, It is seen that in this case, the ratio is not bounded from above. For N F = 4, we will study a larger range, m η /m ϕ = [1,10], noting that a large mass ratio is a natural expectation of models such as [30]. We repeat the calculation described in the previous subsection, for 50 parameter points as before. We choose parameters such that the meson masses are in agreement with experimental constraints, m η /GeV = [10 3 , 2 × 10 -13 Figure 2: Thermal parameters (left) and gravitational wave signatures (right) of a chiral phase transition with N F = 4 and g * = 120, for the benchmarks described in the text. As for N F = 3, it is seen that the latent heat parameter α is typically larger, and the rate parameter β/H is typically smaller for larger ratio m η /m ϕ , associated with a stronger gravitational wave signature (the left and right panels use the same color scaling).
We plot the results in Fig. 2. In contrast to the N F = 3 case, we note that the predicted gravitational wave spectrum of the N F = 4 case may be observed at future interferometer experiments, in line with the naive expectation from the behavior observed in the previous subsection, and the larger ratio m η /m ϕ = [1, 10] 9 . In particular, for ratios m η /m ϕ 7, the gravitational wave signal may be detected by atom-interferometers AION and MAGIS. From the left plane of Fig. 2, we note that the latent heat released by transitions with larger ratio m η /m ϕ is larger, while the predicted nucleation rate is smaller.
The result is interesting, in particular in light of the consideration in subsection 3.1.2. The realization of the large mass ratio relies crucially on the value of the parameter µ Σ at the nucleation temperature, and thus on the explicit breaking of the U (1) A symmetry. As such, further study of the finite temperature behavior of the instanton density are well-motivated.

Discussion
This paper has discussed the gravitational wave signatures of models of dynamical axions and confinement at the TeV scale. The order of the confinement phase transition relies on the number of light fermions at the confinement scale. An analytic argument based on an 9 This m η mϕ upper value estimate comes from the fact that (assuming large separations of scales µSSI and Λ) the η gets its mass from confinement and therefore m 2 η ∼ Λ 4 f 2 , where f is the chiral symmetry breaking scale and is estimated to be Λ ≤ 4πf . Then m 2 η ≤ (4π) 2 Λ 2 and since mϕ ∼ Λ, then m η mϕ ≤ 10. expansion in = 4 − d [11] implies that phase transitions with N F ≥ 3 are first order, and therefore feature a gravitational wave spectrum.
Using the linear sigma model, we studied the cases N F = 3 and N F = 4. The gravitational wave predictions of these models are plotted in Figs. 1 and 2 respectively. We note that the predictions of the first model evade experimental detection at the presently proposed gravitational wave observatories. For N F = 4, however, the signals may be observable at the Big Bang Observer (BBO) [59][60][61], (B-)DECIGO [62,63], and the MAGIS-AION atom interferometers [64][65][66].
An interesting result is that the amplitude of the gravitational wave spectrum depends on the ratio of the mass of the dynamical axion m η , to the mass of the order parameter of the phase transition, the scalar ϕ. In particular, the latent heat released in the transition (conventionally captured in the parameter α), and the nucleation rate (β/H) correlate with m η /m ϕ positively and negatively respectively. In some models, a large ratio may be a natural prediction [30].
The importance of the ratio m η /m ϕ suggests interesting questions for future research, in particular about the origin of the µ Σ parameter in the linear sigma model. This term constitutes an explicit breaking of the global U (1) A symmetry. The determinental interaction originates from instanton dynamics, which is known to have a strong temperature dependence at large temperature. The behavior of this parameter at the confinement scale should be investigated further to allow for more detailed studies of the phase transition.
The analysis in this paper leaves open the question of phase transitions in models with N F > 4, such as the recently proposed high-scale color confinement model [69]. The determinant operator is irrelevant for N F > 4, and the strength of the phase transition in the linear sigma model may instead correlate with a further explicit breaking, or the generation of the η mass in such scenarios. The fact that the mass dimension of µ Σ depends on N F , while the gravitational wave detection prospects also depend strongly on µ Σ means studying chiral phase transitions for different values N F may be interesting. At N F 3, a conformal window is known to exist, though its exact location is the topic of ongoing research [70][71][72].
Our study of gravitational wave signatures shows a sizable hierarchy between the scalar ϕ and pseudo-scalar η is preferred for detection. At the same time, these are TeV-range states which can be searched for at colliders. We used current limits on the colored states to guide our parameter search, but these colored states would also induce loop-level contributions to couplings of the ϕ and η to gluons and hence could be looked for at the LHC via dijet signatures [73,74]. One can envision a future situation where a signature in gravitational waves is found, and that would guide searches for two correlated states in dijets at colliders. Further angular analysis of the dijet final state could also allow us to determine the CP properties of these states, and robustly support the origin of the gravitational wave signature as manifestation of a dynamical axion explanation of the QCD CP problem. This connection between gravitational wave signatures and collider searches is an area that we plan to develop further.
Note added. While this paper was in its final stages, [35] appeared on the arXiv. Al-though the focus of the current work is different, the analysis overlaps partially with the analysis in [35]. To this extent, the works are consistent, although the analysis here applies results from more recent gravitational wave lattice studies, resulting in smaller amplitudes. using dimensional analysis and Σ ∼qq/f 2 .
The topological susceptibility and (A.4) in the DGA approximation are plotted in Fig.  3. It is seen that for high temperatures, the slope approximates T −8 as is known. However, this approximation is no longer trusted in the regime in which the PT takes place. Lattice studies in the regime T < T c are inconclusive, though it is generally expected that at very low temperatures χ(T ) ∼ T −n where 0 n 8.