Emergence of a nematic paramagnet via quantum order-by-disorder and pseudo-Goldstone modes in Kitaev magnets

The appearance of nontrivial phases in Kitaev materials exposed to an external magnetic field has recently been a subject of intensive studies. Here, we elucidate the relation between the field-induced ground states of the classical and quantum spin models proposed for such materials, by using the infinite density matrix renormalization group (iDMRG) and the linear spin wave theory (LSWT). We consider the $K \Gamma \Gamma'$ model, where $\Gamma$ and $\Gamma'$ are off-diagonal spin exchanges on top of the dominant Kitaev interaction $K$. We show that a nematic paramagnet, which breaks the lattice-rotational symmetry, emerges in the quantum model. This phenomenon can be understood as the effect of quantum order-by-disorder in the frustrated ferromagnet discovered in the corresponding classical model. In addition, various classical ordering patterns with large unit cells are replaced by the nematic paramagnet in the quantum model. We compute the dynamical spin structure factors using a matrix operator based time evolution and compare them with the predictions from LSWT. We point out the existence of a pseudo-Goldstone mode, which results from the lifting of a continuous degeneracy via quantum fluctuations. We also discuss these results in the light of inelastic neutron scattering experiments.


I. INTRODUCTION
Recently, there has been a surge of interest in Kitaev materials [1][2][3][4][5][6] due to the promise for the discovery of a Kitaev spin liquid [7]. The presence of strong spin-orbit coupling in these materials gives rise to bond-dependent interactions between effective spin-1/2 moments on the underlying honeycomb lattice [8][9][10]. In addition to a substantial Kitaev interaction, there exist other types of spin-exchange interactions [11] in these materials, which ultimately stabilize a magnetic order instead of the desired quantum spin liquid. A paradigmatic example of Kitaev materials is α-RuCl 3 [12], which has a zigzag (ZZ) ordered ground state [13][14][15][16]. Upon applying a magnetic field, however, the ZZ order vanishes [14,[17][18][19][20][21][22], while a half-quantized thermal Hall conductivity is observed [23]. This would be consistent with the theoretical prediction of the non-abelian chiral spin liquid in the pure Kitaev model under a magnetic field [7]. It leads to the speculation that the field-induced phase is indeed the quantum spin liquid.
The experiment mentioned above [23] and recent inelastic neutron scattering experiments [15,20,[24][25][26] have motivated a number of numerical studies [27][28][29][30][31] in an effort to identify possible nontrivial phases in theoretical spin models proposed for realistic materials. For instance, the KΓΓ model, which is considered as a minimal model describing a broad class of Kitaev materials including α-RuCl 3 , has been studied with an external magnetic field. Exact diagonalization on a 24-site cluster shows that the Kitaev spin liquid (KSL) is stable in a window of fields slightly tilted from the [111] direction [27]. In a field along the [111] axis, however, a tensor net-work study reveals that KSL is confined to the vicinity of the pure Kitaev model, while novel nematic paramagnetic states, which breaks the lattice-rotational symmetry, occupy a significant portion of the phase diagram at intermediate fields [31]. On the other hand, the classical model exhibits a plethora of magnetic orders with large unit cells, and a ferromagnetic phase frustrated by the Γ interaction and the field, before the system becomes fully polarized [30,32,33].
Given the diversity of the numerical methods with their own limitations, the different methods offer complementary information, and it would be important and useful to resolve the differences and establish meaningful connections between them. In this work, we employ the infinite density matrix renormalization group (iDMRG) approach [34][35][36] and linear spin wave theory (LSWT) to investigate the nature of the field-induced phases. Furthermore, we present the dynamical properties of the field-induced phases, unexplored in previous theoretical studies, but essential to understand the experimental signatures of Kitaev magnets under the magnetic field.
We first investigate the phase diagram of the KΓΓ model in a magnetic field along the [111] direction. To understand the effect of the finite sizes and geometries, we use various circumferences and geometries on the cylinder. We find that, independent of the choice of geometry, KSL is confined to a small corner of the phase diagram. A nematic paramagnet (NP), which breaks the rotational symmetry but preserves the translational symmetry, appears at intermediate fields before the system enters the longitudinally-polarized paramagnet (P) at high field that is adiabatically connected to the fullypolarized limit. It is worth emphasising that the NP phase originates in the breaking of lattice-rotation symmetry, and is not the type of quantum spin-nematic associated with the breaking of spin-rotation symmetry [37][38][39][40]. We expect that NP is more likely to be observed than the KSL due to the much wider extent of the former in the parameter space. These results broadly agree with those obtained in the tensor network study [31].
We reveal the origin of NP by demonstrating that it arises from a quantum order-by-disorder effect [41,42] in the frustrated ferromagnet (FF) found in the classical model [30]. In the classical limit, the FF appears in a window of intermediate magnetic fields between a series of magnetically ordered phases and the polarized state at low and high fields respectively. The spin orientation in the FF possesses an azimuthal symmetry, which results in a U (1) degenerate manifold of states. Upon the inclusion of zero point quantum fluctuations via the linear spin wave theory (LSWT) [43], a discrete set of azimuthal angles is selected, lifting the U (1) degeneracy. The selection of azimuthal angles, which depends on the field strength, is consistent with the iDMRG results as well as the magnetizations of the NP states identifed by the tensor network study [31]. Furthermore, the series of large unit cell orders [30] at low fields seen in the classical model is replaced by NP in the quantum model.
We then present the dynamical spin structure factors (DSF) of the NP and P phases of the quantum model, using the matrix product operator based time evolution (tMPO) [44], and compare them to the spin wave dispersions as obtained in the semiclassical LSWT approach. In P at high fields, the DSF clearly shows two magnon bands. The magnon excitation gap shrinks as the field decreases towards the transition into NP. At the critical field, the magnetic excitations become gapless at the Γ and Γ points in the reciprocal space, which can be interpreted as the onset of the pseudo-Goldstone mode associated with the lifted U (1) degeneracy [43] upon entering NP. In NP, while there exists relatively well defined magnon excitations at lower energies, a broad continuum forms at higher energies, which is reminiscent of that as seen in inelastic neutron scattering experiments [15,20,[24][25][26]. Approaching the KSL, the excitation continuum becomes much broader, which may be considered as a signature of the proximate KSL [15,45,46].

II. MODEL
We study spin-1 2 degrees of freedom on a honeycomb lattice with bond-dependent interactions. The Hamilto-nian of interest is given by where i, j γ are neighboring sites connected by a bond with label γ ∈ {x, y, z}. The first term is the Kitaev or bond-dependent Ising exchange, which, in the absence of other interactions, stabilize a quantum spin liquid [7]. The second and third terms are the off-diagonal Γ and Γ exchanges. Classically, the pure Γ model is known to host a spin liquid [47], but its quantum ground state is still under debate [47][48][49]. With dominant K < 0 and Γ > 0 interactions, a finite Γ < 0 exchange stabilizes the long-range ordered zigzag phase at zero field [11,50]. The h term describes the Zeeman coupling of the spins to an external magnetic field along the directionĥ. Here, we consider a trigonometric parametrization of the Kitaev and Γ interactions such that We focus on the range 0 ≤ φ ≤ π/2, fix Γ to be either 0 or −0.06 [51], and apply the field in the [111] direction. The [111] field retains the C 3 symmetry of the KΓΓ model.

III. CLASSICAL U (1) DEGENERACY AND QUANTUM ORDER-BY-DISORDER
We start by considering the classical limit of (1) with O(3) spins. At high fields a narrow window with ferromagnetic order exists where the spins are uniformly aligned but canted away from the [111] field. Such a canting originates from the competition between the field and the Γ interaction [30,32,33], thus motivating the name frustrated ferromagnet (FF) for such a phase.
Assuming a ferromagnetic ansatz, i.e. for all sites i, we write S i = S sin ϑ cos ϕâ + sin ϑ sin ϕb + cos ϑĉ ≡ S, (3) whereâ,b andĉ are the unit vectors along the [112], [110] and [111] directions (or the a, b and c-axes) respectively, ϑ = [0, π] is the polar angle measured from the c-axis, while ϕ = [0, 2π) is the azimuthal angle in the ab-plane measured from the a-axis.
Substituting (3) into (1), the energy turns out to depend only on ϑ but not on ϕ. The critical field, above which the spins are fully polarized, is calculated to be h crit = 3 (Γ + 2Γ ) S.  yielding an emergent U (1) manifold of degenerate ground states. When choosing a particular value of ϕ out of the U (1) degenerate manifold, the spins spontaneously break a continuous symmetry. As a consequence, the spin wave spectrum of the FF order exhibits a Nambu-Goldstone mode [52].
We further investigate whether the classical U (1) degeneracy is lifted by quantum fluctuations, demonstrating a concept known as quantum order-by-disorder [41,42]. Within LSWT, the quantum correction to the energy is given by [43] where E cl is the classical energy and ω kn is the magnon dispersion.
As the field is tuned, our quantum order-by-disorder computation reveals a transition from FF γ to FF αβ , e.g. Fig. 1(b). Fixing Γ = 0, we indicate the location of the FF, as well as its type, in the classical phase diagram of (1) as shown in Fig. 1(c).
We consider now the quantum limit of (1) with spin-1 2 constituents. Fig. 2 shows the phase diagram near the KSL for (a) Γ = 0 and (b) Γ = −0.06, as well as (inset) the entire range 0 < φ/π < 0.5 for Γ = 0. We employ iDMRG [34][35][36] on cylinders with L circ = 10 and twisted boundary conditions (see Appendix B). Firstly, we find the KSL to be confined to a small corner near the Kitaev limit, while a large fraction of the phase diagram is occupied by either NP that spontaneously breaks the lattice-rotational symmetry or ZZ long-range magnetic order that further breaks translational symmetry.
At high fields, the ground state is a longitudinallypolarized paramagnet (P) that is adiabatically connected to the fully-polarized limit. The system enters NP upon lowering the field except at very small values of Γ where a direct transition into the KSL occurs. Similar to the FF phase, NP is characterized by a canting of the fieldinduced magnetic moment away from the [111] field. In contrast to FF, NP extends down to h = 0 if Γ = 0, where the magnetic moment vanishes. Since at zero field NP has already broken the C 3 symmetry, as indicated by anisotropic bond energies, a field-induced magnetic moment is no longer protected to align in the [111] direction.  (3) in the corresponding classical model. A finite Γ < 0, i.e. (b) Γ = −0.06, induces ZZ at low fields such that NP occurs at intermediate fields between ZZ and the longitudinally-polarized paramagnet (P) that is adiabatically connected to the fully-polarized limit. The phase diagram (a,b) is obtained on the rhombic-2 cylinder geometry with Lcirc = 10 (see Appendix B). (c-f) Dynamical spin structure factors S(k, ω) of P and NP near the upper critical field at two opposite limits of the KΓ model as illustrated by the red stars in (a). In both limits, the magnon band at low energies, which is apparent in P, becomes diffuse in NP, whereas the continuum at higher energies persists across the transition P to NP. The characteristic feature of NP is the increased spectral weight at Γ and at low energies, which is indicative of the pseudo-Goldstone mode. Remarkably, the continuum exhibits a very different structure in both limits.
As such, the canting can be regarded as a consequence of the broken C 3 symmetry.
As discussed in the previous section, quantum orderby-disorder selects a discrete subset of states out of the classical U (1) degenerate manifold. Our result strongly suggests that the NP phase originates from the FF phase upon the incorporation of quantum fluctuations. Furthermore, the NP phase is much enhanced in the quantum model as compared to the FF phase in the classical model, which we argue as follows. (a) NP has an increased critical field of the transition into P relative to its classical counterpart. (b) NP supersedes many of the classical long-range ordered states [30,31] in a wide region of the phase diagram with ZZ being the only exception.
The transition from P into NP appears to be continuous with a gap closing and a consecutive range of fields with long correlation length. For Hamiltonians with only local interactions, a long correlation length implies a small spectral gap. [53] Here we expect a small spectral gap from the quantum order-by-disorder scenario. In order to verify our hypothesis, we slightly tilt the magnetic field by 1 • from, and rotate it around, the [111] axis. Within the range of fields with a small gap, we find that the canting direction continuously follows the rotated field with only small and smooth variation in the ground state energy (see Appendix C).
As the field is tuned, the fully quantum model exhibits transitions between the two distinct sets NP γ and NP αβ at certain ratios Γ/|K| similar to the transition between FF γ and FF αβ observed with LSWT in Sec. III. Contrary to LSWT and an earlier tensor network calculation [31], however, iDMRG reveals a reentrance behaviour with multiple transitions between NP γ and NP αβ . Moreover, the order in which the sets appear, as well as the critical fields separating them, depend on the geometry used with iDMRG. This is in stark contrast to the critical fields of transitions from P to NP and from NP to ZZ which are relatively insensitive to the geometry. As a consequence NP γ and NP αβ compete in a wide range in parameter space and, therefore, we do not distinguish between them in the quantum phase diagram Fig. 2. Nonetheless, one can construct the corresponding order parameters for NP γ and NP αβ , see Sec. IV A for details.
A ZZ phase at intermediate fields exists embedded within NP for φ/π 0.25, or Γ/|K| 1. Interest-ingly, the upper transition between NP and ZZ appears to be continuous. Further lowering the field, the transition from ZZ back to NP is likely to be first order. In the zero field limit, the nature of the ground state depends on the system size. As a consequence, it is difficult to draw a firm conclusion on the extent of the ZZ phase. While there is no evidence within iDMRG for a ZZ phase at zero field for all the geometries (up to L circ = 12 sites) explored here, the scenario where the ZZ phase ranges down to zero field as in Ref. 48 cannot be excluded. Remarkably, not only ZZ but also various magnetic orders that appear in the classical model (1) are competing with NP in the quantum model at φ/π = 0.5. We refer interested readers to Sec. IV B for more details.
Upon the inclusion of a small Γ = −0.06, the ZZ phase is stabilized at small fields. The KSL shrinks and borders NP, ZZ and P. The critical field of the transition from P into NP is slightly renormalized to a smaller value due to the Γ interaction, as in the classical model. Beyond that, the phenomenology of the phase transitions remains the same as in the case of Γ = 0, namely the transition between P and NP as well as that between NP and ZZ appear to be continuous.

A. Order Parameter of the Nematic Paramagnet
While the canting of the magnetic moments is an indicator of the broken C 3 symmetry, it is not a good order parameter in the zero field limit as the magnetic moments vanish. A more adequate order parameter can be defined in terms of bond energies E γ = E i,j γ , where γ ∈ {x, y, z} is the label of the bonds. In terms of bond energies, NP γ is characterized by a single preferred bond, E γ < E α E β , while for NP αβ two bonds are (almost) equally preferred, E α E β < E γ . Suppose the bond energies are sorted as above and such that E α ≤ E β , then the order parameter of NP γ and NP αβ can be defined as A finite O NPγ or O NP αβ indicates that the C 3 symmetry is broken, in particular as soon as the NP phase is stabilized (see Fig. 3). The C 3 symmetry remains broken down to the zero field limit. Which of the two sets, NP γ and NP αβ , is selected is, however, dependent on the geometry used in iDMRG. While the rhombic geometry prefers NP αβ near the upper critical field and in the zero field limit, the rectangular geometry prefers the NP γ phase. Both geometries have in common an intermediate transition into a complementary phase, but the corresponding critical fields are distinct. Consequently, NP γ and NP αβ remain competing phases in a wide range of fields.
Order parameters ONP γ and ONP αβ of the nematic paramagnetic phases NPγ and NP αβ at φ/π = 0.1 on (a) the rectangular geometry with Lcirc = 12 and (b) the rhombic geometry with Lcirc = 12. Both exhibit transitions between NPγ and NP αβ as the magnetic field is tuned, but the order of occurrence and the critical fields vary.

B. Zero Field Limit of the KΓ Model
The ground state of the KΓ model (i.e. (1) with Γ = 0) at zero field in the range 0 ≤ φ/π ≤ 0.5 has been a subject of debate recently. Various numerical methods provide different answers [31,48,49,54]. Here, we use different cylinder geometries with circumferences as large as L circ = 12 and vary the bond dimension χ of the matrix product state (MPS), which encodes the quantum wave function, to further support our interpretation that NP extends down to zero field. However, the finite circumference affects the ground state energy significantly and a firm conclusion regarding the ground state in the two-dimensional limit cannot be drawn.
In Fig. 4, we compare the ground state energies E GS of different geometries upon varying χ at φ/π = 0.1, 0.25, 0.5. Narrow cylinders with L circ = 6 (3xL y rho) or L circ = 8 (2xL y rec) shows good convergence with respect to χ and feature an almost flat evolution of E GS . In contrast, cylinders with L circ = 10 or 12 require relatively large χ for convergence. Nonetheless, we can read off the following tendencies. Firstly, the rhombic geometry rho exhibits an NP αβ ground state, while the rectangular geometry rec shows a tendency towards the NP γ at sufficiently large χ. Secondly, E GS depends significantly on L circ where the narrowest cylinder 3xL y rho has the lowest energy, which is about 2 to 3% below E GS of cylinders with L circ = 12. This implies that finite-size effects are significant as will be discussed below. Thirdly, some of the magnetically ordered states observed in the  In the limit h −→ 0, the quantum ground state obtained by iDMRG depends on the cylinder geometry and the MPS bond dimension χ. At (a) φ/π = 0.1 and (b) φ/π = 0.25 the ground state is either NP αβ for rhombic unit cells (rho) or NPγ for rectangular unit cells (rec). The rhombic-2 (rho-2) geometry cannot be associated with either of them, as it does not exhibit the corresponding anisotropy in bond energies, namely Eα E β < Eγ (NP αβ ) or Eγ < Eα E β (NPγ) with a, b, c ∈ {x, y, z}. (c) In the pure Γ > 0 model, i.e. φ/π = 0.5, cylinders with Lcirc = 12 (3xLy rec, and 6xLy rho) exhibit transitions from states with finite magnetic moment and long-range order at small χ into either of the NP phase at large χ. In particular, 6x2 rho exhibits a ZZ order for χ ≤ 450, 6x3 rho exhibits a 12-site order for χ ≤ 256, 3x1 rec and 3x2 rec show a 6-site order for χ ≤ 640 3x3 rec exhibits an 18-site order for χ ≤ 640. This indicates that NPγ, NP αβ and several magnetically ordered states observed in the classical model [30] are close in energy.
classical model (1) [30] appear for smaller χ, in particular at φ/π = 0.5, see Fig. 4 (c). The difference in energy between the magnetically ordered states and the NP phase is much smaller than that between different circumferences. Consequently, finite-size effects make it difficult to conclude the precise nature of the ground state near the Γ limit as h −→ 0.
Furthermore, we find NP αβ to be adiabatically connected to uncoupled KΓ chains along alternating α and β bonds [55]. That is to say, we do not observe a phase transition -other than a change to a different NP αβ orientation -upon adiabatically turning off the couplings on all bonds i, j γ with a label γ ∈ {x, y, z}. In particular, consider the rhombic geometry with bond labels as in Fig. 6(a). The ground state obtained by iDMRG is NP yz with E y = E z < E x . We denote the interactions on the x bond collectively as J x , and multiply it by the factor 1 − η with η ∈ [0, 1]. In doing so, we do not find a phase transition upon tuning η from 0 to 1. On the other hand, if the interactions on the y or z bonds are switched off using a similar parametrization, a different orientation of NP αβ is selected at a small but finite η.
The three possible orientations of the NP αβ state are exactly degenerate in the two-dimensional limit. However, this degeneracy splits on the cylinder due to its finite circumference. This is best understood in the limit of the uncoupled KΓ chains. On the rhombic geometry, NP yz corresponds to a KΓ chain composed of y and z bonds along the circumference. Such a chain is finite with L circ sites and subjected to periodic boundary conditions. The remaining two orientations instead correspond to infinite chains. The ground state energies obey E(L = 6) < E(L = 12) < ... < E(L = ∞). Here E(L = 6) is about 1 to 2% smaller than that of L = ∞, which is of the same energy scale as over which E GS spreads for the different geometries (see Fig. 4). Thus, the significant finite-circumference -or finite-size -effects of the two-dimensional KΓ model can be explained, at least within the NP αβ phase, by its relation to the KΓ chains.

V. DYNAMICS
The dynamical spin structure factor S(k, ω) contains information about the excitation spectrum and can be probed by inelastic neutron scattering experiments. We consider S(k, ω) = γ∈{x,y,z} S γγ (k, ω) with S γγ (k, ω) being the spatio-temporal Fourier transform of the dynamical correlations where γ ∈ x, y, z is the spin component, r a and r b are the spatial positions of the spins, N is a normalization factor defined via dω dk S γγ (k, ω) = dk, and C γγ ab (t) denotes the dynamical spin-spin correlation C γγ ab (t) = ψ 0 |S γ a (t)S γ b (0)|ψ 0 . Given the ground state wave function in MPS form, the time evolution is carried out using the tMPO method [44] until t = 60. The time series is then extended by a linear prediction [56][57][58] and multiplied with a Gaussian distribution. The resulting line broadening amounts to σ ω = 0.036. We restrict ourselves to magnetic fields near the transition between P and NP, where the numerical errors and finite size effects are found to be small. A rhombic geometry with In (a)-(c) P exhibits two magnon bands, with the lower one remaining sharp but the upper one being obscured by the multi-magnon continuum upon decreasing the field. Comparing to the dispersion obtained by LSWT, the lower magnon band (solid line) is shifted downwards in energy and hits zero at the Γ and Γ point at a higher critical field than that predicted classically. Within NP, as in (d), the lower magnon band becomes slightly diffusive. A broad continuum appears, ranging from (almost) zero frequency to about ω = 1.6. The spectral weight is mostly concentrated at low frequencies near the Γ-point, indicating long wavelength excitations corresponding to a pseudo-Goldstone mode. LSWT spectrum in the FF phase, or NP in the quantum limit, is gapless by the Nambu-Goldstone theorem. Some of the spectral weight is redistributed to the M -point, signaling that a small perturbation, e.g. Γ , may be sufficient to drive the system into the zigzag (ZZ) order.
L circ = 6 sites is used such that three separated lines of accessible momenta exist in the first Brillouin zone, as illustrated in Appendix B. We present two sets of S(k, ω), Fig. 2(c-f) and Fig. 5, along the high-symmetry direction Γ − M − Γ in the Brillouin zone. In Fig. 2(c-f) we want to emphasize the similarities and differences of S(k, ω) in the two limits (c,d) Γ −→ 0 and (e,f) K −→ 0. Whereas in Fig. 5 we focus on the evolution of the magnon bands upon reducing the field across the critical field while keeping φ/π = 0.1 fixed. S(k, ω) for additional parameters are presented in Appendix D.
As highlighted before, the NP phase remains stable in a wide range in parameter space. NP extends from the almost (FM) Kitaev limit, φ −→ 0, to the (AF) Γ limit, φ/π −→ 0.5 and beyond when Γ = 0. In the entire range, the lower magnon band which is apparent in P reduces in energy upon approaching the critical field of the transition P to NP, where eventually it condenses at Γ and Γ . Within NP, the main spectral weight is observed at long wavelengths and low frequencies above a small spectral gap. Such excitations are consistent with the slightly lifted accidental U (1) degeneracy and correspond to the slow pseudo-Goldstone mode and is the characteristic feature of NP.
The continuum, however, exhibits a very distinct structure in both limits of φ. Near φ −→ 0, that is near the pure Kitaev model with K < 0, S(k, ω) bears similarities to that of the nearby ferromagnetic KSL. The broad continuum of the KSL is formed by itinerant Majorana fermions subject to a quantum quench caused by the excitation of a flux pair when acting with a spin operator on the ground state [59,60]. Since the flux-pair excitation is gapped, S(k, ω) has a finite gap even though the Majorana fermion spectrum is gapless. The contin-uum above the spectral gap has an almost k-independent width equal to that of the single Majorana fermion dispersion. S(k, ω) of NP and P proximate to the KSL exhibit a similar continuum with the same upper cutoff frequency ω ≈ 1.6. Thus, remnants of the Majorana fermion excitations of the KSL appear to persist at higher energies. Differences between NP and KSL are most obvious at lower energies. While the KSL has a relatively low-lying dispersion and the spectral gap, the NP exhibits a somewhat broadened dispersing low-energy mode associated with the pseudo-Goldstone mode. As a consequence, it may be challenging to discern NP and the nearby KSL using inelastic neutron scattering experiments. Such a phenomenology, which is known as proximate spin liquid [15], was previously observed in the Kitaev-Heisenberg model [45,46].
In the other limit φ/π −→ 0.5, i.e. the pure Γ model with Γ > 0, S(k, ω) of NP at h = 1.4 is dominated by diffuse remnants of the magnon modes that exist in P at h = 1.6, compare Fig. 2(e) and (f). The magnon mode has an enhanced bandwidth and dispersion which remains of equal magnitude upon approaching the transition. This is in contrast to the pure Kitaev model in a [111] field, which exhibits a band that simultaneously flattens and approaches zero in the entire reciprocal space [61]. Moreover, the continuum within NP has more structure and an increased energy scale, ∆ω/Γ ≈ 2.5, compared to that of the KSL, ∆ω/K ≈ 1.6. In comparison to the Γ model (with Γ > 0) at zero field [62], we observe a similar energy scale of the continuum. At low energies, however, the main spectral weight shifts from M at zero field, to Γ at high field, which corresponds to the pseudo-Goldstone mode. Above the pseudo-Goldstone mode, in particular at Γ and Γ , the continuum exhibits a gap similar to the gap observed at zero field. Figure 5 demonstrates the evolution of the magnon modes upon approaching the critical field. Two magnon bands exists due to the unit cell having two sites. The two magnon modes are most pronounced at high fields, e.g. Fig. 5(a) at h = 2.0, where the magnetic moment is close to saturation. Agreement with LSWT is good except from an overall shift of the lower magnon band towards smaller frequencies. Upon reducing the field, e.g. Fig. 5(b) at h = 1.2, both magnon bands move downwards in energy, while their shapes are essentially the same as those at h = 2.0. However, the upper band starts to overlap with the two-magnon continuum resulting in a broader linewidth due to additional decay channels. The bottom of the two-magnon continuum as obtained from LSWT is marked by the dashed line. In the fully quantum computation (tMPO), however, the continuum starts at lower frequencies due to the overall reduction of the lower magnon band as compared to LSWT.
Right before entering NP, e.g. Fig. 5(c) at h = 0.8, the lower magnon band approaches zero frequency at the Γ and Γ points, illustrating a condensation of magnons and the imminent uniform canting of the magnetic moments. LSWT predicts a gap, as the classical model is still in the fully-polarized phase due to a lower h crit . The continuum overlaps with the upper magnon band and obscures it except near the Γ point. Upon the transition from P to NP, the continuum does not change significantly and, thus, confirms that the transition is continuous.
Within NP, e.g. Fig. 5(d) at h = 0.4, a broad continuum appears and extends up to ω ≈ 1.6. Furthermore, a softening of the DSF at the M point signals the development of a significant ZZ correlation when the magnetic field is lowered. Adding a small Γ < 0 indeed leads to the ZZ order as in Fig. 2(b). The ZZ correlation may be related to the presence of an intermediate ZZ order that takes place at larger Γ/|K| (φ/π 0.25) in the quantum phase diagram Fig. 2(a). In the corresponding classical model, the intermediate ZZ order extends down to φ/π = 0.1 [30]. Thus, the ZZ correlation may also be a remnant of the classical ZZ order. The LSWT spectrum is gapless at the Γ point as the classical model is now in the FF phase. Analogous to tMPO, LSWT predicts a condensation of magnons at the M point, which becomes more pronounced upon further lowering the field and approaching the transition into ZZ in the classical model (see Appendix E).

VI. DISCUSSION
We demonstrate that the KΓ and KΓΓ models in a magnetic field along the [111] direction support an NP phase, which is not magnetically ordered, but breaks the lattice-rotation symmetry while preserving translational symmetry. We trace the origin of NP to the FF phase in the corresponding classical model, which has a U (1) degenerate manifold of states. When zero-point quantum fluctuations are incorporated, the U (1) degeneracy is bro-ken down to a discrete subset of ground states related by the C 3 symmetry. This quantum order-by-disorder effect leads to a pseudo-Goldstone mode [43] in the excitation spectrum, as shown in the dynamical spin structure factors of the quantum model.
In realistic Kitaev materials like α-RuCl 3 , the latticerotational symmetry has already been broken by a monoclinic distortion [14,63]. Therefore, it is reasonable to expect that the degeneracy of NP is lifted, i.e. one out of the six canting directions is favored. In fact, a substantial magnetic torque is measured for α-RuCl 3 at high fields [64], indicating that NP may be stabilized in α-RuCl 3 . While monoclinic distortion may have already caused a small canting of the magnetic moment, we argue that the canting of magnetic moments in NP is much more pronounced. It would be necessary to look for a second transition at high fields to verify the existence of the field-induced NP, which is sandwiched between the long-range zigzag order and the longitudinally-polarized paramagnet. The latter only exhibits a negligible torque, while NP has a much larger torque due to the induced canting. The specific nature of the distortion may either enhance or mitigate the extent of NP.
Recently, the iridate K 2 IrO 3 was proposed to exhibit ferromagnetic Kitaev exchange and large offdiagonal exchange while maintaining a C 3 point group symmetry [65]. Considering the wide range in parameter space over which NP occurs, this phase may as well be relevant in K 2 IrO 3 .
An additional implication for experiments is drawn from the dynamical spin structure factor of the NP phase, which exhibits diffusive scattering features masking sharp magnonic excitations. Such a scattering continuum bears some similarities to that of the nearby KSL. With a greater ratio of Γ/|K|, the continuum appears in a wider range of energies. Therefore, the excitation continuum observed in inelastic neutron scattering experiments on α-RuCl 3 [15,20,[24][25][26] could originate from the NP phase. The intriguing question of whether NP is a topologically non-trivial phase is an excellent subject of future study.

VII. ACKNOWLEDGEMENTS
We are grateful to R. Kaneko, K. Penc, J. G. Rau In this section, we work in the crystallographic (abc) coordinates, and write the spin components in (3) as S = S a , S b , S c . For a ferromagnetic order, the Hamiltonian (1) reduces to [30] where N is the total number of unit cells, γ ≡ Γ + 2Γ (assuming nonzero), and the Lagrange multiplier λ has been introduced to constrain the spin magnitude. As any ferromagnetic order saturates the lower bound of the classical energy of the Kitaev interaction with K < 0 [66], we can simply drop it from (A1) in the subsequent analysis. Extremizing H with respect to the variables S a , S b , S c and λ leads to the following equations We study the case where the field is completely aligned in the [111] direction, i.e. h a = 0 and h b = 0, but h c = 0. (A2a) becomes −2γS a + 2λS a = 0. (A2b) yields a similar condition.
(A2c) then implies S c = h c /3γ, which is a physical solution only when h c /3γ ≤ S.
In conclusion, the frustrated ferromagnet (FF) can be realized only when h c < 3γS. As h c ≥ 3γS, the system becomes fully polarized. We thus identify 3γS as the critical field h crit .

Appendix B: Geometries and Accessible Momenta
To enable the use of iDMRG, we wrap the twodimensional honeycomb lattice on the cylinder and wind the one-dimensional MPS structure around the cylinder. Due to the cylindrical geometry, only discrete steps of the wave vectors k 2 /L 2 are available, where k 2 is the reciprocal lattice vector corresponding to n 2 and L 2 is the number of unit cells along the circumference. This results in lines of accessible momenta in the reciprocal space. The way in which the lower and the upper boundaries are connected determines the orientation of these lines, while the circumference determines the spacing between them. Here, we employ three different unit cells, rhombic, rectangular, and rhombic-2. See Fig. 6 for an illustra-  (1), we tilt the magnetic field slightly from the [111] axis by ϑ = 1 • and then rotate it by varying ϕ. The magnetic moment follows the rotating field smoothly with a slight deviation up to 7 • . As discussed in Sec. III, the classical U (1) degeneracy is lifted by quantum fluctuations. The induced gap is of the order of ε ≈ 10 −4 and increases upon lowering the field. The minima are located at ϕ = 0, 2 3 π, and 4 3 π corresponding to FF αβ or NP αβ , respectively. tion of each geometry and their corresponding accessible momenta. For all three geometries, we restrict ourselves to those circumferences L circ that include the K point in the Brillouin zone. In particular, we use (a) rhombic two-site unit cell with L 2 = 3 and 6, which results in L circ = 6 and 12 sites, (b) rectangular four-site unit cell with L 2 = 2 and 3, which results in L circ = 8 and 12 sites, and (c) rhombic-2 two-site unit cell with a shifted boundary condition, such that the vector L 2 n 2 is equivalent to n 1 , and we choose L 2 = 5 resulting in L circ = 10.
The classical model (1) features many different magnetic orders with unit cells containing up to 98 sites [30]. By employing various geometries, in particular (a) and (b) together with L 2 = 1, 2, 3, we allow ordering patterns with up to 18 sublattices per unit cell to fit in at least one of the geometries. As shown in the Sec. III, quantum order-by-disorder lifts the accidental U (1) ground state degeneracy in the classical model (1). In order to see whether this happens in the quantum model using iDMRG, we apply a magnetic field h = h sin ϑ cos ϕâ + sin ϑ sin ϕb + cos ϑĉ (C1) that is slightly tilted away from and rotated around the [111] axis. Here, we choose ϑ = 1 • . Upon rotation of the tilted field around the [111] axis by varying ϕ ∈ [0, 2π), we observe the following. Firstly, the variation in E(ϑ = 1 • , ϕ) is generally small, of the order 10 −4 . The variation is the smallest near the P to NP transition and increases upon lowering the field. Secondly, the tilting angle ϕ m of the magnetic moments continuously follows that of the field ϕ h , though there is a small deviation |ϕ m −ϕ h | 7 • . Thirdly, E(ϑ = 1 • , ϕ) has a period of 2/3π. Deviations arise due to the finite circumference and become more pronounced as the field is decreased. The first observation is consistent with the quantum order-by-disorder scenario since quantum fluctuations only result in small corrections to the energy. As a result, the gap induced by quantum order-by-disorder is small and can be easily overcome by small perturbations, e.g. a slight tilting of the field. While the second observation is a consequence of the first one, it also implies that one can continuously tune from NP αβ to NP γ as well as between different orientations thereof. Moreover, NP γ and NP αβ are related by the same spontaneous symmetry breaking.
The third observation is linked to the broken C 3 symmetry. If the C 3 symmetry is intact, applying a C 3 rotation on the ground state leaves it invariant and does not result in a degeneracy. However, if the C 3 symmetry is spontaneously broken, acting with a C 3 rotation on the ground state transforms the state to a different state with the same energy. This implies at least a threefold degeneracy which is manifest in the 2π/3 periodicity in Fig. 7. Deviations from an exact threefold degeneracy are caused by the cylindrical geometry. They are most pronounced for small L circ and upon reducing the field.
In summary, our observations from tilting and rotating the magnetic field away from and around the [111] axis are consistent with the scenario of NP γ and NP αβ being related to the classical FF phase. Quantum fluctuations induce a small gap and only a discrete subset of states is selected out of the U (1) manifold in the classical model. In the pure Γ model, φ/π = 0.5, the dispersion is still present although reduced. In the entire region φ/π = (0, 0.5] the lower band attains its minimum at the Brillouin zone center, i.e. the Γ and Γ points, and reaches zero frequency at the transition between P and NP. Within NP the main spectral weight is located at low frequencies and long wavelengths, which is consistent with the existence of a pseudo-Goldstone mode. Some spectral weight is distributed around the M point, where the band bends towards lower energies.
Thus, a small perturbation, e.g. Γ < 0, may trigger a transition into the ZZ order.
uniform canting of spins in the NP phase. Regardless of the specific value of φ, S(k, ω) of the NP phase features a nearly vanishing gap at Γ that opens slightly upon further lowering the field. The small gap gives rise to a multi-particle continuum near the zero energy, partially obscuring the single-magnon mode by enabling decay channels. The spectral weight is mostly concentrated just above the gap at the Γ, signifying lowenergy excitations with long wavelengths. Such excitations correspond to the pseudo-Goldstone mode of the slightly lifted U (1) degenerate manifold of states. Moreover, S(k, ω) exhibits a distinctive feature at high energies above Γ , which, as suggested by LSWT, appears to be a remnant of the upper branch of the one-magnon excitation.
Superimposed in Fig. 8 are the two single-magnon bands and the resulting lower edge of the two-magnon continuum calculated by LSWT. While in the P phase the bandwidth and the dispersion generally agree with the DSF computed by tMPO, the lower band of the latter is shifted towards lower energies. Within NP (or FF in the classical model), LSWT exhibits a gapless excitation at Γ and Γ that is the Nambu-Goldstone mode arising from spontaneously breaking a continuous symmetry.
In the limit φ −→ 0 the system enters the KSL phase, which is characterized by the fractionalization of spins into itinerant Majorana fermions with a background of static Z 2 fluxes [7]. The fluxes are gapped excitations of an emergent Z 2 gauge field. S(k, ω) of the KSL features a spectral gap ∆ equal to the energy of a flux pair. A continuum starts at ∆ and has a width equal to that of the single Majorana fermion dispersion [59,60]. The continuum is almost entirely formed by Majorana fermions subject to a quantum quench that caused by exciting a flux pair [59].
In the NP or P phase proximate to the KSL, see Fig. 5(c,d) for φ/π = 0.01 and Fig. 8 for φ/π = 0.05, we observe an equally wide continuum suggesting a likewise description of the high-energy excitations in terms of Majorana fermions. However, neither P nor NP is close to the flux-free state of the KSL and, in fact, fluxes are no longer conserved quantities in the presence of a magnetic field or Γ. Consequently, S(k, ω) at lower energies is distinct from that of the KSL. More precisely, the continuum does not have a constant spectral gap, dispersion evolves, and within P a distinctive magnon band appears.
Upon further increasing φ and thus Γ/|K|, S(k, ω) exhibits enhanced dispersions at both the lower and the upper edges of the continuum in NP and P. Among the values of φ examined here, the magnon bands at φ/π = 0.25 has the largest bandwidth in the P phase. In the pure Γ model (i.e. φ/π = 0.5) the dispersion is still present but its bandwidth is reduced.
Within NP at φ/π = 0.05, 0.15, the excitation gap decreases at the M point together with a significant redistribution of spectral weight there indicating the enhancement of ZZ correlations. Subsequently, a small perturbation like a negative Γ may stabilize a long-range ZZ ordered phase. The LSWT spectrum at φ/π = 0.05 exhibits the same phenomenology due to the nearby transition into the ZZ phase present in the classical model.
Along the lines of accessible momenta not presented here, we observe that the gap between the two magnon bands closes at the K-point at φ/π = 0.15, both in LSWT and tMPO. This is due to a duality transformation that exists in the parameter space of the JKΓΓ in a [111] field mapping the KΓ model at φ = tan −1 (1/2) ≈ 0.148π to the pure FM Heisenberg model at high field. [67] For smaller as well as larger φ, the magnon bands are well separated and known to be topological with non-zero Chern numbers implying edge modes on geometries with open boundaries. [67,68] Moreover, the magnetic moment just above the critical field is expected to be fully polarized only at the φ dual to FM Heisenberg, where the fully polarized stated is in eigenstate. Anywhere else, frustration leads to a reduced magnetic moment which approaches the fully polarized state in the limit h −→ ∞. This motivates to distinguish the longitudinallypolarized paramagnet in the quantum model from the fully polarized phase in the classical model.
At φ/π = 0.5, i.e. the pure Γ > 0 model, the P phase follows the same phenomenology as at smaller φ, see Fig. 5(e,f) and Fig. 8. The dispersing magnon bands shift downwards in energy and the excitation gap closes at the Γ and Γ point upon the transition into the NP phase. As the field is lowered, a multi-magnon continuum starts to form within P and persists across the transition into NP. In contrast, the pure Kitaev model in a [111] field exhibits a band that simultaneously flattens and approaches zero in the entire reciprocal space [61]. At φ/π = 0.5 and h ∈ [1.4, 1.8], the classical ground state is some 6-site magnetic order distinct from P or FF, which is why the LSWT spectra are not plotted.

Appendix E: Additional Plots of Linear Spin Wave Dispersion
The classical limit of the KΓ model exhibits a phase transition from the FF phase to the ZZ long-range order upon lowering the field for a wide range of φ (see Ref. 30 for details). We examine the dispersion obtained via LSWT within the FF phase near this transition. We find that, as the field decreases, the excitation gap at the M point shrinks and approaches zero, as shown in Figs. 9a-9c. This is also seen in the DSF of the quantum model, as discussed in Sec. V and Appendix D. We remark that the DSF at a certain choice of h agrees better with the LSWT spectrum at a lower h than that at the same h.