Perturbative Quantum Monte Carlo Study of LiHoF4 in a Transverse Magnetic Field

0.62 MB·

by  S. M. A Tabei

Perturbative Quantum Monte Carlo Study of LiHoF in a Transverse Magnetic Field
4
S.M.A Tabei,1 M.J.P. Gingras,2,1 Y.-J. Kao,3 and T. Yavors’kii1
1Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada
2Canadian Institute for Advanced Research, 180 Dundas St. W., Toronto, Ontario, M5G 1Z8, Canada
3Department of Physics and Center for Theoretical Sciences,
National Taiwan University, Taipei 10617, Taiwan
Results from a recent quantum Monte Carlo (QMC) study(P.B. Chakraborty et al., Phys. Rev.
B 70, 144411 (2004)) study of the LiHoF Ising magnetic material in an external transverse mag-
8 4
netic field B show a discrepancy with the experimental results, even for small B where quantum
0 x x
0 fluctuationsaresmall. Thisdiscrepancypersistsasymptoticallyclosetotheclassicalferromagnetto
2 paramagnetphasetransition. Inthispaper,wenumericallyreinvestigatethetemperatureT,versus
transverse field phase diagram of LiHoF4 in the regime of weak Bx. In this regime, starting from
n
an effective low-energy spin-1/2 description of LiHoF , we apply a cumulant expansion to derive
4
a
aneffectivetemperature-dependentclassicalHamiltonianthatincorporatesperturbativelythesmall
J
quantum fluctuations in the vicinity of the classical phase transition at Bx = 0. Via this effective
2 classical Hamiltonian, we study the B −T phase diagram via classical Monte Carlo simulations.
x

In particular, we investigate the influence on the phase diagram of various effects that may be at
]
h the source of the discrepancy between the previous QMC results and the experimental ones. For
c example, we consider two different ways of handling the long-range dipole-dipole interactions and
e explore how the Bx−T phase diagram is modified when using different microscopic crystal field
m Hamiltonians. The main conclusion of our work is that we fully reproduce the previous QMC re-
sultsatsmallB . Unfortunately,noneofthemodificationstothemicroscopic Hamiltonianthatwe
- x
t explore are able to provide a B −T phase diagram compatible with the experiments in the small
a x
semi-classical B regime.
t x
s
.
t
a I. INTRODUCTION thesystemfreezesintoan(Ising)spinglassstateataspin
m
glass critical temperature T 6,7. Similarly to the previ-
g
d- A. Transverse Field Ising Model ous example, Tg(Γ) decreases as Γ is increased until, at
Γ=Γ , a quantum phase transition between a quantum
n c
o Phasetransitionsfromordertodisorderaremostcom- paramagnet and a spin glass phase occurs. Extensive
c numerical studies have found the QPT between a quan-
monly driven by thermal fluctuations. However, near
[ tum paramagnet and a spin glass phase8,9,10 to be quite
absolute zero temperature, a system can, via quantum
interestingdue to the occurrenceofGriffiths-McCoysin-
1 fluctuations associated with the Heisenberg uncertainty
v principle,undergoaquantumphasetransition(QPT)1,2. gularities11,12.
3
The transverse field Ising model (TFIM) is perhaps the
4
simplestmodelthatexhibitsaQPT1,3,4. Thismodelwas
4
0 firstproposedbydeGennestodescribeprotontunneling B. LiHoxY1−xF4
. in ferroelectric systems5. The Hamiltonian of the TFIM
1
is given by
0 The magnetic insulator LiHoF , with a magnetic field
4
8 1 B applied perpendicular to the Ising z direction of
H = J σzσz Γ σx , (1) x
0 TFIM −2 ij i j − i the Ho3+ magnetic moments, is a well known exam-
: i,j i
v X X ple of a physical realization of the transverse field Ising
i where σµ (µ = x,y,z) are the Pauli matrices. Since model13,14,15,16. In LiHoF the predominant J interac-
X i 4 ij
σx and σz do not commute, a nonzero field Γ, trans- tion between the Ho3+ ions is the long range interaction
ar veirse to tihe Ising zˆ direction, causes quantum tunnel- between magnetic dipoles which decays as 1/ri3j, where
ingbetweenthespin-upandspin-downeigenstatesofσz, r is the distance between the i and j ions. The sign of
i ij
hence causing quantum spin fluctuations. These fluctu- J depends on the position of j respect to i. The exis-
ij
ations decrease the critical temperature T at which the tence of a large crystal field anisotropy on the magnetic
c
spins develop long-rangeorder. In the simplest scenario, Ho3+ ions16 causes the system to behave as a classical
where J > 0, the ordered phase is ferromagnetic3,4. Ising system with dipolar interactions for zero applied
ij
At a critical field Γ , T vanishes, and a quantum phase magneticfieldB . Thereasonisthatthesingleioncrys-
c c x
transition between the quantum paramagnet (PM) and tal field ground state is an Ising doublet, meaning that
a long-range ordered ferromagnetic state occurs. The the matrix elements of the raising and lowering angular
H canbegeneralizedbyconsideringJ asquenched momentumoperatorJ± vanishwithinthespacespanned
TFIM ij
(frozen) random interactions. Competing ferromagnetic by the two states of the doublet. The Ising direction
J > 0 and antiferromagnetic J < 0 couplings gener- is parallel to the c axis of the body centered tetrago-
ij ij
ates random frustration. For a three dimensional case, nal structure of LiHoF . In zero applied magnetic field
4
2
B ,thesystemiswelldescribedbyalow-energyeffective
x
spin-1/2 classical dipolar Ising model17,18. Because the 3
QMC, A=0
energy gap between the grounddoublet and the first ex-
citedsingletisfairlylargecomparedtotheJ couplings, 2.5 QMC, A=39 mK
ij
thereislittle quantummechanicaladmixingbetweenthe Experiment
ground doublet and the excited state induced by the in- 2
teractions17. However,anonzeroB admixestheground
x
doublet with the excited singlet and splits the ground T) Ferromagnet
doublet. It is this energy splitting which corresponds to B (x1.5 Paramagnet
the effective transverse field Γ in the TFIM description 4
of LiHoF in nonzero B 13,18,19. 1
4 x
The Ho3+ ions may be substituted (i.e. randomly di- 0.5 2 B (T)x
luted) by non-magnetic yttrium (Y3+) ions, with very 0 T (K)
0 0.5 1 1.5
little lattice distortion. This allows one to study the
0
effects of disorder on LiHoxY1−xF4 as an example of 1 1.2 1.4 1.6
a diluted Ising model. Depending on the concentra- T (K)
tion x of magnetic ions, the low temperature phase is
either ferromagnetic13,20 or spin glass21,22,23. Interest- FIG. 1: The discrepancy between the experimental13 phase
ingly, paradoxical behaviors are observed when a trans- diagram of LiHoF4 and quantum Monte Carlo (QMC) sim-
ulations using stochastic series expansion for small B from
verse magnetic field is applied to LiHo Y F , with x
x 1−x 4
Ref. [19]. The whole phase diagram is shown in the inset.
x < 1. In the ferromagnetic regime, (0.25 < x < 1.0),
Atlow temperatureandhigh B ,neglecting thelarge hyper-
when B = 0, a mean-field behavior T (x) x for x
x c ∝ fineinteractionA,generatesasignificantdiscrepancybetween
the paramagnet to ferromagnet temperature transition
theexperimentalquantumcriticalpointandtheoneobtained
is observed. However, in nonzero B , with increasing
x fromsimulation. However,atlowBxandclosetotheclassical
Bx, Tc(Bx) decreases faster than mean field theory pre- criticalpoint,thehyperfineinteractionisnotaquantitatively
dicts24. For Bx = 0, when LiHoxY1−xF4 is diluted be- importantparameter. Otherpossibilitiesfortheoriginofthis
low x 0.25, a conventional spin glass transition is ob- discrepancy haveto be invokedin this regime.
≈
served14,22,23. The signature of the spin glass transition
is the divergence of the nonlinear magnetic susceptibil-
ity χ3 at Tg25. However, surprisingly, χ3(T) becomes shouldhaveTg(x)>0 for allx>0. However,recentnu-
less singular as B is increased from B =0, suggest- merical33,34 andexperimentalworks23 claimthatafinite
x x
ing that no quantum phase transition between a PM temperature paramagnetic to spin glass phase transition
and a SG state exists as T 014,26. Recently, theo- may not occur for x as large as xc 0.2.
retical studies18,27,28,29 have →suggested that for dipole- ≈
coupled Ho3+ in diluted LiHo Y F , nonzero B gen-
x 1−x 4 x
erates longitudinal (along the Ising zˆ direction) random C. LiHoF4 as a TFIM
fieldsthatcoupletothemagneticmomentand(i)leadto
a faster decrease of T (B ) in the ferromagnetic regime In addition to the phenomena arising in the diluted
c x
and (ii) destroy the paramagnet to spin glass transition regime of LiHo Y F , the x = 1 regime also turns
x 1−x 4
inLiHo Y F samplesthatotherwiseshowaSGtran- out to be interesting. There still exist problems for the
x 1−x 4
sition when B =022,23. Recently, for the ferromagnetic pure LiHoF , requiring the properties of this system in
x 4
regime, the influence of these induced random fields on nonzero B to be re-investigated more thoroughly. Per-
x
the behavior of the linear magnetic susceptibility χ in hapssurprisingly,itisjustrecentlythatthepropertiesof
the presence ofan externaltransversemagnetic field has LiHoF inatransverseexternalmagneticfieldhavebeen
4
been experimentally studied20. When LiHo Y F is studiedinquantitativedetailstartingfromatrulymicro-
x 1−x 4
highly diluted (e.g. LiHo Y F ), very interest- scopic spin Hamiltonian19. In Ref. [19], which reported
0.045 0.955 4
ing and peculiar behaviors are observed. AC suscepti- results from a quantum Monte Carlo (QMC) study us-
bility datashowthatthe distributionofrelaxationtimes ing the stochastic series expansion (SSE) technique35, a
narrows upon cooling below 300 mK 21,30,31. This be- general qualitative agreement between the microscopic
havior is quite different from that observed in conven- model and experimental data13 was obtained. However,
tional spin glasses, where the distribution of relaxation as illustrated in Fig. 1, there is significant quantitative
times broadens upon approaching a spin glass transition discrepancybetweenthe Monte Carloresults ofRef. [19]
atT >06,25. Thisso-calledantiglassbehaviorhasbeen and the experimental data of Ref [13]. In particular,
g
interpreted as evidence that the spin glass transition in the discrepancy between experiment and QMC results
LiHo Y F disappears at some nonzero x > 022,23. persists asymptotically close to the classical ferromag-
x 1−x 4 c
This is in contrast with theoretical arguments32 which netictoparamagneticphasetransition,whereB /T and
x c
argue that, because of the long-ranged 1/r3 nature of quantum fluctuations are perturbatively small. For very
dipolar interactions, classical dipolar Ising spin glasses low temperatures and high B , it is crucial to consider
x
3
the hyperfine interaction in order to explain the behav- such systems (e.g. the choice of the boundary geometry,
ior of the phase diagram close to the quantum critical boundaryconditionsandandthe shape ofthe domains.)
point13,19,36. However,for verysmall B /T , the numer- Finite size effects is another issue that needs to be han-
x c
ical results shown in Fig. 1 indicate that the effect of dled quite carefully in systems where ions interact via
the hyperfine interaction is not important close to the long-range interactions.
classical transition at T . There are different ways to incorporate dipolar inter-
c
It was suggested in Ref. [19] that this discrepancy actions in a computationally efficient way. The method
between simulation and experiment, close to the clas- implemented in Ref. [19] is the reaction field method42,
sical transition, may be related to some uncertainty in which truncates the sum of the long-range interactions
the crystal field parameters (CFP) used in the crystal at the boundary of a sphere. The dipoles outside the
field Hamiltonian, which enters in the TFIM description sphere are treated in a mean-field fashion. Due to
of LiHoF , and which is simulated via QMC. Indeed a the semi mean-field nature of this method, the reac-
4
number of CFP sets obtained from different experimen- tionfield methodoverestimatesthe criticaltemperature.
tal works, such as susceptibility measurements16, neu- In the presence of quantum fluctuations, this overesti-
tron scattering15, and electron paramagnetic resonance mation is still at play and can possibly influence the
experiments37, provide rather different values for the Bx-Tc phase diagram as well. The Ewald summation
CFP. Specifically, different CFP would lead to different method43,44,45,46,47 is another method to treat the long-
field(B )dependenteffectivecouplingparametersinthe range dipolar interactions. In the Ewald summation
x
TFIM description of LiHoF , which would result in dif- method, a specified volume is periodically replicated.
4
ferent B vs T phase diagrams. Then, by summing two convergent series effectively rep-
x c
Yet, there are other factors of strictly computational resenting the dipolar interactions between magnetic mo-
nature which may be at the origin of the discrepancy ments i and j, and all the periodically repeated im-
illustrated in Fig. 1. For example, because of the diffi- ages of j, an effective dipole-dipole interaction between
culties associated with dipolar interactions, calculations two arbitrary magnetic moments i and j within the fi-
incorporating long-range dipolar interactions need to be nite size sample to be numerically simulated is derived.
performed quite carefully. Because of the long-range na- Fromageneralperspective,itwouldappearquiteworth-
ture and angular dependence of dipolar interactions, the whiletoinvestigatetheapplicabilityandusefulnessofthe
dipolar sum U(i) = −1/N j(1−3cos2θij)/ri3j is con- EphwaasleddsiuamgrmamatioofnLmiHeothFod. Itnoddeeetde,rmthieneEtwhaeldloswumBmxavtsioTnc
ditionally convergent38,39,40, i.e the value of the sum de- 4
P method, unlike the reaction field one, is less prone to
pends on the shape of the external boundary of the sys-
mean field over-estimations,and can be used as another
tem studied. Here, r is the distance between site i and
ij
methodology to probe the LiHoF problem via simula-
j,andθ istheanglebetweenr andtheIsingspinaxis. 4
ij ij tions34.
The conditional convergence of dipolar sums has been
Another factor whose influence on the B T phase
studied by Luttinger and Tisza40. They performed the x −
diagram that should be studied is the nearest neighbor
dipolar sum for a number of spin structures for systems
exchange interaction J in LiHoF . The strength of
with different external boundary shapes. For example, ex 4
J , which is expected to be comparable to the dipo-
they considered an infinitely large system of dipoles on ex
lar interactions for a 4f ion such as Ho3+, is unknown.
a body centered cubic lattice. They found that when
The strength can be determined such that the classical
the external boundary is spherical, the ground state is
critical temperature matches the experimental value for
antiferromagnetic,while it is ferromagnetic for a needle-
B =0. The estimatedvalue ofJ is highly sensitive to
shaped sample. Later, Griffiths rigorously proved that x ex
the method used to handle the external boundaries and
for zero external field the free energy for a dipolar lat-
finitesizeeffectsinsimulations,bothofwhichhavesignif-
ticesystemhastobeindependentofthesampleshapein
icant effects when using the reaction field (RF) method,
the thermodynamiclimit41. The immediateconsequence
as already found in Ref. [19].
of Griffiths’ theorem is that in zero external field, the
net magnetization of the sample has to be zero. Other-
wise, the field caused by the magnetic moments sitting
ontheboundaryofthe samplewouldcoupletothe dipo- D. Scope of the Paper
larmomentsofthesample,makingthefreeenergyshape
dependent. Therefore,asaresultofGriffiths’theorem41, The above discussion should make it clear that there
domains must form in the sample, such the total mag- aretworatherdistinctavenuestopursueinordertoseek
netization of the sample is zero in the thermodynamic an explanation for the discrepancy between the experi-
limit. Griffiths’ theorem is at variance with Luttinger mental13 B vs T phase diagramofLiHoF andthe one
x c 4
and Tisza40 results because, in their work, the spin con- obtained via QMC19. One avenue, is that the current
figurationswereassumeduniform,anddomainformation microscopic model is incomplete. As mentioned above,
was neglected. This discussionemphasizes the complica- andsuggestedinRef.[19],onepossiblesourceforthisin-
tion of studying systems with dipolar interactions and completeness may be an inaccurate set of CFP. Another
the caution which should be taken while dealing with possiblesourceisthatotherinteractionsotherthanlong-
4
rangemagneticdipolarinteractionsandnearestneighbor agram by performing solely classical Monte Carlo was
exchange may be at play48. Examples of other interac- an original key motivation for the development of the
tionsincludehigherordermultipoleinteractionsandvir- method presented in this paper.
tual phonon exchange48. The other avenue is related to The rest of the paper is organized as follows. In Sec.
the ensemble of computational pitfalls and insuing nu- II,wereviewthecrystalstructureandthephysicalprop-
merical errors that may arise when one deals with long erties of LiHoF in a transverse field B and the effect
4 x
range dipolar interactions through simulations. There- of the choice of crystal field potential on the magnetic
fore, before one delves into exploring a more complex low energy states. In Sec. III, we introduce the full
microscopic Hamiltonian, there is a clear need to re- microscopic Hamiltonian of LiHoF . We discuss how,
4
investigate the “simpler” problem that solely considers for low energies, an effective spin-1/2 Hamiltonian for
long-range dipole-dipole interactions and nearest neigh- LiHoF can be constructed, and explain how one can
4
bor exchange. picture LiHoF in nonzero B as a dipolar TFIM. We
4 x
In this work we aim to scrutinize the individual role then discuss how a semiclassical effective Hamiltonian
of each of the computational issues as potential culprits is derived from the TFIM Hamiltonian by incorporating
for the discrepancy observed in Fig. 1. Because QMC the transverse field term perturbatively via a cumulant
and experiment do not match at B /T 0, we have expansion. In Sec. IV, we employ the semiclassical ef-
x c
developed a tool that allow us to achieve→the goal in an fective Hamiltonian obtained in the previous section in
efficient and computationally simple way. Since this dis- classical Monte Carlo simulations for small Bx. We dis-
crepancyappearsatlowenoughB neartheclassicalT , cuss the results obtained using either the reaction field
x c
where quantum fluctuations are perturbatively small, we or Ewald summation method for the long-range dipole
canexpandthepartitionfunctionZ intermsofthetrans- interactions. We discuss how Jex is estimated and inves-
versemagneticfieldB ,andrecastthepartitionfunction tigate the sensitivity of the determined value upon the
x
as a sum over strictly classical states, using a new ef- choice of the numerical method. Finally, we compare
fective,albeit temperaturedependent, classical Hamilto- the Bx-Tc phase diagrams originating from two different
nian H (T). In H (T), the quantum effects are incor- sets of crystal field parameters. Section V summarizes
eff eff
porated perturbatively, giving us the ability to calculate our results. The paper also contains three appendices.
allthermodynamicalquantitiesinpresenceofsmallquan- Appendix A discusses details pertaining to the crystal
tum fluctuations within a classicalMonte Carlomethod. field Hamiltonian. Appendix B gives some of the inter-
ThereforeclassicalMonteCarlosimulationscanbeeasily mediate steps needed to construct the effective classical
performed using Heff(T) in a very simple way, without Hamiltonian Heff(T). Finally, Appendix C give the for-
the need to perform complicated QMC19,35 simulations mulaeneededtocalculatephysicalthermodynamicquan-
when interested in a regime with weak quantum fluctua- tities when doing classicalMonte Carlosimulations with
tions49. Therefore,wecanfocusontheregionclosetothe eff(T).
H
classical transition and investigate the different possible
origins of the discrepancy in detail.
II. STRUCTURE AND CRYSTAL FIELD
In summary, (i) the complexity of the QMC SSE
method, (ii) the problematic conditional convergence of
dipolarlatticesums,(iii)the questionofcontrolledfinite The magnetic material LiHoF4 undergoes a second-
sizeeffectsanditsroleontheconsistentdeterminationof order phase transition from a paramagnetic to a ferro-
the nearest-neighbor exchange J , and (iv) the possible magnetic state at a critical temperature of 1.53 K13,16.
ex
sensitivityofthe T (B )dependence onthechoiceofthe The critical temperature can be reduced by applying
c x
CFPaltogetherwarrantanewnumericalinvestigationof a magnetic field Bx transverse to the Ising easy-axis
the T (B ) phase diagraminthe LiHoF transversefield direction. The magnetic field induces quantum fluctua-
c x 4
Ising material. Below, we will show that either fortu- tions such that beyond a critical field of Bc 4.9 Tesla,
x ≈
nately or unfortunately, depending on one’s disposition, the system displays a quantum phase transition from a
thefactorsproposedinSectionICasthepossibleorigins ferromagnetic state to a quantum paramagnetic state
of the discrepancy between experiment and simulation at zero temperature13. The magnetic properties of
(see Fig. 1) are apparently not the issue. Therefore, the LiHoF4 are due to Ho3+ rare earth magnetic ions.
originofthe discrepancyremainsunexplained. However, The electronic ground state of Ho3+ is 4f10, which
theperturbativecumulantMonteCarlotoolthatwehave gives small exchange coupling19,50,51, such that the
devised can be used effectively to searchfor the cause of predominant magnetic interaction between the Ho3+
discrepancy. Without it, the discovery of the irrelevance ions are long-range magnetic dipole-dipole interactions.
oftheabovefactorsthroughaclassicalMonteCarlosim- Hund’s rules dictate that the total angular momentum
ulation would have been a more CPU time consuming of a free ion Ho3+, J = 8 (L = 6 and S = 2) and the
burden. Ultimately, the same tool can also be used to electronic ground state configuration is 5I8. LiHoF4 is
explore the role of the small Bx when x = 020,27,28,29. a compound with space-group C46h(I41/a) and lattice
Indeed, constructing the whole x-T (B ) p6hase diagram parametersa=b=5.175˚A,c=10.75˚A,andhas 4 Ho3+
c x
in the “small B ” vicinity of the classical x-T phase di- ions per unit cell positioned at (0,0,1/2), (0,1/2,3/4),
x c
5
(1/2,1/2,0) and (1/2,0,1/4) 50. The crystal has S rive the required effective model in two steps. Firstly,
4
symmetry, which means the lattice is invariant with in LiHoF , in the temperature range that we are inter-
4
respect to a π rotation about the z axis and reflection estedin,whichiscloseorbelowT (B =0)=1.53K,the
2 c x
with respect to the x y plane. highenergyscalesarewellseparatedfromthelowenergy
−
sector. Theenergyscalefordipolarinteractionsbetween
In the crystalstructure, the Ho3+ ions are surrounded nearest-neighborHo3+ ionsisabout0.31K.Thisismuch
by F− ions, which create a strong crystal electric field smaller thanthe energygapbetween the two firstlowest
with S symmetry. This crystalfield lifts the 17-fold de- single ion energy states and the next higher crystal field
4
generacy of the 5I configuration giving a non-Kramers states (>11 K). In this case, one can neglect the higher
8
ground state doublet. The next excited state is a singlet energy states and reduce the full Hamiltonian Hilbert
with an energy gap of 11 K above the ground state space to a smaller subspace spanned by the two lowest
doublet15,16,37,52. The c≈rystal field Hamiltonian and the energystates. Thisenablesustodeduce alowenergyef-
crystal field parametrization is discussed in more detail fective spin-1 Hamiltonian for LiHoF . Secondly, we de-
2 4
inAppendix A. Holmiumis anisotopicallypureelement rive a semi-classical effective Hamiltonian from this low
with nuclear spin I = 7/2, which is coupled to the elec- energy spin-1 Hamiltonian by incorporating the trans-
2
tronicspinJ viathe hyperfine contactinteractionAI J, versefieldtermperturbativelyviaacumulantexpansion.
·
where A 39 mK50,53. We can then perform a simple classical Monte Carlo us-
≈
ingthissemi-classicaleffectiveHamiltoniantoinvestigate
the small B /T regime.
x c
A. Effective Spin-1 Hamiltonian
2
As mentioned in the previous section, there are three
typeofinteractionsthatplayaroleinthemagneticprop-
erties of LiHoF . The main interaction is the long-range
4
dipole-dipoleinteractionbetweentheHo3+magneticions
denoted by
1
H = (g µ )2 LµνJµJν , (2)
dip 2 L B ij i j
i6=j µν
XX
where µ,ν=x,y,z and J is the total angular mo-
i
mentum of Ho3+ ion i. Lµν is the magnetic
ij
dipole interaction written in the form Lµν =
ij
δµν r 2 3(r )µ(r )ν /r 5, where r is the dis-
ij ij ij ij ij
| | − | |
tancebetweenioniandj. g =1.25istheLand´eg-factor
L
(cid:2)of free Ho3+ and µ = 0(cid:3).6717 K/T is the Bohr magne-
B
ton. The dipolar interaction is complemented by a short
range nearest-neighbor Heisenberg exchange interaction
1
FIG. 2: The crystal structure of LiHoF . NN identifies the H = J J J , (3)
4 exch ex i NN
2 ·
first nearest neighbors and NNN identifies the next nearest
i,NN
X
neighbors
where NN denotes the nearest neighbors of site i.
This exchange interaction is considered to be weak and
isotropic19,54. Thethirdinteractionis the hyperfinecou-
pling between the electronic and nuclear magnetic mo-
III. EFFECTIVE THEORY OF LiHoF4 FOR ments
THE LOW Bx/Tc REGIME
H =A (I J ) . (4)
hyp i i
·
Inthissectionwederiveaneffectivemodelsuitablefor i
X
describing LiHoF4 in a small transverse magnetic field The hyperfine constant A 39 mK is anomalously
regime, where Bx/Tc 0 (Tc is the critical temperature largeinHo3+-basedmaterials≈13,19,36. Thus,thecomplete
→
when Bx = 0). The simplicity gained using an effec- Hamiltonian is written as
tive theory gives us the ability to capture the essential
physics, and to easily reinvestigate the influence of the H = VC(Ji)−gLµB BxJxi
different parameters affecting the behavior of the phase i i
X X
diagram of LiHoF in the B /T 0 regime. We de- +H +H +H . (5)
4 x c dip exch hyp
→
6
Projecting the single ion Hamiltonian of Eq. (6) in this
6
two-dimensional subspace for an arbitrary ion i, we get
5 1
H =E(B ) ∆(B )σx, (8)
T x x
− 2
4
elvin)3 wEhβ(eBrex)E−(BExα)(B=x).12(ETαh(eBexn)e+rgyEβd(iffBexr)e)ncaendbet∆w(eBenx)th=e
K twoloweststatescausedbythetransversemagneticfield
(
∆ B can already be interpreted as an effective transverse
2 x
field Γ = ∆(B )/2 acting on S =1 degrees of freedom
x eff 2
at each site. The dependence of ∆(B ) on the magnetic
1 x
transverse field B is plotted in Fig. 3.
x
Since we are henceforth working in a two-dimensional
0
0 1 2 3 4 5 subspace for each ion i, we can write the interactions
Bx (Tesla) between Jµi and Jνj in terms of effective interactions be-
tween Pauli matrices. Indeed, any operator acting in a
FIG. 3: The energy splitting of the ground state doublet, two-dimensional space can be written as a linear combi-
∆(B )≡E (B )−E (B ),inLiHoF asafunctionofB the nationofσµ Paulimatrices plus the unit matrixσ0 11.
transxverse mβagnxetic fiαeldx. The crysta4l field V was obtaxined In order toiexpress Jµ in terms of σµ, we project J≡µ in
c i i i
from Refs. [15,19]. For more details on the crystal field and the subspace spanned by and . Specifically, we
crystal field parametrization (see AppendixA). write the Jµ operator as |↑i |↓i
Jµ =C σ0+ C (B )σν , (9)
µ0 µν x
The first two terms are single ioninteractions, where V
C ν=x,y,z
X
describesthestrongcrystalfieldinteractionsdiscussedin
where
Section II andAppendix A. The secondterm is the Zee-
man interaction. Henceforth, we ignore H since our 1
hyp C = [ Jµ Jµ ] ,
goal, as explained in the Introduction, is to investigate µz 2 h↑| |↑i−h↓| |↓i
the small Bx and small (Tc(0)−Tc(Bx))/Tc(0) regime C = 1[ Jµ + Jµ ] ,
where,asalreadysuggestedbytheresultsofRef.[19]and µ0
2 h↑| |↑i h↓| |↓i
showninFig.1,thehyperfineinteractioneffectsareneg- 1
ligible. The first two single-site (non-interacting) terms Cµx = [ Jµ + Jµ ] and
2 h↑| |↓i h↓| |↑i
in H, denoted as
1
C = [ Jµ Jµ ] .
µy
H =V (J) g µ B Jx , (6) 2i h↑| |↓i−h↓| |↑i
single−site C L B x
−
Based on the crystal field parameters of Refs. [15,19],
can be easily numerically diagonalized for arbitrary the evolution of the various parameters C and C as
µν µ0
transverse field B 19. α(B ) and β(B ) are the two
x x x
| i | i
loweststatesofthesingleionHamiltonian(6)foragiven
B . TheircorrespondingenergiesaredenotedbyE (B ) 6
x α x
andE (B ). AtB =0these twostatesformadoublet, C
β x x zz
but Bx = 0 lifts the degeneracy. The Ising subspace 5 C
6 |↑i 0 xy
and are chosen by performing a unitary rotation on C
the ||α↓(iBx)i and |β(Bx)i states : eters 4 −−00..0042 y0
|↑i = √12(|αi+exp(iθ)|βi) Paramµν 3 −−00..00860 1 Cy2x 3 4 5 Cx0
1 C B (Tesla)
= (α exp(iθ)β ). (7) x
|↓i √2 | i− | i 2
C
yy
The phase θ is chosen such that the matrix elements of
1
the operatorJz between and isrealanddiagonal,
giving for Jz, Jz = C σ|z↑.iSince|↓tihe first excited state, Cxx
i i zz i 0
γ(B ) , above α(B ) and β(B ) , is at an energy at 0 1 2 3 4 5
| x i | x i | x i B (Tesla)
least seven times higher than k T (B ), and is repelled x
B c x
forallB fromthe α(B ) and β(B ) set(seeFig.1of
x x x
Ref. [19]), we hence|forth nieglect| all exicited crystal field FIG.4: TheevolutionoftheCµν parametersusingthecrystal
statesandworkinareducedHilbertspacespannedsolely fieldVcfromRefs.[15,19]. IntheinsetonecanseethatCxy ≈
C . Coefficients that are not plotted are zero.
by α(B ) and β(B ) , or equivalently by and . y0
x x
| i | i |↑i |↓i
7
a function of B is plotted in Fig. 4. We see that C is tions among the effective S = 1 spins. Via Eq. (7), a
x zz eff 2
the largest term compared to all the other C ’s. specific rotated subspace was chosen, such that C = 0
µν zµ
For the Hamiltonian in Eq. (5), the Jµ operators are (µ = x,y,0; σ0 11). As shown in the inset of Fig. 4,
i ≡
substituted by their two dimensional representations in- C ,C ,andC areverysmall,sotheinteractingterms
xy yx y0
troduced in Eq. (9). This leads to a complicated Hamil- containingthesecoefficientscanbeneglected. Therefore,
tonian that acts within the Ising subspace of and neglecting these terms, we obtain
| ↑i
. The projection generates various kinds of interac-
|↓i
1
H = (g µ )2[C2 (B ) Lzzσzσz +2C (B )C (B ) Lzxσzσx
spin−1/2 2 L B zz x ij i j zz x xx x ij i j
i6=j i6=j
X X
+2C (B )C (B ) Lzyσzσy+C2 (B ) Lxxσxσx+C2 (B ) Lyyσyσy ]
zz x yy x ij i j xx x ij i j yy x ij i j
i6=j i6=j i6=j
X X X
1
+ J C2 (B ) σµσµ +(g µ )2C (B )C (B ) Lzxσz
2 ex µµ x i NN L B zz x x0 x ij i
µ i,NN i6=j
X X X
∆(B )
+ C (B )C (B ) 4J +(g µ )2 Lxx x σx . (10)
 x0 x xx x  ex L B ij − 2  i
i j
X X
   
When the external magnetic field B is zero, only
x
C (0)=0andalltheotherC andC vanish. Hence, 0.06
zz µν µ0
6
Fibnoeradtbuessnecanrticebeleyod,fabaynnauemxstibmeerpnrlaeolfcmlianastgsenirceaatclitcdioifipneoltdlae,rrtmIhsseinasgyresmtezomedreocl1ao9nr. 0.05 ||(g(gLL∆µ∆µBB))22CCxy20yCPxxjPLyijjyL|xijx|
caactniobne,nwehgilcehctiesdpwroitphorrteisopneacltttooCthz2ze(Blexa)dingi6=IjsiLnzigjzσiniztσejzr-. cients0.04 |(gL∆µB)2Cx2xPjLxijx|
As we can see from Eq. (10), for pure LiHoF , an effec- ffi
4
tive σxσx and σyσy pair-wise interactionPs as well as a oe0.03
i j i j C
linear transverse field along the x direction are induced of
in the presence of an external magnetic field. As sug- o0.02
gested by Fig. 5, and already assumed in Ref. [19], we ati
R
expect the quantum fluctuations induced by these terms
0.01
viaeitherdipolarorexchangecoupling,tobequitesmall
and negligible comparedto the quantum fluctuations in-
ducedby∆(B ). Forthepure(disorderfree)LiHoF ,the 0
x 4
0 1 2 3 4 5
invariance of the dipolar interactions under lattice mir-
ror symmetries forces Lzx = 0. So the linear term Bx (Tesla)
j ij
with C (B )C (B ) Lzxσz vanishes. Consider-
zz x x0 x Pi6=j ij i FIG. 5: The ratio of the typical value of terms neglected
ing the C (B )C (B ) Lzxσzσx term, because of
zz x xx xP i6=j ij i j in Hamiltonian (11) respect to ∆, using the crystal field Vc
lattice mirror symmetry, one has Lzxσz σx = fromRefs.[15,19]andthedipolarsumisperformedforalong
P i6=j ij i j
0, therefore this term can only contribute via high cylindrical sample.
P (cid:10) (cid:11)
order fluctuation effects beyond the vanishing mean-
field contribution. Since Cx0(Bx) < 1, we expect the
Czz(Bx)
(second order) fluctuation contribution effects from the
above σzσx term to be small. Hence we neglect the
i j
for the generation of the longitudinal random fields in
C (B )C (B ) Lzxσzσx term in the S = 1 ef-
zz x xx x i6=j ij i j eff 2 LiHo Y F whensubjecttononzeroB 18,27,28,asdis-
fective Hamiltonian H . We should emphasize x 1−x 4 x
P spin−1/2 cussed in the Introduction.
that for diluted LiHo Y F , since the lattice mirror
x 1−x 4
symmetriesarebroken,thetwolatterterms,proportional
to i6=jLzijxσiz and i6=jLzijxσiz σjx , can no longer be Hence, the spin-21 Hamiltonian in Eq. (10) can be
further simplified to a familiar looking transverse field
neglected28. Indeed, these are the terms responsible
P P (cid:10) (cid:11) Ising Hamiltonian with a dipolar and nearest-neighbor
8
exchange Ising interaction. close to the classical critical temperature, we are able
to derive an effective classical Hamiltonian analytically,
1 where quantum fluctuations are incorporated perturba-
H = C2 (B ) (g µ )2 Lzzσzσz
spin−1/2 2 zz x  L B ij i j tively. UsingsucheffectiveperturbativeHamiltonian,we
Xi6=j can then perform classical MC simulations. To set the
 stage, we first consider a general transverse field Ising
∆(B )
+ J σzσz x σx.(11) Hamiltonian such as
ex i NN− 2 i
iX,NN Xi = 1 zzσzσz + 1 σzσz
To simplify the calculations, andin order to be consis- H 2 Lij i j 2Jex i NN
i,j i,NN
X X
tent with the notation of Ref. [19] as well as for further
Γ σx h σz. (17)
comparison between our simulation results and those of − i − 0 i
Ref.[19],welumpthewholeBx dependenceinthetrans- Xi Xi
verse field term into a renormalization factor ǫ(B ) is
x Γ is the transverse field in the x direction and h de-
0
defined as
notesanexternallongitudinalfieldalongthez direction.
C (B ) For compactness, note that we passed from dipolar in-
zz x
ǫ(Bx)= Czz(0) . (12) teractions denoted Cz2z(0)(gLµB)2Lzijz to Lzijz and from
exchange interaction C2 (0)J to ( see Eq. (14) ).
We renormalize the Hamiltonian as zz ex Jex
The partition function Z for a system with Hamiltonian
H =[ǫ(B )]2 , (13) (17) is
spin−1/2 x
H
Z = Trace(e−βH)
with, according to Eq. (11), is e
H
= ψ e−βH ψ , (18)
i i
h | | i
= 1C2 (0) (g µ )2 Lezzσzσz +J σzσz {Xψi}
H 2 zz  L B ij i j ex i NN
i6=j i,NN where Z is obtained by tracing over ψ ’s which are, for
X X i
e −gLµBCzz(0)Bx σix , (14) eaxnadmβple,1d/irkecTt.prWodeuccatnofwσriizteeitgheenvHecatmorilsto(n|↑iainan(1d7|)↓ai)s
B
i ≡
X = H +H . H is the classical part of the Hamilto-
where the renormaleized effective transverse magnetic nHian,fo0rwhic1hthe0ψ ’sareeigenvectors. H Γ σx
i 1 ≡− i i
field , is related to the real applied B via is the quantum term, which does not commute with H .
Bx x P 0
The existence of these two non-commuting terms in
e = ∆(Bx) , (15) prevents us from applying classical Monte Carlo techH-
Bx 2gLµBCzz(0) [ǫ(Bx)]2 niques directly to the system. We canderive aneffective
×
consistentweithRef.[19]. IndiscussingMonteCarlosim- classical Hamiltonian as a functional of ψi, such that
ulations below, we also define a renormalized tempera- e−βHeff[ψi] = ψ e−βH ψ . (19)
i i
ture, T, in conjunction with , with T defined as h | | i
H
Referring to the definition above in Eq. (19), since the
T =[ǫ(B )]2T, (16)
e xe e right hand side of Eq. (19) is the matrix element with
respect to ψ , H [ψ ] is a functional depending only
where T is the real physical temperature. i eff i
e on the set o|f σiz eigenvalues. The partition function can
All results presented in the Monte Carlo simulations i
then be written as a classical partition function
section below were obtained by considering the renor-
malized Hamiltonian (14), and performing the simula-
Z = e−βHeff[ψi] . (20)
tions with respect to the renormalizedT and . Before
x
presenting our Monte Carlo simulations of EBq. (14) as {Xψi}
pertain to LiHoF4, we firstdiscuss the teechniqeue we em- By finding an explicit expression for H [ψ ], one can
eff i
ployedtohandlequantumfluctuationsperturbativelyfor
calculate the thermodynamical properties of the system
small B /T .
x c described by by performing classical Monte Carlo
H
simulations using H instead of .
eff
H
B. Effective classical temperature-dependent
Hamiltonian − perturbation expansion To proceed, we write the matrix element hψ|e−βH|ψi
in terms of a cumulant expansion57
In this section, with a focus on the simplified spin 1 ψ e−βH ψ =
2
Hamiltonian of Eq. (14), we aim to implement a cumu- h | | i
lantperturbativeMonteCarlomethodforaspin 12 trans- exp β ψ ψ + ∞ (−β)n ψ ( ψ ψ )n ψ (.21)
verse Ising model55,56. For small quantum fluctuations, "− h |H| i n! h | H−h |H| i | i#
n>1
X
9
To make the notation more compact, by ψ we mean a transverse field Γ/J = 1, we plot the average thermal
| i
typical ψ eigenvector. Using Eq. (21) we can derive energyas a function of temperature obtainedfrom exact
i
| i
the effective Hamiltonian H [ψ ] perturbatively. The diagonalization,timeclusterQMand“perturbativeMC”
eff i
details of the derivation of H [ψ ] are presented in Ap- using the effective perturbative Hamiltonian described
eff i
pendix B. H [ψ ], is to order O(Γ2), given by above. This tests confirms the quantitative correctness
eff i
ofthe perturbative Monte CarloschemeatsmallβΓ2/J.
Heff = H0+βΓ2 {σizF1[2β(hi+h0)] We also computed other thermodynamic quantities (e.g.
Xi hMzi, hMxi ) and these also compared well with QMC
F [2β(h +h )] . (22) and exact diagonalization results.
0 i 0
− }
InEq.(22),h isthetotallocalfieldaffectingthespinat
i
site i caused by all the other spins, and which is
−0.5
h = zzσz σz , (23)
i − Lij j −Jex NN
j6=i NN
X X
andh istheexternallongitudinalfieldinthezdirection.
0 −1
The functions F (x) and F (x) are defined as
0 1
E/J
cosh(x) 1
F (x) − ,
0 ≡ x2
−1.5 Perturbative MC
sinh(x) x
F (x) − . (24) Time Cluster QMC
1 ≡ x2
Exact Diagonalization
In this effective Hamiltonian, the effect of quantum
fluctuationsistakenintoaccountperturbativelytoorder −2
0 1 2 3 4 5
O(βΓ2/[H ]),where[H ]denotestheorderofmagnitude T/J
0 0
of H , the classical part (first two terms) of Eq. (17).
0
To obtain the thermodynamical properties of the sys-
FIG.6: Energyasafunctionoftemperatureforasimpleone
tem for small transverse fields we can therefore perform dimensional nearest-neighbor Ising chain with a transverse
a classical Monte-Carlo on Heff as a classical counter- field of Γ = J and N = 10 spins and periodic boundary
partof the realquantum mechanicalHamiltonian. Since conditions. The energy is obtained by exact diagonalization
we are interested in thermal averages we can calculate of the Hamiltonian, a time-cluster QMC algorithm, and a
thermodynamical quantities by differentiating the parti- classical Monte-Carlo algorithm of the perturbative effective
tion function, which is written in terms of H [ψ ], with Hamiltonian.
eff i
respect to h , Γ or β. The effective Hamiltonian has
0
BeforewepresentourMonteCarloresultsforLiHoF ,
an explicit h and β dependence. For each true ther- 4
0
let us summarize what we have done so far.
modynamicalquantum-mechanicalquantity,weobtaina
pseudo-operator counterpart. For example the pseudo- 1. Since the spin-spin interactions and T (B ) are
c x
operators corresponding to E , M , M , M2 , and small compared to the gap between the low-lying
M4 are calculated in Apphenidixh Cz,iwhherxeiEh, Mzi and states α(B ) and β(B ) with respect to the ex-
h zi z | x i | x i
M are the energy and magnetization operators along cited crystal field state γ(B ) , we can recast the
x x
| i
the z and x direction. ... stands for the Boltzmann full microscopic model of LiHoF in terms of an
4
h i
thermal average. effective transverse field Ising model with effective
Becauseofitsperturbativenaturein(βΓ),thismethod spin-spin interactions and effective transverse field
is not reliable for large transverse fields or low tem- Γ(B ) that depend on the real physical applied
x
peratures. To illustrate the range of validity of this magnetic field B .
x
method we consider a simple one-dimensional nearest-
2. Since we are interested in a regime where B /T is
neighbor transverse-field Ising-model Hamiltonian H = x c
small, we can develop a perturbation expansion of
J σzσz Γ σx with periodic boundary condi-
− i i i+1 − i i thepartitionfunctioninpowersofB /T andrecast
tions. For a one-dimensional chain of 10 ions, we are x
P P thethermalaveragesofrealphysicalobservablesin
able to calculate the exact total energy of the chain by
terms of quantities that can be determined via a
exact diagonalization. To check our perturbative MC
classical Monte Carlo simulation of a further effec-
technique, we calculated the energy of the Ising chain
tive temperature-dependent classical Hamiltonian.
as a function of temperature for a given transverse field.
To make a comparison, we also performed a quantum Having shown that the perturbative cumulant MC can
Monte-Carlo (QMC) simulation on the system. In this quantitativelydescribetheTFIMforsmallβΓ2/[H ],we
0
QMC simulation, we used the Trotter-Suzuki58 formal- proceed in the next section to describe how we use this
ismandappliedacontinuoustimeclusteralgorithmsim- method to study LiHoF at small transverse field B ,
4 x
ilar to the one in Ref. [59]. In Fig. 6, for a quite large B /T 1.
x c
≪
10
IV. PERTURBATIVE MONTE CARLO STUDY and the external magnetic field acting on the domain is
OF LiHoF4 Bzext, then the susceptibility χ of the domain is
In this section we reportresults from the perturbative χ=Mz/Bzext . (25)
Monte Carlo (MC) simulation to study the low trans-
It should be noted that the macroscopic bulk magne-
verse field B properties of LiHoF , using the low field
x 4 tization , is given by = n g µ Jz , where
perturbative effective Hamiltonian in Eq. (22) and using z z 0 L B
M M h i
n = 4/a2c is the number of dipoles per unit of vol-
Eq. (23) for the definition of the local h fields. As dis- 0
i ume and where a2c is the volume of the unit cell. Using
cussed in the Introduction, our primary goal here is to
Jz =C σz,thebulkmagnetization isrelatedtothe
check the quantum Monte Carlo results from stochastic zz z
total moment of the effective Ising sMpins, M = σz,
seriesexpansionofRef.[19],andinvestigatethecontrast- z i i
in the S =1/2 picture by
ing results with the transversefield B phasediagramof eff
x P
Ref. [13] for small B (See Fig. 1). Hence, we are indeed
x 4 g µ C (B )
interested in LiHoF4 in the case of asymptotically small Mz = N L Ba2zcz x hMzi , (26)
B /T . Thetemperatureweuseinoursimulationsisthe
x c
renormalized temperature defined in Eq. (16). Regard-
where N is the total number of dipoles.
ing Eq. (14), the transverse field Γ used in the pertur-
Let us consider consider an imaginary macroscopic
bative effective Hamiltonian (22) is Γ = gLµBCzz(0) x, sphericalcavitydeepinsideaneedle-shapeddomain. The
B
where isdefinedinEq.(15). Forthelocalfieldh ,de- magnetization inside the sphere should be equal to the
x i
fined inBEq. (23), we have zz = C2 (0)(g µ )2Lzz aend uniform bulk magnetization of the long needle-shaped
Lij zz L B ij
=Ce2 (0)J . domain. Apart from the external magnetic field Bext,
Jex zz ex z
In the following subsections, we first discuss the reac- spinsenclosedinthesphereexperienceanadditionalfield
tion field (RF) and the Ewald summation (ES) methods that originates from the spins on the outer boundary
that we use to deal with the long range dipolar inter- of the imaginary sphere embedded in the long needle-
actions, and discuss how the Monte Carlo results in the shaped domain. The magnetic surface charge density
classicalregime,whereB =0,areaffectedbythechoice on the surface of the needle-shaped domain with uni-
x
ofthe methodweuse. Next, wediscussthe sensitivityof form magnetization z produces an internal magnetic
M
the Jex estimates atzeroBx to finite-size effects, bound- field Bneedle = 4πMz. Meanwhile, the magnetic surface
ary conditions and choice of the method to handle the chargedensityonthesurfaceoftheuniformlymagnetized
dipolar lattice sum. We also consider the effect of differ- sphere with magnetization of z induces a (demagneti-
ent Jex on the phase digram, when Bx 6=0 and Bx/T is zation) magnetic field 83πMz iMnside the sphere that is in
small. Finally,weinvestigatetowhatextentthefinalre- the oppositedirectiontothe appliedfieldandtoBneedle.
sultsdependonthesetofcrystalfieldparameterschosen Therefore,the totalfieldBsph inside the sphericalcavity
z
to describe the Ho3+ single ion properties. is62
8π
Bsph =Bext +4π . (27)
A. Reaction Field Method vs Ewald Summation z z − 3 Mz Mz
Method
is uniform for a bulk sample. Now, instead of con-
z
M
sidering a whole needle-shaped bulk, we can also study
Griffiths’theorem41statesthatintheabsenceofanex- anisolatedsphericalsamplewhichaneffective Bsph field
z
ternalfieldthefreeenergyforadipolarlatticesystemhas is applied to it. If we substitute Bext with /χ and
to be independent of the sample shape in the thermody- Bsph with /χ , where χ is zthe susceMptizbility of
namical limit. Therefore, as an immediate consequence, thze sphericMal dzomsapihn, then wespchan write χ as a function
intheabsenceofanexternalfield,thenetmagnetization of χ
sph
ofthesamplehastobezero. Otherwise,forauniform
M = 0, a shape dependent demagnetization field would χsph
M 6 χ= . (28)
coupletothedipolarmomentsofthesample,makingthe 1 4πχ
− 3 sph
free energy shape dependent. Here, the demagnetization
field is the field originating from the magnetic moments If χ is obtained via some calculation procedure for
sph
sitting on the boundary of the sample. Hence, in the a spherical sample, one can use Eq. (28) to determine
thermodynamic limit, domains formin orderfor the sys- the macroscopic susceptibility of the bulk sample within
tem to have a zero magnetization, =0. which the sphere is embedded. Specifically, simulations
M
Experiments on LiHoF show that the results are canbeperformedonafinitesizesphere,andtheeffectof
4
shapeindependent,confirmingGriffithstheoremanddo- the macroscopic bulk surrounding the sphere is incorpo-
main formation60,61. There is evidence that in LiHoF rated in a mean-field manner by considering an effective
4
long needle-shaped domains form along the c axis60,61. field Bsph interacting with the spins inside the spheri-
z
If we assume that there is a uniform macroscopic bulk cal sample. Using this method, called the reaction field
magnetization within a long needle-shaped domain (RF) method, Chakraborty et al. calculated the finite
z
M