Gravitational Waves from First-Order Phase Transition in a Simple Axion-Like Particle Model

TL;DR

This paper explores how a simple axion-like particle model with a scalar field coupled to the Higgs can produce observable gravitational waves during a first-order phase transition, offering a new way to probe ALPs.

Contribution

It demonstrates that certain parameter regions in an ALP model can generate detectable gravitational waves, expanding the methods to explore ALP properties beyond traditional experiments.

Findings

Future gravitational-wave detectors can probe ALP energy scales from 10^3 to 10^8 GeV.

The phase transition can produce observable stochastic gravitational waves.

Detection prospects are complementary to existing ALP searches.

Abstract

We consider a gauge-singlet complex scalar field $\Phi$ with a global $U(1)$ symmetry that is spontaneously broken at some high energy scale $f_a$. As a result, the angular part of the $\Phi$-field becomes an axion-like particle (ALP). We show that if the $\Phi$-field has a non-zero coupling $\kappa$ to the Standard Model Higgs boson, there exists a certain region in the $\left(f_a, \kappa\right)$ parameter space where the global $U(1)$ symmetry-breaking induces a strongly first order phase transition, thereby producing stochastic gravitational waves that are potentially observable in current and future gravitational-wave detectors. In particular, we find that future gravitational-wave experiments such as TianQin, BBO and Cosmic Explorer could probe a broad range of the energy scale $10^3 \, {\rm GeV} \lesssim f_a \lesssim 10^{8} \, {\rm GeV}$, independent of the ALP mass. Since all the ALP couplings to the Standard Model particles are proportional to inverse powers of the energy scale $f_a$ (up to model-dependent ${\cal O}(1)$ coefficients), the gravitational-wave detection prospects are largely complementary to the current laboratory, astrophysical and cosmological probes of the ALP scenarios.

Figures

Equations

$$\Phi (x)=\frac{1}{\sqrt{2}}\left[f_a+\varphi (x)\right]e^{ia(x)/f_a}.$$

$${\mathcal{V}}_0=-{\mu }^2|H{|}^2+\lambda |H{|}^4+\kappa |\Phi {|}^2|H{|}^2+{\lambda }_a{\left(|\Phi {|}^2-\frac{1}{2}{f}_a^2\right)}^2.$$

$${\mathcal{V}}_0$$

$${\mathcal{V}}_0(\varphi ,h)=\frac{{\lambda }_a}{4}{\left({\varphi }^2-f_a^2\right)}^2+\frac{\kappa }{4}{\varphi }^2{h}^2-\frac{{\mu }^2}{2}{h}^2+\frac{\lambda }{4}{h}^4.$$

$${\mu }^2\le \frac{2\lambda {\lambda }_a}{\kappa }{f}_a^2.$$

$${\mu }^2=\frac{2\lambda {\lambda }_a}{\kappa }{f}_a^2.$$

$$\mathcal{V}(\varphi ,T)={\mathcal{V}}_0(\varphi )+{\mathcal{V}}_{\mathrm{CW}}(\varphi )+{\mathcal{V}}_T(\varphi ,T),$$

$${\mathcal{V}}_{\mathrm{CW}}\left(\varphi \right)=\sum _i(-1{)}^F{n}_i\frac{m_i^4\left(\varphi \right)}{64{\pi }^2}\left[\log \frac{m_i^2\left(\varphi \right)}{{\Lambda }^2}-C_i\right].$$

$${\mathcal{V}}_T\left(\varphi ,T\right)=\sum _i{\left(-1\right)}^F{n}_i\frac{T^4}{2{\pi }^2}{J}_{B/F}\left(\frac{m_i^2\left(\varphi \right)}{T^2}\right),$$

$$J_{B/F}\left(y^2\right)={\int }_0^{\infty }\mathrm{d}x{x}^2\log \left[1\mp \exp \left(-\sqrt{x^2+y^2}\right)\right].$$

$${\Pi }_h\left(T\right)={\Pi }_{G_{0,\pm }}\left(T\right)$$

$${\Pi }_{\varphi }\left(T\right)$$

$$\Gamma \left(T\right)$$

$$S_4=\int {\mathrm{d}}^{4}x\left[\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}t}\right)}^2+\frac{1}{2}{\left(\nabla \varphi \right)}^2+\mathcal{V}\left(\varphi ,T\right)\right],$$

$$\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{r}^2}+\frac{3}{r}\frac{\mathrm{d}\varphi }{\mathrm{d}r}={\mathcal{V}}^{\prime }\left(\varphi ,T\right).$$

$$S_3=\int {\mathrm{d}}^{3}x\left[{\left(\nabla \varphi \right)}^2+\mathcal{V}\left(\varphi ,T\right)\right],$$

$$\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{r}^2}+\frac{2}{r}\frac{\mathrm{d}\varphi }{\mathrm{d}r}={\mathcal{V}}^{\prime }\left(\varphi ,T\right).$$

$$h^2{\Omega }_{\text{GW}}\simeq h^2{\Omega }_{\varphi }+h^2{\Omega }_{\text{SW}}+h^2{\Omega }_{\text{MHD}}.$$

$$h^2{\Omega }_{\varphi }\left(f\right)$$

$$S_{\text{env}}\left(f\right)=\frac{3.8{\left(f/f_{\text{env}}\right)}^{2.8}}{1+2.8{\left(f/f_{\text{env}}\right)}^{3.8}}.$$

$$f_{\text{env}}=\left(\frac{f_{\ast }}{\beta }\right)\left(\frac{\beta }{H_{\ast }}\right)h_{\ast }.$$

$$h_{\ast }=16.5\times 10^{-3}\text{mHz}\left(\frac{T_{\ast }}{100\text{ GeV}}\right){\left(\frac{g_{\ast }}{100}\right)}^{1/6}.$$

$$h^2{\Omega }_{\text{SW}}\left(f\right)$$

$$S_{\text{SW}}\left(f\right)={\left(\frac{f}{f_{\text{SW}}}\right)}^3{\left(\frac{7}{4+3{\left(f/f_{\text{SW}}\right)}^2}\right)}^{7/2},$$

$$f_{\text{SW}}=1.008\times \frac{2}{\sqrt{3}{v}_w}\left(\frac{\beta }{H_{\ast }}\right)h_{\ast }.$$

$$h^2{\Omega }_{\text{MHD}}\left(f\right)$$

$$S_{\text{MHD}}\left(f\right)=\frac{{\left(f/f_{\text{MHD}}\right)}^3}{{\left[1+\left(f/f_{\text{MHD}}\right)\right]}^{11/3}\left(1+8\pi f/h_{\ast }\right)},$$

$$f_{\text{MHD}}=0.935\left(\frac{3.5}{2v_w}\right)\left(\frac{\beta }{H_{\ast }}\right)h_{\ast }.$$

$$\beta \equiv -{\frac{\mathrm{d}{S}_E}{\mathrm{d}t}|}_{t=t_{\ast }}=T{H}_{\ast }{\frac{\mathrm{d}{S}_E}{\mathrm{d}T}|}_{T=T_{\ast }},$$

$${\kappa }_{\varphi }=\frac{\alpha -{\alpha }_{\infty }}{\alpha },{\kappa }_v=\frac{{\alpha }_{\infty }}{\alpha }{\kappa }_{\infty },{\kappa }_{\text{MHD}}=ϵ{\kappa }_v,$$

$${\kappa }_{\infty }$$

Full text

1 Introduction

Axion-like particles (ALPs) are light gauge-singlet pseudoscalar bosons that couple weakly to the Standard Model (SM) and generically appear as the pseudo-Nambu-Goldstone boson (pNGB) in theories with a spontaneously broken global $U(1)$ symmetry. ALPs could solve some of the open questions in the SM, such as the strong CP problem via the Peccei-Quinn mechanism [1] and the hierarchy problem via the relaxion mechanism [2]. They could also play an important cosmological role in inflation [3, 4, 5], dark matter (DM) [6, 7, 8], dark energy [9, 10, 11, 12], and baryogenesis [13, 14]. Furthermore, there are recent proposals involving axions to simultaneously address several open issues of the SM in one stroke [15, 16, 17, 18, 19, 20].

A common characteristic among ALPs is that their coupling to SM particles is suppressed by inverse powers of the $U(1)$ symmetry breaking energy scale $f_{a}$. This energy scale can be identified as the vacuum expectation value (VEV) of a SM-singlet complex scalar field $\Phi$, i.e. $\langle\Phi\rangle=f_{a}/\sqrt{2}$, which is assumed to be much larger than the electroweak scale $v_{\rm ew}\simeq 246.2$ GeV to evade current experimental limits [21, 22]. The ALP field $a$ then arises as the massless excitation of the angular part of the $\Phi$-field:

[TABLE]

The particle excitation of the modulus $\phi$ of the $\Phi$-field gets a large mass $m_{\phi}\sim f_{a}$, while the angular part $a$ becomes a pNGB that acquires a much smaller mass $m_{a}$ from explicit low energy $U(1)$-breaking effects. Thus, for the low-energy phenomenology of ALPs, the modulus part $\phi$ can be safely integrated out, and the only experimentally relevant parameters are $m_{a}$ and $f_{a}$.

In this paper, we show that the dynamics of the modulus $\phi$-field around the $f_{a}$ scale can provide complementary constraints on the ALP scenario. In particular, if the parent $\Phi$-field has a non-zero coupling to the SM Higgs doublet $H$, the $U(1)$ symmetry breaking at the $f_{a}$-scale could induce a strongly first-order phase transition (FOPT) [23], giving rise to stochastic gravitational waves (GWs) that are potentially observable in current and future GW detectors (see e.g. [24, 25] for a review on GWs from a FOPT). We find that GW signals of strength up to $h^{2}\Omega_{\rm GW}\sim 10^{-12}$ could be generated, where $\Omega_{\rm GW}$ is the fraction of the total energy density of the universe in the form of GWs today and $h=0.674\pm 0.005$ is the current value of the Hubble parameter in units of 100 km s*-1* Mpc*-1* [26]. Future GW observatories like TianQin [27], Taiji [28], LISA [29, 30], ALIA [31], MAGIS [32], DECIGO [33], BBO [34], Cosmic Explorer (CE) [35] and Einstein Telescope (ET) [36] can probe a broad range $10^{3}\text{ GeV}\lesssim f_{a}\lesssim 10^{8}\text{ GeV}$, independent of the ALP mass. It turns out that the aLIGO [37]  and aLIGO+ [38]  can not probe any of the allowed parameter space for the benchmark configurations considered in this paper (see Fig. 4).

The heavy modulus $\phi$ decouples at low energies, and we are left with the ALP $a$, which has only derivative couplings to the SM particles. These are generated via effective higher-dimensional operators [39], and are proportional to inverse powers of $f_{a}$, up to model-dependent ${\cal O}(1)$ coefficients. The effective ALP couplings to photons, electrons and nucleons are strongly constrained by a number of laboratory, astrophysical and cosmological observables [22]. However, current and future low-energy constraints depend on the ALP mass $m_{a}$, while we find that the GW prospects in the $(m_{a},f_{a})$ plane are largely complementary. For instance, if a stochastic GW signal was found with the frequency dependence predicted by the FOPT111The frequency dependence in this scenario is generically different and can be distinguished from other stochastic GW sources, like inflation [40, 41, 42, 43, 44, 45] or unresolved binary black hole mergers [46]., this would point to a limited range of $f_{a}$ in a given ALP model, which might lead to a positive signal in some of the future laboratory and/or astrophysical searches of ALPs. On the other hand, if we fix the ALP mass $m_{a}$, then current ALP constraints exclude certain ranges of $f_{a}$. If a GW signal is found in the frequency range corresponding to the excluded range of $f_{a}$, then the underlying simple ALP model has to be extended to account for the GW signal.

The rest of the paper is organized as follows: in Section 2 we provide the details of the ALP model, and compute the one-loop effective scalar potential at both zero and finite temperatures. In Section 3 we calculate the GW emission from a strong FOPT at the scale $f_{a}$, including bubble collision, sound wave (SW) and magnetohydrodynamic (MHD) turbulence contributions. The complementary reach of laboratory, astrophysical and cosmological observations is presented in Section 4. The constraints from future precision Higgs data on the ALP model are discussed in Section 5. We summarize and conclude in Section 6. The method used to obtain the power-law integrated sensitivity curves for future GW experiments is described in Appendix A.

2 Scalar Potential in the ALP model

2.1 Tree-level potential

The coupling between the SM Higgs doublet $H$ and the complex field $\Phi$ is described by the tree-level potential (see also Refs. [23, 47])

[TABLE]

The SM Higgs doublet can be parameterized as $H=\left(G_{+},\left(h+iG_{0}\right)/\sqrt{2}\right)$, with $h$ the SM Higgs and $G_{0}$, $G_{+}$ the Goldstone bosons that become the longitudinal components of the $Z$ and $W^{+}$ bosons, respectively. The complex field $\Phi$ can be expressed in the form given by Eq. (1.1). At an energy scale around $f_{a}$, a phase transition (PT) occurs which breaks the global $U(1)$ symmetry. The field $\Phi$ gets a VEV $\langle\Phi\rangle=f_{a}/\sqrt{2}$, and the associated pNGB $a$ is identified as the physical ALP. Depending on the parameters $f_{a}$, $\kappa$ and $\lambda_{a}$, the PT may be strongly first order, in which case it would generate a spectrum of GWs that could be detected in current or future experiments [48], as detailed in Section 3. We note that at this stage the ALP $a$ is neither involved in the scalar potential given by Eq. (2.1), nor in the GW emission from the high-energy scale PT. The low-energy effective couplings of $a$ to SM particles will be discussed in Section 4.1.

In terms of the real scalar field components, the tree-level potential in Eq. (2.1) can be re-written as:

[TABLE]

Setting the field values of the Goldstone modes to zero, we have

[TABLE]

Here $f_{a}$, $\mu$, $\kappa$ and $\lambda_{a}$ are taken as free parameters. By examining the tree-level potential in Eq. (2.3) we have found that a FOPT occurs along the $\phi$ direction, while maintaining the VEV of $h$ equal to zero at the same time, if the parameters satisfy the inequality:

[TABLE]

While this is not the only possible choice that results in a FOPT along the $\phi$ direction, for our present purpose, we will set the $\mu$-parameter to the value that saturates the inequality above, i.e.,

[TABLE]

We have verified numerically that, with this choice of $\mu^{2}$, the VEV of $h$ remains zero up to one-loop level during the phase transition, if there is one. Therefore, in the following analysis, we will ignore the dependence of the effective potential on the $h$-field, and consider $\phi$ as the only dynamical field.

It should be noted that with the specific choice in Eq. (2.5), the scalar $\phi$ will contribute radiatively to the SM Higgs mass, making the latter unacceptably large at the electroweak scale. Nevertheless, this contribution can be cancelled out, e.g. by introducing vector-like fermions at the $\mu$-scale that keep the SM Higgs boson mass at the observed 125 GeV [49]. In the parameter space of interest in this paper, $\lambda_{a}\simeq 10^{-3}$ and $\kappa\simeq 1$ (see Section 3), we have $\mu\sim 10^{-2}f_{a}$. Therefore, the effect of these extra vector-like fermions on the renormalization group (RG) running of the SM couplings, as required for the calculation of the critical temperature (see Section 3.3), is expected to be small and will not be considered here. Since we are focusing on a generic ALP scenario below the $f_{a}$ scale, we will defer a detailed study of ultraviolet-completion involving additional heavy fields to a future work.

2.2 Effective finite-temperature potential

At finite temperature $T\neq 0$, the effective one-loop potential of the scalar fields is [50, 51, 52, 53]:

[TABLE]

where $\mathcal{V}_{\rm CW}$ is the Coleman-Weinberg (CW) potential [54] that contains all the one-loop corrections at zero temperature with vanishing external momenta, and $\mathcal{V}_{T}$ describes the finite-temperature corrections. Working in the Landau gauge to avoid ghost-compensating terms, the CW potential reads:

[TABLE]

The sum runs over all the particles that couple to the $\phi$ field (notice that massless particles do not contribute). In Eq. (2.7), $F=1$ for fermions and [math] for bosons; $n_{i}$ is the number of degrees of freedom of each particle; $C_{i}=3/2$ for scalars and fermions, and $5/6$ for gauge bosons; and $\Lambda$ is the renormalization scale, which will be set to $f_{a}$ throughout this paper.

The finite-temperature corrections are given by:

[TABLE]

where the thermal functions are:

[TABLE]

Here, the minus sign “$-$” is for bosons and the positive sign “$+$” for fermions. We also need to include the resummed daisy corrections, that add a temperature-dependent term $\Pi_{i}(T)$ to the field-dependent mass $m_{i}^{2}$ [53]. To leading order, we have in the ALP model:

[TABLE]

Effectively, the mass terms $m_{i}^{2}$ in Eqs. (2.7) and (2.8) get replaced by $m_{i}^{2}+\Pi_{i}(T)$.

The effective potential in Eq. (2.6) could become complex due to $m_{i}^{2}$ being negative. This is related to particle decay and does not affect the computation of the dynamics of PT (see e.g. Refs. [55, 56, 57] for more details). In the numerical calculations in Section 3 we will always take the real part of the effective potential.

3 First-order phase transition

We fix the decay constant $f_{a}$ at a specific value and scan in $(\kappa,\lambda_{a})$ to find the region in parameter space where a FOPT can take place. In particular, we evaluate the effective potential for given values of $(\kappa,\lambda_{a})$ to look for regions where there is a valid critical temperature $T_{c}$. Here, $T_{c}$ is defined as the temperature at which the two local minima of the effective potential are degenerate.222In our analysis, we have not considered the potential gauge-invariance [58] and RG improvement [47] effects on the calculation of $T_{c}$ and the resultant GW prospects, which are beyond the scope of this paper. We use the package CosmoTransitions [59] for the numerical work and the results are given in the following subsections. We show one example of the effective potential varying temperature $T$ in Fig. 1. In this case, a valid $T_{c}$ can be found, and tunneling can happen for $T<T_{c}$ to give rise to a FOPT.

3.1 Bounce solution

The decay rate of the false vacuum is [60, 61, 62]:

[TABLE]

where the first term corresponds to thermally induced decays and the second term is the quantum-tunneling rate. In Eq. (3.1) $S_{3}$ and $S_{4}$ are the three- and four-dimensional Euclidean actions for the $O(3)$ and $O(4)$-symmetric tunnelling (“bounce”) solutions respectively, and $R_{0}$ is the size of the bubble. For quantum tunneling,

[TABLE]

where $\phi(r)$ is the solution of the $O(4)$-symmetric instanton (with $r=\sqrt{\boldsymbol{r}^{2}+t^{2}}$):

[TABLE]

For thermally-induced decay,

[TABLE]

where $\phi(r)$ is the $O(3)$-symmetric solution of

[TABLE]

3.2 Gravitational waves

The GW signal from a FOPT consists of three main components: the scalar field contribution during the collision of bubble walls [63, 64, 65, 66, 67, 68], the sound wave in the plasma after bubble collisions [69, 70, 71, 72], and the MHD turbulence in the plasma after bubble collisions [73, 74, 75, 76, 77]. Assuming the three components can be linearly superposed, the total strength of GWs produced reads

[TABLE]

Note that the global $U(1)$ symmetry breaking could also generate cosmic strings, which then annihilate to produce GWs [78, 79]. However, this effect turns out to be subdominant for the energy scales under consideration here.

The envelope approximation is often used to calculate the GWs from the scalar $\phi$ contribution, and numerical simulations tracing the envelope of thin-walled bubbles reveal that [65, 80]333For an analytic derivation of GW production in the thin-wall and envelope approximation, see e.g. Refs. [81, 82].

[TABLE]

where $f$ is the frequency; $g_{\ast}$ is the number of relativistic degrees of freedom in the plasma at the temperature $T_{\ast}$ when the GWs are generated; $H_{\ast}$ is the Hubble parameter at $T_{\ast}$; $v_{w}$ is the bubble wall velocity in the rest frame of the fluid; $\alpha\equiv{\rho_{\rm vac}}/{\rho^{\ast}_{\rm rad}}$ is the ratio of the vacuum energy density $\rho_{\rm vac}$ released in the PT to that of the radiation bath $\rho^{\ast}_{\rm rad}=g_{\ast}\pi^{2}T_{\ast}^{4}/30$; $\beta/H_{*}$ measures the rate of the PT; $\kappa_{\phi}$ measures the fraction of vacuum energy that is converted to gradient energy of the $\phi$ field; and $S_{\text{env}}(f)$ parameterizes the spectral shape of the GW radiation,

[TABLE]

The peak frequency $f_{\rm env}$ of the $\phi$ contribution to the spectrum is determined by $\beta$ and by the peak frequency $f_{\ast}=0.62\beta/(1.8-0.1v_{w}+v_{w}^{2})$ [65] at the time of GW production,

[TABLE]

Assuming the Universe is radiation-dominated after the PT and has expanded adiabatically ever since, the inverse Hubble time $h_{\ast}$ at GW production, red-shifted to today, is

[TABLE]

The SW contribution is given by [72]:

[TABLE]

where $\kappa_{v}$ is the fraction of vacuum energy that is converted to bulk motion of the fluid, and the spectral shape

[TABLE]

with the peak frequency

[TABLE]

The MHD turbulence contribution is given by [77, 83]:

[TABLE]

where $\kappa_{\rm MHD}$ is the fraction of vacuum energy that is transformed into MHD turbulence, and the spectral shape can be found analytically:

[TABLE]

where the peak frequency measured today is

[TABLE]

In most of the parameter space of interest, the phase transition occurs in the ‘runaway bubbles in the plasma’ regime, where the $\Omega_{\phi}$ contribution cannot be neglected. In principle, friction from the plasma could stop the bubble wall at some terminal velocity [84], thus rendering the plasma-related GW contributions and their associated uncertainties more important [85, 86]. The friction term, however, turns out to be unimportant in our case, because the $\phi$ field only couples to the Higgs doublet through a scalar quartic interaction, and does not directly couple to the gauge fields. Therefore, the friction from the plasma does not grow with energy [84], thus preserving the runaway behavior. In this limit, the uncertainties in the SW and MHD contributions to $\Omega_{\rm GW}$ due to nonlinearities developing in the plasma do not significantly affect our results.

To calculate the GW signal $\Omega_{\rm GW}$ described above, we need to know the following quantities:

• •

The ratio $\alpha$ of vacuum energy density released in the PT to that of the radiation bath.

• •

The rate of the PT, $\beta/H_{*}$. The smaller $\beta/H_{*}$, the stronger the PT. From the bubble nucleation rate $\Gamma(t)=A(t)e^{-S_{E}(t)}$ [87], with $A(t)$ the amplitude and $S_{E}$ the Euclidean action of a critical bubble, we have:

[TABLE]

where we have assumed the nucleation temperature $T_{n}\simeq T_{\ast}$ (or equivalently $t_{n}\simeq t_{\ast}$ with $t_{n}$ and $t_{\ast}$ respectively the time for bubble nucleation and GW production).

• •

The latent heat fractions $\kappa$ for each of the three processes. For the case of runaway bubbles in a plasma, we have

[TABLE]

where $\epsilon$ is the turbulent fraction of bulk motion, which is found to be at most ($5\%-10\%$) [72]. To be concrete, we choose $\epsilon=0.1$ in this paper. We also have:

[TABLE]

In Eq. (3.17), the sum runs over all particles $i$ that are light in the initial phase and heavy in the final phase; $\Delta m_{i}^{2}$ is the squared mass difference in the two phases; and $c_{i}=n_{i}\,(n_{i}/2)$ for bosons (fermions) with $n_{i}$ the number of degrees of freedom of the particle [88].

• •

The bubble wall velocity $v_{w}$ in the rest frame of the fluid away from the bubble. A conservative estimate for $v_{w}$ is given by [89]:

[TABLE]

• •

The number of relativistic degrees of freedom $g_{*}$ at the time of the PT, which is taken to be the SM contribution of 106.75 plus an additional 1 from the ALP.

3.3 Critical temperature

For the calculation of the GW production associated to the PT, we take $f_{a}$, $\kappa$ and $\lambda_{a}$ as the only free parameters in the scalar potential (2.1). Below the scale $f_{a}$, the scalar sector only contains the SM Higgs and the superlight ALP $a$. The value of $\lambda_{a}$, as well as the other relevant coupling constants appearing in Eqs. (2.10) and (2.11), at a high energy scale $\Lambda<f_{a}$ can be obtained by running the SM RG equations up to the scale $\Lambda$ [90]. We will consider values $f_{a}\leq 10^{8}$ GeV, as the SM vacuum becomes unstable for $\Lambda\gtrsim 10^{8}$ GeV [91] with the current best-fit top quark mass $m_{t}=173.0$ GeV [92]. Similarly, we take $f_{a}\geq 10^{3}$ GeV, because for $f_{a}$ comparable to (or smaller than) the Higgs mass, the LHC Higgs data impose stringent constraints on the coupling $\kappa$ (see Section 5).

Fixing $f_{a}=10^{6}$ GeV, we scan the two parameters $\kappa$ and $\lambda_{a}$. The critical temperature is shown in the left panel of Fig. 2, in units of $f_{a}$. It can be seen that $T_{c}$ gets larger for larger $\kappa$ or larger $\lambda_{a}$, and in the parameter region we focus on, we have $T_{c}\lesssim f_{a}$.

3.4 Bubble nucleation

For the region in parameter space where there is a valid $T_{c}$, bubble nucleation might occur when $T<T_{c}$, i.e., when the two local minima become non-degenerate. The nucleation temperature $T_{n}$ is estimated by the condition:

[TABLE]

where the Hubble constant is given by [85]:

[TABLE]

If a valid solution for $T_{n}$ from Eq. (3.19) can be found, this indicates that the nucleation process will happen, and that the majority of the GW signal is produced at $T_{*}\simeq T_{n}$. $\beta/H_{*}$ is obtained from Eq. (3.14), and the result for the nucleation temperature $T_{n}$ is shown in the right panel of Fig. 2, for the specific value $f_{a}=10^{6}$ GeV. It is clear from Fig. 2 that the region with a viable nucleation temperature $T_{n}$ is a subset of the FOPT region, and that $T_{n}\lesssim T_{c}$.

The two important parameters for computing the GW signal, $\beta/H_{*}$ and $\alpha$, are then evaluated at $T_{n}$. We show the results of these two parameters in the parameter space $(\kappa,\lambda_{a})$ for the case of $f_{a}=10^{6}$ GeV in Fig. 3. $\beta/H_{*}$ can reach values as small as $\sim 10^{2}$; while $\alpha\lesssim 0.05$ for most of the parameter space, but it can reach values of order 1. The region with large $\alpha$ and relatively small $\beta/H_{*}$ is where a relatively large GW signal is expected. This is the region where $0.001\lesssim\lambda_{a}\lesssim 0.1$ and $\kappa\sim{\cal O}(1)$. Note that larger values of $\kappa$ would lead to a breakdown of the perturbation theory, as a Landau pole is developed below the $f_{a}$ scale. Moreover, the one-loop self-energy corrections to the $\phi$ field due to the Higgs and Goldstone modes will be large in this case. In order to avoid these theoretical issues, we will restrict ourselves to $\kappa<6$.

3.5 Detection prospects

We have assumed that the majority of the GW signal is produced at $T_{\ast}\simeq T_{n}$. Since the three GW contributions scale as inverse powers of $\beta/H_{\ast}$,

[TABLE]

we expect that a relatively large GW signal can be generated in the small-$\beta/H_{*}$ region.

The three different components of GW signals from bubble wall collision in Eq. (3.2), SW in the plasma in Eq. (3.2), and MHD turbulence in Eq. (3.2) are added to obtain the total GW emission, given in Eq. (3.6), as a function of frequency $f$. Our numerical simulations reveal that the GW signal obtained in this model can be as large as $h^{2}\Omega_{\text{GW}}\approx 10^{-12}$ for configurations with $\kappa\approx 1$ and $\lambda_{a}\approx 10^{-3}$. For such configurations that can produce considerable GW signals, contributions from both scalar fields and sound waves are comparable, and cannot be neglected, as indicated by the shape of the $h^{2}\Omega_{\text{GW}}$ curves in Fig. 4. In most cases we have $h^{2}\Omega_{\text{SW}}\gtrsim h^{2}\Omega_{\phi}$ and the MHD contribution is much smaller than the other two contributions.

In Fig. 4, we show the detection prospects for the future GW experiments TianQin [27], Taiji [28], LISA [29, 30], ALIA [31], MAGIS [32], DECIGO [33], BBO [34], aLIGO [37], aLIGO+ [38], ET [36] and CE [35]. To see the dependence of the GW signal on the parameters $f_{a}$, $\kappa$ and $\lambda_{a}$, let us first fix $f_{a}=10^{6}$ GeV and $\kappa=1.0$ and vary the quartic coupling $\lambda_{a}$ from $0.001$ to $0.2$. The corresponding GW signal $h^{2}\Omega_{\rm GW}$ is shown in the upper panel of Fig. 4, as a function of the frequency $f$. It is obvious that within the $T_{*}$ region the configurations with smaller $\lambda_{a}$ tend to produce a larger GW signal with a relatively smaller peak frequency, which is preferred by the space-based experiments. Likewise, when we fix $f_{a}=10^{6}$ GeV and $\lambda_{a}=0.001$, the configurations with $\kappa\approx 1$ tend to produce the strongest GW signals at small peak frequencies, as shown in the lower left panel of Fig. 4. When $\kappa$ and $\lambda_{a}$ are fixed, e.g. $\kappa=1.00$ and $\lambda_{a}=0.001$, a larger $f_{a}$ tends to produce a slightly larger GW signal with a larger peak frequency, as seen in the lower right panel of Fig. 4.

The GW detection prospects in the two-dimensional plane of $\kappa$ and $\lambda_{a}$ are shown in Fig. 5 for the benchmark values $f_{a}=10^{3,4,5,6,7,8}$ GeV. For the sake of clarity, we show only the sensitivity regions for three selected GW experiments: TianQin, BBO and CE.

The current and future GW observations are largely complementary to each other. For instance, TianQin, Taiji, LISA  and ALIA are more sensitive to the GWs with a comparatively lower frequency and thus a smaller $f_{a}$; aLIGO, aLIGO+, ET and CE could probe higher frequency GWs, and thus larger $f_{a}$; while MAGIS, DECIGO and BBO are able to cover the frequency range in between. This is explicitly illustrated in Fig. 6, for the benchmark values $\kappa=1.0$ and $\lambda_{a}=0.001$. We use the power-law integrated sensitivity curves for future GW experiments, as described in Appendix A. Although the GW emission in the ALP model does not directly involve the ALP particle $a$, future GW observations could definitely probe a broad range of the decay constant $f_{a}$, which largely complements the low-energy, high-energy, astrophysical and cosmological constraints on the $f_{a}$ parameter, as detailed in Section 4. For related discussions on GW emission from the ALP field itself, see e.g. Refs. [93, 94, 95].

4 Comparison with other ALP constraints

As shown in Figs. 4, 5 of Section 3, current and future GW observations could probe a broad region of the parameter space in the ALP model. In particular, the scale $f_{a}$ in the range $(10^{3}-10^{8})$ GeV can be probed by future GW observatories, as summarized in Fig. 6. At low energies, all the ALP couplings to SM particles are inversely proportional to powers of the decay constant $f_{a}$ (see e.g. Eqs. (4.1) and (4.5)); thus, GW observations are largely complementary to the laboratory, astrophysical and cosmological constraints on the couplings of $a$ to SM particles. For the sake of simplicity, we will focus on the effective CP-conserving ALP couplings to photons ($g_{a\gamma\gamma}$), electrons ($g_{aee}$) and nucleons ($g_{aNN}$). In principle, the ALP could also couple to other SM particles like the muon, tau and other gauge bosons (gluons, $W$ and $Z$ boson), and we could even have CP-violating couplings to SM particles [39, 100, 96] (see e.g. Refs. [96, 22] for more details). In addition, the muon $g-2$ anomaly could be explained by ALP couplings to muons and photons [100, 97, 98, 99] or by flavor violating couplings to muons and taus [101].

4.1 Low-energy effective ALP couplings

Even though the ALP $a$ does not couple directly to the SM Higgs or the real scalar $\phi$ in the potential (2.1), low-energy couplings to SM particles can be induced at dimension-5 or higher. For instance, the effective couplings of $a$ to the SM photon and fermions $f$ can be written as

[TABLE]

Here, $F_{\mu\nu}$ is the electromagnetic field strength tensor and $\tilde{F}^{\mu\nu}$ its dual, $\alpha_{\rm EM}$ is the fine-structure constant, and $C_{a\gamma}$, $C_{af}$ are model-dependent coefficients. Generally speaking, these coefficients are of order one for the QCD axion. Setting the model-dependent coefficients $C_{i}$ to one for simplicity, we rewrite the couplings in Eq. (4.1) as

[TABLE]

where the effective couplings are related to the high scale $f_{a}$ via

[TABLE]

and $m_{f}$ is the corresponding fermion mass. We thus see that the GW limits on $f_{a}$ from Fig. 6 can be used to probe the effective couplings $g_{a\gamma\gamma}$, $g_{aee}$ and $g_{aNN}$.

4.2 Coupling to photons

Following Ref. [22], all the current constraints on the ALP couplings to photons $g_{a\gamma\gamma}$ are collected in the left panel of Fig. 7, while future laboratory and astrophysical prospects are shown in the right panel of Fig. 7. In both panels we also show the parameter space for DFSZ [102, 103] and KSVZ [104, 105] axions, indicated respectively by the yellow region and brown line. The various constraints are explained below:

• •

Given the coupling $g_{a\gamma\gamma}$, the ALP can decay into two photons in the early universe, with a rate depending largely on its mass $m_{a}$ and the magnitude of $g_{a\gamma\gamma}$. If ALPs decay before recombination, the photons produced in the decays would potentially distort the cosmic microwave background (CMB) spectrum and, at earlier times, they would also affect big bang nucleosynthesis (BBN) [106]. The monochromatic photon lines from axion/ALP decays are also constrained by the flux of extragalactic background light (EBL) and direct searches in X-rays and $\gamma$-rays [106]. Furthermore, the photons might also change the evolution of the hydrogen ionisation fraction, $x_{\rm ion}$ [106]. These limits are shown in greenish color in the left panel of Fig. 7.

• •

Assuming ALPs account for all the DM, the regions labelled as “telescopes” in the left panel of Fig. 7 have been excluded by direct decaying DM searches in galaxies [108, 107]. It is promising that future telescopes could probe couplings down to [109]

[TABLE]

Taking the distance to the ALP source to be $d\simeq 2\,{\rm kpc}$, the future sensitivity could reach $g_{a\gamma\gamma}\sim 10^{-13}\,{\rm GeV}^{-1}$. This is shown as the green line in the right panel of Fig. 7.

• •

As a result of the coupling $g_{a\gamma\gamma}$, ALPs can be produced and emitted copiously from dense stellar cores, thus affecting stellar evolution [110, 111, 112]. Large portions of the $\left(m_{a},g_{a\gamma\gamma}\right)$ parameter space have been excluded by measurements of the solar neutrino flux and helioseismology [113], the ratio of horizontal branch (HB) to red giants in globular clusters [114], and SN1987A neutrino data [115]. These limits are labelled respectively as “Sun”, “HB stars” and “SN” in the left panel of Fig. 7.

• •

In the presence of an electromagnetic field, the ALP can be converted to a photon through the $a\gamma\gamma$ coupling [116]. The axion helioscope experiments Brookhaven [117], SUMICO [118, 119, 120] and CAST [121, 122, 123, 124, 125, 126] aim to detect X-rays from $a-\gamma$ conversion in the Sun. The absence of a signal can be used to set the limits on $g_{a\gamma\gamma}$ labelled as “helioscopes” in the left panel of Fig. 7 . Future experiments such as TASTE [127] and IAXO [128, 129] could improve current constraints by over one order of magnitude, also shown by the solid red line in the right panel of Fig. 7.

• •

In a static magnetic field, ALP DM in the $\sim 10^{-6}$ eV mass range can be converted into a microwave photons [116]. Narrow regions around this range have been excluded by the ADMX experiment [130, 131, 132, 133]. They are labelled as “haloscopes” in the left panel of Fig. 7. Although we do not display these limits, let us mention that future stages of ADMX could probe a very narrow range around $m_{a}\sim 10\ \mu{\rm eV}$ [134], while the ABRACADABRA experiment might be sensitive to light ALPs with mass $10^{-14}\,{\rm eV}\lesssim m_{a}\lesssim 10^{-6}\,{\rm eV}$ and couplings down to $g_{a\gamma\gamma}\sim 10^{-19}\,{\rm GeV}^{-1}$ [135].

• •

Light-shining-through-wall (LSW) experiments provide the most stringent laboratory constraints on $g_{a\gamma\gamma}$ for a broad range of ALP mass $m_{a}$. In such experiments, ALPs can be produced from intense photon sources in the presence of magnetic fields and then convert back into photons. The LSW limits from BRFT [136], BMV [137], GammaV [138], LIPPS [139], ALPS [140], OSQAR [141] and CROWS [142] are collectively shown in the left panel of Fig. 7. The LSW limits could be further improved by up to four orders of magnitude by the experiments ALPS II [143], ALPS III [144] and STAX [145], as indicated by the dashed blue, solid blue and solid purple lines in the right panel of Fig. 7. There are also constraints from the polarization experiment PVLAS [146, 136, 147, 148] and from fifth force searches [149, 150, 151, 152, 153, 155, 154], which are however weaker and thus not shown in Fig. 7.

• •

In beam-dump experiments, ALPs can be produced off photons. The limits on $g_{a\gamma\gamma}$ from the experiments CHARM [156], E137 [157], E141 [158] and NuCal [159, 160] are comparatively weaker than those from the astrophysical observations above, excluding a region $g_{a\gamma\gamma}\gtrsim 10^{-7}\,{\rm GeV}^{-1}$ for ALP masses $m_{a}\sim({\rm MeV}-{\rm GeV})$, as shown in Fig. 7. The future experiment SHiP [161, 162] will extend the exclusion regions to higher $m_{a}$, but it will not push to smaller couplings $g_{a\gamma\gamma}$ [163, 164], as indicated in the right panel of Fig. 7. The projected limit from NA62 is expected to be weaker and thus not shown.

For sufficiently small $g_{a\gamma\gamma}$, the ALP might be long-lived and decay outside the detectors in high-energy colliders. There have been searches of single photon plus missing transverse energy $e^{+}e^{-}\to\gamma+\not{E}_{T}$ at LEP [166, 167, 168, 169] and $pp\to\gamma+\not{E}_{T}$ at LHC [170, 171, 172, 173]. Similarly, observations of radiative decays of Upsilon mesons $\Upsilon\to\gamma+\not{E}_{T}$ at CLEO [174] and BaBar [175] can be used to set limits on $g_{a\gamma\gamma}$. If the ALP decays promptly in the detectors, then we have the three photon signature $e^{+}e^{-},\,p\bar{p},\,pp\to\gamma+a\to\gamma\gamma\gamma$ at LEP [176, 177], Tevatron [178] and LHC [179, 180]. These limits could be improved by one to two orders of magnitude at future colliders such as Belle II [181], ILC [182], FCC-ee [183] and at later stages of the LHC [184, 185]. Benefiting from the large proton number in heavy ions, the photon-photon luminosity can be greatly enhanced in heavy-ion collisions compared to proton-proton colliders, and the current LHC bounds on $g_{a\gamma\gamma}$ can be improved by two orders of magnitude with ultra-peripheral heavy-ion collisions [186]. However, even at future colliders, the prospective limits are still too weak, at the level of $g_{a\gamma\gamma}\gtrsim 10^{-5}\,{\rm GeV}^{-1}$ [184, 185, 186], and hence not shown. Additional collider and flavor factory constraints on the coupling $g_{a\gamma\gamma}$ can be found e.g. in Refs. [39, 187, 188].

The following effective operators appear at dimension 6 and 7 [39, 100, 96]:

[TABLE]

with $C_{6,\,7}$ the dimensionless Wilson coefficients and $D_{\mu}$ the covariant derivative for the SM Higgs doublet. These operators generate ALP couplings $haa$ to the SM Higgs and $haZ$ to the $Z$ boson, which induce exotic decays of the SM Higgs, i.e. $h\to aa$ and $h\to aZ$ with the ALPs further decaying into two photons $a\to\gamma\gamma$ [100, 96]. In principle, we can set limits on $g_{a\gamma\gamma}$ from the searches for exotic decays of the SM Higgs $h\to aa\to 4\gamma$ and $h\to aZ\to\gamma\gamma\ell^{+}\ell^{-}$ (with $\ell=e,\,\mu$). However, these limits depend on the coefficients $C_{6,7}$ in Eq. (4.5), and we do not include them in Fig. 7.

Using Eq. (4.3), the GW bounds on $f_{a}$ from TianQin, BBO and CE can be converted to limits on the effective coupling $g_{a\gamma\gamma}$ that do not depend on the ALP mass $m_{a}$, and are shown as the purple, red and orange horizontal bands with dashed border lines in Fig. 7 respectively. When combined, these GW observations are sensitive to the range

[TABLE]

As shown in Fig. 7, some of the regions within this range of $g_{a\gamma\gamma}$ have been excluded by astrophysical and cosmological observations and laboratory experiments, while some are still unconstrained. Assuming that a GW signal from the PT in the ALP model is found in the near future, then we would expect a positive signal in future ALP searches. This can be seen in the right panel of Fig. 7 by the overlap between the regions covered by GW searches and those covered by telescopes [109], helioscopes [127, 128, 129], the LSW experiments ALPS II [143], ALPS III [144] and STAX [145], and beam-dump experiments like SHiP [162].

4.3 Coupling to electrons

Astrophysical and laboratory constraints on the ALP coupling to electrons $g_{aee}$ are collected in Fig. 8. In this figure we also show the parameter space for DFSZ and KSVZ axions. As in the case of $g_{a\gamma\gamma}$, exotic SM Higgs decays $h\to aa$ and $h\to aZ$ can not be used to set robust limits on the coupling $g_{ae}$, since they also depend on the coefficients $C_{6,7}$ in Eq. (4.5). The other constraints are described below:

• •

ALPs can be produced by bremsstrahlung and Compton effects in the Sun. Model-independent constraints on the mass $m_{a}$ and coupling $g_{aee}$ have been imposed by electron recoil searches in the low-background experiments Derbin [189], XMASS [190] and EDELWEISS [191]. The limit from EDELWEISS is the most stringent one, and is shown as the gray region in Fig. 8. There are also constraints from CoGeNT [192] and CDMS [193] on ALPs DM in local galaxies, which exclude however a much narrower region of $g_{aee}$. The limits on $g_{aee}$ from CUORE [194], Derbin [195] and Borexino [196] depend on the effective ALP coupling to nucleons $g_{aNN}^{\rm eff}$, and are not shown in Fig. 8.

• •

If the ALP couples to electrons, it will lead to extra energy losses in astrophysical objects. Constraints from observations of solar neutrinos [197] and Red Giants [198] have excluded a broad region in parameter space. The Red Giant excluded region is shown in green in Fig. 8, while the solar neutrino limits are comparatively much weaker.

• •

ALPs can be produced in beam-dump experiments by bremsstrahlung off an incident electron beam and decay back into electron-positron pairs in the detector [199]. The region excluded by E137 [157] is shown in pink in Fig. 8. The experiment MINOS/MINERvA could improve significantly the current limit [199], as indicated by the dashed gray line.

Given the relation in Eq. (4.3), the GW experiments TianQin, BBO and CE could probe the range

[TABLE]

A sizable fraction of this range for with $m_{a}\lesssim 10$ keV has already been excluded by EDELWEISS [191] and Red Giants [198]. Should the GW experiments TianQin find a GW signal, then the corresponding ALP mass would be expected to be heavier than roughly $10$ keV, which might be tested by the MINOS/MINERvA experiment [199].

4.4 Coupling to nucleons

Limits on the effective coupling $g_{aNN}$ of ALPs to nucleons are collected in Fig. 9, which also displays the parameter space for DFSZ and KSVZ axions.

• •

ALPs can be produced in compact astrophysical objects like neutron stars and supernova cores via nucleon bremsstrahlung $N+N\to N+N+a$, where $N=p,\,n$ represents both protons and neutrons [200, 201]. Neutron star constraints on the coupling $g_{aNN}$ can be found in e.g. [200, 202]. Limits from neutrino bursts from SN1987A are stronger, and exclude couplings $10^{-8}\lesssim g_{aNN}\lesssim 10^{-6}$ for $m_{a}\lesssim 100$ MeV [203, 204, 205, 206, 207, 208, 199], as shown in blue in Fig. 9. Next-generation supernova observations could improve greatly the limits on $g_{aNN}$, depending on how far the next supernova explosion is [209].

• •

The Yukawa couplings of ALPs to nucleons could potentially cause violations of the gravitational inverse-square law, and the effective coupling $g_{aNN}$ is thus constrained by Cavendish-type experiments [150, 210], as shown in brown in Fig. 9. Limits from measurements of Casimir forces are comparatively weaker with $g_{aNN}\lesssim 10^{-2.5}$ [211, 212], and are not shown.

• •

Searches of new long-range spin-dependent forces between nucleons can be used to set limits on the coupling $g_{aNN}$. Magnetometer experiments have excluded couplings $g_{aNN}\gtrsim 10^{-4}$ for ALP masses $m_{a}\lesssim$ meV [213], as shown in pink in Fig. 9.

The GW experiments TianQin, BBO and CE could probe the range

[TABLE]

which is largely complementary to supernova and laboratory constraints.

5 Prospects from precision Higgs data at future colliders

If the scalar $\phi$ in the ALP model resides at the few-TeV scale, it will contribute to the trilinear coupling $\lambda_{3}$ of the SM Higgs through the quartic coupling $\kappa$. This is obtained from the temperature-independent effective potential ${\cal V}$ in Eq. (2.6), after integrating out the $\phi$-field, and reads [214, 215]:

[TABLE]

with the SM contribution $\lambda_{3}^{\rm SM}=m_{h}^{2}/2v_{\rm EW}$. If the quartic coupling $\lambda_{a}\simeq{\cal O}(0.1)-{\cal O}(1)$ as seen in Fig. 5, the mass $m_{\phi}=\sqrt{4\lambda_{a}}f_{a}$ is of the same order as the scale $f_{a}$. Then we can set limits on the $f_{a}$ scale and $\kappa$ by precision measurements of the trilinear SM Higgs coupling at high-energy colliders. Current Higgs pair production data at the LHC lead to the limit $-9\lesssim\lambda_{3}/\lambda_{3}^{\rm SM}\lesssim 15$ [216, 217, 218], which is too weak to exclude any parameter space of the ALP model. Future hadron colliders like the high-luminosity LHC (HL-LHC) and the FCC-hh [219], and lepton colliders such as the ILC [182], will be able to measure the trilinear scalar coupling more precisely and probe the scale $f_{a}$ and the quartic coupling $\kappa$. Indeed, $\lambda_{3}$ can be measured within (30% - 50%) at the $1\sigma$ confidence level by the HL-LHC with an integrated luminosity of 3 ab*-1* [220, 221, 222, 223, 224]. With a larger cross section, the precision can be improved to $\sim 5$% at the future 100 TeV collider FCC-hh with a luminosity of 30 ab*-1* [225], and up to 13% at the 1 TeV ILC with a luminosity of 2.5 ab*-1* [226, 227]. All these sensitivities are shown in Fig. 10, for the benchmark value $\lambda_{a}=0.25$. Future high energy colliders are largely complementary to low energy axion experiments and GW observations for TeV scale $f_{a}$.

6 Conclusion

In this paper we have studied the production of GWs due to a strong FOPT in a generic axion or ALP model, where we extended the SM scalar sector by adding only a complex singlet field $\Phi$. The angular component of $\Phi$ is identified as the axion or ALP field $a$. The original Lagrangian contains only a few free parameters, namely, the ALP mass $m_{a}$, the “axion decay constant” $f_{a}$ and the quartic couplings $\kappa$ and $\lambda_{a}$ in Eq. (2.1). We have explored the prospects for GW emission for $f_{a}$ between $10^{3}$ GeV and $10^{8}$ GeV. Our numerical calculations reveal that in the ALP model we are considering, the GW signal strength could be as large as $h^{2}\Omega_{\text{GW}}\sim 10^{-12}$, which might be detectable at future GW experiments like TianQin, BBO and CE, depending on the GW frequency and on the ALP model parameters (see Figs. 4-6).

At low energies, the ALP couplings to SM particles are universally determined by the decay constant $f_{a}$, up to model-dependent coefficients; in other words, all the effective couplings of ALP depend on inverse powers of $f_{a}$. Therefore, we can convert the GW limits on $f_{a}$ to sensitivities on the effective ALP couplings to SM particles, independent of the ALP mass $m_{a}$. We have considered the CP-conserving couplings of ALP to photons $g_{a\gamma\gamma}$, electrons $g_{aee}$ and nucleons $g_{aNN}$. These couplings are tightly constrained by a large variety of laboratory experiments, and by astrophysical and cosmological observations, which exclude broad regions depending on the ALP mass $m_{a}$. GW experiments would probe sizable regions of the unconstrained parameter space, namely, $10^{-11}\,{\rm GeV}^{-1}\lesssim g_{a\gamma\gamma}\lesssim 10^{-6}\,{\rm GeV}^{-1}$, $10^{-11.5}\lesssim g_{aee}\lesssim 10^{-6.5}$ and $10^{-8}\lesssim g_{aNN}\lesssim 10^{-3}$, which are largely complementary to the laboratory, astrophysical and cosmological constraints. Thus, if a GW signal is found in future GW experiments and interpreted in the framework of axion or ALP models, it can be cross-checked in the upcoming laboratory and/or astrophysical ALP searches. In addition, for $f_{a}$ at the TeV scale, the real component $\phi$ contributes to the trilinear coupling of the SM Higgs. Thus precision Higgs data at future hadron and lepton colliders can be used to probe the $f_{a}$ and $\kappa$ parameters in the ALP model, which is also complementary to low-energy axion/ALP experiments and GW observations.

Acknowledgements

We thank Robert Caldwell, Yanou Cui, Ryusuke Jinno, Arthur Kosowsky, Marek Lewicki, Andrew Long, Alex Pomarol, Michael Ramsey-Musolf and Fabrizio Rompineve for useful discussions and comments on the draft. B.D. also thanks Aniket Joglekar for a discussion on the trilinear Higgs coupling. This work was supported by the US Department of Energy under Grant No. DE-SC0017987. Y.C.Z. is grateful to the Center for High Energy Physics, Peking University where part of the work was done for generous hospitality.

Appendix A Power-law Integrated Sensitivity Curves

Let us briefly outline the procedure used to compute the power-law integrated sensitivity curves for current and future GW experiments. For a detailed description of this method, see Ref. [228].

In the literature, the square root $\sqrt{S_{n}(f)}$ of the strain power spectral density is usually given as a function of frequency, in units of $1/\sqrt{\text{Hz}}$. We first convert it to $\Omega_{n}(f)$:

[TABLE]

Then, given a set of power-law indices $\beta$, e.g. $\beta\in\{-8,-7,...,7,8\}$, we compute for each $\beta$

[TABLE]

where $f_{\text{ref}}$ is some reference frequency. It can be arbitrarily chosen and it will not affect the results. $\rho$ is the integrated signal-to-noise ratio and $T$ is the observation time. Following Ref. [228], for $\rho=1$ and $T=1\text{ year}$ we have:

[TABLE]

The power-law integrated sensitivity curve $\Omega_{\text{PI}}(f)$ is the envelope of all the $\Omega_{\beta}(f)$ curves,

[TABLE]

As an explicit example, the power-law integrated curve $\Omega_{\text{GW}}(f)$ of BBO as well as the series of $\Omega_{\beta}(f)$ are presented in Fig. 11.