{'fraction_non_alphanumeric': 0.03691440760255329, 'fraction_numerical': 0.011031343022237787, 'mean_word_length': 4.206530747555114, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnetocaloric effect in the frustrated square lattice J1-J2 model

We investigate the magnetocaloric properties of the two-dimensional frustrated J1-J2 model on a square lattice. This model describes well the magnetic behavior of two classes of quasi-two-dimensional S = 1/2 vanadates, namely the Li2VOXO4 (X = Si, Ge) and AA'VO(PO4)2 (A, A' = Pb, Zn, Sr, Ba) compounds. The magnetocaloric effect (MCE) consists in the adiabatic temperature change upon changing the external magnetic field. In frustrated systems, the MCE can be enhanced close to the saturation field because of massive degeneracies among low lying excitations. We discuss results for the MCE in the two distinct antiferromagnetic regimes of the phase diagram. Numerical finite temperature Lanczos as well as analytical methods based on the spin wave expansion are employed and results are compared. We give explicit values for the saturation fields of the vanadium compounds. We predict that at subcritical fields there is first a (positive) maximum followed by sign change of the MCE, characteristic of all magnetically ordered phases.


I. INTRODUCTION
Two dimensional (2D) magnetic systems are favorite models to study the influence of quantum fluctuations on magnetic order. Depending on the model they may both prohibit an ordered ground state or select a specific order among classically degenerate states. These phenomena have been studied in great detail for geometrically frustrated systems like trigonal, Kagomé and checkerboard lattice 1 . However they are also present in magnets where the frustration is not the result of lattice geometry but of competition between different (for example nearest-and next-nearest neighbor) magnetic bonds. A prime example is the frustrated J 1 -J 2 model on a square lattice. Its ground state and thermodynamic properties in zero field have been well studied (see Refs. 1,2,3 and references cited therein).
Classically the model predicts three magnetic phases depending on the frustration ratio J 2 /J 1 : The ferromagnet (FM), (π, π) Néel antiferromagnet (NAF) and (π, 0) collinear antiferromagnet (CAF). However it is known that close to the classical CAF/NAF and CAF/FM boundary quantum fluctuations destroy magnetic order and presumably stabilize nonmagnetic order parameters.
The discovery of two classes of layered vanadium oxides Li 2 VOXO 4 (X = Si, Ge) 4,5,6 and AA ′ VO(PO 4 ) 2 (A, A ′ = Pb, Zn, Sr, Ba) 7,8 which are well described by this model has further raised interest in the J 1 -J 2 model. One advantage of the new vanadium compounds is a comparatively low energy scale for the exchange constants of order 10 K. Therefore high field experiments might be a promising way to learn more about their physical properties, indeed the saturation field for these compounds where the fully polarized state is achieved seems within experimental reach.
Therefore in this work we study exhaustively the highfield magnetic and especially the magnetocaloric effects (MCE) in the J 1 -J 2 model. We use a variety of analytical and numerical techniques to investigate the dependence of magnetization, susceptibility, entropy specific heat and adiabatic cooling rate on magnetic field, temperature and frustration ratio. The variation of the saturation field with the frustration ratio is calculated and predictions for the abovementioned compounds are made. We show that the low temperature specific heat is strongly enhanced around the classical phase boundaries where large quantum fluctuations occur.
Our special focus is on the magnetocaloric effect. We will show that the cooling rate, normalized to its paramagnetic value, is strongly enhanced above the saturation field and depends on the frustration angle. We also predict that for subcritical fields the cooling rate is first positive with a maximum at moderate fields and a subsequently changes sign at a larger subcritical field. This behavior is common to all AF phases of the model and can be understood quantitatively from calculations of contours of constant entropy in the (h, T ) plane. The dependence of corresponding characteristic fields on the frustration ratio are also calculated.
In Sec. II we give a brief description of the magnetocaloric effect in magnets. In Sec. III the basic properties and phase diagram of the J 1 -J 2 model are introduced. In Sec. IV we discuss extensively results of the finite temperature Lanczos method (FTLM) for finite 2D J 1 -J 2 clusters. In Sec. V we use analytical methods within mean field or spin wave approximation as an alternative way to study the magnetocaloric properties. In Sec.VI we discuss and compare the results obtained by the various methods. Finally Sec. VII gives the summary and conclusion.

II. THE MAGNETOCALORIC EFFECT IN MAGNETICALLY ORDERED COMPOUNDS
When a crystal containing magnetic ions is placed in a magnetic field the adiabatic or isentropic change of this external parameter causes a temperature change in the sample. This is called the magnetocaloric effect (MCE) which was first discovered by Warburg 9 . It is nowadays interesting in several aspects. Firstly suitable compounds, like paramagnetic salts where demagnetization leads to cooling may be used technically 10,11 . Secondly at high (pulsed) fields the magnetocaloric anomalies at a magnetic phase transition may be used to map out the H-T phase diagrams which are not accessible otherwise. Finally it has recently gained special attention in frustrated magnets. There the behavior around the saturation field may be described by the condensation of a macroscopic number of local magnons 12,13,14 which leads to a giant enhancement of the magnetocaloric cooling rate. The latter is defined as the rate of change of temperature T with magnetic field H at fixed entropy S. Using a Maxwell relation one can write this as where C V is the heat capacity and M the magnetization of the sample. The integrated adiabatic temperature change along an isentropic line with S(T, H) = const which is caused by the variation of magnetic field is then given by Here T 0 and H 0 are the starting values of temperature and field respectively. We take the adiabatic cooling rate for free paramagnetic ions as a reference quantity. As shown later it is simply given by Γ 0 mc = (T /H). The dimensionless magnetocaloric enhancement factor due to interaction effects is then defined by the ratiô In the present work we study the MCE on a Heisenberg square lattice J 1 -J 2 model which incorporates a frustration of nearest-and next-nearest neighbor exchange interactions. It is a suitable model to analyze the magnetic properties of two classes of quasi-two-dimensional vanadates, namely the Li 2 VOXO 4 (X = Si, Ge) 4,5,6 and AA ′ VO(PO 4 ) 2 (A, A ′ = Pb, Zn, Sr, Ba) compounds 7,8 . We investigate several aspects of magnetocaloric properties. We show that indeed it may be used to identify the saturation field and the associated H c (T ) phase boundary between the fully polarized and the AF ordered states. We also discuss whether it may be used as a diagnostic for the appropriate frustration angle (or J 2 /J 1 ratio). Finally we will study whether the frustration effect lead to a visible signature in the anomalies of the magnetocaloric cooling rates, especially close to the saturation fields. We will employ both numerical FTLM methods as well as approximate analytical methods based on mean field and spin wave approximations.

III. THE J1-J2 HEISENBERG MODEL AND EXAMPLES
We first give a brief characterization of the 2D square lattice J 1 -J 2 model in an external field which is defined by the Hamiltonian with the convention that each bond is counted only once. Here J 1 is the nearest-neighbor exchange coupling along the edges and J 2 the next-nearest neighbor exchange coupling along the diagonals of each square. Furthermore h = gµ B H where H is the applied magnetic field. Here g is the gyromagnetic ratio and µ B the Bohr magneton. The zero-field phase diagram may best be characterized by introducing the equivalent parameter set The 'frustration angle' φ is a convenient quantity to characterize the amount of exchange frustration in the model and J c gives the energy scale at which thermodynamic anomalies in specific heat susceptibility etc. are to be expected. For spin-1/2, as function of φ three main phases (FM,NAF,CAF) appear already on the classical level ( Fig. 1) 2 . However for φ ∼ 0.15π (J 2 /J 1 ∼ 0.5) and φ ∼ 0.85π (J 2 /J 1 ∼ −0.5) where CAF meets NAF and FM respectively a large degeneracy of the classical ground state appears and quantum fluctuations lead to non-magnetic "hidden order" phases shown as the shaded sectors in Fig. 1. These phases have been extensively discussed in 1,2,3 and references cited therein.
For the compounds mentioned above generally J c ≃ 10 K. In addition from the high temperature expansion of the susceptibility the Curie-Weiss temperature is obtained as Θ CW = (J 1 + J 2 )/k B . Both quantities can be obtained from experiment and the appropriate pair of exchange constants (J 1 , J 2 ) then has to lie on the intersect of a circle (J c ) and a straight line (Θ CW ) as shown in Fig. 1. Obviously there are always two solutions lying in the NAF (φ − ) and CAF (φ + ) sector. This observation is unchanged by a more detailed analysis of susceptibility and specific heat 2,15 .
Various other methods have been proposed to resolve the ambiguity of frustration angles such as measurement of the spin structure factor 2 , nonlinear susceptibility 16 and the saturation field of the magnetization 17 . Only the former has been tried sofar for Li 2 VO(SiO 4 ) and Pb 2 VO(PO 4 ) 2 18 . In both cases the ground state clearly has CAF order. For this reason we have assigned the other known family members of J 1 -J 2 vanadates to the same sector in Fig. 1 although a confirmation for this conjecture is still lacking. We believe that high field investigations are a further promising method to shed light on these compounds , especially because the comparatively low J c (∼ 10 K) will lead to saturation fields relatively easy to access.
Therefore in this work we study the J 1 -J 2 model in an external field as given by Eq. (4). Thereby we focus on the theory of the saturation field and the magnetocaloric anomalies both around the saturation field and for smaller fields within the ordered phase. We will use both numerical analysis of finite clusters based on the FTLM method as well as analytical methods based on mean field or linear spin wave approximations for comparison. In the latter we focus on the three main phases with magnetic order. The analytical treatment of the magnetocaloric effect in the hidden order phases warrants a separate treatment which takes into account the proper non-magnetic order parameter.

IV. EXACT DIAGONALIZATION FOR FINITE CLUSTERS AT FINITE TEMPERATURES AND FINITE FIELD
We have performed numerical exact-diagonalization calculations for three different clusters; squares with 16 and 20 sites and a 24-site rectangle. All of these tile the lattice in such a way as to be compatible with both the (π, π) NAF and(π, 0) CAF states, once periodic boundary conditions are imposed. Our main focus is on the finite-temperature, finite-field properties of the J 1 -J 2 model. We therefore use the finite-temperature Lanczos method (FTLM) to evaluate the partition function of the model, together with thermodynamic averages of the form Here A is an operator, H is the Hamiltonian of the J 1 -J 2 model including the Zeeman term (Eq. (4)) and Z(T, H) its (field-dependent) partition function. For each frustration angle φ, we perform between 100 and 500 Lanczos iterations with different starting vectors in each symmetry sector of the Hilbert space. We keep between 1 and 100 eigenvalues and eigenvectors of the tridiagonal Lanczos matrix per iteration in order to evaluate the thermodynamic traces discussed in the following. Details of the method can be found in Ref. 19.

A. Level crossings, spinwave instabilities and saturation fields
Before addressing the finite-temperature results, let us examine certain general features of the model at zero temperature. Applying a magnetic field H leads to a Zeeman splitting of the energy levels, and therefore level crossings occur when increasing the field. These level crossings correspond to jumps in the magnetization at zero temperature, until the fully polarized state is reached at a certain critical value of the magnetic field. Figure 2 illustrates this behavior for two different values of the frustration angle φ: As an example for positive (antiferromagnetic) J 1 , we show the field dependence of the energy levels for the pure Néel antiferromagnet (J 2 = 0) on the left side of the figure. In the whole "right half" of the phase diagram (J 1 > 0, J 2 arbitrary), the field dependence of the energy levels is qualitatively similar. For the 24-site cluster considered, 12 spin flips with ∆S z = 1, indicated by the small arrows, occur at the points where the magnetic field is given by Here, E 0 (S z ) denotes the ground-state energy for the subspace with constant S z at zero field where N is the cluster size. The saturation field H sat is reached when S z = N/2. Because the fully polarized state is an eigenstate of the Hamiltonian, the numerical values for H sat from the equation above are exactly identical to what one finds within linear-spinwave theory for the infinite system, with z = 4, S = 1/2, and Q = (π, π) or one of (π, 0), (0, π) is the antiferromagnetic ordering vector. In contrast, for ferromagnetic J 1 < 0 (but φ outside the ferromagnetic regime in the phase diagram), we see a qualitatively different behavior of the energy levels. The right-hand-side of Fig. 2 shows an example for the frustration angle φ/π = 0.84. Instead of 12 level crossings, there are only six, each corresponding to a jump ∆S z = 2 occurring at fields The saturation field is now given by an instability criterion of the fully polarized state towards a two-magnon excitation, (11) For J 1 < 0, this field is larger than the field of the onemagnon instability given by the above equation with k = 1 and therefore determines the predominant instability when lowering the field in the fully polarized state.
A necessary condition for a ∆S z = 1 level crossing to occur is that the lower bound E 0 (S z ) of the energy spectrum at zero field for a fixed value S z is a convex function of S z , i. e., the condition must be fulfilled at H = 0. At the special point J 1 = 0, J 2 > 0 (φ = π/4), the J 1 -J 2 lattice decouples into two independent Néel sublattices, and equality holds above. For the finite size clusters which we consider, enlarging φ further (i. e. making J 1 ferromagnetic) stabilizes two-magnon bound states, and Eq. (12) no longer holds. Level crossings are characterized by ∆S z = 2, and the saturation field H c2 is given by Eq. (11) with k = 2. This is exactly what would be expected where a spin nematic state is selected by quantum fluctuations in applied magnetic field 3 .
However these are finite size results, and must be approached with a little caution. The critical fields associated with one-and two-magnon excitations show quite different finite size scaling properties as a function of φ, as illustrated in Fig. 3. From the three different cluster sizes studied here, the following observations can be made : Firstly, ∆H c is a non-monotonic function of the frustration angle; it has a minimum ∆H c = 0 at the crossover between the NAF and CAF phases for J 1 = 2J 2 (φ/π ≈ 0.15). Secondly, ∆H c changes sign at J 1 = 0 in the CAF phase (φ/π = 1/2) in favor of a ∆S z = 2 instability as above. Thirdly, for −1/2 ≤ φ/π 0.8, |∆H c | is a monotonically decreasing function of 1/N and seems to extrapolate to zero for N → ∞. This means that the one-magnon instability (and conventional canted AF order) is restored in the thermodynamic limit for FM J 1 and all J 2 0.6|J 1 |. Only close to the classical CAF/FM boundary at J 2 = 0.5|J 1 | does a two-magnon instability (with associated nematic order) prevail. These results are in complete agreement with previous exact analytic calculations for two-magnon bound states in the thermodynamic limit, and numerical exact diagonalizations of larger clusters 3 .
Using the two values φ = φ ± for the frustration angle together with the experimental energy scale J c determined from zero-field susceptibility and heat capacity measurements 2,7,8,20 , we can extract the expected values for the two different saturation fields H sat = H ± from the right-hand side plot of Fig. 3. Assuming a value g = 2 for the average gyromagnetic ratio, we arrive at field values between 13 and 24 T, low enough to be reached experimentally. In Fig. 4 we have plotted the predicted values for H sat as a function of the frustration angles φ ± for the known compounds ( Fig. 1), using the corresponding values of J c from zero-field measurements. H sat can be determined for example by a magnetization measurement at sufficiently large fields. For the PO 4 -based compounds, the saturation field together with the zero-field data for the susceptibility and the heat capacity would provide a direct way to determine the exchange constants J 1 and J 2 individually and hence the region of the phase diagram to which the compound belongs, without the need to measure the magnetic ordering vector directly.

B. Magnetization and susceptibility
We have calculated the magnetization m(T, H) and the magnetic susceptibility χ(T, H) = N A µ 0 (∂m(T, H)/∂H) by evaluating the following thermodynamic traces: where we have explicitly included a factor 1/N to account for the volume dependence of these extensive quantities.
In the definition of the susceptibility, we also include the magnetic permeability µ 0 and the Avogadro number N A . In order to make χ(T, H) a dimensionless quantity, we need to multiply it with the characteristic energy scale J c . At T = 0, the magnetization m(T = 0, H) of any finite size system evolves as a series of discrete steps. For a generic AF system with a singlet ground state, m(T = 0, H) takes on all possible (integer) spin values as a function of H, up to the field H c at which the system saturates. Generally, in the thermodynamic limit, m(T = 0, H < H c ) is a smooth curve, and singular features occur only where there is a magnetic phase transition. However in the case of the spin-1/2 J 1 -J 2 model, a "step" at exactly half the saturation magnetization m(T = 0, H) = 1/2 survives in the thermodynamic limit for J 2 ≈ J 1 /2, i. e. where a nonmagnetic ground state separates NAF and CAF order. This half-magnetization "plateau" is believed to be associated with the formation of localized magnon excitations 21 .
Temperature acts to smear jumps in magnetization. For the small clusters which we consider, the steplike behavior in m(T, H) has already disappeared for k B T = 0.2 J c . At the same time, all trace os the halfmagnetization plateau is also lost. The magnetic susceptibility therefore shows a smooth and nearly constant field dependence, see Fig. 5. It drops to zero upon reaching the saturation field. Only at the borders of the ferromagnetic regime for h = 0 do anomalies related to spontaneous magnetization appear.

C. Entropy and heat capacity
The entropy and heat capacity are defined through again using the definition for the thermal averages given in Eqs. (6,7). Z and H are the partition function and the Hamiltonian of Eq. (4), respectively. The right-hand side of Fig. 5 shows a contour plot of the entropy S(T, H) at fixed temperature k B T = 0.2 J c as a function of the frustration angle φ and the magnetic field H. In the ordered phases, the entropy is a smooth and almost constant function of the magnetic field, dropping to 0 for fields higher than the saturation field H sat . Characteristic anomalies can be observed at both edges of the collinear phase, where the entropy crosses a broad maximum as a function of field before again vanishing when crossing the saturation field. Figure 6 shows contour plots of the entropy (left) and heat capacity (right) as a function of magnetic field and temperature for a fixed value of φ = 0.747 π. We have chosen this particular frustration angle for the plots because it is believed to belong to φ = φ + for the compound SrZnVO(PO 4 ) 2 . The wiggly contour lines at low temperatures k B T ≪ J c are finite-size effects, where each temperature minimum corresponds to a Zeeman level crossing at T = 0 as discussed in Section IV A.
The lines of constant entropy (Fig. 6, left) are almost field-independent or even have a slightly negative slope as a function of field for low temperatures T ≪ J c /k B and fields H ≪ H sat . This implies that a sample cools down slightly when increasing H. The behavior of a paramagnet is opposite: Here, isentropic lines are straight lines crossing the origin, and a sample always heats up when increasing the applied field. Of course, for high enough fields and temperatures, the behavior of the entropy of the J 1 -J 2 model is the same as that of a paramagnet. At the saturation field, which is H sat ≈ 1.64 J c /(gµ B ) for φ/π = 0.747, the temperature reaches a minimum when adiabatically changing the field at low temperatures, and rises steeply when increasing the field to higher values H > H sat . In this area of the phase diagram, a J 1 -J 2 compound is a good system for magnetic cooling, especially in view of the low values for the saturation fields (in Tesla) for the experimentally known compounds (see Fig. 4 and its discussion above).
The heat capacity C V (T, H) is characterized by two maxima as shown in the right panel of Fig. 6. One maximum occurs at H = 0 and k B T ≈ 0.5 J c , which is the broad anomaly occurring at the crossover to the 1/T 2temperature dependence for high temperatures: We have This maximum has already been discussed in Ref. 2. A second maximum can be observed at high temperatures T ≫ J c /k B and magnetic fields H ≫ H sat . Here, a fieldinduced gap opens, leading to a Schottky-type anomaly of an effective two-level system, see also Sec. V A and Fig. 12 (center).

D. The magnetocaloric effect
For the numerical calculation of the magnetocaloric effect we express Eq. (1) as the cumulant and normalize the results to the magnetocaloric effect of a paramagnet. The left-hand side of Fig. 7 shows a contour plot ofΓ mc = Γ mc /(T /H) as a function of the frustration angle φ and the magnetic field H. For small fields H ≪ H sat , Γ mc /(T /H) is small, nearly zero or even slightly negative, apart from finite-size effects showing up in particular in the collinear phase. It is only at the sat-uration field where Γ mc /(T /H) develops a large anomaly peaked slightly above H sat (compare with the right plot of Fig. 3). For magnetic fields H ≫ H sat , we eventually reach Γ mc /(T /H) → 1. Apart from a factor T , the magnetocaloric effect is given by the ratio of two quantities, see Eq.  T = 0.2 J c /k B as a function of the frustration angle φ and the magnetic field H. The heat capacity is large in the disordered regions, reflecting the high number of quasidegenerate states. Around J 2 /J 1 = 1/2 (φ/π ≈ 0.148), a two-peak structure evolves when increasing the field, see also Fig 8. Due to the smallness of the saturation field, we currently cannot say whether such a structure also exists at the "mirrored" (J 2 → −J 2 ) position in the phase diagram at J 2 /J 1 = −1/2. When reaching the saturation field, the heat capacity drops and eventually vanishes.
Taken together, it appears naturally that the magnetocaloric effect is peaked around the saturation field. The drop in magnitude inside the nonmagnetic regions can be understood, too, as a consequence of their large specific heat. And since the entropy rises inside these regions when turning on the magnetic field, the magnetocaloric effect must be negative, indicating a cooling of a sample before reaching the entropy maximum. We note that we also observe a change of sign in (∂S/∂H) T as a function of field in the magnetically ordered regions, and return to this point in the context of spin wave theory below.
On the right-hand side of Fig. 9, we have plotted the values of Γ mc (T, H) at the saturation field as a function of the frustration angle φ, again for The open circles denote the absolute values, while the filled circles denote the values relative to a paramagnet. In accordance with the discussion in the previous paragraph, the absolute values of Γ mc (T, H) in the ordered phases are larger than in the nonmagnetic regions. The deviation from the average value is less than a factor two. In contrast, the normalization to the magnetocaloric effect of the paramagnet (filled circles in the right panel of Fig. 9) introduces a strong influence of the saturation field, compare the right plot in Fig. 3. Therefore, the highest enhancement of Γ mc with respect to a paramagnet occurs deep inside the magnetically ordered phases, where the saturation field reaches its maximum values.
The heat capacity C V (T, H sat ) as a function of the saturation field at constant temperature T = 0.2 J c /k B is plotted on the left side of Fig. 9. It is strongly enhanced in the nonmagnetic regions, while roughly constant as a function of the frustration angle φ in the magnetically ordered phases, giving rise to the comparatively weak φ dependence of Γ mc . As discussed above, the enhancement inside the disordered regions is responsible for the suppression of the magnetocaloric effect.
To clarify the field dependence of Γ mc (T, H) further,

V. APPROXIMATE ANALYTICAL TREATMENTS OF THE MODEL
To better understand the exact numerical results for finite clusters it is useful to have approximate analytical results available for comparison. We consider two approaches: Firstly a mean field treatment which provides a reference point for the global behavior of entropy and specific heat in the ordered phase, and the magnetocaloric effect above the saturation field. Secondly we use a linear spin-wave (LSW) approximation to investigate the anomalous enhancement of the magnetocaloric effect around the saturation field, which turns out to be due to a softening of spin excitations at characteristic wave vectors. In this approximation, in contrast to mean field theory, the MCE below the saturation field is nonzero. The spin wave approximation also allows to study subtle effects for subcritical fields which lead to a sign change of the MCE.
The existence of long range order at finite temperatures, i. e. a non-vanishing transition temperature T c associated with magnetic order, is implicit in both treatments. In reality, for the layered vanadates which we wish to describe, T c will be determined by low energy scales which are not present in our model, notably the interlayer magnetic exchange J ⊥ and magnetic anisotropy δJ.
Formally, we cannot break a continuous symmetry such as spin rotation at any finite temperature in 2D, and to be "correct" we should generalize our model to higher dimension and finite anisotropy. However these details make little qualitative (or quantitative) difference for a wide range of temperatures δJ, J ⊥ ≪ T ≪ T c , so we suppress them below. Furthermore, in these calculations we neglect the consequence of interactions between spin waves (see e.g. 22 ) and the break-down of magnetic order at zero temperature on the borders of the CAF phase 3,23 . These effects can be expected to modify the details of critical behavior as a function of magnetic field, but not its broad features, and are left for future investigation.
As discussed in the previous section, the size of the magnetocaloric cooling rate in Eq. (1) is determined by the ratio of the rate of change of entropy with field (∂S/∂H) T,V and its rate of change with temperature (∂S/∂T ) H,V = C V /T . On approaching the critical field both quantities tend to increase sharply, and the resulting increase in Γ mc (H) is a tradeoff between them. It is not immediately obvious for which frustration angle φ the enhancement in Γ mc (H c , φ) should be largest. At modest temperatures, the simple spin wave approximation described below gives considerable insight into this question, and the critical anomalies of the MCE around the saturation field.

A. Calculation of mean field order parameters and thermodynamics
In this section the magnetothermal properties will be investigated in mean field approximation to have a reference for the spin wave and numerical exactdiagonalization methods. The results of the former are, however, not expected to give a realistic description of the MCE. For a unified treatment of AF phases it is advisable to use a four-sublattice description (α, β = A, B, C, D) with each sublattice having N/4 sites for both NAF and CAF. Since we consider only isotropic exchange we may assume without loss of generality that the field is perpendicular to the (xy) plane of the square lattice, i. e., h = hẑ.
In this and the following subsection we refer all extensive quantities like entropy, specific heat etc. to a single site for convenience. The exchange field h ex and total molecular fieldĥ due to Eq. (4) is then given by where the exchange constants J lk αβ are defined per bond. The components of the exchange field h ex and h ex ⊥ which are parallel and perpendicular to the field direction z are related to the respective spin expectation values S and S ⊥ via the equations where the prefactors for the AF and the FM or fully polarized phases (h > h c for any φ) are given by FM: a = z 2 (J 1 + J 2 ); a ⊥ = 0 Then the mean field approximation of the Hamiltonian in Eq. (4) may be written as This also defines selfconsistently the mean field averages via A = Tr{A exp(−βH mf )}/ Tr{exp(−βH mf )} with β = 1/(k B T ). From the above equation the corresponding total mean field internal energy per site U mf = (1/N ) H mf is obtained as Explicitly, using Eq. (21) one obtains NAF, CAF: To calculate thermodynamic quantities the expectation values S and S ⊥ and their temperature derivatives have to be obtained selfconsistently. This is done by diagonalizing H mf which leads to local eigenstates |± = u ± | ↑ + v ± | ↓ where | ↑ , | ↓ are the degenerate free S = 1/2 states. The |± states have energies E ± = E c + ǫ ± given by 1 2 is field independent in the ordered regime. The field dependence of the canting angle θc/2 (counted from the field direction) is also shown. Right panel: Susceptibility χ for three different temperatures. Its value in the ordered regime is T -independent. In all cases the frustration angle is φ = 0.74π corresponding to the CAF choice.
Defining 2ǫ 0 = ∆ = ǫ + − ǫ − as the splitting due to the molecular fieldĥ their coefficients are then obtained as Using these coefficients and the difference of thermal occupation numbers p − − p + = tanh 1 2 β∆ of eigenstates |± one finally obtains selfconsistent equations for the spin expectation values: The selfconsistency is implied via the molecular field expressionsĥ The canting angle θ c /2 of magnetic moments is obtained from minimizing U mf in Eq. (24) The angle is counted from the field-or c-direction and given by tan( where h c =2 S (a +a ⊥ ). Inserting this into Eqs. (27,28) leads to the simple and general result This means that in the ordered phase the molecular field splitting ∆ of spins is field independent up to h c and hence the total moment S is also field independent, i. e. the moment can only be rotated by the field as long as the transverse staggered order exists. This fact has striking consequences for the thermodynamic quantities below h c . For the thermodynamics we also need the temperature derivatives ∂ S i /∂T = −k B β 2 (∂ S i /∂β) (i = , ⊥). They are obtained from Eq. (27) in a straightforward but lengthy calculation and the resulting explicit expressions are given in appendix A. The mean field solution for a CAF value of φ = 0.74 π is shown in Fig. 11. On the left panel the decrease of the staggered OP S x with increasing field and the concomitant increase of the uniform moment S z are shown for small temperature. The total moment S is constant as predicted and practically equal to S=1/2 in the whole field range. This confirms that the moment is simply rotated (canted) by the field without changing its size. The relevant canting angle θ c /2 is also shown in the figure. For moderate temperatures (center panel) the zero field value of S x is already reduced somewhat. As required by Eq. (31) the field still only reorients the moment, i. e., S is a constant less than 1/2 for fields h < h c . Finally for h > h c when the moment is aligned with the field (θ c /2 = 0), the total moment S = S z will be polarized, i. e. increases with field until it reaches the Below the saturation field hc it is field independent because S (see Fig. 11) and also ∆ is field independent in the ordered regime. Center: specific heat dependence on h. For larger temperatures when hc is sufficiently suppressed a Schottky peak evolves above hc. Right: Magnetocaloric enhancement factor. In MF approximation almost no enhancement is visible above hc(T ). Below the critical field the entropy is field independent and hence Γmc drops to zero suddenly according to Eq. (1). asymptotic value of S = 1/2. The right panel shows the susceptibility χ mf for various temperatures. It is constant and T -independent in the ordered regime. These results are qualitatively unchanged for different angles φ.
The desired thermodynamic quantities may now be conveniently obtained from the free energy F mf and internal energy U mf (Eq. 24) of the S = 1/2 system split by the molecular field by an energy ∆(ĥ ,ĥ ⊥ ). The former is given by per site and the entropy S mf = −(∂F mf /∂T ), specific heat C mf V = (∂U mf /∂T ) and susceptibility per site are then obtained as where S ′ , S ⊥ ′ are given in Eq. (A1) of Appendix A. For uncoupled spins a = a ⊥ = 0 and S ′ = 1 4 h cosh −2 1 2 βh which leads to the Schottky specific heat of the two level system. For the magnetocaloric cooling rate Γ mc we need in addition the temperature gradient of the magnetization m mf = gµ B S which is simply given by Then the mean field expression for the cooling rate Γ mc may be obtained from the definition in Eq. (1) using the expressions for C V in Eq. (33) and the temperature gradient in Eq. (33). Thus the solution of the selfconsistent Eqs. (27) for S and S ⊥ and their temperature gradients S ′ and S ⊥ ′ in Eq. (A1) provide all the necessary input for obtaining the magnetocaloric quantities S mf (T, H), C mf V (T, H) and Γ mc (T, H) from the above equations.
In Fig. 12 we show the field dependence of mean field entropy, specific heat and cooling rate as function of field for various temperatures. Obviously S and C V (left and center panel) are field independent below h c caused by the fact that ∆ is constant according to Eq. (31). Consequently, since the cooling rate is proportional to the field-gradient of S (Eq. 1), it suddenly drops to zero below the saturation field h c as seen in the right panel of Fig. 12. In addition it shows that above h c the cooling rate is only slightly enhanced from the paramagnetic value. Although the mean field results are far from realistic, the two main aspects, field independence of S and C V below h c and steplike anomaly in Γ mc at h c still leave their signature in the more advanced spin wave and numerical treatment to be discussed below.
B. The magnetocaloric effect in the linear spin wave approximation

Statistical mechanics and general expressions
In mean field approximation the elementary excitations are local dispersionless spin flips whose energy is determined by the molecular field. This is far from reality especially close to the saturation field when the order parameter breaks down, associated with a softening of the spin excitations at some k-point or even line in the Brillouin zone (BZ). The linear spin wave approximation takes this effect into account and gives a much more realistic description of the MCE. Since spin wave interactions are left out in our approach, however, the singular behavior around h c may be overestimated.
To calculate the magnetocaloric effect of Eq. (1) in spin wave approximation we start from the partition function where the Hamiltonian H is expanded in spin wave coordinates using the Holstein-Primakoff approximation Here E 0 = U mf (T = 0) is the classical (mean field) ground state energy per spin, ǫ λ (h, k) the spin wave dispersion in applied field h, and the sum over k runs over the appropriate magnetic BZ, while λ counts the different spin wave branches within that BZ. The operator α † λk creates magnons with commutation relations In addition to the classical (mean field) ground state energy per spin E 0 there is zero point energy contribution where A(h, k) is the on-sublattice coupling between spins, defined below. For CAF and NAF phases the E zp is always negative; in the FM, where the ground state and spin waves are eigenstates with a single dispersion ǫ(h, k) ≡ A(h, k), E zp vanishes identically. The partition function is essentially that of set of independent simple harmonic oscillators From this we find the internal and free energy per site, using n B (ǫ, T ) = [e ǫ/(kBT ) − 1] −1 for the Bose factor: The entropy per site S = −(∂F/∂T ) follows directly : We can also find the uniform magnetization In the FM phase, where ǫ(h, k) = ω k + h and there is no zero-point term, this simply reduces to with m 0 = ∂E 0 (h)/∂h. Quite generally, we can calculate the MCE as the ratio in Eq. (1) where the magnetization gradient is given by (42) and the specific heat C V = T (∂S/∂T ) = (∂U/∂T ) by From Eqs. (42, 43), Γ mc is obtained using Eq. (1). This reduces the problem to one of evaluating two-dimensional integrals on k for the appropriate spin wave dispersion ǫ λ (h, k). We note that in the special non-interacting case ǫ(h, k) = h, these expressions reduce to those for an ideal quantum paramagnet, with the associated magnetocaloric effect: Incidentally this is a general property of any system for which the partition function Z depends only on H/T , i. e. F = −k B T ln Z(H/T ).

Spin wave dispersion in the ferromagnet/saturated paramagnet
Since spin wave theory assumes broken spin rotation symmetry, we can treat the spontaneously polarized FM phase for h = 0 and the saturated paramagnet for h > h c (φ) on an equal footing. Expanding about the maximally polarized state to O(1/S) we find is the classical ground state energy per spin, which is equal to U mf (T = 0) in Eq. (24). The spin wave dispersion has a single branch given by: where γ(k) = 1 2 (cos k x + cos k y ) and (48) We then have simply i. e. a rigid shift of the entire dispersion with change in magnetic field.
The dispersion will have a single (parabolic) minimum at k = (0, 0) in the FM phase (i. e. for J 1 < 0, J 2 < |J 1 |/2). However in the saturated paramagnetic state, above the critical field h c (φ, J c ), the minimum of the dispersion will be at k = (π, π) where there is a NAF ground state, and at k = (π, 0) (and symmetry points) where there is a CAF ground state. At the classical critical point J 1 = 2J 2 > 0 separating NAF from the CAF, there are line zeros around the zone boundary k x = ±π and k y = ±π. At the classical critical point −J 1 = 2J 2 > 0 separating FM from the CAF, there are line zeros for k x = 0 and k y = 0. Note that these line zeros connect the different minima of the dispersions between which this state must interpolate.

Spin wave dispersion in the canted NAF
Expanding about a canted NAF with ordering vector (π, π) and canting angle θ c we find The classical ground state energy per site is given by which is identical to U mf (T = 0) of Eq. (24). Minimizing this energy fixes the canting angle (Eq. (30)) where (θ c /2) is measured relative to the magnetic field direction z (i. e. the FM has θ c ≡ 0, the NAF θ c ≡ π/2). Exactly the same expression follows from the requirement that the spin wave expansion contains no terms linear in bosons. Eliminating the magnetic field through Eq. (53), we find: These expressions still depend on the applied magnetic field through the canting angle θ c (h). The bilinear form of Eq. (51) (Appendix B) may be diagonalized by a Bogoliubov transformation to give The two-fold degeneracy of the spin wave dispersion of the NAF is lifted by the applied magnetic field. In the (physical) magnetic Brillouin zone centered on k = (π, π), the Goldstone mode is ǫ + (h, k) = 0, while ǫ − (h, k) = h has a finite gap. Since spin waves are not eigenstates, there is now a zero-point energy term in the energy E zp . In order to calculate the rate of change of magnetization with magnetic field in Eq. (42), we also need the field derivative of the dispersion. Considering explicitly ǫ + (h, k), we obtain As far as the periodicity of dispersions is concerned, note that translation by Q=(π, π) leads to therefore the two ǫ ± modes are simply interchanged by translation through the NAF ordering vector Q=(π, π). This means that in simple thermodynamic averages one can work with a single mode in the full square lattice BZ-e. g. ǫ + (h, k)-rather than with the two (physically distinct) modes in the smaller magnetic BZ. The ǫ + (h, k) spin wave dispersion for NAF in the full paramagnetic BZ is shown in Fig. 13 (left panel).

Spin wave dispersion in the canted CAF
The ground state energy per spin is now given by which again is equal to U mf (T = 0) in Eq. (32). Minimizing this energy leads to a canting angle (cf. Eq. 30) Once again we obtain spin waves with a dispersion of the form Eq. (56). After elimination of the field using Eq. (60), the coefficients of the spin wave expansion are given by These coefficients once again satisfy the relation Eq. (58) with Q=(π, 0), i. e. the ǫ ± modes are interchanged under translation through the magnetic ordering vector. Therefore simple averages can again be calculated for a single mode in the full paramagnetic (square lattice) BZ. As in the NAF case, the field gradient of spin wave energies is given by Eq. (57). The ǫ + (h, k) spin wave dispersion for the CAF in the full paramagnetic BZ is shown in Fig. 13 (right panel).

VI. DISCUSSION OF THE ANALYTICAL RESULTS AND COMPARISON WITH THE NUMERICAL FINDINGS
The typical field dependence of entropy, specific heat and MCE in the linear spin wave approximation as calculated from Eqs. (1,39,42,43) are shown in Fig. 14. Entropy, specific heat and MCE are all smooth functions of magnetic field except at the critical field h c at which there is a (2nd-order) phase transition between the paramagnet and canted Néel phases.
The most striking feature of these predictions is the double spike in the MCE at h c . This is a generic feature of a 2nd order phase transition between paramagnetic and ordered phases in applied field 24 and has previously been seen in Monte Carlo simulations of the classical Heisenberg model 12 . This sudden and sharp change in sign of the MCE can easily be understood in terms of the contours of fixed entropy (adiabats), discussed below. It is accompanied by closely related cusps in the entropy and heat capacity, peaked at h c .
In general, as we would expect, entropy and specific heat have a much weaker field dependence in the ordered phase below h c than in the disordered phase above it. The entropy of the (gapped) paramagnetic phase falls FIG. 14: Left panel: Entropy, specific heat and MCE for φ in the CAF regime. The MCE has a maximum at hmax and then changes sign at h0. At the saturation field hc another sign change associated with a spike-like singularity appears due to the field gradient of the entropy. Right panel: Normalized MCE cooling rate Γmc/Γ 0 mc = (H/T )Γmc for five different frustration angles in the NAF (φ = −0.25π, 0) and CAF (φ = 0.25π, 0.5π, 0.75π) regime and shown for moderate fields. The field is normalized to the saturation field hc(φ) given in Fig. 3 (right). Low field maximum of the MCE at hmax(φ) and sign change at h0(φ) are clearly seen to occur for all frustration angles. Pairwise equalities of the MCE are observed due to the symmetries of the spinwave spectrum with respect to φ. rapidly in applied magnetic field, while the Néel phase responds to magnetic field by canting, at nearly constant entropy. As a result the typical (absolute) value of the MCE is much larger above h c than below. These features of the LSW predictions are reminiscent of the mean field theory, as illustrated in Fig. 12.
For sufficiently large h ≫ T, J c we must (and do) recover the response of an isolated paramagnetic spin in either approximation. Where the LSW predictions differ from those of mean field theory is in the singular features at h c , and in the presence of a finite MCE in the ordered phase. This is positive for h → 0, exhibits a shallow maximum at a characteristic field h = h max , changes sign at another characteristic field h = h 0 , before exhibiting the dramatic double spike at h = h c -irrespective of which ordered phase is considered. The absolute values of these critical fields, and the form of the anomalies at h c , do however depend on the structure of the low energy spin spectrum, and therefore on frustration through φ.
The difference between mean field and spin wave results becomes most obvious in a comparison contour plot of the entropy S(T, H) shown in Fig. 15. The former (left panel) has a temperature dependent critical field which vanishes at the (mean field) transition temperature. The entropy does not have any structure below h c due to the field independent molecular field splitting ∆ in Eq.( 31). In spin wave approximation (right panel) the critical field is temperature independent, but the entropy contours show a typical cusp structure around h c with a maximum at h 0 further down which is caused by the excitation of low energy spin waves (cf. Ref. 25,26). According to the definition of Γ mc in Eq. (1) this immediately translates into the sign change of Γ mc at h 0 and its enhancement around h c .
The behavior of the (normalized) MCE in the low to moderate field regime is presented on an enlarged scale in Fig. 14 (right panel) for typical values of φ in the NAF/CAF region. We notice a considerable variation of the characteristic fields h max (φ) and h 0 (φ) with the frustration angle. The maximum MCE, Γ mc (h max ) for k B T /J c = 0.2 is of the order of ten percent of the paramagnetic value Γ 0 mc (h max ). Furthermore a symmetry in the φ dependence is obvious: Firstly the MCE is invariant under reflection at the axis φ = 0.5π or J 1 = 0 when both values of φ lie in the CAF sector. This is obvious from the spin wave dispersion Eq. (56,61) in the CAF regime which is invariant under the simultaneous transformation (J 1 , J 2 ) → (−J 1 , J 2 ) and (k x , k y ) → (k x + π, k y + π). Since the MCE is obtained by integration over the whole BZ, Γ mc is unchanged under sign reversal J 1 → −J 1 . Secondly the MCE for φ = 0 and φ = 0.5π are equal, i. e., it is invariant under the replacement (J 1 , 0) → (0, J 2 ).
Finite size effects prevent the characteristic fields h max (φ) and h 0 (φ) from being identified in FTLM calculations, as shown in Fig. 10. Nevertheless a negative MCE at moderate fields is clearly compatible with the numerical results. In practice, for a finite size cluster, each of the ground state level crossings shown in Fig. 2 shows up as a separate "phase transition" in the FTLM results for the MCE, with associated positive and negative spikes in Fig. 10. For this reason, the sign of the MCE remains ambiguous. However for the saturated paramagnetic phase, where there is no further level crossing in the ground state, the sign of the MCE is correctly resolved and a pronounced enhancement is seen approaching h c (φ) from above. This is in good qualitative agreement with the predictions of LSW theory, where it is clear from Fig. 14 that the largest positive MCE is to be expected just above h c (φ). This maximum arises from the closing of the spin wave gap in the fully polarized phase when h → h + c . It occurs for fields slightly above h c (φ) because temperature acts to "round" the sharp cusp in the entropy contours at the critical field -c.f. Fig. 15. Needless to say, the low energy excitations responsible for the sharply diverging peaks seen in LSW results are not accurately described by a cluster of 24 sites, and so it makes little sense to compare FTLM and LSW predictions at a quantitative level.
While the structure of the MCE is generic to a (second order) phase transition, the details depend strongly on the amount of frustration present. If the frustration angle is deeply within one of the ordered sectors, the softening of the spin waves at h → h + c occurs at the wave vector of the low field AF structure (π, π) or (0, π). However if the frustration angle approaches the transition regions CAF/NAF and CAF/FM, the softening will occur along the whole line in the BZ connecting the wave vectors of the competing structures. This is reminiscent of, but less dramatic than, the situation in certain geometrically frustrated magnets where the spin gap closes simultaneously at h c for an entire branch of excitations across the BZ, leading to the condensation of a macroscopic number of localized magnon modes 12 .
In the present case one should therefore expect a strong enhancement of Γ mc (h → h + c ) for φ close to one of the above boundary regions. The same holds for the specific heat. In Fig. 16 (left panel) we show the peak value C V (h + c , φ) as function of frustration angle. Indeed the specific heat shows a strong enhancement for φ ≃ 0.15π (CAF/NAF) and φ ≃ 0.85π (CAF/FM) of considerable width in φ. This is in good qualitative agreement with the FTLM results for finite clusters presented in Fig. 9 (left panel). The anomaly at the FM/NAF boundary on the other hand is much smaller.
On the right panel of Fig. 16 the corresponding plot for the MCE is shown. The dashed line shows the bare MCE coefficient Γ mc (h → h + c ) as function of φ. It is almost constant except at the phase boundaries where again a large, but much narrower peak appears. This is not immediately obvious since C V enters in the denominator of the expression for Γ mc in Eq. (1). In fact on approaching φ/π ≃ 0.15 the MCE slightly decreases, only very close to the value when the spin wave dispersion softens along the line (π, 0)-(π, π) and equivalent ones in the BZ a very sharp spike appears. This is due to the fact that the spin wave softening leads to a stronger increase of the magnetization gradient (Eq. (42)) as compared to the increase of C V from Eq. (43). The full line in Fig. 16 (right panel) shows the normalized MCE. Aside from the phase boundaries where again sharp peaks appear it is largely determined by the behavior of the saturation field (Fig. 3) Deep within the ordered phases the degree of enhancement of the MCE relative to an ideal paramagnet is chiefly controlled by the saturation field h c (φ). The FTLM and LSW predictions are therefore in excellent qualitative agreement (c.f. Fig. 9 and Fig. 16). However, once again, finite size effects prevent the FTLM method from capturing the full extent of the anomalous enhancement of the MCE in the highly frustrated regions at the borders of the CAF phase. These may in practice be overestimated by LSW theory, since it takes no account of new non-magnetic phases stabilized by fluctuations. None the less we can gain further insight into the strong enhancement of Γ mc (h → h c ) and the singular peak at φ c in Fig. 16 within the LSW approach by expanding the spin wave energies of the fully polarized phase around the incipient ordering vector. Explicitly, for φ in the classical NAF sector we have where q = k − Q is the distance from the NAF vector Q = (π, π) and ∆ h = h−h c is the excitation gap with the NAF critical field h c = 8SJ 1 , The expansion coefficients are given by a 2 = 2a 4 = S(J 1 − 2J 2 ) and a ′ 4 = SJ 2 . This expansion may be inserted into Eqs. (42, 43) and the integration performed approximately analytically. It is assumed that only modes with an energy ǫ q < k B T contribute appreciably to the integral. In performing the integration one has to distinguish two cases: If one is within the NAF sector the second order coefficient a 2 is where Γmc(φ) in Fig. 14 (right panel) is maximal or changes its sign respectively are plotted as function of the frustration angle. In the CAF phase the characteristic fields are symmetric with respect to φ = 0.5π.
nonzero. If one is at the boundary to the CAF regime a 2 = 0 and the dispersion is determined by the mixed fourth oder coefficient a ′ 4 . Performing the integration in this (classical) limit k B T ≫ ∆ h one obtains the approximate expressions for 2J 2 = J 1 (NAF/CAF) In the corresponding quantum limit k B T ≪ ∆ h , both the heat capacity and the rate of change of magnetization with temperature have an activated behavior. However this cancels in between the numerator and denominator of Eq. (1) to give regardless of the degree of frustration present in the model. We note that this is exactly the form predicted at a quantum critical point on the basis of scaling arguments 24 . Returning to Eq. (64,64)-in the first case, inside the NAF sector, the divergence in the MCE for h → h + c is only of a weak logarithmic type. However at the classical boundary to the CAF sector (which widens into the disordered regime due to quantum fluctuations) the singularity becomes a much stronger one essentially of inverse square root type. This is the reason that the MCE at h = h + c shows a large anomalous peak as function of φ when crossing the NAF/CAF boundary. Because of the vanishing of the second order term in Eq. (62) the dispersion has a saddle point at Q leading to a large DOS of low energy spin waves for 2J 2 ≃ J 1 and therefore a stronger algebraic divergence of the MCE at h c appears. The same arguments hold for the CAF/FM boundary.
However a word of caution is appropriate. In our spin wave calculations we assumed that the classical magnetic phases are stable throughout the phase diagram. Strictly speaking this is not true. As indicated in Fig. 1 in the shaded sectors around the CAF/NAF and CAF/FM phase boundary quantum fluctuations lead to instability of magnetic order and select a different nonmagnetic order parameter, presumably staggered dimer 27 and spinnematic 3 . While the broad features of our theory can be trusted, a truly quantitative theory for the nonmagnetic sectors would require to start from the proper order parameter and their associated elementary excitations. In real materials, sufficiently close to h c (T ), the critical behavior of the MCE will also be sensitive to the details of interlayer coupling and magnetic anisotropy. Both refinements remain as an outstanding challenge.
It is also instructive to track the low to moderate field anomalies in Γ mc (h, φ) as function of φ around the phase diagram. In Fig. 17 the maximum field h max (φ) and the field h 0 (φ) at which the MCE changes sign are plotted as function of frustration angle, normalized to the saturation field h c (φ). While h max (φ)/h c (φ) is rather constant throughout most of the range of angles, h 0 (φ)/h c (φ) shows considerable variation in φ. The two characteristic fields are again symmetric around φ=0.5π or J 1 =0 for the same reasons as explained above. Note that the overall double-minimum structure of h 0 (φ)/ h c (φ) as function of φ prevents its use as a criteria to resolve the ambiguity of frustration angles that appears in zero-field thermodynamic considerations mentioned in Sec. III. It is obvious from the FTLM results in Fig. 10 (right panel), that the temperature dependence of the (normalized) MCE above h c is strongly suppressed as h increases. This effect can also be understood from the LSW calculations. Approximating Eqs. (42,43) for small and large temperatures we obtain the ratio of the low and high temperature (normalized) MCE as function of h > h c : As long as the field is not too far above h c there is still a considerable dispersion in ǫ(h, k) and the above ratio is larger than one (Fig. 10), i. e.Γ mc is T -dependent. Once h ≫ h c however the dispersion is negligible compared to the gap energy ∆ h and the ratio in Eq. (66) approaches one, i. e. we recover the behavior of an ideal paramagnet.
A comparison between FTLM and LSW predictions of the temperature dependence of the MCE for fields safely above h c (in order to avoid the logarithmic singularity at h c ) is given in Fig. 18. There is a reasonable agreement in both magnitude and qualitative T -dependence. Note however that the LSW approximation becomes unreliable when T approaches J c /k B and too many spin wave modes are thermally excited.
FIG. 18: Normalized MCE as function of temperature for above-critical field calculated with FTLM an spin wave (LSW) method. Frustration angles are φ = 0 (NAF) and φ/π = 0.74 (CAF) corresponding to the Sr compound. For T ≫ Jc the FTLM result converges rapidly to 1. The LSW result for finite h/hc has an asymptotic large temperature value different form one, only in the limit h ≫ hc it also approaches one.

VII. SUMMARY AND CONCLUSION
We have investigated the magnetocaloric properties of the J 1 -J 2 model using the FTLM method for finite clusters and spin wave analysis starting from the classical magnetic structures. The one-magnon critical field or saturation field obtained from FTLM agrees well with the spin wave result. Finite size scaling results suggest that close to the CAF/FM boundary, the true critical field is determined by a two-magnon instability. This is consistent with the proposed existence of a spin nematic ground state in this parameter range 3 .
Both FTLM and spin wave results predict a strong enhancement in the low temperature specific heat at the saturation field when the frustration angle crosses the phase boundaries. This may be explained by the large degeneracy of low lying states in these regions. At a constant intermediate field the specific heat exhibits a double peak structure around the NAF/CAF boundary. The entropy and specific heat show only moderate field dependence below the saturation field h c . This feature may already be understood in a mean field approach where the entropy is strictly constant for all h < h c .
Likewise the strong enhancement of the MCE just above the saturation field was investigated. In the FTLM results, the MCE was enhanced by up to a factor ten relative to an ideal paramagnetic (at temperatures small compared to the energy scale J c ). Surprisingly, the largest enhancement (from FTLM, relative to an ideal paramagnet) does not occur at the CAF/NAF boundary, but deep within the magnetically ordered sectors. This can be understood in terms of the anomalous enhancement of the specific heat in the frustrated regions, which enters into the denominator of the MCE (c.f. Eq. (1)).
The overall φ dependence of the MCE enhancement ratio reflects that of the saturation field. This is also true for the spin wave results. There, in addition, the enhancement is sharply peaked on the CAF/NAF boundary appears. This is due to the change of the field-scaling behavior above the critical field from logarithmic to inverse square-root when the boundary is crossed. This feature is due to the appearance of Goldstone modes along a line in the BZ when φ has it critical value. The MCE enhancement may also be directly seen from the isentropics or adiabatic temperature curves which exhibit a large slope above the saturation field.
Below the saturation field the results of the FTLM is strongly constrained by finite size effects, however the spin wave analysis provides considerable insight into the systematic properties of the MCE. In absolute terms the MCE is very small-in fact much smaller than the paramagnetic value. This is a direct consequence of magnetic order, and can easily be understood from the mean field picture of a canted Néel state, for which entropy is constant as a function of field. The spin wave theory, however, predicts a flat maximum in the MCE at low fields, followed by a sign change for subcritical fields. These features are present throughout the phase diagram and the characteristic fields are moderately φ-dependent with a double-minimum structure. For this reason it is not likely that these low-field MCE features are useful in the determination of the frustration ratio.
Considerably above the critical field the temperature dependence calculated from FTLM and spin wave theory show reasonable agreement. At temperatures of the order of J c /k B the enhancement of the MCE is substantially reduced and the behavior of the system crosses over to that of an ideal paramagnet.
The most pronounced discrepancies between FTLM and spin wave analysis appear at the classical phase boundaries. As explained above this may be well understood in terms of the absence or presence of low lying collective modes. While the former method underestimates the MCE anomalies at the phase boundaries, the latter overestimates them-in fact it predicts a singular behavior. A more advanced analytical treatment would have to take into account the quantum nature of the ground state around these boundaries, i. e. stacked dimer (J 2 > 0) or spin nematic (J 2 < 0) and the proper associated excitation spectrum.
The present analysis provides some interesting predictions for the experimental investigation of the class of layered perovskites discussed in the introduction. Specifically we give detailed values for the possible saturation fields which should be easily accessible experimentally for SiO 4 and PO 4 vanadates. According to Fig. 4, these critical fields can be used to resolve the ambiguity in parameterizing the model from its low-field susceptibility and heat capacity, far more cheaply than, e.g. neutron scattering. Furthermore we predict a genuine sign change in the MCE for subcritical fields which should be accessible to experiment.
So far as practical applications-for example in cryogen free cooling-are concerned, the goal is to achieve as large a magnetocaloric effect as possible, at as low a field as possible. Here compounds not too far from the phase boundaries FM/NAF and CAF/FM are the most promising because they combine a significant MCE enhancement with very moderate saturation fields. A detailed treatment of entropy as a function of magnetic field in the nematic phase occurring on the CAF/FM border remains an open challenge. However the high density of low energy excitations and low saturation field of this phase means that it looks a priori very promising for magnetothermal applications.