Multicritical hypercubic models

We study renormalization group multicritical fixed points in the $\epsilon$-expansion of scalar field theories characterized by the symmetry of the (hyper)cubic point group $H_N$. After reviewing the algebra of $H_N$-invariant polynomials and arguing that there can be an entire family of multicritical (hyper)cubic solutions with $\phi^{2n}$ interactions in $d=\frac{2n}{n-1}-\epsilon$ dimensions, we use the general multicomponent beta functionals formalism to study the special cases $d = 3-\epsilon$ and $d =\frac{8}{3}-\epsilon$, deriving explicitly the beta functions describing the flow of three- and four-critical (hyper)cubic models. We perform a study of their fixed points, critical exponents and quadratic deformations for various values of $N$, including the limit $N=0$, that was reported in another paper in relation to the randomly diluted single-spin models, and an analysis of the large $N$ limit, which turns out to be particularly interesting since it depends on the specific multicriticality. We see that, in general, the continuation in $N$ of the random solutions is different from the continuation coming from large-$N$, and only the latter interpolates with the physically interesting cases of low-$N$ such as $N=3$. Finally, we also include an analysis of a theory with quintic interactions in $d =\frac{10}{3}-\epsilon$ and, for completeness, the NNLO computations in $d=4-\epsilon$.


Introduction
The original motivations for considering the cubic point group and its generalization, the hypercubic point group H N , in statistical mechanics was to describe the critical properties of classical three-dimensional spin systems in which the interactions are modulated by an underlying cubic lattice [1,2]. In some magnets, like for example Fe or Ni, ions are arranged in crystals with cubic symmetry: non-rotationally invariant single-ion interactions such as i,α (σ α i ) 4 (σ α i is the i-th component of the spin vector σ at the lattice point α) are allowed alongside the "standard" isotropic term i (|σ i | 2 ) 2 . The inclusion of anisotropic terms in the N -component isotropic Ginzburg-Landau (GL) Hamiltonian changes the critical behavior from the O(N )-invariant isotropic one to a distinctive hypercubic critical behavior provided that N > N c 2.9 in d = 3 dimensions. Using the renormalization group (RG) approach to critical phenomena, this can be seen from the fact that the hypercubic critical point is more infrared relevant than the isotropic one for values of N above N c , which therefore includes the physically relevant case of N = 3 [3,4]. Given that N = 3 is just barely above N c , the case of cubic anisotropy is of paradigmatic importance in exemplifying how the symmetry of a system influences the nature of its phase transition.
Another important reason has to do with the inclusion of random disorder in traditional magnetic systems with Z 2 symmetry. For sufficiently low concentrations of non-magnetic impurities, such as Zn atoms in a Fe based Ising-like lattice, a ferromagnetic system still exhibits a second-order phase transition, but its critical properties are different from those of the pure system [5]. To include the effect of disorder in the GL description one usually resorts to the replica approach, in which N identical copies of the system are introduced. Upon integration over the distribution of the disorder, the system has an enhanced H N symmetry and it is possible to show that the critical exponent of the quenched averages of correlators can be obtained from the limit N → 0 [6]. Using again the RG approach, it is possible to see that the replicated theory admits both a decoupled FP, representing N copies of the original pure system, and a new random FP with true H N symmetry, representing the diluted system. Similarly to the previous motivation, the random FP is more infrared relevant than the decoupled one [7].
The above discussions are implicitly based on a GL Hamiltonian with φ 4 -like interactions, but, in principle, higher-order interaction terms of the form φ 2n can also be included in the GL description, both with isotropic and H N symmetry [8,9]. These are generally discarded, at least in qualitative approaches, by following standard considerations based on the canonical dimension of the operators that they would correspond to in the GL Hamiltonian for a three-dimensional system. Higher-order operators, however, play an important role in lower dimensional systems, because they give rise to multicritical behaviors, given that they become increasingly RG relevant when lowering the dimensionality of the system. Rather generally, by means of a standard Hubbard-Stratonovich transformation that introduces a scalar N -component field variable φ i , microscopic terms like i,α (σ α i σ α i ) n , representing a higher-order isotropic interaction, and i,α (σ α i ) 2n , representing a higher-order single-ion interaction, for n ≥ 3 can be expressed as i (φ i φ i ) n and i φ 2n i , respectively. A GL Hamiltonian based on these φ 2n -like interactions can be renormalized perturbatively in d c = 2n n−1 dimensions, and its renormalization gives the critical properties of the multicritical phase transitions in terms of an -expansion in d = d c − dimensions [10]. Notice also how the textbook example of φ 4 is recovered for n = 2.
Multicritical models are important for both theoretical and experimental reasons. Two well-known examples of multicritical behaviors based on the isotropic model with φ 6interaction are provided by the tricritical point of He 3 -He 4 mixtures, which is microscopically described by the so-called Blume-Capel model [11][12][13], and the Flory θ-point of very long polymer chains, which is related to the limit of N → 0 components [14][15][16]. On the more theoretical side, single component's φ 2n models with Z 2 symmetry can be seen as multicritical generalizations of the Ising universality class and are famously connected to the infinite set of unitary minimal models M n+1,n+2 , that are conformal field theories (CFTs) emerging from the representations of the infinite-dimensional conformal group in d = 2 [17]. It is possible, though still unknown for the most part, that more complicated symmetry groups, such as the aforementioned O(N ) and H N , admit the same type of multicritical generalization of the Ising case and are also CFTs with a special meaning in d = 2 [9,18,19].
In this paper, we move towards a better understanding of the last paragraph by analyzing the renormalization group fixed points of multicritical scalar field theories obeying the symmetry of the hypercube. The paper is structured as follows: In Sect. 2, we outline the algebra of hypercubic invariant polynomials, which allows us to determine the basis of invariant polynomials needed in the construction of the GL Hamiltonian. In Sect. 3, we discuss the representation theory of the hypercubic group in relation to quadratic composite operators, because they have a pivotal role in discussing the symmetry breaking patterns connecting various critical points. In Sect. 4, we review the main points behind the perturbative RG method that we use to obtain the -expansion. Sections 5 and 6 contain our results on the multicritical generalizations of the hypercubic model with φ 6 -and φ 8 -like interactions, which are the main results of this paper. Conclusion and perspectives are given in Section 7. We also include appendix A, which summarizes several results for the hypercubic theory in d = 4− dimensions for completeness, and appendix B, which presents a non-unitary generalization of the hypercubic model in d = 10 3 − dimensions.

Algebra of hypercubic invariant polynomials
We begin by reviewing some fundamental facts on the algebra of invariant polynomials of the hypercubic point group H N , which we need to construct a basis for the interaction potentials. For this purpose, let R be the algebra of all polynomials f in the variables {φ 1 , . . . , φ N } with coefficients in C. The hypercubic group acts on {φ 1 , . . . , φ N } as the subgroup of O(N ) that leaves a unit hypercube centered at the origin of R N invariant. The group can be generated from all permutations of the axes, and a parity for each axis, resulting in the isomorphism H N (Z 2 ) N S N , where S N is the permutation group of N objects. We are interested in the subalgebra R H N of all the polynomials f ∈ R satisfying gf = f for each g ∈ H N . It is important to notice two things. First, the algebra R H N admits the decomposition consists of all the invariant homogeneous polynomials of degree k, in the usual sense. Hence, in order to determine R H N , it is sufficient to determine each R H N k . Second, we can introduce a generating function for the dimension of the homogeneous component R H N k . This is a formal power series in a variable t, known as the Hilbert series (also known as the Poincaré series, or the Molien series) of the ring R H N The Hilbert series can be determined explicitly for the ring of any finite group through Molien's theorem. In the case of H N , one finds This result is very useful for us since it determines the size of the ring R H N and, therefore, the number of couplings that we need to parametrize the potential. Once the series is expanded, the coefficient of the monomial t 2n gives the dimension of the corresponding subspace R H N 2n . We summarize the first few values of dim R H N k in a table.
The counting rule above applies, to some extent, rather generally, so it can be used to continue our results analytically in N , and eventually it is responsible for the determination of N -dependent renormalization group functions. It is important to stress, however, that the counting rule of the number of operators applies only for N sufficiently big in comparison to the order of the interaction 2n, that is, N ≥ n. If instead we concentrate on the natural values for which N < n, then the number of independent monomials (and therefore of associated coupling constants) is actually lower, because one or mode interactions of the general basis can be expressed in terms of the others. In these cases, which require separate analyses given in the later parts of Sections 5 and 6, the system of beta functions is reduced in a specific way that accounts for the degeneracy of linear combinations of the operators.
We find most convenient to work always with the general case, and eventually specialize the results later. In this way, one can show that the number of invariant polynomials of degree 2n corresponds to the number p(n) of partitions of the positive integer n. To clarify this relation, we first define generalizations of the Kronecker delta, which are symmetric tensors with 2n indices defined as Homogeneous invariant polynomials of degree 2n in the variables φ i can be obtained by considering all possible ways to contract the product φ i 1 . . . φ i 2n with combinations of the basic tensors δ i 1 ...i 2k with k ≤ n. This can be done exactly in p(n) ways, as expected from the Molien series. A basis made with all possible contractions is manifestly invariant under H N , in fact the hypercubic group acts by reflecting and permuting the variables φ i , namely through the transformations φ i → −φ i and φ i → φ j for all i and j, thanks to the fact that δ i 1 ...i 2n is a symmetric tensor with an even number of indices.

Quadratic deformations
In this section, we consider the study of non-singlet quadratic operators since they turn out to be particularly important in the determination of the so-called crossover exponents [20]. On general grounds, there are physical situations in which the critical behavior depends on the interplay of more than one order parameter, or, alternatively, on different components of a multicomponent order parameter. In these cases, the lattice Hamiltonian typically presents a quadratic term, weighted by some parameter g, which reflects weak anisotropies between such competing order parameters. By approaching the critical temperature m c (g), the isotropic transition crosses over to a transition characterised by a smaller symmetry group. This is the case, for example, of the cubic-to-tetragonal structural phase transitions of anisotropically stressed perovskites like SrTiO 3 , in which the S 4 × Z 2 cubic symmetry for m > m c is reduced to the D 4 × Z 2 tetragonal one for m < m c [21]. The change to the reduced-symmetry form occurs in the vicinity of a crossover "temperature"m = m c (g)+∆m and can be characterized by a crossover exponent Φ according to ∆m ∼ g 1/Φ [22]. In the case of the cubic symmetry, we can introduce the following quadratic spin operator where again σ α i is the i-th component of a spin vector σ at the lattice site α. We then expect two crossover exponents: one Φ axis which characterizes the crossover to a phase where an axial interaction (i 1 = i 2 ) is favoured and another one Φ diag which instead favours the alignment along the diagonal (i 1 = i 2 ) of the cube. In both cases we expect a crossover from an N -component to an Ising-like behavior, however, if the N -component transition is isotropic, then the two crossovers are equivalent since spins can be rotated so that E i 1 i 2 becomes E i 1 i 1 ; if instead the transition is characterized by a truly cubic fixed point (this is the case of the standard cubic transition in d = 3 for which we recall N c 2.9), then the two crossover exponent becomes distinguished. In fact, the separation of these two exponents can be used as an alternative criterion to determine N c , the other being the value of N at which the O(N ) and hypercubic critical points coincide [3].
On the field-theoretical side, the discussion of quadratic operators is naturally encoded in the representation theory of the hypercubic group (a detailed analysis oriented to the construction of the spectrum of composite operators can be found in Ref. [23]). Here we sketch the essentials needed for this paper. For this purpose, it is useful to consider the following expression for the fundamental N -dimensional representation of H N (i.e. the lower dimension non-trivial representation) where the double partition (α, β) represents α objects, even under Z 2 , which transform under an irreducible representation of S α , while the right partition represents β objects, odd under Z 2 , which transform under the trivial representation of S β . One can then construct the tensor product φ i ⊗ φ k , which can be decomposed as follows and we start filling the tableau with indices that are all different. Each Young tableau then corresponds to an irrep of H N . Since the left partition represents objects which are even under Z 2 , we associate an indexed field raised to the zeroth power to the first row and a second field power to the second row. We then symmetrize w.r.t. the rows and antisymmetrize w.r.t. the columns. The same holds for the β partition with increasing odd powers of the field. In terms of these simple rules, we are able to explicitly write the irreps considered as and the last irrep in Eq. (3.3) is anti-symmetric, thus it cannot appear in quadratic oeprators with no derivatives and, consequently, it will be discarded in the following. We understand therefore that the tensor product in Eq. (3.3) is decomposed into the direct sum of the trivial representation S, a traceless diagonal symmetric irrep X, an off-diagonal symmetric Y and finally an anti-symmetric one A. 1 Since these irreducible representations are the representations carried by quadratic scaling operators, we can eventually assign an anomalous dimension γ R for R = S, X, Y to each one of them. For this purpose, one can introduce the following four projectors in the space of quadratic operators P (S) in terms of which we can decompose the general quadratic operator with no derivatives as It could be checked that these projectors satisfy the following relations where dim(R) stands for the dimension of the representation R. In the case at hand we have that dim(S) = 1, dim(X) = N −1 and dim(Y ) = N (N −1)/2. We conclude by stating that the quadratic deformation X can be related to the diagonal crossover exponent Φ diag , while the deformation Y to the axial one Φ axis .

Effective potentials and their renormalization
The decomposition in subalgebras carried out in Section 2 can be rephrased as follows: each of the R H N 2n can be labelled with the critical dimension at which the corresponding polynomials are canonically dimensionless, namely d c = 2n n−1 . In the functional perturbative RG, on the other hand, the upper critical dimension d c (and therefore the label n) is taken as an input for the classification of multicritical universality classes [24,25]. The typical starting point to analyse the critical behavior of a given theory is the following GL action With the only input of d c , one can derive very general expressions for the renormalization group flows β V of the effective potential V (φ) and β Z ij of a general wavefunction, which is related to the anomalous dimension matrix γ ij = − 1 2 Z −1 ik β Z kj (the latter will be given in the form of the singlet γ ij φ i φ j ). In the following, we display these beta functionals using a diagrammatic representation in which k-vertices correspond to k derivatives of the potential w.r.t. the fields, namely in which the gray circle represents the j-vertex in question, hollow circles are the other vertices, and lines organize the contraction of indices. For example, the following diagram represents It is worth sidestepping for a moment and consider one of the most remarkable and striking consequences of the functional formalism, namely the fact that some of the wellknown expressions derived for the beta functionals in the single-component field case [24] can be directly enhanced to their multicomponent counterparts without performing any single additional computation [8]. Few remarks are in order. The basic idea is that the multicomponent (N > 1) beta functionals β V and β Z must reduce to their known single-component cases in the N = 1 limit. In the single-component case, typical monomials contributing to β V and β Z are constructed in terms of derivatives V (k) of the potential, possibly linked by the appropriate powers of V (2) insertions to have proper vacuum diagrams. In promoting these monomials to the multicomponent case, namely rewriting V (2) → V i 1 i 2 and V (k) → V i 1 ...i k , we notice that in dimensional regularization, self-contractions like V i 1 i 1 ...i k or traced mass insertions like V i 1 i 1 (which are generated by the expansion of the propagators in the loop diagrams) do not contribute. Therefore the relative coefficients appearing in front of the monomials in the single-component case are univocally fixed and stay the same in the N > 1 case.
The single-component functional RG flows for the potential and the wavefunction can be found in [24][25][26] and can be used to obtain their generalizations. Let us emphasize that the generalizations to the multicomponent cases of the RG equations for d c = 6, 4, 3 have already appeared in [27], but here we want to make the general point that we can generalize the RG equations for any d c,m to some extent.
For even potentials (apart the trivial LO d c = 4 case), in the NLO beta functional β NLO V , there appears only one term involving V i 1 i 2 and this term is linear in the mass insertion [24]. Because of this fact, there is only one multicomponent diagram generalizing this term, and its relative coefficient can be fixed univocally. A useful example is the d c = 3 case [24]: the N = 1 beta functional is given by for which the NLO monomial V (2) (V (5) ) 2 is the only involving a mass insertion and the V i 1 i 2 can only connect the two five vertices. The vertices in other NLO monomials can be connected only following the single component loop diagram from which they come from. Therefore, the generalization to (4.4) reads (4.5) In the odd case it is possible to derive a general formula for the N = 1 beta functional β V only at LO (see Ref. [25]). Since this formula can be obtained by analytical continuation of the NLO even case, it involves only one V i 1 i 2 in only one monomial and, as before, its multicomponent counterpart can be fixed univocally. The picture is different for the odd NLO case for which, apart the d c = 6 case, no general functionals β V and β Z are known. In the N = 1 single component d c = 6 case, the beta functional is 3 and since V (2) are mass insertions and since the NLO is at two loop (vacuum sunset diagram), there are three propagators where the three mass insertions can fit and thus the knowledge of the N = 1 case only fixes the sum α + β + γ [27]. Finally, in the case of the LO (even and odd) β Z there is only one monomial with no V i 1 i 2 mass insertions and the results of the even case agrees with CFT results [28]. The determination of the NLO contributions in the odd case will complete the knowledge of the universal data of one component scalar field theories. The discussion above provides the rationale behind the analysis of the critical content of scalar field theories endowed with hypercubic symmetry using the functional perturbative RG. In this case, taking advantage of the results of Section 2, the GL effective potential of Eq. (4.1) can be expressed as namely, we assign a coupling constant λ j,n to each of the p(n) invariant polynomials I (2n) j composing the corresponding R H N 2n subspace. For a given d c , dimensionful β-functions for the couplings λ j,n follow by inserting the explicit form of the potential in the corresponding β V and by simplifying the algebra of Kronecker symbols resulting from contracting the fields. In terms of the fixed points of the associated system of dimensionless beta functions, one can finally express universal quantities, such as critical exponents, as power series in = d c − d.

Universal quantities
In dealing with the quadratic operators of Section 3, we can use Eq. (3.8) to source individually each irrep contribution in the path integral and to renormalize them as composite operators. We introduce three sources J (R) , for R = S, X, Y , in the combination Naturally, the indices of the source J (S) of the scalar S operator are redundant as it couples to its operator only through the trace, but contractions are accounted for by the projectors. The computation of the expectation value O qu (φ) requires renormalization of the three sources, corresponding to the above three operators, and consequently gives three independent gamma functions. For almost all values of N , the three operators are renormalized multiplicatively and do not mix, that is, they are already scaling operators (scaling operators thus carry a scaling dimension and a label of the irreps of H N ). From our point of view, the most convenient route for the computation is to determine the renormalization group running of the three composite operators using β V . In fact, the advantage of using a functional formalism is that the running of the potential V (φ) already includes all the information on the running of all relevant operators, including (4.8). To extract this information it is sufficient to replace V (φ) → V (φ) + O qu (φ) in β V and retain only the leading order in the three source couplings J (R) . By construction, this operation provides the gamma functions γ R of the three operators multiplied by the corresponding sources, which can easily be related to either the critical exponents θ R or the scaling dimensions ∆ R . 2 In particular, the critical exponents θ R are obtained as the negatives of eigenvalues of the stability matrix M ij = ∂ j β i . In bridging to CFT, we recall that the eigenvalues θ R are related to the CFT scaling dimension ∆ R of scaling operators by the relation θ R = d − ∆ R . Alternatively, one can retrieve the gamma functions γ R as γ R = θ R − (d + 2 − η)/2, where the anomalous dimension critical exponent η is obtained by diagonalizing the matrix γ ij with an additional factor 2. In particular, the correlation length critical exponent ν is related to the inverse of the quadratic singlet in V (φ) as ν = θ −1 S , while the crossover exponents discussed in Section 3 are found to be Φ axis = θ X /θ S and Finally, we compute the A-function, which is the scalar function from which one can derive the RG flow as a gradient, for all the even models. This can be of particular interest as it entails very important implications on the asymptotic behavior of the renormalization group equations, as well as on the classification of renormalization-group fixed points [27,30,31]. If we call λ a all the couplings of interest and β a the corresponding β-functions, 2 Notice that for almost all values of N the irreducible operators appearing in (4.8) are automatically scaling operators, which means that there is no mixing induced by radiative corrections to our decomposition. In other words, the generator of the dilatations D commutes with the action of the hypercubic group on the irreducible decomposition. This is not true, however, for all values of N . In fact, by changing appropriately N we can end up in the situation in which two or more scaling operators degenerate to the same multiplet because their scaling dimensions coincide. This results in a nontrivial Jordan cell and induces a logarithmic contribution. Interestingly, this happens for the special limits N = 0 (see Ref. [29]) and N = −2 of the hypercubic model, which is obviously reminiscent of the same limit of the O(N ) model that has recently been associated with the loop erased random walk [19].
then A is derived implicitly from and h ab is a suitable positive metric in the space of all couplings. We require the metric to be flat in the space of marginal couplings parametrized as the symmetric tensors V i 1 ,··· ,i 2n , which can be always be done at leading and next-to-leading order [30,31]. Consequently, we express the A-function from contractions of V i 1 ,··· ,i 2n when possible, and it is understood that beta functions are derived from the gradients and there is a flat metric. 3 5 Tricritical models in d = 3 − The discussion of the critical properties in d = 3 − dimensions would require the simultaneous fine-tuning at criticality of other two parameters beside the mass, namely the two quartic non-linearities λ 1 and λ 2 . This being said, in the following we shall still refer to RG fixed points possibly emerging at d c = 3 as tricritical ones, despite their multicritical nature. Critical properties can be accessed in terms of the following potential In this case, LO and NLO beta functionals β V and β Z and the expression for the A-function are given by The β-functions of the couplings, the anomalous dimension matrix γ ij and the A function are then readily obtained by inserting the effective potential of Eq. (5.1) in the βfunctionals above and by simplying the algebra of Kronecker symbols that are responsible for the field contractions in V (φ). The β-function of the quadratic (mass) coupling is The β-functions for quartic couplings are given by and The RG of the marginal (sextic) couplings read and By taking derivatives of A with respect the marginal couplings σ a it is straightforward to determine the metric h ab in the space of the couplings and check that it is positive. We report also the anomalous dimensions associated to the quadratic deformations It is straightforward to check that the gamma function of the singlet S is related to the beta function of m 2 as γ S = 2 − ∂ ∂m 2 β m 2 , but, in general, β m 2 contains more information.
Tri-Ising Tri-O(N ) Table 1. Anomalous dimension η, A-function and critical exponents θ R at the tricritical Ising and the tricritical O(N ) fixed points.
The system of β-functions admits a total of 8 fixed points (both real and complex). However, as a first physical requirement, we restrict ourselves to consider critical theories which are perturbatively unitary: we therefore select only real fixed points with positive anomalous dimension η > 0 and characterized by a real spectrum. A second physical In particular, we highlight the annihilation mechanism (inset) in terms of which we distinguish the solution R N (relevant for the discussion of the tricritical random universality classes in the N → 0 limit) from the standard tricritical hypercubic one H N . In the bottom panel, instead, we present the tricritical hypercubic fixed points relevant to the discussion of the large N limit behavior of the theory. Notice the pair H ± N appearing at N = N (+) ≈ 14000 that prevents the H N solution to approach the decoulped limit at N → ∞. requirement is that the fixed point potential (5.1) is stable, namely, bounded from below: this amounts at considering fixed points for which (σ 1 +σ 2 +σ 3 ) > 0. Under these conditions, beside the trivial Gaussian fixed point, the system of β-functions predicts the tricritical decoupled Ising universality class with coordinates {σ 1 = 0, σ 2 = 0, σ 3 = 3 10 } and the tricritical O(N ) universality class with coordinates {σ 1 = 15 2(3N +22) , σ 2 = 0, σ 3 = 0} . In Table 1 we list the anomalous dimension η, the A function and the critical exponents θ R values at the aforementioned fixed points. Some of these quantities can be found in Refs. [27,32,33].
We also find truly hypercubic fixed points which, due to their complexity, cannot be given in a closed analytical form. However, one can elegantly express the hypercubic fixed point solutions as the roots of the following quintic polynomial there are three possible fixed point solutions relevant to the discussion of the "random" tricritical universality class in the N → 0 limit. In agreement with the logic proposed in Ref. [29] (for which one is led to consider only solutions such that σ 3 > 0 in the limit), we identify R N as the only compatible 'random' tricritical fixed point solution, at which we find η = 0.00137842 2 . While the other two solutions, namely R + N and R − N , distinguished for values N < N (−) , annihilate each other at N (−) , the "random" solution R N ceases to exist at N R 2.015. The mechanism behind this fact is reported in the inset of Fig. 1, left panel. At N H 1.975 a pair of fixed points emerge: one branch (dashed line in the inset) is responsible for annhilating the 'random' solution R N at N R , while the other branch H 0 N , which exists in all the range N > N H is the standard (tricritical) hypercubic one. Concerning the large N limit of the φ 6 -hypercubic theory, we notice that at N (+) 14000, a pair of fixed points emerge, namely H + N and H − N and they become distinguished in the large N limit [27]: while H − N eventually annihilates the H N one, the fixed point H + N tends to the tricritical Ising universality class which is characterized by σ 3 = 3 10 , see the discussion above. This interesting mechanism shows that the large N limit of the hypercubic solution H N does not correspond to the actual large N limit behavior of the hypercubic theory at d c = 3, which is instead carried by the solution H + N . Finally, for small values of N , specifically N = 1, 2, there are less independent monomials in the potential than for the general case N ≥ 3, thus separate discussions are required. For N = 2, the number of independent monomials in (5.1) is reduced from three to two, because, for example, one can write the monomial multiplying σ 2 as a linear combination of the other two. One can thus eliminate the σ 2 in favor of the other couplings, being left with the following system of reduced β-functions This system predicts only the tricritical Ising and the tricritical O(2) universality classes (we could have eliminated any other coupling, but the final RG system would have been the same, as we checked explicitly for this and all cases shown later). This somehow reflects the critical case in d c = 4, for which the N = 2 cubic theory predicts only the Ising and O(2) fixed points (see Appendix A). We deduce that the nontrivial hypercubic fixed point solutions reported in Fig. 1 at N = 2 are simply understood as analytical continuations of the two solutions R N and H N , coming from N = 0 and large-N , respectively. Analogously, in the case N = 1 the number of independent couplings is reduced to one, which we can choose to be σ 1 (thus eliminating σ 2 and σ 3 ). In this case the following reduced β-function predicts only the tricritical Ising universality class (see Table 1 for the associated universal quantities). Concerning the crossover exponents, from Table 1 we see that, as expected, at the tricritical isotropic O(N ) fixed point, the critical exponents associated with the deformations X and Y are identical and so they are the corresponding crossover exponents. At the tricritical Ising fixed point, instead, one can notice that θ X = θ S , while the quadratic deformation associated with Y simply amounts to a shift in the critical temperature. We have also analysed the quadratic deformations at the hypercubic fixed points discussed above, which come from inserting the fixed point values in (5.11). The final results are rather complicate and not particularly illuminating, so we do not report them here.
The anomalous dimension is obtained from diagonalizing γ ij with an additional factor two, as usual, and, in this case, it reads (6.14) The system of β-functions admits a total of 32 fixed points, including real and complex ones. However, as in the tricritical case, we restrict ourselves to consider critical theories which are real and perturbatively unitary, thus we limit our attention only to real fixed points with positive anomalous dimension η > 0 and real spectrum (the spectrum is real for all the operators of the basis of the operators that we included to the potential, we cannot exclude that irrelevant operators come with imaginary parts). We then require the fixed point potential (6.1) to be stable, namely, bounded from below: this amounts at considering fixed points for which the sum of the marginal couplings is positive, (τ 1 +τ 2 +τ 3 +τ 4 +τ 5 ) > 0. These requirements boil down the plethora of fixed points to a manageable subset. Apart the trivial Gaussian fixed point, the system of β-functions admits the tetracritical Ising fixed point with coordinates given by {τ 1 = 0, τ 2 = 0, τ 3 = 0, τ 4 = 0, τ 5 = 3 70 } and LO anomalous dimension η = 9 2 85750 [8]. We then find a fixed point representing the tetracritical O(N ) universality class with coordinates {τ 1 = 105 2(3N 2 +150N +1072) , τ 2 = 0, τ 3 = 0, τ 4 = 0, τ 5 = 0} and LO anomalous dimension given by η = 3(N +2)(N +4)(N +6) 2 2(3N 2 +150N +1072) 2 . The critical exponents of the tetracritical O(N ) universality class, for the specific values N = 2, 3, 4, can be checked against Ref. [8], while the limit N → 0, relevant to the discussion of the tetracritical self-avoiding-walks (SAW) universality class has been given in Ref. [29].
Additionally, we still find several non-trivial fixed points with true hypercubic symmetry, that, because of their complexity, we are not able to express neither in a closed analytical form, nor as roots of a polynomial of the Gröbner basis like in the previous section. However, they can still be accessed numerically, which will be our main tool of investigation from this point forward. We find it useful to report the resulting scenario in terms of the fixed-point coupling τ 5 as a function of the number of components N , similarly to what we have done in the previous section. The solutions are fixed points of the full system of five beta functions continued in N , but, as previously discussed, it is necessary to consider the special cases N = 1, 2, 3 separately, because there are less independent monomials in the respective bases. The result is thus an analytic continuation in N and is shown in Fig. 2.
We immediately recognize a φ 8 -hypercubic solution with H N symmetry that exists in the whole the range of N with the desired physical properties for N > 3, that is denoted H N likewise the group in Fig. 2. In the large N limit, this solution tends asymptotically to the decoupled tetracritical Ising universality class, for which we recall τ 5 = 3 70 , so the hypercubic spins decouple in the limit (which is a rather different behavior in comparison to the tricritical model that did not decouple). At the values N 1 = 4260 and N 2 = 4410, two new pairs of fixed points emerge, respectively 1 H ± N and 2 H ± N . In the large N limit, the upper branches 1 H + N and 2 H + N tend to H N as far as our numerics can tell (asymptotically to the tetracritical decoupled Ising universality class), instead the lower branches 1 H − N and 2 H − N tend to zero. The pairs do not affect the large-N limit of the H N solution because they emerge "below" it, using the coupling τ 5 as reference (see the right panel of Fig. 2).
It is also interesting to discuss the fate of the solution H N in the limit N → 0, in relation to the randomly diluted tetracritical model. Similarly to the tricritical case, we notice that the multicritical solution H N does not correspond to the N → 0 fixed point In particular, we highlight the annihilation mechanism (inset) in terms of which we distinguish the solution R N ceases to exist. In the bottom panel, instead, we present the tetracritical hypercubic fixed points relevant to the discussion of the large N limit behavior of the theory.
that has been argued to describe the "random" tetracritical universality class identified in Ref. [29] on the basis of physical requirements that hold in relation to the replica approach to describe the diluted system. We denote the continuation in N of the random fixed point as R N to make the distinction with H N clear. In the left panel of Fig. 2 one can see that R N ceases to exist at N R 2.2 by annihilating with a not-physical fixed point displayed with a dashed line.
Finally we must discuss separately the cases N = 1, 2, 3 separately, because in these three cases the five hypercubic terms in (6.1) are not independent. For N = 1, the system can be reduced to a single invariant, namely V (φ) = τ 1 φ 8 , which trivially describes only the tetracritical Ising universality class (with fixed point τ 1 = 12 35 ). For N = 2, we find only three independent couplings: one can eliminate τ 2 and τ 3 in favor of the other three couplings being left with the following system of reduced β-functions This system, similarly to the critical and tricritical cases for N = 2, predicts only the tetracritical decoupled Ising fixed point and the tetracritical-O(2) universality classes, the latter being characterised by η = 9 2 59858 (in other words, there does not seem to any multicritical theory of spins on a square, according to our analysis). The last special case in N = 3, in which there are four independent couplings. Eliminating τ 3 in favor of the others, we are led to the following system of reduced β-functions 4 . The fixed point of the general φ 8 -hypercubic solution H N , which is plotted in Fig. 2, is therefore not representative of a truly hypercubic behavior for the special cases N = 1, 2 as the analysis of the reduced systems of β-functions reveals, but is purely understood as an analytical continuation. If one chooses to study the continuation of the solution to the special cases N = 1, 2, it is easy to see that, because of the redundancy of the basis of operators, they are fixed points with pseudo-Gaussian features and their spectrum coincides with either the O(N ) or the decoupled tetracritical case. As far as our analysis can tell, we do not find a genuinely hypercubic fixed point for the interesting value N = 2, and the first natural ones appear at N = 3.

Conclusion
This paper could be considered part of a series of works [8,[34][35][36] that share the common theme of discussing and classifying "discrete" universality classes, namely renormalization group critical points of scalar field theories which are characterised by a finite discrete symmetry group. In particular, here we concentrated on the multicritical behavior of scalar field theories endowed with the symmetry of the hypercubic point group H N . From a group-theoretical point of view, hypercubes (as well as hypertetrahedra) belong to a family of regular polytopes that exist in R N , with N indicating the number of independent components φ i of the associated order parameter. It is also possible to analytically continue N , which makes the results particularly interesting in connection to the study of its critical behavior as a function of N , which makes further universal properties accessible. For example, a notable limit is N → 0 of the hypercubic theory, which is known to be representative of the universality class of randomly diluted spin systems, for which some of us recently investigated multicritical effects of arbitrary distributions of the disorder [29]. In this paper we moved forward to a systematic analysis of the multicritical behavior of hypercubic scalar field theories as a function of N , paying particular attention to some non-trivial annihilation mechanisms involved in the N → 0 and N → ∞ limits.
To assess the continuation in N correctly, we first reviewed the algebra of hypercubic invariant polynomials to determine the most general basis for the interaction potentials. This must be done with care, because, for a given multicritical interaction and finite and small N , the counting of operators is not general, because a limited number of parameters φ i reduces the total number of independent interactions. This analysis is not restricted to a unique quadratic polynomial, but instead we considered the quadratic deformations associated with irreducible representations of the hypercubic group. This is relevant to the discussion of the crossover exponents that characterize the critical behavior of systems which depend on the interplay of more than one order parameter. Besides the general quadratic deformations, we also included all possible relevant H N -singlets in the potentials, which give further informations on the universal properties and can be used to determine some of the conformal data of the underlying CFT [25].
On general grounds, the discussion of the multicritical behavior of scalar field theories with φ 2n interactions involves, perturbatively, the renormalization group analyses and an -expansion from rational space-time dimensions. Within the functional perturbative RG, this analysis can be carried out systematically [24,25]. Most of the RG equations that we adopted can be derived from simpler single component equations, so we provide the theoretical rationale behind the fact that the single-component beta functionals can be directly enhanced to their multicomponent counterpart.
We have focused on two upper critical dimensions, namely d c = 3 and d c = 8/3, corresponding respectively to a φ 6 and a φ 8 scalar field theory, though our analysis could be generalized to higher rank interactions at the cost of numerical complexity. In both cases, our analysis revealed a rich spectrum of fixed points which, in turn, determine the non-trivial mechanism occurring behind the large N as well the N → 0 limit of the theories. In particular, in d = 3 − our analysis showed that the "tricritical" φ 6 -hypercubic fixed point solution, denoted H N as the group, that interpolates the limits N = 2 and N = 3, does not coincide with the continuation of the "random" universality class, denoted R N , that interpolates with the N → 0 limit given in Ref. [29]. A marked difference is that the "random" model fixed point solution exists only for N 2.015. On the opposite side, the analysis of the large N limit reveals that the standard hypercubic solution H N annihilates asymptotically with the lower branch H − N of a pair of fixed point H ± N that emerge from a complex conjugate pair at N 14000. The upper branch H + N instead tends asymptotically to the tricritical decoupled Ising universality class. In other words, the large-N limit of the tricritical hypercubic model does not tend to (N copies of) the decoupled tricritical Ising model, which straightforwardly implies that if one were to construct a large-N expansion in fixed dimension by seeding the decoupled model solution, such expansion would fail. Notice that this does not happen with the standard critical hypercubic model, that in the large-N limit tends to the decoupled Ising.
In the case of d = 8/3 − , we identified a "tetracritical" φ 8 -hypercubic solution, also denoted H N , that in all the whole range N ≥ 0. Differently from the tricritical solution, in the large N limit the tetracritical solution tends to the tetracritical decoupled Ising universality class In the N → 0 limit, instead, it is still distinguished from the "random" tetracritical universality, again denoted R N , that for N = 0 was give in Ref. [29]. The tetracritical random solution exists for N 2.2. The large-N limit of the tetracritical model can reach the decoupled solution asymptotically because there are no pairs of complex conjugate solutions interrupting its path. In fact, we see numerically two pairs of fixed point solutions, 1 H ± N and 2 H ± N , emerging from complex conjugate pairs at some large values of N , but they do not intercept the H N solution. They still participate non-trivially to the large N limit: the upper branches 1 H + N and 2 H + N of these solutions tend to the same large N limit of H N , that is, to the tetracritical Ising universality class, while the two lower branches tend to zero, as far as our numerical analysis can tell.
One speculation that we could make, based on the first few multicritical hypercubic models that we considered in this work, is that not all φ 2n -hypercubic model have the corresponding decoupled model as large-N limit. It could be that only for n even, the hypercubic theories decouple in the limit, while for n odd they do not. The mechanism preventing the decoupling of the spins in the limit would be the birth, at large values of N , of pairs of fixed points, that collide asymptotically with the hypercubic solution. This mechanism is obviously confirmed by our findings for n = 2, 3, 4. One general strategy to try to prove it could follow our procedure: use the Gröbner basis to find a polynomial in the coupling that weights the (sum of the) decoupled interaction and test the behavior of the roots of this polynomial. The application of this idea is essentially what gives us Figs. 1 and 2 for the special cases n = 3, 4. Similar nontrivial behaviors in the large-N limits of O(N ) multicritical models has been already observed, and it also depends on the value of n being either even or odd [9,[37][38][39].
We discussed separately the RG equations of the φ 6 -and φ 8 -hypercubic models for the special values N = 1, 2 and N = 1, 2, 3, respectively. The reason is that for these values the number of independent operators is less than the case of arbitrary N (for example, for N = 1 there is only one scalar, so only one independent marginal interaction can be formed, but similar simplifications occur for N = 2, 3). The presence of a degeneracy in the number of independent invariant monomials, as seen from the point of view of the RG equations for arbitrary N , makes it such that the general N fixed points see pseudo-Gaussian solutions for the lower N cases (the spectrum has a mixture of Gaussian features and nontrivial ones that cancel out if the correct bases are taken into account). To truly explore the cases N = 1, 2, 3, we have found the correct combinations of the general RG equations that gives the reduced (independent) system of β-functions. Our analysis reveals that only in d c = 8 3 and for N = 3 the reduced system of β-functions predicts a true hypercubic fixed-point. All the other cases, instead, revealed no fixed-point solution associated with a truly hypercubic symmetry, but only solutions corresponding to multicritical Ising and O(N ) universality classes. In particular, in the models that we considered there is no fixed point with the symmetries of a square and N = 2.
For completeness, we reviewed the case of d c = 4 (in appendix A) for which, as a new result, we give the NNLO contributions to the A-function in appendix A. We have take the opportunity to discuss the behavior of the quadratic deformations at the hypercubic fixed point, commenting explicitly the physical mechanism occurring when two deformations coincide (which is known to be related to logarithmic CFTs [29,40,41]). Finally, the inclusion of a singlet in the fields' multiplet allows constructing the interactions with an odd number of fields. In this case, we limited ourselves to report the beta functions of the φ 5 -hypercubic theory in d c = 10/3 given in appendix B, which are relevant for the discussion of a generalization of the random universality class presented in Ref. [29].
A future prospect of our approach could be to attempt the construction of multicritical M N -model, that are theories with O(m) n S n symmetry, which can be seen as a natural generalization of hypercubic symmetry [42][43][44]. It is especially interesting to attempt this generalization because M N -models are seeing a resurging interest, mostly coming from the conformal bootstrap literature, which is motivated by the possible presence of two distinct universality classes that could have direct experimental implications [45,46].

A Critical models in d = 4 −
In this appendix we review the known case of the hypercubic model in d = 4 − . The critical properties of the system below d c = 4 can be described in terms of the following potential [1,2,47], which can be obtained following the logic discussed in Section 2. The renormalized model predicts four fixed points: the Gaussian one, the Ising one in which the N components of the field decouple (therefore referred to as decoupled Ising), the O(N )-symmetric and the hypercubic fixed point. The three non-trivial fixed points have been extensively studied in the RG literature: see for example Refs. [3,4,48,49] for present-day six loops order results within the -expansion. There has been special interest in the determination of the stability properties of the aforementioned fixed points: it was shown that while the Gaussian and Ising fixed points are always unstable for arbitrary values of N , the O(N )-symmetric, referred to as isotropic, and the cubic, referred to anisotropic, compete with each other and their respective stability depends on N . For N < N c , where N c is some critical value that depends on d, the isotropic fixed point is stable while for N > N c the opposite is true. The long debated issue is the value of N c for a three-dimensional system, d = 3. To our best knowledge, the reference field-theoretical estimate is still provided by six-loop results of Ref. [3] and recently Ref. [4] which predict N c 2.9 and strengthen the existing arguments in favor of the stability of the cubic fixed point in the physical case N = 3. We mention here that, interestingly, the authors of Refs. [50][51][52] claim evidence of a new 3d CFT characterized by the cubic symmetry group which is different from the one obtainable by standard -expansion starting from d = 4 − . This fact certainly deserves attention and could possibly be investigated by means of the functional non-perturbative RG extending the analysis carried out in Ref. [53]. The determination of N c , as well as the assessment of the existence of this new cubic 3d CFT, are, however, not the scope of the present paper.
In the functional perturbative RG, the LO, NLO and NNLO contributions to the potential's β-functional can be expressed diagrammatically as By plugging the effective potential (A.1) into Eq. (A.2) and switching to dimensionless variables, one can read off the corresponding β-functions β m 2 , β λ 1 and β λ 2 , the renormal-ization group flow of the A function and that of the anomalous dimension matrix γ ij . The procedure is straightforward but requires some computational power.
The corresponding critical exponents and A function read (A. 16) There are two important things to notice. The first one is that the limit N → 0 of (A.15) is divergent, because of inverse powers of N . The reason why this happens is wellknown and it has to do with the fact that for N = 0 the beta functions β λ 1 and β λ 1 have a degenerate solution, so it is necessary to consider a subleading equation when studying the -expansion. This, in turn, makes it such that the expansion is not analytic in , but rather in 1 /2 , and the solution is known as Khmelnitskii fixed point [54]. Obviously, the unpleasant effect is that the expansion converges much less for finite values of , but it is still very important to see the fixed point because it has the physical meaning of describing an Ising model in a lattice that has random impurities or that is randomly diluted, according to some simple distribution [55].
The second important point can be made in relation to the stability and N c at d = 3. Specifically, as seen from Fig. 3 we notice that the critical exponents of the X and Y quadratic operators of the irreps of H N almost coincide for N = 3, implying that axial and quadratic deformations scale very similarly in that configuration. This has to be compared with the analysis of the irreps of the O(N ) model, which contain the same singlet S as the hypercubic model, but combine the X and Y operators in a single traceless tensor with two indices. This implies that if γ X = γ Y then the hypercubic model must coincide with the O(N ) model, so, from our extrapolation to = 1 we deduce that for N = d = 3 the two models must be very close in some sense. In fact, one can check explicitly that the perturbative solution of γ X = γ Y , that gives N c = N c ( ) in d = 4 − , has the same solution as the condition for which the two fixed points (A.13) and (A.15) coincide. It is actually even simpler to determine N c ( ) in this way. . We see that at N = 2, the singlet coincide with the diagonal quadratic deformation. This is expected since at N = 2 the cubic behavior reduces to two decoupled Ising models and the diagonal deformation simply amounts at a shift in the critical temperature (this holds true at all orders in perturbation theory). At N ≈ 3 the axial and diagonal quadratic deformations almost coincide. This is expected too since at N = N c the cubic and isotropic fixed points coincide and the corresponding two crossovers are equivalent, and N c is known to be close to 3. Notice, however, that this extrapolation gives N c 3, while it is believed that N c 2.9 3, and the estimate is fixed by a careful resummation. The analytical expression of the crossover exponents Φ axis and Φ diag can be checked at order O( 2 ) against Ref. [20].