DRalgo: a package for effective field theory approach for thermal phase transitions

DRalgo is an algorithmic implementation that constructs an effective, dimensionally reduced, high-temperature field theory for generic models. The corresponding Mathematica package automatically performs the matching to next-to-leading order. This includes two-loop thermal corrections to scalar and Debye masses as well as one-loop thermal corrections to couplings. DRalgo also allows for integrating out additional heavy scalars. Along the way, the package provides leading-order beta functions for general gauge-charges and fermion-families; both in the fundamental and in the effective theory. Finally, the package computes the finite-temperature effective potential within the effective theory. The article explains the theory of the underlying algorithm while introducing the software on a pedagogical level.


Introduction
The origin of Baryon asymmetry [1] in the universe remains obscure. As such, much powder has been spent throughout past decades to find a sound explanation of this asymmetry [2][3][4][5]. Amongst the suggested mechanisms, the one based on the electroweak phase transitionelectroweak baryogenesis -stands out. While the Standard Model (SM) has no strong firstorder transition on its own [6][7][8], its extensions can contain myriads of them. For new field content to trigger a strong first-order phase transition, their masses have to be around the electroweak scale, and their interaction with the SM Higgs cannot be too feeble. Beyond-the-Standard-Model (BSM) theories that exhibit such transitions provide a direct target for many future-generation colliders [9,10]. Furthermore, thermal phase transitions can, perchance, generate gravitational waves that are observable by next-generation space-based detectors such as LISA [11], DECIGO [12], Taiji [13] and BBO [14]. These waves might open a fresh window into the early universe -and the underlying quantum field theory.
The interplay between BSM phenomenology and gravitational waves is among the most actively studied topics in the high-energy-physics literature [15][16][17][18][19][20][21][22][23][24][25]. This interplay requires solid understanding of thermodynamic properties of different models. It has been long known that determining thermodynamics in theories with non-Abelian gauge fields is challenging because of the Linde problem [26]. In short, there are non-perturbative effects arising from massless vector-bosons: these infrared (IR) effects can only be captured by lattice simulations [27]. Despite this, leading-order perturbation theory is frequently used as an approximation. However, it has been pointed out that such leading-order studies contain large theoretical uncertainties, due to slow convergence of perturbation theory [20,[28][29][30].
Dimensional reduction [31,32] offers a way to overcome these challenges. In this framework, ultraviolet (UV) modes -related to non-zero Matsubara modes in the imaginary time formalism -are integrated out. The resulting effective field theory (EFT) [33,34] describes IR, or long wavelength zero modes, living in three spatial dimensions (3d) (cf. also [35][36][37]). The 3d EFT can be simulated on a lattice and hence by-pass the Linde problem [27,38]. Dimensional reduction can also be viewed as a systematic scheme for thermal resummations, used purely within the realm of perturbation theory [29]. To leading order, dimensional reduction is equivalent to resummation via thermal masses. Concretely, the effective potential, that describes the free energy of thermal plasma, is commonly given in a schematic form [28] V eff ≃ V T =0 + V T + V daisy . (1.1) Here, V T =0 is the Coleman-Weinberg zero-temperature contribution, which is often augmented with the effective potential at one-loop level. Namely, V T is the thermal correction at one-loop (commonly in the high-temperature expansion) and Corrections at NLO are particularly important due to large logarithms. The idea with dimensional reduction is to render some logarithms (lnμ T ) small by matching at a high energy (μ ∼ T ). The remaining logarithms can then be resummed by RG-evolution within the EFT.
Phase transitions often occur below the Debye-mass scale. In such cases temporal scalars can also be integrated out. Thus effectively removing large logarithms of the form ln m D µ , where m D is a Debye mass. This second EFT is said to live at the ultrasoft scale, which is characterised by energies of O(g 2 T ). In summary, 2 the hard scale corresponds to energies E ∼ T ; the soft scale corresponds to energies E ∼ gT ; and the ultrasoft scale corresponds to energies E ∼ g 2 T .
To illustrate next-to-leading order dimensional reduction, we consider a schematic model with scalar mass parameter µ 2 , scalar quartic coupling λ, and gauge coupling g. Given the power counting µ 2 ∼ g 2 T 2 , λ ∼ g 2 , the matching of the mass parameter is where the first line (with even powers of g) results from the first step, and the second line (with odd power of g) from second step of the dimensional reduction. In practice, full O(g 4 ) contributions are included. Going to higher orders, requires a three-loop computation for both steps of the dimensional reduction. The situation is similar for the coupling: In practice for the coupling, O(g 4 ) pieces are neglected since their numerical effect is small (despite being formally of the same order). These contributions arise during the second step of the dimensional reduction at two-loop. Pursuing higher order, requires a two-loop computation for the first step, and a three-loop computation for the second step of the dimensional reduction.
In perturbation theory the definition of 3d EFT parameters are accompanied with a perturbative computation of the effective potential [29]: (1.5) Instead of expanding the result in terms of 4d parameters of the parent theory, resummed couplings and masses are kept along. This improves the overall convergence as the result is less sensitive to the renormalization scale [50,51]. The order O(g 5 ) requires a computation at three-loop level [52][53][54].

Dimensional reduction from a zero-temperature EFT perspective
It is instructive to consider the same physics from a standard EFT perspective. In particular, for equilibrium observables we can, analogous to Kaluza-Klein theories [55], view thermal corrections as an infinite tower of heavy particles. To see the connection with dimensional reduction, we consider a theory with one light and one heavy scalar in Euclidean spacetime: Imagine now that there exists a hierarchy M 2 ≫ m 2 . This is a well-known situation, and when calculating scattering processes one encounters large corrections to m 2 scaling as κM 2 -analogous to the hierarchy problem. As well as large logarithms in the form lnμ 2 M 2 or lnμ 2 m 2 . One of these logarithms will then be large regardless the choice of the RG-scaleμ.
Using thermal masses is equivalent to using a resummation 3 m 2 → m 2 + aκM 2 . This resummation does, however, not solve our problems of large logarithms. Therefore, calculations, even with thermal masses, are sensitive to the RG-scale. The situation is much improved by using an EFT. Therein, the logarithms lnμ 2 M 2 are rendered small by matching at the scalē µ Match ∼ M , and the remaining logarithms lnμ 2 m 2 are taken care of by RG-evolution in the low-energy theory [56].
Once the scalar field Φ has been integrated out, the resulting EFT is of the form To leading order m 2 eff = m 2 + aκM 2 , and λ eff = λ. While at NLO where b, c, d, e are numerical coefficients. It is precisely these kind of corrections that appear in the dimensionally reduced theory. Finally, to do calculations in the IR, we should run the couplings within the EFT from µ Match ∼ M down toμ IR ∼ m eff . In this way, all large logarithms are eliminated and perturbative convergence is improved.
In the case of dimensional reduction, we have an infinite tower of heavy particles with masses M 2 n = (2πnT ) 2 for all integers n = 0. 4 Hence, the high-temperature matching scale is µ Match ∼ T and the low-energy scale is µ 2 IR ∼ g 2 T 2 . Since we are describing equilibrium (static) processes, there is no time dependence -the degrees of freedom of the EFT are static and live in three dimensions.

Installation and running
This section explains how to install DRalgo and presents a tutorial based on the Abelian-Higgs model.

Installation
The current version of DRalgo, 1.0, is installed by placing all the source files either in the applications folder Mathematica/Applications or running Mathematica from the package root directory ./DRalgo.
The required source files are outlined in fig. 1. To create model files, DRalgo uses functions from GroupMath [57]; see ibid. for its installation. This tutorial uses GroupMath Version 1.1.2 but GroupMath is not required for the dimensional reduction itself. If a model file is already available, it can be loaded independently of GroupMath (cf. sec. Q.1). Since GroupMath is an external package, any use of the model-creation features in DRalgo should be accompanied by a citation of [57].
To load DRalgo from the Mathematica applications folder Mathematica/Applications, the following commands need to be executed:

Model implementation
As an example, we consider the Abelian-Higgs model with the Euclidean Lagrangian Here, the field A µ is an U(1) gauge field with gauge coupling g 1 and φ is a complex scalar field charged under U(1) with general charge Y φ . The zero, or temporal, component of the vector field receives a thermally induced mass in the 3d theory. This temporal scalar transforms as a singlet under U(1) and has a Lagrangian where indices i ∈ {1, . . . , d} are spatial. Here, µsqU1 is the temporal-scalar Debye mass (squared), λVLL [1] is the self-interaction coupling, and λVL [1] is the portal coupling to φ. The notation above is aligned with the systematic notation of DRalgo. For earlier studies of this model in the literature, we refer to [35-37, 58, 59].
To create a model in DRalgo, we have to specify the gauge group and the scalar representation. In DRalgo this is done via Group={"U1"}; CouplingName={g1}; RepAdjoint={0}; Higgs={{Yφ},"C"}; RepScalar={Higgs}; RepFermion={}; Here, the adjoint (vector) representation is trivial for the Abelian Higgs model. In the definition of the Higgs scalar {{Yφ},"C<R>"}, Y φ denotes the gauge charge. One additional argument is passed depending on if the scalar has a real (R) or a complex (C) representation. In this simple example, fermions are absent and the corresponding bracket is empty. While not relevant for this model, note that fermions are never part of the 3d EFT Lagrangian as they are integrated out during the first step of the dimensional reduction. The model-input from the gauge sector is generated with the command {gvvv,gvff,gvss,λ1,λ3,λ4,µij,µIJ,µIJC,Ysff,YsffC}= AllocateTensors[Group,RepAdjoint,CouplingName,RepFermion,RepScalar]; The tensors listed above have the following correspondence gvvv structure constants gvff vector-fermion trilinear couplings gvss vector-scalar trilinear couplings λ1 scalar tadpole couplings λ3 cubic couplings λ4 quartic couplings µij scalar-mass matrix µIJ, µIJC fermion-mass matrices Ysff, YsffC Yukawa couplings and their purpose is described in depth in sec. 3. By default all these tensors, except for the gauge ones, are empty. While their order on the left hand side of AllocateTensors is fixed, the above naming of {gvvv,...,YsffC} is arbitrary.
The next task is to specify the scalar potential. Let us start with the mass matrix. There is only one allowed term: φφ † . This term is selected via To understand the syntax, recall that we defined our scalar as RepScalar={Higgs}. Then {1,1} specifies that we are after a term with two (RepScalar[ [1]]) φ fields. For general models, the index 1 can be replaced by any index in RepScalar as defined by the user. In addition, {True,False} specifies that we are after the φφ † term and conversely {True,True} corresponds to a φφ term. The latter term, however, is not allowed in this example due to gauge invariance. Since this example exhibits merely one U(1) invariant combination of φ and φ † , we specify [ [1]] in the CreateInvariant output.
For the quartic tensor only one possible term arises, namely (φφ † ) 2 . We already created a (φφ † ) term above, which can be reused to find the quartic: To see how quartic and mass terms are composed of scalar-field components, one can inspect the variables VMass and VQuartic. This completes the model implementation, and we now have all ingredients to do the dimensional-reduction step.
The actual dimensional reduction is performed with the command This calculates all thermal masses, effective couplings, β-functions and anomalous dimensions. For example, gauge and scalar couplings are given by where the output is a replacement rule. The variables Lb and Lf are dependent on the 4d-renormalization scale (μ) and the temperature (T ) [29,33]. They are where γ E is the Euler-Mascheroni constant. Further information about the functional basis of the effective parameters is collected in sec. Q.11. The corresponding definitions of constants, along with other shorthand notations, are shown with the command The above RG-scaleμ is the hard matching scale, and should be chosen asμ ∼ T to avoid large logarithms. 5 Effective mass parameters are divided into scalar and temporal-scalars. The scalar masses are given by Here, the LO command prints tree-level masses together with one-loop thermal contributions; the NLO command prints two-loop thermal masses. By default these NLO masses contain the running from the matching scale (μ) to an arbitrary 3d scale (μ 3 ). 6 The resulting Debye masses are displayed by Such temporal-scalar couplings are denoted as λVLL[a] for V 4 couplings, λVVSL[a] for V 2 S couplings and λVL[a] for V 2 S 2 couplings, where S represents scalar fields and V temporal scalar fields. The index 'a' systematically labels all linearly independent terms.
The coefficient of the unit operator [34], or the hard-mode contribution to the symmetricphase pressure, is given by 5 In DRalgo the renormalization scales are denoted by µ =μ for the hard scale and µ3 =μ3 for the soft scale. 6 There are different ways to express the NLO masses. We refer to sec. Q.11 for further details.
Here, LO corresponds to one-loop, NLO to two-loop, and NNLO to three-loop level. The LO result describes the pressure of an ideal gas, and is given by P LO (T ) = N π 2 90 T 4 , where N represents the number of degrees of freedom. For the Abelian Higgs model N = 4 with 2 from the photon and 2 from the complex scalar.
DRalgo also provides anomalous dimensions and beta functions in the parent 4d theory; see sec. 3.4. Beta functions are printed with the command Next, to find anomalous dimensions we first need to specify for which particles we want them. As DRalgo stores everything in tensor form, we need to find the positions of all particles
See sec. Q.4 for how to obtain anomalous dimensions for vector-bosons and fermions.

Integrating out temporal scalars
Temporal scalars are often heavy compared to the fields driving the phase transition, and can be integrated out [29,33]. It should be stressed that this is a user-dependent optional step. The resulting EFT is said to describe ultrasoft physics. The command is Couplings at the ultrasoft scale are given by And effective scalar masses are given by In the 3-dimensional theory only scalar masses have non-zero beta functions.
This concludes our tutorial of the dimensional reduction of the Abelian-Higgs model.

Two-loop effective potential
Once the model is loaded, we can calculate the effective potential. This can be done by the user, either in the soft or the ultrasoft theory. Alternatively, DRalgo can calculate the effective potential in the ultrasoft theory. First, we need to create field-dependent masses which requires to specify a vacuum expectation value direction: Here the vacuum expectation value is in the second, imaginary, Higgs component.
For the calculation to proceed, all mass matrices must be diagonal. If not, an error message is printed and it is up for the user to diagonalize the matrices; see section Q.16 for an example. For the model at hand, the mass matrix is diagonal, which the user can confirm by printing the field-dependent masses Here, LO refers to the tree-level, NLO to the one-loop, and NNLO to the two-loop effective potential. Note that all results are given in Landau gauge. In the (two-loop) NNLO part, the renormalization scale is denoted by µ3 =μ 3 .
This completes the first tutorial on quick installation and running. Output from the 3d EFT matching relations can be implemented either to non-perturbative lattice codes, or perturbative analyses in terms of the effective potential. While such implementations are left to the user, appendix A displays how to interface DRalgo output in a Mathematica implementation that determines selected thermodynamic quantities for the Abelian-Higgs model.

Theory behind the scenes: dimensional reduction of a generic model
In this section, we dig deeper into the theory and computations in the back-end of the software. The implementation for dimensional reduction of a generic model is based on tensornotation [60][61][62][63][64] that separates Lorentz algebra from group algebra. The Lorentz algebra is hard-coded, with contractions done by FeynCalc [65,66] and FORM [20,67], while the groupstructure is provided by the user.

Lagrangian in the 4d fundamental theory
Consider a general Lagrangian in Minkowski spacetime with mostly-plus metric: g µν = diag(−1, 1, 1, 1), {σ µ , σ ν } = −2g µν . In terms of the functional-integral measure and the action S, the partition function is Z = D e iS . Correspondingly, the sigma matrices are defined as where σ i are Pauli matrices.
To write down the Lagrangian, in a flattened form, we employ a basis where all scalars and vectors are real. In addition, all fermions are composed of two-component Weyl-spinors [68]. The most general, four dimensional, Lagrangian in Minkowski space is [34,60,61] where L int is the interaction Lagrangian and F µν = i g [D µ , D ν ] is the field strength tensor with the corresponding gauge coupling g. In this notation, the field R i denotes a real scalar-field with scalar-index i. In the Standard Model, i corresponds to the Higgs, neutral Goldstone, and real/imaginary component of the charged Goldstone. Further, the field A a µ corresponds to a real vector-field with vector-index a, the field η is a ghost-field and the field ψ I is a Weylspinor with fermion-index I. In addition, µ ij denote (squared) scalar masses, and M IJ fermion masses. Repeated indices are always summed over irrespective of their vertical placement.
Above we denoted Y iIJ as Yukawa couplings, λ ijkl as scalar quartic, λ ijk as scalar cubic, and λ i as scalar tadpole couplings. The gauge couplings g a jk are (anti-symmetric) representation matrices. For example, for an adjoint scalar-representation in SO(3) g a jk = ǫ ajk . The abbreviation h.c. stands for the hermitian conjugate. The vertical placement of fermion indices is important; e.g. M IJ = M * IJ . We refer to [68] for further details. In the dimensional-reduction step, hard modes with masses ∼ πT are integrated out. This is done by matching the fundamental Lagrangian above to the effective Lagrangian living in three-dimensions. To avoid large logarithms this matching should be performed close toμ = πT whereμ is the RG-scale. This three-dimensional Lagrangian does not contain fermions. Moreover, in the dimensional-reduction step the temporal component of vectorsrepresented by temporal scalar fields in the EFT -obtain Debye masses, as well as thermally generated interactions with other scalars.
The matching is straightforward using Euclidean signature. This is achieved by redefining where quantities are denoted either in Minkowskian (M ) or Euclidean (E) metric. Henceforth, we suppress these subscripts and implicitly assume an Euclidean metric.
With the above redefinitions, the partition function is defined as Z E = D e −S E with the most general tree-level Lagrangian in Euclidean signature All Lorentz indices are contracted using the Euclidean metric δ µν = diag(+1, +1, +1, +1). The Feynman rules for the vertices are 10) where G abcd µνρσ and T µνρ are defined in [61]. Since we only focus on the matching, only the symmetric phase is relevant and vacuum expectation values are not introduced. The propagators are where all momenta are incoming by convention and where ξ is the corresponding gauge parameter which for our investigations in Landau gauge is set to ξ = 0. All the Lorentz structure is contained in the vertices, and the generalized coupling constants take care of the group structure.
To do the matching we need to renormalize our theory. The matching can, and will, introduce kinetic mixing. To allow for this, we express the bare (b) fields and scalar masses as By construction all propagators should be diagonal at tree-level, to wit

Assumed size of masses and couplings
For the current problem there are three energy scales: πT , gT , and g 2 T . Denoted as the hard, soft, and ultrasoft scale, respectively. In DRalgo it is assumed that all scalar (fermion) masses scale as gT or g 2 T . If, on the contrary, some masses would be hard, then they should formally be integrated out together with high-temperature modes. 7 After hard modes are integrated out, it is up to the user whether soft-scale (gT ) fields are integrated out as well. By default all temporal scalars have masses of this order, but some models might also have additional scalars with soft masses. For the couplings, DRalgo assumes that all quartics and cubics scale as g 2 . This scaling is chosen so that scalar masses of O(gT ) can be integrated out safely [23].
In summary, DRalgo assumes the following scaling It is of course possible to consider smaller couplings and masses than those above.

Lagrangian in the 3d effective theory
The three-dimensional theory is of the form where Lorentz indices range between µ, ν = {1, . . . , 3}. In renormalizing the fields, we define Here, the term δZ ab,L 3 corresponds to temporal vectors. The propagators of this theory are as in eq. (3.18) in three dimensions with the absence of the fermionic propagator and the inclusion of the temporal-vector propagator The counterterm Feynman rules are the same as above, barring the new temporal-vector diagram

Matching
Self-energies and general n-point correlations are calculated both in the original 4d theory and the effective 3d theory. The 3-dimensional parameters are chosen, or matched, to give the same correlators as the original theory [29,33].

One-loop thermal scalar masses
We start with the matching for scalar masses. For the matching we assume to be in a regime of soft or ultrasoft external momenta p ∼ gT or p ∼ g 2 T . The corresponding diagrams at one-loop level are illustrated in fig. 2. In their computation, we need the master integrals [69] where ζ s = ζ(s) is the Riemann zeta function. The d-dimensional bosonic integral measure is defined as while for fermions the summation is written as Σ {Q} . We set d = 3 − 2ǫ with Euclidean four-momenta Q 2 = q 2 n + q 2 = (2n + σ)πT 2 + q 2 , and σ = 0(1) for bosons(fermions).
The first diagram in fig. 2 gives where the factor of 1 2 is the symmetry factor, −λ ijnn is the Feynman rule, and The zero Matsubara mode is here ignored. The reason for this is that the zero Matsubara mode corresponds to soft momenta, and should not be integrated out. Conversely, since the n = 0 mode exists both in the EFT and the parent theory, this contribution cancels in the matching. The second diagram is likewise (neglecting higher ǫ terms) where we used the shorthand notation for H ab ij from eq. (3.17). The remaining diagrams give The first diagram contains the external momentum is p, and the last diagram is the contribution from a scalar-mass insertion. With the assumed power-counting (3.22), this contribution is formally of higher order. Finally, there is the fermion loop The second coupling-constant term above is the hermitian conjugate. Also, above we use ǫ b and ǫ f to denote ǫ poles, albeit with some additional factors. They are defined in [29,33], and are where Lb and Lf are given in eq. (2.4).
To perform the matching, we demand At leading order, 3d fields are rescaled by a factor T −1/2 . Next, consider the self-energies. One finds where Π ij (p) = (D 1 + · · · + D 5 ) ij . To avoid confusion we added a subscript 3 to threedimensional quantities.
Here −δZ ij p 2 − δµ ij are 4d counterterms, and cancel all ǫ poles. The matching gives In this step all the Lorentz structure is stripped away, leaving mere group-theory factors.  Figure 3: Diagrams contributing to the thermal corrections to the scalar quartics. Indices correspond to the line notation in fig. 2.

One-loop scalar quartics
Next, we consider scalar quartics. Before calculating the actual diagrams, it is useful to look at what happens on the 3d side. The relevant terms in the Lagrangian are There are two contributions to the matching. First, the coupling constant, h, receives contributions from different loop orders Second, as shown in the previous section, 3d fields are renormalized. In terms of the scalar counterterm, δZ 3,ij , this renormalization is (to leading order) R i(b) = R i + 1 2 δZ 3,ij R j . Combined, these considerations imply that the contribution from the 3d side is In the 4d theory, the first and second diagrams give Finally, the fermion contribution is where the terms (jkl) contain all permutations over indices j, k, l. The matching is performed by demanding In the light of the discussion at the beginning of this section, this implies The renormalized sum of the diagrams above was denoted as Λ ijkl = (D 1 + D 2 + D 3 ) ijkl and one finds where δZ 3,km was given in the previous section.

Higher loop levels
When performing dimensional reduction at NLO, two-loop contributions are required for scalar thermal masses. Here, we do not show details of the two-loop computation, as they are analogous to the one-loop computation presented above. Naturally the number of required diagrams is larger. Due to the use of integration-by-parts identities (IBP) [70][71][72] it can be shown that to the given order, the two-loop master sum-integrals factorize into the one-loop master integrals given in eqs. (3.28) and (3.29). While this streamlines the computation at NLO significantly, contributions at higher orders are non-factorizable. Nonetheless, DRalgo also implements a generic computation of two-loop thermal masses and other 3d parameters at NLO.

Beta functions and anomalous dimensions
DRalgo calculates beta functions and anomalous dimensions by renormalizing 4d parameters (cf. [62][63][64]). As an example, consider the scalar sector. We have renormalized our fields and couplings as On the 4d side, we assume the MS scheme, so δZ ij , δµ ij , and δλ ijkl are chosen to only cancel ǫ poles. We demand that bare parameters are independent of the RG-scaleμ. To leading order in ǫ, we find (t = logμ) (3.52) Using this we can find the anomalous dimensions by demanding that ∂ t R i(b) = 0: where we have introduced the anomalous dimension γ ij . Then, the quartic beta function is determined by Where we have used ∂ t δλ ijkl = −4ǫδλ ijkl ; which follows from counting powers of couplings in δλ ijkl . We identify 2ǫδλ ijkl with the beta function for λ ijkl . This beta function can be written in an equivalent form by writing the counterterm, for λ, as the sum of a bare term and field renormalization: Using this one finds where the last term denotes all possible contractions between one index of λ ijkl and one index of γ ij . Other beta functions are calculated analogously. All gauge beta functions in DRalgo are defined with couplings squared, viz. ∂ t g 2 . Conversely beta functions for masses, scalar couplings, and Yukawa couplings are defined as linear in the parameters, i.e. ∂ t λ and so forth. These conventions are also shown in the DRalgo output.
In the 3-dimensional theory only scalar masses have non-zero beta functions.

Implementing beyond the Standard Model theories
In this section, we demonstrate the use of DRalgo for BSM theories based on the example of the Two-Higgs doublet model (2HDM).

Two-Higgs doublet model with fermions
Consider the 2HDM potential [73,74] The scalars are SU(2) doublets with hypercharge Y φ = 1. This means that the covariant derivative is given by 8 where τ a are Pauli matrices. We further assume that all Standard-Model fermions are present, and consider Yukawa interactions of the form We define the group structure and scalar representations via (4.4) The above notation identifies Rep1 as the left-handed quark doublet, Rep2 as the right-handed up-quark and so forth. The extra factor of 1 2 for hypercharges in the definition of e.g. Rep1 arises from the definition of the covariant derivative similar to eq. (4.2).
Additional fermion-families can be added by grouping multiple instances of RepFermion1Gen together: where the number of generations is chosen to be n f = 3. When creating Yukawa interactions the user decides which index contains which family. For example, above we stacked all generations after each other in RepFermion3Gen. Hence, in RepFermion3Gen indices 1-5 correspond to the first generation, indices 6-10 to the second, and indices 11-15 to the last. Thus the user could use index 1, 6, or 11 as the top quark. An alternative way to add an arbitrary number of fermion families n f is given in sec. Q.9.
Next, to create the tensors, we write The mass terms in the potential in eq. (4.1) are We allowed m 2 12 to be complex, with real part m12R, and imaginary part m12I. For the scalar quartics in eq. (4.1) the corresponding part of the potential is for which we already created all the building blocks. Thus, the quartics can be constructed as Consequently, the quartic tensor itself is defined as For simplicity, we assumed above that λ5H, λ6H, and λ7H are real. One can allow for complex couplings by adding them, and their conjugates, directly in QuarticTerm6 and QuarticTerm7.
To create a complex coupling λ6H, we would write QuarticTerm6=(λ6HR+I*λ6HI)*MassTerm1*MassTerm3+(λ6HR-I*λ6HI)*MassTerm1*MassTerm4; For Yukawa couplings, we only consider the top-quark coupling. If we choose the first family to contain the top-quark, then a Yukawa term ∼ φ † q L u R would involve fermion representations number 1 and 2. When defining a Yukawa interaction, representations should be specified in the order: scalar, first fermion, second fermion. With this in mind, the Yukawa coupling of the first Higgs doublet is 9  9 Here we assume that the top quark resides in the first generation. This has nothing to do with how the top-quark is usually placed in the third generation. Rather, we here place it in the first for simplicity as normally only the top-quark has a sizeable Yukawa coupling. For the sake of generality, we assumed that both Higgs doublets couple to the top-quark. Note that this is not allowed for a realistic model due to flavor-changing-neutral-current constraints.
Finally, assuming that the Yukawa couplings are real, we can define Above we did not consider Yukawa couplings between different generations but such couplings can be added. To wit, when we defined RepFermion3Gen we put the first generation on indices 1-5; the second on indices 6-10; and the third on indices 11-15. To define a coupling between, e.g. the left-handed top quark (assumed to reside in generation 1) and the right-handed charm quark (assumed to reside in generation 2), we would write With the model implementation complete, all dimensional-reduction commands are identical to those of the Abelian-Higgs model. We refer to [77][78][79] for earlier results in the literature.

Integrating out temporal scalars
For the Abelian-Higgs case in sec. 2 we only had one scalar field, whereas with the 2HDM we have two. When integrating out temporal scalars we have two options. First, all scalars are light and we have two active, dynamical doublets. Second, one doublet is heavy and is integrated out when going from the soft to the ultrasoft scale.
For the first case, the corresponding command is For the second case, we assume that the second doublet is heavy. It can then be integrated There is one complication when integrating out one of the Higgs doublets. Namely, generally the two doublets mix through the m 2 12 term. For a complete treatment the mass matrix needs to be diagonalized before the heavy doublet can be integrated out. This lies beyond DRalgo and is an optional step for the user. Instead, the code by default assumes that m 2 12 is smallof O(g 2 T ) in power counting -and performs the diagonalization perturbatively to first order in m 2 12 .

Miscellaneous features
In this section, we discuss specific features of the software on a case-by-case basis.

User-options and features
Lower-order dimensional reduction for speed The user has several options to control what is calculated and how the code operates. For example, if the model has many degrees of freedom, the user might wish to save running-time and only calculate one-loop thermal masses and couplings. This works by specifying Mode->1 when loading the model Group={"U1"}; ImportModelDRalgo[Group,gvvv,gvff,gvss,λ1,λ3,λ4,µij,µIJ,µIJC,Ysff,YsffC,Mode->1]; The default is Mode->2, in which case everything is calculated to NLO. And for the most barebone/fast option select Mode->0; in this case only one-loop thermal masses are calculated.

Model-treatment in the code
DRalgo works by factorizing all group and Lorentz algebra. The Lorentz algebra is hard-coded while the group algebra is supplied by the user. All particles are indexed by the order they appear. If the user has an SU(3) × SU(2) model with gauge bosons, the vectors in the SU (3) group have 8 components, while those of the SU(2) group have 3. The code indexes these components by a = 8 α SU(3) , 3 β SU (2) . Therein, α runs from 1-8, β from 1-3, and a from 1-11. For example, the Debye mass tensor µ ab D stores the information of both groups. For most cases this mass tensor is diagonal.
The scalars are treated differently since DRalgo only deals with real scalar components. For example, a complex scalar Φ is rewritten Φ = 1 √ 2 (φ + iψ). Hence the scalar indices would in this case be i = (φ, ψ). Further, the scalar-mass matrix is stored as µ S,ij . For example, a m 2 φ 2 mass term resides at µ S, 11 . Consider now a complex representation with n scalar components labelled by I, J. In a complex basis vector-scalar-scalar trilinear couplings are of where G a IJ are representation matrices. To convert to a real basis each scalar component is split as Φ I = 1 √ 2 (φ I + iψ I ). The components are then reordered with the real components first: i = (φ I , ψ I ), so i = 1, . . . , 2n. If we call the scalars in the real basis ϕ i , our vector-scalar-scalar trilinear couplings become .

Frequently asked questions
Below, we discuss various questions and problems the user might have encountered.

Q.1. How to save and load my model?
Saving and loading a model is straightforward with DRalgo built-in functions. Once a model is created, and loaded with ImportModelDRalgo (see the .m files in the example folder), it can be saved by writing The loaded model can be your own, or perhaps one of the provided DRalgo model repository. We also encourage you to make your own models available for the community which is possible by submitting the model file via the Issue Tracker on https://github.com/DR-algo/DRalgo. The model will then be verified and added to the model repository. When submitting a model, please refer to a paper or explicitly write out the Lagrangian in an accompanying notebook.

Q.3. How do I check that my model is anomaly free?
Once you have defined your model, the anomaly-free condition is that vanishes identically. Since this condition is not automatically fulfilled for general, non-numeric, U(1) charges, it is the responsibility of the user to choose U(1) charges such that all anomalies cancel.

Q.4. How do I calculate anomalous dimensions?
First find the position of all representations: All anomalous dimensions are then found via For more details regarding anomalous dimensions and beta functions cf. sec. 3.4.

Q.5. Temporal scalar mixing
In most models the temporal-scalar masses are diagonal. However, in cases with multiple U(1) groups there can be mixing. In such a case, the code automatically recognizes this. For example, if the group is Group={"U1","U1","U1"}; the code creates the mixed U(1) masses with the convention that µU1Mix1 is the mixing between groups 1 and 2, µU1Mix2 is the mixing between groups 1 and 3, and µU1Mix3 is the mixing between groups 2 and 3. These masses are displayed with the PrintDebyeMass command. See the ./examples/SMZp.m for an example.
Q.6. What if I miss some couplings?
By default DRalgo assumes that the user has defined all couplings allowed by symmetry in their 4d theory. The code runs even though some couplings are missed. For example, consider the 2HDM with scalar potential as in eq. (4.1) but without λ 6 and λ 7 couplings. For a φ 1 ↔ φ 2 symmetric model no λ 6 /λ 7 -type scalar quartic couplings are induced. However, if the model breaks this symmetry by e.g. a Yukawa sector, then λ 7 and λ 6 couplings are generated at one-loop. In this case DRalgo would still calculate these induced couplings; with the crux that only the λ 1 , . . . , λ 5 couplings are printed by the PrintCouplings[] command, while the induced λ 7 and λ 6 couplings must be found manually via the PrintTensorDRalgo[] command. In addition, DRalgo automatically alerts the user with a message if some new couplings are generated at one-loop. This is a cross-check that no couplings are forgotten. Q.7. Can I run DRalgo without specifying the representation?
For general groups no. However, it is possible for the user to specify arbitrary (non-numeric) U(1) charges. With non-numeric charges the code does not check that various quartic or Yukawa terms are allowed. Hence, it is the responsibility of the user to ensure gauge invariance.

Q.8. How large representations can I use?
In principle any group and representations can be used. In practice the code slows down for huge representations. The code has been tested on general models with ∼100-200 components in a given representation. For example an SO(10) model with a 120-dimensional scalar. Since DRalgo was purposely written to deal with any quartic and Yukawa sector, these are often the bottlenecks. If the user wants to omit two-loop contributions, significantly larger representations are possible. Further still, if the user only wants one-loop thermal masses, almost any model (within reason) can be run in quick order.

Q.9. Can I include an arbitrary number of fermion generations?
Yes, take for example the 2HDM. As described in sec. 4.1, a single family of SM fermions is defined via Defining Yukawa and scalar couplings proceed as before, and the model is loaded with ImportModelDRalgo[Group,gvvv,gvff,gvss,λ1,λ3,λ4,µij,µIJ,µIJC,Ysff,YsffC,Verbose->False]; The user can then define n f = nF fermion families by writing 10 These commands should be used before running Most outputs from DRalgo can directly be compared with the literature. For temporal scalars this requires some care. For example, if the original 4d theory has a SU(2) triplet, then two types of A 2 0 Σ 2 couplings are allowed: DRalgo does not distinguish between these two terms since its output is stored in the form To compare with the two parametrizations, the user can rewrite the κ 1 /κ 2 basis in tensor form. To this end, it is necessary to additionally define the relevant terms of the effective model and compare to its Lagrangian. This procedure is analogous to creating the fundamental model where invariants can be compared using the command The matching of all possible allowed operators in the EFT is automatic in DRalgo. This way also previously disregarded effective coefficients can be determined. One example is a L 3d ⊃ κTr C 3 0 B 0 operator in the Standard Model. Here, C 0 is the gluon (temporal) field, and B 0 is the temporal hypercharge field. The output from DRalgo gives κ ∝ (2Y q + Y d + Y u ). See eq.  Tadpoles are defined analogously. See the 2xsm.m file for a worked example with two real scalars.
Q.13. Do I need to define tadpoles?
If tadpoles are allowed by symmetry we encourage the user to define them -even if they are absent at tree-level.

Q.14. How do I define Dirac and Majorana masses?
Fermion masses are defined analogously to scalar masses. We stress that DRalgo only works with Weyl fermions. For example, assume a theory with two Weyl fermions ψ 1 and ψ 2 . The following terms are then possible Where the first two terms are Majorana masses, and the third is a Dirac-type mass. These masses can be defined in DRalgo via Q. 15. What if a model has Ψ 1 γ 5 Ψ 2 terms?
Since DRalgo only uses Weyl fermions, the inclusion of other fermionic bilinear operators requires some extra work. One example is the operator Ψ 1 γ 5 Ψ 2 which can be expanded in Weyl spinors such that For most cases, scalar quartics can be easily defined using the examples above. For nonstandard representations, the model-building tools in DRalgo might differ by a basis change from other conventions in the literature. Let us take an example, a SU(5) model with an adjoint scalar. Commonly this scalar is written as Φ = Φ a T a where T a are traceless hermitian matrices satisfying Tr T a T b = 1 2 δ ab . The most general (quartic) potential is By default DRalgo uses GroupMath, which defines its invariants differently to those above. In fact, the two scalar quartic operators given in GroupMath can be linear combinations of those above. Thus, the result from DRalgo would be in the form where the invariants inv 1,2 are the output from the CreateInvariant command. Fortunately, it is quite easy to find the relations between the V 1 and the V 2 basis. To do so, first, rewrite everything in tensor form and likewise for the DRalgo output Since λ 1 and λ 2 are defined in different bases, we want to compare invariants. Here we only need three: Comparing the invariants one finds λ H = 1 960 (5.14) The above procedure works for any representation as long as the user can write the quartic operator from DRalgo in their preferred form. See SU5.m for a concrete example.
Q. 19. Can I use dimensionally reduced theory for calculating nucleation rates?
Certainly. The user can use the effective theory and directly find the bounce in the 3d. To first approximation the nucleation rate is e −S 3d , where S 3d is the bounce action in the dimensionally reduced theory. See [80] for the original work, and [58,[81][82][83] for recent examples.
Q.20. Can I study GUT models with DRalgo?
Most GUT models can be readily studied with DRalgo. However, for models with large representations the running-time rises rather rapidly. The bottleneck is not the complexity of e.g. the scalar potential but rather the size of the representations of particles. For example, The required RAM is negligible for most models. However, for models with 50-100 dimensional representations, the requirements start to rise. And for said scenario around 1 GB of RAM is required. If the package drains too much RAM, it is recommended to not include two-loop contributions.
Q.22. How should I report possible bugs?
We are grateful if any bugs are reported. We kindly ask users to supply both a short description and a Mathematica file describing the error. These reports and attached files can provided via the Issue Tracker on https://github.com/DR-algo/DRalgo.

Outlook
High-temperature field theory is pestered by large radiative corrections, which compromises perturbative calculations. In effect, perturbation theory needs to be reorganized in terms of thermal resummations. While at leading order only thermally corrected masses are required, new contributions to couplings arise at higher orders. In addition, there are currently several different schemes for incorporating thermal masses -all depending substantially on the renormalization scale.
As an effective field theory, dimensional reduction by-passes these issues. The ultraviolet sector of the theory is controlled by a matching computation and the infrared sector is resummed to all orders. Indeed, as discussed in this paper, using a three dimensional EFT unambiguously resums masses and couplings. Using the framework is simple from a practical standpoint since three-dimensional theories are super renormalizable. Concretely, only mass parameters are RG-scale dependent and their beta functions are exact at two-loop level. This property not only allows for a straightforward perturbative treatment, but also provides an attractive framework to a non-perturbative lattice study of the 3d EFT since in the effective theory relations between continuum and lattice parameters are exact at two-loop level [84].
Hitherto, dimensional reduction has been used sparsely. With this paper and the associated software, we aim to change this by automating the EFT construction. In summary, DRalgo calculates all effective couplings and masses in the effective theory, the leading-order beta functions both in the 4d parent and 3d effective theory, as well as the effective potential within the EFT. This facilitates studying models with dimensional reduction and requires merely three-dimensional calculations that are analogous to zero-temperature computations. Since the effective theory is fully bosonic, the perturbative calculation can be compared with lattice simulations. This has one clear benefit as a large number of fundamental 4d theories can map into the same effective theory given the mass hierarchy of the additional scalars. Not only does the lattice provide an invaluable cross-check, but pre-existing simulations can be reused and applied to new BSM theories.
In conclusion, gravitational waves have opened up a new gateway to the early universe, and particle physics stands at its threshold. Upcoming experiments are fast approaching both at future particle colliders and gravitational wave observatories. It is therefore important to control theoretical uncertainties at unprecedented precision. In this venture, dimensional reduction is the tool of choice, and DRalgo its harbinger.

A. Interfacing DRalgo output
In this appendix, we collect additional material that is not part of the package itself.
A.1. From 3d effective theory to thermodynamics A possible interface between the output of DRalgo, in the 3-dimensional theory, and the determination of thermodynamic quantities is implemented in ./examples/ah-thermo.m. This implementation comes with a disclaimer: we present a simplified algorithm suitable for the simple case of the Abelian-Higgs model (cf. eq. (2.1)) but this is not part of the package itself. We encourage users to implement their own numerical minimisation routines for thermodynamics optimized for individual models. The algorithm given here is not gauge invariant and should be used with discretion.
The algorithm 2 illustrates our implementation. In the simplified algorithm ah-thermo.m Here, M is the Higgs mass in arbitrary units related to the MS mass parameter by the tree-level relation µ 2 = −M 2 /2. The critical temperature T c is defined from the condition that the effective potential at the broken minimum vanishes since the potential at the symmetric minimum at the origin is normalized to zero. To find this condition, we use an elementary binary search. While the curves are barely discernible, this showcases that the scale-dependence is indeed a higher-order effect once at full NLO dimensional reduction. The additional RG-scale of the 3d EFT is not varied in the plots.
As an indicator of the transition strength, we also output the value of the background field at the critical temperature, φ c /T c , that describes the discontinuity of the minima at T c . Since φ c /T c is gauge-dependent, it should not be given physical meaning but used as a rough indicator that correlates positively with the phase transition strength. This strength can be defined in terms of the released latent heat at the critical temperature T c where p ′ is the derivative of the pressure with respect to temperature and ∆ ≡ (. . . ) low-T − (. . . ) high-T denotes the difference between the broken and symmetric phase. The T -derivative is approximated as a finite difference dT = 0.1. The symmetric-phase part does not contribute here since we normalized the effective potential to zero at all T at the origin. For a commented documentation of the technicalities of this implementation see directly the source ah-thermo.m. The described perturbative determination is done in Landau gauge with ξ = 0 in eq. (3.19) and the results for T c and L/T 4 c are gauge-dependent. Despite this simple user-friendly example, a recipe for a more sophisticated, gauge-invariant determination can be found in e.g. [85] (cf. refs. therein). The final output data for thermodynamics is stored in ./examples/ah-thermo-python-plots/*.dat, plotted in ./examples/ah-thermo-python-plots/plot.py, and shown in fig. 4.
We encourage DRalgo users to develop and optimize their individual implementations for algorithms to determine thermodynamic properties. While including efficient algorithms for determining thermodynamics is conceivable for future versions of DRalgo, in the current version 1.0 these features are not implemented in the package itself.