Dynamical structure factor of magnetic Bloch oscillations at finite temperatures

# Dynamical structure factor of magnetic Bloch oscillations at finite temperatures

Olav F. Syljuåsen Department of Physics, University of Oslo, P. O. Box 1048 Blindern, N-0316 Oslo, Norway
July 15, 2019
###### Abstract

Domain-walls in one-dimensional Ising ferromagnets can undergo Bloch oscillations when subjected to a skew magnetic field. Such oscillations imply finite temperature non-dispersive low-frequency peaks in the dynamical structure factor which can be probed in neutron scattering. We study in detail the spectral weight of these peaks. Using an analytical approach based on an approximate treatment of a gas of spin-cluster excitations we give an explicit expression for the momentum- and temperature-dependence of the spectral weights. Generally the spectral weights increase with temperature and approaches the same order of magnitude as the spin-wave spectral weights at high temperatures. We compare the analytical expression to numerical exact diagonalizations and find that it can, without any adjustable parameters, account for the and momentum-transfer dependence of the numerically obtained spectral weights in the parameter regime where the ratio of magnetic fields and the temperature is . We also carry out numerical calculations pertinent to the material , and find qualitatively similar results.

## I Introduction

Bloch Oscillations (BO)Bloch (); Zener () is a pure quantum mechanical phenomenon where a particle in a periodic potential oscillates when acted upon by a constant external force. Although BO are quite delicate, they have been observed in a range of systems in the last two decadesBOsemiconductors (); BOcondensates (); BOoptics (); BOultrasonics (). In addition BO have also been predictedKyriakidisLoss () to exist in one dimensional magnetic systems, so called magnetic Bloch oscillations (MBO). However, such MBO have so far not been observed experimentally.

One way to search for MBO is to search for its energy spectrum which is the Wannier-Zeeman (WZ) ladderWannier (); Shiba () of equally spaced energy levels. Experimentally one can hope to probe the WZ ladder directly using neutron scattering. At low temperatures neutrons will induce transitions from the ground state to the WZ levels. However, the transition matrix elements of these processes are very small, thus the low temperature inelastic neutron scattering signal will be very weakShinkevich1 (); Rutkevich (). In contrast, at temperatures comparable to the dominant magnetic scale the WZ levels will be thermally populated and neutrons can induce scattering between them. These transitions show up as peaks at low frequenciesKyriakidisLoss () in the neutron scattering dynamical structure factor, and are the topic of this paper.

This paper focuses especially on the spectral weights of the finite temperature low-frequency peaks and how they depend on momentum transfer and temperature. These quantitative details are important in order to distinguish the peaks associated with BO from generic low-frequency peaks in other systems with Zeeman spin-split spectra.

In order to calculate the spectral weights we use both an analytical and a numerical approach. The analytical approach is based on a gas of spin-cluster excitations with collisions treated in an approximate way. This goes beyond the oneKyriakidisLoss () and twoShinkevich1 () domain-wall approximations where collisions are neglected completely. The neglect of collisions can be expected to work well at low temperatures, but is more dubious at high temperatures where one expects frequent collisions between domain-walls. Also the one- and two domain-wall approximations restrict the state space severely to states that constitute the WZ levels. Since the existence of the low-frequency peaks is caused by the thermal population of the WZ levels it is important to verify that thermal population of other excluded states will not destroy them. For these reasons we validate our analytical result by carrying out numerical exact diagonalizations keeping all the states.

While the analytical result presented in this paper is strictly valid for the Ising model in a skew magnetic field we expect it to hold qualitatively for a broader class of Hamiltonians. We substantiate that at the end of the paper by performing numerical exact diagonalizations of a Hamiltonian that also includes extra terms relevant for describing the material .

The outline of the paper is the following. In chapter 2 we present the Hamiltonian and explain qualitatively why it should describe MBO. Chapters 3 and 4 explain the analytical calculation of the dynamical structure factor. The numerical exact diagonalization result is presented in chapter 5, while chapter 6 is devoted to a detailed comparison of the analytical and numerical results. In chapter 7 we present numerical results for the dynamical structure factor for parameter values relevant to the material , and chapter 8 concludes the paper.

## Ii Hamiltonian

The prediction that domain-walls in one-dimensional ferromagnets undergo BO in a magnetic field is rather general and holds for a range of HamiltoniansKyriakidisLoss (). In order to be definite we will focus on the one-dimensional spin-1/2 ferromagnetic Ising Hamiltonian with nearest-neighbor couplings in a uniform skew magnetic field

 H=−∑i(JzSziSzi+1+hxSxi+hzSzi). (1)

We will consider the case where the Ising coupling is the dominant one, and . Whenever numerical energies, frequencies or magnetic fields are quoted, they are in units of . The lattice spacing and are set to unity.

The ground state of the Hamiltonian resembles closely the ferromagnetic all spins up state as quantum fluctuations are small. Excitations can be classified according to how many domain-walls they contain. A domain-wall is an antiferromagnetic arrangement of two neighboring spins and costs an energy . It is these domain-walls that undergo BO in a skew magnetic field. To see why, consider first the transverse field . It couples to the spin operator that flips a spin. When is large the main effect of the transverse field is to flip a spin adjacent to a domain-wall which causes it to move, and will broaden the domain-wall excitations into a band. Thus the role of is to mimic a periodic potential. The longitudinal magnetic field, , plays the role of the external force as it will pull the domain-wall in the direction of the down-spins. The resulting BO of a single domain-wall in this system was recently studied in real-time numerical simulationsCai ().

However, in contrast to other systems with BO, a single domain-wall cannot exist in isolation in the magnetic system. The competition between the longitudinal magnetic field and the Ising coupling means that the system cannot have arbitrarily long spin-down domains. Instead new domain-walls will be created, and the excitations will be pairs of domain-walls, i.e. domain-wall/(anti)-domain-wall bound states. Nevertheless, the magnetic field still causes the domain-walls to oscillate, thus MBO in domain-wall ferromagnets resemble breather oscillations.

## Iii Spin-cluster excitations

At the excitations are domains of consecutive down-spins surrounded by up-spins. For this is slightly altered, and the excitations will be superpositions of such domains with different numbers of down-spinsFogedby78 (). These excitations have been termed multiple-magnon bound statesTorranceTinkham (), but we will call them spin-cluster excitations for simplicity. A spin-cluster state can be described mathematically by labeling the state with consecutive down-spins starting at site as . Defining Fourier-transformed states as , and writing energy eigenstates as superpositions of these, , the energy eigenvalue equation in the sector with two domain-walls takes the form

 hxcos(p/2)[ψp,n(l+1)+ψp,n(l−1)]=(Jz+lhz−Ep,n)ψp,n(l),l=1,2,… (2)

where we have defined , and neglected all terms that change the number of domain-walls. This is expected to be a good approximation for large .

Eq. (2) resembles the recursion relation for Bessel functions and has the solution

 Ep,n=Jz+hzνp,n (3)

and

 ψp,n(l)=CnJl−νp,n(zp) (4)

where is a normalization constant, labels the energy levels and . We have introduced the dimensionless ratio . To get a normalized wave function in the infinite chain limit we have set the coefficient in front of the increasing Neumann function solution to zero. The numbers are obtained by the other boundary condition, Eq. (2) for Fogedby78 ():

 J−νp,n(zp)=0. (5)

For , will depend on for the lowest values of , thus the lowest lying spin-cluster excitations will be dispersiveShinkevich1 (). For , . Thus for , is independent of and the energy levels are equally spaced with a separation . This is the WZL and we expect finite temperature neutron scattering to reveal transitions within the ladder with characteristic frequencies , where . The wave function describing the spin-cluster excitations in this WZL is a superposition of domains of different lengths , each with weight for .

## Iv Dynamical structure factor

We will now use the spin-cluster excitations to calculate the spectral weight of the low frequency peaks analytically. The transverse dynamical structure factor is

 Sxx(q,ω)=1Z∑i∫∞−∞dt2πeiωt⟨i|e−βHSx−q(t)Sxq(0)|i⟩ (6)

where is the partition function and is the inverse temperature. The spin operators in momentum space are and have time dependence . The sum is taken over all energy eigenstates of the system. The time integral can be performed by inserting an extra sum over energy eigenstates

 Sxx(q,ω)=1Z∑ijδ(ω−(Ej−Ei))e−βEi|⟨j|Sxq|i⟩|2. (7)

We will assume that the energy eigenstates can be written as states of almost independent spin-clusters. Thus we specify an energy eigenstate with spin-clusters as where labels the internal quantum number of spin-cluster , and its momentum. we assume that this state has energy . Thus we neglect any interaction energies between different spin-clusters. Such many-spin-cluster states are orthogonal and complete for . For finite we expect orthogonality and completeness to hold approximately when the density of spin-clusters is not too high.

The internal quantum number is an integer, so the dynamical structure factor will have peaks at low energies when for . To get the spectral weight of the peak at we integrate over frequencies in the vicinity of

 Sxxk(q)=1Z∑ije−βEi|⟨j|Sxq|i⟩|2|Ej−Ei=khz (8)

where the sum over states is restricted to states which energies differ from by .

When the operator acts on the state it gives a sum over states which each differ from state by having a single spin overturned. Many of these states do not contribute in Eq. (8) as they have an almost zero overlap with states having energy that differ from by . This is in particular true when creates(destroys) two domain-walls. This increases(decreases) the energy by and will therefore not contribute to low energy peaks at . The only states that contribute are those which is the result of flipping a spin so as to displace an already existing domain-wall. When the spin-clusters are not too densely packed we can treat the spin-clusters independently. In that case the matrix element reduces to the matrix element between two single spin-cluster states as the other spin-clusters are unchanged. However, at the temperatures considered here, one also needs to go beyond the independent approximation and take into account that the movement of a domain-wall can cause it to collide with the neighboring spin-cluster. The result of such a collision is the disappearance of two domain-walls and a merger of two domains. The energy of such a final state differs by the initial state by roughly , and should therefore not be counted as a contribution to the low frequency peaks. We thus need to explicitly exclude final states where a domain-wall is moved such as to merge two domains. We achieve that by including projections which meaning is that the domain in state must be such that its left domain-wall can be moved one lattice spacing to the left without touching another domain-wall for there to be a contribution. We can therefore write the matrix element that contribute to peaks at as

 ∑j|⟨j|Sxq|i⟩|2=di∑r=1|AnrprP←rL+A∗nrprP→rR+Bnrpr+B∗nrpr|2 (9)

where is the number of domains in state and

 Anrpr =12∞∑l=1ψ∗pr+q,nr+k(l+1)ψpr,nr(l)ei(ql−pr)/2, (10) Bnrpr =12∞∑l=2ψ∗pr+q,nr+k(l−1)ψpr,nr(l)ei(ql+pr)/2, (11)

are matrix elements between two single-cluster states, one with quantum numbers and its left(right) domain-wall displaced one lattice spacing to the left compared to the other state which have quantum numbers . The right movement is described by their complex conjugate.

By multiplying out the right hand side of Eq. (9) we get

 Sxxk(q) =1Z∑ie−βEidi∑r=1{12(P←rL+P→rR) ×[|Anrpr|2+(Anrpr+A∗nrpr)(Bnrpr+B∗nrpr)] +P←rLP→rR[A2nrpr+A∗2nrpr]+(Bnrpr+B∗nrpr)212} (12)

where we have set and .

The first term, with , is proportional to the probability to find a spin-cluster with quantum numbers whose left domain-wall can be expanded to the left without touching the neighboring domain (periodic boundary conditions is assumed). This probability does not depend on the spin-cluster index , thus the sum over can be taken into account by replacing the probability with the number of occurrences. We approximate this quantity with the number of occurrences of consecutive down spins surrounded by one up spin to the right and two up spins to the left, , where is the projection operator onto up(down) spins at site i. The brackets denote both quantal and thermodynamical average. The two up-spins to the left allow the displacement of the left domain-wall one lattice spacing to the left without touching a neighboring domain. For this approximation is exact. For small it is valid approximately as the wavefunction is dominated by the term with . Similarly the terms with is approximated by , by and the terms without projections by . This results in

 Sxxk(q) =∑n∑p{12[2|Anp|2+(Anp+A∗np)(Bnp+B∗np)]N↑↑n↑ +[A2np+A∗2np]N↑↑n↑↑+(Bnp+B∗np)2N↑n↑12} (13)

where we have used inversion symmetry .

At high temperatures, , we can approximate the s by setting and use transfer matricesSuzuki (). In the thermodynamic limit the number of down-spins surrounded by up-spins to the left and up-spins to the right is

 N↑…↑↓…↓↑…↑=NP↑Km−1+WKn−1−WKl−1+ (14)

where is the probability of finding a certain spin to be up. is the domain wall factor and is associated with finding two adjacent spins that points up(+) (down(-)). The constant

 λ=cosh(βhz/2)+[sinh2(βhz/2)+e−βJz]1/2 (15)

is the largest eigenvalue of the transfer matrix. Similarly one can find the number of domain-walls that separates exactly up-spins and down-spins. It is in the thermodynamic limit

 N↑…↑↓…↓=NγKn−1+WKm−1− (16)

where .

The expressions for and , Eqs. (10) and (11), involve an infinite sum over . In order to evaluate this we add and subtract terms with negative values of , see Appendix. Substituting the resulting expressions into Eq. (13) gives

 Sxxk(q) =12N↑↑↓J1−k(wq)(J1−k(wq)−2J−1−k(wq)) +12N↑↓J2−1−k(wq) −12δk,0J21(wq)∑n=1(N↑n↑+2N↑↑n↑+N↑↑n↑↑)cos(qn) +Δ(q) (17)

where and is a correction due to extending the sums to negative infinity

 Δ(q) =δk,0x232(3N↑↓↑+8N↑↑↓↑+N↑↑↓↑↑) −δk,0w2q8(2N↑↓↑+6N↑↑↓↑+N↑↑↓↑↑) −δk,1x232(7N↑↑↓↑+2N↑↑↓↑↑) +δk,1w2q16(5N↑↑↓↑+2N↑↑↓↑↑)+O(x4). (18)

## V Numerical results

In order to validate Eq.(17) we will calculate the dynamical structure factor numerically. To calculate at finite temperatures we have diagonalized the Hamiltonian (1) of a chain of spins with periodic boundary conditions. Fig. 1 shows the result as a function of frequency for a fixed value of momentum along the chain at different temperatures. The magnetic fields are and .

At low temperatures the only visible peak occurs at . It can be attributed to the flipping of a single spin in an otherwise all spin-up background. The intensity of this spin-wave excitation diminishes as the temperature is increased, and new low-energy peaks appear at , where is an integer. This is the finite temperature signature of inter-level WZL transitions. At high temperatures the intensity of the peak becomes the same order of magnitude as the spin-wave peak.

Some of the main features of the finite temperature dynamical structure factor in Fig. 1 can be explained qualitatively using what one expects in the limit of vanishing when the Hamiltonian (1) becomes purely classical. In that limit the peak at is associated with flipping a spin at a domain-wall. This process does not create new domain-walls, and the energy cost is therefore only the cost of one extra spin opposing the magnetic field. The peak at can be attributed to flipping a spin inside a domain of spins that are pointing opposite to the magnetic field. This creates two new domain-walls and makes one less spin pointing opposite to the field. These two processes can only occur at finite temperatures as they require the presence of domain-walls. Other features, such as the existence of peaks at cannot be explained from considerations at .

When one expects generally that the location of the peaks in will acquire a dependence on the momentum . This is true for the spin-wave peaks, but not for the low-energy peaks associated with the inter-level WZL transitions. They remain at frequencies for all . However, according to Eq. (17) their spectral weights should depend non-trivially on the momentum . To estimate numerically the spectral weight of each low-energy peak at , , we integrate the dynamical structure factor numerically over a small frequency interval chosen big enough to capture the full peak. We choose , and label the resulting integrated peak intensity as in the analytical calculation. A plot of these as functions of is shown in the inset of Fig. 1 for . The peak has the biggest spectral weight which decreases with increasing , while the spectral weights of the peaks increase with , except for the Bragg contribution at which is due to a finite ground state magnetization in the -direction and is unrelated to BO. Our analysis is restricted to the lowest values of because at higher values of there is an overlap with the peak at . To keep these separate we require .

## Vi Comparisons

To compare the analytical result Eq. (17) with our numerical results we consider each low-frequency peak (k) of the dynamical structure factor separately.

### vi.1 k=1

The expression for simplifies considerably at ,

 Sxxk=1(0) =12N↑↑↓−x232(7N↑↑↓↑+2N↑↑↓↑↑). (19)

The upper panel of Fig. 2 shows how the expression Eq. (19) compares to the exact diagonalization results for . It is rather accurate for , but underestimates the numerical result for .

According to the analytical result the dynamical structure factor for depends on through the variable . Expanding for small we find

 Sxxk=1(q) =Sxx1(0)−w2q16(126N↑↑↓−5N↑↑↓↑−2N↑↑↓↑↑) +O(x4) (20)

The lower panel of Fig. 2 shows the exact diagonalization results for plotted as a function of for different values of and at the same temperature. The data points for fall on the same curve, and an asymptotic linear dependence on , shown as the dotted line, is seen for small . In the inset we have plotted the negative slope of this dotted line, the point at . Repeating the scaling plot for several temperatures we get the remaining points. The solid curve shows the analytical result for the slopes, Eq. (20). For the analytical results agree well with the numerical results, while at higher the analytical result overestimates the slope. The reason for this is probably that interactions between spin-clusters becomes dominant at high temperatures and densities and must be treated more accurately.

### vi.2 k=2

The part of the peak is small and of , whereas depends on through the variable . The leading behavior for small is

 Sxxk=2(q)=Sxx2(0)+N↑↑↓w2q8+O(x4). (21)

In Fig. 3 we have plotted from exact diagonalization vs. for several values of and magnetic fields at a temperature . The points fall on the same curve demonstrating scaling with . The dotted line shows the asymptotic linear behavior of the curve. It was obtained by lumping all the data points into one set, and fitting them to a quadratic function. In the inset we show the slope of this curve together with slopes of scaling functions obtained at other temperatures and compare them with Eqs. (16) and (21). As in the case the comparison is good up to .

### vi.3 k=0

For the zero energy peak, , we need to evaluate the sum over in Eq. (17). Using the classical expressions, Eq. (14), we get

 ∑n=1 (N↑n↑+2N↑↑n↑+N↑↑n↑↑)cos(qn) =(N↑↓↑+2N↑↑↓↑+N↑↑↓↑↑)cos(q)−K−1+K2−−2K−cos(q). (22)

Expanding to order we find

 Sxx0(q) =x232(3N↑↓↑+8N↑↑↓↑+N↑↑↓↑↑) +w2q8(N↑↓+3N↑↑↓) −w2q8(N↑↓↑+2N↑↑↓↑+N↑↑↓↑↑) ×cos(q)−K−1+K2−−2K−cos(q) −w2q8(2N↑↓↑+6N↑↑↓↑+N↑↑↓↑↑). (23)

This expression is a smooth function of . Thus there are zero energy excitations also at a finite . In addition to this expression, at , there is a term coming from the finite ground state magnetization in the -direction, i.e. from intra ground state transitions which cause the discontinuous behavior at seen in the inset of Fig. 1 for . This additional term is proportional to the ground state transverse magnetization squared. Although it is the dominant term, we will ignore it here as it is not related to BO and only occurs at and . Fig. 4 shows a comparison of the numerical results with the expression Eq. (23) for several values of the temperature and ratios of magnetic fields. The agreement is reasonable for low -values and temperatures .

## Vii Cobalt niobate

The material is an experimental realization of a system of weakly coupled one-dimensional Ising-like ferromagnets. It has been studied with neutron scattering and clear evidence of dispersive spin-cluster states have been observed in zero applied magnetic fieldsColdea (). The fact that the excitations are dispersive at zero fields means that they cannot be interpreted as belonging to the WZL. Furthermore, it means that the Hamiltonian of must contain extra terms in addition to the dominant Ising couplingColdea (). A more detailed microscopic model which includes both internal magnetic fields and extra spin-spin couplings has been proposed for this materialKjall (). In particular, the large ratio of the internal fields causes the spin-clusters to be dispersive over a large region of momentum space at low energies, and the WZL-ladder is present only near the zone boundaryRutkevich () which in the neutron scattering experiments is dominated by the “kinetic bound state” caused by the extra spin-spin couplings.

Still, one should expect to see the finite temperature signatures of MBO in if the magnetic fields can be arranged so that . To investigate this closer we have repeated the exact diagonalization of the skew field Ising model including also additional termsKjall (). The extra terms we have added to the Hamiltonian (1) are where and . The magnetic field components, including both internal and external magnetic field contributions, have been set to , , and the temperature is . Fig. 5 shows the resulting dynamical structure factor as a function of frequency for eight different momenta plotted on top of each other.

At high frequencies, , there are dispersive peaks. For low frequencies one can clearly see the two non-dispersive peaks at , associated with transitions between WZL states even in the presence of the additional couplings. The spectral weights of these, and for a smaller peak at , are shown in the inset as a function of momentum. The general trend is the same as without the extra couplings, except that the intensity of the peak shows a larger dependence on and the peaks have in general a larger intensity.

## Viii Conclusions

We have investigated the finite temperature spectral signal of MBO in one-dimensional Ising models in a skew magnetic field using both analytical calculations and numerical exact diagonalizations. The main neutron scattering signature of MBO are finite temperature low-frequency peaks at , where the peak generally has the biggest spectral weight. While the peak frequency does not depend on momentum, its spectral weight does. We have shown that the spectral weight of the peaks can be calculated analytically when and using an approach that treats the excitations of the chain as a gas of spin-cluster excitations.

We have also investigated numerically how these results change when adding other couplings relevant for the material , and conclude that the extra couplings present will not destroy the MBO signatures.

This research was supported by a grant NFR-213606 from the Research Council of Norway, and the numerical calculations were performed on the Abel computer cluster with a Notur grant nn4563k.

## Appendix A Bessel sums

The expressions for the matrix elements needed in the calculation of the dynamical structure factor involve an infinite sum over . Substituting in the Bessel function forms for , Eq. (4), for , and omitting the normalization constant which is very close to unity we get

 Anp =12∞∑l=1Jl+1−n−k(zp+q)Jl−n(zp)ei(ql−p)/2, (24) Bnp =12∞∑l=2Jl−1−n−k(zp+q)Jl−n(zp)ei(ql+p)/2. (25)

In order to evaluate these sums we add and subtract terms with negative values of . For we get

 Anp =12∞∑l=−∞Jl+1−n−k(zp+q)Jl−n(zp)ei(ql−p)/2+ΔAnp =12J1−k(wq)ei(qn+π(1−k)−kp)/2+ΔAnp (26)

where and the sum has been evaluated using Graf’s summation formulaWatson (), and is correcting for the added terms with :

 ΔAnp =−120∑l=−∞Jl+1−n−k(zp+q)Jl−n(zp)ei(ql−p)/2 ≈δn,1(−1)kzkp+qzp2k+2k!e−ip/2 (27)

where we have replaced the sum by its lowest order term in which is . In a similar way

 Bnp =−12J−1−k(wq)ei(qn+π(1−k)−kp)/2+ΔBnp, (28) ΔBnp ≈δn,1(−1)kz1+kp+q2k+2(k+1)!ei(q+p)/2. (29)

## References

• (1) F. Bloch, Z. Phys, 52, 555 (1928).
• (2) C. Zener, Proc. R. Soc. Lond. A 145, 523 (1934).
• (3) E. E. Mendez, F. Agullo-Rueda, and J. M. Hong, Phys. Rev. Lett. 60, 2426 (1988).
• (4) M. Ben Dahan, E. Peik, J. Reichel, Y. Castin, and C. Salomon, Phys. Rev. Lett. 76, 4508 (1996).
• (5) R. Sapienza, P. Costantino, D. Wiersma, M. Ghulinyan, C. J. Oton, and L. Pavesi, Phys. Rev. Lett. 91, 263902 (2003).
• (6) H. Sanchis-Alepuz, Y. A. Kosevich, and J. Sanchez-Dehesa, Phys. Rev. Lett. 98, 134301 (2007).
• (7) J. Kyriakidis and D. Loss, Phys. Rev. B 58, 5568 (1998).
• (8) G. Wannier, Phys. Rev. 117, 432 (1960).
• (9) H. Shiba, Prog. Theo. Phys. 64, 466 (1980).
• (10) S. Shinkevich and O.F. Syljuåsen, Phys. Rev. B 85, 104408 (2012).
• (11) S. B. Rutkevich, J. Stat. Mech. (2010) P07015.
• (12) Z. Cai, L. Wang, X. C. Xie, U. Schollwöck, X. R. Wang, M. Di Ventra, and Y. Wang, Phys. Rev. B 83, 155119 (2011).
• (13) H. C. Fogedby, J. Phys. C: Solid State Phys., 11, 2801 (1978).
• (14) J. B. Torrance Jr. and M. Tinkham, Phys. Rev. 187, 587 (1969), Phys. Rev. 187, 595 (1969).
• (15) M. Suzuki, Phys. Lett. A 301, 398 (2002).
• (16) R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Science 327, 177 (2010).
• (17) Jonas A. Kjäll, Frank Pollmann, and Joel E. Moore, Phys. Rev. B 83, 020407(R) (2011).
• (18) G. N. Watson, A Treatise on the Theory of Bessel Functions, chapter XI, Cambridge Mathematical Library edition, Cambridge, 1995.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters