Anomalous and Anisotropic Nonlinear Susceptibility in the Proximate Kitaev Magnet $\alpha$-RuCl$_3$

The leading order nonlinear (NL) susceptibility, $\chi_3$, in a paramagnet is negative and diverges as $T \rightarrow 0$. This divergence is destroyed when spins correlate and the NL response provides unique insights into magnetic order. Dimensionality, exchange interaction, and preponderance of quantum effects all imprint their signatures in the NL magnetic response. Here, we study the NL susceptibilities in the proximate Kitaev magnet $\alpha$-RuCl$_3$ which differs from the expected antiferromagnetic behavior. For $T<T_c$ = 7.5 K and field $B$ in the ab-plane, we obtain contrasting NL responses in low ($<$ 2 ${T}$) and high field regions. For low fields the NL behavior is dominated by a quadratic response (positive $\chi_2$), which shows a rapid rise below $T_c$. This large $\chi_2>0$ implies a broken sublattice symmetry of magnetic order at low temperatures. Classical Monte Carlo (CMC) simulations in the standard ${K-H-\Gamma}$ model secure such a quadratic ${B}$ dependence of ${M}$, only for ${T}$ $\approx$ ${T}_c$ with ${\chi}_2$ being zero as ${T}$ $\rightarrow$ 0. It is also zero for all temperatures in exact diagonalization calculations. On the other hand, we find an exclusive cubic term (${\chi}_3$) describes the high field NL behavior well. ${\chi}_3$ is large and positive both below and above ${T}_c$ crossing zero only for ${T}$ $>$ 50 K. In contrast, for $B$~$\parallel$~c-axis, no separate low/high field behaviors is measured and only a much smaller ${\chi}_3$ is apparent.


I. INTRODUCTION
Since the demonstration by Kitaev [1] of the existence of a quantum spin liquid state in two dimensions through exact calculations on a honeycomb lattice, there has been an intense experimental search for its realization in real material systems. Aided through insights provided by Jackeli and Khaliullin [2] to engineer (bond dependent) exchange interactions, the materials search has identified several promising systems with prominent attention given thus far to two candidates, the iridium oxides and ruthenium chloride [3,4]. Despite the presence of a very large exchange energy (of the order of 100 K) these systems do not order down to comparatively low temperatures. When they do order magnetically, the features observed in both microscopic probes such as neutron scattering and macroscopic measurements such as magnetometry and thermal response are highly unusual and contrast dramatically with what is expected from conventional magnets. In their thermal transport [5], thermodynamic [6], and microwave response [7], they provide tell-tale signs as to the presence of fractionalized excitations sought after in Kitaev magnets.
In the proximate Kitaev spin liquid candidate, α-RuCl 3 neutron scattering experiments show a low temperature magnetic excitation spectrum consisting of sharp spin wave peaks and a continuum associated with fractional excitations [6,8]. A magnetic transition that sets in at 7.5 K where the spins assume a zig-zag chain pattern located in the ab-plane with two neighboring * https://orcid.org/0000-0002-6792-7463 chains being antiferromagnetically aligned [9] is also established. Raman scattering studies also reveal unconventional magnetic excitations with a broad continuum whose temperature dependence is apparent over a large scale compared to the magnetic ordering temperature [10,11]. The linear susceptibility shows a discontinuity in slope [12] at T N = 7.5 K and exhibits substantial in-plane anisotropy that persists in the normal state [13]. The high temperature (T > 150 K) susceptibility is convincingly Curie like, however, there is an extended intermediate "Kitaev paramagnetic" region [6] beyond T c . The out of plane anisotropies are also significant: the susceptibility parallel to the c-axis is nearly an order of magnitude smaller with only a minor signature at the 7.5 K transition. These magnetic signatures are in stark contrast to what is known in conventional 2D (insulating) antiferromagnets [14].
In this communication we report measurements of the nonlinear DC susceptibilities, χ 2 and χ 3 , in α-RuCl 3 and illustrate that they probe many key aspects of the proximate Kitaev spin liquid state as well as the zigzag antiferromagnetic phase. The equilibrium magnetization in any system can be written in the general form: where the coefficients represent the various order susceptibilities. In particular, the coefficient χ 2 is non zero when time reversable symmetry is broken, while preserving lattice symmetry as for instance [15] in ferromagnets. The nonlinear susceptibility χ 3 in a classical paramagnet is negative at all temperatures and diverges [16] as T → 0 while χ 2 is non existent due to time reversal symmetry.
The negative divergence in χ 3 can be interrupted however if the system develops long range magnetic order. We illustrate in fig. 1(a-c) the three known types of characteristic behavior of both the linear (χ 1 ) and the nonlinear (χ 2 and χ 3 ) susceptibilities for typical ferromagnets, bipartite antiferromagnets and spin glasses, respectively [15,[17][18][19]. In all three cases, significant non zero χ 3 is found only in the vicinity of the critical temperature. It is worth noting that χ 3 assumes only negative values at T > T c for all three types of magnets and tends to zero rather quickly as T increases. The second order susceptibility, χ 2 , is less studied and a non zero χ 2 is possible only when time reversal symmetry is either explicitly or spontaneously broken, as in ferromagnets. Moreover, even for bipartite anti ferromagnets, an effective time reversal symmetry due to the sublattice symmetry results in a vanishing χ 2 even in the ordered state below T c ; see Fig.1

(B).
In α-RuCl 3 we find in contrast that χ 3 assumes significantly positive values over an extended temperature range above the ordering temperature. While there is a divergence of χ 3 at T c as might be expected, its value remains significantly positive down to the lowest temperatures. We also find a quadratic field dependence of M , giving rise to a significantly positive χ 2 which exists only in the ordered state and has a large non-zero value even as T → 0. To our knowledge such behavior of χ 3 and χ 2 have not been observed before. Further, this unexpected quadratic contribution is highly anisotropic and found only when the magnetic field is perpendicular to the high symmetry axis (c-axis).

A. Experimental results
In figure 2, we show the measured magnetization isotherms for B a-axis (i.e. φ = 0 0 as per the definition adopted in [13]) plotted in a manner that facilitates the extraction of the nonlinear susceptibilities. Equation (1) above defines the susceptibility parameters and motivates plots such as Fig. 2. It is apparent from the top nine panels in Fig. 2 that the slope of the lines for B a-axis which are well defined and close to zero at high T , turn significantly positive when lowering T , and increase in magnitude as the ordering temperature T c = 7.5 K is approached. Below T c , the response breaks into two distinct regions, with a crossover threshold field of B* ≈ 2 T . The nonlinear response below B* is considerably larger than at the high field end, 3 T < B < 5 T. Confining ourselves to the quadratic term only, as per Eq. 1 in this low field region, we obtain values of χ 2 as shown in Fig. 3b. The values of χ 2 are positive and large at the lowest temperatures and decrease monotonically towards T c , where it rapidly drops to zero. It is important to note that a non-zero χ 2 is only possible in systems in which time-reversal symmetry is broken [20]. In a strict bi-partite antiferromagnet, however, time reversal symmetry is not broken. The fits described above also yield χ 1 via the intercept which can be used as a consistency check on the linear susceptibility values obtained at constant low field through temperature sweeps. The absence of any feature in χ 1 in the vicinity of 14 K alludes to the high quality of the sample measured.
Another significant feature notable in the nine panels on the top in Fig. 2 is the presence of a clear upward curvature, in the 3 T -5 T region, particularly in the 6 K and 10 K isotherms. This curvature implies the presence of a cubic term which can be extracted by plotting M/B against B 2 on the x-axis as shown in the bottom nine panels of Fig. 2. The values of χ 3 obtained by fitting the linear region in such plots with B between 3 T and 5 T below T c and over the entire field range (0T -5T) above T c , are shown in Fig. 3, panel c. The high linearity of these fits illustrates the absence of a χ 2 term in these regions of the phase diagram. Performing such separate fits in the two regions is the most natural way to analyze our results. Our approach of seperate fits is further motivated by signs of a crossover transition at ≈ 2 T in neutron [12,21], microwave response [7] and differential susceptibility [22] experiments.
As can be seen in Fig. 3c χ 3 is significantly positive near T c and decreases in magnitude as the temperature is increased. A power law fit as shown in the figure provides a reasonable empirical description of the behavior of χ 3 in the (Kitaev) paramagnetic region. Experimentally however, a crossover to negative χ 3 occurs above 50 K (a temperature of the order of the Kitaev-exchange [6]). Below T c the values of the nonlinear susceptibility are obtained through linear fits in two separate regions as explained above. The high field region yields values for χ 3 that exhibit a sharp peak at T = 6.5 K just below the ordering temperature. They remain positive down to the lowest temperature of 2 K with no evidence of a downturn towards zero. While a negative cubic power law dependence is expected from the Curie-law for the third order susceptibility in paramagnets [16], such a positive behavior has not been reported to our knowledge.
We have also measured the nonlinear susceptibility for magnetic fields perpendicular to the ab-plane. Representative magnetization isotherms plotted similar to those in Fig. 2 are shown in supplementary information, Fig. S2. The extracted values of χ 3 are small and positive and increase monotonically as the temperature is lowered as seen in Fig. S3. More importantly, as evident in Fig. S2, the nonlinear response is well described by a single cubic term and there is no evidence for a large quadratic response as found for B a-axis. Given that for B a-axis the deviation from a linear behavior of the magnetization due to a quadratic contribution is FIG. 1. Schematic of the known qualitative behaviors of the linear (χ1) and nonlinear (χ2 and χ3) magnetic susceptibilities in different flavors of magnetic systems: In types A and C, the nonlinear susceptibility, χ3, is negative above the characteristic temperature and has a sharp discontinuity at a characteristic temperature (Tc or Tg). In type B, the sign of χ3 in the paramagnetic regime depends on the coordination number [25]. Below this temperature it is observed to rapidly approach zero as T → 0 for all three phases. χ2 on the other hand exhibits such a discontinuity only for type A (e.g. ferromagnet) where time reversal symmetry is broken [15]. In comparison, in the other two cases it is zero. roughly ten times larger than that arising from the cubic term, such behavior if present for B c-axis -would be easily seen within the resolution of the MPMS SQUID measurements. Considering the anisotropy factor of approximately eight as observed experimentally in linear susceptibility between the ab-plane and the c-axis and assuming that this is entirely due to the g-factor, one can predict an anisotropy in χ 3 of ≈ 8 3 = 512 in reasonable agreement with our experimental results. We also present in the supplementary section, Fig. S4, preliminary data for the in plane nonlinear response when B ⊥ a-axis. While the behavior observed is qualitatively the same as in Fig. 3, there are quantitative differences. Further characterization of this in-plane anisotropy will form part of a separate comprehensive study.

B. Comparison to Traditional Nonlinear Susceptibility Responses
Measurements of the equilibrium third order susceptibility, although rare, have been performed in 2D magnets, frustrated systems, spin glasses and strongly correlated itinerant metamagnets, that is materials which show a rapid rise of magnetization at a critical field B c [23]. In a bipartite antiferromagnet such as FeCl 2 [24], or even in the classic 2D magnets, (C 2 H 5 NH 3 ) 2 CuCl 4 [14], the DC nonlinear susceptibility, χ 3 exhibits a "lambda" like anomaly, just below T c as expected by theory. From below T c the standard response of χ 3 in these systems is a rapid rise to a large positive value at T c , above which it drops with the sign being set by the coordination number [25]. For instance, in the paramagnetic phase of the Kagome system SCGO [26] a large positive value of χ 3 is seen above T c , but it reaches close to zero within 2T c . In many strongly correlated itinerant metamagnets [27,28] in which an order-parameter develops, χ 3 peaks at the ordering temperature, and decreases rapidly as T → 0. Nonlinear susceptibility measurements in such systems typically probe higher order correlations and place strong constraints on the ground state [29,30]. These latter systems are three dimensional electronically but can exhibit a strong magnetic anisotropy due to the g-factor. Nevertheless, in these system, for all directions it is sufficient to include only the cubic term and the possibility of a dual response with a quadratic term has almost never been discussed [31]. Thus, the features reported above in α-RuCl 3 are in contrast to much of what is known about nonlinear susceptibilities in magnets. The unmistakable presence of a large T →0 quadratic term makes α-RuCl 3 a unique 2D quantum antiferromagnet.

C. Comparison to Simulations of the J1-J3-K-Γ Model
In order to understand the nature of the nonlinear susceptibility in α-RuCl 3 it is possible to consider several approaches based on different model Hamiltonians employed thus far [13,32]. Most of these approaches start with the Kitaev-Heisenberg model appended with various choices of off-diagonal terms [33,34]. The correct choice of the Hamiltonian is still very much a matter of debate with the sign of the Kitaev term or even the ne-cessity of the Kitaev term being in question [35][36][37]. We focus on this model and apply two separate calculational tools to study the nonlinear susceptibilities for α-RuCl 3 : (a) we use classical Monte Carlo (CMC) simulations with the model Hamiltonian This so-called J 1 − J 3 − K − Γ model is considered one of the the most successful efforts in modeling α-RuCl 3 [32]. (b) we employ exact diagonalization methods with slight adjustment on parameters as in [38]. For (a) although quantum fluctuations are expected to be strong for such spin-1/2 systems, the fact that the magnet develops a long range antiferromagnetic order justifies the classical spin approximation, as a first step to understanding the magnetic properties of the ordered phase. For (b) the choice is based on the recognition that it reproduces the magnetization isotherms to high fields very well.
With the CMC simulations we reproduce the phase transition to the so called zig zag order at T c ≈ 0.11 K where K is the dominant Kitaev exchange constant, with the other parameter values normalized to K taken as, Γ = 0.5, J 1 = 0.036, J 3 = 0.035, g a = g b = 2.3 and g c = 1.3 [39,40]. Our calculations of the linear susceptibility χ 1 based on the Monte Carlo data, as shown in Fig. 4(a), shows a broad bump near the critical temperature, a feature that is qualitatively consistent with the experiment, Fig. 3(a). However, in stark contrast to the experimental data, our simulations find a positive and significant χ 2 only in the vicinity of T c , Fig. 4(b). Our finding of a vanishing χ 2 at very low temperatures in the zig-zag phase is in fact consistent with an emergent time reversal symmetry of the low T ordered state. To see this emergent symmetry we first note that although there are three equivalent propagation directions for the zig zag order, it has been shown experimentally and theoretically that the ground state of α-RuCl 3 is a single Q zig zag. Spins in such single Q zig zag are collinear, with opposite spins sitting on alternating zig zag chains of the honey comb lattice. The total magnetization of such collinear bipartite antiferromagnetic order is the sum of the two sublattices: M = (M A + M B )n, wheren is the collinear spin direction, and M B = − M A . Their contributions to the second order magnetic susceptibility d (3) cancel each other due to the sublattice symmetry M 3 The persistence of the quadratic χ 2 down to the lowest temperatures in our experiments thus highlights the unusual nature of the zig zag order. This intriguing discrepancy could be due to non trivial quantum fluctuations in the ordered state. Another possible scenario is that the non zero χ 2 results from a highly inhomogenous multidomain zig zag order at low temperatures. Indeed, at the interfaces between different zig zag domains, the two sub lattice collinear spins argument given above no longer applies and a non zero χ 2 could result from the complex non-collinear structure at the domain boundaries. Our experiments also imply that the low temperature χ 2 vanishes at the critical field B * ≈ 2 T. Interestingly, this critical field B * has recently been attributed to the so called Q flop transition identified both in neutron diffraction [12,21] and terahertz spectroscopy [41,42]. Since the presence of the magnetic field in the ab-plane breaks the equivalence of the three zig zag orientations, the Q-flop transition describes a repopulation of zig zag domains in which two energetically unfavorable zig zag domains are replaced by the third one. Such a realignment of the zig zag domains also significantly reduces the non collinear spins residing on the interfaces of different zig zag order, thus giving rise to a reduced χ 2 . Detailed numerical simulations of the Q-flop transition and its effects on χ 2 will be left for a further computational study. The breaking of the sublattice-symmetry, which leads to a nonzero χ 2 in the ground-state, could also come at the dynamical level, for example, due to non trivial quantum fluctuations with non-collinear high-order spin-correlations in the ordered state [27]. Other possibilities such as a stacking of the single-Q collinear zigzag order along the c direction that breaks the sublattice-symmetry cannot be ruled out either. Compared with CMC simulations another intriguing result from our experiments is the persistence of large χ 3 in both high and low temperature regimes. The values of χ 3 from CMC simulations approach zero very rapidly above T c in contrast to the experimental large and positive values that persist for temperatures significantly greater than T c . The low temperature discrepancy of χ 3 can be attributed to the absence of quantum fluctuations in CMC.
However, the importance of such fluctuations is borne out in similar calculations utilizing quantum chemistry methods discussed below. They are also seen in quantum Monte Carlo simulations in the pure Kitaev limit by Kamiya [43] which show a persistent positive χ 3 down to the lowest temperatures. In quantum chemistry methods we use exact diagonalization (ED) of a closely related Hamiltonian employed in [38] which yields accurate results for the magnetization isotherms in good agreement with experimental results at the high field end. The calculated magnetization isotherms in this approach plotted in a manner similar to experiments, are well fit with a single straight line (see Fig. S9 Exact Diagonlization Calculations: Shows the nonlinear susceptibility χ3 in calculations employing exact diagonalization, for B in the ab plane. The results obtained are identical whether B a-axis or B ⊥ a-axis. χ3 rises rapidly to positive values at TN reaches a maximum and saturates to a large positive value at T = 0. However, the zero crossing occurs within 2TN whereas the experimental zero crossing occurs at a much higher temperature. The bottom panel shows the second derivative, d 2 M/dB 2 , at various temperatures. This derivative passes zero at all temperatures -thus implying the absence of a B 2 -term in the magnetization contrary to experimental findings.
zero intercept for all temperatures. This implies that the dual slope response, which we attribute to the presence of the complex multi-domain zigzag order or other sublattice-symmetry breaking mechanisms, is not found in the ED calculations. It is worth noting that the non-zero value of χ 3 for T → 0 is also found in the pure Kitaev model with Quantum Monte Carlo calculations [43]. In such calculations it is possible to secure a crossover of χ 3 at a higher temperature, however only through the antiferromagnetic Kitaev interaction. Many calculations [32] rule out an antiferromagnetic scenario, but our experimental results suggest not to exclude this possibility.
Moving forward, any viable model for α-RuCl 3 must explain the quadratic contribution to the magnetization in the antiferromagnetic zigzag phase evident in our experiments. Further, it also has to account for the persistence of the large positive values of the third order susceptibility for temperatures well above T c . It might be that the details of the material parameters will decide the magnitude and the temperature range over which a positive nonlinearity is stretched.

D. Summary and Outlook
In summary, we have presented measurements of the nonlinear susceptibility in the Kitaev magnet α-RuCl 3 . Most significantly, our work has uncovered an anomalous quadratic response of the magnetization to field which yields a large positive χ 2 in the ordered state as T → 0. This behavior is absent when B c-axis suggesting a strong 2D nature of the order parameter. This anisotropy as well as the measured anisotropy of χ 3 both above and below T c can serve as future characterization tools for pinning down specific models for proximate Kitaev materials. In addition, our observation of extended positive behavior of χ 3 up to 50 K is consistent with previous circumstantial evidence that Kitaev type behavior with associated excitations persist up to fairly large temperatures T ≈ 60 K >> T c [6,10]. The low field crossover i.e. the anomalous response with a quadratic term at low fields does not have to be confined to α-RuCl 3 and the generality of its presence in Kitaev or other spin liquid compounds should be established. Furthermore, our results could place constraints on models which attempt to explain experimental observations in α-RuCl 3 and similar compounds outside the realm of Kitaev physics [35,44]. Is the anomalous χ 2 and the extended positive χ 3 a natural consequence of Kitaev type models when quantum fluctuations are correctly accounted for or is it a peculiarity of α-RuCl 3 [45]? We note that the conditions to approach pure Kitaev without the nuisance of magnetic order experimentally are fairly easy to reachthe zig zag order is destroyed in α-RuCl 3 at relatively low pressures of ≈ 1 GPa [46,47]. Equally important to carry this work forward will be attempts to predict an in plane anisotropy of the nonlinear response in a quantitative manner followed by further detailed experimental work in this regard.

A. Materials and Experiments
Our measurements were performed on high quality single crystals similar to those used for recent linear magnetometry measurements [13]. A Quantum Design Magnetic Property Measurement System SQUID magnetometer capable of reaching 5 T was employed to obtain magnetization isotherms in the temperature range 2 -300 K.

B. Simulations
Classical Monte Carlo as well as quantum chemistry based exact diagonalization (ED) calculations of the lin-ear and nonlinear magnetic response in the K − H − Γ model [32], were performed using available packages on a supercomputer.