Degenerate manifolds, helimagnets, and multi-$\mathbf{Q}$ chiral phases in the classical Heisenberg antiferromagnet on the face-centered-cubic lattice

We present a detailed study of the ground state phase diagram of the classical frustrated Heisenberg model on the face-centered-cubic lattice. By considering exchange interactions up till third nearest neighbors, we find commensurate, helimagnetic, as well as noncollinear multi-$\mathbf{Q}$ orders which include noncoplanar and chiral spin structures. We reveal the presence of subextensively degenerate manifolds that appear at triple points and certain phase boundaries in the phase diagram. We show that within these manifolds, the spin Hamiltonian can be recast as a complete square of spins on finite motifs, permitting us to identify families of exact ground state spin configurations in real space -- these include randomly stacked ferro- or antiferromagnetically ordered planes and interacting ferromagnetic chains, among others. Finally, we critically investigate the ramifications of our findings on the example of the Ising model, where we exactly enumerate all the states numerically for finite clusters.


I. INTRODUCTION
A first acquaintance with the concept of geometric frustration, which has played such a pivotal role in modern condensed matter physics, is generally made in the context of a classic textbook example of the crystal structure of salt (NaCl), namely, the face centered cubic (fcc) lattice. The elementary motif of any covering of the fcc lattice is a triangle, and this feature renders it impossible for antiferromagnetically (AF) interacting spins to simultaneously satisfy all interactions. Indeed, a system of AF interacting spins treated as n-component classical (spin S → ∞) vectors of a fixed length forms an infinitely degenerate one-dimensional ground state manifold at zero temperature [1]. The investigation of the thermodynamic and critical behavior of the fcc O(n) antiferromagnet has had a long history and been the subject of much debate, especially for the n = 1 (Ising) model [2][3][4][5][6] which is now known to undergo a first order transition into a collinear AF state [7][8][9][10][11][12][13], while on the other hand, hardly much is known about the classical n = 2 (XY ) model [1,14]. Concerning the physically realizable case of n = 3 (Heisenberg) spins, after much debate [14][15][16][17], there now appears to be a consensus that the model undergoes a first order phase transition into a collinear AF state [18]. For Heisenberg spins away from the classical limit, i.e., when 1/S = 0, only few studies have addressed the role of quantum fluctuations in the large-S [19] or small-S limits [20], and nearly seven decades after being first attended to [21], the determination of the nature of the ground state of the quantum Heisenberg antiferromagnet on the fcc lattice remains a critically outstanding problem.
Recently, it has been realized that the subextensive degeneracy of the T = 0 ground state manifold is not unique to first neighbor (J 1 ) antiferromagnetic interactions and that upon inclusion of second neighbor interactions (J 2 ) a twodimensional subextensively degenerate ground state manifold in the form of a spiral surface can be stabilized for J 2 /J 1 = 1/2 [22], in addition to three different AF commensurate orders [23]. In this work, we incorporate a third neighbor ex-change coupling J 3 and obtain the T = 0 global phase diagram of the classical J 1 -J 2 -J 3 Heisenberg model for all possible combinations of signs of couplings, which has hitherto never been investigated. Our investigation reveals a rich phase diagram featuring helimagnetic and noncollinear multiple-Q orders which allow for noncoplanar and chiral structures. The most salient finding of our study is that the phase diagram is host to highly (subextensive) degenerate one-and twodimensional ground state manifolds in q-space occurring at triple points (where several phases meet) and phase boundaries. Remarkably, at these triple points and phase boundaries, we are able to express the Hamiltonian as a positive definite sum of complete squares on finite motifs covering the lattice allowing one to understand the origin of the subextensively degenerate manifold of states. This reformulation also permits us to explicitly construct large classes of nontrivial, aperiodic ground states in real space, consisting of randomly stacked ordered planes and frustrated ferromagnetically ordered chains in special crystallographic directions. Considering the case of Ising spins, we are able to completely enumerate, on finite clusters, the type of possible configurations in real and Fourier space, and provide indications of an even richer structure for Heisenberg spins. Our work also provides the basis for understanding the origin of a plethora of fcc magnetic structures in wide a variety of magnetic materials [24]. In particular, a recent investigation of the half-Heusler compound GdPtBi [25] proposed an antiferromagnetic J 1 -J 2 -J 3 Heisenberg model, and argued for the indispensability of a J 3 interaction to match the observed neutron scattering profile.

II. STRUCTURE OF THE PAPER
The remainder of the article is structured as follows: We introduce the fcc lattice, the Heisenberg model and the Luttinger-Tisza method in Sec. III. The phase diagram is derived in Sec. IV. Sec. V presents the commensurate phases and the construction of the possible multiple-Q structures. We shortly discuss the appearing incommensurate phases in Sec. VI. We describe the details of the construction of the spin structures in real and reciprocal space in the ground-state manifolds in Sec. VII. Finally, in Sec. VIII we conclude with a summary of the results. The article ends with the following (mostly technical) Appendices. In Sec. A we define our conventions for the lattice and show the Fourier transform of the interactions. Sec. B contains a table of the phase boundaries. The Ising configurations in the different subextensive manifolds are enumerated in Sec. C for finite clusters and critically compared to results of Sec. VII.

III. THE MODEL AND THE LUTTINGER-TISZA METHOD
The face-centered-cubic (fcc) lattice is an archetypal frustrated lattice, it can be built from (111) triangular planes, in an ABCABC type stacking style, see Fig. 1(b). A triangular lattice is frequently called hexagonal to emphasize its sixfold symmetry. This frustration is even more emphasized if we build the lattice from edge sharing tetrahedra, the two differently oriented tetrahedral building blocks are depicted in Fig. 1(c). Another way of constructing the lattice is an edge sharing octahedral covering, for a picture of the octahedral building blocks see Fig. 1(d).
The Hamiltonian of the classical isotropic Heisenberg describes three dimensional unit vectors |S i | = 1 at the sites R i of the fcc lattice, interacting with first (J 1 ), second (J 2 ), and third neighbor (J 3 ) interactions. The summation indices i, j δ with δ = 1, 2, 3 refer to the δ'th neighbor pairs. There are twelve first, six second and twenty-four third neighbor vectors in the fcc structure. One vector of each neighbor set is presented in Fig. 1(b). We wish to find the ground states of the model Eq. (1) using the method developed by Luttinger and Tisza [26], i.e. by finding the minimum of the exchange interaction in Fourier space J(q). We define the energy in reciprocal space as where the summation runs over the Brillouin-zone (BZ), N is the total number of sites of the lattice with periodic boundary conditions, and is the Fourier transform of the exchange interactions with lattice separation vectors δ = R i − R j , which is presented in Eq. (A4). We used the following convention for the Fourier transform of the spins: In the Luttinger-Tisza method we minimize J (q) with respect to q, i.e. we solve the gradient equation for α = x, y, z. If the equation is satisfied at a point q = Q we check the positive semidefiniteness of the Hessian ∂ 2 J(q) ∂q α ∂q β q=Q 0.
The conditions above are necessary, but not sufficient to have a global minimum: to find a true ground state we have to compare the different local minima and choose the lowest one. We denote the set of points {Q} -the ordering vectorswhere J(q) takes its minimal value by M GS , and call it as the ground state manifold, and choose the Fourier amplitude vectors S Qs to satisfy the spin-length constraint |S i | = 1 for every site. For a given M GS the amplitudes have to be chosen carefully, and we give a detailed analysis of the amplitudes and the corresponding orderings in real space in Sec. V. The ground state energy per site ε(Q) equals to the Fourier transform of the exchange constant evaluated on the M GS . To show this we evaluate Eq. (2) on the M GS :  Fig. 3. About notation: we refer to points in the Brillouin zone either by their names and and coordinates in units of 2π or their respective wave-vector, i.e. QW = (2π, π, 0) ≡ W (1, 1 2 , 0).
The ε(Q Σ = (q Σ , q Σ , 0)) Fourier transform is: where the 0 indices refer to the ground state properties, note that ε(Q) only depends on the ordering vector parametrically. We use the reality of the spin components: S 0 −Q = S 0 * Q , and that J (Q) is constant on M GS and we employ the spin length constraint q |S q | 2 = 1. The latter is true for any state satis- The dimension of the M GS is a crucial ingredient in understanding the physics of these systems. Conventional commensurate ferromagnets or antiferromagnets correspond to zero dimensional manifolds, i.e. a handful of points in highly symmetric positions of the BZ. Incommensurate orders (spin spirals, helices, cycloids) [27] also have zero dimensional M GS , but now at generic points in the BZ. Magnetic Bragg peaks show up at these points in neutron scattering. But in frustrated systems the possibilities are richer: the M GS can be a one-dimensional degenerate manifold as for the J 1 -J 2 model on the square [28,29] and honeycomb lattices [30][31][32] and the J 1 only model on the fcc lattice [19,33] where every point on a line is a possible ordering vector, and with carefully chosen amplitudes we can compose ground states with complicated spatial variation. The situation can be even more complex when the dimension of the M GS is larger: spin spiral surfaces (i.e. two-dimensional M GS -s) were found in J 1 -J 2 Heisenberg model on the diamond [34], fcc [16,22], and body-centered cubic [35] lattices. Furthermore, on the kagome [36,37] and pyrochlore [38][39][40] lattices the whole BZ is the M GS . These extended manifolds give an opportunity to the system to fluctuate between the degenerate ground states making them candidates for classical spin liquids in some temperature range [34,[41][42][43][44].
There is another type of degeneracy, even in unfrustrated spin models without extended M GS -s (i.e. simple cubic ferro-or antiferromagnets of first neighbor couplings) where the interactions are isotropic: breaking the global O(3) rotational symmetry of the model results a family of ground states that can be rotated globally to each other in spin space, resulting a degeneracy parametrized by the three dimensional group of rotations (we will refer to usually called as a trivial degeneracy, since its presence is independent of frustration). Besides the trivial degeneracy, ground states having multiple sublattices can still be indeterminate: we can continuously deform them to each other by a set of local rotations [45][46][47]. This leftover degeneracy can again be characterized by a (continuous or discrete) set of parameters, and we will give a detailed analysis of this scenario in cases of commensurate orders in our model. In the following we construct the ground state phase diagram of the model in exchange parameter space by minimizing the energy ε(q) with respect to the wave-vector.

IV. CLASSICAL PHASE DIAGRAM AND ORDERING VECTORS
In this section we construct the classical, zero temperature, ground state phase diagram of the model Eq. (1) in the J 1 − J 2 − J 3 parameter space. We compare the ε(Q) values for the possible orderings and choose the lowest one for a given set of parameters, these results are collected in Table I, and we present the detailed phase diagram in Fig. 2.
We analyze the phase boundaries by comparing the ground state energies of the neighboring phases summarized in Table VI. The phase boundaries are of second order if the ordering vectors of the two matching phases can be deformed continuously into each other, and of first order if the transition  Table I. Solid black lines mark first order phase transitions, dashed black lines stand for continuous (second order) transitions, and the equations describing these boundaries are collected in Table VI. We label the phases by their ordering vectors given in units of 2π/a, as presented in Fig. 3. The four commensurate phases are the ferromagnet Γ(0, 0, 0), and the three types of antiferromagnets: X(1, 0, 0), L 1 2 , 1 2 , 1 2 and W 1, 1 2 , 0 . The commensurate ordering vectors are depicted in Fig. 3(a), note that all these phases are already present in the absence of J3 [46,48]. Introducing a finite J3 introduces the incommensurate phases ∆(q, 0, 0), Λ(q, q, q) and Σ(q, q, 0), where q has to be optimized according to Eqs. (46a)-(46c), and the possible ordering vectors are depicted in Fig. 3(d)-(f). The phase Σ(q, q, 0) in (b) has a bow-tie shape (enlarged in the inset) with a neck consisting of a single point J2 = J1/2 and J3 = 0 (the green dot), and through this point the L 1 2 , 1 2 , 1 2 and W 1, 1 2 , 0 phases meet. The dark red X − W phase boundary emanating from the first-neighbor antiferromagnetic point J2 = J3 = 0 is degenerate: along this line any of the ground states with ordering vectors residing on the one-dimensional manifold M 1 Z defined by Q = (2π, q, 0) with q ∈ [−π, π] shown in Fig. 3 Fig. 3(e)]. The green dot at J3 = 0, J2 = J1/2 > 0 is a phase of even higher degeneracy: the ground states form the two-dimensional manifold shown in Fig. 3(c). Basic properties of these manifolds are collected in Table II. requires a discontinuous jump of the ordering vector. There are special points of the phase diagram: the triple points and the X(1, 0, 0) − W (1, 1 2 , 0) phase boundary that require particular attention: at these points ground state manifolds extend to lines and a surface, signaling a large but subextensive degeneracy of the ground states. A sidenote about notation: we use two sets of notation for the BZ points: we refer to points in the BZ either by their names and and coordinates in units of 2π or their respective wave-vector, i.e., Q W = (2π, π, 0) = W (1, 1 2 , 0). Basic information about these phases is collected in Table II. In the following we briefly summarize the properties of the phase diagram, but details of the real space picture of the orders are given in Sec. V.
At first glance the spin configurations can be either collinear, coplanar or non-coplanar. Incommensurate spirals are coplanar but not collinear. Commensurate phases can be non-collinear, or even non-coplanar if we use multiple arms of the star of the ordering vector, these are the mutiple-Q phases. Order by disorder (either quantum or thermal) mechanisms tend to select the collinear configurations [46], but disorder and ring-exchange mechanisms favor noncollinearity [47], so the fate of these states depends on further details. In the following sections we give a detailed analysis of the possible configurations. Table II. Multiple points of the phase diagram, corresponding to degenerate manifolds in q-space. In the first column we give the label of the manifold, see Fig. 2. The manifolds themselves are depicted in Figs. 3(b)-(e). In the second column we list phases that meet at the special parameter values given in the third column. In the fourth column the dimension of the manifolds is given, together with their defining equation in q-space, when given in parametric form we mention only one of the crystallographically equivalent directions. The last column gives the energy per site on the manifolds.

Label
Touching phases  Fig. 2(b). Every point on the crisscrosses is energetically degenerate. (c) Two-dimensional energetically degenerate manifold M 2 (a Schwarz P surface [22,51]) corresponding to J3 = 0, J2 = J1/2 > 0 (green dot in Fig. 2(b)). The second row ((d)-(f)) corresponds to incommensurate orderings, in these figures -depending on the exchange parameter values (c.f. Eq. (46a)-(46c))-a single ±Q pair of wave-vectors is chosen as the ordering vector of the developing spin spiral. This row also depicts the one dimensional degenerate manifolds of the special points of the phase diagram. (d) Incommensurate ordering vectors ∆(q, 0, 0) corresponding to spin spirals propagating in the directions of the cubic axes with 6 arms. This is also the manifold M 1 ∆ of the Γ − ∆ − X triple point, see the yellow dot in Fig. 2(a). (e) Incommensurate ordering vectors Λ(q, q, q) corresponding to spin spirals propagating in the directions of the body diagonals of the cubic cell with 8 arms. This is also the manifold M 1 Λ of the Γ − Λ − L triple point, see the orange dot in Fig. 2(a). (f) Incommensurate ordering vectors Σ(q, q, 0) corresponding to spin spirals propagating in the directions of the face diagonals of the cubic cell (only pictured in the horizontal planes for better visibility, this vector has 12 arms). About notation: we refer to points in the Brillouin zone either by their names and coordinates in units of 2π (fractional coordinates) or their respective wave-vector, e.g., W (1, 1 2 , 0) ≡ QW = (2π, π, 0).

V. COMMENSURATE ORDERINGS AND REAL SPACE DESCRIPTION
In this section we describe and analyze the developing orders in detail. For the four commensurate orderings (c.f. Fig. 3(a)) we calculate the Fourier amplitudes, give the constraints on them, count the degeneracies (the number of free parameters describing the order that remain after removing the trivial, global O (3) of the symmetry breaking). We describe the orders in real space and and analyze their symmetry properties. In order to find the number of wave-vectors participating in a given order and to find the constraints on their Fourier amplitudes we must not distinguish between the equivalent q-vectors, where by equivalence we mean differing only in some reciprocal lattice vector G, i.e., q ∼ q if q = q + G. Constraints on the Fourier amplitudes and the number of free parameters can be calculated as follows: we expand the spins S i in Fourier space, keeping only the amplitudes of the arms of the star of the respective ordering vector finite. Afterwards we impose the constraints that for every lattice point the spins have to be real unit vectors, and solve the equations [48] for the Fourier amplitudes: This phase is an ordinary ferromagnet, where all the spins align and only the trivial O(3) degeneracy is present.
B. The X (1, 0, 0) antiferromagnet, Type-I In this phase the nonequivalent Q-vectors form a threearmed star: Since these arms reside on the midpoints of the square-shaped faces of the BZ [52] we can mix them to construct a triple-Q order, provided we choose the Fourier amplitudes appropriately. In order to make the spins real unit vectors we need to consider some constraints on the complex amplitudes S 0 Xα : Since X α and −X α are equivalent (X α ∼ −X α ) and the phase factors e −ıXα·Ri are simply ±1-s, the amplitudes have to be real to ensure the reality of S i : To fix the lengths of the spins we have the following four constraints: The three real amplitudes S 0 Xα mean 9 free parameters. The global O(3) freedom removes 3 of them, and together with the four constraints in Eq. (13) we are left with two free parameters to characterize the degeneracy [46].
The global O(3) freedom of the symmetry breaking and the mutual orthogonality and normalization of the S 0 Xα -s allows as to parametrize them as: where all the parameters are real, and they satisfy the additional constraint: ξ 2 +η 2 +ζ 2 = 1. The ground state manifold thus can be parametrized by a unit vector (ξ, η, ζ). In the following we construct and analyze the developing order in real space.
With the parametrization of Eq. (14) the spin sitting on the lattice point R i = (x, y, z) is (we recall that the coordinates can be integers or half-integers): The superlattice vectors of this order form a simple cubic lattice with the same unit cell as the conventional cell of the original fcc lattice, and the four sublattices form tetrahedra with spins see Fig. 4(a). The spins on every elementary tetrahedron sum up to zero: we refer to this situation as the tetrahedron rule. We can check our Fourier space degeneracy counting in real space: the four sublattice spins mean four unit vectors as eight free parameters, and the global O(3) symmetry and also the tetrahedron rule remove three of them, leaving only two free parameters as expected.
If we use all three arms of the star, creating a triple-Q order (with ξ, η, ζ finite) the spins are noncoplanar: they possess a finite scalar chirality S 1 · (S 2 × S 3 ) on all the triangular faces of each tetrahedron. For -say-ζ = 0, the configuration is coplanar, and if only ξ remains, it is collinear. Thermal or quantum order by disorder effects select a single arm of the star resulting in a collinear structure [46]. There is a possibility, that this type of order is chiral [53] in the sense that applying space inversion I (that leaves the spin-space invariant) centered at some point transforms the spin configuration into another one, that cannot be reached from the original configuration by any lattice translation. This is not the case: the effect of inversion (centered at a lattice point, or at the midpoint of any bond) can be compensated by translating the lattice with a lattice vector: this state is achiral.
In this phase the nonequivalent Q-vectors form a fourarmed star: Since they are on the midpoints of the hexagonal BZ faces we can construct a quadruple-Q order out of them [52]. We expand the spins in Fourier amplitudes: Since L α ∼ −L α and the phase factors are e −iLα·Ri = ±1, the amplitudes have to be real to ensure the reality of S i : We can express the spins in real space as: We can substitute the lattice points in the above equation that yield four independent sublattices on an elementary tetrahedron. The tetrahedron rule does not hold. Shifting this tetrahedron by δr = (1, 0, 0) reverses the directions of the spins, so we have an 8-sublattice antiferromagnet of spin pairs S A , S B , S C and S D , and S A = −SĀ (shifted by (1, 0, 0)), and so on: This type of order is depicted in Fig. 4(a). We can expand the Fourier amplitudes as: The superlattice vectors form an fcc lattice doubled in linear size with respect to the original one, with primitive lattice vectors given by: It is easier to calculate the degeneracies in real space: the four independent unit sublattice spins mean 8 real degrees of freedom, the global O(3) removes 3 of them yielding 5 independent real parameters [46]. In each of the (111) triangular planes only four spins appear, forming a regular 4-sublattice order [54]. These phases are achiral (just like in the Type-I states): i.e., the phases that are related by inversion (either centered on lattice points or on bond midpoints) can be translated to each other by an appropriately chosen lattice vector.
There are 24 symmetry related vectors (the corners of the BZ) belonging to this type of order, and they fall into 6 equivalency classes so the star of W has 6 arms. This can be understood since each corner of the BZ is shared by four truncated octahedra. The arms of the star form three ± pairs: Q α = ±W 1 , ±W 2 , ±W 3 , the classes are: Since W α −W α to ensure reality of the spins in real space we have to combine the ±W α pairs: This type of ordering is also called a triple-Q one. We can decompose the complex amplitudes into real vectors: To fix the lengths of the spins we have the following constraints for the complex amplitudes: or we can express these constraints for the real and imaginary parts of amplitudes: The three pairs of real vectors u α and v α mean 18 free parameters, the equations above give 13 constraints, and we have the global O(3) degrees of freedom (3 free parameters), so we are left with 2 free real parameters for the degeneracy degrees of freedom (in perfect analogy with the Type-I phase).
Using the global O(3) freedom and the orthogonality and normalization of the u α -s we can parametrize them as: where all the parameters are real, and they satisfy the additional constraint: ξ 2 + η 2 + ζ 2 = 1. We use the orthogonality relations between the u α -s and v α -s (see Eq. (29c)-(29e)) to calculate the form of the v α -s (at this point any combination of signs is allowed, resulting 8 possible combinations): With this parametrization the spins become: where each ∓ corresponds to the respective ± in Eq. (31), i.e. if there is ± in front of ξ in Eq. (31), then there is a ∓ in the first component of S i , and so on. Actually these states are not all physically different, there are only two independent phases that form chiral partners (i.e. they are transformed to each other by space inversion). We are going to show this at the end of this subsection.
In order to understand this real space spin pattern we write the lattice point in Cartesian coordinates (R i = (x, y, z), i.e. in our fcc lattice the Cartesian coordinates x, y, and z can either be integers, or some half integer combinations). The spins are: This form is particularly useful for analyzing the properties of the spin structure: we only have to monitor the phase shifts due to transformations (either in real or spin space) and draw the consequences. From this form it can be seen that this pattern is periodic under a translation of 2 along the Cartesian directions (we have to change all x, y, and z by an even number to achieve a 2π phase shift in every component), and this is the smallest possible supercell, with primitive translations: The resulting superlattice is simple cubic, and the edge length of the superlattice cube is 2, and this supercell contains 32 points of the original fcc lattice. There are eight spin directions, that come in four ± pairs (like in the Type-II phase), but they form a complicated 32-sublattice order (we do not even try to visualize this spin pattern here). The four spin directions are not independent: there is still a tetrahedron rule in action (like in Eq. (17)): spins on every elementary tetrahedron sum up to zero. So out of 8 parameters describing the four sublattice spins the tetrahedron rule removes three parameters, and the global O(3) removes another three yielding the correct number of two parameters of the unit vector (ξ, η, ζ). We are left with the task to decide how many physically different phases the discrete ± parameters in Eq. (31) or equivalently the ∓-s in Eq. (33) yield.
We are going to prove that there are only two nonequivalent phases, transformed into each other by space inversion centered on a lattice point (therefore we call them chiral partners or enantiomorphic domains), and those phases cannot be transformed to each other by any lattice translation, i.e. they are really physically different. Let us denote a state in Eq. (33) by their respective sign triplet, i.e. (mmm) and (ppp) have signs (− − −) and (+ + +). We want to study the effect of space inversion (I) and translations t δr by a lattice vector δr = (δx, δy, δz) on the phases of Eq. (33): ϕ = (ϕ x , ϕ y , ϕ z ) , δr = (δx, δy, δz) .
Inversion changes (mmm) → (ppp), so its effect is implemented by a π/2 phase shift in every component of Eq. (33): while a translation acts the on the phases the following way: States that have an odd number of m-s (i.e. (mmm), (mpp), (pmp) and (ppm)), can be translated to each other, so can the ones with an even number of m-s (i.e. (ppp), (pmm), (mpm) and (mmp)) be translated to each other. Applying inversion changes this parity, but translations never. In order to show this on one hand consider the sum and its change under inversion below: on the other hand its change under a translation Eq. (41) is: and this phase change is always an integer multiple of π. So we have showed that states of different parity form disjoint sets (they can never be transformed into each other by translation), we now have to prove, that within each set we can find a translation that transforms the states into each other, i.e. that we only have two physically different states. These translation vectors δr for the (mmm) class states are:

VI. INCOMMENSURATE PHASES
All the commensurate phases were found in the J 1 − J 2 models [23,46,48], see the J 3 = 0 lines in Fig. 2(a) and (b). A novel feature of the J 3 = 0 model is the appearence of incommensurate orderings with propagation vectors along special, highly symmetric directions [55]. For the three incommensurate orderings (c.f. Fig. 3(d)-(e)) we give the dependence of the wave-vectors on the exchange parameters and give their accessible ranges.
In all these phases -since the ordering vectors Q are incommensurate-we have to include both ±Q to ensure reality of the spin components, the resulting spin pattern becomes: where s 1 and s 2 are arbitrary orthogonal unit vectors spanning the plane of spin rotations, and ϕ is an arbitrary phase, and ± accounts for the two possible chiralities. Since our model is isotropic, there is nothing to fix the plane of rotation to the wave-vectors or to the crystallographic axes (another manifestation of the spontaneous breaking of the global O(3) symmetry). Since the wave-vectors are incommensurate we cannot build any multiple-Q ground states of their stars [52]. The three possible incommensurate ordering directions are ∆ (q, 0, 0) with its 6-armed star (see Fig. 3(d)), Λ (q, q, q) with its 8-armed star (see Fig. 3(e)) and Σ (q, q, 0) with its 12-armed star (see Fig. 3(f)).
The vector Q = (π, π, 0), although not a special point in the BZ is compatible with a Nel-type commensurate antiferromegnetic ordering, called the Type-IV phase fcc antiferromagnet, realized in CoN [19,56]. Unfortunately the possible q-values of the Σ (q, q, 0) phase are far from π, and we could not even stabilize this type of Nel-order by introducing quantum fluctuations (in the spirit of [57]).
The optimized q-values of the incommensurate ordering vectors are given by:

VII. GROUND STATES OF THE EXTENDED MANIFOLDS
For the pure fcc model (J 1 > 0, J 2 = J 3 = 0) the degenerate manifold M 1 Z has already been found [33,58], and a class of ground states were constructed from (100)-directed, noninteracting AFM planes. In this section we describe the phases with extended M GS -s of energetically degenerate ordering vectors (see Fig. 3(b)-(e)) that correspond to large ground state degeneracies at special points in the phase diagram in Fig. 2. We explain these degeneracies by a real space construction of covering the lattice with finite motifs (see Fig. 5), and write the Hamiltonian as a positive definite sum over these motifs. Minimizing the Hamiltonian imposes local constraints on the spins on these motifs: any state that satisfies these constraints is an allowed ground state. We extend the construction presented in [33] for the other degenerate manifolds and construct ground states from noninteracting planes, and also find Table III. Motifs used to cover the lattice (see Fig. 5) together with their symbols used in formulas, and their overcounting of the sites and of the J1 first, J2 second, and J3 third neighbor bonds. E.g. in a tetrahedral covering we put two differently oriented tetrahedra on each site, and as consequence each J1 bond is shared by two tetrahedra, and no longer bonds are covered. In the last column we give the reference as a subfigure for the picture of the motif in Fig. 5.

Motif
Symbol Site J1 J2 J3 Subfigure Tetrahedron tetra ground states consisting of ferromagnetic chains (though the chains are now interacting). We have solved the models for the extended manifolds for Ising spins S i ∈ {1, −1}, for finite, symmetric clusters. Details of these results are presented in Appendix C. We also have performed numerical simulations for planar (O(2) or XY ) spins S i = (S x , S y ) i , S 2 x + S 2 y = 1 to guide our intuition about the possible ground states. At the point J 1 = 2J 2 > 0, J 3 = 0 (the green dot in Fig. 2(b)) the ordering vectors of the possible ground states form the two dimensional M 2 manifold, defined by as depicted in Fig. 3(c). This is the only point of the phase diagram with such a large degeneracy [22]. The W 1, 1 2 , 0 and L 1 2 , 1 2 , 1 2 points are parts of this manifold, and this is the point of the phase diagram where the W 1, 1 2 , 0 and L 1 2 , 1 2 , 1 2 phases meet through the neck of the Σ(q, q, 0) phase. At this point the Hamiltonian reads: We can express this Hamiltonian as a sum of complete squares of spins forming edge-sharing octahedra covering the lattice (see Fig. 5(e)): where S 1 , . . . , S 6 refer to the six spins on the sites of an octahedron. Since every first neighbor bond is covered twice, and every second neighbor bond once (see Table III), Eq. (49) exactly reproduces Eq. (48), and this is why we have chosen the octahedral covering for these particular values of exchange parameters. Since J 1 > 0 Eq. (49) is minimized if the spins sum up to zero, on every octahedron -we refer to this rule as the octahedron rule. Every such configuration is a ground state, and every ground state has this property. The additional constant − 3 2 J 1 N gives the ground state energy.
• Substituting the vertices of any octahedron in Eq. (33) of the W 1, 1 2 , 0 phases, and using Eq. (41) to calculate the phase shifts while we walk around the vertices of the octaheron the sum of the spins equals zero.
• A general spiral with Q ∈ M 2 also satisfies the octahedron rule: this can be checked by putting an arbitrary Q in Eq. (45), and summing up the spins on octahedra: the sum vanishes if and only if Q satisfies the defining equation (47) of M 2 .
Order by disorder effects (either thermal or quantum) at the harmonic level select the L 1 2 , 1 2 , 1 2 points on the M 2 surface [22].
At the triple point J 2 = −J 1 > 0, J 3 = 0 (the orange dot in in the phase diagram Fig. 2(a)) the M GS is M 1 Λ . for the degenerate manifold. The possible ordering vectors Λ(q, q, q) smoothly interpolate between Γ(0, 0, 0) and L 1 2 , 1 2 , 1 2 , hence the shape of the manifold, see Fig. 3(e). The Hamiltonian reads We can cover the lattice by signed squares, with signs distributed according to Fig. 5(d): to every site we associate 3 squares, lying in each of the {100} planes. This way every first and second neighbor bond is covered twice, see Table III. This Hamiltonian is minimized if and only if the sums vanish on every square, and the ground state energy per site is given by the additional constant ε = +3J 1 . The three equations on the signed squares are not independent, and instead of them we can use the octahedra containing (a) (b) (c) (d) (e) Figure 5. Finite motifs used to cover the fcc lattice in the exchange parameter regions with high ground state degeneracies. Red, green and blue lines denote first, second, and third neighbor bonds of a motif, respectively. By building up the crystal from these motifs we cover bonds and sites multiple times, this overcounting is summarized for every motif in Table III. (a) Elementary tetrahedra (index "tetra" in formulas).
From each lattice point we draw two differently oriented tetrahedra to cover every first neighbor bond twice. (b) A signed rectangle (index "rect1" in formulas): with 6 differently oriented rectangles put on every site we can cover every first, second and third neighbor bond. "Signed" here means that when writing the complete squares of the spin sums in the Hamiltonian we have to assign a minus sign to the spins sitting in the vertices denoted by white dots, black dots get a plus sign, see Eq. (68). Together with the tetrahedra we use this motif to construct the ground states of the phase corresponding to the manifold M 1 Z , see Fig. 3 60) for the chains are denoted by K1 for the first neighbor red bonds, and by K2 for the second neighbor green bonds. The gray rhombus is the projection of one of the covering signed squares also depicted in Fig. 5(d), minus signs have to be associated to one pair of opposite vertices, say to A and A . (c) Brillouin zone of the lattice depicted in Fig. 6(b), together with the ground state manifold (orange cross) of the Hamiltonian Eq. (60), this manifold is the section of M 1 Λ (see Fig. 3(e)) with the (110) q-plane passing through the origin, this BZ is not a perfect hexagon. Symmetry points of the original three dimensional BZ (see Fig. 3(a)) are shown, together with some less commonly known points K 3 4 , 3 4 , 0 and U 1 4 , 1 4 , 1 .
these squares to cover the lattice (see Fig. 1(d) and Fig. 6(a)). Out of the 3 square equations on orthogonal squares only two are independent per octahedron. Using the notations of Fig. 6(a) for the sites of the octahedron, the ground state spin configuration shall satisfy the equations: where m is proportional to the magnetization of an octahedron. We can solve them introducing the a, b, and c vectors The length constraint |S A | 2 = |S A | 2 = 1 becomes (m ± a) · (m ± a) = 1, and similar equations for b and c hold. Adding and subtracting these equations, we get These equations refer to every single octahedron, but we omit the octahedron index for clarity. For N s component spins the a, b, c, and m counts 4N s degrees of freedom, and there are 6 constraints in Eq. (56), so we expect 4N s − 6 free continuous parameters to describe the ground state of an octahedron. A ferromagnetic order trivially satisfies the rule given by Eq. (54), and this is not surprising: the Γ(0, 0, 0) point is part of this manifold. In the ferromagnet a = b = c = 0. If m = 0 we get an L 1 2 , 1 2 , 1 2 order: all the Type II states described in Sec. V C and shown in Fig. 4(b) can be constructed this way . Among others one can choose the single ordering vector (π, π, π) and get a set of alternating (111) ferromagnetic planes: see Eq.   Table IV). The gray parallelogram shows the projection of a covering square also depicted in Fig. 5 suggests other possible candidate ground states: we can try to construct a family of ground states by stacking ferromagnetic (111) planes, these planes form triangular lattices and they are depicted in Fig. 7(a).
Assuming a state consisting of ferromagnetically ordered (111) planes, and representing a plane by a single effective spin s i of unit length, where now the index "i" enumerates the consecutive planes one can derive an effective onedimensional model: where the effective exchanges can be inferred either from Table IV or from Fig. 7(a), L (111) is the number of (111) planes in the crystal, and the additional constant derives from the in- Table IV. Number of bonds connecting a single point to its neighbors on the nearby (111) planes, see Fig. 1(b) and especially Fig. 7(a). The first column gives the separation of consecutive planes: "0" means the (111) plane containing the chosen point, "1" means the two first neighbor (111) planes, "2" means the two second neighbor and J 3 = 0, we see that the first term disappears, so the planes disentangle, and the second term gives +3J 1 for the correct ground state energy per site of the original model. The resulting ground state is of the form F 1 F 2 F 3 F 4 . . . , where F i denotes the independent ferromagnetic planes. This independence of planes can be further rationalized by noting that only the first neighbor planes are connected by the covering squares (see Fig. 7(a)). Using the notation of Fig. 7(a) one can see that S 1 = S 1 and S 2 = S 2 , since these pairs lie on FM planes. Therefore Eq. (54) is automatically satisfied since S 1 + S 2 = S 1 + S 2 . Such a state can be cooked up by choosing ordering vectors solely from the (q, q, q) line of the M 1 Λ manifold: and depending on the complexity of the real space pattern, any symmetric set of points on the Λ(q, q, q) line can be present in the expansion, as long as we care about the choice of the Fourier amplitudes to produce a real space pattern of unit length spins. Of course we could have chosen any of the symmetry related 111 directions. In the finite cluster Ising solution (see Appendix C for details) we have found the {111} stacking of independent FM planes: these involve only one line of ordering vectors of M 1 Λ . We have found another type of solution where up and down spins form two interpenetrating pyrochlore lattices: the unit cell consists of 8 sites, see In the numerical simulations on planar spins we found another interesting class of ground states: a 3/4 majority fraction of the spins on (111) planes ordered ferromagnetically on a kagome sublattice of the triangular layer (see the orange and green dots in Fig. 7(b)), and the minority spins (purple and blue dots) seemed to be independent of the majority spins, and a similar structure was formed on every (111) plane. In the following we use the notations of Fig. 7(b). We can exploit the octahedral constraint of Eq. (56): we assume the kagome-style ordering described above on consecutive planes indexed by 1, 2, 3,. . . , and monitor how the consequences of the constraints propagate as we step down on octahedra between the planes. We fix the majority spins S 1 and minority spins S 1 on the first plane. For Ising spins fixing the spins on layer "1" determines the spins on all of the consecutive planes, and the resulting pattern is the quadruple-Q order described in the previous paragraph. In the numerical study on the planar spins we have found ground states formed by seemingly independent ferromagnetic chains [59,60] lying in the 111 planes, pointing in one of the 110 directions, a set of such lines are depicted in Fig. 7(c). This numerical finding suggests the following strategy: we assume a FM ordering along the (110) chains (the bonds along the chain are J 1 < 0 ferromagnetic), and derive an effective two dimensional model where we substitute the chains by a single effective spin s i of unit length, where the index "i" refers to points of the lattice formed by the chains, for a picture of the lattice see Fig. 6(b). The effective interactions K δ (δ points to the neighboring chains) are in general very complicated (each point is connected to 16 others, usually by multiple bonds) but the actual exchange parameters (J 2 = −J 1 and J 3 = 0) come to help us and result in a remarkably simple set of nonvanishing effective exchanges: where the indices refer to the bonds depicted in Fig. 6(b). These values can be inferred from the gray rhombus in Fig. 6(b) depicting the projection of one of the covering squares, also shown in Fig. 5(d) (in order to have the correct effective exchanges one needs to take into account all the three differently oriented squares). This lattice is topologically equivalent to a first and second neighbor FM-AFM model with bond strengths K 2 = −K 1 /2 > 0 on the square lattice. The effective two dimensional model reads: where the additional constant derives from the couplings within a chain. This model is strongly frustrated having a codimension-one M GS , depicted in Fig. 6(c): this manifold is nothing but the section of M 1 Λ with the (110) q-plane passing through the origin (the q 110 = 0 plane with notation of Fig. 6(b) and (c)). The large ground state degeneracy can be further rationalized by noting that this Hamiltonian can be  written as a sum of complete squares on signed rhombi (see the gray rhombus in Fig. 6(b)): the resulting rhombus rule S A +S A −S B −S B = 0 is just the signed square rule inherited from the three dimensional problem. Any state obeying the rhombus rule is a ground state for the (110) chains, and this is consistent with the numerical finding of seemingly random chains in the XY -model that actually obey the rhombus rule.
• Almost independent ferromagnetic kagome sublattices in the {111} planes. This type of ordering is realized in the XY and O(3) models. For Ising spins this order reduces to the commensurate quadruple-L structure of intercalating pyrochlore lattices.

(d). The Hamiltonian reads
(1,0,0) (0,1,0)  Fig. 3(b)) with the (100) q-plane passing through the origin. Symmetry points of the original three dimensional BZ (see Fig. 3(a)) are shown. (c) Brillouin zone of the lattice depicted in Fig. 9(a), together with the ground state manifold (dark yellow cross) of the Hamiltonian Eq. (65), at the Γ − ∆ − X triple point. This manifold is the section of M 1 ∆ (see Fig. 3(d)) with the (100) q-plane passing through the origin. Symmetry points of the original three dimensional BZ (see Fig. 3(a)) are shown.
We can cover the lattice by signed rectangles where the signs are distributed according to Fig. 5(c):  Table V. Number of bonds connecting a single point to its neighbors on the nearby (100) planes, see Fig. 1(b) and especially Fig. 8. The first column gives the separation of consecutive planes: "0" means the (100) plane containing the chosen point, "1" means the two first neighbor (100) planes (see Fig. 8(a)), "2" means the two second neighbor (100) planes (see Fig. 8(b)). The second column gives the number of first neighbor bonds (J1) connecting the chosen point to the points on the neighboring planes of the given separation, an so on. Fig. 4(a). Along the same line of reasoning presented in Subsection VII B. we can construct a family of states of ferromagnetic (100) planes, see Fig. 8. Representing a plane by a single effective spin s i of unit length, where now the index "i" enumerates the consecutive planes one can derive an effective onedimensional model: where the effective exchange can be inferred either from Table V or Fig. 8, and L (100) is the number of (100) planes in the crystal. Substituting the actual values J 2 = 2J 1 and J 3 = −J 1 /2, we see that the first two terms disappear, so the planes disentangle, and the last term gives +6J 1 for the correct ground state energy per site of the original model. The resulting ground state is of the form F 1 F 2 F 3 F 4 . . . , where F i denotes the independent ferromagnetic planes. This independence of the planes can be further rationalized by noting that both the first and second neighbor planes are connected by the covering rectangles, and the rectangle rule is satisfied bondwise on every ferromagnetic plane: see the rectangles in Fig. 8, and remember that the signs are distributed according to Fig. 5(c) and the planes are ferromagnetic. Such a state of ferromagnetically aligned independent (100) planes can be Fourier decomposed as and depending on the complexity of the real space pattern, any symmetric set of points on the ∆(q, 0, 0) line can be present in the expansion. Of course we could have chosen any of the symmetry related 100 directions. Stacking antiferromagnetic planes is more restrictive: a rect 2 can connect neighboring planes by J 1 bonds, in this case the J 2 bonds lie in-plane (and connect parallel spins, automatically satisfying the J 2 in-plane bonds), see Fig. 8(a). Another possibility for a rect 2 to connect second neighbor planes by J 2 bonds, and the J 1 bonds lie in-plane (and connect antiparallel spins), see Fig. 8(b). Second neighbor AFM planes are locked: they have to have the same AFM pattern to satisfy the rectangle rule. This leaves us with the two possibilities of stacking: an alternating set of two independent AFM planes A 1 A 2 A 1 A 2 . . . , or we can put independent FM planes between the AFM ones: In the finite cluster Ising solutions (see Appendix C for details) we have found the {100} stacking of independent FM planes: F 1 F 2 F 3 F 4 . . . and the FM stacking with intercalating AFM planes: AF 1 AF 2 A . . . . The alternating AFM stacking is missing here: for Ising spins it is an alternating FM stacking F 1 F 2 F 1 F 2 . . . viewed from a perpendicular direction.
We have performed numerical simulations for planar spins (XY -model): besides the aforementioned planar structures we found (seemingly disordered) ferromagnetic chains along the 100 directions, corresponding to a Fourier pattern of points on two perpendicular lines of M 1 ∆ in q-space. Thus we try to construct a family of states consisting of ferromagnetic (100) directed chains. These chains sit on a square lattice of primitive vectors (0, 1/2, 0) and (0, 0, 1/2), see Fig. 9(a). We represent a chain by a single effective spin s i of unit length, where now the indices i refer to points of the square lattice and δ to the neighbors of the lattice, and we map the system to the effective two dimensional model: For the M 1 ∆ manifold the effective interactions are: K 1 = 2J 1 , K 2 = 0, K 3 = 2J 1 , and K 4 = −J 1 , see Fig. 9(a) for a picture of the generated interactions. Strong, ferromagnetic J 2 = 2J 1 < 0 bonds connect along the chains and N (100) is the number of (100) chains in the crystal. This model has a codimension-one M GS : in the BZ of the square lattice the minima reside on the cross connecting the BZ center to the midpoints of the zone boundary together with the zone corner, see Fig. 9(c), note that this manifold is nothing but the intersection of M 1 ∆ with the (100) q-plane passing through the origin. The energy per site is 6J 1 = −6 |J 1 | (4J 1 comes from the interactions and J 2 = 2J 1 from the additional constant). This Hamiltonian can also be written as a sum of squares on signed rectangles inherited from the rect 2 -s projected to the (100) plane, see Fig. 9(a). This is consistent with the numerical findings of the planar spins: all the configurations found obeyed this projected rectangle rule, but seemed otherwise disordered.
• Stacking of independent ferromagnetic layers separated by the same antiferromagnetic layers on the {100} planes in an AF 1 AF 2 AF 3 . . . style. This type of ordering is realized in all the Ising, XY and O(3) models.
• Stacking of an alternating set of two independent (100) AFM planes A 1 A 2 A 1 A 2 . . . , reaalized in the XY and O(3) models.
• Interacting ferromagnetic chains in the 100 directions, these are absent in the Ising models. On the phase boundary line separating the W 1, 1 2 , 0 and X(1, 0, 0) phases (the red line in Fig. 2(b)) the ordering vectors of the possible ground states form the one dimensional M 1 Z manifold, see Fig. 3(b). Note that this manifold connects the points X(1, 0, 0) and W 1, 1 2 , 0 in the BZ. On this phase boundary given by J 1 > 0, J 3 = J 2 /4, −2 ≤ J 2 ≤ 0 the Hamiltonian reads: Since we have two free parameters (J 1 and J 2 ) expressing this Hamiltonian as the sum of complete squares on finite motifs is a little bit tricky. Here we use the elementary edge sharing tetrahedra of the fcc lattice and signed rectangles: Fig. 5(a) and (b). Two tetrahedra and two rectangles share a nearest neighbor bond, and four rectangles share a second neighbor bond, and each third neighbor bond is covered once by a rectangle, see Table III. "Signed" means that in the complete squares on these rectangles we associate a minus sign to the spins sitting on the sites denoted by white dots in Fig. 5(b), and plus signs to the black dots. The Hamiltonian becomes: Since the prefactors are all positive the Hamiltonian is minimized if and only if the spins sum up to zero on every tetrahedron and on every rectangle (with the appropriate signs), and the additional constant 2(J 2 − J 1 )N gives the ground state energy. The spin sum on rect 1 can be built by subtracting the spin sums of two edge-sharing tetrahedra, so every configuration that satisfies the tetrahedron rule automatically satisfies rectangle rule. At J 2 = J 3 = 0 we do not need the rectangles, and only the tetrahedron rule survives [46]. As an example [33] we can make ground states of (100) independent antiferromagnetically ordered planes in this endpoint: the spins form a checkerboard pattern on the planes of an A 1 A 2 A 3 . . . stacking style, where A i refers to the ith antiferromagnetic plane. This construction extends without modification to the whole X(1, 0, 0) − W 1, 1 2 , 0 boundary. This configuration is indeed a ground state, since both the tetrahedron and the rectangle rules are satisfied bondwise (for the motifs and sign distribution see Fig. 5(a) and (b), for the planes connected by the rectangles see Fig. 8(a) and (b)): spins on first neighbor bonds in a (100) plane are antiparallel and on second neighbor bonds are parallel. Just like in Section VII C the planes disentangle, and the in-plane contribution of interactions gives the correct ground state energy per spin as 2(J 2 − J 1 ).
A configuration of this stacking of AFM (100) planes can be Fourier expanded by combining ordering vectors from the (100) directed lines of the M 1 Z manifold (see Fig. 3(b)) and this nicely explains the shape of M 1 Z . Of course we also could have chosen the stacking direction of planes as (010) or (001). These states appear in the Ising solution, of course there are only two choices of AFM configurations in each plane.
In the simulations of planar spins we found (100)-directed FM chains, in a seemingly disordered distribution. Applying the effective two dimensional model for the chains forming a (100) square lattice one gets Eq. (65) with effective interactions K 1 = 2J 1 , K 2 = J 1 + J 2 /2, K 3 = J 2 , and K 4 = J 2 /2 with J 1 > 0 and −2 < J 2 ≤ 0 (here we exclude the J 2 = −2 Γ − X − W triple point, and discuss it in Subsection VII E), see Fig. 9(a) for a picture of the generated interactions. This model has a codimension-one M GS : in the BZ of the square lattice the minima reside on the BZ boundary, see Fig. 9(c), but be careful: the zone center Γ(0, 0) is not part of the manifold. The Fourier transform of the effective exchange has a local but not global minimum at the BZ center, which gets lower and lower as we move along the X − W line towards the X − W − Γ point, and this minimum becomes degenerate with the M GS on the BZ boundary as we finally reach J 2 = −2J 1 . Note that this manifold is nothing but the intersection of M 1 Z with the (100) q-plane passing through the origin. The ground state energy per site is 2(J 2 −J 1 ) (J 2 −2J 1 comes from the interactions and J 2 from the additional constant in Eq. (65)). This Hamiltonian can also be written as a sum of squares on signed rectangles inherited from the tetrahedra and rect 1 -s projected to the (100) plane, see Fig. 9(a) (the projected tetrahedron rule prohibits Γ(0, 0) being a global minimum). This is consistent with the numerical findings of the planar spins: all the configurations found obeyed this projected rectangle and tetrahedron rule, but seemed otherwise disordered.
To summarize we have found the following candidate ground states for the X(1, 0, 0) − W 1, 1 2 , 0 line: • Stacking of independent antiferromagnetic {100} planes in the style A 1 A 2 A 3 A 4 . . . . This type of ordering is realized in all the Ising, XY and O(3) models.
• Interacting ferromagnetic chains in the 100 directions, these are absent in the Ising models.
At the endpoint J 2 = −2J 1 and J 3 = −J 1 /2 the tetrahedron rule vanishes in Eq. (68), and the only constraint is that spins on rectangles have to satisfy the equation S 1 +S 2 −S 3 − S 4 = 0, this less restrictive condition offers other possibilities (e.g. the appearance of a net magnetization), we devote the next subsection to its analysis.. The Γ(0, 0, 0) − X(1, 0, 0) − W 1, 1 2 , 0 triple point bears striking resemblance to the triple point Γ(0, 0, 0) − ∆(q, 0, 0) − X(1, 0, 0) presented in Subsection VII C, and has a much richer structure than the rest of the X(1, 0, 0) − W 1, 1 2 , 0 line. Here the parameters are J 2 = −2J 1 , J 3 = −J 1 /2, J 1 > 0, and the M GS is M 1 Z ∪ Γ, see Table II. See Fig. 3(b) for the degenerate manifold, and Fig. 2(b) for the point in the phase diagram: the red dot where the X − W boundary line hits the Γ(0, 0, 0) phase. We can cover the lattice by signed rectangles (here the tetrahedron rule does not apply), where the signs are distributed according to Fig. 5(b) now: since J 1 /4 > 0 this Hamiltonian is minimized when S 1 +S 2 − S 3 − S 4 = 0 on every rectangle, and the ground state energy per site is −6J 1 . This rule is compatible with ferromagnetism. We can map the models Γ(0, 0, 0)−X(1, 0, 0)−W 1, 1 2 , 0 and Γ(0, 0, 0)−∆(q, 0, 0)−X(1, 0, 0) to each other by changing the sign of J 1 but keeping the other two interactions intact. In the following we exploit the relationship between the two models, and for the details we refer to VII C.
Stacking of FM and AFM (100) planes works in complete analogy with Subsection VII C, we only need to interchange the words antiferromagnetic and ferromagnetic, and instead of rect 2 we have to use rect 1 . The possible orderings constructed by stacking FM/AFM (100) planes (confirmed by the solution of the Ising model and numerical results on the XY -model) are of the form of an alternate stacking of two independent FM planes: F 1 F 2 F 1 F 2 . . . , of independent AFM layers: A 1 A 2 A 3 A 4 . . . , and of independent AFM layers separated by FM planes of fixed magnetization direction: In the numerical solution of the planar model we find (100) chains again, and we can apply the effective two dimensional model of Eq. (65) on the square lattice, but now with parameters K 1 = 2J 1 , K 2 = 0, K 3 = −2J 1 , and K 4 = −J 1 . Note, that strong ferromagnetic J 2 = −2J 1 < 0 bonds connect along the chains again. This model also has a codimension one M GS : in the BZ of the square lattice the minima reside on the BZ boundary together with the Γ(0, 0) point, see Fig. 9(b). Note that this manifold is nothing but the intersection of M 1 Z with the (100) q-plane passing through the origin extended with the Γ(0, 0, 0) point. The energy per site is −6J 1 . This Hamiltonian can also be written as a sum of squares on signed rectangles inherited from the rect 1 -s, projected to the (100) plane, see Fig. 9(a). This is consistent with the numerical findings of the planar spins: all the configurations found obeyed this projected rectangle rule, but seemed otherwise disordered.
All the ground states found above can be mapped to the ground states of the Γ(0, 0, 0) − ∆(q, 0, 0) − X(1, 0, 0) model by choosing chains along one of the 100 directions and changing the signs of all the spins on every second chain in a checkerboard pattern (i.e. we flip the spins on all the white (100) chains in Fig. 9(a)).
• Interacting ferromagnetic chains in the 100 directions, these are absent in the Ising models.

VIII. CONCLUSIONS
We presented a detailed study of the ground state phase diagram of the classical isotropic J 1 -J 2 -J 3 Heisenberg model on the face-centered-cubic lattice. We found and analyzed in detail -in real and Fourier space-the commensurate Type I, II and III structures (but not the Type IV one), where the multiple-Q orderings allow for noncoplanar and even chiral structures. Besides the commensurate orders, the introduction of the third neighbor coupling resulted in incommensurate spin spiral orders with propagation vectors along high symmetry axes of the crystal.
We also found ground state manifolds in q-space of dimension one and two with subextensive degeneracies at phase boundaries. In all cases, we could express the Hamiltonian as a positive definite sum of complete squares on finite motifs covering the lattice. This reformulation provided us with hints to explicitly construct large classes of nontrivial, aperiodic ground states in real space, consisting of randomly stacked ordered planes and frustrated ferromagnetically ordered chains in special crystallographic directions. We described relations of the real space patterns to the q-space picture.
We thoroughly analyzed the model for Ising spins on finite clusters in the phases with extended manifolds, and determined the number and type of possible configurations in real and Fourier space. We performed numerical simulations on XY models and confirmed the validity of our analytical results. Numerical studies on O(3) spins revealed even richer structures than considered here, these need further investigations.
It is interesting to compare the commensurate orders found here for the fcc lattice to the construction of regular magnetic orders in Ref. [53]. Besides the trivial ferromagnetic order only the Type I antiferromagnet can be regular, and the latter only if we choose |ξ| = |η| = |ζ| = 1/ √ 3 in Eq. (15), i.e. for equally weighted Bragg peaks: this is the three dimensional analogue of the tetrahedral state presented in Ref. [53].
Our work sets the stage for future studies aimed at investigating the finite-temperature classical phase diagram of the J 1 -J 2 -J 3 model including an investigation of its critical phenomenon which has, till date, largely focused only on the nearest neighbor model. In particular, the triple points and phase boundaries which are host to a subextensively degenerate manifold of ground states would provide for a promising route towards potentially realizing classical as well as quantum (for S = 1/2, 1) spin liquids on a three-dimensional lattice, in the scenario when order-by-disorder fails to lift the degeneracy as is known to occur for the pyrochlore [40][41][42] and hyper-hyperkagome lattices [61]. The triple points occurring in the J 1 -J 2 -J 3 Heisenberg model on the simple cubic and body-centered-cubic lattices are known to give way to a quantum paramagnetic phase for S = 1/2 [62][63][64][65]. Given the fact that three of the degenerate manifolds involve a ferromagnetic phase implies that in the scenario that long-range dipolar magnetic orders are absent, multipole orders such as quadrupolar [66][67][68][69][70][71][72][73][74][75][76], and octupolar [37] orders could be stabilized in both classical and quantum models. The role of disorder in stabilizing noncollinear phases will also be an interesting endeavor for future studies [47,77].
is always an integer. The corresponding reciprocal lattice vectors are: −1, 1) , b 3 = 2π (−1, 1, 1) , (A2) we will refer to any point in reciprocal space by its q-triplet, e.g. q = (q x , q y , q z ) = (2π, 2π, −2π) = b 1 . Special points and lines in the Brillouin zone (BZ) have more or less commonly used labels, we will refer to them either by their labels, or the labels with their Cartesian coordinates in parenthesis in units of 2π, e.g. one of the BZ corners of the fcc lattice can be referred to as W , W 1, 1 2 , 0 , or (2π, π, 0). The Fourier transform of the exchange interaction for the fcc lattice with first, second and third neighbor interactions presented in Eq. (1) is defined by and it reads: J(q) = 2J 1 cos q x 2 cos q y 2 + cos q x 2 cos q z 2 + cos q y 2 cos q z 2 + J 2 (cos q x + cos q y + cos q z ) + 4J 3 cos q x cos q y 2 cos q z 2 + cos q x 2 cos q y cos q z 2 + cos q x 2 cos q y 2 cos q z . (A4) Appendix B: To see the way the degeneracy in q-space of the manifolds manifests itself in real space, we have considered Ising spins on finite clusters. The degeneracy of the Ising spins is also the degeneracy for collinear O(3) spin configurations, which is just half of the degeneracy of the Ising manifold if we factor out the trivial O(3) global rotation.
First, we generate the linear set of equations defining the manifold (e.g. the tetrahedron rule) on a finite cluster with periodic boundary conditions. The finite clusters are defined by the superlattice vectors g 1 , g 2 , and g 3 , such that S i+g1 = S i and so on. The clusters of different shape are listed in Table VII, together with the number of planes in different directions. Next, since the set of linear equation is homogeneous, we search for the null space (or kernel) they define. The dimension of the null space D NS depends on the type of the manifold and on the shape of the cluster, and the N − D NS spins in the cluster can be expressed as linear combinations of D NS linearly independent spins. This would suggest that the number of Ising configurations is 2 DNS -however not all of the solutions satisfy the spin length constraint. In order not to miss a configuration, we generate by computer all the 2 DNS linear combinations, and keep only those that give Ising spins on every site. We have collected the numerical results in Table IX and discuss the different manifolds below. The findings for the 1D manifolds are summarized in Table VIII. 1. The M 2 2D manifold with the octahedron rule We discussed M 2 in Sec. VII A. The spins shall obey the octahedron rule (Eq. (50). The numerical results are summarized in the last two columns of Table IX. Seemingly, the dimension of the null space depends randomly on the size of the cluster. Connecting the real space picture to the reciprocal space reveals that the dimension of the null space is equal to the number of discrete q points which satisfy Eq. (47), i.e. which lie on the two-dimensional manifold shown in Fig. 3(c). The number of Ising spin configurations is typically much larger then in 1D manifolds.
The "signed square" constraint (53) provides 3 equations per site, but the number of linearly independent equations grows linearly with the system size, more precisely linearly with the number of the {111} planes, as seen in Table IX and  summarized in Table VIII. This is in perfect agreement with the results of Sec. VII B, that this manifold consists of independent up-or downpointing ferromagnetic triangular {111} planes. Since there are four 111 directions, this leads to 4 × 2 L − 6 configurations in this manifold, the constant 6 is compensating for the multiple counting of the 8 periodic single-Q L 1 2 , 1 2 , 1 2 states. For an even number of planes 8 additional states appear which do not have the layer structure. They are quadruple-Q L 1 2 , 1 2 , 1 2 states, with 4 amplitudes equal in absolute value. The up and down spins form two interpenetrating pyrochlore lattices, the unit cell consists of 8 sites (e.g. the A = B = C = D = 1 andĀ =B =C =D = −1 in Fig. 4(b) represents one of the 8 states).
Let us now see what do Eqs. (56) tell us for Ising spins. In this case a, b, c, m ∈ R and Eq. (56a) becomes e 2 +a 2 = 1 and ea = 0. These two equations are satisfied by either e 2 = 1 and a = 0, or e = 0 and a 2 = 1. The first solution implies b = 0 and c = 0, and the spins in the octahedron are all identical to e = ±1, resulting in a ferromagnetic configuration. When e = 0, there are eight solutions, a = ±1, b = ±1, and c = ±1 corresponding to the choice of the three signs, and the solutions describe structures where on the opposite {111} faces of the octahedra we have opposite spins. We can use these and cos (qΣ/2) = . The third column gives the equations of the phase boundaries, and the last column gives the order of the transition. Compare this table with the phase diagram given in Fig. 2. For pictures of wave-vectors in the Brillouin zone see Fig. 3.
The equation for the phase boundary ε (Q W ) = ε (Q Σ ) is:   The periodic states are the 6 single-X states and the two fully polarized Γ(0, 0, 0) states.
2. AF 1 AF 2 A . . . F L/2 like configuration with alternating ferro-and antiferromagnetic layers, where ferromagnetic layers are independent, but the antiferromagnetic layers are locked with respect to each other. The number of states is 3 × 4 × 2 L/2 − 16. Among these configurations are the 8 quadruple-Q states made from the three X α -s and the Γ(0, 0, 0) Q vector.
Altogether there are 3 × 4 L + 12 × 2 L − 20 configurations. For cluster that are not compatible with the 4-site cubic unit cell, only the F 1 F 2 F 3 F 4 . . . F L states are allowed, so the degeneracy is 3 × 2 L − 4.

The M 1 Z manifold with the tetrahedron rule
The states of the Heisenberg model with only nearest neighbor interactions belong to this manifold. The ground states satisfy the tetrahedron rule -the sum of the spins on every elementary tetrahedron is zero. This constraint gives two linear equations per site, but these equations are not linearly independent. Since the tetrahedra are edge sharing, the num-ber of linearly independent equations is greatly reduced, and scales with the linear size of the cluster, more precisely with the number of {100} parallel planes. The dimension of the null space (shown in Table VIII) is 3L (100) − 3 when there are L (100) parallel planes, the factor of 3 comes from the 3 equivalent directions of the parallel planes in clusters respecting the full cubic symmetry.
We may now count the degeneracy of the manifold assuming Ising spins. Choosing a direction, say [100] and the corresponding set of parallel planes, each of (100) planes is antiferromagnetic with a Z 2 degree of freedom (we can exchange spins on the two sublattices), the total number of Ising configurations is 2 L . Since we have three possible directions, the total number of Ising configurations is 3 × 2 L (100) − 6, 6 compensates for the multiple counting of the periodic configurations consisting of antiferromagnetic planes in two directions and ferromagnetic planes in the third direction -these are the single-X states.
We can extend the covering with tetrahedra by signed rectangles, as discussed in Sec. VII D, to allow for second and third neighbor exchanges. Even though the number of linear equations increases by 6 per site, they do not lower the dimension of the null space, neither do they change the number of Ising configurations.

The M 1 Z ∪ Γ manifold
This is the endpoint of the line M 1 Z in the J 1 − J 2 − J 3 parameter space, where the tetrahedron rule is lost and only the rectangles remain, see Sec. VII E for details. The Γ point appears as an allowed Q vector. The number of equations is 6 per site and the dimension of the null space has increased by one compared to the pure M 1 Z manifold. The allowed Ising configurations are: 1. A 1 A 2 A 3 . . . A L like configurations: these are inherited from the M 1 Z manifold and their number is 3 × 2 L − 6. As we noted, the 6 periodic single-X states belong to this class.
2. F A 1 F A 2 F . . . F A L/2 configurations where ferromagnetic and antiferromagnetic planes alternate along a 100 direction when the number of planes is even. While the rectangles lock the spins in the ferromagnetic layers, the L/2 antiferromagnetic layers retain their Z 2 degree of freedom, and the tetrahedron rule is violated. The number of these states is 12×2 L/2 −16. The factor 12 comes from the 3 choices of the 100 directions, the polarisation of ferromagnetic planes (a factor of 2), and the choice of the first plane to be ferromagnetic or antiferromagnetic (another factor of 2). The periodic states are the 4-sublattice one-half magnetization plateau configurations (8 of them), that consist of the Γ(0, 0, 0) and X quadruple-Q structure, with coefficients equal in absolute value.
3. Pure ferromagnets: 2 Ising degeneracy. Altogether there are 3 × 4 L + 12 × 2 L − 20 configurations, just like for M 1 ∆ . For clusters that are not compatible with the 4 site cubic unit cell, the frustration only allows the 2 FM configurations. Table VIII. Discrete degeneracy of the one dimensional manifolds for Ising spins in finite size clusters respecting the full cubic symmetry of the fcc lattice. The degeneracy depends on the number of the (100) planes (or (111) planes in the case of the M 1 Λ manifold). Depending on the even/odd number of the planes, the frustration reduces the degeneracy.

Manifold
Eq/N Planes Dim. of Number of Ising type No. Parity null space configurations  Table IX. Numerical enumeration of the Ising configurations. The first two columns are the number of (100) and (111) planes in the cluster, the next three columns show the geometry of the clusters with periodic boundary conditions, the number of sites in the cluster is shown in sixth row. The following rows list the dimension of the null space DNS (i.e. the number of linearly independent equations) and the number of Ising configurations.
No. planes