Critical behavior of isotropic systems with strong dipole-dipole interaction: three-loop study

We analyze the critical behavior of isotropic systems with dipole-dipole interaction by renormalization-group methods in fixed space-time dimensions. Working in three-dimensional theory we analytically find three-loop expressions for critical exponents in the limit of dominating dipole-dipole forces. Resummation of the series obtained provides numerical values close to $O(3)$-theory predictions, justifying the applicability of such a simplified model to systems with strong dipole-dipole interaction.


Introduction
The role of the O(n)-symmetric model is difficult to overestimate in studying the critical behavior of a wide range of physical systems. To process this model correctly, Wilson proposed the elegant idea to apply the renormalization group (RG) approach and developed the ε expansion technique as a computational procedure [1,2,3,4]. The five-loop calculation performed almost forty years ago was considered as a record one [5,6,7,8], and the situation remained unchanged for a long time.
Due to the impressive development of numerical methods [9,10,11,12,13,14,15,16,17,18], as well as computational power of modern hardware, the six-and sevenloop results in ε were obtained in Refs. [19] and [20], respectively. The found values of the diagrams allowed to analyze the critical behaviour in the following orders of perturbation theory in the case of different O(n)-symmetric generalizations, which specify symmetries of various physical systems [21,22,23,24].
As known, there are various RG approaches. In addition to the ε expansion, it is possible to construct a theory directly in physical spatial dimensions. To this day, for the three-dimensional (3d) O(n)-symmetric model the record -six-loop -calculation was performed by Nickel et al. in Ref. [25,26]. That made it possible to find the RG functions and accurate values of the critical exponents [27,28].
The real-space RG can be advantageous when the action has a specific term. For example, the three-dimensional approach is extremely natural when the action contains a term that requires the coincidence of order parameter and spatial dimensions. In particular, taking into account dipole-dipole interaction leads to the model containing such a term. We should note that this model has only been studied in the low-order approximation within the ε expansion [29,30] as well as with the help of RG in physical spatial dimension [31,32]. Since calculations for the systems with dipole-dipole interaction is quite complicated, but for a number of physical systems the inclusion of this interaction is important [33,34,35,36], we aim to estimate this effect by considering a simplified model that is a limiting case of complete dominance of dipole forces. To do this, we obtain the three-loop RG series, on their basis we extract the values of critical exponents. Comparison of the obtained numbers with their counterparts in the case of the standard O(n)-symmetric universality class from [27] allows us to evaluate the need to consider the term responsible for dipole-dipole interaction in the analysis of critical behavior.
As a starting point, we reproduce the results for the simple O(n)-symmetric model (SSM), which is presented in Ref. [26,25,37]. We managed to get answers up to three loops completely analytically. The technique for calculating diagrammatic series is further applied to the case of the model with strong dipole-dipole interaction (SDM). All the RG expansions and final numerical estimates of observables are presented both for SSM and SDM.
The paper is organized as follows. In Sec. 2, the model and the renormalization scheme, which we use in this work, are described. After that, the calculation details, including explanation of diagrammatic technique, are presented in Sec. 3. Next, in Sec. 4, the expansions for RG functions and critical exponents as well as their numerical estimates are shown. Finally, in Sec. 5 we draw a conclusion.

Model and renormalization scheme
Thus, we are interested to analyze the critical behaviour of the following model: where ϕ α is d-component bare field, λ is bare coupling constant, m 2 0 is bare mass being proportional to T − T c , where T c is the mean-field critical temperature. Tensor Ω αβ is defined as follows: The modification of tensor Ω αβ in the case of SDM arises from the Fourier transform of dipolar interaction term ∼ (x 2 δ αβ − x α x β )/x 5 . Tensor factor T αβγδ reads as: This model is studied directly in d = 3. The propagator of the theory can be written in the following form: Our calculation scheme is based on the one used in [25] 1 and can be explained on the examples from the SSM model, results of calculation are in one to one correspondence with the three-loop results presented in [25]. For the SDM model we proceed in a similar way up to the redefinition of Ω αβ .
Our final goal is to determine renormalization constants Z i (u) as a function of the renormalized dimensionless coupling u. However, for actual diagrams calculation it is convenient to introduce an intermediate scheme, where the propagator is defined by (4) and perturbative series are organized as expansion in auxiliary parameter W. For dimensionful quantities, such as four-point function in d = 3, result of calculation depends on m 0 and may diverge due to divergent subdiagrams. Such divergences are removed by a specially constructed mass counter-term and after eliminating bare mass with m 2 0 = Z −1 3 m 2 − δm 2 0 the final result becomes finite. The exact form of δm 2 can be determined from the mass derivative of the two-point function and depends on the chosen regularisation prescription. In original work [25] divergent subgraphs were subtracted under the integral sign making sum of the diagram and diagram with counter-term insertion finite. In our work, we make use of dimensional regularization, thus both diagram and its counter-part now have poles in dimensional regularization parameter ε, but the sum is finite and the same as results presented in [25].
The specific choice of connection m 2 0 = Z −1 3 m 2 − δm 2 0 between bare and renormalized mass greatly simplifies calculation allowing independent extraction of renormalization constants in the intermediate scheme as series in parameter W = − √ Z 3 λ/m from particular diagrams. Similar to work [25], we are interested in the set of renormalization constants Z i (W) in the intermediate scheme defined as follows 2 : Here, Especially, in our case we calculate the four-point function with zero external momenta in (5), two-point function with O 2 insertion and zero external momenta in (6) and zero external momentum limit of the momentum derivative of the two-point function with the appropriate projector applied in (7).
Results for renormalization constants as functions of the renormalized dimensionless coupling Z i (u) can be derived from (5), (6), and (7) with: where we used the relation λ = umZ 1 /Z 2 3 , which binds the bare and renormalized couplings.
Having obtained all the necessary Z i (u), we can derive the expressions for all the RG functions. To obtain the series for β-function, we can utilize the following expression: The anomalous dimensions are connected with the renormalization constants in the following way: As well known, the behavior of the system in the vicinity of the critical temperature is determined by infrared-stable fixed point u * , which is determined by zero of the βfunction (β(u * ) = 0). Also, it is convenient to introduce other functions of u, in terms of which one could express the critical exponents: We limit ourselves by consideration of the most traditional critical exponents η, ν, and γ, which are defined by functions from (12) evaluated at u = u * : η = η(u * ), ν = ν(u * ), and γ = γ(u * ). In addition, it is necessary to define the leading corrections to the scaling law, which is governed by the following exponent: However, before we turn to the numerical results, let us tell a few words about the technique for calculating the diagrammatic expansions.

Details of calculation
As was mentioned in Sec. 2, we calculate individual diagrams in the framework of dimensional regularization as expansion in parameter ε = 3 − d. Compared to a more traditional numerical evaluation technique used in [25], the chosen approach at the first sight seems overcomplicated, but it has several attractive features. First of 4 all, working at the three-loop level we are able to perform all calculations analytically, which is especially useful for accurate divergences subtraction and manipulation with the obtained expressions. Secondly, the dimensional regularization framework allows to apply efficient tools for modern multi-loop calculations. Also, the availability of three-loop results for dimensionally regulated massive tadpole [38] integrals justifies our choice of the regularization scheme.
In the three-loop calculation, presented in the paper, a number of diagrams to be considered is not so high, especially compared to six-loop calculation [25]. However, we decided to develop a highly automized setup to simplify tensor manipulations, especially in the case of SDM, and to eliminate possible errors. One more complication comes from the transverse structure of the propagator in SDM (4) leading to the appearance of new classes of diagrams with massless propagators. Calculation of such diagrams cannot be reduced to the set of integrals provided in [25] and they need a separate treatment.
We start our chain of calculations with DIANA [39], which internally calls QGRAF [40] to generate diagrams in both models through the single run. After substitution of Feynman rules for a specific model and application of appropriate projectors, we perform partial fraction decomposition and map the obtained scalar integrals on predefined set of topologies. Inside each topology all integrals are reduced to a small number of master integrals during the reduction step performed with the Laporta algorithm implemented in package FIRE [41] with the inclusion of additional symmetry rules between integrals generated with package LiteRed [42].
After reduction to master integrals, in order to obtain results of the individual diagrams as the series in ε, we substitute the results of the ε expansions for master integrals. It is important to note, that minimal set of master integrals obtained after reduction step does not necessarily match original set of diagrams. Some of needed master integrals may contain poles in ε, or vice versa have divergent coefficients. In the latter case we need to provide higher orders of expansion in ε of these master integrals. For all the most complicated integrals we need only finite parts available from [38,43,44] and for some higher order expansion results can be found in [45,46].
Results for individual diagrams contributions both in the SSM and SDM can be found in Appendix B.1 for the two-point functions and in Appendix B.2 for the fourpoint function and two-point function with the operator insertion. Present tables contain contribution to the expansion in variable W and explicitly reads: where X = O for the case of SSM and X = D for the case of SDM. All presented results are finite, since for the divergent four-point diagram (no.7 in Appendix B.2) we present the result with the subdivergence subtracted. For SSM the contributions to Z 1 from (B.31) and to Z 4 from (B.33) coincide with [25], where a different 5 subtraction scheme was utilized. One more check on the validity of provided results comes from the fact that at least at the two-loop level we can consider a more general type of diagrams depending on additional parameter g responsible for dipole-dipole interaction. The two limiting cases, g → 0 and g → ∞, correspond to the result in SSM and SDM, respectively. As an example for the diagram covering both Γ

Numerical results
Following pioneering works [26,28], we modify the normalization of the coupling constant as v = 11u, which is dictated by the one-loop coefficient of β-function in the case of SSM when n = 3 (O(3)-symmetric model). Moreover, in order to compare the results in the ordinary (labeled as A O ) and dipole-dipole (labeled as A D ) cases, at all stages we present the numerical results for both of them. All the expansions are found analytically. Due to their bulkiness, they are presented in Appendix A. Here, we restrict ourselves to the results through decimal notation. The number of decimal places was chosen to be the same as for Nickel et al. in work [37] despite the fact that they are known with absolute accuracy. Thus, the three-loop expansions for the RG functions have the following form: Let us extract now the numerical values of the critical exponents. For this, as already mentioned, we need to find the value of the coordinate of fixed point v * . Due to the asymptotic nature of renormalization group series, in order to obtain the proper numerical results the various resummation techniques should be applied. The most basic of them is the method of Padé approximants. However, in such an approximation, this method applied to the initial RG expansion does not give any reliable results -all the Padé approximants give drastically different values. The Padé-Borel (PB) method 6  (6) applied to the original RG series does not improve the situation considerably. However, this dramatic picture can be improved in the following way. B. Nickel (see Ref. [19] in Ref. [27]) proposed a method to advance the convergence of series. This method consists of reexpanding of the renormalization group series into alternative power expansions whose coefficients demonstrate a more favorable behavior. The idea is to put a formally small parameter τ at the linear term of the β-function (β(v) = −τv + v 2 + . . . ). Then, similarly to how it is done within the ε expansion, one can iteratively find a fixed point coordinate as series in τ. After further processing with various resummation techniques of new expansion, which already has a more acceptable structure, the formal parameter τ should be equated to unity. Following this recipe for fixed point coordinates in both cases we obtain: Based on these expansions we extract the numerical estimates for coordinates v * O and v * D by means of naive direct summation, Padé approximants, as well as PB technique. These numbers are presented in Table 1. We choose diagonal approximant [2/1] for the Padé and PB estimates, given that approximants should be constructed for series starting with a constant. Noteworthy, as a test, we have the six-loop fixed-point estimate in the SSM case [28,47]: 1.391 (1). Believing in the convergence of estimates with increasing orders of perturbation theory, we can choose its difference from the six-loop value as an error for the final three-loop estimate: 1.394 (4). Following this logic, we can try to interpolate the calculation error per estimate in the case of a model with a dipole-dipole interaction: 2.014 (6). Let us move on to the analysis of the most interesting quantities from the physical point of view -the critical exponents. The situation with exponents is similar. The analysis of the initial RG series leaves much to be desired. To extract some proper estimates, we reexpand the RG expansions for functions η(v), ν(v), γ(v), and ω(v) in Table 2: Comparison of numerical estimates of critical exponents for SSM and SDM. The numbers were extracted based on the τ-series. Padé and PB estimates are obtained according to the near-diagonal 3-loop values ([2/1] and [1/2] since the diagonal is absent in this case) and according to the diagonal two-loop [1/1]. The average value based on the Padé and PB estimates is given as an answer. The dash corresponds to the fact that some of the approximants in the sample are spoiled by the presence of a dangerous pole, which does not allow to make an estimate.

Exponent
Direct terms of τ. They read as follows: ν D (τ) = 0.5 + 0.1225490196τ + 0.0499569216τ 2 + 0.0161209071τ 3 + O(τ 4 ), Once we have the truncated series, we can apply to them the summation techniques, which were used previously for the fixed point coordinates. The corresponding numbers are presented in Table 2. As previously, believing in the succession of the model with the strong dipole-dipole interaction compared to the conventional symmetric model, to estimate the final values of the errors, we resorted to the results obtained within the six-loop calculation of the symmetric theory only with the exchange interaction [28] (ν = 0.7054 (11), η = 0.0340 (25), γ = 1.3866 (12), and ω = 0.779 (6)). Based on the results presented in Table 2, we can conclude that the strong dipole-dipole interaction (g → ∞) does not strongly change the universality class regarding the critical exponents. Only the scaling correction exponent has undergone significant changes. Such numerical results, assuming that we will not find any surprises in higher orders, to some extent justify the constantly neglected term in the action, which is responsible for the dipole-dipole interaction.

Conclusion
In this work, we have calculated the three-loop RG expansions for the isotropic field model in the limiting dipole-dipole region within the renormalization group in the physical spatial dimensionality (d = 3). Such a regime is extremely important to take into account when analyzing the critical behavior in a number of ferromagnets and ferroelectrics. All the results have been found completely analytically, which is of particular value for this area. At each step of our calculations, all computational procedures were performed both for the model with dipole-dipole interaction (SDM) and for the usual isotropic model (SSM). The obtained values for the critical exponents, which are the main measure of the difference between two different universality classes, turn out to be very close between the two models -SSM and SDM, except for the correction for the scaling exponent ω. To some extent, this justifies the fact that, when analyzing the critical behavior of such systems, the term responsible for the dipoledipole interaction is often neglected in the action. The results of the diagrams calculation obtained here can be applied to the cases of different symmetries, and the techniques used here can break into higher orders of perturbation theory. Where we have introduced angle α = arcsin 1 3 according to [43] and Clausen function is defined as Cl 2 (φ) = Im Li 2 e iφ .