Heavy QCD axion model in light of pulsar timing arrays

Recently, pulsar timing array experiments reported the observation of a stochastic gravitational wave background in the nanohertz range frequency band. We show that such a signal can be originated from a cosmological first-order phase transition (PT) within a well-motivated heavy (visible) QCD axion model. Considering the Peccei-Quinn symmetry breaking at the TeV scale in the scenario, we find a supercooled PT, in the parameter space of the model, prolonging the PT with the reheating temperature at the GeV scale.


I. INTRODUCTION
Gravitational wave (GW) experiments have provided new paths to explore the Universe and examine models of physics.Recently, pulsar timing array (PTA) experiments, including the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [1], and the European Pulsar Timing Array (EPTA) [2] along with the Parkes Pulsar Timing Array (PPTA) [3] and Chinese Pulsar Timing Array [4] released their latest results observing significant evidence of a GW following the Hellings-Downs pattern [5] in the angular crosscorrelation of pulsar timing residuals, supporting a stochastic gravitational background (SGWB).
An intriguing candidate that we propose in this paper is a supercooling phase transition (PT) which ends around the GeV scale; such a PT can naturally be accommodated in the context of heavy QCD axion models [50,51].A cosmic FOPT may be accomplished due to a spontaneous broken symmetry and the process is accompanied by the nucleation of bubbles separating true and false vacua.In a supercooled PT, the vacuum energy is dominated and the Universe remains in the false vacuum for a long period and supercools, increasing the PT duration [21,[52][53][54][55][56].
In this work, we consider a scenario within heavy axion models which are very appealing in that not only can the strong CP problem [57] be addressed but also the models have richer phenomenology [58,59] relative to the invisible QCD axion models [60][61][62][63].Moreover, the relation of the axion mass, m a , and its decay constant, f a is modified so that larger m a and lower f a are allowed.These models are also motivated in connection with the Peccei-Quinn (PQ) quality problem [64][65][66] (see Appendix A).
We explore a supercooled PT with a U(1) PQ symmetry breaking at the TeV scale.Based on the nearly conformal dynamics of the PQ scalar field, we analytically find important parameters encoding the GW spectrum and show that the corresponding signals can be well fitted to the PTA data.

II. PQ PHASE TRANSITION IN THE MODEL
The CP-violating θ parameter in the strong interactions is experimentally bounded θ ≲ 10 −10 [67] and the nature of this smallness, known as the strong CP problem, is a big puzzle in particle physics.One of the interesting solutions to this problem is based on a U(1) global symmetry, U(1) PQ , first proposed by Peccei and Quinn [57].At energies higher than the electroweak (EW) scale, the symmetry is spontaneously broken and the axion as the pseudo-Nambu-Goldstone boson is generated.The interaction of the axion with gluons at the QCD scale, (a/ f a + θ )G G, and as a result the vacuum expectation value of the axion can cancel the θ term [68].Considering the astrophysical constraints on the axion decay constant, 10 8 GeV ≲ f a ≲ 10 17 GeV [69][70][71], there are two classes of models, known as invisible axion models [60][61][62][63].However, as mentioned before, f a can be lowered, f a ∼ (1 − 100) TeV in the so-called visible axion models, still addressing the strong CP problem.In this sense, there are diverse scenarios, e.g.enlarging the QCD color gauge group SU(3 + N ) [72] or assuming some hidden sector with a confining scale larger than the QCD cutoff Λ QCD such that additional terms are sup-plemented to the axion potential and change the axion mass [58,73].In the case of adding a hidden sector, one can consider a whole copy of the standard model (SM), with the SM × SM ′ gauge group, and particles of the two sectors can interact gravitationally or by some very feebly couplings.Furthermore, a Z 2 mirror symmetry between particles of the two sectors can be imposed and soft terms breaking the symmetry can induce Λ ′ QCD ≫ Λ QCD [58,73].Here, we consider a DFSZ axion model [60,61], containing a gauge-singlet PQ scalar field, Φ, and two Higgs doublets under SU(2) L , H u and H d , which is supplemented with its hidden copy.Indeed, we consider a theory with the SM × SM ′ gauge group and each sector interacting with the PQ scalar (see Appendix B) and has the same PQ symmetry. 1Therefore, the axion couples to both QCD sectors and its mass would be m a ∼ Λ ′2 QCD / f a [58,73] with f a ∼ TeV.In this work, we do not go through axion interactions and focus on the UV theory and study the PQ PT associated with the PQ symmetry breaking.We explore a supercooled PQ PT along the direction of PQ scalar dynamic, ⟨Φ⟩ = v φ / √ 2. We consider the case where mass parameters are small, µ 2 d , µ 2 u , λ φ f 2 a ≪ f 2 a (see Appendix B) with the similar condition in the hidden sector, and main contributions to the potential arises from the one-loop Coleman-Weinberg (CW) quantum correction [74] and the leading contribution from thermal corrections [75]; thereby masses and the symmetry breaking scale are generated radiatively [76].
At the zero-temperature limit, quantum corrections (see Appendix C) contribute to the potential that is approximately scale invariant as where For the sake of simplicity, we only considered κ u and κ ′ u couplings which are determined at the vacuum and are then evolved to a scale φ through the beta function (see Appendix C).Thus, at the nonzero vacuum, 1 An analogous procedure can also be applied to a KSVZ axion model [62,63].
Considering thermal corrections (see Appendix C) and high temperatures, the origin point would be the minimal, and the potential would be where As temperature goes down, negative values of λ can change the potential slope and induce a bump so that at the critical temperature, T c , two degenerate states, ∆V (T ) = V (0, T ) − V (v φ (T ), T ) = 0, indicating a FOPT, can be generated due to these values of λ , Fig. (1).At temperatures below T c , the vacuum v φ (T ) is the favorable one and at some temperature T n bubbles of the new phase are nucleated.In fact, bubble nucleation occurs when the bubble formation probability per unit Hubble space-time volume, which is proportional to exp (−S 3 (T )/T ) where S 3 (T ) is the bounce action quantifying the tunneling process, would be of the order of 1.As a result, we can find the nucleation temperature T n via the relation [77, 78] where the Hubble parameter H is given by where the reduced Planck mass is M pl ≃ 2.43×10 18 GeV, ξ = π 2 g * /30, and g * ≃ 107 is the effective number of relativistic degrees of freedom, assuming that it does not change during the PT.At low temperatures, T ≪ T c , it is expected v φ (T ) = f a , and ∆V (T ) is well approximated by ∆V .Furthermore, in this limit, the bounce action can be well approximated by [79] and hence we can find the nucleation temperature as 4 ln At the time when the transition completes, from the conservation of energy, one obtains ρ R (T * ) = ρ R (T p ) + ∆V [52] where ρ R (T ) = ξ T 4 is the radiation energy and T p is the percolation temperature.Assuming T p ≃ T n , the reheating temperature can be computed by For the case of fast reheating H * ≡ H(T * ) ≃ H(T n ).Another important quantity in characterizing the generated GWs is the inverse of PT duration calculated by the following relation where From Eq. ( 6) and V (φ 0 , T ) = 0, for T ≪ T c , and Expecting β /H * ∼ O(1 − 10) for the case of supercooling, we fix the parameter γ with the data.The strength of the supercooled PT is also obtained by In the next section, calculating the aforementioned quantities, we obtain the GW energy density spectrum.

III. SGWB SIGNAL
In this section, based on the obtained PT quantities, we calculate the GW spectrum of the supercooled PT.In this sense, sources contributing to the GW production are bubble wall collisions and fluid motions.Relying on the recent numerical computation [80], the present-day GW energy density spectrum is obtained by where the spectral shape of the GW is obtained by and κ is the fraction of the vacuum energy converted to the kinetic energy of bubble walls, κ vac ≃ 1/(1 + 5/(β R eq )), and fluid motions κ f = 1 − κ vac [80].As for the spectral shape, Eq. ( 17), the spectrum's low-and high-frequency power-law tails are controlled by a and b, its transition width from one power law to another is determined by c, and its amplitude depends on Ā.
In our case with α ≫ 1, bubbles collide in the vacuum and the bubble wall velocity can reach the speed of light; thus we consider κ vac ≃ 1 corresponding to β R eq ≫ 1.Therefore, according to Table I of Ref. [80], corresponding fit parameters of the spectral function appropriate for a global broken symmetry, T rr ∝ R −2 (where T rr is the maximum of the stress-energy tensor and R is the bubble radius), would be as Ā ≃ 5.The present redshifted peak frequency is given by where for the case of T rr ∝ R −2 , f p /β ≃ 0.64/(2π) and for the envelope approximation f p /β ≃ 1.33/(2π) [80].
As can be seen from Fig. 2, using datasets of NANOGrav 15 yr, EPTA release 2, and PPTA data release 3 (DR3), for f a = 1 TeV, and some representative values of parameters, κu = κ′ u = 0.0001, κ u = κ ′ u = 0.001 and λ (T ) = −0.0015,which correspond to T n = 100 MeV, T * = 1 GeV, α = 8995, and β /H * = 17, 3.4, 1.7 with the corresponding parameters γ = 1, 5, 10, respectively, the GW signals can be very consistent with the PTA data from the bubble collision spectra for both the spectral shape of T rr ∝ R −2 (solid line) and for the envelope (dashed line).Remarkably, we see that the highfrequency range of the GW signals with the envelope approximation falls within the sensitivity range of future space-based GW experiments such as Laser Interferometer Space Antenna (LISA), Deci-hertz Interferometer Gravitational Wave Observatory (DECIGO), and Big Bang Observer (BBO) and can be probed by these detectors.
Furthermore, using the PTA data, we perform a best-fit analysis over a range of γ parameter (β /H * ) through the χ 2 test based on the relation where Ω th h 2 denotes the GW predicted by the model, Ω exp h 2 represents the observed GW signal by PTA experiments, and Si is the deviation from the midpoint value of each data point in log 10 Ω exp h 2 within the uncertainty range.As shown in Figs. 3, 4 the best-fit point of the GW energy spectrum taking into account all datasets derived from the T rr ∝ R −2 case for the aggregated dataset is γ = 3.28 (β /H * = 5.1) and for the envelope case, the best-fit point would be γ = 13.3 (β /H * = 1.3).DECIGO [82], BBO [83], Einstein Telescope (ET) [84] and Cosmic Explorer (CE) [85].we the find best-fit point in terms of γ parameter, which is 2.52 for EPTA (solid purple), 2.17 for NANOGrav (dashed green), and 1.81 for PPTA (dotted orange).The best-fit value for the aggregated dataset is given by 3.28 (solid blue) corresponding to β /H * = 5.1.
The insert shows GW spectra for the best-fit values of γ in addition to observed signals.

IV. DISCUSSION
Another PQ PT consequence of the model is that the symmetry breaking may induce a similar or different EW symmetry breaking scale in the ordinary and hidden sectors.For instance, as for v ′ EW , this can be realized by these terms of the Lagrangian Based on the GWs with the envelope approximation.wefind the best-fit point in terms of γ parameter, which is 13.55 for EPTA (solid purple), 9.49 for NANOGrav (dashed green), and 7.93 for PPTA (dotted orange).The best-fit value for the aggregated dataset is given by 13.30 (solid blue) corresponding to β /H * = 1.3.The insert shows GW spectra for the best-fit values of γ in addition to observed signals.[73].Consequently, in this case, we obtain the axion mass m a ≳ 1 MeV, while considering v EW ≲ v ′ EW < f a and assuming Λ ′ QCD ≲ v EW , one finds m a ≲ 10 GeV.Heavy axions may interact with visible and hidden (dark) photons so that the axion decay rate would be where ᾱ = 1/137 is the fine structure constant; z = m u /m d ; and E and N are the electromagnetic and QCD anomaly coefficients, respectively, [51].The heavy QCD axions can also couple to SM fermions Decaying the heavy axions to these photons contributes to the effective number of relativistic species which is strongly constrained by CMB measurements, N eff = 2.96 +0.34 −0.33 [86], putting the most stringent constraints on the heavy QCD axion mass.According to the analysis of Ref. [87], for E/N = 8/3 which may be the case of a DFSZ model [88], m a < 100 MeV would be excluded.However, lower masses, 1 MeV < m a < 100 MeV where the dark photon contribution to N eff compensates neutrino dilution from axion decays to photons, are allowed for E/N = 1/3.
We assume a baryon asymmetric dark sector.Thus, considering v EW ≲ v ′ EW < f a , the mass range and relic abundance of stable particles in the hidden sector are approximately the same as those of the ordinary sector.As a result, the relic abundance of these particles cannot give rise to cosmological problems.Moreover, we assume that in the hidden sector there is no seesaw mechanism, which allows the Dirac neutrino mass in this sector as [73] where M R is the right-handed neutrino mass in the ordinary sector.For M R ≳ 10 7 GeV and m a ≳ 1 MeV (Λ ′ QCD ≳ 1 GeV), the hidden neutrinos decay to the hidden electron and pion, v ′ → e ′ + π ′ , and thereby these neutrinos do not contribute to N e f f .
Concerning astrophysical bounds, the model parameters can be constrained from SN 1987A [89], with the supernova core cooling rate due to axion emission.According to Eq. ( 21), one can obtain the axion mean free path as ) where E a stands for typical energies of the axions.Hence, for E a ≳ 200 MeV, f a = 1 TeV and m a ≳ 100 MeV, we find that λ m f p is much smaller than the approximate radius of the supernova core R ∼ 10 km.Thus, because of the scattering and reabsorption, the axions are trapped within the supernova, and cannot take away any energy [90,91].
In addition, the emission of dark photons can lead to a very fast core cooling rate of the supernova and the associated energy loss rate, which is bounded L γ ′ < 10 52 erg/s, and can be obtained as L γ ′ ≃ Γ aγ ′ γ ′ m a (4πR 3 /3)(1.2T 3/π2 ) [58].Thus, taking into account a core temperature T ≃ 10 MeV, f a = 1 TeV, m a = 1 MeV, and z ′ = z ≃ 0.5, the bound can be satisfied, while for m a ≳ 100 MeV, it may require z ′ ≲ 10 −4 .
Considering the axion-electron interactions, the astrophysical bounds on helium burning lifetime of the horizontal branch stars imply g ae ∼ m e / f a < (2.5 − 6) × 10 −13 [51,58], and hence restrict f a to its invisible values.However, this bound does not apply to the heavy axions with m a > 300 keV.In fact, in the stellar cores whose typical temperatures reach 10 keV, the axion production rate is suppressed by the exponential factor exp(−m a /T ).Therefore, astrophysical bounds allow heavy axions with masses above 100 MeV for f a = 1 TeV.
The viable parameter space of the model can also be probed via other experiments such as proton beam dump [92,93], kaon decays [94,95], and colliders [96].For m a ≳ 3m π , the heavy axions may be produced via B → K a, where B and K denote B mesons and kaons, respectively, and then mainly decay into three pions.For m a ≲ 3m π , the dominant decay channels are the axion decay into two photons, electrons, and muons, for m a > 2m e and m a > 2m µ , respectively.In this sense, m a ≲ 400 MeV for f a at the TeV scale is excluded by experiments such as kaon decays and the proton beam dump, while for m a ≳ 400 MeV the search is challenging due to the dominance of the hadronic decay mode [97].
It is noteworthy that the heavy axion cannot be a dark matter candidate in this context.However, some stable particles in the hidden sector may contribute to the dark matter relic abundance. 2.We leave detailed calculations in connection with axion interactions and possible dark matter candidates in the presented scenario to a future work.

V. CONCLUSION
With regard to the recently detected GW signals by PTA experiments, we have proposed a heavy QCD axion scenario, with a PQ symmetry breaking scale around the TeV scale, based on a (DFSZ) axion model which is supplemented with its hidden copy.We investigated a supercooled PQ FOPT which is derived by CW quantum corrections and analytically found important PT quantities, including the reheating temperature, the inverse duration of the PT, and the strength of the PT.We have shown that within the parameter space of the model the generated GWs from the bubble wall collisions can be consistently fitted with the recently observed PTA data.Furthermore, it is shown that the high-frequency range of such GW signals can be probed by future space-based GW detectors.
In addition, Heavy axion models are well-motivated scenarios allowing the (m a , f a ) relation to be relaxed, still addressing the strong CP problem.It is shown that considering the EW symmetry breaking scale in the ordinary and hidden sectors induced after the PQ symmetry breaking, a range of axion masses, 1 MeV ≲ m a ≲ 10 GeV can be obtained.Nevertheless, viable regions of the mass space are limited via CMB observation due to the dark photon contribution to N eff , and also other experimental setups such as the proton beam dump.
As a result, interpreting the recent GW data, the model provides a setup which can be probed by future GW detectors, CMB telescopes and collider experiments.Invisible axion models with f a ≳ 10 8 GeV suffer from the axion quality problem.Indeed, it is expected that the PQ symmetry is explicitly broken by higher-dimensional Plancksuppressed operators which can spoil the axion solution.In our model, the PQ symmetry may be explicitly broken by the following five-dimensional operator

FIG. 1 .
FIG.1.For three different values of temperature as labeled and representative values of A = −0.0053,B = 0038, and D = 0.09, the radiatively induced potential is illustrated.Above T c , the origin is the only minimum.At T c , two degenerate states appear and at lower temperatures the nonzero minimum is favored.The Universe remains in the false vacuum until the tunneling process is well developed at low temperatures due to λ (T ).

FIG. 3 .
FIG.3.Based on the GWs with the spectral shape of T rr ∝ R −2 , we the find best-fit point in terms of γ parameter, which is 2.52 for EPTA (solid purple), 2.17 for NANOGrav (dashed green), and 1.81 for PPTA (dotted orange).The best-fit value for the aggregated dataset is given by 3.28 (solid blue) corresponding to β /H * = 5.1.The insert shows GW spectra for the best-fit values of γ in addition to observed signals.
FIG.4.Based on the GWs with the envelope approximation.wefind the best-fit point in terms of γ parameter, which is 13.55 for EPTA (solid purple), 9.49 for NANOGrav (dashed green), and 7.93 for PPTA (dotted orange).The best-fit value for the aggregated dataset is given by 13.30 (solid blue) corresponding to β /H * = 1.3.The insert shows GW spectra for the best-fit values of γ in addition to observed signals.

ACKNOWLEDGMENTS
SS thanksFazlollah Hajkarim and Seyed Mohammad Mahdi Sanagostar for helpful discussions about data analy-sis.This work is supported in part by the National Key Research and Development Program of China under Grant No. 2021YFC2203004, and the National Natural Science Foundation of China (NSFC) under Grants No. 12075041 and No. 12147102.APPENDIX A: THE PQ QUALITY PROBLEM

2 * /m 2 a[ 51 ]
where M P denotes the Planck mass and φ a is the radial component of the PQ scalar.Such an explicit PQ symmetry breaking term induces a mass term in the axion field, m 2 * ∼ 10 −2 λ a ( f 3 a /M P ), and moves the axion potential minimum away from the CP conserving minimum, ∆ θ ∼ m .Therefore, assuming λ a = O(1), for f a = 1 TeV and m a ≳ 100 MeV, the model satisfies the current bound ∆ θ ≲ 10 −10 .