Unconventional magnetic field response of the hyperhoneycomb Kitaev magnet β −

We present a unified description of the response of the hyperhoneycomb Kitaev magnet β-Li2IrO3 to applied magnetic fields along the orthorhombic directions a, b and c. This description is based on the minimal nearest-neighbor J-K-Γ model and builds on the idea that the incommensurate counter-rotating order observed experimentally at zero field can be treated as a long-distance twisting of a nearby commensurate order with six spin sublattices. The results reveal that the behavior of the system for H ‖ a, H ‖ b and H ‖ c share a number of qualitative features, including: i) a strong intertwining of the modulated, counter-rotating order with a set of uniform orders; ii) the disappearance of the modulated order at a critical field H∗, whose value is strongly anisotropic with H∗ b <H ∗ c H∗ a; iii) the presence of a robust zigzag phase above H∗; and iv) the fulfillment of the Bragg peak intensity sum rule. It is noteworthy that the disappearance of the modulated order for H ‖ c proceeds via a ‘metamagnetic’ first-order transition which does not restore all broken symmetries. This implies the existence of a second finite-T phase transition at higher magnetic fields. We also demonstrate that quantum fluctuations give rise to a significant reduction of the local moments for all directions of the field. The results for the total magnetization for H ‖ b are consistent with available data and confirm a previous assertion that the system is very close to the highly-frustrated K-Γ line in parameter space. Our predictions for the magnetic response for fields along a and c await experimental verification.

Besides the dominant Kitaev anisotropy, the above materials feature additional weaker interactions which generally give rise to a wealth of nontrivial phases competing with the quantum spin liquids [1][2][3][4][5][6][7][8]. A central goal in the field is therefore to map out the various instabilities and identify the distinctive experimental signatures of the most relevant interactions. One of the most promising ways to achieve this goal experimentally is to subject the Kitaev materials in external magnetic fields along different directions. Apart from controlling the interplay of various zero-field competing phases, such external fields can also stabilize new collective spin states. There are, for example, many reports for possible magnetic field-induced quantum spin liquids [37][38][39][40], and a variety of complex multisublattice, single-and multi-Q phases [41][42][43][44][45].
Remarkably, all experimental data reported so far for Kitaev materials show that their response to the magnetic field depends very strongly on its direction. This is true for the layered compounds Na 2 IrO 3 [10], α-Li 2 IrO 3 [46], and α-RuCl 3 [29][30][31]47], as well as for the three-dimensional (3D) iridates β-Li 2 IrO 3 [21,24] and γ-Li 2 IrO 3 [19,48]. Here we revisit the case of the hyper-honeycomb β-Li 2 IrO 3 and show that its strongly anisotropic response signifies a large separation of energy scales between the relevant microscopic interactions, and can thus be used to extract information about the relative strength of these interactions in a direct way.
The main features of β-Li 2 IrO 3 that are known so far are as follows [14,16,21,22,24]. At zero field, the system orders magnetically below T N = 38 K, with the spins forming a noncoplanar, incommensurate (IC) modulation, with propagation wavevector Q = (0.57, 0, 0) in the orthorhombic frame, and two counter-rotating sets of moments [14], similar to those in γ-Li 2 IrO 3 [15] and α-Li 2 IrO 3 [18]. A magnetic field along b destroys the IC order at a characteristic field H * b ∼ 2.8 T, beyond which the spins show a uniform Q = 0 coplanar phase, comprising a ferromagnetic (FM) component along the field and a robust zigzag component along a [21]. These components are also present below H * b , but are too small to be detected at zero field [49,50]. For H a and H c, the system shows a much weaker response, with the IC order remaining robust and the magnetization being linear up to the maximum fields measured (see supplemental material in [21] and also [24]).
On the theory side, it has been established that the magnetism of β-Li 2 IrO 3 can be accurately described by the nearest neighbor (NN) J-K-Γ model [49][50][51][52], where K denotes the Kitaev coupling, J the Heisenberg coupling and Γ the so-called symmetric exchange anisotropy which is present in many Kitaev materials [51][52][53][54][55][56]. In particular, β-Li 2 IrO 3 is believed to be in the regime of large negative K, large negative Γ (with |Γ| < |K|) and small positive J (with J |Γ|), see detailed discussion in [49,50]. Remarkably, in this parameter regime, the critical field H * b depends only on J, specifically [50] µ B H * b ∼ 0.46J (4S /g bb ), where S = 1/2 denotes the classical spin length of the j eff = 1/2 degree of freedom, µ B is the Bohr magneton and g bb is the diagonal element of the electronic gtensor along b. The small value of the experimentally measured H * b is therefore a signature of the smallness of J (J ∼ 4 K).
It has also been shown [49,50] that the IC order of β-Li 2 IrO 3 can be treated as a long-distance twisting of a nearby commensurate period-3 state with Q = 2 3â (in units 2π a ). This state is amenable to a semi-analytical treatment of the problem, with results that are consistent with almost all experimental findings so far, both in zero field and at finite fields along b [49,50]. This analysis explains, for example, the presence of a uniform zigzag component along a on top of the modulated order, and the intensity sum rule of the corresponding Bragg peaks [21].
Here we show that this semi-analytical description can be naturally extended to the cases where H is along a and c. The results, which are cross-checked with classical Monte Carlo simulations, show that the response along a and c directions shares many qualitative features with that along b. Specifically, we find that the period-3 order disappears at a critical field H * , whose value depends strongly on the field direction. Importantly, none of the critical fields depends on the Kitaev interaction K, and moreover H * a are mainly controlled by Γ. A realistic set of coupling parameters, J = 0.4 meV, K = −18 meV and Γ = −10 meV (1) delivers H * b ∼ 2.88 T and T N ∼ 35.5 K (in good agreement with corresponding experimental values of 2.8 T and 38 K [14,21]), and also gives H * a ∼ 102 T and H * c ∼ 13 T. This means that at least the transition at H * c should be accessible experimentally, and the measured value of H * c can provide the value of Γ. The same semi-analytical approach provides a number of additional qualitative findings: (i) The period-3 order is always intertwined with a set of uniform orders, some of which give rise to a finite torque that can be measured experimentally.
(ii) Among these uniform orders, there is always a zigzag component which remains robust above H * and coexists with the FM order along the field. Classically, the zigzag component disappears (but only for g ab = 0, see below) at H * * → ∞ for fields along a and b, but for H c the corresponding field H * * c is finite. In particular, H * * c is governed mostly by Γ, with H * * c ∼ 45 T for the parameters of Eq. (1). (iii) The intensity sum rule between the Bragg peaks at Q = 2 3â and Q = 0, which has been observed experimentally for H b [21], is actually fulfilled for all field directions. As it turns out, this rule is an experimental fingerprint of the spin length constraints. (iv) While the transitions at H * a and H * b are continuous, the transition at H * c is of first order. Moreover, this transition does not restore all broken symmetries, which leads to the prediction of a second thermal phase transition at high enough fields along c. This transition will be demonstrated explicitly by classical Monte Carlo simulations.
One shortcoming of our semi-analytical classical approach is that it overestimates the magnetization at H * b by approximately a factor of two compared to the experimental value. This has led to the assertion [50] that the spin lengths are strongly renormalized by quantum fluctuations due to the close proximity to the special K-Γ line in parameter space, where the system is highly-frustrated [49]. This assertion is now demonstrated explicitly by a semiclassical 1/S expansion. The results confirm that the magnetization correction can be as large as 50%, and a direct comparison with published experimental data shows good agreement for all three field directions.
The rest of the paper is organized as follows. In Sec. II, we recall the structural and symmetry aspects of β-Li 2 IrO 3 that are most relevant for this study. In Sec. III, we review the minimal J-K-Γ model [51,52]  spin Hamiltonian. In Sec. IV, we present the unified semianalytical description for all three orthorhombic directions. This includes the basic spin-sublattice structure of the various configurations (Sec. IV A), their parametrization in terms of Cartesian components and the associated symmetry-resolved static structure factors (Secs. IV B-IV C), the dependence of the critical fields H * on the model parameters (Sec. IV E), and a discussion of the symmetries that are broken in each regime (Sec. IV F). The role of quantum fluctuations is then addressed in Sec. V, along with the direct comparison of the predicted magnetizations with available experimental data. In Sec. VI we discuss how the various transitions can be detected experimentally via measurements of the magnetic torque, for which we provide predictions with and without harmonic spinwave corrections. In Sec. VII, we cross-check our ansätze with classical Monte Carlo simulations and compute the H-T phase diagram for all three orthorhombic directions. Here we also highlight the qualitative difference between the zigzag orders for H a and H b versus the spontaneous high-field zigzag order for H c. A summary and a general discussion is given in Sec. VIII. Auxiliary information and technical details are provided in Appendices A-E.
II. Lattice structure, symmetries & conventions β-Li 2 IrO 3 crystallizes in a hyperhoneycomb structure (shown in Fig. 1) and has the Fddd space group. Its conventional orthorhombic unit cell is set by the crystallographic axes {â,b,ĉ}, which are related to the Cartesian axes {x,ŷ,ẑ} appearing in the spin Hamiltonian below [Eqs. (4)- (5)] bŷ We note here that we stick to the xyz-frame convention of Refs. [49][50][51][52], which is different from the one used in Ref. [14]. The two xyz-frames are related to each other by a two-fold rotation around the x-axis. This is important as the choice of the frame affects the overall sign structure of the Γ interactions.
The orthorhombic unit cell contains four primitive unit cells, of four Ir 4+ ions each (labeled by Ir 1 -Ir 4 in Fig. 1). The Ir 4+ ions form a hyperhoneycomb structure, which can be viewed as a stacking of two types of zigzag chains, which we will denote by xy-and x y -chains. The xy-chains run along the direction a+b and are shown in Fig. 1 by the alternating red and green solid bonds, denoted by x and y respectively. The x y -chains run along a−b and are shown in Fig. 1 by the alternating red and green dashed bonds, denoted by x and y respectively. The two types of chains are interconnected with vertical NN Ir-Ir bonds denoted in Fig. 1 by z (blue solid lines). In total, there are five types of NN Ir-Ir bonds, x, y, x , y and z.
Apart from translations, the crystal structure is invariant under the following point group operations [21]: (i) Inversion I through the center of every x-or yor x -or y -type of bond, such as the center of the Ir 2 -Ir 4 bond of Fig. 1. (ii) Three π-rotations in combined spin-orbit space, C 2a , C 2b , and C 2c , around the axes a, b and c, respectively, passing through the middle of the z bonds, as shown in Fig. 1. In particular, C 2a maps x-bonds to y -bonds and y-bonds to x -bonds in real space, and [S x , S y , Similarly, C 2b maps x-bonds to x -bonds and y-bonds to ybonds in real space, and [S x , S y , S z ] → [−S x , −S y , S z ] in spin space. Finally, C 2c maps x-bonds to y-bonds and x -bonds to y -bonds in real space, and [S x , S y , S z ] → [S y , S x , −S z ] in spin space. (iii) Three glide planes which arise by reflections across the ab-, bcand ac-planes passing through an inversion center, followed by nonprimitive translations by ( 1 and ( 1 4 0 1 4 ), in orthorhombic units, respectively. At this point it is also worth introducing some terminology that we will need later in the analysis of the static structure factors. Following Ref. [14], we define four-component symmetry basis vectors, These vectors represent, respectively, the relative amplitudes of the four sites of the primitive unit cell in the Néel (A), stripy (C), ferromagnetic (F) and zigzag (G) order. Note that, for consistency, our 4-site labeling Ir 1 -Ir 4 of Fig. 1 follows the convention of Fig. 7 of Ref. [14].
For the various components of the static structure factor, we follow the convention of Ref. [49] and denote the modulated components with Q = 2 3â by the letter M and the uniform components with Q = 0 by M . Therefore, M a (A) denotes the modulated Néel (A) component along a, M b (F) denotes the uniform ferromagnetic (F) component along b, and so on. The definitions of these components in terms of the Fourier transform of the spin configuration are given in Appendix A 1.
III. The minimal J-K-Γ model Following earlier works [49][50][51][52], we consider here the minimal microscopic J-K-Γ model mentioned above, supplemented with a Zeeman term H Z to describe the coupling to the external field H. The total Hamiltonian then reads where Here S i denotes the pseudo-spin j eff = 1/2 operator at site i, t ∈ {x, y, z, x , y } labels the five different types of NN Ir-Ir bonds and (α t , β t , γ t ) = (x, y, z), (y, z, x), and (z, x, y) for t ∈ {x, x }, {y, y }, and {z}, respectively. The prefactor σ t equals +1 for t ∈ {x, y , z} and −1 for t ∈ {y, x }, see ± symbols in Figs. 1 and 2. This overall sign structure of the Γ interactions derives from the symmetries mentioned above [51] and our choice of the xyz-frame in Eq. (2). Finally, g i stands for the g-tensor of the i-th Ir ion. As discussed by Ruiz et al [21], these tensors carry a site-dependent, staggered off-diagonal element g ab . Specifically, in the orthorhombic frame, where p i = +1 for spins on the xy chains and −1 for spins on the x y chains. Here we take g aa = g bb = g cc = 2 and g ab = 0.1. We note that, depending on the direction of the field, some of the discrete symmetries mentioned above may or may not be preserved, see Table II. A field along a, for example, breaks both C 2b and C 2c , but still respects C 2a , ΘC 2b and ΘC 2c , where Θ is the time reversal operation. In the following we restrict ourselves to the so-called Kregion of the parameter space with dominant Kitaev interaction, which is believed to be relevant for β-Li 2 IrO 3 [49], and fix the parameters to the representative set given in Eq. (1).
IV. Unified description of β-Li 2 IrO 3 for H along a, b and c axes A. General spin sublattice structure The behavior of β-Li 2 IrO 3 under a magnetic field along the three orthorhombic directions can be described in a unified manner as shown in Fig. 2. For all three directions, a, b and c, the system goes through a low-field phase (0 ≤ H < H * ) with six spin sublattices [A, B and C along the xy-chains, and A , B and C along the x y -chains, see Fig. 2 (a)], followed by a high-field canted phase (H * < H < H * * ) with two spin sublattices [F along the xy-chains and F along the x y chains, see Fig. 2 (b)]. The high-field phase terminates at H * * = ∞ for H a and H b (with a small zigzag component remaining if g ab 0, see Appendix B), whereas H * * c is finite and the classical state reached at H * * c is the fully polarized state. Figure 3 shows a series of representative snapshots, of various ground state configurations for different field directions and strengths (obtained from numerical minimization of the classical ansätze discussed below). As discussed in Ref. [49],  Table I. in the zero-field state, the three sublattices A, B and C along the xy chains form a nearly coplanar 120 • state, and the three sublattices A , B and C along the x y chains form another such nearly 120 • structure, on a different plane, see dotted blue triangles at the top left panel of Fig. 3. Under a magnetic field, the three sublattices of each given chain cant toward each other and eventually get aligned at the characteristic field H * where A = B = C ≡ F and A = B = C ≡ F . For fields along a and b, this intra-chain alignment happens continuously, whereas for fields along c it happens abruptly. Above H * , F and F cant toward the field in a non-uniform way and at a pace that is strongly dependent on the field direction.

B. Basic characterization of the low-field phase (H < H * )
The individual Cartesian spin components of the various configurations are related to each other in a specific way, see the parametrization in Table I. For each given spin sublattice, a spin length constraint must be imposed, for example x 2 1 + y 2 1 + z 2 1 = 1 for the A sublattice, 2x 2 3 + z 2 3 = 1 for the C sublattice of the H a case, etc. The field dependence of the Cartesian components can be obtained by a numerical minimization of the total energy of the system [see Eqs. The low-field phase for H b is described by five Cartesian components (x 1 , y 1 , z 1 , x 2 and z 2 ) or, equivalently, by five structure factors [50]: Three modulated Q = 2 3â components M a (A), M b (C) and M c (F), and two uniform Q = 0 components M a (G) and M b (F). This precise combination is in fact present for all three orthorhombic directions for 0 ≤ H < H * , as it is a property of the zero-field state. In particular, the uniform components M a (G) and M b (F) of the zero-field order reflect the deviation from the perfect 120 • coplanar order mentioned above, see detailed analysis in Ref. [49]. Note further that the modulated components M a (A), M b (C) and M c (F) belong to the Γ 4 irreducible representation, in agreement with experiment [14]. The five structure factors satisfy two constrains which, when normalized appropriately [see Appendix A 1], can be combined to give the Bragg peak intensity sum rule observed experimentally [21]. Namely, where Turning to the low-field phase for H a, here we have ten Cartesian components (see Table I) or, equivalently, ten  structure factors: six modulated components (the three zerofield components plus three induced by the field) and four uniform components (the two zero-field components plus two  Table II of Ref. [14]). Altogether, the ten structure factors satisfy four constraints, and one combination of them gives the Bragg peak intensity sum rule of Eq. (7), where now The low-field phase for H c is described by nine Cartesian components (see Table I) or by nine structure factors: Here the field induces three modulated components (M a (G), M b (F) and M c (C)) and one uniform component M c (F) (which couples directly to the Zeeman field). The modulated components belong to the irreducible representation Γ 3 (see Table II of Ref. [14]), and, as it turns out, they remain at least one order of magnitude smaller than the dominant Γ 4 components, see Fig. 5 (b). Altogether, the nine structure factors satisfy three constraints, and one combination of them gives the Bragg peak intensity sum rule of Eq. (7), where now Let us emphasize that the fulfilment of the intensity sum rule Eq. (7) for all field directions and strengths is a direct fingerprint of the local spin length constraints. The numerical prefactor of 2 in the definition I tot = 2I I +I V reflects the fact that there are twice as many Bragg peaks characterizing the modulated order (Q = ± 2 3â ) compared to the peaks characterizing the uniform order (Q = 0), see detailed analysis and a general proof of Eq. (7) in Appendix A 2.
Note finally that some of the uniform components generated for H along a and c give rise to a finite magnetic torque signal, which will be examined separately in Sec. VI.  Table I. Note also that, for H ≥ H * , the spins lie on the ab-plane for H a and H b, but for H c the spin plane changes continuously. This is related to the fact that the uniform components of the zero-field state all lie in the ab-plane, and so a field applied in this plane will merely reorganize these components and not rotate them out of the plane, unlike what happens for H c.
The zigzag component disappears at a characteristic field H * * . As mentioned above, H * * is infinite for H along a and b but finite for H c, with (see Appendix B 3) which, for J |Γ|, reduces to According to this relation, H * * c depends mostly on Γ, with H * * c ∼ 45 T for the coupling parameters of Eq. (1).

D. Robustness of high-field zigzag orders
We now discuss why the various high-field zigzag orders remain robust up to very high fields, for all three orthorhombic directions. The most direct way to see this is to express the total energies E b , E a and E c in terms of the various static structure factors, see Eqs. (B11), (B23) and (B35), respectively. It turns out that E b and E c contain an explicit cross-coupling term between M a (G) and while E a contains an explicit cross-coupling term between M a (F) and The presence of these terms reveal that the qualitative reason why it is energetically favorable for the system to sustain appreciable zigzag orders up to high fields is the strong Γ interaction. Of course, the actual quantitative details for each field direction derive from the minimization of the total energies under the given constraints. For example, the analytical expression Eq. (13) for H * * c can be derived by minimizing

E. Dependence of H * on microscopic coupling parameters
The characteristic field H * marks the disappearance of the modulated components (and M a (G) and M b (F) for H a). As mentioned earlier, this transition is continuous for H a and H b, but of first-order for H c, see Figs. 4 (a-f). Furthermore, the value of H * depends strongly on the direction of the field. For the coupling parameters of Eq.
This large difference between the critical fields along different directions is related to the strongly anisotropic character of the Hamiltonian, and the different role of the various couplings in each case. For example, as we discussed in Ref. [50], in the parameter regime of interest, H * b depends only on J, which is why H * b is very small. We will now show that H * a and H * c do not depend on K but only on J and Γ, and that the inequality J |Γ| explains why these critical fields are larger compared to H * b . To this end, we will vary the parameters of the model and take a closer look at the evolution of the various contributions to the total energy with the field. Figure 6 (a-c) shows the field-driven evolution of E J , E K , E Γ and E Z , which denote the contributions from J, K and Γ interactions and the Zeeman energy, respectively. The corresponding derivatives of these energies with respect to H are shown in Fig. 6 (d-f). The main finding is that, in the parameter regime of interest, E K remains almost insensitive to H, and this is true for all field directions. This means that the Zeeman field does not act against K, which explains why none of the critical fields H * depends on the dominant coupling of the theory. The results also show that, unlike H * a and H * c , the critical field H * b depends only on J and not on Γ; this is the consequence of the fact that E Γ does not change with H in this direction. These arguments can be formulated mathematically by the following relations that arise from a classical version of Feynman-Hellmann theorem (see Appendix C):  where N is the total number of spins and m (J, K, Γ, H) is the magnetization per site along the field. According to these relations, the fact that ∂E K /∂H ≈ 0 implies that ∂m /∂K ≈ 0, i.e., that the whole magnetization process does not depend on K. Likewise, the fact that ∂E Γ /∂H ≈ 0 for H b implies that ∂m /∂Γ ≈ 0, and therefore the whole magnetization process depends only on J in this field direction. We can go one step further and extract the actual dependence of the critical fields on the relevant couplings by computing these fields for a wider range of parameters. The results are shown in Fig. 7 and demonstrate that the critical fields H * a and H * c depend almost perfectly linearly on J and Γ. Fitting the numerical data for H * gives, in particular, Thus, besides the independence of H * b on K and Γ, [57], we find that H * a is controlled mainly by Γ (given that J |Γ|), whereas H * c is controlled by both J and Γ. Note that the coefficients appearing in Eqs. (18) correspond to the value g ab = 0.1 whose sign and magnitude is chosen arbitrarily here. However, the coefficients do not depend much on this choice. For example, for F. Symmetries Table II shows the symmetry properties of the various fieldinduced configurations for different field directions. The primitive translations (denoted by T ) are broken spontaneously in the low-field phase (0 < H < H * ) due to the modulating components of the order. This symmetry is restored above H * with the disappearance of these components. Furthermore, the lowfield phases preserve the inversion symmetries I around the centers of the FM dimers AA, BB, A A or B B of Fig. 2 (a), while the high-field phases preserve the inversion centers on all x, y, x and y bonds.
Let us now turn to the C 2 -rotation symmetries discussed in Sec. II or their combinations with time reversal Θ. For H b, the symmetries ΘC 2a , C 2b and ΘC 2c of the model are all preserved in both the low-and the high-field phases, emphasizing once again the special role of the b axis [14,21,50].
For H a, on the other hand, among the three symmetries C 2a , ΘC 2b and ΘC 2c , the first two are broken spontaneously in the low-field phase due to M a (G) and M b (F). This symmetry breaking is associated with the choice of the overall sign of M a (G) and M b (F). One can see this more directly from the cross-coupling term of Eq. (B18), according to which the relative signs of M a (G) and M b (F) are fixed by the sign of Γ, but one can still change both signs at the same time without changing the energy. Note that, while a similar cross-coupling term appears between M a (F) and M b (G), see Eq. (16), the individual signs of these two components are fixed by the Zeeman field which couples directly to M a (F), see Appendix B 2. The symmetries C 2a and ΘC 2b are restored at H * a with the disappearance of the M a (G) and M b (F).
The situation for H c has one qualitative difference (besides the abrupt transition at H * c ). Here, among the three symmetries ΘC 2a , ΘC 2b and C 2c of the model, the last two are broken spontaneously in both the low-and the high-field phase, and only get restored at H ≥ H * * c . The symmetry breaking occurs again due to M a (G) and M b (F), which couple via Eq. (15). Hamiltonian As above then, Γ fixes the relative signs of M a (G) and M b (F), but the overall choice of the global sign remains arbitrary.
Altogether, unlike what happens along a, the transition at H * c does not restore all broken symmetries, and one thus expects a second thermal phase transition at high fields, even after the disappearance of the modulated order. This will be shown explicitly in Sec. VII. For completeness, let us recall that the zero-field state breaks C 2a and C 2c , but respects ΘC 2a , ΘC 2c and C 2b [49].

V. Magnetization process & the effect of quantum fluctuations
We now focus on the magnetization per site m, defined as Here N m is the number of spins inside the magnetic unit cell (N m = 48 for H < H * and N m = 2 for H > H * ), µ = 1-N m , S µ is the expectation value of the spin on the µ-th sublattice, and g diag , g off−diag and p µ are defined in Eq. (6). Recalling that p µ = +1 for spins along the xy chains and −1 for spins along the x y chains [see Fig. 1 The magnetizations along the field m (denoted by m a , m b and m c for H along a, b and c, respectively) are given by In agreement with experiment, m b rises much faster than m a and m c . Furthermore, the magnetizations m a and m b first increase monotonously with the field, then show a kink at H * a and H * b , respectively, and then increase at a much slower pace towards a limiting value that is determined by the ratios g ab /g aa and g ab /g bb , respectively (see Appendix B). By contrast, m c shows a finite jump (instead of a kink) at H * c , reflecting the corresponding jumps in Fig. 4 (f). At higher fields, m c shows a kink at H * * c and then saturates. Note that here the exact saturation is only true for classical spins, and the kink in the classical magnetization will be smoothed by quantum fluctuations (as the spin Hamiltonian does not conserve rotations around the field axis, and the fully polarized state is not a true eigenstate).
Let us now compare these classical predictions for m with available experimental data published by Ruiz et al (see supplementing material in [21]), which are shown in Fig. 8 by black lines. Quite generally, while the classical ansätze capture the observed magnetization processes qualitatively, there is a large quantitative discrepancy. For H b, for example, the classical prediction for the magnetization at H * b is about two times larger than the measured value. This deficiency has been recognized previously [50], and has led to the assertion that the system must feature strong quantum fluctuations due to the close proximity to the highly-frustrated K-Γ line [49].
Here we confirm this hypothesis by calculating the leading 1/S corrections to the magnetization from quantum fluctuations. The details of this calculation are provided in Ap- pendix D and the renormalized magnetization curves are shown by the solid blue lines in Fig. 8. The results show that already the leading 1/S corrections reduce the magnetization quite strongly, bringing the curves much closer to the measured data. While subleading higher-order corrections will reduce the magnetization even further, providing a better comparison between theory and experiment, a final quantitative agreement will also require an appropriate re-adjustment of the microscopic couplings. Importantly, our semiclassical results show further that, for H c, the magnitude of the magnetization jump at H * c is significantly reduced by quantum fluctuations, almost to the point that there is no visible change, including the overall slopes of the curves below and above H * c . This renders the detection of this feature in magnetization measurements more challenging and probably explains the absence of the kink in recent measurements [24]. The detection is even more challenging for powder samples given that m c m b . Nevertheless, as we will discuss next [Sec. VI], the transition at H * c should be still visible via the kink in the corresponding magnetic torque.
Interestingly, the expression for ξ that gives the transverse components for H a and H c coincides with the expression for m b , see Eq. (20). Note that, as we discussed in Sec. IV F, the overall signs of M a (G) and M b (F) are chosen spontaneously by the system for both H a and H c, and therefore the sign of the torque (or ξ) is arbitrary for both directions. This aspect has further observable consequences, which will be discussed in Sec. VIII. Figure 9 (b) shows the evolution of τ/H with H for fields along a and c, with and without harmonic spin-wave corrections. First of all, the torque for H a is about 40 times weaker than the torque for H c. This reflects the smallness of M a (G) and M b (F) components for H a, as shown in Fig. 5 (a). Second, the torque for H a remains non-zero up to H * a , whereas the torque for H c remains non-zero up to H * * c . This again stems from the associated behaviors of M a (G) and for H c, in particular, shows a characteristic sharp kink at H * c , reflecting the first order transition between the low-field sixsublattice and the high-field two-sublattice state. Importantly, this kink remains sharp even after we include the leading 1/S spin-wave corrections (blue line, see Appendix D). A measurement of the torque can therefore give direct evidence for the transition at H * c , and thus provide information for the value of Γ via Eq. (18).
Finally, for H c, the torque in the high-field phase scales as Thus, a measurement of the torque at high fields can also be used to extract H * * c and, in turn, an independent constraint on the microscopic parameters J and Γ via Eq. (13).

VII. Effect of thermal fluctuations & classical H-T phase diagram
To cross-check the above zero-temperature results from the classical ansätze and confirm the high-field thermal transition for H c mentioned above, we have performed classical Monte Carlo simulations using the standard Metropolis algorithm combined with the over-relaxation algorithm [58,59]. The simulations were performed on finite-size clusters with a total number of sites N ∈ {48, 96, 144, 192, 240, 288} and periodic boundary conditions. All considered systems, spanned by the unit vectors of the orthorhombic lattice, have at least three periods in the orthorhombic a-direction in order to accommodate Q = 2â/3 order, see more details in Appendix E. The results obtained by a thermal annealing down to T = 5 K show that the total magnetization is almost indistinguishable from the predictions of the semi-analytical approach, lending strong support that the latter delivers quantitatively accurate results for the local physics of the problem.
Let us now turn to the classical H-T phase diagrams, which are shown in Fig. 10 for the three orthorhombic directions. The boundary lines of the counter-rotating order (denoted by 'IC') have been extracted by a finite-size analysis of the so-called Binder cumulant [60] (see Appendix E), of the equally-weighted combination of the three modulated static structure factor components (for all field directions): For H b, the phase diagram contains two distinct phases, the high-T paramagnetic phase and the low-T counter-rotating order, which persists up to H * b ∼ 2.8 T, see Fig. 10 (b). For H a, the counter-rotating order persists up to very high fields (H * a ∼ 102 T), see Fig. 10 (a), and is accompanied by the uniform orders M a (G) and M b (F) which are however extremely weak, see Fig. 5 (a). While these orders onset at the same field H * a as the modulated order at T = 0, it is unclear whether this remains true for finite T . In fact, symmetry considerations alone tell us that the boundaries of the two types of orders can in general be different, as the they break different symmetries (the modulated order breaks translations whereas the uniform orders break C 2a and ΘC 2b , see Table II). Unfortunately, the smallness of M a (G) and M b (F) does not allow for an accurate numerical determination of their transition temperature line.
For H c, there are three distinct phases, see Fig. 10 (c). Apart from the paramagnetic and the modulated phase, there is a robust high-field order associated with M a (G) and M b (F), and the spontaneous breaking of C 2c and ΘC 2b (see Table II). This phase coexists with the modulated order at low H and T , but extends up to very high fields (H * * c ∼ 45 T). Its boundary line has been extracted from the Binder cumulant B O Q=0 associated with Finally, the yellow shading in Figs. 10 (a) and (b) represents the variation of the magnitude of |M b (G)| and |M a (G)|, respectively, from high values (intense yellow) at low T to vanishing values (blue) at higher T .

VIII. Discussion
The study presented here provides a semi-analytical framework for the anisotropic response of β-Li 2 IrO 3 under a magnetic field along the three orthorhombic directions. This framework is based on the minimal nearest-neighbor J-K-Γ model [49][50][51][52] and the hypothesis that the local correlations of the low-field incommensurate order can be captured by its closest commensurate approximant with the right symmetry [49]. The results are in qualitative agreement with almost all experimental facts collected so far, and we have shown how a quantitative agreement can also be reached by including quantum fluctuations.
In addition, our analysis delivers a number of predictions which await experimental verification. First, the critical fields H * that mark the disappearance of the modulated order are highly anisotropic, in particular, Such an anisotropic response, which is also evidenced in susceptibility [21,24], signifies a large separation of energy scales between J and Γ. An explicit dependence of H * on these interactions is derived in this work [Eq. (18)] and can be used to extract the actual strength of Γ (the value of J is estimated ∼ 4 K from the value of H * b [50]). Importantly, the dominant Kitaev coupling K does not affect any of the critical fields, partly because it is ferromagnetic.
Second, for all orthorhombic directions our analysis reveals the presence of various intertwined uniform zigzag and FM orders, some of which remain robust far above H * . The physical origin of this robustness is related to the cross-coupling terms of Eqs. (15)(16). Some of the uniform orders give rise to a finite torque signal and can thus be detected in a direct way. Alternatively, they can also be observed by magnetic X-ray diffraction [21], or by local probes like NMR or µSR.
Third, we have shown that the high-field response for H c is special, in that the disappearance of the modulated order at H * c restores only the translational symmetry and leaves some of the discrete symmetries broken. This implies the presence of a second thermal transition above H * c , which is associated with the onset of the uniform orders M a (G) and M b (F). This transition can then be detected with thermodynamic measurements at high enough fields.
A natural extension of the present study is the investigation of the field-induced behavior of β-Li 2 IrO 3 for general field directions, i.e., away from the orthorhombic axes. As it turns out, a semi-analytical description can be also obtained for fields in the ab-and bc-planes [61]. The emerging picture reveals a remarkable interplay of the various modulated and uniform orders and rich anisotropic phase diagrams, the details of which will be given elsewhere [61]. We can, however, comment on one particular aspect related to the torque signal discussed in Sec. VI. As mentioned there, the predicted torque signals are proportional to the quantity ξ = g bb M b (F) + g ab M a (G), whose sign is chosen spontaneously by the system for H a or H c. However, adding an infinitesimal field along b will actually fix the sign of ξ, since the two are directly coupled to each other. This simple argument shows that, as a function of the angle in the ab-or bc-planes, the torque will show an abrupt reversal when the field passes through the a and c axes, respectively. Such a first-order transition scenario could also be relevant for the explanation of the sawtooth-like torque anomalies observed experimentally in the closely related compound γ-Li 2 IrO 3 [48,62] (see also [63]).
Finally, we would like to touch upon an aspect that may be relevant for the interpretation of the phase transition reported recently around 100 K [64]. As discussed in Ref. [49], the zerofield and zero-temperature configuration contains the uniform orders M a (G) and M b (F), in addition to the modulated order. Given that the two types of order break different symmetries (the modulated order breaks translations whereas the uniform orders break C 2a and C 2c [49]) one generally expects that the two types of order onset at different temperatures. In particular, we have checked numerically (unpublished) that the modulated period-3 six-sublattice order carries a pseudo-Goldstone low-energy mode, similar to other incommensurate phases in related models [65][66][67]. On the other hand, the energy barrier associated with flipping the signs of the uniform orders M a (G) and M b (F) gives rise to a finite energy gap. It is then plausible that the uniform orders onset at a higher temperature T uni compared to T N . While the smallness of the uniform orders does not allow to check this numerically with Monte Carlo, the cross-coupling term of Eq. (15) suggests that T uni could scale with Γ. In such a scenario, a field along the b axis will turn the zero-field line extending from T = 0 up to T = T uni into a line of first order transitions, because the field couples directly to M b (F) (and to M a (G) via g ab ). For very low fields, the proximity to this first-order line would then give rise to hysteresis effects, similar to those observed in Ref. [64]. The actual details of this scenario (in particular, the connection of the measured torque signals with the ones we report here at zero-temperature), remain to be explored. Each orthorhombic unit cell contains four primitive cells (labeled by i = 1-4), and each primitive cell contains four spin sites (labeled by ν = 1-4). Each site can then be labeled by the position R of the orthorhombic unit cell, the position ρ i of the primitive unit cell (relative to R), and the position p ν of the spin sublattice (relative to ρ i ). The physical position of each site can then be written as The Fourier transform of the ν-th spin sublattice is defined as where N is the total number of spins and Q belongs to the reciprocal space of the orthorhombic Bravais lattice. The modulated Q = 2â/3 components of the static structure factor are defined as [see Eq.
Note that the extra prefactors of i in the definitions of M(A), M(C) and M(G) have been inserted to follow the convention of Ref. [14], while the normalization prefactor 1/4 in the right hand side of Eq. (A3) sets the maximum possible magnitude of the various components to S . Similarly, the uniform Q = 0 components of the static structure factor are defined as

Local spin length constraints in terms of structure factors
We now show that the intensity sum rule [Eq. (7)] is a direct consequence of the local spin length constraints. Inverting Eq. (A2) we get where the sum over Q is over the first Brillouin zone of the orthorhombic Bravais lattice. The local spin length constraints then take the form which holds for all R, ν and i if we require The intensity sum rule derives from the q = 0 part, namely To see this let us take the general form of Eqs. (A3) and (A4) for any Q, where f = i or 1 [see Eqs. (A3) and (A4)]. Squaring each row and adding them up gives The only Q vectors inside the first Brillouin zone of the orthorhombic lattice that contribute to this sum are the ones corresponding to Q = ± 2 3â and Q = 0, which leads to the intensity sum rule Eq. (7).
Note that the above analysis can be carried out for quantum spins as well, in which case the various spin-spin correlations, such as S i · S j , must be replaced with the corresponding expectation values S i · S j in the quantum-mechanical ground state of the system, and S i · S i becomes S (S + 1).
According to the above, the Bragg peak intensity sum rule is very general and does not depend on the particular values of the microscopic parameters. This generality was missed in [50], because the components M a , M a , M c and M c defined there differ by a relative prefactor of √ 2 from the ones defined here, while this is not the same for the components M b and M b . As a result, the quantity I tot defined in [50] does not correspond to the intensity defined here, which is why that quantity satisfies the sum rule only for sufficiently small J (compare in particular the two panels of Fig. 4 of [50]).
Appendix B: Auxiliary information for the various ansätze 1. Field along the crystallographic b-axis.
a. Low-field phase for H b According to Table I, the low-field ansatz for H b reads where x 1 , x 2 , y 1 , z 1 and z 2 denote Cartesian components of spins. Due to the spin-length constraints x 2 1 +y 2 1 +z 2 1 = 1 and 2x 2 2 + z 2 2 = 1, only three out of these five parameters are independent. The state can also be parametrized in terms of the five symmetry-resolved static structure factor components M a (A), M b C, M c (F), M a (G) and M b (F), which are related to the Cartesian components by Out of the five structure factor components only three are independent, as there are two spin length constraints. One of them is the Bragg peak intensity sum rule, The second constraint reads This illustrates how the local spin length constraints can lead to effective cross-coupling terms between the modulated and uniform components, i.e., the terms M a (A)M a (G) The total energy per site is given by where N is the total number of spin sites. In terms of the structure factor components, E b takes the form where η aA = −Γ+ 2J +K/2, η bC = −K, η cF = −(Γ+2J +K/2), η aG = 1 /2(Γ+ J +K), η bF = 1 /2(3J +K) .
Note that while there are no cross-coupling terms between the modulated and the uniform components, such terms arise from the spin-length constraints, as shown in Eq. (B4).
b. High-field phase for H b For H ≥ H * b the Cartesian components satisfy the relations and we are left with the two-sublattice ansatz [see Table I] with z 2 1 = 1−2x 2 1 , x 1 > 0 and z 1 < 0. In this phase, the modulated components M a (A), M b (C) and M c (F) vanish identically, and we are left with the two uniform components subject to the constraint M b (F) 2 +M a (G) 2 = S 2 . The total energy Eq. (B6) becomes Minimizing gives the following relation between the magnitude of the field H, and the components x 1 and z 1 : In the limit of very large field, H → ∞, Note that the cross-coupling term − √ 2ΓM a (G)M b (F) in Eq. (B11) favors opposite signs of M a (G) and M b (F), given that Γ < 0. And since the magnetic field favors a positive M b (F), it follows that Γ favors a negative M a (G) (the term ∝ g ab in Eq. (B11) also favors a negative M a (G) if g ab > 0; otherwise M a (G) must turn positive at high enough fields). In other words, the sign of the zigzag component along a is fixed by the field.
2. Field along the crystallographic a-axis.
The state can also be described in terms of ten structure factor components. Among these, the first five are the ones we encounter at zero field (and for finite fields along b). The remaining five include three field-induced modulated components M a (C), M b (A) and M c (G), and two field-induced uniform components M a (F) and M b (G). Their dependence on the Cartesian components is The ten structure factor components obey four constraints. One of them is the Bragg peak intensity sum rule given in Eq. (10). The remaining three constraints involve various types of effective cross-coupling terms, similar to the ones we have seen in Eq. (B4). For example, one of these constraints reads: The total energy of the system reads or, in terms of the static structure factor components, where we have introduced b. High-field phase for H a For H ≥ H * a the Cartesian components satisfy the relations and we are left with the two-sublattice ansatz [see Table I] with 2x 2 1 + z 2 1 = 1, x 1 > 0 and z 1 < 0. The only static structure factor components surviving for H ≥ H * a are the uniform components M b (G) and M a (F), subject to the constraint M a (F) 2 + M b (G) 2 = S 2 , and the total energy Eq. (B18) becomes Minimizing the total energy for H ≥ H * a gives the following relation between H, x 1 and z 1 = − 1 − 2x 2 1 : In the limit of H → ∞, we get Note that the cross-coupling term − √ 2ΓM b (G)M a (F) in Eq. (B23) favors opposite signs of M b (G) and M a (F), since Γ < 0. And given that the magnetic field favors a positive M a (F), it follows that M b (G) is negative (consistent with the term ∝ g ab if g ab > 0). Therefore the sign of the zigzag component along b is fixed by the field. Note that the cross-coupling term − √ 2ΓM a (G)M b (F) in Eq. (B35) favors opposite signs for M a (G) and M b (F), since Γ < 0. However, unlike the cases H a and H b, here none of the signs of M a (G) and M b (F) are fixed by the field, meaning that the system can spontaneously choose either M a (G) > 0 and M b (F) < 0 or M a (G) < 0 and M b (F) > 0. The associated broken symmetries are ΘC 2b and C 2c , see Table II.
Appendix C: Proof of Eqs. (17) Here we show a mathematical proof of Eqs. (17). The proof is based on a classical version of the so-called Feynman-Hellmann theorem known in Quantum Mechanics. We begin by writing the total classical energy of the system as a function of the spherical coordinates {θ i , φ i } of the spins (i = 1-N, the total number of spins) and the free parameters of the model, namely J, K, Γ and H: Let us denote the classical ground state configuration for a given set of J, K, Γ and H by {θ * i , φ * i }, where These angles are found by minimizing the total energy Then the minimum of the classical energy, or the classical ground state energy, is given by where the terms in the second line are the individual contributions to the energy from the J, K and Γ interactions, and the Zeeman field, respectively. We can now formulate the classical version of the Feynman-Hellmann theorem by taking the derivative of the ground state energy with respect to the parameter J, as an example. We have and using Eqs. (C3) we get where in the last step we used the fact that f ({θ * i , φ * i }; J, K, Γ, H) depends linearly on J. Similarly, for the other free parameters we get where m is the magnetization per site along the field. To arrive at Eqs. (17) we need to look at the second derivatives of E. For example, The equality ∂ 2 E ∂J∂H = ∂ 2 E ∂H∂J then gives and similarly for the remaining equations of (17). To reduce the autocorrelation time, we have used the standard Metropolis algorithm combined with the over-relaxation algorithm [58,59]. Namely, one Metropolis sweep was performed after completing ten over-relaxation sweeps where each sweep contains N updates. The over-relaxation process with single spin updates is given by [59], where h R,µ is the local effective field at site i = (R, µ). Compared with other MC updates, over-relaxation usually costs less computing time and has less autocorrelations. However, because the over-relaxation update is a micro-canonical process, we adopt the standard Metropolis algorithm to ensure ergodicity of the simulation. In each Metropolis update, one spin S R,µ is randomly chosen and altered to a new direction confined within a cone defined by dθ ∈ [0, π]. We first rotate the coordinate such that the z-axis coincides with S R,µ , i.e. the center of the cone. Similar to the initialization process, here again we generate two random numbers p 1 and p 2 in the interval (0, 1) and take S z = cos(dθ) + [1 − cos(dθ)]p 1 , at which point we generate a random, uniformly distributed unit vector [S x , S y , S z ] within the cone. Afterwards, we rotate the coordinate back to the original coordinate and compute the change in energy ∆E which is related to the probability of acceptance. In the equilibration process, which is the transient time for the system to reach equilibrium, we gradually adjust the magnitude of dθ such that the acceptance ratio keeps staying within [0.4, 0.6]. In the measurement process, we take one measurement of the observables after every ten Metropolis sweeps. Next, in order to obtain the critical temperatures, we have used the Binder cumulants method. The fourth-order Binder cumulant, U 4 = 1 − O 4 3 O 2 2 , where O denote some long-range order parameter, has a scaling dimension of zero; thus the crossing point of the cumulants for different lattice sizes provides a reliable estimate for the value of the critical temperature T c at which the long range order is destroyed. In Fig. 12, we plot the Binder's cumulants for N ∈ {48, 96, 144, 192, 240, 288} for the case when magnetic field with the magnitude H c = 12T is applied along the c crystallographic axis. The corresponding statistical errors are calculated using a Jackknife binning analysis with ten bins [73].