Gauged Q-ball dark matter through a cosmological first-order phase transition
Siyu Jiang Fa Peng Huang Pyungwon Ko

Revisiting the apparent horizon finding problem with multigrid methods
    Hon-Ka Hui 0000-0002-2815-527X, Lap-Ming Lin 0000-0002-4638-5044
    April 29, 2024
======================================================================



§ INTRODUCTION

Exploring the nature of dark matter (DM) is one of the central issues in (astro)particle physics and cosmology <cit.>. So far, there are no expected signals of conventional DM candidates like Weakly Interacting Massive Particles (WIMPs) in the DM direct detection and collider search experiments <cit.>. Then simple WIMP scenarios are strongly disfavored. 
There are a number of ways to save WIMP scenarios. For example, DM may have substantial couplings only to the 3rd generation fermions <cit.>, or dark sectors may consists 
of two or more stable DM species (see, for example, <cit.>). 
Or one can discard WIMP scenarios and consider other possibilities for DM productions and annihilations or decays in the early Universe. 
This status motivates us to study ultralight or ultra heavy DM candidate (for reviews, see Refs. <cit.>).

Solitons produced in the early Universe are natural candidates of heavy DM 
(see, for example, <cit.> for hidden sector monopole DM accompanied by stable spin-1 vector DM and massless dark radiation, and Ref. <cit.> for hunting for topological DM using atomic clocks).
These solitons are specific field configurations which are classified into two classes, namely, the topological solitons and the non-topological solitions. Recently, as renaissance of the quark nuggets DM proposed by Witten <cit.>, various new ideas on the non-topological solition DM are proposed, where the DM relic density can be produced by the dynamical process of cosmological first-order phase transition (FOPT), such as the Q-ball DM <cit.>. These new mechanisms can naturally avoid the unitarity problem for heavy DM <cit.>.
Dynamical DM mechanisms are specified by the DM penetration behavior into the bubble which depends on the DM mass and bubble wall velocity <cit.>. 

There are extensive discussions on the non-topological solitons in a theory of complex scalar field with global U(1) symmetry, proposed in <cit.> and known as Q-balls <cit.>. And it is natural to study the Q-balls in the gauged case <cit.>, 
by promoting the global U(1) symmetry to the 
local gauge U(1) symmetry. For reviews of U(1) gauged Q-balls, 
see <cit.>. 
Q-balls have been proposed as a potential DM candidate in supersymmetric theories <cit.>. They can also explain the baryon asymmetry of the Universe <cit.>.
The gauged Q-ball DM in supersymmetry model has also
been studied in several papers <cit.>.
It is meaningful to search for other production mechanisms without supersymmetry of Q-ball or gauged Q-ball DM.
In this paper, we study whether the U(1)
gauged Q-balls produced during cosmic phase transition could be a viable DM candidate. If the gauged Q-ball can be stable under certain
circumstances, we still need some mechanism to
(1) produce the charge asymmetry (i.e. locally produce lots of same charge particles to form Q-ball )
(2) and packet the same sign charges in the small size after overcoming the Coulomb repulsive interaction. For the first condition, the primordial 
charge asymmetry could be produced by some early Universe processes such as decays  
of heavier particles. Cosmological FOPT can naturally realize the second condition and can produce phase transition gravitational wave (GW) which can be detected by future GW experiments, such as LISA <cit.>,
TianQin <cit.>, Taiji <cit.>, BBO <cit.>, DECIGO <cit.>, and Ultimate-DECIGO <cit.>.

In this work, for the first time, we study the natural production mechanism of gauged Q-ball DM through a cosmological FOPT. The paper is organised as follows. We describe the basic model that can produce the gauged Q-balls and the numerical solutions of the Q-ball profiles in section <ref>. Basic properties and the stable parameter space of gauged Q-balls are discussed in section <ref>. Thin-wall approximation and the corresponding analytic evaluations are given in section <ref>. Phase transition dynamics in the Standard Model (SM) plus an extra singlet and the relic density of gauged Q-ball DM are elucidated in section <ref>. Signals and constraints of gauged Q-ball DM are given in section <ref>. Concise conclusions and discussions are given in section <ref>.



 


§ GAUGED Q-BALL



 §.§ Friedberg-Lee-Sirlin–Maxwell model

In this work, we adopt the Friedberg-Lee-Sirlin two-component model <cit.>[The Friedberg-Lee-Sirlin two-component model has been reviewed in details in Refs. <cit.>.] plus gauge component, which is called Friedberg-Lee-Sirlin-Maxwell (FLSM) model <cit.>. This model and the corresponding stability of U(1) gauged Q-ball have been discussed in <cit.>.
We begin our discussions with the following Lagrangian density

    ℒ_=(D_μϕ)^†(D^μϕ)+1/2∂_μ h ∂^μ h-1/4Ã_μνÃ^μν-V(ϕ,h)  ,

where the potential V(ϕ,h) reads

    V(ϕ,h) = λ_ϕ h/2 h^2|ϕ|^2+λ_h/4(h^2-v_0^2)^2  .

ϕ and h are the complex scalar field and (real) Higgs field respectively.
D_μ =∂_μ+i g̃Ã_μ and Ã_μν=∂_μÃ_ν-∂_νÃ_μ where Ã_μ is a dark U(1) gauge field and g̃ is the corresponding gauge coupling constant. Ã_μ can be identified as the dark electromagnetic field. We fix the Higgs mass m_h = 125 GeV and vacuum expectation value v_0 = 246  GeV at zero temperature then λ_h = m_h^2/(2v_0^2) ≈ 0.13. 
The complex scalar ϕ gains mass through the portal coupling with the Higgs. 
In the true vacuum, m_ϕ=√(λ_ϕ h/2) v_0.
We assume λ_ϕ h>0, and thus
the Lagrangian density is symmetric under the dark U(1) symmetry which remains unbroken when the Universe goes through the electroweak phase transition.

The local U(1) gauge symmetry leads to the conserved current, 

    J_μ=i(ϕ^†∂_μϕ+2 i g̃Ã_μ|ϕ|^2)  ,

and the corresponding conserved charge,

    Q=∫ d^3 x J^0  .


Once the gauged Q-balls are formed in this FLSM model, one could consider a coherent configuration of ϕ, h, and Ã_μ at a given charge Q. 
The lowest energy state will have no “magnetic field" so the space component Ã_i=0 <cit.>. We assume spherical symmetry for the lowest energy configuration.  Scaling away the physical dimensions, we introduce dimensionless field variables 𝒜, Φ, and ℋ defined in the convention of Ref. <cit.>.

    Ã_t(r)=v_0 g̃/√(2λ_h)𝒜(ρ),   ϕ(r)=v_0/√(2)Φ(ρ) e^-i ω t,    h(r)=v_0 ℋ(ρ)  ,

where ρ≡√(2λ_h)  v_0 r= m_h r. The Lagrangian, with the substitution of the field variables defined above, becomes

    L=     -4 πv_0/√(2λ_h)∫ d ρρ^2[1/2(ν-α^2 𝒜)^2 Φ^2-1/8(ℋ^2-1)^2-k^2/2ℋ^2 Φ^2-1/2(∂_ρΦ)^2. 
       .-1/2(∂_ρℋ)^2+α^2/2(∂_ρ𝒜)^2],

where α≡|g̃|/√(2λ_h), k ≡√(λ_ϕ h)/2√(λ_h)=m_ϕ/m_h, and ν≡ω/√(2λ_h)  v_0.
By varying L with respect to 𝒜, Φ, and ℋ, we find the equations of motion (EoM) for the three fields,

    1/ρ^2∂_ρ(ρ^2 ∂_ρ𝒜)+ (ν-α^2 𝒜) Φ^2=0  ,


    1/ρ^2∂_ρ(ρ^2 ∂_ρΦ)+[(ν-α^2 𝒜)^2-k^2 ℋ^2] Φ=0  ,

and

    1/ρ^2∂_ρ(ρ^2 ∂_ρℋ)-k^2 ℋΦ^2-1/2ℋ(ℋ^2-1)=0  .


The total energy is given by

    E=4 π v_0/√(2λ_h)∫ d ρρ^2 ℰ  ,

where ℰ=α^2/2 (∂_ρ𝒜)^2+1/2(∂_ρΦ)^2+1/2(∂_ρℋ)^2+1/2[(ν-α^2 𝒜)^2+k^2 ℋ^2] Φ^2+1/8(ℋ^2-1)^2. 
And the total charge is given by

    Q=2 π/λ_h∫ d ρρ^2 (ν-α^2 𝒜) Φ^2  .

From Eqs. (<ref>) and (<ref>) we see that

    λ_h Q/2π=-lim _ρ→∞ 4 πρ^2 ∂_ρ𝒜  .

Therefore we get, for large ρ, 𝒜→λ_h Q/2 πρ or Ã_t→g̃Q/4π r. Then Eq. (<ref>) at large ρ becomes  

    1/ρ^2∂_ρ(ρ^2 ∂_ρΦ)-2νg̃^2Q/4πρΦ+(ν^2-k^2) Φ=0  .

It has been shown in Ref. <cit.> for ν<k that this equation at ρ→∞ has the solution of the form,

    Φ(ρ)=C_U e^-√(k^2-ν^2 )ρ  U(1+νg̃^2 Q/4 π√(k^2-ν^2), 2,2 √(k^2-ν^2)ρ)  ,

where C_U is a constant and U(a, b, z)  is the confluent hypergeometric 
function of the second kind. For √(k^2-ν^2)ρ≫ 1 we get

    Φ(ρ) ∼ρ^-1-νg̃^2 Q/4 π√(k^2-ν^2)e^-√(k^2-ν^2)ρ  .

We can see that in the limit g̃→ 0 this form coincides with the nongauged global Q-ball, Φ(ρ) ∼ρ^-1e^-√(k^2-ν^2)ρ. The difference is caused by taking into account the dark electromagnetic potential 𝒜(ρ). For ω=m_ϕ or ν = k, Eq. (<ref>) takes the form of

    1/ρ^2∂_ρ(ρ^2 ∂_ρΦ)-2k g̃^2Q/4πρΦ=0  .

The solution to this equation reads
Φ(ρ)=C_K K_1(√(2 k g̃^2 Q/πρ))/√(ρ),
where C_K is a constant and K_1(z) is the modified Bessel function of the second kind. For large ρ this solution has the form of
Φ(ρ) ∼ρ^-3/4e^-√(2 k g̃^2 Q/πρ) .
This also differs from the global case, in which one expects Φ(ρ) ∼1/ρ for ν=k. For the case ν>k, it can be seen from Eq. (<ref>) that the corresponding solutions for the complex scalar field are oscillatory at ρ→∞, leading to unexpected infinite charge and energy. We think these solutions are unphysical and should be discarded.

It is convenient to write Eq. (<ref>) in the form

    ∂_ρ(ρ^2 ∂_ρ𝒜)=- (ν-α^2 𝒜) Φ^2ρ^2  .

Suppose that 𝒜(0)>ν/α^2. Eq. (<ref>) then implies that ρ^2 ∂_ρ𝒜 is an increasing function of ρ such that ∂_ρ𝒜>0 and 𝒜(ρ)>ν/α^2 for all ρ>0. This possibility is not acceptable, given that ν>0 and 𝒜(ρ) → 0 at ρ→∞. The only acceptable possibility is that 𝒜(0)≤ν/α^2. Then ∂_ρ𝒜<0 and therefore 𝒜(ρ) is a monotonically decreasing function of ρ. We can then say that 𝒜(ρ) obeys the inequalities

    0 ≤α^2 𝒜(∞) ≤α^2 𝒜(ρ) ≤α^2 𝒜(0)≤ν≤  k  .


The energy integral Eq. (<ref>) can be written in a different form. Demanding that the Lagrangian (<ref>) is stationary at ϵ=1 under the rescaling of the form ρ→ϵρ,  leads to the relation dL(ρ→ϵρ)/d ϵ|_ϵ=1=0, this gives the relation:

    3 ∫ d ρρ^2[1/2 (ν-α^2 𝒜)^2 Φ^2-1/8(ℋ^2-1)^2-1/2ℋ^2 Φ^2 k^2] 
    
    	=    ∫ d ρρ^2[1/2(∂_ρΦ)^2+1/2(∂_ρℋ)^2-α^2/2 (∂_ρ𝒜)^2]  .


Substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>), one obtains

    E   = 4 π v_0/√(2λ_h)∫ d ρ[(ν-α^2 𝒜)^2 Φ^2+1/3(∂_ρΦ)^2+1/3(∂_ρℋ)^2+2α^2/3 (∂_ρ𝒜)^2]
       =4 π v_0/√(2λ_h)∫ d ρ[-(ν-α^2 𝒜) ∂_ρ(ρ^2 ∂_ρ𝒜)+ρ^2{1/3(∂_ρΦ)^2+1/3(∂_ρℋ)^2+2α^2/3 (∂_ρ𝒜)^2}] 
       =ω Q+4 π/3v_0/√(2λ_h)∫ d ρρ^2[(∂_ρΦ)^2+(∂_ρℋ)^2-α^2(∂_ρ𝒜)^2]  ,

where in the third line we have integrated by part and use the fact that ρ^2 ∂_ρ𝒜→ -g̃^2 Q/4 πα^2 and 𝒜→ 0 at large ρ. We can see that the energy
for the free field solution when we neglect the variation of Φ and ℋ takes the form

    E_free =m_ϕ Q+“dark electrostatic energy"   .

The “dark electrostatic energy" is roughly proportional to g̃^2 Q^2/R for charges uniformly distributed on scale R.



 §.§ Numerical results of the field configuration, energy, and charge


After qualitative analysis of the gauged Q-ball solution, we begin to numerically solve Eqs. (<ref>), (<ref>), and (<ref>) with the following boundary conditions,

    ∂𝒜/∂ρ=∂Φ/∂ρ=∂ℋ/∂ρ=0  at ρ=0,  𝒜=Φ=0  and ℋ=1  at ρ=∞  .

The first boundary condition is necessary so that the terms 2/ρ(∂𝒜/∂ρ), 2/ρ(∂Φ/∂ρ), and 2/ρ(∂ℋ/∂ρ) do not become singular at ρ=0, and the latter is necessary because the energy density ℰ and charge density (ν-α^2 𝒜) Φ^2 should be integrable 
over the infinite 
spatial volume and the integral should be finite.

It is convenient to introduce the dimensionless quantities for the total energy and charge of the gauged Q-ball,  

    Ẽ   ≡λ_h /2 π m_ϕE=1/k∫_0^∞ d ρρ^2 ℰ  , 
    Q̃   ≡λ_h Q/2 π=∫_0^∞ d ρρ^2 (ν-α^2 𝒜) Φ^2   ,

which can be calculated directly once the numerical solutions of the corresponding differential equations are found. 

We use the relaxation method <cit.> to solve coupled 2nd order ordinary differential equations, Eqs. (<ref>), (<ref>) and (<ref>), with boundary conditions, Eq. (<ref>). The relaxation method solves the boundary value problems by updating the trial functions on the grid in an iterative way. As an example, we fix k=3.5 and we scan over all of the solutions at a given value of α. The results can be seen from figure <ref>. The frequency ν firstly decreases in the direction of the arrow. We call this the “first branch" where the back reaction of gauge field is small. Then the solutions turn on the “second branch" where ν increases and the gauge field dominates. Contrary to the global Q-ball, the parameter ν does not uniquely determine the charge and energy of the gauged Q-ball. For the global case, the energy and charge increase as the ν approaches zero. However, in the case of gauged Q-ball, the ν is replaced by ν-α^2𝒜. Then on the second branch where the gauge field 𝒜 dominates, ν has to increase in order to satisfy ν-α^2 𝒜 > 0.


We choose four specific solutions P1, P2, P3, P4 in figure <ref> for α=1.0 which are marked by the purple triangles and the corresponding numerical profiles of 𝒜, Φ and ℋ are shown in figure <ref>. We can see that as the value of 𝒜(0) of the gauged Q-ball becomes larger, the Higgs field value inside the Q-ball is closer to zero. Actually, when the value of 𝒜(0) becomes larger, the radius, charge and energy also increase, so we can say that the Higgs value is effectively zero inside large gauged Q-balls.

The total charge of gauged Q-ball for different values of gauge coupling α is shown in the left panel of figure <ref>. It can seen that the charge for the gauged Q-ball is also finite at a nonzero α. However, for the global Q-ball the charge is unbounded from above. In order to obtain a gauged Q-ball with relatively large charge, the gauge coupling has to be small enough.





We can define two typical radius for the Higgs field ℋ(ρ) and the Q-ball field Φ(ρ) respectively. The ρ_⋆ is defined by ℋ(ρ_⋆) =1/2 and ρ_b is defined by Φ(ρ_b) = Φ(0)/2. These two radius are shown in the right panel of figure <ref>. It can be seen that ρ_⋆ is generally larger than ρ_b. One may prefer to call ρ_⋆ the radius of the “Higgs ball". Hereafter, we use ρ_⋆ to represent the Q-ball radius which defines the charge and energy of gauged Q-balls.






§ BASIC PROPERTIES OF GAUGED Q-BALLS




 §.§ dE/dQ for gauged Q-balls

It is well known that for non-gauged Q-balls the relations dE/dQ=ω holds. Here we will show that this also holds for gauged Q-balls in FLSM model. From Eq. (<ref>), we have

    dE/dν   = 4 π v_0/√(2λ_h)∫ d ρρ^2 {α^2(∂_ρ𝒜)(∂_ρd 𝒜/dν)+(∂_ρΦ)(∂_ρdΦ/dν)+(∂_ρℋ)(∂_ρd ℋ/dν) .
       +[(1-α^2 d 𝒜/dν)(ν-α^2 𝒜)+k^2 ℋd ℋ/dν] Φ^2+[(ν-α^2 𝒜)^2+k^2 ℋ^2 ] Φd Φ/dν
       . +1/2(ℋ^2-1)ℋd ℋ/dν}  .

After integrating ∂_ρΦ and ∂_ρℋ by part and using Eqs. (<ref>) and (<ref>), we can get

    dE/dν   = 4 π v_0/√(2λ_h)∫ d ρρ^2 [α^2(∂_ρ𝒜)(∂_ρd 𝒜/dν) + (ν-α^2 𝒜)(1-α^2 d 𝒜/dν)Φ^2 +2(ν-α^2 𝒜)^2ΦdΦ/dν] 
       = 4 π v_0/√(2λ_h)∫ d ρρ^2 {(ν-α^2 𝒜)d/dν[(ν-α^2 𝒜)Φ^2]+α^2(∂_ρ𝒜)(∂_ρd𝒜/dν)}
       = √(2λ_h)v_0νdQ/dν +4 π v_0/√(2λ_h)∫ d ρρ^2 {-α^2 𝒜d/dν[(ν-α^2 𝒜)Φ^2]+α^2(∂_ρ𝒜)(∂_ρd𝒜/dν)}_=0  .

In the third line we have integrated ∂_ρ𝒜 by part and used Eq. (<ref>), then the integral vanishes. Finally,

    dE/dω =  dE/dνdν/dω = (√(2λ_h)v_0ωdQ/dω)(1/√(2λ_h)v_0) = ωdQ/dω  ,

which leads to 

    dE/dQ = ω

for dQ/dω≠ 0. The existence of dQ/dω= 0 indicates the locally minimal or locally maximal charge.










 §.§ Stability of gauged Q-balls

The stability of the Q-ball is an important criterion to judge whether it can serve as the DM candidate. Unlike the global Q-balls, the stability of gauged Q-balls is still being discussed. In this subsection, we will systematically analyze four stability criteria of gauged Q-balls and show the viable parameter space of stable gauged Q-balls.


  §.§.§ Quantum mechanical stability


The quantum mechanical stability is satisfied if 

    E<m_ϕ Q	  or  Ẽ/Q̃ < 1  .

This means that the gauged Q-ball is stable against decay to free scalar particles. It should be noted that if the
Q-ball has decay channels into other fundamental scalar particles which have the mass m_i that are smaller than m_ϕ, we need to replace m_ϕ by m_i in Eq. (<ref>). The decay of global Q-balls or gauged Q-balls has been discussed in several works <cit.>. When the effective energy of DM particles inside the gauged Q-balls π/r_⋆∝ ( ω-g̃Ã_t ) is larger than the masses of decay products, the decay process is kinetically allowed. In our case, the gauged Q-balls decay mainly through the process ϕ^†ϕ→ hh. However, as the Q-ball radius r_⋆ is large enough or ( ω-g̃Ã_t ) is small, the gauged Q-balls are stable against decaying into 125 GeV Higgs particles.

The ratio Ẽ/Q̃ is shown in the  figure <ref>. As the frequency ν→ k on the first branch, there is a region of parameter space where Ẽ/Q̃>1. It implies
the existence of a minimal Q-ball charge Q̃_s defined as Ẽ(Q_s)/Q̃_s=1 of the quantum mechanically stable gauged Q-ball.
We can see from figure <ref> that for the global Q-balls the Ẽ/Q̃ decreases with growing Q̃ when Ẽ/Q̃<1. The branch of Ẽ/Q̃>1 corresponds to ν→ k. So the global Q-balls are quantum mechanically stable as ν≪ k and Q̃ > Q̃_s. One would wonder that the gauged Q-ball will destroy the quantum mechanical stability at large charge because the dark electrostatic energy is proportional to Q̃^2 and thus the Ẽ/Q̃ is proportional to Q̃. However, we found the gauged Q-balls is always quantum mechanically stable on the second branch where the gauge field dominates because of the charge of the gauged Q-balls must be finite.



  §.§.§ Stress stability

Now we investigate the effect of the electrostatic repulsion on the stability of the gauged Q-balls. In Ref. <cit.>, the authors pointed out that, just like the hadrons, a necessary condition for stability of the configuration of gauged Q-balls is the balance of the internal forces, called von Laue condition <cit.>

    ∫_0^∞ dr r^2 p(r) = 0   .

Here  p(r) is the radial distribution of the pressure inside the Q-ball, which can be extracted from the energy-momentum tensor by using the following parametrization <cit.>:

    T_ij = (r̂_i r̂_j - 1/3δ_ij)s(r)+δ_ijp(r)  .

s(r) is the traceless part which yields the anisotropy of pressure (shear forces). This kind of stability has also been studied on global Q-balls <cit.>.

One stronger local constraint is that the normal force per unit area
acting on an infinitesimal area element at a distance r, must be directed outward <cit.>,

    F(r) = 2/3s(r)+p(r)>0  .

This is a necessary but not sufficient condition for stability.
By using the rescaled parameters, we have the expressions of p(ρ) and s(ρ) in the FLSM model:

    p(ρ)    = 2λ_hv_0^4[-1/6(∂_ρΦ)^2  -1/6(∂_ρℋ)^2 + α^2/6(∂_ρ𝒜)^2 + 1/2(ν - α^2 𝒜)^2Φ^2 - k^2/2Φ^2ℋ^2 - 1/8(ℋ^2-1)^2 ]  , 
    
    		s(ρ)    =   2λ_hv_0^4[(∂_ρΦ)^2  +(∂_ρℋ)^2 - α^2(∂_ρ𝒜)^2 ]  ,

from which we get 

    F(ρ) = 2λ_hv_0^4[1/2(∂_ρΦ)^2  +1/2(∂_ρℋ)^2 - α^2/2(∂_ρ𝒜)^2 + 1/2(ν - α^2 𝒜)^2Φ^2 - k^2/2Φ^2ℋ^2 - 1/8(ℋ^2-1)^2 ]  .


In figure <ref>, we show the profile of F(ρ) for four 
marked points in figure <ref>. We can see that the F(ρ) has no nodes for P1 and P2 on the first branch. F(ρ) has negative values on the second branch so the gauged Q-balls on the second branch where gauge potential dominates are unstable. It should be noted that the inequality (<ref>) takes a approximation that Q-balls behave as a continuous media , which needs more discussions.

The stress-stability is the most stringent stability criterion in this work. We are not sure if the gauged Q-ball on the second branch would decay into free particles or smaller Q-balls. And it is meaningful to explore that how we can get rid of the strong constraints from stress-stability, in other words, how to get large gauged Q-balls with large gauge coupling. Maybe we can consider an example where two scalar case where the two scalars ϕ and ψ possess opposite charges. If the electrostatic field produced by the two scalars cancels with each other, which implies

    g̃Q_ϕ + g̃'Q_ψ=0  ,

where  g̃ and g̃' are the gauge couplings of ϕ and ψ respectively, this guarantees the electric neutrality of the interior of the gauged Q-balls <cit.>.
The gauged Q-balls will avoid the electric repulsion which leads to the stress-instability even the gauge coupling is large. We expect these Q-balls can also be the DM candidate because they will behave as the ordinary global Q-balls.





  §.§.§ Stability against fission

For global (i.e., nongauged) Q-balls the corresponding stability criterion against fission takes the form

    d^2E/dQ^2<0   .

This clearly leads to E(Q_1)+E(Q_2)>E(Q_1+Q_2) when E(0)=0. However, this may not hold everywhere for the gauged Q-balls due to the presence of the gauge potential. In Ref. <cit.> the authors gave a detailed discussion on the stability against fission of U(1) gauged Q-balls. They pointed out that we could not make any conclusion about the stability
against fission for gauged Q-balls on the second branch. Nevertheless, it has been shown that the gauged Q-balls on the first branch are generally stable because the back-reaction of the gauge field is generally small.



  §.§.§ Classical stability


The problem of classical stability of U(1) gauged Q-balls is discussed in detail in Ref. <cit.>. The classical stability criterion was firstly derived in <cit.> for one-field Q-balls and was discussed for the model with two scalar fields in <cit.>.
The proof of Refs. <cit.> was based on examining the properties of the energy functional of the
system. Instead, the examination of  Ref. <cit.> is based on the Vakhitov-Kolokolov method <cit.> which utilized linearized EoM for the perturbations above the background solution. 

We only consider the spherical perturbations on the gauged Q-ball.
We adopt the following ansatz:

    ϕ(t,r)    = e^-iω tf(r)+e^-iω te^γ t (u(r)+il(r))  , 
    Ã_t(t,r)     = Ã_t(r)+e^γ ta_0(r)  , 
    
    		h(t,r)    = h(r) + e^γ tσ(r)  ,

where f(r), Ã_t(r) and h(r) are the background solutions. u(r), l(r), a_0(r) and σ(r) are the perturbations on the background.
Then we obtain the linearized EoM as below,

    Δ u +(ω -g̃Ã_t)^2u - γ^2 u-2(ω - g̃Ã_t)γ l - 2g̃(ω-g̃Ã_t)fa_0 - Uu -2Su -Yσ = 0  , 
       Δ l+(ω -g̃Ã_t)^2l - γ^2 l+2(ω - g̃Ã_t)γ u -g̃γ fa_0 -Ul =0  , 
       Δ a_0 -2g̃^2f^2a_0 +4g̃(ω-g̃Ã_t)fu-2g̃γ fl=0  , 
       Δσ - γ^2 σ -Wσ - 2Yu=0  .

where Δ = ∑_i=1^3 ∂_i ∂_i is the 3-dim Laplacian operator and

    U(r) = .∂ V/∂ (ϕ^*ϕ)|_ϕ^* ϕ=f^2(r) 
     h=h(r),    S(r) = .∂^2 V/∂ (ϕ^*ϕ)^2|_ϕ^* ϕ=f^2(r) 
     h=h(r) f^2(r)  , 
       W(r)=.1/2∂^2 V/∂ h^2|_ϕ^* ϕ=f^2(r) 
     h=h(r),    Y(r)=.∂^2 V/∂(ϕ^* ϕ) ∂ h|_ϕ^* ϕ=f^2(r) 
     h=h(r) f(r)  .

The boundary conditions are

    u'(0)=l'(0)=a_0'(0)=σ'(0)=0,    u(∞)=l(∞)=a_0(∞)=σ(∞)=0  .


We want to obtain the parameter γ of perturbations u,l,a_0,σ. The classically unstable mode correspond to γ>0. We can do this by using the shooting method in Ref. <cit.>. We introduce four basis solutions Ψ^(i=1,2,3,4) (r)= (u^(i),l^(i),a_0^(i),σ^(i)) which satisfy the Neumann boundary conditions dΨ^(i)/dr|_r=0=0 and Dirichlet boundary conditions Ψ^(1) (0)=(1,0,0,0), Ψ^(2) (0)=(0,1,0,0), Ψ^(3) (0)=(0,0,1,0) and Ψ^(4) (0)=(0,0,0,1). Then we integrate the Eqs. (<ref>) numerically to find the values of Ψ^(i) at large r=r_∞. Now, recall that we are searching for a specific solution
Ψ(r) = c_1 Ψ^(1)(r) +c_2 Ψ^(2)(r)+c_3 Ψ^(3)(r)+c_4 Ψ^(4)(r) which satisfies Ψ(r_∞) = 0. This gives the system of linear equations,

    D̂_c (c_1,c_2,c_3,c_4)^T=0,   D̂_c = (Ψ^(1)T,Ψ^(2)T,Ψ^(3)T,Ψ^(4)T)  .




Equation (<ref>) has nontrivial solutions only if D_c ≡D̂_c = 0. In figure <ref>
we plot Log |D_c| as a function of Log_10γ where we choose α^2=10^-5 and k=3.5. We found that there is no classically unstable mode for ν = 3.45 as there is no solution for γ>0. This also hold for other solutions of 
U(1) gauged Q-balls except for ν→ k on the first branch. We found there is one unstable mode for ν=3.49. Actually, in the limit ν→ k, the contribution of gauged field can be neglected and one expects that the case is similar to the non-gauged Q-ball where dQ/dω>0 indicates there exists classically unstable mode. This is also discussed in  Refs. <cit.>. However, we usually do not have to worry about this because the region already has been excluded by the quantum instability. In Ref. <cit.>, the authors have shown that the gauged Q-balls on the second branch with small gauge coupling are classically unstable with respect to  axisymmetric perturbations. This enhances our confidence that the gauged Q-ball on the second branch is unstable. The region is almost covered by the stress stability criterion.


In summary, the parameter space of gauged Q-balls are shown in figure <ref>. The red line represents the region where the gauged Q-ball is dominated by the gauge field such that it is unstable under stress stability criterion. The green line represents the region where the gauged Q-ball is unstable under quantum stability criterion. We only plot the quantum mechanical stability and the stress stability criterion because they cover the space which is classically unstable and is unstable against fission respectively. The gauged Q-ball is stable only in the region of ν∈ [ν_min,ν_max]. The ν_min corresponds to the maximal charge of gauged Q-balls. In order to form a gauged Q-ball of given charge at given gauge coupling, the charge has to be smaller than the maximal charge. It should also be noted that, if the dark gauge boson of gauged Q-ball kinetically mixes with SM
photons or Z bosons, this would produce distinct experimental signatures which can constrain Q-ball
charge and couplings.




§ THIN-WALL APPROXIMATION

As the radius of gauged Q-ball becomes large, the width of  profile can be neglected and the gauged Q-balls can be depicted by the thin-wall approximation. The Higgs profile can be approximately viewed as a step function, the vacuum value is zero and v_0 inside and outside the Q-ball, respectively. The derivative of the Higgs field only contributes to the surface term of gauged Q-ball which is negligible when the radius of the ball is large.
Then the problems are reduced to those or ones
with two fields, ϕ and Ã_μ. We will discuss the simplified 
piecewise model and show that it behaves closely to the FLSM model. By using the mapping method introduced by Ref. <cit.>, we give some 
semi-analytic results and some analytic evaluations of the maximal charge of 
gauged Q-balls.


 §.§ Piecewise model

If the Higgs field is approximately h(ρ) = v_0Θ(ρ-ρ_⋆), then we can approximately view the complex scalar moving in the piecewise parabolic potential <cit.>. 
The Lagrangian density can be further approximated as

    ℒ_piecewise=(D_μϕ)^†(D^μϕ)-1/4Ã_μνÃ^μν-m_ϕ^2 ϕ^†ϕΘ(1-ϕ^†ϕ/v^2)-m_ϕ^2v^2Θ(ϕ^†ϕ/v^2-1)  ,

where Θ(x) is the Heaviside step function. In our case,

    m_ϕ = v_0 √(λ_ϕ h/2),    v = v_0√(λ_ h/2λ_ϕ h)  .
 
Note that v is chosen so that V(ϕ) is continuous at ϕ^†ϕ = v^2.




This can be understood from the EoM of the Higgs field in the FLSM model,

    h”(r) + 2/rh'(r) + [m_h^2/2 - λ_h h(r)^2 -  λ_ϕ h f(r)^2] h(r) =0  ,

where we use the definition ϕ(r,t) = f(r)e^-iω t and the prime denotes a derivative with respect to r. If the Higgs field 
is approximately a step function, we can neglect the derivatives, then we have

    h^2 ≈m_h^2/2λ_h - λ_ϕ h/λ_hf^2   for   2λ_ϕ h f^2 < m_h^2   , 
    
    		0   for   2λ_ϕ h f^2 > m_h^2   ,

and we could assume f(r) ≈ 0 outside the bubble, for which h^2 ≈m_h^2/2λ_h = v_0^2. Therefore V(ϕ,h) ≈ m_ϕ^2 ϕ^†ϕ outside the Q-ball and V(ϕ,h) ≈λ_h/4v_0^4 inside the Q-ball. The consistency between the piecewise model and the Friedberg-Lee-Sirlin model has been discussed in Ref. <cit.> and the classical stability of gauged Q-balls in piecewise model has been studied in Ref. <cit.>.
The EoM of gauged Q-balls in the piecewise model after the rescaling 
of Eq. (<ref>) are

    1/ρ^2∂_ρ(ρ^2 ∂_ρ𝒜)+ (ν-α^2 𝒜) Φ^2=0  ,


    1/ρ^2∂_ρ(ρ^2 ∂_ρΦ)+[(ν-α^2 𝒜)^2-k^2 Θ(1-4k^2Φ^2)] Φ=0  .

These equations are easier to solve than the FLSM model by using undershooting/overshooting method. The Q-ball energy and charge reads,

    Ẽ =1/k∫_0^∞ d ρρ^2 ℰ,   Q̃ =∫_0^∞ d ρρ^2 (ν-α^2 𝒜) Φ^2  .

Here, ℰ=α^2/2 (∂_ρ𝒜)^2+1/2(∂_ρΦ)^2+1/2[(ν-α^2 𝒜)^2+k^2 Θ(1-4k^2Φ^2)] Φ^2+1/8Θ(4k^2Φ^2-1) which has a different form form the 
FLSM model. 

We solve Eqs. (<ref>) and (<ref>) numerically to get the profiles of 
the complex field and the gauge field. Then we can get the charge and energy of gauged Q-balls by substituting them into Eqs. (<ref>). These numerical results are shown in figure <ref>. We can see that the piecewise model fits 
the FLSM model well when the gauged Q-ball radius is large. The distinctions 
appear at small ρ_⋆ at which the Higgs field value inside is not approximately zero.



 §.§ Mapping gauged Q-balls

We make a further assumption that the complex scalar field is also a step function, Φ(ρ) = Φ_0(1-Θ(ρ-ρ_b)), where we denote Φ_0=Φ(0) and ρ_b is defined by Φ(ρ_b)=Φ_0/2. Then the profile of gauged field is <cit.>

    𝒜(ρ) = ν/α^2
    		1 - sinh(αΦ_0 ρ)/cosh(αΦ_0ρ_b)αΦ_0 ρ,    ρ<ρ_b 
    αΦ_0ρ_b-tanh(αΦ_0 ρ_b)/αΦ_0 ρ,    ρ>ρ_b  .

The Q-ball radius is defined by Φ(ρ_⋆) = 1/2k.



In Ref. <cit.>, the authors propose a mapping between the gauged Q-ball and the global Q-ball. Specifically,

    ν = ν_g αΦ_0 ρ_b (αΦ_0 ρ_b)  ,

where ν_g is the value of frequency for the global Q-ball with the same ρ_b. Then the profile of 𝒜 is given by Eq. (<ref>). This relation holds even for cases beyond the thin-wall approximation.  Interestingly, the global cases in the piecewise model have analytic solutions <cit.>:

    Φ_g(ρ) = 
    	1/2kρ_⋆sin(ν_g ρ)/ρsin(ν_g ρ_⋆),    ρ<ρ_⋆
    1/2kρ_⋆ e^-√(k^2-ν_g^2)ρ/ρ e^-√(k^2-ν_g^2)ρ_⋆,    ρ>ρ_⋆  .

This gives us Φ_0 = 1/2kν_gρ_⋆/sin(ν_g ρ_⋆) and 2sin(ν_g ρ_b) = ν_g ρ_b. Then we have ν_g ρ_b = C_1 where C_1 ≈ 1.89549[The factor C_1 is close to the result of Ref. <cit.> where C_1 ≈ 2.08 by using the definition Φ”(ρ_b)=0.]. The Q-ball radius is defined as

    ρ_⋆(ν_g) = 1/ν_g(π - arctan(ν_g/√(k^2-ν_g^2)))  .


If we use Φ(ρ) = Φ_0(1-Θ(ρ-ρ_b)) and 𝒜(ρ) from Eq. (<ref>), we have further semi-analytic results for charge:

    Q̃ = ν_g ρ_b/α^2(αΦ_0ρ_b(αΦ_0 ρ_b)-1) = ρ_b/α^2(ν-ν_g)  .

In the limit α→ 0 and αΦ_0ρ_b → 0, because x x ∼ 1+x^2/3 +𝒪 (x^3) for  x → 0. Then we have Q̃∼ν_g Φ_0^2 ρ_b^3. And in this case, when ν_g → 0, from Eq. (<ref>) we have ν_g ρ_⋆≃π and sin(ν_g ρ_⋆) ≃ν_g/k, which lead to

    Φ_0 ≃1/2kπ/ν_g/k = π/2ν_g  .

From 2sin(ν_g ρ_b) = ν_g ρ_b, we have 
ν_g ∼C_1/ρ_b, and

    Q̃∝ρ_b^4  ,

which is consistent with the global Q-ball case.
Using the same procedure that derives Eq. (<ref>), we have analytic results for total energy:

    Ẽ   = νQ̃/k + 1/3k∫ d ρρ^2[(∂_ρΦ)^2-α^2(∂_ρ𝒜)^2]  
       ≃νQ̃/k + Φ_0ρ_b^2/12k - 1/3kν^2[αΦ_0 ρ_b(2+sech^2(αΦ_0 ρ_b)) - 3tanh (αΦ_0 ρ_b)] /2α^3Φ_0  .

The second term comes from the integration over discontinuous (∂_ρΦ)^2 by using the approximation of energy conservation <cit.>. 









In the first limit of αΦ_0 ρ_b → 0 and large ρ_b, the last term of Eq. (<ref>) vanishes. We have ν≃ν_g, Φ_0 ≃π/2ν_g≃ρ_⋆/2 and ρ_⋆≃ρ_b, then

    Ẽ≃πQ̃/kρ_b + ρ_b^3/24k  ,

which is just the energy of global Q-ball.

In the opposite limit αΦ_0 ρ_b →∞, because 𝒜→Q̃/ρ, then from Eq. (<ref>),

    ν = α^2 Q̃/ρ_b(1-tanh(αΦ_0 ρ_b)/αΦ_0 ρ_b)^-1≃α^2 Q̃/ρ_b  .

Thus from Eq. (<ref>) the energy of gauged Q-ball is

    Ẽ≃νQ̃/k + Φ_0ρ_b^2/12k - 1/3kν^2 ρ_b/α^2≃2α^2 Q̃^2/3kρ_b +  Φ_0ρ_b^2/12k  .
 
The first term is the Coulomb energy and the second term is the potential energy difference between inside and outside of the Q-ball.
The second term is proportional to ρ_b^2 because in this case the Compton wavelength of the gauge field 1/g̃v_0Φ_0 inside the Q-ball is much smaller than Q-ball radius, r_b = ρ_b/m_h.
So the Q ball is superconducting. The potential
energy is therefore zero inside as well as outside of the Q
ball and is nonzero only in the shell around ρ_b <cit.>.

For a given ρ_⋆, we can solve the Eq. (<ref>) to get the corresponding ν_g. After using Φ_0 = 1/2kν_gρ_⋆/sin(ν_g ρ_⋆), 2sin(ν_g ρ_b) = ν_g ρ_b and Eq. (<ref>) we can get Φ_0, ρ_b and ν of gauged Q-balls. Substitute them into
Eqs. (<ref>) and (<ref>) to get the charge and energy of gauged Q-balls. These semi-analytic results are shown by the red lines in figure <ref>. We find the mapping works well for ρ_⋆≫ 1 on the first branch where the profiles of scalar fields ϕ and h can be safely viewed as step function. The semi-analytic results of energy doesn't work well for the second branch because the Φ_0 ≠1/2kν_gρ_⋆/sin(ν_g ρ_⋆) which is the value in the global case. The value of Φ_0 is lower because the gauge potential dominates. This can be seen from figure <ref>. At the turning point between first branch and second branch, the discrepancy between semi-analytic results and numerical results is about a factor of 𝒪(1).



 §.§ Maximal charge and energy of gauged Q-balls: analytic approximations


We give analytic evaluations of maximal charge and maximal energy of gauged Q-balls for a given gauge coupling. Because the gauged Q-balls should be unstable under stress stability criterion on the second branch where the gauge potential dominates, the maximal charge is approximately defined by . d ν/dρ_⋆|_ρ_⋆ = ρ_max=0. In the limit at large Q-ball with ν_g → 0, we have ρ_b = C_1/ν_g and Φ_0 ≃ρ_⋆/2. 
Then, from (<ref>) we get

    ν = C_1 αρ_⋆/2(C_1 αρ_⋆/2ν_g) = C_1 αρ_⋆/2(C_1 αρ_⋆^2/2π)  ,

where we used ρ_⋆≃π/ν_g when ν_g → 0. Then we have when the charge is maximal . d ν/dρ_⋆|_ρ_⋆ = ρ_max=0,

    sinh(C_1 αρ_max^2/2π) 	cosh(C_1 αρ_max^2/2π) = C_1 αρ_max^2/π  ,

then we get C_1 αρ_max^2/2π = C_2 with C_2 ≈ 1.08866. So in this case αΦ_0 ρ_b ≃ C_2 which is somewhere in the middle of Eq. (<ref>) and Eq. (<ref>). We can also get the minimal frequency 

    ν_min = C_1 αρ_max/2(C_1 αρ_max^2/2π)= √(π C_1 C_2 α/2)(C_2) ≈ 2.26 √(α)  ,

and the corresponding frequency for the global case ν_gmin = π/ρ_max = √(π C_1 α/2C_2). We can see that ν_min/ν_gmin = C_2(C_2) ≈ 1.367 which is independent of the Q-ball size. This implies that as the solutions on the second branch are unstable, the gauged Q-ball lives on the first branch and is close to the global case. 



Finally, we get the maximal charge, 

    Q̃_max = C_1/α^2(C_2 C_2-1) ≈ 0.7 α^-2  ,

from which we can see the charge is unbounded from above in the global case where α=0.
The maximal energy reads 

    Ẽ_max   = ν_minQ̃_max/k + C_1^2 ρ_max^3/24π^2k-1/3kν_min^2[C_2(2+sech^2C_2)-3tanh C_2]/α^3 ρ_max
       = π C_2  C_2 Q̃_max/kρ_max+ C_1^2 ρ_max^3/24π^2k-1/3kπ C_1C_2 ^2 C_2[C_2(2+sech^2C_2)-3tanh C_2]/2α^2 ρ_max  .

By using ν_min= √(π C_1 C_2 α/2) C_2 and ρ_max= √(2π C_2/C_1 α), we have

    Ẽ_max   =1/kα^3/2√(π C_1 C_2/72)[-3C_1  C_2+C_2(1/π+4C_1+3C_1 csch^2 C_2)]
       ≈ 1.51 k^-1α^-3/2  .

This also gives us Ẽ_max∝ρ_max^3 ∝Q̃_max^3/4 which is expected in the global case. 

The analytic results of Q_max = 2π/λ_hQ̃_max and E_max = 2π m_ϕ/λ_hẼ_max are shown in terms of the red lines in figure <ref> and we can see that the maximal charge Q_max fits well with the semi-analytic and numerical results.  However, the analytic and semi-analytic E_max is about 5-6 times larger than the numerical results due to the uncertainties of values of Φ_0.



§ GAUGED Q-BALL DM FROM ELECTROWEAK FOPT

In the above discussions, we have shown that the gauged Q-ball could be stable and hence can make a possible DM candidate. In this section, we begin to discuss the detailed production mechanism of 
the gauged Q-ball DM that is formed during  the electroweak FOPT in the early Universe. We consider the minimal Higss extended model with a singlet scalar field which could trigger a FOPT <cit.>. The phase transition dynamics can also be modified by introducing some new freedoms beyond the standard model.



 §.§ Electroweak FOPT


The electroweak FOPT dynamics is determined by the finite temperature effective potential V_eff (h, T) where h is the real component of the SM Higgs doublet as defined in Eq. (<ref>),


    V_eff (h, T) ≡ V_tree (h)+V_CW(h)+V_T(h, T)  .

The first term V_tree (h)=λ_h(h^2-v_0^2)^2 / 4 is the tree-level SM Higgs potential. V_CW(h) is the one-loop quantum correction to the effective potential, i.e., Coleman-Weinberg potential. Using the on-shell renormalization scheme, we have

    V_CW(h)=∑_i(-1)^F_ig_i/64 π^2[m_i^4(h)(logm_i^2(h)/m_i^2(v_0)-3/2)+2 m_i^2(h) m_i^2(v_0)],

where g_i is the degree of freedom for each particle, F_i=1(0) for fermions(bosons), m_i(h) are masses for =t, W, Z, h, ϕ. The finite-temperature correction term is given by

    V_T(h,T)=∑_i(-1)^F_ig_i T^4/2 π^2∫_0^∞ d x x^2 log[1 ∓ e^(-√(x^2+(m_i^2(h)+Π_i) / T^2))],

where the integral with “-/+ " sign denotes the  contribution of bosons/fermions. Π_i
is thermal masses of species i. Here, we use the daisy resummation scheme proposed by Dolan and Jackiw <cit.>.
It is worth noticing that only the scalar fields and the longitudinal components of the gauge fields have nonzero Π_i. 
For the scalar fields

    Π_h    = (λ_h/2+λ_ϕ h/12 + 3g^2+g'^2/16 + y_t^2/4)T^2,   Π_ϕ   = λ_ϕ h/6T^2  ,

where g and g' are the gauge coupling of SU(2)_L and U(1)_Y, respectively.
For the longitudinal components of the gauge bosons, we have

    Π_W_L = Π_Z_L = 11g^2/6T^2,   Π_B_L =  11g'^2/6T^2  .

Hence, for the longitudinal components of the gauge bosons, their physical masses are eigenvalues of the following matrix

    M_L^2 = ([ m_1^2 + Π_W_L             0             0             0;             0 m_1^2 + Π_W_L             0             0;             0             0 m_1^2 + Π_Z_L        m_12^2;             0             0        m_12^2 m_2^2 + Π_B_L;               ])  ,

where m_1^2=g^2 h^2/4, m_2^2=g'^2 h^2/4 and m_12^2=-gg' h^2/4.


Requiring the ordinary electroweak vacuum with h=v_0=246 GeV as the global vacuum at T=0 or V_eff (v_0, 0)<V_eff (0,0) leads to

    λ_ϕ h≲4 √(2)π m_h/v_0≈ 9.0 .


The phase transition is the process of symmetry breaking in the early Universe. Through a process of bubble nucleation, growth and merger, the Universe transits from a metastable state into a stable state. The critical temperature T_c is defined by the time when the two minima of effective potential is degenerate, V_eff(v(T_c),T_c)= V_eff(0,T_c) with v(T_c) being the vacuum value in the true vacuum at T = T_c. Bubbles begin to nucleate when the temperature drops to the nucleation temperature T_n.
The nucleation rate of bubbles is given by 


    Γ(T) ≈ T^4(S_3(T)/2π T)^3/2 e^-S_3(T) / T  ,

with S_3(T) being the action of the 𝒪(3) symmetric bounce solution <cit.>.

The nucleation temperature T_n is typically defined by

    Γ(T_n) H^-4(T_n) ≈ 1   ,

where H(T) is the Hubble expansion rate,

    H^2(T)=8 π/3 M_pl^2(π^2/30 g_⋆ T^4+Δ V_eff(T))  .

where M_pl=1.22 × 10^19  GeV  is the Planck mass and g_⋆ is the number of relativistic degrees of freedom at temperature T. Δ V_eff(T) is the potential energy difference between false and true vacuum
Δ V_eff(T)=V_eff(0, T)-V_eff(v(T), T).
The potential difference between the inside and outside the bubbles will cause the bubbles expanding in the Universe so that the volume of the false vacuum diminishes with time. The probability of finding a point in the false vacuum reads,

    p(T)=e^-I(T)  ,

where I(T) is the fraction of vacuum converted to the true vacuum <cit.>,

    I(T)=4 π/3∫_T^T_c d T^'Γ(T^')/T^' 4 H(T^')[∫_T^T^' d T̃v_w/H(T̃)]^3  .

In the radiation dominated Universe <cit.>,

    I(T) = 675M_pl^4/4π^5g_⋆^2∫_T^T_cdT^'Γ(T^')/T^' 6(1/T - 1/T^')^3  .

The percolation temperature T_p, is defined by I(T_p) = 0.34 <cit.>. This means that 34% of the false vacuum has been converted to the
true vacuum at T_p. The percolation temperature T_p is also the temperature at which the GW is produced  from the FOPT <cit.>.

We use the following definition of the phase transition strength:

    α(T) ≡(1-T/4∂/∂ T) Δ V_eff/ρ_r,

where ρ_r=π^2g_⋆T^4/30 represents radiation energy density and Δ V_eff is the potential difference between the false and the true vacua. The inverse time duration β at fixed temperature T_i is defined as

    β=H(T) .T d/d T(S_3/T)|_T=T_i  .

We use CosmoTransitions <cit.> to calculate the phase transition dynamics. In the left panel of figure <ref>, we show the three typical temperatures of the FOPT process in the minimal SM plus singlet model. And in the right panel of figure <ref>, we show the wash-out parameter v(T)/T at these different temperatures.





 §.§ Bubble wall filtering during FOPT

As the particles gain mass inside the bubble, due to the energy conservation, only high-energy particles can pass through the bubbles walls and the others are reflected. The condition of penetration in the bubble wall frame reads <cit.>:

    p_z^w>√(Δ m_i^2)  ,

where p_z^w is the particle z-direction momentum in the bubble wall frame, Δ m_i^2=(m_i^in)^2-(m_i^out)^2 is the mass difference between the true and false vacuum where m_i^in is the particle mass in the true vacuum and m_i^out is the mass in the false vacuum <cit.> which we will set to zero in this work.

The particle flux coming from the false vacuum per unit area and unit time can be written as <cit.>

    J_i^w=g_i ∫d^3 p^w/(2 π)^3p_z^w/√((p^w)^2+(m_i^out)^2) f_i^ eqΘ(p_z^w-√(Δ m^2))   ,

where g_i is degrees of freedom of the particle. p^w is the magnitude of the three-momentum of the particles. f_i^ eq is the equilibrium distribution outside the bubble in the bubble wall frame,

    f_i^ eq=1/e^γ_w(√((p^w)^2+(m_i^out)^2)- v_wp_z^w) / T∓ 1  ,

where ∓ is for bosons and fermions respectively. v_w is the bubble wall velocity and γ_w=1/√(1-v_w^2) is its Lorentz boost factor. 
The particle number density inside the bubble n_i^in in the bubble center frame can be written as 

    n_i^in =J_i^w/γ_w   v_w  .

Assuming the particles are massless in the false vacuum, 
 we can integrate Eq. (<ref>) analytically and get <cit.>

    n_i^in ≃g_i T_^3/γ_w v_w(γ_w(1- v_w) m_i^in / T_+1/4 π^2  γ_w^3(1- v_w)^2) e^-γ_w(1-  v_w) m_i^in/T_  ,

where we have used Maxwell-Boltzmann approximation of DM distribution.
One can see that as v_w→ 1,  Eq. (<ref>) approaches n_i^out=g_i T^3 / π^2, which is approximately the equilibrium number density for Boltzmann distribution outside the bubble. 
The fraction of particles i that are trapped into the false vacuum is defined by

    F_i^trap = 1-n_i^in/n_i^out  .

In our case, complex scalars ϕ and ϕ^† are trapped into the false vacuum due to the filtering effect. The symmetric part would annihilate away in terms of the process ϕ+ϕ^†→ h+h, then the asymmetric part survives and composes the charge of gauged Q-balls. 
It can be easily seen that the penetrated particle number density is sensitive to the bubble wall velocity as it appears in the exponent of Eq. (<ref>). The precise calculation of bubble wall velocity <cit.> is beyond the scope of this work and we set it as a free parameter.




 §.§ Charge of Q-ball in the electroweak FOPT


We can
define T_⋆ as the temperature at which the false vacuum or old phase remnants can still form an 
infinitely connected “cluster", just like the definition of percolation temperature <cit.>. The T_⋆ satisfies p(T_⋆)=1-p(T_p)=0.29 which corresponds to I(T_⋆)=1.24.
T_⋆ is the temperature when Q-balls start to form. Below the temperature T_⋆, the false vacuum remnants formed during FOPT may further fragment into smaller pieces. Ultimately, these pieces would shrink into Q-balls if there exists a non-zero primordial charge asymmetry. We can define the critical radius,  r_c, at which the remnant shrinks to an insignificant size before another true vacuum bubble form within it <cit.>. This means <cit.>

    Γ(T_*) (4 π/3 r_c^3) Δ t ∼ 1  .

where Δ t=r_c/v_w is the time cost for shrinking. The number density of the remnants  n_Q^⋆ can be expressed as: 

    n_Q^⋆≃ 0.29 (3/4 π)^1 / 4(Γ(T_⋆)/v_w)^3 / 4  ,

since the condition n_Q^⋆(4 π/3 r_c^3)=p(T_⋆)≃ 0.29. 

The formation of Q-balls requires a nonzero conserved primordial charge which comes from the primordial DM asymmetry η_ϕ=(n_ϕ-n_ϕ^†)/s. The entropy density s_=2π^2g_⋆T_^3/45. 
If the DM asymmetry is produced by thermal freeze-out in the early Universe, it is bounded from above by the equilibrium value,

    η_ϕ≲η_ϕ^eq =  n_ϕ^eq(T)/s(T)≃ 5.1 × 10^-3×(108.75/g_⋆)  ,

where we have used n_ϕ^eq(T) = 2ζ(3) T^3/π^2 with  ζ (3) =1.20206 being the value of Riemann zeta function ζ(s) at s=3. In order to overcome this constraint, we assume the DM is produced by some non-thermal processes like decay. In this work, we do not specify the origin of primordial charge of the complex scalar ϕ.
In the early Universe at higher temperature, new physical processes beyond the standard model may have occurred.  The would-be Q-ball DM particles ϕ
may have their own conserved charge and be created in asymmetric
decays of heavier particles <cit.>. For example, heavy Majorana neutrino could decay into a light fermion and a scalar, like N →χ̅ + ϕ and N →χ + ϕ^† where χ is a fermion <cit.>. Assuming the process is CP-violating, the decay rates of these two channels differentiate from each other at loop level (when the ϕ is the standard model Higgs and χ is leptons, this is the process of leptogenesis.). So the asymmetry between ϕ and ϕ^† appears and is retained until the phase transition in this work. The large DM asymmetry can be discussed in a similar manner to the large lepton asymmetry in the leptogenesis scenarios. Recent  measurement of ^4He  abundance coming from the EMPRESS experiment suggests a large degeneracy parameter of the electron neutrino <cit.>,

    ξ_e=0.05_-0.02^+0.03  .

Since the neutrino oscillations among three flavors, the neutrinos with three flavors have the same amount of asymmetry, and then the total lepton asymmetry reads,

    η_L ≡ 3n_ν_e-n_ν̅_e/s = 3T^3ξ_e/6/2π^2g_BBNT^3/45≃ 5.3 × 10^-3  ,

where g_BBN=10.75 is the relativistic degree of freedom at the epoch of big bang nucleosynthesis. The large lepton asymmetry may comes from the low scale leptogenesis <cit.>, Affleck-Dine mechanism <cit.> or L-ball decay <cit.>.

In a remnant, the trapped Q-charge is given by 
Q_⋆=F_ϕ^trapη_ϕ s_⋆/n_Q^⋆. In figure <ref>, we show the charge of the gauged Q-ball DM formed during the FOPT for different values of bubble wall velocities. We have chosen η_ϕ=η_L.  When λ_ϕ h is larger, both 
the Γ(T_⋆) and the Q-ball number density are suppressed, so that the charge is larger at a given η_ϕ. 



Since n_Q / s and Q does not change in the adiabatic Universe, at present they are

    n_Q=n_Q^⋆/s_⋆ s_0,    Q=Q_⋆  ,

with s_0=2891.2 cm^-3 being the entropy density in current time <cit.>.



We take the approximation T_⋆≈ T_p then the bounce action can be approximated written as <cit.>

    S_3(T_⋆)/T_⋆≃ 131 - 4ln(T_⋆/100 GeV) - 4ln(β/H_⋆/100) + 3ln v_w -2ln(g_⋆/100)  ,

where β/H_⋆=. β/H |_T=T_⋆. By approximately using Γ(T_⋆) ≈ T_⋆^4 e^-S_3(T_⋆)/T_⋆, then we get

    Q_ = 5.35× 10^36 F_ϕ^trap(v_w/0.01)^3(η_ϕ/η_L)(100 GeV/T_⋆)^3(100/β/H_⋆)^3(100/g_⋆)^1/2  .

In order to get stable gauged Q-ball with a given value of charge, we must impose Q_max>Q_, so the gauge coupling g̃ has to satisfy

    g̃<1.28× 10^-18(1/F_ϕ^trap)^1/2(0.01/v_w)^3/2(η_L/η_ϕ)^1/2(T_⋆/100 GeV)^3/2(β/H_⋆/100)^3/2(g_⋆/100)^1/4  .






 §.§ Relic density of gauged Q-ball DM


The DM relic density can also have the contribution from the standard freeze-out, through the process ϕ +ϕ^†↔ h+h. The relic abundance reads,

    Ω_freeze-outh_100^2 ≈2.58× 10^-10 GeV^-2/⟨σ_anni v_rel⟩  ,

where ⟨σ_anni v_rel⟩ is the annihilation cross section. h_100=H_0 /(100  km·s^-1·Mpc^-1)=0.67 where H_0 is the Hubble constant today <cit.>. The cross section of process ϕ +ϕ^†↔ h+h reads λ_ϕ h^2/64π m_ϕ^2(1-m_h^2/m_ϕ^2)^1/2 <cit.>. In our parameter space, where the λ_ϕ h is around 7, the relic abundance from freeze-out is approximately Ω_freeze-outh_100^2 ≈ 2.35× 10^-4.
So we can omit the DM produced from thermal freeze-out.

The DM relic density also receives the contribution from penetrated asymmetric components of DM particles which is given by the excess of ϕ over ϕ^†,

    Ω_asymmetrich_100^2=(1-F_ϕ^trap)η_ϕs_0 m_ϕ  .


The relic density of Q-balls at present is

    Ω_Q h_100^2=n_Q E_Q/ρ_c h_100^2  ,

where ρ_c=3 H_0^2 M_pl^2 /(8 π) is the critical energy density. We have found that, the gauged Q-ball is generally a mixed state of Eq. (<ref>) and Eq. (<ref>) as 0≤α B_0 ρ_b ≤ C_2. So we can write down the energy of gauged Q-ball at a given charge, 

    E_Q ≃π Q/r_⋆ + 4π/3r_⋆^3V_0 + 3/20πg̃^2Q^2/r_⋆  ,

where V_0=λ_h/4v_0^4 is the potential difference between inside and outside of the gauged Q-balls at zero temperature. The first term on the right-side is the zero-point energy of the scalar particles, the second term is the vacuum volume energy inside the Q-ball and the third term is the  electrostatic self energy. By minimizing this expression respect to r_⋆, we obtain,

    E_Q = 4π/3(4V_0)^1/4Q^3/4(1+3g̃^2 Q/20π^2)^3/4  .

We can see that in the limit of zero gauge coupling

    E_Q =4π/3(4V_0)^1/4Q^3/4=4π/3Q^3/4λ_h^1/4v_0  ,

which is just the energy of global Q-ball. By using n_Q=n_Q^⋆s_0/s_⋆ and Q=Q_⋆=F_ϕ^trapη_ϕ s_⋆/n_Q^⋆, we finally arrive at the expression:

    Ω_Q h_100^2 
       ≃ 2.81 ×(s_0h_100^2/ρ_c) (Γ(T_⋆)/v_w)^3/16s_⋆^-1/4(F_ϕ^trapη_ϕ)^3/4λ_h^1/4v_0(1+108^1/4g̃^2F_ϕ^trapη_ϕs_⋆ v_w^3/4/5.4π^7/4Γ(T_⋆)^3/4)  ,

the Γ(T_⋆) can also be expanded by using Γ(T_⋆) ≈ T_⋆^4 e^-S_3(T_⋆)/T_⋆ and Eq. (<ref>), but we keep the expression here to give more accurate results.

Although the expression hitherto is general, we apply these to the minimal SM plus singlet model. We choose λ_ϕ h=6.8 for which T_n=71.65 GeV and T_p= 68.9 GeV. The number density of gauged Q-balls at production is n_Q^⋆≃ 3.0× 10^-31 GeV^-3. In this case, the value of v(T_p)/T_p is approximately 3.5 and there are still 50% of the DM particles trapping inside the false vacuum even at v_w=0.6. But it should be noted that in this case, the contribution from penetrated asymmetric DM will dominate, as can be seen in Eq. (<ref>). This can be avoided in two ways. One way is to increase the phase transition strength and the corresponding v(T)/T so there are little DM particles penetrating into the true vacuum, 
resulting F_ϕ^trap≈ 1. This can be achieved by introducing new freedoms beyond the standard model or considering dark FOPT instead of electroweak FOPT. The other way is model dependent: 
one can introduce new decay channels that penetrated ϕ could decay into 
dark radiation or SM leptons which can account for the lepton asymmetry. The decay process does not destroy the stability of Q-balls as long as π/r_⋆ < m_d where m_d is the mass of decay products. In this work we focus on the gauged Q-ball DM so we do not consider the penetrated asymmetric DM in details.

A strong FOPT requires large portal coupling λ_ϕ h=6.8, and the validity of the perturbative analysis should be roughly smaller than 10 <cit.>. The condition of a strong FOPT leads to large deviation of triple Higgs coupling, which might be detected by the loop-induced e^+ e^-→ hZ process <cit.> at future lepton colliders, such as FCC-ee, CEPC, and ILC.





We show the gauged Q-ball DM relic density in figure <ref>. The colored star represents the value for the gauged Q-ball with maximal charge. The gray region represents the region where the DM is overproduced. We can see that due to the finiteness of the charge of gauged Q-balls, the gauged Q-balls can explain the whole DM at v_w=0.01 only when the rescaled gauge coupling α≲ 10^-16. The DM relic density is slightly enhanced due to the extra electrostatic energy.

In table <ref> we choose four benchmark points that satisfy the correct DM relic density and show the corresponding F_ϕ^trap and η_ϕ/η_L. Besides, we also show the
 fractional change in Zh production relative to the SM prediction at one loop, δσ_Zh. We can see that the value of the required DM asymmetry η_ϕ is close to the lepton asymmetry η_L. This prompts us to speculate that they have the same origin. Actually, if the DM asymmetry comes from the process N →χ̅ + ϕ and N →χ + ϕ^†. We could assume the dark fermion χ is long-lived and decay suddenly into leptons after electroweak phase transition in order to avoid the electroweak sphaleron process. Then we expect the η_ϕ and η_L is at the same order. The detailed model building is beyond the scope of this work and we leave this in our future studies.






§ CONSTRAINTS AND DETECTION OF GAUGED Q-BALL DM



 §.§ Direct detection and astronomical constraints

We keep the gauged Q-ball DM mass and radius as free parameters to discuss its detection potential or constraints. The variation of gauged Q-ball DM mass and radius can be easily realize by varying the primordial DM asymmetry or the phase transition dynamics. The combined constraints of gauged Q-ball are shown in figure <ref>. As we have discussed in the previous sections, the size of the gauged Q-ball is restricted to be finite. The maximal charges or the maximal mass of gauged Q-ball DM at given gauge coupling are marked by the stars in figure <ref>.
The gray region denotes the constraints from cosmic microwave
background (CMB) which is affected by DM-baryon scattering <cit.>. DM cross sections have been probed by a variety of shallow and deep
underground DM or repurposed experiments: XENON1T (orange) <cit.>, Mica (purple) <cit.>, Ohya (green) <cit.>. The Q-ball DM would transfer
energy with other objects primarily through elastic scattering, approximately utilizing their geometric cross-section. When these cross-sections are sufficiently large, the resultant linear energy transfer might trigger observable phenomena. For example, if a
macroscopic DM traverses compact objects such as white dwarfs or neutron stars and triggers
thermonuclear runaway, this could potentially lead to a type IA supernova or superburst, respectively <cit.>. The brown region represents constraints from superbursts in neutron stars <cit.> and the blue region from white dwarf becoming supernovae <cit.>. Combined all constraints in figure <ref>, one can see that the gauged Q-balls could be the DM candidate in a wide region of parameter space.





 §.§ Phase transition GW


The phase transition GW spectra comes from three sources of a strong FOPT: bubble collision, sound wave, and turbulence.

	
  * Bubble collision

	The formula of phase transition GW from bubble collisions <cit.> at the percolation temperature T_p , reads:
	
    Ω_coh_100^2 ≃ 1.67 × 10^-5(H_p/β)^2(κ_h α_p/1+α_p)^2(100/g_⋆)^1 / 30.11 v_w^3/0.42+v_w^23.8(f / f_co)^2.8/1+2.8(f / f_co)^3.8  ,

	where β/H_p and α_p are defined at percolation temperature, respectively. κ_h represents the fraction of vacuum energy converted into the scalar field's gradient energy. f_co is the peak frequency of bubble collision:
	
    f_co≃ 1.65 × 10^-5 Hz(β/H_p )(0.62/1.8-0.1 v_w+v_w^2)(T_p/100 GeV)(g_⋆/100)^1 / 6  .

	
  * Sound wave

The contribution from sound waves could be more significant. The formula GW spectrum from sound waves is <cit.>:
	
    Ω_swh_100^2 ≃    2.65 × 10^-6Υ_sw(H_p/β)(κ_v α_p/1+α_p)^2(100/g_⋆)^1 / 3 
    			  v_w  (f / f_sw)^3  (7/4+3(f / f_sw)^2)^7 / 2  ,

where κ_v reflects the fraction of vacuum energy that transfers into the fluid's bulk motion. The peak frequency of sound waves processes is:

    f_sw≃ 1.9 × 10^-5 Hz1/v_w(β/H_p)(T_p/100 GeV)(g_⋆/100)^1 / 6  .

Υ_sw is the suppression factor of the short period of the sound wave,


    Υ_sw = (1-1/√(1+2τ_swH_p))  ,

where

    τ_swH_p ≈ (8π)^1/3v_w (H_p/β)/√(3κ_vα_p/(4+4α_p))  .


  * Turbulence  

The formula of the GW spectrum from turbulence is <cit.>:

    Ω_turbh_100^2 ≃ 3.35 × 10^-4(H_p v_w/β)(κ_turbα_p/1+α_p)^3 / 2(100/g_⋆)^1 / 3 (f / f_turb )^3/(1+f / f_turb )^11 / 3(1+8 π f / H_p)  ,
	
note that H_p is the Hubble rate at T_p:
 
    H_p=1.65 × 10^-5 Hz(T_p/100 GeV)(g_⋆/100)^1/6  ,

and the peak frequency of turbulence processes is :
	
    f_turb ≃ 2.7 × 10^-5 Hz1/v_w(β/H_p)(T_p/100 GeV)(g_⋆/100)^1 / 6  .

κ_turb represents the efficiency of vacuum energy being converted into turbulent flow:

    κ_turb =ϵ̃κ_v  ,

where the ϵ̃ is set to 0.1.


The total contribution to the GW spectra can be calculated by summing these individual contributions:

    Ω_GWh_100^2=Ω_coh_100^2+ Ω_swh_100^2+ Ω_turbh_100^2  .


We show the GW spectra Ω_GWh_100^2 for four benchmark points in figure <ref>. These four benchmark points are shown in Table <ref>. We choose λ_ϕ h=6.8 and 7.0 and for each value of λ_ϕ h we choose two bubble wall velocities v_w=0.1,0.6. The colored regions represent the sensitivity
curves for future GW detectors LISA <cit.> and TianQin <cit.> with the signal-to-noise ratio (SNR) about 5. We can see that the LISA and
TianQin could detect this new DM mechanism when the bubble wall velocity is relatively large. Taiji <cit.>, BBO <cit.>, DECIGO <cit.>, Ultimate-DECIGO <cit.> could also detect this new DM mechanism by GW signals.








§ DISCUSSION AND CONCLUSION


In this work, we have systematically discussed the gauged Q-ball DM formed 
during a FOPT. We have investigated the stability of gauged Q-balls, including quantum stability, stress stability, fission stability, and classical stability. Different from the global Q-balls, the gauge interaction restricts the size of the stable gauged Q-balls. For a given value of the gauge coupling, the stable
gauged Q-balls can only be realized in the region of charge Q_s<Q<Q_max. The upper limit Q_max and lower limit Q_s mainly comes from the stress stability and  quantum stability criterion respectively. By using the thin-wall approximation, we show that the piecewise model can describe the basic properties of gauged Q-balls in FLSM model well. Based on this, we further give an approximately analytic evaluation of Q_max by using the mapping method. We find the maximal charge is approximately Q_max∝g̃^-2 where g̃ is the gauge coupling of the dark U(1) symmetry. The constraint on the value of gauge coupling g̃ is given by Eq. (<ref>) if the gauged Q-balls are produced by the FOPT. We discuss the relic density of gauged Q-ball DM formed during electroweak FOPT. Even in the minimal SM plus singlet electroweak FOPT model, the gauged Q-balls can comprise all the observed DM. And we have found that in order to satisfy the relic abundance of DM, the original DM asymmetry surprisingly coincides with the observed large lepton number asymmetry. The charge and mass of gauged Q-ball DM can be varied by modifying the phase transition dynamics or the primordial DM asymmetry.

Besides, we give combined constraints on the gauged Q-ball DM from DM direct detection (Mica, XENON1T, Ohya), and astronomical observations (CMB, neutron stars, white dwarfs).
The formation process of gauged Q-ball DM during FOPT also provides phase transition GW signals which could be detected by future GW experiments such as LISA, TianQin, and Taiji. The phase transition dynamics in the early universe provides new formation mechanisms of various soliton DM or dynamical DM. For example, it is reasonable to discuss other species of soliton DM formed during FOPT, such as gauged Fermi-ball DM. The configuration of the Fermi-ball is different from the Q-ball because of the extra Fermi-gas degeneracy pressure. We leave these in our future works.






  
Acknowledgements: 
The authors thank
Emin Nugaev, Julian Heeck, Michael Baker, Kiyoharu Kawana, Mikhail Smolyakov, Mikheil Sokhashvili, Yakov M. Shnir, Dmitry Levkov, Muhammad Fakhri Afif, Chris Verhaaren, Rebecca Riley, Sida Lu, Bingrong Yu, Jiahang Hu for helpful correspondence. This work was supported by the National Natural Science Foundation of China (NNSFC)
under Grant No. 12205387, by Guangdong Major Project of Basic and Applied Basic
Research (Grant No. 2019B030302001), and 
by KIAS Individual Grants under Grants No. PG021403 (PK). 

JHEP












