Six-loop $\varepsilon$ expansion study of three-dimensional $O(n)\times O(m)$ spin models

The Landau-Wilson field theory with $O(n)\times O(m)$ symmetry which describes the critical thermodynamics of frustrated spin systems with noncollinear and noncoplanar ordering is analyzed in $4 - \varepsilon$ dimensions within the minimal subtraction scheme in the six-loop approximation. The $\varepsilon$ expansions for marginal dimensionalities of the order parameter $n^H(m,4-\varepsilon)$, $n^-(m,4-\varepsilon)$, $n^+(m,4-\varepsilon)$ separating different regimes of critical behavior are extended up to $\varepsilon^5$ terms. Concrete series with coefficients in decimals are presented for $m=\{2, \dots, 6\}$. The \textit{diagram of stability} of nontrivial fixed points, including the chiral one, in $(m,n)$ plane is constructed by means of summing up of corresponding $\varepsilon$ expansions using various resummation techniques. Numerical estimates of the chiral critical exponents for several couples $\{m,n\}$ are also found. Comparative analysis of our results with their counterparts obtained earlier within the lower-order approximations and by means of alternative approaches is performed. It is confirmed, in particular, that in physically interesting cases $n=2, m=2$ and $n=2, m=3$ phase transitions into chiral phases should be first-order.


Introduction
Compared to conventional (anti)ferromagnets, where ordering occurs in a collinear manner, phase transitions in systems with frustration are much less studied. The source of such frustration can be the geometric features of crystal lattice as well as the nature of the interactions between nearest and next-nearest spins. In case the number of spin components is greater or equal to two (n 2) the presence of frustration leads to a noncollinear and a noncoplanar spin ordering. As prime examples of systems with mentioned non-trivial spin structures stacked triangular antiferromagnets (STA) and helical magnets (HM) may be considered. In the course of studying the critical properties of such systems, the main question arises: which type of a phase transition takes place in real materials -whether it has to be a continuous one with non-standard set of critical exponents corresponding to so-called chiral universality class or a phase transition should be first-order.
There are many experimental works dedicated to the frustrated spin systems where they were studied by means of optical second-harmonic spectroscopy, measurements of magnetization, susceptibility and specific heat, resistivity measurements, synchrotron X-rays and neutron diffraction, Mössbauer technique etc. [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]. The results of these experiments turn out to be rather contradictory. For materials with the same symmetry the critical exponents measured are noticeably scattered, and even the order of the phase transition was found to be different for the closely related substances.
As for theoretical studies, there are a good number of works investigating the critical properties of frustrated magnets. Here we note only some of key points in the development of the study with emphasis on a field theoretical (FT) approach; those interested in more details and the history of the problem may be addressed to comprehensive reviews [24,25,26,27,28,29].
The first attempts to analyze the critical behavior of the frustrated spin systems by means of renormalization group (RG) approach were undertaken about four decades ago [30,31,32]. These works were intended to explain the critical properties of rare-earth magnets such as Dy, Ho and Tb. It was realized that for proper description of the systems with nontrivial ordering within FT approach an account for the O(n)-symmetric self-interaction is not sufficient. To cope this problem the Landau-Wilson (LW) action with O(m) × O(n) symmetry was suggested which is written down in modern notations, say, in (1). The presence of such a specific symmetry being natural for real systems was shown to lead to the critical behavior which is different from that described by the standard (O(n)-symmetric) universality class [33,34,35].
The reasonableness of studying certain classes of universality depends on whether the corresponding regimes of critical behavior are stable or not. The stability in turn is determined by the structure of the renormalization group (RG) flows. For physically interesting values of m and n relevant phase portraits were analyzed. In the case m 3 the results given by all FT approaches ε expansion [36,37,38], 1/n expansion [37,39], 3D RG technique [40], 2 − ε machinery [41,42], pseudo-ε expansion [38], "exact" (functional, non-perturbative) renormalization group (ERG) [43] -are in favor of existence of some upper bound (marginal) value of n, n + , such that for n > n + the system demonstrates continuous phase transition with specific set of critical exponents while for n < n + the transition is first-order.
For m = 2 the situation is more complicated 1 [37,39,41,42,44,45,46,47,48,49,50,51,52,53,54]. For the same values of n the results obtained by different approaches are controversial [55,56,57]. Some materials with identical symmetry, being expected to belong to the same universality class described by action (1) under physical values of n, demonstrate considerable scattering of critical exponents measured in different samples [58,59,60,61,62]. An interesting picture was obtained within 3D RG analysis: there is the value N c2 < n + such that the systems with n < N c2 have to demonstrate the scaling behavior [49]. According to the six-loop 3D RG analysis the physical values n = 2, 3 are covered by the last case [47,48]. The critical exponents obtained within 3D RG approach [47,63] are in a bad agreement with experimental results. Regarding ERG there were found no stable fixed points within this approach signaling about realization of first-order phase transition [43,55,57]. Contradictory results are observed not only within field-theoretical methods. For STA Monte Carlo calculations performed by different groups argues both in favor of continuous phase transition [64,65,66,67,68] and first-order one [69,70,71,72]. It is worth noting here that systems can undergo first-order phase transition even if some stable fixed points exist on RG flows -this occurs when initial (bare) values of couplings lie outside their regions of attraction. 1 For convenience of reader notations of physically interesting quantities were chosen in a similar way as in [37,38].

2
The realization of this work has become particularly relevant due to the recently obtained sixloop RG expansions for O(n)-symmetric model [73]. These results allow us to extend obtained fifteen years ago five-loop ε expansions for O(m) × O(n)-symmetric model [38] to the six-loop order. It is expected that numerical estimates obtained within this approximation will be a high priority with respect to determination of accurate values of marginal dimensionalities n + , n − and n H .
The paper is organized as follows. In Section 2 the model and renormalization scheme employed are described. In Section 3 all the quantities of interest including RG functions, fixed points and marginal dimensionalities are calculated. The numerical estimates of n + (m, 3), n − (m, 3) and n H (m, 3) for m = {2, . . . , 6} are presented in Section 4. There are also the results concerning the critical exponents for physically interesting couples {m, n}. In Section 5 the numbers obtained in this work, found earlier within the lower-order approximations and by means of alternative approaches are compared and analyzed. In the last section a summary of the main results is presented.

Model
The critical behavior of chiral systems within FT approach is described by the O(n) × O(m)symmetric LW action with two coupling constants [31,33,74,75,76,77,78]: where ϕ 0αi is m-size set of n-component vector fields (α ∈ {1, . . . m}, i ∈ {1, . . . n}), g 01 and g 02 are bare coupling constants, m 0 is a bare mass being proportional to T − T 0 , where T 0 is mean-field transition temperature. The tensor factors T (1) and T (2) entering the O(mn)-symmetric and chiral terms respectively are as follows αi,β j,γk,δl = 1 3 (δ αi,β j δ γk,δl + δ αi,γk δ β j,δl + δ αi,δl δ γk,β j ), δ αi,β j = δ αβ δ i j , In particular, Assuming n m, to provide a noncollinear ordering and positive definiteness of action (1) it is necessary to impose on the coupling constants the following conditions: 0 < g 02 /g 01 < m/(m − 1) [45,37,38]. In the case of negative g 02 the equation (1) describes the magnets with simple unfrustrated ordering and with sinusoidal spin structure [77,78]. For m = 2 the model (1) corresponds to systems with non-collinear but coplanar ordering. Physically most interesting of them are XY (n = 2) and Heisenberg (n = 3) frustrated antiferromagnets. In case of m 3 the action describes the critical behavior of magnets with non-coplanar ordering [36]. 3 In the critical region, the analysis of RG flows demonstrates the presence of four different fixed points (FPs). Two of them always exist -Heisenberg or O(mn)-symmetric (g 02 = 0) and Gaussian (g 01 = g 02 = 0) ones, while appearing of two others depends on the value of n. Upper left picture corresponds to n < n H , upper right one -to n H < n < n − , lower left picture corresponds to n − < n < n + and lower right one -to n > n + . Symbols in boxes mark Gaussian (G), chiral (C + ), antichiral (C − ), sinusoidal (S + ), antisinusoidal (S − ) and Heisenberg (H) FPs. White areas correspond to regions of attraction of stable fixed points where the system demonstrates a scaling behavior.

Renormalization
The model studied is known to be multiplicatively renormalizable 2 . The bare quantities g 10 , g 20 , m 2 0 , ϕ 0 can be expressed via the renormalized ones g 1 , g 2 , m 2 , ϕ by means of the following relations: In terms of renormalized parameters the action (1) acquires the following form αi,β j,γk,δl ϕ αi ϕ β j ϕ γk ϕ δl , where µ is an arbitrary mass scale introduced to make renormalized couplings g 1 and g 2 dimensionless. Use of RG constants eliminates all the divergences, making renormalized Green functions free of them. The multiplicative renormalizability of the model allows to limit ourselves by removal the divergences only in two-and four-point one-particle irreducible Green functions: αi,β j,γk,δl m(m + 2)n(n + 2) αi,β j,γk,δl m(m + 2)n(n + 2) Γ (4) αi,β j,γk,δl .
In the course of renormalization we address here the Minimal Subtraction (MS) scheme. The renormalization constants in this scheme contain only pole contributions in ε and depend only on ε and coupling constants: 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 Bogoliubov-Parasyuk R operation: where R -incomplete Bogoliubov-Parasyuk R-operation, K -projector of the singular part of the diagram andΓ i -normalized Green functions of the basic theory (see e.g. [80,81]) defined by the following relations: The undoubted advantage of the Bogoliubov-Parasyuk approach is that counterterms of the diagrams computed for O(1)-symmetric (scalar) model can be easily generalized to any theory with different symmetries due to the opportunity to factorize the tensor structures (see e.g. [82,83,84]). To calculate tensor factors for particular diagrams of the O(m) × O(n)-symmetric model (1) one should apply projectors (7) to it. To implement these calculations for sufficiently large number of diagrams, we resort to the effective package of tensor algebra FORM [85] and manipulating graphs package Graphine/GraphState [86] while counterterm values can be taken from data obtained in the course of recent 6-loop calculations for O(n)-symmetric model [73,87,88].

RG functions, fixed points, critical exponents and marginal dimensionalities
In this section we present expressions for RG functions, i.e. β functions and anomalous dimensions γ ϕ , γ m 2 , fixed-point coordinates, critical exponents, correction-to-scaling exponents ω 1 and ω 2 and also calculate marginal dimensionalities which determine what regime of critical behavior is realized for particular values (m, n).
With the series for renormalization constants in hand, the RG functions can be calculated by means of the following relations: where Z (1) i -coefficients at first pole in ε from (8). By using these formulas we obtained the RG functions as series in terms of renormalized couplings up to six-loop order. The expansion coefficients were found analytically. Due to the cumbersomeness of six-loop expansions and the fact that they do not represent themselves any physical interest we shall limit ourselves here by presentation of the first terms. So, β functions and anomalous dimensions in the one-loop approximation read: The full expansions up to six-loop order for β functions and anomalous dimensions γ ϕ , γ m 2 are presented in supplementary materials as Mathematica-file (Appendix B). As was said above the critical behavior of the systems described by O(m) × O(n)-symmetric LW action is governed by one of the four fixed points (g * 1 , g * 2 ) of RG equations that are zeroes of β functions: Solving these equations iteratively, we can find solutions as series in powers of ε. In one-loop approximation for arbitrary values of n and m the ε expansions for the coordinates of chiral (c) and antichiral (ac) FPs are as follows: where R = m 2 − 2m(5n + 2) + n 2 − 4n + 52. As in the case of β functions, the fixed-point coordinates themselves are not essential. The purpose to present here the one-loop expressions for FP coordinates is to give an idea about the restrictions which are imposed on the choice of m and n when calculating critical exponents of the chiral class of universality. In particular, under m 1 parameter n has to satisfy the conditions: This fact being a feature of the ε expansion creates some difficulties when considering values of m interesting from physical point of view. To solve this problem in the following sections we address some trick suggested in [37].
Returning to the goals of this section, in order to characterize the chiral class of universality we present here the definitions of critical exponents α, β, γ, η, ν and δ. They can be expressed via anomalous dimensions γ * in the following way: As is well known the critical exponents are related to each other by well-known scaling relations and only two of them may be referred to as independent.
Whether the fixed point is stable or not depends on the eigenvalues ω 1 , ω 2 of the matrix taken at A fixed point is stable only if both eigenvalues are positive. As was already mentioned the numerical values of the marginal dimensionalities n + , n − and n H are of the highest importance for the problem. Let us formulate the conditions to obtain these quantities. Everywhere below m is considered as a fixed parameter. First, to derive series for n ± it is enough to impose the following conditions in all known orders in ε: det As an alternative way one may require instead of (22) the coincidence of coordinates of the chiral and antichiral FPs. In case of n H the conditions being effective from computational point of view have the following form: Thus it suffices to require zeroing of the second coordinate of the antichiral fixed point.

Numerical results
In this section we analyze the series for all marginal dimensionalities under the physically interesting values m = {2, .., 6}. To get the proper numerical estimates from these expansions being in fact asymptotic we have to apply various resummation techniques among which are simple method of Padé approximants and those based upon the Borel transformation. We pay special attention to n + (m = {2, 3}, 4−ε) because these quantities are of prime physical importance.
To give an idea about the numerical structure of corresponding expansions for mentioned values of m we present them with the coefficients in decimals in Table 1. In order to obtain three- dimensional estimates one has to evaluate these series for ε = 1. The dramatic rate of the coefficient's growth makes the direct summation (just equate ε to unity) absolutely fruitless. The strategy of treatment of these divergent series, i. e. the technique of their resummation will be demonstrated in next sections for specific values of m. For m = 2, 3 the coefficients of relevant expansions up to forth order in ε coincide with those found earlier in other works [37,44,38]. Because of the cumbersomeness of analytic expressions for series coefficients, they were put into the Supplementary Materials (Appendix B). As seen from (24), the series for "lower" marginal dimensionalities n − (2, 4 − ε) and n H (2, 4 − ε), being divergent, possess regular structure that may be considered as favorable for their resummation. This allowed us to obtain numerical estimates n − (2, 3) = 1.970(3) and n H (2, 3) = 1.462 (13) stable with respect to addressing different resummation techniques. Regarding n − (2, 3) it is necessary to make some comment. Despite the fact that new estimate is greater than previous one 1.968(1) [38] it can not still be referred to as accurate enough because of the fact that in case n = m = 2 there is a mapping onto the tetragonal model [45] where nontrivial FPs have to exist. Highly likely, this shortcut reflects bad convergence of the ε expansion approach in three dimensions.
Let us consider further the series for n + , which possess the worst structure among others, but corresponding numerical estimates determine the type of phase transition in real systems (n = 2, 3) and therefore play a crucial role. First, we use for resummation simple method based upon Padé approximants. Standard step within this approach is a construction of so-called Padé triangle which is presented in Table 2. Our estimate given by the most reliable approximants is  (8). So strong scattering of Padé estimates is not surprising since the behavior of expansion coefficients is rather irregular. To improve the situation we can use some tricks which are usually employed for acceleration of iteration convergence. As in [38], e. g., we can find expansion for inverse quantity: Constructing Padé approximants for this series we arrive to estimate 6.6(1.8) which also suffers from huge error bar. Another way to improve estimate is to use the insight about the value of searched quantity in different spatial dimensiond = 2 where according to [37] n + (2, 2) = 2. Keeping in mind this point one can reexpand the initial series for n + (2, 4 − ε) Analysis of a + (2, 4 − ε) expansion within Padé approximant approach gives n + (2, 3) = 6.15 (19) while analogous treatment of the series for inverse a + (2, 4 − ε) results in 5.73 (44). The numerical estimates may be considerably improved if we address more powerful Padé-Borel-Leroy (PBL) resummation technique. This approach described in a good number of papers allows one to optimize the resummation procedure by proper choice of the shift parameter b (see, for example, eq. 32 in [79]) and relevant comments). Having performed such an optimization we obtained b opt = 1.02. PBL triangle for n + (2, 3) under this value of b is presented in Table 3 and yields the estimate 5.8 (8). The previous result extracted from the five-loop series is 5.47(7) but its accuracy was claimed by the authors as underestimated [38]. We applied the same approach to the biased expansions (25) and (26) as well. This leads to n + (2, 3) = 6.1(4) and n + (2, 3) = 6.1(1) respectively. Another resummation procedure known to be rather effective is also based on Borel transformation but analytical continuation is constructed by means of some conformal mapping of Borel image. This transformation helps to isolate singularities which exist on ε complex plane. It is important to note that Borel summability (the presence of at most countable number of singularities on negative real axis and factorial or weaker growth of coefficients) is not proved in general for ϕ 4 field theories. For realization of the conform-Borel (CB) procedure it is necessary to know the asymptotic behavior of coefficients of perturbative expansions for Green functions. Such asymptotic can be obtained by means of the steepest descent method applied to action (1) [89,90]. Such an analysis for 3D O(m) × O(n)-symmetric model was done in [40]. We performed analogous analysis for d = 4 − ε and determined corresponding singularity in (g1, g2) plane, with this information it is possible to find coordinate of the singularity of the series in ε (see e.g. [91]). Coordinate of the corresponding singularity presented in Appendix A.
Apart from singularity closest to the origin which defines the convergence radius for Borel image, there are other parameters which can be chosen when following the strategy suggested in [73]. In so doing, we found the estimate n + (2, 3) = 6.0(6) given by CB resummation procedure. To give a general impression of all the numerical results obtained for n + (2, 3) we collect in Figure  2 corresponding estimates as functions of the order of perturbation theory. The final estimate resulting from these data is n + (2, 3) = 5.96 (19).

m = 3: systems with noncoplanar ordering
From the physical point of view the interest of this study is driven by the fact that under m = 3 the model (1) describes the critical behavior in frustrated pyrochlore antiferromagnets [92,93,94,95,96,97,98]. In this case we perform all the same steps as for m = 2. Series coefficients for n +,− (3, 4 − ε) up to five-loop order coincide with those found earlier [37,44,38]. Numerical estimates obtained for n H (3, 3) and n − (3, 3) by means of different resummation procedures turned out to be stable not only with respect to resummation method applied but also when moving from one perturbative order to another. All numbers obtained are presented in Table 10. Our final estimates for n − (3, 3) and n H (3, 3) are 1.408(4) and 0.973 (11) respectively. Note that the estimate for n − is in a very good agreement with its five-loop counterpart 1.409(1) [38].
As for n + (3, 3), the Padé triangle with numerical estimates for initial series is presented in Table 4. Despite the fact that the higher-order estimates are appreciably scattered, they indicate that n + (3, 3) is larger than 9. For improving the convergence of iterations we address as earlier the PBL approach; corresponding triangle is shown in Table 5. Following the chosen strategy of finding the shift parameter b we obtain 9.3(5) for n + (3, 3). Apart from the applying different resummation procedures we use the same tricks as previously to get a series with more friendly structure. First we consider the series for inverse quantity The analysis of corresponding Padé approximants gives 9.2(2.3) while evaluation by means of PBL approach leads to ∼ 10.4. The next step we address employs the idea suggested in [38] to reexpand initial series basing on the assumption that n + (3, 2) = 2 [37]. Thus we find and series for inverse biased part Padé estimates for (28) and (29) are 9.3(9) and 9.3(2) respectively while PBL results are 9.1(4) and 9.6(1.0). The last resummation technique we apply is CB procedure. By means of this approach we found 9.3 (4). As in the case m = 2 we illustrate all the numerical estimates and their accuracy in Figure 3. The final estimate resulting from the data obtained is n + (3, 3) = 9.32 (19). Previous field-theoretical estimates extracted from five-loop ε expansion and found within six-loop 3D RG analysis are 9.5(5) [38] and 11.1(6) [40] respectively. In this subsection we present, without details, numerical estimates of n + (m, 3), n − (m, 3), and n H (m, 3) for all mentioned values of m including already considered. To find these estimates we implemented some of the steps which were applied previously for m = 2, 3. The only things which will be demonstrated here apart from the numbers themselves are Padé and PBL triangles for the series from Table 1 which has the most favorable structure. As an example of a such quantity we choose n H (5, 4 − ε). Corresponding triangles with Padé and PBL estimates are presented in Table 6 and Table 7 respectively. It is remarkable that in case of PBL procedure the numerical estimates given by two best (highest-order and near-diagonal) approximants differ from each another only in the fourth decimals. This tells us once more about the crucial role of a series structure.
Referring to the main goal of this section we collect the whole information concerning numerical estimates for n + , n − and n H for relevant values of m in Table 8. To get the final values 13 we considered all estimates obtained for given quantity within the highest order of perturbation theory as the results of independent measurements. In addition, to clearly illustrate the situation we present in Figure 4 the diagram of stability of different fixed points including the chiral one in axes (m, n). The empty (white) area corresponds to continuous phase transitions into chiral state while the lightest gray marks the region of computational uncertainty of its lower border.

4.2.
Critical exponents 4.2.1. η, γ, and ν: m = 2 In this section we calculate the critical exponents of chiral universality class for m = 2 and some relevant n. As is well known only two of the exponents are independent, others can be found by means of scaling relations. Here we consider only the most popular exponents -η, γ, and ν. In previous section we noted that to get 3D (ε = 1) numerical estimates of critical exponents in case of n 21.798 we face the problem of complexity of fixed-point coordinates due to (18) and consequently complexity of the exponents. To overcome this problem some trick was suggested [37]; its main idea is as follows. For n in range (n + (2, 3), 12 + 4 √ 6) we reexpand the series for fixed-point coordinates and critical exponents after substituting n = n + (2, 4−ε)+∆n, where ∆n is fixed. Values of ∆n are determined by final six-loop estimate of n + (2, 3) and the value of n we want to consider. As was found in previous section n + (2, 3) = 5.96 (19). Thus the interesting values of n for our consideration are n = 6, 7, . . . .
Order   All the coefficients for ν −1 and η up to five-loop contributions are in complete agreement with the results of earlier calculations [38]. Making use of Padé, PBL and CB resummation procedures we obtain numerical estimates for critical exponents. The numbers thus obtained are presented in Table 9. As was mentioned in [38] addressing the trick described (substituting of shifted n + value Table 9: The numerical estimates of critical exponents for m = 2 and n = {n + (2, 3), 6,7,8,16, 32} obtained by means of different resummation strategies applied to six-loop ε expansions. "-" indicates that convergence of corresponding resummation procedure failed. It is remarkable that for η our algorithm did not find any reliable approximant in case of PBL resummation procedure. "?" tells about inconsistency of numerical estimates taking into account error bars. RPresummation procedure.    (4) while PBL one was found to be 0.63 (4). All the above remarks were kept in mind when calculating the error bars for other critical exponents.

4.2.2.
Correction-to-scaling exponents ω 1 and ω 2 : m = 2 Apart from the critical exponents, the correction-to-scaling exponents ω 1 and ω 2 are of prime importance from the physical point of view. As was already said these quantities determine the stability of the fixed point. Basing upon the definition (20) and the trick used above we can calculate corresponding ε expansions for any value of n n + (2, 3). Taking into account the remarks made in previous section we found numerical estimates for ω 1 and ω 2 . They are presented, along with those for critical exponents, in Table 9. It is necessary to give some comment concerning ω 2 (n = 6). Following the adopted strategy we apply all the resummation procedures to expansions calculated at n = 6 itself and at 6.19 in order to take into account the inaccuracy of determination of n + (2, 3). Since the values of ω 2 in corresponding points are very close to zero, even small variation of n can lead to noticeable scattering of numerical estimates. One can see this for ω 2 (n = 6), where we present the range within which the genuine value of ω 2 (n = 6) has to be.

Discussion
Having calculated the numbers of interest we can compare them with the results obtained by means of different field-theoretical approaches and also with estimates found in lower-order perturbative orders. As was already said, among all marginal dimensionalities n + is most interesting from the physical point of view. We collect all the estimates of n + (m, 3) for m = {2, 3, 4, 5} obtained by means of various approaches in Table 10. In fact, this Table is a continuation of Table  1 in [38] supplemented by the six-loop expansion results. Note that our estimate for n + (2, 3) is less than six although corresponding error bar does not exclude the conclusion that systems with O(3) × O(6) symmetry should undergo first-order phase transitions. Analogously, we aggregate the known estimates for the critical exponents in Table 11. The table does not contain the results for n = 6 because in previous works this case was not considered as actual because corresponding estimates of n + (2, 3) were found to exceed 6. However for this value of n, besides estimates obtained in current paper, there are numerical results obtained within ERG and MC approaches. For critical exponent ν the methods of ERG and MC gave 0.707 and 0.700(11) respectively while our estimate is 0.648 (19). For the exponent γ these techniques yielded 1.377 and 1.38(4) whereas our estimate is 1.27 (3).
For n n + (2, 3) our results are close to their counterparts obtained within the lower order in ε. In some cases a noticeable discrepancy with the estimates obtained by means of 3D RG and 1/n expansion is preserved. It is seen that the difference between the 1/n expansion estimates and others decreases with growing n what looks quite natural since the area of applicability of 1/n expansion is large values of n. Concerning the differences between 3D RG and ε expansion estimates it is worthy to recollect that in physical dimensions (3D) both expansion parametersrenormalized coupling and ε -are not small preventing relevant iterations from perfect convergence.

Conclusion
To sum up, we have calculated six-loop RG expansions for O(n) × O(m)-symmetric (chiral) model in 4 − ε dimensions. The series of record length together with various resummation procedures have allowed to obtain advanced numerical results for physically interesting quantities including marginal dimensionalities n + , n − , n H that separate different regimes of critical behavior and determine the order of phase transitions in concrete systems. The estimate for n + (2, 3) found above is 5.96 (19) indicating that transitions into chiral phase in real magnets with noncollinear ordering (n = 2, 3) should be first-order while the chiral class of universality is appropriate to models with n = 6 and bigger. This conclusion is in agreement with the results of MC simulations [46] and ERG analysis [57]. For systems with noncoplanar ordering we have obtained n + (3, 3) = 9.32 (19) enabling one to conclude that chiral phase transitions in such materials should be first-order too. The inaccuracies of the estimates for n + (m, 3) have turned out to be not small even in the highest-order approximation available what is caused mainly by rather unfavorable structure of corresponding series. Nevertheless, these estimates definitely exclude the possibility of continuous transitions into chiral phases in helical magnets and frustrated antiferromagnets. The six-loop ε expansions for chiral critical exponents under n ≥ 6 have been also derived and corresponding numerical estimates have been found.

Acknowledgment
This work has been supported by Foundation for the Advancement of Theoretical Physics "BASIS" (grant 18-1-2-43-1). A.K. is especially grateful to Professor Pasquale Calabrese for his hospitality during the stay at SISSA where the part of the work was done.

Appendix A. Asymptotics of ε expansions for chiral model
In order to resum asymptotic series, one needs to know large order behavior of their coefficients. Having obtained the series for g 1 and g 2 in terms of ε and following suggested in [91] idea about determining the closest to the origin singularity ε b of ε expansions we found for the chiral model: where g * 1,1 and g * 1,1 -first-order coefficients of ε expansions for the fixed-point coordinates. Taking into account suggested in [37] trick (substitution of n = n + (4 − ε, m) + ∆ instead of n), these coefficients are as follows They are also presented in Mathematica file (asymptotics chiral.m) as "g1" and "g2".