Thermodynamics of a Quantum Ising system coupled to a Spin Bath

Thermodynamics of a Quantum Ising system coupled to a Spin Bath

R.D. McKenzie, P.C.E. Stamp Department of Physics and Astronomy, University of British Columbia, Vancouver B.C., Canada V6T 1Z1.
Pacific Institute of Theoretical Physics, University of British Columbia, Vancouver B.C., Canada V6T 1Z1.

We study the effect of coupling a spin bath environment to a system which, at low energies, can be modeled as a quantum Ising system. Results are derived for the thermodynamic properties and response functions, both for a toy model and for the system, in which spin-8 electronic spins couple to a spin- nuclear spin bath: the phase transition then occurs in a system of electronuclear degrees of freedom, coupled by long-range dipolar interactions. The quantum Ising phase transition still exists, and one hybridized mode of the Ising and bath spins always goes soft at the transition.

03.65.Yz, 75.45.+j, 75.50.Xx

In both statistical physics and quantum computation the Quantum Ising model plays a central role QIsing (); QIsing2 (). It is key to understanding quantum phase transitions QPT () (QPTs), and describes a large variety of solid-state and atomic spin systems WolfIsing (); KimSimulations (); BlattSimulations (), as well as many quantum computational systems QIPad1Nishimori (); QIPad2Farhi (); QIPad3SantoroSimulations (); QIPad4SCQubits ().

In all this work a key question stands out: what is the effect of a coupling to a thermal bath? In , often considered the archetypal solid-state Quantum Ising system, some experiments have suggested aeppli05 (); ronnow07 () that the mode softening expected at the QPT is suppressed by coupling to a spin bath environment; very recent experiments ronnow16 () have probed transitions between electronuclear modes moshe05 () in . Spin bath modes (which include nuclear and paramagnetic spins, and various solid-state defects) also cause strong decoherence RPP00 (); sham07 (); decoSQUID (); taka11 (), a key problem for quantum computation. In adiabatic quantum computation QIPad1Nishimori (); QIPad2Farhi (); QIPad3SantoroSimulations (); QIPad4SCQubits (), suppression of this QPT would be expected to radically change the dynamics. Thus a lot turns on the question of how a spin bath affects a Quantum Ising system.

In what follows we address this question, using (i) a toy model and (ii) a realistic model for . We find that even though a gap forms in the Ising mode spectrum, a new hybridized mode between the Ising and bath spins appears, which fully softens at the QPT. In many solid-state spin systems this will be an electronuclear mode; we discuss various experiments for probing this mode, including the NMR used in recent experiments ronnow16 () on .

(i) A Toy Model: The Quantum Ising model on its own is defined by one of the simplest Hamiltonians in physics; in terms of Pauli operators , it is written as


where each two level system feels a transverse tunneling field , and is coupled to its neighbours via the longitudinal interactions . The competition between these two terms causes a QPT QPT () between an ordered phase for (where and in the absence of the spin bath ), and a disordered phase for . In adiabatic quantum computation QIPad1Nishimori (); QIPad2Farhi () one slowly varies either or the in time, to pass through the QPT.

Both “oscillator bath” modes feynV63 () (modelling extended degrees of freedom like phonons or photons) and ”spin bath” modes RPP00 () (modelling spatially localised degrees of freedom like nuclear spins or defects) can couple to the Ising spins. We begin by considering a toy model intended to capture the essential effects of coupling to a spin bath. A lattice of central “Ising” spins couples locally to a set of two level “bath” spins (so that on each lattice site we have spin pair states ). Writing the total Hamiltonian as , we take the spin bath term to be


having both longitudinal and transverse interactions. A mean field (MF) solution is easily found; writing


we drop the fluctuations about the mean polarization determined by , leaving the MF interaction acting on individual sites.

Fig. (1) shows the MF results for the polarizations of both Ising and bath spins at temperature , for the case of an exchange coupled ferromagnetic on a simple cubic lattice (). There is clearly a QPT at the the critical transverse field . Note that increases very rapidly with (becoming infinite when ). This result can be explained by a spin bath “blocking” mechanism moshe05 (): at , with no mechanism for flipping the bath spins, the transverse field at any site is not able to mediate transitions between the degenerate states and . The ordered bath spins act as a longitudinal field destroying the QPT. At finite temperatures, thermal bath spin fluctuations restore the phase transition.

Figure 1: Ground state MF polarizations of the Ising and bath spins for the toy model in 3-dimensions (, with being the exchange coupling), as a function of the normalized transverse field . The thick lines denote longitudinal polarizations, and the thin lines denote transverse polarizations. In blue, we see the polarizations for a system with an isotropic bath (Ising and bath spin polarizations coincide). In red, we see the electronic (solid line) and nuclear (dashed line) polarizations for the toy model with an anisotropic bath.

We can also look at the the dynamic susceptibility of the Ising spins , the Fourier transform of . In Fig. (2), using a random phase approximation (RPA), we show the eigenmodes of for an anisotropic hyperfine coupling with a dominant longitudinal component, for . The key features here are (i) a low-energy collective mode which softens to zero energy at the QPT, with sharply-peaked spectral weight near the QPT, and (ii) a clear effect of the QPT on the higher modes as well.

Figure 2: (a) Energies of the zero temperature RPA modes (in units of the exchange coupling ), at , of the toy model, with an anistropic hyperfine interaction, in 3-dimensions () as a function of the normalized transverse field . (b) Intensity of the RPA modes (arbitrary units), with , associated with the longitudinal Ising spin susceptibility. The RPA modes are colour coordinated with their intensities. (c) Energy of the soft mode in the vicinity of the quantum critical point.

The QPT also shows up in the transverse response functions, and the behaviour of the bath spin response functions is similar. At finite there are six eigenmodes, with the additional three corresponding to transitions between excited states; again, the QPT shows up in the low energy hybridized mode between the Ising and bath spins, with sharply peaked spectral weight when .

These MFT/RPA results for the toy model do not yet account for mode-mode interactions between fluctuations about the MF eigenstates - this is done below.

(ii) Real Quantum Ising systems: We consider the 3-dimensional Quantum Ising magnet . This has subtle (and sometimes controversial) experimental properties LiHo (); BB (); kycia (); luke (), many of which clearly depend on the coupling to its thermal environment aeppli05 (); kycia (); aeppli14 (). It differs in four key ways from the toy model, viz. (i) the Quantum Ising spins result from truncation of spin- ionic spins ; (ii) the bath is now made up of nuclear spins , with spin-, not spin-; (iii) in a transverse applied field, the hyperfine coupling actually generates a transverse term acting directly on the bath spins, absent from the toy model; (iv) the inter-spin Ho-Ho interactions are now long-range dipolar. Clearly any one of these features might render the toy model conclusions invalid; thus we must generalize the previous discussion to include all of them.

Diagonalization of the crystal field Hamiltonian produces a low-energy doublet stateaeppli05 (); girvin (). The electronic spin variables can be written in terms of the in the form , with the being Pauli matrices; the are found numerically from the original crystal field Hamiltonian (see Supplementary Information). At the truncated effective Hamiltonian then reads


where with , where , so that , and where is the dipolar interaction coefficient in terms of the Lande -factor . All other terms , with , can be neglected in this system. There are also much smaller exchange interactions which are included in the numerical calculations (see Supplementary Information).

The nuclear spin bath terms take the form


where, even though we start with an isotropic hyperfine interaction , the final hyperfine coupling in the truncated Ising system is highly anisotropic, with terms like which do not conserve spin. We have dropped quadrupolar couplings since they are small, and include only the nuclei in (5); the other nuclear spins have little influence on the QPT. The hyperfine interactions are large aeppli96 (); aeppli05 (); BB01 (); BB03 (), of the same order as the nearest neighbour dipolar interaction, and strongly influence the phase diagram. The new term is not negligible. It fundamentally changes the nuclear spin dynamics, which no longer depend solely on the Ising spin dynamics.

Figure 3: (a) Energies (in Kelvin) of the zero temperature RPA modes of at , as a function of transverse field (in Tesla). The inset shows the electronuclear soft mode in the vicinity of the quantum critical point. (b) Dominant spectral weights , with , associated with the longitudinal Ising spin susceptibility (arbitrary units). The modes shown in yellow carry negligible spectral weight.

We can again calculate MFT/RPA results for the dynamic susceptibility of the spins; in Fig. 3, we show the low-energy longitudinal electronic collective modes of the system, and their spectral weights, as a function of transverse field (the dependence on the angle which makes in the plane is rather weak). The frequency of the lowest electronuclear mode vanishes at the QPT, at the point where the MF magnetization vanishes.

(iii) Role of quantum fluctuations: It is well known that quantum fluctuations can cause MFT/RPA results to fail near any phase transition QPT (); thus any theory of the QPT must include these. To do so we write the partition function for the system as


where is the MF partition function, the average is the imaginary time ordered average over the electronic and spin bath variables taken over eigenstates of , and is the fluctuation correction. Introducing a Hubbard-Stratonovitch fluctuation field , we decouple in the usual way negele () (see Supplementary Information), to obtain a partition function , with an effective Hamiltonian expanded in powers of the fluctuations:


The coefficients are the effective couplings between -tuplets of fluctuation fields, obtained by performing a cumulant expansion of equation (6). The resulting expressions for the are in terms of time ordered MF cumulants of the Ising spins . The (rather lengthy) expressions for the can be derived in terms of the parameters in ; the results for the and appear in the Supplementary Information. The following results are then important:

(a) the 2nd-order term vanishes at the quantum critical point; this is expected since this term describes the RPA fluctuations.

(b) at T=0, the term is always positive - this means that quartic fluctuations do not destabilize the QPT.

(c) as expected, is only non-zero in the ordered phase. Explicit calculation of is required to find corrections to the mean-field phase diagram of order , where is the effective coordination number. These corrections are incorporated into the results shown in Fig. (4), where we illustrate the effect of quantum fluctuations on the longitudinal Ising spin polarization for both the transverse field Ising model, and the toy model with and . In a dipole-coupled system the effect of quantum fluctuations is quite striking. This is because for any given spin, the dipole interaction favours anti-alignment of all other spins in the transverse plane, leading to an enhancement of the fluctuations. The mode-mode coupling coming from gives small corrections to this, as they are of order .

Figure 4: Ground state longitudinal spin polarization of the transverse field Ising model (blue), and the toy model (red) with and (with being the strength of the coupling), as a function of the normalized transverse field . The MF results (solid line) are for a simple cubic lattice with . The leading order corrections to the MF results are calculated using an expansion in the inverse coordination number for both an exchange coupled system (dot-dashed line), and a dipolar coupled system (dashed line). We take the lattice spacing to be equal to unity.

(iv) Experiments: One would like to know how generally applicable are the results just found, and how to test for them experimentally.

Let us begin with . Unlike the toy model, has correlations between the , and components of the low-energy effective electronic operators. These correlations are introduced by the crystal electric field, and manifest themselves as complex matrix elements, having both real and imaginary components, of the effective Ising spin operators. In Fig. (5), we depict the zero temperature spectral weight of the RPA modes of expected from the response. We see that there is very little absorption due to the low energy modes, with the soft mode dominating any absorption that does occur; the response is similar. At the QPT, the weight of the soft mode seen in and vanishes, and only the higher lying crystal field excitations are able to absorb energy; however, the soft mode should be visible at the QPT in , as illustrated in Fig. 3.

In very recent work ronnow16 (), NMR was used to observe the absorption of the low energy electronuclear modes in . In these experiments, the transverse susceptibility was measured. On the basis of the theory given here, one can then make two remarks:

(a) As noted above, the spectral weight of the soft mode in should vanish near the QPT - thus, to see this mode, measurements should focus on .

(b) Although our results are not directly comparable with experiment, which are performed at finite T, our results include fluctuations about the MF. Note that the fits given by Kovacevic et al. ronnow16 () are to MF theory, ie., they do not include quantum fluctuations, and the soft mode at the phase transition is not apparent.

There are other 3-dimensional quantum Ising systems that can be analyzed in the same way as we have done here. Examples are the molecular magnetic system taka11 (); MST06 (), and a number of -based molecular magnets Mn12 ().

A physical realization of a 1-dimensional quantum Ising system, albeit with weak, frustrated antiferromagnetic couplings between chains, is CoNbColdea (); CoNbSi (); CoNbBalents (); CoNbKinross (). The low energy modes have been probed via NMR CoNbKinross (), and the results used to identify scaling regimes predicted in the 1990s QPT (). The energy spectrum of , measured via neutron scattering CoNbColdea (), is gapped at zero wavevector near the critical point, a fact attributed to the weak interchain couplings that cause the system to order at some . However, we emphasize that hyperfine interactions will certainly lead to low energy electronuclear modes, which will need to be included in any theory of the low-energy spectrum. Given the current experimental energy resolution of neutron scattering, the low energy electronuclear soft mode may be indistinguishable from elastic scattering.

(v) Discussion: We find that Quantum Ising systems coupled to a spin bath, with hybridized modes between the Ising and bath spin variables, must still have a QPT. This is true even when we take account of quantum fluctuations around the MFT results, and even for high-spin electronic Ising systems, or high-spin bath variables. Nor does an independent dynamics for the bath variables (coming from, eg., the extra field in the system) change the result. Although long-range dipolar interactions strongly enhance fluctuations around the MF/RPA results, they do not destroy the QPT.

It is then worth asking - when can the spin bath destroy the QPT? One case is when the bath spins are frozen. In our results, we have assumed thermodynamic equilibrium - but at low , bath spin relaxation times to equilibrium can be much slower than the Ising spin dynamics. In this case the bath may act as a random static potential on the Ising spins, giving more complex effects. In interpreting any experiment where the system is swept through the QPT at a finite rate, due attention will have to be paid to this point.

Figure 5: Zero temperature intensities (arbitrary units), with , of the RPA modes of associated with the transverse susceptibility, as a function of transverse field (in Tesla). (a) Intensities of the electronuclear modes corresponding to the lowest lying crystal field excitation (). (b) The inset shows the intensity of the soft mode (red). The intensities of all other modes (yellow) are negligible.

We thank NSERC in Canada for funding, and G. Aeppli, T. Cox and A. Gomez for discussions.

I Supplementary Information

In this Supplement we give (i) the derivation of an effective low temperature Hamiltonian for the system, starting from the known microscopic Hamiltonian, followed by (ii) the details of the mean field calculations, both for the spin-half toy model, and for the system; and finally (iii) a general field theoretic formalism suitable for describing dynamic fluctuations in quantum Ising systems, and the expressions used in the text to find corrections to the mean field thermodynamics.


The microscopic Hamiltonian for the system involves both the spin- electronic spins , which we will truncate to spin- degrees of freedom , and the spin- nuclear spins ; we will ignore here the nuclear spins on the and sites, since their hyperfine couplings to the electronic spins are too weak to affect the thermodynamic properties.

In what follows we give the details of the truncation to the effective Hamiltonian used in the text.

I (a) Truncation of the Electronic Terms: We consider a single-ion Hamiltonian for the spin- electronic component of the ions of the form


The crystal field Hamiltonian has the form


where we use the standard Stevens’ operators:


with and indicating the Hermitian conjugate. For the we use the estimates of Rønnow et al ronnow07 (). The crystal field eigenstates are mixed and split by an applied transverse field , with a Landé g-factor of , and .

We now diagonalize the electronic single ion Hamiltonian , using a unitary rotation , such that , and . We then truncate the operators down to the two-dimensional subspace involving the two lowest eigenstates of ; the original spin operators may then be expressed in terms of Pauli operators operating on the subspace in the form girvin ()


The lower two electronic eigenstates of are separated from the rest of the electronic eigenstates by a gap of at least . The hyperfine interaction, and the interactions between holmium spins, are too weak to cause significant mixing with the higher lying eigenstates, which justifies the truncation procedure. We apply a second rotation in order to diagonalize the operator in the subspace so that . In terms of the two lowest eigenstates of , and , our basis is , and , where the phase is fixed such that the coefficient of the lowest eigenstate is real and positive. In Fig. (6), we plot the non-zero matrix elements of the effective spin half operators as a function of the transverse field.

Figure 6: The non-zero matrix elements of the effective spin half operators, , for the truncated Hamiltonian, as a function of the applied transverse field (in Tesla). The plot on the left shows the larger matrix elements. The matrix elements shown on the right are significantly smaller.

Apart from truncating the single site terms, we must also truncate the term describing the interaction between the electronic spins . This truncated interaction is given by


where is the shape dependent Fourier transform of the dipolar interaction, and


is the Fourier transform of the exchange interaction, incorporating the four nearest neighbour atoms at and . The electronic dipole-dipole interaction is strongly anisotropic; the physical source of this anisotropy is the deformation of the electronic orbitals due to the crystal electric field. The antiferromagnetic exchange interaction between nearest neighbor sites is much weaker; we use the estimate of Rønnow et al, ronnow07 (). In a long cylindrical sample of the strength of the dipolar interaction at zero wavevector is , which should be compared to the exchange energy .

In terms of the effective spin operators, the total electronic Hamiltonian , may be now written in terms of Pauli operators in the subspace as


The terms neglected in this approximation either vanish due to symmetry considerations, or they are significantly smaller () than the terms given in equation (14) (for a discussion of these terms see Tabei ()). The Ising nature of the system is apparent in the truncated Hamiltonian.

I (b) Truncation of the Nuclear Spin Terms: We now reintroduce the nuclear spins by truncating the hyperfine interaction, , down to the lowest two electronic levels (we replace with the effective spin half operator for the subspace). This is done by by applying the same truncation procedure used above for the electronic spins. Keeping only non-zero terms, the result is then (suppressing the field dependence of all operators):






The effective field is a result of the strong hyperfine interaction in ; the physical transverse field shifts the electronic orbitals, leading to a static effective field. Thus we end up with a nuclear dynamics governed both by this static field, and by the time-varying fields coming from the electronic spins.

We have already seen in Fig. (6) how the matrix elements of the effective spin operators depend on the physical transverse field . In Fig. (7) we show dependence of the effective transverse field mixing the electronic Ising spins, along with the dependence of the parameters in the nuclear spin Hamiltonian . The longitudinal term dominates the hyperfine interactions, with a substantial effective transverse field directly mixing the nuclear spins. The remaining parameters in our model are significantly smaller.

Figure 7: The effective transverse field (in Kelvin), mixing the Ising spins in LiHoF, as a function of the physical transverse field (in Tesla). The upper left inset shows the next largest parameters in the LiHoF Hamiltonian: the longitudinal hyperfine coupling , and the effective transverse field acting directly on the nuclear spins. The lower right inset illustrates the magnitudes of the remaining transverse hyperfine parameters: (uppermost line), (middle line), and the stray field , acting on the nuclear spins in the direction transverse to the easy axis and the direction of the applied transverse field (lower line).

I (c) Full Effective Hamiltonian: The truncated Hamiltonian for the system, including both electronic and nuclear spins, then finally takes the form


The important points to take from the low energy effective Hamiltonian are (i) the Ising nature of the system at low temperatures due to crystal electric field, (ii) the anisotropy of the truncated hyperfine interaction, and (iii) the large effective transverse magnetic field acting directly on the nuclear spins. The effective longitudinal hyperfine interaction is ; the transverse component is over ten times smaller (this anisotropy was noted by Mennenga et al in their specific heat measurements in 1984 Mennenga ()). The effective transverse field acting on the nuclear spins is roughly when the physical transverse field is between 3T and 6T. It is this effective transverse field , rather than the transverse hyperfine interaction , that is mainly responsible for the mixing of the nuclear spin states.


Mean field theory (MFT) starts by approximating the effects of the interactions between the Ising spins by a self-consistent mean field; the Random Phase Approximation (RPA) adds Gaussian fluctuations about this mean field. Both the diagrammatic and functional formulations of this approximation have been known for many years pines (); hubbard59 (); Young ().

In what follows, the details of MFT/RPA calculations for the longitudinal response functions are given (a) for the Ising spins in the spin-half toy model, and (b) for the more complex case of the electronic spins in ; we then (c) give general results for the dynamic susceptibility, both longitudinal and transverse, including the response of the spin bath.

II (a) MFT/RPA for the Spin Half Toy Model (Ising Spins): In the case of the toy model, with a spin half Ising lattice coupled to a set of spin half bath spins, we start from the Hamiltonian


and rewrite this as


where , and denotes a thermodynamic average taken with respect to the mean field Hamiltonian , with single site terms of form


The RPA calculation of response functions then consists in treating the fluctuation term in (20) in a Gaussian approximation.

For the MFT analysis we diagonalize , and make use of single-site Hubbard operators , where , with , are the eigenstates of Hubbard (); Yang (). The mean field Hamiltonian is written in this basis as , with being the eigenvalues of the system. Any one-site operator can be written in terms of the ; in particular we write , and , for the longitudinal spin operators. Making use of the operators, we can calculate finite temperature MFT correlators in terms of the differences between the mean field eigenstates, , and the matrix elements of the relevant operators; for example, the connected imaginary time longitudinal electronic spin correlator , at Matsubara frequency , is given at the inverse temperature by


The are thermal factors, where and is the mean field partition function; the 2nd elastic contribution vanishes in the paramagnetic phase, and in the limit .

The RPA electronic spin correlators are given in terms of in the usual way - in particular, is the connected imaginary time longitudinal spin-spin correlator, where as before ; when evaluated at this is


The dynamic susceptibility follows from analytic continuation back to real frequencies under the prescription , with epsilon being a small positive number later taken to zero. This is the result used to produce the figures in the main text.

II (b) MFT/RPA for the system (Electronic Spins): The full truncated Hamiltonian for the system is given by adding together the electronic and the nuclear terms appearing in equations (14) and (15). The mean field part of this is just the sum , where


in which is just the limit of in (12), and where all the other parameters were discussed in section I above.

The development of the RPA forms for the correlators then proceeds just as for the toy model, but now in an enlarged 16-dimensional Hilbert space incorporating the 2 electronic states and the 8 nuclear states for each site. We get a longitudinal RPA spin correlation function


with . Generalizing the development given for the toy model, we end up with the result


II (c) MFT/RPA for the Dynamic Susceptibility: In Sections II (a) and II (b) we have discussed the longitudinal response of the Ising spins for both the toy model and the system. More generally, the dynamic response of a system will include contributions from the Ising spins and the bath spins,


where and may equal , or , and is the ratio of bath spin and Ising spin gyromagnetic ratios. In we have , and even at energies corresponding to the hyperfine splitting the response is dominated by the electronic contribution.

The dynamic response functions follow from the imaginary time correlation functions , which we write as follows:


In the RPA, the correlation functions are given by


where and may refer to either a bath spin , or an Ising spin . The are connected mean field imaginary time correlation functions; for example,


with the average taken with respect to the single site mean field Hamiltonian .

The RPA spectrum of the a system follows from the zeros of ; the spectral weight carried by these modes depends on the particular response function. Quite generally, we may write the spectral weight carried by the RPA mode as


where is the residue of the pole of . This is the expression used for the the calculations in the main text.


In what follows we begin by (a) formulating the field-theoretic formulation of the Quantum Ising system coupled to a spin bath; then (b) we give the cumulant expansion for the partition function for this system and (c) the fluctuation contributions to the magnetization.

III (a) Form of the partition function: Even in the presence of the spin bath, the partition function of a quantum Ising system may be mapped to a theory. Written in the Matsubara formalism, the partition function is given by


where is the imaginary time ordering operator, , and the average over the fluctuation term is taken with respect to , and


where the are electronic spin operators. We develop the fluctuation theory in terms of spin operators, rather than the Pauli matrices, to avoid confusion between the Pauli matrices and imaginary time indices.

Decoupling the interaction using the Hubbard-Stratonovitch transformation, we introduce a scalar field at each site. The partition function is then


We rewrite this as a functional integral over an effective Hamiltonian ; performing a cumulant expansion in order to write as an expansion in Hubbard-Stratonovitch fields (see III (b)), we find


where the functions are effective coupling constants between the fields defined by


with .

In what follows, we will work up to quartic order in the effective couplings, ie., we will write