The fate of small classically stable Q-balls

The smallest classically stable Q-balls are, in fact, generically metastable: in quantum theory they decay into free particles via collective tunneling. We derive general semiclassical method to calculate the rate of this process in the entire kinematical region of Q-ball metastability. Our method uses Euclidean field-theoretical solutions resembling the Coleman's bounce and fluctuations around them. As an application of the method, we numerically compute the decay rate to the leading semiclassical order in a particular one-field model. We shortly discuss cosmological implications of metastable Q-balls.


Introduction
It is well-known that Q-balls [1][2][3], the solitonic "bags" with attractive globally charged particles inside, are rock stable if their total charges Q are large and the model has a mass gap m. Indeed, the net masses E Q of these objects include negative energy of particle attraction and, as a consequence, grow slowly as the new particles are added, E Q ∝ Q a , a < 1 [3]. This makes sufficiently large Q-balls lighter than the sets of massive particles with the same charge, E Q < mQ, and their fission -forbidden by energy and charge conservation. 1 Exceptional stability singles out Q-balls as the favorite toys for cosmologists [6]. These objects may be generated in the early Universe via first-order phase transition [7][8][9] or fragmentation of the scalar condensate [10][11][12] to form dark matter at the present epoch. They may also participate in the Affleck-Dine baryogenesis [13][14][15] catalyzing the process and hiding the baryon number from the sphaleron transitions [16,17]. In the recklessly exotic scenario Q-balls may even impersonate black holes [18] posing new challenges for the astrophysical observations [19,20]. As a pleasant bonus, the objects of this kind invariably appear [21] in supersymmetric models with flat directions which remain in the club of our favorite extensions of the Standard Model despite obvious lack of support from the experimental data [22]. Weak dependence of Q-ball mass E Q on its charge has another, less widely recognized consequence. At smallest Q this mass exceeds the minimal energy mQ of free particles, which makes small Q-balls generically unstable [1][2][3] in some region Q < Q s , E Qs = mQ s , see Fig. 1a. Still, these objects are classically stable if Q is slightly below Q s . Indeed, the energy deficit E Q − mQ of their decay would be small, as well as the momenta of the final-state particles. But finite-size solitons like that cannot evolve classically into waves/particles of much larger wavelength! This gives a range of charges Q c < Q < Q s where the Q-balls exist, remain classically stable and yet, decay into free particles quantum-mechanically. Being local minima of energy at fixed charge, they cannot "leak" the particles one-by-one. Rather, they remain metastable for a long time and then explode in a spectacular collective tunneling event sketched in Fig. 1b. We stress that the metastability window Q c < Q < Q s of global Q-balls appears in all three-dimensional models with bounded energy considered in literature, see e.g. [1,3,23,24].
In this paper we develop general semiclassical method to compute the decay rate of metastable Q-balls at Q, Q c , Q s 1. To the best of our knowledge, no method of this kind has been ever proposed before, with the closest example describing quantum collapse of a metastable Bose condensate in condensed matter physics [25]. There are two reasons, why. First, the states of Q-balls are not just the local minima of energy like the false vacua [26], but the minima at a fixed charge Q. Second, the charged scalar particles inside the Q-ball are typically described by an oscillating complex scalar field ϕ(x, t) ∝ e iωt which grows without bound in Euclidean time τ ≡ it. These features do not enable one to use powerful techniques developed for false vacuum decay [27][28][29], or even perform Wick rotation at the start of calculations. Our semiclassical method is derived by evaluating path integral for the Q-ball decay rate in the saddle-point approximation [30][31][32][33][34] and going into Euclidean time afterwards. This general procedure is justified at Q 1. We arrive at the semiclassical expression for the rate in terms of a certain classical solution ϕ = ϕ cl (x, τ ) which resembles the Coleman's bounce [27] in Euclidean time τ and yet, differs from it in important details. Our solution interpolates between the Q-ball at τ → −∞ and configuration in the catchment area of the vacuum at τ = 0. Unlike the bounce, it can be continued to Euclidean time only in a very specific way, by turning the charged field ϕ and its complex conjugateφ into independent real functions of x and τ . The real ratio ϕ cl /φ cl satisfies certain boundary conditions at τ = −∞ and τ = 0. Given the solution, one computes the probability of Q-ball decay per unit time, where the suppression exponent F Q is related to Euclidean action on ϕ cl and the prefactor A Q includes fluctuation determinant around this solution. In the semiclassical limit Q ∼ Q c 1 the exponent F Q ∝ Q c becomes large, and A Q ∼ mQ 1/2 c . As a by-product, we obtain expression for the width of metastable Q-balls interacting with neutral finite-temperature bath. In this case the semiclassical solution lives on a finite Euclidean time interval and expressions for the exponent and the prefactor are modified accordingly.
We illustrate the semiclassical method in the model of a complex scalar field ϕ(x, t), with scalar potential shown in Fig. 2a. Overbar in Eq. (1.2) and below denotes complex conjugation, expression for V (ϕφ) will be given in the main text. Note that our model approximates renormalizable Friedberg-Lee-Sirlin two-field model [1] in the limit when the additional field is very massive and can be integrated out. Besides, the potential in Fig. 2a is flat at large fields: V ≈ Λ 4 at |ϕ| > v. In this regard it is similar to potentials appearing in supersymmetric theories with flat directions [21]. The model (1.2) has a family of Q-ball solutions with mass graph E Q = E Q (Q) similar to the one in Fig. 1a. These objects are metastable at charges between Q c ≈ 266 v 4 /Λ 4 and Q s ≈ 1.6 Q c . At v > 5 Λ the values of Q c and Q s are large and semiclassical approximation can be used for all metastable Q-balls. We numerically found the respective semiclassical solutions ϕ cl (x, τ ) and computed the suppression exponent F Q of the decay rate, see Fig. 2b. As expected, the value of F Q tends to zero and infinity in the limits Q → Q c and Q → Q s corresponding to classically decaying and absolutely stable Q-balls. Leaving numerical calculation of the fluctuation determinant to further studies, we roughly estimate the prefactor in Eq. (1.1) as In cosmology metastable Q-balls may form naturally long-living decaying dark matter. In the most exciting scenario the lifetime of these objects is comparable to the age of the Universe, so that their decays can affect structure formation and alleviate [35] emerging tension between the low-redshift [36,37] and high-redshift (CMB) [38] cosmological measurements. This requires, however, moderately small Q-ball charges Q ∼ Q c ∼ F Q : in the above model Γ −1 Q is of order 10 10 years if Q ∼ Q c ∼ 10 2 , see Fig. 2b. Our estimates show that the Q-balls with these charges may be generated in the early Universe via phase transition or fragmentation of the scalar condensate if the generation temperature is raised to T ∼ 10 11 GeV, cf. [39]. Besides, alternative generation mechanisms -say, pair production [40] -may be essential at small Q. This suggests a path to new cosmological scenarios with metastable Q-balls which may be considered separately.
In the introductory Sec. 2 we explain our semiclassical method and show main numerical results. Derivation of the method and details of its numerical implementation are given in Secs. 3 and 4, respectively. We discuss possible cosmological applications of the metastable Q-balls in Sec. 5.

Semiclassical method and its application 2.1 Metastable Q-balls
We start by briefly reviewing the properties of small Q-balls. To be concrete, consider the model (1.2) of complex scalar field ϕ with the potential shown in Fig. 2a. The scalar bosons in this model have mass ≈ m in vacuum ϕ = 0 and become almost massless at large ϕ: V → Λ 4 ≡ m 2 v 2 log(e b + 1)/b as |ϕ| → +∞. This corresponds to short-range attraction between the bosons. Importantly, the potential (2.1) approximately describes decoupling limit of the renormalizable Friedberg-Lee-Sirlin model [1], as we will argue in Sec. 5. The model (1.2) possesses global conserved charge related to phase rotation symmetry ϕ → e −iα ϕ,φ → e iαφ . Conservation of this quantity is vital [2] for the existence and stability of Q-balls -nontopological solitons carrying nonzero Q.
In what follows we use semiclassical approximation which is justified at g ≡ m/v 1. Indeed, change of variables ϕ → vϕ, x µ → x µ /m eliminates all parameters from the classical action (1.2), (2.1) except for the combination g −2 right in front of it. This means that g 2 1 is the semiclassical parameter in our model. Besides, the charge (2.2) also becomes proportional to g −2 after the change of variables suggesting that Q 1 in the semiclassical regime. The simplest way to find nontopological solitons in the model (2.1) is to substitute stationary spherically-symmetric Ansatz into the classical field equations and numerically solve the remaining ordinary differential equation for the real function χ Q (r) with regularity conditions ∂ r χ Q (0) = χ Q (∞) = 0, see Appendix A for details. This gives a family of localized solutions parametrized with the frequency ω, see two of them in Fig. 3a. In Fig. 3b we show the charges Q = Q(ω) of all solutions, Eq. (2.2). Quite surprisingly, no solution exists at Q < Q c ≈ 266/g 2 . At the same time, two solutions, the "Q-ball" and "Q-cloud", are present at each Q > Q c ; the respective subfamilies are marked by solid and thick-dashed lines in Fig. 3b. To explain these features, we recall the main property of nontopological solitons proven in [2]: stationary solutions of the form (2.3) are, in fact, extrema of energy in the subsector of field configurations {ϕ(x), ∂ t ϕ(x)} Q with fixed charge Q. Computing the total masses (2.4) of all solutions, we find that the Q-clouds are always heavier than the Q-balls with the same charge, see Fig. 4a. Besides, the Q-clouds are classically unstable: they do not satisfy the rigorous Vakhitov-Kolokolov criterion dQ/dω ≤ 0 [41][42][43] which is necessary for stability, cf. Fig. 3b. This suggests the structure of the energy functional in the subsector of configurations with fixed Q shown in Fig. 4b. The Q-ball and a configuration describing Q free particles at rest are the two local minima of E at fixed Q separated by an energy barrier. The Q-cloud, to the contrary, is an unstable stationary solution "sitting" precisely on the barrier top [44]; the analogous solutions in scalar theories with false vacuum decay and in gauge theories are called critical bubble [26,27,29] and sphaleron [45], respectively. In Appendix A we confirm the picture in Fig. 4b by demonstrating that all Qballs in our model are classically stable, while the Q-clouds have precisely one unstable mode.
In quantum theory the global minimum of energy in Fig. 4b corresponds to true ground state at a given Q, while the other, local minimum represents metastable state decaying via tunneling through the potential barrier. It is clear from Fig. 4a that the Q-balls are metastable at Q c < Q < Q s when their masses exceed the energies mQ of free particles with the same charge. Below we compute the lifetimes Γ −1 Q of these objects.
Direct inspectation of literature [1,3,23,24,44,46,47] shows that the above picture is typical in weakly interacting theories with bounded energies possessing Q-ball solutions. Namely, nontopological solitons in these models can be sorted into classically stable Q-balls -local minima of energy at fixed charge -and Q-clouds "sitting" on tops of the potential barriers between the Q-balls and free particles in the vacuum. The branches of Q-balls and Q-clouds meet at the critical point Q = Q c , E = E Q (Q c ) corresponding to the smallest soliton at the verge of classical stability. By construction, the Q-cloud masses are higher than the minimal values of energy E Q and mQ. Thus, the critical mass also satisfies E Q (Q c ) > mQ implying that there exists a window of charges Q c < Q < Q s where the Q-balls are metastable, E Q > mQ.

Euclidean solutions for the decay probability
Let us preview our semiclassical recipe to calculate Γ Q postponing its derivation till Sec. 3. For illustrative purposes we will compare this recipe to the powerful Euclidean technique [27][28][29] developed for the decay of metastable (false) vacuum in models of real scalar field(s) φ without conserved charges. The latter process is described by the "bounce" -real Euclidean solution φ cl (x 2 +τ 2 ) interpolating between the false vacuum at τ → −∞ and the catchment area of the true vacuum at τ = 0. After continuation to Minkowski time t ≡ −iτ > 0 this solution describes classical evolution of an expanding true vacuum bubble in the final state. The rate of false vacuum decay is given by the expression similar to Eq. (1.1), where the leading suppression exponent is Euclidean action of the bounce S E [φ cl ] and the prefactor includes fluctuation determinant around this solution.
We expect to find similar, though properly modified, procedure for computing the rate of Q-ball decay. In Sec. 3 we indeed introduce the semiclassical solution ϕ cl (x, τ ), ϕ cl (x, τ ) satisfying Euclidean field equations in the model (1.2), where V is a derivative of V (φϕ) with respect to its argument. Our solution is spherically-symmetric, i.e. depends on r ≡ |x| and τ . It is natural to expect that it coincides with the Q-ball (2.3) at the start of the process, where we introduced the time shift η 0 that cannot be excluded in general. In Sec. 3 we obtain the boundary condition which is consistent with the asymptotics (2.6). Soon we will verify numerically that the solutions of Eqs. (2.5), (2.7) automatically satisfy Eq. (2.6). Surprisingly, ϕ cl andφ cl in Eq. (2.6) are not complex conjugate to each other. This is not a problem, however, for the complex semiclassical method adopted in Sec. 3. Instead, our derivation suggests that ϕ cl andφ cl should be considered as independent real functions of x and τ . Then the value of Euclidean action is real, where β will be sent to infinity in the end of the calculation. Like the bounce, our semiclassical solution arrives to the catchment area of the true vacuum at τ = 0. To quantify this criterion, we use the second boundary condition derived in Sec. 3, This means that the solution is symmetric with respect to time reflections, ϕ cl (x, τ ) =φ cl (x, −τ ). As a consequence, after continuation to Minkowski time t = −iτ > 0 the functions ϕ cl andφ cl become complex conjugate to each other, i.e. represent ordinary classical evolution after the decay. We expect to obtain free waves/particles classically evolving in the vacuum as t → +∞.
Solving Eqs. (2.5) with boundary conditions (2.7), (2.9) at large β, one finds a family of real solutions ϕ cl ,φ cl parametrized with ωβ + η 0 in Eq. (2.7). Apparently, different values of this combination correspond to solutions approaching different Qballs at large negative τ . In what follows we express η 0 = η 0 (Q) using Eq. (2.2) and characterize the solutions with charge Q. Sending β → +∞, we obtain ϕ cl andφ cl interpolating between the Q-balls (2.6) at τ → −∞ and the sector of true vacuum at τ = 0. After that the solutions are continued to τ > 0 using the time reflection symmetry.
Given ϕ cl ,φ cl , η 0 , one computes the suppression exponent of the decay rate (1.1). Note that the difference of Euclidean actions (2.8) on ϕ cl ,φ cl and Q-ball (2.3) in this expression remains finite as β → +∞ because our solutions approach ϕ Q ,φ Q at large |τ |. The third term in Eq. (2.10) is specific to models with conserved charge Q. Notably, it involves the parameter η 0 from Eq. (2.7) at the place typically occupied by the chemical potential of Q in statistical physics. In Sec. 3.3 we derive the Legendre transform formula dF Q /dQ = η 0 supporting this intuition.
Our expression for the prefactor in Eq. (1.1) is way more complicated. It uses small perturbations in the background of the semiclassical solution, where e ±ωτ compensate for exponential behavior of ϕ cl andφ cl as τ → ±∞. The effect of perturbations on the Euclidean action is characterized by its second variation, i.e. the operator Here we introduced the rescaled fields χ cl ≡ e −ωτ ϕ cl andχ cl ≡ e ωτφ cl remaining finite as τ → ±∞ and denoted V cl ≡ V (ϕ clφcl ); the primes are derivatives of this function with respect to its argument. Note thatD cl is Hermitean: it acts on the perturbations (δϕ, δφ) T vanishing at τ → ±∞. The prefactor in Eq. (1.1) involves determinant of this operator, where det means that all zero eigenvalues are excluded from the determinant,D Q is the same asD cl but in the Q-ball background i.e. with χ cl ,χ cl → χ Q (r) in Eq. (2.11), the factors cosh 2 η 0 and B 0 = 2 dτ d 3 x (∂ τ χ cl ) 2 were introduced while deleting zero modes from the determinants. Like in the case of bounce, we expect that det D cl has opposite sign to the same determinant det D Q computed in the background of the metastable solution. Expression (2.12) for the prefactor seems too involved for practical calculations (see, however, [28,34,48,49]). But it is straightforward to estimate A Q by extracting its dependence on the only small parameter in the model -the semiclassical parameter where the mass m is restored on dimensional grounds. start from the stationary Q-ball χ Q (r) (points in Fig. 6a) at τ → −∞ and become lower and wider at τ = 0 where ϕ clφcl < v 2 . To confirm that they indeed arrive to the vacuum sector, we continue ϕ cl andφ cl to real time by numerically solving the field equations from t = τ = 0 to large t, see Fig. 6b. As expected, Minkowski evolution of our solutions describes spreading wave packets corresponding to free particles in the vacuum. Given the solutions, we numerically compute the suppression exponent (2.10) (solid lines in Figs. 2b and 7a). We see that F Q equals zero and infinity at charges Q c and Q s corresponding to classically unstable and absolutely stable Q-balls, respectively. To find practical fitting formula for F Q , let us study its behavior at Q → Q s when the total energy deficit Q p ≡ E Q − mQ of the decay is small. This corresponds to small momenta p ∼ (2m p ) 1/2 ∝ (1 − Q/Q s ) 1/2 of the final-state particles/waves and, as a consequence, slow evolution of the solution near τ = 0:φ cl ϕ cl ∝ e τ p . Thus, large time delay η 0 ∝ −1 p ∝ (1 − Q/Q s ) −1 is needed for ϕ clφcl to approach the Q-ball profile. Integrating the Legendre formula dF Q /dQ = η 0 , we obtain the asymptotic of the suppression exponent

Applying the method
Slightly generalizing Eq. (2.13), we obtain the fitting formula in the entire metastability window Q c < Q < Q s , (2.14) We find that Eq. (2.14) with c 1 = −0.28 and c 2 = −2.6 describes our numerical results with 5% relative precision, see the dashed line in Fig. 7a.

Decay at finite temperature
Consider metastable Q-ball immersed into neutral plasma at temperature T . Thermal fluctuations should decrease the lifetime Γ −1 Q, β of this object kicking it over the potential barrier in Fig. 4b. Fortunately, our semiclassical method is easily generalized to the case of finite temperature T ≡ β −1 : one just considers the above semiclassical solutions in a finite time interval −β/2 < τ < β/2 and imposes Eq. (2.7) at τ = −β/2. This makes the functions ϕ cl andφ cl quasiperiodic: their values at τ = ±β/2 differ by multiplicative factors e ±(ωβ+η 0 ) . The suppression exponent is still given by Eq. (2.10), where S E is the Euclidean action (2.8) in the interval |τ | < β/2. Generalization of the prefactor formula (2.12) is slightly less trivial, see comments in Sec. 3.5.
We numerically obtained the quasiperiodic solutions in the model (2.1). Their suppression exponent F Q, β is represented by the solid ("tunneling") line in Fig. 7b. Notably, the solutions of this kind do not exist at small β and therefore fail to describe thermal transitions in that region. To consider high-temperature processes, we include a trivial semiclassical solution -the Q-cloud -which has the same form (2.3) as the Q-ball of the same charge, but with different frequencyω and profileχ Q (r). This solution satisfies Eqs. (2.5), (2.7), (2.9) at arbitrary β. Instead of interpolating between the sectors of Q-ball and the true vacuum, it simply sits on top of the barrier between them, see Fig. 4b. Computing the exponent (2.10) for the Q-cloud, one obtains linear function F Q, β = β(Ẽ Q − E Q ) shown by the diagonal dashed line in Fig. 7b, whereẼ Q is the mass of this object.
The picture in Fig. 7b is typical [50] for thermal transitions, with two semiclassical contributions to the rate representing direct tunneling through the potential barrier in Fig. 4b and jumps onto its top (activation) induced by thermal fluctuations. The overall suppression exponent corresponds to the least suppressed contribution at each β, given by the "tunneling" and "activation" lines in Fig. 7b at temperatures lower and higher β −1 c , respectively. As a further generalization of the method, one may consider decay of Q-ball in charged plasma at temperature T and chemical potential µ. This process, however, is different from the physical viewpoint [51]. The Q-ball exchanges charge with the medium and reaches equilibrium at Q = Q µ . The respective initial state is described by the grand canonical ensemble giving for the decay rate, where the prefactors are ignored. This integral is saturated either by one of the integration limits or in the vicinity of the saddle point µ = η 0 ; one has to select the least suppressed contribution at every µ and β. We leave this interesting calculation to future studies.

Derivation of the semiclassical method
We are going to derive the semiclassical expression for the Q-ball decay rate in two steps. First, we write this rate in the form of path integral. Second, we evaluate the integral in the saddle-point approximation. We will see that this technique is trustworthy at Q 1.

Quantum states for Q-balls
The first nontrivial step of our program is to define the states of quantum Q-balls. Note that identification of these objects with the classical solutions (2.3) is not satisfactory in quantum theory where the Q-ball state |Q should be explicitly specified in order to compute its decay amplitude. Our definition of Q-balls essentially relies on the presence of a conserved charge in the theory. In the simplest model (1.2) the charge Q is associated with the phase rotation symmetry ϕ → e −iα ϕ. Then in quantum case the operatorQ generates symmetry transformations, i.e. acts as e iαQ |ϕ,φ = |e −iα ϕ, e iαφ (3.1) on the eigenstates |ϕ,φ of field operatorsφ,φ with eigenvalues ϕ(x) andφ(x). Since e 2πiQ = 1 corresponds to full phase rotation, the charge has integer eigenvalues. Then the operatorP projects onto the subspace of states with given Q; note that integration in Eq. (3.2) is performed along imaginary η.
Using the projector (3.2) one can study the subsector of states with fixed Q. It is natural to define quantum Q-ball as a state |Q of minimal energy within this subsector [2]. We obtain the limiting formula whereP Q projects an arbitrary state |i onto the charge-Q subsector, while the operator e −βĤ/2 at large β suppresses all excited states within this subsector leaving only the minimal-energy eigenstate |Q of the HamiltonianĤ. For simplicity we use normalization 2 Q|Q = 1.
In fact, we are interested in small Q-balls representing local (false) minima of energy at fixed Q, see Fig. 4b. One can expect that Eq. (3.3) is still applicable in this case if the initial state |i belongs to the catchment area of the Q-ball and Euclidean time interval β is not exponentially large.

Decay probability in the form of path integral
The total probability of Q-ball decay in time t 0 is obtained by time-evolving the state |Q and projecting it onto the basis of Fock states |f above the vacuum, where Eq. (3.3) was used to obtain the second equality and the limits β → +∞, t 0 → +∞ are assumed from now on. Importantly, all initial and final states |i and |f in Eq. (3.4) belong to the catchment areas of the Q-ball and true vacuum, respectively. In particular, due to unitarity of quantum evolution one would obtain unit probability, if summation over all Hilbert states |Ψ was performed instead of |f .
2 Assuming that a finite spatial box with some boundary conditions is introduced.
To write path integral for the probability (3.4), we use the basis of configuration eigenstates |i ≡ |ϕ i ,φ i , |f ≡ |ϕ f ,φ f and Eqs. (3.1), (3.2) for the projector 3 One can interpret the propagation operators in this expression as describing time evolution from t = iβ/2 to t = t 0 to t = −iβ/2 along the complex time contour C in Fig. 8a. Using path integral representation for the evolution operators, one finally arrives at where the configurations ϕ(x, t),φ(x, t) live on the contour C and satisfy quasi-periodic conditions ϕ(x, t + iβ) = e −η ϕ(x, t) ,φ(x, t + iβ) = e ηφ (x, t) (3.7) at t = −iβ/2. The classical action S[ϕ,φ] is computed along the same contour. Note that the functions ϕ andφ can be continued to the entire Euclidean axis using Eq. (3.7).
Recall that the Euclidean parts of the contour C in Fig. 8a are not related to Wick rotation, they appeared due to the minimization procedure (3.3). Thus, representation (3.6) does not rely on analytic properties of the Q-ball decay amplitude, cf. [29]. Note also that that the "final-state" configurations ϕ(x, t 0 ) ≡ ϕ f (x) in Eq. (3.6) should describe free particles in the vacuum; this makes our integral for the probability essentially different from that for the unity (3.5). Likewise, ϕ(x, ±iβ/2) should belong to the catchment area of the Q-ball.

Saddle-point approximation
In Sec. 2.1 we argued that the action S and charge Q take large values in the semiclassical regime g 2 1. In this case the integral (3.6) of the fast-oscillating exponent e iS can be evaluated in the saddle-point approximation.
We introduce the saddle-point configuration {φ cl , ϕ cl , η cl } as an extremum of the integrand in Eq. (3.6). It will be convenient to consider ϕ cl (x, τ ) andφ cl (x, τ ) as functions of Euclidean time τ ≡ it which takes real and imaginary values on the contour in Fig. 8a. By construction, ϕ cl andφ cl satisfy the field equations (2.5) and the boundary conditions (3.7), Besides, extremization with respect to η gives equation where Eqs. (2.5), (3.8) were used in the differentiation. Note that Eq. (3.9) coincides with the classical expression (2.2) for the charge. Importantly, ϕ cl andφ cl are not necessarily complex conjugate to each other: integrations in Eq. (3.6) can be continued to independent ϕ-andφ-contours in the functional space. Solving the field equations (2.5) with the quasiperiodicity conditions (3.8), one obtains the saddle-point configuration {ϕ cl ,φ cl } for every η cl . The latter parameter is then related to the conserved charge by Eq. (3.9). Now, we compute the integral (3.6) by noting that at g 2 1 (S 1) its integrand is sharply peaked in small vicinity of the saddle point, ϕ = ϕ cl + e ωτ δϕ ,φ =φ cl + e −ωτ δφ , η = η cl + δη . Here the factors e ±ωτ in front of field perturbations compensate for exponential growth of ϕ cl andφ cl as τ → ±∞, cf. Eq. (2.6). Substituting Eq. (3.10) into Eq. (3.6) and taking the Gaussian integrals over small δϕ(x, τ ), δφ(x, τ ) and δη, one obtains the standard saddle-point formula, Here N is the unknown normalization factor from the functional measure, the second multiplier is the saddle-point value of the integrand, while the determinant ofD cl defined in Eq. (2.11) and the factor dQ/dη cl = id 2 S/dη 2 account for fluctuations of {δϕ, δφ} and δη, respectively. Now, let us us simplify the above semiclassical recipe. First, the double-bent Minkowskian part of the contour C is redundant: analytic functions ϕ cl andφ cl can be considered on the Euclidean time axis in Fig. 8b. There remains, however, an important selection rule: when continued to real time, these functions should describe the final state with free particles at t = t 0 → +∞.
Second, the saddle-point equations have two discrete symmetries at τ ∈ R: simultaneous complex conjugation of all fields ϕ cl ,φ cl , η cl and time reflection In general case this gives four equally suppressed saddle-point solutions, each producing a term with complex exponent in Eq. (3.11). We expect, however, to find a single dominant saddle-point configuration, just like for other tunneling processes. Then this saddle point is symmetric, where Eqs. (2.4), (2.8), (3.9) were used. We arrive at the convenient expression (2.10) for F Q . Fourth, it is natural to expect that the semiclassical solution is spherically-symmetric like the Q-ball, i.e. ϕ cl andφ cl depend on r ≡ |x| and τ .
Finally, we derive the Legendre transformation formula that was used in Sec. 2. To this end we explicitly differentiate the suppression exponent (2.10) with respect to Q at fixed β, use Eqs. (3.9), (3.13), (3.14) and the property of Q-ball dE Q /dQ = ω [1]. We obtain dF Q /dQ = η 0 which gives physical interpretation of η 0 . This relation can be also used as a cross-check of numerical calculations.

Expression for the prefactor
Let us further simplify Eq. (3.11). To cancel the unknown constant N , we divide this expression by unity (3.5) which is given by the same path integral as in Eq. (3.6) but with an important distinction: now the integration runs over all final-state configurations ϕ(x, t 0 ),φ(x, t 0 ), not just the ones from the vacuum sector. As a consequence, the dominant saddle-point configuration for that integral is the Q-ball {ϕ Q ,φ Q , η 0 = 0}, Eq. (2.3). We obtain the saddle-point result similar to Eq. (3.11), whereD Q is the same fluctuation operator as in Eq. (3.11) but in the Q-ball background. Using Eq. (3.13) one finds that dQ/dη → β −1 dQ/dω + O(β −2 ). This factor cancels in the ratio of Eqs. (3.11) and (3.15), where F Q is the suppression exponent from Eq. (3.11). The operatorD cl introduced in Eq. (2.11) is Hermitean; recall that it acts on perturbations ψ ≡ (δϕ, δφ) T with quasiperiodic boundary conditions δϕ(x, τ + β) = e η 0 δϕ(x, τ ) , δφ(x, τ + β) = e −η 0 δφ(x, τ ) , (3.17) cf. Eqs. (3.7) and (3.10). As a consequence, the eigenmodes ψ k (x, τ ) of this operator form orthonormal basis 4 and its eigenvalues λ k are real, We immediately run into a problem: detD cl = 0 due to zero eigenmodes 4 We still assume that the finite-size spatial box is introduced.
It is clear that the functional integration in the directions of zero modes cannot be performed using the saddle-point method. Indeed, in the standard case one represents and takes the Gaussian integrals over c k . On the other hand, parameter c 0 in front of ψ . The respective integral is not Gaussian; rather, it gives the full period of the process, Thus, we should exclude zero eigenvalue of ψ The modes due to phase rotations and spatial translations are treated in the same way as ψ where the asymptotics (2.6) were used in Eq. (3.18). Collecting all factors, we obtain expression (2.12) for the prefactor in Eq. (1.1), where det is the product of all nonzero eigenvalues of the operator and we computed the norm B 0 of ψ (0) 0 using Eq. (3.18). By definition, the determinants in Eq. (2.12) involve quasiperiodic boundary conditions (3.17) with β → +∞. However, once all zero modes are excluded, any kind of boundary conditions with finite δϕ can be imposed at τ → ±∞. Indeed, the eigenvalue problem forD cl has the form of a stationary Schrödinger equation for the four-dimensional particle with a spin. The regions τ → ±∞ are classically forbidden for this particle: the background solution at large |τ | coincides with the Q-ball which is classically stable, i.e. has only exponentially growing or decaying with τ perturbations ψ k . It is clear that the values of λ k do not depend on the boundary conditions imposed deep inside the classically forbidden regions. Thus, one can consider all determinants in Eq. (2.12) with vanishing boundary conditions δϕ, δφ → 0 as τ → ±∞ instead of the quasiperiodic ones.

Finite-temperature solutions
Now, consider Q-ball decay at a finite temperature T ≡ β −1 . The initial state of this process is specified by the density matrixρ β = Z β e −βĤP Q , whereP Q restricts the statistical ensemble to states with charge Q, while Z β is a normalization constant ensuring trρ β = 1. As before, we tacitly assume thatρ β is projected onto the states |i from the Q-ball sector. Note, however, that the quantity (3.4) used in the above calculation involves the same density matrixρ β ≡ e β(E Q −Ĥ)P Q but with different normalization. Moreover, Eq. (3.5) at finite β gives, where the trace is taken in the Q-ball sector. Thus, the ratio of Eqs. (3.4) and (3.5) at finite β that was computed above gives the probability of Q-ball decay at temperature T = β −1 .
To generalize Eq. (2.12) for the prefactor, one should avoid simplifications related to the limit β → +∞. At finite temperature the quantities dQ/dη cl in Eqs. (3.11) and (3.15) do not cancel, the original quasiperiodic conditions (3.17) should be used for all operators, and orthogonal family of zero modes with norms B α should be properly defined. We leave derivation of finite-temperature prefactor to interested readers.

Numerical implementation
In this Section we solve the boundary value problem for the semiclassical solution ϕ cl (r, τ ),φ cl (r, τ ). To this end we adopt units with m = 1 and set v = 1 by field rescaling.

Framework
It is convenient to consider functions with finite asymptotics, y cl (r, τ ) = r e −ωτ ϕ cl ,ȳ cl (r, τ ) = r e ωτφ cl , where ω is the frequency of the decaying Q-ball. We introduce 5 N r × N τ = 150 × 300 lattice with sites {r j , τ i } = {∆r(j + 1), −∆τ (i + 1/2)}, 0 ≤ j ≤ N r , 0 ≤ i ≤ N τ , and uniform spacings ∆r, ∆τ . The spatial and temporal extents of our lattice r j < L r and −β/2 < τ i < 0 are large enough for the solutions to approach their asymptotics: L r = 20, β = 40 ÷ 200. Our unknowns are the values of y j, i ≡ y cl (r j , τ i ) andȳ j, i on the lattice sites. We discretize the field equations (2.5) in the standard second-order manner using and similar expressions for the derivatives with respect to r. Boundary conditions (2.7), (2.9) have natural lattice generalization, where τ 0 = −∆τ /2 and τ Nτ −1 = −β/2 + ∆τ /2 are the first and the last lattice sites. Finally, we impose Dirichlet and Neumann conditions at r = 0 and L r , respectively, cf. Eq. (4.1). After discretization Eqs. (2.5) with boundary conditions (4.2), (4.3) constitute a sparse system of 2N r N τ nonlinear algebraic equations for the same number of unknowns ϕ j, i ,φ j, i . A convenient technique to solve large systems of this kind is provided by the Newton-Raphson method, see e.g. [52]. One starts with the initial guess y for the solution. Substituting y = y (0) + δy,ȳ =ȳ (0) + δȳ into the lattice equations and ignoring nonlinear terms in δy j, i , δȳ j, i , one arrives at the sparse linear system that gives δy and δȳ. Then one refines the guess by setting y (0) → y (0) + δy,ȳ (0) → y (0) + δȳ and repeats the procedure. The iterations converge to the exact solution if the original guess was good enough. Finally, one adjusts the value of η 0 in Eq. (4.2) to fix the charge of the solution to that of the Q-ball with frequency ω. Note that the above algorithm should be supplemented with fast linear solver, cf. [33].

Obtaining the solutions
A drawback of the Newton-Raphson method is its high sensitivity to the initial guess: if y (0) ,ȳ (0) are not close enough to the semiclassical solution, the method may diverge or produce useless trivial solutions -the Q-ball, Q-cloud, or charged condensate in the vacuum. This poses a problem of finding an approximate solution at some β and η 0 .
To construct such a solution, we use a nontrivial observation [30] that the (quasi)periodic finite-β solutions describe, besides the thermal decays, transitions between the sectors of Q-ball and vacuum at fixed energies E. Indeed, thermal density matrix is diagonal in the basis of energy eigenstates, which implies that the semiclassical solutions describing fixed-energy processes contribute with factors e −βE to the thermal rates. This gives physical interpretation to all our quasiperiodic solutions, even the ones producing subleading contributions in Fig. 7b.
A particular semiclassical solution with energy E close to the Q-cloud massẼ Q describes transmission through the infinitesimally small classically forbidden region near the barrier top in Fig. 4b. It should coincide with the Q-cloud, a configuration "sitting" on top of the barrier, up to small corrections. We demonstrated in Appendix A, however, that almost all small perturbations of the Q-cloud grow exponentially with τ and therefore cannot satisfy the quasiperiodic conditions (2.7). The notable exception is the negative (unstable) mode which periodically depends on τ via e iγτ . In Fig. 9a we plot the profile ξ − (r) of this mode obtained in Appendix A. Adding it to the Q-cloud y Q (r) we obtain the approximate solution with charge Q and period β = 2π/γ, where we used Eq. (A.2) from Appendix A and recalled that ξ * − (r) e −iγτ also satisfies the linearized field equations in the background of Q-cloud. Besides, we introduced the Q-cloud frequencyω and small amplitude A − of the perturbation. It is straightforward to check that the configuration (4.4) is real and satisfies the semiclassical equations (2.5), (2.7), (2.9) with η 0 = β(ω − ω).
This suggests the following strategy of numerical calculations at a given chargesay, at Q ≈ 1.33 Q c . We compute the configuration (4.4) with small A − = 10 −3 and use it as an initial guess for the solution with period β = 2π/γ − δβ, δβ = 0.2 and η 0 = 5.7. After several Newton-Raphson iterations the new solution is obtained with acceptable precision. Then we use this solution, in turn, as an initial guess at even smaller period β. Continuing to change the value of β in small steps, we obtain a complete 6 branch of solutions shown by points in Fig. 9b, see examples in Fig. 10. At last, we arrive at the solution with large β (Fig. 5b)  Once the solution with large β is found, we start changing ω and Q(ω) in small steps. As before, at each step we use the solution with charge Q as an initial guess for the one with charge Q ± δQ. This gives the semiclassical solutions in the entire metastability region Q c < Q < Q s . We compute their suppression exponents (2.10) and obtain Figs. 2b and 7a.
We supported the above numerical procedure with several tests. First, we checked that all numerical solutions with large β are close to the Q-balls at τ ≈ −β/2, see Fig. 6a, while their energies E = E(Q) coincide with E Q . Second, the suppression exponent F Q computed via Eq. (2.10) has the expected properties: it equals zero at Q = Q c and grows to infinity as Q → Q s . Third, the numerical values of F Q and η 0 are consistent with the Legendre relation dF Q /dQ = η 0 . Fourth, we checked that the semiclassical solutions and the values of F Q are not sensitive to the lattice parameters L r , ∆r, ∆τ within relative precision of 2%. Besides, we monitored the energy E and charge Q of the solutions which were conserved with accuracy better than 10 −3 . These tests convinced us that the lattice solutions give correct suppression exponent at the precision level of 2%.

Discussion: Prospects for cosmology
We developed general semiclassical method to calculate the decay rate Γ Q = A Q · e −F Q of metastable Q-balls at Q 1. The method can be applied in arbitrary models at the cost of numerically obtaining certain Euclidean solutions ϕ cl (x, τ ) that enter the semiclassical expressions for exponent F Q and prefactor A Q of the rate. We generalized the method to finite-temperature processes.
To illustrate the method, we performed explicit numerical calculations in the model (2.1). In particular, we computed the exponent F Q of the Q-ball decay rate, see Eq. (2.14), and estimated 7 the prefactor A Q ∼ mQ 1/2 c . Notably, the model we use is close to a certain limit of the celebrated Friedberg-Lee-Sirlin (FLS) model [1] which describes complex and real fields ϕ and ζ with renormalizable potential Here m ≡ Λ 2 /v and m ζ ∝ λ 1/2 Λ are masses of the fields, λ is the coupling constant. At m ζ m the field ζ is very massive and cannot be excited in processes with relatively low momenta p ∼ m. Minimizing Eq. (5.1) with respect to ζ at fixed ϕ, one obtains low-energy potential for the remaining field, Note that numerical evaluation of the fluctuation determinant in the prefactor involves application of two-dimensional Gelfand-Yaglom theorem and renormalization, which are interesting open tasks beyond the scope of this paper, see [28,34,48,49].
This potential is close to the smooth one, Eq. (2.1) that was used in our calculations, see Fig. 11a. Thus, our numerical data estimate decay rate of the FLS Q-balls. Although we did not directly address supersymmetric models with flat directions [21], the decay rate of their Q-balls is still expected to exhibit the behavior similar to that in Fig. 2b. In particular, Γ Q should be unsuppressed at the critical charge Q c corresponding to classically unstable Q-ball, it should sharply decrease with charge at Q > Q c and reach zero at the point of absolute stability Q = Q s . Deep inside the metastability window Q c < Q < Q s the suppression exponent of the rate should be large, F Q ∼ Q c 1. Let us discuss possible cosmological application of metastable Q-balls with Q ∼ Q c . These objects may form dark matter with naturally large lifetime, so that their decays in the late-time Universe produce warm dark matter component -energetic ϕ-particles leaving the galaxies. The latter process may change structure formation or explain [35] appearing tension between the low-redshift [36,37] and CMB [38] measurements of the Hubble constant.
Note that the decays of dark matter Q-balls are observable only if their lifetime Γ −1 Q is comparable to the age of the Universe or smaller. This constrains the Q-ball charge Q ∼ Q c 150, where we used estimates F Q ∼ Q c , A Q ∼ mQ 1/2 c in the mid of the metastability window and assumed wide interval of masses m = 1 eV ÷ 10 19 GeV in the prefactor. Thus, we need a mechanism to generate small Q-balls in the early Universe.
The two standard mechanisms of forming dark matter Q-balls involve phase transition [7][8][9] or fragmentation of the scalar condensate [10][11][12] at high temperature T . Typically, these mechanisms act in a similar way, by collecting charge Q from the cosmological horizon of size H −1 into the compact objects which further relax into Q-balls. Assuming initial charge asymmetry 8 ∆ Q = n Q /s, where n Q and s are the charge and entropy densities at temperature T , one estimates the charges of formed Q-balls, where in the second equality we introduced the number of relativistic degrees of freedom g * ∼ 10 2 , used Friedmann equation H ≈ g 1/2 * T 2 /M pl and expression for the entropy density s ∼ g * T 3 . After that we recalled that Q ∼ 10 2 . Equation (5.3) relates the generation temperature T to the asymmetry ∆ Q .
The present-day mass density of Q-balls should coincide with that of dark matter ρ DM ∼ 10 −6 GeV/cm 3 . Immediately after generation the concentration of these objects is n Q−balls = n Q /Q = s∆ Q /Q. Since expansion of the Universe conserves the ratio n Q−balls /s, the mass density of Q-balls at the present epoch is where we introduced the entropy density now s 0 ∼ 10 3 cm −3 and estimated the mass of a small Q-balls as E Q ∼ mQ. Equations (5.3) and (5.4) relate parameters m, T , and ∆ Q of the cosmological model. The standard generation mechanisms naturally work at comparable values of m and T . Taking T ∼ 10 −2 m to prevent immediate decay of Q-balls due to the temperature fluctuations (see Fig. 7b), we obtain ∆ Q ∼ 10 −22 and m ∼ 10 2 T ∼ 10 13 GeV. Note that this generation temperature is much higher than the ones in typical models with large Q-balls, cf. [7][8][9][10][11][12]. Similar small-mass Q-balls were considered as natural candidates for self-interacting dark matter [39].
Since the value of ∆ Q is way smaller than the observed baryon asymmetry ∆ B ∼ 10 −10 , we cannot interpret the charge Q as the baryon number thus relating small Q-balls to the Affleck-Dine baryogenesis. Indeed, bold substitution of ∆ Q = ∆ B into Eqs. (5.3), (5.4) gives essentially different values of mass m ∼ GeV and temperature T ∼ 10 15 GeV which is problematic for the generation mechanisms. Hence, in order to pack the baryon number inside metastable Q-balls one has to consider special mechanisms of their formation or at least some modification of the standard ones. Alternatively, one can avoid identification of the charge Q with the baryon number, since cosmology of small metastable Q-balls is interesting enough by itself.
We want to obtain the parameters γ of regular perturbations ξ R , ξ I . To this end we appeal to the shooting method, again. We discretize Eqs. (A.3) using the lattice {r j } from Sec. 4 and impose Neumann condition ξ R, I (r Nr−1 ) = ξ R, I (r Nr ) at the rightmost link of the lattice. After that we introduce two basis solutions Ψ (A) (r j ) ≡ (ξ  Fig. 11b we plot log |D p | as a function of γ for the Q-cloud with ω ≈ 0.99. Parameter γ of the unstable (negative) mode corresponds to a sharp dip in this graph, D p = 0. The mode Ψ = (ξ R , ξ I ) in Fig. 9a is then obtained using Eq. (A.4) and normalization condition c 2 A + c 2 B = 1. Performing the above procedure at different ω, we learn that all Q-clouds in Fig. 3b have precisely one unstable mode with γ > 0, while the Q-balls have none. At the critical point dQ/dω = 0 we obtain γ = 0.