Hidden Photon and Axion Dark Matter from Symmetry Breaking

A light hidden photon or axion-like particle is a good dark matter candidate and they are often associated with the spontaneous breaking of dark global or gauged U(1) symmetry. We consider the dark Higgs dynamics around the phase transition in detail taking account of the portal coupling between the dark Higgs and the Standard Model Higgs as well as various thermal effects. We show that the (would-be) Nambu-Goldstone bosons are efficiently produced via a parametric resonance with the resonance parameter $q\sim 1$ at the hidden symmetry breaking. In the simplest setup, which predicts a second order phase transition, this can explain the dark matter abundance for the axion or hidden photon as light as sub eV. Even lighter mass, as predicted by the QCD axion model, can be consistent with dark matter abundance in the case of first order phase transition, in which case the gravitational wave signals may be detectable by future experiments such as LISA and DECIGO.


Introduction
Spontaneous symmetry breaking usually happens in the thermal history of the Universe. Within the standard cosmology, ΛCDM model, there are electroweak symmetry breaking (EWSB) and chiral symmetry breaking. The electroweak and chiral symmetries are restored at a high temperature due to the thermal effect but breaks when the temperature sufficiently redshift due to the expansion of the Universe. After the symmetry breaking, the standard model (SM) weak bosons and pions naturally appear with non-vanishing masses.
A clear evidence of new physics beyond the SM is the presence of dark matter (DM), the origin of which is a mystery of particle theory and cosmology. Except for its longevity, abundance, and coldness, most of properties, such as the mass, spin, interactions, are not known. The mass of the DM may be so small and the interaction between the DM and SM particles is likely to be so weak that it is consistent with the longevity and the nondetections in various experiments e.g. [1][2][3][4][5]. A simple possibility to realize both the lightness and weakness is that the DM is associated with a symmetry breaking at a high energy scale.
Such setups are naturally realized if the DM is an axion, axion-like particle (ALP) or hidden photon similar to the pion or weak bosons (see reviews [6][7][8][9][10][11][12]). Then, the interaction rates are suppressed by positive powers of the mass to the symmetry breaking scale. The question is how to produce these DM candidates in the early Universe.
In this paper, we discuss the possibility that a light axion (ALP) or hidden photon DM is produced via a hidden symmetry breaking. The simplest UV completion model to give ALP or hidden photon DM a mass is to introduce a dark Higgs field. The dark Higgs field is assumed to spontaneously break (approximate) global symmetry in the case of ALP and gauge symmetry in the case of hidden photon. This dark Higgs field is interacting with the SM particles and thus in the early Universe the hidden symmetry is restored due to the thermal effect. The symmetry breaking occurs when the temperature becomes small enough. We point out that the symmetry breaking is necessary followed by a parametric resonance production of the (would-be) Nambu-Goldstone (NG) bosons if the dark Higgs is not thermalized at the moment. In particular, we focus on a minimal setup where the dark sector and the SM sector is communicated only through the portal coupling between the SM Higgs and dark Higgs. A large fraction of the energy of the cold dark Higgs condensate is transferred into that of the NG bosons. Consequently, the cold DM abundance is explained. 2 A relevant topic may be the DM production from topological defects (see Refs. [40][41][42][43][44] for ALPs and Ref. [45] for hidden photons.) These defects may appear via a symmetry breaking. However, depending on the symmetry group or the breaking patterns, the defects may not appear like the case of the EWSB. In this case, our mechanism is more important. Even in a hidden U (1) symmetry breaking, which will be our concrete example, and which are also studied in the context of the topological defects, our mechanism provides complementary 1 The most misalignment production mechanisms of the hidden photon suffer from theoretical inconsistency or observational constraints [19]. 2 Ref. [32] also considered hidden photon production from parametric resonance effect induced by the dark Higgs dynamics. While Ref. [32] mainly focused on the regime of broad resonance, which corresponds to the case of large initial value of the dark Higgs field, we consider the small field regime corresponding to thermal phase transition. We will show that it generically leads to marginally broad or narrow resonance and also we take into account effects of the Higgs portal coupling.
parameter regions. Other than the ALP or hidden photon DM, heavy DM production is discussed relevant to the bubble wall dynamics in a first order phase transition (PT) by coupling the DM to certain Higgs fields [46][47][48][49][50]. Compared with those studies, our DM is much lighter than the symmetry breaking scale, and our mechanism works not very relevant to the bubble dynamics. In particular, our mechanism also works in a 2nd order PT or cross-over, where the bubbles are not created. This paper is organized as follows. In the next section 2 we will discuss the NG boson production at the symmetry breaking with a simple Higgs portal potential, which leads to a second order PT. In section 3 we show how the DM mass is generated and how the abundance can be explained via our mechanism. In section 4 we discuss the case for a 1st order PT and the gravitational wave. The last section 5 is devoted to conclusions and discussion.

NG boson production at symmetry breaking
Let us consider the the spontaneous symmetry breaking of a hidden global continuous symmetry in the early Universe. Later we will gauge or explicitly break this group to give mass to the (would-be) NG boson. In this part, we show that the NG boson can be efficiently produced soon after the symmetry breaking or PT if the dark Higgs is not thermalized at the moment.

Zero-temperature potential of dark and SM Higgs fields
To be concrete, we consider a minimal dark sector in which there is one dark Higgs field which spontaneously breaks the hidden global U (1) symmetry. In this minimal setup, the only renormalizable interaction between the SM and dark sector is the portal coupling between the SM and dark Higgs fields. The most general dark and SM Higgs potential is given by Here Φ (H) is the hidden (SM) Higgs field (doublet) which will break the U (1) (SU(2) L × U (1) Y ) symmetry, λ P (> 0), λ(> 0) and λ H (> 0) are coupling constants, µ 2 H (125 GeV) 2 /2 is the bare Higgs mass term in the SM, and m 2 Φ (> 0) is the dark Higgs mass squared parameter. Λ 4 is needed to cancel the cosmological constant. Here is the dark Higgs vacuum expectation value (VEV) by introducing which in the last term we have cancelled the contribution to the SM Higgs boson mass. We will discuss the tuning to the SM Higgs boson mass, later. The interaction between the dark sector and the SM sector is controlled by the portal coupling constant, λ P . Note that the portal coupling in (1) may ensure the absolute stability of the electroweak vacuum [51,52]. Assuming that Φ is much heavier than the electroweak scale, we can integrate out it below the scale ∼ m Φ to obtain the effective four-point coupling constant of the SM Higgs as λ eff λ H − λ 2 P /(2λ). If m Φ 10 10 GeV and the following condition is satisfied, it is shown that the quantum-corrected effective potential never becomes negative. If the condition (3) is not satisfied, either low-scale inflation or high-scale inflation with some additional SM Higgs interaction is required in order to avoid the collapse of the vaccum [53][54][55][56][57][58][59]. On the other hand, if the condition (3) is satisfied, we need not to worry about such details of the inflaton and Higgs dynamics.

Phase transition and dynamics of dark Higgs
In the early epoch, the Universe is filled by hot and dense plasma. Here we assume that the reheating occurs in the SM sector and thus the dense plasma, characterized by the temperature T , is composed by the SM particles. We assume that Φ is not fully thermalized before the symmetry breaking, i.e.
with g being the effective relativistic degrees of freedom, and M pl ≈ 2.4 × 10 18 GeV being the reduced Planck mass. Here the thermalization rate of the dark Higgs field is given by We will show that this condition before the PT is important for our DM production mechanism to work. At high temperature, the radial component, s = √ 2|Φ|, gets a thermally corrected effective potential as [60,61] where · · · represents irrelevant terms, including the Coleman-Weinberg corrections as well as the higher order terms. 3 Here is one remark. In this case we do not have a cubic term of s from the ∼ −m H (s, T ) 3 T term in the free energy density. This is because the SM Higgs mass is m 2 H ∼ T 2 + λ P (−v 2 Φ + s 2 /2) + 2µ 2 H , the first term of which comes from the daisy resummation at s ∼ 0 [60]. Due to the contribution −λ P v 2 Φ , which is required to cancel the SM Higgs mass at the vacuum s ∼ √ 2v Φ , m 2 H cannot be approximated as m 2 H ∼ λ P s 2 for s √ 2v Φ . Thus there is no parameter region for s 3 term to appear. This is a peculiar feature of the SM Higgs contribution to the thermally corrected potential: the SM Higgs is (almost) massless at the finite VEV of s = √ 2v Φ . Thus we expect that the PT of Φ is the second order.
As the Universe expands, the temperature T decreases. One can easily see that the symmetry is broken at the temperature As we have explained above, there is no cubic term or potential barrier in the thermal effective potential and hence the PT is expected to be the 2nd order. A similar discussion can be also made to the SM Higgs potential, which is broken while U (1) is symmetric if We can easily find that the symmetry breaking of U (1) occurs prior to the electroweak symmetry breaking if This condition will be assumed, so that we can safely neglect the dynamics in the H direction. Now let us see the s dynamics around the PT. When T T crit , the VEV of Φ is temper- with being the temperature dependent effective mass. The radial component of the dark Higgs s may follow the potential minimum just at around the transition, but after the transition s 3 Here it is assumed that λ P λ. As will be explained later, we will mainly consider the phenomenologically preferred case of λ 2 P ∼ λ. Then, as far as both λ P and λ are smaller than unity, this assumption is justified.
starts to oscillate around the temperature dependent minimum if The energy density of the coherent oscillation comes from part of the potential energy Λ 4 . This can be seen by solving the equation of motion in a simplified setup by neglecting the contribution from the NG bosons:s We will solve this equation by taking the initial conditions s(0) = 0, anḋ We setṡ with a tiny non-vanishing value initially, because we expect that there is a thermal fluctuation, which kicks s at random. If Φ is initially thermalized, we expectṡ ∼ T 2 . We take into account of the thermal fluctuation by this initial condition. As we will see soon that this initial condition is insensitive to the result as long as (4) is satisfied. In Fig. 1, we show the numerical result of s/( 14 GeV, the oscillation takes place within a few 1/m Φ which is much shorter than one Hubble time (m Φ /H ubble (T crit ) 2338 in this case). After the onset of oscillation, the oscillating amplitude decreases in time. We also notice that v T Φ is settled into v Φ within O(1) Hubble time. This is the case Γ Φ th /H ubble 2.58 × 10 −22 . When m Φ H ubble , and Γ Φ th /H ubble are larger, the transition is faster, as shown in the right panel. Here m Φ = 0.1 GeV, v Φ = 10 8 GeV and Γ Φ th /H ubble 2.58 × 10 −7 . In this case the oscillation amplitude is much smaller To understand this behavior, let us consider two time periods for the s evolution. Soon after the PT, s is placed at the hilltop with a non-vanishing negative mass squared. Then s starts to slow-roll towards the bottom of the potential as long as there is a tinyṡ or s at T T crit . The slow-roll lasts much longer than 1/m Φ because of the vanishingly small V T at around the hilltop. The time, ∆t osc , for the onset of the oscillation (measured from the instant of T = T crit ) can be obtained by solving the equation of motion with neglecting the Hubble friction. We can take s ∝ exp t 0 m T Φ dt . When the exponent becomes larger than O(1), the slope becomes so steep and evetually s starts to oscillate. Since the evolution is in exponential, the result will not be very sensitive to the initial condition. By expanding m T Φ around T ∼ T crit , we obtain the time scale for the onset of oscillation   Figure 1: The time dependence of s and v T Φ around the 2nd order PT given in red and blue solid lines, respectively. Effects of particle production and dissipation/decay are neglected. In the left panel we have taken m Φ = 0.1 GeV, v Φ = 10 14 GeV, and in the right panel after the PT. The coefficient,C, which is not very different from O(1), logarithmically depends on the initial condition (13) when it is small enough. After the onset of oscillation, the number density is an adiabatic invariant, Then we can define for later convenience. It represents an effective suppression factor of the s coherent oscillation abundance. It is evaluated as C eff ∼ 10C 3/2 m Φ /(M pl λ P ) by inserting v T Φ at the timing of onset of s oscillation (14).
The contours of C eff by numerically solving the equation of motion is given in the (m Φ , v Φ ) plane in Fig. 2. For simplicity we take λ = λ 2 P , in which case C eff scales as 10(C 3/2 v φ /M pl ). We find a non-negligible number of n s produced due to the oscillation of s after the PT. We emphasize that the discussion is based on the condition that Φ is not thermalized or Γ Φ th H ubble at the PT. When Γ th and H ubble are close to each other soon after T = T crit , s starts to oscillate without being trapped at around the potential top. This will lead to a different conclusion with a much suppressed NG boson production.
Here we emphasize that we neglect the particle production of the (would-be) NG boson for illustrative purpose. In practice, we cannot neglect the interaction with the NG boson. 16 In the grey shaded region the dark Higgs is (close to be) thermalized at the phase transition or our estimation here is invalid.
As we will show soon, the NG boson production happens at a similar time scale of the oscillation frequency, m T Φ .

Particle production at the phase transition
Let us consider the particle production at around the onset of oscillation, i.e. at around the symmetry breaking. By The decay rate to the SM Higgs multiplets is given by 4 The results are the same for the non-linear parametrization Φ = v Φ e (s+ia)/ by neglecting the SM Higgs boson mass. However, when T m Φ , this process is kinematically forbidden due to the heavy thermal mass of H. Instead, as we shall see later, a thermal dissipation effect may work.
Parametric resonance production of a-condensate The decay rate to the NG mode pair, on the other hand, is given by By taking into account the thermal corrections, m Φ in the r.h.s should be replaced to be m T Φ . With oscillating s, we should take into account the parametric resonance effect for s → a's since there could be a Bose-enhancement effect [37,38,[62][63][64][65][66][67][68][69][70][71][72][73][74]. The resonance q-parameter [67] can be obtained for the s-a system soon after the PT as where s amp denotes the oscillation amplitude of s and we have used s amp v T Φ . In particular, this is maximized at the onset of oscillation, where s amp ∼ v T Φ . During the parametric resonance the phase space distribution function of a in the resonance band increases expo- Thus, the growth rate q a m T Φ is as fast as the oscillation m T Φ at the onset of oscillation. This is a generic prediction associated with the Higgs dynamics around the symmetry breaking in the early Universe. Since m T Φ is much larger than the Hubble scale (Eq. (11)), the parametric resonance effect transfers the most of the s oscillation energy density to the NG mode a after O(1) oscillations of the s field. The resulting number density of a from this effect is estimated as with S being the entropy density of the SM plasma, g ,s being the relativistic degrees of freedom for the entropy density and n (para) a the produced number density of a. Here we have taken the temperature at the onset of oscillation to be T crit .
The parametric resonance effect for the s-H system, on the other hand, is suppressed again due to the large thermal mass of H. Although s cannot decay into H, s can scatter with the thermal plasma at a rate [75][76][77][78][79][80][81][82] where we have considered the production of a top quark pair with y t ∼ 1 being the top Yukawa coupling by taking account of the Higgs thermal mass. Here we have neglected the temperature dependence on v Φ for the time scale to consider. This can be also found from the equation of motion of s in an effective action (see appendix A). Through this process, the energy stored in s may be dissipated into the SM plasma. In order for the parametric resonance effect not to be blocked by such a scattering process, we need [67] q a m This is easily satisfied soon after the onset of oscillation (see Eqs. (8) and (10)). Consequently, soon after the PT, s has a good environment for producing a via parametric resonance.
Dissipation of a-condensate We also need to consider the dissipation of produced a soon after the PT. The dissipation rate of the k-mode of the NG boson is found to be (see the Appendix B for derivations): This can be seen as the a-a annihilation and thus it is proportional to the number density of the k-mode NG boson. This effect is most important at the NG boson production k ∼ m Φ /2 ∼ √ λ P T crit /2. 6 We can neglect the annihilation effect if soon after the production, so that a is kept intact. 5 We note this is different from some of the result in the references due to the time scales where we consider m Φ Γ th , T , with Γ th being the thermalization rate of the Higgs, Γ th ∼ y 2 t T . For example, Eq. (27) of Ref. [75] implicitly assumes Γ th m Φ T and it is different from our expression. On the other hand, Eq. (3.34) of Ref. [78] is derived under the same assumption as ours and hence the result is consistent. 6 Strictly speaking, k ∼ C 1/3 eff m φ /2 in the second order phase transition. This is because at the production of the NG boson, the mass of s is smaller due to the aforementioned thermal correction. As we will see, the annihilation effect, even that it is overestimated, is not important in the second order case. On the other hand, in the first order PT case in Sec. 4, the production of a is delayed and the momentum is not suppressed by C If this is not satisfied, some of the NG bosons annihilate into the SM plasma, while some NG bosons remain. The remnant n a can be estimated from Γ like the WIMP scenario. Therefore soon after the PT we get which originates from the PT. Much after the PT, the annihilation effect on a-condensate is suppressed and we will neglect it.
Before moving to the next section, let us discuss some other components produced relevant to the PT.
Remnant of s-condensate The efficiency of the parametric resonance for producing the NG mode a decreases as the amplitude of s decreases due to the energy transfer into a. Eventually the parametric resonance stops and a tiny component of the s-condensate should remain. It is estimated as follows. The resonance peak is at k m T Φ /2 while the width of the resonance band in the momentum space is given by ∼ q a m T Φ . The redshift of the momentum and the enhancement of m T Φ due to the Hubble expansion takes the produced a away from the resonance band within a time scale ∆t ∼ q a /H ubble . The resonant enhancement stops when the exponential growth factor q a m Φ ∆t ∼ q 2 a m Φ /H ubble ∼ 1. From this we can estimate s amp at which the parametric resonance stops and hence the remnant of the energy density of s as [38,67] One can also estimate the order of it from the Boltzmann equation by taking account the Bose-enhancement effect [37]. 7 As far as the narrow resonance with q a 1 is concerned, we expect that the back-reaction such as aa → ss, or sa → sa is not important because it is either kinematically invalid or the rate is suppressed by q 4 a . On the other hand, due to the tachyonic instability the fluctuation of s itself may also develop within a few oscillations [83,84], which may tend to stop the resonant enhancement of a. Still, however, the conclusion that the most energy of s is transferred to a should remain valid. The condition (22) may not be satisfied for a smaller q a . In such a case, before the Hubble expansion becomes important in preventing the production of a, the dissipation may be more important. In this case, the s condensate is easier to thermalize than our estimation which does not change our conclusion.
Topological defects After the PT, topological defects may be formed. In our U (1) case, there are cosmic strings produced after the PT. Cosmological effects of cosmic strings in our scenario will be briefly discussed in the next section.
3 Dark matter production at second order phase transition Let us apply the mechanism of NG boson production discussed in Sec.2 to the DM production. The DM, if dominant, must be cold and thus we should somehow give mass to the NG boson to make it non-relativistic around and after the galaxy formation era. In Secs. 3.2 and 3.3, we will provide two possibilities for generating masses of the DM: explicitly breaking the global U (1) and gauging the U (1), in which case the DM becomes an axion(-like particle) and hidden photon, respectively. In the latter case, we can discuss the most properties of DM by looking at the NG boson Lagrangian according to the equivalence theorem since we are interested in the light hidden photon DM and it is highly relativistic at the production. Thus, in the Sec. 3.1, we first discuss general model-independent features by assuming that the NG boson acquires a mass term of m a and discuss the thermal history after the PT.

Dark components after phase transition
(Not much) After the PT, we have five kinds of cosmic components other than the SM particle plasma: the a condensate from parametric resonance, the remnant s condensate, the topological defects, a produced from thermal scattering, and s produced from thermal scattering. We use the "condensate" to distinguish the "cold" component, whose typical momentum is much smaller than the cosmic temperature T , which is the typical momentum of a "particle" from the thermal scattering. They will be discussed separately.
As we will discuss soon, the remnant of s may be dissipated away due to (21), decay via Eqs. (17) or via the mixing with the SM Higgs boson when the Universe cools down. For simplicity, let us assume that the remnant of s condensate does not dominate the Universe and its subsequent interactions do not play important role on cosmology. 8 Thus, the s will neither contribute to nor dilute the DM abundance. Due to this assumption, we can first calculate the DM. We will check that this assumption is satisfied in the parameter region of interest. We will also come back to the case that s once dominates the Universe in the last section, by considering s as an inflaton.
a-condensate as dominant dark matter component The produced a condensate later composes the DM when it acquires the mass m a and becomes non-relativistic. We can calculate the abundance of the (would-be) NG boson a from Ω a ∼ m a n a S S 0 ρ c , where S 0 (ρ c ) is the entropy density (critical density) today. This explains the observed DM abundance if [86] Ω a h 2 = Ω DM h 2 ∼ 0.12, with h 0.67 being the present Hubble parameter in unit of 100 km/s/Mpc. Also, to explain the coldness of the DM we use the conservative bound calculated in [37] (see also Refs. [87,88]), which gives a lower bound on the DM mass. Interestingly, since this becomes by using Eq. (28), the coldness bound is automatically satisfied from Eq. (8). Thus the DM from the symmetry breaking is naturally cold. For explanation of the effects of the following constraints from the thermal history, we first show the contour plot of the DM mass given in Fig. 3. Again here we take which is the largest λ 2 P satisfying (8), and corresponds to (almost) the lightest DM according to (30). Note that this choice is consistent with the condition for the absolute stability of the electroweak vacuum (Eq. (3)).
Thermal history for the remnant s-condensate and s particle To discuss the evolution of other components, let us introduce T s→HH and T s→aa which are defined by −1/4 M pl Γ i . As we have explained that the decay s → HH is thermally blocked and dissipation is important. T s→HH should not be considered as the decay temperature.
The dissipation rate (21) is smaller than Γ s→HH with T for T m Φ . Therefore the dissipation can remove s condensate away if and only if If (33) is satisfied, in the figure, T th is always greater than m Φ and the electroweak scale. Thus s-condensate evaporates.
We must also consider the thermal production of s particles since the production rate, which is dominated by the inverse decay, is given as Γ HH→s 4πT . The production rate via tt → sH has a similar form. This is comparable to Γ sH→tt . Then at T ∼ T th , the scondensate disappears, but, instead, s-particles are thermalized. The thermalized s mostly interacts with the SM particles if T s→aa < m Φ , until s becomes non-relativistic. On the other hand, if T s→aa > m Φ , a particles are produced via the s decay and a are also thermalized. In the end, s would decay to SM thermal plasma. Since we focus λ ∼ λ 2 p , the decay rate to aa is smaller than the decay rate to HH if λ 1. 9 Since Eq. (33) in the parameter region of interest, the components of s condensate and particles disappear from the Universe not much later than T ∼ m Φ . The decay of s should not cause cosmological problems as long as they happen at a high enough temperature. In particular, we take from the viewpoint of the big-bang nucleosynthesis [89][90][91][92][93][94][95][96][97][98]. 10 This is the lower limit of the horizontal axis of the figure. After the decoupling/decay of s, the Universe is composed by three components: a-condensate, a-particles, and topological defects. We also mention that s in the sub GeV mass range, which mixes with the SM Higgs with a mixing angle θ H ∼ λ P v Φ /m h , θ H 10 −3 , may be excluded by the accelerator bounds or BBN constraint. A large fraction of the allowed range may be tested in the SHiP experiment [117].
Freeze-in production of a Although s dominantly decays into SM particles (via mixing with the Higgs if it is lighter than 2m h ), the rare decay into aa provides a freeze-in production of DM. The produced abundance of a can be estimated as 9 We note that we may also consider the s decay to aa when λ 2 λ 2 P in general. 10 Since the mass is relatively heavy, we neglect bounds on s particles from stellar cooling arguments.
This explains the DM abundance, Ω th a ∼ Ω DM with m a shown on the contours below the red solid line. However it is subdominant above the red solid line. Notice that the produced DM tends to be warm and is intension with the Ly-α data for m a O(10) keV. Therefore, the freeze-in region is disfavored.
Constraints from topological defects/coherent oscillation The topological defects or coherent oscillation contributes to the DM abundance depending on the nature of the DM mass. When v Φ is sufficiently large, these contribution cannot be neglected. Therefore we do not consider the region above the blue dashed band. These production will be discussed in more details in later in this section.
Irrelevant constraints and consistency Before ending this section let us mention some constraints that are irrelevant and not shown in this figure. a is in kinetic equilibrium with the thermal plasma, if the scattering of a process aH → aH is too fast. The scattering rate is given by This form is justified when E a T m 2 Φ , where E a ∼ m Φ (T /T crit ) is the energy of the produced NG boson energy. When E a T m 2 Φ , it is much slower. Above the red line for the freeze-in, this process is always slower than the Hubble expansion.
Generally, there is another contribution to the freeze-in production of a from direct thermal scattering. The production rate via the portal coupling is given as When T m Φ , this is suppressed since the NG boson-Higgs interaction comes from the higher dimensional term λ P 2m 2 Φ |H| 2 (∂a) 2 , which is generated by integrating out s. When T m h , it is suppressed by a Boltzmann factor. 11 In the parameter region of focus, this production is subdominant compared with the a production from the decay of the thermally produced s.
Since we assumed that s never dominates the Universe to estimate the DM abundance, i.e. the remnant of s does not dominate the Universe at T th T T crit , we need to check whether this is the case. In fact, this condition gives an upper bound of v Φ which is much higher than the bound from topological defects/coherent oscillation. At the PT s should oscillate, i.e. m Φ H ubble (T crit ), which is also satisfied in the shown region. Lastly let us mention the fine-tuning on the SM Higgs boson mass. The dark Higgs field acquires a large vacuum expectation value which contributes to the SM Higgs boson mass via the portal coupling. It may be one of the sources of the fine-tuning problem of the SM Higgs mass, if this contribution is much larger than the electroweak scale. Interestingly, in the viable parameter region the portal coupling contribution is negligible compared with the SM Higgs boson mass. In this sense, it may be viewed as a natural parameter region.

Axion production via Peccei-Quinn symmetry breaking
Having discussed generic feature of the NG mode production at the symmetry breaking, now we look into more details of the case of axion DM. Suppose that the global U(1) symmetry is explicitly broken by a small amount, which gives a potential for the NG mode, axion. The axion potential is assumed to be of the form This can be either made if the "Peccei-Quinn" (PQ) field Φ [103][104][105][106], which takes a role of dark Higgs field discussed so far, anomalously couples to some non-abelian gauge fields, which generate the axion potential due to non-perturbative dynamics or with some explicit breaking term like δL ∝ Φ + Φ † . Here f a = √ 2v Φ /N DW with N DW being the domain wall number, and we take N DW = 1 to evade the cosmological domain wall problem. In this case, domain walls are temporary formed at the onset of the oscillation of a, i.e. at m a ∼ H. However, each domain wall is bounded by a string. Soon after the domain wall formation the wall tension dominates the dynamics of the string-wall system and the domain walls collapse.
In the axion model, we have two additional sources of the DM production other than that we have discussed so far, i.e. production at the symmetry breaking. One comes from the misalignment mechanism [13][14][15] i.e. from the axion coherent oscillation. The abundance is estimated as where we have taken the misalignment angle θ a = π/ √ 3, and we have assumed a temperature independent potential of V a . This contribution can explain the DM around the blue dashed line, above which it dominates over our PT production. The other is the ALP radiation when the domain walls collapse [40][41][42][43]. This contribution is more or less comparable to the misalignment one. Since there is a theoretical uncertainty on the numerical estimation of this contribution, we simply assume that these two contributions are the same order and we just use (39) as a representative one. These contributions would dominate over our production mechanism at high v Φ (much) above the blue dashed line in Fig. 3. This turns out to be subdominant due to the small decay constant in the region of interest (i.e. below the blue dashed line).
So far we have implicitly assumed that m a and f a are independent. In the case of the QCD axion, which is well motivated from the viewpoint of strong CP problem in QCD, the potential is generated via the non-perturbative dynamics of the QCD and hence m a and f a are related. In this case we have m a f a = √ χ 0 , where χ 0 is the topological susceptibility, which we adopt χ 0 ≈ (0.0756 GeV) 4 [107] (See also Refs. [108][109][110][111][112]). The region compatible with this relation cannot be found in this figure because the DM is too heavy to be the axion. Strictly speaking, in the case of QCD axion, we need to take care of the existence of additional particles and topological defects, depending on the concrete UV completion model. In the KSVZ model [113,114], we may have thermalized light PQ quarks in the symmetric phase. 12 There is no domain wall problem in the KSVZ model since N DW = 1 with a minimal number of PQ quarks. On the other hand, there is a domain wall problem in the DFSZ scenario [115,116] in which N DW = 6. To solve the problem we may introduce a tiny PQ breaking term in order to let the domain walls collapse soon after the onset of the coherent oscillation of the axion. Note that in the DFSZ model there is an additional Higgs doublet coupled to the PQ field, and hence the thermal potential discussed so far may be different, which may lead to a first order PT. We will come back to the possibility of the first order PT in Sec. 4.

Hidden photon production via hidden U (1) breaking
Next we discuss the case of gauging the hidden U (1) symmetry in order to make the NG mode massive. The Lagrangian of the hidden sector, including the Higgs-portal coupling, is given as where F is the field strength of the Hidden photon, g the gauge coupling, D µ Φ ≡ (∂ µ + igA µ )Φ is the covariant derivative of the dark Higgs. In this case, we can still calculate the (longitudinal component of the) hidden photon DM abundance from (28) thanks to the equivalence theorem, by taking m a = m A = √ 2gv Φ with m A being the mass of A. Since we are interested in the case of very small g, and since the interaction of the dark Higgs to the transverse gauge boson is suppressed by the coupling g, we can safely neglect the production of transverse mode.
In principle we can write down the kinetic mixing term, F µν F µν Y , between gauge fields of the SM U (1) Y and hidden U (1). However this can be neglected if we take g small enough or assume a charge conjugation symmetry in the hidden sector, A → −A, Φ → Φ † to forbid the kinetic mixing. In either case we do not need to care the thermal production of (transverse components of) the hidden photon via the gauge interaction.
In the hidden photon model, there is an additional contribution to the hidden photon abundance from the cosmic string network formed during the symmetry breaking. As shown in Ref. [45], cosmic string loops emit (longitudinal component of) the hidden photon as far as the loop size is smaller than m −1 A . This production is dominant above the blue dashed line, which is taken from Ref. [45]. Another contribution may be from the inflationary period or (pre)heating [22][23][24][25]. This component, however, is sensitive to the inflation scale and the reheating dynamics, and is subdominant if the inflation scale is not very high and not shown here.
An important difference between the axion case and hidden photon case is that cosmic string networks remain until present day in the latter case. The cosmic string tension is constrained by several observations. A robust constraint comes from the CMB observation, which indicates v Φ 2 × 10 15 GeV [118]. The cosmic string networks necessarily produce string loops in order to maintain the scaling solution, and string loops emit gravitational waves [119,120]. In the present case, because of the smallness of the hidden photon mass, loops lose their energy dominantly through the emission of longitudinal vector boson if the loop size is smaller than m −1 A and through the gravitational waves if the loop size is larger [45]. There are orders-of-magnitude uncertainties of the typical loop size, but for wide range of parameters the string loops contribute to the stochastic gravitational waves at the nano-frequency range, at which pulsar timing arrays have a good sensitivity. The recent NANOGrav result [121] gives an upper bound on the symmetry breaking scale as v Φ 5 × 10 13 GeV if the loop size is about one-tenth of the Hubble horizon scale, but it is relaxed as v Φ 10 15 GeV if the loop size is smaller [122]. If the symmetry breaking scale is close to this upper bound, it is possible to explain the NANOGrav evidence of the gravitational waves. 13 4 Light dark matter from first order phase transition So far we have considered a simple setup of the second order phase transition of a hidden U (1) global or gauge symmetry in the early universe. This is true if the dark Higgs only has a portal coupling to the SM Higgs boson. On the other hand, the dark Higgs field may also 13 Since the dominant contribution to the gravitational waves at the NANOGrav frequency range comes from loops that is going to decay at present, size of such loops is large enough to forbid the emission into the hidden photon. However, on the very high frequency range, at which laser interferometer gravitational wave detectors are sensitive, the signal may be greatly reduced due to the emission into the hidden photon. Such correlations between the low and high frequency gravitational wave signals may be a smoking-gun of this scenario.  have other couplings in general. In particular if the NG boson is the QCD axion, it should either coupled to heavy Higgs boson in the DFSZ model or PQ quarks in the KSVZ model. The inclusion of the new thermal and Coleman-Weinberg contributions to the potential may lead to a first order phase transition. In the first order PT, s stays longer around the hilltop of the potential, then undergoes a tunneling and starts to oscillate. Therefore the suppression factor C eff (16) tends to be close to one.
Strictly speaking, in such a Higgs potential that leads to a first order PT, the bubble wall may take away a fraction of the energy stored in the potential in the symmetric phase, like the well-known reheating problem in the old inflation. Then the Higgs oscillation amplitude in the broken phase should be suppressed according to energy conservation. This bubble wall expansion, however, gets a friction due to the pressure induced by the interactions the outside thermal plasma and wall, and reach a terminal velocity [123][124][125][126] (see also Ref. [127].) in which case, we expect that the Higgs field in the broken phase exhibits a coherent oscillation.
In this section, let us simplify the discussion with the assumption that the PT takes place not too later than T ∼ T crit , and the number density of the s oscillation is given by Eq. (20) with C eff taken as a free parameter to take account of the model-dependence and the uncertainty due to the bubble wall dynamics. Moreover we neglect the effect on the NG boson production due to the bubble wall dynamics. By these assumptions, our previous discussions remain intact.
In Fig.4 we show the parameter region with C eff ≈ 1 [left panel] and C eff = 0.1 [right panel], by particularly focusing on the QCD axion range. Below the blue solid line, the number of the NG boson is produced too much initially, and the annihilation takes place promptly. The abundance is given by Eq. (25). We find that 0.01 C eff 1, the QCD axion produced by the PT can explain the present DM abundance.
In general, during the first order PT, gravitational waves are produced via the bubble collisions or plasma sound wave [128,129]. The typical frequency of the gravitational wave is determined by the bubble size, denoted by β −1 , at the collision. It depends on the details of the dark Higgs interactions, and it is estimated as where T reh is the thermal temperature at the completion of the PT. In our scenario, T reh ∼ T crit . In Fig.4 contours of the critical temperature ∼ T reh are shown by the green dotted lines. This does not depend on C eff . Interestingly, when the QCD axion DM is successfully produced during the PT, we obtain 1 mHz f GW 10 Hz taking account of the model dependence of β/H ubble ∼ 10-10 3 , which may be within the sensitive range of LISA [130] and DECIGO [131].

Conclusions and discussion
In this paper we have proposed the hidden photon or axion-like particle DM production via continuous symmetry breaking with a dark Higgs, taking account of the interaction between the dark Higgs and the SM Higgs. It is a minimum setup that accounts for the dark global or gauged U(1) symmetry breaking. Even in this simple setup, the dark Higgs dynamics and its consequence for the DM production are complicated partly due to thermal effects. We found parameter regions that are consistent with present DM abundance. In our scenario the DM mass can be as light as 1 eV. The light DM may be warm and can be tested in the future observations of the 21cm line [133]. On the other hand, given a setup that the PT is the first order, the DM can be much lighter and there is a possibility that the QCD axion produced by the dark Higgs dynamics takes a role of DM. In this case, the gravitational waves from the PT may be tested in the future.
Some additional comments are in order. In the main part we have considered the case where the dark Higgs does not dominate the energy density of the Universe before the PT. An interesting alternative possibility may be that s is the inflaton, which means s dominates the Universe and must reheat the Universe later. Let us suppose that U (1) is gauged and g is chosen so that Then the curvature at the hilltop of the potential may satisfy according to the weak gravity conjecture (WGC) [134,135]. 14 The least tuned region saturates the WGC [135,136], g 2 M 2 pl = V Φ . Thus we may have a local tiny minimum at the potential maximum. An old inflation takes place there (for the e-folds, see e.g. Refs. [137,138]) and later s tunnels through the potential barrier. After the tunneling, still the curvature of the potential may be suppressed enough and then the quartic hilltop inflation happens there [139][140][141] if the quartic coupling is negative. Note that the WGC required the potential to be flat and the slow-roll condition is satisfied. The potential can have a minimum stabilized by the quartic term and a higher dimensional term [21,138,[142][143][144][145]. 15 Soon after the slow-roll inflation ends, the hidden photon DM is produced via the parametric resonance as we have discussed in the main part. However, as we have also mentioned, the s-condensate may not completely disappear and remain slightly. This may dominate the Universe again at the later stage, and then its decay reheats the Universe again. The detail of the reheating is complicated due to thermal corrections and we leave it for our future study. In any case, irreverent to the detailed thermal history, an unique prediction is the relation between m A and H inf , In particular, if the higher dimensional terms are suppressed by M pl , v Φ ∼ 10 12 GeV and H inf ∼ 1 GeV [21]. This predicts m A ∼ keV and g ∼ 10 −18 .

A Parametric resonance in thermal environment
Here let us discuss whether significant particle production via parametric resonance may occur in a thermal environment. This is important since in a large parameter region the 14 We are not sure if the WGC can work in a false vacuum. However, this issue also exists in the original paper explaining the hierarchy problem of the SM since the electroweak vacuum is essentially false vacuum in the SM [135]. 15 The inflaton, s, can also be stabilized by two or more higher dimensional terms. In such a case the VEV and the mass of s are the typical scales of the higher dimensional terms [146].
resonance parameter for the s-H system is larger than unity. As a toy model, we consider where g is a portal coupling, H is a scalar complex field and is massless at the vacuum s = 0. In this Appendix we neglect the expansion of the Universe for simplicity. Later we will discuss the case where H is the SM Higgs, but for a while we keep H just as a general complex scalar field. As is well known, particle production of H happens due to the s coherent oscillation. When H is coupled to thermal bath, this process becomes more involved.
Let us first consider the case H does not interact with any other particle and it is only produced by s. Then a broad parametric resonance occur if q ∼ gs 2 amp /m 2 Φ 1. The number density of H, n H , increases exponentially. We emphasize that s is now wave-like, and we cannot describe the evolution from perturbation theory for particles. For instance one can easily find the s n → HH process has a rate proportional to q n m Φ , and the particle picture is highly non-perturbative. We can describe the evolution of s by solving the equation of motion (in the 1PI effective theory). Following [67], we can write down the equation of motion (EOM) of s ass where H 2 ≡ d 3 k (2π) 3 |H k | 2 , and · · · represents terms with higher order in s, s k . These neglected terms would be important when resonance lasts long enough and the processes known as re-scattering would occur, which, however, is not our focus. We emphasize that this EOM should include all the effects (of non-perturbative series in the particle picture) which only involve the s zero modes. By With this, the EOM looks like that s moves in an effective potential of V eff ∼ m 2 Φ s 2 /2 + √ gn H |s|. If n H increases s amp should decreases so that V eff [s amp ] is kept. Since n H increases exponentially due to the parametric resonance, the amplitude of s is decreased. Now let us consider the parametric resonance in thermal environment where H is thermalized with a temperature T and suppose that the s oscillation time scale is much longer than the thermalization time scale: m Φ T so that the s oscillation is nearly adiabatic with respect to the thermal bath. There may be a coupling of H to other fields like gauge bosons, with typical coupling of g ∼ 1, which may induce the H's thermal mass of ∼ g T . Then the total mass of H is expressed as The energy density of H is expressed as Since we assume that ρ s is smaller than ρ H , ρ H is kept almost constant during the oscillation of s. This means that the two point function, satisfying |H| 2 E 2 H ∼ ρ H , is bounded by In the left hand side we divide ρ H by maximal Higgs energy T 2 while in the right hand side we divide it by the minimal one, m 2 H . This implies that |H| 2 cannot change much for g ∼ 1 and hence s amp does not decrease much.
Let us estimate the dissipation rate of s in this setup following the arguments in Refs. [147][148][149]. We introduce a time scale ∆t(k), which represents a typical time scale for the thermal distribution of H, i.e., the time scale for a distribution function n k ≡ |H k | 2 w k reaches to the equilibrium distribution n eq k = (−1 + exp (− k 2 + gs 2 /T )) −1 . with w k = k 2 + g 2 s 2 . The k dependence of ∆t[k] ∝ g 4 is model dependent and we here assume that ∆t decreases fast enough if k is smaller than T , i.e. the scatterings of IR modes are efficient. Since s is time varying and the thermalization time scale is finite, n k at the time t exhibits the equilibrium distribution at slightly earlier time t − ∆t[k], Thus we can estimate with C being an O(1) numerical coefficient. We note that s is time-dependent and the time derivative in Eq. (52) is non-vanishing. Eq. (52) represents the deviation from the thermal equilibrium, and we approximated the dominant contribution from the integrant around k ∼ T since when k T (k T ) it is suppressed by ∆t (Boltzmann suppressed).
By inserting Eq. (52) into the EOM (45) and multiply both sides by 2ṡ/m Φ , we obtain the evolution equation for the number density n s . Then it is found that the term proportional toṡ in Eq. (52) leads to the effective friction of s and leads to the dissipation of s energy density. We note that by taking a time average (over a few 2π/m Φ ), terms without time derivatives in the EOM (45) are cancelled out. We then arrive aṫ This is smaller than the naïve estimation of s annihilation contribution ∼ n 2 s × σ ss ∼ g 2 n 2 s /(4πm 2 Φ ) since ∆t[T ]m Φ 1. By counting the number of s one may identify the process corresponding to the particle picture (at least in q 1 limit). For example, s 2nṡ term in the EOM should correspond to the scattering of n zero modes of s. As we can see it is suppressed by (gs 2 /T 2 ) n .
For the symmetry breaking system discussed in the main part of this paper, we similarly obtainṅ The leading term, by noting ∆t[T ] ∼ (y 2 t T ) −1 , corresponds to the dissipation term (21).

B Dissipation of (would-be) NG boson condensate
To discuss the dissipation of the produced NG boson or would-be NG boson, whose momentum is much smaller than the temperature and the occupation number is extremely large, we may also apply a similar method to the case of dissipation of s given in App. A. 16 The Lagrangian under consideration is By assuming spherical symmetric distribution of a k and assuming that only | k| = k modes dominate, the Hartree approximation reads 16 The dissipation of QCD axion, which is coupled to the gluon though the anomaly, has been discussed in Ref. [82]. In our present model, we do not necessarily assume such interactions and the dominant source of axion dissipation comes from the interaction with the SM Higgs (55).
where · · · includes a k =k modes or higher order in 1/m 2 Φ . We neglect the other modes again due to the small occupation (note that in our scenario q a m Φ ∼ m Φ H ubble , which means that the NG bosons are soon produced and the spectrum is nearly monochromatic). Then we obtain d dt with n a k [t] = a amp k [t] 2 k/2, t ave is a time scale much longer than k but so short that a amp k can be taken as constant. We note the contribution from the first term of (57) is negligible with large enough t ave since it includes terms of even number ofȧ k . In this case, the integral consists only total derivatives by usingä k [t] ≈ −k 2 a k [t]. The non-vanishing contribution comes from the second term of (57), which has a leading contribution of By assuming again that ∆t[k] is larger at larger t we obtain the integral dominates at around k ∼ T, and d dt Substituting this expression into the equation of motion (56) and performing the t integral by using the explicit form of a k [t] and n a k , we arrive aṫ This gives the dominant contribution to the k modes scattering of the NG bosons. This is equivalent to have a dissipation rate of where we have used ∆t[T ] ∼ (y 2 t T ) −1 .
The same result is also obtained from a diagrammatic approach. The optical theorem tells us that the annihilation cross section of the NG bosons is given by σv aa ∼ (k 0 ) −2 ImM(aa → aa), where k 0 denotes the total incoming NG boson energy and M(aa → aa) the amplitude. Since the SM Higgs obtains large thermal mass, we need to take account of its thermal width for the Higgs propagating in the loop. At the one-loop level, it is evaluated as [78] ImM(aa → aa) ∼ λ 2 where f B (q 0 ) = (e q 0 /T − 1) −1 is the Bose-Einstein distribution and we take the Breit-Wigner form for the spectral density with Γ th being the thermal width of the SM Higgs, which is expected to be (∆t[T ]) −1 . Assuming Γ th m Φ (∼ k 0 ), we obtain the dissipation rate for the NG boson through Γ (a) dis ∼ σv aa n a k and the result is the same as (62).