Quantum and Classical Phases of the Pyrochlore Heisenberg Model with Competing Interactions

We investigate the quantum Heisenberg model on the pyrochlore lattice for a generic spin $S$ in the presence of nearest-neighbor $J_{1}$ and second-nearest-neighbor $J_{2}$ exchange interactions. By employing the pseudofermion functional renormalization group method, we find, for $S=1/2$ and $S=1$, an extended quantum-spin-liquid phase centered around $J_{2}=0$, which is shown to be robust against the introduction of breathing anisotropy. The effects of temperature, quantum fluctuations, breathing anisotropies, and a $J_{2}$ coupling on the nature of the scattering profile, and the pinch points, in particular, are studied. For the magnetic phases of the $J_{1}$-$J_{2}$ model, quantum fluctuations are shown to renormalize phase boundaries compared to the classical model and to modify the ordering wave vectors of spiral magnetic states, while no new magnetic orders are stabilized.


I. INTRODUCTION
The classical nearest-neighbor Heisenberg antiferromagnet on the pyrochlore lattice stands as an epitome of geometric frustration in three dimensions as shown by its failure to develop magnetic long-range order down to absolute zero temperature, realizing what has been dubbed a "cooperative paramagnet" [1]. This failure is a consequence of the extensive classical ground-state degeneracy [1][2][3][4] which proves severe enough to prevent a thermal "order-by-disorder" mechanism [5][6][7] from selecting a unique ground-state ordering pattern [3,4,8,9]. In contrast to thermal fluctuations, the impact of quantum fluctuations remains much less understood and constitutes a critically outstanding problem. In the regime of large spin S, using an effective Hamiltonian approach [10], it is known that at harmonic order in 1/S, the extensive classical ground-state degeneracy exp[O(L 3 )] (L is the linear dimension of the system) is partly lifted, yielding a subset of collinear states with a massive, albeit subextensive, degeneracy exp[O(L)] [11][12][13][14]. It turns out that the consideration of higher-order terms in a 1/S expansion also fails to select a unique ground state [15]. Indeed, while quartic corrections in boson operators do break the degeneracy of * yiqbal@physics.iitm.ac.in the harmonic ground states, there still remains a family of (almost) degenerate (exp[O(L)]) states [16]. Thus, the fate of the semiclassical (1/S) approach remains unsettled due to weak selection effects at the anharmonic level. In the opposite extreme quantum limit of small S, there is reasonably strong evidence for a quantum paramagnetic ground state. Investigations of the S = 1/2 antiferromagnet claim for either a valence-bond crystal [17][18][19][20][21][22][23][24] or a quantum-spin-liquid [25][26][27][28][29][30][31] ground state. We note that a J 1 -J 2 -J 3 S = 1/2 model derived from a strong-coupling expansion of a one-band half-filled Hubbard model on the pyrochlore lattice has been proposed to host a quantum spin liquid [32,33]. In the much-less-investigated case of S = 1 [19,[34][35][36], there have been suggestions of a ground state with tetrahedral symmetry breaking [37].
The "cooperative paramagnet" ground state of the classical nearest-neighbor Heisenberg antiferromagnet is known to be extremely fragile, in that magnetic longrange order is induced upon the inclusion of various perturbations, such as further neighbor Heisenberg interactions [2,[38][39][40][41], dipole interactions [42], Dzyaloshinsky-Moriya anisotropy [43,44], single-ion anisotropy [45,46], lattice distortions [47][48][49][50][51][52], and bond disorder [53][54][55]. In particular, further neighbor Heisenberg interactions are found to stabilize a plethora of intricate classical magnetic orders [56,57]. However, in the low-spin-S regime, where the strong possibility of a quantum para- [111] FIG. 1. The nearest-neighbor (J1) and next-nearest-neighbor (J2) bonds in the pyrochlore lattice. magnetic ground state for the nearest-neighbor quantum Heisenberg antiferromagnet exists, the impact of the above-mentioned perturbations on the paramagnet remains largely unexplored. This topic is of high significance and importance when considering the behavior of real materials. In this paper, we carry out a broad investigation of the J 1 -J 2 Heisenberg model for a generic spin S on the pyrochlore lattice: whereŜ i is a quantum spin-S operator at a pyrochlore lattice site i. The indices i, j and i, j denote sums over nearest-neighbor and second-nearest-neighbor pairs of sites, respectively [see Fig. 1]. The investigation of the low-temperature properties of this Hamiltonian in the small-S regime is notoriously difficult. This is a methodological challenge for which numerically exact and unbiased methods are not yet available. Indeed, traditional quantum many-body numerical methods such as density-matrix renormalization group and tensor network approaches [58,59], while successful in one and two dimensions, become unfeasible in three dimensions due to entanglement scaling and system-size limitations. Quantum Monte Carlo methods [60,61], while able to reach sufficiently large system sizes, are, in principle, restricted to unfrustrated systems, while variational Monte Carlo approaches [62,63], which are shown to be extremely successful in two dimensions [64][65][66], require very large correlation volumes to extract reliable estimates in the thermodynamic limit. Finally, the bold diagrammatic Monte Carlo method can reach down only to moderately low temperatures [31]. Thus, one is essentially left with only mean-field approaches based on Schwinger bosons [67], semiclassical analysis based on spin waves, or linkedcluster expansion methods [68], which capture magnetic order accurately but are unsuitable for studying paramagnetic behavior deep in the collective paramagnetic (spin-liquid) regime. In this respect, the pseudofermion functional renormalization group (PFFRG) framework has an important feature in the form of a built-in balance towards the treatment of ordering and disordering tendencies for three-dimensional frustrated magnets [69].
By employing PFFRG for the spin-S J 1 -J 2 Heisenberg model, we find for S = 1/2 an extended quantum-spinliquid regime centered around J 2 = 0, with an extent of −0.25 (3) J 2 /J 1 0.22 (3) while, for S = 1, its span is reduced by approximately a factor of 2, −0.11(2) J 2 /J 1 0.09 (2). For S = 1/2 and S = 1, the spin susceptibility profile of the nearest-neighbor antiferromagnet in the [hhl] plane features a bow-tie pattern, characteristic of the well-known Coulomb spin-liquid phase [70]. The bow ties are found to be robust up to temperatures T /J 1 ∼ 1. However, the inclusion of even a small J 2 coupling is shown to shift the spectral weight away from the pinch points, causing the bow ties to rapidly disappear upon cooling, similar to the findings for the corresponding classical model [71]. In the opposite limit of large S, quantum fluctuations lift the extensive degeneracy of the classical ground-state manifold either only partially to a subextensive one or completely (which would then potentially induce long-range magnetic ordering). The J 1 -J 2 parameter space is known to host seven different classical magnetic orders [56], which we also find in the S = 1/2 model. Moreover, we show that quantum fluctuations do not stabilize any new phases, such as long-range dipolar or quadrupolar magnetic orders, and valence-bondcrystal states.
The paper is organized as follows: In Sec. II, we describe the PFFRG method (Sec. II A) employed for the quantum treatment of the model, starting with a description of its formalism (Sec. II A 1) followed by details of its numerical implementation in Sec. II A 2. In Secs. II B and II C, we discuss schemes used to obtain the ground state of classical spin models, namely, the Luttinger-Tisza method [Sec. II B] and the iterative minimization of the energy [Sec. II C] (the reader interested mainly in the results can directly jump to Secs. III and IV). Employing these methods, we begin with a treatment of the ground-state and low-energy physics of the nearest-neighbor Heisenberg antiferromagnet in Sec. III, starting first with a classical analysis [Sec. III A 1] of the isotropic and breathing lattices and then moving on to the quantum treatment of the S = 1/2 [Sec. III B] and S = 1 [Sec. III C] models for both isotropic and breathing lattices. Finally, the section ends by addressing the problem of the ground state of the large-S quantum Heisenberg antiferromagnet [Sec. III D]. Next, in Sec. IV, we deal with the J 1 -J 2 Heisenberg model, by first revisiting the classical phase diagram [Sec. IV A], and subsequently present the results for the quantum model in Sec. IV B. We also discuss the impacts of quantum fluctuations on the nature of phases and phase boundaries. We end the paper with a summary of the results in Sec. V, followed by an outlook and discussion of future directions in Sec. VI.

II. METHODS
A. Pseudofermion functional renormalization group method

Formalism
The key idea of the PFFRG method [72] is to express the spin-1/2 operators in terms of pseudofermions [73], where σ µ αβ are Pauli matrices (µ ∈ {x, y, z}) andf iα (f † iα ) denote spin-α fermionic annihilation (creation) operators. For the implementation for spin systems with local S > 1/2 spins, we adopt the approach of Ref. [74], where multiple copies of spin-1/2 degrees of freedom are introduced at each lattice site; i.e., the local spin operators are replaced byŜ while the couplings J ij remain independent of the fermion "flavor" κ. If all individualŜ iκ "spins" ( κ ∈ {1, . . . , M }) align ferromagnetically (see below for details), they realize the largest possible magnitude S = M/2 on each site, thus implementing the desired effective magnetic moment. In terms of pseudofermions, the substitution in Eq. (3) amounts to equipping the fermion operators with an additional index κ: Pseudofermionic representations for spin operators generally require some caution, since they introduce additional spurious states with zero (Q i ≡ f † i↑ f i↑ +f † i↓ f i↓ = 0) or two (Q i = 2) fermions at a site i. Such states carry no spin (S = 0), and the physical spin-1/2 degrees of freedom are realized in the singly occupied subspace with Q i = 1. The pseudofermionic approach is guaranteed to be faithful only if the contribution from the S = 0 states is negated. For a proper implementation of spins S > 1/2, one additionally needs to ensure that the spin flavors κ combine to the largest local moment S = M/2 while smaller spins with S = M/2 − 1, . . . are eliminated from the Hilbert space. A convenient approach that simultaneously fulfills both constraints is to add an on-site local level repulsion term A( M κ=1Ŝ iκ ) 2 to the Hamiltonian. For negative A, this term reduces the energies of all levels with finite magnetic moments, where the largest reduction occurs in the sector with the highest spin. An |A| chosen sufficiently large guarantees that the low-energy subspace of the Hamiltonian is the one without any nonoccupied or doubly occupied states for each κ. Furthermore, the M spin-1/2 copies combine into an effective spin S = M/2. We emphasize, however, that, for the ground states of generic Heisenberg spin models (such as the pyrochlore systems studied here), a vanishing level repulsion term A = 0 turns out to be sufficient to fulfill both pseudoparticle constraints. This simplification is because, for two-body spin interactions, the energy naturally scales with the spin length squared such that the largest local moment is energetically favored even for A = 0 (note, however, that counterexamples can be constructed [75]).
Rewriting the spin Hamiltonian in terms of Eq. (4), the resulting fermionic model is treated within the standard FRG framework for interacting fermion systems [76][77][78]. A somewhat unusual situation arises here: the system is purely quartic in the fermions without any quadratic kinetic terms that could be used as a noninteracting starting point in a perturbative expansion. Within FRG, this situation is addressed by summing up infinite-order diagrammatic contributions in different interaction channels as well as accounting for vertex corrections between them. Particularly, as explained in more detail below, the summation is such that, in the large-S and the large-N limits, where N generalizes the spin symmetry group from SU(2) to SU(N ), the leading diagrammatic contributions in 1/S and 1/N are both treated exactly [79]. As a consequence, classical magnetically ordered states (typically favored at large S) and nonmagnetic spin liquids or dimerized states (as obtained at large N ) [80] may both be described within the same methodological framework.
Because of the absence of fermion kinetic hopping terms, the bare fermionic propagator is strictly local and takes the simple spin-independent form where iω denotes a frequency on the imaginary Matsubara axis. Within the standard PFFRG scheme [72], this propagator is dressed with an infrared steplike regulator function: which interpolates between the limits Λ → ∞ (where the fermionic propagation is completely suppressed) and the original cutoff-free theory at Λ = 0. This modification generates a Λ dependence of all one-particle irreducible m-particle vertex functions as described by the FRG flow equations. For the self-energy Σ Λ (iω) and the two-particle vertex Γ Λ (1 , 2 ; 1, 2) (the label "X" stands for site, frequency, and spin variables, respectively, i.e., X ≡ {i, iω, α}). A diagrammatic version of these equations is illustrated in Fig. 2, where the arrows denote dressed and Λ-dependent propagators FIG. 2. Diagrammatic representation of the PFFRG equations for (a) the self-energy Σ Λ (iω) (gray disk) and (b) the twoparticle vertex Γ Λ (1 , 2 ; 1, 2) (gray squares). Arrows denote the fully dressed propagator G Λ (iω), and slashed arrows denote the single-scale propagator S Λ (iω). The gray hexagon in (b) is the three-particle vertex. Note that the right-hand side of (b) contains additional terms where the slashes in the first to fifth terms appear in the respective other propagator. For a spin-S generalization, the first term on the right-hand side of (a) and the second term in (b) are multiplied with a factor of 2S.
and slashed lines denote the single-scale propagator Because of the locality of fermion propagators, the twoparticle vertex Γ Λ (1 , 2 ; 1, 2) effectively depends on two site indices only, i.e., Γ Λ (1 , 2 ; 1, 2) ∼ δ i1i 1 δ i2i 2 . As illustrated in Fig. 2, this restriction allows one to connect incoming and outgoing arrows of Γ Λ (1 , 2 ; 1, 2) in a way that on-site variables remain constant along fermion lines. The FRG equations in Fig. 2 show a systematic interplay between the RG flows of different vertex functions where the Λ derivative of each m-particle vertex couples to all m -particle vertices with m m + 1. To reduce this infinite hierarchy of intertwined equations to a finite and numerically solvable set, we neglect the threeparticle vertex in Fig. 2(b) albeit not in entirety, as certain three-loop terms obtained from the Katanin truncation scheme are included and which amount to self-energy corrections [81], as described below; however, all higher vertices are completely discarded. However, this approximation effectively amounts to discarding three-body spin correlations such that the description of spin phases with chiral order parameters Ŝ i ·(Ŝ j ×Ŝ k ) is not possible [82]. Still, parts of the three-particle vertex can be included by applying the so-called Katanin truncation [81], which replaces the single scale propagator by the full Λ derivative of the dressed propagator While the additional Katanin terms formally have the structure of the three-particle term [the last term in Fig. 2(b)], they should rather be understood as selfenergy corrections [81]. Indeed, the Katanin truncation ensures full self-consistency at the two-particle level in the sense that the self-energy is completely fed back into the flow of Γ Λ . This feedback is particularly important for the description of strongly fluctuating spins which requires two-particle vertex renormalizations beyond the bare ladder summations. Together with the initial conditions defined in the limit Λ → ∞ (where the self-energy vanishes and the two-particle vertex reduces to the bare couplings J ij ), the closed set of differential equations is now amenable to numerical treatment.
According to standard diagrammatic Feynman rules, the implementation of spins S > 1/2 via the local replication of S = 1/2 degrees of freedom [see Eq. (3)] introduces additional sums over flavor indices κ for all closed fermion loops in the PFFRG equations. Since the bare couplings J ij are independent of κ, this summation simply leads to an extra factor M = 2S in the Hartree contribution for the self-energy [the first term on the right-hand side of Fig. 2(a)] and in the RPA contribution for the two-particle vertex [the second term on the right-hand side of Fig. 2(b)]. Increasing S consequently strengthens the RPA term with respect to the other terms, indicating that these diagrams are responsible for the formation of classical magnetic long-range order. Indeed, one can show that, in the absence of finite-temperature divergencies of subleading 1/S diagrams, the bare RPA channel (which accounts for only leading 1/S diagrams) correctly reproduces the classical limit S → ∞ where the PFFRG becomes identical to the Luttinger-Tisza method [74]. We mention that a correct treatment of the classical nearest-neighbor Heisenberg antiferromagnet indeed requires accounting for the effects of subleading 1/S diagrams as discussed in Appendix A. In a similar way, the PFFRG method can be generalized to treat SU(N ) spins with N > 2. In such a scheme, the ladder channels [first and fifth terms on the right-hand side of Fig. 2(b)] contribute with an additional factor of approximately N , indicating that these terms describe nonmagnetic spin liquids or dimerized states. In analogy to a large S generalization, they become exact in the limit N → ∞. This built-in balance between large-S and large-N terms represents the key property of the PFFRG that allows one to study magnetic order and disorder tendencies on fair footing. The PFFRG was initially developed in two dimensions [72]; however, subsequent refinements have made it capable of handling a wide spectrum of frustrated magnetic Hamiltonians for multilayer systems and in three dimensions [69,[83][84][85][86][87][88][89][90][91][92][93][94][95][96][97][98][99][100][101][102].
2. Numerical solution of PFFRG flow equations and probing the nature of the ground state To solve the PFFRG equations numerically, we approximate the spatial dependence of Γ Λ (1 , 2 ; 1, 2) by discarding all vertices with a distance between sites i 1 and i 2 greater than some maximal value. In our calculations, we use a distance of approximately 11.5 nearest-neighbor lattice spacings, which corresponds to a total volume of 2315 correlated spins. Likewise, the continuous frequency arguments of the vertices are approximated by discrete meshes, for which we typically use a combination of linear and logarithmic grids consisting of 64 discrete frequency points.
By fusing the external legs (1, 1 ) and (2, 2 ) of the two-particle vertex Γ Λ (1 , 2 ; 1, 2), one can calculate the static spin-spin correlator where T τ (with τ being the imaginary time) is the imaginary time-ordering operator. Transforming χ zz ij into k space yields the wave-vectorresolved susceptibility χ(k): which is the central outcome of the PFFRG to probe the system's magnetic properties. Note that, since in the Heisenberg case the susceptibility is always isotropic, we omit the component indices xx/yy/zz in the susceptibility χ(k). Here, the first summation is carried out over the four sites of a given primitive unit cell, and the prefactor of 1/4 is the inverse of the total number of sites in the unit cell. This quantity has the periodicity of the extended Brillouin zone but not of the first Brillouin zone, and thus the susceptibilities are always presented in the former. Henceforth, all wave vectors k are expressed in units where the edge length of the pyrochlore cubic unit cell is one. The onset of long-range dipolar magnetic order is signaled by a divergence in the Λ flow of the susceptibility as observed in the thermodynamic limit. This divergence is a manifestation of the fact that the spin-spin correlations do not decay in the limit of long distances, which would ultimately cause the Fourier transform χ(k) to diverge. However, in the numerical calculations, we employ a frequency discretization and keep only a limited spatial range of the two-particle vertices; hence, the Fourier transform amounts to a finite site summation that no longer diverges. Thus, these divergences end up being regularized, manifesting themselves as kinks or cusps at some critical Λ c in the Λ evolution of the susceptibility (henceforth referred to as "breakdown of the RG flow") [see Appendix B for a discussion on the detection of magnetic instabilities in the RG flow].
The type of magnetic order is characterized by the wave vector at which the breakdown of the RG flow occurs. In 3D, the PFFRG ordering scales, i.e., Λ c , are directly related to the ordering temperatures T c via Tc J = 2πS(S+1) 3 Λc J [69]. The conversion factor 2πS(S + 1)/3 between the RG scale Λ and the temperature T can be obtained by comparing the limit of PFFRG where only the RPA diagrams contribute [74], i.e., a mean-field description, and the conventional spin mean-field theory which is formulated in terms of the temperature instead of Λ [103]. On the other hand, nonmagnetic (absence of dipolar magnetic order) ground states are signaled by a susceptibility flow that continues to evolve smoothly down to the (numerical) limit Λ → 0. Even in the absence of long-range dipolar magnetic order, the momentum profile of χ(k) at Λ 1 allows one to determine the dominant types of short-range spin correlations or to identify competing ordering tendencies.
In the absence of long-range dipolar magnetic order in the ground state, we can further probe for possible spin-nematic [3,4,104] and valence-bond-crystal orders [17][18][19][20][21][22][23][24] by computing the corresponding nematic and dimer response functions. Here, we are particularly interested in studying the tendency of the quantum paramagnet towards spontaneous breaking of either spin rotation symmetry, i.e., nematic order, or translational symmetry, i.e., dimer order. The onset of these orders is marked by the divergence of the corresponding orderparameter susceptibility, which is given by a four-spin correlator. For spin-nematic order, this correlator is the standard nematic correlation function [105,106] (with µ, ν = x, y, z denoting the three directions in spin space and i, j representing the lattice sites) is a symmetric traceless tensor. For dimer order, it is the singlet-singlet correla- In PF-FRG, such correlators are represented by the fermionic four-particle vertex, and, while the PFFRG formalism could, in principle, be straightforwardly extended to obtain the RG flow equation for the four -particle vertex, their numerical solution is, at present, not feasible due to limitations posed by computational complexity limitations and memory requirements. The fact that the four-particle vertex is a priori excluded from the RG equations implies that the RG flow of the spin susceptibility [Eq. (11)] is unaffected by the possible presence of competing nematic and dimer orders. Hence, we adopt a simple recipe within the PFFRG framework to calculate the nematic (dimer) response function η SN (η VBC ) which measures the propensity of the system to sup-port nematic (valence-bond-crystal) order. It amounts to adding a small perturbation to the bare Hamiltonian which enters the flow equations as the initial condition for the two-particle vertex. The perturbing term for probing spin-nematic order iŝ (12) which strengthens (weakens) the xx and yy (zz) component of the couplings J ij on all nearest-neighbor bonds and where 0 < |δ| J. This term induces a small bias towards the lowering of spin-rotational symmetry in such a way that spin isotropy is always retained for spin rotations in the xy plane; i.e., the spin-rotational symmetry is broken down from SU(2) to U(1). Similarly, the perturbing term for probing dimer order iŝ which strengthens the couplings J ij on all bonds in S [J ij → J ij + δ for i, j ∈ S] and weakens the couplings The bond pattern P ≡ {S p , W p } (the subscript "p" labels the strong and weak bonds corresponding to a pattern "P ") employed here specifies the spatial pattern of symmetry breaking one wishes to probe. These modifications amount to changing the initial conditions of the RG flow at large cutoff scales Λ. As Λ is lowered, we keep track of the evolution of all nearestneighbor spin susceptibilities χ ij . We then define the nematic response function for a given pair of nearestneighbor sites by where χ xx ij (χ zz ij ) are the correlators on the strengthened (weakened) bonds. Similarly, the dimer response function is given by where, χ Sp (χ Wp ) denotes χ ij ∈ S p (χ ij ∈ W p ). The normalization factor J/δ ensures that the RG flow starts with an initial value of η SN/VBC = 1. If the absolute value η SN/VBC decreases or remains small under the RG flow, the system tends to equalize, i.e., to reject the perturbation on that link, while, if η SN/VBC develops a large value under the RG flow, it indicates that the system is tending to develop an instability towards the probed nematic or valence-bond-crystal order.

B. Luttinger-Tisza method
The classical limit of a system of n quantum spins described by a Heisenberg model is achieved by first normalizing the spin operators by dividing them by their angular momentum S and then taking the limit S → ∞ [107,108]. This procedure yields the corresponding classical spin system wherein the spin operators in Eq. (1) are replaced by ordinary vectors of unit length at each lattice site i. For general interactions, the classical Hamiltonian to be minimized reads as where by i/j we denote the primitive lattice site separated by the lattice translation vectors R ij and α/β denotes the sublattice site index. The underlying primitive lattice of the pyrochlore lattice is the face-centered cubic lattice, and the pyrochlore structure is composed of four interpenetrating face-centered cubic lattices. The Luttinger-Tisza method [109][110][111] attempts to find a ground state of Eq. (16) by enforcing the spin-length constraint only globally, i |S 2 i | = S 2 n, where n is the total number of lattice sites, which is termed the weak constraint. This relaxed constraint implies that sitedependent average local moments are now permissible, which, strictly speaking, take us beyond the classical limit by approximately incorporating some aspects of quantum fluctuations [112].
To solve this relaxed problem, we decompose the spin configuration into its Fourier modesS α (k) on the four sublattices of the pyrochlore lattice Inserting this equation into Eq. (16) results in with the interaction matrix given bỹ The optimal modes satisfying the weak constraint are then given by the wave vector k, for which the lowest eigenvalue of Eq. (19) has its minimum. The eigenvector corresponding to this eigenvalue gives the relative weight of the modes on the sublattices [113], which means that the optimal modes do not fulfill the strong constraint (|S 2 i | = S 2 , i.e., fixed spin-length constraint on every site) if the components of the eigenvector do not have the same magnitude. If, however, this condition is met, the true ground state of the classical model is a coplanar spiral determined by the optimal Luttinger-Tisza wave vector [114]. There are also cases where one can construct an explicit parametrization of the ground state purely from the optimal modes in the pyrochlore lattice, as is the case with the cuboctahedral stack state described in Sec. III A 1.

C. Iterative minimization of the classical Hamiltonian
To find the ground state of the classical Heisenberg Hamiltonian in parameter regions where the Luttinger-Tisza method is not exact-i.e., a state constructed solely from the optimal modes does not fulfill the strong constraint-we employ an iterative minimization scheme which preserves the fixed spin-length (strong) constraint at every site [56]. Starting from a random spin configuration on a lattice with periodic boundary conditions, we choose a random lattice point and rotate its spin to point antiparallel to its local field defined by This rotation results in the energy being minimized for every spin update and thereby converging to a local minimum. We choose a lattice with L = 32 cubic unit cells in each direction, and thus a single iteration consists of 16L 3 sequential single-spin updates. One can therefore view this scheme as a variant of classical Monte Carlo with Metropolis updates at zero temperature, where we accept only optimal updates. This iterative scheme is carried out starting from ten up to 50 different random initial configurations per parameter set to maximize the likelihood of having found a global energy minimum. The exact number depends on convergence of the resulting energies. From the minimal energy spin configuration, the spin structure factor is computed, which is, up to a normalization constant, the same as the susceptibility defined in Eq. (11), but now for a finite system. Although it is not guaranteed that this scheme ends up in the global energy minimum, we find that, in all cases where an exact ground state is known, the iterative minimization scheme recovers the ground state, even when there exist nonoptimal states corresponding to local energy minima and having the same wave-vector content as the true ground state. This scheme also provides us with the opportunity to use spin configurations built from various (which can be arbitrarily chosen) parametrizations as a starting point of the minimization to check the quality of these parametrizations and also compare the competition between two states directly at a phase boundary.
As the iterative minimization works in direct space, we naturally see lattice symmetry breaking inherent to the ordered ground state, which cannot be captured by symmetry-preserving Fourier-space-based methods such as Luttinger-Tisza.
In the following section, we investigate the ground state of the general J 1 -J 2 Heisenberg model, both in the small spin-S regime (employing PFFRG) as well as the corresponding classical model using a combination of the Luttinger-Tisza method and iterative energy minimization schemes. We first begin with a discussion of the nearest-neighbor Heisenberg antiferromagnet.

III. THE NEAREST-NEIGHBOR HEISENBERG ANTIFERROMAGNET
We begin by investigating the ground state and behavior of the spin-spin correlation functions of the Heisenberg model with only a nearest-neighbor antiferromagnetic interaction and for a general pyrochlore lattice with nonzero breathing anisotropŷ where J up > 0 and J down > 0 are two different antiferromagnetic couplings on the nearest-neighbor bonds within the up and down tetrahedra, i.e., i, j up and i, j down , respectively. Hereafter, we parametrize these couplings in terms of a single angle ϕ and an overall energy scalẽ J: From a material perspective, the isotropic version of the model, i.e., ϕ = π/4, proves to be of relevance in understanding the low-temperature dynamics in chromium spinels [57,115]. On the other hand, the spatially anisotropic version of the model, wherein the up and down tetrahedra feature different exchange couplings, i.e., J down /J up = 1, the so-called breathing pyrochlore is realized in the recently synthesized spinels LiGaCr 4 O 8 , LiInCr 4 O 8 , LiInCr 4 S 8 , LiGaCr 4 S 8 , CuInCr 4 S 8 , and CuInCr 4 Se 8 [116][117][118][119][120][121][122][123][124][125][126][127][128][129][130] and in a pseudospin S = 1/2 Yb-based compound Ba 3 Yb 2 Zn 5 O 11 [131][132][133]. In these compounds, the magnetic Cr 3+ (Yb 3+ ) ions, which carry S = 3/2 (S = 1/2), form an alternating array of small and large tetrahedra, resulting in different exchange couplings for the two sets of tetrahedra. We begin by reviewing the established results for the classical Heisenberg antiferromagnet on the isotropic and breathing [130] pyrochlore lattices. While a number of the results given below have previously been published in the literature, reestablishing them here sets the stage for our own original results.
Henceforth, all temperatures for the isotropic classical and quantum models are expressed [00l] in units of J 1 S(S + 1) and J 1 , respectively (and we omit the factor of √ 2), while for the breathing model they are expressed in units ofJS(S + 1) andJ for the classical and quantum models, respectively. In the classical limit of Eq. (22), the Heisenberg spin operatorsŜ i reduce to standard three-component vectors S i . In the ensuing analysis, it proves convenient to introduce the magnetization M T of the T th tetrahedron, where the index α = 1, 2, 3, and 4 labels the four spins within the T th tetrahedron. In terms of M T , the Heisenberg Hamiltonian can be recast as a disjoint sum of the square of the magnetizations M T over the "up" and "down" tetrahedra, From Eq. (25), it follows that any state which satisfies the condition M T = 0 on each tetrahedron T is a classical ground state. The dimension of the ground-state manifold turns out to be countably infinite, which is best illustrated via a "Maxwellian counting argument" [3,4], which proceeds as follows: For a system of N s classical Heisenberg spins, we have the number of degrees of freedom F = 2N s (three degrees of freedom with one spin-length normalization constraint). In the ground state, all three components of M T should be zero on every tetrahedron, which gives the number of constraints K = 3N c , where N c is the number of tetrahedral clusters, and N s = 2N c (each tetrahedron has four spins, but each spin is shared between two tetrahedra). Hence, under the assumption that all constraints can be satisfied simultaneously and are all linearly independent, we arrive at the number of ground-state degrees of freedom D = F − K = 4N c − 3N c = N c which is an extensive quantity. If the constraints are not all linearly independent, then one underestimates D; however, for the pyrochlore Heisenberg antiferromagnet, it is known [3,4]  plane susceptibility shown in the inset is evaluated. Above this temperature, the width of the pinch points is seen to reproduce the T 1/2 behavior [4]. In Appendix A, we show how the inclusion of higher-order diagrammatic contributions in 1/S cure this spurious divergence.
that the corrections to the estimate for D are at most subextensive. The extensive (exp[O(L 3 )]) degeneracy of the ground-state manifold proves severe enough to preclude a finite-temperature phase transition, thus realizing a zero-temperature "cooperative paramagnet" [1] with nonzero entropy [134], referred to as a "classical spin liquid" [3,4,8,9,135]. Indeed, at low temperatures, the Heisenberg model not only fails to develop long-range dipolar magnetic order of the Néel type but also does not have conventional nematic order [3,4] of the type characterized by an order parameter which takes on its maximal value in a perfectly collinear state [104]. At T = 0, the classical spin liquid features critical, i.e., algebraic, spinspin correlations of dipolar character [136], which is a consequence of the local constraint that the magnetization M T on each tetrahedron is identically zero for any ground state [137][138][139][140][141]. These dipolar correlations most visibly show up in the Fourier transform of the two-spin correlator, where they form a pattern of bow ties [see Fig. 3(a)] with sharp singularities termed pinch points [see the encircled point in Fig. 3(a)] [3,9,27,70,142]. The dipolar nature of the correlations in the T → 0 regime is, in fact, a common feature of all classical O(N ) nearest-neighbor antiferromagnets for which the system remains paramagnetic down to T = 0 [140]. This fea-ture excludes the N = 2 (XY -spins) case, as this case is known to show a thermal order-by-disorder transition to collinear ordering for spins which have a global easy plane [3,4] as well as those with local sublatticedependent easy planes which are perpendicular to the local 111 axes [45,[143][144][145][146][147]. The limit N = 1 (Ising spins) is realized in various spin-ice materials A 2 B 2 O 7 (A ≡ Dy, Ho and B ≡ Ti, Sn) which, at low but nonzero temperatures, host a classical spin liquid featuring dipolar correlations and the associated pinch points [148]. Coming back to the case of N = 3 (Heisenberg spins) at finite temperatures, we note that thermal fluctuations lead to violations of the M T = 0 constraint and generate a finite correlation length ξ which, at low temperatures, diverges as T −1/2 [4]. At distances r ξ, the algebraic nature of the real-space spin-spin correlations changes into an exponential. Consequently, at finite temperatures the pinch points acquire a finite width ∼ 1/ξ [71] [see Figs. 3(a)-3(c)] which, at low temperatures, goes to zero as T 1/2 [4] [see Fig. 4].

Breathing case
As in the case of the isotropic pyrochlore lattice, the Heisenberg Hamiltonian in the presence of breathing anisotropy [Eq. (22)] can be straightforwardly recast as a disjoint sum of terms, each involving the magnetization M T [Eq. (24)] of a tetrahedron T : (26) It is clear that when J up and J down are both antiferromagnetic, any state in which M T = 0 on every up and down tetrahedron T is a classical ground state. Thus, in the presence of a breathing anisotropy, the extensive degeneracy of the isotropic model remains intact, and, consequently, the ground state at low temperatures remains a classical spin liquid [130]. However, as one moves away from the isotropic point ϕ = π/4, the appearance of the bow-tie pattern with a decreasing temperature, and the development of the pinch-point singularities in the limit T → 0, becomes progressively slower on approaching the decoupled tetrahedron limit, which is because the correlation length is proportional to the product J up J down /J 2 = cos φ sin φ [130], and, hence, the development of the correlations is slower when closer to the decoupled tetrahedron limit. In Fig. 3, we show the spin susceptibility profile for two values of the breathing anisotropy, ϕ = 3π/16 and ϕ = π/16, to enable a comparison with Fig. 8 in Ref. [130]. As expected, the development of the bow-tie pattern of scattering with sharp singularities as T → 0 becomes progressively slower as one moves towards the decoupled tetrahedron limit.
In the following section, we consider the regime of small spin S where strong quantum fluctuations are expected to significantly alter the ground state and nature of the spin-spin correlations.

Isotropic case
The investigation of the low-temperature (T J 1 ) physics of Eq. (22) in the small spin-S regime proves to be of utmost physical interest by virtue of the fact that in this limit the model harbors strong correlations which conspire with amplified quantum fluctuations to set the stage for a potential realization of a quantum spin liquid. However, it is precisely in this regime that the model acquires a notorious reputation for difficulties due to its nonperturbative character which makes the conclusions obtained from perturbative approaches unreliable [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]149]. Herein, we address this problem within the PFFRG framework, which is particularly suited for addressing this regime due to its nonperturbative character.
To probe the propensity of the system towards developing long-range magnetic order at any wave vector k, we track the evolution of the susceptibility χ(k) with Λ for all wave vectors k in the extended Brillouin zone (EBZ) of the pyrochlore lattice. As discussed in Sec. II A, the onset of magnetic long-range order at a particular k is signaled by the presence of kinks or cusps in the Λ flow of χ(k), whereas a smooth monotonically increasing behavior of χ(k) down to Λ → 0 points to a quantum-disordered ground state. For S = 1/2, we observe that the Λ evolution of the susceptibility χ(k) ∀ k ∈ EBZ [see Fig. 5(a) for the EBZ] is smooth and displays a monotonically increasing behavior down to Λ → 0 with no detectable signatures of an instability or a kink [see also Appendix B]. A numerical maximization of the susceptibility function in the EBZ finds feeble maxima at the high-symmetry W points, i.e., at k = 2π(2, 1, 0) [see Fig. 5(a)]. The RG flow of the susceptibility evaluated at the W point is shown in Fig. 5(b), wherein the smooth nature of the flow gives strong evidence in favor of a quantum paramagnetic ground state of the S = 1/2 quantum Heisenberg antiferromagnet on the pyrochlore lattice, in agreement with previous works [17-27, 29, 30].
The corresponding reciprocal space spin susceptibility profile in the EBZ evaluated at the lowest simulated temperature T /J 1 = 1/100 is shown in Fig. 5(c). The profile appears to be of a highly diffusive character along the edges and surfaces of the EBZ. So as to reveal the nature of the correlations, we plot χ(k) in the [hhl] plane (i.e., , wherein one clearly sees the characteristic bow-tie pattern, albeit with a softening and broadening of the pinch points due to quantum fluctuations [31,[149][150][151][152][153]. Indeed, in the small spin-S regime, the spin-flip exchange processes in the Heisenberg Hamiltonian become important and generate quantum fluctuations which dynamically violate the zero magnetization per tetrahedron constraint. Since it is this constraint which is ultimately responsible for the singular and perfectly sharp pinch points observed in the classical model, its violation in the quantum spin-1/2 model leads to a regularization or a softening of the pinch-point amplitude as their singular character disappears. In addition, quantum fluctuations also generate a finite correlation length ξ for the direct-space spin-spin correlations, such that at distances r ξ the dipolar nature of the correlations changes into an exponential. Consequently, the pinch points undergo "broadening," which can be quantified by their FWHM. Indeed, the FWHM is determined by the inverse of this correlation length, i.e., FWHM ∼ 1/ξ. In Fig. 7, we show the variation of χ(k) along the width of the pinch point, i.e., along the white vertical line in Fig. 6(a), and for S = 1/2 the FWHM of the pinch point is determined to be 1.6π at the lowest simulated temperature T /J 1 = 1/100.
Our finding of relatively rounded pinch points is in agreement with the results of Refs. [25,26,31], which also observe pinch points of a similar nature. The fact that the overall bow-tie pattern of susceptibility appears rather intact (despite relatively rounded pinch points) lends support to the view that the low-temperature (T /J 1 = 1/100) paramagnetic phase of the S = 1/2 nearest-neighbor Heisenberg antiferromagnet respects to a good degree the zero net magnetic moment per tetrahedron constraint, i.e., the "ice rules"-as also found in Ref. [31]. The temperature evolution of the susceptibility in the [hhl] plane is shown in Fig. 6. To obtain a quantitative picture, we plot in Fig. 8 the susceptibility along a 1D cut (white line in the T /J 1 = 1/100 plot of Fig. 6) across the width of the pinch point in the bowtie structure. On increasing the temperature by even an order of magnitude, i.e., up to T /J 1 = 1/10, it is found that the susceptibility profile and the width of the pinch points remain essentially unchanged. In the temperature range T /J 1 = 1/10 till T /J 1 ∼ 1, the pinch-point width is seen to increase (approximately) linearly (see the inset of Fig. 8) in contrast to the T 1/2 behavior expected classically (see Fig. 4). However, the fact that the overall bow-tie structure remains relatively intact up till T ∼ J 1 seems to suggest that the ice rules govern the physics (to a good degree of accuracy) over a surprisingly large temperature range as also found in Ref. [31]. We also study the behavior of the direct-space spin-spin correlations with the temperature and find that, for any given distance, it is only their amplitude that varies with the  temperature, while their signs remain constant over the entire temperature range, in agreement with the findings of Ref. [26]. Also, the signs of all correlators up to the 16th neighbor as obtained from PFFRG agree with those obtained in Table I of Ref. [26]. This agreement is interesting in light of the fact that Ref. [26] evaluates the equal-time spin-spin correlators, i.e., S(q, ω) integrated over the frequency, whereas we compute only the ω = 0 correlator, which implies that an integration over frequencies does not change the sign.
Early investigations into the nature of the ground state of the S = 1/2 nearest-neighbor Heisenberg antiferromagnet, predominantly based on perturbative approaches in the intertetrahedra coupling, found the ground-state to be a valence-bond crystal [17][18][19][20][21][22][23][24]. Using PFFRG, we probe for possible instabilities of the quantum paramagnet towards valence-bond-crystal formation. We consider three simple dimerization patterns which, respectively, break the translational symmetry along (i) all three tetrahedral axis directions (VBC 3D ), (ii) two tetrahedral axis directions (VBC 2D ), and (iii) one tetrahedral axis direction (VBC 1D ). The dimer response functions η P VBC [Eq. (15)] of all three VBCs are found to decrease under the RG flow [see Fig. 9(a) for the RG flow of η P VBC ] which lends support towards the scenario of a symmetric quantum-spin-liquid ground state as opposed to the previously proposed scenario of a VBC ground state. The disagreement between our findings and those of previous studies [17][18][19][20][21][22][23][24], which argue for a VBC ground state, is likely explained by the fact that a common thread of these approaches is the inherent symmetry breaking already built in to the scheme considered therein, which then biases the conclusion towards a VBC ground state. That being said, here we investigate VBCs only up to an eight-site unit cell, and the possibility of VBCs with larger unit cells cannot, in principle, be ruled out.
The possibility of the occurrence of spin-nematic order in the classical nearest-neighbor Heisenberg antiferromagnet is discussed in Refs. [3,4], wherein it is found that the system evades such nematic order [104]. Here, we investigate for the possibility of nematic order [see Sec. II A 2] in the S = 1/2 nearest-neighbor Heisenberg antiferromagnetic model. We plot the RG flow of the nematic response function η SN [Eq. (14)] in Fig. 10, wherein one observes that η SN remains less than one throughout the RG flow (albeit displaying nonmonotonic behavior) and sharply decreases at low temperatures (T J 1 ). This result indicates that the system tends to reject spontaneous breaking of SU(2) spin rotational symmetry via a quadrupolar order parameter in the ground state of the S = 1/2 nearest-neighbor isotropic Heisenberg antiferromagnet. Though our results are at variance with Ref. [154], which argues for a nematic quantum spin liquid featuring spin-nematic order in the S = 1/2 nearestneighbor Heisenberg antiferromagnetic model, we mention that, since we a priori exclude the fermionic fourparticle vertex from the RG equations and hence we cannot calculate the nematic susceptibility, our calculation of the nematic response function by applying symmetry breaking is approximative in character. Thus, we do not definitively exclude the possibility of the realization of a nematic quantum-spin-liquid ground state.

Breathing case
In a breathing pyrochlore system, the ratio of the interto intratetrahedra coupling J down /J up provides a convenient interpolation parameter which connects the decoupled tetrahedron and the isotropic limits. It is of interest to investigate the stability of the isotropic model ground state and the evolution of the spin-spin correlations as a function of J down /J up . The RG flow of the dominant susceptibility for different values of the breathing anisotropy is shown in Fig. 11(a), wherein we observe a smooth flow down to Λ → 0, in similarity with the finding for the isotropic model [see Fig. 5  to an extended region of parameter space (accessible by tuning J down /J up ) over which a quantum paramagnetic phase is stabilized. We also assess the stability of the paramagnetic phase against dimerization into the type of VBC orders considered for the isotropic model and find that the system rejects the applied symmetry breaking under the RG flow, hinting at a possible quantumspin-liquid state. In the strongly anisotropic limit, we cannot totally exclude the possible scenario of a ground state with more involved patterns of symmetry breaking, e.g., lattice nematic order or VBC with a larger unit cell. Indeed, in the S = 1/2 breathing kagome Heisenberg antiferromagnet, the situation is contentious: with one work finding VBC [155] while the other finds lattice nematic order [156]. So, further work on the (strongly) anisotropic breathing pyrochlore is probably warranted to ascertain whether it remains without VBC or lattice nematic order down to the limit of the decoupled tetrahedron. Furthermore, we find that the bow-tie pattern of scattering seen in the [hhl] plane is remarkably robust with regard to the introduction of breathing anisotropy, and the width of the bow tie increases only marginally even for strong values of anisotropy [see Fig. 12(a)]. This result shows that in the quantum paramagnetic ground state the low-energy physics is approximately governed by the ice rules.

Isotropic case
Increasing the spin S from S = 1/2 to S = 1 renders the effects of quantum fluctuations less pronounced, thus favoring conditions amenable for stabilizing long-range magnetic order. Previous investigations of the S = 1 Heisenberg antiferromagnet have not been able to reach an unambiguous conclusion regarding the presence or absence of magnetic order [19,157]. The Λ evolution of the susceptibility at the k vector where it has its maximum value, i.e., the high-symmetry W point, is shown in Fig. 13(a). The RG flow is not seen to exhibit any instabilities as would be signaled by the presence of kinks and, on the contrary, appears to be of a smooth character [see Appendix B for an analysis on the detection of possible magnetic instabilities in the S = 1 RG flow]. Similar flow behaviors of the susceptibility are exhibited for all wave vectors k ∈ EBZ. These observations lead us to the interesting conclusion that in increased spatial dimensionality (here, 3D) if geometric frustration is severe enough, such as on the pyrochlore lattice, then even for S = 1 quantum fluctuations are able to prevent the onset of long-range magnetic order in the Heisenberg antiferromagnet, thereby stabilizing a quantum paramagnetic ground state. The susceptibility profile in the [hhl] plane is qualitatively similar to the one obtained for S = 1/2; however, the pinch points become slightly sharper as reflected by the decrease in FWHM to 1.42π compared to 1.6π for S = 1/2, evaluated at the lowest simulated temperature T /J 1 = 1/100 [see Fig. 7].
To assess the stability of this paramagnetic phase against spontaneous dimerization, we study the dimer response functions of three candidate VBC states described in Sec. III B. The Λ evolution of the dimer response functions for the three VBCs [see Fig. 9(b)] shows that, similar to the S = 1/2 case, the system strongly rejects the corresponding applied symmetry breaking. With the present data, we cannot, as in the S = 1/2 case, rule out the possibility of VBCs with larger unit cells and more complicated patterns of symmetry breaking being stabilized. Nonetheless, from the current PFFRG results, the predicted ground state would be a quantum spin liquid.

Breathing case
Upon tuning a breathing anisotropy, i.e., J down /J up = 1, we observe that the RG flows [see Fig. 11(b)] do not develop any signatures of a kink or an instability [as inferred from an analysis based on the method of detection of instabilities as explained in Appendix B] down to the strongly anisotropic limit and remain smooth as Λ → 0, pointing to the absence of magnetic long-range order. Thus, our results show that even for S = 1, where quantum fluctuations are expected to be less pronounced, there exists an extended region in parameter space hosting a quantum paramagnet which can be accessed from the isotropic point (J down /J up = 1) by tuning the breathing anisotropy. We probe this paramagnetic phase for possible VBC instabilities, and find that the system rejects the applied symmetry breaking; however, as in the case of S = 1/2, we do not exclude the possibility of a ground state featuring a more elaborate pattern of symmetry breaking [37]. We also observe that the bow-tie pattern and the pinch-point width remain essentially unchanged compared to the isotropic model [see Fig. 12(b)], indicating that the ice rules continue to dictate the low-energy physics of the quantum paramagnetic ground state even for strong breathing anisotropy.

D. Large spin-S regime
As quantum fluctuations decrease in strength with increasing spin S, magnetic long-range order might be expected to ultimately prevail. Indeed, we find that, for S = 3/2, the RG flow of the dominant susceptibility [see Fig. 13 Based on this analysis [see Fig. 25], we conclude that for S = 3/2 and beyond there is an onset of magnetic long-range order in the nearest-neighbor isotropic Heisenberg antiferromagnet. It is worth emphasizing that, for the finite S values studied in our manuscript, the correct balance between leading 1/S terms and subleading contributions is already incorporated in the PFFRG [see Sec. II A]. For this reason, the PFFRG at any finite S is still well justified even if plain RPA in the large-S limit, i.e., treating only leading 1/S diagrams, produces the aforementioned artifact of finite-temperature divergence of the susceptibility [see Fig. 4]. However, with increasing S, the PFFRG becomes numerically more challenging (and also more sensitive to errors), because it becomes progressively difficult to account for the proper interplay between (large) leading 1/S and (much smaller but still important) subleading terms in our numerical algorithm. For this reason, we applied the PFFRG only to "moderate" spin magnitudes smaller than eight and use plain RPA in the infinite-S limit [see Appendix A]. Therefore, we are unable to comment on the long-standing issue of the presence or absence of long-range magnetic order in the large-S quantum Heisenberg antiferromagnet.
Determining the precise nature of the magnetic order (if any) for intermediate values of S constitutes an intriguing and challenging question which has remained unanswered to date. The problem of the ground state of the large-S quantum antiferromagnet on the pyrochlore lattice is addressed extensively using effective Hamiltonian approaches [10][11][12][13][14][15][16]. However, due to the weak selection effects operating at both the harmonic and anharmonic level, no definitive conclusion on the nature of the ground state has yet been reached. Addressing this problem within the PFFRG scheme, we study the evolution of the spin susceptibility profile with increasing values of S in order to figure out whether quantum fluctuations are successful in distilling a unique (magnetically ordered) ground state with a given wave vector k ∈ EBZ out of the extensively degenerate classical ground-state manifold. In Fig. 14, we show the variation in the susceptibility along a path passing through the high-symmetry points [see Fig. 5(a)] for increasing S values. One observes that, while the susceptibility increases with increasing S, there is no clear enhancement at any given wave vector, and the susceptibility profile evaluated at and above the critical breakdown temperature in Figs. 14(b)-(d) remains essentially unchanged compared to that of the S = 1/2 and S = 1 paramagnetic phase, with just an overall enhancement. The absence of pronounced Lorentzian peaks points to the fact that the quantum order-by-disorder selection effects as captured by one-loop PFFRG [72] may be extremely feeble down to the lowest cutoff or temperature considered, even upon the inclusion of higher orders in 1/S embedded within the PFFRG calculation framework [74]. It will be interesting to investigate the large-S limit beyond one loop formulations of PFFRG, e.g., by employing the recently formulated multiloop PFFRG which sums up all parquet diagrams to arbitrary order in the interaction [158][159][160].

A. Classical phase diagram
Given the absence of long-range order at a nonzero temperature in the classical nearest-neighbor Heisenberg pyrochlore antiferromagnet, any weak perturbations to that model have strong effects on the thermodynamic and magnetic properties of the system that may result in, e.g., magnetic long-range ordering. Indeed, the inclusion of a second-nearest-neighbor Heisenberg coupling J 2 to the classical nearest-neighbor Heisenberg model on the pyrochlore lattice is known to stabilize a plethora of intricate magnetic orders [see Table I and Fig. 15], part of which is investigated in Refs. [38][39][40][41], with a full exploration of the J 1 -J 2 parameter space reported in Ref. [56]. Despite the fair number of results available in the literature for this classical J 1 -J 2 model, we find and report below some corrections and/or amendments to the current knowledge about the classical phases of this system.
We find the J 1 -J 2 model to host seven different classical magnetic orders, in addition to a classical spin-liquid (cooperative paramagnetic) phase found for the nearestneighbor antiferromagnetic model. Employing an approach which combines a Luttinger-Tisza analysis with an iterative energy minimization on large system sizes of  Table I for a description of the phases and the location of the phase boundaries. 32×32×32 cubic unit cells (i.e., 524 288 spins), we present a refined analysis of the classical phase diagram and the nature of its magnetic orders. The principal differences in our findings compared to those presented in Ref. [56] can be attributed to the substantially reduced finite-size effects in our calculations compared to those of Ref. [56], which are based on a 4 × 4 × 4 cubic unit cell (1024 sites) system. In addition, we identify within the EBZ of the pyrochlore lattice the ordering wave vectors of the classical magnetic orders [see Table I] as would be determined in neutron-scattering experiments. It is important to discuss these states in detail here, since, as we will see in the next section, the quantum (S = 1/2 and S = 1) models harbor the same long-range ordered states.
The pure nearest-neighbor Heisenberg antiferromagnet (J 2 = 0) features an extensively degenerate manifold of classical ground states whose sole shared feature is that the sum of the spins on every tetrahedron is identically zero [see Sec. III A 1]. It is shown in Ref. [2] that an infinitesimal amount of antiferromagnetic second-nearestneighbor coupling J 2 > 0 proves sufficient to partially lift this degeneracy by selecting a nonextensive subset of the ground states of the pure nearest-neighbor antiferromagnet. These states are such that the spins within each of the four face-centered cubic (fcc) sublattices of the pyrochlore lattice order ferromagnetically, and therefore this state is dubbed k = 0. However, the sublattices are not aligned parallel to each other, but the state preserves the constraint of zero spin sum per tetrahedron, resulting in an ordering wave vector at k = 2π(2, 0, 0) and symmetry-related points in the EBZ. This result can perhaps be most easily understood by noting that a secondnearest-neighbor interaction J 2 is equivalent to a thirdnearest-neighbor interaction J 3 of the opposite sign, i.e., J 3 = −J 2 , as long as every tetrahedron satisfies the zero spin sum ("ice rule") constraint [40]. Since J 3 couples only spins on the same sublattice, it is straightforwardly optimized by selecting states with ferromagnetic ordering within each sublattice. This state turns out to be an exact Luttinger-Tisza eigenstate of theJ k αβ matrix in Eq. (19) with an energy per spin E = −2J 1 − 4J 2 . Given that the ordering is fixed only within each sublattice separately, there remains the freedom of choosing the relative orientation of the individual ferromagnetically aligned sublattices while respecting the zero spin sum per tetrahedron constraint. Hence, at T = 0 there exists a ground-state degeneracy characterized by three angular degrees of freedom. Therefore, the distribution of spectral weight between the dominant k = 2π(2, 0, 0)type vectors is not fixed. At T = 0, the breaking of the cubic pyrochlore symmetry is not energetically determined by the interactions; however, for finite temperatures entropic effects could select a unique ground state. The relative weights of the dominant peaks in the structure factor then serve as a measure of the collinearity of the sublattices, with the case of only one of them being present corresponding to a fully collinear state. Irrespective of the relative orientation of the sublattices, the ferromagnetic correlations within each of these manifest themselves in the spin structure factor by subdominant peaks of equal intensity at all of the 2π(1, 1, 1) points at the edge of the EBZ. The spectral weight of any one of the given subdominant peaks is exactly one-eighth of the total weight of the dominant peaks. The aforementioned k = 0 state minimizes the energy only in the regime where antiferromagnetic J 1 > 0 is dominant over sufficiently weak antiferromagnetic J 2 . Since the J 2 bonds are twice as many as the J 1 bonds, the J 2 interaction becomes dominant when J 2 /J 1 > 1/2 (θ 26.56°), resulting in a phase transition to a planar spiral ground state with one of the symmetry-related k = 2π(k, 0, 0)-type wave vector as the ordering wave vector. This state is also an eigenstate of the Luttinger-Tisza matrix Eq. (19), thus giving the exact expression k = (2/π) arccos[−J 1 /(4J 2 )−1/2] for the wave vector and an energy per spin of E = −J 2 1 /(2J 2 ) − 6J 2 . This wave vector differs from the one given in Ref. [56] by a factor of 2, which is due to the fact that the transformation done on this state to map it into an equivalent spin-chain model [161] was apparently not performed correctly. The pure second-nearest-neighbor antiferromagnet (J 1 = 0, J 2 = 1) also falls into this region and has a 120°spiral structure on each fcc sublattice. Taking into account the relative phases of the spirals between the sublattices, we find a resulting ordering wave vector k = 2π(4/3, 0, 0) in the EBZ of the pyrochlore lattice. In the planar spiral, and corresponding to the aforementioned dominant peaks at k = 2π(k, 0, 0)-type ordering wave vectors, there also exist subdominant peaks at ordering wave vectors of the k = 2π(3 − k, 1, 1) type in the EBZ. The k and 3 − k entries of the dominant and subdominant ordering wave vectors, respectively, always appear in the same component for each of these wave-vector pairs. The subdominant peaks are a signature of the correlations within the fcc sublattices of the pyrochlore lattice and have a fixed relative amplitude of one-quarter of the dominant peak.
The planar spiral order is stable against J 1 < 0 now becoming ferromagnetic (keeping J 2 > 0 antiferromagnetic), up to J 2 /J 1 = −0.68 (θ ≈ 145.78°). Beyond that point, the ground state changes to a noncoplanar structure, the so-called double-twist (DT) state, first uncovered in a frustrated antiferromagnet on an octahedral lattice [161]. Its name derives from the fact that the spins form two different kinds of spirals in two perpendicular directions but both governed by the same type of wave vector. In reciprocal space, this state features two pairs of k = 2π(3/4, 3/4, 0)-type wave vectors on different reciprocal space planes; the first pair, e.g., could be located in the k x -k y plane with k = 2π(3/4, 3/4, 0) and k = 2π(3/4, −3/4, 0), while the second pair, e.g., could be located in the k y -k z plane with k = 2π(0, 3/4, 3/4) and k = 2π(0, 3/4, −3/4). In the first plane, e.g., the k x -k y plane, two dominant peaks in the structure factor are located at the aforementioned wave vectors and have identical spectral weight. In the second plane, e.g., the k y -k z plane, subdominant peaks with approximately 59% of the spectral weight of the dominant ones are located at the aforementioned wave vectors. An approximate parametrization of such a state is given in Ref. [56]. Both pairs of wave vectors control the ordering on the individual fcc sublattices. The relative orientations of the spins on the sublattices lead to the appearance of additional subdominant peaks at k = 2π(5/4, 5/4, 0)-type wave vectors. For example, corresponding to the pair of dominant peaks in the k x -k y plane, there appear a pair of subdominant peaks at wave vectors k = 2π(5/4, 5/4, 0) and k = 2π(5/4, −5/4, 0) carrying approximately 29% of the amplitude of the dominant peaks. Similarly, corresponding to the pair of subdominant peaks in the k y -k z plane, there appear a pair of weaker peaks at wave vectors k = 2π(0, 5/4, 5/4) and k = 2π(0, 5/4, −5/4) carrying approximately 13% of the amplitude of the dominant peaks (in the k x -k y plane). The particular choice of planes chosen for the dominant and subdominant planes is not fixed by the Heisenberg model, but is determined by the spatial symmetry breaking when entering this phase.
Decreasing antiferromagnetic J 2 > 0 further, we encounter a phase transition at J 2 /J 1 ≈ −0.475(5) (θ ≈ 154.59°) to a state which is similar to the multiply mod-ulated commensurate spiral of Ref. [56], for which the transition point is estimated to be J 2 /J 1 ≈ −0.43. In reciprocal space, this state is characterized by the presence of four dominant commensurate ordering wave vectors of the k = 2π(3/4, 1/2, 1/4) type in the EBZ, for all of which the 1/2 component is in a common direction. We also find subdominant ordering vectors of the k = 2π(3/4, 0, −3/4) type; the zero component is the one which is 1/2 in the dominant k = 2π(3/4, 1/2, 1/4) wave vectors. This result is a consequence of a magnetic structure wherein the spins trace out multiple spirals in different directions in direct space which are controlled by the above wave vectors. Our refined analysis reveals that the observed commensurability of the wave vectors found in Ref. [56] is an artifact of large finite-size effects at play in that work. The imposition of periodic boundary conditions in the simulation of a L × L × L cubic unit cell system allows only those k vectors whose components are integer multiples of 2π/L. This restriction implies that an incommensurate ordering wave vector which is proximate to a commensurate one leads to an observed peak at the commensurate position. Indeed, we find that, for J 2 /J 1 ≈ −0.47, the four incommensurate ordering wave vectors of k = 2π[0.81(2), 0.50(2), 0.19(2)] type evolve continuously (at least within the used k-space numerical resolution of 2π/32) towards the commensurate values which are taken on at the transition point to the cuboctohedral stack (CS) state in Fig. 15. At the same time, the subdominant ordering vector stays unchanged, but its weight relative to the weight of the dominant peak varies from approximately 26% at its border with the DT state to approximately 32% at its border to the CS state. Our calculations show that, while the manner in which dominant and subdominant wave vectors control this state does not change, the dominant wave vector it is composed of does evolve as a function of J 2 /J 1 . Our findings are also supported by a Luttinger-Tisza analysis, which shows that there are incommensurate wave vectors with slightly lower energy close to the commensurate point. In this parameter regime, the Luttinger-Tisza state does not fulfill the strong spin-length constraint [see Sec. II B] but needs to be supported by the subdominant wave vectors we find, in order to be able to construct a normalized state. Because of the incommensurability of the dominant wave vector, we simply refer to this state as a multiply modulated spiral (MMS).
At J 2 /J 1 = −0.3965(5) (θ ≈ 158.37°), the MMS state evolves into the CS state [56,161]. Its name derives from the fact that, in a construction of the pyrochlore lattice as a stacking of alternating kagome lattice and triangular lattice layers in a [111] direction, the spins in each kagome layer are arranged such that they point towards the 12 vertices of a cuboctahedron, forming a 12sublattice magnetic structure first found on the kagome lattice [162,163]. At the same time, the spins on the triangular layers point to the eight midpoints of the triangular faces of the same cuboctahedron. This noncoplanar state is built up from any three wave vectors of the k = 2π(1/2, 1/2, 1/2) type, e.g., k = 2π(1/2, 1/2, 1/2), k = 2π(−1/2, 1/2, 1/2), and k = 2π(1/2, −1/2, 1/2) with identical spectral weight, and is stacked along the [111] direction parallel to the fourth wave vector of this type, e.g., k = 2π(1/2, 1/2, −1/2). The spin configuration in this state can be expressed analytically (see Ref. [56]). Each of the dominant ordering vectors is accompanied by a subdominant wave vector of k = 2π(1/2, 1/2, 3/2) type with approximately 18% of the spectral weight of the dominant vectors. From the parametrization, it follows that the average energy per spin, E = J 1 (3/4 + √ 6/2), is independent of J 2 (an extensive discussion how this originates from the state can be found in Ref. [56]). Thus, decreasing J 2 further does not change the energy of this state but, rather, lowers the energy of competing states.
At J 2 /J 1 = (−3/8 + √ 6/12) (θ ≈ 170.30°), the energy of the ferromagnet becomes lower than that of the CS state and occupies the largest extent of the J 1 -J 2 parameter space. Just as for the k = 0 state, the ferromagnetic ordering within the sublattices features subdominant ordering wave vectors at all the k = 2π(1, 1, 1)-type points in the EBZ, which have a spectral weight of one-quarter of the dominant k = 2π(0, 0, 0) vector. The pure J 2 ferromagnet proves to be fairly robust against moderately strong antiferromagnetic J 1 coupling.
For J 2 /J 1 −1.09 (θ ≈ 312.53°), the antiferromagnetic J 1 exchange destroys the ferromagnetic order, and a phase transition occurs to a family of states dubbed the Kawamura states after the group which investigated them in great detail [38]. This phase is made up of a family of degenerate ground states with dominant incommensurate wave vectors around the k = (k, k, 0) points with k ≈ 2π(5/4) and subdominant ones at k ≈ 2π(3/4) having approximately 22% of the spectral weight of the dominant vectors. In addition, we find stronger subdominant ordering at k ≈ 2π(1, 1/4, 7/4)-type vectors with approximately 55% spectral weight. There are two classes of ground states, composed of either four or all six of the ordering wave vectors, the latter therefore respecting the cubic symmetry of the pyrochlore lattice. In the case of a ground state composed of four of the six wave vectors, the Heisenberg model a priori does not determine which four are selected. A common feature of both these states is that they are superpositions of spirals with the pertinent wave vectors which, when combined, realize a noncoplanar state. The parameter k for the dominant ordering starts with a value k ≈ 2π(1.31) at the phase boundary to the ferromagnetic state J 2 /J 1 = −1.09 and approaches k = 2π(5/4) as J 2 → 0. The Kawamura states also approximately fulfill the zero spin sum per tetrahedron constraint, so they can likewise be considered as perturbed eigenstates of the pure J 1 -only antiferromagnetic model.

B. Quantum Phase Diagram
The regime of small spin S in highly frustrated magnets harbors strong quantum fluctuations which display intriguing effects such as (i) melting magnetic orders to potentially realize a quantum spin liquid, (ii) fostering the birth of new kinds of magnetic orders, (iii) shifting the pitch vector of spiral magnetic states, and (iv) shifting the phase boundaries relative to that found for the same Hamiltonian in its classical S → ∞ limit. With the aim of investigating these possibilities, we carry out a study of the quantum phase diagram of the J 1 -J 2 Heisenberg pyrochlore model for low values of spin S, which, to the best of our knowledge, had not been performed before the present work. We first address the important question concerning the possibility of stabilizing a quantum paramagnetic phase in the presence of a J 2 coupling. At the classical level, and as discussed in the previous section, it is shown [2] that the presence of an infinitesimal further neighbor J 2 coupling induces long-range magnetic order at low temperatures. However, strong quantum fluctuations in the small-S regime may destabilize those classical magnetic orders. Therefore, the question arises, in what range of J 2 /|J 1 |, with either antiferromagnetic J 1 > 0 or possibly even ferromagnetic J 1 < 0, may a quantum-spin-liquid phase be potentially realized.
By employing PFFRG, we map out the full J 1 -J 2 quantum phase diagram for S = 1/2, S = 1, and S = 3/2, which is shown in Fig. 16. Our most important finding, which is the main result of our work, is the presence of an extended quantum paramagnetic phase for the S = 1/2 model [see Fig. 16(a)] and, perhaps surprisingly, also for the S = 1 model [see Fig. 16 , while, for S = 1, its span is reduced by half to −0.11(2) J 2 /J 1 0.09(2) but remains nonetheless appreciable. For S = 1/2, we show the representative RG flows within the paramagnetic regime for a point in the antiferromagnetic J 2 regime [ Fig. 17(a)] and one in the ferromagnetic J 2 regime [ Fig. 17(b)]. These display a smooth and monotonically increasing behavior with no signatures of a kink, pointing to the absence of magnetic long-range order. The paramagnetic character of the ground state also shows up in the spin susceptibility profile in the form of an absence of sharp maxima in the EBZ which would be a signature of incipient Bragg peaks (IBPs) marking the onset of magnetic long-range order, along with a diffuse spectral weight caused by quantum fluctuations. Indeed, the antiferromagnetic J 2 spin susceptibility profile [see Fig. 17(c) for the S = 1/2 result] displays weak maxima at k = 2π(2, 0, 0) (and symmetry-related points), which correspond to the dominant Bragg peak wave vectors of the underlying k = 0 parent classical magnetic order (see Sec. IV A). Similarly, the spin susceptibility profile for ferromagnetic J 2 [see Fig. 17(d) for the S = 1/2 result] features a smeared distribution of spectral weight forming homogeneous ringlike features on the surface of the Brillouin zone (see Fig. 18 for the [hhl] plane scattering profiles). Classically, this parameter regime hosts the Kawamura magnetic order with dominant and subdominant Bragg peaks at k ≈ 2π(5/4, 5/4, 0) and k ≈ 2π(3/4, 3/4, 0) (and symmetry-related points). A comparison of the S = 1/2 paramagnetic spin susceptibility profiles, i.e., Fig. 17(c) for J 2 antiferromagnetic and Fig. 17(d) for J 2 ferromagnetic, with those of the respective parent classical magnetic orders, i.e., k = 0 [ Fig. 19(b)] and Kawamura [ Fig. 19(h)] states, lends support to the view that the quantum paramagnetic ground state may be viewed as a molten version of the parent magnetic orders under the action of quantum fluctuations.
The inclusion of a J 2 coupling also substantially modifies the nature of the paramagnetic scattering profile at low temperatures (see Fig. 18 for the [hhl] plane scattering). We find that for antiferromagnetic J 2 > 0 there is an enhancement of the pinch-point scattering amplitude as found in the corresponding classical model [71], while for ferromagnetic J 2 the scattering intensity at the pinch points is strongly suppressed and instead redistributes to form a hexagonal cluster pattern of scattering [71]. In Fig. 20, we plot the relative weight of the susceptibility (at T /J = 1/100) with respect to its value at the pinch point, i.e., (χ/χ pinch point ) along a 1D cut (marked by a white line in Fig. 18). This plot clearly reveals the degree of enhancement at the pinch point as an antiferromagnetic J 2 coupling is cranked up, while, for ferromagnetic J 2 , we see clearly the drifting of the maxima of susceptibility away from the pinch point and its enhancement at the wave vectors of the Kawamura state. The overall structure of the paramagnetic scattering profile is seen to be robust up to high temperatures T /J ∼ 1 [see Fig. 18]. Although the above results and discussions are for the quantum paramagnet in the S = 1/2 model, the findings for the S = 1 model differ only quantitatively, and the entire discussion for S = 1/2 holds true for S = 1, albeit for the smaller collective paramagnetic regime of the S = 1 model.
We now move on to the discussion of the magnetically ordered phases in the low-spin regime of the J 1 -J 2 model. A comparison of the classical and quantum phase diagrams in Fig. 16 shows that all the classical magnetic orders are present in the low-spin regime of the model and that no new magnetic orders are found to be stabilized by quantum fluctuations, as is found for the Heisenberg model on the square lattice [164]. Starting our discussion with the k = 0 order, we find that its span is considerably diminished for the S = 1/2 model [see Table I for phase boundaries], due to the fact that it gives way to an extended spin liquid phase around the J 2 = 0 point. The RG flow of the dominant susceptibility evaluated in the middle of the k = 0 phase [J 2 /J 1 ≈ 0.36, marked by a black disk in Fig. 16(a)] clearly shows signature of an instability [see Fig. 21], indicating the onset of k = 0 mag-  Fig. 16(a). Also shown, the Brillouin zone, a "truncated octahedron," with the high-symmetry points labeled.
netic order with a Néel temperature of T c /J ≈ 0.39 (2) which is given by the position of the instability, marked by an arrow in Fig. 21. The spin susceptibility profile evaluated for J 2 /J 1 = 0.36 at the instability point is shown in Fig. 19(b), wherein one observes the dominant IBP at the high-symmetry X points [ Fig. 19(a)], i.e., k = 2π(2, 0, 0) (and symmetry-related points), and the subdominant peaks at the L points [ Fig. 19(a)], i.e., k = 2π(1, 1, 1), and symmetry-related points, are also seen to be clearly resolved. Although both thermal and quantum order from fluctuation effects (order by disorder) are in principle captured in our simulations [102], we cannot make a statement about the collinearity of the ground state, as the PFFRG in its current formulation does not allow for lattice symmetry breaking; i.e., all symmetry-related IBPs have the same height. As discussed in Sec. IV A, classically, the collinear k = 0 state is selected by thermal fluctuations [40], and quantum fluctuations are likely to select the same state [7].
The k = 0 state undergoes a phase transition at J 2 /J 1 = 1/2 to an incommensurate planar spiral magnetic order. The RG flow of the dominant susceptibility 21. RG flows of the spin susceptibility for S = 1/2 at the ordering wave vectors of the seven magnetically ordered phases evaluated at the data points marked by black disks in Fig. 16(a). The points at which the solid lines become dashed (marked by arrows) indicate an instability in the flow, indicating an onset of magnetic order. evaluated deep inside the spiral ordered phase [J 2 = 1, marked by a black disk in Fig. 16(a)] features an instability at T c /J 2 ≈ 0. 73(3) [marked by an arrow in Fig. 21] pointing to the onset of magnetic order at this temperature. The corresponding spin susceptibility profile evaluated at the instability point is shown in Fig. 19(c). We find that the effect of quantum fluctuations on the planar spiral order is twofold: (i) it leads to a shift of the spiral wave vector compared to its classical value [165] and (ii) is found to increase the region of stability of the planar spiral beyond its classical domain. First, concerning the shift in the spiral wave vector, we show in Fig. 22 its evolution across its domain of existence for the classical and the quantum models. The wave vector is found to decrease monotonically as one traverses the spiral domain starting from its boundary with the k = 0 to the DT magnetic order. Meanwhile, the shift |δk| ≡ |k qu | − |k cl | from the classical k cl wave vector to the quantum k qu wave vector changes nonmonotonically across the domain of the planar spiral ordered phase [see the inset in Fig. 22]. For the most part of the spiral ordered regime, we find that quantum fluctuations always increase the wave-vector value, leading to more antiferromagnetic types of order. The shift δk achieves a maximal value of approximately 4% of the classical value near the boundary to the k=0 order. Second, concerning the increase in the region of stability of the planar spiral order, we find that there is a strong renormalization of the phase boundary of the planar spiral with the DT order, which gets shifted from its classical value of J 2 /J 1 ≈ −0.68 to J 2 /J 1 ≈ −0.537(6) for the S = 1/2 model [see Fig. 16(a) and Table I], implying a significant enhancement of the domain of existence of the planar spiral order.
At J 2 /J 1 = −0.537 (6), the planar spiral gives way to the DT magnetic order, whose RG flow evaluated at J 2 /J 1 ≈ −0.43 [marked by a black disk in Fig. 16(a)] and tracked at the dominant wave vector becomes unstable at T /J ≈ 0.39 (2) [marked by an arrow in Fig. 21]. The corresponding spin susceptibility profile is shown in Fig. 19(d), wherein, besides the dominant one, the subdominant peaks are also clearly resolved. We find that the DT phase in the S = 1/2 model occupies a similar extent in parameter space as in the classical model, albeit with displaced phase boundaries. As the ratio J 2 /J 1 is lowered, we find that at J 2 /J 1 = −0.347(2) the susceptibility at the ordering wave vectors of the MMS phase becomes stronger compared to that at the DT ordering wave vectors, and the MMS order is stabilized. However, the extent of the MMS phase in the S = 1/2 model is reduced to approximately one-third of its classical extent and thus now occupies only a tiny sliver in parameter space. Just as in the classical model, the IBPs of the quantum model are still located at incommensurate wave vectors, which are, however, shifted compared to those of the classical model. In Fig. 21, we show the RG flow evaluated at the optimal quantum wave vectors for J 2 /J 1 ≈ −0.335 [marked by a black disk in Fig. 16(a)], which reveals the onset of magnetic order at a Néel temperature of T c /J = 0.39 (2). The associated spin susceptibility profile is shown in Fig. 19(e). At J 2 /J 1 = −0.326(2), the MMS phase ends and the susceptibility at the CS order wave vectors becomes dominant. The CS phase for S = 1/2 has an appreciable extent in parameter space comparable to the classical model but with shifted phase boundaries. The spin susceptibility profile evaluated for J 2 /J 1 = −0.24 [see Fig. 19(f)] shows that the dominant IBP is located along the line joining the origin and the high-symmetry L point, and that the peak undergoes substantial smearing due to quantum fluctuations. The instability feature at T /J = 0.39 (2) in the RG flow [ Fig. 21] appears feeble, possibly hinting at the "weakness" of the CS magnetic order. It is of interest to note that the analogous cuboctohedral kagome orders [162,163] found in Heisenberg models with longrange interactions also display an extremely feeble signal of an instability in their RG flow [89,90].
Finally, as we lower J 2 /J 1 further, the ferromagnetic J 1 coupling becomes dominant enough to drive the system into a ferromagnetic ordered state which onsets at J 2 /J 1 = −0.153 (5). On comparison with the classical transition boundary at J 2 /J 1 ≈ −0.171, we see that the antiferromagnetic CS order intrudes into a portion of the phase diagram occupied by the ferromagnetic order at the classical level, as expected from general considerations [93,166]. For the J 1 = −1-only model [marked by a black disk in Fig. 16(a)], we show the RG flow of the k = (0, 0, 0) susceptibility in Fig. 21, wherein we observe a strong signal of an instability. We obtain an estimate of the critical (Curie) temperature T c /|J 1 | = 0.77(4), which is equal within two error bars to the Quantum Monte Carlo value of T /|J 1 | = 0.718 [167] [see Table II for a comparison with other methods]. In Table II, we also provide for a comparison the Curie temperatures of the simple cubic lattice which has the same coordination TABLE II. The ordering (Curie) temperatures for the S = 1/2 nearest-neighbor quantum Heisenberg ferromagnet (in units Tc/|J1|) (columns 2 and 3) and its corresponding classical (S → ∞) model (in units of Tc/[|J1|S(S + 1)]) (columns 5 and 6) on the pyrochlore and simple cubic lattices as obtained by PFFRG and compared with estimates obtained from quantum Monte Carlo (QMC), classical Monte Carlo (CMC), high-temperature expansion (HTE), rotation-invariant Green's function method (RGM), and random-phase-approximation (RPA). The fact that T pyro c /T SC c < 1 can be attributed to finite-temperature frustration effects [167]. We also quote the result in the mean-field approximation (MFA), which is insensitive to the difference between the pyrochlore and simple-cubic lattice, since it depends only on the coordination number.

Method
Pyrochlore number z = 6 as the pyrochlore lattice but is bipartite. It is of interest to observe that, for both the S = 1/2 and classical (S → ∞) models, the Curie temperature of the pyrochlore lattice is lower compared to the simple cubic lattice, a fact which can be attributed to finite temperature frustration effects [167][168][169][170][171].
The spin susceptibility profile [see Fig. 19(g)] also reveals the presence of subdominant IBPs at the L point besides the dominant peak at the Γ point. As expected, the ferromagnetic phase occupies an entire quadrant of the phase diagram spanning from the limit J 1 = −1 till J 2 = −1 and gets destabilized only when a significant antiferromagnetic J 1 coupling is added to the J 2 = −1 ferromagnetic model. Our PFFRG calculations identify the value of J 2 /J 1 = −1.252(5) when the ferromagnetic order gives way to the antiferromagnetic Kawamura state, whereas classically the transition occurs at J 2 /J 1 −1.09. Herein, similar to the CS state, we observe that quantum fluctuations extend the region of stability of the antiferromagnetic Kawamura order at the cost of the ferromagnetic state [93,166]. The optimal wave vectors of the Kawamura state evolve within the region it occupies in the phase diagram; however, their value remains close to 2π(5/4, 5/4, 0). In Fig. 21, we show the RG flow of the susceptibility evaluated at the optimal wave vectors for J 2 /J 1 = −0.634(4) [marked by a black disk in Fig. 16(a)]. The signature of an instability is not very pronounced and appears to be located around T c /J = 0.54 (2). The corresponding spin susceptibility profile is shown in Fig. 19(h), wherein one observes that quantum fluctuations cause a significant diffusing of the spectral weight for both the dominant and subdominant IBPs [99].
The quantum phase diagram for the S = 1 model [see Fig. 16(b)] appears qualitatively similar to the one for S = 1/2, with the only differences being quantitative ones, such as the location of the phase boundaries, value of optimal wave vectors, etc. As we gradually increase the value of the spin S, we see that the quantum phase diagram starts going over into the classical one, as is already manifestly apparent for S = 3/2 [see Fig. 16(c)].

V. SUMMARY
In this paper, we employed the PFFRG method to investigate the long-standing problem of the effects of quantum fluctuations on the pyrochlore lattice for generic spin S in a Heisenberg model with nearest-neighbor J 1 and second-nearest-neighbor J 2 couplings. For the spin S = 1/2 nearest-neighbor Heisenberg antiferromagnetic model with spatially isotropic couplings, we find a quantum paramagnetic ground state [Sec. III B 1]. The paramagnet appears robust against potential instabilities towards the formation of either a valence-bond crystal [ Fig. 9(a)] or spin-nematic order [ Fig. 10], thus providing evidence in support of a quantum-spin-liquid ground state. The reciprocal space susceptibility plotted in the [hhl] plane displays the characteristic bow-tie pattern [Fig. 5(d)]. However, the dynamic violation of the zero magnetization per tetrahedron constraint due to quantum fluctuations manifests itself as (i) a regularization or softening of the pinch-point amplitude which loses its singular character and (ii) the generation of a finitecorrelation length ξ which endows the pinch points with a finite width ∼ 1/ξ [ Fig. 7]. The fact that the bowtie structure of susceptibility appears intact indicates that the low-temperature phase of the S = 1/2 nearestneighbor Heisenberg antiferromagnet respects the ice rules to a good degree of accuracy. An increase in temperature is seen to be associated with an overall decrease in the scattering intensity, while the bow-tie pattern appears to be remarkably robust up till T ∼ J 1 [Figs. 6 and 8], suggesting that the ice rules govern the physics over a surprisingly large temperature range. We find that, within a significant segment of this temperature range up till T ∼ J 1 , the width of the bow tie as measured by its full width at half maximum increases (approximately) linearly [ Fig. 8].
For the spin S = 1 nearest-neighbor Heisen-berg antiferromagnet with spatially isotropic couplings [Sec. III C 1], we find that, strikingly, the ground state remains magnetically disordered [ Fig. 13(a)] with no instability towards dimerizing into a valence-bond-crystal structure [ Fig. 9(b)], pointing to the realization of a rare scenario of a S = 1 quantum spin liquid in three dimensions. The formation of the bow-tie pattern of scattering now features relatively sharper pinch points, as seen by a decrease in their full width at half maximum compared to S = 1/2 [ Fig. 7]. This decrease is as expected, since with increasing spin, quantum fluctuations decrease in strength, and the ice rules are better fulfilled. We find that the bow-tie structure remains robust up till T ∼ J 1 , similar to what is observed for S = 1/2.
In the presence of breathing anisotropy (of arbitrary strength) in the nearest-neighbor Heisenberg antiferromagnet, we find that, for both S = 1/2 [Sec. III B 2] and S = 1 [Sec. III C 2], the quantum paramagnetic nature of the ground state remains intact [ Fig. 11]. The reciprocal space spin susceptibility profile is still characterized by bow ties and the associated "rounded" pinch points, whose width is found to remain essentially unchanged from the isotropic point down to the strongly anisotropic limit [ Fig. 12]. Our results thus point to the presence of an enlarged region in parameter space over which the low-temperature physics is approximately governed by the ice rules.
For the nearest-neighbor isotropic Heisenberg antiferromagnetic model with spin S > 1 [Sec. III D], we find that for S = 3/2 and beyond long-range dipolar magnetic order finally sets in [see Figs. 13 and 25]. We mention that, for the finite S values studied in our manuscript, the correct balance between leading 1/S terms and subleading contributions is already incorporated in the PFFRG [see Sec. II A]. However, with increasing S, the PFFRG becomes numerically more challenging (and also more sensitive to errors), because it becomes progressively difficult to account for the proper interplay between (large) leading 1/S and (much smaller but still important) subleading terms in our numerical algorithm. For this reason, we applied the PFFRG only to moderate spin magnitudes smaller than eight and use plain RPA in the infinite S limit [see Appendix A]. Therefore, we are unable to unambiguously address the question of the nature of the ground state (presence or absence of long-range magnetic order) in the large-S nearest-neighbor quantum Heisenberg antiferromagnet.
Upon inclusion of a J 2 coupling [Sec. IV], the complete parameter space of the J 1 -J 2 Heisenberg model is shown to host seven different kinds of magnetic orders in the classical model [ Fig. 15]. We have reported some corrections and/or amendments to previously known results [56] concerning the nature of the magnetic orders and the classical phase diagram [ Table I]. For low values of spin, i.e., S = 1/2 and S = 1, quantum fluctuations are shown to stabilize an extended domain of quantumspin-liquid behavior centered around the point J 1 > 0 and J 2 = 0, i.e., the nearest-neighbor Heisenberg antifer-romagnet [ Fig. 16]. For S = 1/2, the quantum spin liquid ranges from −0.25(3) J 2 /J 1 0.22 (3), while for S = 1, its span is reduced by half to −0.11(2) J 2 /J 1 0.09 (2) but remains nonetheless appreciable. The introduction of even a small J 2 coupling is seen to substantially modify the reciprocal space scattering profile at low temperatures such that the bow-tie structure becomes quickly obliviated accompanied by an enhancement (decrement) for antiferromagnetic (ferromagnetic) J 2 in the spectral weight at the wave vector (k = (0, 0, 4π)) where the pinch-point did exist [see Fig. 18]. Indeed, we find that for antiferromagnetic J 2 > 0 there is an enhancement of the pinch-point scattering amplitude as found in the corresponding classical model [71] [ Fig. 18 (first row) and Fig. 20], while for ferromagnetic J 2 the scattering intensity at the pinch points is strongly suppressed and instead redistributes to form a hexagonal cluster pattern of scattering [ Fig. 18 (second row) and Fig. 20] [71]. Interestingly, we do not observe the stabilization of a paramagnetic phase by frustrating the nearest-neighbor Heisenberg ferromagnet, i.e., in the regime J 1 < 0 (FM) and J 2 > 0 (AFM). The phase boundaries between magnetically ordered phases get significantly modified compared to the classical model [ Fig. 16], and the wave vectors of spiral orders get shifted by quantum fluctuations [Fig. 22]. Finally, we provide the Néel and Curie temperatures for different magnetically ordered phases, and for the S = 1/2 nearest-neighbor Heisenberg ferromagnet we benchmark our PFFRG results with available numerically exact quantum Monte Carlo and other methods [ Table II].

VI. OUTLOOK AND FUTURE DIRECTIONS
Our analysis of quantum effects on the pyrochlore lattice lays new avenues towards further exploration in search of novel quantum phases in a more generic symmetry-allowed Hamiltonian [152,153,[179][180][181][182] relevant for a large class of materials. Indeed, it has been shown at the classical level that anisotropic nearestneighbor spin interactions can stabilize novel phases such as spin liquids and spin nematics and a plethora of intricate magnetic orders [106,149,152,181,183,184]. The simplest extension to an XXZ model has been argued to serve as a minimal model of quantum spin ice [180] and has recently been shown to host spin nematic order and a variety of spin-liquid phases, albeit considered only at the classical level [106,184]. Surprisingly, little is known about the role of quantum fluctuations beyond a perturbative treatment [152,153,157,[185][186][187][188][189][190]. In particular, the nature of the competing ordered or disordered quantum phases in the low spin-S regime of the XXZ model remain open questions, and it will be interesting to investigate if, and to what extent, the quantum-spin-liquid phase of the isotropic model [31] found in this work remains stable in the presence of XXZ anisotropy.
Our identification of extended regimes of quantum spin liquid and, in general, quantum paramagnetic behavior in the S = 1/2 and S = 1 models in the presence of breathing anisotropy or J 2 coupling sets the stage for future theoretical and numerical studies aimed at identifying the precise nature of the quantum-spin-liquid phase, e.g., gapped or gapless spin liquid, and its associated gauge structure, SU(2), U(1), Z 2 , etc. One promising approach would be to carry out a fermionic projective symmetry group (PSG) classification [191][192][193] of the meanfield spin-liquid states on the pyrochlore lattice for both symmetric [30] and chiral spin liquids [194] similar to what has been accomplished on other lattices [195][196][197][198].
The ground-state energies of the corresponding projected variational wave functions could then be calculated from variational Monte Carlo methods [199,200], enabling one to identify the most competitive variational ground state, which could then be improved by a subsequent application of Lanczos steps to obtain an estimate of the true ground-state energy [64,66,92,201,202]. Recently, the PFFRG method has been successfully combined with a self-consistent Fock-like mean-field scheme to calculate low-energy effective theories for emergent spinon excitations in spin-1/2 quantum spin liquids [203]. In this approach, the two particle vertices, i.e., the effective spin interactions from PFFRG, are taken as an input for the Fock equation yielding a self-consistent scheme to determine spinon band structures beyond mean field. The precise forms of such free spinon Ansätze are dictated by a PSG classification of quantum spin liquids [191], allowing for a systematic investigation of kinetic spinon properties. It would be of interest and importance to apply this scheme to the pyrochlore Heisenberg antiferromagnet and compare the findings with those of variational Monte Carlo calculations. To address the issue of the nature of the elementary excitations and, in particular, to reveal the possible presence of a spinon continuum which is a manifestation of fractionalization and a hallmark of a quantum-spin-liquid phase, one needs a knowledge of the dynamical structure factor S(q, ω). The PFFRG framework can also be formulated directly in the real frequency domain employing the Keldysh formalism, which would allow one to obtain the complete S(q, ω). We leave the treatment of the Keldysh formalism and its application to the pyrochlore Heisenberg antiferromagnet as an important and exciting future endeavor.
From a materials perspective, a fascinating class of transition-metal-based fluorides with the pyrochlore structure have recently come into the limelight. This family of materials is at the boundary between quantum spin liquid, magnetic order, and magnetic freezing (or glassy regime). Their importance stems from the availability of large high-quality single crystals. Prominent candidate spin-liquid examples include the S = 1 NaCaNi 2 F 7 [204], which may be a first realization of a S = 1 quantum spin liquid in three dimensions [205], and the related higher-spin fluoride compounds featuring a high frustration index (f = Θ CW /T c ), such as NaCaCo 2 F 7 [206][207][208], NaCaFe 2 F 7 , NaSrFe 2 F 7 , and NaSrMn 2 F 7 [209], which, nonetheless, either show signs of long-range magnetic order at low temperatures or undergo spin freezing [210]. With the PFFRG formalism in place, it would be useful in such a material context to extend the mapping of the quantum phase diagram in the presence of longer-ranged Heisenberg couplings which will most likely give rise to additional novel phases compared to the seven phases of the classical J 1 -J 2 Heisenberg model, as, for instance, shown in Ref. [57] for classical spins. It would seem likely that most of the abovementioned materials could be placed to a good degree of approximation in the extended phase diagram so determined.
Given that frustrated quantum spin systems are challenging to deal with theoretically and, in three dimensions, pose a formidable barrier to most quantum manybody numerical methods, PFFRG is one of the very few methods that can be used to shed light on the physics at play in these systems, with the field now poised to benefit from the arrival of more materials. It is in this broader context that we investigated and presented in this paper the rich example of the J 1 -J 2 Heisenberg model on the pyrochlore lattice.  In this Appendix, we present further details about how the susceptibility of the nearest-neighbor pyrochlore Heisenberg model in the infinite-S limit as depicted in Fig. 4 is computed. Particularly, we explain why the simple RPA-type summation which we use to obtain these results reproduces the correct pinch-point singularity but also results in a spurious divergence of the susceptibility at a finite temperature T which is not expected from the exact (numerical) solution [3,8]. We further present analytical arguments why the summation of further diagram classes can cure this artifact by regularizing the divergence.
We begin by reviewing the PFFRG scheme in the large-S limit and explain that, to leading order when S → ∞, the PFFRG becomes identical to a simple RPA-type approximation (for further details, see Ref. [74]). As briefly mentioned in Sec. II A 1, the generalization of the PF-FRG for arbitrary spin S amounts to introducing fermion flavors f i↑κ , f i↓κ with κ = {1, . . . , 2S} on each lattice site i which add up to a total spin S. Furthermore, to avoid diverging energy scales in the large-S limit, it is convenient to renormalize all interactions via J ij → J ij /(2S). As a consequence of the additional flavor index κ, the Feynman diagrams acquire an extra factor of 2S for each closed fermion loop. Hence, when formulating the PF-FRG equations for arbitrary S, the second term on the right-hand side in Fig. 2(b) (the so-called RPA channel) acquires a prefactor of 2S, indicating that, among all interaction channels in Fig. 2(b), this term is singled out at large S. The flow equation for the two-particle vertex at S → ∞, where only the RPA term contributes on the right-hand side, can be readily solved [74] and leads to the RPA-type diagram series shown in Fig. 23(a). These two-particle vertex diagrams are precisely the ones, and no others, of leading order in 1/S. This result is evident from the fact that, for a given number of interaction lines, they each maximize the number of loops. Specifically, each term of the series has n bare interaction lines and n − 1 fermion loops, resulting in an overall order of 1/S. Having established that, to leading order in 1/S, the PFFRG generically reduces to an RPA-type approximation, we now study the structure of this approximation in the context of the nearest-neighbor pyrochlore Heisenberg model. In the following, we are interested only in the static frequency components (ω = 0) of the two-particle vertex Γ Λ (1 , 2 ; 1, 2). We thus omit the arguments 1 , 2 , . . . and write Γ Λ (1 , 2 ; 1, 2) → Γ Λ ij with the site indices i, j as subscripts. Furthermore, the propagators considered are the bare (i.e., without self-energy corrections) and Λ-regularized ones from Eq. (6). The RPA diagram series may be expressed in a self-consistent form [second FIG. 23. RPA-type approximations for the two-particle vertex in the large-S limit. Dashed lines are the bare interactions Jij, and lines with an arrow are the bare and Λ-regularized pseudofermion propagators from Eq. (6). Gray boxes denote the two-particle vertex in different approximations. (a) Plain RPA scheme summing up diagrammatic terms of the order of 1/S. (b) An example of a contribution to the two-particle vertex of the order of (1/S) 2 . (c) Improved RPA scheme, RPA', regularizing the divergence of the two-particle vertex occurring in plain RPA. See the text for details. line of Fig. 23(a)], leading to Here, Π Λ is the ω = 0 component of the bare fermion loop and is given by Π Λ = S/(πΛ). The solution of Eq. (A1) can be obtained via a Fourier transform, giving whereJ(k) is the interaction matrix in sublattice space as given by Eq. (19).Γ Λ (k) is also analogously defined in sublattice space, and 1 denotes the identity matrix in the same space. To better understand the physical implications of Eq. (A2), we diagonalizeJ(k) via is a diagonal matrix whose elements are the eigenvalues ofJ(k). It follows that For the nearest-neighbor pyrochlore Heisenberg antiferromagnetic model, the lowest bands ofJ d (k) take the form of two degenerate flat modes with an energy −2J 1 [2,113]. As a result of these flat modes, the matrix Π Λ 1 + 2SJ −1 d (k) in Eq. (A3) becomes singular at Λ = J 1 /π for all wave vectors k, which leads to a diverging susceptibility at the corresponding (finite) temperature T = 2π 3 S(S + 1)Λ. However, as explained further below, this divergence is a methodological artifact of the plain RPA treatment within which only the leading 1/S diagrammatic contributions are considered.
The flat modes inJ d (k) are also responsible for the pinch-point singularities in the susceptibility [140]. To see this, we first note that (up to irrelevant overall factors from fusing external fermion lines) the susceptibility χ Λ (k) of Eqs. (10) and (11), rewritten in sublattice coordinates, is related to the two-particle vertexΓ Λ (k) via Here, α, β are sublattice indices and ξ α denote the sublattice displacements, i.e., site coordinates r i , unit cell coordinates R i , and displacements ξ α fulfilling r i = R i + ξ α . Since the lowest (flat) modes give the dominant contribution to the susceptibility and also describe the physics of pinch points we are interested in, we may approximate Eq. (A3) by neglecting higher-energy bands inJ d (k). Using Eqs. (A3) and (A4), one then obtains where γ = fm only sums over the flat modes (fm). The numerator in this expression (which is used to plot the inset in Fig. 4) contains the pinch-point pattern, while the denominator produces the aforementioned singularity at finite Λ. This analysis shows that in plain RPA, as obtained from PFFRG in leading order in 1/S, the pinch points are correctly reproduced. However, their manifestation within this plain RPA scheme is implicitly tied with a divergence of the k-dependent susceptibility for all k that define the flat modes. Thus, the physically correct paramagnetic (broadened) pinch points observed in plain RPA exist only above the instability, and so their discussion in plain RPA is bounded from below by the instability at Λ = J 1 /π. We now investigate how Eq. (A5) is modified when adding diagrams of order higher than 1/S. Within PFFRG, such higher orders are generally described by the other interaction channels on the right-hand side in Fig. 2(b), i.e., those corrections to RPA which do not contain a fermion loop. In contrast to the leading order in 1/S discussed above, where all diagrammatic contributions to the two-particle vertex are exactly included in the PFFRG, higher orders are treated only approximately. A thorough analytical discussion of all subleading diagrams implicitly included within the PFFRG computational scheme is, admittedly, very challenging, because, already to the order of (1/S) 2 , they may not be represented by a simple series of diagrams such as the one shown in Fig. 23(a). Furthermore, from a more technical perspective, it is a rather involved task to apply the PFFRG at large but finite S and systematically explore the effects of different diagrammatic orders in 1/S. This hurdle arises because of numerical difficulties in capturing the subtle competition between large leading 1/S and much smaller, but still important and possibly singular subleading terms, when the frequency dependence of the vertex functions is approximated by a finite grid (which is a computational necessity within PFFRG).
To still be able to investigate general properties of higher diagrammatic orders in 1/S, we, therefore, use a different strategy. We take as a starting point the S → ∞ limit (as described above) and then incorporate "by hand" subleading diagrams to study their effects on the spurious divergence encountered in a plain RPA treatment. Subleading diagrams of the order of (1/S) 2 are obtained by feeding back the RPA two-particle vertex into a fermion loop of the RPA series as shown in Fig. 23(b). In the following, we discuss a generalization of such terms (dubbed RPA ) where (i) the feedback of the RPA takes place in every fermion loop and (ii) the insertion is performed self-consistently as shown in Fig. 23(c). The resummation of such diagram classes also involves contributions from orders higher than (1/S) 2 . This type of approximation first amounts to replacing the bare fermion loop Π Λ by Π Λ + Π Λ , where Π Λ is the loop diagram with the RPA series reinserted as depicted in Fig. 23(c). Using the fact that only the local two-particle vertex Γ Λ ii contributes to this diagram, one finds Without the loss of generality, we choose the "11"sublattice component of the two-particle vertex, since all sublattices are equivalent in the paramagnetic regime. Also note that, in order for the calculation to be analytically tractable, we perform a static approximation where the two-particle vertex is assumed to be ω independent. The self-consistency for Π Λ is closed using Eq. (A3) and replacing Π Λ → Π Λ + Π Λ , giving Γ Λ 11 (k) = −M (k) Π Λ 1 + Π Λ 1 + 2SJ −1 d (k) −1 M † (k) 11 .
(A7) Here again, we consider only the contribution from the flat modes in J d (k) and neglect higher-energy bands. Furthermore, we write the momentum integral (which is a positive dimensionless number) as and again use Π Λ = S/(πΛ), leading to This is a quadratic equation for Π Λ which can be solved to yield the susceptibility where, when compared to Eq. (A5), an additional contribution from Π Λ appears in the denominator. From the two solutions for Π Λ following from Eq. (A9), the correct one is identified by the condition that the leading order at large S must be a contribution ∼ 1/S as is the case for the bare RPA. One then obtains Most importantly, this expression no longer has a divergence in Λ while the pinch-point pattern given by the k-dependent term [first line of Eq. (A11) and numerator of Eq. (A5) which generate the pinch points] persists. The Λ-dependent second and third line of Eq. (A11) is plotted in Fig. 24 for S = 1000 and S = 10 000. It can be seen that the diverging susceptibility of the RPA scheme is regularized by the higher-order terms such that χ Λ (k) becomes bounded in the vicinity of the singularity. Yet, certain artifacts still remain in the RPA scheme such as a steplike behavior of the susceptibility and a finite interval where χ Λ (k) becomes imaginary (the size of this interval shrinks with increasing S). We expect that such spurious behavior would become further regularized upon including more diagrammatic contributions.
In summary, even though this analysis is based on an approximate resummation of a certain class of diagrams, it demonstrates that higher-order terms have a significant effect even in the large-S limit and may counteract the diverging susceptibility observed in the bare RPA calculation leading to Eq. (A5). This calculation also shows that-even though counterintuitive at first sight-leading 1/S diagrams are not sufficient to treat the classical limit S → ∞ exactly. One may, therefore, conclude that, while the spatial structure of the spin correlations at large S is already correctly described by plain RPA, thermal fluctuations are much more intricate in pseudofermionic formulation. This conclusion may possibly indicate that pseudofermions are not ideally suited to describe the thermodynamics of spin systems in the classical large-S limit. We also emphasize, however, that such methodological subtleties do not affect the PFFRG at finite (but not too large) S, where the correct balance between classical magnetic phenomena and quantum fluctuations is captured by the interplay between leading 1/S and leading 1/N diagrammatic contributions [where N generalizes the spin symmetry group from SU(2) to SU(N ); see Sec. II A 1 for details]. Here, we present the details of the numerical procedure [98] used to detect the onset of long-range magnetic order in the RG flow. The expected divergence of the spin susceptibility [Eq. (11)] at a critical Λ which would signal the spontaneous breaking of SU(2) spin-rotation symmetry towards long-range dipolar magnetic order is, in practice, regularized due to two numerical approximations in the PFFRG method: (i) the discretization of the frequencies in the arguments of the vertex functions and (ii) the finite spatial extent of the two-particle vertex function. Both these approximations regularize the divergence to a finite maximum, or a feeble kinklike feature when the ordered magnetic moment is small. In addition, the discretization of the frequencies induces the artifact of oscillations in the susceptibility flow, especially at small Λ. The distinct advantage of the method presented here lies in its ability to detect such kinks even in the presence of pronounced frequency oscillations and a small ordered magnetic moment. To illustrate the method, we focus on the transition with increasing spin S, from the paramagnetic into the magnetically ordered phase, for the nearest-neighbor pyrochlore Heisenberg antiferromagnet.
The appearance of a finite maxima or a kinklike feature in the susceptibility evolution with decreasing Λ is marked by a change in the slope of the RG flow. However, as the susceptibility flow is plagued by oscillations due to frequency discretization, one encounters a difficulty in defining the slope. As each discrete frequency grid point produces a small peak or an upturn in the susceptibility flow, we compute the slope in a manner that averages out these oscillations. To this effect, one connects two adjacent peaks via a straight line which represents a tangent of the susceptibility and approximates χ(k) between the two peaks. A kink in the RG flow now manifests as a change in the slope, i.e., a finite-angle α, between the two neighboring tangents, as shown in Fig. 25(a), which then serves as a measure of the size of the kink. We One divides a fixed Λ interval into two regions, I and II, which by construction lie between two adjacent pairs of susceptibility kinks or peaks. Within each region, χ(k) is approximated by a tangent connecting the neighboring kinks. We find that, while the angle between the two tangents is negligible for S = 1, it acquires a sizable finite value for S = 3, implying that two curves represent different phases. To obtain a more quantitatively robust measure for the size of the kink, we consider additional pairs of adjacent peaks, compute the angles between their tangents, and calculate the average angleᾱ. (b) Averaged angleᾱ as a function of spin S. Whileᾱ is negligible and almost constant for S = 1/2 and S = 1 (labeled by filled circles), there is a pronounced increase for higher values of spin S 3/2 (empty circles). Based on this behavior, we estimate the phase transition of the spin-S nearest-neighbor Heisenberg antiferromagnet to occur for S = 3/2. first choose a fixed Λ interval wherein multiple kinks, potentially representing magnetic instabilities, appear to be located. We then consider tangents between different pairs of adjacent peaks and calculate the averageᾱ of the absolute value of these angles within a given Λ interval. The angleᾱ then serves as a relatively robust quantitative measure of the change in slope (i.e., the size of the kink) within this Λ interval; a largerᾱ implying a more pronounced kink. To locate the phase transition, we plotᾱ as a function of the spin S [see Fig. 25(b)]. For S = 1/2 and S = 1, we observe a small and constant value ofᾱ ≈ 1.6°, followed by a sudden increase at S = 3/2 indicating a transition point to magnetic longrange ordered state.