Quantum and Classical Phases of the Pyrochlore Heisenberg Model with Competing Interactions
Yasir Iqbal Tobias Müller Pratyay Ghosh Michel J. P. Gingras Harald O. Jeschke Stephan Rachel Johannes Reuther Ronny Thomale





yiqbal@physics.iitm.ac.in
Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India

Institute for Theoretical Physics and Astrophysics, Julius-Maximilian's
University of Würzburg, Am Hubland, D-97074 Würzburg, Germany

Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India

Perimeter Institute for Theoretical Physics, Waterloo, Ontario, Canada N2L 5G7
Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1
Quantum Materials Program, Canadian Institute for Advanced Research, MaRS Centre, West Tower 661 University Avenue, Suite 505, Toronto, Ontario, M5G 1M1, Canada

Research Institute for Interdisciplinary Science, Okayama University, 3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530, Japan

School of Physics, The University of Melbourne, Parkville, Victoria 3010, Australia
Institut für Theoretische Physik, Technische Universität Dresden, D-01069 Dresden, Germany

Dahlem Center for Complex Quantum Systems and Fachbereich Physik, Freie Universität Berlin, D-14195 Berlin, Germany
Helmholtz-Zentrum Berlin für Materialien und Energie, D-14109 Berlin, Germany

Institute for Theoretical Physics and Astrophysics, Julius-Maximilian's
University of Würzburg, Am Hubland, D-97074 Würzburg, Germany




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. 


Quantum and Classical Phases of the Pyrochlore Heisenberg Model with Competing Interactions
    Ronny Thomale
    December 30, 2023
===========================================================================================






§ 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” <cit.>. This failure is a consequence of the extensive classical ground-state degeneracy <cit.> which proves severe enough to prevent a thermal “order-by-disorder” mechanism <cit.> from selecting a unique ground-state ordering pattern <cit.>. 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 <cit.>, it is known that at harmonic order in 1/S, the extensive classical ground-state degeneracy exp[𝒪(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[𝒪(L)] <cit.>. It turns out that the consideration of higher-order  terms in a  1/S expansion also fails to select a unique ground state <cit.>. Indeed, while quartic corrections in boson operators do break the degeneracy of the harmonic ground states, there still remains a family of (almost) degenerate (exp[𝒪(L)]) states <cit.>. 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 <cit.> or a quantum-spin-liquid <cit.> 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 <cit.>. In the much-less-investigated case of S=1 <cit.>, there have been suggestions of a ground state with tetrahedral symmetry breaking <cit.>.  

The “cooperative paramagnet” ground state of the classical nearest-neighbor Heisenberg antiferromagnet is known to be extremely fragile, in that magnetic long-range order is induced upon the inclusion of various perturbations, such as further neighbor Heisenberg interactions <cit.>, dipole interactions <cit.>, Dzyaloshinsky-Moriya anisotropy <cit.>, single-ion anisotropy <cit.>, lattice distortions <cit.>, and bond disorder <cit.>. In particular, further neighbor Heisenberg interactions are found to stabilize a plethora of intricate classical magnetic orders <cit.>. However, in the low-spin-S regime, where the strong possibility of a quantum paramagnetic 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:




    Ĥ=J_1 ∑_⟨ i,j ⟩Ŝ_i·Ŝ_j
    +J_2 ∑_⟨⟨ i,j ⟩⟩Ŝ_i·Ŝ_j,

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. <ref>]. 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 <cit.>, while successful in one and two dimensions, become unfeasible in three dimensions due to entanglement scaling and system-size limitations. Quantum Monte Carlo methods <cit.>, while able to reach sufficiently large system sizes, are, in principle, restricted to unfrustrated systems, while variational Monte Carlo approaches <cit.>, which are shown to be extremely successful in two dimensions <cit.>, 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 <cit.>. Thus, one is essentially left with only mean-field approaches based on Schwinger bosons <cit.>, semiclassical analysis based on spin waves, or linked-cluster expansion methods <cit.>, 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 <cit.>. 

By employing PFFRG for the spin-S J_1-J_2 Heisenberg model, we find for S=1/2 an extended quantum-spin-liquid 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 <cit.>. 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 <cit.>. 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 <cit.>, 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-bond-crystal states.
  
The paper is organized as follows: In Sec. <ref>, we describe the PFFRG method (Sec. <ref>) employed for the quantum treatment of the model, starting with a description of its formalism (Sec. <ref>) followed by details of its numerical implementation in Sec. <ref>. In Secs. <ref> and <ref>, we discuss schemes used to obtain the ground state of classical spin models, namely, the Luttinger-Tisza method [Sec. <ref>] and the iterative minimization of the energy [Sec. <ref>] (the reader interested mainly in the results can directly jump to Secs. <ref> and <ref>). Employing these methods, we begin with a treatment of the ground-state and low-energy physics of the nearest-neighbor Heisenberg antiferromagnet in Sec. <ref>, starting first with a classical analysis [Sec. <ref>] of the isotropic and breathing lattices and then moving on to the quantum treatment of the S=1/2 [Sec. <ref>] and S=1 [Sec. <ref>] 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. <ref>]. Next, in Sec. <ref>, we deal with the J_1-J_2 Heisenberg model, by first revisiting the classical phase diagram [Sec. <ref>], and subsequently present the results for the quantum model in Sec. <ref>. 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. <ref>, followed by an outlook and discussion of future directions in Sec. <ref>.
  


§ METHODS





 §.§ Pseudofermion functional renormalization group method




  §.§.§ Formalism




The key idea of the PFFRG method <cit.> is to express the spin-1/2 operators in terms of pseudofermions <cit.>,

    Ŝ^μ_i=1/2∑_α,βf̂^†_iασ^μ_αβf̂_iβ,

where σ^μ_αβ are Pauli matrices (μ∈{x,y,z}) and f̂_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. <cit.>, 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

    Ŝ_i→∑_κ=1^M Ŝ_iκ ,

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. (<ref>) amounts to equipping the fermion operators with an additional index κ:

    Ŝ_iκ^μ=1/2∑_αβf̂_iακ^†σ_αβ^μf̂_iβκ.

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 (∑_κ=1^MŜ_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 [ Heisenberg systems (on any lattice) with S=1 and single-ion anisotropies Δ∑_i (S_i^z)^2 provide a simple exception wherein if Δ is positive (and
sufficiently large) this term would always energetically prefer the
unphysical spin sector S=0 over all other sectors.]). 

Rewriting the spin Hamiltonian in terms of Eq. (<ref>), the resulting fermionic model is treated within the standard FRG framework for interacting fermion systems <cit.>. 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 [The class of diagrams representing the leading order in 1/S contributions, which are thus of random phase approximation (RPA)-type and responsible for the formation of classical magnetic order, are summed up exactly. Similarly, the class of diagrams capturing contributions to leading order in 1/N, and thus responsible for the formation of nonmagnetic states, are also summed up exactly. However, an accurate treatment of the S→∞ limit may require a consideration of subleading terms in 1/S, thus going beyond a bare RPA treatment [see Appendix <ref>]]. As a consequence, classical magnetically ordered states (typically favored at large S) and nonmagnetic spin liquids or dimerized states (as obtained at large N) <cit.> 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

    G_0(iω)=1/iω ,

where iω denotes a frequency on the imaginary Matsubara axis. Within the standard PFFRG scheme <cit.>, this propagator is dressed with an infrared steplike regulator function:

    G_0(iω)⟶ G_0^Λ(iω)=Θ (  | ω | - Λ) /iω,

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. <ref>, where the arrows denote dressed and Λ-dependent propagators

    G^Λ(iω)=Θ (  | ω | - Λ)/iω-Σ^Λ(iω)

and slashed lines denote the single-scale propagator

    S^Λ(iω)=δ (  | ω | - Λ)/iω-Σ^Λ(iω) .

Because of the locality of fermion propagators, the two-particle vertex Γ^Λ(1',2';1,2) effectively depends on two site indices only, i.e., Γ^Λ(1',2';1,2)∼δ_i_1 i_1'δ_i_2 i_2'. As illustrated in Fig. <ref>, 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. <ref> 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 three-particle vertex in Fig. <ref>(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 <cit.>, 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 <cit.>. Still, parts of the three-particle vertex can be included by applying the so-called Katanin truncation <cit.>, which replaces the single scale propagator by the full Λ derivative of the dressed propagator

    S^Λ⟶ -d/dΛG^Λ=S^Λ-(G^Λ)^2d/dΛΣ^Λ .

While the additional Katanin terms formally have the structure of the three-particle term [the last term in Fig. <ref>(b)], they should rather be understood as self-energy corrections <cit.>. 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. (<ref>)] 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. <ref>(a)] and in the RPA contribution for the two-particle vertex [the second term on the right-hand side of Fig. <ref>(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 <cit.>. 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 <ref>. 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. <ref>(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 <cit.>; however, subsequent refinements have made it capable of handling a wide spectrum of frustrated magnetic Hamiltonians for multilayer systems and in three dimensions <cit.>.



  §.§.§ 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

    χ_ij^zz=∫_0^∞ dτ < T_τS^z_i(τ)S^z_j(0)  > ,

where T_τ (with τ being the imaginary time) is the imaginary time-ordering operator. 

Transforming χ_ij^zz into 𝐤 space yields the wave-vector-resolved susceptibility χ(𝐤):


    χ(𝐤)=1/4∑_i=1^4∑_jχ_ij^zze^i𝐤·(𝐫_i-𝐫_j),

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 χ(𝐤). 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 𝐤 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 χ(𝐤) 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 <ref> 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 T_c/J=(2π S(S+1)/3)Λ_c/J <cit.>. 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 <cit.>, i.e., a mean-field description, and the conventional spin mean-field theory which is formulated in terms of the temperature instead of Λ <cit.>. 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 χ(𝐤) 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 <cit.> and valence-bond-crystal orders <cit.> 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 order-parameter susceptibility, which is given by a four-spin correlator. For spin-nematic order, this correlator is the standard nematic correlation function ∑_μ,ν⟨𝒪_ij^μν𝒪_kl^μν⟩, where 𝒪_ij^μν=Ŝ_i^μŜ_j^ν-(δ_μν/3)Ŝ_i·Ŝ_j <cit.> (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 correlation function D_ijkl=⟨(Ŝ_i·Ŝ_j)(Ŝ_k·Ŝ_l)⟩-⟨Ŝ_i·Ŝ_j⟩^2. In PFFRG, 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. (<ref>)] 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 support 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 is

    ℋ̂_SN=δ∑_⟨ ij⟩(Ŝ_i^xŜ_j^x+Ŝ_i^yŜ_j^y)-δ∑_⟨ ij⟩Ŝ_i^zŜ_j^z

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 is 

    Ĥ_ VBC=δ∑_⟨ i,j⟩∈SŜ_i·Ŝ_j- 
                       δ∑_⟨ i,j⟩∈WŜ_i·Ŝ_j,

which strengthens the couplings J_ij on all bonds in S [J_ij→  J_ij+δ for ⟨ i,j⟩∈ S] and weakens 
the couplings in W [J_ij→  J_ij-δ for ⟨ i,j⟩∈ W]. 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 nearest-neighbor spin susceptibilities χ_ij. We then define the nematic response function for a given pair of nearest-neighbor sites by 

    η_SN=J/δ(χ_ij^xx)_Λ - (χ_ij^zz)_Λ/(χ_ij^xx)_Λ + (χ_ij^zz)_Λ ,

where χ_ij^xx (χ_ij^zz) are the correlators on the strengthened (weakened) bonds. Similarly, the dimer response function is given by

    η^P_ VBC=J/δ(χ_S_P)_Λ-(χ_W_P)_Λ/(χ_S_P)_Λ+(χ_W_P)_Λ,

where, χ_S_p (χ_W_p) 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. 



 §.§ 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 →∞ <cit.>. This procedure yields the corresponding classical spin system wherein the spin operators in Eq. (<ref>) are replaced by ordinary vectors of unit length at each lattice site i. For general interactions, the classical Hamiltonian to be minimized reads as

    H = ∑_i,j, α, β J_αβ(𝐑_ij) 𝐒_i, α·𝐒_j, β,

where by i/j we denote the primitive lattice site separated by the lattice translation vectors 𝐑_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 <cit.> attempts to find a ground state of Eq. (<ref>) by enforcing the spin-length constraint only globally, ∑_i|𝐒_i^2|=S^2n, where n is the total number of lattice sites, which is termed the weak constraint. This relaxed constraint implies that site-dependent average local moments are now permissible, which, strictly speaking, take us beyond the classical limit by approximately incorporating some aspects of quantum fluctuations <cit.>.

To solve this relaxed problem, we decompose the spin configuration into its Fourier modes 𝐒̃_α(𝐤) on the four sublattices of the pyrochlore lattice

    𝐒_i, α = 1/√(N/4)∑_𝐤𝐒̃_α(𝐤) e^𝐤·𝐫_i,α.

Inserting this equation into Eq. (<ref>) results in

    H = ∑_𝐤∑_α, βJ̃_αβ(𝐤) 𝐒̃_α(𝐤)·𝐒̃_β(-𝐤),

with the interaction matrix given by

    J̃_αβ(𝐤) = ∑_i,j J_αβ(𝐑_ij) e^𝐤·𝐑_ij.


The optimal modes satisfying the weak constraint are then given by the wave vector 𝐤, for which the lowest eigenvalue of Eq. (<ref>) has its minimum. The eigenvector corresponding to this eigenvalue gives the relative weight of the modes on the sublattices <cit.>, which means that the optimal modes do not fulfill the strong constraint (|𝐒_i^2|=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 <cit.>. 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. <ref>.



 §.§ 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 <cit.>. 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

    𝐡_i = ∂ H/∂𝐒_i = ∑_j J_ij𝐒_j.

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 16 L^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 

    ℱ(𝐤) = 1/16 L^3|∑_i 𝐒_i e^𝐤·𝐫_i|^2

is computed, which is, up to a normalization constant, the same as the susceptibility defined in Eq. (<ref>), 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.
 
 
 


§ 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 anisotropy

    Ĥ=J_ up∑_⟨ i,j ⟩_ upŜ_i·Ŝ_j + J_ down∑_⟨ i,j ⟩_ downŜ_i·Ŝ_j,

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 scale J̃:  


    J_ up=J̃cos(φ),     J_ down=J̃sin(φ).


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 <cit.>. 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_4O_8, LiInCr_4O_8, LiInCr_4S_8, LiGaCr_4S_8, CuInCr_4S_8, and CuInCr_4Se_8 <cit.> and in a pseudospin S=1/2 Yb-based compound Ba_3Yb_2Zn_5O_11 <cit.>. 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 <cit.> 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.



 §.§ Classical model



  §.§.§ Isotropic case


At the isotropic point of Eq. (<ref>), we have J_ up=J_ down=J̃/√(2)≡ J_1. Henceforth, all temperatures for the isotropic classical and quantum models are expressed in units of J_1S(S+1) and  J_1, respectively (and we omit the factor of √(2)), while for the breathing model they are expressed in units of J̃S(S+1) and J̃ for the classical and quantum models, respectively. In the classical limit of Eq. (<ref>), the Heisenberg spin operators Ŝ_i reduce to standard three-component vectors 𝐒_i. In the ensuing analysis, it proves convenient to introduce the magnetization 𝐌_𝒯 of the 𝒯th tetrahedron,

    𝐌_𝒯=∑_α=1^4𝐒_𝒯,α,
 
where the index α=1,2,3, and 4 labels the four spins within the 𝒯th tetrahedron. In terms of 𝐌_𝒯, the Heisenberg Hamiltonian can be recast as a disjoint sum of the square of the magnetizations 𝐌_𝒯 over the “up” and “down” tetrahedra,

    ℋ_ isotropic=J_1/2∑_𝒯𝐌_𝒯^2 -  const.
 



From Eq. (<ref>), it follows that any state which satisfies the condition 𝐌_𝒯=0 on each tetrahedron 𝒯 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” <cit.>, 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 𝐌_𝒯 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 <cit.> that the corrections to the estimate for D are at most subextensive. The extensive (exp[𝒪(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” <cit.> with nonzero entropy <cit.>, referred to as a “classical spin liquid” <cit.>. 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 <cit.> of the type characterized by an order parameter which takes on its maximal value in a perfectly collinear state <cit.>. At T=0, the classical spin liquid features critical, i.e., algebraic, spin-spin correlations of dipolar character <cit.>, which is a consequence of the local constraint that the magnetization 𝐌_𝒯 on each tetrahedron is identically zero for any ground state <cit.>. 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. <ref>(a)] with sharp singularities termed pinch points [see the  encircled point in Fig. <ref>(a)] <cit.>. 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 <cit.>. This feature 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 <cit.> as well as those with local sublattice-dependent easy planes which are perpendicular to the local ⟨111⟩ axes <cit.>. The limit N=1 (Ising spins) is realized in various spin-ice materials A_2B_2O_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 <cit.>. Coming back to the case of N=3 (Heisenberg spins) at finite temperatures, we note that thermal fluctuations lead to violations of the 𝐌_𝒯=0 constraint and generate a finite correlation length ξ which, at low temperatures, diverges as T^-1/2 <cit.>. 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/ξ <cit.> [see Figs. <ref>(a)–<ref>(c)] which, at low temperatures, goes to zero as T^1/2 <cit.> [see Fig. <ref>].



  §.§.§ Breathing case



As in the case of the isotropic pyrochlore lattice, the Heisenberg Hamiltonian in the presence of breathing anisotropy [Eq. (<ref>)] can be straightforwardly recast as a disjoint sum of terms, each involving the magnetization 𝐌_𝒯 [Eq. (<ref>)] of a tetrahedron 𝒯:


    ℋ_ breathing=J_ up∑_𝒯∈ up𝐌_𝒯^2+
    J_ down∑_𝒯∈ down𝐌_𝒯^2-  const.
 

It is clear that when J_ up and J_ down are both antiferromagnetic, any state in which 𝐌_𝒯=0 on every up and down tetrahedron 𝒯 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 <cit.>. 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_ upJ_ down/J̃^2 =cosϕsinϕ <cit.>, and, hence, the development of the correlations is slower when closer to the decoupled tetrahedron limit. In Fig. <ref>, 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. <cit.>. 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.







 §.§ Spin-1/2 model




  §.§.§ Isotropic case


The investigation of the low-temperature (T≪ J_1) physics of Eq. (<ref>) 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 <cit.>. 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 𝐤, we track the evolution of the susceptibility χ(𝐤) with Λ for all wave vectors 𝐤 in the extended Brillouin zone (EBZ) of the pyrochlore lattice. As discussed in Sec. <ref>, the onset of magnetic long-range order at a particular 𝐤 is signaled by the presence of kinks or cusps in the Λ flow of χ(𝐤), whereas a smooth monotonically increasing behavior of χ(𝐤) down to Λ→ 0 points to a quantum-disordered ground state. For S=1/2, we observe that the Λ evolution of the susceptibility χ(𝐤) ∀ 𝐤 ∈  EBZ [see Fig. <ref>(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 <ref>]. A numerical maximization of the susceptibility function in the EBZ finds feeble maxima at the high-symmetry W points, i.e., at 𝐤=2π(2,1,0) [see Fig. <ref>(a)]. The RG flow of the susceptibility evaluated at the W point is shown in Fig. <ref>(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 <cit.>. 



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. <ref>(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 χ(𝐤) in the [hhl] plane (i.e., k_x=k_y plane) [see Fig. <ref>(d)], wherein one clearly sees the characteristic bow-tie pattern, albeit with a softening and broadening of the pinch points due to quantum fluctuations <cit.>. 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. <ref>, we show the variation of χ(𝐤) along the width of the pinch point, i.e., along the white vertical line in Fig. <ref>(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. <cit.>, 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. <cit.>. The temperature evolution of the susceptibility in the [hhl] plane is shown in Fig. <ref>. To obtain a quantitative picture, we plot in Fig. <ref> the susceptibility along a 1D cut (white line in the T/J_1=1/100 plot of Fig. <ref>) across the width of the pinch point in the bow-tie 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. <ref>) in contrast to the T^1/2 behavior expected classically (see Fig. <ref>). 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. <cit.>. 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. <cit.>. Also, the signs of all correlators up to the 16th neighbor as obtained from PFFRG agree with those obtained in Table I of Ref. <cit.>. This agreement is interesting in light of the fact that Ref. <cit.> 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 <cit.>. 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_3 D), (ii) two tetrahedral axis directions (VBC_2 D), and (iii) one tetrahedral axis direction (VBC_1 D). The dimer response functions η_ VBC^P [Eq. (<ref>)] of all three VBCs are found to decrease under the RG flow [see Fig. <ref>(a) for the RG flow of η_ VBC^P] 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 <cit.>, 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. <cit.>, wherein it is found that the system evades such nematic order <cit.>. Here, we investigate for the possibility of nematic order [see Sec. <ref>] in the S=1/2 nearest-neighbor Heisenberg antiferromagnetic model. We plot the RG flow of the nematic response function η_ SN [Eq. (<ref>)] in Fig. <ref>, 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. <cit.>, which argues for a nematic quantum spin liquid featuring spin-nematic order in the S=1/2 nearest-neighbor Heisenberg antiferromagnetic model, we mention that, since we a priori exclude the fermionic four-particle 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 inter- to 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. <ref>(a), wherein we observe a smooth flow down to Λ→0, in similarity with the finding for the isotropic model [see Fig. <ref>(b)]. Our results thus point 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 quantum-spin-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 <cit.> while the other finds lattice nematic order <cit.>. 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. <ref>(a)]. This result  shows that in the quantum paramagnetic ground state the low-energy physics is approximately governed by the ice rules.  





 §.§ Spin-1 model




  §.§.§ 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 <cit.>. The Λ evolution of the susceptibility at the 𝐤 vector where it has its maximum value, i.e., the high-symmetry W point, is shown in Fig. <ref>(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 <ref> 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 𝐤∈ 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. <ref>].

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. <ref>. The Λ evolution of the dimer response functions for the three VBCs [see Fig. <ref>(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. <ref>(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 <ref>] 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 <cit.>. We also observe that the bow-tie pattern and the pinch-point width remain essentially unchanged compared to the isotropic model [see Fig. <ref>(b)], indicating that the ice rules continue to dictate the low-energy physics of the quantum paramagnetic ground state even for strong breathing anisotropy. 



 §.§ 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. <ref>(b)] shows feeble signatures of the development of an instability or kink at the point marked by an arrow. This faint feature, appearing in the S=3/2 RG flow, develops into a pronounced kink (marking the breakdown of the RG flow) for increasing values of S [see Figs. <ref>(c) and <ref>(d)]. The details of the scheme employed to detect the instability or kink are given in Appendix <ref>. Based on this analysis [see Fig. <ref>], 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. <ref>]. 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. <ref>]. 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 <ref>]. 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 <cit.>. 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 𝐤∈ EBZ out of the extensively degenerate classical ground-state manifold. In Fig. <ref>, we show the variation in the susceptibility along a path passing through the high-symmetry points [see Fig. <ref>(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. <ref>(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 <cit.> 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 <cit.>. 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 <cit.>.





§ J_1-J_2 HEISENBERG MODEL



 §.§ Classical phase diagram



[table]capposition=bottom



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 <ref> and Fig. <ref>], part of which is investigated in Refs. <cit.>, with a full exploration of the J_1-J_2 parameter space reported in Ref. <cit.>. 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 nearest-neighbor antiferromagnetic model. Employing an approach which combines a Luttinger-Tisza analysis with an iterative energy minimization on large system sizes of 32×32×32 cubic unit cells (i.e., 524288 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. <cit.> can be attributed to the substantially reduced finite-size effects in our calculations compared to those of Ref. <cit.>, 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 <ref>] 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. <ref>]. It is shown in Ref. <cit.> that an infinitesimal amount of antiferromagnetic second-nearest-neighbor 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 𝐤=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 𝐤=2π(2,0,0) and symmetry-related points in the EBZ. This result can perhaps be most easily understood by noting that a second-nearest-neighbor interaction J_2 is equivalent to a third-nearest-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 <cit.>. 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 the J̃_αβ^𝐤 matrix in Eq. (<ref>) 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 𝐤=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 𝐤=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 𝐤=2π(k,0,0)-type wave vector as the ordering wave vector. This state is also an eigenstate of the Luttinger-Tisza matrix Eq. (<ref>), thus giving the exact expression k=(2/π) arccos[-J_1/(4 J_2) - 1/2] for the wave vector and an energy per spin of E=-J_1^2/(2 J_2) - 6 J_2. This wave vector differs from the one given in Ref. <cit.> 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 <cit.> 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 𝐤=2π(4/3,0,0) in the EBZ of the pyrochlore lattice. In the planar spiral, and corresponding to the aforementioned dominant peaks at 𝐤=2π(k,0,0)-type ordering wave vectors, there also exist subdominant peaks at ordering wave vectors of the 𝐤=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 <cit.>. 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 𝐤=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 𝐤=2π(3/4,3/4,0) and 𝐤=2π(3/4,-3/4,0), while the second pair, e.g., could be located in the k_y-k_z plane with 𝐤=2π(0,3/4,3/4) and 𝐤=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. <cit.>. 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  𝐤=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 𝐤=2π(5/4,5/4,0) and 𝐤=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 𝐤=2π(0,5/4,5/4) and 𝐤=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 modulated commensurate spiral of Ref. <cit.>, 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 𝐤=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 𝐤=2π(3/4,0,-3/4) type; the zero component is the one which is 1/2 in the dominant 𝐤=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. <cit.> 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 𝐤 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 𝐤=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. <ref>. 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. <ref>] 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 <cit.>. 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 12-sublattice magnetic structure first found on the kagome lattice <cit.>. 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 𝐤=2π(1/2,1/2,1/2) type, e.g., 𝐤=2π(1/2,1/2,1/2), 𝐤=2π(-1/2,1/2,1/2), and 𝐤=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., 𝐤=2π(1/2,1/2,-1/2). The spin configuration in this state can be expressed analytically (see Ref. <cit.>). Each of the dominant ordering vectors is accompanied by a subdominant wave vector of 𝐤=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. <cit.>). 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 𝐤=0 state, the ferromagnetic ordering within the sublattices features subdominant ordering wave vectors at all the 𝐤=2π(1,1,1)-type points in the EBZ, which have a spectral weight of one-quarter of the dominant 𝐤=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 <cit.>. This phase is made up of a family of degenerate ground states with dominant incommensurate wave vectors around the 𝐤=(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 𝐤≈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.



 §.§ 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 <cit.> 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. <ref>. 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. <ref>(a)] and, perhaps surprisingly, also for the S=1 model [see Fig. <ref>(b)]. In Figs. <ref>(a) and <ref>(b), quantum fluctuations are seen to melt away a significant portion (around J_2=0) of the classical domain of existence of the 𝐤=0 and Kawamura magnetic orders. For S=1/2, the paramagnet 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. For S=1/2, we show the representative RG flows within the paramagnetic regime for a point in the antiferromagnetic J_2 regime [Fig. <ref>(a)] and one in the ferromagnetic J_2 regime [Fig. <ref>(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. <ref>(c) for the S=1/2 result] displays weak maxima at 𝐤=2π(2,0,0) (and symmetry-related points), which correspond to the dominant Bragg peak wave vectors of the underlying 𝐤=0 parent classical magnetic order (see Sec. <ref>). Similarly, the spin susceptibility profile for ferromagnetic J_2 [see Fig. <ref>(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. <ref> for the [hhl] plane scattering profiles). Classically, this parameter regime hosts the Kawamura magnetic order with dominant and subdominant Bragg peaks at 𝐤≈2π(5/4,5/4,0) and 𝐤≈2π(3/4,3/4,0) (and symmetry-related points). A comparison of the S=1/2 paramagnetic spin susceptibility profiles, i.e., Fig. <ref>(c) for J_2 antiferromagnetic and Fig. <ref>(d) for J_2 ferromagnetic, with those of the respective parent classical magnetic orders, i.e., 𝐤=0 [Fig. <ref>(b)] and Kawamura [Fig. <ref>(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. <ref> 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 <cit.>, 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 <cit.>. In Fig. <ref>, 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. <ref>). 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. <ref>]. 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. <ref> 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 <cit.>. Starting our discussion with the 𝐤=0 order, we find that its span is considerably diminished for the S=1/2 model [see Table <ref> 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 𝐤=0 phase [J_2/J_1≈ 0.36, marked by a black disk in Fig. <ref>(a)] clearly shows signature of an instability [see Fig. <ref>], indicating the onset of 𝐤=0 magnetic 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. <ref>. The spin susceptibility profile evaluated for J_2/J_1=0.36 at the instability point is shown in Fig. <ref>(b), wherein one observes the dominant IBP at the high-symmetry X points [Fig. <ref>(a)], i.e., 𝐤=2π(2,0,0) (and symmetry-related points), and the subdominant peaks at the L points [Fig. <ref>(a)], i.e., 𝐤=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 <cit.>, 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. <ref>, classically, the collinear 𝐤=0 state is selected by thermal fluctuations <cit.>, and quantum fluctuations are likely to select the same state <cit.>.

The 𝐤=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 evaluated deep inside the spiral ordered phase [J_2=1, marked by a black disk in Fig. <ref>(a)] features an instability at T_c/J_2≈ 0.73(3) [marked by an arrow in Fig. <ref>] pointing to the onset of magnetic order at this temperature. The corresponding spin susceptibility profile evaluated at the instability point is shown in Fig. <ref>(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 <cit.> 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. <ref> 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 𝐤=0 to the DT magnetic order. Meanwhile, the shift 
|δ𝐤|≡|𝐤_ qu|-|𝐤_ cl| from the classical
𝐤_ cl wave vector to the quantum 𝐤_ qu wave vector changes nonmonotonically across the domain of the planar spiral ordered phase [see the inset in Fig. <ref>]. 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 δ𝐤 achieves a maximal value of approximately 4% of the classical value near the boundary to the 𝐤=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. <ref>(a) and Table <ref>], 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. <ref>(a)] and tracked at the dominant wave vector becomes unstable at T/J ≈ 0.39(2) [marked by an arrow in Fig. <ref>]. The corresponding spin susceptibility profile is shown in Fig. <ref>(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. <ref>, 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. <ref>(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. <ref>(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. <ref>(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. <ref>] appears feeble, possibly hinting at the “weakness” of the CS magnetic order. It is of interest to note that the analogous cuboctohedral kagome orders <cit.> found in Heisenberg models with long-range interactions also display an extremely feeble signal of an instability in their RG flow <cit.>. 


[table]capposition=top
