Six-loop $\varepsilon$ expansion study of three-dimensional $n$-vector model with cubic anisotropy

TL;DR

This paper performs a six-loop epsilon expansion analysis of the three-dimensional n-vector model with cubic anisotropy, providing refined estimates of critical properties and confirming the stability of the cubic fixed point for n=3.

Contribution

The study extends previous analyses by calculating six-loop expansions and applying resummation techniques to improve estimates of critical exponents and fixed point stability.

Findings

Cubic fixed point is stable for n=3.

Refined critical exponents for the cubic universality class.

Comparison with previous theoretical and lattice results.

Abstract

The six-loop expansions of the renormalization-group functions of $\varphi^4$ $n$-vector model with cubic anisotropy are calculated within the minimal subtraction (MS) scheme in $4 - \varepsilon$ dimensions. The $\varepsilon$ expansions for the cubic fixed point coordinates, critical exponents corresponding to the cubic universality class and marginal order parameter dimensionality $n_c$ separating different regimes of critical behavior are presented. Since the $\varepsilon$ expansions are divergent numerical estimates of the quantities of interest are obtained employing proper resummation techniques. The numbers found are compared with their counterparts obtained earlier within various field-theoretical approaches and by lattice calculations. In particular, our analysis of $n_c$ strengthens the existing arguments in favor of stability of the cubic fixed point in the physical case $n = 3$.

Tables

Table 1. Table 1: Padé triangle for the ε 𝜀 \varepsilon expansion of n c subscript 𝑛 𝑐 n_{c} . Here Padé estimate of k 𝑘 k -th order (lower line, RoC) is the number given by corresponding diagonal approximant [L/L] or by a half of the sum of the values given by approximants [L/L − - 1] and [L − - 1/L] when a diagonal approximant does not exist. Three estimates are absent because corresponding Padé approximants have poles close to the physical value ε = 1 𝜀 1 \varepsilon=1 .

$M\setminus L$	0	1	2	3	4	5
0	4	2	4.5885	$-$1.2858	15.5412	$-$41.0807
1	2.6667	3.1283	2.7917	3.0684	2.5692
2	-	2.8930	2.9576	2.8828
3	1.9518	-	2.9138
4	-	2.7887
5	0.4549
RoC	4	2.3333	3.1283	2.8424	2.9576	2.8983

Table 2. Table 2: Padé-Borel-Leroy estimates of n c subscript 𝑛 𝑐 n_{c} obtained from ε 𝜀 \varepsilon expansion ( 25 ) under the optimal value of the shift parameter b o ​ p ​ t = 1.845 subscript 𝑏 𝑜 𝑝 𝑡 1.845 b_{opt}=1.845 . The estimate of k 𝑘 k -th order (lower line, RoC) is the number given by corresponding diagonal approximant [L/L] or by a half of the sum of the values given by approximants [L/L − - 1] and [L − - 1/L] when a diagonal approximant does not exist. Two estimates are absent because corresponding Padé approximants turn out to be spoiled by dangerous poles.

$M\setminus L$	0	1	2	3	4	5
0	4	2	4.58848	-1.28584	15.5412	-41.0807
1	2.75996	3.05988	2.87042	2.92283	2.91341
2	-	2.93394	2.91132	2.91499
3	2.57775	2.91419	2.91416
4	-	2.91416
5	2.39138
RoC	4	2.3800	3.0599	2.9022	2.9113	2.9146

Table 3. Table 3: Padé triangle for the ε 𝜀 \varepsilon expansion of β 𝛽 \beta . Five estimates are absent because corresponding Padé approximants have poles lying between ε = 0 𝜀 0 \varepsilon=0 and ε = 2 𝜀 2 \varepsilon=2 .

$M\setminus L$	0	1	2	3	4	5	6
0	0.5	0.3611	0.3792	0.3421	0.4301	0.1779	1.059
1	0.3913	0.3771	0.3670	0.3682	0.3648	0.3740
2	0.3791	-	0.3681	0.3673	0.3674
3	0.3586	0.3715	-	0.3674
4	-	0.3674	0.3693
5	0.2983	-
6	-

Table 4. Table 4: Padé-Borel-Leroy estimates of β 𝛽 \beta obtained from corresponding ε 𝜀 \varepsilon expansion under the optimal value of the shift parameter b o ​ p ​ t = 3.460 subscript 𝑏 𝑜 𝑝 𝑡 3.460 b_{opt}=3.460 . Several boxes are empty because of dangerous poles spoiling corresponding Padé-Borel-Leroy approximants.

$M\setminus L$	0	1	2	3	4	5	6
0	0.5	0.3611	0.3792	0.3421	0.4301	0.1779	1.059
1	0.3952	0.3768	0.3674	0.3674	0.3664	0.3703
2	0.3808	-	0.3674	0.3674	0.3672
3	0.3653	0.3728	-	-
4	-	0.3697	0.3690
5	0.3474	0.3691
6	-

Table 5. Table 5: Padé triangle for the ε 𝜀 \varepsilon expansion of γ 𝛾 \gamma . Five estimates are absent because corresponding Padé approximants have poles close to the physical value ε = 1 𝜀 1 \varepsilon=1 .

$M\setminus L$	0	1	2	3	4	5	6
0	1	1.2222	1.3489	1.3198	1.4885	1.0451	2.651
1	1.2857	-	1.3252	1.3446	1.3663	1.3925
2	1.4275	1.3543	1.3939	1.3733	-
3	1.2905	1.3848	1.3770	1.3879
4	-	1.3754	1.3832
5	0.9069	-
6	-

Table 6. Table 6: Padé-Borel-Leroy estimates of γ 𝛾 \gamma obtained from six-loop ε 𝜀 \varepsilon expansion under the optimal value of the shift parameter b o ​ p ​ t = 0.090 subscript 𝑏 𝑜 𝑝 𝑡 0.090 b_{opt}=0.090 . Empty boxes correspond to the approximants spoiled by dangerous poles.

$M\setminus L$	0	1	2	3	4	5	6
0	1	1.2222	1.3489	1.3198	1.4885	1.0451	2.651
1	-	-	1.3263	1.3438	1.3711	1.3825
2	-	-	-	-	-
3	1.1893	1.3631	1.3604	-
4	-	1.3605	1.3629
5	1.142	-
6	-

Table 7. Table 7: The values of critical exponents for the cubic class of universality obtained by means of Padé-Borel-Leroy resummation of the six-loop ε 𝜀 \varepsilon expansions. Corresponding Padé estimates and the differences between Padé-Borel-Leroy estimates and their Padé counterparts are also presented.

n = 3	$\alpha$	$\beta$	$\gamma$	$\delta$	$\eta$	$\nu$
PBL resum.	$-$0.09(9)	0.3684(13)	1.368(12)	4.733(4)	0.036(3)	0.700(8)
Pade resum.	$-$0.11(6)	0.368(3)	1.379(8)	4.772(17)	0.038(4)	0.703(5)
Difference	0.02(11)	0.0004(33)	$-0.011\left(44\right)$	$-0.039\left(17\right)$	$-0.002\left(5\right)$	0.003(9)

Table 8. Table 8: Six-loop estimates of critical exponents versus scaling relations

Scaling relation:	1	2	3	4
Deviation from zero	-0.005(14)	0.016(13)	0.0121(36)	0.007(45)

Table 9. Table 9: The values of correction-to-scaling exponents ω 1 subscript 𝜔 1 \omega_{1} and ω 2 subscript 𝜔 2 \omega_{2} for the cubic class of universality obtained by means of Padé-Borel-Leroy resummation of the six-loop ε 𝜀 \varepsilon expansions. Corresponding Padé estimates and the differences between Padé-Borel-Leroy estimates and their Padé counterparts are also presented.

n = 3	${\omega }_1$	${\omega }_2$
PBL resum.	0.799(4)	0.005(5)
Pade resum.	0.78(11)	0.008(38)
Difference	0.02(11)	$-0.003\left(38\right)$

Table 10. Table 10: Marginal order parameter dimensionality n c subscript 𝑛 𝑐 n_{c} given by the ε 𝜀 \varepsilon expansion technique, 3D RG approach, Monte Carlo simulations and the pseudo- ε 𝜀 \varepsilon expansion machinery. By the number of loops we mean the order of approximation.

Number
of loops	$n_c$	Paper	$n_c$	Paper	$n_c$	Paper
	$\epsilon$ expansion		3D RG		Others
1	4	[5]-1974			Monte Carlo
2	2.333	[5]-1974	2.0114	[44]-1983	3	[23]-1998
3	3.128	[5]-1974	3.003	[45]-1984
4	2.918	[21]-1997	2.9	[17]-1989
	2.96(11)	This work-2019	2.89(2)	[25] -2000
	2.958	[20]-1995
5	3	[22]-1997	2.89-2.92	[24]-2000
	2.855	[21]-1997
	2.87(5)	[26]-2000
	2.91(3)	This work-2019			Pseudo-$\epsilon$ expansion
6	2.915(3)	This work-2019	2.89(4)	[26]-2000	2.86(1)	[32]-2016
					2.862(5)	[28]-2000

Table 11. Table 11: Critical exponent values given by multi-loop ε 𝜀 \varepsilon expansion calculations versus those resulting from 3D RG analysis. Error bar for four loop estimate of η 𝜂 \eta is absent because it can not be evaluated within approach described in Sec.  4.1 .

Number
of loops	$\eta$	$\nu$	Paper	$\eta$	$\nu$	Paper
	$\epsilon$ expansion			3D RG
3				–	0.700	[45]-1984
4	0.034	0.68(3)	This work-2019	0.0331	0.6944	[17]-1989
				0.0332	0.6996	[25] -2000
5	0.0375(5)	0.6997(24)	[46]-1998	0.025(10)	0.671(5)	[24]-2000
	0.0374(22)	0.701(4)	[26]-2000
	0.0353(21)	0.686(13)	This work-2019
6	0.036(3)	0.700(8)	This work-2019	0.0333(26)	0.706(6)	[26]-2000

Table 12. Table 12: Comparison of critical exponents for cubic (this work) and Heisenberg ( [ 35 ] ) classes of universality for n = 3 𝑛 3 n=3 . The numbers with asterisk were obtained from six-loop ε 𝜀 \varepsilon expansion estimates for η 𝜂 \eta and ν 𝜈 \nu via scaling relations.

n = 3	$\alpha$	$\beta$	$\gamma$	$\delta$	$\eta$	$\nu$
Cubic	$-$0.09(9)	0.3684(13)	1.368(12)	4.733(4)	0.036(3)	0.700(8)
Heisenberg	$-$0.118(6)*	0.3663(12)*	1.385(4)*	4.781(3)*	0.0378(5)	0.7059(20)

Equations

$$S=\int d^{D}{x}\Biggl{\{}\frac{1}{2}\left[(\partial\varphi_{0\alpha})^{2}+m_{0}^{2}\varphi_{0\alpha}^{2}\right]+\frac{1}{4!}\left[g_{01}T^{(1)}_{\alpha\beta\gamma\delta}+g_{02}T^{(2)}_{\alpha\beta\gamma\delta}\right]\varphi_{0\alpha}\varphi_{0\beta}\varphi_{0\gamma}\varphi_{0\delta}\Biggr{\}},$$

$$\begin{array}{r}{T}_{\alpha \beta \gamma \delta }^{(1)}=\frac{1}{3}({\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{\delta }_{\alpha \delta }{\delta }_{\gamma \beta }),\\ T_{\alpha \beta \gamma \delta }^{(2)}={\delta }_{\alpha \beta \gamma \delta },{\delta }_{{\alpha }_1\dots {\alpha }_n}=\{\begin{array}{ll}1,& {\alpha }_1={\alpha }_2=\dots ={\alpha }_n\\ 0,& \text{otherwise}.\end{array}\end{array}$$

$$T_{\alpha \beta \gamma \delta }^{(1)}{T}_{\alpha \beta \gamma \delta }^{(1)}=\frac{n(n+2)}{3},T_{\alpha \beta \gamma \delta }^{(1)}{T}_{\alpha \beta \gamma \delta }^{(2)}=n,T_{\alpha \beta \gamma \delta }^{(2)}{T}_{\alpha \beta \gamma \delta }^{(2)}=n.$$

$$\begin{array}{r}{m}_0^2=m^2{Z}_{m^2},g_{01}=g_1{\mu }^{\epsilon }{Z}_{g_1},g_{02}=g_2{\mu }^{\epsilon }{Z}_{g_2},{\phi }_0=\phi Z_{\phi },\\ Z_1=Z_{\phi }^2,Z_2=Z_{m^2}{Z}_{\phi }^2,Z_3=Z_{g_1}{Z}_{\phi }^4,Z_4=Z_{g_2}{Z}_{\phi }^4.\end{array}$$

$$S^{R}=\int d^{D}{x}\;\Biggl{\{}\frac{1}{2}\left[Z_{1}(\partial\varphi_{\alpha})^{2}+Z_{2}m^{2}\varphi_{\alpha}^{2}\right]+\frac{1}{4!}\left[Z_{3}g_{1}\mu^{\varepsilon}\;T^{(1)}_{\alpha\beta\gamma\delta}+Z_{4}g_{2}\mu^{\varepsilon}\;T^{(2)}_{\alpha\beta\gamma\delta}\right]\varphi_{\alpha}\varphi_{\beta}\varphi_{\gamma}\varphi_{\delta}\Biggr{\}},$$

$${\Gamma }_{\alpha \beta }^{(2)}={\Gamma }^{(2)}{\delta }_{\alpha \beta },{\Gamma }_{\alpha \beta \gamma \delta }^{(4)}={\Gamma }_1^{(4)}{T}_{\alpha \beta \gamma \delta }^{(1)}+{\Gamma }_2^{(4)}{T}_{\alpha \beta \gamma \delta }^{(2)},$$

$${\Gamma }_1^{(4)}=\frac{3(T_{\alpha \beta \gamma \delta }^{(1)}-T_{\alpha \beta \gamma \delta }^{(2)})}{n(n-1)}{\Gamma }_{\alpha \beta \gamma \delta }^{(4)},{\Gamma }_2^{(4)}=\frac{(n+2)T_{\alpha \beta \gamma \delta }^{(2)}-3T_{\alpha \beta \gamma \delta }^{(1)}}{n(n-1)}{\Gamma }_{\alpha \beta \gamma \delta }^{(4)}.$$

$$Z_i(g_1,g_2,\epsilon )=1+\sum _{k=1}^{\infty }{Z}_i^{(k)}(g_1,g_2){\epsilon }^{-k}.$$

$$Z_i=1+K{R}^{\prime }{\bar{\Gamma }}_i,$$

$${\bar{\Gamma }}_1=\frac{\partial }{\partial m^2}{\Gamma }^{(2)}{\mid }_{p=0},{\bar{\Gamma }}_2=\frac{1}{2}{\left(\frac{\partial }{\partial p}\right)}^2{\Gamma }^{(2)}{\mid }_{p=0}{\bar{\Gamma }}_3=\frac{1}{g_1{\mu }^{\epsilon }}{\Gamma }_1^{(4)}{\mid }_{p=0},{\bar{\Gamma }}_4=\frac{1}{g_2{\mu }^{\epsilon }}{\Gamma }_2^{(4)}{\mid }_{p=0}.$$

$$\begin{array}{r}{\beta }_i(g_1,g_2,\epsilon )=\mu \frac{\partial g_i}{\partial \mu }{\mid }_{g_{01},g_{02}}=-g_i\left[\epsilon -g_1\frac{\partial Z_{g_i}^{(1)}}{\partial g_1}-g_2\frac{\partial Z_{g_i}^{(1)}}{\partial g_2}\right],i=1,2,\\ {\gamma }_j(g_1,g_2)=\mu \frac{\partial \log Z_j}{\partial \mu }{\mid }_{g_{01},g_{02}}=-g_1\frac{\partial Z_j^{(1)}}{\partial g_1}-g_2\frac{\partial Z_j^{(1)}}{\partial g_2},j=\phi ,m^2,\end{array}$$

$${\beta }_i=g_i\left[-\epsilon +\sum _{l=1}^6\sum _{k=0}^l{C}_{{\beta }_i}^{k,(l-k)}{g}_1^k{g}_2^{l-k}\right],i=1,2,$$

$${\gamma }_j=\sum _{l=1}^6\sum _{k=0}^l{C}_{{\gamma }_j}^{k,(l-k)}{g}_1^k{g}_2^{l-k},j=\phi ,m^2.$$

$${\beta }_1(g_1^{\ast },g_2^{\ast },\epsilon )=0,{\beta }_2(g_1^{\ast },g_2^{\ast },\epsilon )=0.$$

$$\begin{split}g_{1}^{*}=\frac{\varepsilon}{n}+\varepsilon^{2}\biggl{(}-\frac{106}{27n^{3}}+\frac{125}{27n^{2}}-\frac{19}{27n}\biggr{)}+\sum_{k=3}^{6}C^{(k)}_{g_{1}}\varepsilon^{k}+\oldcal{O}\left(\varepsilon^{7}\right),\\ g_{2}^{*}=\frac{\varepsilon(n-4)}{3n}+\varepsilon^{2}\biggl{(}\frac{424}{81n^{3}}-\frac{178}{27n^{2}}+\frac{31}{27n}+\frac{17}{81}\biggr{)}+\sum_{k=3}^{6}C^{(k)}_{g_{2}}\varepsilon^{k}+\oldcal{O}\left(\varepsilon^{7}\right),\end{split}$$

$$\alpha =2-\frac{D}{2+{\gamma }_{m^2}^{\ast }},\beta =\frac{D/2-1+{\gamma }_{\phi }^{\ast }}{2+{\gamma }_{m^2}^{\ast }},\gamma =\frac{2-2{\gamma }_{\phi }^{\ast }}{2+{\gamma }_{m^2}^{\ast }},\eta =2{\gamma }_{\phi }^{\ast },$$

$$\nu =\frac{1}{2+{\gamma }_{m^2}^{\ast }},\delta =\frac{D+2-2{\gamma }_{\phi }^{\ast }}{D-2+2{\gamma }_{\phi }^{\ast }}.$$

$$g_1^{\ast }$$

$$g_2^{\ast }$$

$$\zeta (3,5)=\sum _{0<n<m}\frac{1}{n^3{m}^5}\simeq \mathrm{0.037707672985.}$$

$$g_1^{\ast }=0.33333\epsilon +0.13443{\epsilon }^2-0.13363{\epsilon }^3+0.16124{\epsilon }^4-0.43104{\epsilon }^5+1.3278{\epsilon }^6+\text{oldcal}O\left({\epsilon }^7\right),$$

$$g_2^{\ast }=-0.11111\epsilon +0.053955{\epsilon }^2+0.061933{\epsilon }^3+0.050592{\epsilon }^4-0.18841{\epsilon }^5+0.95219{\epsilon }^6+\text{oldcal}O\left({\epsilon }^7\right).$$

$$\Omega =\left(\begin{array}{cc}\frac{\partial {\beta }_1(g_1,g_2)}{\partial g_1}& \frac{\partial {\beta }_1(g_1,g_2)}{\partial g_2}\\ \frac{\partial {\beta }_2(g_1,g_2)}{\partial g_1}& \frac{\partial {\beta }_2(g_1,g_2)}{\partial g_2}\end{array}\right)$$

$${\omega }_2=\epsilon \frac{n-4}{3n}+{\epsilon }^2\frac{(n-1)(-848+660n+72n^2-19n^3)}{81n^3(n+2)}+\sum _{k=3}^6{C}_{{\omega }_2}^{(k)}{\epsilon }^k+\text{oldcal}O\left({\epsilon }^7\right),$$

$${\omega }_2(n_c,\epsilon )=0.$$

$$n_c$$

$$\eta$$

$${\nu }^{-1}$$

$$\gamma$$

$${\phi }_1\to \frac{{\phi }_1+{\phi }_2}{\sqrt{2}},{\phi }_2\to \frac{{\phi }_1-{\phi }_2}{\sqrt{2}},$$

Full text

††journal: Nuclear Physics B

1 Introduction

As is well known, the systems undergoing continuous phase transitions demonstrate the universal critical behavior. This leads to the concept of classes of universality introduced decades ago. They are determined by the general properties of the system such as spatial dimensionality, symmetry, and the number of order parameter components, thereby its microscopic nature does not play any role in the vicinity of phase transition temperature. There is a set of universal parameters such as critical exponents, critical amplitude ratios, etc. that characterize the critical behavior of the systems belonging to the same universality class.

The analysis of critical phenomena in a broad variety of materials can be performed on the base of three-dimensional $O(n)$-symmetric $\varphi^{4}$ field model. In case of one-component – scalar – order parameter ($n=1$) one deals with the Ising model describing phase transitions in uniaxial ferromagnets, simple fluids, binary mixtures, and many other systems. There is also a great numbers of substances with the vector ordering, e.g. easy-plane ferromagnets, superconductors and superfluid helium-4 ($n=2$), Heisenberg ferromagnets ($n=3$), quark-gluon plasma in some models of quantum chromodynamics $(n=4)$, superfluid helium-3 ($n=18$) and the neutron star matter ($n=10$). On the other hand, if we consider real materials with more or less complex structure, some anisotropy of the order parameter often exists. Perhaps, simplest example of such a material is a cubic ferromagnet.

Initially, to describe its thermodynamics near Curie point the $O(3)$-symmetric theory neglecting crystal anisotropy has been used. The detailed analysis performed later within the renormalization-group (RG) approach has shown, however, that for proper description of the critical behavior of real cubic crystals one should take into account the presence of the anisotropy, i. e. add to the Landau-Wilson Hamiltonian an extra term invariant with respect to the cubic group of transformations. It looks as $g_{2}\sum_{\alpha=1}^{n}\varphi^{4}_{\alpha}$, where $\varphi_{\alpha}$ is $n$-vector ordering field and $g_{2}$ – anisotropic coupling constant. This new quartic coupling, in particular, accounts for the fact that in real ferromagnets ($n=3$) the vector of magnetization ”feels” the crystal anisotropy and can lie only along the axes or spatial diagonals of cubic unit cell in the ordered phase.

This model with two coupling constants – $g_{1}$ (isotropic) and $g_{2}$ – was carefully examined since 1972 [1] by many researches. As was found, its RG equations describing evolution of quartic couplings under $T\to T_{c}$ possess four fixed points: Gaussian $(0,0)$, Ising $(0,g_{I}^{*})$, Heisenberg $(g_{H}^{*},0)$ and cubic$(g_{1}^{*},g_{2}^{*})$. One of the most important issues involved in the study is the determination of the stability of these fixed points or, in other words, what critical regime takes place in real ferromagnets. Analyzing the RG flows it was shown that the first two points are always unstable for arbitrary values of order parameter dimensionality $n$ whereas the last two of them corresponding to the Heisenberg (isotropic) and cubic (anisotropic) modes of critical behavior compete with each other. Which regime turns out to be stable depends on $n$. For $n<n_{c}$, where $n_{c}$ is some marginal value of spin dimensionality, the isotropic (Heisenberg) critical regime is stable while for $n>n_{c}$ the cubic critical behavior is realized. If initial (”bare”) values of coupling constants lie outside the regions of fixed points attraction critical fluctuations strongly modify the behavior of the system converting the second-order phase transition into the first-order one. Figure 1 illustrates the situation.

Thus, in the case $n>n_{c}$ the cubic quartic term is certainly relevant and has to be taken into account. This results in the emergence of new class of universality corresponding to the anisotropic – cubic – critical behavior. So, the value of $n_{c}$ becomes of prime physical importance since it determines the true regime of the critical behavior in real cubic ferromagnets and of some other systems of interest.

Detailed study of the $n$-vector cubic model including evaluation of critical exponents and $n_{c}$ was carried out by many groups [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] having used both field-theoretical methods and lattice calculations. Early numerical estimates of $n_{c}$ obtained in the lower-order approximations within the $\varepsilon$ expansion approach [2, 4, 5, 6] and in the frame of 3D RG machinery [12, 15, 16] turned out to be in favor of the conclusion that $n_{c}>3$, while lattice calculations implied $n_{c}$ is practically equal to 3 [13]. This made the study of the cubic class of universality less interesting from the physical point of view. Later, however, the higher-order analysis including resummation of RG perturbative series was performed and shown that numerical value of $n_{c}$ falls below 3 [17, 18, 19, 20, 21, 22, 24, 25, 26, 28, 32]. To date, the most advanced estimates of $n_{c}$ obtained within the $\varepsilon$ expansion, 3D RG and pseudo-$\varepsilon$ expansion approaches are $n_{c}=2.855,2.87$ [21, 26], $n_{c}=2.89,2.91$ [24, 26] and $n_{c}=2.86$ [28, 32], respectively.

These numbers differ from each other appreciably what may be considered as a stimulus to find the value of $n_{c}$ with higher accuracy. On the other hand, recently the $\varepsilon$ expansions of record length – six-loop – for $O(n)$-symmetric $\varphi^{4}$ field theory [33, 34, 35] were calculated. This paves the way to analysis of the critical behavior of the cubic model within the highest-order $\varepsilon$ approximation including getting precise numerical estimates for critical exponents and $n_{c}$. Such an analysis is the aim of this work.

The paper is organized as follows. In Sec. 2 we write down the fluctuation Hamiltonian (Landau-Wilson action) of $n$-vector cubic model and describe the renormalization procedure. In Sec. 3 the six-loop $\varepsilon$ expansions for $\beta$ functions, critical exponents and $n_{c}$ are calculated. The six-loop $\varepsilon$ series for cubic fixed point coordinates and critical exponents are also presented here for the physically interesting case $n=3$. In Sec. 4 the $\varepsilon$ expansions for ”observables” – $n_{c}$ and critical exponents – are resummed and corresponding numerical estimates are found. In Sec. 5 the numbers obtained are discussed and compared with their counterparts given by alternative field-theoretical approaches and extracted from the lower-order approximations. Sec. 6 contains the summary of main results and concluding remarks.

2 Model and renormalization

In this work we address the field-theoretical RG approach in spatial dimensionality $D=4-\varepsilon$.111Original six-loop calculations [33, 34, 35] were performed in space dimension $D=4-2\varepsilon$ which is more common for high energy physics. The critical behavior of the cubic model is governed by the well-known Landau-Wilson action with two coupling constants

[TABLE]

where $\varphi_{0\alpha}$ is $n$-component bare field, $g_{01}$ and $g_{02}$ being the bare coupling constants. The tensor factors $T^{(1)}$ and $T^{(2)}$ entering the $O(n)$-invariant and cubic terms respectively are as follows

[TABLE]

In particular,

[TABLE]

The action (1) is seen to be physical (positively defined) if $g_{02}>-g_{01}$ for $g_{01}>0$ and $g_{02}>-ng_{01}$ for negative $g_{01}$.

The model is known to be multiplicatively renormalizable. The bare parameters $g_{10},g_{20},m_{0}^{2},\varphi_{0}$ can be expressed via the renormalized ones $g_{1},g_{2},m^{2},\varphi$ by means of the following relations

[TABLE]

Using these relations we arrive to the renormalized action

[TABLE]

where $\mu$ is an arbitrary mass scale introduced to make couplings $g_{1}$ and $g_{2}$ dimensionless. Renormalization constants are defined in a way enabling to absorb divergences from all Green functions, so that renormalized Green functions are free of divergences. Due to multiplicative renormalizability of the model it is enough to remove divergences in two- and four-point one-particle irreducible Green functions:

[TABLE]

[TABLE]

In this paper we employ the Minimal Subtraction (MS) scheme where renormalization constants acquire only pole contributions in $\varepsilon$ and depend only on $\varepsilon$ and coupling constants:

[TABLE]

Renormalization constants can be found from the requirement of the finiteness of renormalized two- and four-point one-particle irreducible Green functions. Another way to calculate renormalization constants is use of Bogolubov-Parasiuk $R^{\prime}$ operation:

[TABLE]

where $R^{\prime}$ – incomplete Bogoludov-Parasiuk $R$-operation, $K$ – projector of the singular part of the diagram and $\bar{\Gamma}_{i}$ – normalized Green functions of the basic theory (see e.g. [36, 37]) defined by the following relations:

[TABLE]

One of the most important advantages of the Bogoludov-Parasiuk approach is that counterterms of the diagrams computed for O(1)-symmetric (scalar) model can be easily generalized to any theory with non-trivial symmetry due to the factorization of the tensor structures (see e.g. [38, 39, 40]). To calculate tensor factors for particular diagrams of the cubic model (1) one should apply projectors (7) to it. Such an operation can be automated with FORM [41] and GraphState [42] while counterterm values can be taken from data obtained in the course of recent 6-loop calculations for $O(n)$-symmetric model [35].

3 Six-loop expansions for RG functions, cubic fixed point coordinates, critical exponents and $n_{c}$

The RG functions, i. e. $\beta$ functions and anomalous dimensions $\gamma_{\varphi}$, $\gamma_{m^{2}}$ are related to renormalization constants $Z_{i}$ by the following relations:

[TABLE]

where $Z^{(1)}_{i}$ – coefficients at first pole in $\varepsilon$ from (8).

We calculated the RG functions as series in renormalized coupling constants up to six-loop order. They are found analytically and presented in Tables 1, 2, 3 and 4 of Supplementary materials (see A) in the form

[TABLE]

[TABLE]

The critical regimes of the system are controlled by the fixed points $(g_{1}^{*},g_{2}^{*})$ of RG equations that are zeroes of $\beta$ functions:

[TABLE]

As was already mentioned, for the model under consideration there are four fixed points: Gaussian $(0,0)$, Ising $(0,g_{I}^{*})$, Heisenberg $(g_{H}^{*},0)$ and cubic $(g_{1}^{*},g_{2}^{*})$. Since six-loop $\varepsilon$ expansions analysis of Ising and Heisenberg models have been performed earlier [33, 34, 35] we concentrate on the cubic critical behavior. To calculate $\varepsilon$ expansions for critical exponents we have to find those for coordinates of the cubic fixed point. Solving (14) by means of iterations in $\varepsilon$ for the cubic fixed point we find:

[TABLE]

where higher-order coefficients $C^{(k)}_{g_{1}}$, $C^{(k)}_{g_{2}}$ are presented in Tables 5 and 6 of Supplementary materials (see A).

To fully characterize the cubic class of universality, we need to calculate the critical exponents $\alpha$, $\beta$, $\gamma$, $\eta$, $\nu$ and $\delta$. They can be expressed via $\gamma_{m^{2}}^{*}\equiv\gamma_{m^{2}}(g_{1}^{*},g_{2}^{*})$ and $\gamma_{\varphi}^{*}\equiv\gamma_{\varphi}(g_{1}^{*},g_{2}^{*})$ in the following way:

[TABLE]

The critical exponents are related to each other by well-known scaling relations and only two of them may be referred to as independent.

It is instructive to present $\varepsilon$ expansions of cubic fixed point coordinates for physically important case $n=3$. They are as follows:

[TABLE]

[TABLE]

where $\zeta(3,5)$ is double zeta value [35]:

[TABLE]

To give an idea about the numerical structure of these expansions we present them also with the coefficients in decimals:

[TABLE]

The character of a fixed point and, in particular, its stability is determined by the eigenvalues $\omega_{1}$, $\omega_{2}$ of the matrix

[TABLE]

taken at $g_{1}=g_{1}^{*}$, $g_{2}=g_{2}^{*}$. If both eigenvalues are positive the fixed point is stable and describes true critical behavior. At the same time, the roles of $\omega_{1}$ and $\omega_{2}$ in governing the cubic critical behavior are quite different. The eigenvalue $\omega_{1}$ determines the rate of flow to the cubic fixed point along the radial direction in the plane $(g_{1},g_{2})$, while $\omega_{2}$ controls approaching this point normally to the radial ray. In particular, when $n\to n_{c}$ the cubic fixed point tends to coincide with Heisenberg one and $\omega_{2}$ goes to zero. So, the dependence of $\omega_{2}$ on $n$ and its numerical value at $n=3$ are essential in the problem we study. That is why here we write down the $\varepsilon$ expansion for $\omega_{2}$ only. It reads:

[TABLE]

where coefficients $C^{(k)}_{\omega_{2}}$, along with those for $\omega_{1}$, are presented in Tables 7 and 8 of Supplementary materials (see A).

With $\varepsilon$ expansion for $\omega_{2}$ in hand we can find $\varepsilon$ series for the marginal dimensionality of the fluctuating field $n_{c}$. It may be extracted from the equation

[TABLE]

Solving it by iterations in $\varepsilon$ we obtain:

[TABLE]

or, in decimals,

[TABLE]

Six-loop $\varepsilon$ expansions for critical exponents $\eta$ and $\nu$ corresponding to the cubic class of universality result directly from those for anomalous dimensions and scaling relations (16). In its turn, six-loop $\varepsilon$ expansions for $\gamma_{\varphi}$ and $\gamma_{m^{2}}$ originate from RG series (13) and $\varepsilon$ expansions for the cubic fixed point coordinates. Since $\varepsilon$ expansions for the critical exponents under arbitrary $n$ are extremely lengthy they are presented in Tables 9 and 10 of Supplementary materials (see A). Here we write down them only for physically interesting case $n=3$:

[TABLE]

[TABLE]

Of significant interest is also the critical exponent of susceptibility $\gamma$ which is usually measured in experiments and extracted from lattice calculations. Coefficients of its $\varepsilon$ expansion at the cubic fixed point under arbitrary $n$ are presented in Table 11 of Supplementary materials (see A). For $n=3$ this expansion is as follows:

[TABLE]

All calculated $\varepsilon$ expansions are rather complicated and need to be checked up. We compared them with known five-loop series [20] and found complete agreement. In the Ising ($g_{1}\to 0$) and Heisenberg ($g_{2}\to 0$) limits our $\varepsilon$ expansions are found to reduce to their counterparts for $O(n)$-symmetric model [35] under $n=1$ and arbitrary $n$ respectively. Our $\varepsilon$ expansions should also obey some exact relations appropriate to the cubic model with $n=2$. Such a system possesses a specific symmetry: if the field $\varphi_{\alpha}$ undergoes the transformation

[TABLE]

the coupling constants are also transformed:

[TABLE]

but the structure of the action itself remains unchanged [1]. Since the RG functions are completely determined by the structure of the action, the RG equations should be invariant with respect to any transformation conserving this structure [10]. It means that under the transformation (30) the $\beta$ functions should transform in an analogous way while all the observables including critical exponents should be invariant with respect to above replacement (see [10, 29, 43] for details and extra examples). The expansions (12) and (13) do satisfy these symmetry requirements. Moreover, transformation (30) converts the Ising fixed point into cubic one and vice versa making them dual under $n=2$. Six-loop $\varepsilon$ expansions (15) reproduce this duality.

4 Resummation and numerical estimates

With six-loop $\varepsilon$ expansions in hand we can obtain advanced numerical estimates for all the quantities of interest. It is well known that $\varepsilon$ expansions as other field-theoretical perturbative series are divergent and for getting proper numerical results some resummation procedures have to be applied. In this paper we address the methods of resummation based upon Padé approximants [L/M] which are the ratios of polynomials of orders $L$ (numerator) and $M$ (denominator) and Borel-Leroy transformation. The Padé-Borel-Leroy technique enables one to optimize the resummation procedure by tuning the shift parameter $b$ and proved to yield accurate numerical estimates for basic models of phase transitions. Much simpler Padé technique that is certainly less powerful will be also used, mainly in order to clear up to what extent the numerical results depend on the resummation procedure. Note that both approaches do not require a knowledge of higher-order (Lipatov’s) asymptotics of the $\varepsilon$ expansions coefficients finding of which is a separate non-trivial problem.

4.1 Resummation strategy and error estimation

Application of Padé approximants and use of Padé-Borel-Leroy resummation technique are rather straightforward and were described in detail in a good number of papers and books. At the same time, the determination of the final estimate of the quantity to be found and evaluation of corresponding error bar (apparent accuracy) are somewhat ambiguous procedures. The point is that the choice of a subset of approximants which can be accepted as working and used to get the asymptotic or averaged estimate of a given order usually may be tuned within a very wide range what may lead to unreliable (unstable) results and overestimation of the accuracy.

Here we suggest clear and consistent strategy for calculating estimates with Padé approximants and Padé-Borel-Leroy technique which is aimed to yield the stable results and reasonable error estimates from order to order. While finding numerical values of physical quantities with Padé approximants we use the following procedure. To estimate the value in $k$-th order of perturbation theory we take into consideration approximants of $k$ and $k-1$ orders (particular values of $[L/M]$ depend on the observable). The reason of accounting for such a subset is to provide the results stable from order to order while keeping the contribution from $k$-th order dominant. From this set of approximants we exclude ”maximally off-diagonal” ones, in particular $[0/M]$ and $[L/0]$ as they are known to possess bad approximating properties. We exclude also approximants which have poles in the interval $\varepsilon\in[0,2\varepsilon_{phys}]$ (in our case $\varepsilon_{phys}=1$). The reason for this is as follows: if there is a pole in $\varepsilon\in[0,\varepsilon_{phys}]$ the approximant simply cannot be used to estimate the value at $\varepsilon_{phys}=1$, but even if the pole lying outside this area is still close to $\varepsilon_{phys}=1$ such an approximant cannot give reliable estimate as unphysical pole contribution dominates in this case. Particular choice of the upper bound ($2\varepsilon_{phys}$), namely multiplier 2 is based on our experience and tries to keep a balance between dropping out unsuitable approximants and keeping a total number of working approximants as large as possible.

To estimate the error bar (apparent accuracy) we consider values given by different approximants as ”independent measurements” of the quantity and use $t$-distribution $t_{p,n}$ with $p=0.95$ confidence level, i. e. estimates for the value itself and its error are calculated with the following formulas:

[TABLE]

In the case of Padé-Borel-Leroy resummation the procedure is almost the same except the fact that we have an additional – tuning – parameter $b$. For each particular value of $b$ we perform Borel-Leroy transformation of the original series, construct Padé approximants of $k$ and $k-1$ orders for Borel-Leroy transform and drop out approximants $[0/M]$, $[L/0]$ and those spoiled by pole(s) on positive real axis. To find the optimal value of $b$ we perform discrete scan over $b\in[0,20]$ with $\Delta b=0.01$ and search for the value of $b$ which minimizes the standard deviation. The final estimate and error bar are then calculated with (31) for this value of $b$.

4.2 Marginal field dimensionality $n_{c}$

Let us start from the estimation of the fluctuating field marginal dimensionality $n_{c}$. As seen from (25) $\varepsilon$ expansion for $n_{c}$ is alternating and its coefficients rapidly grow in modulo. The former property makes employing Padé approximants not meaningless. The results of Padé resummation of the series (25) under the physical value $\varepsilon=1$ are shown in Table 1.

Applying the procedure described in section 4.1 to the data collected in Table 1 we obtain $n_{c}^{(4)}=2.9\pm 0.4$, $n_{c}^{(5)}=2.94\pm 0.12$ and $n_{c}^{(6)}=2.89\pm 0.14$ as the four-loop, five-loop and six-loop estimates respectively. These estimates are seen to converge to the value close to 2.9 but the rate of convergence and the accuracy are certainly very low.

Since higher-order coefficients of the $\varepsilon$ expansion for $n_{c}$ are big and rapidly grow use of Borel-Leroy transformation that factorially weakens such a growth should significantly accelerate the convergence and refine the estimate itself. This transformation looks as follows

[TABLE]

Padé-Borel-Leroy resummation procedure consists of transformation (32) and analytical continuation of the Borel transform $F(y)$ by means of Padé approximants. It includes also the choice (tuning) of the shift parameter $b$ enabling one to achieve the fastest convergence of the iteration scheme. The results of the Padé-Borel-Leroy resummation of the six-loop series for $n_{c}$ are presented in Fig. 2 and Table 2. The figure shows the behavior of relevant six-loop and five-loop Padé-Borel-Leroy estimates as functions of the parameter $b$ and illustrates, in particular, the emergence of the optimal value $b_{opt}$. Note that the curves in Fig. 2 are drawn only within the regions where Padé approximants of the Borel-Leroy transform have no positive axis poles. Padé-Borel-Leroy estimates of various approximants obtained under the optimal value of $b$ which was found to be $b_{opt}=1.845$ are collected in Table 2.

As is seen the application of Padé-Borel-Leroy machinery indeed makes the iteration faster convergent and corresponding estimates much less oscillating. Being processed according to our strategy (Section 4.1) the numbers presented in Table 2 give us $n_{c}^{(4)}=2.96\pm 0.11$, $n_{c}^{(5)}=2.91\pm 0.03$ and $n_{c}^{(6)}=2.915\pm 0.003$ at the four-, five- and six-loop levels. The last, highest-order value

[TABLE]

we accept as a final result of our calculations.

4.3 Critical exponents

Since the coordinates of the fixed points depend on the normalization conditions adopted their numerical values being non-universal are not interesting from the physical point of view. That is why further we proceed directly to evaluation of critical exponents characterizing the cubic class of universality at $n=3$. Starting from the six-loop $\varepsilon$ expansions for $\eta$ and $\nu^{-1}$ and using well-known scaling relation we obtain $\varepsilon$ expansions for exponents $\alpha$, $\beta$, $\gamma$, $\nu$ and $\delta$. Then we perform Padé and Padé-Borel-Leroy resummation of all the series in hand. As the Padé-Borel-Leroy resummation procedure turns out to be most effective (regular and fast convergent) for $\beta$ and $\gamma$ we present here details of evaluation of these two exponents. Numerical values of $\beta$ and $\gamma$ obtained within Padé and Padé-Borel-Leroy resummation approaches are collected in Tables 3, 4, 5 and 6.

Similar tables were calculated for the exponents $\alpha$, $\delta$, $\eta$ and $\nu$. All the final estimates and error bars obtained according to the resummation procedure suggested in Section 4.1 are presented in Table 7.

What is the accuracy of numerical results just found? Some idea on the point may be obtained looking at the differences between the Padé-Borel-Leroy and Padé estimates presented in Table 7. However, much more definite conclusions concerning an actual accuracy of our calculations can be made on the base of the analysis to what extent the numbers obtained obey exact scaling relations between the critical exponents. One can combine critical exponents in different ways. We choose the next set of independent relations:

[TABLE]

that are ”normalized to unity” to get the estimates of accuracy more uniform. Since the calculated values of critical exponents are approximate they can not meet the scaling relations precisely and emerging discrepancies may be considered as a measure of achieved accuracy. The discrepancies relevant to scaling relations (34) along with their error bars originating from the estimates of the critical exponents themselves (Table 7, upper line) are presented in Table 8.

As is seen the deviations from exact scaling relations are small demonstrating the consistency of our approach and indicating that actual computational uncertainty of found numerical estimates is of order of 0.01.

To finalize this section, in Table 9 we present, for completeness, the values of correction-to-scaling exponents $\omega_{1}$ and $\omega_{2}$ obtained by resummation of corresponding $\varepsilon$ expansions. Despite the fact that zero lies inside the error bar for $\omega_{2}$ the median value of this exponent, being very small, turns out to be positive. Moreover, keeping in mind the results of independent evaluation of $n_{c}$ we may state that the value of $\omega_{2}$ given by six-loop $\varepsilon$ expansion analysis is certainly positive. More accurate estimates for $\omega_{2}$ can be obtained within the higher-order (seven-loop, etc.) approximations or by means of more sophisticated resummation procedure such as Borel transformation combined with conformal mapping which will be a subject of a separate paper.

5 Discussion

In this section we will compare our results with those obtained earlier within the lower-order approximations and by alternative methods.

The first quantity of interest is the marginal spin dimensionality for which we get the value $n_{c}=2.915(3)$. It is worthy to note that the $\varepsilon$ expansion for this quantity has rapidly growing coefficients (see eq. (25)) what prevents Padé approximants from giving accurate enough numerical results while Padé-Borel-Leroy approach yields stable estimates with an accuracy increasing from order to order. The results of previous studies performed within the $\varepsilon$ expansion approach and RG machinery in fixed dimensions (3D RG) as well as the numbers extracted from the Monte Carlo simulations and the six-loop pseudo-$\varepsilon$ expansion are aggregated in the Table 10.

In addition, the values of $n_{c}$ collected in Table 10 are depicted at Fig. 3 to visualize the trend these values demonstrate under increasing order of approximation. This trend enable us to conclude that $n_{c}$ is certainly less than 3 for the 3D cubic model that justifies the significance of studying the cubic class of universality.

The other quantities of prime physical importance are critical exponents of the cubic universality class. We should stress that to get estimates for critical exponents we perform resummation of the series for each exponent separately and afterwards checked a validity of several scaling relations (34). Despite the fact that sometimes the relations are satisfied with inaccuracies exceeding corresponding error bar estimates, these deviations are not too large lying within 3$\sigma$ interval. This may be considered as a proof of the consistency of the results obtained and a demonstration of the numerical power of the $\varepsilon$ expansion approach.

It is worthy to compare our estimates with their analogs given by the lower-order approximations and with the results of multi-loop 3D RG analysis. The data enabling one to do such a comparison are collected in Table 11. The numbers presented in both columns are seen to rapidly converge to the asymptotic values that differ from each other only tiny coinciding in fact within the declared error bars. It confirms the conclusion that the field theory is a powerful instrument enabling one to get precise numerical results provided the calculations are performed in high enough pertubative order. On the other hand, addressing the six-loop $\varepsilon$ approximation shifts the estimates only slightly indicating that they should be very close to the exact values still unknown.

Another point to be discussed is to what extent – quantitatively – the critical exponents of the cubic class of universality differ from those of the 3D Heisenberg model. Since for $n=3$ the cubic fixed point lies near the Heisenberg one corresponding differences are known to be rather small. In Table 12 we present the estimates of critical exponents for cubic and Heisenberg classes of universality obtained in the six-loop approximation. As expected, the differences between numerical values of critical exponents for these two classes are really small. So, it is hardly believed that measuring critical exponents in physical or computer experiments one can distinguish between cubic and Heisenberg critical behaviors.

Conclusion

To summarize, we performed six-loop RG analysis of the critical behavior of $n$-vector $\varphi^{4}$ model with cubic anisotropy in the framework of $\varepsilon$ expansion approach employing the minimal subtraction scheme. We calculated $\varepsilon$ expansions for marginal spin dimensionality $n_{c}$ and critical exponents $\alpha$, $\beta$, $\gamma$, $\delta$, $\eta$, $\nu$, $\omega_{1}$, $\omega_{2}$ for the cubic class of universality. We resummed these diverging series with Padé approximants and using Padé-Borel-Leroy technique. Obtained numerical estimates for critical exponents turn out to be self-consistent in the sense that they are in accord, within the computational uncertainties, with the scaling relations. Six-loop contributions are found to shift five-loop estimates only slightly but they improve numerical results considerably diminishing their error bars. Our results confirm and strengthen the conclusion that cubic ferromagnets (D = 3, n = 3) belong to cubic class of universality and their critical behavior is described by critical exponents differing from those of 3D Heisenberg model. At the same time, the critical exponents of 3D cubic and Heisenberg models are numerically so close to each other that it makes their behaviors practically indistinguishable if one limits himself by measuring critical exponents only.

Acknowledgment

It is a pleasure to thank Professor M. Hnatič and Professor M.Yu. Nalimov for fruitful discussions. E.I. and A.K. are especially grateful to the Professor Hnatič for support and hospitality during their stay in Slovakia. This work has been supported by Foundation for the Advancement of Theoretical Physics ”BASIS” (grant 18-1-2-43-1).

Appendix A Supplementary materials

In Supplementary materials we present expansions of RG functions and critical exponents for arbitrary $n$. In rg_expansion_coefficients.pdf we list coefficients $C_{k}^{i,\,j}$ of the expansions of beta functions $\beta_{1}$, $\beta_{2}$ (12), anomalous dimensions $\gamma_{\phi}$, $\gamma_{m^{2}}$ (13) and $\varepsilon$ expansions of coordinates of the cubic fixed point (15), correction-to-scaling exponents $\omega_{1}$, $\omega_{2}$ (22) and critical exponents $\eta$ (26), $1/\nu$ (27) and $\gamma$ (28) corresponding to cubic universality class.

Additionally, for RG functions ($\beta_{1}(g_{1},g_{2})$, $\beta_{2}(g_{1},g_{2})$, $\gamma_{\phi}(g_{1},g_{2})$, $\gamma_{m^{2}}(g_{1},g_{2})$) we provide Mathematica file with their expansions (rg_expansion.m). For critical exponents we present Mathematica files for all non-trivial fixed points: cubic (cubic_crit_exp.m), Ising (ising_crit_exp.m) and Heisenberg (heisenberg_crit_exp.m). Each file contains $\varepsilon$ expansion for exponents $\alpha$, $\beta$, $\gamma$, $\delta$, $\eta$, $\nu$ as well as for $1/\nu$, correction-to-scaling exponents $\omega_{1}$, $\omega_{2}$ and coordinates of fixed points $g_{1}^{*}$, $g_{2}^{*}$. In the file corresponding to cubic fixed point (cubic_crit_exp.m) we also present expansion for marginal spin dimensionality $n_{c}$.