QCD-Collapsed Domain Walls: QCD Phase Transition and Gravitational Wave Spectroscopy

For a discrete symmetry that is anomalous under QCD, the domain walls produced in the early universe from its spontaneous breaking can naturally annihilate due to QCD instanton effects. The gravitational waves generated from wall annihilation have their amplitude and frequency determined by both the discrete symmetry breaking scale and the QCD scale. The evidence of stochastic gravitational waves at nanohertz observed by pulsar timing array experiments suggests that the discrete-symmetry-breaking scale is around 100 TeV, assuming the domain-wall explanation. The annihilation temperature is about 100 MeV, which could naturally be below the QCD phase transition temperature. We point out that the QCD phase transition within some domains with an effective large QCD $\theta$ angle could be a first-order one. To derive the phase diagram in $\theta$ and temperature, we adopt a phenomenological linear sigma model with three quark flavors. The domain-wall explanation for the NANOGrav, EPTA, PPTA and CPTA results hints at a first-order QCD phase transition, which predicts additional gravitational waves at higher frequencies. If the initial formation of domain walls is also a first-order process, this class of domain-wall models predicts an interesting gravitational wave spectroscopy with frequencies spanning more than ten orders of magnitude, from nanohertz to 100 Hz.


Introduction
It has been pointed out a half century ago that spontaneous breaking of discrete symmetries can lead to productions of domain walls in the early universe and overclose the universe for a high symmetry-breaking scale [1].Additional ad-hoc explicitly symmetry breaking operators are usually introduced to bias the potential energy in different domains and collapse walls.One more natural way to collapse the domain walls is to use the known QCD instanton effects in the Standard Model (SM), which was pointed out in Ref. [2].The minimal assumption is that the discrete symmetry is anomalous under the QCD interactions such that the QCD instanton effects generate an effective potential for the discrete-symmetry order parameter, explicitly break the symmetry and collapse the walls at around the temperature of QCD phase transition.
Discrete symmetries are ubiquitous among particle physics models.For instance, the spontaneously breaking of time-reversal was studied by T. D. Lee a long time ago [3].More recently, discrete matter symmetries are introduced in the Nelson-Barr mechanism [4][5][6] to solve the strong CP problem [7,8].To explain quark and lepton masses and mixings, many discrete Abelian and non-Abelian flavor symmetries have also been introduced to achieve certain matrix structures (see [9] for instance).Some of those discrete symmetries could be anomalous under the QCD interactions.Early literature about general discrete symmetry anomalies can be found in Refs.[10,11], while more specific studies related to flavor symmetries can be found in Refs.[12,13].
In this paper, we want to point out an interesting feature about the QCD-anomalous discrete symmetries.The QCD dynamics inside different domains can be thought of as different QCD's with different θ angles (the strong CP angles).The domains with θ = 0 are energetically preferred if a Nelson-Barr-like mechanism exists to enforce zero θ angle before the discrete symmetry breaking (the QCD-anomalous discrete symmetry could be embedded in the Nelson-Barr mechanism).The finite-temperature QCD phase transition for θ = 0 is a crossover one [14,15].However, the QCD phase transition with a θ angle close to π could be a first-order one.One hint about this possibility is the first-order phase transition along the θ direction at θ = π [16][17][18][19].Another hint is based on the phenomenological linear sigma model coupled to quarks (LSMq) [20], in which the finite-temperature QCD phase transition is shown to be a first-order one at θ = π [21,22].Although a robust answer for the phase diagram in θ − T requires some non-perturbative tool, we will extend the two-quark-flavor study in Ref. [21] to the more realistic three-quark-flavor case and demonstrate a region of θ in [θ c , 2π − θ c ] centered at π to have a first-order QCD phase transition.The existence of first-order QCD phase transition could have many phenomenological consequences including the formation of quark nuggets [23][24][25][26].
On the gravitational wave (GW) side, the annihilation of domain walls can generate stochastic gravitational wave background (SGWB) with its amplitude and frequency determined by both the discrete-symmetry-breaking scale f or the wall tension σ ∼ f 3 as well as the socalled potential bias parameter that measures the potential energy differences between different domains [27].The models with QCD-anomalous discrete symmetries are more economical or predictive because the QCD intanton effects are the source of the potential bias parameter and have the known contributions at the QCD scale O(100 MeV).In our study here, we will use Z N as an example to study the domain-wall evolution as well as the annihilation-generated GW.Other than studying the scaling case for domain wall evolution in the radiation-dominated universe [28,29], we also consider the case with a domain-wall-dominated universe and the corresponding GW.In our study, we pay some special attention to the cluster of potential bias values when N > 2 such that different domain walls annihilate at different temperatures and could generate a GW spectrum different from the minimal one with N = 2 and one single annihilation temperature.Some of the GW spectrum features can be used to identify the group theory properties of the discrete symmetry and the detailed effective potential in terms of the domain-wall order parameter.
On the GW experimental side and in the last few years, three of the current pulsar timing array (PTA) experiments NANOGrav [30], EPTA [31], and PPTA [32] had reported strong evidence for a common-spectrum red process across pulsars in their data.Their results were also confirmed by the IPTA collaboration [33], combining the data sets of the three collaborations.More recently, an evidence of the key Hellings-Downs correlation [34] to confirm the gravitational wave origin of the red spectrum has been shown at around 3σ confidence level by the NANOGrav collaboration [35], which is also supported by EPTA [36], PPTA [37] and CPTA [38].Other than the explanation of astrophysical super massive black-hole binaries, many new physics model explanations have also been proposed including domain-wall collapse [39][40][41][42][43][44][45][46].
The NANOGrav collaboration has also presented their results for the domain-wall explanation of GW sources [47].For the domain walls with the discrete symmetry anomalous under QCD, the preferred symmetry breaking scale f is around 100 TeV with the annihilation temperature around 100 MeV, which is close to but below the QCD phase transition temperature.According to the phase diagram calculated in this paper based on the LSMq model, some order-one fraction of domains have a first-order QCD finite-temperature phase transition.This means that the PTA data implies an important consequence of the strong dynamics during early-universe evolution, which could also have many other phenomenological consequences.
Our paper is organized as follows.In Section 2, we introduce a discrete Z N symmetry that is anomalous under QCD.Section 3 contains the domain wall evolution of both the scaling behavior and domain-wall-dominated cases.The QCD phase diagram in θ − T is presented in Section 4 with the detailed calculation in Appendix A. The gravitational wave signatures from domain-wall annihilation, QCD and discrete-symmetry phase transitions are calculated in Section 5. Some formulas for GW spectra from phase transition are kept in Appendix B.

QCD-anomalous discrete symmetry
For simplicity, we consider the domain walls related to the spontaneous breaking of a Z N symmetry.Introducing a complex scalar field S transforming as S → e i 2π/N S under Z N , a Z N -invariant (non-)renormalizable potential, V (S), exists to determine the N -fold vacua after the symmetry breaking The discrete symmetry breaking scale f will be assumed to be much higher than the QCD or the electroweak scale in this study.For the simplest Z 2 case, the order-parameter field could be a real scalar field with a simple renormalizable potential V (S) = λ 4 (S 2 − f 2 ) 2 .For N > 2, the single-field potential could be V (S) = −m 2 SS † + λ(SS † ) 2 − µ(S N + S †N ) with the mass dimension of µ as 4 − N .For N > 4, this potential contains non-renormalizable terms and could be replaced by a renormalizable model with additional fields charged under Z N .
For the class of models with the Z N symmetry anomalous under the QCD interaction, the effective interaction below the scale f but above the QCD scale is Here, G µν =1 2 ϵ µναβ G αβ with G αβ as the gluon field tensor; θ 0 is the UV QCD θ angle that will be set to be zero later; 1 C(r ψ ) is the Dynkin index for the representation of a chiral fermion ψ under SU (3) c with C(3) = 1/2; q ψ is a (mod N ) integer and is the Z N charge of the heavy chiral fermion ψ that obtains a mass after Z N breaking.Note that, in different domains with different ⟨S⟩ j , the effective QCD θ angle is θ j = θ 0 + ψ 2 q ψ C(r ψ )2πj/N .For θ 0 = 0 and n f number of ψ fermions with q ψ = 1 and in 3 of SU (3) c , one has θ j = 2π j n f /N .To have all domain walls collapsed by the QCD instanton effects or distinct θ angles for different domain numbers j's, a necessary condition is gcd (n f , N ) = 1.Otherwise, the discrete symmetry is only broken to Z gcd (n f ,N ) by QCD with remaining uncollapsed domain walls affecting BBN observables.In the following, we will simply take gcd (n f , N ) = 1 with QCD breaking all Z N symmetry such that different domains have θ angles as After QCD phase transition, the QCD instanton effects generate an effective potential for arg (⟨S⟩ j ) or a biased potential for different domains.At T = 0 and the leading order in chiral expansion, the two-flavor potential is [50] Note that the above formula is valid in the small effective θ angle limit or when j ≪ N .For a large θ angle or j/N ∼ 1, the QCD vacuum deviates dramatically from the one used by the chiral Lagrangian (for instance, the operator G µν G µν may develop a large vacuum expectation value, which is absent in the ordinary chiral Lagrangian vacuum).Although a reliable effective potential in θ j requires a non-perturbative derivation and is absent at the current moment, we will use the above potential to guide us through the qualitative evolution of domain walls.As a comparison, we also present the effective potential in the LSMq model in Appendix A.
Taking the quark mass ratio z ≡ m u /m d = 0.49 [51] in the MS scheme at 2 GeV, m π = 135 MeV and f π = 92 MeV, the maximum potential difference among different domains is and a smaller N -dependent V max bias when N = odd.Note that the above effective potential provides the leading-order value of the topological susceptibility χ 1/4 top = 76.4MeV and is close to the NLO [52] and Lattice QCD results [53].
Because the j = 0 or θ j = 0 domain has the lowest effective potential, other domains with a nonzero j will eventually disappear with the corresponding walls collapsing at a temperature (potentially) below the QCD phase transition temperature.Different walls could have different biased potential on their two sides and therefore could collapse at slightly different temperatures, which depend on N and the detailed effective potential.

Domain wall evolution
The time-independent domain wall solutions are related to the topological structure of the symmetry breaking.For the simplest Z 2 case, an analytic solution exists for the domain wall

√
λf 3 , which we will treat as a model parameter.In the following study, we will consider a high inflation scale with the reheating temperature T 0 rh higher than the symmetry breaking scale or T 0 rh ≫ f .A thermal phase transition for the discrete symmetry happens at a temperature T form ∼ f to form the domain walls via the Kibble-Zurek mechanism [54,55].Therefore, we take σ and T form as two model parameters for both the Z 2 and the general Z N cases.

Overview of domain wall evolution
Before we discuss the domain wall evolution, we first list a few critical early-universe moments that are labeled by the corresponding temperatures.
• T 0 rh : the reheating temperature immediately after inflation.• T form : the initial domain wall formation temperature.Depending on the dynamics related to the order-parameter field S, a nontrivial first-order or second-order phase transition could happen at this time.It is around the discrete symmetry breaking scale or T form ∼ f .
• T scal : the time when domain walls reach the scaling behavior during a radiation-dominated universe.After this time, the domain walls have a fixed co-moving number density.The domain wall curvature or the averaged separating distance L scales linearly in time t.
• T dom : the time when the domain walls dominate the total energy of the universe.Note that if the domain walls collapse before they dominate the universe, there is no domainwall-dominated period.
• T caus : the time when the wall separation speed reaches the causal limit or the speed of light.Any walls survive after this time will not collapse.
• T j QCD : the QCD phase transition temperature.Note that different domains labelled by j have different θ j 's and hence different T j QCD 's as well as different orders of QCD phase transition.
• T ann : the domain wall annihilation time.Due to the biased potential among different domains from QCD instanton effects, the vacuum energy drives the domain walls to collapse.After this time, only the j = 0 or θ = 0 domain survives to the current universe.For N > 2, a set of T ann 's are anticipated because of different values of biased potential.
• T w rh : the "wall reheating temperature" after the domain walls collapse and convert their energy to radiation.
Depending on the relation between T dom and T ann , two distinct situations can happen.When T ann > T dom (denoted as Case I), the domain walls disappear before they dominate the universe.The whole domain wall evolution is approximately in the radiation-dominated universe and can be described by the scaling behavior till the annihilation time.The evolution of temperature as a function of time is plotted in the left panel of Fig. 1, where the rough scales for different temperatures are provided with f ∼ 10 5 GeV.For this example, the annihilation temperature is comfortably higher than the BBN temperature, so the BBN observables will not be affected by the early existence of domain walls.Also note that the wall reheating temperature is only slightly higher than T ann because the energy contained in the walls are subdominant compared to the main radiation energy.
For the other case with T ann < T dom (denoted as Case II and illustrated in the right panel of Fig. 1), the temperature changes its dependence on time from the early T ∝ t −1/2 in the radiation-dominated universe to T ∝ t −2 after domain-wall dominance.After a period of a quick drop of temperature that reaches T ann , the collapse of domain walls reheats the universe to T w rh , which could be much higher than T ann .For this case, T ann could be lower than T BBN , but with T w rh above T BBN to guarantee a radiation-dominated universe before the BBN time.
Figure 1: Schematic plots of early universe temperature evolution as a function of time.The left panel depicts purely radiation-dominated universe, while the right panel contains a short period of domain-wall dominance.The QCD phase transition time is labelled by the red star with T j QCD denoting a set of QCD phase transition temperatures of different θ domains.The purple cross in the right panel indicates the possible causal-limit temperature T caus , which is not reached in this example.

Details of domain wall evolution
In this section, we detail the domain wall evolution through cosmic history.We consider T 0 rh > T form , so that domain walls form and survive in the radiation-dominated universe after inflation and reheating.For T < T 0 rh , we have the radiation energy density ρ R (T ) = (π 2 /30) g * (T )T 4 and the Hubble scale H(T ) = (π 2 /90) 1/2 g * (T ) 1/2 T 2 /M pl , where g * (T ) denotes the total radiation degrees of freedom and M pl = 1/ √ 8πG = 2.43 × 10 18 GeV is the reduced Planck mass with G as the Newton constant.The scale factor during the radiation-dominated universe scales with time as a ∝ t 1/2 , and thus H(t) = 1/(2t) and T ∝ t −1/2 .We will remain agnostic regarding the order of the phase transition for the discrete symmetry breaking, as the information about the initial size of the domain wall population determined by the Kibble-Zurek mechanism [54,55] will be erased by the subsequent evolution.
The evolution of domain walls is governed by two forces: the tension force p T = ρ w = σ/L, where L is the curvature radius of the walls and controls the average distance between domain walls, and the friction force due to the reflection of plasma particles off domain walls.The wall network evolution, taking into account effects of the Hubble damping and the friction force, is determined by velocity-dependent one-scale (VOS) model [56] which agrees with the numerical simulation.In this model the root mean square velocity of the wall v and L are related by the coupled differential equations as follows [56]: where c w and k w are constant phenomenological parameters determined by the numerical simulation and l d is the damping scale given by l −1 d = 3H + l −1 f .Here l f is the friction length.If the friction force is negligible, which would be the case if the relativistic plasma particles perfectly transmit through the walls, then l −1 d = 3H (see [57] for the case with dominant friction term and its effect on domain wall evolution).In such cases the scaling solution with L = L 0 t and v = v 0 is reached where L 0 and v 0 are constants determined by c w and k w .For the case of a radiation-dominated universe, numerical simulations suggest that the scaling solution is reached at the temperature T scal ≈ T form /30, with L 0 = 1.2 and v 0 = 0.42 (which give k w = 0.66 and c w = 0.81) [56].While this result is only valid for Z 2 , for general Z N the scaling solution is also seen in numerical simulations with ρ w = A σ/t and A ≈ 0.4N [58,59].
As the domain wall energy density in the scaling regime goes as ρ w ∝ t −1 , whereas the radiation energy density during the radiation-dominated era goes as ρ R ∝ t −2 , domain walls will eventually dominate the total energy density of the universe.This happens at the temperature T dom that is determined by ρ w = ρ R and given by Parametrically one has t dom = 3M 2 pl /(4Aσ).As domain walls cannot dominate the energy density of universe till today or even BBN time, they would have to annihilate.There are two cases: I) walls annihilate before they dominate the total energy density (T dom < T ann ); II) walls annihilate after they dominate the energy density of the universe (T dom > T ann ).We will consider both cases below.

Case I: T dom < T ann
In this case the radiation energy dominates the universer until the walls annihilate.Thus, the scaling solution described above, with ρ w ∝ t −1 , is valid till the annihilation time.As different domains, with different θ j 's, have different potential energies as given by Eq. ( 4), this leads to a vacuum pressure force p V = V bias , with V bias being the vacuum energy difference between adjacent domains.Note that depending on the form of the effective potential, there could be a series of V ij bias labelled by two adjacent domain numbers.In this section, we keep V bias as a general parameter and will come back to the indexed V ij bias later.As p T tries to stretch the walls as depicted by the scaling solution, p V acts to collapse the domains with higher vacuum energies.The collapse starts when p T = p V , which corresponds to the temperature Parametrically one has t ann = A σ/V bias .The larger the surface tension, the longer it takes to collapse the walls, while the larger V bias , the earlier annihilation happens.For a higher symmetry breaking scale or larger σ, it takes longer to annihilate the walls and thus leads to a lower T ann .
After annihilation the energy contained in the domain walls gets transferred to the radiation, reheating the universe to T w rh given by where F is given by This implies that for most of the parameter space T ann ≈ T w rh for Case I. Given the tight constraints from BBN observables, we impose a constraint of T ann ≈ T w rh ≳ 3 MeV [60,61].Using Eq. ( 8), this translates into an upper bound on σ as Later we will show that this bound is practically irrelevant because domain-wall dominance and/or the causal-wall-separation limit will be reached before one reaches the above upper limit.
Note that we have used temperature-independent V bias so far.In fact, for the potential in Eq. ( 4) the thermal corrections are known for T < T j QCD [52,62,63] based on the chiral Lagrangian of the θ = 0 vacuum (for the GW spectra presented in Sec. 5, we will check the finite-temperature effects).For T > T j QCD the θ-dependent potential is not known reliably.The dilute instanton gas calculation [64], relying on the finite-temperature perturbative QCD, is only valid for T > 10 6 GeV [52] [see [2,65] for the domain wall annihilation assuming the validity of dilute instanton gas calculation up to O(GeV)].In this paper we will consider the case with T ann < T j QCD as we will be immune to the above uncertainties and it also opens up a new phenomenological study of QCD phase transition at finite QCD theta angle (see Sec. 4 for details).For the case with T ann > T j QCD , we would have θ = 0 in all Hubble patches at T QCD and the QCD phase transition would be a "boring" crossover one at T QCD ≈ 170 MeV [66].If we demand that T ann < T j QCD ≈ 125 MeV (see Fig. 3), we obtain a lower bound on σ σ ≥ 9.2 × 10 15 GeV 3 A 0.8 V bias (100 MeV) 4   g * (T ann ) 10 Finally, we demand T dom < T ann for this case so that domain walls annihilate under radiation domination.This imposes an upper bound on σ given by For σ larger than the above value, the universe will enter the domain-wall-dominated era before wall annihilation.
Comparing the two constraints from opposite ends in Eq. ( 12) and Eq. ( 13), one can see that for Case I there is only a small parameter space in σ that allows for potential non-trivial QCD phase transitions in different domains unless a smaller value of V bias is given.

Case II: T dom > T ann
In this case when t > t dom , the scaling solution given in the previous subsection with L ∝ t is no longer valid.In the domain-wall-dominated universe, we have H = ρ w /(3M 2 pl ) = σ/(3LM 2 pl ).Substituting it into Eq.( 5), with the boundary conditions L(t dom ) ≡ L dom = t dom /A and v ≃ 0, where we have defined the causal wall separation distance L caus = 3M 2 pl /σ, which is the average domain wall distance at the causality limit L = H −1 = L caus .The time at which L → L caus is given by t caus = t dom (8A − 3).For t > t caus , we have L > H −1 , i.e., domain walls are separated by a scale larger than the horizon size.In such a case the assumptions of Friedmann-Robertson-Walker (FRW) cosmology, homogeneity and isotropy, break down.Moreover, the vacuum energy inside the domain would likely determine the Hubble evolution and guide the universe into a period of inflation.Accordingly, the domain walls cannot collapse, and thus this scenario with L > H −1 is clearly ruled out.Hence, we demand that domain walls annihilate before the causality limit is reached, i.e., t ann < t caus .
In the domain-wall-dominated scenario, we can derive from the annihilation condition ρ w = V bias Requiring t ann < t caus then imposes an upper bound on σ Combined with Eq. ( 13), the domain-wall-dominated universe happens for σ satisfying 2.6 × 10 16 GeV 3 A 0.8 After domain walls annihilate, the energy contained in the wall network is transferred to radiation and reheats the universe to a temperature T w rh of which is obviously higher than the BBN temperature of O(1 MeV) and thus safe from the BBN constraints.
To derive the values of other characteristic temperatures, we use Eq. ( 14) and H = ȧ/a = 2/(t + 3 t dom ), which give At t = t ann , the domain wall annihilation temperature is given by Together with the condition given in Eq. ( 17), one can see that T ann < T j QCD and therefore non-trivial QCD phase transitions could happen in Case II.
At t = t caus , one has Demanding T ann > T caus gives the same upper bound on σ as in Eq. ( 16).Note that one can have T caus or T ann less than T BBN , since the universe that is reheated after wall annihilation has T w rh > T BBN .The lowest value of T caus can be calculated by taking the smallest possible value of σ [from Eq. ( 13)] that would cause domain wall domination, which is For a large A, i.e., for a large N of the Z N discrete symmetry, one can have a much smaller T low caus .
We summarize the different situations for different values of σ in Fig. 2. In this figure, we also show the preferred range of σ by the PTA data, which turns out to sit near the boundary between radiation dominance and wall annihilation before QCD phase transition.The PTA GW results prefer a value of σ that is in the radiation dominance case with the annihilation temperature comparable to the QCD phase transition temperature (see Sec. 5.4 for details).Here, V max bias is assumed to be valid up to T ann = 148 MeV.

QCD phase transition with non-zero θ
The QCD with θ = 0 in the SM has been shown to have a crossover when it transits from the high-temperature quark-gluon plasma phase to the low-temperature hadronic phase [14,15].For a non-zero θ, there is no robust calculation for the strength of finite-temperature phase transitions.On the other hand, it has been pointed out a long time ago by Dashen [16] and later by Witten [17] using Large-N c expansion that CP is spontaneously broken for θ = π with a first-order transition along the θ direction.More recently, this result is obtained and confirmed by using the generalized symmetry tool and 't Hooft anomaly matching [18,19].This result also suggests QCD works quite differently for θ = π or close to π.Without a non-perturbative tool to analyze the QCD phase diagram with a non-zero θ, one could use some phenomenological model to gain some insights about the QCD finite-temperature phase transition.In Appendix A, we adopt the so-called linear sigma model coupled to quarks (LSMq) introduced in Ref. [20].For two quark flavors, the finite-temperature phase transition with a nonzero θ has been studied in Ref. [21].A comparison between the LSMq and Nambu-Jona-Lasinio (NJL) model has been studied in Ref. [22].In Appendix A, we perform a detailed study about finite-temperature phase transitions for the LSMq model with three quark flavors and matched meson spectra.
The results in Appendix A are summarized in the phase diagram of θ − T in Fig. 3, where we show the phase transition temperatures for different θ's.For θ = 0, the crossover transition is recovered.The phase transition temperature is around 146 MeV, which is close to the Lattice QCD result (≈ 170 MeV) [66].For θ close to π, first-order phase transitions happen with a lower phase-transition temperature as θ gets closer to π.The critical points (labelled by the blue dots) are located at (θ, T ) = (0.7π, 136 MeV) and (1.3π, 136 MeV).Therefore, for θ ∈ [θ c , 2π − θ c ] with θ c ≈ 0.7π, the QCD phase transition is first-order based on the LSMq model.We want to emphasize again that the LSMq model is just a phenomenological model to provide the hint of a first-order QCD phase transition.The true phase diagram could be similar to the one in Fig. 3, although the precise locations of critical points or the value of θ c are subject to additional theoretical and numerical studies.
If T ann < T j QCD , different domains with different θ j 's could undergo phase transitions before the domain walls are annihilated due to the biased potential.For those domains with θ j ∈ [θ c , 2π − θ c ], the QCD phase transitions are first-order and could generate additional GW's beyond the one generated by domain wall collapse.Later, we will show that the PTA observed GW prefers a range of σ with T ann < T j QCD , or a non-trivial QCD phase transition.

Gravitational wave signatures
In this section we will consider GW emissions from three different sources 1) collapse of the domain wall network, 2) QCD phase transition, and 3) phase transition from the discrete symmetry breaking at T form .While the generation of GW's from the last two sources depend on the model parameters and could be absent, GW emission from domain wall annihilation is guaranteed if the PTA observed GW is due to domain-wall annihilation.As we will show, the GW's from above cases will span more than ten orders of magnitude in frequency, and could be probed with future GW experiments.

GW from domain wall collapse
As we pointed out in Sec. 3, after formation domain walls evolve to the scaling regime and annihilate due to the vacuum energy difference across the wall; the walls could annihilate under radiation domination T ann > T dom , or during the domain-wall-dominated era T ann < T dom .In this paper we will mainly focus on the former case but will also estimate the frequency and amplitude of emitted GW's in the latter case.The estimation of GW amplitude could be performed using the Einstein's quadrupole formula [67] with the gravitational radiation power as Here, Q is the transverse-traceless part of the quadrupole moment of matter.Two time-derivatives of this formula come from using the tensor virial theorem to convert the volume integration of the spatial part of energy-momentum tensor into double time-derivatives of the quadrupole moment.The third time-derivative is simply from calculating the propagating graviton field.For the domain wall case [68], one has Q ∼ M wall L(t) 2 with L being the curvature radius of the wall and the mass of the wall M wall ∼ σL 2 .The energy density of GW's released per unit Hubble time ρ GW ∼ P GW H −1 /L 3 , where the factor L −3 can be thought of as the number density of GW sources.We will evaluate ρ GW separately for the radiation domination and the domain-wall domination cases.The peak-frequency of the GW at the annihilation time, for both cases, is given by f (t ann ) ∼ H(t ann ).

Case I: collapse during radiation-dominated era
In the case of domain wall annihilation during the radiation-dominated era, we can use L(t) ∝ t to obtain P GW ∝ G σ 2 t 2 , which leads to ρ GW ∼ G σ 2 .The overall constant can be fixed by comparing to the spectrum obtained through numerical simulations [29,58,68] (see Ref. [27] for a review).In particular, the spectrum of GW's per unit logarithmic frequency interval has a peak value Here, k is the comoving wave number.The peak amplitude at t ann is then given by with ρ c (t ann ) as the critical energy density at t ann .Here, εGW = 0.7 ± 0.4 is a phenomenological parameter determined by numerical simulations [29].The peak frequency of the spectrum has f peak (t ann ) ≈ H(t ann ) based on the simulation up to N = 6 [29].
After GW productions at the domain wall annihilation time, the amplitude and frequency of the GW are red-shifted till today.The peak amplitude today (t 0 ) is given by [27] Ω GW h 2 (t 0 ) and the peak frequency is given by f peak = 1.1 × 10 −8 Hz g * (T ann ) 10 1/2 g * s (T ann ) 10 where g * s ≈ g * stands for the effective relativistic degrees of freedom for the entropy density.The numerical simulations also find that the GW spectrum scales as f 3 for f < f peak , as expected from the causality argument [69,70], and as f −1 for f > f peak [29] (the numerical simulation also shows a harder spectrum for a larger N [58]).Note that we are using T ann and σ as model parameters to describe the GW amplitude and frequency.While T ann is determined using V bias and σ (see Eq. 8), V bias (T ) could be a function of temperature.In fact, the finite temperature corrections to the potential in Eq. ( 3) are evaluated for T < T QCD and are given by [52] V ⟨S⟩ j ; T where and m 2 π (⟨S⟩ j ) = −V ⟨S⟩ j /f 2 π with m π = 135 MeV being the pion mass in the normal (j = 0) vacuum.The finite temperature correction changes V bias by around 10% for T ≈ 125 MeV, above which we expect that the QCD phase transition could change the potential drastically [52].Using Eq. ( 8), we expect a max 5% change in T ann .This will lead to a maximum of 5% and 20% change in the peak frequency and amplitude, respectively.Given the larger uncertainty in other parameters (ε GW and A), we will neglect the uncertainty from finite temperature corrections.As a result, we will use the T = 0 potential in Eq. (3) to calculate V bias .
Note that while the spectral characteristics of GW from domain wall annihilation mentioned above were obtained for the Z 2 case, for Z N the results should be obtained with some slight modifications.For the general Z N case, the energy difference V ij bias between two adjacent domains (labelled by θ i and θ j ) could be the same or different.Similar situations are also pointed out in [71] and [72] under the context of U (1) N PQ and Z N symmetries, respectively, in the latter of which the Z 3 case is studied in detail.If V ij bias are all the same, all domain walls collapse at around the same time with a common T ann .Then, there is a single peak in gravitational wave spectrum.In this case, according to numerical simulations [58], the peak amplitude and frequency match with Z 2 .While the spectrum at low f < f peak is same as the Z 2 case, for f > f peak the amplitude is slightly enhanced relative to the Z 2 case [58].This enhancement for a large N at higher frequencies is caused by the production of many configurations with sub-Hubble sizes, which after collapse lead to higher frequency GW's [27].
For the case with various V ij bias values, which could be the case for non-linear potentials such as the one in Eq. ( 3), one has sequential annihilation of domain walls: domain walls with a larger V ij bias collapse earlier with a higher T ann .Since different T ann 's lead to different GW peak frequencies, a multi-peak GW spectrum is anticipated from the sequential annihilation of walls.The amplitude of GW produced at some T ann is determined by the fraction of the total number of domain walls undergoing collapse at that temperature.This fraction is determined by the statistical distribution of V ij bias for a given discrete symmetry Z N .The situation for Z 2 is very simple with V bias = (100.4MeV) 4 .For Z 4 , one has domains with θ j = 0, π/2, π, 3π/2 which we denote by j = 0, 1, 2, 3. Using the potential in Eq. ( 3), one has V 01 bias = V 03 bias ≈ (79.4 MeV) 4 , V 12 bias = V 23 bias ≈ (88.8 MeV) 4 , and V 02 bias ≈ (100.4MeV) 4 , where As initial domains, i.e., j = 0, 1, 2, 3, are created with equal probabilities, the statistical distribution of V ij bias for the N = 4 case can be obtained, as shown in the left panel of Fig. 4 and labelled by "χPT".For comparison, we also show the probability distributions using the effective potential constructed from the cos function: bias ≈ (100.4MeV) 4 , and the one in the LSMq model.In the right panel of Fig. 4, we show the probability distributions for the N = 6 case.For this case, there are potentially 15 V ij bias 's in total.Two of them, V 15 bias and V 24 bias , are zero and thus not included in the right panel of Fig. 4. Depending on the surrounding wall evolution, those domain walls have a more complicated collapsing possibility and are not included in our later analysis.A similar treatment is also used for the N = 4 case by not including the V 13  bias wall in the later analysis.After knowing the probability distributions of V ij bias , we can then estimate the GW spectrum.For simplicity, we just sum the GW's generated from collapsed domain walls at different temperatures according to the specific V ij bias and Eq. ( 8) with A = 1.6 for N = 4.With this assumption, we ignore the potential additional wall evolution after some fraction of walls are collapsed.Using the χPT and N = 6 model as an example and with σ = 1 × 10 15 GeV 3 , we have 1/13 of walls collapse at around T ann ≈ 187 MeV, 2/13 of walls collapse at T ann ≈ 169 MeV, 150 MeV, 112 MeV, 79.6 MeV, and 4/13 of walls collapse at T ann ≈ 127 MeV.For each fraction, the GW amplitude is calculated using Eq. ( 24) and the specific T ann including the spectrum feature of f 3 for f < f peak and as f −1 for f > f peak .Following this treatment, we show both the total and sub-components of GW spectra in Fig. 5 for N = 2 (left panel) and N = 6 (right panel) and for both χPT and Cos potentials for comparison.For the N = 6 case, one can see that the GW spectrum is dominated by the smallest V ij bias with the smallest T ann due to the larger inverse power dependence of the GW amplitude in T ann [see Eq. ( 25)].Furthermore, this spectrum also has an interesting plateau-like feature that could be used to identify the specific discrete symmetry group with a more precise measurement of the GW spectrum.
Figure 5: The GW spectra induced by the (sequential) domain-wall collapses based on the "χPT" (red) and "Cos" (green) potentials for N = 2 (left) and N = 6 (right).The dashed lines in the right panel denote the spectra induced by the collapses of individual walls based on the "χPT" potential, and the red solid line denotes the total spectrum after taking into account the statistical distribution.

Case II: collapse during domain-wall-dominated era
The GW spectrum analysis till now only considered the wall annihilation during the radiationdominated period.Now, we turn to the GW generation from wall collapses during the walldominated period.Using L(t) from Eq. ( 14), one has where κ is a constant.The energy density of GW's released per unit Hubble time is If the GW was emitted at t = t dom , then we can compare the above equation to the corresponding one from the radiation-dominated regime in Eq. ( 23) to fix κ = π 9 εGW A 7 .Similar to the case of radiation domination, the peak frequency of GW is at where we have used the fact that for instantaneous reheating H(t ann ) = H(T w rh ) with T w rh given by Eq. (18).Taking into account the redshift factor from T w rh to today, the peak amplitude is In above equations we have taken g * s (T w rh ) = g * (T w rh ).

GW from QCD phase transition
As mentioned in Sec. 4, the LSMq model suggests a possibility of QCD first-order phase transition (PT) in the regions with non-zero θ values, in particular for θ ∈ [θ c , 2π − θ c ] with θ c ≈ 0.7π (see Fig. 3) with the phase transition temperature of ≈ 125 MeV.If T ann ≲ 125 MeV, then we expect to have QCD first-order PT for the domains with θ ∈ [θ c , 2π − θ c ] and GW productions.
As we do not have a trustworthy model to describe the QCD phase transition dynamics, we will follow a model-independent approach to describe the GW spectrum with the nucleation temperature fixed at ≈ 125 MeV.Fig. 3 suggests that only in the domains with θ ∈ [θ c , 2π − θ c ] are the FOPT and thus generation of GW's enabled.For a general Z N case, as we expect to have domains with θ = 2πj/N with 0 ≤ j ≤ N − 1, the domains with ⌊N − N θ c /2π⌋ ≥ j ≥ ⌈N θ c /2π⌉ will undergo QCD FOPT.Since at the time of phase transition we expect each θ j to cover 1/N fraction of the universe, we expect FOPT in only a fraction of universe, denoted by ζ, populated by domains satisfying the above condition on j.For θ c ≈ 0.7π, one has ζ = 1/2, 0, 1/4, 2/5, 1/6 for N = 2, 3, 4, 5, 6, respectively.Those fractions will be the multiplication factors for the standard GW amplitudes from FOPT, where one has the whole universe undergoing the phase transition.The GW's from FOPT are mainly described by two parameters: α GW ≈ ∆V (T n )/ρ R (T n ) denotes the strength of the PT, and β GW denotes the rate of bubble nucleation with β GW /H(T n ) being the commonly used parameter.Here ∆V (T n ) denotes the vacuum energy difference between the false and true vacua at the nucleation temperature T n .The SGWB from FOPT can come from three different processes: collision of bubbles, sound waves generated from bubble expansion, and turbulence in plasma.The GW spectra from these three processes have been evaluated numerically (see Appendix B for formulas).Note that the ζ factor will act as a multiplicative suppression factor for the standard formulas from all three sources (see [69] for a recent review), as we have only ζ fraction of the universe undergoing FOPT.Since determining the precise values of α GW and β GW for QCD FOPT with a non-zero θ is challenging, we take model-independent approach and show the spectra for a range of parameters.In Fig. 7 we show the spectra for (α GW , β GW /H(T QCD n )) = (0.5, 10 4 ), (0.5, 10 5 ) with ζ = 0.5 and compare them with the sensitivity curves of GW experiments.

GW from potential phase transition at T form
In addition to QCD phase transition, we also expect a possible phase transition at T form leading to discrete Z N symmetry breaking and the formations of domain walls.To describe the phase transition dynamics we need to know the complete UV model leading to the effective operators in Eq. ( 2), as new particles and interactions modify the tree-level potential [73].Depending on those UV model parameters, the thermal effective potential for the field S, V (S, T ), could allow nearly degenerate minima and thus a FOPT at a temperature of T form ∼ f ∼ σ 1/3 .As discussed for the QCD case, we can then expect another source of GW's from this FOPT.In Fig. 7, we show the expected GW spectra for T form = 2 × 10 5 GeV, (α GW , β GW /H(T form )) = (0.5, 10 4 ), (0.5, 10 5 ) and compare them with the sensitivities of GW experiments.
Figure 6: The preferred parameter space in T ann and α DW based on the NANOGrav 15-year data for the domain-wall model [35,47].The two contours correspond to one-and two-sigma confidence levels, respectively.The brown (N = 2) and orange (N = 6) curves are for the QCDannihilated domain wall models with a relation between α DW and T ann after fixing V bias .The N = 2 curve intersects with the two-sigma contour at points with the σ values of (0.66, 1.2) × 10 16 GeV 3 .For N = 6, the intersections have σ = (0.9, 1.3) × 10 15 GeV 3 .The green band indicates the range of the QCD phase transition temperatures for different θ angles (see Fig. 3 based on the LSMq model).

Hints from pulsar time array experiments
The current and future pulsar timing array experiments including NANOGrav [74,75], EPTA [76], PPTA [77], MeerTime [78], CHIME [79], and SKA [80] are sensitive to gravitational waves with frequencies in the range 10 −9 Hz − 10 −7 Hz.Since collapsing domain walls with T ann ≈ 100 MeV lead to GW's with a peak frequency within the above range [see Eq. ( 26)], one can probe the QCD-collapsed domain-wall scenario with the PTAs.A few years back, three of the current pulsar timing array experiments NANOGrav [30], EPTA [31], and PPTA [32] had reported strong evidence for a common-spectrum red process across pulsars in their data.The result was also confirmed by the IPTA collaboration [33], combining the data sets of the three collaborations.While the common-red spectrum was observed in all of the data sets, the Hellings-Downs correlations [34] required to confirm the observation of GW background was lacking back then.
Recently, with a larger pulsar data set and longer experimental run-time, the Hellings-Downs correlation curve has been observed by NANOGrav at 3 sigma [35], EPTA at 3 sigma [36], PPTA at 2 sigma [37] and CPTA at 4.6 sigma [38].These results suggest the first observation of SGWB.The domain wall collapse during the radiation-domination era as a new-physics explanation for the observed GW background has been considered by the NANOGrav collaboration [47].In Fig. 6, we quote NANOGrav's one-and two-sigma contours in the T ann and α DW plane that is favored by the NANOGrav 15-year data-set.Where α DW is defined as the fraction of domain-wall energy density over the critical energy density at the time of annihilation, which  ) = (0.5, 10 5 ) (blue).Dot-dashed lines are from possible first-order phase transition of the discrete symmetry breaking with (α GW , β GW /H(T form )) = (0.5, 10 4 ) (yellow) and (α GW , β GW /H(T form )) = (0.5, 10 5 ) (purple).The gray band denotes the rough range of GW spectrum observed by NANOGrav [35]. is given by α DW = F/(1 + F), where F is given by Eq. (10).Using Eq. ( 8) and Eq. ( 10) with V bias ≈ (100 MeV) 4 for N = 2, and V bias ≈ (67 MeV) 4 for N = 6, one can obtain a relation between α DW and T ann , which is shown in Fig. 6 by the brown (N = 2) and orange (N = 6) curves.For the N = 6 case, we expect to have a sequential collapse of domain walls.The GW spectrum shown in Fig. 5 is dominated by the smallest T ann or smallest V ij bias , which from Fig. 4 is given by V 1/4 bias ≈ 67 MeV.

GW spectroscopy
In Fig. 7, we show the predicted GW spectra from domain-wall annihilation, QCD and discretesymmetry breaking phase transitions, assuming that both QCD and discrete-symmetry phase transitions are first-order.The detailed formulas for GW spectra from phase transition can be found in Appendix B. The GW spectroscopy based on the QCD-anomalous discrete symmetries span more than 10 orders of magnitude in frequencies from 10 −9 Hz to 100 Hz.For comparison, we also show the sensitivities from the current and future GW experiments: ET [81], Ad-vLIGO [82], DECIGO [83], TianQin [84], Taiji [85], LISA [86], SKA [87], IPTA [88], EPTA [89], CE [90], BBO [91].In this plot, the preferred GW spectrum range from NANOGrav results is shown in the gray band.

Discussion and conclusions
The domain-wall collapse, in addition to the production of gravitational waves, could also form primordial black holes (PBH) [92][93][94][95].The relevant quantity which determines whether PBHs are produced is p(t) = 2GM DW (t)/L(t), the ratio of Schwarzschild radius of the domain wall to the size of the wall.When p(t) ≥ 1, domain walls could collapse into black holes.Including both the vacuum energy and surface energy, the mass of the domain wall M DW ≈ 4π 3 V bias L 3 + 4πσL 2 .It is clear that a maximum of p(t) is reached when L(t) is maximum.As walls decelerate due to the potential bias term with a deceleration of V bias /σ, the maximum wall size is reached when the wall velocity becomes zero and the wall starts to shrink.The maximum wall size is estimated to be where v 0 is the initial wall velocity.The maximum value of p(t), p max , is then given by Here, we have used v 0 = 0.42 and A = 0.8, and α DW is the parameter used in the PTA result (see Fig. 6).As long as α DW < 0.33, we have p max < 1 and PBH's are not produced from domain wall collapse.For the case with α DW > 0.33, PBH's might be formed with their abundance determined by the details of the collapse process [59] including sphericity [92].
We also briefly comment on implications on quark nuggets for the QCD-collapsed domain wall scenario.The possibility of first-order phase transitions for some domains with a nonzero θ angle could fulfil the mechanism of "Cosmic Separation of Phases", as pointed out in Ref. [23].Interestingly, collapsed domain walls can serve an alternative way to form quark nuggets.As the walls collapse, the baryon number is accumulated into a small pocket with the quark degenerate Fermi pressure to withhold the vacuum pressure.For both formation mechanisms, an initial nonzero baryon-number chemical potential is required to exist from early universe.The usual calculations of quark nugget properties are based on the θ = 0 domain [96] (see also [97] for updates).The scenario studied in this paper could suggest other types of quark nuggets with a nontrivial QCD vacuum inside (the one with a nonzero θ).For the phenomenological approach by introducing a "bag parameter", one can define a generalized bag parameter: B (had) θ=0 ,(q/g) θ , which are related to the vacuum energy differences between the hadronic phase with θ = 0 and quark-gluon plasma phase with a generic θ.Based on the LSMq model, the bag parameter is calculated to be a monotonic decreasing function of θ from 0 to π with B (had) θ=0 ,(q/g) θ=0 = (222 MeV) 4 and B (had) θ=0 ,(q/g) θ=π = (219 MeV) 4 .The stability of quark nuggets with a nonzero θ requires its mass per baryon to be lighter than the nucleon mass with the same θ (or even different θ if the quantum decay into nucleons with a different θ angle is allowed; Ref. [98] has suggested that the nucleon mass decreases monotonically as θ increases).Quark nuggets with a nonzero θ could be more stable than the ordinary quark nuggets because the nontrivial topological angle can further protect them from destructions.The properties, formation and evolution of this new type of quark nuggets deserve further investigation.
In this paper, we have considered only the case with the discrete-symmetry-breaking scale f after inflation or reheating.For the opposite case, our universe could end up with a QCD with θ ̸ = 0 without a domain wall in the visible universe.Since the current universe is consistent with a QCD with (approximately) zero θ, quantum or thermal tunneling are needed to transit the vacuum to the QCD with zero θ.For a high scale f much above the QCD scale ∼ 100 MeV, those tunneling rates are exponentially suppressed and a viable scenario is hardly anticipated.For the case with f close to the QCD scale, the tunneling rates can be large enough, although additional QCD-charged fermions to mediate the QCD-anomalous discrete symmetry could have a too low mass to be phenomenologically viable.
The scenario with QCD-anomalous discrete symmetries has the interesting interplay between the effective θ angle and QCD phase transition.One may wonder whether a similar interplay exists for QCD axion models, where different Hubble patches could have different axion field values or also different effective θ angles.The existence of an interesting interplay relies on the relation of the QCD phase transition temperature T QCD and the axion oscillation temperature T osc , when the axion particle mass is comparable to the Hubble scale and the axion field starts to oscillate around the effective θ = 0 value.Based on the dilute instanton gas calculation [99,100], the oscillation temperature T osc is qualitatively higher than T QCD .Effectively, one has a zero θ during the QCD phase transition time with a crossover phase transition.Some interesting consequences for axion properties based on first-order QCD phase transitions (studied in Refs.[101][102][103]) may not happen.
In conclusion, we have demonstrated that the PTA-data preferred domain-wall models could have an interesting interplay with the QCD dynamics with a nonzero θ angle, assuming an underlying QCD-anomalous discrete symmetry.The discrete-symmetry-breaking scale is around 100 TeV with the domain-wall annihilation temperature of around 100 MeV.Some domains with an effective large θ angle could undergo a first-order QCD phase transition, which can lead to other phenomenological consequences including GW's at a higher frequency.Many future GW experiments could test the QCD-collapsed domain wall models studied here.the eigenvalues are Using the partially conserved axial current relation, ⟨σ x,y ⟩ and h x,y are given by [104] ⟨σ where f π,K are the pion and kaon decay constants, respectively.Once we have the relation between potential parameters and experimental observables, we can fix the parameters using the observed values.We choose experimental values of f π , f K , m π , m K , m 2 η + m 2 η ′ , and m σ = 600 MeV [107] to fix the six potential parameters λ 1,2 , h x,y , µ 2 , κ.In this way the remaining meson masses, such as the masses of η and η ′ , could be checked against the measured values, and are in fact found to match with the experimental values.
Here, the first two vacua with θ = 0 are CP conserving with the first one being the global vacuum.The last two vacua with θ = π are degenerate and CP-violating.
Next we want to derive the vacuum structure at finite temperature to determine the nature of finite-temperature QCD phase transition.In the LSMq [20], in addition to the potential V vac in Eq. (38), one also has the Yukawa interactions between the quarks and the meson fields given by L Yukawa = q i / ∂ − gT a σ a + iγ 5 π a q, q = u, d, s .
In fact at non-zero temperature the quark degrees of freedom are the only one giving rise the thermal potential where ν q = 24, ν s = 12, E q = p 2 + M 2 q , and E s = p 2 + M 2 s , with M q = g ⟨σ x ⟩ 2 + ⟨π x ⟩ 2 /2 and M s = g ⟨σ y ⟩ 2 + ⟨π y ⟩ 2 / √ 2. Because of the isospin symmetry, we require that the constituent light quark mass M q = gf π /2 should be roughly 1/3 of a nucleon mass, which gives g = 6.6 and accordingly the constituent s-quark mass M s ∼ 422 MeV.After fixing V T , one has the total finite temperature potential as V T + V vac .We can minimize this potential for any given T and θ value to obtain the vacuum expectation values ⟨σ x,y (T )⟩, ⟨π x,y (T )⟩ as a function of temperature to study phase transition.
In Fig. 8 we plot the thermal evolution of the condensates with θ = 0, π, respectively.One can see from the left panel that when θ = 0 there is a crossover phase transition at T ≈ 146 MeV as expected from the lattice QCD simulations.On the contrary, in the right panel with θ = π there is a first-order phase transition at T ≈ 129 MeV.The procedure can be carried out for different θ values to figure out the nature of phase transition.We use CosmoTransitions [108] to identify the transitions and show the resulted θ − T phase diagram in Fig. 3.

B Gravitational wave spectrum from phase transition
In this appendix we list the GW spectrum formulas for three processes: collision of bubbles, sound waves generated from bubble expansion, and turbulence in plasma [109][110][111][112][113]. We also take into account the ζ factor which denotes the fraction of the universe undergoing FOPT.The spectrum is given by: Ω col (f )h 2 ≈ 1.4 × 10 −13 ζ 0.5 g * (T n ) 20 where v w is the wall velocity and κ col, sw, turb denotes the fraction of vacuum energy that is converted into kinetic energy, bulk motion of the fluid, and of the turbulence, respectively.They are given by κ col = 1 1 + 0.715 α GW 0.715 α GW + 4 27 where ξ turb ∼ 0.1 is the fraction of turbulent bulk motion.The red-shifted peak frequencies of GW spectra are The red-shifted Hubble parameter to today is given by h * ≈ 1.6 × 10 −8 Hz T n 125 MeV g * (T n ) 20

. 92 ……T
Case II: domain-wall dominanceT caus < T ann < T dom < T j QCD Case I: radiation dominance T BBN < T dom < T ann < T j QCD 0BBN , T dom < T j QCD < T ann ……

Figure 2 :
Figure 2: Different situations for different values of σ with A = 0.8 and V bias = (100 MeV) 4 .The PTA GW results prefer a value of σ that is in the radiation dominance case with the annihilation temperature comparable to the QCD phase transition temperature (see Sec. 5.4 for details).Here, V max bias is assumed to be valid up to T ann = 148 MeV.

Figure 8 :
Figure 8: Thermal evolution of various condensates with θ = 0 (left) and θ = π (right).In the left panel one has a crossover phase transition at T ≈ 146 MeV, while in the right panel one has a first-order phase transition at T ≈ 129 MeV.