PhaseTracer: tracing cosmological phases and calculating transition properties

We present a C++ software package called PhaseTracer for mapping out cosmological phases, and potential transitions between them, for Standard Model extensions with any number of scalar fields. PhaseTracer traces the minima of effective potential as the temperature changes, and then calculates the critical temperatures, at which the minima are degenerate. PhaseTracer is constructed with modularity, flexibility and practicality in mind. It is fast and stable, and can receive potentials provided by other packages such as FlexibleSUSY. PhaseTracer can be useful analysing cosmological phase transitions which played an important role in the very early evolution of the Universe. If they were first order they could generate detectable gravitational waves and/or trigger electroweak baryogenesis to generate the observed matter anti-matter asymmetry of the Universe. The code can be obtained from https://github.com/PhaseTracer/PhaseTracer.


Introduction
Since some fundamental symmetries are expected to break in the early Universe the accompanying phase transitions are critical to our understanding of the phenomenon (see e.g. Ref. [1] for a review). These transitions occur when the system of scalar fields transits between two distinct minima, also called vacua, of the effective potential. If the vacua are separated by a barrier the transition is first order, and if the deeper minimum is charged a symmetry breaks spontaneously.
Finite temperature corrections to effective potentials result in important modifications to the free energy typically restoring spontaneously broken symmetries at high temperatures [2,3]. Consequently, if a symmetry is broken at zero temperature, a phase transition probably occurred as the Universe cooled. We therefore expect that there was an electroweak phase transition associated with electroweak symmetry breaking. There may have been other phase transitions in the early Universe, such as one that breaks symmetries of a Grand Unified Theory (GUT) that embeds the Standard Model (SM) gauge groups into a unified gauge group [4][5][6][7][8].
Alternatively, other gauge groups may have broken at inter-mediate scales, for example, extra U (1) gauge groups that are fairly generic predictions of string theory [9][10][11][12][13][14] and can solve the μ-problem of the Minimal Supersymmetric Standard Model (MSSM) [15][16][17]. These phase transitions played an important role in the evolution of the Universe and it is vital to understand their detailed mechanisms.
First order electroweak phase transitions are particularly interesting as they could help to satisfy Sakharov's third condition for baryogenesis [18] -a departure from thermal equilibrium. Provided there is sufficient CP violation, they could trigger electroweak baryogenesis and explain the observed baryon asymmetry of the Universe (see e.g., Refs. [19][20][21][22] for reviews). Determining if an electroweak baryogenesis mechanism in a particular extension of the Standard Model can successfully predict the observed baryon asymmetry of the Universe is rather involved, as described by Refs. [19][20][21][22]. Nonetheless finding the critical temperature of the phase transition is a very important step in this calculation and many studies have focused on this and on the calculation of the order parameter of the first order phase transition (see for example Refs. [23][24][25][26]).
This would provide new tests of models beyond the SM (BSM), including electroweak baryogenesis, that are complementary to collider physics experiments and measurements of electric dipole moments (EDMs). Even if the energy scale associated with the phase transition is far higher than that which could be probed in any foreseeable collider experiments, gravitational wave detection is still a possibility. For example, it was shown in Ref. [57] that gravitational waves from the breakdown of a Pati-Salam group [4] at about 10 5 GeV can give rise to detectable signatures at the Einstein Telescope.
The gravitational wave spectrum is determined by thermal quantities such as the nucleation temperature (and rate) of the bubbles, T n , the latent heat released during the phase transition (relative to the radiation energy density of the plasma), α, the (inverse) duration of the transition, β, and the velocity of the expanding bubble walls, v w . There are significant subtleties involved in calculating these and a variety of approaches have been taken in the literature [58][59][60][61][62]. Nonetheless, finding the first order phase transitions and their associated critical temperatures is an important step towards testable gravitational wave predictions.
To help study phase transitions and their associated phenomenology, we present PhaseTracer , a code to analyse the possible phases and phase transitions as the Universe cooled for any scalar potential. There are two other public codes with similar goals: CosmoTransitions [63] and BSMPT [64]. 1 Our code is most similar to the former, which facilitated over a hundred studies. Compared to Cos-moTransitions, the advantages of PhaseTracer are speed, as it is written in C++ rather than Python; flexibility, as it can be linked to models in BSMPT and Flexible-SUSY [66][67][68]; robustness, as it includes a test-suite of known analytic and numerical results and correctly handles discrete symmetries; and lastly, active maintenance and development by several authors. Our code is public and we encourage contributions (such as repository pull requests) from the community.
BSMPT and PhaseTracer have different strategies for finding and tracing the temperature dependence of vacua.
In the cases that we tested, however, PhaseTracer was able to identify more transitions than BSMPT. The public code Vevacious uses a numerical polyhedral homotopy continuation method, a powerful way to find all the roots of a system of polynomial equations, to identify vacua at zero temperature. This technique is more robust than the one used presently in PhaseTracer . However it requires the potential to be written in a special symbolic format and can only be used directly to find the tree-level minimum at zero temperature. After this, adjustments to the minima locations from one-loop Coleman-Weinberg corrections are then found using the MINUIT algorithms [69]. Minima at finite temperature are also found via MINUIT, using each minimum at zero temperature as a starting point and assuming that these two minima belong to the same phase, which is not always true.
With PhaseTracer our aim is to improve the calculation of phase transitions on several fronts. We intend to locate all minima at zero temperature and a high temperature and find the thermal phases of the potential by following the temperature dependence of the minima. We attempt to reliably distinguish phases during the evolution of the minima with temperature. We try to properly map the temperature dependence of all phases being aware that minima may appear and disappear with the temperature evolution. Finally, we aim to 1 A C++ version of CosmoTransitions is adopted by Vevacious [65] but the latter is mainly set up to test whether or not the observed electroweak vacuum is unstable and should have transitioned to a deeper underlying minimum.
properly identify the potential transitions between the identified thermal phases including cases when the potential may have discrete symmetries.
The structure of our paper is as follows. In Sect. 2 we provide quickstart instructions for installing and running PhaseTracer . Then in Sect. 3 we describe the physics problem and our numerical methods to solve it, as well as presenting the zero-and finite-temperature corrections that we include in our calculations. After this we briefly describe the structure of the code in Sect. 4 before providing detailed instructions and examples on how to implement new models in Sect. 5. In Sect. 6 we provide an extensive set of example results and comparisons between PhaseTracer and other codes or analytic results. Finally in Sect. 7 we present our conclusions.
• The NLopt library, 3  Furthermore, one of our example programs (THDMIISNMSSMBC) uses a potential built with Flexible-SUSY. Our build script attempts to build FlexibleSUSY, but there are further dependencies; see the FlexibleSUSY manuals [66,68] or the FlexibleSUSYREADME file for further details. 4 See http://eigen.tuxfamily.org. 5 See http://www.boost.org.
To visualize our results, our package includes plot scripts that require: • Python and the external modules scipy, numpy, matplotlib and pandas. • gnuplot On Ubuntu, they can be installed by $ sudo apt install gnuplot python-numpy python-scipy python-matplotlib python-pandas Finally, for testing, we provide BSMPT programs that require • libCMAES. If this is absent, it is built automatically but this requires the autoconf and libtool build tools.

Downloading and running PhaseTracer
PhaseTracer can be obtained by doing, git clone https :// github . com / PhaseTracer / PhaseTracer where the master branch will have the latest stable version, and the development branch will have the latest indevelopment features.
To build PhaseTracer on a UNIX-like system with the Make build system installed as the default build tool, run at the command line: The resulting library and executables are located in the lib/ and bin/ subdirectories of the main package directory, respectively. The built-in examples can be executed by running $ ./ bin / r u n _ 1 D _ t e s t _ m o d e l -d $ ./ bin / r u n _ 2 D _ t e s t _ m o d e l -d $ ./ bin / s c a n _ Z 2 _ s c a l a r _ s i n g l e t _ m o d e l -d The results described in Sect. 6 will be generated, including Fig. 3. By default, the example THDMIISNMSSMBC is not compiled, as it needs FlexibleSUSY. To build the THDMIISN MSSMBC model, run at the command line: Finally, by default comparisons with BSMPT are not compiled, but may be enabled by the cmake -D BUILD_WITH_ BSMPT=ON flag.

The physics problem
Our aim is to trace the minima of the effective potential for some arbitrary BSM physics theory, and find the critical temperatures for possible phase transitions between them. The first task is to construct the effective potential.

Constructing the effective potential
The effective potential typically includes the following contributions, where V tree is the tree-level potential expressed in terms of MS parameters. The contribution V 1-loop stands for the one-loop zero temperature corrections to the effective potential, which in the R ξ gauge are given by [70], where the sums over φ, V and f represent sums over the scalar, vector and fermion fields respectively and we have introduced the renormalization scale, Q, the field dependent MS mass eigenstates of the considered BSM theory m i , where there i stands for the BSM fields which are summed over, and finally n i , which gives the numbers of degrees of freedom for field i. With the potential constructed in this way we are assuming that the MS parameters have been chosen or adjusted to fulfill the electroweak symmetry breaking con-ditions at the one-loop order, as would automatically be the case if one is obtaining them from a mass spectrum generator. Please also note that if the parameters are input in a different scheme then one may need to modify these corrections, see e.g. section 3 of Ref. [71] for a review of the known one-loop corrections in the MS, DR [72,73] and DR [74] schemes in the Landau gauge. One could also add counterterms to the potential 6 in (1), i.e. change V eff → V eff + V c.t. modifying the MS scheme. This can be used to e.g. ensure that the one-loop zero temperature minimum matches the tree-level minimum and that the pole masses correspond to tree-level mass eigenstates.
The contribution V

1-loop T
is the one-loop finite-temperature correction to the effective potential and in the R ξ gauge it is given by, where J B and J F are the thermal functions, The upper (lower) sign above applies for bosons (fermions). Strong Boltzmann suppression occurs in these thermal functions for m 2 T 2 . Finally if the Debye masses are provided then the daisy corrections can be included. Following the Arnold-Espinosa approach we include these in an additive potential of the form, where in the first two terms we sum over the Higgs fields (including Goldstone bosons) and in the second two terms we sum over the massive gauge bosons. We use m 2 to denote field-dependent mass eigenvalues that include Debye corrections. If there are no Debye corrections to the masses this contribution will vanish. We always take the real part of the potential, which deals with the possibility of tachyonic masses in (5). If Debye masses are not supplied by the user then this contribution will be neglected.
The effective potential presented is not gauge invariant [70,75]. In our implementation, as described in Sect. 5, we allow the user to choose any R ξ gauge for the calculation, allowing them to test variation in a variety of gauges. In future releases we plan to add additional methods to allow an alternative gauge-independent calculation based on the procedure in Ref. [70]. We also note that there are well-known problems associated with perturbative approaches near phase transitions [76,77]. Nevertheless, lattice approaches are not currently a tractable alternative for phenomenological studies in arbitrary models of new physics.

Critical temperature and transition strengths
Having constructed a temperature dependent effective potential, we wish to trace its possible phases as the temperature changes. Tracing phases means tracking the temperature dependence of the minima of the effective potential. We will describe our method to do this in a separate section below. The result of this tracing are the phases of the potential, that is the locations of all minima as functions of temperature.
After we identified the phases of the potential we also wish to find critical temperatures at which the potential has degenerate minima. We define the critical temperature, T C , to be the temperature at which the potential has two degenerate minima, Here x min and x min are two distinct minima in the field space, that is they are separated by a barrier. Without loss of generality we can assume that the deepest vacuum of the system was x min above T C . In most cases, below T C the minima x min becomes the deepest vacuum. Just below the critical temperature, however, the system may remain at x min , i.e., in a false vacuum. Somewhat below T C the system may fluctuate over or tunnel through the barrier to the lower minimum, x min , called the true vacuum [78][79][80].
We also define the quantity which is useful for characterising the strength of the transition, and N is the number of scalar fields considered. By default all scalar fields in the potential are considered in the quadratic sum defining v(T ), but as described later this may also be restricted, allowing, for example, one to output γ EW ≡ (v EW (T c ) − v EW (T C ))/T C which may be more relevant for electroweak baryogenesis.
The strength of a first order phase transition can be a useful heuristic for assessing whether the phase transition is strong enough for a successful electroweak baryogenesis mechanism, with γ EW > 1 being a commonly applied rule of thumb. Similarly the strength of a first order phase transition also has an impact on the detectability of gravitational waves that are generated from the transition.
However it is important to note that for a full calculation of the baryon asymmetry of the Universe or gravitational wave spectra many other steps are required, which Phase-Tracer alone does not take care of. For example another very important quantity, the so-called nucleation temperature can only be calculated by using PhaseTracer in combination with a cosmological bounce solver, such as Bubble-Profiler [81], which calculates the bounce action for the potential at a given temperature. A first order transition proceeds via bubbles of broken phase nucleating in space-time. If these bubbles prove to be stable and grow, somewhat below T C there will be a nucleation temperature, T n , at which the mean number of bubbles nucleated within the relevant spacetime volume is one. For the early Universe the time of nucleation can be written as where H is the Hubble parameter and (T ) is the nucleation rate per unit volume, which may be expressed as, where S E is the bounce action that can be obtained from a bounce solver. Ultimately, the physical problem becomes computing observable quantities -such as gravitational wave signatures or the baryon asymmetry of the Universe -using some of the above thermal properties of the transition. For reviews on how to do such calculations see e.g. Refs. [1,[20][21][22]82,83]. While PhaseTracer does not provide all of these quantities, the outputs of Phase-Tracer , namely the phase structure and its temperature evolution, the critical temperatures and the strengths of the transitions, are important ingredients in such calculations and have quite broad utility.

Locating minima of the potential
We locate all minima of the potential at T = 0 and T = 1000 GeV (which may be changed by PhaseTracer:: set_t_low and PhaseTracer::set_t_high). To do so, we first generate a set of guesses by uniform sampling. We polish the guesses using local minimization and keep the unique minima. There are obviously more sophisticated techniques for global optimization that we do not implement. Nonetheless our method is fast and we have tested that it gives reliable results (see Sect. 6). If the performance of our implementation is a bottleneck or fails, we advise that a user implements techniques tailored to their particular physics problem or supplies guesses for the locations of the minima (PhaseTracer ::set_guess_points). After locating the minima at low and high temperature, we trace the high-temperature minima down in temperature, and the low-temperature minima up.

Tracing a minima
After identifying the minima, we trace them with temperature using PhaseFinder::trace_minimum with an estimate of the derivative of the fields with respect to temperature, Note that the minima, x min , carries an implicit dependence on temperature, i.e., by x min we denote the trajectory of a minimum with temperature. We find the derivative by noting that must vanish when evaluated along the trajectory of a minima, x min , since by definition of the trajectory for all temperatures, T , and all fields, x i . Thus we find the derivative in Eq. 9 by solving This is implemented in PhaseFinder::dx_min_dt. The derivatives of the potential with respect to the fields and temperature are found by a finite difference method. We use it to predict the minima at temperature T 1 = T 0 + T , We polish this prediction with the PhaseFinder::find_min method to find x 1min . Let us furthermore define a postdiction for x 0min in a similar manner, Thus we trace a minimum in steps of T by guessing the minimum at temperature T 1 = T 0 + T , polishing that guess, and iterating. We interrupt tracing if the step in temperature T is smaller than a user-specified level (changed by PhaseTracer::set_dt_min) and if any of the following occurs: • We encounter a discontinuous jump in the minimum.
We define a discontinuous jump by a significant difference between the minima, x 1min and x 0min , and a significant difference between a minima and its expected location, that is, between x 1min and x predict or between x 0min and x postdict . A significant difference is defined by the absolute and relative tolerances, x_abs_jump and x_rel_jump, that we use in Eq. 16. Their default values are shown in Table 1. • The minimum goes out of bounds or into a forbidden region of field space. • The Hessian is not positive semi-definite at what should be a minima. There is a possibility that the PhaseFinder ::find_min method may return a saddle point, instead of a local minimum. Thus if the Hessian at x 1 has any negative eigenvalues then the phase ends. • The determinant of the Hessian matrix is zero, as this indicates that there is probably a transition. Note that this criteria (disabled by PhaseTracer::set_check_hessian_ singular(false)) can separate phases by first-and second-order transitions (without it, they are separated only by first-order ones).
If tracing isn't interrupted, we stop once we reach the desired temperature. If a jump or the Hessian indicated that a phase ended, we check where the minimum went by finding and tracing the new minimum. For performance, before tracing a minimum, we check that it doesn't belong to an existing known phase. This enables our method to efficiently trace intermediate phases that do not exist at the low or high temperature.
At each iteration, we re-evaluate the step size. If we find reason for interrupting tracing, we reduce the step size by a factor of two. If the guess of the minima, x guess lay far away from the real minima, we reduce the step size by a factor of two; if not, we increase the step size by a factor of two. We define far away by the relevant absolute and relative tolerances, x_abs_identical and x_rel_identical, that we use in Eq. 16. Their default values are shown in Table 1.
Our strategy for tracing a minima relies on local optimization of a reasonable guess, found from Eq. 9. By default, we use a Nelder-Mead variant called subplex [84] implemented in Nlopt [85] as LN_SBPLX. This can fail. If we are tracing a second-order phase transition, the minimum won't change smoothly and thus our guess, based upon a first order Taylor expansion, may be poor. If the potential is pathological in the vicinity of the local minimum -e.g., particularly flat -the local optimization may fail even with a reasonable guess. Unfortunately, second-order transitions can present both problems simultaneously. In these cases, a phase may be incorrectly broken into two pieces that are discovered separately by the algorithm, since due to the discontinuity, tracing up and down can yield different results. The field values at the joint may be unreliable and the phases could slightly overlap.
To remedy this, note that we can predict first-and secondorder transitions by singularities in the Hessian matrix in Eq. 12. For second-order transitions, a singular Hessian matrix means that there may be multiple solutions for the change in minimum with respect to temperature. The presence of multiple solutions for dx min /dT indicates a possible jump discontinuity in dx min /dT at that temperature and thus a second-order transition. In cases in which a phase ends by a first order transition, on the other hand, two extrema -a minima and a barrier -merge into a double root. The Hessian matrix must vanish at a double root. Let us denote the magnitude of the smallest eigenvalue of the Hessian matrix at x min at temperature T by λ T min and at zero temperature by λ 0min . We judge the Hessian matrix to be singular if where hess is a relative tolerance for this, stored in a member of the PhaseFinder class of PhaseTracer , double hessian_singular_rel_tol. This has a default value of 0.01 which we recommend for most cases, but this may be reset by the user, via the method set_hessian_ singular_rel_tol(double tol). We use the smallest eigenvalue of the Hessian at zero temperature as an appropriate numerical scale, as if Eq. 15 is satisfied it means that there must be cancellations in the Hessian matrix caused by the finite-temperature corrections at that particular temperature. This partially alleviates the potential problems discussed in the previous paragraph.

Dealing with discrete symmetries
Discrete symmetries are common in scalar potentials. In a model with n discrete symmetries, each phase is potentially duplicated 2 n times (though it may transform trivially under some transformations). The duplicates, however, cannot be ignored. Let us denote two phases by P and Q, which transform to P and Q by a discrete symmetry of the model. Suppose that at some critical temperature a transition from P → Q is possible. The discrete symmetry means that the transitions P → Q and P → Q are equivalent, and that the transitions P → Q and P → Q are equivalent, but it does not mean that e.g., P → Q and P → Q are equivalent. In fact, in the presence of n discrete symmetries, each transition belongs to a set of up to 2 n inequivalent "cousin" transitions.
Thus, to account for discrete symmetries efficiently, we allow a user to specify them if they want in the virtual std::vector<Eigen::VectorXd> apply_symmetry (Eigen::VectorXd phi) function. This insures that only one copy of a phase is traced, but all copies related by discrete symmetries are taken into account when calculating possible transitions. The function should return the result of each discrete symmetry on the set of fields, i.e., for n symmetries the returned std::vector should have size n. We do not address global or local continuous symmetries, leaving it to a user to gauge fix them appropriately in their potential if they wish. See Sect. 6.2 for an example.

Finding possible first order transitions
Having identified the phases, we check for possible first order transitions (FOPTs) between them. For every pair of phases, we find the temperature interval in which they coexist. If the interval is non-vanishing and if the difference in potential between the two phases, V , changes sign over the interval, we look for a critical temperature with a root-finding algorithm. By default, we assume there can be no more than one critical temperature between any two phases. This can be switched off by TransitionFinder:: set_assume_only_one_transition(false) and specifying the minimum separation between critical temperatures by e.g., TransitionFinder::set_separation(1.5). With these settings, we would look for critical temperatures in every 1.5 GeV interval in temperature in which the phases coexisted. The current default value for this setting is 1 GeV, based on experience with phase transitions that we have tested PhaseTracer on, i.e. the examples discussed in Sect. 6.
To check the validity of critical temperatures by eye, in potential_line_ plotter.hpp we provide a plotting routine that shows the potential along a straight line between the true and false minima at the critical temperature. We suggest that it is called as e.g., PhaseTracer:: potential_line_plotter (EffectivePotential:: Potential &P, tf.get_transitions()). For a genuine critical temperature, we should see degenerate minima separated by a barrier. We warn readers that we encountered cases in which numerical artefacts in the potential -e.g., small jump discontinuities in potentials that were constructed in a piece-wise manner -were located by our code and mistaken for minima. Thus, we advise checking a few transitions by eye if the potential contains any possible numerical pathologies.

PhaseTracer structure
The problem is divided into three steps: constructing the potential, finding phases and checking for critical temperatures between them. They are performed by the Effective-Potential library, and the PhaseFinder and Transition Finder classes, respectively. The latter are implemented in the phase_finder.{hpp,cpp} and transition_finder .{hpp,cpp} files in the source code, respectively. The structure is illustrated in Fig. 1.

The EffectivePotential library
New models are implemented in our EffectivePotential library, which provides pure virtual classes that represent a tree-level (EffectivePotential::Potential) and a one-loop effective potential (EffectivePotential:: OneLoopPotential). The latter automatically includes oneloop zero-and finite temperature corrections, as well as daisy corrections.
PhaseTracer comes with example models that publicly inherit from the pure virtual EffectivePotential:: Potential class. The pure virtual classes make it easy to implement a new model. For example, we may implement the simple 1D model in Sect. 6.1 as, We have overriden two pure virtual methods: V, the scalar potential as a function of the fields and the temperature, and get_n_scalars, the number of scalar fields in this problem. We have furthermore overridden the virtual method forbidden, which ensures that our field always has nonnegative values. In the implementation of this model packaged with PhaseTracer in models/1D_analytic_test_ model.hpp we furthermore implement analytic derivatives of the potential. In this simple example we are not including any additional perturbative corrections.
To include such perturbative corrections one can instead implement one-loop effective potentials through the pure vir-tual EffectivePotential::OneLoopPotential class. For example, we may implement the 2D example in Sect. 6.2 as  This time we have overridden five methods of the pure virtual OneLoopPotential class to define our model: V0, the tree-level potential; get_scalar_masses_sq, the field-dependent scalar squared masses; get_scalar_dof , the numbers of degrees of freedom for the scalars; get_n_scalars, the number of scalar fields (2) in this problem; and lastly, apply_symmetry, which returns the result of the model's discrete symmetry, φ → −φ, on the fields.
The OneLoopPotential class itself implements methods for the one-loop zero-temperature and finite-temperature corrections to the potential. Fermion and vector contributions are calculated if the methods get_{fermion/vector} _masses_sq and get_{fermion/vector}_dofs are implemented. By default, it works in the ξ = 0 (Landau) gauge, but this may be changed to any R ξ gauge by the OneLoopPotential::set_xi method. Note, though, that it is up to the user to correctly implement the ξ -dependence of the scalar masses if they depart from ξ = 0. Daisy contributions are added using the Arnold-Espinosa [86] method if Debye masses are supplied by implementing get_{fermion /vector}_debye_sq. Furthermore, counter-terms can be added to the potential by overriding the virtual method counter_term.
The only methods that must be overriden (i.e., that are pure virtual), however, are the tree-level potential and the number of scalar fields. The get_scalar_dof and get_scalar_mass_sq methods aren't compulsory -they default to one degree of freedom per scalar field and numerical estimates of the eigenvalues of the Hessian matrix of the tree-level potential, respectively. If you change the numbers of degrees of freedom, it is vital to change get_scalar_mass_sq too, as the order of the eigenvalues when found numerically won't be stable and won't corre-spond correctly to the intended degrees of freedom per scalar field. For this reason and for accuracy, wherever possible it is wise to implement derivatives and eigenvalues analytically. By default, there are no fermions or vectors in the model and the apply_symmetry method trivially returns the coordinates with no changes.
The OneLoopPotential class automatically includes one-loop thermal corrections to the potential in (3). This requires an implementation of the thermal functions in (4). We interpolated them from a set of look-up tables generated from CosmoTransitions.

Running PhaseTracer
PhaseTracer may be called in an example program as follows (using the EffectivePotential::OneDimModel as an example), # include < iostream > # include " models /1 D_test_model . hpp " # include " phase_finder . hpp " # include " t r a n s i t i o n _ f i n d e r . We see the number of phases found, and for each phase, the minimum and maximum temperature, the corresponding field values and potential, and an explanation about why the phase ended. In this case, the first phase ended at T = 1000 GeV because that was the highest temperature under consideration (see PhaseFinder::set_t_high) and at T ≈ 33 GeV because the fields made a discontinuous jump. Each phase is numbered by a key, e.g. === phase key = 0 ===. The line std::cout << tf; produces output about the transitions with the format: We see the number of potential first order phase transitions, and for each one, the keys of the false and true phases (i.e., === transition from phase 0 to phase 1 ===), followed by the true and false vacua at the critical temperature, information about which elements of the field changed in the transition, the critical temperature, transition strength γ ( or γ EW if TransitionFinder::set_n_ew_scalars is used) and difference in potential between the true and false vacua at the critical temperature. The latter serves as a sanity check: it should of course be close to zero.

The Phase and Transition objects
The objects containing this information may be accessed though PhaseFinder::get_phases() and Transition Finder::get_transitions(), which return a std::vector of Phase and Transition objects, respectively. The Phase object is a struct containing, amongst other things, the attributes: • std::vector<Eigen::VectorXd> X -The field values, x min , through the phase • std::vector<double> T-The temperature, T , through the phase • std::vector<Eigen::VectorXd> dXdT -The change in minima with respect to temperature, dx min /dT , through the phase • std::vector<double> V -The potential, V (x min , T ), through the phase • phase_end_descriptor end_low -The reason why the phase ended at the lowest temperature • phase_end_descriptor end_high -The reason why the phase ended at the highest temperature // Summary information about the transition which would tell us the critical temperature of the first transition, and then print a summary of information about the true phase, the key corresponding to the true phase, and finally a summary of information about the transition.

Settings for finding phases and transitions
There are many adjustable settings that control the behavior of PhaseFinder and TransitionFinder objects, listed in Tables 1 and 2, respectively. They are altered by setters, i.e., set_{setting_name} methods, and read by get_ {setting_name} methods. The detailed usage of settings is introduced in Sect. 3 when we describe our algorithms of tracing phases and finding critical temperatures.
The attribute n_ew_scalars in both PhaseFinder and TransitionFinder is the number of scalar fields charged under electroweak symmetry, which we assume are the first n_ew_scalars scalar fields. By default PhaseFinder checks that the zero-temperature vacuum agrees with v = 246 GeV only if a user sets n_ew_scalars to a non-zero value. By default TransitionFinder calculates the transition strength γ using all scalar fields; it calculates γ EW only if a user sets n_ew_scalars.
When comparing floating point numbers, we typically check relative and absolute differences, i.e., our checks are of the form, When tracing phases in temperature, we restrict the largest and smallest possible step sizes. The largest permissible stepsize dt_max is the minimum of dt_max_abs and the temperature interval multiplied by dt_max_rel. Similarly, the smallest permissible step-size dt_min is the maximum one found from the relative and absolute settings.

Plotting results
The header include/phase_plotter.hpp provides a function phase_plotter(PhaseTracer::TransitionFinder tf, std::string prefix = "model") that takes a TransitionFinder object and plots the phases, using prefix to construct file names for the resulting plots and data files: • {prefix}.dat -Data file containing phases and transitions in a particular format • phi_1_phi_2_{prefix}.pdf etc -The phases and transitions on plane of the first and second field, e.g., the left panel of Fig. 3 • phi_T_{prefix}.pdf -The minima, x min , versus temperature, T , for every phase, e.g., the right panel of Fig.  3 • V_T_{prefix}.pdf -The potential, V (x min , T ), versus temperature, T , for every phase.
The plots are produced by the script make_plots/phase_ plotter.py. To make plots with L A T E X fonts, export MATPLOTLIB_LATEX. Note that in models with discrete symmetries supplied by a user, the plots do not show cousin transitions (see Sect. 3.5).
We furthermore provide include/potential_line_plotter .hpp to check critical temperatures by plotting the minima between the true and false vacuum at the critical temperature -see Sect. 3.6 for further details.

One-dimensional test model
To exhibit the usage of PhaseTracer , we consider a onedimensional potential, for φ ≥ 0 and κ < 0 as a simple dummy model, without physical meaning. We consider the specific numerical constants, The simplicity of the model means that it can be solved analytically, giving T C = κ 2 + 4λm 2 4cλ = 59.1608 (19) with minima at 0 and −κ/(2λ) = 50 separated by a barrier at −κ/(4λ) = 25. We show the potential and the phase transition in Fig. 2. We find a numerical result of T C = 59.1608 with minima at 4.92767e−6 and 50.0003, i.e. the numerical result matches the exact analytic result up to at least six sig-   7 This result may be reproduced by 7 The agreement on T C is often much better than the nominal relative precision, TC_rel_tol = 1e-4, as the toms748_solve rootfinding algorithm converges rapidly in many cases. However in principle it can be worse than the nominal precision, because there is uncertainty on the two minima.

$ ./ bin / r u n _ 1 D _ t e s t _ m o d e l -d
The argument --d -indicates that we want to produce debugging information including plots named *1D_test_model.
pdf. Relative tolerance for judging whether a transition was trivial as fields did not change trivial_abs_tol 1e−3 GeV Absolute tolerance for judging whether a transition was trivial as fields did not change assume_only_one_transition true Assume at most one critical temperature between two phases separation 1 GeV Minimum separation between critical temperatures Fig. 2 The one-dimensional potential at different temperatures (left) and the subsequent phase structure of the model (right). On the right panel, the lines show the field values at a particular minimum as a function of temperature. The arrows indicate that at that temperature the two phases linked by the arrows are degenerate and thus that a FOPT could occur in the direction of the arrow. The label x was autogenerated by our plotting code. In this example it stands for φ

Two-dimensional test model
To compare the performance of PhaseTracer with Cos-moTransitions, we adopt the test model with two scalar fields in CosmoTransitions. The tree-level potential is The resulting field-dependent mass matrix for the pair of scalar fields is 8 8 We correct a typo for the off-diagonal elements in Ref. [63]. The typo was not present in the CosmoTransitions-2.0.3 source code.
Here m 1 , m 2 , v and μ take the values listed in Table 1 of Ref. [63], i.e. 120 GeV, 50 GeV, 246 GeV and 25 GeV, respectively. In addition to the scalar fields in the tree-level potential, we add a scalar boson X with n X = 30 degrees of freedom and the field-dependent mass where Y 2 A = 0.1 and Y 2 B = 0.15. We set the renormalization scale in the one-loop potential to 246 GeV.
We compare our results against CosmoTransitions-2.0.3 (at present the latest version) using its default settings. 9 After checking our implementation of the potential, by making sure we obtain exactly the same potential for given field values and temperature, we test the critical temperature and true and false vacuums found by our code. Our results may be reproduced by

$ ./ bin / r u n _ 2 D _ t e s t _ m o d e l -d
The argument --d -indicates that we want to produce debugging information including plots named *2D_test_model.

pdf.
We show in Table 3 that our results agree extremely well 10 with CosmoTransitions for the critical temperature and field values for the first possible transition; however, our code is about 100 times faster for this problem. 11 The phase transitions in this problem are shown in Fig. 3. We see in the left panel a parametric plot of the two fields as functions of temperature. The distinct phases are shown by blue, orange and green colors, and the FOPT is shown by a black arrow. In the right panel, we see each field as a function of temperature.
Note, however, that PhaseTracer found two transitions whereas CosmoTransitions only found one. In the Cosmo-Transitions example, the duplicated phase and thus the cousin transition is removed by forbidding φ 1 < 0 and it returns only the transition with the smaller action. 12 However, this appears to be a coincidence: forbidding φ 2 < 0 instead of φ 1 < 0 would result in CosmoTransitions finding only the transition with the larger action 13 . Note that in the automatically generated plots, such as Fig. 3, the symmetric cousin transitions are not shown, as their inclusion makes the plots hard to read. 10 In fact in this case they agree up to 9 significant figures, with default settings. 11 We computed timings of PhaseTracer in C++ using the average of 1000 repeats in example/time.cpp and example/time .hpp with a desktop with an Intel Core i7-6700 CPU @ 3.40GHz processor. For CosmoTransitions, we averaged the run time of "for i in range(1000): model1().calcTcTrans() " inside examples/testModel1.py of CosmoTransitions-2.0.3 on the same machine. 12 In the code, it is φ 1 < −5 to allow for slight inaccuracy in the location of the phase. 13 In fact, if φ 2 < 0 is forbidden, CosmoTransitions-2.0.3 would only find Phase 2 shown in Fig. 3.

Z 2 Scalar Singlet Model
As a more realistic two-dimensional example, we consider the Z 2 scalar singlet extension of the SM. Unlike the SM, this model can accommodate a 125 GeV Higgs and a first order phase transition, as well as possibly providing a viable dark matter candidate. The potential consists of the ordinary SM Higgs potential and gauge-and Z 2 invariant interactions involving a real singlet, After EWSB, the real scalar obtains a tree-level mass m 2 S = μ 2 s + 1 2 λ Hs v 2 . We implement one-loop and daisy corrections to this tree-level potential. Under certain circumstances, a first order transition between an EW preserving and an EW breaking vacuum may be possible and the critical temperature can be found analytically; see, e.g., Ref. [36] for an analysis of the PTs in this model.
We reproduce the behaviour of the critical temperature as a function of λ Hs and m S that was found in Ref. [36] (see, e.g., Figure 1 of Ref. [36]) in Fig. 4. On the right panel, we show that the differences between numerical results from PhaseTracer and the analytic formulae are at most about 0.01%. To reproduce it, $ ./ bin / s c a n _ Z 2 _ s c a l a r _ s i n g l e t _ m o d e l $ python make_plots / c o m p a r e _ a g a i n s t _ a n a l y t i c _z2 . py The first command scans the relevant parameter space and writes a data file named Z2ScalarSingletModel_Results .txt that the second command plots in ms_lambda_Z2_SSM .pdf.

Two-Higgs-Doublet Models
We next consider three Two-Higgs-Doublet Models (2HDMs) which are implemented in BSMPT. We do not reimplement their potentials; instead, we link PhaseTracer directly to the models implemented in BSMPT and directly call the BSMPT potentials. The general 2HDM tree-level scalar potential is The two Higgs doublets take the form  The star represents a phase that does not change with temperature. The meaning of the right panel is same as the right panel in Fig. 2. The labels x 1 and x 2 were autogenerated by our plotting code. In this example they stand for φ 1 and φ 2 , respectively Second, we consider the complex, CP-violating 2HDM (C2HDM) with the parameters Type = 1 t a n β = 4.64487 Note that in this case we consider a complex phase for λ 5 . Lastly, we consider the Next-to-2HDM (N2HDM) with the parameters Type = 1 t a n β = 5.91129 Compared to the 2HDM models, this model contains an extra complex singlet with Z 2 symmetric interactions, The singlet field S is expanded as S = v S + ζ S with v S denoting singlet VEV. To run our code on these models, $ ./ bin / run_R2HDM $ ./ bin / run_C2HDM $ ./ bin / run_N2HDM Note that these examples require BSMPT; see Sect. 2.2 for installation instructions. Regarding finding the critical temperature, our code is slower than BSMPT as we perform more checks and thoroughly map out the whole phase structure. We show results from all three models with our code and BSMPT in Table 4. In BSMPT, the tolerance of the bisection method used to locate the T C is 0.01 GeV , while our default relative precision TC_rel_tol is 0.01%. Thus our results for R2HDM and C2HDM benchmark points are in agreement with that from BSMPT in the range of these errors. In addition, we find an extra group of transitions at T C = 161.205 GeV for R2HDM, which is missed by BSMPT because it focuses only on the transition that starts in the symmetric vacuum. For the N2HDM benchmark point, the discrepancy of T C is larger than the expected precision. This is related to the fact that BSMPT assumes that the false vacuum lies at the origin and with its default settings BSMPT misses the minimum around at (0, 0, 0, 0, 301.0) at T 121.00 GeV, mistakenly treating the minimum around at (0, 0, 32.5, 176.5, 297.1) as the deepest minima at that temperature. 14 By increasing the number of random starting points for finding multiple local minima in BSMPT, it finds this minima and agrees with our result, T C = 120.73 GeV, though the agreement may be partially an accident since the false vacuum does not lie at the origin, contrary to the assumption in BSMPT.

Next-to-Minimal Supersymmetric Standard Model
In Ref. [26] we explored the Next-to-Minimal Supersymmetric Standard Model (NMSSM) with a preliminary version of PhaseTracer . We include our NMSSM model (in which we match the NMSSM to a model with two-Higgs doublets and a singlet; see Ref. [26] for further details) as an example. The effective potential depends upon FlexibleSUSY (see Sect. 2.2 for the relevant installation instructions) for the matching conditions, tadpole equations and field-dependent masses. We show a benchmark point from this model in Because this model point is particularly pathological, there is a chance that PhaseTracer may miss Phase 4 and part of Phase 3 shown in Fig. 5. This is because at T = 127.1 GeV, there exists a saddle point at (h u = 0 GeV, h d = 211.0 GeV, s = 0 GeV) near a minimum. When PhaseTracer traces the phases near this saddle point, it may correctly find the real minimum near the saddle point or mistake the saddle point for a minimum. In the latter case, when PhaseTracer discovers that the Hessian isn't positive semi-definite, it stops tracing Phase 3 at the saddle point, and then will miss Phase 4. We will address this problem by improving the PhaseFinder ::find_min method in a future version so that it cannot mistake saddle points for minima. The behaviour of PhaseTracer in these cases isn't deterministic since it depends on random number generation for locating minima.  The labels x 1 , x 2 and x 3 were autogenerated by our plotting code. In this example they stand for h u , h d and s, respectively

Conclusions
In this paper we have presented PhaseTracer , a fast and reliable C++ software package for finding the cosmological phases and the critical temperatures for phase transitions. For any user supplied potential, PhaseTracer first maps out all of the phases over the relevant temperature range. 15 Once the phases have been identified PhaseTracer checks each pair of phases to see if there may be a first order phase transition between them and calculates the critical temperature.
To do this PhaseTracer uses a similar algorithm to that of CosmoTransitions, but our implementation has a number of advantages: (i) PhaseTracer is orders of magnitude faster, as we have demonstrated in Sect. 6.2; (ii) PhaseTracer has been carefully designed to provide clear error information for cases where a solution cannot be found and has many adjustable settings which may be varied to find solutions in such cases; (iii) PhaseTracer will be maintained by an active team of developers (the authors of this manual) and is distributed with a suite of unit tests to validate the code and ensure results do not change unintentionally as the code evolves. In addition Phase-Tracer has been designed so that it can be easily linked with FlexibleSUSY spectrum generators, making it easy to embed it in a much wider tool set for investigating the phenomenology of an arbitrary BSM extension. Furthermore we have also made it possible to link to potentials implemented in BSMPT, another tool for finding critical temperatures. This makes it possible for users to compare results in both BSMPT and PhaseTracer , as we have demonstrated for two Higgs doublet and two Higgs plus singlet extensions of the Standard Model in Sect. 6.4. This is particularly useful since BSMPT has a complementary approach, implementing a method that is simpler and faster, but finds exactly one transition that starts from the symmetric vacuum. Therefore BSMPT will miss other transitions and the transition it does find may not be in the cosmological history, while PhaseTracer will find all potential transitions.
With PhaseTracer it is therefore very easy to find the phases and critical temperatures of arbitrary Standard Model extensions. This is an important and useful tool for investigating both electroweak baryogenesis and gravitational waves, and can be used in thorough phenomenological investigations of realistic Standard Model extensions with many parameters, as we have demonstrated in Ref. [26] where an early prototype of PhaseTracer was used. Furthermore Phase-Tracer is part of a wider goal to develop a set of tools 15 The upper and lower bound on the temperature is chosen by the user, though by default these are taken to be T = 0 and T = 1 TeV, which are the values relevant for electroweak phase transitions. that can be used to automatically calculate the baryon asymmetry of the Universe, predict the stochastic gravitational wave background from first order phase transitions and test this against observations from gravitational wave experiments. A. also acknowledges the hospitality of Nanjing Normal University while working on this manuscript. The work of C.B. and Y.Z. was supported by the Australian Research Council through the ARC Centre of Excellence for Particle Physics at the Tera-scale CE110001104. The work of P.A. and C.B. is also supported with the Australian Research Council Discovery Project Grant DP180102209.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Our paper describes a code that is publicly available. There is no dataset associated with our paper.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .