Table Of Content

Theory of triplon dynamics in the quantum magnet BiCu PO 2 6 Kyusung Hwang1 and Yong Baek Kim1,2,3 1Department of Physics and Centre for Quantum Materials, University of Toronto, Toronto, Ontario M5S 1A7, Canada 2Canadian Institute for Advanced Research/Quantum Materials Program, Toronto, Ontario MSG 1Z8, Canada 3School of Physics, Korea Institute for Advanced Study, Seoul 130-722, Korea (Dated: June 20, 2016) We provide a theory of triplon dynamics in the valence bond solid ground state of the coupled spin-ladders modelled for BiCu PO . Utilizing the recent high quality neutron scattering data 2 6 [NaturePhysics12,224(2016)]asguidesandatheoryofinteractingtriplonsviathebondoperator formulation,wedetermineaminimalspinHamiltonianforthissystem. Itisshownthatthesplitting 6 ofthelowenergytriplonmodesandthepeculiarmagneticfielddependenceofthetriplondispersions 1 can be explained by including substantial Dzyaloshinskii-Moriya and symmetric anisotropic spin 0 interactions. Takingintoaccounttheinteractionsbetweentriplonsandthedecayofthetriplonsto 2 the two-triplon continuum via anisotropic spin interactions, we provide a theoretical picture that n can be used to understand the main features of the recent neutron scattering experimental data. u J 7 I. INTRODUCTION interactions. It is found that the low energy proper- 1 ties of the triplons can only be explained in the pres- There has been considerable interest in the emergence ence of substantial Dzyaloshinskii-Moriya and symmet- ] of non-trivial paramagnetic ground states of quantum ric anisotropic spin interactions.16 It is shown that the l e magnets with interacting spin S = 1/2 local moments. interactions between the triplons renormalize the triplon r- Such non-trivial quantum paramagnets would occur in dispersionsandmoreimportantlytheanisotropicspinin- st low dimensional systems or on geometrically frustrated teractions are primarily responsible for the triplon decay . lattices. Twoprominentexamplesofsuchquantumpara- intotwo-triploncontinuum. Ourstudyofthetriplondy- t a magneticgroundstatesarequantumspinliquidsandva- namicsinthissystemprovidesausefulframeworktoun- m lence bond solids.1 The determination of the underlying derstandtherolesofvariousanisotropicspininteractions - Hamiltonian in such systems, however, has been a chal- andpresentsanopportunitytodeterminethespinHamil- d lengeasthisoftenrequiresthedetailedinformationabout tonian of the valence bond solid (VBS) ground states in n o the spectra of the elementary excitations as well as the considerable details. c ground state. The rest of the paper is organized as follows. We in- [ Morespecifically,thedispersionoftheS =1/2charge- troduce our model Hamiltonian for BiCu2PO6 and the neutral spinon excitations or more accurately the two- VBS ground state in Sec. II. The VBS state and its 3 v spinon continuum would be an essential information to triplon excitations are described using the bond oper- 7 determine the Hamiltonian for quantum spin liquids. In ator formulation in Sec. III, where the importance of 4 the case of the valence bond solids, the triplon disper- anisotropic spin interactions is discussed in comparison 7 sionsanddynamicswouldplaysimilarroles. Intwo-and withtheexperimentallymeasuredtriplondispersions. In 3 three-dimensionalquantummagnets,suchinformationis Sec. IV,itisshownthattheanisotropicspininteractions 0 quite scarce and this has remained as one of the main are responsible for the triplon decay into the two-triplon . 1 challenges both in experimental and theoretical studies continuum. Finally our results are summarized with an 0 of these systems. This is in contrast to magnetically or- outlook in Sec. V. 6 deredsystemswhereaccuratedeterminationofspinwave 1 spectrumhasbeenaroundforalongtime,whichisoften : v used to infer the minimal spin Hamiltonian. II. MODEL HAMILTONIAN AND VBS STATE i X In this work, we present a theoretical study of triplon dynamics in the valence bond solid ground state of the r We start by describing our model Hamiltonian for a coupledspin-laddersystem,designedtoexplainmagnetic BiCu PO and the VBS ground state. properties of BiCu PO .2–12 In particular, we investi- 2 6 2 6 gate a possible minimal spin Hamiltonian that is consistent with previous experimental results. The recent high quality neutron scattering experiments reported in A. Hamiltonian Ref. [12]enableustoconstructthemodelusingthevalu- able information on triplon dispersions and two-triplon We consider the lattice structure in Fig. 1 and in- continuum. Combining the theoretical results obtained troduce the spin model described in Eq. (1) with two in the bond operator formulation13–15 of the spin model types of anisotropic spin interactions, which are the and the spectra of the collective modes measured in Dzyaloshinskii-Moriya (D ) and anisotropic & symmet- ij the experiments, we determine various anisotropic spin ric(Γµν)interactions, aswellastheisotropicHeisenberg ij 2 (i,j) J Da Db Dc ij ij ij ij (3,1) J Da Db Dc 1 1 1 1 (1(cid:48),3) J Da −Db Dc 1 1 1 1 (4,2(cid:48)) J Da −Db Dc 1 1 1 1 (2,4) J Da Db Dc 1 1 1 1 (1,1(cid:48)) J Da 0 Dc 2 2 2 (4(cid:48),4) J Da 0 Dc 2 2 2 (2(cid:48),2) J(cid:48) D(cid:48)a 0 D(cid:48)c 2 2 2 (3,3(cid:48)) J(cid:48) D(cid:48)a 0 D(cid:48)c 2 2 2 (1(cid:48)(cid:48),2) J 0 Db 0 3 3 (4,3(cid:48)(cid:48)) J 0 Db 0 3 3 (1,2) J 0 Db 0 4 4 (4,3) J 0 Db 0 4 4 FIG. 1. (Color online) Lattice structure of the spin-1/2 moments in BiCu PO . 2 6 TABLE I. Coupling constants {J ,D } determined by the ij ij crystal symmetry Pnma of BiCu PO . The table lists the 2 6 interactions (Jij). coupling constants at twelve interaction links in a unit cell. DM vectors are decomposed along the a,b,c axes, i.e. D = H=(cid:88)(cid:0)JijSi·Sj +Dij ·Si×Sj +ΓµijνSiµSjν(cid:1), (1) Diajaˆ+Dibjˆb+Dicjcˆwithaˆ,ˆb,cˆbeingorthonormalvectorsaliojng i>j the crystallographic axes. Listed sites, i,j, are denoted with numbers in Fig. 1. The other links on the lattice can be where Si is a S = 1/2 moment at site i, µ,ν ∈ {x,y,z}, generatedbyactinglatticetranslationsonthetwelvelinksin and summation convention for repeated Greek indices is the table. assumed. We set the x,y,z directions along the crystallographic a,b,c axes (Fig. 1), respectively. Based on the crystal symmetry Pnma of BiCu2PO6,2,3 there are five We constrain the coupling constant matrix Γµijν of the anisotropic&symmetricinteractionbyrequiringthefol- independent links as denoted with different colors in the lowing condition: figure. Accordingly, there are five independent Heisen- berg interactions: J1 along the zigzag legs (magenta), J2 DµDν δµνD2 and J(cid:48) along the straight legs (green and orange, respec- Γµν = ij ij − ij . (2) 2 ij 2J 4J tively), J on the rungs of the ladders (blue), and J as ij ij 4 3 theinter-laddercouplings(gray). DistinctionbetweenJ2 Hereunderlyingassumptionisthatthespinexchangein- and J2(cid:48) arises from the existence of two inequivalent Cu teraction is generated by the antiferromagnetic superex- sites (denoted with different colored balls). The antifer- change mechanism. Then, the condition comes from the romagnetic couplings J1, J2, J2(cid:48) form triangular struc- factthatbothofD andΓinteractionsoriginatefromthe tures on the lattice and generate magnetic frustration in spin-orbit coupling in the microscopic Hubbard model the system. Superexchange pathways are given by Cu- (seeRef. [16]orAppendixA).Asmentionedearlier,J is 3 O-Cu for J1 and J3, and Cu-O-O-Cu for J2, J2(cid:48), and not dominated by the antiferromagnetic superexchange. J4.2,3,5 As pointed out in Refs. [3] and [5], the Cu-O-Cu Hence, Eq. (2) is not applied to the J3 links (gray in bondanglefortheJ3 exchangeiscloseto90◦ incontrast Fig. 1). We find that D3 and Γ3 are not essential for to the other exchange interactions, implying that J3 can describing overall magnetic anisotropies in the system so beweakantiferromagneticorferromagneticaccordingto thatwewillignoreD andΓ afterwards(D =Γ =0). 3 3 3 3 the Goodenough-Kanamori rule.17 We assume that the We will investigate the magnetic field response of the coupling J3 is weak in magnitude compared to the other system later. In this case, we consider the Zeeman inter- exchange interactions, J1, J2, J2(cid:48), J4, dominated by the action: antiferromagnetic superexchange. (cid:88) The Dzyaloshinskii-Moriya (DM) vectors {Dij} can HZ =−gµBH· Si, (3) alsobedeterminedbythecrystalsymmetry.16Welistthe i symmetry-constrained DM vectors in Table I for twelve with H being the magnetic field. In principle, we could links in a unit cell. In the table, DM vectors are decom- considersymmetry-allowedg-tensorsfortwoinequivalent posed along the a,b,c axes: Dij = Diajaˆ+Dibjˆb+Dicjcˆ Cu2+ ions to allow anisotropy in the Zeeman interaction with aˆ,ˆb,cˆbeing orthonormal vectors along the crystal- as well. However, this will introduce more parameters lographic axes. The minus signs and zero values in the neededtobedetermined,renderingthetheorymorecom- vector components arise from the pseudo-vector nature plicated. For simplicity, we assume an isotropic g-factor ofD ,andmirror&inversionsymmetriesinthesystem. with g =2. ij 3 phase. III. BOND OPERATOR FORMULATION Bondoperatortheory13–15isausefulframeworkforde- scribingavalencebondsolidphaseanditsspin-1triplon excitations. The theory is built upon the bond operator representation of the spin operators S forming a L,R valence bond. (cid:16) (cid:17) Sα = 1 s†t +t†s−i(cid:15) t†t , L 2 α α αβγ β γ (cid:16) (cid:17) (4) Sα = 1 −s†t −t†s−i(cid:15) t†t , R 2 α α αβγ β γ whereα,β,γ ∈{x,y,z}, and(cid:15) isthetotallyantisym- αβγ FIG.2. (Coloronline)Dimercoveringoftherung-VBSphase. metric tensor. The bond operators s† and t† create the α Dimers (thick black lines) denote the valence bonds formed spin-singlet and spin-triplet states of S , respectively, L,R along the link. The figure also shows the convention for the andfollowthebosonstatistics. Inordertokeepthephys- dimerindex(1,2)andthespinindex(L,R)withinadimerin icalHilbertspaceconsistingofthefourstates(thesinglet the bond operator theory. and triplet), the hardcore constraint should be imposed: s†s+t†t =1. α α We express the original spin Hamiltonian Eq. (1) in B. Valence bond solid terms of the bond operators by using the representation Eq. (4) with the dimer covering for the rung-VBS state Now we discuss the valence bond solid phase as the in Fig. 2. In the resulting bond operator Hamiltonian, ground state of BiCu PO . The VBS phase is depicted the rung-VBS state can be described by condensing the 2 6 in Fig. 2. Here, the valence bonds are formed at the s-bosons at all the dimers: (cid:104)s(cid:105) = (cid:104)s†(cid:105) = s¯. The hard- J -links or the rungs of the spin ladders (denoted with core constraint is incorporated at each dimer with the 4 thick black lines). We call this phase the rung-VBS in Lagrange multiplier µ as −µ(s†s+t†αtα−1). Exploiting this paper. the crystal symmetry of BiCu2PO6, we set the parameters {s¯,µ} to be uniform across all the dimers. Then, The existence of a VBS phase in BiCu PO has been 2 6 we end up with the following form of the bond operator hinted through earlier studies. A finite spin gap and Hamiltonian. elementary spin-1 excitations in the compound are evi- dencesforaVBSstate. Thefinitespingapwasobserved H+H =H +H +H . (5) Z quad cubic quartic in various experiments such as the heat capacity, mag- Here,theHamiltonianisarrangedaccordingtotheorder netic susceptibility, NMR, and neutron scattering.2–5,11 of the t-boson operator. H consists of the quadratic The spin-1 excitations have been detected in recent in- quad terms like t†t, t†t†, and their Hermitian conjugates. elastic neutron scattering experiments.11,12 The spin-1 H contains the cubic terms t†t†t and t†tt. H character of the excitations was inferred by investigat- cubic quartic has the quartic terms of the form t†tt†t. In the above ing their behaviours under external magnetic fields. The expression, we also included the Zeeman interaction [Eq. structure of the VBS state was investigated in the elab- (3)],whichonlyhasquadratictermssincepossiblelinear orate work by Tsirlin et al.5 They constructed a Heisen- terms are cancelled. berg spin ladder model and studied it by using various Below, we develop a simple quadratic bond operator numerical techniques and experimental informations. It theory by keeping only the quadratic part H in the was shown that the ground state of the Heisenberg spin quad Hamiltonian. Via this quadratic theory, we describe the ladder model has strong bond strength (cid:104)S ·S (cid:105) at the i j low energy triplon excitations around the spin gap ob- rungs in the exact diagonalization studies. This numeri- served in experiments. We will consider higher order in- calresultsuggeststhattherung-VBSphasearisesintheir teractions later in this paper. model, wheretheydidn’tconsideranisotropicspininter- actions. The spin model in Eq. (1) is a generalization of the Heisenberg spin ladder model with anisotropic spin A. Quadratic Hamiltonian interactions, which turns out to be extremely important to describe the neutron scattering data. The quadratic bond operator Hamiltonian has the fol- In the next section, we show that the model in Eq. lowing form in the momentum space after the Fourier (1) provides an excellent description of the triplon exci- transformation. tationsseeninthescatteringexperimentsonBiCu PO . 2 6 In the following, we will describe the triplon excitations H =E + 1(cid:88)Λ†(k)M(k)Λ(k), (6) in the bond operator theory developed for the rung-VBS quad o 2 k 4 where Eo is a function of s¯and µ, and 35  t1(k)  30 H=0 t (k)  2  Λ(k)=t†1(−k) (7) ] 25 t†(−k) V 2 me 20 [ with q) 15 ( t (k) t† (−k) ωn 1x 1x 10 t1(k)= t1y(k) , t†1(−k)= t†1y(−k) , (8) t1z(k) t† (−k) 5 1z andsimilarlyfort2(k)andt†2(−k). Here,thesubscripts, 0 0 0.5 1 1.5 2 1 and 2, in the t-operators indicate the two dimers in q=(h=0 and 3, k, l=1) r.l.u. a unit cell (shown in Fig. 2). As mentioned earlier, the x,y,zdirectionsaretakenparalleltothecrystallographic FIG. 3. (Color online) Triplon energy dispersions, ω (q), in n a,b,c axes. The J, D, Γ, and H interactions in the thequadraticbondoperatortheory. Thedispersionsarecom- original spin Hamiltonian are transformed to the triplon paredwiththeneutronscatteringresultsinRef. [12](denoted hopping and pairing amplitudes contained in the 12×12 wthcolordots). Thespingappositions,±q(cid:63) =±(0,0.425,0), matrix M(k). The quadratic Hamiltonian has two no- are marked by the arrows. The gray region indicates the table features: (i) its dependences on J and J(cid:48) appear multi-triploncontinuumcomputedwithEq. (13). Errorbars 2 2 only through the sum J +J(cid:48), and (ii) the D interac- intheneutronscatteringdatarepresenttheenergywidths(or 2 2 1 decay rates) of the measured quasiparticle peaks. tionscanceleachotheratthequadraticlevelwithoutany contribution to H . These features will be discussed quad again later. The detailed form of H is provided in quad B. Triplon dispersions Appendix B. ThequadraticHamiltonianisdiagonalizedviatheBo- The triplon dispersions obtained in the quadratic the- goliubov transformation leading to the following form: ory can be fitted with the experimental results12 by 6 controlling the coupling constants {J,D,Γ}. Figure 3 (cid:88)(cid:88) Hquad =Egs+ ωn(k)γn†(k)γn(k). (9) shows the best fit (black lines) with the neutron scat- k n=1 tering data (color dots) along the momentum direction q=(h=0 and 3, k, l=1), yielding the following set of Here, γ†(k) is the Bogoliubov quasiparticle operator n the coupling constants. or the triplon with the excitation energy ω (k). The n constant term corresponds to the ground state energy, J =J =J(cid:48) =J =8 meV, 1 2 2 4 Egs =(cid:104)Hquad(cid:105). With this setting, the parameters s¯and J3 =0.2J1, µforthegroundstatearedeterminedbythesaddlepoint Γaa =−Γbb =0.039J , (12) 1 1 1 equations: Γab =Γba =0.135J , 1 1 1 (Da =0.6J , Db =0.45J ). ∂(cid:104)H (cid:105) ∂(cid:104)H (cid:105) 1 1 1 1 quad =0, quad =0. (10) ∂s¯ ∂µ The other coupling constants not shown here are set to zero since they are found to be irrelevant for describing In the ground state, we also compute the magnetization, essentialfeaturesoftheneutronscatteringdata. Asmen- M, under nonzero magnetic fields as follows. tioned earlier, the D couplings cancel each other in the 1 1 (cid:88) quadratic Hamiltonian H whereas the Γ couplings M= (cid:104)gµ S (cid:105), (11) quad 1 N B i survive at the quadratic level. Hence, we control the Γ1 i couplings instead of the D couplings. The D couplings 1 1 with N being the number of the spin moments. in Eq. (12) are the values obtained from the Γ cou- 1 Before discussing the results of the quadratic theory, plings and the relationship Eq. (2). As will be shown wecommentonourconventionsabouttheBrillouinzone later, the D interactions appear in the cubic terms of 1 and momentum vectors. For direct comparisons of our H . Their effects on the triplon excitations will be cubic theoretical computations with experimental results, we investigatedinalaterpartofthispaper. Furtherdiscus- use the extended zone scheme for the Brillouin zone and sions on the parameter regime of Eq. (12) are provided denote momentum vectors in reciprocal lattice unit as in Appendix C. q = (h,k,l), which means q = hG +kG +lG with The quadratic theory with the coupling constants in a b c G being the reciprocal lattice vectors corresponding Eq. (12) yields six nondegenerate, triplon dispersions a,b,c tothelatticevectorsR alongthea,b,caxes,repsec- (see Fig. 3). It is due to the fact that the anisotropic a,b,c tively. andsymmetriccouplings,Γ ,completelybreaktheSO(3) 1 5 spinrotationsymmetryexistingatthelevelofHeisenberg (green and blue dots in Fig. 4). Indeed the magnetic model, and there are two dimers in a unit cell. The spin field dependence of the two higher energy modes follows gap (minimum excitation energy) occurs at the incom- c +c H+c H2 behaviour. This is due to the fact that 0 1 2 mensurate momentum positions: ±q(cid:63) = ±(0,0.425,0) the triplon modes do not possess well-defined spin quan- and their equivalent momenta translated by reciprocal tum number (S·Hˆ, spin component along the field di- lattice vectors (denoted with arrows in the figure). We rection) as a result of magnetic anisotropies. Instead the find that the lowest three dispersions are in good agree- usual spin quantum numbers (S·Hˆ = +1,0,−1) in the mentwiththeneutronscatteringresultsaroundthespin spin isotropic case are mixed in the triplon modes. The gap. Among the other higher-energy three dispersions, non-linear magnetic field dependence can also be found only one of them is experimentally observed and qual- in other field directions. itatively consistent with the theoretical result (see the It can be seen that certain triplon modes have almost highest line of blue dots). constantenergyinmagneticfield. Forinstance,whenthe In regions away from the spin gap, however, the magnetic field is applied along the a axis, the energy of quadratic theory cannot fully explain the results from the lowest energy triplon mode is basically constant (red the experimental measurements. To be specific, inside dots in Fig. 4), implying that the spin character of this the gray region, the lowest dispersion bends downward mode is dominated by the quantum number S·aˆ = 0. (red), and the lowest two triplon modes decay at certain Similar behavior is also found in the second/third low- momenta(redandgreen). Thesefeaturesarebelievedto est mode under the magnetic field along b/c axis (see betheeffectsoftriploninteractionscomingfrom, for ex- the middle and right panels in Fig. 4). Accordingly, ample,thecubictermsgeneratedbytheD1 interactions. we can characterize three triplon modes approximately These effects will be discussed later. At the moment, we as the S·aˆ = 0, S·ˆb = 0, S·cˆ = 0 states from the focus on the low energy part around the spin gap (below lowest to the highest energy modes (the spin characters the gray region) and see if the quadratic theory provides can be directly identified by taking the spin projections a satisfactory description of the low energy spin dynam- of the triplon mode eigenvectors). In other words, three ics in BiCu2PO6. As we have already seen, the theory is triplon modes have their own approximate spin quanti- in good agreement with the neutron scattering results in zationaxes. Iftheappliedmagneticfieldisnotalongthe the low energy region. This fact supports the idea that quantization axis, triplon modes follow non-linear mag- the spin excitations observed in the neutron scattering neticfielddependenceandcanbecharacterizedbymixed aretriplonexcitations,confirmingtherung-VBSstatein spin states. BiCu PO . Below we provide more evidences for this 2 6 conclusion. D. Magnetization and critical field C. Magnetic field response We now discuss the magnetization process of the system. Recent high magnetic field experiments in Ref. Asanothermeantocheckthattheobservedexcitation [7] show that BiCu2PO6 undergoes a cascade of field- modesarethetriplonsintherung-VBSphase,theirfield induced phase transitions with anisotropic magnetic re- responses can be examined in the theory and compared sponses to different field directions. Upon increasing with the experimental data. Turning on the Zeeman in- thefield,themagnetizationmonotonicallyincreaseswith teraction in the quadratic bond operator Hamiltonian, different slopes depending on the field direction until wecomputethetriplonexcitationspectrumasafunction the system reaches the transition where the spin gap is of the magnetic field. The obtained spectra are plotted closed. The average slope or susceptibility, χavg = ∆∆MH, in Fig. 4 at the spin gap positions, ±q(cid:63). Existing ex- has the sequence of χc > χb > χa and concomi- avg avg avg perimental data12 is only for the fields along the a axis, tantly the critical field, H , has the opposite sequence of c which is also denoted in the figure with color dots. One Ha (= 23T) > Hb (= 21T) > Hc (= 20T), where the c c c can find that the theoretical result is consistent with the superscripts denote the applied field direction. experimental data both qualitatively and quantitatively The above experimental features can be explained in (leftpanel). Forexample,itwasobservedinexperiments the quadratic bond operator theory. Figure 5 shows the that the lowest mode (red) is almost not reacting to the magnetizationspredictedfromthetheory,whicharecon- magneticfieldwhereastheothertwohigherenergymodes sistent with the susceptibility sequence pattern observed aremovingdownward(green)andupward(blue),respec- intheexperiments. Ontheotherhand,thecriticalfields tively. These behaviours are well captured in the theo- canbereadfromtheprevioustriplonenergyplotsinFig. retical result. The observed magnetic field response also 4. Extrapolating the triplon spectra (gray dashed lines), tellsusthatthetriplonmodesdonotpossessanydegen- we find that Ha (= 27T) > Hb (= 20T) > Hc (= 19T) c c c eracy as predicted in the theory. in the quadratic theory. Although the value of Ha is a c Notice that, in general, the magnetic field depen- bit larger than the measured value, the right trend in dence of the triplon dispersion is not linear in magnetic the sequence of the critical fields is well captured in the field, especially for the two higher energy triplon modes theory. 6 7 7 7 6 6 6 V] 5 V] 5 V] 5 e e e m 4 m 4 m 4 ★) [ ★) [ ★) [ q 3 q 3 q 3 ± ± ± (n (n (n ω 2 ω 2 ω 2 1 1 1 0 0 0 0 5 10 15 20 25 30 0 5 10 15 20 0 5 10 15 20 H [T] // a H [T] // b H [T] // c FIG. 4. (Color online) Magnetic field dependence of the triplon excitations in the quadratic bond operator theory. The three plots show the triplon excitation energies at the spin gap positions, ω (±q(cid:63)), as a function of the magnetic field, H, for the n field directions along the a,b,c axes, respectively. In the case of the field along the a axis (the left), the theoretical results are compared with the neutron scattering data in Ref. [12] (colored dots). For an estimation of the critical field, H , in each plot c the theoretical results are extrapolated with a functional form (c +c H+c H2) as depicted with a dashed line.18 0 1 2 E. Importance of triplon interactions (red and green), and (iii) broadened energy width or increased decay rate in the third lowest dispersion (blue). Notice that the gray region in the figure denotes multi- In the previous sections, we have observed remarkable triploncontinuumwiththelowerboundarycomputedas consistency between the quadratic bond operator theory the minimum energy of two-triplon excitations: and experiment in the low energy regions. This fact im- pliesthatBiCu2PO6hastherung-VBSgroundstatewith E2(q)= min[ωm(q−k)+ωn(k)], (13) thetriplonsastheelementaryexcitations. Moreover,the k,m,n theorytellsusthatthecouplingsinEq. (12)areminimal where ω (k) is the single-triplon dispersion in the n interactions. quadratictheorywithm,nbeingbandindices. Thisfact On the other hand, the experimental data deviates proposesapicturethatthetriplonsarestronglyinteract- from the theoretical calculations in the gray region of ingwithinthetriploncontinuumregionsothattheylose Fig. 3 with the following features: (i) downward bend- their prominence as quasiparticle modes inside the re- ing of the lowest triplon dispersion (red), (ii) abrupt de- gion. Thus,weneedtoconsidertheinteractionsbetween cays of the lowest twotriplon modes at certain momenta triplons to capture triplon decay processes. IV. EFFECTS OF TRIPLON INTERACTIONS 6 H//a H//b 5 Now we consider the influence of triplon interactions H//c on the triplon dynamics in BiCu PO . Two main ef- ] 2 6 Cu 4 fects are expected from the triplon interactions: triplon /B energy renormalization and decay.19,20 These two effects µ 3 3 are closely related with the triplon phenomenology ob- -0 servedintheneutronscatteringexperiments(substantial [1 2 M energysplittingsaroundthespingapanddecayphenom- ena inside the triplon continuum). 1 The triplon interactions are taken into account by ex- 0 tendingthepreviousquadraticbondoperatortheory. For 0 2 4 6 8 10 12 14 the spin Hamiltonian [Eq. (1)] with the couplings in Eq. (12),wearrangethecorrespondingbondoperatorHamil- H [T] tonian in the following way: FIG. 5. (Color online) Magnetization, M, as a function of H=H [Γ ]+H [D ]+H [D ,Γ ], (14) quad 1 cubic 1 quartic 1 1 magneticfield,H,obtainedinthetheory. Themagnetization is computed with Eq. (11) for the magnetic fields along the wherewehavedenotedthedependenceontheanisotropic a,b,c axes. The figure shows the the magnetization compo- couplingsinthesquarebrackets. Wewilltakethemean- nents along the field directions. Other components perpen- fieldapproximationsforH andusethecubicinter- quartic dicular to the field are zero. actions, H , to describe the triplon decay processes. cubic 7 In this approach, the triplon modes and their decays are nels. Details of the decoupling scheme are explained in identified via the peak structures in the spectral weight Appendix D. function of the triplon Green’s function. This approach reveals that the D couplings have noticeable contribu- 1 tionsontheenergysplittingsaroundthespingapaswell as the decay phenomena inside the triplon continuum. B. Triplon decay channels of H cubic Readersinterestedinthe resultsrather thanthecalcula- tional details are advised to directly move to Sec. IVF. The cubic terms in H provide decay and fusion cubic processes of the t-bosons (t† (cid:10) t†t†). The processes are inducedbythecubictermsfromtheD andJ couplings A. Mean-field approximations for H 1 3 quartic whereastheothercouplings,J ,J (=J(cid:48)),J ,andΓ ,do 1 2 2 4 1 not have such cubic terms because of symmetry reasons. The quartic terms in H provide two-body scat- quartic More details about the existence of the cubic terms and terings of the t-bosons. We include the two-body inter- the associated symmetry are provided in Appendix E. actions in the mean-field approximations, leading to the following mean-field Hamiltonian: Wewillinvestigatetheeffectsofthecubicinteractions viatheGreen’sfunctionapproach.19–21 Inthisapproach, H +H →H . (15) thetriplonsdeterminedfromthemean-fieldHamitonian, quad quartic MF H , are taken as the bare particles. We express the MF For the mean-field decouplings, we employ particle- interaction part H in terms of the bare triplons (γ) cubic particle (Q = (cid:104)tt(cid:105)) and particle-hole (P = (cid:104)tt†(cid:105)) chan- from H , which leads to the following form: MF 6 (cid:88) (cid:88) 1 H = Y (k ,p ;q )γ†(k )γ†(p )γ(q )+H.c. cubic 2! 2 l m n l m n l,m,n=1k+p−q=0 6 (cid:88) (cid:88) 1 + Y (k ,p ,q )γ†(k )γ†(p )γ†(q )+H.c. (16) 3! 3 l m n l m n l,m,n=1k+p+q=0 Here, we use the shorthand notations: k = (k,l), and where the average (cid:104)···(cid:105) is taken for the ground state l so forth for p , q . In addition to the decay and fusion of H with the time-ordering operator T. Here, we m n terms for the triplons (Y and Y∗), we have the source set (cid:126) = 1, and Γ(k,t) = eiHtΓ(k)e−iHt with Γ(k) = 2 2 andsinktermsfortheγ-triplons(Y3andY3∗)intheabove [γ1(k),··· ,γ6(k)|γ1†(−k),··· ,γ6†(−k)]T. The Green’s expression. The vertex functions, Y2(kl,pm;qn) and function is written as the following 12×12 matrix: Y (k ,p ,q ), are functions of the singlet-condensation 3 l m n (s¯), the Bogoliubov transformation matrix of H , the (cid:20)G11(k,t) G12(k,t)(cid:21) MF G(k,t)= , (19) coupling constants (J’s, D , Γ ), and the lattice vec- G21(k,t) G22(k,t) 1 1 tors (R ). The vertex function Y (k ,p ,q ) is totally b,c 3 l m n symmetric whereas the other Y (k ,p ;q ) is symmet- with the normal Green’s function parts G11,22 and the 2 l m n ric only for the first two arguments: Y (k ,p ;q ) = anomalousfunctionpartsG12,21 asthe6×6submatrices. 2 l m n Y (p ,k ;q ). The triplon self energy has the same matrix structure: 2 m l n (cid:20)Σ11(k,t) Σ12(k,t)(cid:21) Σ(k,t)= . (20) C. Green’s function formalism Σ21(k,t) Σ22(k,t) The Green’s function and the self energy are related It is convenient to recast the total Hamiltonian as fol- by the Dyson equation: lows. G(k,ω)=G (k,ω)+G (k,ω)Σ(k,ω)G(k,ω). (21) 0 0 H=H +V, (17) 0 In the momentum-frequency space, the bare Green’s withthekineticpartH0 =HMF (fromHquad+Hquartic) function G0(k,ω) for the Hamiltonian H0 is a diagonal and the interaction part V =Hcubic. We now define the matrix with the matrix elements triplon Green’s function: G(k,t)=−i(cid:104)T (cid:2)Γ(k,t)Γ†(k,0)(cid:3)(cid:105), (18) (cid:2)G101(k,ω)(cid:3)mn = ω−ωδm(kn)+iη, (22) m 8 the spectral function is defined as 1 A(k,ω >0)=− Im{tr[G(k,ω)]}. (25) π As a simple example, one can check that A (k,ω) = 0 (cid:80)6 δ[ω − ω (k)] for the noninteracting Hamiltonian n=1 n H . With the triplon interaction V, the delta-function 0 peaks representing the bare triplon modes are modified into finite-size peaks having renormalized energy and FIG. 6. One-loop diagrams for the self energy Σ11(k). nonzerowidth. Ifthepeakisstillwell-definedwithlarge height and narrow width, the associated triplon mode survives as a quasiparticle with a finite life time. In the where ω (k) (m = 1,··· ,6) are the triplon eigenval- m next section, we analyze the triplon modes and their de- ues obtained from H and η = 0+. The other piece MF cay processes by looking at the peak structures of the of the diagonal elements is given by the relationship spectral function. G22(k,ω) = G11(−k,−ω). Using the Dyson equation, 0 0 we can express the full Green’s function in terms of the bare Green’s function and the self energy. For example, F. Results one can obtain (G11)-1 =(G11)-1−Σ11−Σ12[(G22)-1−Σ22]-1Σ21, Here we present the results of the interacting triplon 0 0 theory. We control the coupling constants such that G21 =[(G22)-1−Σ22]-1Σ21G11, thequasiparticlepeakstructuresinthespectralfunction 0 (23) A(q,ω) match the neutron scattering data. The result- and similarly for G22,12. ingrenormalizedcouplingconstantsarethenobtainedas follows. J =J =J(cid:48) =J =10 meV, D. One-loop self energy 1 2 2 4 J =0.2J , 3 1 (26) Da =Db =0.3J , Let us consider one-loop self-energy correction. One- 1 1 1 Γab =Γba =0.045J loop diagrams are drawn in Fig. 6 for the part Σ11(k) 1 1 1 of the self energy. The diagrams are translated into the Compared to the previous estimation in Eq. (12), following equation: the Heisenberg couplings have been increased and the [Σ11(k)] anisotropic couplings, D1 and Γ1, have been somewhat mn reduced with the inclusion of the triplon interactions. = 1(cid:88)(cid:88)Y2∗(pa,(k-p)b;km)Y2(pa,(k-p)b;kn) Such reductions of the anisotropic couplings reflect the 2 p a,b k0−ωb(−k+p)−ωa(p)+iη fact that the D1 couplings have sizeable energy renormalization effects contributing to the energy splittings + 1(cid:88)(cid:88)Y3(pa,(-k-p)b,km)Y3∗(pa,(-k-p)b,kn). (24) around the spin gap. 2 −k −ω (k+p)−ω (p)+iη 0 b a The spectral function, A(q,ω), calculated for the cou- p a,b pling constants [Eq. (26)], is plotted in Fig. 7 with the Here, we are using the abbreviated notation k = (k0,k) twodifferentstyles, (a)lineand(b)colormap. InFig. 7 with k0 and k being the frequency and the momentum, (a), we find that the triplon phenomenology mentioned respectively, and m,n,a,b (= 1,··· ,6) are triplon band before is well captured in the spectral function. First of indices. Besides the two diagrams in the figure, there all, the spin gap is found at k = 0.575 with the energy, is one more possible diagram having the vertices Y2 and ∆th (cid:39) 1.1meV, which is comparable to the measured Y∗. However,itvanisheswithnocontributiontotheself value,∆ =1.67meV. Aroundthespingap,threesepa- 2 ex energy. Other parts of the self energy can be calculated rated triplon modes are found with the energy splittings in similar ways. of ∆ω ∼2 meV consistently with the experimental re- th sult, ∆ω = 2 ∼ 3meV. The energy splittings around ex the spin gap are mainly the energy renormalization ef- E. Spectral function fectsoftheD couplingsamongtheanisotropiccouplings 1 in Eq. (26). One can notice this from the comparable TheGreen’sfunctioncanbecalculatedbyinsertingthe sizes of ∆ω and Da,b(= 3meV), and the small ratio ex 1 one-loop self energy into the Dyson equation [Eq. (21) Γ /D =0.15. 1 1 or (23)]. We will extract information about the triplon Next, moving inside the triplon continuum region modes by computing the spectral weight function of the (gray), thetriplonmodesundergosubstantialdecaypro- Green’sfunction. Forpositivefrequenciesofourinterest, cesses as shown in significantly broadened quasiparticle 9 (a) (b) 20 20 101 18 18 16 16 0 10 14 14 ] 12 ] 12 V V me 10 me 10 10-1 [ 8 [ 8 ω ω 6 6 -2 10 4 4 2 2 0 0 10-3 0.5 0.6 0.7 0.8 0.9 1.0 0.5 1.0 1.5 q=(h=0 and 3, k, l=1) r.l.u. q=(h=0 and 3, k, l=1) r.l.u. FIG. 7. (Color online) Spectral function, A(q,ω), in the interacting triplon theory. The spectral function is computed for the coupling constants in Eq. (26) and plotted in the two different styles, (a) with lines and (b) with a color map. It is compared with the neutron scattering data of Ref. [12] (denoted with color dots). The gray shaded region in (a) and the solid gray line in (b) denote the multi-triplon continuum computed with Eq. (13) and the bare triplon dispersions from H (= H ). The 0 MF dashedgraylinerepresentsthelowerboundaryofthecontinuumobtainedinthequadratictheory. Intheplot(b),thetriplon continuum is partially revealed by two-triplon states generated by the triplon fusion channels in V(=H ). The two-triplon cubic states are spread down to the solid gray line with weak intensities in the plot. peakstructures. Remarkably,thetriploncontinuumpre- spinstructurefactorhavebeenobservedinthecompound dicted in the bond operator theory matches well the re- SrCu (BO ) .22,23 2 3 2 gion where the decay processes are found in the neutron On the other hand, the bended dispersion could be an scattering experiments. This can be seen in comparison effectofthelevelrepulsionbythecontinuumasproposed with the experimental data in Fig. 7 (b). The second inRef. [12]. Inourtheory,theeffectsofthecouplingsbe- (green)andthird(blue)lowesttriplonmodesobservedin tweenthesingle-andtwo-triplonexcitationswereimple- the experiments disappear or lose its prominence above mented via the one-loop self energy in the single-triplon the solid gray line that represents the lower boundary of propagator;thelevelrepulsioneffectappearstoberather the triplon continuum calculated with HMF. For com- small at the one-loop level. To investigate the full level parison,wealsoshowthelowerboundaryobtainedinthe repulsioneffect, itmaybenecessarytotakeintoaccount quadratictheory(dashedgrayline)inthesameplot. The higher order corrections or consider the single- and two- decay processes are mainly caused by the D1 couplings triplon sectors on equal footing. with a minor contribution from the J coupling (see Sec. 3 In both cases, these considerations would provide IVB). It is confirmed by conducting computations with a unique opportunity to investigate interesting multi- one of the two couplings being turned off. particledynamicscausedbyanisotropicspininteractions. We leave this problem as an interesting subject of future study. G. Discussions Now we comment onthe couplingconstants estimated from the triplon theory. Due to many independent ex- Amongseveralexperimentalfeaturesinthetriplondy- change couplings in the spin model, we considered the namics, the downward bending of the lowest triplon en- simplest parameter regime [Eq. (26)] that captures es- ergy dispersion [red dots in Fig. 7 (b)] is not clearly sential features of the neutron scattering data (see Ap- explained in the current theory. For the bended part pendix C). Nonetheless, we find that the overall energy of the dispersion, we can think of two possibilities: (i) scale of the parameter regime is consistent with a previ- two-triplon bound state and (ii) the level repulsion by ous estimation in Ref. [5]. Table II shows our results in the triplon-continuum. The former idea naturally arises comparison with theirs. Despite various differences be- by noting that the bended part has a similar shape to tween the two works, both results have an almost same that of the gray line. Although the bended part disap- averagevalueoftheHeisenbergcouplingsperaunitcell: pears at certain momentum, it could be a matrix ele- (4J +2J +2J(cid:48) +2J +2J )/12 (cid:39) 8.6 meV. A major 1 2 2 3 4 ment effect. Such two-triplon state contributions to the difference is that in Ref. [5] J(cid:48) is evaluated to be much 2 10 Ref. [5] This work helpful discussions. This work was supported by the Magnetic susceptibility Neutron scattering NSERC of Canada and the Canadian Institute for Ad- (poly-crystal) (single-crystal) vanced Research. This research was also supported in Exact diagonalization Bond operator theory part by Perimeter Institute for Theoretical Physics. Re- (Heisenberg model) (generic model in Eq. (1)) search at Perimeter Institute is supported by the Gov- J1=140 K(cid:39)12 meV J1 =J2 =J2(cid:48) =J4 =10 meV ernmentofCanadathroughIndustryCanadaandbythe J2 =J1 J3 =0.2J1 ProvinceofOntariothroughtheministryofResearchand J2(cid:48) =0.5J1 D1a =D1b =0.3J1 Innovation. Computationswereperformedonthegpcsu- J4 =0.75J1 Γa1b =Γb1a =0.045J1 percomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation un- dertheauspicesofComputeCanada;theGovernmentof TABLEII.ComparisonofourworkwithRef. [5]. Thesecond and third rows indicate experimental results and theoretical Ontario; Ontario Research Fund - Research Excellence; approaches employed in the two works. and the University of Toronto. smallerthanJ whileinourtriplontheorysuchadissim- 2 ilaritybetweenJ andJ(cid:48) isnotcrucialfordescribingthe Appendix A: Anisotropic spin interactions 2 2 neutron data. We hope this point is clarified in future studies. ThelowenergyspinHamiltonianinEq. (1)canbede- The DM interactions responsible for magnetic rived from the microscopic Hubbard model consisting of anisotropies in BiCu PO were estimated in our study. theelectronhoppinghandtheon-siteCoulombrepulsion 2 6 Although various experimental results could be well de- U: scribedandunderstoodbyourtheory,themagnitudesof (cid:88) (cid:88) theestimatedDMinteractionsarequitelarge(D1a,b/J1 = H = c†iαhij,αβcjβ +U ni↑ni↓. (A1) 0.3) compared to the values usually found in copper ox- i>j i ides. This may change once higher order contributions are taken into account in the triplon self-energy (beyond The electron hopping amplitude generically consists of the one-loop level). For a more accurate estimation of the spin-independent (t) and spin-dependent (v) parts: the coupling constants, the microscopic spin model may be directly studied with numerical techniques. Further hij,αβ =tijδαβ +ivij ·σαβ, (A2) experimentalinformationsuchaselectronspinresonance (ESR)measurementswillbealsohelpfulfordetermining where σ are the Pauli matrices and α,β ∈{↑,↓} are the the DM interactions more precisely.24–26 spinindices. Thespin-dependenthoppingshavetheirori- gin in the atomic spin-orbit coupling. The corresponding hopping amplitude v is a three-component pseudo- ij V. SUMMARY AND OUTLOOK vector satisfying v = −v . Taking the large Coulomb ji ij interaction limit (U/h(cid:29)1) with the half electron filling In this paper, we provided theoretical analysis of the (one electron per site) and developing a degenerate per- rung-VBS ground state in BiCu PO by constructing a turbationtheory,onecanobtainthespinHamiltonianin 2 6 minimal spin Hamiltonian and developing a comprehen- Eq. (1)asalowenergyeffectivemodel.16,27Thecoupling sive theory of triplon dynamics. In comparison with the constants are defined in the following way. neutron scattering experiment data, it is shown that the anisotropic spin interactions (D1 and Γ1) are crucial to Jij =4t2ij/U, explain the unusual quantum numbers carried by the triplons and the decay processes of the triplons to the D =8t v /U, (A3) ij ij ij multi-triplon continuum. Ourresultswouldprovideessentialinformationforvar- Γµν =(cid:0)8vµvν −4δµνv2(cid:1)/U. ij ij ij ij ious ongoing studies of BiCu PO . In particular, the re- 2 6 centhigh-fieldexperimentsfoundaseriesoffield-induced It is clear from the expressions that the Dzyaloshinskii- quantum phase transitions.7 Nature of the field-induced Moriya and anisotropic & symmetric interactions share phaseshasnotbeenclearlyunderstood,inpartduetothe thesameorigininthespin-orbitcoupling. Theirrelation- lack of the spin Hamiltonian incorporating anisotropic ship in Eq. (2) is obtained from the above microscopic spin interactions. We believe that the spin Hamiltonian expressions for the coupling constants. presented here, together with the information about the spin content of the triplons, is a good starting point for the study of the field-induced phases. Appendix B: Quadratic Hamiltonian ACKNOWLEDGMENTS The quadratic Hamiltonian takes the following form. WethankKempPlumbandYoung-JuneKimforshar- ingtheirneutronscatteringdataonBiCu PO andmany 2 6