Grüneisen parameter of quantum magnets with spin gap

# Grüneisen parameter of quantum magnets with spin gap

Abdulla Rakhimov    Zabardast Narzikulov    Andreas Schilling National University of Uzbekistan, Tashkent 100174, Uzbekistan
Institute of Nuclear Physics, Tashkent 100214, Uzbekistan
Physik-Institut, University of Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland
###### Abstract

Within a path-integral formalism (variational Gaussian approximation) we obtained analytical expressions for thermodynamic quantities of the system of triplons in spin gapped quantum magnets such as magnetization, heat capacity and the magnetic Grüneisen parameter . Near the critical temperature, is discontinuous and changes its sign upon the Bose-Einstein condensation (BEC) of triplons. We predict that in the low-temperature limit and near the critical magnetic field , diverges as , while it scales as as the magnetic field approaches the quantum critical point at .

###### pacs:
75.45+j, 03.75.Hh, 75.30.D

## I Introduction

The properties of condensed matter at low temperatures have always been of high interest. Phenomena such as novel types of superconductivity/superfluidity, quantum phase transitions or different types of topological order still fascinate a growing community of researchers. For condensed matter systems, P. Debye and W. F. Giauque independently suggested in 1926 to use the magnetocaloric effect (MCE) of paramagnetic materials to reach temperatures significantly below 1 K. This effect, which describes the temperature changes of a magnetic material in response to an adiabatic variation of the magnetic field, forms the basis of magnetic refrigeration. The observation of a giant MCE even around room temperature, indicating the potential of the MCE for an environment-friendly room-temperature refrigeration, has stimulated additional work (see recent review by Wolf et al. wolf ()).

The magnetocaloric effect (MCE) and the related magnetic Grüneisen parameter,

 ΓH=1T(∂T∂H)S, (1)

quantify the cooling or heating of a material when an applied magnetic field is changed under adiabatic conditions with constant entropy . In such a process the exchanged heat is zero,

 δQ=TdS=T(∂S∂T)HdT+T(∂S∂H)TdH=0, (2)

and hence

 ΓH=−1CH(∂S∂H)T, (3)

where is the heat capacity at constant magnetic field . Experimentally can be directly accessed by measuring the change in temperature at constant entropy upon magnetic field variation using Eq. (1). Mathematically corresponds to the gradient of the temperature in the landscape along an isoentropic line. Another equivalent expression for the Grüneisen parameter using the magnetization ,

 ΓH=−1CH(∂M∂T)H, (4)

can be derived from the grand thermodynamic potential and suitable Maxwell relations using , with the chemical potential, the number of particles, and and pressure and volume, respectively, which we assume in the following to be constant as we restrict ourselves only to the magnetic subsystem.

The Grüneisen parameter is usually discussed in terms of the quantum critical point111Below we consider only the magnetic Grüneisen parameter, where quantum fluctuations play the major role. Although here we concentrate on a mean-field analysis, where these fluctuations are not taken into account, for the sake of comparison, we briefly discuss the quantum critical point properties. First, if the transition occurs at a given , then , where in our notation is a universal prefactor Zhu (); garst (). For example, for a dilute Bose gas in the symmetry-broken state garst (). The temperature dependence of in the critical regime also shows a divergence as with a certain critical index. It is predicted that has a different sign on each side of the quantum phase transition garst (). These divergences and the sign change of are the hallmarks to identify quantum critical points. These properties have been experimentally confirmed by Gegenwart et al., who developed a low-frequency alternating-field technique to measure down to low temperatures gegenw1 (), in order to classify a number of magnetic systems ranging from heavy-fermion compounds to frustrated magnets gegenw2 (); gegenrev ().

There is class of quantum magnets referred to as zero field gap quantum magnets zapf (); giam (). In a subclass of these materials containing dimers of two entities, the spin gap between exited triplet () and singlet ground () states closes beyond a critical magnetic field due to the Zeeman effect. As a result, bosonic quasiparticles (”triplons”) arise, which may undergo a BEC below a critical temperature . Although experimental data on thermodynamic properties are available for many of such systems (for a review, see zapf (); giam ()), quantitative measurements of the MCE and the associated Grüneisen parameter are rare zapf (); dtn2012 (); zvyagin (). Experimentally it is very difficult to explore the behavior of in the zero-temperature limit. This topic has not yet been systematically addressed for these materials, to the best of our knowledge, neither theoretically nor experimentally, with, perhaps, only a single exception dtn2012 (). From simple arguments, the property can be easily considered for non-interacting Bose systems using and . From Eq. (4), all materials belonging to the non-interacting BEC universality class should therefore obey i.e. dtn2012 (). Nevertheless, as the triplon bosonic quasiparticles in the magnetic insulators to be considered here are known to constitute an Bose gas zapf (); nikuni (); giam (), a consideration of the effects of interaction on the Grüneisen parameter is of utmost interest. The aim of the present work is to investigate the properties of for such magnets within a mean field approximation. We will show that and , and demonstrate that changes its sign at the transition.

## Ii The free energy and entropy of the triplon gas

For the thermodynamics of a dimerized quantum magnet is determined by the system of triplon quasiparticles with integer spin if we neglect the phonon contribution for the moment. In a constant external magnetic field, the number of triplons is conserved in the thermodynamic limit, and they can experience a Bose-Einstein condensation (BEC)zapf (); nikuni (); giam (). Although the critical temperature or the density of triplons of the BEC may be obtained within Hamiltonian formalism ourannals (); ourprb (); ournjp1 (); ournjp2 (), the thermodynamic potential and, in particular, the entropy can be also derived by using a Gaussian functional approximation andersen (), which is in fact, equivalent to the Hartree-Fock-Bogoliubov (HFB) approach.

In this formalism one starts with the action

 A[ψ†,ψ]=∫β0dτ∫d3r{ψ†[∂∂τ−^K−μ]ψ+U2(ψ†ψ)2}, (5)

where , is the chemical potential, here given as with the Lande -factor matsubara (); zapf (); nikuni (); giam (), is the operator of kinetic energy, represents a constant for repulsive triplon -triplon interaction, and is the Bohr magneton. The complex fields, and satisfy the standard bosonic periodicity conditions in that and are periodic in with period . The operator gives rise to the bare dispersion of triplons as defined, for example, in the bond operator representation matsumoto (). The integration in coordinate space may be taken in the first Brillouin zone with the volume , which we set here ourlatt (). Then the thermodynamical potential can be obtained from

 Ω=−TlnZ (6)

where the grand-canonical partition function is given by the path integral bellac ()

 Z=∫Dψ†Dψe−A[ψ†,ψ]. (7)

Due to the complications related to the term in (5), the path integral cannot be evaluated exactly. In the present work we shall use a variational perturbation theory klbook () as outlined in Refs. ourjstatmeh (); ourpra77 () for finite systems. Referring the reader to the Appendix A for the calculation details, we obtain for

 Ω=Ωcl+Ω2+Ω4, Ωcl=−μρ0+Uρ202+12∑k(Ek−εk)+T∑kln(1−e−βEk), (8) Ω2=12[A1(Uρ0−X2−μ1)+A2(3Uρ0−X1−μ1)], Ω4=U8(3A21+2A1A2+3A22),

and for = (1,2) or (2,1), respectively,

 Ai=Gjj(τ,r,τ′,r′)∣∣r→r′,τ→τ′=T∑n∑kεk+Xiω2n+E2k=∑kεk+XiEkWk, (9)

where

 Wk=12coth(βEk2)=12+nB(Ek), nB(x)=1ex−1, (10)

with

 Ek=√εk+X1√εk+X2 (11)

being the dispersion relation of the quasiparticles. Here and are variational parameters defined from the principle of minimal sensitivity andersen () as

 ∂Ω(X1,X2,ρ0)∂X1=0,∂Ω(X1,X2,ρ0)∂X2=0.

The normal and the anomalous densities become

 ρ1=∫⟨˜ψ†˜ψ⟩d3r=A2+A12, σ=∫⟨˜ψ˜ψ⟩d3r=A2−A12, (12)

respectively. From (II), (8) and (12) one obtains for and

 X1=−μ+U(2ρ1+3ρ0+σ), (13) X2=−μ+U(2ρ1+ρ0−σ). (14)

The stability condition yields

 μ−Uρ0−2Uρ1−Uσ=0, (15)

where is the condensed fraction summing up to the total density . In general, explicit expressions for all thermodynamic quantities can be inferred from given in (8). In particular, differentiating with respect to temperature yields the entropy

 S=−(∂Ω∂T)H=−∑kln[1−exp(−βEk)]+β∑kEk(eβEk−1), (16)

while the heat capacity in constant magnetic field becomes

 CH=T(∂S∂T)H=β2∑kEk(Ek−TE′k,T)eβEk(eβEk−1)2. (17)

The resulting magnetic Grüneisen parameter is

 ΓH=−gμBCH(∂S∂μ)T=μBgβ2CH∑kEkE′k,μeβEk(eβEk−1)2, (18)

where and , which are given explicitly in the Appendix B.

As we noted above, the present approximation is equivalent to the HFB approximation. Another similar approach, the Hartree - Fock - Popov (HFP) approximation which is widely used in the literature zapf (); yamada (); nikuni (), can be formally obtained from the HFB relations by neglecting the anomalous density, i.e. by setting in the above equations.

For further considerations, we have to discuss the normal () and the condensed phase () of the system separately.

### ii.1 Normal phase, T≥Tc

When the temperature exceeds a critical temperature , the condensate fraction as well as the anomalous density vanish, i.e., , and . In this normal phase both approximations, HFB and HFP, coincide.

The basic equations (13) and (14) have the same trivial solutions as

 X1=X2=2Uρ−μ. (19)

Inserting this into Eq. (11) gives

 Ek(T≥Tc)≡ωk=εk−(μ−2Uρ)≡εk−μeff, (20)

defining the effective chemical potential . Differentiating both sides of Eq. (20) with respect to and using Eq. (17) gives the following expression for the heat capacity:

 CH(T≥Tc)=β2∑kωkeβωk(ωk−2Uρ′T)(eβωk−1)2. (21)

The triplon density, which defines the longitudinal magnetization (i.e., the component parallel to ) via

 M=−∂Ω∂H=−∂Ω∂μ∂μ∂H=μBgρ, (22)

is given by the solution of the nonlinear equation

 ρ(T)=ρ1=A1+A22=∑k1eβωk−1=∑k1e(εk−μ+2Uρ)β−1, (23)

where we used Equations (9), (12) and (19). Note that in this phase, the staggered magnetization , which is a hallmark for the BEC state in dimerized spin systems, vanishes.

For the Grüneisen parameter we have from Eqs. (4) and (22)

 ΓH(T>Tc)=−gμBCHρ′T, (24)

where may be obtained from Eq. (23) (see Appendix B).

The critical density , i.e. the density of quasiparticles at the critical temperature , is reached as soon the effective chemical potential vanishes, and hence

 ρc=ρ(Tc)=μ2U. (25)

With this condition we may obtain the critical temperature as the solution of the equation

 μ2U=∑k1eεk/Tc−1, (26)

which will later be used to optimize the input parameters of the model by comparing experimental data with the calculated dependence.

### ii.2 Condensed phase, T<Tc

In the condensed phase where the symmetry is spontaneously broken, one has to implement the Hugenholtz - Pines pines () theorem relating the normal and the anomalous self energies and to each other, i.e.

 Σn−Σan=μ. (27)

In our notation this leads to the equation ournjp2 ()

 X2=Σn−Σan−μ=0, (28)

or

 μ−U(2ρ1+ρ0−σ)=0, (29)

where we have used Eq. (14). Due to Hugenholtz- Pines theorem, the excitation energy becomes gapless,

 Ek(T

where is the velocity of the first sound for the quasiparticles with effective mass . Eliminating from Eqs. (13) and (29) one obtains the basic equation

 Δ=X12=μ+2U(σ−ρ1), (31)

where 222see ref. redbook () for the origin of the term in (33)

 σ = −Δ∑kWkEk, (32) ρ1 = ∑k[Wk(εk+Δ)Ek−12], (33)

and

 Ek=√εk√εk+2Δ. (34)

Equation (29) with gives the same expression for the critical density as in Eq. (25), which proves the self consistency of this approach. Taking from Eq. (34) with into Eq. (17) gives

 CH(T

where is given in the Appendix B.

For practical calculations Eq. (31) can be rewritten as

 Z=1+˜σ(Z)−˜ρ1(Z), (36)

where , , and . After solving the equation (36), the longitudinal and the staggered magnetizations and in the condensed phase, respectively, become

 M(T≤Tc)=gμBρ=gμBρc(Z+1), (37) M2⊥(T≤Tc)=12g2μ2Bρ0=12g2μ2Bρc(2Z−˜σ),

where we used

 ρ=Δ+μ2U,ρ0=ΔU−σ. (38)

The Grüneisen parameter is with Eqs. (37) and (38)

 ΓH(T≤Tc)=−gμBρ′TCH=−gμBΔ′T2UCH, (39)

where is given in Eq. (35).

The main difference between the HFB and HFP approximations manifests itself in the condensed phase. In particular, the basic equation (31) simplifies to

 ΔHFP=μ−2Uρ1=Uρ0, (40)

where is formally the same as in Eq. (33).

## Iii Low temperature expansion

In the present section we will derive analytical expressions in the limit. We shall perform the low-temperature expansion as a function of the dimensionless parameter . In fact, for the majority of spin gap quantum magnets, the effective mass is small, e.g. for TlCuCl misgu (), so that any power series in the small parameter should quickly converge.

In general, three dimensional momentum integrals e.g in Eq. (33) can not be taken analytically. So, to overcome this difficulty we use Debye - like approximation yuklaser () . To this end the inegration over momentum in Brillouene zone is replaced by the Debye sphere, whose radius is chosen such that to retain the normalization condition

which gives the dimensionless Debye radius , i.e . In practical calculations we use dimensionless momentum variable . Next, we replace the simple symmetric three-dimensional bare dispersion

 εk=J0(3−coskxa−coskya−coskza), (42)

which is frequently used as a model dispersion relation in gapped quantum magnets giam (), by . Then the momentum integration may be approximated as

 ∑kf(εk)=V(2π)3∫π/a−π/af(εk)dkxdkydkz=∫10f(εq)dqxdqydqz≈π2∫Q00q2dqf(εq) (43)

where . As to the phonon dispersion, similarly to the case of optical lattices, one may use long- wave approximation yuklaser ():

 Eq=√εq√εq+2Δ≈cπq (44)

with the sound velocity at zero temperature . With these approximations for the low-temperature limit, most of the integrals can be evaluated explicitly in terms of logarithmic and polylogarithmic functions of the argument , i.e., as a function robinson () . Since decreases quickly with increasing we may expand in powers of to extract a leading term. On the other hand one may also introduce the Debye temperature and make an expansion in powers of as in solid state physics.

We refer the reader to the Appendix B for the further details of the calculation. The final result for the entropy becomes

 S=2π2(˜T)345γ3+O(˜T5) (45)

where and . The derivative of (45) with respect to gives the heat capacity

 CH=TdSdT≈2π2(˜T)315γ3=2π2T315c3, (46)

which is common for interacting BEC systems since its measurement in superfluid helium helium (). Note that for an ideal Bose gas, i.e., for a system of noninteracting particles, the dispersion is not linear but quadratic, and ouriman ().

To find an expression for the Grüneisen parameter, we use Eq. (37) with the relations

 ΓH=−1CH(dMdT)H=−gμBCH(dρdT)H=−gμB2UCHΔ′T. (47)

The expansion for becomes

 Δ′T=−α1˜T−α3˜T3+O(˜T5) (48)

where

 α1=83UUQ20+4c, (49) α3=8U45γ23π2UQ20+12π2c+10Uγ2UQ20+4c2. (50)

Inserting from (46) we find

 ΓH=15gμBα1γ24π2U1˜T2+15gμB(2γα3−α21)8Uπ2γ+O(˜T2). (51)

This is one of the central results of our paper. We will further simplify and discuss it later in the Discussion section (see Eqs. (72) to (74)).

Using Eqs. (37), (B.13) and (B.14), the low-temperature expansions for the magnetizations become

 M=gμBρ(T)=M(0)−gμBα14Uγm˜T2+O(˜T4) (52)

and

 M2⊥=M2⊥(0)−g2μ2Bc(3α1+Um)24Uγ2˜T2+O(˜T4). (53)

Both quantities vary as in the low-temperature limit while for a non-interacting Bose Einstein condensate, one has the dependence.

Finally, we mention that the above relations for thermodynamic quantities given in Eqs. (45)-(53) are also valid in the HFP approximation, but with slightly modified

 α1|HFP=83UUQ20+8c, (54) α3|HFP=16U45γ23π2UQ20+24π2c+5Uγ2UQ20+8c2. (55)

## Iv Properties near Tc

The behavior of thermodynamic quantities in the temperature region is crucial for the nature of a phase transition. According to the Ehrenfest classification, a discontinuity in a second derivative of at with a continuous first derivative indicates that the transition is of second order huangbook (). In the present section we will study , and .

### iv.1 T→Tc+0 region

Here and . From Eqs. (21) and (B.1)-(B.3) we have

 C(+)H=−S3+2US1ρ′T (56)

where with

 S3=−β2c∑kε2keβcεk(eβcεk−1)2,S1=−βc∑kεkeβcεk(eβcεk−1)2,ρ′T=βcS12S2−1,S2=−Uβc∑keβcεk(eβcεk−1)2. (57)

It can be easily shown that

 limT→Tc+0ρ′T=0 (58)

since in this limit, in the denominator of Eq. (57) at small momentum has an infrared divergence, since the integrand behaves as in this limit. while the numerator is finite. Thus, Eq. (56) becomes

 C(+)H=−S3=β2c∑kε2keβcεk(eβcεk−1)2. (59)

From Eq. (58) we may immediately conclude that the Grüneisen parameter at vanishes,

 Γ(+)H=−1C(+)HlimT→Tc+0(dMdT)H=−gμBC(+)HlimT→Tc+0(dρdT)H=0, (60)

in agreement with the prediction of Garst et al. garst (). The entropy of Eq. (16) is with

 S(+)=−∑kln[1−exp(−βcεk)]+βc∑kεk(eβcεk−1). (61)

In the HFP approximation, the equations (56) - (61) remain unchanged, since it coincides for with the HFB approximation.

### iv.2 T→Tc−0 region

Here , and and hence again, i.e., the dispersion is the same on both sides of the critical temperature. For this reason the entropy is continuous at , . The heat capacity and the Grüneisen parameter are

 C(−)H=β2c∑kεk(εk−TcΔ′T)eβcεk(eβcεk−1)2, (62)
 Γ(−)H=−gμBΔ′T2UC(−)H, (63)

where we used the relation . The defined in (Appendix B) for the HFB approximation may be rewritten as

 Δ′T|HFB=βcUS42(2S5+1),S5=−Uβc∑kTc+εkeβcεk−Tceβcεkεk(eβcεk−1)2,S4=−4βc∑kεkeβcεk(eβcεk−1)2. (64)

In the HFP approach,

 Δ′T|HFP=2βcUS12S2+1, (65)

where and are the same as in Eqs. (57).

#### iv.2.1 HFB approximation

From Eqs. (59), (60) and (62) we can express the discontinuities in and as

 ΔCH=C(−)H−C(+)H=−βc∑kεkΔ′Teβcεk(eβcεk−1)2>0, (66)
 ΔΓH=Γ(−)H−Γ(+)H=−gμBΔ′T2UC(−)H>0, (67)

where is given in (64) and in (62).

From Eqs. (LABEL:jumpch) and (67) it is clear that not only but also the Grüneisen parameter has a finite jump near the critical temperature, and therefore the transition is of second order according to the Ehrenfest classification.

It is easy to show that our results satisfy self-consistently the Ehrenfest relation

 ΔCH=−Tc(dHcdT)T=TcΔ(dMdT)T=Tc. (68)

Using equations (26), (37), (58) and (LABEL:jumpch) leads to a modified Ehrenfest relation for the discontinuity in the Grüneisen parameter in triplon systems,

 ΔΓH=ΔCHTcC(−)H(dHc/dT), (69)

which can be easily derived from Eqs. (26), (62), (LABEL:jumpch) and (67).

#### iv.2.2 HFP approximation

Here is given by Eq. (65), where is finite but diverges as it has been shown in the previous section. Thus , and hence

 ρ′T|HFP=(Δ+μ2U)′=0. (70)

Therefore, we may conclude from Eqs. (59), (62) and (70)

 C(−)H|HFP=C(+)H|HFPΓ(−)H|HFP=Γ(+)H|HFP. (71)

In other words, there is no discontinuity in the HFP approximation, neither in the heat capacity nor in the Grüneisen parameter, which is in sharp contrast to the HFB approximation used here, and to experimental heat-capacity measurements.

## V Discussion

### v.1 Sign change of ΓH at Tc

In the previous section, we have shown that at the critical temperature . It is also easy to show that must change its sign there. Using Eqs. (60) and (B.3) we have with a for . Approaching the critical temperature from below where (see Eq. (63)), for .

### v.2 Divergence of ΓH near the transition

Rewriting Eq. (51) around the QCP in the limit in a compact form (see Appendix C), we obtain

 ΓH≈Gt(H−Hc)T2+GrH−Hc, (72)

with

 Gt=10g2μ2BUmQ20π2=5g2μ2Bπ2Gr (73)

and

 Gr=2UmQ20, (74)

where the next higher-order terms are and , respectively. Here we used the relation

 U=4πasm, (75)

where is the s- wave scattering length. The first term in Eq. (72) dominates in a fixed magnetic field for temperatures , with , while the second term dominates in the opposite limit when approaches the QCP from above at a fixed low temperature .

The fact that diverges as at low enough temperatures is one of our main results. Remarkably, the classification of a number of magnetic systems ranging from heavy-fermion compounds to frustrated magnets done by Gegenwart et al. gegenw2 () reveals that the majority of considered systems shows indeed a similar behavior, with some exceptions like dtn2012 (), however.

The phase boundary between the condensed and the uncondensed states in spin gapped quantum magnets, respectively, can be expressed by a power law of the form . As experimental data on insulating spin systems often show (see Table 1 in the next section), and our results are in line with experimental observations.

The behavior being well established gegenw2 (); garst () in the QCP systems, is obviously also realized in the systems discussed in the present work (see Eq. (72)). We note, however, that this relation cannot be directly applied to the continuous systems such as atomic gases where . In this case, renormalization procedures may lead to different dependences.

### v.3 Universality of Gr

We now discuss whether the value of the dimensionless parameter is universal in spin-gapped triplon systems or not. We make use of our result (74), and state that only depends on the product of the material parameters and . We will now show that within our assumptions, and are not really independent of each other. The full width of the model dispersion relation (42) is . The lower and the upper bounds of the gapped lowest magnon band with bandwidth determine, in a crude approximation, the width of the magnetic phase (i.e., the values of and , respectively) by the Zeeman shift of the corresponding lowest and highest lying triplet states, respectively giam (). Although our approach is only valid in the dilute limit near , we can formally extrapolate it to the fully polarized state with at , and then . The state at corresponds with Eq. (37) to at . With Eq. (29) we find , and therefore .

In real systems, however, the dispersion relation will deviate from Eq. (42), and the bandwidth then differs from the value containing the low-energy effective mass . Moreover, the magnon bands are not rigidly shifted by the Zeeman effect in magnetic fields between and matsumoto (), so that can hold only by up to a factor of unity. In this sense, relations (73) and (74) are not strictly universal, but depend on the details of the magnon spectrum. In Table 1 to be presented below, we have indeed a ratio between 3.2 - 6.3, deviating from the value of 3.

### v.4 Numerical results for real systems

In the previous sections, we have given general expressions for , , , and , and elaborated the limiting cases and . We can use these results to numerically evaluate these quantities over the full range of temperatures. In the following we will restrict ourselves to and . To do this, we have to assume a set of realistic material parameters , , and which we take from experimental data for BaCrO, SrCrO and TlCuCl misgu (); nikuni (); kofuprl (); aczelSr (); wangprl () (see Table 1).

To begin with, we show in Fig. 1 the phase diagrams as calculated from Eq. (26) for BaCrO and SrCrO, together with experimental data taken from Refs. kofuprl (); aczelSr (); wangprl (). For our calculation, we fixed , , and , and fitted according to Eq. (26).

In Table 1, we compare the exponent as obtained from a power-law fit according to to our numerically obtained data, with corresponding fits to the experimental data in the same temperature range () and to a range of values for TlCuCl from the literature nikuni (); tanaka (). These exponents are in fair agreement with our expectation .