Observational prospects for gravitational waves from hidden or dark chiral phase transitions

TL;DR

This paper investigates gravitational wave signals from first-order chiral phase transitions in hidden sectors, finding that upcoming detectors like LISA, DECIGO, and BBO could observe these signals for transitions at specific energy scales, with transition dynamics being more rapid than previously assumed.

Contribution

The study provides the first detailed calculations of gravitational wave spectra from hidden chiral phase transitions using effective models, revealing significantly faster transition dynamics than prior estimates.

Findings

Transition inverse duration $eta/H$ is at least 10,000, much larger than previous assumptions.

LISA can potentially detect signals from ~100 MeV transitions.

DECIGO and BBO could observe backgrounds from 1 GeV to 10 TeV transitions.

Abstract

We study the gravitational wave (GW) signature of first-order chiral phase transitions ($\chi$PT) in strongly interacting hidden or dark sectors. We do so using several effective models in order to reliably capture the relevant non-perturbative dynamics. This approach allows us to explicitly calculate key quantities characterizing the $\chi$PT without having to resort to rough estimates. Most importantly, we find that the transition's inverse duration $\beta$ normalized to the Hubble parameter $H$ is at least two orders of magnitude larger than typically assumed in comparable scenarios, namely $\beta/H\gtrsim\mathcal{O}(10^4)$. The obtained GW spectra then suggest that signals from hidden $\chi$PTs occurring at around 100 MeV can be in reach of LISA, while DECIGO and BBO may detect a stochastic GW background associated with transitions between roughly 1 GeV and 10 TeV. Signatures of transitions at higher temperatures are found to be outside the range of any currently proposed experiment. Even though predictions from different effective models are qualitatively similar, we find that they may vary considerably from a quantitative point of view, which highlights the need for true first-principle calculations such as lattice simulations.

Tables

Table 1. Table 1: Benchmark points used throughout the paper. The model parameters which correspond to the shown mass spectra and pion decay constant are compiled in Appendix   D .

benchmark point	$f_{\pi }$ [MeV]	$m_{\pi }$ [MeV]	$m_{\sigma }$ [MeV]	$m_{{\eta }^{\prime }}$ [MeV]	$m_a$ [MeV]
A	$72.3$	$0$	$248$	$458$	$491$
B	$90$	$0$	$400$	$671.5$	$696.8$
C	$73.6$	$0$	$291$	$328.353$	$430.8$
D	$107.5$	$0$	$693.5$	$535.251$	$791.7$

Table 2. Table 2: Parameters characterizing the (hidden) chiral phase transition as predicted by the considered effective models for the benchmark points from Table   1 . All calculations were performed assuming g ⋆ = 47.5 subscript 𝑔 ⋆ 47.5 g_{\star}=$47.5$ which corresponds to n f = 3 subscript 𝑛 𝑓 3 n_{f}=3 light hidden fermions and n c 2 − 1 = 8 superscript subscript 𝑛 𝑐 2 1 8 n_{c}^{2}-1=8 hidden gluons. The contents of the last two columns were obtained by least-square fits of the function b ​ ( 1 − T / T c ) − γ 𝑏 superscript 1 𝑇 subscript 𝑇 𝑐 𝛾 b(1-T/T_{c})^{-\gamma} to the results of our explicit calculations of S 3 / T subscript 𝑆 3 𝑇 S_{3}/T .

benchmark point	effective model	$T_c$ [MeV]	$T_n$ [MeV]	$g_{\star }\alpha$	$\beta /H$ [${10}^4$]	$\gamma$	$b$
A	NJL	$71.71$	$70.52$	$\mathrm{3.396\hspace{0.17em}25}$	$1.79$	$1.76$	$1.25\text{⋅}{10}^{-1}$
	PNJL	$121.76$	$121.36$	$1.083$	$9.43$	$1.82$	$5.04\text{⋅}{10}^{-3}$
	LSM	$101.76$	$101.02$	$\mathrm{0.769\hspace{0.17em}405}$	$4.41$	$1.8641$	$1.8480\text{⋅}{10}^{-2}$
B	NJL	$107.14$	$106.38$	$2.565$	$4.27$	$1.80$	$2.33\text{⋅}{10}^{-2}$
	PNJL	$140.49$	$140.17$	$\mathrm{2.018\hspace{0.17em}75}$	$13.75$	$1.87$	$2.04\text{⋅}{10}^{-3}$
	LSM	$145.83$	$145.27$	$\mathrm{0.692\hspace{0.17em}788}$	$8.59$	$1.8944$	$4.5587\text{⋅}{10}^{-03}$
C	NJL	$90.84$	$90.58$	$\mathrm{1.220\hspace{0.17em}75}$	$11.09$	$1.81$	$3.95\text{⋅}{10}^{-3}$
	PNJL	$131.33$	$131.12$	$0.9025$	$45.66$	$1.85$	$2.40\text{⋅}{10}^{-4}$
	LSM	$100.52$	$99.94$	$\mathrm{1.081\hspace{0.17em}91}$	$5.6649$	$1.8705$	$1.1314\text{⋅}{10}^{-02}$
D	NJL	$180.31$	$180.27$	$\mathrm{0.403\hspace{0.17em}75}$	$162.59$	$1.92$	$1.39\text{⋅}{10}^{-5}$
	PNJL	$198.27$	$198.25$	$0.3325$	$244.94$	$1.86$	$9.7\text{⋅}{10}^{-6}$
	LSM	$175.25$	$174.51$	$\mathrm{1.211\hspace{0.17em}77}$	$7.8341$	$1.9089$	$5.0214\text{⋅}{10}^{-03}$

Table 3. Table 3: Model parameters corresponding to the mass spectra listed in Table   1 . The UV cutoff scale of the (P)NJL model is fixed to 𝚲 = 𝟗𝟑𝟎   𝐌𝐞𝐕 𝚲 times 930 𝐌𝐞𝐕 \Lambda=$930\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ for all benchmark points.

benchmark point	(P)NJL model		Linear Sigma Model
	$𝑮$ [GeV-2]	${𝑮}_{𝑫}$ [GeV-5]	${𝒎}^{\mathrm{𝟐}}$ [${\mathrm{𝐌𝐞𝐕}}^{\mathrm{𝟐}}$]	${𝝀}_{𝝈}$	${𝝀}_{𝒂}$	$𝒄$ [MeV]
A	$3.84$	$-90.65$	$-\mathrm{𝟒𝟐𝟎𝟗}$	$16.76$	$12.91$	$\mathrm{𝟐𝟑𝟔𝟗}$
B	$3.988$	$-106.445$	$\mathrm{𝟒𝟖𝟒𝟖}$	$25.54$	$15.22$	$\mathrm{𝟒𝟎𝟗𝟏}$
C	$4.0$	$-60.0$	$\mathrm{24\hspace{0.17em}371.2}$	$14.84$	$13.99$	$1196.075$
D	$5.0$	$-60.0$	$\mathrm{192\hspace{0.17em}722.186}$	$33.25$	$25.14$	$2176.0$

Equations

$$\frac{\beta }{H}\approx 4\log \frac{M_{\text{Pl}}}{T_n}\approx \mathcal{O}(100),$$

$$\displaystyle\beta/H\simeq\mathcal{O}(${10}^{4}$)\;{.}$$

$${\mathcal{L}}_m=\sum _{i=1}^{n_f}{\bar{q}}_i(\mathrm{i}/D-m_i)q_i=\sum _{i=1}^{n_f}\left[{\bar{q}}_{Li}\mathrm{i}/D{q}_{Li}+{\bar{q}}_{Ri}\mathrm{i}/D{q}_{Ri}-m_i({\bar{q}}_{Li}{q}_{Ri}+\text{h.c.})\right],$$

$${\mathcal{G}}^{\prime }$$

$$\mathcal{G}$$

$${\mathcal{L}}_{\text{NJL}}=\mathrm{Tr}\bar{q}\mathrm{i}/\partial q+2G\mathrm{Tr}({\Psi }^{†}\Psi )+G_D(\det \Psi +\text{h.c.})\text{with}{\Psi }_{ij}={\bar{q}}_j(1-{\gamma }_5)q_i,$$

$$q$$

$$q_L$$

$$\Psi$$

$$-4G⟨\Psi ⟩=(\sigma +\mathrm{i}{\eta }^{\prime })\text{mathds}1+2(a_a+\mathrm{i}{\pi }_a)T^a.$$

$$[T^a,T^b]=\mathrm{i}{f}_{abc}{T}^c\text{and}\{T^a,T^b\}=d_{abc}{T}^c$$

$$V_{\text{eff}}^{\text{NJL}}(\bar{\sigma },T)=V_0^{\text{NJL}}(\bar{\sigma })+V_{\text{CW}}^{\text{NJL}}(\bar{\sigma })+V_{\text{FT}}^{\text{NJL}}(\bar{\sigma },T),$$

$$V_0^{\text{NJL}}(\bar{\sigma })$$

$$V_{\text{CW}}^{\text{NJL}}(\bar{\sigma })$$

$$V_{\text{FT}}^{\text{NJL}}(\bar{\sigma },T)$$

$$M$$

$$Z_{\sigma }^{-1}(\sigma )={\frac{\text{mbox}d{\Gamma }_{\sigma \sigma }(p^2,\sigma )}{\text{mbox}d{p}^2}|}_{p^2=0},$$

$$f_{\pi }^2=4n_c{M}_c^2{I}_0(v_{\sigma })$$

$${\mathcal{L}}_{\text{PNJL}}^{\text{MFA}}={\mathcal{L}}_{\text{NJL}}^{\text{MFA}}-V_{\text{glue}}(L,T),$$

$$T^{-4}V_{\text{glue}}(L,T)=-\tfrac{1}{2}a(T)L\bar{L}+b(T)\log\bigl{[}1-6L\bar{L}-3(L\bar{L})^{2}+4(L^{3}+\bar{L}^{3})\bigr{]}\;{,}$$

$$a(T)=a_0+a_1\frac{T_{\text{glue}}}{T}+a_2{\left(\frac{T_{\text{glue}}}{T}\right)}^2\text{and}b(T)=b_3{\left(\frac{T_{\text{glue}}}{T}\right)}^3.$$

$$a_0=3.51,a_1=-2.47,a_2=15.2,b_3=-1.75.$$

$$V_{\text{eff}}^{\text{PNJL}}(\bar{\sigma },L,T)=V_0^{\text{PNJL}}(\bar{\sigma })+V_{\text{CW}}^{\text{PNJL}}(\bar{\sigma })+V_{\text{FT}}^{\text{PNJL}}(\bar{\sigma },L,T)$$

$$V_0^{\text{PNJL}}(\bar{\sigma })=$$

$$V_{\text{CW}}^{\text{PNJL}}(\bar{\sigma })=$$

$$V_{\text{FT}}^{\text{PNJL}}(\bar{\sigma },L,T)=$$

$$\displaystyle+T^{4}\Bigl{(}-\tfrac{1}{2}a(T)L^{2}+b(T)\log\bigl{[}1-6L^{2}-3L^{4}+8L^{3}\bigr{]}\Bigr{)}\;{,}$$

$${\mathcal{L}}_{\text{LSM}}=\mathrm{Tr}{\partial }_{\mu }{\Phi }^{†}{\partial }^{\mu }\Phi -V_{\text{tree}}^{\text{LSM}}(\Phi ),$$

$$\displaystyle\begin{split}V_{\text{tree}}^{\text{\tiny{LSM}}}(\Phi)={}&-m^{2}\operatorname{Tr}(\Phi^{\dagger}\Phi)+\tfrac{1}{2}(\lambda_{\sigma}-\lambda_{a})\bigl{(}\operatorname{Tr}(\Phi^{\dagger}\Phi)\bigr{)}^{2}+\tfrac{3}{2}\lambda_{a}\operatorname{Tr}\bigl{(}(\Phi^{\dagger}\Phi)^{2}\bigr{)}\\ &+\sqrt{\tfrac{2}{3}}\,c\bigl{(}\det\Phi+\det\Phi^{\dagger}\bigr{)}\;{.}\end{split}$$

$$\Phi =\frac{1}{\sqrt{6}}(\sigma +\mathrm{i}{\eta }^{\prime })\text{mathds}1+(a_a+\mathrm{i}{\pi }_a)T^a$$

$$\Phi \to U_L\Phi U_R^{†}\text{with}{U}_{L,R}\in \text{U(}3{\text{)}}_{L,R},$$

$$⟨\Phi ⟩=\frac{v_{\sigma }}{\sqrt{6}}\text{mathds}1\text{with}{v}_{\sigma }=\sqrt{\frac{3}{2}}{f}_{\pi },$$

$$V_{\text{eff}}^{\text{LSM}}(\bar{\sigma },T)=V_0^{\text{LSM}}(\bar{\sigma })+V_{\text{FT}}^{\text{LSM}}(\bar{\sigma },T),$$

Full text

1. I Introduction

2. II Low-energy effective models of QCD-like theories

3. II.1 The Nambu–Jona-Lasinio (NJL) model

4. II.2 The Polyakov loop enhanced NJL model

5. II.3 The linear sigma model

6. III Chiral phase transition

7. IV Stochastic gravitational wave signal

8. V Results

9. V.1 Hidden chiral phase transitions at $\mathcal{O}(100\,\mathrm{MeV})$

10. V.2 Hidden chiral phase transitions at higher temperatures

11. VI Summary and conclusions

12. A Collection of integrals

13. B Details of the (P)NJL model

14. C Basics of the CJT formalism

15. D Model parameters

I Introduction

Since their first direct detection on Earth [1, 2, 3], gravitational waves (GWs) have provided a new and complementary approach towards observational astrophysics. Not only do they allow to study violent transient phenomena like black hole or neutron star mergers, but they also offer a unique way to probe the physics of the early Universe. In particular, it is well known that first-order cosmic phase transitions (PTs) can create a stochastic background of gravitational waves [4, 5, 6, 7, 8, 9]. Even though the Standard Model (SM) of particle physics does not feature such a PT [10, 11, 12, 13, 14], they are predicted by a plethora of beyond-the-SM (BSM) theories. For instance, successful electroweak baryogenesis typically requires a strong first-order electroweak phase transition (EWPT), which can e.g. be obtained in models with an augmented scalar particle content [15]. Of course, gravitational waves may also be produced during phase transitions in dark or – more generally – hidden sectors. Correspondingly, studies investigating gravitational wave signatures of a wide range of new physics scenarios have received a lot of attention in the recent years, see e.g. [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36].

In the present article, we focus on BSM theories with a new confining gauge force $\text{SU(}n_{c}\text{)}$ acting on $n_{f}$ flavors of hidden fermions in the fundamental representation. Similar to the situation in quantum chromodynamics (QCD), strong (i.e. non-perturbative) dynamics are anticipated to spontaneously break the hidden counterpart of chiral symmetry in such models. Unlike in real-world QCD,111However, the crossover transition predicted for real-world QCD can influence the spectrum of inflationary gravitational waves [37, 38, 39, 40, 41]. the corresponding chiral phase transition can, however, be strong and of first order. It is the gravitational wave background associated with such a transition that we will be interested in, cf. [16].

Hidden sectors of the above described type are particularly attractive for several reasons. They are, for instance, featured in various models of dark matter (DM) like strongly interacting massive particle (SIMP) DM [42, 43], or composite DM as discussed in e.g. Refs. [44, 45, 46, 47, 48, 49, 50]. Similarly, they are introduced in Twin Higgs models [51, 52], Hidden Valley theories [53] or vector-like confinement models [54]. Lastly, strongly coupled hidden sectors commonly occur in theories based on classical scale invariance, since they offer one of the two222The other one being the (perturbative) Coleman-Weinberg mechanism [55]. known possibilities to dynamically generate a scale in a previously scaleless setup. Several realistic classically conformal SM extensions along those lines were suggested e.g. in Refs. [56, 57, 58, 59, 60, 61]. Reflecting this rich variety of applications, well-motivated scenarios exist at vastly different scales ranging from $100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ to above $10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, see Ref. [16] for an overview.

Given the presented appealing physics case for BSM theories with new confining gauge forces, it is desirable to quantify gravitational wave signatures from a potential hidden chiral phase transition, in order to assess the prospects for their observation. Due to the problem’s non-perturbative nature, the corresponding analyses generally prove difficult. Accordingly, no true first-principle calculations exist yet for the setup under consideration.333The bubble nucleation rate in a first-order EWPT was computed using first-principle lattice methods in Ref. [62]. Similar calculations were performed in purely scalar theories [63]. Based on some of the aforementioned lattice results and employing dimensionally reduced effective field theory, the authors of Ref. [64] recently presented a non-perturbative analysis of the GW spectrum associated with a first-order EWPT. There are, however, various ways to proceed. For a start, one may simply estimate the relevant quantities based on certain assumptions and/or general arguments. For instance, the inverse duration of the phase transition $\beta$, which is a key parameter in understanding the associated GW spectrum, is commonly estimated by [5, 6]

[TABLE]

where $H$ is the Hubble parameter at the time of the transition, $M_{\text{\tiny{Pl}}}$ is the (reduced) Planck mass and the second relation holds for a wide range of transition temperatures $T_{n}$. Analyses of GW signals based on the estimate of Eq. 1 can e.g. be found in Refs. [65, 66, 16].

Here, we will follow a different approach. Namely, we will model the dynamics in the strongly coupled hidden sector by employing various effective theories, which are known to work well in capturing non-perturbative aspects of the low-energy regime of real-world QCD. Thus we will be able to explicitly calculate relevant quantities characterizing the hidden chiral phase transition such as the bubble nucleation rate, the transition temperature or $\beta/H$. Interestingly, we find that all our calculations result in values for $\beta/H$ which are much larger as compared to the rough estimate of Eq. 1, namely

[TABLE]

Since the energy density of the GW background decreases – depending on the nature of its source – linearly or even quadratically with $\beta/H$, and since the spectrum’s peak frequency increases linearly with $\beta/H$ (see e.g. [67]), conclusions about the observability of the GW signal change considerably, as we will demonstrate in the main text.

The article is structured as follows. In Section II, we briefly introduce the low-energy effective models that we employ throughout the paper. Those are the Nambu–Jona-Lasinio (NJL) model [68, 69] and its Polyakov-loop enhanced variant (PNJL) [70], as well as the linear sigma model (LSM) [71]. For definiteness, we will consistently focus on the case of QCD-like hidden sectors with $n_{c}=3$ colors and $n_{f}=3$ hidden fermion flavors, which we assume to be massless for simplicity. In Sections III and IV we review important aspects of first-order phase transitions and of the associated GW background spectra, respectively. Section V then contains this work’s main results. Namely, we use the aforementioned effective models to explicitly calculate different parameters that are crucial in characterizing the hidden PT. Furthermore, we determine the corresponding GW signals and evaluate the prospects for their detectability with future observatories by calculating appropriate signal-to-noise ratios. We present our findings in such a way that they can be applied to hidden or dark sectors over a wide range of inherent scales, i.e. , transition temperatures. Finally, we conclude in Section VI.

II Low-energy effective models of QCD-like theories

In the present work we are interested in strongly interacting hidden sectors which we assume to be describable by a theory similar to QCD. In particular, the hidden analogue of the chiral phase transition ($\chi$PT) will be triggered by non-perturbative effects, requiring a treatment beyond standard perturbation theory. To that end, we here employ different effective models for the purpose of reliably capturing the phase transition’s dynamics. In order for such an approach to be successful, all of these models necessarily have to share the relevant symmetry properties of the underlying hidden sector, which we briefly outline in the following.

The matter part of a QCD-like theory with $n_{f}$ fermion flavors $q$ transforming in the fundamental representation of a hidden $\text{SU(}n_{c}\text{)}$ gauge group reads

[TABLE]

where $i$ is a flavor index. Color indices ranging from one to $n_{c}=3$ are suppressed for readability, but are to be understood as being summed over as well. In the chiral limit of vanishing fermion masses, $m_{i}=0$, the above Lagrangian classically exhibits a global chiral symmetry $\mathcal{G}^{\prime}:=\text{U(}n_{f}\text{)}_{L}\times\text{U(}n_{f}\text{)}_{R}$. Separating the transformations which treat the chiral components $q_{L}$ and $q_{R}$ alike (vector subgroup) from those which transform them differently (axial-vector subgroup), as well as using the Lie algebra isomorphism $\text{U(}N\text{)}\cong\text{SU(}N\text{)}\times\text{U(}1\text{)}$, the group $\mathcal{G}^{\prime}$ can be conveniently written as

[TABLE]

Nevertheless, the theory’s vacuum does not respect the full group $\mathcal{G}$, which is spontaneously broken down to $\text{SU(}n_{f}\text{)}_{V}\times\text{U(}1\text{)}_{V}$ by a finite expectation value of the fermion condensate, $\langle\bar{q}_{i}q_{i}\rangle\neq 0$. However, finite-temperature effects are anticipated to restore the chiral symmetry group $\mathcal{G}$ at sufficiently high temperatures [74, 75, 76]. The properties of the associated chiral phase transition depend on both the number $n_{f}$ of fermion flavors and their masses, as summarized in the famous Columbia plot [77]. Since we are interested in the gravitational wave signatures from a first-order hidden $\chi$PT, we will focus on the scenario of $n_{f}=3$ massless flavors in the following. This choice corresponds to the minimal number of flavors needed to obtain a first-order $\chi$PT [78, 77].

Apart from the correct symmetry structure, a successful low-energy effective model must, of course, also incorporate all degrees of freedom most relevant to the physics process of interest. In the case of the chiral phase transition in QCD-like theories, those will primarily be the hidden counterparts of the well-known mesons from real-world QCD, see e.g. [79]. In particular, the presence of a $\sigma$-like scalar flavor singlet is crucial since it carries the same quantum numbers as the aforementioned fermion condensate $\langle\bar{q}_{i}q_{i}\rangle$. Thus, by acquiring a finite vacuum expectation value (vev) $v_{\sigma}$ it triggers the spontaneous breakdown of chiral symmetry as indicated above. In the course of the present work we will additionally take into account the pseudoscalar octet $\vec{\pi}$, the pseudoscalar singlet ${\eta^{\prime}}$ and the scalar octet $\vec{a}$. Depending on the concrete model, further degrees of freedom may be taken into consideration.

The remainder of the present section briefly introduces the three effective theories considered throughout this article. In order to meaningfully compare predictions about gravitational waves as derived from the individual models, we choose their free parameters such that all models produce an identical meson spectrum, i.e. the same $\sigma$, ${\eta^{\prime}}$ and $\vec{a}$ masses, as well as the same pion decay constant $f_{\pi}$. Since we exclusively discuss the chiral limit, the pion octet is an exact Nambu-Goldstone boson and thus massless. The benchmark meson spectra used throughout the paper are listed in Table 1.

II.1 The Nambu–Jona-Lasinio (NJL) model

The NJL model was originally proposed to describe the nucleus [68, 69]. Nowadays, it is primarily used as a low-energy effective model of strongly coupled theories like QCD. The NJL model only describes the quarks, as its Lagrangian does not contain any of the mesons which are the physical degrees of freedom at low energies. As we will explicitly see the Lagrangian contains four- and six-fermion interactions which render the model perturbatively non-renormalizable. The theory therefore necessarily contains a UV cutoff $\Lambda$. Following Ref. [58], we employ a four-dimensional momentum cutoff scheme and fix $\Lambda=$0.93\text{,}\mathrm{G}\mathrm{e}\mathrm{V}$$ throughout the paper. An extensive review of the NJL model is given in, for example, Ref. [80].

The Lagrangian of the three-flavor NJL model in the chiral limit is given by [58]

[TABLE]

where the $q_{i}$ are Dirac fields corresponding to the (hidden) quarks, and $\Psi$ is a fermion bi-linear with flavor indices $i,j$ running from $1$ to $n_{f}=3$. Spinor and color indices are implicit. As we are only interested in the chiral limit, the above Lagrangian does not contain an explicit quark mass term. Importantly, the NJL model thus exhibits the global symmetries of massless QCD, see e.g. [81]. To be more precise, the first two terms are invariant under $\mathcal{G}^{\prime}=\text{U(}3\text{)}_{L}\times\text{U(}3\text{)}_{R}$, whereas the last term is the ’t Hooft determinant which mimics the anomalous $\text{U(}1\text{)}_{A}$ breaking mentioned below Eq. 4 [82, 73]. One can explicitly demonstrate that these symmetries are indeed present by realizing that the Dirac spinors $q$ transform under $\text{U(}1\text{)}_{V}$ and $\text{U(}1\text{)}_{A}$ as

[TABLE]

respectively, where $q=q_{L}+q_{R}$. The operators $U_{L}$ and $U_{R}$ correspond to unitary matrices with determinant equal to one. The above transformations clearly show the $\mathcal{G}^{\prime}$ invariance of the kinetic term. Establishing the symmetry properties of the other two terms is slightly more involved. As a start, note that the bi-linear $\Psi$ as defined in Eq. 6 can also be rewritten in terms of the chiral fermion fields, namely $\Psi_{ij}=\bar{q}_{Rj}q_{Li}$. Employing Eqs. 8 and 7, one then finds that $\Psi$ transforms as follows under the (global) symmetries of massless QCD:

[TABLE]

These transformations show that the second term of the Lagrangian is also invariant under the full group $\mathcal{G}^{\prime}$. The last term, the ’t Hooft determinant, is trivially $\text{U(}1\text{)}_{V}$ invariant. It also exhibits an $\text{SU(}3\text{)}_{L}\times\text{SU(}3\text{)}_{R}$ symmetry, because the $\text{SU(}3\text{)}$ operators $U_{L}$ and $U_{R}$ have unit determinant. In contrast, Eq. 9 reveals that $\Psi$ transforms non-trivially under $\text{U(}1\text{)}_{A}$, which is therefore broken by the ’t Hooft determinant. Thus, the NJL Lagrangian of Eq. 6 has all of the symmetry properties we require for the NJL model to be a suitable effective theory describing the low-energy regime of QCD.

To be able to perform explicit calculations in the NJL model we use the (self-consistent) mean field approximation (MFA) [83, 84, 85, 81]. Within the MFA one defines the expectation value of the fermion bi-linear $\langle\Psi\rangle$ as a sum over effective meson fields, namely

[TABLE]

The $\text{SU(}3\text{)}$ generators in the fundamental representation $T^{a}$ are normalized so as to satisfy $\operatorname{Tr}(T^{a}T^{b})=\delta_{ab}/2$ and fulfill the appropriate $\text{SU(}3\text{)}$ algebra relations

[TABLE]

with the usual antisymmetric and symmetric $\text{SU(}3\text{)}$ structure constants $f_{abc}$ and $d_{abc}$. Correspondingly, the indices $a$, $b$ and $c$ run from $1$ to $n_{f}^{2}-1=8$. The meson fields introduced into the theory via Eq. 10 are to be understood as bosonic auxiliary fields meaning that they do not have tree-level kinetic terms and are thus classically non-propagating. However, fermionic quantum fluctuations will eventually induce kinetic terms and hence render the fields dynamical. More details regarding the MFA can be found in Appendix B. In particular, the final form of the mean field Lagrangian in terms of the fermion and effective meson fields is given in Eq. 61.

In order to determine the dynamics of the chiral phase transition the model’s effective potential is required. Note that only the $\sigma$ field acquires a finite vev during the $\chi$PT.444The effective scalar meson fields are defined in terms of their constituent fermions as $\sigma\sim\langle\bar{q}\mathds{1}q\rangle$ and $a_{a}\sim\langle\bar{q}T_{a}q\rangle$. The isospin symmetry $\text{SU(}3\text{)}_{V}$ present in the Lagrangian is not spontaneously broken, so that all mesons whose definition contains generators with off-diagonal elements, i.e. , terms which couple fermions of different flavor, cannot acquire a finite vev. The only mesons whose vev can be non-zero are therefore $\sigma$, $a_{3}$ and $a_{8}$. Furthermore, we assume all quark masses to be equal. The vacuum should also respect this symmetry, so that $\langle\bar{u}u\rangle=\langle\bar{d}d\rangle=\langle\bar{s}s\rangle$ must hold. Since $\sigma$ is the only meson proportional to the identity matrix, it is the only one which can obtain a finite vev. Note that pseudoscalars cannot acquire a non-zero vev because the vacuum is parity even. The effective potential will therefore only be a function of the corresponding classical field $\bar{\sigma}$, which is defined as the expectation value of $\sigma$ in the presence of some external source. The tree-level potential $V_{0}^{\text{\tiny{}}}$ can be found in Appendix B. Apart from the tree-level terms the one-loop vacuum and thermal contributions, $V_{\text{\tiny{CW}}}^{\text{\tiny{}}}$ and $V_{\text{\tiny{FT}}}^{\text{\tiny{}}}$, are required and are readily obtained by integrating out the fermions. Since all mesons are non-propagating at tree-level they do not contribute to the one-loop effective potential: only the fermions contribute. The finite-temperature effective potential of the NJL model at one loop can thus be written as

[TABLE]

with

[TABLE]

Using the finite-temperature effective potential one can determine the nature of the $\chi$PT by computing the potential’s global minimum $\bar{\sigma}_{\text{min}}(T)$ for each temperature, where $v_{\sigma}=\bar{\sigma}_{\text{min}}(T=0)$. It can be shown that for suitable $G$ and $G_{D}$ the transition is indeed of first order.

As we will see in Section III, studying the dynamics of the $\chi$PT necessitates the calculation of the probability that the scalar field which drives the transition, $\sigma$, tunnels from the theory’s false vacuum to its true vacuum. Recall, however, that all meson fields are classically non-propagating in the NJL model. Since tunneling processes are inherently dynamical, tunneling of the $\sigma$ field can only be described in the current setup, once we compute the $\sigma$ field’s kinetic term as it is induced via fermionic quantum fluctuations. At this stage of the approximation, the auxiliary field $\sigma$ is thus promoted to a propagating quantum field. In doing so, the crucial quantity is the (finite-temperature) wave-function renormalization, which is defined as

[TABLE]

where $\Gamma_{\sigma\sigma}(p^{2},\sigma)$ is the $\sigma$ field’s one-loop propagator, cf. Eq. 64a. Note that in order to determine $Z_{\sigma}^{-1}(\sigma)$ this propagator is to be evaluated at the given field value $\sigma$, but not necessarily at $v_{\sigma}$. Additionally, the loop integrals in the expression for the propagator quoted in Appendix B are generally evaluated at finite temperature when calculating $Z_{\sigma}^{-1}$, and at zero temperature to compute the meson pole masses (see below). The full expression for $Z_{\sigma}^{-1}$ can be found in Eq. 65.

Eventually, we want to compare results about gravitational waves as obtained from the different considered effective models, provided that they predict the same meson mass spectra as well as the same pion decay constant $f_{\pi}$. In the NJL model, the masses of the mesons – in the broken phase – are derived from the roots of their propagators, which, in turn, are determined from the mean field Lagrangian $\mathcal{L}_{\text{\tiny{MFA}}}$, see Appendix B. The one-loop meson propagators can be found in Eqs. 64a to 64d. The pion decay constant on the other hand is given by (see e.g. [80, 86])

[TABLE]

with the constituent quark mass $M_{c}=M(v_{\sigma})$, which corresponds to the effective fermion mass defined in Eq. 13d evaluated at the global minimum of the effective potential at zero temperature. The integral $I_{0}(v_{\sigma})$ can be found in Appendix A.

Finally, note that the NJL model only has three free parameters in the chiral limit, namely $\Lambda$, $G$ and $G_{D}$. Since we always fix $\Lambda=$0.93\text{,}\mathrm{G}\mathrm{e}\mathrm{V}$$ only two free parameters remain. Therefore, not all of the physical parameters $m_{\sigma}$, $m_{a}$, $m_{\eta^{\prime}}$ and $f_{\pi}$ are independent of each other. Using $f_{\pi}$ and $m_{\sigma}$ as input to determine $G$ and $G_{D}$, the other meson masses ($m_{a}$ and $m_{\eta^{\prime}}$) become predictions of the model. The pion mass $m_{\pi}$ is zero since we are working in the chiral limit.

II.2 The Polyakov loop enhanced NJL model

Real QCD exhibits not only a chiral PT, but also a (de)confinement PT, whose dynamics is governed by the gluons. Both of the aforementioned phase transitions are known to occur at similar temperatures from lattice QCD, see e.g. Ref. [87]. Therefore, gluon dynamics might also be relevant to the physics of the chiral phase transition. The NJL model discussed above does not include any gauge dynamics and is consequently unable to describe confinement. This issue was remedied in 2004 by Fukushima [70]. He incorporated gluon dynamics into effective QCD models at finite temperature by adding the Polyakov loop. This resulted in the so-called Polyakov-enhanced NJL (PNJL) model. A recent review of the Polyakov loop can be found, for instance, in Ref. [87]. The Lagrangian of the PNJL model (again in the mean field approximation) is given by

[TABLE]

where $\mathcal{L}_{\text{\tiny{NJL}}}^{\text{\tiny{MFA}}}$ can be found in Eq. 61 and $V_{\text{glue}}$ is the background gluon potential in terms of the deconfinement transition’s (pseudo) order parameter $L$, which denotes the expectation value of the traced Polyakov loop [87].555Our notation is slightly different from the one typically used in the literature (e.g. Ref. [87]), where the expectation value of the traced Polyakov loop is denoted by $\Phi$. Here, we use $L$ to avoid confusion with the meson field $\Phi$. Unfortunately, the gluon potential cannot yet be directly calculated from first principles. Still, the impact of including gauge dynamics can be effectively captured by appropriately parametrizing $V_{\text{glue}}$. Data from, for example, lattice QCD can then be used to fix the values of the thus introduced additional free parameters. A discussion on different parametrizations can be found in Ref. [88], where it is also shown that the possible choices for the gluon potential are equivalent for temperatures below two or three times the critical temperature of the deconfinement transition. The dynamics of the $\chi$PT is hence not expected to depend strongly on the choice of potential. In this paper we use the gluon potential first introduced in Ref. [89], namely

[TABLE]

with

[TABLE]

The potential is chosen such that the value of $L$ is constrained to lie between [math] and $1$, with the deconfined and confined phases being reached for $L\to 1$ and $L\to 0$, respectively. Furthermore, since we are always working at zero chemical potential, $L$ is real and thus $L=\bar{L}$ [87]. The critical temperature $T_{\text{glue}}$ of the QCD deconfinement PT in the pure gauge limit, i.e. assuming infinitely heavy quarks, is known to be $T_{\text{glue}}=$270\text{,}\mathrm{M}\mathrm{e}\mathrm{V}$$ from lattice simulations [90]. However, when taking into account three mass-degenerate quark degrees of freedom as in our model, $T_{\text{glue}}$ is reduced to $178\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ [90]. Lastly, the dimensionless parameters $a_{i}$ and $b_{3}$ in Eq. 18 are determined from lattice QCD calculations as well [89], namely,

[TABLE]

Including the Polyakov loop modifies the thermal part of the NJL model’s effective potential from Eq. 13c in two ways. First, the background gluon potential $V_{\text{glue}}$ of Eq. 17 is to be added. Second, the quark coupling to gluons must be taken into account. In contrast, the vacuum part of the effective potential is unchanged with respect to the NJL model. In particular, the meson propagators and thus also the meson masses are therefore the same to before. Combining all of the above, the final expression for the PNJL model’s thermal one-loop effective potential is

[TABLE]

with

[TABLE]

where $r\equiv r(\bar{\sigma},T)=M(\bar{\sigma})/T$ with $M$ denoting the effective fermion mass of the NJL model as introduced in Eq. 13d. The MFA was applied to arrive at the quoted expression for $V_{\text{\tiny{FT}}}^{\text{\tiny{PNJL}}}$ [87].

Using the above potential, we find that the deconfinement transition described by the Polyakov loop $L$ is a crossover, while the $\chi$PT with order parameter $\bar{\sigma}_{\text{min}}$ is of first order for all considered benchmark points. To proceed we again need to determine how the $\sigma$ field666The wave-function renormalization $Z_{\sigma}$ is the same in the NJL and the PNJL model. tunnels from the theory’s false vacuum to its true vacuum. In doing so, we choose $L$ to minimize the effective potential for any given point $(\bar{\sigma},T)$, i.e. $L:=L_{\text{min}}(\bar{\sigma},T)$. This will reduce $V_{\text{eff}}^{\text{\tiny{PNJL}}}$ to be only a function of $\bar{\sigma}$ and $T$, thus allowing us to use the standard formalism for tunneling in one field dimension, which will be outlined in Section III. Note that our choice of $L$ is similar to the one made in Ref. [91].

II.3 The linear sigma model

Just as the Nambu–Jona-Lasinio (NJL) model, the linear sigma model (LSM), as first introduced by Gell-Mann and Levy [71], is a phenomenological theory to effectively describe certain non-perturbative aspects of the low-energy dynamics of QCD or related strongly coupled theories. Unlike the NJL model the LSM is a renormalizable, purely scalar theory with the scalar fields being identified with the well-known mesons of QCD. There exist several variants of the LSM which differ by the global internal symmetry that they are based on. Since we are interested in a QCD-like theory with $n_{f}=3$ massless fermion flavors, we will here discuss the chiral limit of the $\text{U(}3\text{)}\times\text{U(}3\text{)}$ LSM, whose Lagrangian reads as follows:

[TABLE]

where the tree-level potential $V_{\text{tree}}$ is given by (see e.g. [92])

[TABLE]

The complex-valued bosonic field $\Phi$ is a $3\times 3$ matrix defined as

[TABLE]

and thus collects both the scalar meson fields $(\sigma,\vec{a})$ and the pseudoscalar ones $({\eta^{\prime}},\vec{\pi})$. The $\text{SU(}3\text{)}$ generators $T^{a}$ were introduced below Eq. 10, and the index $a$ runs again from 1 to $n_{f}^{2}-1=8$.

Under a general chiral $\mathcal{G}^{\prime}=\text{U(}3\text{)}_{L}\times\text{U(}3\text{)}_{R}$ transformation, the bosonic field $\Phi$ behaves as

[TABLE]

so that $\Phi_{ij}$ has the same quantum numbers as the fermion bilinear $\bar{q}_{Rj}q_{Li}$, where $q$ are the matter fields introduced in Eq. 3. It is straightforward to derive from Eqs. 24 and 25 that $\sigma$ and ${\eta^{\prime}}$ transform as singlets under the vector subgroup $\text{U(}3\text{)}_{V}$ of $\mathcal{G}^{\prime}$, whereas $\vec{a}$ and $\vec{\pi}$ are in its adjoint representation. Using Eq. 25, one furthermore finds that the kinetic term in Eq. 22 and all operators in the first line of Eq. 23 are invariant under $\mathcal{G}^{\prime}$. In fact, with the exception of the $\lambda_{a}$-term, all of the aforementioned terms are even symmetric under the larger group $\text{O(}2n_{f}^{2}\text{)}$. In contrast, the last term in Eq. 23 is only invariant under $\mathcal{G}=\text{SU(}3\text{)}_{V}\times\text{SU(}3\text{)}_{A}\times\text{U(}1\text{)}_{V}$ and thus accounts for anomalous $\text{U(}1\text{)}_{A}$ breaking [82, 73]. Spontaneous chiral symmetry breaking $\mathcal{G}\to\text{SU(}3\text{)}_{V}\times\text{U(}1\text{)}_{V}$ is finally triggered by the scalar singlet $\sigma$ acquiring a finite vacuum expectation value, i.e. , $\sigma\to v_{\sigma}+\sigma$, or in terms of the $\Phi$ field from Eq. 24

[TABLE]

where the latter relation is e.g. demonstrated in Ref. [92].

It is well-known that the LSM at finite temperature is plagued by infrared divergences which lead to the breakdown of standard perturbation theory [93]. Many-body resummation schemes are therefore necessary for the purpose of obtaining robust results. In the following, we will employ the composite operator or 2PI formalism due to Cornwall, Jackiw and Tomboulis (CJT) [94] to implement the so-called Hartree-Fock approximation [95, 96]. Details of the corresponding calculations are relegated to Appendix C and can e.g. also be found in Refs. [97, 98, 99, 100].

Here, we only quote the final result for the $\text{U(}3\text{)}\times\text{U(}3\text{)}$ LSM finite-temperature effective potential, which reads

[TABLE]

where $\bar{\sigma}$ denotes the expectation value of the $\sigma$ field in the presence of an external source, such that $\bar{\sigma}\to v_{\sigma}$ as the source vanishes. The individual components in Eq. 26 are given by

[TABLE]

Note that here we neglect vacuum contributions, which was previously demonstrated to not change the outcome qualitatively [98, 99]. A few further comments regarding Eq. 27 are in order. First, the sum in Eq. 27b runs over all meson species introduced in Eq. 24, i.e. , $i\in\{\sigma,{\eta^{\prime}},\pi,a\}$. The corresponding statistical weights $g_{i}$ are $g_{\sigma}=g_{\eta^{\prime}}=1$ and $g_{\pi}=g_{a}=n_{f}^{2}-1=8$. Moreover, we have defined $r_{i}\equiv r_{i}(\bar{\sigma},T)=m_{i}(\bar{\sigma})/T$ and $R_{i}\equiv R_{i}(\bar{\sigma},T)=M_{i}(\bar{\sigma},T)/T$. Here, the $m_{i}(\bar{\sigma})$ denote the field-dependent tree-level meson masses [92],

[TABLE]

while the $M_{i}(\bar{\sigma},T)$ constitute effective meson masses, which are self-consistently calculated within the CJT formalism for given $T$ and $\bar{\sigma}$ as described in Appendix C. Lastly, the thermal integrals $J_{\text{\tiny{B}}}$ and $I_{\text{\tiny{B}}}$ are defined in Eq. 57 of Appendix A. Formally, the first term in Eq. 27b reproduces the one-loop term of the conventional perturbative series for the effective potential at finite temperature, see e.g. [101]. However, the tree-level masses $m_{i}$ usually appearing as an argument of the $J_{\text{\tiny{B}}}$ are replaced by the effective masses $M_{i}$. The second term goes beyond standard perturbation theory, as well. In total, it can be shown that the CJT effective potential in the Hartree-Fock approximation accounts for the resummation of all daisy and super-daisy diagrams [102].

Let us finally mention that for the $\text{U(}3\text{)}\times\text{U(}3\text{)}$ LSM in the chiral limit, i.e. $m_{\pi}=0$, the number of free Lagrangian parameters $(m^{2},\lambda_{\sigma},\lambda_{a},c)$ equals the number of observable vacuum quantities $(f_{\pi},m_{\sigma},m_{\eta^{\prime}},m_{a})$. Hence, evaluating Eq. 28 at the physical vacuum, i.e. for $\bar{\sigma}=\sqrt{3/2}\,f_{\pi}$, yields a system of equations that can be uniquely solved for the former set of parameters given the pion decay constant and a full set of meson masses.

We will use the finite-temperature effective potential of Eq. 26 for all further calculations related to the linear sigma model.

III Chiral phase transition

As we have mentioned in the beginning of Section II, the chiral phase transition (PT) in QCD-like theories with three (or more) massless fermion flavors is predicted to be of first order on general grounds [78]. For the low-energy effective models introduced in the previous section, we can explicitly demonstrate the existence of such a first-order chiral PT by using the appropriate finite-temperature effective potentials in Eqs. 12, 20 and 26, respectively. To that end, we examine how the potential’s global minimum $\bar{\sigma}_{\text{min}}$, which, in turn, determines the given theory’s true groundstate, changes with temperature. Fig. 1 shows the function $\bar{\sigma}_{\text{min}}(T)$ as calculated within the investigated low-energy models for one of the benchmark points collected in Table 1. The appearance of two degenerate minima at some critical temperature $T_{c}$ clearly signals the occurrence of a first-order phase transition. As an aside, note that the figure also explicitly demonstrates the equivalence of the NJL and the PNJL model at zero temperature.

From a physics point of view, it is well known that a first-order PT proceeds via the nucleation and subsequent growth of bubbles of the true vacuum inside an expanding universe, which is still in the metastable symmetric phase. Hence, the transition’s properties crucially depend on both the nucleation rate $\Gamma$ of the aforementioned bubbles and the Hubble parameter $H$. In what follows, we will briefly introduce these two quantities in turn.

First, the temperature-dependent bubble nucleation rate $\Gamma$, which is also referred to as the decay rate of the false vacuum, can be calculated as [103, 5, 4]

[TABLE]

In the above equation, the theory’s three-dimensional Euclidean action $S_{3}$ is to be understood as being evaluated for the $\text{O(}3\text{)}$-symmetric bounce or tunneling solution $\bar{\sigma}_{b}(r)$, i.e. $S_{3}\equiv S_{3}[\bar{\sigma}_{b}(r)]$. The field configuration $\bar{\sigma}_{b}(r)$, in turn, is determined as the solution of the scalar background field’s equation of motion subject to the boundary conditions $\bar{\sigma}\to 0$ as $r\to\infty$ and $\mbox{d}\bar{\sigma}/\mbox{d}r=0$ at $r=0$, where $r$ is the radial coordinate of three-dimensional space. If the $\bar{\sigma}$ field represents a fundamental degree of freedom, as is e.g. the case in the LSM, the theory’s Euclidean action can be written as

[TABLE]

so that the background field’s equation of motion reads

[TABLE]

In contrast, the $\bar{\sigma}$ field is composite in the (P)NJL model (i.e. non-propagating at tree-level). The Euclidean action therefore has to be slightly modified with respect to Eq. 30 and is now given by

[TABLE]

where $Z_{\sigma}\equiv Z_{\sigma}(\bar{\sigma})$ is the wave-function renormalization introduced in Eq. 14. The background field’s equation of motion is then found to be

[TABLE]

Throughout the present work, we solve Eqs. 31 and 33 using the CosmoTransitions package [104] and an appropriately customized version thereof, respectively. The same code is employed to compute the quantity $S_{3}[\bar{\sigma}_{b}(r)]$ based on Eqs. 30 and 32.

Next, let us discuss the expansion of the universe as it is governed by the temperature-dependent Hubble parameter $H$. In the absence of significant supercooling, the universe will be radiation dominated during the whole phase transition so that $H(T)$ is given by

[TABLE]

where $M_{\text{\tiny{Pl}}}$ is the reduced Planck mass, $M_{\text{\tiny{Pl}}}=$2.435\text{\cdot}{10}^{21}\text{,}\mathrm{M}\mathrm{e}\mathrm{V}$$, while the effective number of relativistic degrees of freedom in the thermal plasma is denoted as $g_{\star}$.

Given the bubble nucleation rate and the Hubble parameter of Eqs. 29 and 34, the nucleation temperature $T_{n}$ is defined as the temperature, for which one bubble per Hubble volume is produced on average, i.e. ,

[TABLE]

Typically, the above integral is dominated by temperatures very close to $T_{n}$. Hence, the defining condition in Eq. 35 is approximately equivalent to

[TABLE]

where we used Eqs. 29 and 34 to arrive at the second relation, ignoring the slowly varying factor $(S_{3}/(2\pi T))^{3/2}$ in the expression for the nucleation rate $\Gamma$. Notably, the right-hand side of Eq. 36 depends on the temperature only logarithmically, so that the nucleation condition roughly translates to $S_{3}/T\approx\mathcal{O}(100)$ for a wide range of temperatures [5, 4].

The gravitational wave spectrum associated with a first-order PT, however, not only depends on the nucleation temperature, but also on various other quantities that characterize the transition’s properties. First, the phase transition’s inverse duration $\beta$ compared to the expansion rate of the Universe at the time of the transition is a crucial parameter that can be calculated as777Strictly speaking, Eqs. (37) to (39) are to be evaluated at the transition temperature – usually denoted as $T_{*}$ in the literature – rather than at the nucleation temperature $T_{n}$. However, as we will explicitly see in Section V, the chiral phase transitions under consideration exhibit virtually no supercooling, so that the thermal plasma is practically not reheated after the PT has completed and $T_{n}\simeq T_{*}$ holds to good approximation. On a related note, $H$ will throughout the paper refer to the Hubble parameter at the time of the transition unless stated otherwise.

[TABLE]

In the absence of a calculable model Eq. 37 cannot be evaluated without further assumptions on $S_{3}/T$. One way is then to follow Refs. [5, 4] and approximate $\mbox{d}S_{3}/\mbox{d}T\gtrsim S_{3}/T$. In that case $\beta/H$ can be estimated as

[TABLE]

where the last relation holds for a wide range of transition temperatures $T_{n}$. Since first-principle calculations of $\beta/H$ are lacking for theories exhibiting a chiral phase transition in a strongly coupled hidden sector, the rough estimate in Eq. 38 was employed for a long time in the corresponding literature, see e.g. [65, 66, 16]. In contrast, the use of effective models in the present work permits us to determine the function $S_{3}(T)/T$, so that we can evaluate Eq. 37 explicitly.

A second important quantity in the calculation of gravitational wave spectra is usually denoted as $\alpha$ and encodes the energy released during the phase transition. Correspondingly, it can be regarded as a measure of the PT’s strength. Several definitions of $\alpha$ can be found in the literature, including one in terms of the transition’s latent heat normalized to the radiation energy density (see e.g. [26]),

[TABLE]

Notably, $\alpha_{L}$ and $\alpha_{T}$ are equivalent for strong phase transitions with large amounts of supercooling, i.e. $\Delta V_{\text{eff}}\gg T_{n}\partial\Delta V_{\text{eff}}(T)/\partial T$. Since the considered $\chi$PTs will, however, turn out to be weakly first order, a definition for $\alpha$ needs to be chosen. In the remainder of the paper, we will always set $\alpha:=\alpha_{T}$ as required by our use of Eq. 42 from Ref. [105].

IV Stochastic gravitational wave signal

The production of gravitational waves (GW) during violent astrophysical processes like black hole mergers is a well-known phenomenon, which is predicted by Einstein’s theory of General Relativity [106]. The first GW signals of this kind were recently detected by the LIGO and VIRGO collaborations [1, 2, 3]. Apart from these transient signals due to astrophysical sources, a stochastic background of GWs may exist as a result of, for instance, first-order cosmic phase transitions (PT) or an inflationary phase in the early universe [107]. Observing such a background would thus provide a unique opportunity to investigate primordial physics which would otherwise be impossible to study.

As discussed in Section III, a first-order PT proceeds via the nucleation and subsequent growth of bubbles of the true groundstate. Gravitational waves are produced when the aforementioned bubbles collide via several different processes, namely, collisions of bubble walls or, equivalently, scalar field shells [108], as well as sound waves [109] and magnetohydrodynamic turbulence [110] in the thermal plasma. The relative importance of these contributions is determined by the dynamical properties of the phase transition and hence depends on the underlying particle physics model. Still, it is possible to broadly distinguish between two different transition scenarios which are usually referred to as runaway and non-runaway, see e.g. [67].

In the former case the friction exerted by the thermal plasma on the nucleated bubbles is too small to slow down their accelerated expansion, which, in turn, is driven by the released latent heat. Consequently, the bubble wall velocity $v_{b}$ will continue to increase until it eventually reaches the speed of light $c$. In the non-runaway scenario on the other hand, the friction between bubbles and the surrounding plasma is large enough to counteract the bubbles’ accelerated expansion. The bubble wall velocity hence approaches a terminal value, which is bounded by $c$ but may still be relativistic.

In order to determine which of the above cases applies to the QCD-like hidden sector models discussed in the present paper, the following considerations are helpful: as outlined in Section II, the chiral phase transition in those models is driven by the $\bar{\sigma}$ field. The latter is expected to have sizable interactions with the theory’s remaining degrees of freedom. The friction exerted by the plasma on the bubbles defined by $\bar{\sigma}$ is therefore anticipated to be non-negligible, which already strongly suggests a non-runaway scenario (i.e. Case I of Ref. [67]). Furthermore, we will see in the next section that our calculations predict rather small values for the phase transition strength $\alpha$ defined in Eq. 39b, which also indicates a transition of this type. Accordingly, we will assume in the following that the $\chi$PT proceeds via non-runaway bubbles in all of the considered models and for all of our benchmark points.888Note that, if the $\chi$PT proceeded via runaway bubbles but with the same $\alpha$ and $\beta/H$, the total energy density of the GW background would decrease with respect to the non-runaway case. The reason is that part of the energy which is transformed into bulk motion for non-runaway bubbles will be deposited in the scalar field gradient in a runaway scenario. The associated contribution to the GW spectrum is, however, suppressed by an additional factor of $(\beta/H)^{-1}$ and is thus irrelevant for $\beta\gg H$ [67]. This should be kept in mind in the following.

In non-runaway scenarios, the fraction of the PT’s latent heat that goes into the kinetic energy of the scalar field is negligible, so that the contribution to gravitational wave production from scalar field shell collisions becomes irrelevant [67]. In contrast, and as a result of the large friction between bubbles and the surrounding thermal plasma, a substantial portion of the available energy is converted into bulk motion in the form of sound waves and magnetohydrodynamic turbulence. In the following, we will briefly discuss both of the aforementioned contributions to GW production in turn, starting with that from sound waves. The corresponding energy density is [67]

[TABLE]

where the spectral shape $S_{\text{sw}}(f)$ and the spectrum’s peak frequency $f_{\text{sw}}$ are given by

[TABLE]

respectively. As indicated above, $\alpha$ and $\beta/H$ characterize the strength and inverse duration of the PT (cf. Section III), while $g_{\star}$ is the effective number of relativistic degrees of freedom in the plasma and $\kappa_{v}$ denotes the fraction of latent heat converted into bulk motion. The actual terminal wall velocity $v_{b}$ of non-runaway bubbles is determined by the details of the underlying particle physics model. Calculating its specific value is, however, beyond the scope of the present work. In the following we will instead assume a certain range of different $v_{b}$ values and investigate how the resulting GW spectra change. In doing so we will concentrate on highly relativistic bubble wall velocities, $v_{b}\geq$0.75$$, since those lead to the strongest GW signals and are thus most interesting from an observational point of view. In the considered case of non-runaway bubbles with large wall velocity $v_{b}\lesssim c$ an estimate for the efficiency factor $\kappa_{v}$ is then given by [105]

[TABLE]

That being said, one should keep in mind that the spectrum in Eq. 40 was derived from simulations assuming that GW production via sound waves is sufficiently long-lasting. To be more precise, Eq. 40 is only reliable if $\tau_{\text{sw}}H>1$ [111], where $\tau_{\text{sw}}$ is related to the duration of the sound wave period and reads

[TABLE]

with the plasma’s root-mean-square four-velocity $\bar{U}_{f}$. As we will explicitly see in Section V, the $\chi$PTs under consideration are throughout predicted to proceed very fast, $\beta/H\gtrsim\mathcal{O}(${10}^{4}$)$, so that the condition $\tau_{\text{sw}}H>1$ is unlikely to be satisfied. In order to check this explicitly, we follow Refs. [105, 26] to first calculate $\bar{U}_{f}$. Subsequently, we use Eq. 43 to determine the product $\tau_{\text{sw}}H$ assuming $v_{b}=1$. Indeed, we find $\tau_{\text{sw}}H=\mathcal{O}(10^{-3})$ across all benchmark points and models. From a physics point of view, $\tau_{\text{sw}}H<1$ means that sound waves can efficiently source gravitational waves only over a period shorter than a Hubble time, since the transition proceeds too fast. Correspondingly, Eq. 40 is expected to overestimate the energy density of the produced GWs [26]. The authors of Refs. [26, 36] therefore suggested accounting for a potentially shortened sound wave period in the case of fast PTs by multiplying the amplitude of Eq. 40 by $\tau_{\text{sw}}H$, namely

[TABLE]

Note that although the use of this reduction factor is perfectly reasonable from a physics perspective, it is not yet backed up by dedicated numerical simulations, but is rather subject of current research.

The second relevant source of stochastic gravitational waves in the scenario with non-runaway bubbles is the plasma’s turbulent motion. Its contribution to the total GW amplitude is [67]

[TABLE]

where the spectral shape $S_{\text{tb}}(f)$ and the spectrum’s peak frequency are given by

[TABLE]

with $h_{\star}$ being the Hubble rate at the time of GW production, redshifted to today, namely $h_{\star}=$16.5\text{\cdot}{10}^{-6}\text{,}\mathrm{m}\mathrm{H}\mathrm{z}$\cdot(T_{n}/$100\text{,}\mathrm{M}\mathrm{e}\mathrm{V}$)\cdot(g_{\star}/100)^{\nicefrac{{1}}{{6}}}$. Typically, one assumes that only a rather small part of the bulk motion is turbulent, namely $\kappa_{\text{tb}}=\epsilon\kappa_{v}$ with $\epsilon\approx$5\text{,}%$$, see e.g. [67, 26]. If this is indeed true, the corresponding contribution to GW production will be irrelevant and the stochastic background will arise virtually exclusively from sound waves.

However, as it was pointed out e.g. in Refs. [26, 36], the above assumption on $\epsilon$ is expected to become invalid in the case of very fast transitions: After a shortened sound wave period the plasma is anticipated to enter a regime of non-linear dynamics, so that a relatively large amount of the available energy can be transferred to turbulent motion. The choice $\epsilon\approx$5\text{,}%$$ will therefore probably underestimate the associated GW signal if $\beta/H\gg 1$. To account for the described effect, the authors of Ref. [36] proposed setting $\kappa_{\text{tb}}=\kappa_{v}$ and multiplying the amplitude of Eq. 45 by $(1-\tau_{\text{sw}}H)$,

[TABLE]

Notably, the reduction of the sound wave amplitude quoted in Eq. 44 together with the above-described enhancement of the GW signal due to turbulence results in both contributions being of almost equal order of magnitude.

As already mentioned before, neither Eq. 44 nor Eq. 47 have been derived from numerical simulations. Rather, they are the result of heuristic arguments and thus may change in the future. In the following, we will therefore compute the predicted GW spectra both for the conventional case, i.e. solely based on Eq. 40, and explicitly taking into account the (short) duration of the transition via Eqs. 44 and 47.

V Results

Having summarized the basics of first-order (chiral) phase transitions and the associated gravitational wave signal in the previous two sections, we will now apply the described formalism to derive predictions for the gravitational wave spectrum within the low-energy effective models introduced in Section II. Let us once more stress that using the aforementioned models puts us in the position to explicitly calculate the parameters characterizing the chiral phase transition (e.g. $\beta/H$) without having to resort to heuristic arguments or rough estimates. We will do so in the following, starting with the benchmark points introduced in Table 1, for which the transition is expected to occur around the QCD scale. Building on these results we will, in a second step, also determine the gravitational wave signal originating from hidden chiral phase transitions at (much) higher temperatures. For the reasons outlined at the end of the last section, the gravitational wave spectra obtained from Eq. 40 will throughout be compared to those that are based on Eqs. 44 and 47.

V.1 Hidden chiral phase transitions at $\mathcal{O}(100\,\mathrm{MeV})$

Employing the finite-temperature effective potentials of Eqs. 12, 20 and 26, we can directly compute the bounce solution $\bar{\sigma}_{b}(r)$ for a given temperature $T$. Eqs. 30 and 32 then allow us to determine $S_{3}/T$ as a function of the dimensionless ratio $T/T_{c}$, where the critical temperatures $T_{c}$ for each benchmark point and model are listed in the first column of Table 2. The result of such a calculation is shown in Fig. 2 for one of our benchmark points. Interestingly, we observe that $S_{3}/T$ quite precisely follows a function of the simple form

[TABLE]

which was already anticipated by Hogan in Ref. [5]. Fitting our data points to the function in Eq. 48, the parameters $\gamma$ and $b$ can be determined, see the last two columns of Table 2. Intriguingly, the exponent $\gamma$ is very similar among all benchmark points and effective models. In contrast, the coefficient $b$ is found to vary by several orders of magnitude. Using the best-fit values for $\gamma$ and $b$, we plot the function of Eq. 48 alongside the data in Fig. 2 (solid lines). For comparison, we also show natural cubic splines interpolating between the data points (dashed lines). All further calculations involving $S_{3}/T$ will, however, make use of the above described fit, since this is less sensitive to numerical errors.

Given $S_{3}/T$ as a function of $T$ we can now straightforwardly compute the nucleation temperature $T_{n}$ by solving Eq. 36. Most importantly, we always find that $T_{n}\lesssim T_{c}$ implying that there exists virtually no supercooling, which a posteriori justifies our choice of the Hubble parameter in Eq. 34, as well as our assumption that the nucleation temperature approximately coincides with the transition temperature (cf. footnote 7 on page 7).

Next, the phase transition’s inverse duration $\beta$ normalized to the expansion rate of the Universe at the time of the transition can be determined from Eq. 37. The outcome is again compiled in Table 2. Although the precise values for $\beta/H$ vary between different models and benchmark points, we generally find them to be of order ${10}^{4}$ or even larger. Note that these results are in stark contrast to the usual assumption999See also Eq. 1 as well as our discussion in Section III after Eq. 37. $\beta/H\approx\mathcal{O}(100)$, which is being made in discussions of (hidden) chiral phase transitions throughout the literature, see e.g. [16, 66].101010There have been earlier indications that $\beta/H$ can attain much larger values than $\mathcal{O}(100)$, e.g. in Refs. [21, 22, 28]. As we will see in more detail below, this has far-reaching consequences on the observational prospects for the associated gravitational wave signal. Employing the functional form for $S_{3}/T$ in Eq. 48, we can even learn why the conventionally used approximation for $\beta/H$ fails. For that purpose, we compute the derivative

[TABLE]

As stated earlier, we always find $T_{n}$ to be very close to $T_{c}$ so that the second term in the above expression dominates and $\mbox{d}S_{3}/\mbox{d}T$ is much larger than $S_{3}/T$. Thus, the central assumption in arriving at $\beta/H\approx S_{3}/T\approx\mathcal{O}(100)$ via Eq. 38 is no longer justified. Alternatively, we can also use Eqs. 37 and 48 to calculate $\beta/H$ directly, giving

[TABLE]

which exhibits the same enhancement with respect to $S_{3}/T$ for $T_{n}$ near $T_{c}$. Since we generally expect Eq. 48 to be valid in the vicinity of the critical temperature, we suspect large $\beta/H\gg S_{3}/T$ to be a generic feature of models without large supercooling. As the quantity $S_{3}/T$ evaluated at $T_{n}$ varies only very little for a wide range of nucleation temperatures (cf. Eq. 36), the exact value of $\beta/H$ then chiefly depends on the ratio $T_{n}/T_{c}$, as well as on the exponent $\gamma$. This behavior is clearly visible in Table 2, in particular for our benchmark point D. Here, the (P)NJL model predicts practically degenerate $T_{n}$ and $T_{c}$, which leads to huge values for $\beta/H$ of order ${10}^{6}$.

Lastly, the transition strength $\alpha$ introduced in Section III is obtained using Eq. 39b, see Table 2 for our results. Let us note that, whereas $T_{n}$ and $\beta/H$ depend only very mildly on the effective number of relativistic degrees of freedom, $\alpha$ is by definition rather sensitive to the precise value of $g_{\star}$. By Eqs. 40 and 42 the same is then true for the energy density $\Omega h^{2}$ of the GW background (but not for its peak frequency, cf. Eq. 41). To be more precise, it is straightforward to see that increasing $g_{\star}$ effectively decreases $\Omega h^{2}$. For definiteness and unless explicitly stated otherwise, we have fixed $g_{\star}=$47.5$$ in all calculations, which corresponds to $n_{f}=3$ light hidden fermions and $n_{c}^{2}-1=8$ hidden gluons. For models with different hidden sectors the presented results must be appropriately adapted.

With the quantities $T_{n}$, $\beta/H$ and $\alpha$ at hand, and assuming that the chiral phase transition proceeds via non-runaway bubbles with some terminal wall velocity $v_{b}$, we are now able to compute the predicted gravitational wave spectra $\Omega h^{2}$ as described in Section IV. Our findings for the benchmark points from Table 1 are displayed in Fig. 3. Recall that both the position and height of the spectrum’s peak are solely determined by the sound wave contribution to $\Omega h^{2}$.

Importantly, Eq. 41 reveals that values of $\beta/H\approx\mathcal{O}(${10}^{4}$)$ imply a peak frequency at approximately $1\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$ for bubble nucleation at roughly $100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$. A first-order chiral phase transition at such temperatures is consequently predicted to produce a stochastic GW background within the LISA frequency band. At the same time, a large $\beta/H$ suppresses the peak height according to Eqs. 40, 44 and 45. Unfortunately, Fig. 3 seems to suggest that, indeed, the GW signal is too weak to be detected by LISA. However, comparing the predictions of different effective models for a given benchmark point (i.e. a fixed mass spectrum), we find that the parameters characterizing the chiral phase transition as well as the resulting gravitational wave spectrum may vary considerably. Also the variance across our benchmark points is sizable. We therefore conclude that our effective model analysis ultimately fails to conclusively answer whether the GW signal can be expected to be observed by LISA. Consequently, first-principle calculations, like lattice simulations, are required in order to get quantitatively robust predictions.

V.2 Hidden chiral phase transitions at higher temperatures

In the previous section, we exclusively studied benchmark points with an inherent scale of order $100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$. However, as we have argued in the introduction, there are equally well motivated BSM scenarios, where a hidden chiral phase transition is anticipated to occur at higher temperatures. In order to investigate those cases as well, we will now simply consider scaled-up versions of the benchmark points in Table 1, obtained by multiplying all meson masses etc. by a common dimensionless factor $\xi>1$, i.e. ,

[TABLE]

Since it is fully determined by model parameters, the critical temperature $T_{c}$ will then scale with the same factor. For the same reason the action functional as a function of $T/T_{c}$ will remain unaltered, i.e. ,

[TABLE]

In contrast, the nucleation temperature $T_{n}$ only approximately scales with $\xi$. An exact scaling is violated by the fact that Eq. 36 contains the Planck mass as an absolute energy scale, which is kept constant. Accordingly, the ratio $T_{n}/T_{c}$ is observed to decrease mildly for growing $\xi$, while it always stays close to one. A small change in $T_{n}/T_{c}$ can, however, have a sizable effect on $\beta/H$ for the rescaled point since $S_{3}/T$ has a pole at $T=T_{c}$, close to which its derivative varies rapidly, see Fig. 2 and Eq. 50. The resulting scale dependence of $\beta/H$ is displayed in Fig. 4. Note that $\beta/H$ stays of order ${10}^{4}$ for a wide range of critical temperatures in all considered models. Lastly, in line with its definition in Eq. 39b the parameter $\alpha$ is observed to not depend strongly on the precise value of $T_{n}/T_{c}$, so that it can be safely assumed as constant in $\xi$ for any fixed choice of $g_{\star}$.

Following the above discussion and based on the results compiled in Table 2, we can now compute the GW signal for any given $\xi$. Importantly, the spectrum’s peak frequency is proportional to the nucleation temperature and will thus approximately scale linearly with $\xi$, cf. Eq. 41. For hidden chiral phase transitions occurring between roughly $1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and $10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ (corresponding to $\xi$ ranging from $10$ to ${10}^{5}$) the associated GW backgrounds are therefore anticipated to lie within the DECIGO/BBO frequency band. If we, for the moment, ignore any potential additional suppression of the GW amplitude due to a shortened sound wave period in very fast transitions, our effective model study strongly suggests that the considered experiments are, in fact, sufficiently sensitive for the stochastic background to be observed, cf. the solid bands in Figs. 5 and 6. Conversely, the inclusion of the aforementioned effect via the rough estimates in Eqs. 44 and 47 drastically reduces the GW signal, so that a detection becomes less likely.

The above described analysis can, of course, be repeated for larger scaling factors $\xi$ in order to obtain results appropriate to hidden chiral phase transitions at even higher temperatures. Crucially, we find that for $T_{n}\gtrsim$10\text{,}\mathrm{T}\mathrm{e}\mathrm{V} (or, equivalently, $\xi\gtrsim${10}^{5}) the considered models unanimously predict GW spectra with peak frequencies outside of the sensitivity bands of all planned space-based gravitational wave observatories, cf. Fig. 7. A stochastic GW background originating from such high-temperature $\chi$PTs is therefore unlikely to be detectable.

In order to make all of the above qualitative statements slightly more precise, let us finally quantify the discovery potential of the considered experiments by calculating so-called signal-to-noise ratios (SNR), which are defined as (see e.g. Ref. [112])

[TABLE]

A few comments on the quantities appearing in Eq. 53 are in order. First, $\Omega_{\text{\tiny{GW}}}(f)h^{2}$ represents a gravitational wave signal defined by the stochastic background spectra introduced in Section IV and shown in Figs. 3 and 5 to 7. On the other hand, $\Omega_{\text{noise}}(f)h^{2}$ denotes a given observatory’s effective strain noise power spectral density, expressed as energy density parameter [117]. As already mentioned in the caption of Figs. 3 and 5, we adopt strain noise power spectral densities from Refs. [113, 115, 114]. Frequency integration in Eq. 53 is performed over the respective experiment’s bandwidth defined by $f_{\text{min}}$ and $f_{\text{max}}$. The duration of an observation in seconds is denoted as $t_{\text{obs}}$ and will be set to five years in the following. Let us finally remark that a more complete treatment should also take into account the effect of unresolvable astrophysical foregrounds from black hole, neutron star, and white dwarf mergers on the signal significance, which is, however, beyond the scope of the current analysis. The resulting signal-to-noise ratios for one of our benchmark points are displayed in Fig. 8 and confirm our qualitative discussion from before.

VI Summary and conclusions

In the present paper, we investigated the prospects for observing a stochastic gravitational wave background associated with a first-order chiral phase transition (PT) in a strongly coupled hidden or dark sector. For definiteness we studied a new confining $\text{SU(}n_{c}\text{)}$ gauge force with $n_{c}=3$ colors acting on $n_{f}=3$ flavors of massless fermions. In order to reliably model the strong dynamics that drive the aforementioned transition, we employed various effective theories known to capture non-perturbative aspects of the low-energy regime of real-world QCD: the Nambu–Jona-Lasinio model, its Polyakov-loop enhanced variant and the linear sigma model.

This approach allowed us to explicitly calculate the temperature-dependent action functional $S_{3}/T$, which is, in turn, crucial to determine key characteristics of the hidden PT like its duration or the nucleation rate of critical bubbles of the true vacuum. We found that, in the vicinity of its pole at the critical temperature $T_{c}$, the function $S_{3}(T)/T$ is of the simple form $b(1-T/T_{c})^{-\gamma}$ for some $b,\gamma>0$. Intriguingly, the exponent $\gamma$ is very similar among all considered benchmark points and effective models, namely $\gamma\simeq$1.81.9$$.111111Interestingly, the found values are very close to the critical exponent $\gamma$ of the mean cluster size in standard percolation theory in three dimensions [118, 119, 120]. Even though the discussion of a possible connection is beyond the scope of the present article, it might be worthwhile to investigate it in the future.

As a next step, we then computed the phase transition’s nucleation temperature $T_{n}$, as well as its inverse duration normalized to the Hubble parameter $\beta/H$. Our results consistently predict the transition to proceed very fast with virtually no supercooling, $T_{n}\lesssim T_{c}$. Crucially, our findings of $\beta/H\gtrsim\mathcal{O}(${10}^{4}$)$ are in stark contrast to $\beta/H\approx$100$$ which is typically assumed in the literature whenever explicit calculations are unavailable. Lastly, we demonstrated that the observed enhancement of $\beta/H$ by two orders of magnitude has significant impact on the predicted gravitational wave spectrum. The latter was determined for various transition temperatures appropriate for different well-motivated beyond-the-Standard Model scenarios featuring a new strongly coupled sector. Our study suggests that a gravitational wave background associated with a hidden chiral PT at temperatures of order $100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ ($1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ to $10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$) may be detected by LISA (DECIGO/BBO). In contrast, if the transition occurs at temperatures above roughly $10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ the corresponding gravitational wave signal is unlikely to be observed by any currently proposed experiment. When looking at these findings, it has to be kept in mind that very fast transitions are generally anticipated to come along with a shortened period of GW production via sound waves and an appropriately suppressed spectrum. However, there are to date no results from dedicated numerical simulations to incorporate such a suppression in a quantitatively robust way, and thus it is hard to make reliable predictions.

Let us further remark that although all discussed low-energy effective models tend to give qualitatively similar results, the predicted gravitational wave spectra may vary considerably. This highlights the necessity of true first-principle calculations, like lattice simulations, for the purpose of reaching also quantitative precision.

While the present work focused on the the case of a new $\text{SU(}3\text{)}$ gauge sector with three massless fermion flavors, our analysis can, in principle, likewise be applied to different combinations of $(n_{c},n_{f})$. Those are then, however, not automatically guaranteed to give a first-order PT (see Ref. [16] for an overview). We also do not expect our results to change qualitatively if the hidden fermions acquire a (sufficiently small) current mass. Thus, our study and analogous studies can be employed to investigate the prospects for observing potential gravitational wave signals from a large variety of well-motivated new physics scenarios featuring strongly coupled hidden or dark sectors.

Acknowledgements.

We would like to thank Manfred Lindner, Jan M. Pawlowski, Dirk H. Rischke, Shinji Takeda, Hong Mao, Masatoshi Yamada, Christian Döring, Thomas Hugle, and Vedran Brdar for valuable discussions, as well as Carroll L. Wainwright for correspondence regarding the CosmoTransitions package. The work of J.K. is partially supported by the Grant-in-Aid for Scientific Research (C) from the Japan Society for Promotion of Science (Grant No.19K03844).

Appendix A Collection of integrals

In this appendix, we collect integrals used in various expressions throughout the paper. For instance, the one-loop meson propagators of the (P)NJL model as given in Appendix B can be brought into a compact form by defining

[TABLE]

where $M_{c}=M(v_{\sigma})$, with $M(\bar{\sigma})$ given by Eq. 13d, and $v_{\sigma}$ denotes the zero-temperature vacuum expectation value of the $\sigma$ field. All of the above integrals correspond to diagrams which arise from integrating out the fermion fields. Similarly, the expression for the pion decay constant in the (P)NJL model contains

[TABLE]

Note that the integrals in Eqs. 54 and 55 are divergent and must therefore be regularized. We employ a hard four-dimensional Euclidean cutoff $\Lambda$ to do so.

Next, the integrals relevant for determining the wave-function renormalization $Z_{\sigma}$ in the (P)NJL model (cf. Eq. 65) are given by

[TABLE]

Lastly, the bosonic (B) and fermionic (F) thermal integrals needed to compute the one-loop thermal contributions to the effective potentials in Section II are defined as

[TABLE]

Appendix B Details of the (P)NJL model

This appendix elaborates on certain aspects of the (P)NJL model and is meant to complement the corresponding discussions in Section II. Let us start with a slightly more detailed account of the self-consistent mean field approximation (MFA) [83, 84, 85, 81]. For convenience, we repeat here the Lagrangian for the three-flavor NJL model which is given by

[TABLE]

Using the MFA, the above Lagrangian can be written in terms of the fermion fields $q_{i}$ and auxiliary meson fields, cf. also Ref. [58]. This is done by going from the perturbative to the Bogoliubov-Valatin (BV) vacuum and using Wick’s theorem for operator products normal ordered with respect to the BV vacuum. Consequently, one can split up the Lagrangian into an interacting part $\mathcal{L}_{\text{int}}$ and a mean field term $\mathcal{L}_{\text{\tiny{MFA}}}$, such that $\mathcal{L}_{\text{\tiny{MFA}}}$ contains all terms which are at most quadratic in the fermion fields. One can then rewrite $\mathcal{L}_{\text{\tiny{MFA}}}$ in terms of $\Psi$ and its average $\langle\Psi\rangle$, which is defined as a sum over mesonic auxiliary fields

[TABLE]

or, equivalently,

[TABLE]

Working along the steps outlined above one can thus write the mean field Lagrangian in terms of the effective meson fields $\sigma$, $a_{a}$, $\pi_{a}$, ${\eta^{\prime}}$ and the fermion field $q$. The NJL Lagrangian in the MFA is given by

[TABLE]

where repeated flavor indices $a$, running from 1 to $n_{f}^{2}-1=8$, are summed over and we have defined $\pi:=2\pi_{a}T_{a}$ and $a:=2a_{a}T_{a}$, which are matrices in flavor space. The effective fermion mass $M$ was already introduced in Eq. 13d and is given by

[TABLE]

The tree-level potential is

[TABLE]

Next, we briefly discuss the zero-temperature one-loop meson propagators, whose roots define the meson’s effective pole masses. The propagators can be determined based on the mean field Lagrangian of Eq. 61 by calculating all one-loop 1PI diagrams with two external $\phi$ lines, where $\phi$ stands for one of the meson fields. The correlator $\mathrm{i}\Gamma_{\phi\phi}$ is then given by the sum of all of these diagrams. To find the relevant diagrams one has to keep in mind that only the fermions can run in the loop, because the mesons are represented by auxiliary fields in the NJL model and are thus non-propagating at tree-level. Furthermore, note that at zero temperature $\sigma$ has a non-zero vev, $v_{\sigma}$, which implies that diagrams with $v_{\sigma}$ as an external source also contribute to the meson propagators. Taking all these considerations into account, the one-loop meson propagators are computed to be

[TABLE]

Explicit formulas for $I_{V}$, $I_{S}$ and $I_{P}$ can be found in Appendix A. Equivalent results are given, for example, in Ref. [121], where the meson propagators were determined for non-zero quark masses. Our results agree with theirs when taking the chiral limit of the latter.

Finally, we provide the full expression for the $\sigma$ field’s wave-function renormalization $Z_{\sigma}$ as obtained within the MFA from its definition in Eq. 14:

[TABLE]

with $r\equiv r(\sigma)=|M(\sigma)|/T$ and

[TABLE]

where $M\equiv M(\sigma)$ is again the effective quark mass from Eq. 62 and $\Lambda$ denotes the four-dimensional hard momentum cutoff used to regularize the divergent vacuum parts of the occurring loop integrals. Expressions for the thermal parts $\ell_{i}$ can be found in Eq. 56 of Appendix A.

Appendix C Basics of the CJT formalism

The results of the linear sigma model presented in this work’s main part were derived within the composite operator formalism due to Cornwall, Jackiw and Tomboulis (CJT) [94]. Their formalism extends the approach to QFT based on the conventional effective action $\Gamma[\bar{\phi}]$ by employing a generalized functional $\Gamma[\bar{\phi},\bar{G}]$, the so-called two-particle irreducible (2PI) action. Here, $\bar{\phi}$ and $\bar{G}$ are expectation values121212To be more precise, the mentioned quantities are expectation values in the presence of external sources. of some quantum field $\Phi$ and of the corresponding two-point function, respectively. While the original CJT paper [94] only considered the vacuum case $T=0$, the composite operator method was subsequently extended to finite-temperature field theory problems in Refs. [102, 122], see also [123]. It was then employed to calculate the chiral phase transition in different versions of the linear sigma model e.g. in Refs. [97, 98, 99, 100].131313Furthermore, a variant of the CJT formalism, the so-called symmetry-improved CJT method [124], was applied to the $\text{O(}4\text{)}$ linear sigma model in Ref. [125]. In the present appendix, we will only give a brief overview of the CJT formalism’s very basics. For more details we refer the interested reader to the given references.

As mentioned above, the 2PI effective action is a functional of the generally spacetime-dependent expectation values $\bar{\phi}(x)$ and $\bar{G}(x,y)$, with a given theory’s physical solutions satisfying the following stationarity conditions:

[TABLE]

In situations with translational invariance, however, $\bar{\phi}$ becomes constant and $\bar{G}$ can be taken as a function of $(x-y)$ only. All relevant information is then encoded in the 2PI effective potential $V_{\text{\tiny{CJT}}}^{\text{\tiny{}}}[\bar{\phi},\bar{G}]$. In the chiral limit of the $\text{U(}3\text{)}\times\text{U(}3\text{)}$ linear sigma model (LSM) discussed in Section II.3, where mixing between the different meson species is absent, the aforementioned potential can be written as [99, 100]

[TABLE]

Here, the sum runs over all meson species and the corresponding multiplicities $g_{i}$ are given by $g_{\sigma}=g_{\eta^{\prime}}=1$ and $g_{\pi}=g_{a}=8$ for $n_{f}=3$. Finite-temperature momentum space integration is abbreviated as

[TABLE]

with $\omega_{n}=2\pi nT$ being the bosonic Matsubara frequencies. Furthermore, $V_{0}$ is the model’s tree-level potential given in Eq. 27a, while the $\Delta_{i}$ represent the tree-level meson propagators, i.e.

[TABLE]

with the field-dependent tree-level meson masses $m_{i}$ of Eq. 28. Finally, the quantity $V_{2}$ is computed as the sum of all two-particle irreducible vacuum graphs with each line corresponding to a full propagator $\bar{G}$, as well as generally $\bar{\sigma}$-dependent vertices.

Importantly, the conventional effective potential $V_{\text{eff}}$ as a function of the classical field $\bar{\sigma}$ is then obtained from Eq. 67 as

[TABLE]

where the full, field-dependent propagators $\hat{G}_{i}$ are determined via appropriate gap equations, which derive from the second stationarity condition in Eq. 66, and read

[TABLE]

Using the CJT effective potential from Eq. 67, the above gap equations can be equivalently written as

[TABLE]

being the scalar self-energy. Note that the well-known form of the one-loop one-particle irreducible effective potential $V_{\text{eff}}^{(1)}$ will be exactly reproduced from Eq. 70, if $V_{2}\equiv 0$ is assumed. Consequently, any non-trivial improvement over $V_{\text{eff}}^{(1)}$ requires a finite $V_{2}$. On the other hand, an exact determination of $V_{2}$ is in general not feasible since it involves the calculation of infinitely many diagrams. In practice, one must therefore resort to truncations of $V_{2}$.

Following Refs. [97, 98, 99, 100], we use the so-called Hartree-Fock approximation [95, 96] throughout the present work. At two-loop level, this corresponds to only retaining contributions from double-bubble graphs (cf. Fig. 9a) in computing $V_{2}$ [94]. In contrast, sunset diagrams as the one displayed in Fig. 9b are consistently ignored. According to its definition in Eq. 72 the scalar self-energy thus exclusively obtains contributions from graphs of the tadpole topology shown in Fig. 9c.

In order to solve the system of gap equations in Eq. 72 it is now generally convenient to make an ansatz for the full field-dependent propagator, namely

[TABLE]

In the Hartree-Fock approximation, the effective meson masses $M_{i}$ are the tree-level masses $m_{i}$ dressed by the previously mentioned tadpole contributions and can therefore be assumed to be momentum independent, i.e. $M_{i}\equiv M_{i}(\bar{\sigma})$, see e.g. [102]. Eq. 72 consequently reduces from a set of integral equations in $\hat{G}_{i}(k;\bar{\sigma})$ to a system of fixed-point equations in the effective masses $M_{i}(\bar{\sigma})$. Explicit calculations eventually yield [99, 100]

[TABLE]

where we have defined $R_{i}\equiv R_{i}(\bar{\sigma}):=M_{i}(\bar{\sigma})/T$, while the field-dependent tree-level meson masses $m_{i}(\bar{\sigma})$ can be found in Eq. 28. Note that, following Refs. [99, 100, 21], we consistently ignore vacuum contributions in evaluating loop integrals throughout all LSM-related calculations. The remaining thermal parts, $J_{\text{\tiny{B}}}$ and $I_{\text{\tiny{B}}}$, of the relevant functions can be found in Eq. 57.

For a given field value and at a given temperature, the system in Eq. 74 can be solved for the four effective masses $M_{i}$, which, in turn, define the full, field-dependent meson propagators according to the ansatz in Eq. 73. The sought-after finite-temperature effective potential can then be computed via Eq. 70, which eventually leads to the final expression quoted in Eq. 26.

Appendix D Model parameters

In Table 3, we collect our specific choices for the free parameters of the (P)NJL and linear sigma models, which correspond to the benchmark mass spectra listed in Table 1 and used throughout the paper. How to calculate the meson masses as well as the pion decay constant for a given set of parameters is explained in the main text, namely around Eq. 15 for the (P)NJL model and at the very end of Section II.3 for the linear sigma model.

Note that the scaled-up benchmark points discussed in Section V are obtained by rescaling the model parameters of Table 3 with $\xi$ according to their mass dimension, e.g. $G\to\xi^{-2}G$ etc., where the dimensionless factor $\xi$ was introduced in Eq. 51.