{'fraction_non_alphanumeric': 0.031065455161715053, 'fraction_numerical': 0.015395846641202466, 'mean_word_length': 4.282694029336823, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating the Electroweak Phase Transition in the SU(2) Higgs Model

Numerical simulations are performed to study the finite temperature phase transition in the SU(2) Higgs model on the lattice. In the presently investigated range of the Higgs boson mass, below 50 GeV, the phase transition turns out to be of first order and its strength is rapidly decreasing with increasing Higgs boson mass. In order to control the systematic errors, we also perform studies of scaling violations and of finite volume effects.


Introduction
The masses of elementary particles in the Standard Model are generated via the Higgs mechanism by the non-zero vacuum expectation value of the scalar Higgs field. At high temperatures, above the scale of the vacuum expectation value, the Higgs mechanism is not operative, the symmetry of the vacuum gets restored [1]. In fact, in the early universe, according to the big bang cosmology, matter first existed in the symmetry restored phase. As a consequence of expansion and cooling, a non-zero vacuum expectation value of the scalar field was developed in the phase transition between the symmetric phase at high temperatures and the Higgs phase at lower temperatures. The properties of this electroweak phase transition might have a substantial influence on the later history of the Universe. For instance, since the number of baryons is not conserved in the minimal standard model [2], the small baryon asymmetry of the Universe could perhaps be created in non-equilibrium processes during a strong enough first order electroweak phase transition [3,4]. This offers the possibility that the baryon asymmetry can be explained within the minimal standard model. The resolution of this question is therefore a major challenge for elementary particle physics.
The standard calculational method for the study of the symmetry restoring electroweak phase transition is resummed perturbation theory [5,6,7,8]. In the Higgs phase perturbation theory is expected to work well for not very high Higgs boson masses, since the couplings are small. In the high temperature symmetric phase, however, the situation is similar to high temperature QCD: irreparable infrared singularities occur which prevent a quantitative control of graph resummation [9]. Since the calculation of physical quantities characterizing the phase transition requires the knowledge of both phases, there is a priory no reason why perturbation theory could provide a quantitative treatment of the electroweak phase transition. Indeed, the results of perturbation theory show bad convergence.
For a better understanding several non-perturbative methods have also been tried. For simplicity, fermions and the U(1) gauge field are often omitted. This can be expected on general grounds to be a reasonable first approximation. In this way one is left with the SU(2) Higgs model describing the interaction of a four-component Higgs scalar field with the SU(2) gauge field. Possible non-perturbative approaches include a block spin procedure leading to evolution equations for average actions [10], the ǫ-expansion at 4−ǫ spatial dimensions [11] and, of course, numerical simulations. After pioneering works [12,13], recent numerical simulations concentrated on the understanding of the finite temperature behaviour of the SU(2) Higgs model at large Higgs boson masses near and above the W-boson mass [14]. Another non-perturbative approach is based on dimensional reduction, studying the three-dimensional effective Higgs theory, which is obtained in the high-temperature limit [15,16,17]. A further simplification leads to an effective scalar theory [18], which has also been studied numerically in the reduced model [19].
The non-perturbative investigations of the electroweak phase transition did not yet lead to a convincing unique picture. Therefore, we decided to perform a large scale numerical simulation of the symmetry restoring phase transition in the SU(2) Higgs model. We stay in the original four-dimensional theory without reduction. This has the advantage of keeping the number of bare parameters small and not introducing any further approximations beyond the lattice regularization. First results have been published in a recent letter [20]. Here we give a detailed description of the techniques used and include additional results. As it is known from previous studies [14], for Higgs boson masses near and above the W-boson mass the numerical simulations in the original four-dimensional model are technically difficult. Therefore we restrict the present calculations to smaller Higgs boson masses below 50 GeV. Since this region of parameters of the minimal standard model is already excluded by experiments, our present scope is merely theoretical because we would like to check the validity of some other theoretical approximation schemes, e.g. resummed perturbation theory. We plan to extend this investigation to heavier Higgs boson masses in future papers.

Lattice action
The lattice action of the SU(2) Higgs model is conventionally written as Here U xµ denotes the SU(2) gauge link variable, U pl is the product of four U's around a plaquette and ϕ x is a complex 2 ⊗2 matrix in isospin space describing the Higgs scalar field and satisfying The bare parameters in the action are β ≡ 4/g 2 for the gauge coupling, λ for the scalar quartic coupling and κ for the scalar hopping parameter related to the bare mass square µ 2 0 by µ 2 0 = (1 − 2λ)κ −1 − 8. Throughout this paper we set the lattice spacing to one (a = 1), therefore all the masses and correlation lengths etc. will always be given in lattice units, unless otherwise stated.
In order to fix the physical parameters in a numerical simulation one has to define and compute some suitable renormalized quantities at zero temperature. The physical Higgs mass M H and W-mass M W can be extracted from correlation functions of different quantities (see section 5). The renormalized gauge coupling can be determined from the static potential of an external SU(2) charge pair, measured by Wilson loops (see section 6).
Since we are interested in the study of the symmetry restoring phase transition as a function of temperature, we use asymmetric lattices: the small temporal extensions L t = 2, 3, . . . represent the discretized inverse temperature L t = 1/(aT ). The other three (spatial) extensions of the lattice have to be much larger, for reaching the thermodynamical limit. In order to fix the physical parameters at the phase transition we have to determine the zero temperature renormalized parameters at the phase transition points for the L t = 2, 3, . . . lattices. A renormalized gauge coupling near the physical value g 2 R ≃ 0.5 can be obtained near β = 8 [21]. As stated before, we would like to have lighter Higgs boson masses than studied in [14]. Hence for the bare quartic coupling we have chosen values near λ = 0.0001 (in this paper referred to as low) and near λ = 0.0005 (referred to as high). In the present paper the inverse temperature in lattice units will be restricted to L t = 2 and L t = 3. Therefore the indices of the four sets of numerical simulations will be • l2 for low λ and L t = 2; • l3 for low λ and L t = 3; • h2 for high λ and L t = 2; • h3 for high λ and L t = 3.
In the next section the numerical simulation methods will be discussed. An important tool for the orientation in bare parameter space will be the invariant effective potential introduced in section 3. Then different groups of numerical simulation results will be discussed: the location of the phase transition points in section 4, masses and correlation lengths in section 5, the renormalized gauge coupling and the renormalization group trajectories in section 6, the latent heat in section 7 and finally the interface tension in section 8. The last section is devoted to the discussion of results and to a summary.

Monte Carlo simulation
In this section some aspects of the applied Monte Carlo simulation techniques are discussed. This can be skipped by readers not interested in technical details.
The simulations have been performed on the Alenia Quadrics computers of DESY. The Quadrics Q16 is a massive parallel machine with SIMD 1 architecture which consists of 128 processors (nodes). Depending on the goals and features of the respective simulation, we use different strategies: • A lattice is assigned to each node. No time is wasted for the communications between the nodes. Limitations of memory allow this only for small enough lattice extensions.
• The Q16 may be switched to consist of 16 independent 2 3 tori.
• The whole machine is arranged as a three-dimensional torus.
In this way the lattice is distributed over 1, 8 or 128 nodes and one obtains 128, 16 or 1 independent data sets from one run, respectively. Of course, the lattice extensions have to be multiples of the corresponding tori. The Quadrics offers 32 bit floating point arithmetics. This is sufficient for most of our purposes, except for building global averages, when we use a simple variant of software based double precision arithmetics for summation.

Updating
In a Monte Carlo simulation the autocorrelation of subsequent configurations is one of the main problems one has to deal with. The autocorrelation is usually worse at phase transitions. This phenomenon is called critical slowing down. Unfortunately, the continuum limit has to be performed in this region of parameter space (in our case κ, β and λ).
Due to the small λ-values we use, large fluctuations of the squared Higgs field length ρ 2 x ≡ 1 2 Tr(ϕ + x ϕ x ) occur in the Higgs phase and this expectation value shows the largest autocorrelation of all investigated quantities. In order to reduce the autocorrelation, the updating has to propose a big change of ρ x , with a high probability to accept it. This is achieved by an overrelaxation algorithm for ρ x , developed by two of us [22]. As usual, this overrelaxation algorithm is non-ergodic, thus we have to mix it with some ergodic updating steps.
A further improvement of the autocorrelation was achieved by the replacement of the ergodic Metropolis step considered in [22] with a heatbath algorithm proposed to us by Burkhard Bunk [23]. This algorithm works on the four real components of the ϕ x -field. As far as we know, this algorithm has not appeared in a publication, therefore we would like to describe it in some detail.
The idea is to divide the action in two parts, a quadratic one and the rest proportional to λ. In a heatbath step a new ϕ x is proposed according to the quadratic part of the action and the remaining parts are taken into account by an additional accept-reject procedure. Starting from the lattice action (1), the ϕ x -dependent part can be written in the following form: Here we used the notations In this form of the action a term proportional to ρ 2 x can be shifted from the quadratic part to the quartic part. We introduce a factor ζ x to express this freedom: For a good acceptance the minima of quadratic and quartic parts in (4) should coincide. So we get for ζ x the equation This cubic equation cannot be solved fast enough in an updating. Starting from the observation ζ x = 1 for |b x | = 1, we split ζ x in ζ x = 1 + ε x and get the following approximate expression for the optimal ζ x : This approximation works very well in a sufficient range of |b x |. In practice the average acceptance of this algorithm turned out to be larger than 98%. For the SU(2)-variables U xµ and α x ≡ ϕ x /ρ x we use standard overrelaxation methods [22,24]. Because the above described heatbath algorithm for ϕ x offers also an ergodic update for the angular part α x of the scalar field ϕ x , we need an ergodic update for U xµ only. For this purpose we use the heatbath algorithm described in [25,26].
In all updatings, the random number generator proposed and implemented by Martin Lüscher [27] is applied. It is based on an algorithm of Marsaglia and Zaman [28]. The latter algorithm is known by the name RCARRY, if the parameters have been chosen appropriately. RCARRY offers an extrem long period > 10 171 , but unfortunately it owns some short range correlations. As it has been shown [27], on long range, a chaotic nature of the algorithm comes to light. Skipping from time to time some hundreds of numbers in the sequence the correlation is practically eliminated. Due to the skip this random number generator is relatively slow. In order to benefit from the parallel architecture, the random number generator has to be initialized independently on the nodes of the machine.
For updating the field configurations a combination of the above described five algorithms is used. We choose some basic sequence of elementary updatings for the different sets of field variables, which is repeated periodically many times. The whole sequence, which visits every variable at least once but usually many times, will be called sweep. The optimization of this basic sequence making up a sweep is a difficult but important task, which will be discussed in the next subsection.

Autocorrelations
Let us consider the autocorrelation of a quantity Q[U xµ , ϕ x ] measured on a sequence of field configurations. Q n is the value of Q[U xµ , ϕ x ] measured on the n-th configuration and Q ≡ 1 N N n=1 Q n is the average over the configurations. The autocorrelation function for this quantity is defined as: By the use of the integrated autocorrelation time the statistical error can be estimated from the variance: The curves for Γ ρ 2 x (t) and Γ Cf (t) coincide.
We investigated the integrated autocorrelation time for four characteristic quantities: ρ 2 x , ρ 2 (x+Lt/2) ρ 2 x , U pl and the largest calculated Wilson loop. The second quantity characterizes the correlation function of ρ 2 at the largest distance. In zero temperature simulations the largest Wilson loop had the size L s /2 ⊗ L t /2, with L s and L t denoting the extensions of the lattice in time and space direction, respectively. At finite temperature only the 1 ⊗ 1 Wilson loop was considered. For Wilson loops not every orientation was taken: two sides were always in the direction of the largest lattice extension. Further on, we refer to the correlation function as Cf and to the Wilson loop as Wl. If equation (8) is evaluated on a finite sequence of configurations, one has to decide where to truncate the sum over t. As mentioned before, the largest τ int -values were found for ρ 2 x , so we truncated at the first zero of Γ ρ 2 x (t). On larger lattices we usually had 16 independent configurations in the computer. They were evaluated separately and an estimate for the statistical error of the integrated autocorrelation time was obtained from the variance.
We made some investigations how to optimize the autocorrelation by changing the number of calls of the various updatings in the complete sweeps but we did not try to optimize the ordering of the algorithms. Of course, we optimized the autocorrelation in CPU time since the above changes affect the time requirements of the complete sweeps. It should be mentioned that the optimal number of calls depends on lattice size and parameter range. The measurement routines were called after each sweep. Figure 2: Autocorrelation functions at the phase transition line on a 3 · 24 2 · 96 lattice. The left picture refers to the symmetric phase, the right one to the Higgs phase. Evaluated sweeps: 32000 sweeps in the symmetric phase and 80000 sweeps in the Higgs phase. In both cases the same updating scheme was applied. In the symmetric phase τ int,ρ 2 x = 21 ± 4 sweeps and in the Higgs phase τ int,ρ 2 x = 217 ± 66 sweeps.
The autocorrelation function was investigated in the Higgs phase both for T = 0 and for finite values of T . In the symmetric phase only finite T was considered. A typical example of autocorrelation in the Higgs phase is shown by fig. 1. The autocorrelation function Γ ρ 2 x (t) is to a good approximation a single exponential. There is no significant difference between Γ Cf (t) and Γ ρ 2 x (t). The autocorrelation functions for quantities depending only on U xµ show a fast fall off for t values very small compared to τ int,ρ 2 . For larger values of t the exponential descent of Γ U pl (t) and Γ W l (t) is the same as of Γ ρ 2 x . We always found the initial fall off to be larger for Γ W l (t) than for Γ U pl (t).
A comparison of autocorrelations in the two phases is given in fig. 2. In the symmetric phase at finite T we found the autocorrelation time for the quantities U pl and Wl to be less than 1. The left hand side of figure 2 displays the extremely fast descent of these autocorrelation functions. The behaviour of Γ ρ 2 x (t) differs from the behaviour in the Higgs phase, because there is a strong curvature in the logarithmic plot. This could be a signal for a dense spectrum of states. The integrated autocorrelation time in the symmetric phase is usually much smaller than the one in the Higgs phase at the same parameters and lattice extensions. Integrated autocorrelation time τ int,ρ 2 x of ρ 2 x for 5 different updating schemes. The lattice size is 16 3 · 32, the parameters are β = 8 and λ = 0.0005. The first item is at κ = 0.12885 and all the others are at κ = 0.1289. For each item more than 32000 sweeps were evaluated. On a 16 3 · 32 lattice with parameters β = 8, κ = 0.1289 and λ = 0.0005 we compared four different compositions of the complete sweep. These parameters give a point in the Higgs phase. The results for the largest autocorrelation time τ int,ρ 2 x are given in table 1 in sweeps and in CPU seconds of Q16, assuming that the lattice is distributed on the whole machine. A comparison of the third and fourth rows shows that more work on the SU(2) variables has no influence on the autocorrelation: the autocorrelation time in sweeps is about the same for both updating schemes. The fact that more overrelaxation for ρ x does not lead to a better autocorrelation is plausible [22].
The autocorrelation measured in sweeps decreases significantly if more heatbath is called but, due to the time needed for the heatbath algorithms, there is no significant difference between the second and the third row of table 1, if the autocorrelation time is measured in CPU time.
The updating scheme in the first row of table 1 is the best combination we found. It was therefore used in many points. Because of the large autocorrelations in the Higgs phase at finite temperatures, which were typically about 10 times longer in sweeps than at T = 0, we could not compare different updating schemes there. Another comparison of updating schemes was performed on a 18 3 · 36 lattice with parameters β = 8.15, κ = 0.1281 and λ = 0.00011. The conclusions were very similar.

Multicanonical simulation
An important problem of Monte Carlo simulations of a system with first order phase transition is the supercritical slowing down. At the transition point the tunneling rate between the two phases is exponentially suppressed for any local update algorithm (e.g. overrelaxation, heatbath). To overcome this problem the multicanonical algorithm was developed [29]. The basic idea is an enhancement of the mixed states, which are suppressed due to the additional free energy of the interfaces. This enhancement is reached by an extra term in the action, i.e. S → S + f (O). This term can be a function of any order parameter O. The easiest way is to use the action and a continuous function f (S) = β k S + α k with constant β k , α k for S in the interval I k = (S k , S k+1 ]. In fact, instead of the lattice action in eq. (1), we used as an order parameter the modified action which is natural to take if ρ x is used as an integration variable in the path integral. This choice is particularly convenient, since all the overrelaxation algorithms (ρ x , α x , U xµ ) can be used without changes. The intervals I k and the parameters α k and β k are chosen in such a way that the multicanonical probability distribution P mc L is nearly flat. This is achieved if f (S log ) ≈ log(P L ) between the two maxima and is constant elsewhere. Here P L is the canonical probability distribution of the action S log . The distribution P L is obtained in a multicanonical simulation by reweighting P mc L with exp(β k S log + α k ). In practice a first choice for the multicanonical parameters is made and they are optimized afterwards. If necessary, the procedure is repeated until P mc L becomes flat. A first guess can be obtained from smaller lattices. In this way the distribution of the action S log and the link variable L ϕ , defined in eq. (24), was measured on 2 · 4 2 · 64 and 2 · 4 2 · 128 lattices at the "low" value of the quartic coupling (λ = 0.0001).
For larger lattices two problems arise. The parameters have to be tuned very precisely and the autocorrelation times become even for optimally tuned values very large, of the order of O(10000) sweeps. To solve these problems we combined the multicanonical method with the constrained simulation method [30]. In what follows we call this way of simulation constrainedmulticanonical method.
We divide the interval between the two maxima of P L into subintervals. These are chosen to have an overlap with their neighbours. Starting in one phase we tune the multicanonical parameters such that P mc L is flat in a given subinterval and suppressed elsewhere. This means that f (S log ) is approximately equal to log(P L ) in this subinterval and increases rapidly beyond the boundaries. By moving the subinterval one is going from one maximum to the other. At the end every set is reweighted. In case of large overlaps between neighbouring intervals the absolute normalization can be obtained with small errors.
To ensure that this method yields the same result as the pure multicanonical one, we performed simulations using both methods on 2 · 4 2 · 128 lattice. The results coincide within statistical errors. Only the constrained-multicanonical algorithm has been used for simulations on 2 · 8 2 · 128 lattice.
The technical realization of the multicanonical approach by the Metropolis algorithm is straightforward. As it has been emphasized above, due to our special choice of the order parameter, the overrelaxation algorithms can be used without changes. The modifications for the heatbath algorithms are more involved (see e.g. [31]). A description of our implementation of heatbath algorithms for U xµ and ϕ x is given in the appendix. Since the heatbath algorithms are more efficient, we always used them instead of the Metropolis algorithms.
In our simulations the acceptance rate for the multicanonical heatbath algorithms was very good: for the modified gauge algorithm at least 99% and for the ϕ-algorithm at least 96%. The overlaps for the neighbouring intervals were chosen to be approximately 40%. The number of subintervals was 5 for the 2 · 4 2 · 128 lattice and 13 for the 2 · 8 2 · 128 lattice. The autocorrelation times were on average about 500 sweeps for the constrained-multicanonical simulations. For the two smaller lattices we measured about 2000 and 7000 sweeps as autocorrelation times with the pure multicanonical algorithm.
The easiest way to parallelize the multicanonical algorithm on the Quadrics Q16 machine is to simulate several lattices independently on each node. For any other implementation there is a need for communication between the different nodes for each updating step, since f is a function of the global action. Another disadvantage of partitioning the lattice would be a decrease in the acceptance rate due to simultaneous change of several variables.

Invariant effective potential
In the perturbative approach to the electroweak phase transition the most important quantity to compute is the effective potential. Interesting physical observables like latent heat, surface tension or masses can be extracted from it. Of course, these latter quantities can also be obtained from the non-perturbative approach of numerical lattice simulations by measuring suitable observables. Nevertheless, a direct comparison of the effective potential itself from both methods would obviously be desirable. On the lattice, however, the action is gauge invariant and so are the observables, as demanded by Elitzur's theorem. In perturbation theory the effective potential is calculated in different gauges. (The most popular one is the Landau gauge.) In order to compare this with lattice results one ought to fix the gauge on the lattice, a notoriously difficult task in particular for non-abelian gauge groups.
A way out is the study of the gauge invariant effective potential that has been initiated recently [32,33] 2 . In this approach one considers composite gauge invariant operators in the standard Legendre transformation framework. The obvious advantage of this approach is that the potential can be evaluated perturbatively and, since it is gauge invariant, it can be directly compared to lattice simulations. It offers therefore a conceptually clean and directly accessible tool of confronting results obtained in perturbation theory with numerical data. In this paper we want to report about our first experiences with the gauge invariant potential. We will use it mainly for the determination of the transition points. We postpone a discussion of its renormalization and the extraction of physical quantities to a future publication.
The starting point for the gauge invariant effective potential for the length square of the Higgs field is the free energy F (J) in the presence of a constant external source J with S the action eq. (1) and Ω the lattice volume. From this the effective potential is obtained by a Legendre transformation The perturbative evaluation of the gauge invariant effective potential is done in the standard loopwise semiclassical expansion. To proceed we choose the unitary gauge and set the angle of the Higgs field to α x = 1. Since the free energy is gauge invariant this has no effect on the effective potential. The action becomes with S g the plaquette action for the gauge field defined in (1). A stationary point is obtained for and U x,µ the unit matrix asρ The last equation is easily inverted for J(ρ 2 ) and the effective potential to tree level is To get the one-loop effective potential we have to consider fluctuations around the stationary point (16). The fluctuations to one-loop consist of a gauge part and a Higgs part and at this level no mixing appears. One obtains where the masses are related to the parameters in the lattice action (14) andρ 2 by and the momenta in the lattice integrals The solution of the effective potential given above corresponds to the broken phase of the SU(2)-Higgs model. In [33] it was emphazised that there exists another stationary point which belongs to the symmetric phase of the model and which is given byρ = 0. In this case the tree level potential is trivially zero and we have to start with the one-loop formula for the free energy with To obtain the effective potential one has to solve eq. (13) for J(ρ 2 ). In [33] the solution has been given for the three dimensional Higgs model in a closed form. A description of the effective potential in the symmetric phase has been found with a quite characteristic asymmetric shape.
In our case we have to work with lattice integrals or finite lattice sums. Then the solution can no longer be given in a closed form. However, one can perform the Legendre transformation and solve eq. (13) numerically. The result of this procedure for the points where our simulations are performed confirm the general shape of the potential in the symmetric phase and are qualitatively in agreement with the potential extracted from the distributions of ρ 2 values from the simulations. However, we do not have a quantitative understanding of the symmetric phase yet. We hope to come back to this question in a future publication.
A final remark concerns the lattice simulations where the gauge invariant effective potential is obtained from a distribution of the operator under consideration. The potential computed in this way is the so-called constraint effective potential [35]. In the infinite volume limit this potential coincides with the one defined by means of the Legendre transformation above. In perturbation theory both approaches differ in the treatment of the zero modes. For the oneloop result (18) this amounts to leaving out the k = 0 mode in the finite lattice sums for the constraint effective potential.

Phase transition points
A numerical simulation of the SU(2)-Higgs model should start by first determining the phase transition points. Physically the transition is triggered by a temperature change. Keeping all other parameters fixed, this could only be achieved on the lattice by asymmetric couplings, which is possible but cumbersome. It will become clear later (see sections 5 and 6) that in the parameter range we are interested in β and λ are fixing, to a good approximation, the renormalized parameters. A change of κ is reflected mainly in a change of the lattice spacing a. Therefore if one crosses the transition at fixed β, λ by changing κ, the essential change is in the physical temperature T = 1/(aL t ). (The physical volume is assumed to be large enough such that its change with a 3 is not important.) Thus we are looking for the phase transition in the hopping parameter at κ = κ c , for fixed β, λ.

One-loop effective potential and transition points
We found that for searching the transition point the gauge invariant effective potential can be very helpful. It can serve as a tool to provide quite accurate information about κ c which helps to select the κ-values where simulations are then performed. We define the transition point κ c as the κ-value where the symmetric and broken minima of the gauge invariant effective potential are degenerate. For the computation we used the one-loop formula eq. (18) for the broken phase and the trivial, ρ 2 = 0, minimum for the symmetric phase. Although this is certainly not the exact value for the symmetric minimum, we will see in the following that for small λ the transition κ ′ s are in very good agreement with numerical data.
For the computation of the gauge invariant effective potential on a finite lattice of size L x · L y · L z · L t the integrals in (18) have been replaced by the corresponding lattice sums. Following the experience in QCD [36], we also used the mean field improved gauge coupling g → g/ U pl , with U pl the measured plaquette value. This choice for the gauge coupling appeared to be very useful to achieve better agreement with simulation results.
At the values of the Higgs mass where our simulations are performed the electroweak phasetransition is certainly of first order. The effective potential is such that besides the absolute minimum a secondary minimum exists. In the simulations this phenomenon leads to hystere-sis effects when the system gets stuck in the wrong minimum. Therefore hysteresis effects in thermal cycles are often taken as an indication for a first order phase transition. The effective potential allows to compute the values of ρ 2 also in the false vacuum and should hence reproduce the hysteresis. In fig. 3 we show a thermal cycle at β = 8, λ = 0.0001 on a 2 · 4 2 · 32 lattice. The data points, connected by a long dashed line to guide the eye, show a clear hysteresis. The solid line represents the values of ρ 2 in the absolut minimum obtained from the gauge invariant effective potential. The short dashed line indicate the values of ρ 2 in the second minimum. The figure demonstrates that the numerical data are very well described by the one-loop gauge invariant effective potential. The agreement gets worse for the high point at λ = 0.0005. Here the transition points from the perturbatively evaluated potential and the simulations show some discrepancy, see table 2.
In section 7 we will compute the latent heat. For this we need transition κ ′ s for L t = 2, ..., 5. For the higher L t = 4, 5-values numerical simulations are very demanding as one would have to scale the other extensions of the lattice accordingly. Therefore we will resort there to the values of κ c as obtained from the effective potential. For this purpose we performed a finite size scaling analysis of κ L c on various size lattices where V = L x · L y · L z and L t is kept fixed. We computed κ L c for various V from the gauge invariant effective potential and fitted it to (21). In all the fits performed, we found a value of ν = 1.00 (2). This again confirms the first order nature of the electroweak phase transition for Higgs masses below 50 GeV. The obtained results for κ ∞ c for various L t , β and λ are plotted in fig. 12 where we discuss the lines of constant physics. The numerical values in our four basic points are: l2 :

Two-coupling method and transition points
Provided hysteresis effects in thermal cycles are seen at a first order phase transition, the twocoupling method is useful for a precise determination of the position of the phase transition point. Let us consider the largest extension of the lattice to be the z-direction. In this direction the lattice is divided into two halves. The idea is to choose different coupling parameters in both halves, to enforce one part to stay in the symmetric phase and the other one in the Higgs phase. We assume that the z-direction is long enough for accomodating a pair of interfaces between the phases. (If the z-direction is not too long, the appearance of more interface pairs has a negligible probability.) By warming up the configuration far enough in both directions from the transition point one can achieve that at the beginning of the simulation at the envisaged pair of parameters both phases and the interface pair are present. In this configuration the system can sensitively react to free energy differences by shifts of the interface positions in the Table 2: Hopping parameter κ c at the transition point, obtained with the two coupling method. As lattice size, the total size of both parts is given. For comparison, the estimates obtained on lattices with half the size from the one-loop invariant effective potential are also given. (2) 0.12826(1) z-direction. If the configuration stays for many autocorrelation times in the mixed state the free energies of the two phases at the chosen parameters are such that this situation is stable against transitions to a unique phase. Thus the parameter sets of the two phases give a lower and an upper bound on the transition parameters. In fact, due to the additional free energy associated with the interfaces, the phase transition from a two-phase situation to a unique phase occurs somewhat earlier than the equality of the free energies of the two phases. This parameter shift goes, however, exponentially to zero for increasing lattice volumes.
Of course, the two-coupling method for the determination of the transition point works only if one is able to tell one phase from the other. A possible way is to perform hysteresis runs on lattices with the same extensions as half of the lattice for the two-coupling method. The hysteresis plots are used for the distinction of the two phases at a given κ.
As stated in the introduction of this section, we are interested in the position of the phase transition in the hopping parameter at κ = κ c for fixed β and λ. Thus the only parameter chosen differently in both parts of the lattice is the hopping parameter. At each pair of κ the system was observed for at least 10 autocorrelation times. κ c was defined as the mean value of the best lower and upper bounds. The best estimates for the transition point κ c obtained with this method are given in table 2.

Multicanonical method and transition points
The most precise way to determine the transition point has been provided by a combination of the multicanonical method and κ-reweighting [37]. In this method first an order parameter distribution is generated through a Monte Carlo simulation. This distribution is then extrapolated to nearby κ's to find the transition point. The action with the additional logarithmic term (S log in eq. (10)) and the link-variable (L ϕ in eq. (24)) have been considered as order Figure 4: The distribution of the action density obtained from a constrainedmulticanonical simulation on 2 · 8 2 · 128 lattice. Note that the left hand peak corresponds to the Higgs phase and the right hand one to the symmetric phase.
parameters. The transition point is determined by the equal height signal: κ c is the hopping parameter value for which the heights of the two peaks in the distribution are equal. The systematic uncertainty of κ c can be estimated by comparing its values obtained by equal height signal and equal area signal for different order parameters.
The multicanonical simulation has been performed at κ = κ v , close to the real transition point (|κ v − κ c | = O(10 −5 )). The approximate transition point κ v was determined by the oneloop invariant effective potential. The distribution of an order parameter at a nearby κ ′ , with all other couplings fixed, can be obtained by attaching a weight to the different configurations. The action density distribution at the low point for κ v = 0.1283 on a lattice with Ω = 2·8 2 ·128 sites is shown in fig. 4. The distribution in neighbouring points is obtained by eq. (22). The two peaks are of equal height at κ c = 0.128307 ( fig. 5). As it can be seen, not only the heights are equal but also the widths of the two peaks are quite similar, thus the equal height condition for κ c is roughly equivalent to the equal area condition. At the same time the flat regime between the peaks is almost constant. This means that in a multicanonical simulation the two phases can mix with each other with an arbitrary mixing ratio. The supression is due to the interfaces between the phases. Taking L ϕ as an order parameter, the transition point corresponding to the equal height signal is given by a slightly higher κ=0.128308. This leads to a transition point at κ c = 0.128307(2 + 1), where the first number in brackets denotes the systematic the second one the statistical error estimate. Similarly for Ω = 2 · 4 2 · 128 and for Ω = 2 · 4 2 · 64 the transition points are at κ c = 0.128366(3 + 1) and κ c = 0.128367(3 + 1), respectively.

Masses and correlation lengths
Important characteristic features of any statistical physical system are the correlation lengths. At zero temperature their inverses give the masses of the low lying particles. Particularly interesting are also the (inverse) correlation lengths on the two sides of the phase transition.

Zero temperature masses
The physical Higgs mass M H can be extracted from correlators of quantities as the site variable or, using ϕ x ≡ ρ x α x , the link variables     (16)  48 The W-boson mass M W can be obtained similarly from the composite link fields (r, k = 1, 2, 3) Due to lattice symmetry only the diagonal correlators in (r, k) are non-zero and they can be averaged before extracting masses. Numerical simulations to determine the zero temperature masses were performed on 12 3 · 24 and 16 3 · 32 lattices at the phase transition points of L t = 2 lattices. The former were scaled up to 18 3 · 36 at the corresponding L t = 3 points. Since the gauge-Higgs system is in the Higgs phase at low temperatures and in the symmetric phase at high temperatures, the T = 0 points at the transition points are always in the Higgs phase. We did neither attempt to investigate the particle spectrum at T = 0 in the symmetric phase nor did we determine the κ-dependence of the masses in the Higgs phase in a wider range of κ. These questions could be addressed in future studies. The only information on κ-dependence is obtained from two points on 16 3  Statistical errors on masses were always estimated by subdividing the data sample into subsamples. Performing the fits in subsamples gives estimates of standard deviations of fit parameters. This is particularly straightforward on the Quadrics Q16, because these lattices can be simulated on 8 nodes which is repeated 16 times with independent sequences of random numbers. The 16 parallel sets of statistically independent results can be used for the estimate of fit parameter errors. The error estimates in last digits are always given in parentheses. The obtained results are shown in table 4.
The masses are in every case well determined. The correlators are always dominated by the lowest state. The statistical errors of the ratios of Higgs-and W-boson masses show that these two masses are not strongly correlated in the data samples. The values of T c /M W in table 4 are obtained under the assumption that the points are at the finite temperature transition points of L t = 2 and L t = 3 lattices, respectively. Therefore T c = 1/2 and T c = 1/3, respectively.
Points h2[16/85] and h2 [16/90] are at slightly different κ values: κ = 0.12885 and κ = 0.12890, respectively, but otherwise identical. Their comparison shows that the ratio of Higgsto W-boson masses remains almost constant in this κ-range. The same holds also for the renormalized gauge coupling (see next section). The only significant difference is a small change in the overall scale shown by an about 5% decrease of the masses in lattice units. This feature is rather helpful because it makes fine tuning of the scalar hopping parameter κ less relevant.
Comparing points h2 [12] with h2[16/90] we see that the finite volume effects are not large. The only statistically significant deviation is seen in the W-boson mass, which is about 5% smaller on the 12 3 · 24 lattice than on 16 3 · 32 (for the change of the static potential see section 6). This implies that within our statistical errors the 16 3 · 32 lattice can probably be considered to represent the infinite volume limit. Nevertheless a systematic study of finite volume dependences in future studies is desirable.

Correlation lengths at transition points
An important effect of finite temperatures is the change in correlation lengths. This is particularly interesting at a first order phase transition, where correlation lengths stay finite and may eventually differ in the two metastable states. The inverse correlation lengths at finite temperatures may also be called "temperature dependent masses". In order to distinguish them from the real masses we shall denote them by small letters: m W for the inverse correlation length in the W-boson channel, m H the one in the Higgs-boson channel.
In our numerical simulations we make use of the fact that on a large enough lattice the metastability at the first order phase transition becomes so strong that in a finite amount of computer time the system stays in the phase it started from. At the first sight this seems to lead to mathematically ill-defined averages, because the transition speed of the lattice from one phase to the other does certainly depend on the particular updating algorithm. However, as it is well known (and is discussed in detail in section 4) the two phases can be well distinguished in the distribution of some order parameter which shows two well separated peaks on a large enough lattice. The averages in the two metastable states can be defined to belong to the two peaks in the order parameter distribution. The advantage of using large lattices is that one can sit right at the phase transition, hence no extrapolation is necessary (see also section 7). This is more important than saving a little on lattice sizes and staying in an unclear situation with an underdeveloped two peak structure, when no clear separation of the two phases is possible.
The lattice shapes we consider are typically of the type L t · L 2 xy · L z with L t = 2, 3 the inverse temperature direction, L xy ≫ L t the "small" spatial directions and L z ≫ L xy the "long" spatial direction. The finite volume effects are mainly determined by L xy . The long direction L z is chosen to force the interfaces between the two phases to be built perpendicular to this direction. This is useful for the study of the interface tension (see section 8). The parameters of our numerical simulations are collected in table 5. The results are contained in table 6.
The behaviour of the correlators in the Higgs phase (points with index a) is in general similar to the zero temperature points discussed in the previous subsection. There is, however, a substantial difference in the symmetric phase: in general the determination of the inverse correlation lengths is much more difficult. The main reason is that the correlators are not dominated by the contribution of the lowest state. There is obviously a dense spectrum contributing in     (see table 5). Even in this case the lowest mass in the W-boson channel is still rather unprecise. For an illustration of the errors and fits see figures 8 and 9.
There is a clear signal for a discontinuity of m W and m H between the two sides of the phase transition, whereby the most significant change is an increase of m H when going from the Higgs to the symmetric phase. The values of m W,H are rather small, therefore the lattice volumes may be too small, especially at the l2, l3 points. There the infinite volume values may still be somewhat different. From this point of view the situation is much better at the "high" points, where the volumes had to be chosen much larger due to the smaller value of the interface tension.

Renormalized couplings
In the SU(2) Higgs model there are two independent dimensionless renormalized couplings which are conventionally defined at zero temperature. Corresponding to the bare gauge coupling g and bare quartic coupling λ one can introduce the renormalized couplings g R , respectively, λ R . One combination is conventionally fixed in the Higgs phase by the ratio of masses This corresponds to the relations of the masses to the renormalized vacuum expectation value: For the correct normalization of λ R one has to remember that it corresponds to the bare coupling in perturbation theory λ 0 , which is related to λ by and λ 0 is often multiplied by 24 = 4!. Note that for fixed gauge coupling g R , because of the Weinberg-Linde bound on the Higgs boson mass [38], R HW and λ R defined by (26) have positive lower bounds.
The mass ratio R HW has been determined in the previous section (see table 4). Now we shall introduce the renormalized gauge coupling g R in a way which is convenient for numerical simulations.

Renormalized gauge coupling
The static potential at distance R can be obtained from the Wilson loops W (R, T ) by We determined in the points defined by table 3 the on-axis Wilson loops of size R ⊗ T with 1 ≤ R ≤ L s /2 and 1 ≤ T ≤ L t /2 on L 3 s · L t lattices. Every rectangular Wilson loop with two sides in time direction was included in the statistics. These Wilson loops were evaluated after transforming the gauge configuration to temporal gauge. The time dependence has been fitted by three exponentials in order to determine the large time asymptotics. We checked that by leaving out one exponential and going to larger time distances the lowest energy V 0 does not change. In figure 10 we show a typical example for the behaviour of the logarithm of the ratio W (R, T − 1)/W (R, T ). As it can been seen from the given fit parameters, the ground state and the excited states are separated by a large energy gap. The situation for all other points is quite similar. Due to this, there is no need for an optimization of the source to enhance the ground state signal. The static potential can be fitted well by a Yukawa-term with lattice corrections, as discussed in [21]. It takes the form The parameter M is the screening mass which is closely related but not exactly equal to the physical W-boson mass M W determined from the correlation functions in section 5.1. G(M, R, L s ) describes the difference between the continuum potential and the finite lattice version to lowest order. For on-axis R values it is given by Figure 11: The continuum potential extracted (left picture) and the lattice correction term (right picture) in points l2 and l3 as a function of distance. M is the screening mass. The dotted curve is the best fit for l3. where In figure 11 we show the continuum potential V cont (R) = V (R) − C − D G(M, R, L s ) and the correction term δV (R) = D G(M, R, L s ) for the two points l2 and l3. One can see that on the L t = 3 lattice (point l3) the lattice artifacts are at larger physical distances smaller but at R = 1, 2 larger than on the L t = 2 lattice (point l2).
The renormalized gauge coupling can be defined from the parameters of the static potential as g 2 R = 16πA/3. In table 7 we collected the best fit parameters together with this global definition of g 2 R . Another local definition of g 2 R , which is shown in the last column, will be described below. Obviously the renormalization effect on bare coupling g 2 = 0.5 is moderate. The nearby equality of A and D shows that the potential is dominated by the one vector boson exchange [39]. One can see also that finite volume effects are small (points h2 [12] and h2 [16/90]) and that a small change in κ (points h2 [16/85] and h2 [16/90]) leaves g 2 R almost constant. Another way to define the renormalized gauge coupling was discussed in [21]. We use a slightly different version which is similar to the procedure recently proposed for pure gauge Table 7: Summary of the fit parameters for the static potential and the renormalized gauge coupling for both definitions.  (7) theory [40]. Given the potential, the renormalized coupling at distance R I is defined by This can be interpolated to a physical distance, for instance r = M −1 . R I is defined by the force as the solution of the equation The value of M in eq. (33) can be taken from the fit (29) restricted to large distances. The reason of introducing a local definition of g R is that the short distance potential is not purely Yukawa like. Specifying the distance is therefore in principle important. However, our resolution is up to now not fine enough to see the logarithmic dependence on distance. The time dependence of the Wilson loops as well as the distance dependence of the static potential were described very well by the ansätze, i.e. χ 2 was of order one in both cases. The statistical errors were determined from independent subsamples (see previous section). The error for the renormalized gauge coupling squared g 2 R (M −1 ) quoted in table 7 is the sum of the statistical errors for potential and mass, and the systematical error. The interpolation to the distance M −1 was done by a linear approximation of the two neighbouring points. The systematical errors are estimated from the change arising by the inclusion of a third point. The given errors are usually dominated by the statistical errors of M. Note that the value of g 2 R (M −1 ) in point l2 is corrected as compared to [20].

Renormalization group trajectories
The continuum limit of quantum field theories on the lattice is taken along renormalization group trajectories, also called lines of constant physics (LCP's). Going to the continuum limit along such lines in bare parameter space the same physical theory is reproduced with increasing precision. This means that the physical results of numerical simulations should be the same, apart from "scale breaking lattice artifacts". In the SU(2) Higgs model, in order to define the LCP's, one has to keep fixed the values of two independent renormalized couplings, say the above discussed renormalized quartic (λ R ) and gauge (g R ) couplings. The third parameter characterizing the points along LCP's can be chosen as (Remember that in this paper the lattice spacing is set to a = 1, therefore M W is measured here in lattice units.) Since our numerical simulations are performed in the weak coupling region of parameter space, the change of the bare couplings g 2 ≡ 4/β and λ 0 defined by eq. (27) along LCP's can be well approximated by the solutions of the one-loop perturbative renormalization group equations: The integration of these equations can be started at the transition points of the L t = 2 lattices. At a distance ∆τ = log(L t /2) (L t = 3, 4, 5, . . .) in parameter τ one obtains the corresponding values of (g 2 , λ 0 ) for the transition points of lattices with temporal extension L t . In order to have a perturbative prediction also for the change of the third bare parameter κ, one can make use of the one-loop perturbative invariant effective potential. The infinite volume predictions for our "low" and "high" values of the quartic coupling in the range 2 ≤ L t ≤ 5 are shown by fig. 12. As we discussed in section 4, the measured values of κ c are very well reproduced by the one-loop invariant effective potential for L t = 2, 3 at our "low" value of the quartic coupling. At the "high" quartic coupling there is a parallel shift between prediction and measurements: the prediction is too low by ∆κ c = 0.00025(1). As a consequence, the derivatives of κ along these LCP's can be well determined: The errors given here include an estimate of the systematic error by comparing different polynomial interpolations of the points 2 ≤ L t ≤ 5 (first entry in parentheses) plus the error coming from the uncertainties of κ c obtained in the numerical simulations (second entry).
Although the estimates of the LCP's by one-loop perturbation theory are very useful for guiding the numerical simulations, finally one has to confront them with numerical simulation data. Having this in mind we have chosen the L t = 3 points with index l3 and h3 on the trajectories of the solutions of eq. (36) which start at l2 and h2, respectively. (For β and λ we kept only the first few digits which are sufficient within our typical statistical errors.) Comparing the values of R HW and g 2 R in the points l2 and l3, respectively, h2 and h3 (see tables 4 and 7) one can conclude that the one-loop predictions from eq. (36) are very good: the measured values of both renormalized couplings in corresponding points are very close to each other, in most cases equal within errors.

Latent heat
An important characteristic feature of first order phase transitions is the latent heat, i.e. the discontinuity of the energy density ǫ. The pressure P is continuous, therefore ∆ǫ can be obtained from the discontinuity of δ ≡ ǫ/3 − P . This latter quantity is related to the trace of the energy-momentum tensor and is theoretically simple, because it can be obtained from the derivative of the action density with respect to the lattice spacing a. This, in turn, implies that it is closely related to the LCP's for the continuum limit discussed in the previous section.
As it has been shown in [20], the latent heat ∆ǫ in the SU(2) Higgs model is given by The obtained values of global averages "above" and "below" phase transition points, i.e. in the Higgs, respectively, symmetric phase.  (11) 6.517 (7) 6.61214 (4) Here, besides the notations in (24), we used The sign of the discontinuities like ∆ǫ etc. will be defined in such a way that they are differences of values in the symmetric phase minus Higgs phase. The obtained average values of the global quantities in eqs. (23), (24), (39) and of the lattice action per point are shown in table 8. The indices of points where the numerical simulations were performed are defined in table 5. As it has been discussed above (see e.g. section 5.2), the small errors in table 8 could be achieved by avoiding errors of extrapolation. This is made possible by performing the simulations on large enough lattices, where the strong metastability of phases can be exploited. In this way the contribution of the statistical errors in table 8 to the errors of ∆ǫ is negligibly small compared to the errors of the derivatives. In particular, in our case the errors in eq. (41) are dominated by the errors coming from eq. (37).
Besides the discontinuities of global quantities, for ∆ǫ in eq. (38) the derivatives of bare parameters along the LCP's are needed. These have been discussed in the previous section. ∂β/∂τ and ∂λ/∂τ can be obtained to a good approximation from eq. (36). The values of ∂κ/∂τ have also been evaluated from the one-loop invariant effective potential and from numerical simulation data. They are given by eq.
We see from here that within errors the results on L t = 2 and L t = 3 lattices coincide.
Comparing the contributions of different terms to ∆ǫ in (38) one sees that the term with L ϕ dominates. This is in accordance with the observation that the lattice spacing is mainly determined by the hopping parameter κ (see sections 4 and 5).

Interface tension
At the transition point of the electroweak phase transition mixed states can appear, where different bulk phases are separated by interfaces. The interface tension, σ, is the free energy density associated with these interfaces. The dynamics of a first order phase transition is to a large extent determined by the interface tension. In particular, the nucleation of bubbles or droplets is essentially influenced by it [41,42]. In perturbation theory the determination of σ is a task which is far from being trivial. The main problem is the bad infrared behaviour of the finite temperature effective action in the SU(2)-Higgs model. Therefore, non-perturbative lattice simulations are very useful.
In this section we use two different methods, namely, the two-coupling method and Binder's histogram method combined with multicanonical updating, to determine σ. In both cases we use elongated lattices. This geometric choice of the lattices results in interfaces which are perpendicular to the long direction. Other configurations (e.g. multiple walls, large spherical bubbles) are highly suppressed. The inverse temperature will be restricted in this section to L t = 2.

Two-coupling method and interface tension
One way to determine the interface tension in the SU(2) Higgs model is based on a modified version of the Potvin-Rebbi two-coupling method [43]. Similarly to subsection 4.2, a lattice is considered with a long extension in the z-direction. The long direction with periodic boundary conditions is divided into two equal halves with different couplings, in such a way that an interface pair is created near the hypersurfaces where the couplings are changing.
Since we have three bare parameters, in principle one can choose any of them (or some combination) to be different in the two halves of the lattice. A physically good choice would be to change the temperature only. One can move, for instance, along an LCP (see section 6.2) when only the lattice spacing a is changing. Assuming that the volume is very large, the only relevant change is in the temperature T = 1/(aL t ). A simpler way, which comes close to this, is to change only κ. The reason is that in the latent heat the contribution of the ϕ-link L ϕ dominates which is conjugate to the hopping parameter κ.
Therefore we decided to split the value of κ. On an L t ·L x ·L y ·L z lattice with L z ≫ L x,y,t , and in our case L t = 2, L x = L y ≡ L xy , the two halves in the z-direction have hopping parameters κ 1,2 which are slightly below and above the transition point: The interface tension is defined by the "free energy" (more precisely "free action") where S is the lattice action (1), Ω ≡ L x L y L z L t is the number of lattice points and d(s) is the spectral density of states as a function of the action density s ≡ S/Ω: This implies that the probability distribution of s is w(s) ∝ exp[−Ωf (s)]. The interface tension between the states with κ 1 and κ 2 is given in lattice units by Let us write the action on the lattice with two halves as where S 0 is the piece independent from κ. Then in the thermodynamical limit we have for if . . . κ 1 ,κ 2 denotes expectation values for a lattice with two halves at κ 1,2 . Therefore we obtain from (44) Here the notation σ(κ 1 , κ 2 ) emphasizes the dependence on κ 1,2 . Finally, since S 1,2 is given by the ϕ-link L ϕ , we have where L (1,2) ϕ (κ, κ ′ ) denote the expectation values of L ϕ averaged in the two halves, if the hopping parameters are κ and κ ′ , respectively.
In order to obtain the physically interesting interface tension one has, of course, to perform the non-commutative limits L x,y,z → ∞ and ∆κ ≡ κ 2 − κ c = κ c − κ 1 → 0. For a given lattice extension ∆κ cannot be arbitrarily small, because if the difference in free energy density becomes small tunneling into the other phase can occur and the interfaces disappear. The presence of the interfaces can, however, be monitored to ensure the applicability of (48).
For small ∆κ the integrals in (48) can be well approximated by the average of the integrand at the two end points. Besides, for equal arguments we obviously have L (1) ϕ (κ, κ) = L (2) ϕ (κ, κ). This gives We are interested in the limit ∆κ → 0. Therefore in the square bracket only contributions of order 1/∆κ matter. Such contributions cannot come from terms with equal hopping parameters in the two halves. Therefore let us introduce the parametrizations Then for ∆κ → 0 a finite volume estimator of the interface tension iŝ In view of this formula the question arises about the origin of the contributions of order O(1/∆κ). For small ∆κ the difference in free energies of the two phases is of the order O(∆κ). Therefore the interfaces are not fixed at the hypersurfaces where κ changes, but are penetrating into the neighbouring regions with constant κ. If the κ change is at z = z 0 , the probability of the interface position at z > z 0 is ∝ exp[−const. ∆κ(z −z 0 )]. Integrating over z gives O(1/∆κ). In fact, as the numerical simulation data show, the distributions of z-slices of L ϕ as a function of z can be well fitted by exponentials.
The above formulas also show that the extrapolation to ∆κ = 0 is delicate. The intervals in κ 1,2 for the fits in eq. (50) have to be carefully chosen. These forms are, in fact, replacing (48) in the region where σ(κ 1 , κ 2 ) depends linearly on κ 1,2 . For given lattice extensions one cannot take intervals too close to κ c . This is a kind of round-off effects, which usually appear at first order phase transitions in finite volumes but get smaller for increasing lattice sizes. First, according to eq. (51) one has to choose L z large enough in order that the two interfaces do not interact with each other. In this case the left hand side becomes independent from L z . Second, for increasing L x , L y the free energy differences grow proportional to (L x L y ∆κ). The probability of penetration into the regions with constant κ is exponentially decreasing. This allows to consider small ∆κ's of the order of O(1/(L x L y )).
In the numerical simulations with two κ's the lattice sizes extended up to 2 · 16 2 · 128 in the l2 point and 2 · 32 2 · 256 in the h2 point. The numerical simulation data on these lattices are collected in table 9, respectively, table 10. Fits with the parameters a, b, c in eq. (50) were performed, omitting quadratic and higher order terms in ∆κ. For the transition points κ c = 0.12830 and κ c = 0.12887 were assumed in points l2 and h2, respectively. The data entries in the tables with italic were not taken into account in the fits, in order to obtain acceptable Table 9: The average values of L (1,2) ϕ obtained in two-κ simulations with κ = κ 1,2 in point l2 on 2 · 16 2 · 128 lattice.  = 0.008 (2) .
The errors here also contain some subjective estimates of the systematic errors coming from the choice of fit intervals. The statistical part of the errors was estimated by repeating the fits with normally distributed input data.
Concerning the volume dependence ofσ we also collected data in point l2 on 2 · 8 2 · 128 lattice. Within statistical errors of about 10% no deviation from the result given in (52) could be detected.

Multicanonical method and interface tension
In this subsection we present our results for the interface tension obtained by the histogram method [44] at the low point, L t = 2.
At κ c the probability distribution of an order parameter (e.g. action density s log ≡ S log /Ω defined by (10), or link variable L ϕ in (24)) develops two peaks. They correspond to pure phases and the suppressed configurations between the peaks are dominated by mixed states where the phases are separated by interfaces. Defining κ c by the equal height signal the suppression at infinite volume is given by the interface tension where p max corresponds to the heights of the peaks and p min to the minimum in between.
Recently intensive studies have been carried out in order to understand the finite size corrections to σ [45,46,47,48,49]. We follow [48] and for our simulations at the "low" value of the quartic coupling λ = 0.0001 use elongated lattices with extensions Ω = 2 · 4 2 · 64, Ω = 2 · 4 2 · 128 and Ω = 2 · 8 2 · 128. In this case the finite size corrections are particularly simple where L x = L y ≡ L xy , L z is the longest extension and c is an unknown constant. The obvious advantage of this choice is that practically all mixed configurations contain two surfaces perpendicular to the z direction. More than two surfaces or surfaces in any other directions are suppressed by many orders of magnitude. (Corrections beyond the gaussian approximation [49] are negligible in our case.) The basic picture behind the above formulas assumes that the two interfaces are thin and do not interact. Fig. 13 is a typical picture of the average L ϕ values for different z-slices in a mixed configuration. As it can be seen, in the case of Ω = 2 · 4 2 · 64 the interfaces are very close to each other and they cannot be considered as thin (the spatial width in lattice units  is approximately 30). In this case there is no flat region between the two peaks of the S log distribution. The situation is much better for the even longer lattices (e.g. Ω = 2 · 8 2 · 128, see fig. 14). In the latter case the distribution of the action density s log at κ = κ c is shown in fig. 5. It is remarkably flat between the peaks. This is another sign that the interfaces are well separated, thus mixed configurations can appear with any bulk phase ratio.
In order to suppress statistical noise we have used a third order polynomial fit to the histograms near the extrema and calculated the extremum values from the fits. (However, due to our good statistics the simple use of the histogram-bins has resulted in almost the same results.) Since in case of multicanonical simulations one lattice has been attached to each node, we had 128 independent simulations. Our error estimates have been obtained by inspecting these independent samples.
The finite volume interface tensions turn out to be σ Ω /T 3 c = 0.63 (3) for Ω = 2 · 4 2 · 128 and σ Ω /T 3 c = 0.80 (2) for Ω = 2 · 8 2 · 128, respectively. Combining these two values one gets from This result shows that the infinite volume limit is not far away from our largest volume result and it is in very good agreement with the results of the previously discussed two-κ method in (52).

Summary and discussion
The main conclusion of this paper and ref. [20] is that the electroweak phase transition in the SU(2) Higgs model is of first order for Higgs boson masses below 50 GeV . The strongest indication for this is the two-peak structure of order parameter distributions (see e.g. fig. 5) showing a strongly suppressed flat region between the peaks which corresponds to mixed phases. For M H ≃ 50 GeV the phase transition became rather weakly first order. Still we could simulate large enough lattices to reach long living metastability. In fact, we used this phenomenon to study the properties of both metastable phases separately (see e.g. section 5.2). Investigations of the correlation lengths in both phases showed that they behave discontinuously at the phase transition. The strength of the first order phase transition can be characterized by the latent heat (section 7) and by the interface tension (section 8). Both these quantities play a central role in the dynamics of first order phase transitions and can be determined in numerical simulations. Our results show that at small Higgs boson masses (M H ≃ 20 GeV ) the transition is very strong, with a large latent heat (eq. (41)) and interface tension (eqs. (52), (55)). At Higgs boson masses near M H ≃ 50 GeV , however, the phase transition becomes weak: the dimensionless latent heat decreases by about one order of magnitude, the dimensionless interface tension by about two orders of magnitude (see the same equations). The rate of decrease is qualitatively the same as given by two-loop resummed perturbation theory [8]. For this we refer to the figures of ref.
[20] which are very similar after the inclusion of the new results here. Since this is our first large scale numerical simulation of the electroweak phase transition, we could not yet achieve a full quantitative control of the errors. This definitely requires further essential work. Nevertheless, our lattices are large and the statistical errors in most cases small. In the present paper first attempts were also made to estimate the systematic errors due to finite lattice spacing and finite volume. The comparison of the results on L t = 2 and L t = 3 lattices shows remarkable consistency at the 10 % level. This means that lattice artifacts are not large. In particular, the latent heat shows a remarkably good scaling between L t = 2 and L t = 3 lattices. The apparent general tendency of small scaling violations already at such small temporal lattice extensions is consistent with the expected dominance of the lowest Matsubara modes motivating the dimensional reduction [15,16,17].
Comparing the results for T c from the lattice simulations and those from perturbation theory, the values obtained from the simulations are always smaller. However, the proper comparison should be done not with continuum but with lattice perturbation theory, for relatively small L t 's. We have carried out this analysis (see fig. 15). In order to determine the critical points and masses the gauge invariant effective potential of section 3 has been used (again, the role of the symmetric phase and effects due to higher order corrections were neglected). One can see that for m H = 18 GeV the perturbative result is in complete agreement with the lattice one. A scaling violation can be observed and it goes approximately like 1/L 2 t ∝ a 2 . (See the inlet of fig. 15, where L t ≈ 5 − 7 practically reaches 2.47, indicated as a dashed line. This value corresponds to L t → ∞, which can be also obtained by the use of the one-loop effective potential without high-temperature expansion.) In the case of m H ≈ 49 GeV the situation is somewhat different. The perturbative values are larger than those obtained by the lattice simulations. The scaling violation is much smaller and within the errors it is reproduced by the perturbative treatment.
Up to now we did only a few checks of finite volume effects. These showed that the results are only changing at the few percent level. Further systematic tests have to be done in the future. Nevertheless, the lattices are large enough to support strong metastability of the two phases. Therefore the conclusions concerning the strength of the first order phase transitions are firm.
In general, numerical simulations of the electroweak phase transition in the four dimensional SU(2) Higgs model in the range of Higgs boson masses below M H ≃ 50 GeV turned out to be feasible and powerful. For these relatively light Higgs masses there is no need to go to the three dimensional reduced model. The situation could, however, be different for heavier Higgs bosons. In fact, recent numerical studies of the reduced SU(2) Higgs model are concentrating on M H ≃ M W [50].
In future numerical simulations in the four dimensional model one has to consider also L t = 4 lattices for confirming the smallness of scaling violations. A systematic study of finite volume effects is desirable. Of course, an extension towards heavier Higgs bosons would be very interesting.

A Appendix: multicanonical heatbath algorithms
Here we describe the modifications for the heatbath algorithms in multicanonical simulations. For sake of simplicity we first explain the general strategy in case of the gauge field and then the realization of the algorithms for the fields U and ϕ. We use the notations introduced in section 2.3. In general, the indices old and new will denote the value before the update and the proposal, respectively.
We distinguish two cases. If all possible values of the new action S new log remain in the interval I k = (S k , S k+1 ], one has to replace S by (1 + β k )S and the generation of the distribution is as usual (see section 2.1). If S new log crosses an interval boundary, the two distributions which are different for the two intervals have to be matched. This is done in such a way that for the interval I k the distribution becomes W k ∼ exp[−(1 + β k )S] and for I k+1 it is W k+1 ∼ exp[−(1 + β k+1 )S]. In addition the distribution for the combined interval I k and I k+1 have to be continuous. The idea is to construct a distribution which is a hull of W k and W k+1 . To obtain the desired distribution in each interval one has to correct the hull distribution by imposing an acceptreject step. The correction factor (denoted by Θ) has to be positive and not larger than one. Clearly if the hull distribution is close to W k and W k+1 the correction factor is about one and the acceptance rate is high.
This is done by the algorithm described in [26]. In the multicanonical case the first step is to check, as discussed above, whether the global action S new log can cross an interval boundary. This depends on a 0 through S new log (a new 0 ) = S old log + ξ(a old 0 − a new 0 ). If the minimal action S new log (1) and the maximal action S new log (−1) are in the same interval I k one replaces ξ by ξ(1 + β k ) in (56) and proceeds in the usual way. In case that the action crosses an interval boundary S k+1 , which defines the corresponding a bound If a 0 ≥ a bound 0 (S new log ∈ I k ) no correction factor is needed. For a 0 < a bound 0 (S new log ∈ I k+1 ) one has to impose Θ = exp[ξ(β k+1 − β k )(a 0 − a bound 0 )] as an additional accept-reject factor (see figure 16).
To avoid that the action S new log crosses more than one interval boundary, the intervals must be large enough.
In practice the values of β k 's are small (O(10 −2 −10 −3 )). Therefore the two curves are closer together than in figure 16 and the acceptance rate is high. Also the values of ξ are higher, i.e. there is a sharp peak near a 0 = 1.

A.2 Scalar field
The algorithm which is used for ϕ-updating was described in section 2.1. We write the action of eq. (4) as a sum of the quadratic and quartic part and a constant S(ϕ x ) = S (2) (ϕ x ) + S (4) (ϕ x ) + const. ′ .
For the multicanonical simulations one has to generate the distribution for the four components of ϕ x . The first step is to generate the distribution as a hull analytically by replacing S (2) (ϕ x ) → (1 + β min )S (2) (ϕ x ). β min denotes the minimal β k . Let I k denote the interval in which the minimum S 0 log of the action can be found. If S new log is in this interval I k , the correction factor is where the meaning of the constant c is described below. In case of crossing an interval border, i.e. S new log = S − 3 x log(ρ x ) > S k+1 , one needs an additional factor The total correction factor is in this case Θ 1 Θ 2 . This procedure generates the action distribution in I k and in I k+1 correctly. Due to the definition of I k , S new log cannot be in I j (j < k). Since our intervals are large enough the probability that S new log is in I j (j > k + 1) is negligible. In practice it is quite difficult to find S 0 log , but a slightly smallerS log < S 0 log can be used, which could lead to a somewhat smaller acceptance rate. In addition it is necessary that the whole correction factor (which includes all terms except the analytically generated quadratic term) has to be smaller or equal one. To be sure that no problem arises we subtracted a small positive extra term c in the exponent of Θ 1 . This term reduces the acceptance. It must be tuned in such a way that the correction factor is never greater than one and the acceptance rate is high.
Of course, there is always only one accept-reject step in both updating algorithms. This means that the analytically generated distribution will be corrected with the product of all factors.