Three-loop effective potential of general scalar theory via differential equations

We consider the scalar sector of a general renormalizable theory and evaluate the effective potential through three loops analytically. We encounter three-loop vacuum bubble diagrams with up to two masses and six lines, which we solve using differential equations transformed into the favorable $\epsilon$ form of dimensional regularization. The master integrals of the canonical basis thus obtained are expressed in terms of cyclotomic polylogarithms up to weight four. We also introduce an algorithm for the numerical evaluation of cyclotomic polylogarithms with multiple-precision arithmetic, which is implemented in the Mathematica package cyclogpl.m supplied here.


Introduction
The effective potential [1,2] plays a very important role in investigating spontaneous symmetry breaking. In the Standard Model, it has been a topic of numerous studies over many years [3][4][5][6][7], with partial results at the three- [8] and four-loop orders [9]. The analysis of the effective potential at the two-loop order leads to the conclusion that the electroweak vacuum may be stable, critical, or slightly meatastable up to very high energies of the order of the Planck scale [10][11][12]. Such an analysis at the three-loop order would require, apart from four-loop renormalization group functions [13][14][15][16][17][18], also the matching conditions, which are presently known through the two-loop order only [19]. An important step has recently been taken in Ref. [20], where the three-loop potential has been studied in a general renormalizable theory evaluating the loop integrals numerically [21][22][23].
In this work, we consider the purely scalar sector of a general renormalizable theory and evaluate the effective potential through three loops analytically. The integrals that appear in our calculation can have up to two different mass scales. In the case of O(n) symmetry, we reproduce the known result for the scalar ϕ 4 theory with spontaneous symmetry breaking. This theory is a matter of interest for the study of phase transitions, and its effective potential has been calculated in a series of papers [24,25]. The results presented here also reproduce the contribution of the scalar sector to the effective potential in the Standard Model.
This paper is organized as follows. In Section 2, we introduce the Lagrangian of the scalar sector and parametrize the three-loop effective potential in terms of three-loop master integrals. In Section 3, we discuss the evaluation of the master integrals with two different mass scales with the help of differential equations.
In Section 4, we present our results. The numerical evaluation of the cyclotomic polylogarithms that appear in our results is discussed in the appendix.

Effective potential in the scalar theory
Let us consider the scalar theory described by the Lagrangian where τ i , λ i , and λ ij are couplings. In the Standard Model, the scalar part of the potential reads In the broken phase, it can be parametrized as where φ is the vacuum expectation value, H is the Higgs field, and G 0 and G ± are the scalar would-be Goldstone bosons. This parametrization corresponds to the following choice of parameters in Eq. (1): The three-loop contribution to the effective potential of the theory in Eq. (1) can be schematically represented as the following sum of Feynman diagrams: Here, use the following notation: in each term, the symmetry and color factors are shown explicitly and the part corresponding to the scalar loop integral with all divergences properly subtracted is depicted as a figure.
We use loop functions with subtracted divergences as introduced in Ref. [20] to separate the calculation of the complicated three-loop integrals from the renormalization and the calculation of the simple lower-loop integrals.

Evaluation of the three-loop integrals
Let us introduce some general notations. All the three-loop integrals which appear in this calculation can be mapped, depending on the masses attached to the internal lines, on one of the three topologies A, B, or C shown in Fig. 1. We set m 2 = 1 and define the ratio x = m 2 1 /m 2 2 . After this rescaling, the respective master integrals read where k ab = k a − k b and the loop integration measure is defined as d[k i ] = d d k/π d/2 e εγ E , with d = 4 − 2ε being the dimension of space-time and γ E being Euler's constant. The integrals of topology A are single-scale integrals and have been known for a long time [26]. Recently, they have been evaluated through weight six [27]. The integrals of topologies B and C, appearing in lines (5)-(10), have two mass scales and depend on the ratio x. For each of these two topologies, we have constructed a set of reduction rules with the help of the LiteRed package [28] and identified the set of the master integrals. Differentiating all the master integrals with respect to x and reducing the right-hand sides to the set of master integrals, we obtain the following systems of differential equations: The elements of the matrices M b and M c are rational functions of the mass ratio x and the space-time dimension d. It is convenient to switch from the original basis {J b , J c } to the new, canonical basis {g b , g c } of master integrals introduced in Ref. [29]. The integrals in the canonical basis possess the property that the coefficients of their expansions in ε have uniform transcendentality weight, and the matrix of the differential equations has a special, so-called ε form. An early discussion of the uniform transcendentality weight of the ε expansion may be found in Refs. [30,31]. For topologies B and C, we choose the following sets of master integrals to be transformed to the canonical bases by suitable transformation matrices: To find the explicit forms of the transformation matrices T b and T c and of the matrix of the differential equations in the ε form, we use the public tools Fuchsia [32] and CANONICA [33], which are implementations of the algorithms of Refs. [34] and [35], respectively. Both algorithms heavily rely on the polynomial, rational form of the transformation matrices T b and T c . To fulfill this condition, we first make the variable transformation illustrated in Fig. 2, upon which we can successfully find the transformation matrices T b and T c . The corresponding systems of differential equations then take the form Now, the matrices B and C are independent of the space-time dimension d. We can proceed further and decompose them into sums of constant matrices with all dependence on the kinematic variable y factorized out. These decompositions read: The matrices B a,b and C a,b are constant, with rational entries, and all dependence on y is contained in the functions f b a (y), defined as where Φ n (y) is the n-th cyclotomic polynomial relevant for our calculation. These are given by As a nice property, these functions have transparent integration rules, which lead to a special type of functions, namely iterated integrals which generalize harmonic polylogarithms. These so-called cyclotomic harmonic polylogarithms were introduced in Ref. [36]. The systems of differential equations in ε form in Eq. (15) have the nice property that, after the expansion of all the master integrals in ε, the differential equations for the series coefficients completely decouple and can be written in the form and similarly for topology C, where g b {ε n }(y) denotes the O(ε n ) coefficient of the ε expansion of g b (y).
Integration of the decompositions in Eq. (16) leads to an iterative integration of the functions in Eq. (17). Starting from the lowest order of the expansion in ε, which is just a constant, we can integrate order by order introducing the definition of the cyclotomic harmonic polylogarithm of weight zero, H (y) = 1, and using the rule At each step of integration, we need to fix the integration constant. This is done by using the expansions of the integrals in the small-x limit. These expansions are not naive, but contain, on top of the power-like terms, also logarithms due to subgraphs.

Results
The results for the integrals in the canonical basis up to transcendental weight four can be found in the ancillary files of our submission to the arXiv. To test the validity of the obtained results, we perform a number of comparisons with already known results and carry out numerous expansions in different limits. Two-scale integrals with topologies B and C reduce in the limit x → e iπ 3 to one-scale integrals with topology A and can be calculated with the help of the MATAD package. The four-line [37] and five-line [22] integrals have already been known before. Therefore, we only present here the expressions for the most complicated integrals with six lines.
After calculating all the needed integrals, we are ready to present our final results for the loop functions in lines (5)-(10). We do this for each diagram separately. Defining L g = ln(m 2 g /µ 2 ) and L h = ln(m 2 h /µ 2 ), we have for the diagrams with three one-loop subgraphs: For the diagrams with convolutions of two-loop and one-loop subgraphs, we have: For the three-loop four-line diagrams, we have: For the three-loop five-line diagrams, we have: For the three-loop six-line diagrams of ladder type, we have: For the three-loop six-line diagram of Mercedes Benz type, we have: The symbols S2, T1ep, E3, DM, DN, and D6 denote the finite and O(ε) parts of single-scale integrals and are defined in Ref. [26].
Here, we only present the most complicated six-line integrals: In particular, the solution of Eq. (A.2) around x 0 = 0 is given by the product where c is the constant vector (0, 0, . . . , 0, 1) T , which is determined by the boundary conditions. In order to analytically continue h(x) away from x = 0, we follow the idea of Refs. [42,43]. We can match the solutions in different regions at some particular point that belongs to both regions. We find that, to cover the unit circle completely, we can evaluate W (x) at x 0 = 0 and at six other points placed symmetrically on the unit circle, x k = e iπk/6 , k = 0, 1, . . . , 5. We fix the boundary conditions at x 0 = 0 as described above. The matchings to the six other expansions can be taken at the pointsx k = 1 2 e iπk/6 , i.e. in the middle of the straight lines connecting x 0 and x k .
The above algorithm has been realized in the Mathematica package cyclogpl.m. It allows one to evaluate, with multiple-precision arithmetic, the cyclotomic polylogarithms in Eq. (A.1) with b j ≤ 12 and arbitrary weights. An alternative numerical implementation of cyclotomic harmonic polylogarithms is described in Refs. [44,45].