Bubble nucleation and growth in slow cosmological phase transitions

We study the dynamics of cosmological phase transitions in the case of small velocities of bubble walls, $v_w<0.1$. We discuss the conditions in which this scenario arises in a physical model, and we compute the development of the phase transition. We consider different kinds of approximations and refinements for relevant aspects of the dynamics, such as the dependence of the wall velocity on hydrodynamics, the distribution of the latent heat, and the variation of the nucleation rate. Although in this case the common simplifications of a constant wall velocity and an exponential nucleation rate break down due to reheating, we show that a delta-function rate and a velocity which depends linearly on the temperature give a good description of the dynamics and allow to solve the evolution analytically. We also consider a Gaussian nucleation rate, which gives a more precise result for the bubble size distribution. We discuss the implications for the computation of cosmic remnants.


Introduction
A first-order phase transition in the early universe causes disturbances in the plasma which may result in the production of cosmic relics such as topological defects [1], gravitational waves [2], the baryon asymmetry of the universe [3], or baryon inhomogeneities [4,5]. In a first-order phase transition, the system is initially supercooled in a phase which is metastable below the critical temperature T = T c . The transition to the stable phase proceeds through the nucleation and expansion of bubbles, and the properties of the generated cosmic relics depend strongly on this dynamics. The relevant quantities are the bubble nucleation rate per unit volume per unit time, Γ, and the velocity of bubble walls, v w .
The nucleation rate vanishes for T ≥ T c and, below T c , it grows very rapidly as the temperature decreases. Bubble nucleation becomes appreciable when Γ becomes comparable to the Hubble rate H =ȧ/a [6], where a is the scale factor and a dot means a derivative with respect to time. Since Γ is extremely sensitive to the difference T c − T , the most reliable approach is to compute it numerically. However, in order to simplify the treatment of the dynamics, approximations are often used. Since Γ has the well-known form Γ = Ae −S , where S (the instanton action) is a rapidly varying function of the temperature, the most common approximation is to assume a constant factor A and linearize S(T (t)) around a certain time t * This gives an exponential rate Γ(t) = Γ * e β * (t−t * ) .
The bubble wall velocity also vanishes at T = T c , since the pressure is the same in the two phases. At lower temperatures, the pressure is higher in the stable phase, and the walls move. Their velocity v w depends on the friction with the plasma, and is also affected by non-trivial hydrodynamics. Besides, the interactions between bubbles must be taken into account in the global dynamics of the phase transition. To simplify the treatment, a common approximation is to assume that v w remains constant during the phase transition. Since the pressure difference goes roughly with the difference T c − T , this will be a reasonable approximation as long as the temperature has a relatively small variation, δT ≪ T c − T . In principle, the temperature decreases with time due to the adiabatic expansion of the universe, and this condition translates into a condition for the time intervals, δt ≪ t − t c . For an exponential nucleation rate, this requires β −1 * ≪ H −1 , which is satisfied in general.
In some calculations, more drastic simplifications are needed, such as considering a constant nucleation rate or assuming that all bubbles nucleate at the same time. These approximations are quite different from the exponential growth, and are in general regarded as rough approximations. However, the exponential approximation is not always valid either. For instance, for very strong phase transitions, the function S(T ) may have a minimum at a certain temperature T m , and the supercooling may be such that this temperature is reached. In such a case S(T ) cannot be linearized, and a Gaussian approximation for Γ(t) is more suitable [7] (besides, the phase transition becomes slow [8]). Nevertheless, in most cases the temperature is not close to the minimum of S(T ), and the latter may be assumed to be a monotonically increasing function of T . For decreasing T we have, around a given time t * , an exponentially growing rate.
However, these arguments do not take into account the fact that latent heat is released as bubbles expand. This energy reheats the plasma, causing the temperature to reapproach T c . As a consequence, both Γ and v w may decrease, slowing down the phase transition. The latent heat is released at the bubble walls, and the reheating is in general inhomogeneous. This fact makes the general treatment of the phase transition difficult, except in some special cases. One of them is the case of detonation bubbles [9], in which the velocity of the walls is so high that the fluid in front of them remains unperturbed. In this case, the reheating occurs only inside the bubbles, so in the old phase the nucleation rate grows according to the adiabatic cooling. Thus, for detonations, an exponential nucleation rate is generally a good approximation, and the phase transition is quick enough to assume a constant wall velocity.
In contrast, for smaller velocities, the wall propagates as a deflagration front [10] and is preceded by a shock wave which carries away the released energy. In this case, the temperature outside the bubbles may be very inhomogeneous, and the nucleation rate may vary significantly from region to region. Nevertheless, the shock waves propagate supersonically and the latent heat quickly distributes throughout space. For small enough wall velocities, we may assume a homogeneous temperature as bubbles nucleate and expand [5]. Once the rate at which latent heat is released exceeds the energy decrease due to the adiabatic expansion, the temperature will stop decreasing and start to grow. This dynamics has been discussed, to different extents, in Refs. [5,11,12,13,14,15]. The treatment in these works is either numerical or involves rough approximations. This is because, even assuming a homogeneous temperature, its variation is related to that of the wall velocity and the nucleation rate through non-trivial integro-differential equations. In contrast, in the detonation case, the assumptions of a constant wall velocity and an exponential nucleation rate simplify considerably the computations. Depending on other approximations, it is even possible to obtain the complete development of the phase transition analytically. The results are functions of a few free parameters, such as v w and β * , which can be computed numerically for specific models, and one obtains a very precise description of the phase transition (see [7] for a recent discussion).
Notice that, in the slow deflagration case, the temperature will generally have a minimum T m at a time t m separating the supercooling and reheating stages. Around this time, the exponent S(T (t)) is quadratic in t − t m , and Γ(t) should be well approximated by a Gaussian function. Moreover, it was pointed out in Refs. [11,12] that the time during which Γ is effectively active can be quite smaller than the total duration of the phase transition. In such a case, even a delta function may be a suitable approximation for Γ(t). On the other hand, slow deflagrations will occur in general in phase transitions with little supercooling. In such a case, it is a good approximation to consider v w to first order in T c − T . In this paper we shall investigate this scenario. We shall discuss the validity of these approximations for the nucleation rate and the wall velocity, and we will show that they simplify considerably the numerical treatment, even allowing to obtain analytic solutions. These approximations depend on a few parameters, which, for a given model, can be estimated numerically with simple computations.
The plan is the following. In the next section we study the general dynamics of a phase transition mediated by slow deflagrations, and we discuss the possible approximations for the distribution of the latent heat and the dependence of the wall velocity on the temperature. In Sec. 3 we consider several approximations for the nucleation rate. In particular, we show that assuming an exponential rate until the time t m gives a good estimation of the temperature T m . However, assuming a Gaussian nucleation rate gives a better approximation for the final bubble size distribution. On the other hand, assuming a simultaneous nucleation at t = t m is a good approximation for computing the evolution of T and v w . We study this evolution for a delta-function rate in Sec. 4. In Sec. 5 we discuss on cosmological applications of our treatment, and in Sec. 6 we summarize our conclusions.
2 Development of a slow phase transition

Effective potential and model parameters
To describe a first-order phase transition, we shall consider a model consisting of a scalar field φ with a spontaneous symmetry-breaking tree-level potential of the form V (φ) = −(m 2 /2)φ 2 + (λ 0 /4)φ 4 , and particles which acquire their masses m i from the vacuum expectation value the field φ. For m i /T 1, the one-loop finite-temperature effective potential V (φ, T ) has an expansion in powers of φ/T (see, e.g., [16]), The coefficients in Eq. (1) depend on the parameters m and λ 0 as well as on the particle content. These coefficients are constant, except for λ, which has a logarithmic dependence on T . Since the temperature will not depart significantly from the critical temperature, we shall neglect this dependence. The complete free energy density is given by where the constant ρ V is the false-vacuum energy density and the second term is the radiation component. This term is proportional to the total number of degrees of freedom of the species in the plasma, g * , while only a part ∆g of these degrees of freedom have strong enough couplings to φ and contribute to the term V (φ, T ). At high enough temperatures, V (φ, T ) has an absolute minimum at φ = 0 ≡ φ + , while at low enough temperatures the absolute minimum is at There is a temperature range in which these two minima coexist. The lower temperature in this range is T = T 0 , below which φ = 0 becomes a maximum. The critical temperature, defined by the equation F (φ + , T c ) = F (φ − , T c ), is given by In this model, the dimensional parameter T 0 determines the temperature scale of the phase transition. This parameter can be related to the zero-temperature minimum v, which is usually considered as the energy scale of the theory, 2DT 2 0 = λv 2 . The two phases are characterized by the functions F ± (T ) = F (φ ± , T ), from which we can derive thermodynamic quantities such as the pressure p ± = −F ± , the entropy density s ± = −dF ± /dT , and the energy density ρ ± = F ± − T dF ± /dT . Thus, in the high-temperature phase, we have ρ + = ρ V + ρ R , where ρ R is the radiation energy density, and the energy difference between the two phases is given by ρ + − ρ − = T ∂V /∂T − V . Its value at the critical temperature is the latent heat L. In this model we have where φ c ≡ φ − (T c ) is the jump of φ at the critical temperature, which is given by For E = 0 we have T c = T 0 and φ c = 0, which means that the phase transition is second order. The strength of the phase transition is usually measured by the parameter φ c /T c . For φ c /T c 1 the phase transition is said to be strongly first order, while for φ c /T c ≪ 1 it is said to be weakly first order. Using dimensionless quantities such as T /T c and φ(T )/T , we see that the dynamics will depend mainly on the dimensionless constants D, E, λ which determine the shape of the effective potential. In general, for φ c /T c ∼ 1, most of the relevant quantities, such as L/T 4 c , will be of order 1. However, the dynamics of reheating will depend on the ratio L/ρ Rc , where ρ Rc = ρ R (T c ) (a value L ≪ ρ Rc will not cause a significant temperature change in the plasma). Since only a few of the degrees of freedom contribute to the free energy difference V (φ, T ), we expect in general L ∼ T 4 c , while we have ρ Rc ∼ g * T 4 c . Hence, we will have, typically, L/ρ Rc ∼ 1/g * . On the other hand, the dynamics will have a dependency on the temperature scale through the Friedmann equation where M P = 1.22 × 10 19 GeV is the Plank mass. The false-vacuum energy density in this model is given by ρ V ≃ λv 4 /4, and for T ≃ T c ∼ v we have ρ V /ρ R ∼ λ/g * . For many physical models there is a large number of degrees of freedom in the plasma, and we have For our general treatment we may regard the constants v, g * , D, E, and λ as free parameters. For specific computations we shall consider electroweak-scale values for v and g * , namely, v = 250GeV, g * = 100, and we shall set the value of λ by demanding a natural value for the scalar mass m φ = √ 2λ 0 v 2 . Assuming λ 0 ≃ λ, we choose λ = 0.125, corresponding to m φ ≃ v/2. On the other hand, the parameter E is generally smaller, since it depends cubically on the couplings of φ with gauge fields. Hence, according to Eq. (7), this model does not naturally give very strong phase transitions 1 . Nevertheless, for E ≥ λ/2 we have φ c /T c ≥ 1. We shall consider the value E = 0.075, corresponding to φ c /T c = 1.2. The coefficient D, in contrast, is quadratic in the couplings and involves a sum over all particle species, and we may have D ∼ 1. Typically, we have E < λ < D, and Eq. (4) gives T c − T 0 ≪ T 0 . Notice that, according to Eq. (6), for φ c /T c ≃ 1 and E ≪ D we have L/T 4 c ≃ 2D ∼ 1, as expected. For these natural values, we expect a reheating ∆T /T ∼ L/ρ R ∼ 10 −2 . This competes with the amount of supercooling, since the phase transition will take place in the range T 0 < T < T c , and we have (T c − T 0 )/T c ≃ E 2 /(2λD) ∼ 10 −2 . We shall consider a couple of values of the ratio L/ρ Rc ; namely, L/ρ Rc = 0.025 (corresponding to L ≃ 0.82T 4 c and D ≃ 0.33), and L/ρ Rc = 0.05 (corresponding to L ≃ 1.64T 4 c and D ≃ 0.62).

Initial nucleation and growth
The temperature variation is governed by the adiabatic-expansion equation ds/dt = −3Hs. In our model, for the high-temperature phase we have s + ∝ T 3 . Hence, be-fore the phase transition the temperature decreases with a rate Between T = T c and T = T 0 the effective potential has a barrier between the minima, and the phase transition may proceed by bubble nucleation. For T ≤ T 0 the barrier separating the minima disappears, and the phase transition will proceed by spinodal decomposition.
In the bubble-nucleation range, there is a probability of nucleating a bubble per unit volume per unit time given by [17,18] where S(T ) = S 3 (T )/T , A(T ) = T 4 [S 3 (T )/(2πT )] 3/2 , and S 3 is the spherically symmetric extremum of the three-dimensional instanton action [18]. This extremum also gives the configuration of the nucleated bubble. We will solve numerically the equation for the instanton configuration by the undershoot-overshoot method (see [7] for details).
It is well known that the action S 3 diverges at T = T c and vanishes at T = T 0 . Hence, the nucleation rate vanishes at the critical temperature and reaches values Γ ∼ T 4 as T approaches the value T 0 . The latter is, relatively, a huge rate, since the phase transition will occur, roughly, for Γ ∼ H 4 , and we have T 4 /H 4 ∼ M 4 P /T 4 (which is generally large unless the scale v is very close to the Planck scale). Hence, the phase transition will generally complete before reaching the spinodal decomposition temperature. We shall take the time t H defined by the equality Γ = H 4 as a reference time. The corresponding temperature T H is given by the equation S − 3/2 log(S/2π) = 4 log(T /H).
The time t N at which bubble nucleation effectively begins is usually defined by the condition that there is one bubble in a Hubble volume, and generally we have t N > t H (see [7] for a recent discussion). During this supercooling stage, the number density of bubbles is given by with the time-temperature relation given by Eqs. (9) and (8), with ρ = ρ + . This expression does not take into account the fraction of volume occupied by bubbles, which at this time is negligible, and we do not include the dilution of the number density as a −3 either, since we have ∆a/a ∼ ∆T /T ≪ 1. Thus, the time t N is determined from the condition nH −3 = 1.
A nucleated bubble grows due to the pressure difference between the two phases. The motion of a bubble wall causes bulk fluid motions as well as temperature gradients in the plasma, which affect the propagation of the phase transition front. It is well known that the treatment of hydrodynamics is quite simple for the bag equation of state (EOS) [9]. In our case, the free energy (2) has exactly the bag form in the phase +, where ε + = ρ V and a + = π 2 g * /30. On the other hand, since we have little supercooling, for the phase − we may use the approximation which also has the bag form. Taking into account the relations F − (T c ) = F + (T c ) and L = 4T 4 c (dF − /dT 4 − dF + /dT 4 )| Tc , we may write with ε − = ε + − L/4 and a − = a + − 3L/(4T 4 c ). With this approximation, several results will depend only on the variable This quantity is larger for stronger phase transitions and smaller for weaker ones 2 . The force which drives the motion of a bubble wall is not simply given by the pressure difference p − (T ) − p + (T ) since, in the first place, the temperature varies across the wall. Using a linear approximation for the temperature variation inside the wall, the driving force can be written as [19] where T + and T − are the values of T in front and behind the wall, respectively. The relation between T + and T − is obtained by energy-momentum conservation [20]. This relation involves also the values of the fluid velocity on each side of the wall. In the wall reference frame, we denote the magnitude of the incoming fluid velocity by v + and that of the outgoing fluid velocity by v − . For non-relativistic velocities we have and the bag EOS gives Notice, also, that we have the relation a − /a + = 1 − 3α(T c ). There is also a friction force, which arises as a consequence of the departures of the particles distributions from their equilibrium values inside the wall (see, e.g., [22]). A phenomenological approach to this force is often used, and will be sufficient for our purposes. In the wall reference frame, the friction is modeled by a function of the average fluid velocityv, and the latter is approximated byv = (v − + v + )/2 [23]. In the non-relativistic case we have where the friction coefficient η is a free parameter which may be inferred by matching to the full microphysics computation. The wall velocity is obtained from the steady state condition F dr + F fr = 0. To solve this equation, we need additional relations between the fluid velocities v ± and v w . These relations are obtained from the fluid equations and the boundary conditions. For an isolated bubble, the fluid is at rest far in front of the wall as well as far behind it (at the bubble center). For the small wall velocities we are interested in, which are subsonic with respect to the bubble center, it turns out that we have a vanishing fluid velocity inside the bubble, i.e., v − = v w . Hence, Eqs. (16)(17)(18)(19) give where α + = α(T + ). Inside the bubble we also have a constant temperature.
Since the wall is subsonic with respect to the fluid in front of it, this fluid is affected by the wall motion. The fluid equations relate the variables next to the wall, T + , v + , to the boundary conditions far in front. It turns out that the fluid profile ends in a discontinuity, a shock front which propagates supersonically in the phase +. The temperature falls from the value T + in front of the wall to a value T sh at the shock front. Beyond this discontinuity we have an unperturbed fluid at temperature T . The matching conditions here are similar to those for the phase transition front. For small fluid velocities, the profile of the shock wave can be solved analytically. In this case, the temperature variation is small, (T + − T sh )/T sh ∼ v 2 w . Besides, the velocity of the shock front is v sh ≃ c s , and the shock discontinuity becomes very weak, T sh ≃ T (the error of this approximations is of order e −1/v 3 w [28]). Hence, we have We shall use the complete expression (20) for the numerical computations, but we will see that Eq. (21) is a good approximation in our case. The friction coefficient η is not directly related to the parameters of the effective potential, and we shall regard it as an independent free parameter. According to Eq. (21), for dimensionally natural values η ∼ T 4 c ∼ L we will have v w ∼ (T c − T 0 )/T c ∼ 10 −2 for our model parameters. For specific computations we will use the value η/T 4 c = 0.5. For η/T 4 c ≪ 1, the non-relativistic approximations in Eqs. (17)(18)(19) will no longer be valid. In the general case, Eqs. (18) become where γ ± = 1/ 1 − v 2 ± . Several extrapolations of the friction force (19) to the relativistic case have been proposed. Some of them [24,25] take into account the fact that the friction may saturate in the ultra-relativistic limit and the wall may run away [26]. However, for a potential of the form (1) the bubble wall cannot run away, and a reasonable extrapolation for the friction force is given by where γv = (v + γ + + v − γ − )/2 (for a recent discussion, see [7]).
We have different kinds of hydrodynamic solutions, corresponding to different branches of the v + -v − relation. If the velocity of the incoming flow is lower than the speed of sound in the plasma, c s = 1/ √ 3, we have a deflagration, while if the incoming flow is supersonic, we have a detonation (see e.g. [21] for details). Although we are interested in the case of very slow deflagrations, for comparison we shall also consider the detonation case. A detonation wall moves supersonically with respect to the fluid in front of it. As a consequence, this fluid is not affected by the wall, and we have v + = v w and T + = T , where T is the value of the temperature in the absence of any bubbles. Using these conditions in Eqs. (16), (22)(23) (with the + sign), and (24), we may solve for v w as a function of η and α (for more details, see, e.g., [27]). If α is too small or the friction is too large, it turns out that there is no detonation solution. In this case we have to consider a deflagration, which is compatible with smaller velocities. To obtain a detonation, we shall consider the value η/T 4 c = 10 −3 .

Global dynamics
For isolated bubbles, the temperature T is determined by the adiabatic expansion, where a c ≡ a(T c ). We may assume that, as temperature varies, the walls instantly take the terminal velocity v w (T ) (see, e.g., [5]). Once bubbles begin to meet each other, the growth of a cluster of bubbles will depend on the motion of the uncollided walls. In the case of detonations, the conditions in front of a wall are always the same as for an isolated bubble (namely, v + = v w , T + = T ), and we expect the results for the wall velocity to remain essentially unchanged. Also, the nucleation dynamics in the phase + is not affected by the presence of previously nucleated bubbles. Therefore, v w and Γ are functions of the outside temperature T , which is given by Eq. (25). In contrast, for deflagrations, the fluid in which a bubble expands is affected by shock fronts coming from other bubbles. The main effect is an increase of the temperature, which will decrease the wall velocity. Besides, the reheating will diminish the nucleation of new bubbles. Initially, the isolated deflagration bubbles are contained inside larger spheres whose surfaces are the shock fronts. Each of these "shock bubbles" has a radius R sh ≃ (c s /v w )R, where R is the radius of the bubble. To characterize the moment after which bubbles cannot be regarded as isolated, let us consider the time t 0 at which, in average, neighboring shock bubbles have just met each other. We may estimate this time by the condition that the average shock radius matches the average bubble separation d (which depends on the bubble number density). At this time, the fraction of volume which is in the new phase is For slow walls, we will have f − ≪ 1, which means that most of the phase transition will occur with interacting bubbles. Nevertheless, in this case we have extremely weak shocks which softly change the boundary conditions for the wall motion. Hence, we may assume that Eq. (20) still applies locally, with an inhomogeneous temperature T + .
Once the shock waves of a bubble have reached several other bubbles and bounced several times between two neighboring bubbles 3 , we may assume that the released energy is homogeneously distributed. It will take the shocks a few times t 0 to complete this homogenization. For v w ≪ c s this will happen when the fraction of volume in the new phase is still very small. Hence, during this homogenization process the reheating will still be insignificant, and we may just assume a homogeneous distribution of the latent heat from the beginning (we have checked these features with our numerical computations). This approximately instant spreading of the released energy will continue during the whole phase transition, due to the velocity difference between walls and shocks.
Notice that the temperature cannot be homogeneous everywhere, since we have discontinuities at the phase transition fronts, given by Eqs. (17)(18). For the bag EOS, the matching condition , the energy density has the same variation on both sides of the walls. This is consistent with the assumption of a homogeneous distribution of the released energy, and we may assume homogeneous temperatures T ± in each phase. To study the temperature variation, the most straightforward way is to consider the conservation of entropy, where s c = s + (T c ), and f ± is the fraction of volume occupied by each phase. Since the system is out of equilibrium, the entropy is not exactly conserved. In the appendix we estimate the error of this approximation, which for the present case is negligible. Like in the detonation case, we shall denote the homogeneous temperature outside the bubbles by and T − (T ) is given by Eq. (18). The wall velocity is still given by Eq. (20), with T + = T and α + = α. Thus, for very slow deflagrations the wall velocity and the nucleation rate depend on a single variable T , like in the detonation case. The main qualitative difference in the dynamics of these two cases is the reheating term proportional to f − in Eq. (27), which is absent in Eq. (25). For T ≃ T c we have s + − s − ≃ L/T c and s c ∝ ρ Rc , so this term is proportional to the ratio L/ρ Rc . Thus, Eq. (27) takes into account the released energy as well as the heat capacity of the plasma. Indeed, ignoring the last term, this equation gives a temperature increase ∆T ≃ (dT /dρ R )Lf − , as expected.
To compute the fraction of volume (either in the detonation or the deflagration case), we shall assume that the new phase is composed of spherical bubbles which may overlap. At time t, the radius of a bubble which nucleated at time t ′ is given by where we neglected the initial radius. This is generally valid unless T is very close to the Planck scale. The fraction of volume in the old phase is given by [29,30] f where Here, t c is the time at which T = T c . In these equations we have ignored the effect of the scale factor on physical lengths (namely, the stretching of R and the dilution of the density of nucleated bubbles). We shall take into account this effect in the numerical computations, but it is negligible for the cases with little supercooling we consider. Indeed, although for deflagrations the phase transition may last quite longer than for detonations, we will have a duration of order 10 −2 H −1 and, thus, a/a c ≃ 1. Notice, on the other hand, that the variation of a cannot be ignored in Eqs. (25) or (27), since a small change in T will cause a large change in Γ. A measure of progress of the phase transition is given by the fraction of volume f + (t). As in our work [7], we shall define a few reference points in this evolution. The first of them corresponds to the "initial" moment t I at which f + (t I ) = 0.99. The second one is the percolation time t P , which is approximately given by f + (t P ) = 0.71. Another reference time, which is often considered, is the time t E at which the fraction of volume has fallen to f + (t E ) = 1/e. Finally, we define the "final" time t F by f + (t F ) = 0.01.

Numerical results
As discussed in previous subsections, we shall fix the potential parameters to typical values, we shall consider η ∼ T 4 c , which gives slow deflagrations, and we shall consider a couple of representative values for the ratio L/ρ Rc . In order to compare with the detonation case, we shall also consider the case η ≪ T 4 c . In the left panel of Fig. 1 we show the evolution of the temperature. The time is normalized to the Hubble time H −1 c , where H c ≡ H(T c ). The solid curve corresponds to the case L/ρ Rc = 0.025. The initial part of the curve corresponds both to deflagrations (η/T 4 c = 0.5) and detonations (η/T 4 c = 0.001), since the supercooling stage does not depend on the friction. The evolution of the temperature is initially determined by Eq. (25), and the curves would separate once the reheating becomes appreciable in the deflagration case. However, for the detonation the phase transition completes sooner due to the higher wall velocity.
The reference times t H , t N , t I , t P , t E t F are indicated on the curve by dots for the deflagration and by ticks for the detonation. We see that the first two dots (t H and t N ) coincide with the first two ticks. In contrast, the rightmost (purple) tick, indicating the time t F for the detonation case, is to the left of the blue dot, which indicates the time t I for the deflagration case. This means that, for detonations, a 99% of space is in the new phase even before a 1% is reached for deflagrations. In the deflagration case, we observe that, as soon as the fraction of volume occupied by bubbles begins to be noticeable (i.e., at t ≃ t I ) the reheating becomes noticeable too. Once bubbles have filled most of space (t ≃ t F , purple dot), the release of latent heat ceases and the temperature decreases again. The solid line in the right panel corresponds to the wall velocity. We see that the velocity decreases during reheating 4 , as expected from Eq. (21).
The dashed curves correspond to the case L/ρ Rc = 0.05. Notice that in this case we have less supercooling that in the previous one. This is because larger L roughly corresponds to a larger parameter D, and this gives a smaller value of (T c − T 0 )/T c [see Eqs. (6) and (4)], which is reflected in a smaller value of (T c − T N )/T c . In spite of the smaller supercooling, the maximum velocity is similar to that of the previous case (v w ∼ 10 −2 ). This is because the pressure difference is roughly proportional to L, which is reflected in v w , as can be seen in the approximation (21). On the other hand, in this case there are no detonation solutions, no matter how small the friction coefficient.
Like in the previous case, the temperature begins to grow at t ≃ t I and decreases again for t ≃ t F . In this case, though, this time interval is longer, and we have a stage of approximately constant temperature T r very close to T c . This occurs due to the larger latent heat, which would actually be enough to reheat the system to a temperature T > T c . Nevertheless, the backreaction of reheating on bubble growth prevents this to happen. In our case, the released energy gets quickly distributed and the increase in energy density is given by Lf − , which is initially small. As the temperature gets close to T c , the wall velocity decreases significantly, which can be appreciated in the right panel of Fig. 1. The phase transition slows down, preventing further release of latent heat. The reheating temperature T r , as well as the wall velocity during this stage, are determined by the balance between the rate at which energy is injected, which is roughly given by Ldf − /dt, and that at which the adiabatic expansion takes energy from the plasma, which is roughly given by 4ρ R H. This gives the equation (L/ρ R )df /dt = H. We discuss this approximation further in Sec. 4.
These results are in qualitative agreement with previous works (see, e.g., [5,12,14,15]). The effects of a significant velocity slow-down have been investigated in some detail in Refs. [5,12,14]. In Refs. [4,13,31,32], the limit of a long phase-coexistence stage at T ≃ T c has been investigated. Here, we shall consider the general case, and we shall focus on the dynamics of nucleation. Some quantities of interest are the average nucleation rate Γ (t) = f + (t)Γ(t), the number density of bubbles, the average distance between centers of nucleation, and the distribution of bubble sizes, where t R is the time at which the bubble of radius R was nucleated, which is obtained by inverting Eq. (28) for t ′ as a function of R.
In Fig. 2 we show the evolution of some of these quantities for the detonation and the two deflagration cases. In the detonation case (upper plots), the development of the phase transition is determined by the extremely quick growth of the nucleation rate (see the left panel). The main features are the following. Once bubble nucleation effectively begins, the phase transition completes in a relatively short time, i.e., t F − t N ≪ (t N − t c ). Moreover, the variation of the fraction of volume occurs in an even shorter time t F − t I . Most bubbles nucleate in this short interval near t F . In the right panel we see that the average bubble size grows very slowly during most of the phase transition. This is due to the constant nucleation of very small bubbles. Only when the volume in the old phase becomes small and the nucleation of bubbles turns off, the average bubble radius begins to grow with velocity v w . On the other hand, the average separation between centers of nucleation inherits initially the rapid variation of Γ, and becomes constant whenΓ vanishes. At t ≃ t F we have d ∼R, as expected.
The second row of Fig. 2 corresponds to the deflagration case L/ρ Rc = 0.025. The model parameters are exactly the same as in the previous case, except for the friction. We see that the evolution of the various quantities is quite different. In particular, since Γ is very sensitive to the temperature, the nucleation rate turns off as soon as the reheating begins. At this moment we still have f + ≃ 1, soΓ(t) almost coincides with Γ(t). Since bubble nucleation ceases so soon, the average distance d reaches its final value at t ≃ t I (in contrast, in the detonation case this happens at t ≃ t F ). Similarly, the transition from an approximately constant average radius to the behavior dR/dt = v w occurs at t ≃ t I and not at t ≃ t F . Notice that in the deflagration case the final bubble separation is smaller, so the final number of bubbles is higher than in the detonation case. This is because, in this slower phase transition, the time at which the nucleation stops is later than for the detonation case (see the left panels). In this small time difference, the nucleation rate reaches quite higher values.
The effect of reheating is more marked in the third row of plots, corresponding to the deflagration case with L/ρ Rc = 0.05. Since in this case the interval t F − t N is longer (due to the slow-down after reheating) the nucleation rate looks like a sharp peak at t ≃ t I . Like in the previous case, at t ≃ t I the distance d becomes constant and the average   radiusR begins to grow with velocity v w . The subsequent change of slope ofR and f + corresponds to the decrease of v w . In Fig. 3 we plot the distribution of bubble sizes for the three cases at the times t I , t P , t E and t F . We also indicate, for comparison, the size (at each of these times) of a bubble which nucleated at t = t N . In the fast wall case, we see that the size distribution maximizes at R = 0 until the phase transition is quite advanced (more precisely, until f + = e −1 ). This is again due to the rapidly increasing nucleation of vanishingly small bubbles until the time t E . In contrast, in the slow wall case, for t t I no new bubbles are nucleated, and the size distribution just shifts to larger radius. Hence, the maximum separates from R = 0. For the same reason, the ratio of the average radiusR(t) to that of the largest bubbles, R(t N , t) is smaller in the detonation case than in the deflagration case.
For some applications, such as gravitational-wave generation, it is more appropriate to consider the volume-weighted average radius. However, in the case of slow deflagrations, there will not be a large difference between the weighted and unweighted averages, since all bubbles nucleate around the same time t ≃ t I and, hence, have similar sizes. To illustrate this, in Fig. 4 we show the radius distribution together with the volume-weighted distribution at the time t = t E . We see that only in the detonation case the difference may be relevant.

Approximations for the nucleation rate
In the computations of the previous section we have made several approximations in order to simplify the treatment. However, much simpler approximations are often used, such as assuming a constant wall velocity and even a constant nucleation rate. As we have seen, none of these is a good approximation in our case. For the wall velocity we may use the relatively simple form (21). On the other hand, for the nucleation rate, it is not easy to find a simple approximation. For T close to T c , the thin wall approximation can be used for the instanton action, and we have S(T ) ∝ 1/V (φ − , T ) 2 [16]. Linearizing V as in Eq. (13), we obtain a nucleation rate of the form where the constants A and B depend on the potential parameters D, E, λ (see, e.g., [11]). This expression shows that Γ is a very rapidly varying function of T , and assuming Γ = constant will generally be a bad approximation. From Eq. (34) we may obtain Γ(t) using the appropriate time-temperature relation, such as Eq. (25) or Eq. (27). However, analytic approximations for the nucleation rate may introduce large errors. Therefore, it is usual to consider a semi-analytic approach, which consists in linearizing the exponent S(t) in Eq. (10) around a certain time t * . This procedure only requires to compute numerically S and its derivative at T * = T (t * ).

Exponential rate
As pointed out in Ref. [7], the linearization of S(t) actually involves linearizing both functions S(T ) and T (t), and any of these approximations may fail. Nevertheless, except in special cases of very strong phase transitions, the temperature T * will not be close to a minimum of S(T ), and we may expand S to first order, provided that the temperature variation is small enough. Assuming that (in the radiationdominated era) the temperature decreases with time as dT /dt = −HT , then for a short enough time we may write with H * = H(T * ). Hence, we have an exponentially growing nucleation rate where Γ * = Γ(T * ) and β * = (HT dS/dT )| T * . This rate is simpler than Eq. (34), and may be a much better approximation than the latter if t * is conveniently chosen, so that the values of interest of t are close to it 5 . The exponential rate is a good approximation for a phase transition mediated by detonations (see [7] for a recent discussion). On the other hand, as can be appreciated in Fig. 2, this is not a good approximation for deflagrations. This is because Eq. (36) will no longer be valid in the presence of reheating.
Nevertheless, this approximation will always be valid in the initial stages of the phase transition. For the initial number density (11), the nucleation rate (37) gives n(t) ≃ Γ(t)/β * , and the condition for the onset of nucleation is The quantity β * is quite sensitive to the temperature, which implies that the exponential rate is a good approximation only in a small interval around t * . Nevertheless, since Γ grows quickly with time, the interval of interest is generally small. For instance, in Eq. (11) it is convenient to choose t * = t N , since Γ decreases rapidly for smaller times. Thus, Eq. (38) gives the equation (in the last term we have used the approximation T N ≃ T c ). This equation gives a very good estimate for the temperature T N [7]. For the detonation case, we may use the exponential rate beyond t = t N . For a small time interval we may also assume a constant wall velocity, in which case Eqs. (28)(29)(30) can be integrated analytically [6]. We have From Eq. (40) we may obtain analytic expressions for quantities such as n or dn/dR, as well as analytic estimations of time intervals. These analytic results can be applied to physical models by computing the parameters Γ * , β * and v w (T * ). Moreover, the parameter Γ * is not too relevant, since the dynamics of nucleation depends essentially onΓ/Γ ≃ β * . The parameter β −1 * . Therefore, the time β −1 * characterizing the dynamics will be generally much smaller than t * − t c . This is also in agreement with the top panels of Fig. 2.
In the deflagration case, these results will be valid as long as the reheating is not appreciable. According to Eq. (27), the reheating rate will be roughly proportional to df − /dt. From Eqs. (38)(39)(40)

Sudden reheating
For t > t N , the reheating will be noticeable at some point, and the exponential rate (37) will break down. Since the temperature variations are relatively small, we may still use the linear approximation (35) for S(T ). However, we must replace the linear timetemperature relation (36) with the relation (27), which takes into account the reheating. In order to obtain analytic results, we need to simplify further this equation. We shall accomplish this by assuming a completely homogeneous temperature, T − ≃ T + ≡ T . This is valid since T + − T − ∼ (L/ρ Rc )(T c − T + ) (see the appendix). Using this approximation in Eq. (27), we obtain the simple relation 6 where ∆s(T ) = s + (T ) − s − (T ). If we differentiate Eq. (41), to lowest order in T − T c and L/ρ Rc we obtain (see the appendix for details) where a dot indicates a derivative with respect to time, and we have defined the parameter In terms of bag parameters, we have r = 3α(T c ). As we have seen, for t ≃ t N we haveḟ − ≃ 0. In the approximation (42), at this stage the temperature decreases linearly,Ṫ ≃ −T c H c . This behavior is observed in Fig. 1. In the units of this figure, the slope of the curve is −1. In our numerical examples, the reheating becomes noticeable for t ≃ t I , i.e., for f − ≃ I ∼ 10 −2 . In the general case, this will happen when the reheating term in Eq. (42), rḟ − , becomes comparable to the adiabatic cooling term 3H c . Assuming that this happens for small f − , we have f − ≃ I. Assuming also that the dynamics is still characterized by the time β −1 * , we haveḟ − ∼ β * I. We thus obtain the condition f − ∼ (3/r)(H c /β * ) for the reheating to become noticeable. Notice that, even though r is a small number, H c /β * is even smaller, so the approximation I ≪ 1 is generally consistent. For our example cases, this estimation gives f − ∼ 10 −2 .
More precisely, the equality of the cooling and reheating terms in Eq. (42), gives the condition for the minimum temperature T m . In order to estimate the time t m , we notice thatḟ − (t) has an extremely rapid growth at this stage, so an instant before t m the reheating term in (42) was negligible. Therefore, we have adiabatic cooling until almost t = t m . This can be appreciated in Fig. 1. Then, to solve the condition (44) we may use the exponential-rate approximation. Assuming small I, we haveḟ − ≃İ, and using the result (40) we have 8πv 3 w Γ(t m )/β 3 * = 3H c /r. We may also choose t * ≃ t m , and we obtain where Γ m = Γ(T m ), v m = v w (t m ), and This gives a semi-analytic equation for T m , The wall velocity in the last term can be replaced by the estimation (21) as a function of T m . Moreover, due to the logarithmic dependence, it can be even replaced by the constant v w (T N ) with no significant error. With the latter simplification, the approximation (47) gives the value of T m with a relative error of order 10 −4 with respect to the numerical computations of Sec. 2. Using the result in the linear approximation (36) (with t * = t m ), we obtain t m − t c with a relative error of order 10 −3 . This estimation of t m is so good because, in the first place, the reheating is almost instantaneous, which justifies taking the limit t * → t m . In the second place, the use of Eq. (40) involves a time interval of order β −1 m . We may obtain estimations for longer time intervals using this approximation; however, the errors will be larger. For instance, the time t m − t H is obtained by comparing the result (45) with Γ(t H ) = H 4 . This gives Similarly, comparing (45) These analytic relations assume an exponential rate, as well as a constant wall velocity, for time intervals which are several times β −1 m (e.g., t m − t H ∼ 10β −1 m ). For our examples, the agreement with the numerical computation is around a 25%.
Since f − grows so rapidly at t = t m , we may assume a sudden reheating at this point, such that Γ(t) vanishes for t > t m . This corresponds to a nucleation rate of the form Since we have f + ≃ 1 for t < t m , we can make the replacementΓ(t) ≃ Γ(t) in Eqs. (31)(32)(33), which simplifies the computation of the quantities n,R and dn/dR. Moreover, the final number density of bubbles is given by n f = n(t m ), and we have so the final average distance between centers of nucleation is given by This estimation gives the value of d f to a 6% of the numerical computation. On the other hand, assuming a constant wall velocity for t < t m , the exponential rate gives a constant value for the average bubble radius,R This is in agreement with the initial behavior ofR in the plots on the right of Fig. 2. The small slope in the plots is due to the fact that v w is not constant and the nucleation rate is not exactly exponential. For t > t m , the computation ofR(t) is involved even with these simplifications, and we shall consider a different approximation for Γ(t) below. Nevertheless, notice that the value ofR at t = t F is approximately given by the average bubble separation,R(t F ) ≃ d f . In our examples, the value (52) is a factor ∼ 10 larger than the initial value (53).

Gaussian nucleation rate
The approximation (50) gives a function Γ(t) with a sharp peak at t = t m , while the actual nucleation rate has a differentiable maximum. As we have seen, for the computation of several quantities this qualitative difference is not important. However, it will be reflected, for instance, in the shape of the bubble size distribution. To obtain a smooth function, we may use the quadratic expansion of T (t) around its minimum at t m . Together with the linear approximation (35) for S(T ) around T * = T m , this gives a Gaussian nucleation rate, where This approximation will be valid only for t close enough to t m . Nevertheless, away from t = t m the nucleation rate decreases rapidly and its value is less relevant. The first factor in Eq. (55) is related to the constant β m defined in Eq. (46) 7 . On the other hand, from Eq. (42) we obtain the last factor,T = rT cf− /3. For I ≪ 1 we havë f − =Ï, and we may write We have seen that Eq. (47) provides a good estimation for the values of T m and t m .
To computeÏ(t m ) from Eqs. (28)(29)(30), we must take into account that the velocity has a maximum at t = t m , as can be seen either from Eq. (20) or Eq. (21), and also in Fig. 1. We thus haveÏ The right hand side does not involve any derivatives at t = t m , so we can use again the approximation (50) to estimate the integral. Notice that the latter is given by n(t m )R(t m ) = Γ m v m /β 2 m . Finally, using Eq. (45), we obtain 7 Notice that β m is defined as a function of the temperature T m and does not depend on the dynamics. In particular, β m does not coincide with dS/dt| tm , as it would in the absence of reheating. In the present case the latter derivative vanishes, since S has a minimum at t m .
The characteristic time for the dynamics of bubble nucleation is thus α −1 m ∼ β −1 m . For the detonation case, the time β −1 * roughly characterizes the duration of the whole phase transition. In contrast, for the deflagration case, β −1 m only gives the time in which bubble nucleation occurs. The transition is longer, due to the decrease of both the nucleation rate and the wall velocity.
For a Gaussian nucleation rate and a constant wall velocity, the evolution of I(t) has been obtained analytically in Ref. [7] (for a very strong phase transition with v w ≃ 1). In the present case, the approximation v w ≃ v m is valid only for t ≃ t m . We shall not be interested in the evolution of f − (t) during a small interval in which f − ≃ 0. Nevertheless, in this short time the bubble size distribution is formed. In Eq. (33), t R is the nucleation time of a bubble which has radius R at time t. Since most bubbles nucleate around t = t m (where f + ≃ 1), we haveΓ(t R ) ≃ Γ(t R ) and v w (t R ) ≃ v m . Now, we must invert the relation R(t ′ , t) to obtain t ′ = t R (t), but we are only interested in times t ′ within a short interval around t m . Therefore, we may use the approximation , in which all the bubble sizes have approximately the same evolution R(t m , t), except for the initial dispersion, v m (t m − t ′ ). Inverting this relation, For any of our approximations (50) or (54), we see that dn/dR will be a function of and we see that for this distribution we haveR(t) = R(t m , t). For the distribution given by Eq. (50), the latter equality will be a good approximation at late times, due to the relatively small dispersion ∆R ∼ v m β −1 m . If a particular moment in the development of the phase transition is characterized by the value ofR, we may evaluate this distribution without even solving the evolution. For instance, at t = t F we haveR ≃ d (below we quantify the error of this approximation), where d takes its final value d f given by Eq. (52). Thus, we have a fully analytic approximation for the size distribution at the end of the transition, A similar expression is obtained for the sudden-reheating approximation. In Fig. 5 we compare these approximations and the numerical computation. We consider the normalized distribution n −1 dn/dR for a better comparison (so that order 1 errors in the analytic determination of Γ m cancel out). We see that the Gaussian nucleation rate gives a better approximation. We remark, though, that in both analytic approximations the position of the peak, R = d f , was computed from the sudden-reheating result (52).

Delta function rate
The evolution of the phase transition for t > t m is difficult to describe analytically. The integrals (28)- (29), involving R(t ′ , t) and v w (t ′′ ), are not easy to solve, except in the simple case of a constant v w , which gives R(t ′ , t) = v w (t − t ′ ). Even for the simple approximation (21) for v w (T ) and using the simple relation (42) between T , t and f − , we have an integrodifferential equation for R and v w . The main problem is that we have to deal, even after Γ turns off, with bubbles which nucleated at different times t ′ . Therefore, a considerable simplification is achieved by assuming that all bubbles nucleate at t = t m (we partially used this approximation in the previous subsection). In this case, all bubbles have the same radius, and we deal with a single variable R(t m , t). Thus, we consider a nucleation rate of the form where n f is the final number density of bubbles, which can be estimated from Eq. (51). Integrals like (32) and (30) are now trivial. We haveR(t) = 0 for t < t m andR(t) = R(t m , t) for t > t m . Thus, we may obtain a relatively simple equation forR(t) since we have for t > t m , where the temperature depends on t and f − . We have I = (4π/3)n fR 3 , so (with d f = n −1/3 f ). Even without solving the equation forR (which we do in the next section), we may check that this approximation is suitable to estimate quantities at later times. According to Eq. (63), the equalityR = d f occurs for f + = exp(−4π/3) ≃ 0.015. This happens before the time t F , which is defined by the condition f + = 0.01. In this approximation, the latter corresponds toR ≃ 1.032d f . Thus, we see that the equalityR = d occurs very close to the time t F . This is in agreement with the numerical computations of Sec. 2 and justifies the approximationR(t F ) ≃ d f in Eq. (60).

Analytic calculation of the evolution
We shall now compute the evolution of the phase transition for t > t m with the approximation of a delta-function nucleation rate. For this aim we consider Eq. (62) with the approximation (21) for the function v w (T ), which we write in the form On the other hand, in the present approximation, the linearized equation (42) gives the and where we have defined the quantity which parametrizes the reheating. Thus, we have

General behavior
From Eq. (68) we see that the evolution of the phase transition depends on a few parameters, namely, the wall velocity v m , the time t m − t c , the bubble separation d f , and the parameter q. Under the present approximations, this parameter gives the ratio of the released energy to the energy which is needed to reheat the system from T = T m back to Therefore, we may expect qualitative differences for q < 1 and q > 1. Indeed, in our two numerical examples, we have q ≃ 0.48 for the case r = 0.025 and q ≃ 1.89 for the case r = 0.05, and the difference is clear in Figs. 1 and 2.
For q < 1, we see from Eq. (65) that the temperature reached during reheating is bounded by T < T m + q(T c − T m ). This maximum value can only be reached if the second term on the right-hand side of Eq. (65) is negligible, for which the variation of f − must be very rapid. Initially, this is the case. For t close enough to t m , Eqs. (62) and (63) give which has a variation of order 1 in a time t − t m ∼ d f /v m . From our previous results, this is typically ∼ 10β −1 m ≪ t m − t c , so the last term in Eq. (65) can indeed be neglected during a time of this order. The approximation (70) will break down only if v w decreases significantly from its maximum v m . According to Eq. (67), the wall velocity will decrease at most by a factor 1 − q. Except in the limit q ≃ 1, this will be an order 1 factor, and the variation of f + will occur in a time At t = t F , the temperature will reach the maximum T r ≃ T m + q(T c − T m ) ≃ T m + (r/3)T c , and the wall velocity will reach a minimum v r ≃ v m (1 − q). For t > t F , the temperature is given by On the other hand, for the case q > 1, neglecting the term proportional to t − t m in (65) gives a reheating T r > T c . Nevertheless, this term cannot be neglected in this case, since, as T gets close to T c , the wall velocity (64) decreases significantly, and the growth of f − slows down. Although we will still have a variation f − ∼ 1 in a time of order d f /v m , according to Eq. (67), when the fraction of volume reaches a value f − ≃ 1/q we will have v w /v m ≪ 1. At this point the phase transition enters a phase-equilibrium stage at T ≃ T c . In this stage we have, from Eq. (65), the evolution Hence, the condition f − = 1 is reached for t F − t m ≃ (q − 1)(t m − t c ). For t > t F , the temperature is given by

Semianalytic computations
We have reduced the calculation of the phase transition dynamics to a relatively simple equation for the average bubble radiusR(t), Eq. (68). Although this equation cannot be integrated analytically in the general case, it represents a considerable simplification for the numerical computation. We shall now compare its solution with the results of the more complete treatment of Sec. 2. All the parameters appearing in Eq. (68), as well as the initial condition forR, can be estimated with the analytic and semi-analytic equations derived in Sec. 3. We may use the initial conditionR(t m ) = 0, which is consistent with the approximations (68). We have checked that this is in very good agreement with the results of Sec. 2. However, an even better agreement is obtained if we use as initial condition for the average radius the value given by Eq. (53),R(t m ) = v m /β m . Similarly, from Eqs. (40) and (45) we may obtain an initial value for the fraction of volume, These values ofR and I are consistent for an exponential nucleation rate but not for a delta-function nucleation rate 8 and, hence, they will give a slightly different evolution.
Here, we shall show the results for f + corresponding to the initial condition (72) and the results forR corresponding to the initial condition (53). In Figs. 6 and 7 we plot the fraction of volume and the temperature. The solid curves are those of Figs. 1 and 2, while the dashed curves correspond to the solution of Eq. (68). The other curves correspond to analytic approximations described below. We see that the simultaneous nucleation is a very good approximation in these cases. We remark that we have used also the analytic approximations of Sec. 3 for the parameters in the equation and the initial conditions, which introduce errors as well. As expected, the maximum error occurs at the end of the phase transition. This is better appreciated in the temperature curves.
In Fig. 8 we show the evolution ofR. The numerical computation described in Sec. 2 (solid line) is plotted from t = t N to t = t F . The dashed line is the solution of Eq. (68) and is plotted from t = t m to t = t F . The time t F is different in each curve, since it is defined by f + (t F ) = 0.01 for each calculation.

Analytic approximations
For t < t m , we have f − = 0, and Eq (65) is equivalent to T = T c − T c H c (t − t c ). For the average radius we haveR ≃ 0, and we may also use the valueR(t) ≃R(t m ), which is given by Eq. (53) and is indicated in Fig. 8 by a dash-dot-dot line.
To obtain analytic results for t > t m , we must make further approximations to Eq. (68). As already discussed, for the case q < 1, neglecting the difference t − t m could be a good approximation until the end of the phase transition. In this case, we have v m dt/dR = [1 − q(1 − e −I(R) )] −1 and, for small I, we obtain v m dt/dR = (1 + qI), which can be readily integrated. We have whereR m =R(t m ). The functionR(t) is thus obtained by inverting this quartic polynomial. The approximation (73) will break down for I 1. Nevertheless, in this case we have f + ≃ 0, so the error will be irrelevant for several quantities, such as the temperature or the velocity. Indeed, in Eqs. (65) and (67), at late times we have f − ≃ 1 and the evolution is given by the last terms. The curves are shown with a dotted line in Fig. 6. We see that, indeed, this analytic approximation is very close to the numerical result. The value ofR for this case is shown in the left panel of Fig. 8 (dotted line). As expected, it departs from the numerical result only at the end of the phase transition. Notice that, usingR = ( 3 4π I) 1/3 d f , analytic approximations can be readily obtained for the reference times t I , t P , t E and t F defined in Sec. 2.
On the other hand, for q > 1, we should not expect this approximation to remain valid until the end of the phase transition, since the approximation of neglecting the difference t − t m in Eq. (68) breaks down. Indeed, in Fig. 7 this analytic approximation (dotted lines) departs from the numerical computation as soon as the phase transition slows down.
For the subsequent slow stage, we may use the rough approximation T = T c , which gives the linear function of Eq. (71). This approximation, which is indicated with dashdotted cyan lines in Figs. 7 and 8, can be used until f + vanishes, where it can be matched to the approximation f + = 0. Although this rough approximation reproduces quite well the behavior of f + and T , it is not useful for the estimation of the wall velocity, since it corresponds to v w = 0. A better approximation can be obtained with a recursive trick. Notice that Eq. (71) is equivalent to f − = 1/q + (3H c /r)(t − t m ), which could have been obtained directly from Eq. (42) with the conditionṪ ≃ 0, which gives rḟ − = 3H c (the balance between the injection and extraction of energy from the plasma). In terms ofR, this condition is 4πṘR 2 f + /d 3 f = 3H c /r. Inserting it on the left-hand side of Eq. (68), we obtain the analytic relation This relation gives the black dash-dotted curves in Fig. 7 and in the right panel of Fig. 8.
Notice that the approximations are very good, except near the end of the phase transition.
In particular, the value of t F for the analytic approximation has a relatively large error.

Implications for cosmic relics
Although it is out of the scope of this paper to compute the possible relics from a phase transition, we wish to discuss on the implications of the dynamics we have just studied for their formation mechanisms. Several simplifications are often used in the literature.
As already mentioned, the most common approximations are a constant wall velocity and an exponential nucleation rate, and the results for the remnants of the phase transition depend on the free parameters v w and β * . In the case of small v w , the dynamics also depends on a few parameters such as v m and β m . In this case, however, the computation of cosmic remnants will be more involved, due to the non-trivial variation of the quantities. Below we consider two of the possible relics and discuss on the computation of the relevant quantities which are involved in their formation.

Topological defects
An important possible consequence of a cosmological phase transition is the formation of topological defects (see [34] for reviews). We shall consider for concreteness the case of cosmic strings. The simplest scenario in which these objects may arise is that in which a global U(1) symmetry is spontaneously broken at the phase transition. Inside each bubble, the phase angle θ of the Higgs field takes different values, and, when the walls of two bubbles collide, θ interpolates smoothly between these values [1]. As three or more bubbles meet, a total phase equilibration may not be possible due to topological obstruction. In such a case, the phase will change by ∆θ = 2π around the line at which the bubbles meet. At this line the Higgs field vanishes, and a cosmic string is formed.
In the case of a gauge symmetry, the phase difference is gauge dependent. One way of dealing with this issue is to use a gauge invariant phase ∆θ, defined as the line integral of the covariant derivative D µ θ = ∂ µ θ + eA µ [35]. In this case, during phase equilibration a magnetic flux is generated in the false vacuum region near the intersection of two colliding bubble walls. When a third bubble arrives, the fluxes corresponding to each pair of bubbles combine and, if there is a total phase change of 2π, a flux quantum is trapped inside the string. An additional mechanism of string formation is due to the presence of magnetic fields before the phase transition, which can be produced by thermal fluctuations [36]. After the phase transition, this magnetic field will be trapped in quantized flux tubes. If the phase transition is quick enough, this mechanism may produce a larger density of strings [37], or strings with higher winding number [38].
In either case, by the end of the phase transition, a random network of cosmic strings with some characteristic length scale ξ is expected to be formed. The statistical properties of such a network were studied in numerical simulations with cells of size ξ [39,40,41]. These calculations give, for instance, the proportions of closed loops and infinite strings. However, the characteristic length ξ is in principle given by the separation between nucleation centers, which is not a constant. In the first place, bubbles nucleate at random points, so the separation between neighboring bubbles has a dispersion ∆d around its averaged, with ∆d ∼d. In the second place, bubbles nucleate at different times, and those which were nucleated at the beginning of the phase transition are larger than those which were nucleated near the end. As a consequence, we will have inhomogeneities in the average separationd. In particular, for an exponentially growing nucleation rate, regions which were converted later to the broken-symmetry phase contain a much larger number density of bubbles. In contrast, for a slow phase transition, all bubbles nucleate in a relatively short time, and we expect a homogeneous average separation.
Similarly, for the bubble size distribution, in the exponential nucleation case we have at the end of the phase transition, while for the case of slow deflagrations the dispersion is quite smaller. Indeed, according to Even for a single scale d, the final characteristic length depends on the dynamics of phase equilibration or magnetic field diffusion, and hence on that of bubble nucleation and growth. If the phase angle in each bubble is uniformly distributed between 0 and 2π, the probability of trapping a string between three bubbles is 1/4. This gives a string length density of order 1/(4d 2 ). If bubbles expand at approximately the speed of light, this is a good approximation, but for v w ≪ 1, phase equilibration between two collided bubbles may complete before the wall of a third bubble reaches the meeting point, thus reducing the probability of trapping a string. For a gauge theory, the evolution of the phase difference is related to the spreading of magnetic flux [35], and the conductivity plays a role in the process. When two bubble walls collide, the magnetic flux generated at their intersection will spread in the symmetric phase. If part of the flux escapes to distances greater than the bubble radius before a third bubble arrives, the probability of defect trapping will be suppressed. Regarding the mechanism of flux trapping of already existing magnetic fields, it has not been much investigated for first-order phase transitions. Nevertheless, it is clear that the density of defects will be smaller for slower phase transitions.
Different kinds of simulations have been performed (mainly in 2+1 dimensions; see e.g. [42,43,44]) to study the dependence of defect formation on the dynamics of the phase transition. In these simulations, a constant wall velocity as well as a constant nucleation rate were assumed. As already discussed, the latter is generally a bad approximation. If we assume that such a constant rate turns on at a certain time t 0 and then takes a value Γ(t) = Γ 0 , we obtain a fraction of volume f . This result is qualitatively different from both the detonation and the slow deflagration cases.
For the deflagration case, a simultaneous nucleation at t = t m ≃ t I is a good approximation and is simpler than a constant rate. Unfortunately, in this case v w changes during the phase transition. Nevertheless, the dynamics is simplified by the fact that all the bubbles have similar sizes, and our analytic approximations may be useful in the calculation. Without entering into the details of the formation mechanisms, we notice that, although a common feature seems to be that smaller wall velocities reduce the probability of trapping defects in bubble collisions, the bubble separation is also smaller for lower velocities, d ∼ v m β −1 m , so the characteristic time between successive collisions, δt ∼ β −1 m , is rather independent of the wall velocity.

Electroweak baryogenesis and baryon inhomogeneities
The generation of the baryon asymmetry of the universe (BAU) may occur in the electroweak phase transition (see [45] for a review). The mechanism requires a first-order phase transition. In front of the walls of expanding bubbles, chiral asymmetries in particle number densities are generated due to C-and CP-violating scattering processes at the interfaces. These asymmetries bias the baryon-number violating processes (the sphalerons) in the symmetric phase. A net baryon number density is thus formed and enters the bubbles, where baryon number violation is turned off. A successful electroweak baryogenesis requires sufficient CP violation as well as a strong enough phase transition. The latter requirement is expressed quantitatively by the condition φ − (T )/T 1, which guarantees that sphaleron processes are suppressed in the broken-symmetry phase, thus avoiding the washing out of the generated BAU.
This mechanism has also an important dependence on the wall velocity. The chiral densities formed in front of the walls will be larger for higher velocities. However, the walls must also be slow enough for the sphalerons to have enough time to produce baryons. As a result, the generated BAU peaks for a certain wall velocity v w = v peak , which depends on the interaction rates and diffusion constants, and is generally in the range 10 −2 v peak 10 −1 (see, e.g., [46,47,48,49,50]).
In general, computations of electroweak baryogenesis for specific models focus on the sources of CP violation and on the condition φ/T > 1, and assume some (fixed) value for the wall velocity. On the other hand, the velocity of the electroweak bubble wall has been investigated for several models (see, e.g., [22,51,52,53,54]). Such computations generally focus on the determination of the friction of the wall with the plasma, taking into account the hydrodynamics of an isolated bubble. The resulting v w depends on the temperature T outside the bubble, which for application to specific models is usually evaluated at the onset of nucleation, T = T N .
As we have seen, the wall velocity may vary significantly after the time t N , especially if its initial value is in the range which is favorable for baryogenesis (as in our numerical examples). Depending on the model, the effect of this velocity decrease may be either an enhancement [5] or a suppression [14] of the generated baryon number density n B . Indeed, for v w > v peak we have roughly n B ∝ v −1 w , while for v w < v peak we have roughly n B ∝ v w . Thus, if the initial velocity is lower than v peak , the decrease of v w will cause a suppression, while if the initial velocity is (sufficiently) larger than v peak , we will have an enhancement.
In either case, a consequence of the velocity variation during baryogenesis is the formation of baryon inhomogeneities [5,12], due to a varying baryon number density which is left behind by the moving walls. A spherically-symmetric density profile is formed inside each bubble (at least, until bubbles meet each other). Since all bubbles nucleate almost simultaneously (at t = t m ), the inhomogeneities will have similar sizes and profiles. A density n B (v w (t)) is generated at a distance r = R(t m , t) from the bubble center. At the bubble center the value is approximately given by n B (v m ), and at a distance R(t m , t) there will be either an enhancement or a suppression by a factor v w (t)/v m (if v m is far enough from v peak ). Hence, since R(t m , t) ≃R(t), the profile n B (r) inside the bubble is essentially given by the parametric curve of v w (t)/v m vs.R(t).
The effect of reheating on the BAU and the formation of baryon inhomogeneities were investigated numerically in Refs. [5,14,12,15]. In Fig. 9 we plot the wall velocity vs. the average bubble radius, which gives an idea of the inhomogeneity profile, for our numerical, semi-analytic, and analytic results. Solid lines correspond to the complete numerical computation, dashed lines correspond to the delta-function rate, and dotted lines correspond to the analytic approximations. In the latter case, the curves are easily  Fig. 8, we observe that the analytic approximations depart from the numerical results near t = t F . We remark that the significant error observed in the right panel of Fig. 9 occurs inR(t) and not in v w (t), since Eq. (74) breaks down for small values of f + , and therefore does not have an impact on Eq. (67).
It is useful to find also simple formulas for the essential features, such as the amplitude of the inhomogeneities and their characteristic size. In Ref. [12], these basic characteristics were obtained as functions of model parameters, using rough approximations such as the analytic nucleation rate (34) and the estimation T m ≃ T N . We shall now obtain expressions in terms of the quantities v m , β m , and q.
The maximum size scale of the inhomogeneities is the final bubble size d f , which is given by Eq. (52). This value is indicated by the solid vertical lines in Fig. 9. However, the profile may have a shorter variation, depending essentially on the ratio q of Eq. (69). For q < 1, the maximum temperature T r is reached at the end of the phase transition, so the minimum velocity v r is taken at the boundaries of the inhomogeneities. Therefore, the characteristic length of the profile is the whole length d f . In this case, we have a variation v r /v m ≃ 1 − q ∼ 1. This value of v r /v m is indicated by the horizontal green line in the left panel of Fig. 9.
In the case q > 1, the wall moves with velocity v w ∼ v m until the fraction of volume becomes f − ∼ 1/q. After that, we have a much smaller velocity v r until the end of the phase transition. Thus, the baryon profile is composed of two main parts; namely, we have roughly a value n B (v m ) inside a sphere of radiusR r , while between this radius and d f , we have roughly a value n B (v r ). An estimate of the minimum velocity is obtained by setting f − = 1/q in Eq. (67), which gives v r /v m = (t r − t m )/(t m − t c ). The corresponding time t r − t m is given by Eq. (74) for I = log[q/(q − 1)]. We thus obtain 9 v r ≃ 3 4π This value is indicated by a green horizontal line in the right panel of Fig 9. The valuē R r corresponding to this value of I is indicated by a green dashed vertical line in Fig. 9. In this approximation, the reheating temperature T r is given by Taking into account that r is proportional to q, if we vary this parameter we see that v r decreases roughly as q −1 .
On the other hand, the sizeR r decreases more slowly.

Conclusions
The dynamics of a phase transition which proceeds by the growth of deflagration bubbles is quite different from that of a phase transition mediated by detonation bubbles. The former is much more difficult to study, due to the reheating caused by the release of latent heat. Nevertheless, the limit of very small velocities, v w < 10 −1 , is relatively simple, as the quantities depend on a homogeneous temperature T (t). Thus, the distinctive characteristic of the dynamics is a homogeneous reheating during the phase transition, which causes the nucleation rate to turn off and the wall velocity to decrease. In spite of the aforementioned simplification, this scenario is still more complex than the detonation case. In the latter, an exponential nucleation rate and a constant wall velocity can be assumed, and the evolution of the phase transition can be solved analytically. In contrast, in the slow-deflagration case, the nucleation rate, the wall velocity, and the temperature are linked through non-trivial equations. This kind of phase transition has been considered in a number of works [5,11,12,13,14,15]. In Refs. [5,12,14,15] the development of the phase transition was computed numerically, while in Refs. [11,13] some features, such as the minimum and maximum temperatures reached during the transition, were estimated analytically. Only in the limit of a phase transition in equilibrium at T = T c the development can be solved exactly [31]. If the latent heat is large enough, the temperature eventually gets very close to T c , and the subsequent stage can be approximated by that limiting case. However, this phaseequilibrium stage is only a part of the evolution, and does not always occur. The main aim of the present paper was to find analytic approximations for the nucleation rate, which allow to solve analytically the development of the phase transition in the general case, and to contrast the results with a complete numerical computation for a physical model.
During supercooling, the nucleation rate can be approximated by an exponential, like in the detonation case. However, when reheating begins bubble nucleation quickly turns off. With the approximation of a sudden reheating, we have found a simple semianalytic equation for the minimum temperature T m . Bubble nucleation occurs in a short interval around the corresponding time t m . Assuming a constant wall velocity v w = v w (T m ) in this interval, we obtained analytic approximations for several quantities, such as the final average bubble separation (which is constant for t t m ). We have found that the quantities estimated at t = t m are in very good agreement with the numerical computation. However, for the calculation of the bubble size distribution, a Gaussian nucleation rate (which gives a Gaussian distribution) is more appropriate.
In the evolution for t > t m , the wall velocity can no longer be taken as a constant, and the simplest realistic approximation is a velocity of the form v w ∝ T c −T , which is valid for a phase transition with little supercooling. Since the nucleation of bubbles is concentrated around t = t m , a great simplification is achieved by considering a nucleation rate of the form Γ ∝ δ(t − t m ). These approximations allowed us to obtain analytic solutions for the development of the phase transition after t = t m . The effects of reheating are encoded in the parameter q defined in Eq. (66). For q < 1, the temperature reaches a maximum T r ≃ T m + q(T c − T m ), and the velocity decreases only by a factor 1 − q. In this case, we have solved analytically the complete evolution. For the case q > 1, this analytic solution only describes the reheating stage, after which the temperature gets very close to T c and the velocity can decrease by a few orders of magnitude. For this longer phaseequilibrium stage, we have found a refinement of the usual approximation in which the phase transition develops at T = T c .
These analytic solutions depend on a few parameters which can be evaluated at t = t m . We have verified that these approximations describe remarkably well the evolution of the quantities which are relevant for the generation of cosmic remnants. This agreement with the numerical calculation shows in particular that a delta-function nucleation rate is a good approximation. This is interesting since, in numerical simulations, the bubbles are sometimes nucleated simultaneously for simplicity. For this kind of phase transition, this approximation is more appropriate than using an exponential nucleation rate. time (i.e., we consider the "envelopes" of bubble clusters). We assume that every surface element δA moves with a velocity v w perpendicular to the surface.
In the reference frame of a surface element, we assume that the fluid velocity is perpendicular to the wall. We assume deflagration conditions, in which the outgoing flow velocity v − is greater than the incoming velocity v + . Therefore, a portion of fluid which passes through the surface has a smaller entropy density but a larger volume in the − phase. In a time ∆t the entropy changes by where we have assumed non-relativistic velocities, and we have used the deflagration relation v − = v w . Integrating over the uncollided wall area, we obtain where V is the total volume 10 . On the other hand, we have and Eq. (78) givesṡ If we neglect the entropy increase, Eq. (80) gives the right-hand side of Eq. (26), and Eq. (79) is the left-hand side. In order to compare the size of the different contributions to the temperature variation, it is convenient to differentiate (79). We obtaiṅ To convert this equation into an equation for the temperature, notice thatṡ + /s + = 3Ṫ + /T + and, from Eqs. (18), we havė Using also the relation T − s − v − = T + s + v + in the last term of Eq. (81), we obtain We shall show that the two terms proportional to (T + −T − )/T − can be generally neglected.
We have also checked in our numerical computations that Eqs. (27) and (83) do not give appreciable differences. The former corresponds to neglecting the last term in (83) (the entropy-production term).
In the first place, we have (s + − s − )/s + ∼ L/ρ Rc . More precisely, from Eq. (18) we with α c ≡ α(T c ) = L/(4ρ Rc ) and α + ≡ α(T + ) = α c T 4 c /T 4 + . Therefore, for small L/ρ Rc and T + ≃ T c , we may expand (84) in powers of α, and we obtain plus terms of order (L/ρ Rc ) 2 and (L/ρ Rc )(T c − T + )/T + . In our numerical examples, these terms are of order 10 −4 , since we have L/ρ Rc ∼ (T c − T + )/T + ∼ 10 −2 . In the second term on the right-hand side of Eq. (83) we may neglect the part of order L/ρ Rc . In contrast, in the first term we cannot do so, sinceḟ − becomes much larger than H.
On the other hand, we know that (T + − T − )/T − is small, and we may neglect it in the left-hand side of Eq. (83). In contrast, in the right-hand side, the last term could be comparable to the first one. From Eq. (18) we have Hence, to lowest order we obtain and we may neglect the entropy-production term. With these approximations, Eq. (83) becomesṪ which is equivalent to Eq. (42).