Emergence of a nematic paramagnet via quantum order-by-disorder and pseudo-Goldstone modes in Kitaev magnets
Matthias Gohlke Li Ern Chern Hae-Young Kee Yong Baek Kim
cond-mat.str-el

	


Theory of Quantum Matter Unit, Okinawa Institute of Science and Technology Graduate University,
  Onna-son, Okinawa 904-0495, Japan

Department of Physics, University of Toronto, Toronto,
  Ontario M5S 1A7, Canada

Department of Physics, University of Toronto, Toronto,
  Ontario M5S 1A7, Canada
Canadian Institute for Advanced Research/Quantum
  Materials Program, Toronto, Ontario MSG 1Z8, Canada

Department of Physics, University of Toronto, Toronto,
  Ontario M5S 1A7, Canada
School of Physics, Korea Institute for Advanced Study,
  Seoul 130-722, Korea
     
	


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 ΓΓ' model, where Γ and Γ' are off-diagonal
spin exchanges on top of the dominant Kitaev interaction K.
Focusing on the magnetic field along the [111] direction,
we explain the origin of the nematic paramagnet,
which breaks the lattice-rotational symmetry and exists in an extended window of magnetic field,
in the quantum model.
This phenomenon can be understood as the effect of quantum order-by-disorder
in the frustrated ferromagnet with a continuous manifold of degenerate ground states
discovered in the corresponding classical model.

We compute the dynamical spin structure factors using a matrix operator based time evolution
and compare them with the predictions from LSWT.
We, thus, provide predictions for future inelastic neutron scattering experiments on Kitaev
materials in an external magnetic field along the [111] direction.
In particular, the nematic paramagnet exhibits a characteristic pseudo-Goldstone mode
which results from the lifting of a continuous degeneracy via quantum fluctuations.



	Emergence of nematic paramagnet via quantum order-by-disorder 

and pseudo-Goldstone modes in Kitaev magnets
    Yong Baek Kim
    December 30, 2023
============================================================================================================








§ INTRODUCTION
 

Recently, there has been a surge of interest in Kitaev materials <cit.>
due to the promise for the discovery of a Kitaev spin liquid <cit.>.
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 <cit.>.
In addition to a substantial Kitaev interaction, there exist other types of spin-exchange interactions <cit.> 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 <cit.>,
which has a zigzag (ZZ) ordered ground state <cit.>.
Upon applying a magnetic field, however, the ZZ order vanishes
<cit.>,
while a half-quantized thermal Hall conductivity is observed <cit.>.
This would be consistent with the theoretical prediction of the non-abelian chiral spin liquid
in the pure Kitaev model under a magnetic field <cit.>.
It leads to the speculation that the field-induced phase is indeed the quantum spin liquid.

The experiment mentioned above <cit.>
and recent inelastic neutron scattering experiments <cit.>
have motivated a number of numerical studies <cit.>
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 <cit.>.
In a field along the [111] axis, however, a tensor network 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 <cit.>.
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 <cit.>.



In this work, we investigate the origin of the nematic paramagnet in the 
quantum K ΓΓ' model and its relation to the classical ground states,
employing the infinite density matrix renormalization group (iDMRG) approach
<cit.>
and linear spin wave theory (LSWT).
Our main findings are:

    
  * The nematic paramagnetic state arises due to quantum order-by-disorder effect
          in the continuous manifold of degenerate classical ground states,
          dubbed frustrated ferromagnet, in the corresponding classical model. 
    
  * The dynamical structure factors for the nematic paramagnetic state
          reveal the existence of a pseudo-Goldstone mode
          resulting from the quantum order-by-disorder effect.
          Our predictions for the dynamical structure factors can be compared to
          future neutron scattering experiments to confirm the nematic paramagnetic state.
          The dynamical structure factors in the corresponding regions in the phase diagram
          are shown in Fig. <ref>.


Furthermore, in the process of achieving these goals,
we attempt to resolve the differences and establish meaningful connections between
the results of different numerical studies,
with the understanding that different numerical methods possess their own limitations
and offer complementary information.
We elaborate all of these findings in a concise fashion in the rest of the introduction.
Further and complete details can be found in the main text.

We first investigate the phase diagram of the KΓΓ' model in a magnetic field along the [111] direction.
In agreement with Ref. [lee_magneticfield_2019], we find the KSL to be confined to a small corner of the phase diagram,
while a nematic paramagnet (NP) appears at intermediate fields before the system enters the longitudinally-polarized paramagnet (P) at high field
that is adiabatically connected to the fully polarized limit.
The NP phase breaks lattice-rotational symmetry but preserves translational symmetry.
In that context, it is worth emphasizing that NP is not the type of quantum spin-nematic associated with the breaking of spin-rotation symmetry
<cit.>.
We expect that NP is more likely to be observed than the KSL due to the much wider extent of the former in parameter space. 

Next, we reveal the origin of NP by demonstrating that it arises from 
a quantum order-by-disorder effect <cit.> 
in the frustrated ferromagnet (FF) found in the classical model <cit.>.
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) <cit.>,
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 our iDMRG results as well as the magnetizations of the NP states identifed by the tensor network study <cit.>.
Furthermore, the series of large unit cell orders <cit.> observed 
at low fields 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) <cit.>,
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 <cit.> 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
<cit.>.
Approaching the KSL, the excitation continuum becomes much broader, which may be considered
as a signature of the proximate KSL <cit.>.
We suggest to perform neutron scattering experiments with the [111] magnetic field to test
our predictions for the pseudo-Goldstone mode and the nematic paramagnetic state. 






§ MODEL
 
We study spin-1/2 degrees of freedom on a honeycomb lattice with bond-dependent interactions. The Hamiltonian of interest is given by

    ℋ =      ∑_⟨ i,j ⟩_γ K_γ S^γ_i S^γ_j 
    		 + ∑_⟨ i,j ⟩_γ, α,β≠γΓ_γ( S^α_i S^β_j + S^β_i S^α_j) 
       + ∑_⟨ i,j ⟩_γ, α,β≠γΓ'_γ( S^α_i S^γ_j + S^β_i S^γ_j + S^γ_i 
    S^α_j + S^γ_i S^β_j ) 
       - h ∑_i 𝐡̂·𝐒_i ,

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 <cit.>. The 
second and third terms are the off-diagonal Γ and Γ' exchanges. 
Classically, the pure Γ model is known to host a spin liquid 
<cit.>, but its quantum ground state is still
under debate <cit.>.
With dominant K<0 and Γ>0 interactions, a finite Γ'<0 exchange stabilizes the
long-range ordered zigzag phase at zero field
<cit.>.
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 

    K= -cosϕ    and    Γ = sinϕ .

We focus on the range 0 ≤ϕ≤π/2,
fix Γ' to be either 0 or -0.06
[The choice of Γ'=-0.06 is based on
    (a) Γ' being sufficiently large to induce ZZ order and 
    (b) Γ' being sufficiently small to keep a small region of KSL in the phase diagram when using iDMRG.
    Consequently, our choice for Γ' differs slightly from the ones used in other works <cit.>.],
and apply the field in the [111] direction.
The [111] field retains the C_3 symmetry of the KΓΓ' model.






§ CLASSICAL U(1) DEGENERACY AND QUANTUM ORDER-BY-DISORDER
 



We start by considering the classical limit 
of (<ref>) 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
<cit.>,
thus motivating the name frustrated ferromagnet (FF) for such a phase.

Assuming a ferromagnetic ansatz, i.e. for all sites i, we write

    𝐒_i = S ( sinϑcosφ 𝐚̂
                     + sinϑsinφ 𝐛̂
                     + cosϑ 𝐜̂)
                     ≡𝐒 ,

where 𝐚̂, 𝐛̂ and 𝐜̂ are 
the unit vectors along the [112̅], [1̅10] 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 (<ref>) into (<ref>),
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.
Below h_crit, the spins are canted away from the [111] direction by ϑ = cos^-1( h/h_crit),
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 <cit.>.

We further investigate whether the classical U(1) degeneracy is lifted 
by quantum fluctuations, demonstrating a concept known as quantum order-by-disorder <cit.>.
Within LSWT, the quantum correction to the energy is given by <cit.>


    E = S^2 [ ( 1+1/S) E_cl + 1/S1/2∑_𝐤nω_𝐤n + O ( 1/S^2) ]  ,


where E_cl is the classical energy and ω_𝐤n is the magnon dispersion.

Our calculation shows that either of the following sets of azimuthal angles,
{π/3, π, 5π/3 } and { 0, 2π/3, 4π/3 },
are more energetically favorable than others (see Fig. <ref>(a) and (b)).
Within each set, the three angles yield the same energy,
i.e. the U(1) degeneracy is broken down to a C_3 degeneracy.
For φ=π/3, π and 5π/3, the spins cant towards the cubic 
axes [010], [001] and [100] respectively,
leading to S^α = S^β < S^γ, where α, β, γ∈{x,y,z}
denote the spin components and are also related to the labels of bonds. 
As we will show later, the selected state is related to the NP phase occurring in the quantum model.
We therefore denote this set as NP_γ.
On the other hand, for φ=0, 2π/3, and 4π/3,
the spins cant towards the [110], [011] and [101] directions respectively,
leading to S^α = S^β > S^γ, which we denote as NP_αβ.
The magnetization of NP_γ (NP_αβ) coincide with that of the NP1 (NP2) phase reported in Ref. [lee_magneticfield_2019]. 

As the field is tuned, our quantum order-by-disorder computation reveals a transition from NP_γ to NP_αβ, e.g. Fig. <ref>(b).
Fixing Γ'=0, we indicate the location of the FF, as well as the corresponding NP states selected by quantum fluctuations,
in the classical phase diagram of (<ref>) as shown in Fig. <ref>(c).






§ QUANTUM SPIN-1/2 GROUND STATE
  



We consider now the quantum limit of (<ref>) with spin-1/2 constituents. 
Fig. <ref> 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 <cit.>
on cylinders with various geometries and circumferences as illustrated in Fig. <ref>.
The phase diagram is obtained on the rhombic-2 geometry with L_circ=10 sites
and twisted boundary conditions.
The other geometries with upto L_circ=12 are used checking against possible classical ordering patterns,
see Sec. <ref> for more details.

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 longitudinally-polarized 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 field-induced 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.
As such, the canting can be regarded as a consequence of the broken C_3 symmetry.

As discussed in the previous section, quantum order-by-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
<cit.> 
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 <cit.>.
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 <ref>). The magnitude of the variation is of similar order as in LSWT. 

As the field is tuned, the fully quantum model exhibits transitions between different NP states at certain ratios Γ/|K|
similar to the transition between NP_γ and NP_αβ observed with LSWT in Sec. <ref>.
Contrary to LSWT and an earlier tensor network calculation <cit.>, however,
iDMRG reveal multiple transitions between NP_γ, NP_αβ, and other NP states.
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 different NP states compete in a wide range in parameter space
and, therefore, we do not distinguish between them in the quantum phase diagram Fig. <ref>.
In Sec. <ref> we introduce an order parameter for NP and discuss the different NP states in more detail.

A ZZ phase at intermediate fields exists embedded within NP for ϕ/π≳ 0.25, or 
Γ/| K |≳ 1.
The upper transition between NP and ZZ is likely continuous.
Further lowering the field, the transition from ZZ back to NP appears to be first order 
and its critical field depends on the geometry used.
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 long-range magnetically ordered phase
at zero field for all the geometries (up to L_circ=12 sites) explored here,
scenarios for the two-dimensional limit in which either the ZZ phase ranges down to zero field <cit.>
or other magnetic orderings <cit.> are stabilized cannot be excluded.
In fact, various magnetic orders that appear in the classical model (<ref>)
are obtained as meta-stable states at small χ and, thus,
these orderings are competing with NP in the quantum model.
We refer interested readers to Sec. <ref> 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.



 §.§ 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 suitable 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_γ = ⟨ℋ⟩_⟨ i,j ⟩_γ, where γ∈{x,y,z} is the label of the bonds.

The broken C_3 symmetry manifests in different E_γ on each bond,
while if the C_3 symmetry is not broken, e.g. at high fields, the bond energies satisfy 
E_x = E_y = E_z. 
This motivates to define an order parameter 𝒪^γ_C_3 on each bond γ
by the difference of E_γ to the average bond energy as

    𝒪_C_3^x    = E_x - 1/3(E_x + E_y + E_z ), 
    𝒪_C_3^y    = E_y - 1/3(E_x + E_y + E_z ),  
    𝒪_C_3^z    = E_z - 1/3(E_x + E_y + E_z ).


In terms of bond energies, NP_γ is characterized by a single preferred bond,
E_γ < E_α≃ E_β and hence
𝒪^γ_C_3 < 0 < 𝒪^α_C_3 = 𝒪^β_C_3,
while for NP_αβ two bonds are equally preferred,
E_α≃ E_β < E_γ and hence
𝒪^α_C_3 = 𝒪^β_C_3 < 0 < 𝒪^γ_C_3.
An arbitrary NP state has three different bond energies, E_α≠ E_β≠ E_γ, and hence 
𝒪^α_C_3≠𝒪^β_C_3≠𝒪^γ_C_3.



In any case, a finite 𝒪^γ_C_3 indicates 
that the C_3 symmetry is broken, in particular as soon as the NP phase is stabilized (see Fig. <ref>).
The C_3 symmetry remains broken down to the zero field limit.
Whether NP_γ, NP_αβ, or a different NP state gets 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.
Moreover, both geometries have in common an intermediate transition into a phase,
which exhibits the homogeneous canting of the magnetic moments characteristic for NP, 
but the azimuthal angle φ neither belongs to NP_γ nor NP_αβ.
Instead, degenerate states with φ = n π/3 ±δ occur,
where n is an odd (even) integer corresponding to NP_γ (NP_αβ).
In the two-dimensional limit, this φ implies a six-fold degeneracy.  

In conclusion, the different sets of NP states compete over a wide range of fields,
and their difference in energy appears to be of the order of,
or smaller than, finite-size corrections. 
Whether the intermediate NP state is an effect of quantum order-by-disorder
or of the finite circumference has to remain an open question.
In that context one should recall that the quantum order-by-disorder computation
which we present in Sec. <ref> is semi-classical. 
Non-linear terms, which were not included, may be relevant to stabilize an intermediate NP phase
different from NP_γ or NP_αβ and with a six-fold degeneracy.



 §.§ Zero Field Limit of the K Γ Model
 



The ground state of the K Γ model (i.e. (<ref>) with Γ'=0) at zero field
in the range 0 ≤ϕ/π≤ 0.5 has been a subject of debate recently.
Various numerical methods provide different answers
<cit.>.
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, as we discuss in the following, 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. <ref>, we compare the ground state energies E_GS of different geometries upon varying χ at ϕ/π = 0.1, 0.25, and 0.5.
Narrow cylinders with L_circ = 6 (3xL_y rho) or L_circ=8 (2xL_y rec)
show 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) exhibits 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 corresponding classical model <cit.> appear as meta-stable ground states
if an MPS ansatz with small χ is used,
in particular at ϕ/π = 0.5, see Fig. <ref> (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.


Remarkably, we find NP_αβ to be adiabatically connected to uncoupled K Γ chains <cit.>
along alternating α and β bonds
[Likewise, we find NP_γ to be adiabatically connected to uncoupled K Γ dimers on bonds with label γ. This is consistent with a recent series expansion on a dimer ansatz <cit.>, which finds the dimerized ansatz to be stable up to almost isotropic couplings.].
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. <ref>(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, NP_xz or NP_xy respectively,
is selected at a small but finite η.


The three possible orientations of the NP_αβ are exactly degenerate
in the two-dimensional limit.
This degeneracy, however, is lifted on the cylinder due to its finite circumference,
which is best understood in the limit of the uncoupled K Γ chains.
On the rhombic geometry, see Fig. <ref>(a), NP_yz corresponds to a K Γ chain
composed of alternating 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. <ref>.
Thus, the significant finite-circumference — or finite-size — effects of the two-dimensional K Γ model can be explained,
at least within the NP_αβ set of the NP phase,
by its relation to the K Γ chains. 
We expect this type of finite-size effect to be quite ubiquitous in any numerical study
of an extended spin-1/2 Kitaev model with the Kitaev and other anisotropic exchange interactions 
using a finite geometry
<cit.>.






§ DYNAMICS
 



The dynamical spin structure factor 𝒮(𝐤,ω)
contains information about the excitation spectrum and can be probed by inelastic neutron scattering experiments.
We consider 
𝒮(𝐤,ω) = ∑_γ∈{x,y,z}𝒮^γγ (𝐤,ω)
with
𝒮^γγ (𝐤,ω) being the spatio-temporal Fourier transform
of the dynamical correlations 

    𝒮^γγ (𝐤,ω) = N ∫dt   e^i ω t∑_a,b e^i 𝐤· (𝐫_b - 𝐫_a)  C^γγ_ab(t)  ,

where γ∈x,y,z is the spin component,
𝐫_a and 𝐫_b are the spatial positions of the spins,
N is a normalization factor defined via 
∫ dω∫d𝐤 𝒮^γγ ( 𝐤,ω) = ∫d𝐤,
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 <cit.>
until t = 60.
The time series is then extended by a linear prediction <cit.> 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 L_circ = 6 sites is used such that
three separated lines of accessible momenta exist in the first Brillouin zone,
as illustrated in Fig. <ref>.

We present two sets of 𝒮(𝐤,ω),
Fig. <ref>(c-f) and Fig. <ref>,
along the high-symmetry direction Γ-M-Γ' in the Brillouin zone.
In Fig. <ref>(c-f) we want to emphasize the similarities and differences of
𝒮(𝐤,ω) in the two limits (c,d) Γ⟶ 0
and (e,f) K ⟶ 0.
Whereas in Fig. <ref> we focus on the evolution of the magnon bands
upon reducing the field across the critical field while keeping ϕ/π=0.1 fixed.
𝒮(𝐤,ω) for additional parameters are presented in Appendix <ref>.

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.

On the contrary, the continuum exhibits a very distinct structure in both limits of ϕ.
Near ϕ⟶ 0, that is near the pure Kitaev model with K<0,
𝒮(𝐤,ω) 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
<cit.>.
Since the flux-pair excitation is gapped, 𝒮(𝐤, ω) has a finite gap
even though the Majorana fermion spectrum is gapless.
The continuum above the spectral gap has an almost 𝐤-independent width equal to that
of the single Majorana fermion dispersion.
𝒮(𝐤,ω) 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,
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 <cit.>,
was previously observed in the Kitaev-Heisenberg model <cit.>.

In the other limit ϕ/π⟶ 0.5, i.e. the pure Γ model with Γ>0,
𝒮(𝐤, ω) 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. <ref>(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 <cit.>.
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 <cit.>, 
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 <ref> 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. <ref>(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. <ref>(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. <ref>(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. <ref>(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. <ref>(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. <ref>(a).
In the corresponding classical model, the intermediate ZZ order extends down to
ϕ/π = 0.1 <cit.>.
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 <ref>).







§ 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 broken 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 <cit.>
in the excitation spectrum, as shown in the dynamical spin structure factors of the quantum model.

In realistic Kitaev materials like , the lattice-rotational symmetry
has already been broken by a monoclinic distortion <cit.>.
Therefore, it is reasonable to expect that the degeneracy of NP is lifted,
i.e. one or two out of the six canting directions is favored.
In fact, a substantial magnetic torque is measured for  at high fields <cit.>,
indicating that NP may be stabilized in .
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.

Considering the wide range in parameter space over which NP occurs,
this phase may as well be relevant to other Kitaev materials.
In the context of the NP phase, with its spontaenously broken C_3 symmetry,
a particularly interesting example is the iridate K_2IrO_3.
Recently, K_2IrO_3 was proposed to exhibit ferromagnetic Kitaev exchange and large off-diagonal exchange while maintaining a C_3 point group symmetry <cit.>.

An important 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 <cit.>
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.




    

§ ACKNOWLEDGEMENTS
 
    
    We are grateful to R. Kaneko, K. Penc, J. G. Rau, J. Romhanyi, and N. Shannon for insightful discussions.
    M.G. acknowledges support by the Theory of Quantum Matter Unit of the Okinawa Institute of Science and Technology Graduate University (OIST)
    and by the Scientific Computing section of the Research Support Division at OIST for providing the HPC ressources.
    L.E.C. and Y.B.K. were supported by the Killam Research Fellowship from the Canada Council for the Arts, the NSERC of Canada, and the Center for Quantum Materials at the University of Toronto.
    HYK is supported by NSERC of Canada and the Center for Quantum Materials at the University of Toronto.
    This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.






§ MINIMIZATION OF THE CLASSICAL ENERGY IN THE FERROMAGNETIC PHASE
 
In this section, we work in the crystallographic (abc) coordinates, and write the spin 
components in (<ref>) as 𝐒=(S^a,S^b,S^c). 
For a ferromagnetic order, the Hamiltonian (<ref>) reduces to 
<cit.>

    [b]
    H =     N { K |𝐒|^2 + γ[ - (S^a)^2 - (S^b)^2 + 2 
    (S^c)^2 ] 
       - 2 (h^a S^a + h^b S^b + h^c S^c) } + λ( |𝐒|^2 - 
    S^2 ) ,

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 <cit.>, we can simply drop it 
from (<ref>) in the subsequent analysis. Extremizing H with respect to 
the variables S^a, S^b, S^c and λ leads to the following equations


    ∂ H/∂ S^a = -2 γ S^a - 2 h^a + 2 λ S^a=0, 
    
       ∂ H/∂ S^b = -2 γ S^b - 2 h^b + 2 λ S^b=0, 
    
       ∂ H/∂ S^c = 4 γ S^c - 2 h^c + 2 λ S^c=0, 
    
       ∂ H/∂λ = (S^a)^2+ (S^b)^2 + (S^c)^2 - S^2=0.


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. (<ref>) becomes -2 γ S^a + 2 
λ S^a=0. (<ref>) yields a similar condition. 
 

i. Suppose that S^a ≠ 0 or S^b ≠ 0. Then λ=γ. 
(<ref>) then implies S^c=h^c/3γ, which is a physical solution only when 
h^c/3γ≤ S. 
 

ii. Suppose that S^a=0 and S^b=0. Then, (<ref>) and 
(<ref>) are satisfied. Futhermore, λ≠ -2 γ, otherwise 
(<ref>) would imply h^c=0. Therefore, S^c=h^c/(λ+2 γ), with 
λ chosen to satisfy the normalization (<ref>). 


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.



§ SIGNATURES OF THE ACCIDENTAL U(1) DEGENERACY
 



As shown in Sec. <ref>,
quantum order-by-disorder lifts the accidental U(1) ground state degeneracy
in the classical model (<ref>).
In order to see whether this happens in the quantum model using iDMRG, we apply a magnetic field

    𝐡 = h ( sinϑcosφ 𝐚̂ + sinϑsinφ 𝐛̂ + cosϑ 𝐜̂)

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

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 ground state degeneracy.
In NP the C_3 symmetry is broken spontaneously,
which is apparent from the anisotropic bond energies and the canting of the induced magnetic moment.
As a consequence, 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. <ref>.
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.



§ ADDITIONAL PLOTS OF DYNAMICAL SPIN STRUCTURE FACTOR
 



We complement the discussion on dynamical spin structure factors 𝒮(𝐤,ω) by presenting the data for additional sets of parameters, as shown in Fig. <ref>.
We restrict ourselves to the phases close to the transition between P and NP,
where the numerical errors and finite size effects are found to be small.
Here we employ the rhombic geometry with L_circ = 6 sites
that has three separated lines of accessible momenta in the first Brillouin zone,
as shown in Fig. <ref> (a,b).

The P phase clearly exhibits two magnon bands, where the lower band attains its minimum at the Brillouin zone center, i.e. the Γ and Γ' points, for any ϕ/π∈(0,0.5]. Upon approaching the transition, the lower band touches zero frequency at Γ and Γ', which is consistent with the uniform canting of spins in the NP phase.

Regardless of the specific value of ϕ, 𝒮(𝐤,ω) 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 low-energy excitations with long wavelengths.
Such excitations correspond to the pseudo-Goldstone mode of the slightly lifted U(1) degenerate manifold of states.
Moreover, 𝒮(𝐤,ω) 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. <ref> 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 ℤ_2 fluxes <cit.>.
The fluxes are gapped excitations of an emergent ℤ_2 gauge field.
𝒮(𝐤,ω) 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 <cit.>.
The continuum is almost entirely formed by Majorana fermions subject to a quantum quench 
that caused by exciting a flux pair <cit.>.

In the NP or P phase proximate to the KSL, see Fig. <ref>(c,d) for ϕ/π=0.01 and Fig. <ref> 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, 𝒮(𝐤,ω) 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|,
𝒮(𝐤,ω) 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 J K ΓΓ' in a [111] field
mapping the K Γ model at ϕ = tan^-1 (1/2) ≈ 0.148π
to the pure FM Heisenberg model at high field <cit.>.
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 <cit.>.
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 longitudinally-polarized 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. <ref>(e,f) and Fig. <ref>.
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 <cit.>.
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.



§ 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. [chern_magnetic_2020] 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. <ref>-<ref>. This is also seen in the DSF of the quantum model, as discussed in Sec. <ref> and Appendix <ref>.
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.






