Thermodynamics of anisotropic triangular magnets with ferro- and antiferromagnetic exchange

We investigate thermodynamic properties like specific heat $c_{V}$ and susceptibility $\chi$ in anisotropic $J_1$-$J_2$ triangular quantum spin systems ($S=1/2$). As a universal tool we apply the finite temperature Lanczos method (FTLM) based on exact diagonalization of finite clusters with periodic boundary conditions. We use clusters up to $N=28$ sites where the thermodynamic limit behavior is already stably reproduced. As a reference we also present the full diagonalization of a small eight-site cluster. After introducing model and method we discuss our main results on $c_V(T)$ and $\chi(T)$. We show the variation of peak position and peak height of these quantities as function of control parameter $J_2/J_1$. We demonstrate that maximum peak positions and heights in N\'eel phase and spiral phases are strongly asymmetric, much more than in the square lattice $J_1$-$J_2$ model. Our results also suggest a tendency to a second side maximum or shoulder formation at lower temperature for certain ranges of the control parameter. We finally explicitly determine the exchange model of the prominent triangular magnets Cs$_2$CuCl$_4$ and Cs$_{2}$CuBr$_{4}$ from our FTLM results.


I. INTRODUCTION
The 2D triangular S = 1/2 Heisenberg antiferromagnet (HAF) is the archetypical geometrically frustrated quantum magnet [1]. Therefore it has been studied in countless investigations. Although there is no pure physical realization due to in-plane symmetry breaking and inter-plane coupling several compounds fall into the wider class of anisotropic triangular magnets where the exchange coupling J 2 along one of the triangle sides is different from J 1 along the other two sides [2]. In fact this generalization is not a drawback but opens interesting possibilities. It allows to consider a smooth interpolation between square lattice HAF to isotropic triangular magnet to quasi-1D chain compounds as function of a single tuning parameter J 2 /J 1 .
Most of the work on this magnetic system was focused on the zero temperature phase diagram as function of anisotropy ratio, spin wave excitations and influence of quantum fluctuations as well as effects of external fields like plateau formation in the magnetization. We will not touch these fundamental topics but refer to Ref. [2] and references cited therein. Here we will focus on a more practical issue, namely the determination of the exchange constant or at least their anisotropy ratio from thermodynamic quantities of triangular magnets. This can be done by analyzing the experimental temperature dependence of specific heat and magnetic susceptibility where the former is less favorable due to lattice contributions [3]. A simplified but not unambiguous determination of exchange constants is possible by identifying position and height of the maximum in thermodynamic quantities and comparing with theoretical predictions. Further useful tools of diagnosing the exchange model are magnetothermal, e.g. magnetocaloric properties as well as high field magnetization and saturation field. This program has sofar mainly been carried out for frustrated J 1 -J 2 square lattice magnets [4][5][6].
Here we focus on the numerical investigation of much less studied finite temperature properties of anisotropic triangular quantum spin systems. Previous investigations used variants of Monte Carlo methods for the quantum case [7][8][9], also for the classical case [10,11] or analytical methods for the low temperature regime [12] and high temperature series expansion [13]. In the present work we use the finite temperature Lanczos method (FTLM) based on the exact diagonalization (ED) of finite lattice tiles for this purpose. The method is based on evaluating the partition function by averaging over random starting vectors in the ED procedure. Using periodic boundary conditions and discrete symmetries we can go up to largest cluster size of N = 28 sites for thermodynamic averages. We also consider smaller clusters of size N = 16, 20, 24 by FTLM. Furthermore we present the full analytical solution of the spectrum for the smallest N = 8 cluster as a reference point. However, we note that unlike zero temperature ED results for finite clusters the FTL method cannot be used for finite size scaling analysis of thermodynamic quantities. Nonetheless in the special case of the 1D chain system (J 1 = 0) we can compare to exact results which are in excellent agreement with the N = 28 cluster results. This indicates that our FTLM results are trustworthy to use as a representation for the thermodynamic limit.
Our main emphasis is the calculation of peak position T max and height C V (T max ) and χ(T max ) in the temperature dependence of specific heat and susceptibility, respectively. We investigate its systematic variation with anisotropy control parameter J 2 /J 1 or φ = tan −1 (J 2 /J 1 ), in particular when it is tuned between the simple special cases mentioned above. This information is of great importance for the analysis of thermodynamic data of triangular magnets. The T max (φ) dependence turns out to be considerably more asymmetric than for the related square lattice model [4]. We show that in certain ranges of the control parameter an indication of second maximum or shoulder in these quantities appears. Furthermore as an example we discuss the analysis of the full temperature arXiv:1506.09104v1 [cond-mat.str-el] 30 Jun 2015 1. Schematic view of anisotropic triangular exchange models. For J1 = 0 (φ = ±π/2) this reduces to 1D J2-spin chains ( ), for J2 = 0 (φ = 0) to the square lattice J1-HAF ( ) and for J1 = J2 (φ = π/4) to the isotropic triangular system ( ).
dependence of χ(T ) for the most prominent triangular quantum spin systems Cs 2 CuCl 4 and the isostructural Cs 2 CuBr 4 using FTLM results. We demonstrate that the derived exchange model is in excellent agreement with the one obtained from direct spectroscopic methods like inelastic neutron scattering (INS, Cs 2 CuCl 4 only) spin wave results and electron spin resonance (ESR) experiments, both in high fields. In Sec. II we define the model and its parametrization. The reference of the full solution for the eight-site cluster is discussed in Sec. III. The FTL method and its technical implementation are discussed to some extent in Sec. IV. Our main results on specific heat and susceptibility are presented in Secs. V and VI, respectively and the explicit comparison to Cs 2 CuCl 4 and Cs 2 CuBr 4 is discussed in Sec. VII. Finally Sec. VIII gives the conclusion and outlook.

II. ANISOTROPIC TRIANGULAR EXCHANGE MODEL AND ITS PARAMETRIZATION
The anisotropic J 1 − J 2 exchange model on the triangular lattice ( Fig. 1) is given by with where e x , e y are unit vectors along cartesian directions. Figure 1 may also be interpreted as tilted square lattice model with one diagonal J 2 bond cut out [2]. For the anisotropic triangular-lattice model, different parametrizations of the exchange energies are customary. As introduced in Ref. [4] and used subsequently we prefer the polar representation defined by where J c gives the overall energy scale and φ is the anisotropy control parameter. This parametrization allows for an easy interpolation between important geometrical limiting cases, namely the square-lattice Néel antiferromagnet with J 2 = 0 (φ = 0), the isotropic 120 • triangular antiferromagnet with J 2 = J 1 (φ = π/4), the antiferromagnetic chain with J 1 = 0 (φ = π/2), and their ferromagnetic counterparts. However, there are alternative possibilities in the literature which we briefly mention here for ease of comparison (Fig. 2). Many authors regard the model as an extension to the one-dimensional spin chain, therefore they use the exchange J 2 along these chains as the overall energy unit and parametrize their results in terms of α := J 1 /J 2 with J 1 being the interchain coupling. From a square-lattice perspective in turn, the exchange parameter J 1 appears to be the natural energy unit, and correspondingly 1/α := J 2 /J 1 is used as well. Both model parameters α and 1/α however are problematic when trying to describe the full phase diagram, in particular the interpolation between the square-lattice antiferromagnet (α → ∞) and the one-dimensional chain (1/α → ∞). To overcome this, the function f := J 2 /(J 1 +J 2 ) = 1/(1+α) has been introduced in Ref. [14], which remains finite in both limits and between these.
However the latter parametrization also does not cover the full phase diagram unambiguously, which is why we have introduced the universal energy scale J c and the anisotropy angle φ. In the square-lattice case, this has the additional advantage that J c , apart from a global factor √ 2, denotes the magnetocaloric energy scale J mc as well [4]. To facilitate the comparison of our results with the literature, Fig. 2 displays the dependence of the quantities introduced here on the anisotropy control parameter φ.
The nearest neighbor (n.n.) HAF model on the triangular lattice is the generic 'geometrically frustrated' spin system where the exchange energy of n.n. bonds cannot be minimized simultaneously for all bonds (see upper corner of Fig. 1). It is worthwhile to quantify this intuitive notion of 'frustration'. A measure for it is the total loss of exchange energy due to frustration relative to the exchange energy without it. For the basic triangular (three-site) plaquette in Fig. 1 the degree of frustration is then given by Here E is the ground state energy of the frustrated triangle. E t := E (J 2 = 0) and E d := E (J 1 = 0) are the ground state energies of its unfrustrated trimer and dimer parts, respectively with where the minimum is taken from the three different energy eigenvalues of the triangle. Using these expressions in Eq. (4) the degree of frustration (in per cent) in the triangular HAF as function of control parameter is shown in Fig. 3. There the dotted lines indicate the classical boundaries between FM, AF and spiral phases as discussed in Ref. [2]. κ(φ) vanishes identically in the unfrustrated (J 2 < 0 or φ < 0) case. In the frustrated regime (J 2 > 0 or φ > 0) it achieves its maximum values of κ = 4/7 corresponding to 57 % frustration at the isotropic point ( ) and the Spiral/FM phase boundary while it vanishes for the 1D chain case ( ) where J 1 = 0. We note that the reduction in the ordered ground state moment in the frustrated regime does not directly follow the degree of frustration but is the result of the subtle interplay of the latter with quantum fluctuations [2].

III. PRECURSOR: CLASSICAL PHASE DIAGRAM AND EIGHT-SITE FULL SOLUTION
Within the finite temperature Lanczos method the thermodynamic quantities will be computed directly us- ing an averaging procedure over the low energy spectrum only. To obtain a better intuition for the triangular Heisenberg model it is useful to calculate also the full spectrum by direct diagonalization of small clusters. We will demonstrate that it develops characteristic signatures at the classical phase boundaries and special symmetry points as function of control parameter using the eight-site cluster.
A. Spectrum and classical phase diagram As discussed in Ref. [4] for the square-lattice model, we can express the spectrum of the triangular eight-site cluster on tile 8:2-2 with periodic boundary conditions as a sum over Hamiltonians on complete graphs. We define where L denotes an ordered list of lattice sites. In our case it has either two elements (a bond), four elements (a square or tetrahedron), or eight elements (a cube). We then can write for the Hamiltonian. The tile and its labelling is illustrated in Fig. 4. From Eq. (7) with S = 1/2 being the maximum expectation value for each local spin S i , we get the corresponding list of eigenvalues by hierarchically constructing all possible multiple-spin configurations starting with the basic two-spin singlets and triplets on the pairs {1, 5}, {3, 4}, {2, 8}, and {6, 7}: From all pairs of two-spin states, we can construct the four-spin states with total spin S = 0, 1, 2, each possible pair of these four-spin states in turn can be combined to an eight-spin state with total spin S = 0, 1, 2, 3, 4 and spin degeneracy (2S + 1) . Because the total spin may be composed in several ways by the spins of sub-clusters there exist additional degeneracies. In this way, we can TABLE I. All energy levels of the eight-site cluster. S1−8 is the total cluster spin and the next six columns give the subcluster spins. The last column denotes the additional degeneracy of levels on top of the (2S1−8 + 1)-fold spin degeneracy. S1−8 S1345 S15 S34 S2678 S28 S67 energy deg.
construct a total of N N/2 = 70 states for the eight-site cluster. Table I displays the complete set of eigenvalues, together with their (additional) degeneracies, total spins, and total spins on the sub-lattices. Of particular interest are those states corresponding to classical ground states (marked by an asterisk in the last column): the columnar antiferromagnet (first state in the table, energy −6J 2 ), the Néel antiferromagnet (seventh state, energy 2J 2 − 6J 1 ), and the ferromagnet (last state, energy 2(2J 1 + J 2 ). These three states replace each other as ground states as a function of the anisotropy ratio φ, see  Table I. Solid lines denote singlet states, the long-dashed line denotes the ferromagnet, for further line types see the figure caption. We overlay this spectrum to the classical phase diagram of the model discussed in detail in Ref. [2]. The classical FM, AF and spiral phases are separated by the thin vertical lines. The thin dotted lines indicate as special cases the square-lattice Néel antiferromagnet ( , φ = 0 or J 1 > 0, J 2 = 0), isotropic antiferromagnetic triangular lattice ( , φ = π/4 or J 2 ≡ J 1 > 0) corresponding to the non-collinear 120 • commensurate spiral structure and one-dimensional antiferromagnetic chains ( , φ = π/2 or J 1 = 0, J 2 > 0).
Characteristic for the spectrum at these six special points is the high number of degenerate excited states with different total spin. At the borders of the classical ferromagnetic phase also the ground state changes from the fully polarized state to a singlet. However the classical phase boundary between the antiferromagnet and the spiral phase does not correspond to an equivalent change of (in this case nonmagnetic) ground states. This is due to the fact that with our eight-site tile we cannot approximate any incommensurate state at all, therefore the "best" approximation to the "true" ground state remains the Néel state for 1/2 ≤ J 2 /J 1 ≤ 3/4 (equivalent to 0.148 ≤ φ/π ≤ 0.205), and the columnar antiferromagnetic state for larger values of J 2 /J 1 .
Similar to the square-lattice case [4], we observe level crossings near the classical phase boundary between AF and spiral phase. For the square lattice, we took this as an indication for the appearance of a nonmagnetic phase, coinciding with the suppression of the ordered moment in a region around the classical transition observed within linear spin-wave theory [2]. Qualitatively the same happens here, however this analogy should not be taken too far, because no indication of that kind exists for the classically invisible crossover to one-dimensional chains at φ = π/2, which is in reality surrounded by a large nonmagnetic region. And empirically we know that linear spin-wave theory rather underestimates the size of these nonmagnetic regions. At the crossover from the spiral to the ferromagnetic phase for φ/π ≈ 0.852, a similar type of pattern of level crossings exists. Linear spin-wave theory gives inconsistent results in this case.

B. Thermodynamics of the eight-site cluster
Given the eigenvalue E I with degeneracy d I and total spin S I for each state I listed in Table I, we can easily evaluate the partition function Z = I d I (2S I + 1) e −βE I to calculate the magnetic susceptibility χ and the specific heat c V according to where β = 1/(k B T ) with k B being the Boltzmann constant. Here and in the following, we express the molar susceptibility in units of N L µ 0 (gµ B ) 2 J −1 c , where N L is the Loschmid constant, µ 0 the magnetic permeability constant, g the gyromagnetic ratio, and µ B the Bohr magneton. For the dimension of the specific heat c V , we use the universal gas constant R = N L k B .

IV. GENERAL REMARKS ON FINITE TEMPERATURE METHODS
The finite temperature properties of spin systems can be treated with analytical high temperature expansions [13,[15][16][17][18] or with the numerical FTL method [19] which will be employed in this work. As a reference we first discuss briefly the single first order term of the expansion method.

A. High temperature approximation
The high-temperature behavior of the susceptibility χ and the specific heat c V to quadratic order in β is determined by the Curie-Weiss energy Θ and the magne-tocaloric energy scale J mc , which are defined through The sums run over all bonds n connecting an arbitrary but fixed site i with its (not necessarily nearest) neighbors at sites {i + n}. We then have [18] In principle we should be able to determine the exchange parameters J 1 and J 2 already from high-temperature fits of the experimental results to the expressions above. However having fixed J 2 mc and Θ determines 2J 1 +J 2 and |J 1 − 4J 2 | but leaves the sign of the latter undetermined. When expressed with the parameters J c and φ, this is equivalent to the fact that there are, apart from special cases, always two possible values φ ± for the anisotropy parameter. These two values φ ± , however, can lie in two different thermodynamic phases with completely different properties. This ambiguity is similar to the one in the square lattice J 1 − J 2 model and its implications there were discussed in Refs. 4 and 20.
Although the coefficients of the high-temperature expansions for χ(T ) and c V (T ), being polynomial functions of J 1 and J 2 , are known up to at least eighth order [18], this situation essentially will not change by including further higher-order terms in a high-temperature expansion [20], and it remains difficult to determine J 1 and J 2 unambiguously solely from fits to the high temperature dependence of χ and c V . One powerful further diagnostic is the investigation of saturation fields [6,21] provided that they are in an accessible range.

B. Finite-temperature Lanczos method
To overcome this ambiguity, we use the finitetemperature Lanczos method [5,19,22] to evaluate the thermodynamic functions directly and compare specific heat and susceptibility temperature dependence over the whole temperature range above the finite size gap region. The method is based on the evaluation of thermodynamic traces using the eigenvalues and many-particle wave functions determined by numerical exact diagonalization of the Hamiltonian matrix on finite tiles. After mapping the Hamiltonian onto a sparse matrix, we use the iterative Lanczos algorithm [23] to generate the first few (between 1 and 100) extremal eigenvalues and the corresponding wave functions. Due to the Boltzmann weight, these are the most important eigenvalues contributing to the partition function. We classify our wave functions according to the expectation value of the z component Ω z = N i=1 S z i of the total spin Ω, the crystal momenta k and the point group symmetries of the tile. This brings the Hamiltonian matrix into block diagonal form, allowing us to go to tile sizes up to N = 28 in our finite-temperature calculations. In this way we can evaluate thermodynamic traces on industry standard computer hardware.
The tiles we are using for the diagonalization procedures are illustrated in Fig. 6. We apply the same labelling as described in Ref. [24] for the square-lattice case, adapted to the triangular lattice. There are two important differences: At finite temperatures we cannot easily perform any kind of finite-size scaling analysis of the partition function, therefore we directly take the results of the different lattice tilings. In addition there exists at least one phase corresponding to a classical phase with an incommensurate ordering vector. Due to the finiteness of our k space grid, we cannot model this closely. However for thermodynamic properties like susceptibility and heat capacity all thermally populated states contribute, and the exact modeling of the ground state as a function of our parameters is of secondary importance. For zero temperature ED approach the introduction of twisted boundary conditions may provide a way to circumvent this problem [25].
We attempt to determine the thermal expectation value of an arbitrary static operator A by the fundamental traces over the statistical operator, To sample an as large as possible part of the Hilbert space, we start N R = O(100) iterations with different random wave functions or starting vectors |r , such that eventually we use no more than O 10 4 eigenvalues and wave functions per Hamiltonian block. It may be shown [19] that this procedure yields an asymptotically exact result. In summary, the thermal expectation value of A is approximated by where the index s denotes the summation over all symmetry sectors of the Hilbert space with dimension N s st . If the operator A is a conserved quantity with [H, A] = 0, we can replace A by its quantum numbers A r,s j , and Eq. (16) further simplifies to The general definitions of the volume magnetic susceptibility and the heat capacity are where H = B a /µ 0 is the external magnetic field defining the z direction and F = (−1/β) ln Z is the canonical free energy. Written in terms of the expectation values defined above, we get the cumulants for the molar susceptibility and the specific heat of a tile with N sites, respectively. In our case, both Ω z and trivially H commute with H, allowing us to use Eq. (18) to determine the thermal expectation values. Furthermore spontaneous magnetic order is absent because our tiles are finite and we do not include an external magnetic field, therefore we can safely set Ω z β = 0. As mentioned before, in the comparison of calculated values and experimental results two strategies are possible. One can identify the maximum position and value of thermodynamic quantities or perform a fit over the whole temperature range to extract the exchange model parameters. We will discuss both approaches.

V. HEAT CAPACITY
First we discuss the heat capacity which requires only the cluster eigenvalues for its evaluation. A comparison of the theoretical temperature dependence to experiments is, however, not straightforward due to the lattice contribution to the heat capacity [3].

A. Temperature dependence
To demonstrate characteristic features, Figs. 7 and 8 show the temperature dependence of the specific heat c V (T ) according to Eq. (21) for a selected range of anisotropy parameters φ in the antiferromagnetic (Fig. 7) and the spiral phase (Fig. 8). The FTLM calculations are trustworthy only down to a temperature range of T ≈ J c /(N k B ). For lower temperatures they are dominated by the artificial finite size gaps. This excluded region is indicated by grey bars in the figures.
For the AF phase with values −0.4 ≤ φ/π ≤ −0.21 a single peak indicated by dots that shifts continuously to higher values with increasing φ is observed. This is qualitatively clear already from the eight-site spectrum (Fig.5) which shows an increasing average excitation gap from the ground state in that range of φ. Furthermore a plateau (or very flat second maximum) at the lowest temperatures close to the finite size gap region is visible.
For the spiral phase the values 0.25 ≤ φ/π ≤ 0.44, correspond to the region between the isotropic triangular lattice with J 2 /J 1 = 1 and deep inside the disordered phase with quasi-one-dimensional chains (J 2 /J 1 ≈ 5). The tile 28:2-4 was used for the numerical calculations in both cases.
In Fig. 8, the characteristic maxima of c V (T ), indicated by the small black dots in the figure, have positions which fall into two different temperature ranges: For the isotropic triangular case with φ/π = 0.25, the maximum sits at T ≈ 0.18J c /k B , reducing in magnitude and shifting slightly to lower temperatures T ≈ 0.15J c /k B upon increasing φ from its isotropic value to φ/π ≈ 0.31 or J 2 /J 1 ≈ 1.5. At this point, a second maximum develops starting at T ≈ 0.32J c /k B , shifting to higher temperatures as φ is increased towards the disordered region. Furthermore, the maximum T characteristic for the triangular lattice rapidly disappears with increasing φ. Only the C V (T ) curves for φ/π = 0.31 and φ/π = 0.32 show both maxima simultaneously. At high temperatures T J c /k B , all the curves show the (1/T ) 2 temperature dependence expected from Eq. (13).  The eight-site cluster, introduced as an illustration of the model and its overall spectrum, clearly is not very useful quantitatively to discuss the heat capacity in the thermodynamic limit. Only the overall behavior of the maximum temperature is qualitatively similar to our findings for the larger tilings, however the two minima in the spiral phase are shifted towards the classical phase borders. The values c V (T max ) for the eight-site cluster are monotonically increasing in the antiferromagnetic phase, followed by a minimum at φ/π ≈ 0.1 and a maximum at the isotropic point, φ/π = 0.25. On the ferromagnetic-J 1 side for 1/2 < φ/π ≤ 1, c V (T max ) in- creases again to a second maximum, followed by a minimum at the crossover to the ferromagnet at J 2 /|J 1 | = 1/2. At the one-dimensional point J 1 = 0, naturally both T max and c V (T max ) differ strongly from the exact values.

����-� ��
Also the 16-site tiling shows similar features. For the larger clusters of size N = 20 and more, both maximum temperature positions and maximum values of the specific heat clearly show a double-peak structure as function of φ. The positions T max decrease with increasing cluster size, apart from the regions around the isotropic point and in the spiral phase near the ferromagnet. At the one-dimensional point, agreement with the infinite-chain result [26] (white circles) is achieved already with N = 20 tiling. And therefore results for the N = 28 tile represent well the thermodynamic limit as far as peak position and height is concerned. De-noting the largest value in each sector of Fig. 9 by T * max we observe that T * max (AF)/T * max (spiral) = 1.58. This asymmetry in T max (φ) is considerably larger than the corresponding one in the square lattice Inside the AF and spiral sectors the lowest T max is reached close to the isotropic triangular point with T max ≈ 0.18J c /k B because there exchange frustration is most pronounced (Fig. 3). This leads to a high density of low energy excitations (c.f. Fig. 5) and therefore a low T max .
Furthermore, in parts of the antiferromagnetic as well as the spiral phase, a second maximum c V (T max2 ) at very low temperatures appears. T max2 depends roughly linearly on φ while c V (T max2 ) remains constant and small, see also Fig. 7. With the 28-site tiling being the largest possible, we cannot judge whether this second-low temperature maximum in the antiferromagnetic phase merely is a finite-size effect. This is different from the situation in the spiral phase, where the partially observed second maximum is at temperatures T max2 far larger than the finite-size gap, see also Fig. 8. Similar to the behavior around the isotropic point, at the crossover to the ferromagnetic region a second maximum gradually appears in c V (T ) which then evolves to the "ferromagnetic" maximum while the "spiral" maximum disappears. This also happens at comparatively low temperatures with an irregular behavior of the maximum values as function of φ. We attribute this irregularity to the fact that in particular near the borders of the spiral phase the ground state and possible low- lying excited states have incommensurate ordering vectors Q = Q AF + δQ and Q FM + δQ respectively with |δQ| 1. These are not contained in the coarse grid of crystal momenta of our finite tiles.

VI. SUSCEPTIBILITY
Evaluation of the susceptibility with FTLM is more involved because it requires the calculation of magnetic moment matrix elements. On the other hand this quantity can be more easily compared to experimental results. Most known exchange parameters for spin systems are due to the application of this method.

A. Temperature dependence
Similar to the specific heat, we have calculated the temperature dependence of the static magnetic susceptibility χ(T ) according to Eq. (20) with the tilings shown in Fig. 6. Fig. 11 shows the results in a selected range of the anisotropy parameter, 0.57 ≤ φ/π ≤ 0.86, calculated with tiling 28:2-4. This parameter range corresponds to the interpolation between ferromagnetically coupled (J 1 < 0) quasi-one-dimensional antiferromagnetic chains with J 2 /J 1 ≈ −4.5 and the crossover to the ferromagnet at the classical boundary J 2 /J 1 = −1/2. For the former, a very broad maximum is characteristic which gradually evolves into the T = 0 divergence of χ(T ) in the ferromagnetic phase, which is the expected behavior.

B. Exchange anisotropy dependence
In the same way as for the specific heat, we follow the positions and values of the characteristic maxima of χ(T ) as a function of the anisotropy parameter, using the tiles displayed in Fig. 6 again. This is shown in Fig. 10. The solid lines denote the eight-site results according to Eq. (8). For the same reasons as discussed for the heat capacity, strong deviations from the larger tilings occur in particular in the spiral phase. In contrast to c V (T ), the maximum positions and values of χ(T ) are already converged for tilings of size N = 20 and larger, apart from the crossover to the ferromagnetic phase, where a tiledependence of the maximum temperatures clearly can be observed.
Similar as for the specific heat the T max (φ) dependence has a pronounced asymmetry in AF and spiral sectors. Denoting again the largest value in each sector by T * max we obtain for the triangular lattice T * max (AF)/T * max (spiral) = 1.58 which is much larger than the asymmetry value T * max (AF)/T * max (CAF) = 1.11 for the square lattice case. Similar to specific heat behavior the T max (φ) minimum inside AF and spiral sectors is reached around the most frustrated isotropic triangular position T max ≈ 0.35J c /k B .
At the one-dimensional point, as for the heat capacity, the results from Ref. [26], k B T max /J c = 0.641 and χ(T max )/(N L µ 0 (gµ B ) 2 /J c ) = 0.147 are accurately reproduced with our method. In general we think that in particular our results on the magnetic susceptibility can be used to accurately determine both exchange constants J 1 and J 2 individually, which we show in the following section.

VII. APPLICATION TO Cs2CuCl4 AND Cs2CuBr4
As a demonstration of the usefulness of the method, we apply our findings to the compounds Cs 2 CuCl 4 and Cs 2 CuBr 4 . Fig. 12 shows the temperature dependence of the magnetic susceptibility. The dots and the open circles denote the experimental data taken from Ref. [27], the solid and dashed lines denote the corresponding fits with our FTLM data. In order to avoid misunderstandings, we have plotted the data and the curves in two different ways: The top plot in Fig. 12 displays χ(T ) in electromagnetic (CGS) units from where we have determined the values for our model parameters reproduced in Table II. Using these fitted values for J c , φ, and g, we have plotted the same data again in dimensionless units in the bottom plot of Fig. 12. The main effect of replacing Cl with Br is an increase of the overall energy scale J c by a factor 3.4, whereas the anisotropy angle φ changes only by about 5 %. Thus only this large change in J c is responsible for the decrease of the maxi-   Table II). The bottom plot displays exactly the same data, this time in dimensionless units using Jc, φ, and g for the two compounds from Table II. mum in χ(T ), the broadening of the maximum, and the shift of the maximum towards higher temperatures. The anisotropy ratio and therefore the position of Cs 2 CuBr 4 in the phase diagram remains essentially the same as for Cs 2 CuCl 4 , it is only slightly moved away from the quasione-dimensional regime. Historically the first measurement of χ(T ) of Cs 2 CuCl 4 yielded slightly different values [30], however the data analysis was performed over a limited temperature range assuming weakly coupled onedimensional chains. As pointed out in Sec. V A, the steep decrease of χ(T ) for T → 0 in both cases is an artifact of the finiteness of the tiling used for the FTLM calculations. Assuming a finite-size gap ∆ ≈ J c /N , the corresponding temperatures where we expect finite-size effects to dominate are T ∆ ≈ 0.2 K for Cs 2 CuCl 4 and T ∆ ≈ 0.6 K for Cs 2 CuBr 4 . Experimentally, the lowest temperature was T min ≈ 2 K, well above the finite-size gaps.
Our result from the FTLM fit to thermodynamic data are compared in Table II to results from direct spectroscopic methods: Inelastic neutron scattering (INS) [28] (Cs 2 CuCl 4 only) and electron spin resonance (ESR) [29], both in fields above the saturation field [2] which is µ 0 H sat ≈ 8.4 T for Cs 2 CuCl 4 and µ 0 H sat ≈ 30 T for Cs 2 CuBr 4 . Taking the spectroscopic data at high fields has the advantage that the ground state is fully polarized, corresponding to a single "all spins up" spinor product state |FM . The excitations on top of this are the N orthonormal single-particle excitations |ψ i = (1/ √ 2S)S − i |FM , i = 1 . . . N with N = O(N L ). Because [H, Ω z ] = 0, the Hamiltonian can not generate more than these one-spin-flip states, and its spectrum can be determined exactly by Fourier transform.
The agreement between the different methods is almost perfect, for exchange parameters as well as g-factors. As already mentioned in Ref. [2] if Cs 2 CuCl 4 is interpreted as a purely 2D system this set of exchange parameters puts the compound very close to the quasi-1D spin liquid regime centered at φ = 0.5π. In reality, however, the magnetic order is stabilized below T N = 0.62 K by a finite inter-plane coupling of the order J ⊥ ≈ 0.017 meV along the crystallographic a direction [31].
For Cs 2 CuBr 4 , we see deviations of the experimental data from our result at the lowest temperatures. These are not due to impurities [27,32], but can be regarded as an indication for a tendency towards magnetic order in this compound. The overall energy scale J c is more than three times larger than for Cs 2 CuCl 4 , however the estimate for the anisotropy angle φ is essentially the same. We therefore expect that Cs 2 CuBr 4 orders magnetically as well, possibly at a higher temperature than its Cl counterpart.
The agreement of our thermodynamic analysis with direct spectroscopic results gives us confidence that the FTL method may be applied to the analysis of the whole Cs 2 CuCl 4−x Br x substitutional series [27,32].

VIII. CONCLUSION AND OUTLOOK
In this work we have shown that finite temperature Lanczos method for finite clusters is a versatile tool to investigate the little known finite temperature properties of frustrated triangular quantum magnets, in contrast to QMC approach which is not suitable in this case. The FTL method complements the analytical spin wave approach which is useful only for the very low temperature regime and is a more straightforward alternative to the high temperate series expansion method.
For this quantum spin model one can use a single control parameter and tune the system through the phase diagram where several special cases with very different frustration degree like frustrated Néel, unfrustrated HAF, frustrated spiral and isotropic triangular 120 • phase as well as unfrustrated spin chains and frustrated FM may be realized.
For the largest investigated cluster with N = 28 sites we obtain a trustworthy representation of the thermodynamic limit behavior since the exact known result for the spin chain case is produced very well and shows little difference to the N = 24 cluster. As a main result we gave the systematic variation of peak position and peak height of specific heat and susceptibility as function of J 2 /J 1 . Both values are strongly suppressed for the most frustrated isotropic triangular magnet (Figs. 9 and 10).
Furthermore a surprisingly large asymmetry in particular for the maximum position exists between the AF and spiral phase region. It is considerably larger than between the AF and CAF regions in the square lattice J 1 -J 2 model. In addition we found that the simple single peak shape of c V (T ) and χ(T ) may be distorted due to smaller side maxima or shoulders at lower temperature for some ranges of the control parameter.
The most reliable method to extract the exchange parameters is the fitting of χ(T ) over the whole temperature range (above the finite size gap) using FTLM results. We have demonstrated this for two of the most typical 2D triangular magnets, Cs 2 CuCl 4 and Cs 2 CuBr 4 . We have obtained excellent agreement with results from the INS investigation and with ESR results in the fully polarized state. Our method has the additional advantage that it can easily be extended to the whole substitutional series Cs 2 CuCl 4−x Br x (0 ≤ x ≤ 4) [27] to extract the systematic variation of frustration anisotropy control parameter φ and energy scale J c as a function of Br concentration [33]. We note that our method can also be applied to magnetocaloric measurements of the organic charge-transfer salts where localized S = 1/2 magnetic moments are residing on a possibly distorted triangular lattice as well. This includes the Et x Me 4−x Z[Pd(dmit) 2 ] 2 (x = 0, 1, 2) family of compounds [34,35] as well as κ-(ET) 2 B(CN) 4 and κ-(ET) 2 Cu 2 (CN) 3 [36] and κ-(BEDT-TTF) 2 Cu[N(CN) 2 ]Cl [37].