How arbitrary are perturbative calculations of the electroweak phase transition?

We investigate the extent to which perturbative calculations of the electroweak phase transition are arbitrary and uncertain, owing to their gauge, renormalisation scale and scheme dependence, as well as treatments of the Goldstone catastrophe and daisy diagrams. Using the complete parameter space of the Standard Model extended by a real scalar singlet with a $\mathbb{Z}_2$ symmetry as a test, we explore the properties of the electroweak phase transition in general $R_\xi$ and covariant gauges, OS and $\overline{\text{MS}}$ renormalisation schemes, and for common treatments of the Goldstone catastrophe and daisy diagrams. Reassuringly, we find that different renormalisation schemes and different treatments of the Goldstone catastrophe and daisy diagrams typically lead to only modest changes in predictions for the critical temperature and strength of the phase transition. On the other hand, the gauge and renormalisation scale dependence may be significant, and often impact the existence of the phase transition altogether.


Introduction
A central motivation for the ambitious program of the next generation experiments that probe physics beyond the Standard Model (SM) is to understand aspects of electroweak symmetry breaking in the cosmological context [1][2][3]. If the electroweak phase transition (EWPT) was strongly first order, it would satisfy Sakharov's third condition for baryogenesis [4] and the collision of sound shells could create a stochastic gravitational wave background that is expected to peak at a frequency range visible to LISA [1]. In the SM, however, the EWPT is predicted to be a crossover [5][6][7][8][9] and new states near the electroweak scale would be needed to change it to a strong first-order transition . Either a 27 TeV upgrade to the LHC [48] or the ambitious 100 TeV proposal [49] could potentially confirm the SM prediction [5][6][7][8][9].
Phenomenological studies usually assess the impact of new states on the EWPT through the perturbative effective potential. As has long been recognised, the perturbative effective potential depends on a choice of gauge and gauge fixing parameter ξ [50], a choice of renormalisation scheme [51,52] and a choice of renormalisation scale [53]. Besides these choices, there are different ways to remedy infrared (IR) divergences caused by Goldstone boson loops and different treatments of daisy diagrams [54]. These subjective factors lead to theoretical uncertainties and a degree of arbitrariness associated with the perturbative EWPT predictions. For a recent update on non-perturbative approaches, see Ref. [55,56]. An effective field theory approach is also being developed that promises a gauge-independent evaluation of the effective potential free of some of the problems of the traditional approaches [55][56][57][58][59].
In this work we are predominantly motivated by the possibility of a first-order EWPT. In that case, the ground state makes a discontinuous jump while electroweak symmetry breaks. As is customary, we characterise a transition by a critical temperature, below which the transition is possible, and a transition strength. More detailed descriptions could include, for example, the latent heat, which characterises the energy available to phenomena such as gravitational waves, and the nucleation temperature. A critical temperature is defined as any temperature at which the effective potential contains two degenerate minima, i.e., where Φ and Θ are distinct locations in the field space where the potential has a minimum at the temperature T c . For non-manifestly gauge covariant calculations the strength of the transition is typically quantified by where the sum only includes components of the fields that break EW symmetry. Within electroweak baryogenesis, the condition that γ EW 1 approximately corresponds to avoiding washing-out baryon asymmetries by sphaleron processes [50]. The rate of sphaleron processes is governed by the spaleron action, which in the context of baryogenesis could be a more appropriate characterization of a phase transition than γ EW . The sphaleron action is gauge independent if computed carefully [60][61][62]. However, computing it requires solving a system of coupled differential equations, and hence it tends to be calculated for benchmarks rather than scans. Such computationally expensive calculations are not currently feasible for this study because we wish to explore the arbitrariness of perturbative calculations across the parameter space of a model. So we leave that for future work. We compute T c and γ EW numerically from the effective potential using PhaseTracer; see the manual [63] for an overview of the numerical methods used. For the purposes of the work described here we extended 1 PhaseTracer to include other approaches such as the so called expansion approach of Ref. [50]. These features should be available in a future release of PhaseTracer and are not to our knowledge implemented in any other public software.
In this work, for the first time we present a comprehensive study of uncertainties that arise in the simplest widely used techniques in theoretical calculations. 2 Such techniques include gauge-independent methods that lack resummation at leading order [50] and gauge-dependent methods that include leading-order resummation. Additionally, we examine numerical differences arising from different renormalisation schemes and gauge sensitivity arising from R ξ and covariant gauges and gauge-dependent tadpole constraints. Quantifying gauge dependence is difficult, as although there are indications from perturbativity, it is unclear how much to change the gauge parameter. Our investigation of renormalisation scheme and scale dependence, and the two gauge-independent approaches, on the other hand, do not depend on these choices. Approximate as they may be, such techniques are fast and simple and have been used in phenomenological analyses in many models, including large multi-parameter scans. In this paper we are agnostic about the theoretical properties, assumptions and shortcomings in the methods, and instead focus on their numerical differences and abitrariness, and leave it to the reader to decide the most appropriate method for their analysis. As our benchmark model, we consider the SM augmented by a real scalar singlet (the SSM) and a Z 2 discrete symmetry imposed on the potential. This model is ideal in its simplicity for highlighting key features of the parameter space as it includes both tree-level and thermal barriers depending on whether the discrete symmetry is broken at non-zero temperatures (but is enforced at zero temperature). It also includes cases where the gauge boson thermal loops are sufficiently important that gauge dependence might become crucial, as it does in the SM with a light Higgs boson [50]. We build on recent progress in examining the impact of theoretical uncertainties on the EWPT. The gauge dependence of tunneling rates and the implications for gravitational wave signals in gauge-dependent techniques were examined in Ref. [60,64,65] and the scale dependence in multiple models was investigated in Ref. [66][67][68]. Here, we consider the combined impacts of choices of scale, gauge and gauge fixing parameter, treatment of higher-order daisies, and the Goldstone catastrophe (GC).
The paper is organised as follows. In Sec. 2 we define our model and discuss various treatments of the effective potential, including choices of renormalisation scheme and gauge. In Sec. 3, we describe modifications to the effective potential to account for daisy diagrams and avoid the GC. In Sec. 4, we show our numerical results on the impact of the choice of renormalisation scheme, scale, resummation methods, and gauge. Finally, the conclusions and discussions are given in Sec. 5.

The effective potential
Evidently, analysing the dynamics of the phase transition requires an accurate treatment of the effective potential. In the following subsections, we discuss the methods considered in this paper. Unfortunately, all methods suffer from some of the following potential theoretical ambiguities: 1. Gauge dependence -the effective potential is fundamentally gauge dependent as it only contains one-particle irreducible loop diagrams. This doesn't include external leg corrections that would otherwise be included in scattering amplitudes to cancel the gauge dependence with other diagrams and maintain gauge independence [69]. The choice of gauge can strongly impact the EWPT predictions; for example, in the SM the strength of the EWPT at leading order is proportional to terms involving the gauge parameter, ξ 3/2 [50]. 3 2. Scale dependence -the one-loop potential and resulting phenomenological predictions could strongly depend on the scale Q and choice of renormalisation scheme. The dependence originates from regularizing the loop corrections in a particular scheme, as usual.
3. The IR problem [70][71][72] -the perturbative expansion may break down near the critical temperature. The loop expansion involves the mode occupation g → gn B ∼ gT /m for gauge bosons and λ → λT 2 /m 2 for scalars as shown in the daisy diagram in Fig. 1. 4 This expansion diverges when T m. Specifically the Goldstone contributions can get large deep in the broken phase and the gauge boson contributions can get large for small field values when m 0. There exists two prescriptions that alleviate the IR problem by resumming the most dangerous daisy diagrams to all orders; see Sec. 3.2.
4. Goldstone catastrophe (GC) [74][75][76] -Goldstone bosons could be massless, e.g., in the ξ = 0 gauge in the tree-level vacuum, and lead to divergences in the secondderivatives of the potential. Although terms in one-loop corrections to the potential of the form x 2 ln x → 0 vanish as x → 0, the second derivatives diverge, meaning that we could not compute e.g., the Higgs masses at one-loop. Whilst the photon is also massless, its mass isn't field-dependent, and doesn't lead to these problems. See Sec. 3.1 for a discussion of the treatments of this problem considered in this work. The occurrences of these problems in the treatments in this paper are summarised in Tab. 1. In addition, uncertainties on the SM nuisance parameters, such as the top mass, result in parametric uncertainties. The benchmark model (the SSM) used in our calculations is the SM extended by a real scalar singlet. The tree-level Higgs potential for the Higgs doublet and the real scalar field respects a Z 2 symmetry, S → −S, Note the convention of the negative sign in front of the µ 2 h and µ 2 s terms. Denoting the fields as where we treat φ h and φ s as real background fields, the tree level effective potential takes the form We impose constraints on the free parameters in the potential to fix the masses of the scalar particles and the vacuum. First, we require that at zero temperature an extremum occurs at φ h = v h 246 GeV and φ s = v s = 0, where v ≡ (v h , v s ) and the potential is to be evaluated at T = 0. We could impose this constraint on the tree-level or one-loop potential. The second condition in eq. (6) is satisfied trivially with a vanishing singlet VEV, while, at tree-level, the first condition for the Higgs requires, We only fix the zero-temperature vacuum in our model -at finite temperature the vacuum needn't satisfy eq. (6). Typically, at high-temperature the vacuum lies at the origin, φ = (0, 0). As the Universe cools, there is a smooth transition in the singlet direction, (0, 0) → (0, v s (T )). From that state, there is the possibility of a FOPT that breaks EW symmetry, (0, v s (T )) → (v h (T ), 0). Subsequently, the ground state smoothly evolves to the desired (v h , 0) state at zero temperature. This was demonstrated to be the only possible pathway for a FOPT in a high-temperature approximation [77]. For further discussion see appendix A. Second, we require that the vacuum reproduces the experimentally measured Higgs mass, m h = 125.25 GeV [78], and the desired scalar singlet mass, which we treat as an input in this work. If we neglect the momentum in the self energy (see remarks in Sec. 3.3 detailing the impact of this approximation) this corresponds to, To ensure the potential has the correct form for this we eliminate the λ h quartic by fixing m h and exchange µ 2 s for the singlet mass m s . This leaves only m s , λ s and λ hs as free parameters in our model.
The desired squared masses are positive, which ensures the zero-temperature VEVs correspond to minima of the potential. After fixing the Lagrangian parameters, we check that the global minimum of the potential at zero temperature is as desired. We may impose eq. (8) on the tree-level or one-loop potential. The tree-level scalar masses are where we use m h,s to emphasise that these are scheme-dependent tree-level masses. When solving this at higher order with the one-loop Higgs pole mass and one-loop tadpoles, we numerically solve for λ h including the one-loop tadpoles iteratively. We detail treatments of the one-loop corrections to this potential in the following subsections.

MS + R ξ : The MS scheme in the standard R ξ gauge
This is a standard treatment of the effective potential. The one-loop corrections to the potential in the standard R ξ [79] gauges 5 and MS scheme are given by [50,53] where Q is the renormalisation scale and m i are field-dependent MS masses (for their specific forms in our model see appendix C), and we have suppressed the scale dependence of the masses. The sum over h represents a sum over the scalar fields, which at the electroweak breaking minimum can be separated into physical Higgs and Goldstone bosons. The sum over V includes transverse and longitudinal massive gauge bosons. The third term is a gauge dependent sum over Faddeev-Popov (FP) ghosts that contribute to the effective potential in ξ = 0 gauges. There is one FP ghost corresponding to each vector boson with mass Lastly, the sum over f includes relevant fermion fields. The one-loop thermal corrections in the R ξ gauge are Note that in the tree-level minimum, the Goldstone boson masses are proportional to the corresponding gauge boson masses and ξ and thus equal the FP ghost masses in eq. (13), and that the degrees of freedom match, n G = n FP . As a result, the ξ-dependence from Goldstone and FP ghost contributions cancel in eqs. (12) and (14), as expected from Nielsen's identity. 5 Although we refer to this as the standard R ξ gauges, it is important to note that we do not use fixed VEVs in the gauge fixing Lagrangian, as is standard when calculating most observables. However since the field in the effective potential should be allowed to vary, the use of a fixed VEV would spoil cancellations of the mixing between the scalar (goldstones) and the gauge sectors, so instead the field is allowed to fluctuate away from the minimum. Although we still refer to this as the standard R ξ gauge here, it is sometimes referred to as the background field R ξ gauge in the literature, see e.g. [80], to distinguish it from the case with a fixed VEV.
We construct the total effective potential by adding the one-loop terms to the treelevel potential, Additionally, there are further possible modifications to the potential to treat the GC and resum daisy diagrams; see Sec. 3.

MS + cov: The MS scheme in the covariant gauge
The previous MS treatment considered the standard R ξ gauge. The standard R ξ gauge used for effective potentials includes scalar fields that we vary in the gauge fixing Lagrangian. However this means that a different gauge is used for calculating the potential at each field value, which may be a cause for some concern [54,81,82]. 6 To remedy this, the effective potential of the SSM using the covariant gauge (sometimes called the Fermi gauges in the literature) [81][82][83][84] in the MS scheme has been used in recent work [48]. The gauge fixing Lagrangian in this gauge is, This avoids the issue of a fixed VEV at the expense of mixing between the Goldstones and the gauge sector. We here use the same potential as derived in Ref. [48], but with a Z 2 discrete symmetry enforced, In the above, the field-dependent fermion, Higgs and gauge boson masses take their usual form. There are two gauge-fixing parameters, ξ W and ξ B , that appear only implicitly through the masses of the mixed Goldstone-ghost scalars, 7 where and the gauge dependence appears only in the terms Because the Υ terms are proportional to the derivative of the tree-level potential, the ξ-dependence vanishes in the tree-level minima, as expected from Nielsen's identity. If Υ > 0, by increasing ξ sufficiently it may be possible to make the square roots in eqs. (20) and (21) imaginary. As we take the real part, in this regime the potential would become independent of ξ.
The multiplicities for the Goldstone-like modes are (n 1 , n 2 ) = (2, 1) respectively to account for six terms in the effective potential, that is double the number of Goldstone modes. This is because the Goldstone and ghost contributions mix in the covariant gauge, so the total number of terms in the gauge-dependent effective potential is the same as in the R ξ gauge. The positive and negative roots in eqs. (20) and (21) may both produce zero eigenvalues, potentially resulting in GC. The zeros occur in the tree-level vacuum, where Υ = χ = 0 for all ξ. For the negative roots, however, when ξ = 0 the zero eigenvalues are field-independent and so cannot contribute to the GC.

HT: The high-temperature approximation
The high-temperature (HT) approximation of the total effective potential is where the Debye coefficients c h and c s are given in eqs. (47) and (48). At zero temperature, this reduces to the tree-level potential, V HT (φ, T = 0) = V 0 (φ). Equation (24) is obtained by considering only the leading terms of a high-temperature approximation of the finite-temperature correction in eq. (14). Namely, for y ≡ m/T 1, the thermal functions are approximated by the leading terms in the expansions [70] The zero-temperature Coleman-Weinberg corrections are neglected. Importantly, the gauge parameter ξ and renormalisation scale Q do not appear in this high-temperature expansion, leaving an explicitly gauge-and scale-independent high-temperature approximation to the potential. However, the fact that Q does not explicitly appear in the potential may in fact worsen the scale dependence, as explicit and implicit scale dependence are expected to partially cancel. The tadpole constraints are applied at tree-level, which avoids a potential source of gauge dependence. In this scheme, the potential contains only quadratic and quartic terms, and the critical temperature and transition strength can be solved analytically,

PRM: The expansion
This is a technique for computing the critical temperature and transition strength in a gauge-independent way [50]. The method exploits the Nielsen identity [85], which, when truncating the effective potential to one-loop, means that the one-loop effective potential is gauge independent at tree-level turning points. To maintain gauge independence, we ultimately use the tree-level, one-loop and high-temperature expansion of the effective potential in Sec. 2.3. As usual, the renormalisation scale enters explicitly in the Coleman-Weinberg potential and implicitly in the Lagrangian parameters. We first find the minima of the tree-level potential, φ tree . It is possible that the structure of minima in the tree-level potential does not match that of the one-loop potential. Therefore if there is only one minimum at tree-level, we further consider a saddle point. If there are no saddles, we consider a maximum. To ensure gauge independence of the critical temperature, the two points must be turning points. 8 To ensure that no gauge dependence creeps into the tree-level potential, we fix Lagrangian parameters only using tree-level tadpoles and tree-level constraints on masses. To ensure the validity of Nielsen's identity, we cannot modify the field-dependent masses in the Coleman-Weinberg potential, even by terms that are formally order O( 2 ) or resum daisy terms by either the Parwani or Arnold-Espinosa methods.
We then turn to the one-loop potential with no daisy resummation. We find the temperature at which the two points considered are degenerate in the one-loop potential. This is the critical temperature, T c . Lastly, at that critical temperature, we find the minima of the gauge-invariant high-temperature expansion of the potential, φ HT . We use this in our estimate of the transition strength. For the strength of a transition involving one field that breaks electroweak symmetry, In multi-field cases beyond our model, we would use the distance from φ HT to the h-axis. Importantly, because we use tree-level constraints and tadpoles to maintain gauge independence, we predict v 246.22 GeV and m h 125.25 GeV only at tree-level. The most common approach for calculating the masses and observables in BSM theories to fix v = 246.22 GeV and use the loop corrected effective potential to fix parameters in the potential such that the minimum lies at v = 246.22 GeV. In this way, one-loop tadpole diagrams vanish and do not need to be included in calculating observables. This is the approach implemented in most public tools, and therefore one needs to proceed with care when using the PRM method. The one-loop predictions for them could be quite different and disagree with measurements, and, depending on how they are computed, be gauge dependent. As we shall discuss, there is a risk that we are merely trading sensitivity to the choice of gauge for imprecision in agreement with experimental measurements; neither are desirable. Alternatives approaches, such as the Fleischer-Jegerlehner treatment of tadpoles [86] may provide an approach to ensure that λ is extracted from the Higgs mass at the order without re-introducing gauge dependence into to effective potential (see also Ref. [87] for a recent discussion about the differences in these approaches). Later we study the impact of extracting λ h at one loop by reintroducing gauge dependence via one-loop EWSB conditions and label this PRM + 1LHiggs 1LTad.
Lastly, we note that whilst the daisy treatments are an all-order resummation, and thus do not naturally fit into an -expansion scheme, Ref. [50] proposed a gauge-invariant modification to eq. (34) that could be incorporated into the PRM scheme. This method requires the effective three-dimensional theory and, as far as we know, is only rigorously proven to capture the O( 2 T 3 ) terms correctly. The conflict between IR divergences, which appear to require an all-order resummation, and the expansion is further discussed in Ref. [88].

OS-like: The on-shell-like scheme
The on-shell-like (OS-like) scheme is often used in studies of first-order phase transitions (see e.g. Ref. [39,66,[89][90][91][92][93][94][95][96][97]). This scheme has the practical advantage that the treelevel VEVs and mass eigenstates are preserved at the one-loop level. That is, the masses extracted from the one-loop potential are the same as those extracted from the treelevel potential, and both potentials have a minimum at the same location in field space. Therefore in this work we compare results from the MS scheme and this OS-like scheme to indicate the differences in results that may be expected between the two most commonly used schemes.
Uncertainties from zero-temperature higher-order corrections are often estimated by comparing results in different renormalisation schemes. While both schemes can provide controlled approximations allowing predictions that can be tested in experiments, the predictions will not in general be the same, with the difference coming from theoretical uncertainties arising from missing higher-order corrections in the respective schemes. Therefore, comparing the calculation in different schemes provides sensitivity to missing higher-order corrections and can be used as an estimate of this theoretical uncertainty.
Ref. [89] demonstrated that for field-dependent masses of the form we may pick counter-terms that satisfy the OS conditions and bring the zero-temperature one-loop corrections in the Landau gauge into the simple form where i ranges over the entire field-dependent mass spectrum (including scalars, gauge bosons and fermions), and we have definedñ i ≡ (−1) 2s i n i to absorb the sign for particle species of spin s i . The last term is field-independent and so is often ignored. The finitetemperature corrections are identical to those in Sec. 2.1. We make a common simplification: we ignore the limitation eq. (30) and use eq. (31) in our SSM model with the field-dependent masses in appendix C. The OS conditions remain satisfied; however, eqs. (12) and (31) differ by terms that do not appear in the treelevel Lagrangian and cannot be interpreted as counter-terms. However, we predominantly consider transitions between ground states in which mixing between the Higgs and singlet vanishes, that is, between (v h (T ), 0) and (0, v s (T )) ground states. This potential is correct along these axes as the field-dependent masses are in the form eq. (30). See e.g. Ref. [32] for further discussion about implementing an OS-like scheme in scalar extensions of the SM.

Modifications to the effective potential
The treatments in Sec. 2 may be modified to alleviate the Goldstone catastrophe and resum daisy diagrams.

Goldstone catastrophe
The Goldstone catastrophe (GC) occurs when a massless goldstone boson entering loop functions leads to an infrared divergence. The GC should appear directly in the MS ξ = 0 potential at three-loop order and higher and in the electroweak symmetry breaking conditions at two-loop order and higher [75]. Furthermore it is very common in BSM physics investigations for the effective potential to be calculated using the effective potential approximation or OS-like schemes, and as discussed in Sec. 2 this causes the second derivative of the one-loop potential to diverge. For example, in the MS case this happens when the tree-level masses that appear in loop functions of the Coleman-Weinberg potential are fixed using a tree-level EWSB condition 9 , such that for ξ = 0 we have massless Goldstone bosons. This particular issue could be avoided simply by relating the parameters to the Higgs pole masses with self energies evaluated with the momentum set equal to the Higgs mass [11] (see also Sec. 3.3 for comments on the impact of including the momentum), but here we consider instead two more direct alternatives that avoid the massless Goldstone and could be used to avoid the GC appearing directly in the potential and electroweak symmetry breaking conditions: 1. GC SelfEnergy Sol: Adding the tree-level Goldstone boson mass with a one-loop self-energy correction to all Goldstone masses [74,75] where the self energy contribution for the Goldstone can be given by, where V CW does not include Goldstone boson contributions. This correction depends on the renormalisation scale and gauge. 10 2. GC Tadpole Sol: Applying EWSB constraints eq. (6) at one-loop when calculating the masses that enter the one-loop corrections to the potential.
In both cases this effectively adds a correction to the potential which is formally beyond the one-loop precision we aim for, while shifting the Goldstone masses away from the treelevel value. The former is formally re-summing a sub-class of corrections to all-orders, providing a principled solution (please see the extensive discussion in Refs. [74,75] for details), while the latter is simply a convenient way to avoid the massless tree-level Goldstone that is easy to apply and done in the literature [66], while only introducing spurious two-loop pieces to the effective potential.

Daisy diagrams
There are two popular methods for resumming daisy diagrams: 1. Parwani [98]: This resums the finite-temperature daisy diagrams. To use the Parwani method, the mass eigenvalues appearing throughout the one-loop potential are replaced by thermal mass eigenvalues (see appendix C) and one forgoes the inclusion of any explicit daisy term. With this method it is important to only substitute in the thermal masses with lowest order from the high temperature expansion to avoid a introducing spurious linear terms [99,100].
2. Arnold-Espinosa [101]: This resums only the zeroth Matsubara modes in the daisy diagrams. We add the daisy term wherem 2 i (φ, T ) are thermal masses including Debye corrections (see appendix C), and we have suppressed their gauge and scale dependence.
Both methods introduce terms to the effective potential that are cubic in the coupling g for gauge bosons and to the 3/2 power in scalar couplings. The thermal mass is, by contrast, quadratic in g and linear in scalar couplings. Neither treatment changes the zero-temperature potential. See Ref. [102,103] for recent progresses about the two-loop order computations.

Renormalisation group equations
We take into account both the implicit and explicit scale dependence of the effective potential by using renormalisation group equations (RGEs) to run the MS parameters and compute the full MS renormalised effective potential. To do this we implement a FlexibleSUSY 2.6.1 [104,105] 11 version of the effective potential in PhaseTracer. The FlexibleSUSY spectrum generator has two-loop RGEs, one-loop threshold corrections for the extraction of the gauge and Yukawa couplings and it also extracts the Higgs potential parameters using one-loop electroweak symmetry breaking conditions and oneloop self energies. The latter include the p 2 contributions, which are neglected when masses are calculated using the one-loop effective potential approximation. We have checked that it will lead to 1.3% difference on extracted λ h and 1.0% difference on output

Results
We now investigate different sources of uncertainty in the analysis of a PT, including gauge dependence, renormalisation scale dependence, treatment of the Goldstone catastrophe (GC), resummation of daisy diagrams, and renormalisation schemes. We do so by applying the treatments discussed in Secs. 2 and 3 to the SSM. After considering them separately, we compare their typical sizes across the available parameter space. The numerical uncertainties associated with our computational methods are discussed in appendix D.

Benchmark point
First, in Fig. 2 we show the ξ-dependence of the critical temperature T c and transition strength γ EW in the R ξ and covariant gauges for different methods, using a benchmark This benchmark lies near the center of the region of parameter space in which a FOPT could occur. We vary the gauge parameter between ξ = 0 and 60. Note that the upper limit may lie beyond constraints from perturbativity [60,81], as increasing the gauge parameter aggravates IR divergences that spoil perturbation theory. The legend follows the notation introduced in Secs. 2 and 3. In the three MS methods, we used the Arnold-Espinosa technique to resum daisy terms, and set Q = m t . The gauge dependence of the OS-like scheme is not presented as eq. (31) is only valid at ξ = 0. As a reference, in Fig. 2 we show the PRM and HT methods which, of course, are independent of the gauge parameter ξ. The left panels of Fig. 2 show results for the R ξ gauge. In the MS scheme, with ξ increasing, the critical temperature T c increases slightly for small values of ξ, and then decreases significantly for large values of ξ. With ξ increasing even further, there is no FOPT for this benchmark point anymore in the gauge-dependent schemes. The choice of treatment of the GC has limited impact.
In Fig. 2 we investigated ξ ≤ 60. However, there may be constraints on perturbativity that impose ξ 1 [60,81]. We see in Fig. 3 that varying the gauge between ξ = 0 and 1 results in only slight changes in the properties of the FOPT. We furthermore see that some methods do not depend smoothly on ξ. For example, the T c of MS blows up at ξ = 0. This originates from the GC: due to the catastrophe, the second-derivatives of the oneloop potential are singular, and so the Lagrangian parameters required in the tree-level potential to maintain finite scalar masses diverge. Although using the modified µ from one-loop tadpole cures the GC at ξ = 0, the Goldstone may be massless at tree-level in In contrast, adding the one-loop self-energy to the Goldstone masses avoids a catastrophe at any ξ as seen in the blue curve for MS + GC SelfEnergy Sol.
In PRM + 1LHiggs 1LTad, we use one-loop tadpole conditions to set the Lagrangian parameter µ h and extract λ h and µ s at the one-loop level from the measured Higgs pole mass and the input singlet pole mass, respectively. The one-loop tadpoles are required if one wants to extract λ h from the Higgs mass at the one-loop level using the standard approach we applied in the MS calculation, where the tadpoles vanish in the mass corrections and do not need to be explicitly included 12 . A one-loop extraction of λ h is required for a consistent order calculation of the effective potential consistent with the measured Higgs pole mass. As a result, the T c and γ EW calculations in this method are closer in spirit to those in the MS method and numerically are closer when far away from the ξ values where the MS and pure PRM results match. This, however, introduces gauge dependence, which contradicts the purpose of the PRM method. This gauge dependence originates from the gauge dependence of the MS scalar masses (which are extracted at the one-loop level) used in the one-loop tadpole conditions. Just as in the MS scheme, however, T c significantly changes with increasing ξ for this case at this benchmark point.
We show results in the covariant gauge in the right panels of Figs. 2 and 3, where ξ ≡ ξ B = ξ W . As anticipated, results from the covariant and R ξ gauges are identical at ξ = 0. The results for the MS calculation without any modification to solve the GC cannot be seen as the GC exists for all ξ in the covariant gauge.
Different to the R ξ gauge, in the covariant gauge the T c for MS + GC SelfEnergy Sol and + GC Tadpole Sol respectively increase and decrease almost linearly with ξ increasing. Meanwhile, there are some discontinuities in the lines of the covariant gauge when ξ is small, while the lines then become continuous when ξ is large. This is because we take the real part of the square root in eqs. (20) and (21) and for sufficiently large ξ, the ξ-dependent parts are imaginary. This behavior can also be seen in the results of PRM + 1LHiggs 1LTad in the covariant gauge, i.e. there are two steps in the results at small ξ and the results are independent of ξ for large ξ. This gauge independence, however, does not necessarily mean that this method is a better choice, because it is simply due to the fact that we ignored the imaginary gauge-dependent parts of the squared masses. There are no results when ξ > 46.5 for MS + GC Tadpole Sol because our iterative solver for λ h using one-loop EWSB did not converge. The gauge dependence of the MS + GC SelfEnergy Sol in the covariant gauge is also fairly strong, though it shows a different tendency.
The gauge dependence in the calculations of T c and γ EW originates from two sources: indirectly from the Lagrangian parameters chosen by solving one-loop tadpole conditions and explicitly in the one-loop potential. The indirect dependence causes T c to decrease with increasing ξ, as shown by the PRM + 1LHiggs 1LTad results in R ξ gauge, which only suffers from this source. The explicit dependence, on the other hand, increases T c , as shown by the MS + GC SelfEnergy Sol results in the covariant gauge with large ξ. In this case, the ξ dependence vanishes in the zero-temperature vacuum, as the square roots in eqs. (20) and (21) are imaginary, such that the tadpoles are gauge independent. However, the potential remains gauge dependent at other field values. When both contributions appear at the same time, such as in the MS scheme with the R ξ gauge and the MS + GC Tadpole Sol in the covariant gauge, the indirect dependence through tadpoles overpowers the explicit dependence.

Two-dimensional scans
The above discussion is based on a single benchmark point. The magnitude of the gauge dependence, however, may depend strongly on the Lagrangian parameters. In Fig. 4, we show the uncertainty from the choice of gauge parameter in the MS + GC SelfEnergy Sol prediction for the critical temperature T c with three two-dimensional scans of the Lagrangian parameters. We quantify this uncertainty through where T c (ξ = 0, 25) is the critical temperature with ξ = 0, 25 in the MS + GC SelfEnergy Sol scheme. As discussed in the introduction, this quantification is somewhat arbitrary.
We attempted to pick an upper limit, ξ = 25, that strikes a balance between understating the uncertainty and constraints from perturbativity. If we were to increase the upper limit further, eventually there would be no points with FOPT simultaneously at ξ = 0 and the upper limit. We scan over two of the three free parameters, while fixing m s , λ hs and λ s to their values in eq. (35) in the left, middle and right panels, respectively. To help understand the gauge dependence further, we show T c (ξ = 0) in the lower panels. We found no valid points in the blank regions of Fig. 4. The causes for this and the phase structures of points in those regions are discussed in appendix A. In the colored regions, the critical temperature changes with the varying input parameters. The trends can be derived from the expression of T c in the high-temperature approximation, eq. (27). Alternatively, consider the two minima at temperature . We may Taylor expand the difference in potential around T = 0 [23,112] where Since this difference vanishes at the critical temperature and ∆V (T = 0) < 0, we find Accordingly, a smaller gap at zero temperature may imply a lower T c . In the SSM at tree level we have where 1 2 λ hs v 2 h > m 2 s . We can approximately infer that with increasing λ s and m s , and decreasing λ hs , the gap ∆V (T = 0) and thus T c must increase. This agrees with the lower panels of Fig. 4.
In the upper panels of Fig. 4, we show the uncertainty in T c caused by changing the gauge parameter ξ in the regions in which a FOPT could occur for both ξ = 0 and ξ = 25. The change ∆ ξ T c varies from −73 GeV to +55 GeV, and increases smoothly with λ hs decreasing, m s increasing and λ s increasing. In Appendix B we also provide Tab. 2 which shows precise input parameters and the output critical temperatures of the benchmark point in eq. (35) and points around it for ξ = 0 and ξ = 25, to give quantitative results that also show the trends discussed here.

Renormalisation scale dependence
The renormalisation scale dependence of the PT properties originates from the scaledependence of the one-loop corrections to the effective potential. If RGEs are used, this may be partly balanced by renormalisation-scale dependence of the Lagrangian parameters. In this work, we run the Lagrangian parameters to different scales using Flexible-SUSY, which uses the MS scheme with ξ = 1 by default. We do not consider the OS-like scheme in this section, since there is no free renormalisation scale to vary for that method.

Benchmark point
In Fig. 5, we show the renormalisation scale dependence in the MS scheme and PRM methods, by varying the renormalisation scale from Q = 1 2 m t to Q = 2m t . We choose this interval because it is expected to capture uncertainty from the most important missing high-order corrections from the top quark at T = 0. Here we use the one-loop resummed Goldstone mass to treat the GC. In each of the panels, we vary one input parameter around the benchmark point in eq. (35). The top and bottom two rows of panels in Fig. 5 show the uncertainties on T c and γ EW , respectively. The lower edge, center, and upper edge of the bands in the first and third rows indicate results at Q = 1 2 m t , m t , and 2m t , respectively. To highlight the uncertainties, we separately show the fractional changes, (T c (Q = 1 2 m t ) − T c (Q = 2m t ))/T c (Q = m t ) and (γ(Q = 1 2 m t ) − γ(Q = 2m t ))/γ(Q = m t ), directly by colored lines in the second and fourth rows of panels.
We see almost no parameter dependence in the results at λ hs 0.2 and m s 90 GeV. In these regions the transition is between (φ h = 0, φ s = 0) and (φ h = 0, φ s = 0) and thus is independent of parameters related to the singlet. Unlike Ref. [113], we show PRM results even if the high-temperature approximation of the effective potential contained a single minimum at the critical temperature. For λ s ≈ 0.15, that single minima was at the origin for Q = 1 2 m t leading to γ EW = 0. The change in T c for the PRM methods is typically mild and less than about 10%, though γ EW shows more sensitivity, typically less than about 20%.
A common and important feature in all the results, however, is that the uncertainties rise sharply when the critical temperature is sensitive to the Lagrangian parameters. This is not only true for the uncertainty caused by scale dependence, but also true for the uncertainty caused by gauge dependence discussed in the above subsection, and the uncertainties presented in the following subsection. As the variation of renormalisation scale or gauge parameter slightly changes the Lagrangian parameters, it is obvious that when the critical temperature is sensitive to the parameters, the critical temperature and transition strength become sensitive to the renormalisation scale or gauge parameter. This happens in the region of large λ hs , small λ s and small m s , where the critical temperature is relatively small. There are no results for T c < 30 GeV in Fig. 5 because we find points in the plots by a grid scan of the parameters, and the grid misses fine-tuned cases in which T c < 30 GeV. To find results in that region we performing a more detailed scan in a smaller parameter range where we expect to find smaller T c . We show the MS result for the resulting low temperature range in Fig. 6. To obtain T c < 5 GeV, the parameters have to be tuned to one part in 10 7 . The reason can be seen from eq. (39). A zero critical temperature means no energy gap between the EW VEV in the φ h direction and a local minimum on φ s direction, so the parameters should cancel each other.
With RGE running and the Parwani method, the scale uncertainties in the MS result for T c are about 16 GeV for T c ≈ 150 GeV. As T c decreases, the relative uncertainty changes from approximately constant to sharply increasing in magnitude at T c < 40 GeV. Using the Arnold-Espinosa method or eliminating daisy resummation results in similar uncertainties with the same tendency. When T c 20 GeV, the daisy corrections, which are proportional to the temperature T 4 (see eq. (34)), are negligible. Consequently, in this regime the scale dependence of different treatments of daisy resummation are approximately the same.
We see that ∆ Q T c of the MS scheme without RGE running is smaller than the one with RGE running, which is counter-intuitive. For example at leading order there are cancellations such that 13 [68] µ d dµ and there can be further cancellations of scale dependence in the temperature-dependent parts of the potential. In the Parwani method, scale dependence in finite-temperature bosonic contributions cancels scale dependence in temperature-dependent terms introduced in the CW potential by the Parwani method. The cancellation eq. (40) is, however, spoilt by next-to-leading order terms, such as two-loop running of parameters in the tree-level potential and one-loop running of parameters in the CW potential, and in general requires that we consider the β-functions for all Lagrangian parameters and the anomalous dimensions of the fields. In our calculation, we apply two-loop RGEs and omit the anomalous dimension, such that eq. (40) no longer holds exactly. Furthermore, we include finite-temperature corrections from the top quark, and the scale cancellations in the Parwani method only occur for bosonic corrections. Consequently, results that include RGE running suffer from greater scale uncertainties than results without RGE running, and the Parwani method makes the situation even worse. This is somewhat counter-intuitive, as the derivative of the Parwani improved Coleman Weinberg potential with respect to the RG scale has terms that are quadratic in temperature and negative in sign, so one might expect they naively cancel the scale dependence of the quadratic contribution to the thermal functions. However, the combined β function for the thermal mass turns out to be negative as it is dominated by top quark interactions. So the new scale-dependent contribution in Parwani adds to, rather than cancels, the scale dependence of the un-resummed theory. In Fig. 7, we display the zero-temperature potentials of the MS scheme with the Arnold-Espinosa method for the benchmark point of λ s = 0.1, λ hs = 0.3 and m s = 65 GeV. We can see clearly that eq. (40) does not hold with RGE running in our calculation, and that the potential suffers from greater scale dependence when we include RGE running. We also mark the tree-level minimum since the PRM method evaluates the one-loop potential at the minima of the tree-level potential. Fig. 7 shows that changes in the one-loop potential at the tree-level minima with RGE running are smaller than that without RGE running, as with RGE running the tree-level minima approximately lie on a horizontal line of constant potential. Therefore, the ∆ Q T c for the PRM method with RGE running is smaller than ∆ Q T c without RGE running. Note that for the PRM method, T c (Q = 1 2 m t ) > T c (Q = 2m t ). Fig. 5 shows that ∆ Q T c of the PRM method with RGE running is rather small and roughly independent of Lagrangian parameters, typically about 5% for T c , though it can exceed about 20% for very small critical temperatures. Without RGE running, ∆ Q T c can be larger than 15%, and increases rapidly with T c decreasing. This implies that proper RGE running is quite essential in the PRM method to reduce the renormalisation scale dependence. Looking only at ∆ Q T c , the MS scheme has larger scale uncertainties than the PRM method, but considering ∆ Q γ, the MS scheme has smaller scale uncertainties than the PRM method. The increased uncertainty in γ in the PRM method may be connected to the fact that the numerator, ∆φ h , and denominator, T c , are obtained from two different potentials in this method.
Finally, we consider a broader variation of the renormalisation scale from Q = T /2 to 2πT . This interval was chosen to capture uncertainty from missing higher-order corrections at finite temperature, specifically involving the soft scale ∼ gT and the first Matsubara mode as discussed in Ref. [67]. Here, the renormalisation scale changes throughout a calculation, as whenever the potential is evaluated at a temperature T , the renormalisation scale is set to a scale Q ∝ T . We perform calculations at Q = T , Q = T /2 and Q = 2πT by solving the following equations iteratively for T , By doing so, we ensure that at temperature T = T c the potential was computed at the renormalisation scales Q = T /2, T and 2πT , as required. We show the variation in the critical temperature in Fig. 8. We see that the uncertainty bands are broader compared to Fig. 5, especially for the PRM method for which the relative uncertainty was always

Two-dimensional scans
We now turn to two-dimensional planes of the SSM parameters to check the renormalisation scale dependence throughout the model's parameter space. Here we change the scale using RGE running, use the MS scheme and the R ξ gauge, and include daisy corrections. In Fig. 9, we show the difference of critical temperature ∆ Q T c = |T c (Q = 1 2 m t ) − T c (Q = 2m t )| on the λ s -λ hs , m s -λ s and m s -λ hs planes. The colored regions are obtained by two-dimensional scans with m s = 62.5 GeV fixed in the left panel, λ hs = 0.3 fixed in the middle panel, and λ s = 0.1 fixed in the right panel. In the blank regions there were no points with a FOPT at both Q = 1 2 m t and Q = 2m t . The tendency exhibited in the two-dimensional planes in Fig. 9 is consistent with that in the one-dimensional results. The uncertainties at first decrease with λ hs increasing, λ s decreasing and m s decreasing, but subsequently rise sharply. In most of the parameter space, the uncertainties lie around 8 GeV -16 GeV. The magnitude of the uncertainties can reach more than 50 GeV when T c < 60 GeV (see Fig. 6). The maximum value of ∆T c is dependent on the resolution of the scan. With two-dimensional grid scans of 200 points per dimension in the ranges 0.1 < λ hs < 0.5, 0.01 < λ hs < 0.3 and 10 GeV < m s < 120 GeV, we found a maximum ∆T c of 40.

Daisy resummation dependence
In Fig. 10, we present differences in T c and γ EW with different treatments of daisy diagrams -the Parwani method and the Arnold-Espinosa method described in Sec. 3.2, and discarding resummation of daisy diagrams. We can see that in the SSM the effects of daisy resummation on transition properties are relatively small (though this effect may be bigger in other models; see e.g. Ref. [11,32,114]).
The changes in T c caused by switching between the Parwani and the Arnold-Espinosa methods are less than 3% and the changes caused by switching off daisy resummation altogether were less than 6%. The changes decrease with λ hs increasing or λ s and m s decreasing, which means T c decreasing (see top panels of Fig. 11). Unlike changes caused by varying the gauge parameter or renormalisation scale, varying the treatment for daisy resummation does not alter the Lagrangian parameters. Therefore, here the changes in T c do not rise sharply at low T c . On the other hand, the relative changes in γ EW are typically even smaller, but increase rapidly for large T c , simply because that γ EW are very small. The absolute changes in γ EW are less than 0.06 when T c > 100 GeV and γ EW < 0.4.

Comparison of renormalisation schemes and other methods
We now turn to the renormalisation schemes. In Fig. 11, we compare results of the MS and OS-like scheme, together with the high temperature approximation which is independent of renormalisation scheme, and the PRM method which depends on the MS scheme. In both the MS scheme and the OS-like scheme, we use ξ = 0, the Arnold-Espinosa method for daisy resummation, and resum the Goldstone boson mass to solve the GC. The renormalisation scale in the MS scheme is set to m t . For the MS potential used in the PRM method we adopt the same settings as in the MS scheme, except that we do not resum daisies.
Since the zero-temperature one-loop corrections depend on the choice of renormalisation scheme, the tadpole conditions result in different Lagrangian parameters, λ h , µ h and µ s , in each scheme. Therefore, every part of the effective potential can be different in the MS and the OS-like schemes. However, the differences in T c and γ EW are typically mild with the OS and MS schemes giving surprisingly similar results for the choice of renormalisation scale we make for the MS scheme (Q = m t ). It was also found in previous work that the OS-like result lies within the scale variation uncertainty [66]. This suggests that the differences to the OS-like scheme is not really accessing any higher-order corrections beyond the MS scale variation though the very close agreement we see must be coincidental in nature.
Except for the flat region λ hs < 0.17 in the left panels and m s > 88 GeV in the right panels of Fig. 11, T c in the MS scheme is larger than T c in the OS-like scheme, and vice-versa for the transition strengths in the bottom panels. With critical temperature decreasing, the difference ∆T c = T c (MS)−T c (OS-like) increases monotonically from about 1 GeV to more than 10 GeV. This occurs for the whole parameter space, which can be seen from the results of the two-dimensional scans in Fig. 12.
The T c of the high-temperature approximation and the PRM method are far lower than that of the OS-like scheme and the MS scheme. There is no result for small λ hs , large λ s and m s , because the expression of T c of HT, eq. (27), assumes a minimum of the form (0, v high s ) exists at T c . The differences between results from the HT approximation and the OS-like scheme are easy to understand. They both use tree-level tadpole conditions to solve the Lagrangian parameters, λ h , µ 2 h and µ 2 s , which results in exactly the same tree-level potential. The one-loop finite-temperature corrections are also supposed to be similar. Since the zero-temperature one-loop correction in the OS-like scheme is zero at the EW VEV, the only obvious difference is the zero-temperature one-loop correction of the OS-like scheme around the minimum of φ h = 0. In this minimum, the fermions and vectors are massless, and the only non-vanishing contribution to the zero-temperature one-loop correction comes from scalars and is always positive. Therefore, to compensate this increase in the potential, the OS-like scheme needs larger thermal corrections, i.e. higher temperature than in the HT approximation, to achieve degenerate minima. As for the PRM method, it also uses tree-level tadpole conditions, and the zero-temperature one-loop corrections in the MS scheme reduce the vacuum energy gap, so it leads to lower T c than the HT approximation. However, in the HT approximation, both the vacua tend towards the origin as temperature increases, as can be seen in appendix A, which also reduces the vacuum energy gap, while the vacua in the PRM method are fixed. These two impacts cancel each other to some extent, leading to similar results in these two schemes. In contrast the results can also diverge more at higher critical temperatures. In that case the HT expansion is, as you would expect, becoming a better approximation of the MS calculation and gets closer to the MS result.

Comparison in whole parameter space
Here we compare all the possible changes investigated above in the whole parameter space. As shown in Fig. 13  With these choices, about 20% of points had a FOPT at ξ = 0 and Q = m t . Of those points, about 50% didn't have a FOPT at ξ = 25 and about 10% didn't have a FOPT at either Q = 1 2 m t or 2m t (or both). In about 80% of cases gauge dependence was the greatest source of uncertainty. The top panels show the samples on planes of Lagrangian parameters with points with the smallest scale uncertainty shown on top. There are no strong trends on these planes, though samples of small uncertainty tend to lie at large m s . In the bottom panels, the uncertainties are shown versus T c and γ EW , colored by the source of the maximal uncertainty. When T c < 100 GeV, varying the gauge parameter between 0 -25 always results in a greater change in the critical temperature than varying the renormalisation scale between 1 2 m t -2m t . When T c > 100 GeV, on the other hand, varying the renormalisation scale can have the greatest impact, and the lower bound of maximal uncertainty on T c is given by points for which the maximum uncertainty comes from changing the renormalisation scale. The relative maximal uncertainty on T c derived from gauge dependence decreases with T c increasing, while for scale dependence the relative maximal uncertainty always lies around 10%. Conversely, varying the renormalisation scale can have the greatest impact on γ EW for small γ EW , and max |∆γ EW |/γ EW derived from gauge dependence increase with γ EW increasing. The relative maximal uncertainty on γ EW can scarcely reach about 500%, because large γ EW corresponds small T c , which become sensitive to the Lagrangian parameters. The relative uncertainties can be as large as 100% when one of the T c 0 GeV, but it is rare to encounter such cases in general scans because the parameters need to be finetuned to achieve low T c . In Fig. 15, we display the distribution of relative and absolute uncertainties of T c and γ of the 3-dimensional random scan. We find that absolute (relative) uncertainties of T c and γ for most of the samples are smaller than 13 GeV (10%) and 0.4 (14%), respectively, but that no samples have uncertainties smaller than 10 GeV (9%) and 0.19 (11%). In extreme and rare cases, max |∆T c | and max |∆γ EW | reached 77 GeV and 0.49 respectively, and in an even more detailed scan the maximum uncertainty would be even larger. We see from the colors of the bars that for points with moderate or large uncertainties, the greatest source of uncertainty is almost always gauge dependence. This final remark is subject to the choice of scales and gauge parameters that we compare. With a smaller range of gauge parameters, ξ = 0 to ξ = 3, it was instead found that scale dependence was the greatest source of uncertainty for the majority of the parameter space.

Conclusions
First-order phase transitions could have occurred in the history of our Universe, helped lead to the abundance of matter over anti-matter, and left behind observable gravitational wave signatures. In perturbative calculations in models of particle physics, however, the properties and occurrence of first-order phase transitions are sensitive to somewhat arbitrary choices, including choices of renormalisation scale and gauge. We investigated the impact of those choices in detail, using the electroweak phase transition in the SM extended by a real scalar singlet as a benchmark. We explored the three-dimensional parameter space of this model through one-and two-dimensional slices as well as threedimensional scans.
We refrain from making recommendations about the best choices, and instead wish to emphasise the sizes of the uncertainties and which choices influence the results the most. The scale dependence in the MS scheme and gauge dependence were the most significant sources of arbitrariness. Regarding the former, we found that • Changing Q from 1 2 m t to 2m t induced moderate changes in the PT properties, typically less than about 10%.
• Scale dependence was particularly severe in the PRM scheme, and only somewhat alleviated by using RGE running.
• Counter-intuitively, applying two-loop RGE running in the MS scheme worsened scale dependence, possibly because this was not a strict fixed-order calculation and because we neglected the anomalous dimension of the fields.
• Similarly, the Parwani method worsened scale dependence, despite the possibility of cancellations in the scale dependence between the CW and finite-temperature parts of the potential.
• Whenever the PT properties strongly depend on the choices of Lagrangian parameters, they inevitably strongly depend on the choice of scale.
• Only about 10% of our points with a FOPT at Q = m t didn't have a FOPT at both Q = 1 2 m t and 2m t . In an OS-like scheme, on the other hand, we found quite similar results to the MS scheme at Q = m t , though in extreme cases the differences in critical temperature reached about 10 GeV. Regarding gauge dependence, we found that • In the MS scheme the gauge dependence of the critical temperature and transition strength could be O(100%) in extreme cases when ξ was varied from 0 to 25 for both the R ξ and covariant gauges, though was typically milder.
• When varying between ξ = 0 and 1, with an upper bound consistent with an O(1) bound from perturbativity [60,81], the gauge dependence was extremely mild.
• The gauge-independent PRM and HT predictions for PT properties were similar to each other and quite different from those from the MS method.
• Despite being gauge independent, gauge dependence creeps into the PRM and HT methods when Lagrangian parameters are determined from tadpole conditions or mass constraints in a gauge-dependent way. In these cases the gauge dependence was similar to that in the MS scheme.
• In the R ξ gauge, we can avoid the GC at ξ = 0 by using the Higgs mass squared parameter obtained from one-loop tadpole conditions in the calculation of the fielddependent Goldstone masses. However, even with this approach, Goldstone catastrophes still occurred for other choices of ξ for which Goldstones were massless in the EW vacuum.
Finally, the impact of the different treatments of daisy diagrams and the GC were relatively small, always changing the critical temperature by less than about 10 GeV. In summary, we investigated the major sources of arbitrariness in perturbative treatments of FOPTs, including renormalisation scales and schemes, and gauges and gauge fixing parameters. We found that moderate variations in these choices, especially the choice of gauge parameter, may lead to significant changes in the predictions for the FOPTs.
In the SSM, the FOPT almost always happens between a minimum of the form φ s ≡ (v high  Fig. 16. Fig. 16 shows the evolution of minima of the effective potential with temperature for different cases. In all the cases, there is only one minimum at high temperature, T > 400 GeV. In the top left panel, a second-order phase transition occurs at T 300 GeV from the origin to a minimum at (φ h = 0, φ s = 0), indicated by orange lines. With temperature decreasing, a FOPT becomes possible below T c = 114 GeV, after which the minimum at (φ h = 0, φ s = 0) becomes the deepest minimum of the effective potential. There are parameter points for which no FOPT occurs during the evolution of the potential with temperature. They may be split into three cases: • The minimum at (φ h = 0, φ s = 0) is never the deepest minimum, so there is no transition from the minimum at (φ h = 0, φ s = 0) to it, as shown in the top right panel of Fig. 16. This happens for large λ hs , small λ s and m S , i.e. the top left corners of the λ s -λ hs and m s -λ hs planes in Fig. 4.
• There is never a minimum at (φ h = 0, φ s = 0) during the evolution, and the phase transition from the origin to (φ h = 0, φ s = 0) is not a strong first-order transition. This is similar to the situation in the SM. This case is displayed in the bottom left panel of Fig. 16, and corresponds to small λ hs and large λ s and M S . This kind of transition is ignored in 2-dimensional and 3-dimensional scans.
• The phase transition from the minimum at (φ h = 0, φ s = 0) to the minimum at (φ h = 0, φ s = 0) is not first-order, which is presented in the bottom right panel of Fig. 16. It is caused by large λ s and small m s .
In our calculations, there are more complicated phase structures than above cases, caused by numerical problems, such as the fact that the thermal functions are not smooth.
To maintain consistency, unless otherwise stated we only compare and show the type of FOPT shown in the top left panel of Fig. 16 Table 2: Values of Lagrangian parameters, critical temperature, γ EW , and VEVs for points around the benchmark point in the MS + GC SelfEnergy Sol scheme. All dimensionful quantities are in GeV or GeV 2 . The blue numbers show difference between quantities evaluated at ξ = 0 and 25.
In Tab. 2, we present input parameters, parameters extracted from tadpole conditions, critical temperature and transition strength of the benchmark point used in Sec. 4, for different gauge parameters. Here we used the MS scheme in R ξ gauges, the Arnold-Espinosa technique to resum daisy terms, GC SelfEnergy Sol to solve the GC and set Q = m t .
longitudinal component of the W boson receives the Debye correction 11 transition using the HT method and the SSM as an example. The differences between the analytic formulae and the numerical results from PhaseTracer are less than about 0.01%. In this section, we further discuss numerical uncertainties and issues we encountered in this work.
The SM parameters involved in the effective potential models include mass of the SM-like Higgs m h , SM EW VEV at zero-temperature v, the gauge coupling constants g and g , and the Yukawa coupling y t , y b and y τ . Since we use FlexibleSUSY for our RGE running, we use the default values in FlexibleSUSY in all the models to be consistent, m h = 125.25 GeV, v = 247.4554 GeV, g = 0.6477, g = 0.3586, y t = 0.9341, y b = 0.01547, y τ = 0.010014.
Modifying these values will of course affect the values of T c ad γ EW . However, by regenerating all the results with m h = 126.0 GeV and v = 246 GeV we have checked that it will not qualitatively change our above discussions and conclusions.
There is no closed-form solution for calculation of the thermal functions, J B (y 2 ) and J F (y 2 ), in eq. (25). The precision of several numerical methods to do this have been fully discussed in [116]. PhaseTracer adopts the same method as CosmoTransitions [117], calculating a series of points by quadrature integration in SciPy library for Python and then using B-spline interpolation to get thermal functions of any point. This method is very fast and pretty accurate except around y 2 = 0. We find that the second-order derivative of the thermal function obtained in this way is not continuous at y 2 = 0. Fortunately, it is very rare that the minima at T c encounter y 2 = 0, so the T c of the FOPT is not ruined.
Nevertheless, the discontinuity of the second-order derivatives of the thermal function around zero mass may introduce an additional artificial minimum in the effective potential, which will form an artificial phase. In some cases, there is a transition between this artificial phase and an ordinary phase, especially if the Arnold-Espinosa method is used for daisy resummation. The Parwani method can partially avoid this situation, because the Debye corrections to the field-dependent masses are positive such that we do not encounter zero masses at the minimum. In any case, the artificial phase cannot affect our results because we focus on the FOPT described in appendix A, and ignore other transitions.
As already discussed, the number density of parameter points sampled in the scan will impact the maximum change in the critical temperature and transition strength found. Since the maximum change occurs when T c is very sensitive to the input parameters, the more detailed scan we have, the larger the maximum uncertainty we are likely to obtain. Besides, the local minimum finding algorithm adopted in PhaseTracer has large numerical uncertainty when the potential around the minimum tends to be flat, and has a chance of misidentifying a saddle point as a minimum. This causes noise in our plots. The more detailed scan we have, the more noisy points will appear.
These numerical problems do not only exist in the SSM or only when using Phase-Tracer. As the SSM is quite a simple model, we anticipate that further issues and larger uncertainties could be found in more complicated models. Therefore, they should be carefully considered in any perturbative calculations of phase transitions.