text:
# Transport and the Order Parameter of Superconducting Sr2RuO4

\[

## Abstract

Recent experiments make it appear more likely that the order parameter of the unconventional superconductor Sr2RuO4 has a spin-triplet $`f`$-wave symmetry. We study ultrasonic absorption and thermal conductivity of superconducting Sr2RuO4 and fit to the recent data for various $`f`$-wave candidates. It is shown that only $`f_{x^2y^2}`$-wave symmetry can account qualitatively for the transport data.

PACS numbers: 74.70.Pq, 74.20.Rp, 74.25.Ld

\]

The high-$`T_c`$-analog compound Sr2RuO4 poses a fundamental challenge to the study of unconventional superconductivity . Clean samples become superconducting at about $`T_c=1.5K`$ and the thermodynamics show that it is unconventional. The precise nature of its superconducting order parameter remains elusive, however. The simplest spin triplet case is the time-reversal-symmetry-breaking odd-parity $`p`$-wave state which has an order parameter $`𝐝=\widehat{𝐳}(k_x+ik_y)`$ . Spin-triplet character is consistent with the Knight shift  and neutron scattering  measurements which indicate that the spin susceptibility shows no obvious suppression below $`T_c`$. Moreover, $`\mu `$SR experiment  confirms that the order parameter does indeed break time-reversal symmetry. The central feature of this $`p`$-wave order parameter is that it exhibits finite gap all over the Fermi surface (FS). However, more recent thermodynamic and transport measurements at low temperatures have found powers laws in the temperature dependences. This includes specific heat $`C_VT^2`$ , penetration depth $`\mathrm{\Delta }\lambda T^2`$ , Ru spin-lattice relaxation rate $`1/T_1T^3`$ , ultrasonic attenuation $`\alpha T^3`$ , and thermal conductivity $`\kappa T^2`$ . These powers laws in $`T`$ are consistent with line nodes in the order parameter.

The combination of all experiments therefore favors the spin-triplet odd-parity f-wave symmetries which feature both time-reversal symmetry breaking and line nodes. Maki and Yang  have proposed an $`f`$-wave model with a superconducting order parameter $`𝐝(T)=\mathrm{\Delta }(T)\widehat{𝐳}k_z(k_x+ik_y)^2`$ (the so-called “3D” model) for Sr2RuO4. This 3D gap has line nodes where the basal plane ($`k_z=0`$) intersects the FS. It is similar to one of the leading candidates for the order parameter of UPt3 (known as the $`E_{2u}`$ model in that hexagonal system ). Graf and Balatsky  have proposed another candidate $`f`$-wave order parameter $`𝐝=\mathrm{\Delta }(T)\widehat{𝐳}(k_x+ik_y)k_xk_y`$, which we will call the $`f_{xy}`$ model. This order parameter has line nodes perpendicular to the basal plane at where the planes $`k_x=0`$ and $`k_y=0`$ intersect the FS. This nodal structure is analogous to the $`p`$-wave model $`𝐝=\mathrm{\Delta }(T)\widehat{𝐳}(\mathrm{sin}k_x+i\mathrm{sin}k_y)`$ proposed by Miyake and Narikiyo . It is suggested that $`f_{xy}`$-wave order parameter could arise from a pairing interaction mediated by antiferromagnetic spin fluctuations with peak intensities at the incommensurate wavevectors $`𝐪=(\pm 2/3\pi ,\pm 2/3\pi )`$. The latter has been confirmed by inelastic neutron scattering . Finally there is the $`f_{x^2y^2}`$-wave order parameter $`𝐝=\mathrm{\Delta }(T)\widehat{𝐳}(k_x+ik_y)(k_x^2k_y^2)`$ with line nodes where the planes $`k_x=\pm k_y`$ intersect the FS. This has been suggested  partly by analogy to high $`T_c`$ cuprates where a similar $`x^2y^2`$ (but with spin-singlet) gap structure has been conclusively established. It is worth noting that the above three models all account reasonably well for the specific heat and superfluid density data .

The structure of all of these gaps in momentum space is rather complicated. Since all thermodynamic and transport quantities involve Fermi surface averages, it is not easy to see how they may be distinguished experimentally. In the study of heavy-fermion superconductors, it has been found that thermal conductivity and ultrasonic attenuation can be key tools in this situation. These quantities can be studied as a function of propagation direction and, in the case of the attenuation, also of polarization dependence. The key point is that this directional sensitivity can not only test for the presence of point or line nodes, but can also determine the positions of nodes and the orientations of nodal lines. These experiments  offer the best chance to distinguish between order parameters that have different nodal structures but share similar thermodynamic properties.

In this paper, we calculate the thermal conductivity and the transverse and longitudinal ultrasonic attenuation of Sr2RuO4. We consider the three major candidates of $`f`$-wave symmetry mentioned above, namely the 3D, $`f_{xy}`$, and $`f_{x^2y^2}`$ models. It will be shown that only the $`f_{x^2y^2}`$ model can account for the data even qualitatively. The in-plane $`k_x^2k_y^2`$ symmetry offers an intriguing similarity to the high-Tc order parameter. If this gap structure is indeed manifested in Sr2RuO4 it would suggest the presence of anti-ferromagnetic (AF) fluctuations in the ruthenate superconductors, perhaps peaked at or near $`𝐪=(\pi ,\pi )`$ (such as one finds in the cuprates). The transport data thus suggests that AF fluctuation may also exist in ruthenates.

Ultrasonic attenuation experiments in Sr2RuO4 have been done in the hydrodynamic regime where the sound frequency $`\omega `$ satisfies $`\omega <<1/\tau `$, where $`\tau `$ is the electronic relaxation time . In this case the Boltzmann equation is appropriate, and the attenuation $`\alpha _{ij}`$ for sound propagation along the direction $`\widehat{q}\widehat{i}`$ with polarization $`\widehat{\epsilon }\widehat{j}`$ is :

$`\alpha _{ij}(T)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑\omega \left[{\displaystyle \frac{F(\omega )}{\omega }}\right]\tau _s(\omega ,T)`$ (1)
$`\times `$ $`\mathrm{Re}\left({\displaystyle \frac{\sqrt{\omega ^2|𝐝_𝐤|^2}}{\omega }}[(\widehat{i}\widehat{𝐤})(\widehat{j}\widehat{𝐤}){\displaystyle \frac{1}{d}}\delta _{ij}]^2\right).`$ (2)
Here $`d`$ denotes the dimensionality, $`F`$ is the Fermi function, and the angle brackets indicate an FS average. Note the second term inside the square bracket vanishes for the transverse case. Similarly the thermal conductivity in the $`i`$th direction when a temperature gradient is imposed along the $`i`$th direction is computed using:

$`\kappa _i(T)`$ $``$ $`{\displaystyle \frac{1}{T}}{\displaystyle _0^{\mathrm{}}}𝑑\omega \omega ^2\left[{\displaystyle \frac{F(\omega )}{\omega }}\right]\tau _s(\omega ,T)`$ (3)
$`\times `$ $`\mathrm{Re}\left({\displaystyle \frac{\sqrt{\omega ^2|𝐝_𝐤|^2}}{\omega }}(\widehat{ı}\widehat{𝐤})^2\right).`$ (4)
The FS of Sr2RuO4 has been determined by de Haas-van Alphen experiments. It has three sheets denoted by $`\alpha `$, $`\beta `$, and $`\gamma `$ . These are in good agreement with LDA calculations . The bands are two-dimensional, a feature similar to cuprates. Among these three bands, the $`\gamma `$ band has the largest FS and has the highest density of states (DOS) at the Fermi energy. For simplicity, we use a single cylindrical FS for $`f_{xy}`$ and $`f_{x^2y^2}`$ models and a three-dimensional spherical FS for the 3D model. We shall comment on the effects of FS anisotropy below. We will assume that the momentum dependence of the gap is independent of temperature, which leaves only the overall magnitude $`\mathrm{\Delta }(T).`$ The temperature dependence of the gap in Sr2RuO4 has not been directly measured. However, one may deduce the power law at low temperatures from the pattern of nodes and near $`T_{c\text{ }}`$the power is given by mean field theory. For line nodes, a reasonable interpolation for the temperature dependence is $`\mathrm{\Delta }(T)=\mathrm{\Delta }(0)[1(T/T_c)^3]^{1/2}`$. This form has only one parameter $`\mathrm{\Delta }(0)`$, reducing the arbitrariness in the fitting procedure.

Inspection of Eqs. (2) and (4) shows that the only remaining unknown is the relaxation time $`\tau _s(T)`$. Our procedure will be to use the thermal conductivity to determine the relaxation time and the magnitude of the ground state gap $`\mathrm{\Delta }(0)`$, then use these to constrain the fits of the attenuation.

In Fig. 1, we first show the best fits against the thermal conductivity $`\kappa _x`$ along the  direction in the normal and superconducting states. The normal state is achieved experimentally by applying a sufficiently strong external magnetic field. The data are those of Tanatar et al. . In the normal state at $`0.5T_cTT_c`$, the scattering rate can be well described by an empirical form: $`1/\tau _n=\mathrm{\Gamma }_0(1+aT^2/T_c^2)`$ with $`\mathrm{\Gamma }_00.01\pi T_c`$ and $`a=0.15`$. That is, the inelastic scattering rate is observed to be roughly $`15\%`$ of the elastic part at $`T_c`$. For the effective scattering rate in the superconducting state, we take the same temperature dependence as in the normal state. This is clearly appropriate for $`TT_c`$. At lower temperatures $`TT_c`$, due to the opening of the gap, the inelastic rate could decrease faster than $`T^2`$. However, since the elastic scattering rate dominates over the inelastic part, the exact temperature dependence of inelastic rate is not so important. (This is one big difference between this system and the heavy fermion and high-$`T_c`$ materials.) Fitting to the superconducting thermal conductivity data then fixes the other parameter $`\mathrm{\Delta }(0)`$. It turns out that all three theoretical candidates can fit the data for $`\kappa _x`$ in the superconducting state. There is no noticeable difference in goodness-of-fit for the different models. The parameters determined by the fits are $`\mathrm{\Delta }(0)=1.65T_c`$ for $`f_{x^2y^2}`$ model, $`\mathrm{\Delta }(0)=3.3T_c`$ for $`f_{xy}`$ model, and $`\mathrm{\Delta }(0)=4.05T_c`$ for 3D model. The difference in size between these various values is not significant, since each gap form has a different relation between $`\mathrm{\Delta }(0)`$ and the root-mean-square FS average gap. It is the latter quantity that is more important in determining the overall scattering strength.

The parameters determined by the fit in Fig. 1 are now used to calculate the transverse and longitudinal ultrasonic attenuation without any adjustable parameters (Fig. 2) for all three models.

As seen in Fig. 2, the experiment exhibits the very clear feature that the transverse T100 $`\alpha `$ is systematically higher than the longitudinal L100 $`\alpha `$ when scaled to the normal state values at $`T_c`$. Note that the polarization vector for the transverse sound is also in the $`xy`$ plane. This sets a crucial constraint on theoretical models.

Of the three candidate gaps, only the $`f_{x^2y^2}`$-wave has this property. The reason for this can be seen semi-quantitatively as follows. Inspection of Eq. (2) shows that one big difference between the transverse and longitudinal cases is a factor of $`\widehat{k}_x^2\widehat{k}_y^2`$, and $`(\widehat{k}_x^21/d)^2`$, respectively. The other factors restrict the integral to the FS, and, as a function of temperature, weight the points on the FS according to how small the gap is. At the lowest temperatures, only the nodal regions contribute. Thus we may think of the integral as being over an effective FS located near the nodes, with a polarization-dependent weighting factor. The 3D gap has an effective FS that is a circle at the equator $`k_z=0.`$ For this geometrical object, we have $`\widehat{k}_x^2\widehat{k}_y^2=(\widehat{k}_x^21/d)^21/8`$ with $`d=2`$. The angle brackets indicate an average over the effective FS. Thus the average for longitudinal polarizations is about the same as that for transverse. For the $`f_{xy}`$ state, the effective FS is a set of two lines at $`k_x=0`$ and $`k_y=0`$, and we have $`\widehat{k}_x^2\widehat{k}_y^20`$, while $`(\widehat{k}_x^21/d)^2`$ is finite ($`d=2`$), favoring longitudinal polarization. In contrast, for the $`f_{x^2y^2}`$ state, the effective FS is at $`k_x=\pm k_y`$. The averages therefore lead to $`\widehat{k}_x^2\widehat{k}_y^2`$ finite, while $`(\widehat{k}_x^21/d)^20`$ ($`d=2`$ again), favoring transverse attenuation. These averages show that the $`f_{x^2y^2}`$ will naturally have more transverse attenuation than its competitors. Although we have plotted only one set of curves for each of the states, this is a general feature. The fact that the attenuation data favor $`f_{x^2y^2}`$ over other candidates was already noted by Lupien et al.  without calculation.

Since the data in Figs. 1 and 2 are taken on different samples ($`T_c=1.44`$K and 1.37K, respectively), one does not expect a quantitatively good fit in Fig. 2 using the same parameters in Fig. 1. At first glance, our results of $`\alpha `$ for $`f_{x^2y^2}`$ and $`f_{xy}`$ waves are qualitatively similar to Graf and Balatsky’s  when a resonant scattering limit ($`\delta =90^\mathrm{o}`$) is taken under the self-consistent $`T`$-matrix approximation. Lupien et al.  pointed out that while the sound attenuation data favors $`f_{x^2y^2}`$ wave, the low $`T`$ power law is much lower than the calculated result of Graf and Balatsky for the longitudinal case. In contrast, as shown in Fig. 2, our calculation based on a phenomenological input of scattering time gives satisfactory good fit for the longitudinal case including the low $`T`$ power law. It also leads to better fit for the transverse case especially at $`TT_c`$, compared to Graf and Balatsky’s almost linear in $`T`$ result. The still rather large discrepancy between the theoretical fits and the experimental data (especially for the T100 mode) could be due to the order parameter which is more involved than having the pure $`f_{x^2y^2}`$ symmetry. This remains to be studied in more details.

It is important to note one major difference between our phenomenological treatment and the self-consistent $`T`$-matrix approximation used by Graf and Balatsky. In the unitary limit, the usual $`T`$-matrix approximation for the impurity scattering will assume a strongly frequency dependent scattering rate $`1/\tau _{\mathrm{imp}}1/\omega `$ in the superconducting state . In contrast, we have adopted a constant elastic scattering rate (i.e., $`\mathrm{\Gamma }_0`$). Our phenomenological input relies truly on the good fit to the thermal conductivity data in the superconducting state (see Fig. 1).

In-plane thermal conductivity measured by Izawa et al.  with a rotating in-plane magnetic field has shown a four-fold pattern with a minimum when $`𝐇[110]`$. This angle dependence is consistent with the $`f_{x^2y^2}`$ and 3D models but rules out the $`f_{xy}`$ model . In terms of variation of the amplitude, Izawa et al. claimed that the state $`𝐝=\mathrm{\Delta }(T)\widehat{𝐳}(k_x+ik_y)(\mathrm{cos}ck_z+|b|)`$, originally proposed by Hasegawa et al. , is more consistent with the data based on the approach in Refs. . Here $`c`$ is the spacing between two adjacent Ru layers and $`|b|1`$. This model has a circular-like line node parallel to the basal plane with any $`k_z`$ satisfying $`1k_zc1`$. This gap is symmetrical for any given parallel $`𝐤`$ plane and will lead to a sound attenuation very similar to that of the 3D model above. Thus it also does not produce the proper polarization dependence.

Finally we comment on the effect of Fermi surface anisotropy. While in Sr2RuO4 the FS has multiple sheets and there seems to be a coupling between the $`\alpha `$ and $`\beta `$ bands, the analysis of transport data given here will not be significantly modified providing that the $`D_{4h}`$ symmetry of the crystal is preserved. Moreover, our result for the 3D model stands if the spherical FS is replaced by an ellipsoidal one that has circular projection parallel to the basal plane. As long as transport data is concerned, it is likely that interband effects are unimportant.

In conclusion, the transport data of Sr2RuO4 presently give strong preference to the $`f_{x^2y^2}`$-wave order parameters with perpendicular line nodes at $`k_x=\pm k_y`$. A sharper comparison of theory and experiment will be possible if the data for $`\alpha `$ and $`\kappa `$ is taken on the same sample. Furthermore, we suggest that the out-of-plane transverse sound attenuation measurement should be performed under the polarization $`\widehat{q}\widehat{x}`$ and $`\widehat{\epsilon }\widehat{z}`$. This will shed more light on the information of the location of line nodes and the overall nodal structure.

Financial support from NSC of Taiwan (Grant No. 89-2112-M-003-027) and the NSF under the Materials Theory program, Grant No. DMR-0081039 is acknowledged.


format: text

text:
# Quasiparticle excitation in and around the vortex core of underdoped YBa2Cu4O8 studied by site-selective NMR

\[

## Abstract

We report a site-selective 17O spin-lattice relaxation rate $`T_1^1`$ in the vortex state of underdoped YBa2Cu4O8. We found that $`T_1^1`$ at the planar sites exhibits an unusual nonmonotonic NMR frequency dependence. In the region well outside the vortex core, $`T_1^1`$ cannot be simply explained by the density of states of the Doppler-shifted quasiparticles in the $`d`$-wave superconductor. Based on $`T_1^1`$ in the vortex core region, we establish strong evidence that the local density of states within the vortex core is strongly reduced.

\]

The electronic structure in the vortex state is a fundamental problem in the physics of superconductors with unconventional symmetries. Recent experimental results have revealed that the quasiparticle (QP) excitations in the vortex state of high-$`T_c`$ cuprates are very different from those of conventional $`s`$-wave superconductors. Despite extensive studies on this subject, the question of the QP excitations both inside and outside the vortex core is still far from being settled. In fact, STM experiments have revealed many unexpected properties in the spectrum of the vortex core . For example, a zero energy peak in the tunnel conductance due to the extended zero energy QP states, which is expected in the conventional theory of the $`d`$-wave superconductor, is not observed at the center of the vortex core. Although the novel electronic structure of the vortex core has been discussed in terms of several intriguing models, it remains controversial . On the other hand, well outside the vortex core, the QPs are considered to be in extended states due to the presence of gap nodes. The dominant effect of these delocalized states in the magnetic field is the Doppler shift of the QP energy spectrum, which occurs due to the presence of a supercurrent flow around each vortex and generates a nonzero density of states (DOS) at the Fermi surface . The $`H`$\- and $`T`$-dependencies of the specific heat , thermal conductivity , and NMR Knight shift  at low temperatures have generally been accepted as strong evidence for this effect. However, rather surprisingly, neither the local DOS with a $`1/r`$ asymptotic tail expected from the Doppler shift nor the 4-fold symmetry of the extended QPs has been observed in the STM spectrum . Thus, a discrepancy still exists between the STM findings and the results of other experiments including thermodynamic measurement. The situation therefore strongly confronts us with the need for a powerful probe of the QP excitation with spatial resolution.

The measurements of the NMR spin-lattice relaxation rate $`T_1^1`$ provide microscopic information on the low energy QP excitations in the superconducting state. Recent experimental and theoretical studies have established that the frequency dependence of $`T_1^1`$ measurements serves as a site selective probe for the electronic spectrum in the vortex state of the high-$`T_c`$ cuprates . Although the site-selective NMR, compared to STM, does not possess an atomic resolution , it has two advantages. First, it is a bulk probe, free from surface effects. Second, as pointed out by several authors , the tunneling conductance obtained by STM is not simply related to the local DOS owing to the very anisotropic tunneling matrix element along the $`c`$-axis. The site selective NMR, on the other hand, especially at the oxygen sites in the plane where the antiferromagnetic spin fluctuation is cancelled, provides direct information on the local DOS . Very recently, the frequency dependence of $`T_1^1`$ at the 17O sites (denoted as $`{}_{}{}^{17}T_{1}^{1}`$) in the vortex state of optimally doped YBa2Cu3O7 has been reported . A key feature of the data is that $`{}_{}{}^{17}T_{1}^{1}`$ at the planar sites increases monotonically with NMR frequency well outside the vortex core, which is consistent with the delocalized QP DOS induced by the Doppler shift. In this Letter, we attempted to measure the site-selective NMR in underdoped YBa2Cu4O8 extending the measurement to the vortex core region which has not yet been reported. We have found an unusual NMR frequency dependence of $`{}_{}{}^{17}T_{1}^{1}`$ at the planar sites in the vortex lattice state, which is markedly different from that of YBa2Cu3O7. On the basis of these results, we discuss the QP structure in and around the vortex core of YBa2Cu4O8.

High quality underdoped YBa2Cu4O8 with $`T_c`$=81 K has been synthesized under a high-oxygen-pressure atmosphere. Previous NMR studies have shown that the normal state spin fluctuations at the planar Cu and O sites for YBa2Cu4O8 are similar in behavior to those for the 60 K phase of YBa2Cu3O6+x . The NMR data were obtained on an aligned powder of underdoped YBa2Cu4O8, which was isotopically enriched with 17O. The 17O spin echo signals were obtained by a conventional pulse spectrometer under a field cooling condition (FCC) in a constant field ($`H_0`$=9.4 T) applied parallel to the $`c`$-axis by using a highly homogeneous superconducting magnet. We stress here that the NMR measurements under the FCC are important because the Bean critical current associated with the sweeping magnetic field not only produces the field gradient in the crystal but also seriously influences $`{}_{}{}^{17}T_{1}^{1}`$ by producing an additional Doppler shift to the QPs.

The final NMR spectra were obtained by convolution of the respective FT-spectra of the spin echo signals measured by a step of 50 kHz. Given that the spectra of the central transition of the chain O(1), the planar O(2,3), and the apical O(4) overlap at low temperature, $`{}_{}{}^{17}T_{1}^{}`$ is determined by the recovery curves of the quadrupole split O(2,3) satellite. The inset of Fig. 1(a) shows the planar 17O (O(2) and O(3)) spectrum of the satellite transition at 100 K. The slightly-split between the satellite spectra of the oxygen sites of O(2) and O(3) is due to a small difference in electric field gradient. The procedure to determine $`T_1`$ is as follows. The spin echo intensities are measured as a function of time after saturation pulses. Then the FT spectrum of each echo is obtained. Each frequency component of the FT spectra is fitted individually by $`R(t)=(M(\mathrm{})M(t))/M(\mathrm{})=1/35\mathrm{exp}(2Wt)+3/56\mathrm{exp}(6Wt)+1/40\mathrm{exp}(12Wt)+25/56\mathrm{exp}(20Wt)+25/56\mathrm{exp}(30Wt)`$ for the (3/2$``$1/2) transition. Here, we calculated $`{}_{}{}^{17}T_{1}^{}`$ (=1/2$`W`$) for each frequency point at the 3 kHz interval with a gaussian weight function of $`\sigma `$=3 kHz.

The local field profile in the vortex state is obtained by approximating $`H_{loc}(𝒓)`$ with the London result,

$$H_{loc}(𝒓)=H_0\underset{G}{}\mathrm{exp}(i𝑮\mathbf{}𝒓)\frac{\mathrm{exp}(\xi _{ab}^2𝑮^2/2)}{1+𝑮^2\lambda _{ab}^2},$$
(1)
where $`𝑮`$ is a reciprocal vector of the vortex lattice, $`𝒓`$ is the distance from the center of the core, $`\xi _{ab}`$ is the in-plane coherence length, and $`\lambda _{ab}`$ is the in-plane penetration length. The thin solid lines in Figs. 1 (a) and (b) depict the histogram at particular local fields which are given by the local field distribution $`f(H_{loc})=_\mathrm{\Omega }\delta [H_{loc}(𝒓)H_{loc}]d^2𝒓`$, where $`\mathrm{\Omega }`$ is the magnetic unit cell. In the calculation we used $`\xi `$=15Å and $`\lambda _{ab}=1600`$ Å reported in Ref. and assumed the square vortex lattice. The magnetic field is the lowest at the center of the vortex square lattice (C-point in the inset of Fig. 1(b)), and it is the highest at the center of the vortex core (A-point). The intensity of the histogram shows a peak at the field corresponding to the saddle point of the vortices (B-point). Because the local NMR frequency is proportional to the local magnetic field, we are able to produce a spatial image of the low-lying QP excitation spectrum.

The solid lines in Figs. 1 (a) and (b) represent the NMR spectrum at the planar and the apical sites in the superconducting state, respectively. A clear asymmetric pattern, which originates from the local field distribution associated with the Abrikosov vortex lattice (the

Redfield’s pattern), is observed well below $`T_c`$ in both sites. The broken lines in Figs. 1 (a) and (b) depict the spectrum convoluted with the Gaussian broadening function, $`f(H_{loc})=e^{2H_{loc}^2/\sigma ^2}`$. (The broken line in Fig. 1(a) is obtained by the sum of the spectrum at the O(2) and (3) sites.) Using $`\sigma `$=24 kHz, the theoretical curves well reproduce the data in both planar and apical sites. The open circles in Figs. 1 (a) and (b) display the NMR frequency dependence of $`{}_{}{}^{17}T_{1}^{1}`$ at the planar (denoted as $`{}_{}{}^{17}T_{1,pl}^{1}`$) and at the apical ($`{}_{}{}^{17}T_{1,ap}^{1}`$) sites, respectively. A clear nonmonotonic frequency dependence is observed at the planar sites. With increasing frequency, $`{}_{}{}^{17}T_{1,pl}^{1}`$ first decreases, shows a minimum near the B-point, then increases. Finally, at high frequency when approaching

the vortex core regime (A-point), $`{}_{}{}^{17}T_{1,pl}^{1}`$ again decreases. This frequency dependence is in contrast to the result of Ref., in which $`{}_{}{}^{17}T_{1,pl}^{1}`$ shows monotonic increases with frequency. In contrast to $`{}_{}{}^{17}T_{1,pl}^{1}`$, $`{}_{}{}^{17}T_{1,ap}^{1}`$ shows little frequency dependence.

Figure 2 shows the frequency dependence of $`{}_{}{}^{17}T_{1,pl}^{1}`$ at various temperatures. The inset shows the half width at half intensity of the spectrum of 17O(2) as a function of temperature. One can clearly see a slight increase of the line width below $`T_c`$ and the large broadening due to the vortex lattice below $`T_{melt}45`$ K. $`{}_{}{}^{17}T_{1,pl}^{1}`$ becomes nearly independent of frequency at $`45`$ K, above which the vortex lattice melts into disordered vortex liquid. The frequency dependences of $`{}_{}{}^{17}T_{1,pl}^{1}`$ are remarkable below $`45`$ K. This result reinforces the fact that the observed frequency dependence of $`{}_{}{}^{17}T_{1,pl}^{1}`$ at lower temperatures originates from the field distribution associated with the vortex lattice. Figure 3 displays the $`T`$-dependencies of $`{}_{}{}^{17}T_{1,pl}^{1}`$ and $`{}_{}{}^{17}T_{1,ap}^{1}`$. For the planar sites, we show the largest and the smallest values of $`{}_{}{}^{17}T_{1,pl}^{1}`$ for corresponding temperatures. The $`T`$-dependence of the

smallest values is very similar to results reported previously.  Because $`T_1`$ is determined from averaged recovery curves over all frequency components by a conventional method, $`T_1`$ is attributed to the value from the most dominant components of the signals, namely, in the superconducting region far from vortex cores (around the B point). With decreasing $`T`$, $`{}_{}{}^{17}T_{1,pl}^{1}`$ decreases faster than $`T`$ above $`T_c`$ due to the pseudogap, and does not show an anomaly at $`T_c`$ . On the other hand, $`1/^{17}T_{1,ap}`$ is nearly proportional to $`T`$ except at $``$45 K, where a peak in $`{}_{}{}^{17}T_{1,pl}^{1}`$ associated with the vortex lattice melting is clearly observed.

We first discuss the relaxation mechanism. It should be noted that the spin diffusion effect from the core region is unimportant in the present case because of the large broadening of the spectra below $`T_M`$ . Generally, $`{}_{}{}^{17}T_{1}^{1}`$ is composed of a contribution from vortex vibration $`{}_{}{}^{17}T_{1,VV}^{1}`$ and a contribution from the QP scattering $`{}_{}{}^{17}T_{1,QP}^{1}`$ :

$$\frac{1}{{}_{}{}^{17}T_{1}^{}}=\frac{1}{{}_{}{}^{17}T_{1,VV}^{}}+\frac{1}{{}_{}{}^{17}T_{1,QP}^{}}.$$
(2)
The vortex vibration influences $`{}_{}{}^{17}T_{1}^{}`$ through the transverse fluctuation of the local magnetic field. It has been suggested that the vortex vibration plays an important role in $`T_1`$ for very anisotropic Tl2Ba2CuO6 . We consider that in YBa2Cu4O8, the QP scattering predominates $`{}_{}{}^{17}T_{1}^{}`$ at the planar sites, while the vortex vibration seemingly influences $`{}_{}{}^{17}T_{1}^{}`$ at the apical sites, based on the following explanation. First, as shown in Figs. 1 (a) and (b), $`{}_{}{}^{17}T_{1,ap}^{1}`$ shows much smaller frequency dependence than $`{}_{}{}^{17}T_{1,pl}^{1}`$. Second, if the contribution from the vortex vibration dominates, then $`{}_{}{}^{17}T_{1,ap}^{1}`$ should be the same as $`{}_{}{}^{17}T_{1,pl}^{1}`$. However, as shown in Fig. 3, $`{}_{}{}^{17}T_{1,ap}^{1}`$ are smaller than $`{}_{}{}^{17}T_{1,pl}^{1}`$ by a factor of $`5`$. Third, the vortex lattice melting seriously influences $`{}_{}{}^{17}T_{1,ap}^{1}`$, while its influence on $`{}_{}{}^{17}T_{1,pl}^{1}`$ is fairly small.

Having established that the $`{}_{}{}^{17}T_{1,pl}^{1}`$ is dominated by the QP excitation, the next question is the frequency dependence. We first discuss $`{}_{}{}^{17}T_{1}^{1}`$ in the low frequency region well outside the core. Within the framework of the semi-classical approximation, it has been pointed out that the local DOS $`N(E=0,r)`$ well outside the core is dominantly produced by the Doppler shift of the delocalized QP spectrum, $`E(𝒑)E(𝒑)𝒗_s𝒑`$, in the presence of the supercurrent flow with velocity $`𝒗_s=\mathrm{}/2mr`$. As a result, there is a finite $`N(E=0,r)(𝒗_s)`$, which decays as $`1/r`$ for a single vortex and varies slightly more modestly than $`1/r`$ in the vortex lattice state. According to Ref., this local DOS gives rise to a monotonous increase of $`{}_{}{}^{17}T_{1,pl}^{1}`$ with frequency.

In the intermediate frequency regime where $`{}_{}{}^{17}T_{1,pl}^{1}`$ increases, the Doppler shift appears to be important for the local DOS. However, the observed decrease of $`{}_{}{}^{17}T_{1,pl}^{1}`$ when going from C- to B-point is apparently inconsistent with the Doppler shift. There are two possible origins for this unusual behavior in $`{}_{}{}^{17}T_{1,pl}^{1}`$. One is the competition between the thermally excited QPs and the Doppler-shifted QPs in the presence of the strong antiferromagnetic fluctuation which is strongly suppressed by the Doppler shift, as pointed out in Ref.. The other is the vortex lattice effect, which is not taken into account in the semi-classical theory of Doppler shift. According to the numerical calculation, the localized QPs in the Andreev bound state around the vortex core extend to $`(\pm \pi ,\pm \pi )`$ directions, and can overlap the QPs of the neighboring vortices . This effect can increase the QP DOS at the C-point, depending on the vortex lattice structure. At present it is not clear which effect is relevant to $`{}_{}{}^{17}T_{1,pl}^{1}`$ at low frequency. The fact that the low frequency behavior of $`{}_{}{}^{17}T_{1,pl}^{1}`$ cannot be simply explained by the Doppler shifted QPs implies some key features of the QP excitation outside the vortex core.

Finally we discuss $`{}_{}{}^{17}T_{1,pl}^{1}`$ at high frequencies, at which $`{}_{}{}^{17}T_{1,pl}^{1}`$ begins to decrease after reaching maximum. It is apparent that this behavior is again incompatible with the Doppler shift. As pointed out in Ref., the semi-classical description of the QPs breaks down in the vicinity of the core region ($`r<2\xi `$) where the superconducting order parameter varies rapidly. Therefore it is natural to consider that the reduction of $`{}_{}{}^{17}T_{1,pl}^{1}`$ occurs in the vortex core regime. It should be noted that the antiferromagnetic fluctuation, if it exists, always increases $`{}_{}{}^{17}T_{1,pl}^{1}`$. These considerations lead us conclude that the decrease of $`{}_{}{}^{17}T_{1,pl}^{1}`$ is caused by the suppression of the QP DOS within the vortex core. This is the first experimental evidence of the reduction of the local DOS within the vortex core, which was probed by a bulk measurement. In the pure $`d`$-wave superconductors, the local DOS at the center of the vortex core should be enhanced, as it is in $`s`$-wave superconductors, owing to the interference of the Andreef scattered QPs at the both sides of the pair potential. It would follow that $`{}_{}{}^{17}T_{1,pl}^{1}`$ should be enhanced in the vortex core regime, which is apparently inconsistent with the present results. The present results are seemingly consistent with the STM measurements, which show that only one QP energy level is present within the core for the optimally doped YBa2Cu3O7 and Bi2Sr2CaCu2O8+δ . Several possible origins for the reduction of DOS have been proposed as follows: (1) different symmetry, such as $`d_{x^2y^2}+id_{xy}`$, is induced within the core ; (2) the Andreef bound state is split due to the large $`\mathrm{\Delta }^2/\epsilon _F`$ value where $`\mathrm{\Delta }`$ is the superconducting gap and $`\epsilon _F`$ is the Fermi energy; (3) the strong antiferromagnetic interaction changes the core state dramatically, as discussed by the calculation based on the $`tJ`$ model  and the SO(5) theory . To clarify the origin of the reduction of the local DOS within the vortex core, more detailed information of the QP structure, especially in the underdoped regime, is strongly required.

In summary, we measured the site-selective spin-lattice NMR relaxation rate in and around the vortex core of underdoped YBa2Cu4O8. Nonmonotonic NMR frequency dependence of $`T_{1,pl}^1`$ indicates that $`T_{1,pl}^1`$ well outside the vortex core regime cannot be simply explained by the Doppler-shifted delocalized quasiparticles of the $`d`$-wave superconductor and that the local DOS in the vortex core regime is strongly reduced.

We thank T. Isobe and S. Shamoto for assisting with the crystal preparation. We also acknowledge helpful discussions with M. Ichioka, Y. Kato, K. Machida, K. Maki, D.K. Morr, M. Ogata, M. Takigawa, A. Tanaka, and I. Vekhter. This work has been supported by a Grant-in-Aid from Ministry of Education, Culture, Sports, Science and Technology.

Note added in proof. After submission (at April 12, 2001) of this manuscript, we are aware of an article (in Nature 413 (October 4, 2001) 501) of similar NMR study on the spatial dependence of $`T_1^1`$ in near-optimally doped YBa2Cu3O7.


text:
# Asymmetric Anomalous Diffusion: an Efficient Way to Detect Memory in Time Series

## 1 Introduction

One of the basic expectations of the statistical analysis of time series of biological, sociological and financial interest is that the process under study might depart from fully random behavior. Following the prescriptions of the authors of Refs. we use the time series to create a sort of random walk trajectory whose statistical properties are then carefully studied to detect the deviations from ordinary Brownian motion. To do so a given site of the time series is interpreted as the “time” at which the walker makes a jump. The intensity of the jump is determined by the value of the time series at that given site, and the walker jumps either forward or backward according to whether the sign of the value is positive or negative. The most natural way to establish a comparison between this walker and the ordinary random walker of statistical mechanics would be to compare the two resulting diffusion processes to one another, but this would imply that infinitely identical walkers are available. Unfortunately a time series means only one trajectory, thereby generating the problem of how to create a number of “independent” trajectories large enough to produce a diffusion process. The authors of Refs. addressed this problem with an ingenious technique called *Detrended Fluctuation Analysis* (DFA) which allows us to bypass this difficulty. The whole pseudo-random trajectory is divided into non-overlapping intervals of length $`l`$, each of which can be imagined as being a trajectory running for a ”time” $`l`$. Thus we have a collection of a relatively large number of ”independent” trajectories that can be used to produce the wanted diffusion process. These trajectories might be affected by a local trend and this local trend might result in an artificial deviation from ordinary Brownian motion. For this reasons the authors of Refs. for any interval evaluate the local bias and then record only the deviation from the local trend. The square root of the resulting variance, denoted by them with the symbol $`F(l)`$, is expected to yield for very large $`l`$’s

$$F(l)l^\delta ,$$
(1)
where $`\delta `$ denotes the scaling parameter, a crucial property to assess the departure from ordinary Brownian motion. In fact, ordinary Brownian motion should yield $`\delta _B=0.5`$ and the deviation of the observed value of $`\delta `$ from $`\delta _B=0.5`$ is considered to be a measure of the deviation from a totally random behavior.

It has to be pointed out that scaling is a property implying that

$$p(x,l)=\frac{1}{l^\delta }F(\frac{x}{l^\delta })$$
(2)
and that in some special cases, discussed in this paper, the scaling value detected by means of Eq.(1) does not afford a satisfactory account of the property of Eq.(2). Furthermore, as we shall see, the scaling $`\delta `$, even if it reflects correctly the property of Eq.(2), might not disclose an exhaustive information on the dynamic source of the deviation from ordinary Brownian motion.

Here we illustrate a technique of analysis, called Diffusion Entropy (DE) method, originally applied to detect memory effects in time series of sociological interest. We shall show that this technique can, in part, solve the earlier mentioned problems. We use a sequence of the same kind as that adopted in earlier work to mimic the statistical properties of DNA sequences. This sequence consists of strings of extended length, filled with only $`+`$’s or $``$’s. We consider as significant event only the change of sign. Thus, this kind of sequence is changed into a new sequence realized by assigning the vanishing value to the sites with a sign identical to that of the following one. Only the sites with a sign opposite to that of the next site are assigned a non vanishing value, which is always the same positive quantity, $`W`$. Thus we obtain a sequence of rare events. When one of these events occurs, the random walker makes a jump of constant intensity, always in the same direction. We prove that the resulting diffusion process makes much more ostensible the correlated nature of the original time series.

We also show that this procedure makes faster the detection of the proper scaling parameter $`\delta `$, namely, a scaling value which is a fair reflection of the property of Eq.(2). As earlier mentioned, the methods determining scaling through variances might yield scaling parameters significantly departing from the ideal scaling value. This is the scaling that an artificial time series, whose asymptotic behavior is dictated by compelling theoretical arguments, is expected to reach in the long-time limit.

The outline of the paper is as follows. In Section 2 we describe the main properties of the DE method, and with the help of an artificial sequence we show that the observation of a symmetric diffusion process results in a slow convergence towards the emergence of the expected scaling parameter. Section 3 illustrates a theoretical approach to the scaling of the asymmetric diffusion process using two distinct procedures, the Continuous Time Random Walk (CTRW) and the Generalized Central Limit Theorem (GCLT), yielding the same result. This section also shows, with the help of the DE method applied to some artificial sequences, that the detection of the correct scaling is faster than in the case of symmetric diffusion. Finally, Section 4 shows our detection technique at work on real data. The illustrative example of Section 4 is based on a time series resulting from the record of human heart beats. A balance on the results obtained in this paper is made in Section 5.

## 2 Diffusion Entropy

The main idea of this approach to scaling is remarkably simple. Let us consider a sequence of $`M`$ numbers $`\xi _i`$ , with $`i=1,\mathrm{},M`$. The purpose of the DE algorithm is to establish the possible existence of a scaling, either normal or anomalous, in the most efficient way as possible, without altering the data with any form of detrending. Let us select first of all an integer number $`l`$, fitting the condition $`1lM`$. This integer number will be referred to by us as “time”. For any given time $`l`$ we can find $`Ml+1`$ sub-sequences defined by

$$\xi _i^{(s)}\xi _{i+s},s=0,\mathrm{},Ml.$$
(3)
For any of these sub-sequences we build up a diffusion trajectory, labelled with the index $`s`$, defined by the position

$$x^{(s)}(l)=\underset{i=1}{\overset{l}{}}\xi _i^{(s)}=\underset{i=1}{\overset{l}{}}\xi _{i+s}.$$
(4)
Let us imagine this position as referring to a Brownian particle that at regular intervals of time has been jumping forward or backward according to the prescription of the corresponding sub-sequence of Eq.(3). This means that the particle before reaching the position that it holds at time $`l`$ has been making $`l`$ jumps. The jump made at the $`i`$-th step has the intensity $`|\xi _i^{(s)}|`$ and is forward or backward according to whether the number $`\xi _i^{(s)}`$ is positive or negative.

We are now ready to evaluate the entropy of this diffusion process. To do that we have to partition the $`x`$-axis into cells of size $`ϵ(l)`$. When this partition is made we have to label the cells. We count how many particles are found in the same cell at a given time $`l`$. We denote this number by $`N_i(l)`$. Then we use this number to determine the probability that a particle can be found in the $`i`$-th cell at time $`l`$, $`p_i(l)`$, by means of

$$p_i(l)\frac{N_i(l)}{(Ml+1)}.$$
(5)
At this stage the entropy of the diffusion process at the time $`l`$ is determined and reads

$$S_d(l)=\underset{i}{}p_i(l)ln[p_i(l)].$$
(6)
The easiest way to proceed with the choice of the cell size, $`ϵ(l)`$, is to assume it independent of $`l`$ and determined by a suitable fraction of the square root of the variance of the fluctuation $`\xi (i)`$.

Before proceeding with the illustration of how the DE method works, it is worth making a comment on the way we use to define the trajectories. The method we are adopting is based on the idea of a moving window of size $`l`$ that makes the $`sth`$ trajectory closely correlated to the next, the $`(s+1)th`$ trajectory. The two trajectories have $`l1`$ values in common. The DFA, on the contrary, is based on non-overlapping windows, and consequently trajectories with different labels are totally independent the one from the other. A motivation for our choice is that we are having in mind a possible connection with the Kolmogorov Sinai (KS) entropy. The KS entropy of a symbolic sequence is evaluated by moving a window of size $`l`$ along the sequence. Any window position corresponds to a given combination of symbols, and, from the frequency of each combination, it is possible to derive the Shannon entropy $`S(l)`$. The KS entropy is given by the asymptotic limit $`lim_l\mathrm{}S(l)/l`$. We believe that the same sequence, analyzed with the DE method, at the large values of $`l`$ where $`S(l)/l`$ approaches the KS value, must yield a well defined scaling $`\delta `$. To realize this correspondence we carry out the determination of the DE using the same criterion of overlapping windows as that behind the KS entropy.

Details on how to deal with the transition from the short-time regime, sensitive to the discrete nature of the process under study, to the long-time limit where both space an time can be perceived as continuous, are given in Ref.. Here we make the simplifying assumption of considering so large times as to make the continuous assumption valid. In this case the trajectories, built up with the above illustrated procedure, correspond to the following equation of motion:

$$\frac{dx}{dt}=\xi (t),$$
(7)
where $`\xi (t)`$ denotes the value that the time series under study gets at the $`łth`$ site. This means that the function $`\xi (l)`$ is depicted as a function of $`t`$, thought of as a continuous time $`t=l`$. In this case the Shannon entropy reads

$$S(t)=_{\mathrm{}}^{\mathrm{}}𝑑xp(x,t)ln[p(x,t)].$$
(8)
We assume

$$p(x,t)=\frac{1}{t^{\delta (t)}}F\left(\frac{x}{t^{\delta (t)}}\right).$$
(9)
This is a generalization of the ordinary scaling assumption of Eq.(2), which, in fact, can be recovered by assuming $`\delta (t)`$ to be time independent so that we can denote it as $`\delta `$. For the sake of simplicity we make the simplifying assumption that $`F(y)`$ maintains its form, namely, that the statistics of the process are unchanged. This is questionable in the case where the scaling is a function of time. However, making this assumption we obtain the following results. Let us plug Eq.(9) into Eq.(8). After a simple algebra, we get:

$$S(\tau )=A+\delta (\tau )\tau ,$$
(10)
where

$$A_{\mathrm{}}^{\mathrm{}}𝑑yF(y)ln[F(y)]$$
(11)
and

$$\tau ln(t).$$
(12)
It is evident that this kind of technique to detect scaling does not imply any form of detrending, and this is one of the reasons why some attention should be devoted to it. An interesting way to check the efficiency of this technique is realized by the study of the artificial sequence of Ref. . This sequence is built up in such a way as to realize long sequences of either $`+`$’s or $``$’s. The probability of finding a sequence of only $`+`$’s or only $``$’s of length $`t`$ is given by

$$\psi (t)=(\mu 1)\frac{T^{\mu 1}}{(t+T)^\mu }.$$
(13)
Here we focus our attention on the condition $`\mu <3`$ and we raise the reader’s attention on the interval $`[2,3]`$. In fact, this kind of sequence is the same as that adopted in earlier work  for a dynamic derivation of Lévy diffusion, which shows up when the condition $`2<\mu <3`$ applies. It corresponds to a particle travelling with constant velocity throughout the whole time interval corresponding to either only +’s or only -’s, and changing direction with no rest, at the end of any string with the same symbols.

We will refer to this model as Symmetric Velocity Model (SVM). We know from the theory of Ref. that the scaling of the resulting diffusion process when $`2<\mu <3`$ is

$$\delta =\frac{1}{\mu 1}.$$
(14)
Note, however, that this diffusion process has a finite propagation front, with ballistic peaks showing up at both $`x=t`$ and $`x=t`$. The intensity of these peaks is proportional to the correlation function

$$\mathrm{\Phi }_\xi (t)=\left(\frac{T}{T+t}\right)^{\mu 2}.$$
(15)
As a consequence of this fact, the whole distribution does not have a single rescaling. In fact, the distribution enclosed between the two peaks rescales with $`\delta `$ of Eq.(14) while the peaks are associated to $`\delta =1`$. Furthermore, it is well known that the scaling of the second moment is given by

$$\delta =\frac{4\mu }{2}.$$
(16)
Thus, it is expected that the scaling detected by the DE method might not coincide with the prediction of Eq.(14) for the whole period of time corresponding to the presence of peaks of significant intensity. We think that the Lévy scaling of Eq.(14) will show up at long times, when the peak intensity is significantly reduced. This conjecture seems to be supported by the numerical results illustrated in Fig.1. We see in fact that the scaling predicted by Eq.(14) is reached after an extended transient, of the order of about 20,000 in the scale of Fig.1. This time interval is about 2000 larger than the value assigned to the parameter T, of Eq.(13), which is, in fact, in the case of Fig.1, T= 10.

In conclusion, this section proves that the DE method applied to the SVM yields, for the scaling parameter $`\delta `$, the correct value of Eq.(14), rather than the value that would be obtained measuring the variance of the diffusion process, Eq.(16). However, the time necessary to make this correct value emerge is very large. Furthermore, as proved by the theory of Ref., the adoption of SVM would make the scaling parameter $`\delta `$ insensitive to $`\mu `$ in the whole interval . This means that the adoption of the DE method would not allow us to distinguish a process with $`\mu `$ very close to 1 from one with $`\mu `$ very close to 2. In Section 3 we shall see how to overcome these limitations.

## 3 Asymmetric Diffusion

Let us consider again the kind of artificial sequence that we analyzed in Section 2. Let us change perspective. Let us explicity focus on the events corresponding to changing sign. For this reason we build up a new sequence where the sites of the laminar region are given the vanishing value. This means that the states of constant velocity are here perceived as states of rest. The random walker can make a jump only at the moment when the time series analyzed in Section II changes signs. If at the moment of changing sign the walker makes a jump, of the same intensity, forward or backward, according to whether the sign change is positive or negative, we shall refer to this process as Symmetric Jump Model (SJM). The scaling corresponding to the SJM has been studied by Shlesinger in the pioneer paper of Ref., which yields the following prescription when the condition $`1<\mu <2`$ applies,

$$\delta =\frac{1}{2}(\mu 1).$$
(17)
When $`\mu >2`$, the theory of Ref.  predicts

$$\delta =1/2.$$
(18)
We see that the whole region $`\mu >2`$, including that with $`\mu <3`$, would be indistinguishable from an ordinary Brownian diffusion. For this reason, in addition to focusing our attention on the events corresponding to the sign change, we decide to disregard whether the change of sign is from - to + or vice-versa. This generates a kind of asymmetric diffusion process that we term Asymmetric Jump Model (AJM). This section is devoted to the discussion of the AJM by means of two distinct theoretical approaches, both resulting in the same prediction, an interesting form of asymmetric Levy diffusion. Then we generate two artificial sequences corresponding to $`\mu >2`$ and $`\mu <2`$, respectively, and we analyze them by means of the DE method.

### 3.1 Continuous Time Random Walk

The first theoretical approach rests on the Continuous Time Random Walk (CTRW) of Montroll and Weiss . This means that the distance between two nearest neighbor events, $`\tau _it_{i+1}t_i`$, is considered to be a random number, described by the probability density $`\psi (\tau )`$. Note that we are here using the notation $`\tau `$ rather than $`t`$ to make clear to the reader the difference between the absolute time t and the time distance between two nearest neighbor events. As earlier said, we assume that the random walker always moves in the same direction. Moreover, we also set the condition that the steps have constant length. Thus, the probability that a step of length $`x`$ is taken is given by

$$\mathrm{\Pi }(x)=\delta (xW).$$
(19)
Applying to this case the CTRW formalism of Montroll and Weiss, we write for $`p(x,t)`$, the probability of the walker being at position $`x`$ at time $`t`$,

$$p(x,t)=F^1L^1\widehat{p}(k,s)F^1L^1[\frac{1\widehat{\psi }}{s}\frac{1}{1\lambda (k)\widehat{\psi }(s)}],$$
(20)
where $`\lambda (k)`$ and $`\widehat{\psi }(s)`$ are the Fourier and Laplace transform of $`\mathrm{\Pi }(x)`$ and $`\psi (t)`$, respectively. Note that

$$\lambda (k)=exp(ikW).$$
(21)
By Taylor series expansion of Eq.(20) we get:

$$\widehat{p}(k,s)=\frac{1\widehat{\psi }}{s}\underset{n=0}{\overset{\mathrm{}}{}}\widehat{\psi }^ne^{inkW}$$
(22)
By evaluating the inverse Fourier transform of Eq.(22) we arrive at:

$$\widehat{p}(x,s)=\frac{1\widehat{\psi }}{s}\underset{n=0}{\overset{\mathrm{}}{}}\widehat{\psi }^n\delta (xnW).$$
(23)
This means that for large distances we can adopt the following expression

$$p(x,s)=[\widehat{\psi }(s)]^{\frac{x}{W}}\frac{1\widehat{\psi }(s)}{s}.$$
(24)
We shed light on the meaning of Eq.(24) with the following remarks. Let us denote by $`\psi _x(t)`$ the probability of making $`x`$ jumps up to time $`t`$, with the last jump occurring at $`t`$. Let us denote by $`\mathrm{\Phi }_n(t)`$ the probability of making $`n`$ jumps at earlier times, $`t^{}<t`$. The two functions are connected by the relation

$$\mathrm{\Phi }_n(t)=_0^t\psi _n(t^{})\mathrm{\Psi }(t^{}t)𝑑t^{},$$
(25)
where

$$\mathrm{\Psi }(t)_t^{\mathrm{}}\psi (t^{})𝑑t^{}.$$
(26)
Using the properties $`\widehat{\psi }_n(s)=[\widehat{\psi }(s)]^n`$ and $`\widehat{\mathrm{\Psi }}(s)=(1\widehat{\psi }(s))/s`$, and Laplace transforming Eq.(25), we get

$$\widehat{\mathrm{\Phi }}_n(s)=\widehat{\psi }_n(s)\widehat{\mathrm{\Phi }}(s)=[\widehat{\psi }(s)]^n\frac{1\widehat{\psi }(s)}{s}.$$
(27)
We see that the result of Eq.(24) can be recovered by identifying $`n`$ with $`x/W`$, as it is legitimate to do when the number of jumps done is so large as to perceive $`n`$ as a continuous number.

Let us focus now our attention on the case where the distribution of waiting times, $`\psi (t)`$, has the following asymptotic behavior

$$\psi (t)\frac{const}{t^\mu },$$
(28)
with $`\mu >1`$ so as to insure the normalization condition. Let us focus our attention on the case where

$$1<\mu <2.$$
(29)
This is the region where jumps of intensity $`x`$ occurring with equal probability in both the positive and the negative direction of the $`x`$-axis, i.e. the SJM, would result in a diffusion process with the scaling $`\delta `$ of Eq.(17). We now show that, as a benefit of adopting the AJM rather than the SJM, we derive an enhanced scaling parameter. In fact, for very small values of $`s`$,

$$\widehat{\psi }(s)1cs^\gamma ,$$
(30)
where

$$\gamma =\mu 1.$$
(31)
This is the same property as that used by Zumofen and Klafter. By plugging Eq.(30) into Eq.(24) and by inverting the Laplace transform with the nice method of ref., we obtain

$$p(x,t)\frac{t}{x^{\frac{1}{\gamma }+1}}L(\frac{t}{x^{\frac{1}{\gamma }}},\gamma ,\gamma ),$$
(32)
where $`L(y,\gamma ,\gamma )`$ denotes a fully asymmetric Lévy stable law of index $`\gamma `$. It is straightforward to express this important result in the form of Eq.(9) with the time independent scaling

$$\delta =\mu 1.$$
(33)
By comparing Eq.(33) to Eq.(17) it becomes evident the scaling enhancement produced by making the diffusion process asymmetric. In fact, at $`\mu =3/2`$ a transition is made from subdiffusion to superdiffusion.

### 3.2 Generalized Central Limit Theorem

The result of Section 3.1 can be recovered from within the attractive perspective of the Generalized Central Limit Theorem (GCLT) in a form derived from properly adopting the procedure of Feller  to our purposes. To make closer the connection with the formalism of Feller let us assume $`W=1`$. The function $`p(x,t)`$ in this case is nothing but the probability for the walker to make $`x`$ steps up to time $`t`$, or, equivalently, the probability of occurrence of $`x`$ events up to time $`t`$. It has to be stressed that according to the perspective expressed by Refs., in accordance with the arguments of Gaspard and Wang, this way of proceeding can also be used to measure the complexity or randomness of the process under study, in the special case where the time series under study is totally memory-less. In fact, in this specific case the occurrence of any sporadic event does not have any connection with the earlier or subsequent events. It has to be pointed out that the kind of time series that we are having in mind in this paper is characterized by events of the same intensity, or made so by ignoring the real intensity of each event. If there is no correlation between intensity and time occurrence of the event, the simplifying assumption that the events have the same intensity does not have any significant consequence on the scaling of the resulting diffusion.

If the simplifying assumption applies, that no form of memory exists, in addition to the inverse power law nature of the waiting time distribution $`\psi (t)`$, we can safely apply the Feller arguments. Feller  studied the case where the survival probability $`\mathrm{\Psi }(t)`$ has the following time asymptotic expression

$$\mathrm{\Psi }(t)\frac{1}{t^\gamma }.$$
(34)
The Feller notation is related to our notation by $`\gamma =\mu 1`$. In the article of Ref., Feller studies bot the case $`0<\gamma <1`$ and the case $`1<\gamma <2`$, which refers to our two regimes $`1<\mu <2`$ and $`2<\gamma <3`$, respectively. In the former case, Feller finds the important property

$$Pr\left(N_twt^\gamma \right)=G_\gamma (w^{1/\gamma }),$$
(35)
where, according to the notation used by Feller  $`Pr(N_t>y)`$ denotes the probability that the number of events occurring at times $`t^{}t`$ is greater than $`y`$. The function $`G_\gamma (x)`$ denotes the stable law whose characteristics function is given by

$$\widehat{G}_\gamma (z)=\mathrm{exp}\left\{|z|^\gamma \left[\mathrm{cos}\left(\frac{\pi \gamma }{2}\right)i\mathrm{sin}\left(\frac{\pi \gamma }{2}\right)sgn(z)\right]\mathrm{\Gamma }(1\gamma )\right\}.$$
(36)
The symbol $`\mathrm{\Gamma }`$ denotes the Euler function. The meaning of the function $`G_\gamma (w)`$ is made transparent by the case where the walker makes jumps in the same direction, with the probability $`\mathrm{\Pi }(x)1/x^\gamma `$, at regular interval of times. In this specific case, which is the case of an asymmetric Levy flight, the position of the random walker at time $`t`$, $`X_t`$, fulfills the property

$$Pr\left(X_twt^{1/\gamma }\right)=G_\gamma (w).$$
(37)
The function on the right hand side of Eq. (37) is the same as that appearing in the right hand side of Eq. (35). However, we know that it has to do with an asymmetric Levy process. Consequently, it is a stable law that the theory of  makes it possible to write under the form

$$G_\gamma (w)=\underset{0}{\overset{w}{}}L(y,\gamma ,\gamma )dy,$$
(38)
where $`L(y,\gamma ,\gamma )`$ is the function defined in Section 3.1. We note that Eq. (35) can also be written as

$$_{wt^\gamma }^{\mathrm{}}dN_tp(N_t)=G_\gamma (w^\gamma ).$$
(39)
By differentiating Eq. (39) with respect to $`w`$ and setting$`N_t=x`$ (the position of the walker at time $`t`$, labeled as $`x`$, is equal to the number of steps taken up to time $`t`$, namely, $`N_t`$) we get

$$p(x,t)\frac{t}{x^{\frac{1}{\gamma }+1}}L(\frac{t}{x^{\frac{1}{\gamma }}},\gamma ,\gamma ),$$
(40)
which agrees with the scaling predicted by the CTRW arguments used in section 3.1 since it leads to the scaling parameter $`\delta =\mu 1`$.

In the case where $`1<\gamma <2`$ we are led by the corresponding Feller prescription  to express the probability distribution in the reference translating with the mean velocity (the inverse of the mean waiting time). With an argument similar to the one earlier used we obtain

$$p(x,t)\frac{1}{t^{\frac{1}{\gamma }}}L(\frac{x}{t^{\frac{1}{\gamma }}},\gamma ,\gamma ),$$
(41)
and it is now easy to prove that this yields, for $`2<\mu <3`$, the scaling prescription

$$\delta =\frac{1}{\mu 1}.$$
(42)
Results are summarized in Tab. 1.

### 3.3 Numerical Calculations

This section is devoted to illustrating the results of numerical calculations that confirm the benefits of the strategy suggested in this paper. First of all with Fig.2 we show that the AJM makes it possible to reach much earlier the correct Lévy scaling of Eq.(14). The numerical work shows that the time scale at which the DE method yields the correct scaling is, in this case, equal to about $`2,000`$ ($`T=10`$).

We also prove that the resulting scaling value is a genuine reflection of the property of Eq.(2). To do so, we express the diffusion process in the reference framework moving with the velocity $`TW/(\mu 2)`$ corresponding to the theoretical prediction of section 3.2. Then we consider a reference time $`t_{n_0}`$ and later times $`t_n>t_{n_0}`$. For any time $`t_n`$ we consider the squeezing parameter $`R(n)=(t_n/t_{n_0})^\delta `$, with $`\delta `$ given by the scaling detected by the DE method. The distance scale is squeezed by the quantity $`R(n)`$ and the ordinate scale is enhanced by the quantity $`1/R(n)`$. If the scaling detected by the DE is correct and if this procedure is applied to the distribution corresponding to the time $`t_n`$ this scaling procedure should make it identical to the distribution corresponding to the time $`t_{n_0}`$. Fig. 3 shows that this is really so, since all the distributions at times $`t_n`$ are made to virtually coincide with the distribution at time $`t_{n_0}`$.

We have explored also the case $`\mu <2`$. Figure 4 shows the result of the DE method in this case, and Fig.5 shows the result of scaling back to earlier times the distributions corresponding to later times, in this case. The result is not as accurate as that illustrated by Fig.3. However, we judge it to be satisfactory, especially if we take into account the fact that the region $`\mu <2`$ corresponds to the case of extremely large time correlations, with an infinite sojourn time.

## 4 The Method in Action in the Case of Real Sequences

As an example of how this technique can be applied to real data we considered time series concerning human heartbeats. These series were taken from Physionet  and are a record of time intervals between successive beats. To apply our method we adopted a coarse graining procedure: times between beats were first multiplied by $`10`$, with only the integer part of the resulting number taken into account. This sequence of integer numbers was then transformed into a new one according to the following prescription: we wrote $`0`$ if the number was equal to the preceeding one, $`1`$ if it was different. We note that this would have resulted almost in a sequence of only 1’s, without any 0, had we not adopted the coarse graining procedure. From this sequence we generated an asymmetric diffusion process (AJM) which we then studied by means of DE method described in the previous sections to derive the scaling parameter $`\delta `$. This corresponds to having a walker taking a step in a fixed direction when the time interval between successive beats, as seen on a coarse grained scale, differs from the preceding one. We also analyzed the probability distribution $`\psi (t)`$ of waiting times between the $`1^{}`$s, which shows an inverse power-law behavior with index $`\mu >2`$. We thus expect the correct scaling to be predicted by Eq.(42) rather than by Eq.(33). Results obtained from different databases are reported in table 2, whose third column represents the value of $`\mu `$ we derive from $`\delta `$. An example for a particular database (Ltdb15814) is shown in Fig.6 where we derive the scaling $`\delta `$ using DE, and in Fig.7, where we show how the inverse power law of index $`\mu =1+1/\delta `$ indeed fits the data used to obtain $`\psi (t)`$.

## 5 Concluding Remarks

This paper is a sequel to that of Ref., where the method of DE was applied for the first time. However, the emphasis here is on the detection of the scaling more than on the property that the authors of Ref. referred to as memory. On the basis of the results of the present paper, the memory of Ref. seems to have to do with the transition from the short-time to the long-time regime. The main result here reached, however, has to do with proving the convenience of recording the times at which events occur, an event being a sign change, regardless of whether it implies an increase or a decrease of the intensity signal. If we record the time at which the event occurs, and we interpret it as a jump of the random walker always in the same direction, with a fixed intensity, then the diffusion process becomes asymmetric, with the big benefit of enhancing the value of the scaling parameter in the interval $`1<\mu <2`$, and of making it sensitive to the power law index $`\mu `$ in the whole interval $`1<\mu <3`$. The AJM has to be compared to the SJM, a model that the theoretical analysis of the 1974 pioneer work of Shlesinger  shows to be insensitive to $`\mu `$ in the range $`2<\mu <3`$. The value of the scaling parameter $`\delta `$, stemming from the AJM, as earlier repeatedly pointed out, is enhanced with respect to that of the SJM, in the whole interval $`2<\mu <3`$. If we compare the AJM to the SVM of Ref., we find that it provides information on $`\mu `$ even in the interval $`1<\mu <2`$, where the SVM would steadily result in the ballistic prescription $`\delta =1`$.

It is interesting to notice that the AJM proposed by this paper emphasizes the phase transition character of a transformation moving from $`\mu >2`$ to $`\mu <2`$. As pointed out in Ref., this is a transition from a stationary to a non-stationary regime. In fact, Ref.  shows that $`\mu <2`$ implies the lack of an invariant distribution, thereby making any observation carried out for a finite time equivalent to monitoring an out of equilibrium process. We are convinced that the systematic application of the prescriptions of this paper to studying time series generated by complex dynamic processes, might lead to discover these non-stationary cases.

The application to the real sequences of physiological interest of Section 4 shows that the detection of the scaling parameter $`\delta `$ is much more accurate than the direct evaluation of $`\mu `$ from the experimental data, the only ambiguity being the fact that the value of $`\delta `$ is compatible with both $`\mu (1)=1+1/\delta `$ and $`\mu (2)=1+\delta `$. So, the experimental error affecting the direct determination of delta must only be smaller than the difference $`\mu (1)\mu (2)`$. Only when the value of $`\mu `$ is very close to the transition value $`\mu =2`$, the adoption of the DE method to detect $`\mu `$ becomes ambiguous, and additional criteria must be found.

A possible efficient criterion could be given by the joint use of the SVM and AJM method. In fact, as clearly illustrated by Table 1, if the SVM yields $`\delta =1`$, this means that $`\mu <2`$. If, on the contrary, the SVM, within the statistical accuracy of the numerical procedure, yields the same scaling as the SVM, we have to select $`\mu >2`$.

Another benefit of the adoption of the AJM is the faster attainment of the correct scaling. It is not yet quite clear which is the true reason why the correct scaling of (14), in the region $`2<\mu <3`$, is significantly faster than with the SVM. One reason might be that the probability distribution of the AJM diffusion process is sharper than that of the SVM. Some more research work is required to assess this point. We hope that the results already obtained are interesting enough as to draw the attention of the researches of this field of investigation on the benefit stemming from converting the experimental data into asymmetric diffusion processes. As the last but not least remark, we want to point out that the DE method is probably the most convenient way to detect the scaling of an asymmetric diffusion processes, since this technique is not affected by the presence of biases, and does not require the adoption of any form of de-trending.


text:
# Combinatorial nature of ground state vector of O(1) loop model

## Abstract

Hanging about a hypothetical connections between the ground state vector for some special spin systems and the alternating-sign matrices, we have found a numerical evidence for the fact that the numbers of the states of the fully packed loop model with fixed link-patterns coincide with the components of the ground state vector of the dense O$`(1)`$ loop model considered by Batchelor, de Gier and Nienhuis. Our conjecture generalizes in a sense the conjecture of Bosley and Fidkowski, refined by Cohn and Propp, and proved by Wieland.

I would say that imagination is a form of memory.

(V. Nabokov)

In paper  we made some conjectures related to combinatorial properties of the ground state vector of the XXZ spin chain for the asymmetry parameter $`\mathrm{\Delta }=1/2`$ and an odd number of cites. In the subsequent paper  Batchelor, de Gier and Nienhuis considered two variations of this model along with the corresponding O$`(n)`$ loop model at $`n=1`$ and notably increased the number of models and related combinatorial objects (see also ).

During last months we have accumulated a lot of data on relation between the spin systems and alternating sign matrices (ASMs). In this note we limit ourselves to one qualitative conjecture which relates some classes of the states of the fully packed loop (FPL) model with the ground state vector of the dense O$`(n)`$ loop model (see  and references therein). The states the FPL model is in bijective correspondence with the ASMs. Therefore if our conjecture is true then one can investigate the FPL model, and the ASMs, using the methods elaborated in mathematical physics for integrable model (see, for example, ).

The background information on the ASMs and their different combinatorial forms can be found in the recent review by Propp  and in references therein. For more details on enumerative problems related to the states of the FPL model see the article by Wieland .

Following the review paper by Propp  we define the “generalized tic-tac-toe” graph as the graph formed by $`n`$ horizontal lines and $`n`$ vertical lines meeting $`n^2`$ intersections of degree $`4`$, with $`4n`$ vertices of degree $`1`$ at the boundary. Then we number the vertices of degree $`1`$. We start with the left top vertex and number clockwise every other vertex. Now consider subgraphs of the underlying tic-tac-toe graph such that each of the $`n^2`$ internal vertices lies on exactly two of the selected edges and each numbered external vertex lies on a selected edge, while each unnumbered external vertex does not lie on a selected edge (see, for example, figure 2).

These subgraphs are the states of the FPL model. The states the FPL model are in bijective correspondence with the states of the square-ice model with the domain-wall boundary conditions , and with the ASM matrices. In particular, the number of such states is equal to the number of ASMs, usually denoted by $`A_n`$. Each state of the FPL model define a so-called link-pattern describing the pairings of the external vertices. We depict such a pattern as a circle with $`2n`$ vertices placed on it and connected pairwise inside the circle without intersection (see figure 2). Having in mind the relation to ASMs we denote the number of the states of the FPL model corresponding to the link-pattern $`\pi `$ by $`A_n(\pi )`$.

Consider as an example the case of $`n=4`$. In this case the FPL model has 42 states and 14 different link-patterns which are given in table 1.

The coincidence of the numbers $`A_n(\pi )`$ for different link-patterns is explained by the conjecture of Bosley and Fidkowski, refined by Cohn and Propp, and proved by Wieland . In accordance with this conjecture, if two link-patterns $`\pi `$ and $`\pi ^{}`$ are connected by a rotations or a reflection then $`A_n(\pi )=A_n(\pi ^{})`$.

Let us now formulate a conjecture on the numbers $`A_n(\pi )`$. To this end define $`2n`$ operations $`h_i`$, $`i=1,\mathrm{}2n`$, on link-patterns in the following way. For a fixed $`i`$ consider a general link-pattern. Let this link-pattern has the $`i`$-th vertex linked to the $`j`$-th one and the $`(i+1)`$-th vertex linked to the $`k`$-th one. Two different case are possible, the first is when $`j=i+1`$ and $`k=i`$. In this case the operation $`h_i`$ left the link-pattern unchanged. In the second case $`ji+1`$ and $`ki`$. In this case the operation $`h_i`$ replaces the two links under consideration by the link of $`i`$-th vertex and the $`(i+1)`$-th vertex and the link of the $`j`$-th vertex and the $`k`$-th vertex. To be understandable by the widest audience of readers, we decided to formulate our conjecture in three different ways.

The first formulation is in terms of the game theory. Let two players A and B play a game of chance. Both persons randomly choose a state of the FPL model. The player A checks a link-pattern of his state. The player B has to choose randomly one of the operations $`h_i`$ and to apply it to the link-pattern of his state. Only the resulting link-pattern is of importance for him. Players aim at obtaining some fixed link-pattern.

###### Conjecture 1

Both players have equal chances to win.

Consider, for example, the case of $`n=4`$ with the wining link-pattern be the very first one of table 1. The chances of the player A are $`𝖯_A=7/42=1/6`$. Let us evaluate the chances of the player B.

$`𝖯_\mathrm{B}`$ $`=`$ $`𝖯\left[\text{}\right]𝖯[\{h_1,h_3,h_5,h_7\}]+𝖯\left[\text{}\right]𝖯[\{h_5,h_7\}]`$
$`+`$ $`𝖯\left[\text{}\right]𝖯[\{h_1,h_7\}]+𝖯\left[\text{}\right]𝖯[\{h_1,h_3\}]`$
$`+`$ $`𝖯\left[\text{}\right]𝖯[\{h_3,h_5\}]+𝖯\left[\text{}\right]𝖯[\{h_3,h_7\}]+𝖯\left[\text{}\right]𝖯[\{h_1,h_5\}].`$
Using the data from table 1 one obtains

$$\frac{7}{42}\times \frac{4}{8}+\frac{3}{42}\times \frac{2}{8}+\frac{3}{42}\times \frac{2}{8}+\frac{3}{42}\times \frac{2}{8}+\frac{3}{42}\times \frac{2}{8}+\frac{1}{42}\times \frac{2}{8}+\frac{1}{42}\times \frac{2}{8}=\frac{1}{6}.$$

Thus our conjecture in this case is true.

The second formulation is in terms of the combinatorics.

###### Conjecture 2

For any $`n=1,2,\mathrm{}`$, one has

$$\underset{i=1}{\overset{2n}{}}\underset{\pi ^{}\mathrm{\Pi }_i(\pi )}{}A_n(\pi ^{})=2nA_n(\pi ),$$

where $`\mathrm{\Pi }_i(\pi )`$ is the set of link-patterns $`\pi ^{}`$ such that $`h_i(\pi ^{})=\pi `$.

The last formulation is in terms of the dense periodic O$`(1)`$ model. The state space for this model is the set of all formal linear combinations of the link-patterns. The operations $`h_i`$ extended to this space by linearity become linear operators. The Hamiltonian of the model is the sum of the operators $`h_i`$:

$$H=\underset{i=1}{\overset{2n}{}}h_i.$$
(1)
Consider the action of this operator on the vector

$$\mathrm{\Psi }=\underset{\pi }{}\pi A_n(\pi ).$$
(2)
We have

$$H\mathrm{\Psi }=\underset{i=1}{\overset{2n}{}}\underset{\pi }{}h_i(\pi )A_n(\pi )=\underset{\pi }{}\pi \underset{i=1}{\overset{2n}{}}\underset{\pi ^{}\mathrm{\Pi }_i(\pi )}{}A_n(\pi ^{})=2n\underset{\pi }{}\pi A_n(\pi ),$$

where the last equality follows from conjecture 2. Hence we can write

$$H\mathrm{\Psi }=2n\mathrm{\Psi }.$$

Thus, the vector $`\mathrm{\Psi }`$ is an eigenvector of the Hamiltonian of the dense periodic O$`(1)`$ loop model with an odd number of sites.

Let us consider our usual case of $`n=4`$. The matrix of the Hamiltonian $`H`$ and the component form of the vector $`\mathrm{\Psi }`$ have in this case the form

$$H=\left(\begin{array}{cccccccccccccc}4& 0& 2& 0& 2& 0& 2& 0& 2& 0& 2& 0& 2& 0\\ 0& 4& 0& 2& 0& 2& 0& 2& 0& 2& 0& 2& 0& 2\\ 1& 0& 3& 0& 0& 1& 0& 1& 0& 0& 0& 2& 0& 0\\ 0& 1& 0& 3& 0& 0& 1& 0& 1& 0& 0& 0& 2& 0\\ 1& 0& 0& 0& 3& 0& 0& 1& 0& 1& 0& 0& 0& 2\\ 0& 1& 1& 0& 0& 3& 0& 0& 1& 0& 2& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 3& 0& 0& 1& 0& 2& 0& 0\\ 0& 1& 1& 0& 1& 0& 0& 3& 0& 0& 0& 0& 2& 0\\ 1& 0& 0& 1& 0& 1& 0& 0& 3& 0& 0& 0& 0& 2\\ 0& 1& 0& 0& 1& 0& 1& 0& 0& 3& 2& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 1& 2& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 1& 0& 0& 0& 0& 2& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 1& 0& 0& 0& 0& 2& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0& 0& 0& 2\end{array}\right),\mathrm{\Psi }=\left(\begin{array}{c}7\\ 7\\ 3\\ 3\\ 3\\ 3\\ 3\\ 3\\ 3\\ 3\\ 1\\ 1\\ 1\\ 1\end{array}\right).$$

Here we take the basis vectors in the order given in table 1. Note that the sum of the matrix elements of $`H`$ belonging to any its fixed column is equal to 8 which coincides with $`2n`$. This is a general property of the Hamiltonian $`H`$ valid for any $`n`$ which implies that the spectral radius of $`H`$ is equal to $`2n`$ (see, for example, \[9, section 9.3\]). Therefore, the vector $`\mathrm{\Psi }`$ corresponds to the largest eigenvalue of the Hamiltonian $`H`$. Actually this eigenvalue is nondegenerate that can be proved using the positivity of the corresponding transfer matrix. Taking all this into account we give the last formulation of our conjecture.

###### Conjecture 3

The vector $`(\text{2})`$ is the unique eigenvector of the Hamiltonian $`(\text{1})`$ corresponding to its largest eigenvalue $`2n`$.

In paper  Batchelor, de Gier and Nienhuis formulated two conjectures on the eigenvector under consideration111Actually Batchelor, de Gier and Nienhuis take the Hamiltonian with the minus sign. In this case the eigenvector under consideration is the ground state of the model.. They conjectured that the sum of its components is equal to $`A_n`$ and its largest component is equal to $`A_{n1}`$. Our conjecture reveals the combinatorial nature of all components of the eigenvector with the largest eigenvalue.

Note that the Hamiltonian $`H`$ has the rotational and reflection invariance. It can be easily shown that the eigenvector with the maximal eigenvalue must have rotational and reflection invariance. Hence, our conjecture generalizes partially the conjecture proved by Wieland in paper  and mentioned above.

In conclusion note that we have verified our conjecture up to $`n=7`$.

Acknowledgments The authors would like to thank Bernard Nienhuis for his clear explanation of the Hamiltonian for the dense O$`(1)`$ loop model. We have to acknowledge the beautiful survey of James Propp which inspired us to investigate the FPL model. We are also thankful to him for his interest to our work and discussions.

The work was supported in part by the Russian Foundation for Basic Research under grant #01–01–00201.


text:
# On the Power Spectrum Density of Gamma Ray Bursts

## 1  INTRODUCTION

The origin and the fundamental mechanism of the gamma-ray bursts (GRBs) have not been revealed despite the fact that more than thirty years have passed since the first discovery of them. In the numerous luminosity profiles of GRBs, we notice that most GRBs have very rapid milli-second time variability in their intensity (Schaefer & Walker 1999) though there is vast variety in luminosity profile itself for each GRB. We believe that this characteristic variability must have important information to reveal the GRB mechanism. Therefore, in this letter, we would like to concentrate on the analysis of this variability.

In studying time sequence of objects, the Fourier transform technique is useful and often yields indispensable information on the scaling properties and characteristic time scale. For example, Beloborodov et al. (1998) used this method for 214 light curves of long GRBs ($`T_{90}>20\mathrm{sec}`$) and reported that the averaged power spectrum density (PSD) of GRBs shows the power-law with the index $`1.67\pm 0.02`$ over two decades in frequency range. This value is very closed to $`5/3`$ which suggests the Kolmogorov spectrum of velocity fluctuations in turbulent medium.

On the other hand, the power-law in PSD does not necessarily specify the whole mechanism as we know in various examples. Actually the power-law in PSD can be derived by many reasons; some specific burst profile in luminosity, Levy-type random noise in the background, and the superposition of similar shots and so on. For example, the power-law in PSD of X-ray emission from the Active Galactic Nuclei (AGN) is considered to be originated from the superposition of many similar shots (Mineshige & Yonehara 2001). Chang & Yi (2000) tried to reproduce the power-law index, -5/3, and the individual PSD distribution of GRBs by the superposition of many decaying pulse shots.

Here in this paper, we would like to determine the origin of the power-law in PSD of GRBs by analyzing the detail of PSD for each burst data.

First, we demonstrate that the individual PSD of GRBs does not exactly have the power-law index $`1.67`$ but has wide variation which is well approximated by the Gaussian distribution. Second, we argue that power-law index of PSD is determined by the profiles of an individual shot in the light curve and not by the superposition of many similar shots nor by Levy-type random noise. Then we discuss the general nature of PSD of superposed shots and apply this method to the actual data of GRB light curves. Finally we observe how the averaged PSD shows clear power-law index $`1.67`$.

## 2 POWER LAW IN PSD OF GRB

In our analysis, we use the data detected by the Burst and Transient Source Experiment (BATSE) on Compton Gamma Ray Observatory (CGRO) with 64ms resolution 111 ftp://cossc.gsfc.nasa.gov/pub/data/batse/ascii\_data/64ms/. We use the light curves in the energy band $`50<h\nu <100`$ keV, excluding unclear high and low energy bands, and subtract the background by the linear fitting. Then for excluding the noise contamination, we select the data with peak count rate larger than 250 counts per 64ms bin; finally 297 data set remain.

In order to see the generality of the power index $`5/3`$, we first calculate the power spectra, which is the absolute square of the Fourier transform of the time sequence, for individual light curves. We assume the following fitting function for the power spectrum $`P(f)`$, expecting the coexistence of the power-law component and the thermal white-noise component.

$$P(f)=Af^\alpha +B.$$
(1)
Figure 1 shows the distribution of power-law index $`\alpha `$ thus obtained. The distribution is well fitted by the Gaussian and the mean value is $`1.76`$, which is close to the value $`1.67`$ reported in the reference Beloborodov et al. (1998). We emphasize here that the individual data has the power-law PSD though the power index is widely distributed with the variance 0.65.

## 3 WHAT DETERMINES THE INDIVIDUAL POWER LAW?

Now let us consider the meaning of the power-law in PSD. As is well known, the power-law itself in the PSD can be realized in various origin. The specifications of the origin is significant for the analysis of the central engine of GRBs.

Mineshige & Yonehara (2001) shows at least the following three ways to realize the power-law in PSD.

1. Specific global profile of a burst.
2. Levy-type random process.
3. Superposition of self-similar shots.

We examine the above possibilities in the following in this order.

Examining the first possibility, we realize that each GRB data generally has multiple shots and rapid fluctuations as well as noise whose origin is not seem to be GRB. We try to identify these components in each GRB data step by step.

First we discuss the light curve with a single shot. This simple type of light curve shows a clear power-law behavior in its PSD. We consider the following two model cases; a) the decay type, and b) the grow-and-decay type.

a) The decay type

The decay type shot has the flux $`x`$ at time $`t`$ as

$$x(t)=ht^p\mathrm{exp}(f_0t)\theta (t),$$
(2)
where the positive parameters $`h`$, $`p`$ and $`f_0`$ respectively represent “intensity”, “sharpness” and “inverse of the duration” of the shot. $`\theta (t)`$ is the step function. Power spectrum $`P`$ for this shot is written as

$$P(f)=h^2[\mathrm{\Gamma }(1p)]^2\left[f_0^2+f^2\right]^{p1},$$
(3)
where $`\mathrm{\Gamma }(x)`$ is the Gamma function.

b) The grow-and-decay type

The grow-and-decay type shot generally has asymmetric profile in the burst; we separate the grow (left side) and decay (right side) of the shot. Each side of the shot is written as

$`x(t)`$ $`=`$ $`h_L(t)^{p_L}\mathrm{exp}(tf_{0L})\theta (t)+h_R(t)^{p_R}\mathrm{exp}(tf_{0R})\theta (t)`$ (4)
Power spectrum of this type is written as

$`P(f)`$ $`=`$ $`h_L^2\{\mathrm{\Gamma }(1p_L)\}^2\left(f_{0L}^2+f^2\right)^{p_L1}+h_R^2\{\mathrm{\Gamma }(1p_R)\}^2\left(f_{0R}^2+f^2\right)^{p_R1}`$ (5)
$`+2h_Lh_R\mathrm{\Gamma }(1p_L)\mathrm{\Gamma }(1p_R)`$
$`\times \left(f_{0L}^2+f^2\right)^{\frac{p_L1}{2}}\left(f_{0R}^2+f^2\right)^{\frac{p_R1}{2}}\mathrm{cos}\left[(p_L1)\theta _L+(p_R1)\theta _R\right],`$
where $`\theta _{\mathrm{L},\mathrm{R}}=\mathrm{arctan}\left(f/f_{0\mathrm{L},\mathrm{R}}\right)`$.

Since the GRB data generally has multiple shots, the above single shot profiles must be superposed before we use. For simplicity for fitting procedure, we assume that all shots in the individual light curve has the same profile except the overall amplitude 222 When we fit actual data of GRBs, we use inverse of duration $`f_0`$ as a fitting parameter, but it does not change the index of PSD at all. ; i.e. the $`k`$-th shot in the light curve $`x_k(t)`$ is written as $`A_kx(ta_k)`$, where $`A_k`$ is the relative amplitude and $`a_k`$ is the location of the $`k`$-th shot with $`x(t)`$ being the fixed function for each GRB data. Then we can easily calculate the PSD for multiple shots since the Fourier transforms of the first and the k-th shots are simply related with each other as $`\stackrel{~}{x}_k=A_ke^{ifa_k}\stackrel{~}{x}_1`$.

The light curve $`x_T(t)`$ which has $`n`$ shots

$$x_T(t)=\underset{k=1}{\overset{n}{}}x_k(t),$$
(6)
has the PSD

$`P_T(f)`$ $`=`$ $`P(f)\left[{\displaystyle \underset{j=1}{\overset{n}{}}}A_j^2+{\displaystyle \underset{j>k}{}}2A_jA_k\mathrm{cos}f(a_ja_k)\right],`$ (7)
where $`P(f)`$ is PSD for $`x(t)`$.

It is apparent that the last expression consists of the “constant” term and the “oscillating” term; the former term determines the global profile of PSD and the latter fluctuations around it.

We can actually see this structure in PSD of GRB data. We first identify the maximum shot in the GRB data and locally fit this shot by the single-shot form argued in the above. Then we subtract this first fit from the original data yielding one-shot subtracted data. Second, we identify the maximum shot in this subtracted data and locally fit this shot by the single-shot form with the same parameter $`p`$. Then we subtract this second fit from the one-shot subtracted data yielding two-shots subtracted data. After repeating this process several times, the reduced data becomes almost flat.

We applied this method for many individual GRB data. One typical example for GRB910602(#257) is shown in Fig.2. Here we use the decay type shot for fitting identifying two shots in the data. The full PSD of this light curve is plotted in Fig.3 with the solid line. The PSD of the first fit is plotted with the chain line which already shows smooth power-law. The PSD of the superposition of the first and the second fits is plotted with the broken line which shows the oscillation around the chain line. It is almost clear that the PSD of the superposition of first few fits are sufficient to faithfully reproduce the main part of the original PSD. Other GRB data sets show similar behavior as this one.

Thus we have confirmed the first possibility: Specific global profile of a burst determines the power-law behavior in PSD.

In order to exclude the second possibility, the Levy-type random process, we have examined the peak distribution analysis. We calculate the distribution function of all peak-intensity in the data of artificially produced Levy-type light curve and that in the actual data of GRBs. The former shows a clear power law however the latter does not. We can now claim that the power-law in PSD of GRBs is not originated from the Levy-type random process.

In order to exclude the third possibility i.e. superposition of self-similar shots, we calculate the PSD of the few-shots subtracted data. For the previous data set, we calculate the PSD of the two-shots subtracted data in Fig.3 with the dotted line. This residual component is 10-100 times smaller than the main component and is almost flat. If GRBs have self-similar structure in their light curves, the few-shots subtracted data should also show power law (but its range must be narrower than that of original data). Thus we can now claim that the power-law in PSD of GRBs is not originated from the superposition of self-similar shots

We have also applied the same analysis to the data set of the soft gamma-ray repeaters (SGR). We use the data of SGR1900$`+`$14 333 http://ssl.berkeley.edu/ipn3/sgr1900$+$14.lightcurve. The first one-shot fit has been enough to reproduce the power-law behavior in its original PSD 444Fitting paramerer $`p`$ is $`0.045`$ and power index of original PSD is $`1.8`$. and the one-shot subtracted data shows a flat PSD except several very sharp peaks which are considered to be originated from the central pulsar of SGR.

## 4 WHAT DETERMINES THE AVERAGED POWER LAW?

We have averaged all the 297 PSD of GRB data and obtain the Figure 4. Before taking the average, we have normalized each GRB luminosity data so that the maximum of the count rate be unity. This is because the original light curves have two to three orders of difference in maximum count rates and therefore the naively averaged PSD is determined by few GRB data. Since the magnitude of the count rate itself is not physical (we don’t know the distance!), such discrimination is not a proper manipulation. We have also tried the total-count-rate (fluence) normalization as well as maximum-count-rate normalization; yielding no significant difference between them.

Thus averaged PSD shows much clear power-law (Fig.4) than the individual PSD of GRB data. From Fig.4 we again observe that the individual PSD shows power-law and the power index simply fluctuates variously. After the superposition of many PSD data, there appears smooth power-law behavior with the power index of the central value. It is important to observe in Fig. 4 that there is no systematic correlation in each PSD to yield the global smooth power-law behavior. Especially there is no systematic distribution of time scales nor correlation of turning points which, if any, would have yielded clear envelope when many PSD are superposed.

The above reasoning for the appearance of the smooth power-law in averaged PSD is also supported by the fact that the power index $`1.67`$ which is close enough to the mean value $`1.76`$ of the index distribution within 14 percent of the variance. If the clear power-law were realized by the envelope of many systematic distribution of burst time scales, the averaged power index would be significantly smaller than the mean power index.

## 5 CONCLUSIONS AND DISCUSSIONS

Analyzing 297 power spectrum density (PSD) of GRBs, we obtain the following results in this paper. a) Individual GRB data shows power law behavior in the PSD. The power index ($`\alpha `$) distribution is well approximated by the Gaussian form with the mean $`1.76`$ and the variance $`0.65`$. b) Power law behavior in PSD for individual GRB data is determined by the shot profile and not by the Levy-type noise nor the superposition of many self-similar shots. c) Power law index in averaged PSD is simply determined by the mean value of the index ($`\alpha `$) distribution. We found no correlation mechanism for producing clear power law in averaged PSD.

From the result b), we doubt the turbulence (Beloborodov et al. 1998, 2000) as the origin of the power law in PSD of GRBs. This is because the self-similar cascade of eddy, which would naturally yield many self-similar shots in the emission, yields the Kolmogorov power law spectrum of velocity fluctuations in turbulent medium.

The result b) provides sharp contrast with the case of power law in the X-ray emission from AGN. Many authors have demonstrated that some self-similar cascade model, which yields many self-similar shots, successfully describes the emission profile and the power law in PSD of AGN. Thus we naturally expect that the emission mechanism of AGN and that of GRB are quite different with each other despite the similar power law behavior in their PSD.

As we have seen in our analysis, some special form of shots, which yield the power law in PSD, is very characteristic and would be a good criterion for restricting various models of GRB generation mechanisms. The special shot profile is more severe checking point of the validity of the model than simply the power-law behavior of PSD (Panaitescu et al. 1999; Spada et al. 2000). An urgent interest then would be the question whether the popular ”Internal shock model” can explain these special shot profiles.

There are still important questions to be answered in the near future. 1) The shots in the light curves of GRBs have distinctly different two types; the decay type and the grow-and-decay type. This property is different from that of blazars, even if their similarity in emission process is often suggested. 2) What is the distribution of the sharpness parameter $`p`$ and the ratio of power-law component and noise $`A/B`$ in equation (1) ? Is it also Gaussian or any other form from which we can extract the information of the distance or classification of GRB?

We all would like to thank Takashi Okamura for many useful suggestion. I. J. would also like to thank Prof. T. Yokobori and Prof. A. T. Yokobori, Jr for their hospitality.