Bubble nucleation and growth in very strong cosmological phase transitions

Strongly first-order phase transitions, i.e., those with a large order parameter, are characterized by a considerable supercooling and high velocities of phase transition fronts. A very strong phase transition may have important cosmological consequences due to the departures from equilibrium caused in the plasma. In general, there is a limit to the strength, since the metastability of the old phase may prevent the transition to complete. Near this limit, the bubble nucleation rate achieves a maximum and thus departs from the widely assumed behavior in which it grows exponentially with time. We study the dynamics of this kind of phase transitions. We show that in some cases a gaussian approximation for the nucleation rate is more suitable, and in such a case we solve analytically the evolution of the phase transition. We compare the gaussian and exponential approximations with realistic cases and we determine their ranges of validity. We also discuss the implications for cosmic remnants such as gravitational waves.


Introduction
A cosmological first-order phase transition may have several observable consequences, such as the generation of the baryon asymmetry of the universe (for a review, see [1]), the formation of topological defects [2], or the production of gravitational waves (see [3] for a review). The estimation of these relics involves computing the development of the phase transition in order to determine quantities such as the bubble wall velocity, the number density of nucleated bubbles, and the distribution of bubble sizes. Depending on the computational method, several approximations are usually necessary. In particular, it is customary to consider a constant wall velocity and specific analytical forms for the nucleation rate as a function of time. Among the latter, the most common assumptions are an exponentially growing rate, a constant rate, or just a simultaneous nucleation. The validity of these approximations depends on the strength of the phase transition.
For a first-order phase transition, the free energy density F is a function of an order parameter φ and has two minima separated by a barrier. The stable phase corresponds to the absolute minimum, while the metastable phase corresponds to a local minimum. For concreteness, we shall consider the case in which at high temperature T we have a stable minimum φ + = 0, corresponding to a symmetric phase. At low temperatures, this minimum becomes metastable, while a broken-symmetry minimum φ − becomes the stable one (we shall use the subscripts ± for quantities corresponding to each of these phases). At the critical temperature T c the two minima of F (φ, T ) are degenerate. Several quantities, such as the energy density, are different in each phase. Therefore, if we assume that the phase transition takes place at T = T c , these quantities are discontinuous functions of T . The discontinuities depend on the jump of φ. Thus, the value of φ − /T at T = T c can be regarded as an order parameter. The phase transition is usually said to be weakly first order if φ − /T ≪ 1 and strongly first order if φ − /T 1. We shall be mainly interested in the case of very strong phase transitions, for which φ − /T ≫ 1.
At T = T c or higher, the bubble nucleation rate per unit volume, Γ, vanishes. Therefore, the phase transition does not occur exactly at T = T c . Below T c , the nucleation rate grows continuously from Γ = 0, and may reach, in principle, a value Γ ∼ T 4 c ∼ v 4 , where v is the energy scale characterizing the phase transition. In many cases there is a temperature T 0 < T c at which the barrier separating the two minima disappears and the local minimum φ + = 0 becomes a maximum (see [4] for a review). At this temperature, the phase transition would proceed through spinodal decomposition rather than by bubble nucleation. Nevertheless, before reaching T = T 0 we would have Γ ∼ v 4 , which is extremely large. Indeed, the development of the phase transition depends crucially on the relation between Γ and the expansion rate H [6], since the number of bubbles N nucleated in a causal volume V ∼ H −3 in a cosmological time t ∼ H −1 will be N ∼ ΓH −4 . While Γ varies significantly, H is roughly given by H ∼ T 2 /M P ∼ v 2 /M P , where M p is the Planck mass. Thus, for T close enough to T c we will have Γ ≪ H 4 and the nucleation of bubbles will be too slow, while for Γ ∼ v 4 we have ΓH −4 ∼ (M P /v) 4 ≫ 1, since in most cases v is several orders of magnitude below M P . Therefore, the phase transition will generally occur at an intermediate temperature between T c and T 0 , such that ΓH −4 ∼ 1, and, due to the rapid growth of Γ, it will end in a time ∆t ≪ t.
For a weakly first-order phase transition, the barrier between minima is relatively small and the temperature T 0 is very close to T c . On the other hand, for stronger phase transitions we have a larger barrier which persists at smaller temperatures or even at T = 0. In such a case there is no temperature T 0 (see [5] for examples of such models), and the phase transition can last longer.
Since in many cases the duration of the phase transition is short, the nucleation rate is sometimes approximated by a constant. This is convenient, e.g., for involved computations such as numerical simulations, in which sometimes it is even necessary to nucleate all the bubbles simultaneously. However, it is the variation of Γ with time what actually determines the dynamics of the transition. In general, the nucleation rate is of the form Γ = A exp(−S), where the quantities S and A are functions of temperature and, hence, of time. Nevertheless, the variation of Γ is dominated by the exponent S. A common approximation is to assume that the duration of the phase transition is short enough to justify the expansion of S to linear order about a suitable time t * . This gives a nucleation rate of the form where Γ * = Γ(t * ) and β = −(dS/dt)| t * = H(T dS/dT )| t * .
In the last equality, the temperature variation dT /T = −da/a is assumed 1 , where a is the scale factor. The quantity β is usually positive. This is because, as the temperature decreases from T c to T 0 , the barrier between minima decreases until disappearing. Hence, the nucleation rate monotonically grows with time and the exponent S decreases. With the approximation (1), the parameter β determines the time scale for the phase transition. Thus, the duration is given by ∆t ∼ β −1 and the typical bubble size is given by where v w is the bubble wall velocity. The latter can be assumed to be a constant if the duration of the phase transition is short enough. The expression (1) is widely used as an approximation for the nucleation rate, either to simplify numerical computations or to obtain analytic results. One might expect that the essential conditions on which it is based, namely, that S(t) can be linearized and that β is positive, are quite generally valid. However, both of these conditions may break down. One case in which this happens is that of a very strong phase transition, such that the barrier persists at T = 0. In such a case, as T descends from T c the nucleation rate will initially grow (as the barrier decreases), but for low enough temperatures Γ will stop growing and will begin to decrease, since the presence of a barrier at small temperatures will prevent the nucleation of bubbles. Thus, in this case Γ has a maximum and S has a minimum at a certain temperature T m . Assuming that the system reaches such amounts of supercooling, the linear approximation for S(t) will not be valid at the corresponding time t = t m . In such a case, the sign of β will be different on each side of t m , and Γ will clearly depart from the form (1).
This approximation may also break down due to reheating. Indeed, the energy density difference between the two phases (latent heat) is released at the bubble walls and reheats the surrounding plasma. Such a temperature increase may cause a decrease of Γ. Nevertheless, the effect will depend on the wall velocity. Indeed, if the phase transition fronts propagate as detonations with supersonic velocities, the fluid in front of them remains unperturbed. Thus, the reheating occurs only inside the bubbles, while in the metastable phase the nucleation rate continues growing due to the adiabatic cooling. The same holds for a runaway wall, i.e., a wall which has not reached a steady state and propagates at almost the speed of light. For smaller wall velocities the bubbles expand as deflagrations. A deflagration wall is preceded by a shock wave which carries the released energy outside the bubble. In this case, the plasma in the metastable phase will be reheated, and the nucleation rate will be suppressed. This is a complex scenario, since the reheating is inhomogeneous, causing an inhomogeneous Γ. Besides, the growth of a bubble will be affected by the latent heat released by other bubbles, which complicates further the dynamics. We shall address this case elsewhere and focus here on very strong phase transitions, for which we generally have detonations or runaway walls.
We shall study the dynamics of such very strong transitions. The main aim of this paper is to assess the limit of validity of the approximation (1), and to give an alternative analytic description of the transition, based on a simple approximation for the nucleation rate beyond this limit. The plan of the paper is the following. In the next section we identify the main aspects of the dynamics which are relevant for the possible cosmic remnants. In Sec. 3 we use a realistic toy model to study the dynamics of the nucleation and growth of bubbles in very strong phase transitions. In Sec. 4 we consider analytic approximations for the nucleation rate. We show that, in the very strong limit, the exponential approximation breaks down and a gaussian approximation is more appropriate. The latter allows to obtain analytical results for the dynamics as much as the former. We analyse the ranges of validity of these approximations. In section 5 we compare the analytic results with the numerical computation, and we discuss the effect of using different approximations on the formation of cosmic remnants. Our conclusions are summarized in Sec. 6.

Cosmological consequences of a very strong phase transition
The most important consequence of a very strong phase transition is probably the generation of sizeable gravitational radiation [7]. Gravitational waves (GWs) are produced by the collisions of bubble walls [8] as well as by the motions caused in the plasma [9] (see also [10,11,12]). The process of bubble collisions is usually computed with the envelope approximation [13], in which bubbles are modeled by thin spherical shells which expand with a constant velocity. This approximation has been used in numerical simulations [13,14,15] in which the bubbles are nucleated at random positions with a rate given by Eq. (1). More recently, an analytic approach using the same approximations was considered [16]. The general result is a GW spectrum which depends on the wall velocity and the constant β. The quantity to be computed is the normalized energy density per logarithmic frequency of GWs, Ω GW (f ) = (1/ρ tot )dρ GW /d(log f ), where ρ tot is the total energy density in the plasma. For v w ≃ 1, the peak frequency and amplitude from the envelope approximation are [14] f where ρ wall is the average energy density which accumulates in the thin shells. This corresponds to the gradient energy of the Higgs field, which is a free parameter in the simulation. Although fluid motions generally produce a larger signal of GWs, bubble collisions become important for runaway walls, since in this case a significant fraction of the energy goes to the wall rather than to the fluid. The effects of the fluid motions could in principle be calculated using the envelope approximation, if one assumes that the energy is concentrated in a thin region [9]. However, this is generally not the case (for a recent discussion, see [15]). Moreover, the fluid motions remain after the bubbles have filled all the space and their walls have already disappeared. Lattice simulations [17] of the field-fluid system suggest that the strongest source of GWs are the sound waves generated in the plasma. A fit to these results is given in Ref. [18]. For the peak of the spectrum we havẽ where ρ fl is the energy density in bulk fluid motions, and d is the average distance between bubble centers. It is important to mention that in Ref. [18]f andΩ are actually given in terms of β rather than d, although in the simulations all the bubbles were nucleated simultaneously and thus the scale is set by the bubble separation. The dependence on β was imposed by the replacement d = (8π) 1/3 v w /β, which corresponds to the analytic approximation we discuss in Sec. 4. Another possible remnant of a cosmological phase transition is a network of cosmic strings or other topological defects. Moreover, these defects can only be formed in a phase transition of the universe. In the simplest case in which a global U(1) symmetry is spontaneously broken, the phase angle θ of the Higgs field takes different values in uncorrelated regions. For a first-order phase transition, these regions are the nucleated bubbles. As bubbles meet and the variation of the phase from one domain to another is smoothed out, a vortex or string may be trapped in the meeting point of several bubbles [19]. This determines a segment of string whose length depends on the bubble size. Therefore, the strings will have the shape of random walks, where the step size depends on the size distribution of bubbles at the end of the phase transition.
The statistical properties of the system of strings were studied using Monte Carlo simulations, in which a large volume is divided in cubic cells of fixed size, corresponding to the correlation length ξ [20]. The result is that about 20% of the total length of strings is contributed by closed loops, while the rest are infinite strings. The number density of loops per loop length is given by dn/dl ∼ ξ −3/2 l −5/2 . For a first-order phase transition, the characteristic length ξ is roughly given by the average distance d between bubble centers. More realistic simulations of the phase transition were performed in 2+1 dimensions, nucleating bubbles at random positions and at random times [21,22]. The resulting number of defects per unit volume depends on the number density of bubbles and the wall velocity. Numerical simulations for gauge symmetry braking were also performed (see e.g. [23]). The probability for trapping a defect actually depends on the dynamics of the Higgs phase θ. This dynamics is important for slow bubble walls, since the phase difference between two bubbles could equilibrate before a third bubble meets them. This results in the formation of a smaller number of defects for smaller velocities. On the other hand, for very strong phase transitions we generally have v w ≃ 1, in which case the dynamics of phase equilibration can be neglected (see, e.g., [24]).
Other cosmological consequences of the phase transition depend on similar aspects of the dynamics. For instance, the formation of the baryon asymmetry of the universe in the electroweak phase transition (electroweak baryogenesis) depends strongly on the velocity of bubble walls. Furthermore, in this case, variations of the wall velocity may lead to baryon inhomogeneities [25,26]. The shape of these inhomogeneities will be similar inside each bubble, and their size will be determined by the bubble radius. We shall not discuss this possibility here, since the mechanism of electroweak baryogenesis generally requires small velocities, v w ∼ 10 −2 -10 −1 [27], which are not likely in the case of very strong phase transitions (see, however, [28,29,30,31]). In any case, all these examples show that two aspects of the phase transition dynamics are relevant for the formation of cosmic remnants, namely, the wall velocity and the distribution of bubble sizes.
A complete calculation of the remnants involves the whole evolution of the phase transition. Notice that, in order to take into account the global dynamics (i.e., the expansion and collisions of many bubbles) in a realistic way, the dynamics of bubble nucleation is generally simplified. This is accomplished by considering specific forms for the nucleation rate, which are not necessarily realistic. Thus, for lattice simulations it is usual to either nucleate all the bubbles simultaneously or at a constant rate, while in less time-consuming simulations the exponential approximation (1) is often used.

Phase transition dynamics
In order to investigate the general features of the phase transition dynamics, we shall consider a simple free energy depending on a few parameters. For several considerations it actually suffices to use simple functions F + (T ) and F − (T ) to describe the free energy density of each phase, which in a real model are given by the corresponding minima, F ± (T ) = F (φ ± , T ). The pressure is given by p ± = −F ± and the energy density by ρ ± = T s ± − p ± , with s ± = dp ± /dT . Thus, these functions define the equation of state (EOS) in each phase. However, a realistic treatment of bubble nucleation requires specifying the whole dependence of F on the order-parameter field φ. Therefore, we will consider a free energy density of the form where the first term is the field-dependent part, which vanishes for φ = 0, the second term corresponds to massless particles which do not interact with φ (with a total number of effective degrees of freedom g * ), and the constant ρ V is the vacuum energy density for φ = 0. Therefore, the total energy density of the phase + is given by is the radiation energy density. The Hubble rate will be essentially determined by the quantity ρ + , while the nucleation rate will depend only on the φ-dependent part V (φ, T ).

Effective potential and nucleation rate
The simplest physical model which may provide very strong phase transitions is given by a finite-temperature effective potential of the form For A = 0, Eq. (7) has the well-known form of the high-temperature expansion of the one-loop effective potential, for a tree-level potential of the form V 0 = −m 2 φ 2 + λ 0 φ 4 (see, e.g., [32,33]). It is well known that tree-level modifications of this potential can easily give strong phase transitions, and a simple example is the cubic term −Aφ 3 . The high-temperature minimum of V is φ + = 0, the low-temperature minimum is given by and the critical temperature is given by The strength of the phase transition is usually measured by the quantity φ c /T c , with A further quantity of interest will be the latent heat The strength of the phase transition is directly related to the cubic term of the potential. Indeed, for A = 0 and E = 0 we have T c = T 0 , φ c = 0, and the phase transition is second order. In this case, for T < T c the minimum moves continuously to φ − (T ) > 0, while the point φ = 0 becomes a maximum. For A > 0 or E > 0 the phase transition is first order. In this case, at T < T 0 the point φ = 0 is a maximum and the only minimum is φ − . At T > T 0 the maximum becomes a minimum, and the two minima coexist up to a temperature T 1 , where the minimum φ − disappears. The critical temperature T c lies between T 0 and T 1 . As can be seen in Eq. (10), for A = 0 the quantity φ c /T c is proportional to the parameter E. This means that in this case the first-order nature of the phase transition is due to a fluctuation-induced cubic term, which is not present in the tree-level potential. In a realistic model, E depends cubically on the couplings of bosons with φ, and it is difficult to obtain a relatively strong phase transition (i.e., one in which φ c /T c > 1). In contrast, Eq. (10) shows that adding a tree-level cubic term increases the order parameter, and we obtain a strongly first-order phase transition for A/T c ∼ λ. Furthermore, runaway walls are only possible if the first-order nature of the phase transition is due to zero-temperature terms [34].
It is convenient to relate the parameters of the potential to physical quantities. In the first place, the zero-temperature minimum v = φ − (T = 0) is given as a function of the parameters by Eq. (8). If we want to vary the shape of the potential without changing the characteristic scale of the theory, then we must fix the value of v. This sets the value of the spinodal-decomposition temperature, 2 Another useful relation is φ 2 This temperature decreases for increasing A/v, indicating that the barrier between minima persists at lower temperatures. Notice that for A/v > λ/3 we have T 2 0 < 0, which indicates that in this case the potential has a barrier even at T = 0. Hence, large values of A/v will imply large amounts of supercooling. On the other hand, the value of λ should be such that the curvature of the potential is naturally given by the scale v. Thus, we may choose as a physical parameter the scalar mass and set m φ ∼ v. The value of the dimensionless parameter E, as well as the ratio A/v may be determined by choosing a definite value for φ c /T c (e.g., demanding a strongly first-order phase transition). The value of A/v has a stronger effect on the dynamics, as it affects also the spinodal temperature T 0 , which vanishes for A/v = (m φ /v) 2 /3. Finally, the parameter D can be set as a function the latent heat using Eq. (11). The constant ρ V in Eq. (5) is determined by the condition that the vacuum energy vanishes in the true vacuum at T = 0, i.e., ρ V = −V (v, 0). This gives We see that for A/v > λ/2 the false vacuum energy becomes negative. This actually means that for large A/v the minimum φ = 0 becomes stable even at zero temperature, while the minimum v becomes metastable. Indeed, for A/v = λ/2 the critical temperature vanishes. This sets an upper limit The probability of bubble nucleation per unit volume per unit time is of the form [35] where S(T ) is an instanton action and A(T ) is a dimensional determinant. In the limit of very low temperatures the nucleation rate can be approximated by the false-vacuum decay (quantum tunneling) result, in which S is the four-dimensional instanton action S 4 and the prefactor is given by A ≃ v 4 [S 4 /(2π)] 2 [36]. The action is given by where an O(4)-symmetric configuration is assumed. The instanton corresponds to the extremum of this action and is given by the equation with the boundary conditions dφ/dr| r=0 = 0, lim r→∞ φ(r) = 0. At higher temperatures, the tunneling is affected by thermal fluctuations. In the high-temperature limit we have S = S 3 (T )/T and the prefactor in Eq. (15) is given by A ≃ T 4 [S 3 (T )/(2πT )] 3/2 , where S 3 is the three-dimensional instanton action [37]. Assuming spherical configurations, we have and the configuration of the nucleated bubble is given by the extremum of S 3 . Its equation is similar to Eq. (17), and the same boundary conditions apply. It is usual to consider the smaller of S 4 and S 3 /T as an approximation for S(T ). Analytic approximations for S are sometimes considered (see, e.g., [32]). These are useful to examine the parametric dependence as well as the behaviour as the temperature varies. It can thus be seen that S diverges at T = T c , where the minima of the potential are degenerate, and therefore the nucleation rate vanishes. Furthermore, S vanishes at T = T 0 , where the potential barrier disappears, and hence near this limit we have Γ ∼ T 4 ∼ v 4 . Due to the rapid variation of S with T and the exponential dependence of Γ, analytic approximations generally introduce large errors in the nucleation rate. A better approximation is to consider an expansion of S around a certain temperature as in Eqs. (1)(2), computing numerically S and its derivatives at that point. We shall compute S by solving Eqs. (17) and (19) with the undershoot-overshoot method and then integrating numerically the solution to obtain S 3 and S 4 .
The dimensionless quantity S does not depend on the scale but only on the form of the potential. For the model (7) this can be seen by making the change of variables R = vr, Φ = φ/v in Eqs. (16)(17)(18)(19). With this rescaling, it is easy to see that S depends only on the dimensionless parameters D, E, λ and the ratios A/v, T /v. For concreteness, in this paper we will fix the values of m φ /v, E and D, and we will vary only the parameter A/v, which plays a more determining role for the strength of the transition. We will set m φ = v/2 (similar to the electroweak relation) and we will choose a value of E such that we already have φ c /T c = 1 for A = 0. This gives λ = 1/8 and E = 1/16. Finally, we set D ≃ 0.44, which corresponds to a natural value for the latent heat 3 , L ∼ T 4 c . In Fig. 1 we plot S as a function of T /T c . The rightmost curve corresponds to A = 0, and we see that for this kind of phase transition (without tree-level cubic term) the temperature T 0 is very close to the critical temperature. Thus, we have an extremely rapid variation. The rest of the curves correspond to increasingly stronger phase transitions. The spinodal temperature T 0 decreases as A is increased, and for A = v/12 (dashed line) we have T 0 = 0. For higher values of A/v the potential barrier persists at zero temperature, and therefore the action S 3 never vanishes. In this case, S 3 /T diverges at T = 0 (black solid curve). This indicates that, for very low temperatures, bubble nucleation by thermal activation becomes impossible. In that case, the nucleation may occur by quantum tunneling, and we show also the action S 4 (blue solid line).
Although S(T ) does not depend on v, the dynamics of nucleation does depend on the scale through the expansion rate, which determines the rate at which the universe cools where M P = 1.22 × 10 19 GeV is the Planck mass. As already discussed, bubble nucleation will become appreciable when Γ ∼ H 4 . For T ∼ v, the value Γ = H 4 corresponds to S = S H ≃ 4 log(M P /v). This value is generally large, provided that v is a few orders of magnitude below the Planck scale. Due to the logarithmic dependence, for the wide range of scales 10 8 GeV > v > 1MeV we have 100 S H 200 (see the horizontal green lines in Fig. 1). For a given scale, the value S = S H is reached at smaller temperatures for stronger phase transitions. If T 2 0 ≥ 0 (dashed-dotted, dotted, and dashed curves), then S 3 /T becomes arbitrarily small as T decreases, and the condition Γ ∼ H 4 will eventually be reached. In contrast, when S 3 /T has a minimum (as in the black solid line), if the condition Γ ∼ H 4 has not been fulfilled by the time this minimum is reached, then it will never be fulfilled (Γ never gets high enough). Indeed, notice that in such a case S 4 (blue solid line) is higher than S 3 /T until the latter has become very high, and therefore the quantum tunnelling is also very suppressed 4 . The universe will thus remain stuck in the false vacuum and enter an inflationary stage. We are not interested in such a scenario 5 , and from now on we shall only consider cases in which the nucleation action is given by

Bubble growth and global dynamics
Once nucleated, a bubble begins to grow due to the pressure difference between the two phases, p − (T ) −p + (T ) = F + (T ) −F − (T ). This initial driving force is affected by the local reheating (see, e.g., [40,41,42,43,44,45,46,47]), and thus depends on the values T + and T − of the temperature on each side of the wall (the wall width is in general microscopic in comparison with fluid profiles). Analytic expressions can be obtained by using the bag EOS, as an approximation for the free energy density in each phase. The model (5) has the form of the bag EOS for the + phase, with ε + = ρ V and a + = π 2 g * /90. For the − phase, the approximation (21) can be used by defining temperature-dependent parameters ). Using a linear approximation for the temperature inside the wall, we obtain the driving force [48] where ∆ε = ε + − ε − . The temperatures T ± are related to the speeds of the incoming and outgoing flows in the reference frame of the wall, which we denote v + and v − , respectively. We have (see, e.g., [49]) and where α T = ∆ε/(a + T 4 + ). For 0 < v ± < 1, the relation (24) has four branches. In the first place, the signs in front of the square root give two values of v − , either supersonic or subsonic, for a given v + . Similarly, there are two possible values of v + for a given v − . Notice that the argument of the square root is positive only for large enough or small enough v + . Since we are interested in very strong phase transitions, we need to consider the branch which is compatible with fast walls, namely, that with large v + and large v − . These hydrodynamic solutions are called weak detonations. In this case, the fluid outside the bubble is not affected by the motion of the wall. Thus, the incoming flow velocity v + is given by the wall velocity v w , and the temperature T + is given by the initial condition outside the bubble. These fluid equations are valid for a steady state wall 6 propagating at constant velocity, as well as for a runaway wall which has reached the ultra-relativistic limit v w ≃ 1.
Besides the macroscopic effects of reheating, the wall motion is affected by microphysics inside the wall [50,51,52,53,54,55,56,57,58,59]. The departures from equilibrium of the particle distributions result in a friction force. This force is well understood in the non-relativistic (NR) and ultra-relativistic (UR) limits. In the NR case, the friction grows linearly with the velocity, while for v w → 1 it approaches a constant value [34]. We shall use the approximation [60] which interpolates between these regimes and depends on two free parameters η N R and η U R . Here,v is the average fluid velocity in the wall frame,v = (v − + v + )/2. For small velocities Eq. (25) gives The NR friction coefficient η N R depends on the diffusion of particles which interact with the wall. It is not directly related to the effective potential, and we shall consider it as a free parameter. On the other hand, η U R can be obtained from is the UR limit of the driving force (22) and F U R net is the UR limit of the net force. The latter is given by the mean field value of the effective potential [34]. For the model (7) we have Notice that only for A > 0 we can have F U R net > 0. If F U R net is negative, it means that in the UR limit the friction is higher than the driving force. In that case the UR regime cannot be reached.
The steady state wall velocity is given by the equation F dr = F fr , and can be computed using Eqs. (22)(23)(24)(25)(26). Like for S(T ), the result does not depend on the scale. The wall velocity is a function of the dimensionless ratios η N R /T 4 and η U R /T 4 , as well as of the parameter combination α T . The latter is essentially given by the ratio of latent heat to radiation energy density, α T ∼ L/ρ R , and contains the information on hydrodynamics. For F U R net > 0, these equations may give v w > 1, indicating that the steady state condition cannot be fulfilled and we have a runaway wall. In that case, we must set v w = 1, and the value of F dr − F fr will match the UR limit (26).
This calculation of the wall velocity is valid for an isolated bubble. When bubbles merge to form larger domains, the shapes of the walls, as well as the fluid profiles, change. Although the above equations do not depend on the wall shape 7 , the fluid profiles do [49], and this may also change the boundary conditions. Nevertheless, for detonations the temperature T + is not affected by the wall motion (in contrast, for deflagrations the + phase would be reheated). Inside the bubbles we will have inhomogeneous reheating, but we do not expect this to have a significant effect until the end of the phase transition. Thus, T + is only determined by the adiabatic expansion and we have T + ≡ T , with T (t) given by Eq. (20). Notice that the scale factor is given by a similar equation, since the square root in Eq. (20) is just the expansion rate H ≡ȧ/a. Actually, in the Friedmann equation, the average energy densityρ for a large volume containing the two phases should be used. Nevertheless, the energy difference between the two phases goes into reheating and fluid motions in the − phase, and the average energy density in that phase is approximately given by ρ + (T ).
At time t, the radius of a bubble which nucleated at time t ′ is given by where we have neglected the initial radius, which is a good approximation if the phase transition takes place a few orders of magnitude below the Planck scale. The fraction of volume occupied by the + phase at time t is given by [61,62] f with where t c is the time at which T = T c . The exponential in Eq. (28) takes into account bubble overlapping, and the scale factors in Eq. (29) take into account the dilution of nucleated bubbles. Since bubbles only nucleate in the + phase, the number of bubbles per unit time per unit volume which are nucleated in a large volume containing regions of both phases is given by the average nucleation ratē Thus, as bubbles fill all space, their nucleation turns off. The number density of bubbles at time t is given by We shall estimate the average distance between centers of nucleation by while the mean radius is given bȳ The distribution of bubble sizes is given by where t R is the time at which the bubble of radius R was nucleated, which is obtained by inverting Eq. (27) for t ′ as a function of R and t.
The temperature at which the phase transition "effectively" begins should be close to the temperature T H at which the equality Γ = H 4 is fulfilled, which, as already mentioned, is roughly given by the value S H ≃ 4 log(M P /v). More precisely, T H and S H are determined by the equation (with S = S 3 /T ). However, a sensible definition of the beginning of the transition should involve the dynamics of nucleation. It is usual to define the "onset of nucleation" by demanding that there is already one bubble in a Hubble volume. We thus define the corresponding time t N and temperature T N by the condition nH −3 = 1. Using instead the fraction of volume in the + phase as a measure of the progress of the phase transition, we may define the initial time as that t I at which the variation of this quantity becomes noticeable, e.g., by the condition f + (t I ) = 0.99. Below we shall compare these definitions.
With the quantity f + (t) we may define a few more reference points in the evolution of the phase transition. A representative moment at which there are already many bubbles in contact is when the percolation threshold is reached. This is the time at which there is an infinite cluster of connected bubbles. For spherical bubbles of equal size R and number density n, this happens when 1 − exp(−I) ≃ 0.29, with I = (4π/3)R 3 n [63]. Therefore, we define the percolation time t P by f + (t P ) = 0.71. Another characteristic time which is often considered is the time t E at which the fraction of volume in the false-vacuum phase has fallen to f + (t E ) = 1/e. Equivalently, we have I(t E ) = 1. Notice that the value f + (t) = 0, corresponding to the completion of the phase transition, is approached asymptotically for t → ∞ in the approximation (28)(29). However, in general f + becomes very small in a time of the order of t E − t I . Physically, this indicates that the phase transition actually completes. We shall define a representative time t F at which f + has become small by f + (t F ) = 0.01. We shall denote the values of the various quantities at all these reference points with corresponding subscripts H, I, N, P , E and F . The dynamics of the phase transition will depend mainly on the parameters of the effective potential 8 V (φ, T ). The Hubble rate and the wall velocity depend also on the radiation energy density, and hence on the number of degrees of freedom of the plasma. Besides, H depends on the scale. For concreteness, we shall set g * = 100 and v = 250GeV (i.e., electroweak scale values; as discussed above, for other scales there will not be qualitative differences). For the wall velocity we need to consider a specific value for the non-relativistic friction parameter, and we shall set η N R = 0.25T 4 c . The precise value of η N R is not very important for the strong phase transitions we are interested in, since the friction will be dominated by the UR coefficient, which is determined by Eq. (26).
In the left panel of Fig. 2 we show the order parameter at some reference temperatures, for a range of A/v in which the phase transition is very strong. We have φ c /T c 3, and we see that by the onset of nucleation the order parameter is significantly larger. This is due to a combination of two effects. Firstly, the minimum φ − is closer to the zerotemperature value (which is larger than φ c ), and secondly the temperature has decreased from T c . We also see that the strength of the transition increases very rapidly with A/v, so for A/v ≃ 0.13 the transition becomes so strong that T N vanishes. In the right panel of Fig. 2 we show the wall velocity. In this parameter range we have detonations and runaway walls, while for lower values of A/v we have deflagrations. Thus, there are detonations for only about 15% of the total range of A/v and runaway walls for only 8% of the range. These proportions are typical (see, e.g. [5]), although they depend of course on the rest of the parameters. Notice that for this parameter range the wall velocity does not change significantly during the development of the transition. This happens because the minimum of S is surpassed very soon after t = t H , and the nucleation rate starts to decrease even before there is one bubble per Hubble volume. In spite of this, we see that the various reference points are reached (the point F lies out of the range of the figure, T F ≃ 0.1).
In Fig. 4 we show the evolution of some quantities for a few of the cases considered in Fig. 3. We also indicate the reference times. The plots on the top correspond to A/v = 0.11. The significant amount of supercooling can be appreciated in the fact that the time from t c to t N is quite longer than the time from t N to t F , in which the phase transition actually develops. We also see (left panel) that the interval [t I , t F ] in which the fraction of volume varies in also quite shorter than the latter. Furthermore, most bubbles nucleate in this short interval near the end of the transition. In the top right panel of Fig. 4, we see that the average distance between bubble centers and the average bubble size become comparable by the end of the transition, as expected. The bubble separation has initially a very rapid variation, since Γ is increasingly high, but becomes constant at t = t F , as bubbles cease to nucleate. We also see that the average radius grows very slowly during most of the phase transition, even though all bubbles grow with a velocity v w ≃ 1. The small average is caused by the constant nucleation of vanishingly small bubbles. This effect disappears for t ≃ t F , as the space where bubbles nucleate is considerably reduced, and we have dR/dt ≃ v w . Although we will be mostly interested in higher values of A/v, it is worth mentioning that the evolution will be qualitatively similar for weaker phase transitions, provided that bubbles expand as detonations. Quantitatively, the variation of f + andΓ will be even more concentrated at the end of the phase transition, since the function S(T ) will have a steeper slope (see Figs. 1 and 3). Thus, if T N is close to the critical temperature, the nucleation rate will grow at a dramatically increasing rate. In such a case, the dynamics is dominated by the nucleation of bubbles, rather than by their growth. In contrast, for stronger phase transitions we have more supercooling, the temperature gets closer to the minimum of S, and the nucleation rate varies more slowly.
In the left panel of Fig. 5 we plot the distribution of bubble sizes for the same case A/v = 0.11, at the times t I , t P , t E and t F . We see that, for t ≤ t E , vanishingly small bubbles are more abundant. This is due to the rapid growth of the nucleation rate. On the other hand, for t > t E the maximum of the distribution separates from R = 0, indicating that the average nucleation rate has decreased. Notice that the width of this distribution is essentially the same as that ofΓ, as expected from Eqs. (30) and (34). We also indicate in Fig. 5 the radius of the "largest" bubbles, R(t N , t i ), at each reference time t i . We see that, for this case, the average size is generally quite smaller than the largest bubble.
The second row of plots in Fig. 4 corresponds to the case A/v = 0.126. In this case we have an important amount of supercooling (S is quite close to the minimum, as can be seen in Fig. 3). Therefore, we have a larger time t H − t c . Besides, due to the slower growth of Γ the development of the phase transition is slower, and at the end of it we have fewer and larger bubbles. On the other hand, the variation of f + begins relatively sooner (i.e., closer to t N in the interval [t N , t F ]). As a consequence, the curve ofΓ is less concentrated. This is reflected in the bubble size distribution (see the central panel of Fig. 5), which is wider in relation to the largest bubble size R(t N , t). Due to the longer phase transition, the effect of the scale factor begins to be noticeable in the curve of d as well as in the evolution of dn/dR. The third row of Fig. 4 and the right panel of Fig. 5 correspond to the case A/v = 0.128. In this case the minimum of S is crossed during the phase transition, and the nucleation rate turns off even in the + phase. Thus, the effective rateΓ drops to zero more quickly than f + . Besides, since Γ reaches a maximum, it varies much more slowly than in the previous cases. This can be seen e.g. in the curves ofΓ and the bubble size distribution, which are quite wider. The dynamics is no longer dominated by bubble nucleation, but also depends on the bubble growth and on the expansion rate of the universe. The effect of bubble growth can be seen in the curve ofR, which reaches the behavior ∝ v w t sooner than in the previous cases (i.e., closer to t N ). This can be appreciated also in the evolution of the bubble size distribution, which at t = t P has already a non-vanishing maximum. The effect of the expansion of the universe is reflected, e.g., in the distance d between bubble centers. As soon as bubble nucleation slows down, d stops decreasing and begins to grow. This shows that the variation of the scale factor during the phase transition is not negligible. The effect is more evident in the distribution of bubble sizes, where the height of the curve of dn/dR decreases with t due to the dilution of the number density as n ∝ a −3 .
Although these three examples correspond to very strong phase transitions, we see that the dynamics becomes qualitatively different near the minimum of S. The difference is particularly evident if we compare the leftmost and rightmost panels in Fig. 5. The bubble size distribution is not only quantitatively and qualitatively different, but it also evolves quite differently.
In discussing GW generation it is more appropriate to look at the volume-weighted bubble distribution, R 3 dn/dR, since the energy injected into the walls of a bubble and into the surrounding fluid is proportional to the volume of the bubble. In Fig. 6 we show the volume-weighted distribution together with the radius distribution, for the cases A/v = 0.11 and A/v = 0.128. For a better comparison, we have divided the distributions by their integrals. The volume-weighted distribution vanishes at R = 0, and therefore is peaked at a larger radius. For this distribution we still observe a different behavior in the

Analytic approximations
Analytic approximations are useful for a better understanding of the physics as well as for specific applications. In particular, analytic approximations for the nucleation rate are often used in simulations of the phase transition. If the duration is short enough, then the linearization of S(t) is a good approximation and we may use a nucleation rate of the form (1) [64,6]. In this approximation, the constant β −1 will set the time scale for the dynamics. One expects that the approximation will be good if this time scale is shorter than that for the temperature variation, which is given by H −1 . Therefore, the approximation will be valid for β/H ≫ 1. More concretely, consider the second-order expansion of S(t) about a certain t * , and all the quantities are evaluated at t = t * . The linear approximation will be valid if we can neglect the last term in (36), which implies the conditionṡ The first condition requires a small variation of H in the time (t − t * ), while the second condition is just the requirement of a short time compared to H −1 . Since in general H/H ∼ H, these are equivalent. Assuming that they are fulfilled, the third condition requires that the coefficient β = −dS/dt does not change too much if it is evaluated at the time t instead of t * . Indeed, in a time ∆t we have a variation ∆β ≃ 2α 2 ∆t. Hence, the last condition in (38) just demands that ∆β ≪ β.
For the exponential rate to be a good approximation, Eqs. (38) should be valid at least for the characteristic time t − t * ∼ β −1 . Therefore, we have The second condition does not depend on H and imposes a relation between the temperature derivatives of S(T ). As we have seen, the phase transition in general takes place for high values of S, namely, S ∼ 4 log(M P /v) ≫ 1. Hence, for natural values of the derivatives we have T dS/dT ∼ T 2 d 2 S/dT 2 ∼ S, which give β/H ∼ S, α/H ∼ √ S. Both quantities are generally large, and so is the ratio β/α ∼ √ S. Therefore, the two conditions (39) will be fulfilled in most cases. Of course, there is a special case in which this argument breaks down, namely, when the temperature gets too close to the value T m corresponding to the minimum of S(T ), at which β vanishes. For the case considered in Figs. 3 to 6 (v = 250GeV) we have S ∼ 150. Table 1 shows the values of β/H at several times t * for the specific cases considered in Figs. 4 to 6. The value of α/H at the minimum is also shown. In all these cases we have α/H ∼ 10, while β/H is roughly ∼ 10 2 . However, the latter will be higher for weaker phase transitions, and it becomes small and negative for t * very close to the minimum (cf. Fig. 3).  Table 1: The values of β/H at the reference points, and the value of α/H at the minimum of S, for the cases considered in Fig. 4.

A/v β/H| T H β/H| T N β/H| T I β/H| T P β/H| T E β/H| T F α/H|
As T approaches T m , the second condition in Eq. (39) will be violated before the first one (since in general we have α ≫ H). Thus, the linear approximation will break down when the value of β(T * ) becomes comparable to that of α(T * ). This establishes a rough criterion for the validity of this approximation, namely, β(T * ) α(T * ), and this condition will only break down near the minimum. Thus, we expect the exponential rate to be a good approximation for β(T * ) a few α(T m ).
If T gets too close to T m , we may set t * = t m in Eq. (36), and we obtain a gaussian approximation for the nucleation rate, with Γ m = exp[−S(T m )]. This quadratic approximation for S(t) will also break down for t far enough from t m . In particular, we expect that it will break down for |t − t m | H −1 . Therefore, if the temperature T m is not reached during the phase transition, in order to use the approximation (40) we must have at least S(T ) close enough to S(T m ). On the other hand, if temperatures below T m are reached, then Eq. (40) can be safely used, even if the phase transition takes a long time to complete, since Γ will be negligible after a time ∼ α −1 , and, as we have seen, α −1 ∼ ( √ SH) −1 ≪ H −1 .

Exponential nucleation rate
As can be appreciated in Fig. 3, in many cases the function S(T ) is approximately linear during the phase transition or at least in the relevant interval [t I , t F ] (where f + varies and the size distribution is formed). This can also be seen in the first row of table 1, where β varies by a 10% between t I and t F . In such cases, the exponential nucleation rate (1) should be a good approximation. Besides, as we have seen, we may also assume a constant wall velocity. With these assumptions Eqs. (27)- (34) can be solved analytically [6]. We shall also neglect the evolution of the scale factor for simplicity, although this is not necessary to obtain analytic results. Besides, when comparing a time interval with the characteristic cosmic time H −1 , we shall assume a constant H in that interval, which is consistent with a constant β. We shall later compare these approximations with the numerical results obtained above. With these approximations the radius is just given by R(t ′ , t) = v w (t − t ′ ), and taking, with little error, t c → −∞ in Eq. (29), we obtain [6] where I * = 8πv 3 w Γ * /β 4 . Thus, we have I(t) = (8πv 3 w /β 4 )Γ(t). With f + = e −I , this gives which is easily integrated to obtain the number density, Actually, at the beginning of the transition we may ignore the factor e −I in Eq. (42), and we obtain (for f + ≃ 1) n = Γ(t)/β.
The distribution of bubble sizes is given by and for the average radius we have the expression In the last equality we have used the relationΓ = dn/dt ′ . To see the qualitative behavior of this approximation, it is convenient to set v w = 1 and choose the value Γ * = H 4 (corresponding to t * = t H ), so that, if time and distance are in units of H −1 , all the above expressions depend on the single parameter β/H. In Fig. 7 we plot some of these quantities for β/H = 100. As expected, these plots show the qualitative behavior of the first row of Fig. 4, even though the present value of β(t H ) is closer to that corresponding to the third row (see table 1). There is a qualitative difference, though, between Fig. 4 and Fig. 7, namely, that for most of the phase transition the average radius remains constant in the latter, while in the former there is a small slope. This difference is mainly due to the variation of β, which is neglected in the present approximation.
More quantitatively, at the beginning of the phase transition we haveΓ ≃ Γ and the first integral in Eq. (46) givesR(t) = v w /β. Moreover, this result does not change to first order in I. Indeed, from Eq. (43), the second integral in (46) . This givesR(t) = v w /β + O(1 − f + ) 2 , which will remain constant until f + departs appreciably from f + = 1. At first sight, such a constant average seems in contradiction with the fact that, according to Eq. (45), the bubble size distribution moves with velocity v w without changing its shape. However, part of this curve falls into negative sizes. It is easy to see that the peak is at R p (t) = v w (t − t E ). Therefore, we have R p (t) < 0 for t < t E . This means that before t = t E the physical distribution actually decreases from R = 0 (in agreement with the left panel of Fig. 5). As a consequence,R will remain small until t = t E (and by this time the phase transition is fairly advanced). On the other hand, for t t F we haveΓ ≃ 0 and the first integral in (46) givesR ≃R(t F ) + v w (t − t F ), as observed in the figures.
In contrast, since the volume-weighted distribution vanishes at R = 0, its average radius is close to its maximum R ′ p . The latter is given by R ′ p = x(t) v w /β, where x is the solution of the equation (x−3)e x = xI(t). This peak is initially at R ′ p = 3 v w /β (for I = 0) and moves to R ′ p ≃ 3.14 v w /β at t = t E , while at t = t F we have R ′ p ≃ 3.49 v w /β. This small variation is in agreement with the left panel of Fig. 6. Another important distance scale is the bubble separation d = n −1/3 . According to Eq. (44), d initially reflects the rapid variation of Γ. However, near the end of the transition n reaches a constant value n = β 3 /(8πv 3 w ), and we have d = (8π) 1/3 v w /β. For the exponential rate it is straightforward to obtain the relations among the reference times defined in the previous section. At t = t H we have f + = exp(−8πv 3 w H 4 /β 4 ), which in the usual case H/β ≪ 1 may be written as f + (t H ) ≃ 1 − 8πv 3 w (H/β) 4 (in most cases the second term is actually negligible). Then, we have n(t H ) ≃ H 4 /β, which indicates that the number of bubbles per Hubble volume is very small, n(t H )H −3 ≃ H/β ≪ 1. Therefore, the time t H is always earlier than the onset of nucleation. At the time t N we may use Eq. (44), which gives Comparing with Γ(t H )/H 4 = 1, we obtain t N − t H ≃ β −1 log(β/H). (48) Notice that at the time t N we have I(t N ) = 8πv 3 w (H/β) 3 ≪ 1, i.e., f + ≃ 1. On the other hand, the subsequent times t i , i = I, P, E, F , correspond to specific (larger) values I(t i ) = I i , with I i = 0.01, 0.34, 1, and 4.6. We thus have Γ(t i ) = I i β 4 /(8πv 3 w ). Comparing with (47), we have, for example, For the last four points the time differences are given by t i − t j = β −1 log(I i /I j ). For instance, for the interval in which f + varies we have If we take the width ofΓ(t) as t F − t I , then the width of the distribution of bubble sizes is ∆R ≃ 6v w β −1 . We see that all the time intervals are of order β −1 , but for the time from t H to t I we have log(β/H) enhancements. In contrast, the time in which the fraction of volume varies is fixed in units of β −1 . This difference becomes important for β/H ≫ 1, i.e., for weaker phase transitions. In the interval [t I , t F ], the time t E is particularly interesting, since it corresponds to the moment of maximum bubble nucleation, as can be easily seen by differentiating Eq. (42). This agrees with the first case considered in Fig. 4, and is also in good agreement with the second case.

Gaussian nucleation rate
The gaussian approximation (40) will be more appropriate to describe a phase transition like the one shown in the last row of Fig. 4. Neglecting the variation of the wall velocity, the evolution of the phase transition can be solved again analytically, although the expressions are more involved. For the moment we will also ignore the variation of the scale factor. At the beginning of the phase transition the number density of bubbles is just given by the integral of Γ, and we have (for f + ≃ 1) with The function (1 + erf)/2 varies between 0 and 1, hence n max is the maximum density of bubbles that can be nucleated. For n max H −3 < 1 there will always be, in average, less than a bubble per Hubble volume, and the time t N will not exist. Hence, the existence of t N gives the condition where x = α(t − t m ) and For t t m − α −1 we have I(t) ≃ 0, while at t = t m we have I = I m and for t t m + α −1 we have . This cubic growth contrasts with the exponential behavior in Eq. (41). In Fig. 8  Although the latter reproduce better the behavior of the curves in Fig. 4, we see that the quantitative difference between the two approximations is not significant. The curves in Fig. 8 are plotted from t = t H , which is readily obtained from the nucleation rate (40), This time is shown in Fig. 9 (black solid line). The black dashed line corresponds to the time at which Γ reaches the value H 4 again, after crossing its maximum, and is obtained by changing the sign in the right hand side of Eq. transition takes place around the minimum of S, we will have t m − t H = a few α −1 . This can be appreciated in Fig. 9 (where α −1 ∼ 0.1H −1 ). The time t N is obtained from Eq. (51), and is given by the inverse of the error function, This time is indicated by a red dotted line in Fig. 9, which coincides with the red solid line in most of the range. The latter takes into account the variation of the scale factor (see below). The two approximations only separate for Γ m /H 4 very close to the limit in which t N ceases to exist. Even this lower bound is very similar for the two approximations. If we ignore the scale factor, this limit is given by the condition (53) The average nucleation rateΓ(t) = Γ(t) exp[−I(t)] is an important quantity since it determines the bubble number density as well as the bubble size distribution and, hence, the various size scales. Obtaining, e.g., the bubble separation or the average radius from Eq. (54) requires a numerical computation. Nevertheless, assuming that the nucleation rate turns off before f + decreases, we may just setΓ(t) = Γ(t). This rough approximation is actually quite good in the case of Fig. 8 (as can be seen also in the left bottom panel of Fig. 4), and it is even better for stronger phase transitions. In this approximation all the bubbles nucleate for f + ≃ 1, and the number density is given by Eq. (51). Thus, the average bubble separation achieves a final value d ≃ n −1/3 max ≃ π −1/6 (α/Γ m ) 1/3 . On the other hand, the radius distribution is given by This distribution has a width ∆R ∼ v w α −1 and its peak is given by (which is negative for t < t m ). On the other hand, the peak of the volume-weighted distribution is given by Therefore, at t = t m we already have a positive peak R ′ p ≃ 1.2v w α −1 . This gives a good description of the behavior of the peaks in the right panels of Figs. 5 and 6. As t increases R p and R ′ p approach each other, since for t − t m ≫ α −1 we have R ′ p ≃ v w (t − t m ). This behavior can also be appreciated in Fig. 6.

Gaussian nucleation rate with exponential expansion
If we take into account the variation of the scale factor, we may still obtain analytic expressions, as long as we assume a simple form for a(t) or, equivalently, for H =ȧ/a. For a short enough time ∆t ≪ H −1 this is not relevant, and we may set H(t) ≃ H(t m ). On the other hand, longer times ∆t H −1 only matter in cases in which there is strong supercooling. In such a case the energy density begins to be dominated by the false vacuum (since ρ R /ρ V ∼ T 4 /v 4 ), and the expansion rate will again be approximately constant. Therefore, a constant H should be a good approximation in general for very strong phase transitions, and we have In Ref. [6], such an exponential expansion was considered with an exponential nucleation rate. However, as we have seen, for the cases in which the exponential rate is a good approximation, the duration of the phase transition is generally short enough, so that the scale factor can be safely ignored, i.e., ∆t ∼ β −1 ≪ H −1 . In contrast, if the minimum of S is reached we may have ∆t > α −1 . We thus use the relation (61) together with the gaussian nucleation rate (40) in Eqs. (27)- (34). The initial number density of bubbles (i.e., for f + ≃ 1) is given by with n max defined in Eq. (52). This result is similar to Eq. (51), with the additional exponential factor which takes into account the dilution n ∼ a −3 . For t − t m ∼ α −1 , we have in general H(t − t m ) ∼ H/α ≪ 1 and we recover Eq. (51). Due to this suppression, the maximum number of bubbles which can be nucleated in a Hubble volume is now smaller than n max H −3 . As a consequence, we have a stronger condition than Eq. (52) for the existence of the time t N . Nevertheless, for the typical values H/α ∼ 0.1 the difference is small, as can be seen in Fig. 9 (solid and dotted red lines), and the bound on Γ m /H 4 is still roughly given by Eq. (53). Thus, according to Eq. (62), a total bubble density n ≃ n max is nucleated in a time t − t m ∼ α −1 , and for t − t m α −1 the error function saturates and we have (neglecting the term H 2 /α 2 ) For t − t m ∼ H −1 the exponential suppression becomes important, and eventually the condition nH −3 = 1 will be achieved again (red dashed line in Fig. 9). This second solution only makes sense if it is achieved during the phase transition; once the universe is in the true vacuum the assumption of an exponential expansion will not be a good approximation. For a constant wall velocity, Eq. (27) gives R(t ′ , t) = v w H −1 (e H(t−t ′ ) − 1). Then Eq. (29) gives, defining τ = H(t − t m ), where The result (54) is recovered if we assume t − t m α −1 and expand Eq. (64) to third order in H/α. On the other hand, for t − t m α −1 the expression becomes even simpler, since the error functions saturate and we obtain (neglecting the terms H 2 /α 2 ), A numerical comparison of Eqs. (54) and (64) is given by the orange and purple lines in Fig. 9, corresponding to the times at which the values I = 1 (t = t E ) and I = log 100 (t = t F ) are reached, respectively. The solid lines are obtained from Eq. (64), while the dotted lines ignore the variation of the scale factor. The latter is a very good approximation for t − t m α −1 . On the other hand, for t − t m α −1 we may use the approximation (66), and the equation I(t i ) = I i gives For α/H ∼ 10, this becomes a very good approximation for t−t m > 2α −1 (Γ m < 0.02I i α 4 ). In the range α −1 t − t m H −1 we may use the approximation I i /I max ≪ 1 in Eq. (67), which gives As can be seen by inspection of Eq. (64), and more directly in the approximation (66), for t − t m ≫ H −1 we have a constant I(t) = I max . Formally, this implies that the phase transition never ends, since f + (t) approaches the asymptotic value e −Imax > 0. Nevertheless, notice that even when I(t) grows exponentially, as in Eq. (41), the limit f + = 0 is reached in an infinite time. This is a characteristic of the approximation (28)(29), which lead us to define the time t F by f + (t F ) = 0.01 in Sec. 3. In the present case, for large enough I max the fraction of volume in the false vacuum will reach a very small value within a Hubble time, and this means that the phase transition is actually completed. On the other hand, for smaller values of I max the situation changes. Reaching the value I = I F requires I max > I F , which gives the restriction As this bound is approached, the time needed to reach the value I = I F becomes infinite (and we have a similar bound for any value I = I i ). This is apparent in Eq. (67) and can be appreciated in Fig. 9, since the purple solid line has an asymptote at this value of Γ m /H 4 (and similarly for the orange solid line). It is worth mentioning that, in contrast to the gaussian rate, an exponential nucleation rate Γ ∝ e β(t−t * ) gives an exponentially growing I(t) ∝ Γ(t), even in the case of an exponential scale factor. In spite of this, for small β the completion of the transition will take many Hubble times (see Ref. [6] for details). In this scenario, it is convenient to consider the physical volume in the false vacuum, V phys (t) ∝ a 3 (t)f + (t). This volume always grows at the beginning of the phase transition, and in the case of small β, it will grow during most of the phase transition. Hence, there is a new "milestone", namely, the time t e at which the derivative vanishes. For the exponential nucleation the time t e always exists, but for β H this time is later than t F [6]. In such a case, at t = t F the physical volume in the false vacuum has never stopped growing, and we cannot say that the phase transition is coming to an end. Nevertheless, as we have seen, in general we have β ≫ H.
For the gaussian nucleation rate, even though we have in general α −1 ≪ H −1 , the phase transition can take a longer time. For t − t m α −1 we may use the approximation (66), and we see that dI/dt reaches a maximum at (t − t m ) = log(3)H −1 . The maximum value is (4/9)HI max , so the existence of the time t e requires I max > 27/4, which gives another bound on Γ m , For smaller Γ m , the volume in the false vacuum always grows and the phase transition never completes. Notice that this bound is slightly more restrictive than Eq. (69) (since 27/4 > I F ). For larger Γ m the time t e exists, and it can be checked that it is earlier than t F . Therefore, in this case dV phys /dt becomes negative before t = t F . However, since dI/dt decreases for t − t m > log(3)H −1 , there is a second time t ′ e at which dV phys /dt becomes positive again. Requiring that this does not happen before the time t F gives the condition which is a little more restrictive than (69) and (71). In Fig. 9 the times t e and t ′ e are indicated by solid and dashed cyan lines, respectively. We shall take the condition (72) as the bound for the completion of the phase transition. For larger Γ m /H 4 , when the value f + = 0.01 is reached the physical volume in the false vacuum is decreasing, and the phase transition will soon complete.
For the case of For this case we cannot say that the phase transition will ever complete, even though the value f + = 0.01 is reached. We also see that this value of f + is achieved in a relatively long time, t F − t H ≃ 2.5H −1 . This case (Γ m /H 4 ≃ 11) corresponds to the case A/v = 0.129 of the physical model, and was considered in Fig. 3 (second curve from the top). In that figure, the long duration of the phase transition is apparent from the low temperatures at which the points P and E are reached (the value T F ≃ 0.1 lies outside the range of the figure). Notice also that in Fig. 9 the vertical line barely touches the solid red line. Accordingly, in the numerical computation of Sec. 3 the value nH −3 = 1 is exceeded only during a short time (the maximum number of bubbles achieved is nH −3 ≃ 1.2).
To estimate the bubble size distribution we must now take into account the scale factors in Eq. (34). At a given time t, the time t R at which a bubble of radius R was nucleated is given by t R = t − H −1 log(HR/v w + 1), and the scale factors at these times are related by a/a R = 1 + HR/v w . As in the previous subsection, we shall use the approximationΓ(t) ≃ Γ(t), which is equivalent to assuming that f + (t R ) ≃ 1 at the time in which most bubbles were nucleated. Thus, the size distribution is given by This reduces to Eq. (58) for small HR. The main difference is in the denominator of Eq. (73), which takes into account the dilution of the nucleated bubbles. Thus, for increasing time the peak, R p = v w H −1 (e H(t−tm) −1), moves to higher R, and the amplitude at R = R p decreases as (1 + HR p /v w ) −4 . This is in agreement with the right panel of Fig. 5. We see that for large t the radius growth will be dominated by the exponential expansion, as any physical length. This effect will be noticeable if the completion of the phase transition takes longer than H −1 . For t − t m α −1 but still smaller than i.e., the result (59) has a small modification. Similarly, for the peak of the volume-weighted distribution, to linear order in H(t−t m ) we obtain the result (60). Using the approximation (68), we have On the other hand, the average bubble separation d = n −1/3 reaches a minimum value d ≃ n −1/3 , which is very similar to the values of R p and R ′ p .
5 Application to specific models

Computing the parameters
In order to apply the analytical approximations obtained in the previous section, it is only necessary to compute, for a given model, either the quantities Γ * , β, v w at a suitable temperature T * or the quantities Γ m , α, v w at T = T m . Since T m is just the minimum of S(T ), its computation is relatively simple, as it does not involve the dynamics of the phase transition. On the other hand, there is some freedom in choosing the value of T * . In principle, we might take, e.g., any value between T H and T F , provided that the linearization of S is valid in the whole range. Notice that T * = T H is the simplest choice, since T H is obtained by solving Eq. (35), which does not involve the dynamics of the phase transition (and therefore is as simple as obtaining T m ). However, for a real model the linear approximation for S is valid only in a small time interval, and the parameter β may vary significantly in the interval [t H , t F ] (see Table 1). Therefore, the computations will be more precise if β is evaluated at a later time, such as t * = t E , where most bubbles are nucleated. It is relatively easy to follow the phase transition up to the time t N , since it only requires to evaluate the integral 10 in Eq. (31) withΓ = Γ. Continuing to later times, instead, involves solving numerically Eqs. (27)- (29), which is precisely the calculation one is trying to avoid by using approximations such as (1) or (40). Therefore, it is very common to evaluate β, as well as v w , at T = T N . A relatively simple and accurate approximation for T N can be obtained by using the exponential rate approximation with t * = t N . Although this will not be a good approximation for computing the whole development of the phase transition (which involves times t − t N ∼ 10β −1 ), we can use it to compute the number of bubbles nucleated before t N . Indeed, since Γ falls quickly for t < t N , the integral (31) effectively involves times The result (assuming a constant scale factor and wall velocity in the time ∼ β −1 ) is given by Eq. (47), with β evaluated at t * = t N . In terms of T N we have The same reasoning can be applied for the computation of I(t) (see, e.g., [43]). As long as the time t is far enough from the minimum, t m − t α −1 , the nucleation rate falls quickly for t ′ < t in the integral (29). Then, we may use the exponential rate, and the best choice for t * is t * = t. For instance, to estimate the temperature T E we may set t = t * = t E in Eq. (41), which gives 8πv 3 w Γ(T E )/β 4 = 1, with β evaluated at T E . We thus have the equation The wall velocity in the last term can be ignored for very strong phase transitions, i.e., setting v w = 1 will not introduce a significant error in the value of T E . These equations are very similar to that for T H , and give T N and T E without considering details of the dynamics of the phase transition. The main difference with Eq. (35) is the appearance of a term involving the derivative of S, and the approximations will break down if S ′ (T ) becomes too small. Nevertheless, they will provide good estimations for T N and T E quite close to T m . We have checked this by comparing with the values obtained from the numerical evolution. For instance, in the very strong cases A/v = 0.126 and A/v = 0.128, the difference in T N as given by Eq. (74) or by the numerical evolution is on the order of a 0.1%. On the other hand, the agreement for T E is on the order of a 1% in the case A/v = 0.126, while for A/v = 0.128 the approximation (75) breaks down since S ′ becomes negative. Some applications require to compute certain quantities at a given reference temperature, for which solving Eq. (74) or (75) may be enough. For other applications the time intervals will also be important, and the time-temperature relation may be needed. Notice that estimations of time intervals such as Eqs. (48)(49)(50), which assume a constant β, cannot be used to determine the time elapsed from the critical temperature. If the equation of state is simple enough, the Friedmann equation will give an analytic expression for the time-temperature relation. For instance, for the bag EOS the expansion rate is of the form H = √ a + bT 4 . In this case Eq. (20) can be solved analytically. Using the condition T = ∞ at t = 0, we have [42] T 2 = a/b/ sinh(2 √ at). Our model (5) has the bag form in the + phase 11 , and we have a = (8π/3)ρ V /M 2 P and b = (4π 3 g * /45)M −2 P . Thus, we may write This equation is valid while the system is in the + phase, and thus it can be inverted to obtain the time t c as well as the time t N . Moreover, as already discussed, for a phase transition mediated by detonations or runaway walls the average energy density is given by ρ + (T ) until the end of the transition. Therefore, Eq. (76) is still valid at later times, and we have We remark that Eq. (77) is independent of the nucleation rate and can be used, e.g., to compute the time t m if the temperature T m is reached during the phase transition. After the phase transition there will be some reheating, and Eq. (76) will no longer hold.

Comparison with the numerical computation
We shall now compare the different approximations with the numerical computation. For that aim we shall consider the evolution of the fraction of volume, as well as the normalized distribution of bubble sizes at the time t = t E . Thus, for instance, for the exponential rate we have f Notice that the latter equation does not depend explicitly on t * . Different choices of this time will be reflected in Eq. (78) through different values of β and v w . The parameter I * appearing in f + is also a function of β and v w . For instance, for t * = t N we have I * = 8πv 3 w (H/β) 3 , and for t * = t E we just have I * = 1. On the other hand, the value of t * is related to the corresponding temperature T * by the analytic expression (77). We shall consider the temperatures T * = T N and T * = T E , which can be obtained from the semi-analytic approximations (74) and (75), respectively. We remark that, in these approximations, the only parameters which we shall compute numerically are T * , β(T * ) and v w (T * ). Similarly, for the gaussian rate the results depend only on the numerical values of α(T m ) and Γ(T m ). The function f + (t) depends on the parameter t m , which is again obtained from Eq. (77), while the size distribution depends only on the difference t E − t m , which is obtained from Eq. (54) as a function of α and Γ m .
In Fig. 10 we consider the cases discussed in Sec. 3. The solid lines correspond to the numerical computation. In the top panels we consider the case A/v = 0.11, and we show the results of the exponential rate, with t * = t E (dashed lines) and t * = t N (dotted lines). As expected, using t * = t E gives an excellent approximation, while for t * = t N the approximation is not so good. We do not include the gaussian approximation, since this case is far from the minimum of S.
The central panels of Fig. 10 correspond to the case A/v = 0.126. In the left panel it can be appreciated that, at the beginning of the phase transition, the evolution of f + is better approximated by the curve corresponding to t * = t N , while, as expected, that corresponding to t * = t E is closer to the numerical result in the last stages. On the other hand, the gaussian rate (dashed-dotted line) gives also a reasonably good approximation for f + (t). Moreover, the shape of this curve describes the evolution of the fraction of volume better than the curve of t * = t N . Indeed, for the latter the evolution is too quick due to a large value of β(T N ). In the right panel it can be seen that the case t * = t E and the gaussian rate give quantitatively similar results for the size distribution. However, the latter reproduces more correctly the shape of the curve, with a non-vanishing maximum. Notice that this is a borderline case; for lower values of A/v the exponential rate gives a better approximation, while for higher values the gaussian approximation will be better. As can be seen in the figure, in this limit the error of both approximations is smaller than a 10%. The plots at the bottom of Fig. 10 correspond to the case A/v = 0.128. In this case the gaussian rate is clearly a good approximation, while the exponential approximation breaks down. Indeed, for t * = t E we have β < 0 and the approximation cannot even be used. On the other hand, we see in the left panel that the case t * = t N is a good approximation only for the first stages of the evolution.

Implications for cosmic remnants
As discussed in Sec. 2, numerical simulations for computing certain cosmological consequences of the phase transition require several simplifications on the dynamics. One usual approximation is to assume a constant wall velocity. In the case of a very strong phase transition this is generally a good approximation since we have v w ≃ 1. On the other hand, in simulations bubbles are often nucleated simultaneously, which corresponds to a delta-function nucleation rate, or at random times, corresponding to a constant nucleation rate. The most realistic approximation is an exponential rate of the form (1), whereby the value of Γ * is generally chosen so that the desired number of bubbles is nucleated in a sample volume. This number is restricted by the limitations of the simulation. Nevertheless, the dynamics of the phase transition is determined by the parameter β rather than by Γ * . Thus, the results of the simulation will generally be functions of β, as in Eq. (3). As we have seen, this exponential approximation for Γ breaks down around the minimum of S(T ), and the gaussian approximation is more appropriate. Our numerical and analytic results show that several aspects of the dynamics, such as, e.g., the bubble size distribution, cannot be directly extrapolated from weaker phase transitions to very strong cases. To our knowledge, a gaussian nucleation rate has not been considered so far in numerical simulations (and it is certainly out of the scope of this work). In the gaussian case, the dynamics will depend not only on the parameter α but also on Γ m .
In principle, one expects that the cosmic remnants depend particularly on the bubble size distribution and, hence, on the nucleation rate. For instance, the GWs inherit the wavenumber spectrum of the source [65]. However, recent simulations in the envelope approximation [15] seem to contradict this expectation. In Ref. [15] bubbles were nucleated either using the exponential rate or simultaneously, and the resulting GW spectrum is described in both cases by the same power laws at high and low frequencies. Besides the similarity in the shape of the spectra, the peak frequency and amplitude agree within the order of magnitude, although in the exponential-rate case the frequency is lower and the amplitude higher (see Fig. 3 in [15]). This similarity is somewhat unexpected. Notice in particular that for a simultaneous nucleation the bubble size distribution will peak around the average bubble separation, while for an exponential rate smaller bubbles will be much more abundant 12 . In any case, we remark that for GW generation one expects the volume-weighted bubble distribution to be more relevant. For the exponential rate, the peak of the volume-weighted distribution has a value R ′ p ≃ 3.49 v w /β at the end of the transition, which is not too different from the final bubble separation, d ≃ 2.93 v w /β. This may explain the similarities. Moreover, if we assume that R ′ p is the relevant length in the exponential nucleation while d is the relevant length in the simultaneous nucleation, their numerical difference is in agreement with the difference in the peaks of the corresponding GW spectra. Notice also that even assuming a direct relationf = 1/R ′ p gives (for v w ≃ 1) f = 0.29 β, which roughly agrees with the result of the envelope approximation, Eq. (3) above.
Regarding GWs generated by sound waves, the actual free parameter in the simulations of Ref. [17] is the average separation d, so this is in principle the relevant quantity. On the other hand, the peak frequency in Eq. (4) suggests a relevant length scale ≃ 0.3d. In this case, though, the width of the sound shell plays a relevant role [66]. In any case, the results of the simulations of Ref. [17] must be interpreted with care, since relatively weak phase transitions were considered [corresponding to the potential (7) with A = 0], most of them with deflagration walls, and it is not clear whether these results extrapolate directly to stronger phase transitions. Beyond that, there is no obstacle in formally extrapolating the approximations of Eq. (4), since d is a well-defined physical quantity. In contrast, the approximations (3) for the collision mechanism depend on the parameter β which becomes negative for strong phase transitions.
To illustrate this problem, we will apply Eqs. (3) and (4) to the model (7). We only consider the runaway regime. For a runaway wall the quantities ρ wall and ρ fl appearing in Eqs. (3) and (4) are given by [48] ρ wall = F U R net , ρ fl ≃ 4(ρ R + 3∆ε − 3F U R net )(0.15v 2 − − 0.132v 3 − + 0.065v 4 − ), where the quantities F U R net , ρ R , ∆ε and v − are defined in Sec. 3. We evaluated these quantities, as well as d and β, at t = t F . The result is shown in Fig. 11. The solid lines correspond to the GWs generated by bubble collisions and the dashed lines correspond to the result from sound waves. The latter have in general a higher intensity, except for very strong phase transitions. This is expected, since for stronger supercooling more energy goes to the bubble walls and less energy goes to the fluid. This is why the GW intensity from sound waves suddenly falls near the value A/v = 0.129. For this model, this is approximately the limit beyond which the phase transition is stuck in the false vacuum while the temperature of the plasma decreases. On the other hand, it is clear that the abrupt growth inΩ coll and the abrupt decrease inf coll are due to the vanishing of β for A/v ≃ 0.127. Beyond this value, β becomes negative and we cannot apply the approximation (3).
In order to extrapolate the results of Ref. [14], where the exponential-rate approximation was used, to very strong phase transitions, we may first write Eqs. (3) in terms of a physical length. Assuming that the actual relevant scale in the collision mechanism is given by the peak of the volume-weighted bubble distribution, as discussed above, we set instanton action S(T ) which determines the nucleation rate is a dimensionless quantity which does not depend on the scale of the theory but only on the shape of the effective potential. Using a simple yet physical model, we have considered model parameters for which S(T ) has a minimum, in which case the phase transition can be very strong.
As long as the minimum is not reached, the main differences with weaker phase transitions are only quantitative. On the other hand, if S gets too close to its minimum, the dynamics becomes qualitatively different. Besides, in this case the usual approximation of linearizing S around a time t * near the completion of the transition, which leads to an exponentially growing nucleation rate, breaks down. In particular, the parameter β = d log Γ/dt, which characterizes the duration of the phase transition and the typical bubble size, will vanish or become negative, which poses a problem for applying results obtained with this approximation to physical models. Nevertheless, around its minimum, S is approximately quadratic, and the nucleation rate can be approximated by a gaussian function of time. In this case, the relevant dimensional parameter α is given by α 2 = (d 2 Γ/dt 2 )/(2Γ). We have solved analytically the development of the phase transition for the gaussian approximation, both for a static universe and for an exponentially growing scale factor. Ignoring the scale factor is justified if the duration of the transition is short enough, which in this case is determined by the time scale α −1 . For cases in which the duration is longer than this time, an exponential scale factor is generally a good approximation. In the general case, the time and length scales are determined by the parameter α as well as by the value of the nucleation rate at the minimum, Γ m . Thus, for instance, the average bubble size at the end of the transition isR ∼ (α/Γ m ) 1/3 . The gaussian approximation also allowed us to determine analytically the limit in which the system remains stuck in the metastable phase, which gives the bound Γ m /H 4 α/H.
We have compared the analytic results for the gaussian and for the exponential rates with those obtained with the fully numerical computation. As we have seen, the exponential-rate approximation is valid for β ≫ α. For β α the gaussian approximation can be used, and in the limit between the two regimes, namely, β = a few α, both approximations are quite good.
For physical models, the parameter region for which the phase transition occurs around the minimum of S is relatively small. In our example model, in which we varied a single parameter, only a 3% of its range corresponds to phase transitions of this kind. In any case, such strong phase transitions are possible and may have important consequences, particularly for the generation of gravitational waves. As we have shown, the dynamics of a very strong phase transition cannot be obtained by direct extrapolation from that of weaker phase transitions. In order to consider this possibility in a simulation (and take into account the dynamics in a realistic way), a gaussian nucleation rate should be considered instead of the usual exponential approximation. In the case of gravitational waves, a possible extrapolation for results which are functions of β is to write β in terms of a physical length such as the peak of the volume-weighted bubble distribution. For a very strong phase transition, this length can be obtained analytically from the gaussian approximation. In this case, the extrapolation consists in the replacement β → 4.1(Γ m /α) 1/3 .