Reentrant incommensurate order and anomalous magnetic torque in the Kitaev magnet $\beta$-$\text{Li}_2\text{IrO}_3$

We present a theoretical study of the response of $\beta$-$\text{Li}_2\text{IrO}_3$ under external magnetic fields in the $ab$, $bc$ and $ac$ crystallographic planes. The results are based on the minimal nearest-neighbor $J$-$K$-$\Gamma$ model and reveal a rich intertwining of field-induced phases and magnetic phase transitions with distinctive signatures that can be probed directly via torque magnetometry. Most saliently, we observe: i) an unusual reentrance of the incommensurate counter-rotating order for fields in the $ab$-plane, and ii) a set of abrupt torque discontinuities which are particularly large for fields rotating in the $bc$ plane, and whose characteristic shape resembles closely the ones observed in the 3D Kitaev magnet $\gamma$-$\text{Li}_2\text{IrO}_3$. An experimental confirmation of these predictions will pave the way for an accurate determination of all relevant microscopic parameters of this 3D Kitaev magnet.

To understand the behavior of the torque one first needs to elucidate the ground state phase diagram for different field orientations and strengths.To this end, we generalize the semianalytical framework developed previously in [30,31,40] to study the field-induced phase diagram for fields in the ab, bc and ac crystallographic planes.This framework builds on the idea that the incommensurate (IC) counter-rotating order observed experimentally at zero field can be treated as a longdistance twisting of a nearby commensurate order with six spin sublattices [40].This approach has already been applied successfully for fields along the orthorhombic directions a, b and c [30,31].It has been shown in particular that the three phase diagrams share a number of qualitative features, including: i) a strong intertwining of the modulated, counterrotating order with a set of uniform orders, some of which give rise to a finite torque signal, ii) the disappearance of the modulated order at a critical field H * , whose value is strongly anisotropic with H * b < H * c H * a , and iii) the presence of a robust zigzag phase above H * .It was further shown that, for H c, the disappearance of the modulated order proceeds via a first-order transition which does not restore all broken symmetries, which in turn gives rise to a second transition at H * * c that survives at finite temperatures [31].The extension of this approach for fields applied in the ab, bc and ac planes shows that the interpolation of the phase diagram from one crystallographic axis to another comes with a number of unexpected results, including: i) a reentrant IC phase that survives in a finite region of field strengths and orientations for fields in the ab plane (and to a lesser extent in the ac plane), ii) a region featuring a two-step disappearance of the IC order, proceeding via two consecutive metamagnetic jumps, and iii) large discontinuous reversals of the magnetic torque, with characteristic sawtooth-like angular profiles that resemble closely the ones reported for γ-Li 2 IrO 3 .These results will be elucidated in Secs.III and IV, but we first recap the main aspects of the microscopic model in Sec.II.

II. Microscopic model
The Hamiltonian of the J-K-Γ model reads [30,31,[38][39][40] Here t ∈ {x, x , y, y , z} labels the different types of NN bonds (see Fig. 1), H t i j includes the couplings between NN sites i and j of a type-t bond, α t is the Cartesian component x, y, or z associated with t, and (β t , γ t ) label the remaining two Cartesian components coupled by the Γ interaction.The prefactor σ t alternates between +1 and −1 for t ∈ {x, y , z} and t ∈ {y, x }, respectively.The term H Z stands for the Zeeman interaction, where the g-tensors carry a staggered, site-dependent off-diagonal element g ab [25] (in the orthorhombic frame), where p i = 1 for spins on the xy chains and −1 for spins on the x y chains (see Fig. 1).More specifically, for fields in the bc, ab and ac planes, the Zeeman terms read (with h = µ B H) where θ bc and θ ab denote the angle between H and b, whereas θ ac is the angle between H and a, see insets of Fig. 2. The coupling parameters are taken as in [31], namely J = 0.4 meV, K = −18 meV, Γ = −10 meV, g aa = g bb = g cc = 2 and g ab = 0.1.
To find the classical ground states we follow the strategy developed in [30,31] and use a general six-sublattice ansatz which is also cross-checked by independent unconstrained minimizations of Eq. (1) on finite-size clusters.As it turns out, the classical ground state for H in the ab plane can be obtained by using the ansatz reported in [30,31] for H a. This is due to the fact that the ansatz for H b is a special case of the more general ansatz for H a. The same situation takes place for H in the bc plane, where the classical ground state is obtained by using more general ansatz for H c. By contrast the classical ground states in the ac plane cannot be obtained from either the ansatz for H a or the ansatz for H c, but by a more general unconstrained six-sublattice ansatz.

III. Phase diagram
Figure 2 summarizes the emerging picture for the phase diagram as we rotate from one crystallographic axis to another.The blue solid line shows the boundary of the IC order, i.e., the critical field H * , as obtained by tracking the associated structure factor I I (see definition in App.A and [31]), whose field evolution is shown in the insets of Fig. 2 for a number of representative field orientations.Altogether, the interpolation of the phase diagram between different crystallographic axes comes with a number of unexpected results.Most notably, we find that while the critical field H * grows very fast as we approach the a axis, reflecting the inequalities H * a found in [31], the way H * grows in the ab-plane gives rise to a reentrant IC order in a finite range of field orientations (θ ab ∈ [62 • , 69 • ] and θ ab ∈ [111 • , 118 • ]), see inset of Fig. 2 (b).In this range, the magnetic field drives the system through four different phases: the low-field six-sublattice phase (including the IC plus the uniform orders), followed by a two-sublattice order (including only the uniform orders), followed by the reentrant six-sublattice order, and then back to the two-sublattice order (whose robust zigzag component diminishes only at H → ∞).
Our classical Monte Carlo (MC) simulations shown in Fig. 3 for θ ab = 69 • demonstrate that this quite unusual reentrant behavior remains robust in a finite temperature range (up to 18 K for the parameters chosen in the present study) and should therefore be observable experimentally, given in addition that the reentrance can occur at a field as low as 22 T (for θ ab = 69 • ).The critical temperatures shown by red circles in Fig. 3 are computed using the crossing of the Binder cumulant [31,41] where O Q=2â/3 is the modulated order (see App.A and [31]).
A similar reentrant behavior is also observed for fields in the ac plane, for θ ac ∈ [34.8 • , 35.1 • ] (see e.g. the dashed vertical line in Fig. 4), although this region is likely too narrow to survive at high enough temperatures to be observed.
Figure 4 shows in addition another interesting feature, which will show up in our discussion of the torque signal below.Namely that, in the region extending for θ ac ∈ [17.0 • , 34.8 • ], the disappearance of the IC order proceeds with two consecutive metamagnetic jumps.Such jumps can be seen in the inset of Fig. 2 (c).The characteristic fields at which these jumps take place are denoted by H * ac,1 and H * ac,2 in Fig. 4. The two-step disappearance of the IC order is accompanied by a similar two-step behaviour of the uniform (zigzag and ferromagnetic) orders, as seen in the evolution of the associated structure factor I V (see Fig. 9 (c)).This arises from the fact that the total intensity from all intertwined orders obeys the intensity sum rule 2I I +I V = S 2 [30,31].
To conclude this section, we comment on the red dashed lines shown on top of the boundaries of the IC phase in Fig. 2.These lines correspond to the following approximations for the critical fields These expressions arise from simple geometrical considerations, i.e., by associating the onset of the transition with the field at which its projection along the axes b or c reaches the critical fields H * b or H * c , respectively.The expressions work surprisingly well in a a wide range of angles, especially for fields in the bc and ab planes.Importantly, the above linear relationship between H * bc and 1/ cos θ bc has been reported by Modic et al for γ-Li 2 IrO 3 [12].

IV. Magnetic torque
Having established the general form of the phase diagram we now turn to the magnetic torque, whose behavior as a function of field strength and orientation is summarized in Fig. 5.The torque τ = (τ a , τ b , τ c ) is given as τ = m×H, where   m = (m a , m b , m c ) is the magnetization per site.The data shown in Fig. 5 reveal a wealth of complex features that are distinctive for each crystallographic plane.Some of these features reflect naturally the phase transitions discussed above, but there are again some additional striking features of experimental relevance.Among these, the most notable is the large discontinuous reversal of the torque at H c as the field rotates in the bc plane, and its characteristic 'sawtooth' shape (see Fig. 5 (d)) that resembles closely the one reported in γ-Li 2 IrO 3 [12,13].This discontinuity was already predicted in [31] and is attributed to the fact that the uniform (zigzag and ferromagnetic) orders that drive the torque are chosen spontaneously inside the region 0 < H ≤ H * * c for H c. A torque reversal of similar origin occurs at H a when H rotates in the ab plane (and for 0 < H ≤ H * a ), but the height of this reversal is likely too small to be observed.
Before we discuss these features in detail let us first briefly recall the main aspects of the torques for H a, H b and H c [31], which are included in the insets of Fig. 5.

A. Magnetic torque for fields along the crystallographic axes
The torque for H b is zero for all the fields.For H a and H c, the torques are, correspondingly, along c and a.The torque for H a remains non-zero up to the critical field H * a , whereas the torque for H c remains non-zero up to the second critical field H * * c , which marks the transition into the fully polarized state [31].In both cases, the sign of the torque is chosen spontaneously and its magnitude is proportional to the m b -component of the magnetization, τ c /H = m b and τ a /H = m b .
Note, however, that m b is significantly different in the two cases, leading to a much larger |τ a | for fields along c compared to |τ c | for fields along a [31].Note also that for H c, the zigzag and ferromagnetic orders are comparable to each other, but the dominant contribution to τ a comes from the ferromagnetic order, since g bb g ab .Moreover, the uniform orders have a non-monotonic behavior as a function of H and exhibit a finite jump at the critical field H * c , leading to a sharp kink in τ a , see dark blue line in Fig. 5 (a) associated with θ bc = 90 • .

B. Magnetic torque for fields away from the crystallographic axes
Fields within the bc plane.This is the setup explored experimentally in γ-Li 2 IrO 3 [12,13].Here the torque is along the a axis, with The sign of τ a for H c (θ bc = 90 • ) is chosen spontaneously below H * * c , is zero for H b (θ bc = 0 • ) [31] and is chosen by the field for any intermediate angles.
The corresponding torque data are shown in Figs. 5 (a, d).At low θ bc and H we find τ a ∝ H 2 sin θ bc , which is the expected behavior in the linear magnetization regime, where m b = χ b H cos θ bc , m c = χ c H sin θ bc and τ a = 1 2 (χ b − χ c )H 2 sin 2θ bc , where χ b and χ c are the magnetic susceptibilities.This behavior ends with a kink at the critical field H * bc associated with the disappearance of the IC order, which happens via a second order phase transition.
With increasing θ bc , the low-field behavior of τ a /H turns from linear to nonlinear (see Fig. 5 (a)), and for θ bc close to 90 • the kink at H * turns to a discontinuity, consistent with previous results for a metamagnetic transition at H c [31].The size of the discontinuity depends non-monotonously on θ bc , with the largest jump taking place at θ bc = 90 • .At this angle, there is a second phase transition at H * * c , above which the classical magnetization saturates and the torque vanishes identically (dark blue line in Fig. 5 (a)).Recall that the presence of the second transition at high field is due to the fact that the zigzag order for H c can take two possible orientations, and the system chooses one of them spontaneously.For all other fields in the bc-plane, the orientation of the zigzag order is determined by the field and the second phase transition is absent.For all θ bc 90 • , the torque is slowly approaching zero with increasing field, and becomes identically zero only at infinite field.
Turning to the angular dependence of the torque (see Fig. 5  (d)), the low-field sinusoidal behavior turns into a characteristic sawtooth shape at high fields, featuring a large torque reversal at θ bc = 90 • .Such a reversal is also present at lower fields, but becomes much more pronounced at higher fields.The reversal is due to the explicit breaking of the discrete symmetries ΘC 2b and C 2c by the component of the field along b [31].Note that the sawtooth torque profile will give rise to a pronounced anomaly in the magnetotropic coefficient k = (1/H)∂τ a /∂θ bc (see Fig. 6).This thermodynamic quantity characterizes the curvature of the free energy with respect to the field orientation, and can be measured directly via torsion magnetometry [21].
Field within the ab plane.Because of a very strong anisotropy of the critical fields, H * a /H * b 50, a rather different picture is observed for H in the ab-plane.In this setup the torque is along the c axis, with The corresponding data are shown in Fig. 5 (b, e).Here the low-field behavior τ c ∝ H 2 sin θ ab persists up to a critical field H * ab H * b .At this field the torque shows a kink which persists for a large range of angles, θ ab 50 • and θ ab 130 • .This   behavior is similar to the case of the field sweeping in the bc plane except that at H > H * ab the torque shows almost fieldindependent behavior up to high fields, H ∼ H * a , and then slowly decreases.This behavior reflects the robustness of the uniform orders contributing to m a and m b .
As the field gets closer to a, the contribution from the acomponent of the field increases.For the ranges θ ab ∈ [50 • , 62 • ] and θ ab ∈ [118 • , 130 • ], τ c /H then exhibits two sharp kinks appearing due to the competition between the two contributions in the Zeeman term in Eq. ( 4).In particular, Fig. 5  Here, the competition between the two contributions in Eq. ( 4) leads to the unusual reappearance of the modulated order at intermediate fields, as discussed above.As the field further approaches the a-axis (70 • θ ab 110 • ), these two kinks disappear but another interesting feature is observed.At the angles and magnitudes of the field at which the first and the second term in Eq. ( 11) become equal, the magnetic torque simply vanishes.In Fig. 5 (e) this situation is shown for θ ab = 70 • (110 • ) and θ ab = 80 • (100 • ).Finally, when θ ab = 90 • (H a), the torque signal is barely visible since the contribution of the first term in Eq. ( 11) becomes vanishingly small for all values of H. Also, for H a the sign of the ferromagnetic uniform order ( b) is spontaneously chosen, which results in a small jump in the torque as shown in the inset of Fig. 5(e).
Field within the ac plane.Unlike the two previous cases, when we sweep the field in the ac plane the torque is no longer perpendicular to the plane of the field.In particular, all three components τ a , τ b and τ c remain finite for all angles θ ac {0 • , 90 • , 180 • }.This aspect leads to a complex non-monotonic evolution of the torque.Here, the torque components are related to the magnetization components via  Two critical fields at some field directions are also clearly seen in τ a /H and τ c /H (see Fig. 7).At angles close to the caxis, both τ a /H and τ c /H show a jump at H * ac H * c and then go to zero above the critical field H * * ac H * * c .When the orientation of the field is closer to the a-axis, both τ a /H and τ c /H remain negligibly small, and the total torque is oriented almost along the b-axis.
Finally, let us discuss what happens with the torque in the region depicted in Fig. 4, which includes a narrow range (around θ ac = 35 • ) featuring a reentrant IC phase and a more extended range (θ ac ∈ [17.0 • , 34.8 • ]) featuring a two-step disappearance of the IC order.Figure 8 shows what happens at two representative angles, θ ac = 35 • and θ ac = 34.5 • , one from each range.As expected, the torque components show three abrupt discontinuities at θ ac = 35 • and two discontinuities at θ ac = 34.5 • , reflecting the respective evolution of the IC order.

V. Discussion
The above numerical results demonstrate a rather complex interplay between the modulated and uniform (zigzag and ferromagnetic) orders of β-Li 2 IrO 3 as we rotate the field in different crystallographic planes.The resulting phase diagram shows a number of unexpected results.Most saliently, we find a rare example of a reentrant IC phase for fields rotating in the ab plane.The stability region of this reentrant phase extends in a finite range of field orientations (θ ab ∈ [62 • , 69 • ] and θ ab ∈ [111 • , 118 • ]), field strengths and temperatures.This unusual behavior is yet another manifestation of the strongly anisotropic character of the underlying microscopic interactions in β-Li 2 IrO 3 , and we have provided numerical estimates for its experimental verification.A similar reentrant behavior is also observed for fields sweeping in the ac-plane but only for a very narrow range of angles.
Our numerical results for the behavior of the magnetic torque offers further insights and reveals a number of distinctive signatures that are unique for each crystallographic plane.Quite generally, the sequence of field-induced phase transitions give rise to torque anomalies (jumps and kinks) which can be directly seen with torque magnetometry.Among the distinctive aspects of the torque is the observation of a large discontinuous reversal as the field rotates in the bc plane and passes through the c axis.The field orientation profile of this sharp discontinuity resembles closely the sawtooth-like behaviour reported by Modic et al in γ-Li 2 IrO 3 [12,13].In that compound the discontinuity extends well beyond H * c 16 T and persists up to 55 T, which likely corresponds to the specific value of H * * c in γ-Li 2 IrO 3 [13].While Modic et al [42] attributed this torque anomaly of γ-Li 2 IrO 3 to the presence of a chiral spinorder, the more conventional symmetry breaking mechanism presented here for β-Li 2 IrO 3 appears more likely, given the close similarity of the two compounds.
The experimental observation of the above distinctive features can be used to map out a quantitative phase diagram and extract accurate values of the microscopic exchange couplings of β-Li 2 IrO 3 .This compound is rather unique in this respect, because its microscopic model has been better understood [30,31,[38][39][40] compared to its 2D and 3D analogues, αand γ-Li 2 IrO 3 .
Given the static structure factor components, we can compute the corresponding contribution to the Bragg peak intensities I tot = 2I I + I V .We first examine the modulated order contribution which is given by The uniform orders can have six different components, which correspond to the zigzag (ZZ) order or ferromagnetic (FM) order along a, b, or c direction, i.e., M a (G), M a (F), M b (G), M b (F), M c (G), and M c (F).As discussed in Ref. [31], when the field is applied along one of the orthorhombic axes, the uniform components are: (A5) Note that for H c, M a (G) and M b (F) only remain finite for H < H * * c .Correspondingly, when the field is rotated on the crystallographic planes, there is always a competition between those components in Eq. (A5). Figure 9 shows the individual non-vanishing contributions I V,i = |M i | 2 of each uniform component to the total Bragg peak intensity I V = 6 i=1 I V,i .This complicated evolution clearly signifies the anisotropic character of the Hamiltonian in Eq. ( 1).

FIG. 1 .
FIG.1.Orthorhombic unit cell of the hyperhoneycomb lattice.The five NN bonds of the J-K-Γ model are marked in red for t ∈ {x, x }, green for t ∈ {y, y }, and blue for t ∈ {z}.

FIG. 2 .
FIG. 2. Global phase diagram as a function of magnetic field strength H (vertical axes) and angular orientation (horizontal axes) as we rotate in the bc plane (a), ab plane (b), and ac plane (c).The color coding tracks the evolution of the Bragg peak intensity I I of the modulated IC order, with the thick blue solid line showing its boundary, i.e., the critical field H * .The red vertical line at H c depicts the fully polarized state at H H * * c .The insets show the field dependence of I I at some representative angles.

FIG. 3 .
FIG. 3. The field-temperature phase diagram for H in the ab plane with θ ab = 69 • obtained from classical MC simulations.The critical temperatures at which the modulated order disappears are shown by red dots and the black line is a guide to the eye.

1 FIG. 4 .
FIG. 4. Evolution of the Bragg peak intensities, I I , when the field direction in the ac-plane is within the range of angles, θ ac ∈ [34 • , 36 • ].

0 /FIG. 5 .
FIG.5.Field (a-c) and angular (d-f) dependence of the ratio of the torque to the magnetic field: (a,d) τ a /H for the field applied in the bc plane; (b,e) τ c /H for the field applied in the ab plane; (c,f) τ b /H for the field applied in the ac plane.The insets in (d, e) zoom-in the jumps of τ a /H at H c and τ c /H at H a, respectively.These jumps originate from the explicit breaking of some of the discrete symmetries by the field component along b[31].

FIG. 6 .
FIG.6.Magnetotropic coefficient k = (1/H)∂τ a /∂θ bc as a function of θ bc for several field strengths.The magnetotropic coefficient at θ bc = 90 • is not well-defined due to the discontinuity of τ a .

FIG. 7 .
FIG.7.Evolution of (a) τ a /H and (b) τ c /H as a function of H for fields rotating in the ac plane, for several field orientations θ ac .To each angle there correspond two curves with opposite torque signs, reflecting the fact that this sign is chosen spontaneously when H ⊥ b.
(b)  shows that τ c /H for θ ab = 60 • (120 • ) has two kinks, one at low field close to H * b and the other at intermediate field close to H * a .Especially interesting is the behavior of the torque for θ ab ∈ [62 • , 69 • ] and θ ab ∈ [111 • , 118 • ].

FIG. 8 .
FIG. 8. Non-monotonous evolution of |τ a |/H (a), |τ b |/H (b) and |τ c |/H (c) as a function of H at θ ac = 34.5 • and θ ac = 35 • , where the destruction of the IC order happens via two-step process.At θ ac = 35 • the reentrant behavior is observed in addition to the two-step process.

FIG. 9 .
FIG. 9. Evolution of the individual Bragg peak intensities from the uniform orders with field in (a) bc-plane, (b) ab-plane, and (c) ac-plane.
H a : M a (G), M a (F), M b (G), M b (F), H b : M a (G), M b (F), H c : M a (G), M b (F), M c (F).