Effects of Dzyaloshinskii-Moriya interactions in volborthite: Magnetic orders and thermal Hall effect

Volborthite offers an interesting example of a highly frustrated quantum magnet in which ferromagnetic and antiferromagnetic interactions compete on anisotropic kagome lattices. A recent density functional theory calculation has provided a magnetic model based on coupled trimers, which is consistent with a broad 1 3 -magnetization plateau observed experimentally. Here we study the effects of Dzyaloshinskii-Moriya (DM) interactions in volborthite. We derive an effective model in which pseudospin1 2 moments emerging on trimers form a network of an anisotropic triangular lattice. Using the effective model, we show that for a magnetic field perpendicular to the kagome layer, magnon excitations from the 1 3 -plateau feel a Berry curvature due to the DM interactions, giving rise to a thermal Hall effect. Our magnon Bose gas theory can explain qualitative features of the magnetization and the thermal Hall conductivity measured experimentally. A further quantitative comparison with experiment poses constraints on the coupling constants in the effective model, promoting a quasi-one-dimensional picture. Based on this picture, we analyze low-temperature magnetic phase diagrams using effective field theory, and point out their crucial dependence on the field direction.


I. INTRODUCTION
The last two decades have witnessed increasing interest in highly frustrated quantum magnetism, boosted by the development in theoretical concepts and simulation methods as well as a huge variety of material realizations [1][2][3]. On one hand, such magnets do not easily find energetically stable states, which results in a wide range of behavior under perturbations. On the other hand, they offer an attractive possibility of a quantum spin liquid (QSL), which evades any ordering in terms of a conventional order parameter down to zero temperature [4][5][6]. Antiferromagnets on triangular and kagome lattices are typical examples of geometrical frustration in two dimensions. Numerical studies on the spin- 1 2 kagome antiferromagnetic Heisenberg model have provided indications of a QSL ground state with either gapped [7][8][9][10] or gapless [11][12][13][14] low-energy excitations although its precise nature is still under active debate. The enigmatic nature of the kagome antiferromagnet has stimulated experimental studies on a number of copper minerals [15] such as herbertsmithite [16,17], volborthite [18,19], and vesignieite [20], which host layers of spin- 1 2 moments arranged in a kagome pattern. Thermodynamic and neutron scattering measurements for herbertsmithite have identified the gapless QSL behavior with fractionalized excitations [16,21] while the NMR experiment has detected a small intrinsic excitation gap [22]. As for the scenario of rich behavior due to frustration, the spatially anisotropic triangular antiferromagnet Cs 2 CuCl 4 offers a particularly interesting example. Experiments have revealed continuum of excitations indicative of spin fractionalization [23,24], and rich magnetic phase diagrams that crucially depend on the field direction [23,25]. A theoretical understanding of these results has been made by view-ing the system as weakly coupled Heisenberg chains [26][27][28][29][30], and thereby an extreme sensitivity to weak magnetic anisotropy and inter-layer couplings has been pointed out [27].
In this context, volborthite Cu 3 V 2 O 7 (OH) 2 ·2H 2 O is a fascinating material for which a wealth of field-induced phenomena have been observed in powder [31][32][33][34] and single-crystal [35][36][37][38][39][40] samples. The material was originally considered as a candidate for a spin- 1 2 kagome antiferromagnet [18], and a sign of strong frustration has been found in the low magnetic transition temperature around 1 K in comparison with the Curie-Weiss temperature −155 K [32,[41][42][43]. X-ray diffraction measurements for single crystals have, however, suggested highly anisotropic arrangements of magnetically active orbitals at the crystallographically distinct Cu sites, indicating a strong spatial anisotropy in magnetic interactions [19,35,43]. Furthermore, magnetization measurements for single crystals have revealed a wide 1 3magnetization plateau which starts at H = 27.5 T [35]; according to recent Faraday rotation measurements, this plateau continues as high as 100 T or even above 160 T, depending on the sample setting condition [39]. Such an extremely wide 1 3 -plateau is in sharp contrast with the relatively short plateau in the kagome antiferromagnetic model [44][45][46]. Below this plateau, NMR measurements have identified three distinct phases [35,36], as schematically shown in Fig. 1(a). Phase II shows doublehorn NMR spectra indicative of an incommensurate spindensity-wave (SDW) order while indications of bimagnon condensation upon entering Phase N from the plateau have been found in Refs. [36,37]. Thermal measurements [37] have indicated that Phase N is divided into two phases, N 1 and N 2 .
In Ref. [47], density functional theory (DFT) calcu-lations have been performed on the basis of the singlecrystal structural data [35] to determine the microscopic spin model of volborthite. Using DFT+U , four leading exchanges have been identified, as schematically shown in Fig. 2(a): antiferromagnetic J and J 2 as well as ferromagnetic J and J 1 with a distinctive hierarchy J > |J 1 | > J 2 , |J |. The dominance of the J coupling naturally leads to a coupled-trimer picture: on each trimer formed by the J coupling, the spin states are restricted to the lowest-energy doublet, which can be viewed as pseudospin-1 2 states, at zero field and such pseudospins interact with each other through the inter-trimer couplings. This sharply contrasts with the coupled frustrated chain model with J = J [48], which was obtained previously based on powder structural data. In the coupled-trimer model, the 1 3 -plateau state can be interpreted as a product of polarized trimers, and a wide plateau extending to H = 225 T has been predicted. By means of a strong-coupling expansion, an effective model has been derived for the pseudospin- 1 2 degrees of freedom living on an anisotropic triangular lattice, as shown in Fig. 2(b,c). This model shows a tendency towards condensation of magnon bound states preceding the plateau, indicating the emergence of a bond nematic order  Schematic magnetic phase diagrams at low temperatures based on experiments (single crystals) and theory. See Fig. 2 for the definitions of directions. In (a) the experimental results for H z, Phase II shows double-horn NMR spectra indicative of an incommensurate SDW order [35,36]. Phase N can be divided into two phases, N1 and N2 [37]. The theoretical results in (b) and (c) based on effective field theory for a quasi-one-dimensional regime exhibit a marked sensitivity to the field direction due to the DM interactions. In (b), we assume that the magnitude of the strongly relevant interaction γ between the second-neighbor chains [defined in Eq. (58) later] is suppressed (|γ | 0.1 K). The nature of the obtained phases are described in detail in Sec. IV. Shaded areas indicate the regimes where the field-theoretical approach is less effective (owing to spatially oscillating interactions or a vanishing velocity). However, the analysis of the multimagnon spectra indicates that the bond nematic state due to condensation of bimagnons appears just below the plateau for a certain parameter range [47]. 54]. This can provide a scenario for Phase N (or N 2 ) observed experimentally [55]. The coupled-trimer model has also stimulated theoretical studies of possible QSLs above magnetic ordering temperatures [56,57].
The theoretical analyses of Ref. [47] have mostly been based on isotropic Heisenberg interactions. Although Ref. [47] has also given estimates of the Dzyaloshinskii-Moriya (DM) interactions for the leading couplings J and J 1 , their effects on the magnetic properties have not been analyzed in detail. Given a high sensitivity of frustrated magnets, the DM interactions can significantly influence the low-temperature magnetic orderings. Furthermore, these interactions can give rise to nontrivial transport properties such as a magnon thermal Hall effect [58][59][60][61]. In fact, a thermal Hall effect has been observed in volborthite in magnetic fields up to 15 T perpendicular to the magnetic layer [38]. The dependence of the observed transverse thermal response on the magnetic field indicates that the effect is due to spin excitations. Interestingly, the effect has been observed even above the magnetic transition temperatures, where the system may be described as a cooperative paramagnet.
In this paper, we theoretically study the effects of DM interactions in volborthite on the basis of the coupledtrimer model of Ref. [47]. By incorporating the effects of the DM interactions, we derive an effective pseudospin-1 2 model on an anisotropic triangular lattice. In the resulting model, the magnetic anisotropy is characterized by a single effective DM vector (see D z 1 in Fig. 2), which leads to a significant simplification of the analysis. We show that for a magnetic field perpendicular to the kagome layer, magnon excitations from the 1 3 -plateau (fully polarized pseudosopins) feel a nonzero Berry curvature due to the effective DM interaction, leading to a thermal Hall effect. This effect disappears when the field is changed to the direction of the screw axis (the x direction of Fig.  2), reflecting a symmetry. Although our analysis of the thermal Hall effect is based on the spin wave picture valid just below the plateau, we expect that the magnitude and qualitative features of the thermal Hall conductivity do not change abruptly as we lower the magnetic field to the regime |H| 15 T investigated experimentally [38]. Comparison of the theory with experimental data of the magnetization process and the thermal Hall conductivity poses constraints on the coupling constants in the effective model, promoting a quasi-one-dimensional picture. Based on this picture, we further analyze magnetic orders at low temperatures using effective field theory (in close parallel with the theory for Cs 2 CuCl 4 [26,27]), and point out that the magnetic phase diagram can sensitively depend on the field direction owing to the DM interaction, as shown in Fig. 1(b,c). We thus find a striking similarity to the physics of Cs 2 CuCl 4 .
The rest of this paper is organized as follows. In Sec. II, we describe the microscopic spin-1 2 model obtained by DFT+U for volborthite in Ref. [47]. We then perform a strong-coupling expansion to derive an effective pseudospin-1 2 model on an anisotropic triangular lat-  (3) and (b,c) the effective pseudospin-1 2 model (17) for a single magnetic layer of volborthite. Here, (b) shows only the nearest-neighbor interactions which appear in the first-order strong-coupling expansion, and (c) shows further-neighbor interactions which appear in the second-order expansion. Some of the interactions are shown only for representative bonds for brevity. See Eqs. (4), (5), and (16) and Table I for the estimates of the (effective) coupling constants based on DFT+U [47]. For simplicity, we work with the approximate lattice in which all the spins on the layer reside in the xy plane, and the central sites of neighboring trimers are connected by the vectors u = (b, −a, 0)/2, v = (b, a, 0)/2, and b = (b, 0, 0). In (a), three spins on each trimer are labeled by j = 1, 2, 3 in ascending order from left to right. In defining DM vectors, every bond is oriented from left to right, and an arrow (a dot or a cross) indicates the in-plane (out-of-plane) component. The (effective) DM vectors are subject to constraints due to crystal symmetries such as the inversion I about the center of each trimer and the two-fold screw (21) axis along each J1-J2 chain (see Appendix A). A bar on a DM vector indicates π rotation around the x axis.
tice. In Sec. III, we perform a spin wave analysis starting from the 1 3 -plateau state, and analyze the field-and temperature-dependences of the magnetization and the thermal Hall conductivity. In Sec. IV, we analyze lowtemperature magnetic orders in light of an effective field theory for a quasi-one-dimensional regime. In Sec. V, we present a summary and an outlook for future studies.

II. MODEL
In this section, we describe the microscopic spin-1 2 model obtained by DFT+U [47], and then perform a strong-coupling expansion to derive an effective pseudospin-1 2 model on an anisotropic triangular lattice. Before presenting the theoretical models, let us briefly review the low-temperature crystal structure of volborthite. Below 155 K, volborthite shows the P2 1 /a structure (space group No. 14) with the lattice constants [19,35] We introduce the Cartesian coordinate (x, y, z) in such a way that the x and y axes are along the crystallographic b and a directions, respectively. In this coordinate, the primitive vectors of the lattice are represented as b = (b, 0, 0), a = (0, a, 0), c = c(0, cos β, − sin β). (2) A single kagome layer extends in the xy or ab plane as shown in Fig. 2, and two kinds of kagome layers are stacked alternately in the c direction. The layer dependence of exchange couplings is estimated to be small [47], and neglected in the subsequent analysis.
A. Microscopic spin-1 2 model The microscopic spin-1 2 model determined by DFT+U [47] is shown in Fig. 2(a). We label each trimer by the location r of its central site. Let S r,j (j = 1, 2, 3) be the three spin-1 2 operators on the trimer. There are two types of trimers, A and B; the sets of the central sites of those trimers are denoted by A and B as well.
The Hamiltonian is given by where the vectors u, v, and b connect between neighboring trimers as shown in Fig. 2. For a trimer at r ∈ X = A, B, the intra-trimer interactions are given by Here, a Zeeman energy in an external magnetic field H has been included. Throughout this paper, a bar on a vector indicates π rotation around the x axis. The inter-trimer interactions along u and v are given for r ∈ A by D = (D x , D y , D z ) = (12, 7, −17) K, The DM vectors D and D 2 on the J and J 2 bonds have not been estimated. We note that the presence of the ferromagnetic J 1 < 0 coupling is consistent with negative magnetostriction along the x axis observed recently [40].
For the direction of the external field H, we mainly consider the following two cases: Here, H and h point in mutually opposite directions because of the negative g-factor −g < 0 of an electron. Experiments have mainly been performed for H ±ẑ.

B. Isolated trimer
We first consider the Hamiltonian for an isolated trimer of type A, which is given by We note that the Hamiltonian for a trimer of type B can be treated by replacing D byD in the above. When D = 0 and h = 0, because of the SU(2) symmetry and the inversion symmetry around the site 2, the eigenstates of Eq. (7) are classified into the quadruplet {|q µ }, the even-parity doublet {|d µ }, and the odd-parity doublet {|d µ }, where µ is the eigenvalue of 3 j=1 S z j . Their wave functions are given by . The Zeeman term in Eq. (7) commutes with the Heisenberg terms, and only shifts the eigenenergies by −hµ when h = (0, 0, h). The eigenenergies in this case are calculated as At h = 0, the ground states are the even-parity doublet |d ± 1 2 ; they are split for h = 0. For low temperatures and low fields (T, |h| J), it is therefore legitimate to focus on the pseudospin-1 2 subspace Span{|d ± 1 2 }, as has been done in Ref. [47].
When D = (D x , D y , D z ) = 0, however, the evenparity doublet {|d ± 1 2 } is modified. By assuming |D| J, the new eigenstates are given to first order in D by We adopt the lowest-energy doublet {|d ± 1 2 r } at h = 0 as the local basis on each trimer r ∈ A ∪ B, and perform a strong coupling expansion [62] to derive the effective Hamiltonian (see Refs. [63,64] for related calculations for distorted diamond chains). Using this doublet, we introduce a pseudospin-1 2 operator where σ = (σ x , σ y , σ z ) are Pauli matrices. We also introduce a local projection operator P r and a global one P as P r = |d + 1 2 r r d The first-order effective Hamiltonian is derived by projecting the inter-trimer interactions onto the degenerate manifold V 0 = r Span({|d ± 1 2 r }). Using Eq. (10), the projection of the spin operators gives the followings to first order in D: where X r = A, B indicates the set of trimers which r belongs to. Using these relations, the projection of the inter-trimer interactions gives  (17). The vectors u and v are defined in Fig. 2. Only for J3 and J 3 , the relative vector ∆r depends on the sublattice X, as indicated in the corresponding rows. The first-and second-order perturbative estimates (third and fourth columns) are calculated from Eqs. (14) and (18), respectively, for the parameter set (4) obtained by DFT+U .
relative vectors ∆r 1st-order 2nd-order [(∆r, X) for J3, Notably, the effects of the three DM interactions D, D 1 , and D are now combined into the single effective DM interaction D 1 . Furthermore, D 2 has no contribution to first order (as the contributions from the two DM interactions ±D 2X between two trimers cancel out). By using Eq. (5) and setting D = 0 for simplicity, we obtain the estimate The second-order contributions of the exchange couplings to the effective Hamiltonian have been calculated in Ref. [47], and found to have the magnitudes of several Kelvin, which are comparable to |D 1 | above. Taking into account these contributions also, we obtain the full effective Hamiltonian The coupling constants J ∆r,X show seven nonzero different values as listed in Table I, and their second-order expressions are given by Most of these interactions do not depend on the sublattice X, and have the translational invariance of the anisotropic triangular lattice. Only the interactions J 3 , J 3 , and D 1 depend on the sublattice X, and double the unit cell of the effective model; the primitive vectors of the system then change to ±u + v. By using Eqs. (14) and (18) for the parameter set (4), the effective coupling constants are estimated in Table I.
Although the effective coupling constants other than J 1 and J 2 in Table I have comparatively small magnitudes, they can significantly influence the lowtemperature properties of physical quantities and the phase diagram as we discuss in Sections III and IV. To examine such sensitivity to the small effective coupling constants, we allow a slight modification of the J 3 coupling from the second-order expression as This modification is also useful in finding a parameter set consistent with the experimental data of the magnetization and the thermal Hall conductivity, as we discuss in Sec. III E. Microscopically, δJ 3 can arise from the second-neighbor interaction J 2 along the direction of the J bond in the original model-this interaction can be written down as and the first-order perturbation theory gives δJ 3 = (4/9)J 2 . In this section, we analyze the effective pseudospin-1 2 Hamiltonian (17) at high fields by using the spin wave theory. We mainly consider the case in which the field is applied in the −z direction [see Eq. (6a)]. In experiments for single crystals, a broad magnetization plateau at 1 3 of the saturation has been found to appear for H > H c1 27.5 T [35,39]. In the effective model, this corresponds to the fully polarized state of pseudospins. Near but below the 1 3 -plateau phase, the system can be viewed as a finite-density gas of magnons. We consider the temperature regime where these magnons are not condensed or forming bimagnon bound states. We show that magnon Bloch states acquire a nonzero Berry curvature through a combination of the effective Heisenberg and DM interactions, giving rise to a thermal Hall effect [58][59][60][61]65]. We calculate the field-and temperature-dependences of the magnetization and the thermal Hall conductivity. Comparison with experimental data [35,38] indicates the necessity to modify the effective coupling constants significantly from the DFTbased estimates (Table I and Eq. (16)), and we discuss some constraints imposed on these constants. We also discuss the cases of other field directions briefly. We note that our analysis of the thermal Hall effect based on the effective model sharply contrasts with those based on isotropic kagome models [66,67].

A. Spin wave Hamiltonian
We consider the effective pseudospin- 1 2 Hamiltonian (17) with h = (0, 0, h) (h > 0). The plateau state at high fields h corresponds to the fully polarized state of pseudospins r |d + 1 2 r , which we view as the magnon vacuum in the following. Then, T ± r = T x r ± iT y r play the role of magnon annihilation and creation operators, and the magnon occupation number at the site r is given by n r = T − r T + r = 1 2 − T z r . In the following, we replace T + r by the bosonic annihilation operator a r satisfying the commutation relation [a r , a † r ] = δ rr , and introduce an infinite on-site repulsion U 0 → ∞ to impose the hardcore constraint. Furthermore, concerning the DM interactions, we only take into account the z component of the DM vectors to simplify the analysis. This treatment is justified below the saturation field and at low temperatures as the x and y components of the DM vectors vanish in the low-energy effective field theory as we see in Sec. IV.
Our effective Hamiltonian (17) can thus be rewritten in terms of the bosonic operators as where N t is the number of trimers in the system and The first term in Eq. (21) is the energy of the magnon vacuum. The second term is the Zeeman term, which plays the role of a magnon chemical potential. The third and fourth terms, H kin and H int , are magnon kinetic and interaction energies, which do not depend on h. We note that the total number of magnons, r n r , is conserved in this Hamiltonian since we have neglected the x and y components of the DM vectors. We comment that the J 1 and D z 1 terms in H kin can be combined into the form where ν := arg(−J 1 + iD z 1 ). When J 4 = 0, the phase factor e −i Xr ν in Eq. (25) can be removed by performing the gauge transformation a r = e iν a r for r ∈ B while keeping the other terms in the Hamiltonian unchanged. However, such removal of phase factors from H kin is not possible for J 4 = 0 as the J 4 term acquires phase factors by the above gauge transformation. Therefore, the presence of J 4 is crucial for finding nontrivial effects of the complex hopping amplitude of magnons induced by D z 1 .

B. Magnon Bloch states
We first analyze the kinetic part H kin [Eq. (23)] of the spin wave Hamiltonian and determine the Bloch states of magnons. By performing the Fourier expansion where the sum is over the discrete wave vectors in the first Brillouin zone, H kin can be rewritten as Here, the 2 × 2 matrix M (k) is given by where I is the identity matrix, σ = (σ x , σ y , σ z ) are the Pauli matrices, and  (40) and (E1)] for the lower magnon band over the first Brillouin zone. These are calculated for the second-order model (see Table I) with the value of D z 1 in Eq. (16).
with k u = k · u and k v = k · v. The two Bloch energy bands are calculated as It is easy to show that the two bands touch along the entire boundary of the Brillouin zone (|k x | = π/b or |k y | = π/a), which can be interpreted as the Kramers degeneracy due to certain antiunitary symmetries as explained in Appendix B. By parameterizing the vector J (k) in terms of the polar coordinates as the eigenstates corresponding to Eq. (30) are given by Taking account of the second term (Zeeman energy) in Eq. (21) also, we find that the single-magnon Bloch state has an excitation energy E ± (k) + h measured from the vacuum (the 1 3 -plateau state). Introducing E 0 := min k E − (k), we find that the magnon excitation energy is gapped over the entire Brillouin zone for h > h c1 := −E 0 . At h = h c1 , the bottom of the lower magnon band touches zero, leading to a Bose-Einstein condensation of magnons [68,69] at zero temperature. For h < h c1 , repulsion between magnons stabilizes a finite-density condensate of magnons. Since the interaction part H int [Eq. (24)] of the spin wave Hamiltonian has no contribution in the single-magnon problem, h c1 gives the exact single-magnon condensation point.
In systems with competing ferromagnetic and antiferromagnetic interactions, however, multi-magnon bound states are formed under certain conditions and can condense before the single-magnon states do with lowering the field h. A condensation of two-magnon bound states (bimagnons) gives rise to a bond nematic order [49][50][51][52][53][54]. In Ref. [47], it has been found that while the first-and second-order models in Table I show a conventional single-magnon condensation, slightly modified models with reduced J 4 show a condensation of bimagnons. Here we do not address the nature of the low-temperature phase below the 1 3 -plateau phase in further detail in this section. We are instead interested in the thermodynamic behavior of magnons at finite temperatures well above the regime where a single-or twomagnon condensation occurs or magnon bound states are formed.
The lower energy band E − (k) is plotted in Fig. 3(a) for the second-order model (see Table I) with the value of D z 1 in Eq. (16). It shows the minimum value The single-magnon condensation point is therefore given by h c1 = 35.2 K = (gµ B /k B ) × 24.4 T, which reasonably agrees with the low-field end H c1 27.5 T of the 1 3plateau observed experimentally [35][36][37].
For a later purpose, we further consider the expansion of E − (k) around the minima k = (±Q/b, 0): where the coefficients are given by (C x , C y ) = (42.8, 5.12) K for the case of Fig. 3(a). This leads to a constant density of states G = 1/(π C x C y ) = 0.0215 K −1 in units of N t /2 at low energies.

C. Magnetization
To study a finite-density gas of magnons, it is important to treat the magnon interaction part H int in Eq. (24) properly. To this end, we perform the mean-field decoupling of this part, which results in the effective chemical potential µ = −(h − h c1 ) − 2U n for magnons [69]. Here, n is the magnon number per trimer, and U is the effective interaction parameter which encompasses the effects of all the interaction terms in Eq. (24). Furthermore, µ is measured relative to E 0 = −h c1 so that the condensation occurs at µ = 0. Since there is a subtlety in the mean-field treatment of the infinite on-site interaction U 0 , it is challenging to determine U through a microscopic calculation. We instead determine U later in such a way that consistency with the experimental magnetization data [35] is achieved. The magnon density n is obtained as a function of µ and T as where the division by 3 comes from the fact that n is defined per trimer. The magnetic field H is related with the chemical potential µ and the magnon density n as For given T and U , we can thus obtain the theoretical magnetization curve around H = H c1 by calculating n(µ, T ) in Eq. (34) as a function of µ < 0 and plotting the relation between Eqs. (35) and (36). Although the values of H c1 are slightly different between theory and experiment (as noted in Sec. III B), the magnetization curves as a function of H − H c1 can be compared in a quantitative manner. In Fig. 4(a), we compare the theoretical and experimental magnetization curves for T = 1.4 K and 4.2 K. In the absence of an interaction (U = 0), the magnon density n ∝ −∆M diverges with lowering the field H to H c1 . A finite magnon density n ∝ −∆M for H < H c1 in the experimental data can therefore be interpreted as a result of magnon repulsion. We determine the effective interaction parameter U in such a way that the slope of the T = 1.4 K experimental magnetization curve slightly below the 1 3 -plateau is reproduced (see Appendix C). With U = 14.0 K, we indeed find good agreement between theoretical and experimental curves for −4 T < H − H c1 < −2 T and T = 1.4 K. For −2 T < H − H c1 or T = 4.2 K, however, the magnon density n ∝ −∆M tends to be larger in experiment than in theory. This indicates that the density of states of magnons from the 1 3 -plateau state in volborthite is larger than that expected in the effective model with the present parameter values. This simply implies that volborthite is more frustrated than our present model.
We can indeed estimate the density of states of magnons by fitting the experimental data using the asymptotic form which is derived in Appendix C. Here, G is the density of states of magnons in units of N t /2 in the low-energy limit.
The best fit as shown in Fig. 4(b) gives G = 0.045 K −1 , which is roughly twice as large as the value G = 0.0215 K −1 for the model used in Fig. 4(a). This indicates the necessity to modify the effective coupling constants {J ∆r,X } to have consistency with experiment, as we discuss in more detail in Sec. III E.

D. Thermal Hall conductivity
The thermal Hall conductivity in the clean limit is given by [59][60][61] where 2N c is the number of layers in the system, and V = (N t /2)ab · N c c sin β is the volume of the system (we remind the reader that two inequivalent layers are alternately stacked in the c direction in volborthite). Here, we assume that inter-layer interactions can be neglected and the contributions from different layers can simply be summed up, hence the multiplication by 2N c in Eq. (38). The function c 2 (ρ) is given by where Li m (z) = ∞ n=1 z n /n m is the polylogarithm function; see the inset of Fig. 5(d) for a plot of c 2 (ρ). The Berry curvature Ω α (k) (α = ±) is defined as where ∂ i = ∂ ∂ki (i = x, y) and ij is an antisymmetric tensor with xy = − yx = 1.
The dimensionless Berry curvature Ω − (k)/(ba) for the lower band is plotted in Fig. 3(b). It shows the maximal value Ω − (k)/(ba) 3 × 10 −4 near k = ±k * := ±(π/(2b), 0), and vanishes at the Brillouin zone boundary. At k = k * , the following simple expression of the Berry curvature is available (see Appendix E for the derivation): This expression suggests that the sign and the magnitude of the Berry curvature are controlled not only by the DM interaction D z 1 but also by the long-range effective couplings J 3 − J 3 and J 4 which arise from the second-order strong-coupling expansion. Furthermore, the magnitude of the Berry curvature in Eq. (41) depends significantly on the nearest-neighbor effective coupling J 1 .
The importance of the effective coupling constants discussed above for the magnitude of the Berry curvature can be understood as follows. In order to have nonzero Ω ± (k), the vector J (k) in Eq. (29) must form a solid angle as we vary k in different directions. When D z 1 = 0 or J 3 − J 3 = 0, one component of J (k) vanishes and J (k) is thus constrained to a 2D plane. When J 4 = 0, we have J x (k) ∝ J y (k), and J (k) is again constrained to a 2D plane. Therefore, D z 1 , J 3 − J 3 , and J 4 must all be nonzero to have nonzero Ω ± (k) and κ xy . The importance of these coupling constants can also be understood from the following viewpoints: (i) When J 3 −J 3 = D z 1 = 0, the effective model acquires the invariance under the translations by u and v, and no longer satisfies minimal requirement of a two-sublattice structure for a finite Berry curvature; (ii) When J 4 = 0, the complex hopping amplitude of magnons due to D z 1 can be transformed into a real one by a gauge transformation as discussed in Sec. III A. We note that the dimensionless Berry curvature Ω − (k)/(ba) in Fig. 3(b) takes only small values of order 10 −4 over the entire Brillouin zone because the vector J (k) changes only around the −x direction owing to −J 1 |J 3 − J 3 |, |J 4 |, |D z 1 | in the present model. To evaluate κ xy , it is useful to rewrite Eq. (38) as where the coefficient is evaluated as We note that κ xy /T is an antisymmetric function of the magnetic field H: κ Let us now compare our results with the experimental data of Watanabe et al. [38] for |H| 15 T. Although the experimental field range is away from the 1 3plateau, it is expected that magnitudes and qualitative features of the thermal Hall conductivity do not change abruptly as we vary |H| at sufficiently high temperatures T 2 K where there is no phase transition. The appearance of a peak in the T dependence of |κ xy (−H)|/T in Fig. 5(d) is qualitatively consistent with the experimental data for H = 15 T. However, the experimental peak value κ xy (−H)/T ∼ 4×10 −5 W·K −2 ·m −1 is two orders of magnitudes larger than our data in Fig. 5(c,d). Furthermore, our data in Fig. 5(c,d) always show κ xy (−H) < 0 while the experimental data show κ xy (−H) > 0 for T 4 K. We note that the sign and the magnitude of the thermal Hall conductivity depend crucially on various effective coupling constants as seen in the representative value of the Berry curvature in Eq. (41). Specifically, the magnitude of the right-hand side of Eq. (41) significantly increases as we weaken J 1 from the present estimate. We discuss this issue further in Sec. III E.
In passing, we note that in general, the integral of the Berry curvature over the Brillouin zone for an isolated band, N α := BZ 2π Ω α (k), is topologically quantized to integers [71]. In the present case, however, the two bands touch at the Brillouin zone boundary, and thus the quantization rule does not apply individually for each band. Indeed, in Fig. 3(b), Ω − (k)/(ba) shows positive but small values over the Brillouin zone, and the full integration of it gives only a tiny non-integer value N − = 7.5 × 10 −4 . This contrasts with the case of SrCu 2 (BO 3 ) 2 , for which a theory predicts the emergence of topologically nontrivial bands and an associated large thermal Hall conductivity through a combination of DM interactions and a magnetic field [72].

E. Modification of the coupling constants
We have seen that the magnetization and the thermal Hall conductivity calculated from our effective model can capture some qualitative features of the experimental data [35,38], but there are some quantitative discrepancies. This indicates the necessity to modify the effective coupling constants, {J ∆r,X } and D 1 , from the DFTbased parameter set in Table I and Eq. (16). In tuning these coupling constants, the following observations provide a useful guidance. Firstly, the magnetization data have indicated that the low-energy density of states of magnons, G, must be roughly twice as large as the value for the DFT-based parameter set (see Fig. 4). Secondly, the experimental data of κ xy for H = 15 T [38] are two orders of magnitude larger than our data for H ≈ 20-35 T in Fig. 5; if we assume a smooth change of κ xy between these values of H, the Berry curvature Ω − (k) must be two orders of magnitude larger than our calculated data   Fig. 3(b). Lastly, the experimental value H c1 = 27.5 T [35,37] of the low-field end of the 1 3 -plateau poses another constraint. We can therefore summarize the quantitative requirements for the effective model as follows: Here, for the Berry curvature Ω − (k), we take the representative value at k = k * , for which a simple analytical expression (41) is available. The low-energy density of states, G, can be determined by the curvature (C x and C y in the x and y directions as in Eq. (33)) of the dispersion relation around the minima. We note that for a given parameter set, accurate values of h c1 , C x , and C y can be obtained by numerically minimizing the dispersion relation in Eq. (30); however, to see how they depend on the effective coupling constants, the analytical (yet approximate) expressions in Eq. (D3) in Appendix D (valid for 2J 2 > −J 1 others) are useful. To satisfy the requirements in Eq. (44), we modify the coupling constants in the following way. To enlarge Ω − (k * )/(ba) by two orders of magnitude, one must reduce the magnitude of the effective coupling J 1 consider-  (46). For the effective interaction parameter U = 14.0 K tuned appropriately [see Appendix C], the theoretical curves show a good agreement with the experimental data [35]. The curves for the non-interacting case U = 0 are also presented for comparison.
ably, which is achieved by reducing |J 1 | or enlarging |J | from the DFT+U estimates in Eq. (4). Since the reduction in |J 1 | leads to enlargement in h c1 , we must reduce J 2 (and thus J 2 ) at the same time to satisfy Eq. (44a). We also introduce a modification δJ 3 to the J 3 coupling as in Eq. (19) to change the sign of the Berry curvature to negative as required in Eq. (44c). After some examination, we have arrived at the following modified parameter set:  Fig. 6(a) shows minima at k = ±(Q/b, 0) with Q/(2π) = 0.473. The Berry curvature Ω − (k * )/(ba) in Fig. 6(b) is negative over the entire Brillouin zone, and have much larger amplitudes than in Fig. 3(b). Reflecting this, the thermal Hall conductivity in Fig. 8(c,d) is positive and have much larger amplitudes than in Fig. 5(c,d). Furthermore, the calculated magnetization curves show

F. Other field directions
We have so far considered the case in which the field H is applied in the −ẑ direction [Eq. (6a)]. Here we briefly discuss the case of other field directions.
When the field H is changed to the −ŷ direction, the calculation of κ xy can be done in parallel with the case of H −ẑ discussed above by just replacing D z 1 by D y 1 . Therefore, by measuring κ xy for two different directions of the field (−ŷ and −ẑ), one can determine the relative sign and magnitude of D y 1 in comparison with D z 1 . When the field H is changed to the −x direction, the system restores the 2 1 screw axis symmetry. In this case, we have κ xy = 0 because the heat current in the y direction changes its sign under the 2 1 operation.

IV. EFFECTIVE FIELD THEORY FOR A QUASI-ONE-DIMENSIONAL REGIME
In this section, we analyze the low-temperature phases of the effective pseudospin- 1 2 Hamiltonian (17) by means of a field-theoretical method. The key point of our approach is the anisotropic triangular lattice structure of the effective model as seen in Fig. 2(b). Among various effective couplings in Table I, J 2 on the horizontal bonds has the largest magnitude. If we look only at the J 2 coupling and the Zeeman term, the system can be viewed as an array of decoupled Heisenberg chains in a magnetic field. In this case, the low-energy description of each chain is given by the Tomonaga-Luttinger liquid (TLL) theory unless the pseudospins are fully polarized by the field. We can then include the other couplings in the effective Hamiltonian, which couple different chains and give rise to a variety of magnetic orders. We analyze the competition among those inter-chain couplings through a perturbative renormalization group (RG) method and the chain mean field theory.
To facilitate the quasi-one-dimensional viewpoint, it is useful to represent the positions r ∈ A ∪ B using the x and y coordinates. The J 2 and Zeeman parts of the effective Hamiltonian can then be rewritten as Here, the coordinates y ∈ Z and x ∈ Z + y/2 are measured in units of a/2 and b, respectively. Even (odd) y's correspond to the sites in A (B). Interactions between nearest-neighbor chains are given by the J 1 and J 4 couplings and the effective DM interaction, which are expressed as Interactions between next-nearest-neighbor chains are given by the J 2 , J 3 , and J 3 couplings, which are combined into the form Throughout this section, we do not consider the J 5 coupling as it only leads to a slight modification of the field-theoretical parameters (such as the velocity v and the compactification radius R) of each chain. Our analysis in this section is closely analogous to those in Refs. [26,27,73], in which spatially anisotropic triangular antiferromagnets have been studied in relation with Cs 2 CuCl 4 . Therefore, we take some notations similar to these references, refer to some of their results, and adapt them to the present problem. Some doubt may be cast on the applicability of the coupled-chain approach to volborthite since J 1 has a magnitude comparable to J 2 in the estimates in Table I. However, numerical studies have indicated that predictions of the coupled-chain approach can qualitatively continue up to rather large values of inter-chain couplings. Examples include the SDW and cone phases in spatially anisotropic triangular antiferromagnet [73] and the dimer and vector chiral phases in the J 1 -J 2 XXZ chain [74][75][76][77][78]. Furthermore, in the present problem, |J 1 | can potentially be modified into a much smaller value as we have discussed in Sec. III E.

A. Field theory for a Heisenberg chain
Here we briefly summarize the field-theoretical (bosonized) description of a spin-1 2 antiferromagnetic Heisenberg chain in a magnetic field, which corresponds to a part of Eq. (47) with fixed y. The Hamiltonian is given by where x runs over integers or half-integers, and S x is the spin-1 2 operator at the position x. The magnetic field h is chosen to be in the z direction. The magnetization m = (1/L) x S z x = S z x , with L being the number of spins in the chain, is conserved in this Hamiltonian.
For any m less than the saturation, i.e., m ∈ (−1/2, 1/2), the low-energy description of Eq. (50) is given by the TLL theory with the Hamiltonian Here, the bosonic fields φ(x) and θ(x) satisfy the commu- is the Heaviside step function (with Θ(0) = 1/2). The field φ(x) is compactified on a circle of radius R, and θ(x) is analogously compactified with the radius 1/(2πR); the vertex operators appearing in Eq. (53) below are consistent with this compactification. The velocity v and the compactification radius R are smooth functions of m. The magnetization curve m(h) and the dependences of v/J and R on m can be determined by numerically solving the integral equations obtained from the Bethe ansatz [79][80][81]. In particular, the exponent η(m) = 2πR(m) 2 monotonically decreases from η(0) = 1 to η(1/2) = 1/2 with the increase in m; see Fig. 10(c).
At a fixed magnetization m, the low-energy fluctuations of spins occur around particular wave vectors k. Such wave vectors are k = 0 and π ± 2δ with δ = πm for the "longitudinal" spin component S z x along the field direction, and k = ±2δ and π for the "transverse" spin components S ± x = S x x ± iS y x perpendicular to the field. The spin operators can thus be decomposed as Here, the operators S µ k (x) (µ = z, ±) describe slowly varying components, and are expressed in terms of the bosonic fields as The coefficients A 1 , A 2 , and A 3 are also smooth functions of m, and have been determined numerically in Ref. [82]. Among the operators in Eq. (53), S z ∓π±2δ and S ± π with comparatively smaller scaling dimensions ∆ SDW = 1/(2η) and ∆ ± = η/2, respectively, can play particularly important roles in the low-energy physics. The other operators, S z 0 and S ± 2δ , have scaling dimensions 1 and η/2 + 1/(2η), respectively.

B. Field theory for coupled chains
Let us now consider the coupled-chain problem with Eqs. (47), (48), and (49). Here we present the fieldtheoretical expressions of various inter-chain couplings, and discuss their relevances to the low-energy physics on the basis of their scaling dimensions. More detailed discussions on the competition among those couplings are given in Sec. IV C.

Magnetic field in the z direction
We first consider the case in which the field h is applied in the z direction. Using Eq. (52) for each chain labeled by y, the pseudospin operators are expressed as T z x,y = m + S z y;0 (x) + e i(−π+2δ)x S z y;−π+2δ (x) + e i(π−2δ)x S z y;π−2δ (x), T ± x,y = e i2δx S ± y;2δ (x) + e −i2δx S ± y;−2δ (x) + e ±iπx S ± y;π (x) (54) with m = (1/N t ) x,y T z x,y and δ = πm. The operators S µ y;k (x) (µ = z, ±) can be expressed in terms of the bosonic fields φ y (x) and θ y (x) defined on each chain, as in Eq. (53). The J 1 and J 4 couplings between nearest-neighbor chains in Eq. (48) are then expressed as with As discussed by Starykh et al. [26,27,83] (and numerically demonstrated in Ref. [73]), the interactions in Eq. (55) lead to competition between SDW and cone orders. Namely, the γ SDW term with the scaling dimension 2∆ SDW = 1/η induces the incommensurate SDW order at low fields while the γ cone term with the scaling dimension 1 + 2∆ ± = 1 + η induces the incommensurate cone order at high fields (below the saturation); see Fig. 10(c) for the plots of these scaling dimensions. The interaction S z y;0 S z y+1;0 with the scaling dimension 2 is marginal for any m ∈ (−1/2, 1/2). The ellipsis in Eq. (55) indicates other terms which have larger scaling dimensions and are less important in the low-energy physics.
The simple scenario of the SDW-cone competition in Eq. (55) may break down if we also consider other couplings in the effective Hamiltonian. Specifically, the J 2 , J 3 , and J 3 couplings between next-nearest-neighbor chains in Eq. (49) are expressed as with The γ SDW and γ SDW terms have the same scaling dimension as the γ SDW term while the γ cone term has the same scaling dimension as the γ cone term. The γ term has a smaller scaling dimension 2∆ ± = η, and grow much faster than the other terms in the RG flow; when γ > 0 (γ < 0), it has the effect of stabilizing (destabilizing) the cone order induced by the γ cone term. However, whether this term dominates the low-energy physics depends on the initial (bare) value of γ , and this issue is analyzed in more detail in Sec. IV C.
For m = 0, the effective DM interaction in Eq. (48) is expressed as We note that the x and y components of the DM interactions disappear in a perturbative treatment because T ± x,y and T z x,y have Fourier components with separated wave vectors for m = 0 as seen in Eq. (54). The interaction in Eq. (59) has the form similar to the γ cone term in Eq. (55). In fact, the two terms can be combined as where we definē S ± y;π (x) := exp ±i γ cone := γ 2 cone + (D z 1 /2) 2 , ν := arg (γ cone + iD z 1 /2) .
We note that the other terms in Eqs. (55) and (57) remain unchanged under the "gauge transformation" of S ± y:π (x) done here. Therefore, the effects of the D z 1 term are to enlarge the amplitudeγ SDW of the cone-inducing term and to modify the resulting cone structure slightly via the gauge transformation. The gauge transformation introduced here is analogous to the one discussed in the spin wave analysis; see the last paragraph of Sec. III A. Contrary to that case, the condition J 4 = 0 is not required in the present discussion. This implies that in the low-energy theory, the thermal Hall conductivity κ xy vanishes even for J 4 = 0.
At sufficiently low fields, the x and y components of the DM interaction can also play certain roles in the lowenergy physics because the shift δ = πm of momenta in Eq. (54) vanishes as m → 0. For a better understanding of this regime, it is useful to consider the case of precisely zero field, i.e., the case of m = 0, as we do next.

Zero magnetic field
For m = 0, Eq. (54) can simply be written as which is based on the following mapping: S z y;0 → M z y , e −i π 2 y S z y;−π+2δ + e i π 2 y S z y;π−2δ → N z y , S ± y;2δ + S ± y;−2δ → M ± y ≡ M x y ± iM y y , e ±i π 2 y S ± y;π → N ± y ≡ N x y ± iN y y .
Here, the uniform and staggered components, M y (x) and N y (x), have the scaling dimensions 1 and 1/2, respectively. The effective DM interactions can then be expressed as The term N b y N c y+1 has the smallest scaling dimension 1, and grows fastest in the RG flow. If this term dominates the low-energy physics, an "orthogonal" order in which spins rotate by ±90 • in the yz plane (Fig. 9) appears. Once this order appears at zero field, it is expected to persist in the low-field regime.

Magnetic field in the x direction
When the field is applied in the x direction, we can perform the same line of analysis as in Sec. IV B 1 by expressing T x x,y andT ± x,y := T y x,y ± iT z x,y in terms ofS z y;k (x) and S ± y;k (x), respectively, in a way analogous to Eq. (54). The only difference occurs in the expression of the effective DM interaction H DM . Specifically, for m := S x x = 0, only the x component of the DM interaction remains in a perturbative treatment, and it is expressed as This interaction is a finite-field version of the N b y N c y+1 term in Eq. (64), and has a small scaling dimension 2∆ ± = η. It thus grows as fast as the γ term in Eq. (57) along the RG flow, and potentially dominates the low-energy physics over the entire range of the magnetic field below the 1 3 -plateau. If this happens, a canted orthogonal order in which pseudospins T x,y have the constant magnetization m in the x direction and rotate by ±90 • in the yz plane (as in Eq. (82) below) appears up to the 1 3 -plateau. Unfortunately, in 51 V NMR measurements as were done in Refs. [35,36], this order would show no direct signal because of the cancellation of the internal fields at the V site as seen in Fig. 9.  13) and (82) with Θ0 = 0. This order is induced by the D x 1 interaction in Eq. (64) and formed in the yz plane. When a magnetic field is applied in the x direction, the spins further acquire nonzero averages, S x r,1 = S x r,3 ≈ 2 3 m and S x r,2 ≈ − 1 3 m in the x direction, leading to a canted order. In 51 V NMR measurements as were done in Refs. [35,36], this order would show no direct signal because the internal fields at the V site from the surrounding spins on the hexagon cancel out except the uniform component in the x direction.

C. Chain mean field theory
We now quantitatively analyze the competition among the inter-chain couplings which are described in Sec. IV B. Specifically, following Ref. [27], we calculate the critical temperatures associated with different magnetic orders using the chain mean field theory. The order with the highest critical temperature is expected to be selected among the competition. We first summarize our results in Sec. IV C 1, and then describe the details of the analysis in the subsequent sections. The processes of calculations go essentially the same way as in Appendix D of Ref. [27], and we roughly describe the ideas and adapt their results to the present model.

Summary of the results
We consider the problem using the modified parameter set in Eqs. (45) and (46). If we directly use this parameter set, we have a rather large value (γ = 2.5 K) of the relevant γ coupling, and the γ term dominates the lowenergy physics over the entire range of the magnetic field below the 1 3 -plateau for h ẑ. To obtain a rich phase diagram as observed experimentally [35][36][37] (and shown (Color online) Critical temperatures associated with the SDW, incommensurate cone, and canted orthogonal orders, determined by the chain mean field theory for (a) h ẑ and (b) h x. We use the modified parameter set in Eqs. (45) and (46), but change the value of δJ3 to δJ3 ∈ {3.6, 3.8, 4.0} K, which correspond to γ ∈ {0.1, 0.3, 0.5} K, respectively. In (a), the critical temperature for the orthogonal order (dashed horizontal lines at T ∼ 3 K) is estimated from the data for smallest m in (b); the stability of this order against the change in m is unfortunately beyond the scope of the present analysis. (c) Scaling dimensions of the inter-chain couplings that induce the three types of orders. The γSDW and γcone terms in Eq. (55) have the scaling dimensions 2∆SDW and 1 + 2∆±, respectively, whose crossing at m ≈ 0.3 leads to the scenario of the SDW-cone competition. The γ and D x 1 terms in Eqs. (57) and (65) have the smaller scaling dimension 2∆±, and substantially change the magnetic phase diagrams unless their magnitudes are suppressed.
The calculated critical temperatures are displayed in Fig. 10(a,b); see also Fig. 1 for the resulting phase diagrams at zero temperature. We first look at the result for (a) h ẑ. For γ = 0.1 K, which is sufficiently small, we find that the scenario of the SDW-cone competition in Eq. (55) essentially holds: as we lower the temperature, the SDW and cone orders first set in for m 0.3 and m 0.3, respectively. For γ = 0.3 and 0.5 K, in contrast, the cone order is stabilized significantly by the γ coupling, and it wins against the SDW order over the full range of the magnetization m. At sufficiently small m, the orthogonal order induced by the D x 1 interaction sets in even before the cone order does, as indicated by dashed horizontal lines at T ∼ 3 K in Fig. 10(a). We also note that close to the saturation m = 1/2, the bond nematic order appears due to condensation of bimagnons for a certain range of parameters [47]; unfortunately, we are not aware of an appropriate method for describing this order within bosonization for the present anisotropic triangular system.
Our result indicates that in order to obtain the SDW order over an extended range of the field as observed experimentally [35,36], the following constraint is required on the value of γ : Unfortunately, we have not been able to find a parameter set which simultaneously satisfies the requirements in Eqs. (44) and (66). We further note that the range of γ in Eq. (66) is too narrow to be satisfied in a realistic system. Since our analysis is performed in the limit of weak inter-chain couplings, we expect that the constraint in Eq. (66) is loosened with the increase in the magnitudes of the inter-chain couplings in order to be consistent with experiment. This indicates that a nontrivial stabilization mechanism of the SDW order exists beyond the scope of the perturbative RG approach. We also note that defects present in volborthite crystals [19] can have nontrivial effects on the stabilization of the SDW order-as the SDW state breaks only the translational symmetry, such defects can act as random fields on the SDW order parameter, as is known in the context of collective pinning [84] (for further discussion, see Sec. V A of Ref. [83]).
We next look at the result for (b) h x. In this case, the canted orthogonal order first sets in owing to relevant D x 1 over the full range of m below the saturation.
Our results indicate that the magnetic phase diagrams depend sensitively on the direction of the magnetic field, as summarized in Fig. 1(b,c). The experimental investigations of the phase diagram [35][36][37] have mainly been conducted for magnetic fields perpendicular to the kagome plane, i.e., h ẑ. It will be interesting if these experiments are extended to other directions of the magnetic field to uncover nontrivial roles of the DM interactions as predicted here.

SDW phase
The inter-chain couplings related to a magnetic ordering in the longitudinal component T z x,y are summarized as For γ SDW > 0 and γ SDW = γ SDW = 0, the ground state of this Hamiltonian is clearly given by φ y (x) = φ 0 (constant). Small γ SDW and γ SDW would not modify this ground state since the expansion of (67) around this state does not produce any term linear in φ y 's. By setting e iφy/R =ψ = |ψ|e iΦ0 and performing the mean-field decoupling of the inter-chain couplings, we obtain The resulting state is the incommensurate SDW order in the longitudinal component along the magnetic field: Combining Eq. (68) with the unperturbed Hamiltonian, which is the TLL theory (51) for each chain, we obtain a collection of sine-Gordon models decoupled into different chains. We can then calculate the finite-temperature average e iφy/R perturbatively in powers ofψ, obtaining the self-consistent equation. The leading-order result is

Cone phase
We next discuss magnetic orderings in the transverse components T ± x,y , first focusing on the case of h ẑ. In this case, the effective DM interaction can be treated by the suitable gauge transformation of S ± y;π (x) in Eq. (61). Reflecting this, we introduce the shifted fieldθ y (x) as The inter-chain couplings related to a magnetic ordering in the transverse component are then summarized as As first pointed out by Nersesyan et al. [75], the interaction like theγ cone term here leads to an incommensurate transverse order. To describe such an order, it is useful to set 2πRθ y = −q 0 x + 2πRθ y , where the first term on the r.h.s. describes the incommensurate rotation of spins and the second the slowly varying component. The inter-chain couplings in Eq. (73) are then rewritten as − γ cos 2πR θ y −θ y+2 + 2(−1) y q 0 γ cone sin 2πR θ y −θ y+2 .
By further setting e i2πRθy =ψ = |ψ|e iΘ0 and performing the mean-field decoupling, we obtain The resulting state is the incommensurate transverse order (the cone order) with We can determine the associated critical temperature T cone in a way similar to Sec. IV C 2. However, reflecting the transformation fromθ y toθ y , the susceptibility should be evaluated at the wave vector q 0 . The condition for the critical temperature T cone is then given by The wave vector q 0 is determined in such a way as to maximize T cone . In this way, a set of implicit equations for determining T cone and q 0 are obtained as (see Eq. (D23) in Ref. [27]) 4Im Ψ ∆ ± 2 + ir = 2π sinh(2πr) cosh(2πr) − cos(π∆ ± ) + λ cone s λ cone sr + λ , where we introduce When D x 1 is sufficiently weak, a similar cone order can also appear for h x. Such a case can be analyzed by setting D z 1 → 0 in the above argument.
The inter-chain coupling related to this order is then summarized as Here, we did not include the γ cone and γ cone terms as they vanish after the mean-field treatment. By setting e i2πRθy =ψ = |ψ|e iΘ0 and performing the mean-field decoupling, we obtain The resulting state is the commensurate transverse order (the canted orthogonal order) with T ± x,y = T y x,y ± iT z x,y = A 3 |ψ|e ±i[π(n+x)+Θ0] . (82) The condition for the critical temperature T orth is given by which is independent of Θ 0 . By solving this equation, the transition temperature is obtained as This is similar to the second equation in Eq. (77) but with r = 0 because of the commensurate nature. We note that Θ 0 in Eq. (82) should in the end be fixed at a certain value as the effective spin model (17) does not possess a spin rotational symmetry around any axis owing to the DM interactions; unfortunately, the value of Θ 0 cannot be determined by the present mean field approach.

V. SUMMARY AND OUTLOOK
In this paper, on the basis of the coupled-trimer model of Ref. [47], we have investigated the effects of DM interactions on the magnetic properties of volborthite. By means of a strong-coupling expansion, we have derived an effective pseudospin-1 2 model on an anisotropic triangular lattice. In the effective model, the magnetic anisotropy is characterized by a single effective DM vector D 1 (in contrast to four vectors in the original model), which leads to a significant simplification of our analysis. We have performed a spin wave analysis starting from the 1 3 -plateau state for the case of magnetic fields perpendicular to the kagome layer. The magnon Bloch states have been found to acquire a nonzero Berry curvature, which gives rise to a thermal Hall effect. Our magnon Bose gas theory can explain qualitative features of the magnetization and the thermal Hall conductivity measured experimentally. Through a further quantitative comparison with experiment, we have derived some constraints on the effective model as in Eq. (44). In particular, the requirement of enlarging the Berry curvature by two orders of magnitude leads to a much smaller magnitude of the J 1 coupling, promoting a quasi-one-dimensional picture. Based on this picture, we have analyzed magnetic orders at low temperatures using effective field theory. The requirement that the SDW phase appear for an extended range of the magnetic field poses the constraint (66) on the magnitude of the relevant γ coupling between the second-neighbor chains. Assuming this, we have predicted the magnetic phase diagrams as schematically shown in Fig. 1(b,c), which sensitively depend on the field direction. Unfortunately, we have not been able to find a parameter set which simultaneously satisfy the constraints in Eqs. (44) and (66), leaving open the issue of more precise determination of the microscopic spin model of volborthite.
Our analysis of the thermal Hall effect has been focused on the regime just below the 1 3 -plateau, where the system can be described as a low-density gas of magnons. This approach is less effective with lowering the magnetic field as the mutual interactions between magnons become more significant. It is worth noting that the pseudospin-1 2 effective model on an anisotropic triangular lattice has a structure similar to that of Cs 2 CuCl 4 , and may support fractionalized excitations such as spinons, psinons, and antipsinons in such an intermediate-field regime [28,29]. It is an interesting theoretical challenge to calculate the thermal Hall conductivity based on those fractionalized excitations. Such a calculation can be directly compared with the thermal Hall conductivity data up to 15 T of Watanabe et al. [38]. If experimental measurements can be extended to higher fields, it would provide an exciting possibility of observing the crossover from fractionalized excitations to magnons through transport properties. It would also be interesting to investigate the role of magnon bound states, which appear below 1 K around the low-field end of the 1 3 -plateau, on transport properties.
The experimental investigations of the magnetic phase diagram [35][36][37] have mostly been conducted for the case of H z as shown in Fig. 1(a). Furthermore, the nature of Phase I has yet to be explored in single crystals. We expect that the prediction of a crucial dependence on the field direction and the characterization of different phases in Fig. 1(b,c) in this work stimulate further experimental studies. The nature of the two-step transition to Phase I with decrease in temperature [36,37,43] also merits further investigation in both theory and experiment. Here we discuss how the symmetry of the space group P2 1 /a (No. 14) imposes constraints on the DM vectors as shown in Fig. 2.
We therefore find the appearance of the DM interactions −D and −D 1 on the respective bonds.
In these ways, we have constraints on the relative signs of the DM vectors as in Fig. 2(a).
Similar symmetry consideration also applies to the DM vectors in the effective model. It leads to the relative signs of the DM vectors on the J 1 bonds as shown in Fig. 2(b). Furthermore, one can show that the DM interactions on the J 2 bonds strictly vanish. This is because such a DM interaction, if present, is mapped onto the DM interaction with the reversed sign on the same bond under the site-centered inversion followed by translation.

Appendix B: Magnon band touching
Here we argue that the touching of the two magnon bands (30) at the Brillouin zone boundary can be understood as the Kramers degeneracy due to certain antiunitary symmetries (see Refs. [85][86][87][88] for related arguments for other space groups).
We place the origin of the coordinate (x, y, z) at the center of a trimer of type A. We introduce the 2 1 screw axis operation C 2x about the axis (y, z) = (−a/4, 0), the inversion I about the origin, and time reversal Θ. Under these operations, the coordinate, the momentum, and the spins are transformed as follows: k →k, S r,j →S r ,j ; I : r → −r, k → −k, S r,j → S −r,4−j ; Θ : r → r, k → −k, S r.j → −S r,j .
We note that C 2x and I are unitary while Θ is antiunitary.
In the absence of a magnetic field H, the Hamiltonian has the symmetries under all of the three operations C 2x , I, and Θ. In the presence of a magnetic field H in the z direction, the symmetries under C 2x and Θ are lost while that under I is retained. Yet, the Hamiltonian is still symmetric under the following product of operations: It is also useful to consider the following product, which also leaves the Hamiltonian invariant: k →k, S r,j → −S −r ,4−j .
Since (C 2x Θ) 2 and (IC 2x Θ) 2 are equal to the translations by b and a, respectively, we have (C 2x Θ) 2 = e ik·b and (IC 2x Θ) 2 = e ik·a in the subspace with the wave vector k. Thus, the Kramers degeneracy due to the antiunitary symmetry C 2x Θ occurs when e ik·b = 1 and k is invariant under this operation-this explains the band touching for k x = π/b. Similarly, the Kramers degeneracy due to the antiunitary symmetry IC 2x Θ occurs when e ik·a = 1 and k is invariant under this operation, i.e., when (k·a, k·c) = (π, 0) and (π, π); since inter-layer couplings are neglected in our present model, this degeneracy occurs for arbitrary k · c, explaining the band touching for k y = π/a. Henceforth, we assume 0 < −J 1 < 2J 2 and J 2 1 > 4J 2 J 2 . When J 3 +J 3 = J 4 = J 5 = 0, E − (k) above is minimized at k = (±Q 0 /b, 0) with When J 3 + J 3 , J 4 , and J 5 are finite but their magnitudes are sufficiently smaller than |J 1 |, we can expect that these minimum points change only perturbatively. We can therefore search for the minima of Eq. (D1) by expanding it around k = (±Q 0 /b, 0). The resulting expression is Eq. (33) with Q = Q 0 +δQ, and the first-order perturbative estimates of h c1 = −E 0 , δQ, and (C x , C y ) are obtained as For the second-order model (see Table I), we have h est c1 = 34.3 K, Q est /(2π) = 0.382, (C est x , C est y ) = (41.2, 5.76) K; these agree reasonably with the accurate values h c1 = 35.2 K, Q/(2π) = 0.369, (C x , C y ) = (42.8, 5.12) K given in Sec. III B.