Peccei-Quinn Phase Transition at LIGO

The LIGO observatories can potentially detect stochastic gravitational waves arising from phase transitions which happened in the early universe at temperatures around $T\sim 10^{8}$ GeV. This provides an extraordinary opportunity for discovering the phase transition associated with the breaking of the Peccei-Quinn symmetry, required in QCD axion models. Here we consider the simplest Peccei-Quinn models and study under which conditions a strong first-order phase transition can occur, analyzing its associated gravitational wave signal. To be detectable at LIGO, we show that some supercooling is needed, which can arise either in Coleman-Weinberg-type symmetry breaking or in strongly-coupled models. We also investigate phase transitions that interestingly proceed by first breaking the electroweak symmetry at large scales before tunneling to the Peccei-Quinn breaking vacuum. In this case, the associated gravitational wave signal is more likely to be probed at the proposed Einstein Telescope.


Introduction
The recent detection of gravitational waves (GWs) by LIGO [1] represents the beginning of a new era in the exploration of the universe. In a few years LIGO-VIRGO has compiled a sizable catalogue of detected binary merger events [2], and the prospects to further increase the sensitivity and even to build more observatories look promising.
In the zoo of candidates for GW signals there is one that stands out from the point of view of high-energy physics: the stochastic GW backgrounds originating from cosmological first-order phase transitions in the early universe. First-order phase transitions develop by the formation of bubbles that expand, collide and percolate. The bubble wall collisions are violent events that occur everywhere in space at a given cosmological time, leading to sizable stochastic signals that remain as a relic cosmological background analogous to the cosmic microwave background, but in GWs. Since GWs are a form of radiation, after their production the fraction of the energy density that they carry keeps constant in the radiation dominated epoch, thereby giving a relic background that can be detected now, no matter how early they were produced and how high the temperature of the universe was. The temperature of the phase transition is encoded directly into the power spectrum of the signal, mainly in the peak frequency that scales as f peak ∝ T . A first-order phase transition at T ∼ TeV peaks in the frequency sensitivity band of LISA, while GW observatories with higher frequency sensitivity bands can probe even higher energies [3].
The main motivation for this work is that the LIGO frequency band corresponds to first-order phase transitions which could have happened when the early universe was at a temperature around 10 7 − 10 8 GeV. This roughly coincides with the lowest possible energy scale where the Peccei-Quinn (PQ) symmetry U (1) P Q had to be broken in QCD axion models which solve the strong CP problem of the SM [4,5]. In other words, the axion solution to the strong CP problem predicts a phase transition that can occur around this scale. Then, LIGO-VIRGO has the chance to discover this PQ phase transition if it was of first-order and "strong enough".
The purpose of this work, then, is to browse through the simplest incarnations of the PQ mechanism and see in which cases a detectable first-order phase transition is obtained. We will focus on the minimal KSVZ [6,7] and DFSZ [8,9] models, as well as supersymmetric and strongly coupled versions of them. As we will see, the most important requirement is that the models manage to give a strong enough (and long enough) transition. We will show that this is the case for certain regions of the parameter space of the DFSZ model, and is more favorable, when the PQ breaking is driven by the Coleman-Weinberg mechanism [10]. Also, we will show that strongly-coupled models of PQ breaking lead to long periods of supercooling which end with strong GW signals detectable at LIGO. Furthermore, we investigate the occurrence of a two-step PQ phase transition in DFSZ constructions, with an intermediate second-order electroweak phase transition at very high scales, before ending in the PQ broken minimum with a first-order phase transition. Crucially, we show that it is possible to obtain a significant amount of "cooling" in these cases, albeit much milder than in the aforementioned supercooled scenarios.
We must remark that astrophysical bounds on the PQ scale F a require F a 10 8 GeV, that is slightly above the scales at which LIGO is most sensitive. Nevertheless, as we will see, the temperature of the phase transition can be actually slightly smaller than F a (by up to a factor ∼ 10), which in the end allows LIGO to probe these scenarios. The capability of LIGO to probe PQ phase transitions has also been pointed out and partially discussed in [11,12].
We will also include in our analysis the projected sensitivity for the Einstein Telescope (ET) [13]. The enhanced sensitivity with respect to LIGO offers the opportunity to probe a much larger area of the parameter space. Therefore ET holds a great promise to probe axion physics.
The article is organized as follows. In Section 2 we show the sensitivity of Advanced LIGO and the proposed ET to the parameters α, β/H * and T * of the phase transition. In Section 3 we present the simplest models of PQ breaking, the KSVZ and DFSZ models, analyze their type of phase transitions, and study their GW signals. Section 4 is for conclusions.
2 Sensitivity of Advanced LIGO and ET to first-order phase transitions In this section we show that if the spontaneous breaking of the PQ symmetry occurred via a first-order cosmological phase transition, then this would have left a stochastic GW signal potentially detectable by LIGO as well as future GW observatories. Indeed, in this case the transition proceeds by bubble nucleation and the collisions between the bubbles as well as the motion of the thermal plasma which surrounds them are sufficiently violent events to generate significant GWs. Let us start by reviewing some basic notions that characterize first-order phase transitions and how they can source GW backgrounds. First-order phase transitions occur when there are at least two minima in the scalar potential (which generically depends on the temperature T ) and the universe, initially trapped in the minimum with higher energy at high T , transits to the minimum with lower energy either by thermal fluctuations or quantum tunneling. In both cases this proceeds at a certain 'nucleation' temperature T = T n through the formation of bubbles of a critical radius R which then expand and percolate. The rate at which bubbles are produced per unit volume is given by Γ = A e −S B where S B is the action of the critical bubble, or bounce, and A is a prefactor that is usually of order 1/R 4 . In order for the phase transition to be completed in an expanding universe, we must have Γ H 4 where H is the Hubble rate. The nucleation temperature T n is therefore determined by Γ ∼ H 4 which leads to where we have taken the approximation A ∼ T 4 . The calculation of S B depends on the details of the potential and has to be performed case by case. The parameters which characterize the first-order phase transition, and which are relevant for the GW signal, are the following: 1. The temperature T * at the time when the phase transition completes. It can be estimated from energy conservation by equating the latent heat ∆V (the difference of the potential between the false and true vacuum) plus the energy density in the thermal bath at the nucleation temperature to the energy density of a thermalized plasma, ρ γ (T * ) = ρ γ (T n ) + ∆V with ρ γ (T ) = π 2 g * /30 T 4 . Assuming that the number of relativistic degrees of freedom, g * , does not change much between T * and T n , one gets 2. The strength of the first-order phase transition α, characterized by the energy density going into the bubbles over the thermal energy density of the surrounding plasma: 3. The inverse of the duration of the phase transition β = [(dΓ/dt)/Γ] Tn [14], which can be approximately determined as where we have assumed fast reheating so that H * ≡ H(T * ) H(T n ), and the −4 arises from A ∝ T 4 .
4. The bubble wall velocity v w , which is determined by the interaction of the bubble walls with the surrounding plasma. The latter exerts a friction force on the propagation of the walls. In very strong phase transitions (α 1) one expects that the pressure difference across the bubble walls dominates over the friction of the plasma and bubbles run away, thus v w → 1, except in certain cases [15,16]. In weaker phase transitions (α 1), we will take the estimate that v w is expected to be close to the speed of sound in the plasma [17].
The collisions of bubbles during the phase transition can source GWs of a sizable amplitude. Production of GWs in a first-order phase transition has been much discussed previously -see e.g. [14,18,19] for recent reviews. The generated GW signal represents a stochastic background and as such it is best characterized by its power spectrum. It is customary to express it in terms of the fraction of the present energy density in GWs per unit decade in frequency, This signal can be separated into three distinct contributions, arising from the collision of the scalar wall profiles, the sound waves in the plasma and from turbulence, respectively. The shape and size of each contribution can be estimated separately as reviewed in [14,18,19]. In all cases the power spectrum has a maximum at a characteristic frequency basically determined by the inverse duration β, and deviates from the maximum by two different power laws. In this work we will simply assume the following expressions for the GW spectra as functions of the parameters of the phase transition, quoted in [14,19]: • From bubble wall collisions, with h the dimensionless Hubble parameter, κ φ an efficiency parameter which can suppress the contribution from bubble collisions when the effects of the thermal plasma cannot be neglected, and the peak frequency today given by GeV with α ≈ 3.5, and F a = 10 9 GeV with α ∼ 10 6 . We have set v w = 1 and show the signals which arise from only bubble wall collisions with κ φ = 1 (blue lines) and from only sound waves in the plasma with κ sw = 1 (red lines).
• From sound waves in the plasma, with the peak frequency today given by T * 10 7 GeV g * (T * ) 100 1 6 . (2.10) The efficiency parameter κ sw ≤ 1 quantifies the fraction of the latent heat which goes into bulk motion. Here we shall assume the expression obtained e.g. in [20,21] (see [14] for a recent discussion), which holds for v w ∼ 1, • The contribution from turbulence Ω t is suppressed (while also being more uncertain), and we will set it to zero for our estimates.
A convenient way to know whether a signal is detectable by a given GW observatory is to compare the power spectrum to the so-called power-law integrated curves [22], which express the sensitivity as the minimal Ω needed for detection as a function of f . In this work we will be interested in the frequency range which can be probed by ground-based interferometers. We show in Fig. 1 the sensitivities of the current Advanced LIGO (with O1 and O2 data [23]), as well as the projected sensitivities of the design Advanced LIGO and ET [13]. For illustration, we also include in the figure some representative power spectra which arise in the PQ model discussed in Section 3.2.1. We also include the indirect limits resulting from CMB data [24]. The CMB bound is on the integral df Ω GW /f , so how this translates to a bound on the spectral density depends on the shape of the assumed spectrum. This is why we show this bound as a thick line in Fig. 1. One can then easily obtain which part of the parameter space (α, β, T * , v w ) corresponds to detectable signals by simply checking whether the power spectrum overlaps with the instrument sensitivities. In Fig. 2, we show this in the T * − β/H * plane for the two representative values α = 0.1 and α = 3 (the detectable region is to the left of the lines). Clearly, for strong first-order phase transitions with α 3, LIGO at design sensitivity can detect signals that fall into the relevant range for the PQ models, T * ∼ 10 7 − 10 8 GeV, and it can reach values of β/H * as large as 10 2 − 10 3 . Interestingly, even the current O1 and O2 runs of LIGO are capable of ruling out phase transitions with β/H * 10. As can be seen in Fig. 2, the improvement on these figures by ET would be rather impressive.
On the other hand, for small α the possibility to detect a first-order phase transition almost completely fades away at LIGO. This is illustrated in Fig. 3 where we show the LIGO sensitivity in the T * − α plane. By taking reasonable values of β/H * O(10), one clearly sees that in order to possibly detect a signal at LIGO the transition needs to be strong, that is, with α 1. The situation could be slightly improved with ET which could reach down to α ∼ 0.1.
One must be aware that the collection of unresolvable black hole and binary neutron star mergers creates an additional stochastic GW background [23], the so-called 'popcorn' background. Given the event rates of these mergers, the magnitude of the popcorn in the LIGO frequency band is around h 2 Ω ∼ 10 −9 , which enters in the detectability range for ET and marginally so for LIGO at design sensitivity. This signal represents a 'foreground' for α T * /GeV Figure 3: Sensitivity lines for current LIGO (solid) and LIGO at design sensitivity (dashed) in the T * − α plane assuming fixed β/H * equal to 10 and 100 as indicated. Blue or red colors refer respectively to whether the signal is assumed of bubble wall collisions type, (2.7) with κ φ = 1 and v w = 1, or of sound waves type, (2.9) with κ sw set by (2.11) and v w = 1/ √ 3. The points to the right of the curves represent detectable signals.
the cosmic GW backgrounds, and it should be subtracted away in order to be able to detect a possible background from cosmological phase transitions. This seems in principle feasible since the power spectra from popcorn and phase transitions differ significantly [23].
It is interesting to note that that PQ models predict actually two more stochastic sources of GWs in addition to the possible one from the PQ phase transition. Indeed, since the PQ symmetry is a global U (1) symmetry, it is granted that global cosmic strings will form at the symmetry breaking scale F a . Cosmic string networks radiate GWs, but this is negligible for global strings. Also, at temperatures of order GeV, QCD effects further break U (1) P Q and lead to domain walls, attached to the global strings. The string-wall network then disappears around the QCD scale via rather violent processes where large topological defects collapse and collide. This string-wall network anihilation is similar to a cosmological phase transition and it may give a larger signal. For QCD axion models the peak frequency of this signal must be around f ∼ 10 −10 − 10 −7 Hz, which is in the sensitivity range of Pulsar Timing Array observatories. Unfortunately, the recent numerical simulations of these networks [25] suggest that the spectrum of this signal falls a bit short to be detectable.

Peccei-Quinn Phase Transition and its GW signal
Having seen the current and future reach of GW interferometers, we now move to the particle physics motivation of this work: the QCD axion solution to the strong CP problem. We start by providing a lightning description of axion physics to set notations, then we investigate the occurrence of a first-order phase transition in the simplest PQ constructions.
Axion models are characterized by having a global U (1) P Q symmetry with a U (1) P Q − SU (3) c − SU (3) c anomaly. The U (1) P Q is assumed to be spontaneously broken by the vacuum expectation value (VEV) of a scalar Φ at some scale F a . The axion a(x) then is the Nambu-Goldstone boson that arises from this breaking, Φ = e ia(x)/Fa F a / √ 2 + · · · . Due to the U (1) P Q − SU (3) c − SU (3) c anomaly, the axion couples to gluons as which leads to a potential for the axion through QCD instantons. This gives a = 0, solving the strong CP problem, and an axion mass We can categorize PQ models into two different types, depending on the origin of (3.1). Those referred to as KSVZ models [6,7] contain extra quarks which are responsible for the anomaly and which generate the term (3.1). On the other hand, those referred to as DFSZ models [8,9] contain extra scalars which, after being integrated out, generate the coupling where q refers to SM quarks. By a chiral rotation of q, the axion can be moved from (3.3) to (3.1). Below we discuss the minimal versions of these types of models, their phase transitions and potential GW signals.

KSVZ axion models
The minimal model of this type consists of a scalar Φ and an extra quark Q L , Q R with PQ charges q Φ , 0 and −1, respectively. The interactions, dictated by the PQ symmetry, are given by where we have set q Φ = 1/n. Unfortunately, in this model in which Φ only interacts with itself and an extra fermion, the phase transition is second order, and no significant GWs are expected to be produced from the phase transition. We could couple Φ to the SM Higgs, e.g., |H| 2 (κ|Φ| 2 − µ 2 ). However, in order to achieve a viable electroweak (EW) symmetry breaking, we need to tune κ Φ 2 ≈ µ 2 . This constraint has not allowed us to find a region of the parameter space where the PQ phase transition is strongly first-order (see also [12]).
We will see later that supersymmetric versions of the KSVZ model can however have a strong first-order phase transition.

DFSZ axion models
This type of models instead consist of the PQ scalar Φ and one extra scalar SU (2) L doublet beyond the one in the SM. We denote the two doublets as H 1 and H 2 . Their hypercharges are Y = 1 and Y = −1 and we choose their PQ charges as 0 and −1, respectively, while the PQ charge of Φ is q Φ . The model should also contain at least one SM quark charged under PQ. A minimal option is that only u R is charged under PQ, with PQ charge 1. The interactions are then fixed by the U (1) P Q symmetry to be for the quarks in the first family, while the rest of the SM fermions couple only to H 1 . The scalar potential is given by , and all couplings are real (κ 3 can be made real by a field redefinition). We will for definiteness fix all couplings to be positive in this section. For the real parts of the U (1) EM -neutral components, where λ 12 = λ 3 + λ 4 . The mass matrix of h 1,2 at the PQ-breaking minimum φ = f is given by In order to obtain the observed electroweak scale, the determinant of the mass matrix has to be tuned such that This is the hierarchy problem which we do not address here but which will be considered below. The SM Higgs is given by the linear combination H = cos θ H 1 +sin θH 2 (H 2 = iσ 2 H * 2 ) which diagonalizes M 2 H and whose mass squared is of order m 2 W . Notice that the mixing angle θ enters into the expressions for the SM fermion masses: m d = y d cos θ v/ √ 2 and m u = y u sin θ v/ √ 2. By integrating out the heavy Higgs doublet, one gets the coupling (3.3). The original DFSZ proposal [8,9] has n = 2 (q Φ = 1/2) and all three SM up-type quarks are charged under PQ. This choice leads to a cosmological problem [26] after the QCD phase transition, since the domain wall number parameter N DW is larger than one (in particular N DW = 6 in the original DFSZ proposal). This can be evaded by the introduction of a further small source of explicit breaking of the PQ symmetry [27]. Here instead we make a different choice and focus on n = 1 (q Φ = 1). In this case, if only the first-family u R is charged under the PQ symmetry, we have N DW = 1 and we avoid the domain wall problem. Other choices for the PQ charges and for n will, however, not substantially change our results on phase transitions in these models.
From now on, since m W f , we drop the EW scale in our computations. Thus, the tuning (3.9) reduces to In our study of the DFSZ model, we use (3.10) to fix the parameter κ 3 . The potential (3.7) is then characterized by nine parameters: the scale f , the mass parameters µ 2 1 , µ 2 2 , the selfcouplings λ 1 , λ 2 and λ φ , the quartic couplings κ 1 , κ 2 , λ 12 . Furthermore, the potential (3.7) is a function of the three scalar fields h 1 , h 2 and φ. Nevertheless, we will focus on cases where h 2 either vanishes or can be assumed to quickly track its minimum during the phase transition. We will therefore not study its dynamics during the phase transition and only consider its loop effects on the potential for h 1 and φ.
It is thus only in the two-dimensional field space of h 1 and φ that we will look for a first-order phase transition. In this field space, the potential with signs as chosen in (3.7) can have two minima away from the origin O, which we denote with A and B, located along the φ and h 1 direction respectively (see Fig. 4): Our universe will correspond to the PQ-broken minimum A. Therefore, in order to avoid any danger of having an energetically more favorable vacuum at B, we require V (A) < V (B). This implies the following lower bound on λ φ : This lower bound is only valid at tree level and can be modified by loop corrections. The point B can either be a local minimum or a saddle point of the potential. It is a local minimum (the mass of φ is positive at B) if the following upper bound on λ φ is satisfied: We then find two possibilities for a strong first-order phase transition in the DFSZ construction (shown in Fig. 4): The barrier can be either induced by thermal corrections, mainly thanks to the cubic term T φ 3 , or by one-loop corrections of Coleman-Weinberg type (see (A.1)). The latter is more promising for a strong first-order phase transition, but it requires the mass parameters to be very small compared to F a .
II. O → B → A, first along the h 1 direction and later along some φ − h 1 trajectory. If (3.13) is satisfied, a tree-level zero-temperature barrier separates the minima A and B which can lead to a first-order phase transition in the second step. In this case the universe goes through an intermediate phase with a large EW symmetry breaking scale.
We explore the two possibilities above in the following subsections Sec. 3.2.1 and Sec. 3.2.2, respectively.

Thermal and Coleman-Weinberg driven first-order phase transition
Let us consider the phase transition in the direction of φ (trajectory I in Fig. 4). To ensure that h 1 stays zero during the phase transition, we roughly need T n µ 1 for signs as chosen in (3.7) (otherwise one first rolls/tunnels towards the h 1 -direction, leading to a trajectory like II in Fig. 4). This limits the smallest T n that is achievable. The smaller T n , however, the stronger is the GW signal as we will see below. Another option is to flip the signs of both κ 1 and µ 2 1 in (3.7). One can show that if µ 2 1 is chosen sufficiently large, a tachyonic direction in h 1 and h 2 only develops for φ very close to its minimum. Both fields can therefore be consistently set to zero and their dynamics ignored during the phase transition. 2 We will further assume that all couplings to h 1 are sufficiently small and it is sufficiently heavy that we can also ignore its loop-corrections to φ and h 2 . We are then left with φ and h 2 , where the latter affects the dynamics of the phase transition only via loop corrections. The resulting potential for φ at loop-level and for finite temperatures is discussed in Appendix A.
A first-order phase transition can occur due to a thermal barrier generated by the cubic term ∼ T φ 3 , mostly arising from loops of h 2 . Nevertheless, when the daisy masses are included (see Appendix A), this cubic term is diminished and the barrier is usually small (see e.g. [28] and [29] for a recent discussion). The resulting values of α are then small and those of β/H * large, leading only to a weak GW signal. This can be seen for example for the point marked by p 1 in Fig. 5, calculated with κ 2 = 2, λ φ ∼ 10 −2 and λ φ f ∼ 10 6 GeV. As can be seen in the plot, the phase transition for this case has β/H * ∼ 100, while α ∼ 0.2. Note that the barrier in this case already has a contribution from the Coleman-Weinberg corrections which we discuss below. A purely thermal barrier would have even larger β/H * and smaller α.
A second more promising possibility for a strong first-order phase transition arises in the limit in which the mass parameters are small, µ 2 2 , λ φ f 2 f 2 . In this case, the T = 0 potential for φ becomes almost scale invariant and can be written as (3.14) Due to one-loop corrections, λ φ (φ) depends logarithmically on φ and the potential is thus of Coleman-Weinberg type [10] (when the logs are large, this potential must be RG-improved).
If λ φ (φ) is negative for small φ and turns positive for large φ, a minimum develops close to where the coupling crosses zero. More precisely, the minimum is determined by where β λ φ = dλ φ /d ln φ. Notice that now F a ≡ φ = f . Considering only the couplings λ φ and κ 2 , we have From (3.15), we can fix one parameter, say λ φ , and therefore we are left with only one free coupling, κ 2 . Using (3.15), we obtain at the minimum The phase transition of Coleman-Weinberg models with a potential given by (3.14) was first studied in [30]. Let us sketch here how this proceeds. When non-zero temperature effects are included, the potential at small φ is always dominated by thermal corrections which lead to where D φ is given in (A.7). Therefore at any non-vanishing temperature, the curvature of the potential is always positive near φ = 0 and this point is a (local) minimum. In fact, at very high temperatures, the thermal corrections are so large that the minimum (3.15) at φ = F a is lifted, and the point φ = 0 is the only minimum of the potential. This implies that at a certain temperature T c , the two minima are degenerate, and it becomes favorable to tunnel from φ = 0 to φ = F a . Notice that the barrier is generated thanks to λ φ being negative for φ ≤ φ .
As was discussed in [30], O(3)-symmetric bubbles dominate tunneling in this case and in the limit of small temperatures T their action is well approximated by From this, we see that S B can slowly evolve from large values to small values, since −λ φ (T ) grows as T decreases. This can eventually allow the criterion in (2.1) to be satisfied and thus the phase transition to happen at some temperature T n F a . While trapped in the false vacuum, the universe inflates with H 2 = ∆V /(3M 2 P ) and supercools. We can calculate T n using (2.1) where now (3.20) From (3.16), we see that the smaller κ 2 , the slower does −λ φ (T ) grow with decreasing T and therefore the more supercooling we have. Notice that there is a lower bound for T n , since S n also decreases with T n and at some point becomes too small and S B can never reach its value. Due to this (long) period of inflation, where the temperature drops exponentially, the thermal plasma is diluted, and we have α 1. Furthermore, from (2.4) and (3.19), we obtain We can now see under which conditions a slow transition can be achieved. In principle, since β λ φ in (3.21) is one-loop suppressed, one would expect that β/H * ∼ 1 can be easily achieved. However, also −λ φ (T n ) is one-loop suppressed near the minimum as follows from (3.15). In order to make it larger than that, one needs T n F a . To be more explicit, let us consider the one-loop coupling λ φ (φ) ∼ −β λ φ ln F a /φ. We then roughly obtain from (3.21) which reaches values of order one at T n F a . Having α 1 and the possibility of β = O(1), this scenario then can lead to a maximal signal in GWs, which we expect to be mainly sourced by the collision of runaway bubble walls themselves since supercooling exponentially dilutes the thermal plasma around them. From (2.2) together with (3.17), we can relate T * to F a : This predicts T * to be slightly below F a , making LIGO and ET (see Fig. 2) quite suited to test the interesting region F a ∼ 10 8 − 10 10 GeV. We have calculated the properties of the phase transition by numerically solving the bounce equation for the potential (3.7) plus its thermal and Coleman-Weinberg corrections as discussed in Appendix A. We have performed this calculation for both O(3)-and O(4)symmetric bubbles and have confirmed that the former indeed dominate. The resulting values of T * and β/H * for F a = 10 8 , 10 9 and 10 10 GeV are shown in Fig. 5. We have chosen the two relevant mass parameters, µ 2 and λ φ f (where λ φ is the tree-level coupling), equal for concreteness and hierarchically smaller than F a in order to be in the Coleman-Weinberg regime. Furthermore, we have fixed λ φ as discussed above and scanned over different values of κ 2 ≤ 2. This gives rise to the solid lines in Fig. 5 which for each F a from top to bottom correspond to µ 2 = λ φ f = (10 −2 ,10 −3 ,10 −4 ) F a , respectively. By decreasing κ 2 , one moves along these lines towards smaller T * (as is expected from (3.23)). As follows from the discussion above, as long as T n µ 2 , λ φ f , we have that T n and β/H * decrease if one lowers κ 2 . This regime corresponds to the parts of the lines in Fig. 5 with positive slope. Eventually, however, one reaches T n ∼ µ 2 , λ φ f . Since we have chosen the mass of φ to be tachyonic (cf. (3.7)), this mass compensates the thermal barrier (cf. (3.18)) at lower temperatures and the phase transition thus always happens at T n ∼ λ φ f if one lowers κ 2 further. Since this removal of the barrier happens rapidly at around T n ∼ λ φ f , β/H * then begins to grow again for decreasing κ 2 . This regime corresponds to the parts of the lines in Fig. 5 with negative slope. We thus find that for every given hierarchy between F a and µ 2 , λ φ f , there is a minimal β/H * that can be reached. Furthermore, the dash-dotted lines in Fig. 5 show results of an analytical approximation following (3.19), (3.21) and (3.23) for the case µ 2 = λ φ f = 0. As expected, this case allows to reach much lower values of β/H * . Note also that the solid lines only delimit points with α ≥ 3, while some representative points with α < 3 are shown in red. The values of α always increase on the parts of the lines with positive slope, while they eventually decrease again on the parts with negative slope. The restriction to α ≥ 3 was made since the amplitudes of the GWs becomes independent of this parameter in the limit of large α (see (2.7) and (2.9)). The current and expected reaches of the GW observatories are then shown in Fig. 5 for α = 3 (and setting v w = 1 and κ φ = 1 or κ sw = 1). Since the amplitudes of the GWs increase by about 40% when going from α = 3 to very large α, the true reaches in the very supercooled regime are slightly higher than what is shown in Fig. 5. We see from Fig. 5 that part of the parameter space could be already detected at LIGO, while other parts will have to wait for ET.

Cooled two-step phase transition
We now focus on the case of a two-step first-order PQ phase transition, along the trajectory II in Fig. 4. Let us first understand under what circumstances the two-step phase transition can occur and be strong enough to source a detectable GW signal. From Fig. 3, it is clear that LIGO can probe only transitions with α > 1. In the case of a standard two-step transition, where the minimum B develops at a temperature T h 1 which is higher than the temperature T φ at which the PQ minimum appears, these values of α are difficult to obtain. Indeed in this situation the universe cannot cool much if the transition is to be completed, since the barrier between the two minima is already present at tree level. This is in contrast with the previously discussed Coleman-Weinberg driven scenario.
However, in the DFSZ scenario a new possibility arises: namely, that T φ > T h 1 , but that below T φ the universe is stuck for a while at the origin, due to a loop-induced barrier which opposes rolling/tunneling along the φ direction. In this case, a two-step transition can occur, as below T h 1 the universe tracks the local minimum in the h 1 direction (second order/crossover phase transition). If T h 1 is sufficiently small, large values of α are obtained whenever the transition can complete. For this reason, here we focus on this cosmological history.
We already know of one way to realize this: that is, to make use of the Coleman-Weinberg induced barrier in the φ direction. Alternatively, a barrier induced by φ 3 T terms arising from thermal loops may also suppress tunneling, although it requires large values of κ 2 . In both cases, the crucial ingredient which is peculiar to the DFSZ scenario is the presence of extra bosonic fields coupled to φ, beyond the content of the doublet H 1 . For concreteness, here we focus on the case in which tunneling along φ is suppressed because of the barrier induced by Coleman-Weinberg corrections due to h 2 loops. We then discuss the values of λ φ and κ 1 which allow for this scenario to occur, while we keep the rest of the parameters fixed as follows.
f . Also we take µ 2 = 0.1f and κ 2 ∼ 1. Furthermore, λ 1 is related to the SM Higgs quartic coupling, 3 which at the energies we consider is of order 0.01. For this reason we take λ 1 0.01.
A local minimum in the h 1 direction occurs if the upper bound (3.13) on λ φ is respected. For λ φ 10 −3 , this is easily satisfied and the potential in the φ direction is dominated by Coleman-Weinberg corrections due to h 2 . This also ensures that the tree-level lower bound on λ φ is relaxed, as the minimum A is always the global minimum of the potential. Interestingly, completion of the transition from B to A is facilitated in this case, since the minima are always significantly non-degenerate.
We then proceed to a numerical investigation of the parameter space for this type of two-step transition. As mentioned above, even though the potential is a function of three fields, we can focus on the dynamics of φ and h 1 only. The rest of the fields of the DFSZ model will only affect the potential of φ and h 1 at the loop level. These are all components of the doublets H 2 and H 1 , the imaginary part of Φ, the EW gauge bosons and the top quark. We fix µ 1 = 0.09f, µ 2 = 0.1f, λ 1 = 0.05, λ 2 = 0.01, λ 12 = 10 −3 and the gauge couplings as well as the top Yukawa coupling to 0.6, as appropriate for f ∼ 10 8 − 10 10 GeV. Finally, in order to consider interesting frequencies of the GW signal, we fix f = 10 8 GeV. For f 10 9 GeV, the transition necessarily requires very small values of β/H * to be detectable by LIGO and/or ET.
We vary λ φ and κ 1 while requiring that tunneling along the φ direction does not occur until at least T h 1 . We find that this condition is respected for any value of κ 1 , as long as λ φ 0.002. For values of κ 1 close to the lower bound κ c = 2µ 2 1 /f 2 0.02, the local minimum B appears at T h 1 2 · 10 7 GeV, while T φ ∼ 5 · 10 7 GeV.
We show the evolution of the latent heat parameter α for temperatures below T h 1 in Fig. 6 for representative choices of parameters λ φ and κ 1 . It is clear that α 1 can be obtained with these choices of parameters if there is just a mild cooling of ∼ 20 %, i.e., if T n 0.8 T h 1 . Alternatively, one can consider smaller values of µ 1 , µ 2 and λ φ , according to (3.13). In this way T h 1 can be made smaller, therefore ensuring that values of α above one are obtained even when the universe immediately tunnels below T h 1 .
Tunneling from B to A is numerically investigated by means of the multi-field tunneling package AnyBubble [31]. We find, as expected, that O(3) bubbles only provide a closed window for tunneling to occur: namely, the tunneling action S 3 /T initially decreases as the difference in vacuum energy of the two minima slightly increases (because the PQ minimum becomes deeper), then reaches a minimum value after which it grows again rapidly (because ∆V T remains constant (and then S 3 ≈ constant), while the temperature keeps decreasing (and then S 3 /T becomes larger)). For values of λ φ and κ 1 close to the line determined by the upper bound (3.13), we find that tunneling occurs very rapidly below T h 1 , with α 0.2 and β/H * 10 2 , as expected since in this region the tree-level barrier is small. However, as we move away from this limit, we find points in parameter space where T n 1.5 · 10 7 and α 1. For these points, we also find β/H * 100, since the transition occurs only after some cooling. These values are enough to make the associated GW signal detectable at ET independently of the main source of GWs and even at design LIGO, if sound waves are the dominant source of GWs. While we leave a detailed numerical scan of the values of β/H * in the parameter space of the model for future work, we expect that small regions with β/H * 10 should arise as we move further away from the upper bound (3.13), close to the region in which the universe remains stuck in B forever. 4 This would open up the possibility to detect the signal at LIGO, independently of the specific source of GWs.
In this latter respect, our two-step PQ phase transition may be characterized by a further peculiarity. Indeed, for α 1 it is not clear whether bubbles can achieve a runaway regime, nor whether the main source of GWs is the collisions of the walls or the sound waves in the thermal plasma, or in fact an admixture of both. Since our transition involves the EW gauge bosons, one should consider the implications of transition radiation [16] as these particles change mass across the bubble walls. However, in our case the EW symmetry is initially broken at B, with gauge bosons receiving masses m W ∼ µ 1 in the second-order transition from the origin to B. In the first-order transition from B to A the gauge bosons become light, which is the opposite of the case discussed in [16]. Therefore, in our case it should be possible for bubbles to run away even if they are surrounded by a thermal plasma, which would lead to v w 1 and a GW signal sourced by both sound waves and bubble collisions. Having an early phase of broken EW symmetry, with very massive gauge bosons at high energies, may also lead to interesting possibilities for baryogenesis at high scales. We leave the interesting questions above for future work.

Supersymmetric versions
A possibility to have the EW scale naturally smaller than F a without fine-tuning (and also F a M P ) is to supersymmetrize the above models. For the KSVZ models this implies that the interactions of Φ with the quarks Q L,R must arise from the superpotential term (for n = 1) while for DFSZ models Notice that in this latter case, when Φ gets a VEV, (3.25) generates a supersymmetric mass for the Higgs doublets. Since this mass must be of order the EW scale, this requires κ ∼ TeV/F a , making this term irrelevant in the scalar potential. The above superpotentials, however, leave the VEV of Φ undetermined. The latter can be generated once we add soft supersymmetry breaking (SSB) terms, which are also required to get realistic models for the EW scale. The relevant potential for φ is then simply given by 5 where m 2 φ (φ) is the SSB mass of φ and its dependence on φ arises from loop effects. The potential (3.26) can lead to a nonzero minimum for φ, similar to the Coleman-Weinberg model, by demanding that m 2 φ is positive at large φ but "runs" towards negative values as φ decreases. The VEV of φ then occurs at around m 2 φ ( φ ) ∼ 0, or, more precisely, at where β m 2 φ = dm 2 φ /d ln φ arises at the quantum level and it is then one-loop suppressed. For example, from the interaction (3.24), we have where m Q L,R and A y Q are respectively the SSB mass of the scalar component of Q L,R and the trilinear SSB term. It is easy to choose the SSB parameters such that the minimum of the potential (3.27) occurs at the desired value φ = F a .
Let us consider the phase transition of this model. At high temperatures the potential is given by where D φ is defined in (A.6). 6 The critical temperature is at where m 2 φ,min corresponds to the minimal value of m 2 (φ). As long as this minimal value is negative and occurs at φ > 0, as we will assume from now on, the potential at T c will have a thermal barrier, and a first-order phase transition will be possible. We can estimate the bounce action of a thermal O(3)-symmetric bubble as [32] where the minimization is over the tunneling point φ tun . The latter in this case corresponds to the smallest possible φ tun , determined by V (φ tun ) ≈ V (0): Since we have assumed that |m 2 φ (φ)| decreases with φ after it has reached |m 2 φ,min |, φ tun also decreases as T drops. Therefore S B decreases till it reaches S n where bubbles form and complete the phase transition. We can estimate the resulting value of α as 33) and the value of β/H * as β H * which lies close to the LIGO and ET range for interesting values of F a . Nevertheless, the predicted values of β/H * from (3.34) are quite large, 100, which makes it impossible to be seen at LIGO, since bubble collisions would be the main source of GWs in this case, and only ET could be able to detect this type of phase transition -see Fig. 2

Strongly-coupled PQ models
After discussing the possibility of a first-order phase transition in the KSVZ and DFSZ models, let us now move to a different class of realizations of the PQ mechanism. We consider the case in which the PQ symmetry arises as an accidental global symmetry of a new strong sector that, similarly to the U (1) A in QCD, is broken at the scale where condensates are formed. This scale can be chosen to be of order F a .
GWs can arise in this case from the deconfined-to-confined phase transition which proceeds in the following way. At high temperatures (T F a ) the strong sector is expected to be in a deconfined phase, where the constituents are not confined into hadrons. As the temperature drops below T c ∼ F a , the confined phase becomes energetically favorable, and the model can go through a phase transition. For a gauge theory with a large number of colors N , this phase transition is expected to be of first order, and indeed this can be proven to be the case for holographic models [33,34]. To address this phase transition quantitatively, we will follow the strongly coupled models studied in Ref. [35][36][37] which have a weakly-coupled five-dimensional version via holography. This helps to reduce the number of parameters, although the conclusions can be extended to models without holographic versions [36]. The requirements for the strongly-coupled PQ model are the following. We assume that the strong sector has a global U (1) P Q ⊗SU (3) c symmetry with an U (1) P Q −SU (3) c −SU (3) c anomaly (this means that its constituents must be colored under SU (3) c ). We also assume that the confinement scale Λ c is determined by a potential for the dilaton µ given by where the dependence of the quartic coupling λ(µ) on µ is dictated by the explicit breaking of scale invariance (several examples are given in [36]). We identify the mass gap Λ c with the dilaton VEV, µ = Λ c . We further assume that confinement also leads to the spontaneous breaking of U (1) P Q . The axion is then the corresponding (composite meson) Nambu-Goldstone boson. 7 The U (1) P Q − SU (3) c − SU (3) c anomaly guarantees the coupling (3.1), with an axion decay constant 37) where N 1 plays the role of the number of "colors" of the strong sector. The free-energy of the unconfined phase is given by F dec −π 2 N 2 T 4 /8, while in the confined phase F conf = V eff ( µ ). Thus, the critical temperature at which the confined phase is energetically favorable follows as [36] T c 0.3 × 10 10 GeV (Λ c m dil ) 1/2 10 10 GeV , (3.38) where m dil is the dilaton mass. The rate of the phase transition from the unconfined to the confined phase is in most of the cases dominated by vacuum tunneling whose bounce action is roughly given by [36] S B ∼ 24N 2 |λ(µ tun )| , (3.39) where µ tun T Λ c /T c . We are interested in phase transitions with large values of α and small values of β/H * , as this maximizes the GW strength. As in the case studied in Sec. 3.2.1, this arises when there is a period of supercooling, which in this case happens when the universe stays for a while in the unconfined phase before the phase transition takes place. In order to achieve that, |λ(T )| must slowly increase as T decreases, so that S B slowly approaches S n .
In this case we have α 1 while where β λ = dλ/d ln µ. From this, we see that long periods of supercooling, where S B evolves slowly towards S n , can give rise to small values of β/H * . This can be appreciated in Fig. 7, where we consider λ(µ) = b 0 (ln(Λ c /µ) − 1/4) and vary b 0 , or equivalently, T n . Starting at T n = 0.02 Λ c and going to smaller values, we move from the right to the left along the black solid lines of Fig. 7 (taking N = 3 and choosing different values of F a ). The value of T * is the reheating temperature after the phase transition is completed which is found to be T * 1.8 √ N /g 1/4 * T c [36]. Using this and (3.37), we obtain the relation We see from Fig. 7 that the phase transition of this model can be detected by LIGO (current and design sensitivity) if there is enough supercooling. The smaller is F a , the more likely is the detection of the GWs. Even though this scenario realizes supercooling, which strongly dilutes the thermal plasma around the bubbles and leads to v w 1, it is possible that sound waves and turbulence are still the main source of GWs. This is important because in this case detection could be easier, as can be appreciated in Fig. 2. The reason for this is that the deconfined-to-confined phase transition involves gauge bosons (dark gluons) which receive a mass across the bubble walls. As pointed out in [16], these can be radiated off as particles cross the bubble walls. This so-called transition radiation generates friction on the motion of the bubble walls and can halt their acceleration. More concretely, transition radiation leads to an upper bound on the γ factor of the bubble walls, given by [36] γ c ∼ Λ c T n 3 .

(3.42)
If bubbles collide significantly after reaching γ c , then most of the energy available in the phase transition goes to the thermal plasma, since the bubbles are not in the runaway regime even if v w is very close to one. However, bubbles can also collide before they have time to reach γ c . In this case, bubble collisions are the dominant source of GWs. Let us then estimate the amount of supercooling required to be in this latter regime. Following [36], the maximal γ factor achieved before collision is Matching the equation above to (3.42) we obtain Thus we see that for F a ∼ 10 8 − 10 10 GeV, sound waves and turbulent motion in the plasma are expected to be the dominant source of GWs when T n 10 −2 − 10 −3 F a . For longer supercooling, bubble collisions are the main source instead.
Finally, let us conclude this subsection by noting that in principle an alternative option for a long period of supercooling is to have λ(T ) evolving too slow (for a holographic example see [41])) such that the condition Γ H 4 is not met and the universe gets trapped in the unconfined phase. As discussed in [35,36], the universe could still exit supercooling at the QCD scale, where a new contribution to the dilaton potential arises. In order for this to happen, we need the strong sector to have constituents which are charged under SU (3) c . This is indeed the case for the axion models discussed here, since, as we have mentioned, the strong sector must have an SU (3) c symmetry in order for the axion to couple to GG. Nevertheless, exit due to QCD effects is not possible here since F a is much larger than the scale where QCD becomes strong, and to exit supercooling at such low temperatures, S B would need to be of order one.

Conclusions
We have shown that LIGO has the possibility to detect GWs arising from a phase transition which occurs in the early universe at temperatures around 10 8 GeV. As shown in Fig. 3, however, detection requires the phase transition to be strong enough with values of the latent heat parameter α > 1. For these types of phase transitions LIGO will be able to detect GWs for values of the inverse transition time β/H * up to ∼ 10 3 . On the other hand, the proposed ET observatory will be able to access phase transitions with slightly smaller values of α but with much larger β/H * . In particular, as shown in Fig. 2, ET will access phase transitions with α 0.1, and β/H * 10 6 .
The breaking of the PQ symmetry, required in QCD axion models, is a particularly well motivated example of such a phase transition. Indeed, the PQ phase transition would have to occur at temperatures T ∼ 10 8 −10 12 GeV, if the initial axion misalignment is not tuned to small values. The main message of this work is that LIGO, at current and design sensitivity, will be able to probe some of the simplest realizations of the PQ mechanism.
In particular, we have shown that DFSZ realizations have the right ingredients to generate a GW signal, which is in the reach of LIGO. This occurs when the PQ symmetry breaking is of Coleman-Weinberg type, that is when the mass parameters of the model are small and the minimum is generated by quantum effects. Our key results are presented in Fig. 5, which shows that PQ scales up to F a 10 11 GeV can be probed by LIGO and even more by ET. Furthermore, we have discussed an alternative type of phase transition in the DFSZ model, which is due to a zero-temperature tree-level barrier. This would proceed via an intermediate step where the EW symmetry is broken at high scales, before tunneling from this to the PQ broken phase. We have shown that this case can exhibit α 1, while the typical values of β/H * make its GW signal suited for detection at ET. A more detailed investigation of the parameter space which allows for a detectable two-step PQ transition is left for future work, as are also the phenomenological implications of the associated high-scale breaking of the EW symmetry.
For KSVZ realizations, we have shown that the simplest model does not lead to a strong first-order phase transition. However, supersymmetric KSVZ and DFSZ models can exhibit a first-order phase transition, with naturally small mass scales. We have found that the PQ symmetry breaking can be driven by supersymmetry-breaking effects, giving a first-order phase transition with α 1 and β/H * 100. We have continued our exploration of PQ phase transitions by considering models where the symmetry is broken by strong dynamics. In this case supercooling arises rather generically. The transition from the unconfined to the confined phase in these realizations can be strong enough to give a GW signal detectable at LIGO. Our key results for this type of phase transition are presented in Fig. 7.
Interestingly, other proposed observatories, like DECIGO [42] and BBO [43], would be able to probe the small frequency tails of the broad GW spectra generated by the strongest first-order phase transitions which we have discussed in this work. Looking further into the future, GW detectors with sensitivity at higher frequencies than LIGO and ET, such as [44], will open the possibility to discover phase transitions from QCD axion models with F a up to 10 11 GeV and weaker than the ones that we considered here.
Overall, as laboratory experiments progress in their search for the QCD axion at low energies, we have shown that LIGO can already join this effort by hearing the axion 'birth' at the very high PQ scale.
Note added: While preparing this manuscript we became aware of the work of [45] which also considers models with a PQ phase transition detectable at LIGO. phase transitions is the dependence of the field masses on the values of the scalar fields {φ i }, which is usually of the form m 2 a ∼ c + bφ 2 i , with c and b constants. For the scalar fields, the masses m 2 a are to be taken in the mass eigenstate basis, i.e. they are the eigenvalues of the i × i-dimensional mass matrix obtained from the tree-level scalar potential.
The tree-level zero-temperature potential receives the following corrections: 1. Coleman-Weinberg: at zero temperature, the one-loop correction to V 0 ({φ i }) using dimensional regularization and the MS renormalization scheme is given by: Here F = 1 for fermions and F = 0 for bosons. Similarly, c a = 3/2 for scalars and fermions and c a = 5/2 for vectors.

2.
Thermal: at finite temperature T , the one-loop thermal correction to V 0 ({φ i }) is given by: Here the functions J B/F are defined as Eqs. (A.4) and (A.5) deliver an important message for phase transitions driven by thermal corrections: since m 2 a ∼ c + bφ 2 , the leading thermal corrections due to bosons take the form Both fermions and bosons can contribute to D φ . On the other hand, only bosons can contribute to E φ and induce a cubic term in φ. This latter term is important, since it can induce a barrier separating two minima in field space. A further more subtle point is related to the infrared singularity in the high temperature limit of V T [46][47][48], as defined in (A.2). The standard strategy to avoid this problem is to replace the bosonic squared masses m 2 i with the dressed squared masses m 2 i (φ j ) + 2D φ i T 2 (before diagonalization of the scalar mass matrix), where D φ i = 2[∂ 2 φ i V T /T 2 ] φ i ,T =0 . This replacement is done everywhere in V T as well as in V CW . These so-called daisy corrections generically weaken the strength of a phase transition, since at high temperatures T 2 m 2 a , they screen the field dependence of the leading order cubic terms in the bosonic thermal potential.
For reference, let us conclude this section by providing the expressions for the daisy masses of the real, U (1) EM -neutral components of Φ, H 1 and H 2 which we have used in our work (we do not list the daisy masses of the imaginary and charged components, while those of the EW gauge bosons can be found in [28]):