text:
# Multi-order TODCOR: application to observations taken with the CORALIE echelle spectrograph

## 1 Introduction

TODCOR (TwO-Dimensional CORrelation) is a two-dimensional correlation technique to derive the radial velocities of both components of double-lined binary spectra (Zucker & Mazeh 1994). It was introduced as a generalization of the (one-dimensional) cross-correlation technique (Simkin 1974; Tonry & Davis 1979) to deal with the difficulties encountered in double-lined spectra, when the lines of the two components cannot be easily resolved. Assuming the observed spectrum is a combination of two known spectra shifted by the radial velocities of the two components, TODCOR calculates the correlation of the observed spectrum against a combination of two templates with different shifts. The result is a two-dimensional correlation function, whose peak simultaneously identifies the radial velocities of both the primary and the secondary.

One of the advantages of TODCOR is its ability to use different templates for the primary and the secondary. When the primary and the secondary are of different spectral types, the simultaneous use of two different templates utilizes all the spectral information contained in the observed spectrum of the combined system. This feature of TODCOR is most important when deriving the radial velocities of faint secondaries (Mazeh & Zucker 1994). However, the ability to detect faint secondaries using TODCOR is still limited, depending mainly on the signal-to-noise ratio and the number of spectral lines. Therefore, Mazeh et al. (2002) and Prato et al. (2002), who searched for faint companions, had to apply TODCOR on IR spectra, where the flux ratio is favorable for detecting cooler stars.

The modern spectrographs offer another path to enhance the ability to detect faint companions. Due to the progress in detector technology, many of the modern spectrographs produce multi-order spectra. In order to enhance our ability to detect a faint companion, while maximizing the precision of the measured velocities, we need to combine the spectral information in all the relevant orders. Originally, TODCOR was devised to analyse only single-order spectra, and further generalization is therefore needed in order to use the information in multiple orders.

When applying TODCOR to each order separately, the weak signal of the secondary may produce a certain local peak of the correlation function at the correct secondary velocity. However, this peak can be easily topped by spurious random peaks. This rules out, for example, calculation of the radial velocity independently for each order and then averaging the velocities, since many of the secondary velocities would have resulted from wrong peaks. Concatenation of the spectral orders to one single spectrum would require special treatment to the gaps and overlaps between adjacent orders. Co-adding overlapping regions and bridging the gaps both require interpolation, which would introduce artificial noise into the analysed spectrum.

The approach we suggest here is to calculate the correlation function for each order separately, and then combine the correlation functions of all the orders. The combination emphasizes the relevant correlation peak, and averages out the spurious random peaks. A simple average of the correlation functions may not be efficient enough, since the combination scheme has to consider the different spectral information in the different orders, and weigh them accordingly. Zucker (2003) introduces such a scheme, based on a few plausible statistical assumptions.

Using this scheme, multi-order TODCOR is applied in this work to the CORALIE spectra of HD 41004. As Santos et al. (2002) have shown, these spectra are composed of the spectra of the two visual components of the system, that are separated by $`0.5\mathrm{}`$. Santos et al. found that the radial velocities measured by CORALIE for these spectra showed a minute periodical variability. They suggested that this apparent variation was actually related to a much larger variation in the velocity of the faint B component. The presence of the spectrum of HD 41004 B caused a variable asymmetry in the line shape of the composite spectra, which was reflected in a minute variability of the measured radial velocity.

TODCOR, when applied to the multi-order CORALIE spectra of HD 41004, derived the radial velocities of both components without any assumptions regarding the orbit. The derived velocities confirm the conjecture of Santos et al.: the B component indeed moves periodically with a $`1.3`$-day period, and the A component shows a long-term trend. We present here an orbital solution for the B component and derive the slope of the long-term trend of the A component, based on the precise velocities derived by TODCOR.

The next section briefly reviews the previous results concerning HD 41004. The analysis and its results are presented in Sections 3 and 4. Section 5 discusses the results and their implications. Section 6 concludes the paper with a few remarks.

## 2 Characteristics of HD 41004

HD 41004 is a visual double system, consisting of a K1V–M2V pair. According to the Hipparcos catalogue, the pair is a common-proper-motion pair, with a V-magnitude difference of $`3.68`$ and a separation of $`0.541\mathrm{}\pm 0.033\mathrm{}`$. The basic stellar parameters of HD 41004 A are summarized in Table 1, reproducing Table 1 of Santos et al. (2002). The effective temperature, surface gravity and metallicity were calculated by Santos et al. through Strömgren photometry, which also showed that the star is photometrically stable within the instrumental precision. An independent estimate of the metallicity was obtained from analysis of the Cross-Correlation Function of CORALIE. The two independent estimates of the metallicity are quoted in the table.

The principal result of Santos et al. (2002) is the detection of a radial-velocity periodical variability in the CORALIE spectra of HD 41004. The variability pattern presented in their paper is consistent with the presence of a planet orbiting HD 41004 A with a $`1.3`$-day period. However, Santos et al. rejected this possibility based on the bisector shape analysis of the CORALIE spectra. This analysis revealed a periodic line shape variation, having the same period as the radial-velocity variation.

The interpretation Santos et al. suggested is the presence of an object orbiting the M2V star HD 41004 Ba. In this model A (the K star) and B are orbiting each other in a wide orbit, while B is further composed of two objects, Ba (the M star) and Bb, in a very close orbit, with a period of $`1.3`$ days. According to this interpretation, the CORALIE spectra of HD 41004 are a combination of the A and B components, because their separation is much smaller than the diameter of the CORALIE fiber. Since the flux of the B component in the relevant wavelength range is only about $`3\%`$ of the flux of the primary, large radial-velocity variations of B cause only small variations of the measured radial velocity, but their effect is manifested in the bisector shape.

Santos et al. performed some simulations to test their interpretation and concluded that their results are consistent with the presence of a brown dwarf orbiting HD 41004 Ba, with a radial-velocity amplitude of the order of $`5\mathrm{km}\mathrm{s}^1`$. They also found a long-term linear trend in the radial velocities, which did not agree with the older measurements. They ascribed it to the motion of HD 41004 A around HD 41004 B, but did not rule out an additional component that may be involved in this variation.

## 3 Analysis

The data we analysed comprised the 86 CORALIE spectra of HD 41004 used by Santos et al. (2002), which were obtained between November 2001 and February 2002, except for one spectrum obtained in December 2000. The system was further monitored and in this paper we add 17 spectra obtained as of March 2002.

CORALIE is a fiber-fed, cross-dispersed echelle spectrograph, mounted on the $`1.2`$-m Leonard Euler telescope at La Silla (Queloz et al. 1999). With a resolution of $`\lambda /\mathrm{\Delta }\lambda =\mathrm{50\hspace{0.17em}000}`$, it covers the wavelength range $`3800`$$`6900`$ Å, with $`68`$ echelle orders. The CORALIE system uses a software code that produces automatically the radial velocities, but for this work we used the reduced spectral orders.

The analysis used $`32`$ CORALIE orders within the spectral range $`4780`$$`6820`$ Å, after having excluded the orders that are heavily polluted by telluric lines. Bluer orders were excluded because the secondary signal is expected to be too weak in these orders. Multi-order TODCOR was used to combine the correlation functions from the $`32`$ orders. Note that the pixel-to-radial-velocity scale is different for each order. Therefore, before combining the correlation functions we had to interpolate and re-sample them into a pre-determined scale.

TODCOR requires two spectral templates as similar as possible to the expected primary and secondary spectra, and an assumed value for the flux ratio between the two spectra. We searched for the best templates in two datasets of spectra obtained by ELODIE and CORALIE. The CORALIE templates, appropriate for solar-type stars, were obtained as part of a program to derive precise abundances of planet-hosting and non-planet-hosting stars (Santos et al. 2000, 2001). The templates for the later stars were obtained using ELODIE as part of a program studying the binarity of close M dwarfs (Delfosse et al. 1998, 1999). The flux-ratio information was taken from Pickles (1998), using his measurements of typical stellar spectral energy distribution (SED), normalized according to the V-magnitude difference of HD 41004. By convolving the templates with a rotational broadening profile (e.g., Gray 1976), we created an additional degree of freedom which expanded our template library, allowing a better fit of the templates to the observed spectra.

In all the configurations of two templates that had spectral types similar to K1V and M2V, the secondary radial velocities exhibited a clear nearly circular orbit with an amplitude of around $`6\mathrm{km}\mathrm{s}^1`$, and a linear trend of the primary. We finally chose the configuration which yielded the smallest residuals in the secondary velocities, relative to the best-fit orbital solution. Table 2 lists the known stellar characteristics of HD 52698 and GJ 393 – the two templates finally chosen. The stellar parameters of HD 52698 (effective temperature, surface gravity and metallicity) were derived as in Santos et al. (2001), whereas an estimate of the $`v\mathrm{sin}i`$ was obtained using the calibration of the CORALIE Cross-Correlation Function width (see the Appendix in Santos et al. (2002)).

Table 2 quotes the original equatorial rotational velocities for both templates. Even mild broadening of the primary template ($`v\mathrm{sin}i=1\mathrm{km}\mathrm{s}^1`$), on top of its original $`v\mathrm{sin}i`$ of $`2\mathrm{km}\mathrm{s}^1`$, degraded the solution considerably. On the other hand, the secondary template could be broadened by $`v\mathrm{sin}i`$ up to $`5\mathrm{km}\mathrm{s}^1`$ without significantly affecting the solution. Finally we used a broadening of $`v\mathrm{sin}i=1\mathrm{km}\mathrm{s}^1`$, which yielded the smallest residuals of the secondary velocities.

Figure 1 shows the flux ratio as a function of wavelength, corresponding to the chosen spectral types of K1V and M2.5V, using the SED library of Pickles (1998). The dashed line in the figure shows the flux ratio obtained by assuming a black-body radiation law for the two spectra, with a temperature of 3500 K for the secondary. It is clear that while the black-body model roughly fits the detailed flux-ratio, it is not accurate enough and there are large differences, e.g., around $`6500`$ Å. In any case, using a constant flux ratio as in the single-order application of TODCOR is clearly not sufficient for such a wide spectral range.

Figure 2 demonstrates the way multi-order TODCOR measures the radial velocity of HD 41004 B. The Figure presents one-dimensional “cuts” of the two-dimensional correlation function. In each “cut” the correlation is shown as a function of the secondary velocity, while the primary velocity is fixed according to the location of the two-dimensional maximum. Because of the small flux-ratio between the two templates, we expect the correlation value to change very little when the secondary template alone is shifted. However, in each “cut” we still expect to see a relative peak around the correct secondary velocity. The Figure is divided into two columns, corresponding to two different exposures ($`\mathrm{JD}=\mathrm{2\hspace{0.17em}452\hspace{0.17em}312.722}`$ and $`\mathrm{JD}=\mathrm{2\hspace{0.17em}452\hspace{0.17em}343.590}`$). Using the best-fit orbit, the expected secondary velocity for the spectrum used in the left column of the Figure is $`35.2\mathrm{km}\mathrm{s}^1`$, while for the right column it is $`42.2\mathrm{km}\mathrm{s}^1`$. These velocities are marked on the figure by dashed lines. The primary velocity is represented by a dotted line in the left column, while in the right column it almost coincides with the secondary velocity. For this graphic demonstration only, a third-order best-fit polynomial was subtracted from all the plotted functions in order to accentuate the local peak.

The three upper panels in each column demonstrate the problems in the single-order TODCOR. On the left, the first upper panel exhibits a very prominent peak in a wrong velocity. The second and third panels show a moderate peak at about the expected velocity In almost all $`32`$ orders, some local peaks appeared around the correct velocity. The lower panel shows the result of combining the correlation functions of all the analysed orders. The correct peak is clearly emphasized relative to the spurious peaks.

Comparing the two columns of Figure 2, we see that we cannot know in advance which orders present the correct peak. Thus, we have to combine all $`32`$ orders in order to have the correct peak emphasized. The right column of the Figure demonstrates another advantage of TODCOR: in this spectrum the secondary velocity is almost identical to the primary velocity. In the conventional one-dimensional cross-correlation, there is no way to measure the two velocities, due to the blending of the two correlation peaks. TODCOR, by using two different templates, allows the measurement of both velocities.

The listing of the radial velocities of A and B and the corresponding times, can be obtained at the CDS.

## 4 Results

Figure 3 shows the resulting orbit of HD 41004 B, while the orbital elements are summarized in Table 3. Applying multi-order TODCOR to many orders ($`32`$), combined with the large number of measurements ($`103`$) yielded a very precise orbital solution. Thus, the radial-velocity amplitude ($`K`$) was found with a precision of about $`1\%`$. This fine precision allowed also a very accurate estimate of the orbital eccentricity ($`e`$). Although very small, the eccentricity is still non-vanishing, with a significance level of $`210^5`$ according to the Lucy & Sweeney (1971) test. Assuming a mass of $`0.4M_{\mathrm{}}`$ for HD 41004 Ba, the companion minimum mass is $`19M_\mathrm{J}`$. The uncertainty of $`0.25M_\mathrm{J}`$ does not take into account the uncertainty in the mass of HD 41004 Ba$`M_{\mathrm{Ba}}`$. A $`20\%`$ uncertainty in $`M_{\mathrm{Ba}}`$ would result in a $`2.3M_\mathrm{J}`$ uncertainty in $`M_{\mathrm{Bb},\mathrm{min}}`$.

Figure 4 shows the velocities of HD 41004 A. The velocities measured after $`\mathrm{JD}=2,452,200`$ show a clear linear trend. The best-fit line has a slope of $`+105\pm \mathrm{\hspace{0.17em}10}\mathrm{m}\mathrm{s}^1\mathrm{year}^1`$, and is also shown in the Figure. As Santos et al. (2002) have already noticed, the first isolated measurement does not agree with the linear trend implied by the other measurements. A Lomb-Scargle periodogram of the de-trended velocities (Figure 5) shows no hint of the $`1.3`$-day periodicity, again proving that only the B component participates in the periodic motion, as suggested by Santos et al. (2002). In light of the trend in the velocities of component A, we tried also to fit the B velocities with an additional trend but the resulting trend was not statistically significant.

The mean radial velocity of A during the linear part is $`42.5768\pm 0.0009\mathrm{km}\mathrm{s}^1`$, which is very close to the center-of-mass velocity of B – the difference is only $`\mathrm{\Delta }RV=0.95\pm 0.06\mathrm{km}\mathrm{s}^1`$. However, this difference depends also on the estimated velocities we used for the templates, and therefore we adopt a conservative error estimate of $`0.1\mathrm{km}\mathrm{s}^1`$ for $`\mathrm{\Delta }RV`$.

## 5 Discussion

The radial velocities derived in the previous section support the triple-system model suggested by Santos et al. for HD 41004. Our results allow, nevertheless, a somewhat more detailed study of the system.

### 5.1 The long-term trend of the radial velocity of HD 41004 A

As we have demonstrated in the previous section (see Figure 4), the radial velocities of HD 41004 A show a clear long-term variation over time. For the period November 2001 to February 2002 this variation can be nicely fitted by a linear increase. Note, however, that one single point, at $`\mathrm{JD}=\mathrm{2\hspace{0.17em}451\hspace{0.17em}902.774}`$, which was observed one year before all the others – on December 2000, deviates substantially from this linear fit. This point may indicate that the extrapolation of the linear approximation is no longer valid for the time of that point. Thus, we may actually be seeing the curved part of a long-term orbit, caused by a fourth component in the system. Such an object would have an orbital motion with a much shorter period than that of the AB system, which is at least 100 years, as inferred from the observed separation. In fact, follow-up measurements we are currently performing hint that this may indeed be the case, but a detailed solution of the orbit would still be premature.

### 5.2 The close orbit and the rotation of HD 41004 Ba

As mentioned in Section 3, the process by which we chose the templates allows a very crude estimate of the equatorial rotational velocity. The best-fit orbit was attained when the M2 template was convolved with a rotation profile of $`v\mathrm{sin}i=1\mathrm{km}\mathrm{s}^1`$. Our procedure is not accurate enough to use this value as a true measurement of $`v\mathrm{sin}i`$. However, trying large values for $`v\mathrm{sin}i`$ degraded the quality of the orbital solution. For example, while the quality of the solution remained almost unchanged for $`v\mathrm{sin}i`$ as large as $`5\mathrm{km}\mathrm{s}^1`$, the $`OC`$ RMS was doubled at $`v\mathrm{sin}i=10\mathrm{km}\mathrm{s}^1`$. Thus, we believe that the equatorial velocity of Ba is of the order of a few km s-1, and probably less than 10 km s-1. Even this crude estimate is sufficient to conclude that $`v\mathrm{sin}i`$ is substantially lower than the expected equatorial velocity of $`20\mathrm{km}\mathrm{s}^1`$, assuming synchronization (Santos et al. 2002).

At such a close orbit, with a period of $`1.33`$ days, it is usually assumed that synchronization and alignment of the spin with the orbital motion had been established (Santos et al. 2002). An obvious explanation of the apparently small $`v\mathrm{sin}i`$ would be small orbital inclination. Assuming alignment, small inclination implies that the mass of the unseen companion is considerably larger than its derived minimum value.

### 5.3 The small eccentricity of the close orbit

A binary or a star-planet system with a period as short as $`1.33`$ days is naïvely believed to have been circularized. However, our radial-velocity solution shows that the close orbit has a finite non-zero eccentricity – $`e=0.065\pm 0.014`$, which has to be explained.

Let us first estimate the timescale of circularization due to processes occurring in Ba. The relevant process is dissipation of the equilibrium tide through interaction with the convective envelope (Zahn 1989). We follow Rasio et al. (1996) and write:

$$\tau _{\mathrm{circ}}=\frac{\tau _\mathrm{c}}{f}\frac{M}{M_{\mathrm{env}}}\frac{1}{q(1+q)}(\frac{a}{R})^8,$$

where the parameter $`\tau _\mathrm{c}(MR^2/L)^{\frac{1}{3}}`$ is the eddy turnover timescale, and its numerical value is about $`0.5`$ yr. The numerical value of $`f`$ is obtained by integrating the viscous dissipation of the tidal energy through the convective zone and is of order unity (Zahn 1977). In our case the tidal pumping period – $`P/2`$ – is much smaller than $`\tau _\mathrm{c}`$, and therefore only convective eddies with turn-over time less than $`P/2`$ contribute to the dissipation. The value of $`f`$ is then reduced to $`(P/(2\tau _\mathrm{c}))^\alpha `$ (Zahn 1989). The correct value of $`\alpha `$ is debated but generally assumed to be either $`1`$ (Zahn 1992) or $`2`$ (Goldreich & Keeley 1977; Goldman & Mazeh 1991). We assume that the mass of the convective envelope – $`M_{\mathrm{env}}`$ – is very close to the stellar mass – $`M`$ (e.g., Chabrier & Baraffe 2000). The resulting timescale is about $`10^{10}`$ yr, assuming $`\alpha =1`$, and certainly much larger if $`\alpha =2`$. At the age of $`1.5`$ Gyr, it therefore seems that circularization through dissipation in Ba has not taken place.

Now let us examine the possible involvement of Bb in circularization processes. If the object is substellar or maybe even a very late M-dwarf, radiative zones in the atmosphere may form either due to the internal physics (Burrows et al. 1997) or due to the external heating by the primary (Guillot et al. 1996). These would necessitate another approach to calculate the circularization timescale, probably based on dissipation of the dynamical tide (Zahn 1977). The theory of orbital evolution through tidal dissipation is still debated, specially the mechanisms for dissipation of the dynamical tide (e.g., Claret et al. 1995), and therefore it is not clear whether it contributes significantly in the case of HD 41004 B.

## 6 Conclusion

We have presented in this work the first application of multi-order TODCOR to echelle spectra. The case of HD 41004 demonstrates the unique capabilities of this technique. It utilizes all the available prior knowledge regarding both spectral components (the different templates), it is not limited to a fixed flux ratio (we used the SED to calculate it) and it incorporates optimally the data from all the relevant spectral orders. Based on the best-fit solutions of the present case, we estimate the precision of the secondary velocity of HD 41004 to be $`0.56\mathrm{km}\mathrm{s}^1`$ and that of the primary velocities to be $`10\mathrm{m}\mathrm{s}^1`$. The radial velocities yielded accurate orbital elements of the unseen companion of HD 41004 B, and an accurate measurement of the radial acceleration of HD 41004 A. The latter suggested an additional component may be present in the system. Follow-up observations, currently underway, tend to confirm this hypothesis.

The combination of multi-order TODCOR together with the high signal-to-noise and high resolution of the CORALIE spectra render this analysis a very promising path toward expanding the database of spectroscopic binaries and multiple systems (like HD 41004). It may also facilitate the detection of planets in binary stellar systems, which are lately the focus of an increasing interest.

###### Acknowledgements.

This research was supported by the Israeli Science Foundation (grant no. 40/00). Support from Fundação para a Ciência e Tecnologia, Portugal, to N.C.S. in the form of a scholarship is gratefully acknowledged.


format: text

text:
# Doping dependence of the mass enhancement in (Pb,Bi)2Sr2CaCu2O8 at the antinodal point in the superconducting and normal state

(March 19, 2024)

## Abstract

Angle-resolved photoemission spectroscopy (ARPES) is used to study the mass renormalization of the charge carriers in the high-Tc superconductor (Pb,Bi)2Sr2CaCu2O8 in the vicinity of the $`(\pi ,0)`$ point in the superconducting and the normal state. Using matrix element effects at different photon energies and due to a high momentum and energy resolution the bonding and the antibonding bands could be separated in the whole dopant range. A huge anisotropic coupling to a bosonic collective mode is observed below Tc for both bands in particular for the underdoped case. Above Tc, the more isotropic coupling to a continuum or a mode at much higher energy is significantly weaker.

Sixteen years after the discovery of high-Tc superconductors (HTSCs) by Bednorz and M ller Muller  the understanding of the electronic transport in these systems is still full of unresolved issues. This holds both for the normal and for the superconducting state. In particular, it is not clear whether in these compounds which are close to an antiferromagnetic insulator the interaction between the electrons is large enough to destroy coherence, i.e., the existence of well defined quasi-particles, and to result in a new state of matter Anderson . On the other hand, one can possibly analyze the available experimental data in terms of a strong coupling to bosons such as phonons or spinfluctuations, leading to coherent and incoherent states. Moreover, there are intermediate scenarios Varma . Therefore it is interesting to investigate in which dopant and temperature range coherent states are present.

Traditional angle-resolved photoemission spectroscopy (ARPES) has been an important experimental probe to map out the bandstructure of solids. Recent improvements of the energy and the momentum resolution opened the possibility to study the many-body properties of solids. The deviation from the independent particle dispersion close to the Fermi level gives information on the renormalization of bands or on the enhanced effective mass. The width of the peaks reflects the scattering rate or the inverse lifetime of the photohole which should behave in a similar way as the charge carriers in a p-type doped HTSC. Both, the mass renormalization and the finite width of the spectral weight are normally described in terms of the real and the imaginary part of the self energy $`\mathrm{\Sigma }`$, respectively.

Numerous ARPES studies on HTSCs Damascelli ; Campuz  have shown that the renormalization is strongly momentum dependent, i.e., it is a function of the position on the Fermi surface. Up to now there are many ARPES studies along the nodal direction (the $`\mathrm{\Gamma }(\pi ,\pi )`$ line) where in the superconducting state no gap is opened. All these studies indicate that there is a mass renormalization m\*/m=1+$`\lambda `$ at the Fermi level by a factor of about 2 corresponding to a coupling constant $`\lambda `$ of about 1. The renormalization extends to an energy of about 70 meV below the Fermi level indicated by a kink in the dispersion. Currently there is no consensus on the origin of this renormalization. Coupling to collective bosonic modes such as phonons or spin fluctuations are discussed.

At the antinodal point (the $`(\pi ,0)(\pi ,\pi )`$ Fermi surface crossing) where the superconducting gap has a maximum, only a few experimental studies exist Damascelli ; Campuz ; Kaminski ; Dessau . Moreover, the situation is more complicated by the fact that in the drosophila of the ARPES studies on HTSCs, Bi2Sr2CaCu2O8 (Bi2212), the so-called bilayer splitting into a bonding (B) and an antibonding (A) band appears in this region of the Brillouin zone due to an interaction between two adjacent CuO2 layers. This bilayer splitting could not be resolved in many previous low resolution experiments and therefore the existence of coherent states could be overlooked. On the other hand, the investigation of bilayer systems offers the unique possibility to study the interaction of CuO2 layers which is of extreme importance for the understanding of the mechanism of high-Tc superconductivity, since Tc increases with the number of interacting layers up to 3 and then decreases again.

In this Letter we present high resolution ARPES data of Bi2212 around the $`(\pi ,0)`$ point for the whole dopant range from the overdoped (OD) via optimally doped (OP) to underdoped (UD) samples. Contrary to previous experiments Kaminski , we find dispersing renormalized bands not only in the superconducting state but also in the normal state in the whole dopant range, even for UD samples. Since we used the photon energy dependence of the excitations from the B and A bands we could determine the respective renormalization of both bands. It turned out that above Tc the renormalization is rather isotropic and weaker than in the superconducting state. No kink and no specific energy scale could be detected. In the superconducting state the renormalization due to a coupling to a collective bosonic mode is very strong and increases with decreasing dopant concentration. The results are in contrast to the traditional view, in which due to correlation effects in the normal state of UD samples no dispersing states exist and in the superconducting state the number of coherent states are expected to increase. Thus, the present results indicate a continuous transition from a more conventional metal to a strange metal with a small number of coherent states and not a transition into a new state of matter when going from the OD to the UD region.

The ARPES experiments were carried out at the BESSY synchrotron radiation facility using the U125/1-PGM beam line and a SCIENTA SES100 analyser. Spectra were taken with various photon energies ranging from 17 to 65 eV. The total energy resolution ranged from 8 meV (FWHM) at photon energies h$`\nu `$=17-25 eV to 22.5 meV at h$`\nu `$=65 eV. The momentum resolution was set to 0.01 Å-1 parallel to the $`(\pi ,0)(\pi ,\pi )`$ direction and 0.02 Å-1 parallel to the $`\mathrm{\Gamma }(\pi ,0)`$ direction. Here we focus on spectra taken with photon energies of 38 eV and 50 (or 55 eV) to discriminate between B and A bands. The polarization of the radiation was along the $`\mathrm{\Gamma }(\pi ,0)`$ direction. Measurements have been performed on (1$`\times `$5) superstructure-free Pb-Bi2212 single crystals ranging from OD (OD69, Tc= 69 K) via OP (OP89K, Tc= 89 K) to UD(UD77, Tc= 77 K) samples. The dopant concentration was calculated from the empirical equation of Tc vs. the dopant concentration Tallon .

In Fig. 1 we show a collection of data for wave vectors close to the $`(\pi ,\pi )(\pi ,\pi )`$ line, centered around the $`(\pi ,0)`$ point. As shown previously Kordyuk ; Borisenko  and supported by theoretical calculations Bansil , due to matrix element effects the data taken with h$`\nu `$= 38 eV represent mainly the B band with some contributions from the antibonding band while the data taken with h$`\nu `$= 50 (or 55) eV have almost pure A character. The subtraction of the latter from the h$`\nu `$= 38 eV data represents almost the pure B bands. This was performed by scaling the h$`\nu `$=50 eV spectral weight in such a way that at the $`(\pi ,0)`$ point no negative intensity occurs. Using this procedure one clearly recognizes that in the whole doping range and even in the UD samples the B and the A band can be well separated.

In the superconducting state these data show strong changes upon reducing the dopant concentration: the B and most clearly seen the A band moves further and further below the Fermi level and the renormalization strongly increases as indicated by the appearance of a flat dispersion below the gap energy. At $``$70 meV for the B band a kink and a reduction of spectral weight is observed, both increasing with decreasing dopant concentration.

In Fig. 2 we show on an enlarged scale data for the B and A band of the UD77 sample. Also included are dispersions derived from constant energy scans or momentum distribution curves (MDCs) in red for the B band and from constant momentum scans or energy distribution curves (EDCs) in blue for both bands. As previously discussed Campuz ; Dessau  the EDC and the MDC dispersions are quite different for weakly dispersing states, i.e., in the superconducting state, while they are similar in the normal state. Below Tc the EDC derived dispersion shows the typical BCS-type back-dispersion due to particle-hole mixing. At 120 K a pseudogap is observed and there is no observable back-dispersion. This could explain the absence of a coherent peak in the tunnelling spectroscopy data in the pseudogap regime Fischer .

The EDC depicted in Fig. 2, which was measured in the superconducting state at kF of the B band, is typical of a strong coupling of the charge carriers to a bosonic mode Damascelli ; Campuz ; Dessau . For this sample, the dip-energy, $`E_D=\mathrm{\Delta }+E_M`$, occurs at 70 meV which with a superconducting gap, $`\mathrm{\Delta }`$= 30 meV, corresponds to a mode energy $`E_M`$ of 40 meV Eschrig . Then the dispersion E(k) in the vicinity of the gap can be approximated by the expression E(k) =$`(ϵ_k^2+\mathrm{\Delta }^2)^{1/2}`$ where $`ϵ_k`$ is a renormalized band dispersion which is approximated by a linear dispersion $`ϵ_k=v_F^b(kk_F)/(1+\lambda )`$ where $`v_F^b`$ is the Fermi velocity of the bare particles and $`k_F`$ is the Fermi wave vector Schrieffer . The bare Fermi velocities were taken from a previous evaluation of the Fermi surface and data along the nodal direction in terms of a tight-binding bandstructure (also shown in Fig. 2) Kordyuk\_gap\_ARPES . This bandstructure is close to that derived from an analysis of the anisotropic plasmon dispersion Nucker  and to that obtained from LDA bandstructure calculations Andersen . The values for $`\lambda `$, presented in Fig. 3, and $`\mathrm{\Delta }`$ values (UD: 30 meV; OP: 22 meV; OD: 16 meV) were obtained by fits of E(k) to the EDC dispersions. In OP and OD samples, where steeper dispersions are observed, similar values of $`\lambda `$ were also derived from MDC dispersions. For the A band reliable values for $`\lambda `$ could be only obtained for an UD sample since this band is just crossing the Fermi level.

With decreasing dopant concentration there is a strong increase of the coupling strength from $`\lambda `$ of about 3 to huge values of about 8. We emphasize that the increase seems to be continuous and there is no break between the OD and the UD region. For the A band in the UD77 sample a similar huge mass renormalization has been detected (see Fig. 3). The results for the A and the B band in the superconducting state can also be summarized in a different way. With decreasing dopant concentration the coherence factor Z=1/(1+$`\lambda `$) decreases from 0.25 to 0.1. This means that for UD samples below Tc only about 1/9 of the spectral weight near EF represent coherent states.

Above Tc we see dispersive states even in UD samples and the renormalization of both bands is strongly reduced (see Fig. 2). There is also no kink in the MDC dispersion at $``$70 meV due to a coupling to a bosonic mode. The mode is either at much higher energies or there is a continuum which leads to the renormalization. The coupling constant as derived from the differences in the slopes of the measured dispersion and the bare particle dispersion is close to one, independent of the dopant concentration. These values are also depicted in Fig. 3. Comparing these values with similar $`\lambda `$ values at the nodal point indicates that the renormalization in the normal state is rather isotropic along the Fermi surface.

To obtain further information on the renormalization of the dynamics of the charge carriers in the normal state, we have evaluated the MDC width of the A and the B bands which multiplied by $`v_F^b`$, is a measure of the imaginary part of $`\mathrm{\Sigma }`$ plus contributions due to the finite energy and momentum resolution. Typical MDC cuts measured with a photon energy of 38 eV are shown in Fig. 4 for UD and OD samples. These cuts can be well fitted by 4 Lorentzians, 2 corresponding to the B and 2 others corresponding to the A band. For the OD and the UD samples the derived scattering rates at EF are almost the same for both the B and the A bands. Values between 80 and 160 meV are derived for the antinodal point at T= 120 K. These values are not far from our values ($``$ 100 meV) derived at the same temperature and at the same energy at the nodal point. Moreover, at energies smaller than about 100 meV the scattering rates are only slightly dependent on energy (see Fig. 4). This is an indication that the modes to which the charge carriers are coupled, are at higher energies or form a continuous spectrum.

The results obtained can be summarized in the following way: in the normal state a rather isotropic mass renormalization of the electronic states without a clear energy scale is observed. The strength of the coupling corresponds to $`\lambda `$ $``$1. Below Tc and in particular in the UD samples a much higher strongly anisotropic coupling to a collective bosonic mode could be detected. We emphasize that the difference in the mass renormalization at the nodal point between the superconducting and the normal state is rather small (less than 5 %) .

It is difficult to interpret the entirety of our results in terms of a conventional isotropic coupling to phonons. A probably more promising model would be a coupling to spin fluctuations Chubukov . Although this scenario has been applied in many previous ARPES studies the present work provides a much more detailed picture. Above Tc the system starts from charge carriers which are coupled to a continuum of spin fluctuations. The strength of the coupling is not strongly dopant dependent. Both the A and the B band feel a similar coupling as derived from the mass renormalization and from the scattering rates. Previous ARPES studies at the $`(\pi ,0)`$ point could not resolve the bilayer splitting and could not follow the flat dispersion of bands and therefore came to the conclusion that in the UD region above Tc there are only incoherent states. In the present study, renormalized dispersive and possibly coherent states are even detected above Tc in UD samples.

In the spin fluctuation scenario, below Tc, the opening of the gap leads via a feed-back process to a magnetic resonance mode at EM detected by inelastic neutron scattering neutron  to which the charge carriers couple. Due to the fact that the bilayer splitting of the band could be resolved, a quantitative analysis of the coupling strength could be performed for the B band in the entire doping range and for the A band for an UD sample. In a recent theoretical work Eschrig  it was pointed out that according to magnetic susceptibility measurements using inelastic neutron scattering the magnetic resonance mode couples the A band only to the B band and vice versa. There is no coupling via the resonance mode within a band. In that paper Eschrig  the authors came to the conclusion that in the OD region the mass enhancement of the B band should be large since the van Hove singularity of the A band comes close to the Fermi level and leads to a large susceptibility $`\chi _{AB}`$. This is at variance with our experimental results where $`\lambda `$ for the B band does not show a maximum in the OD region. On the other hand, it is remarkable that the coupling of the B band to the resonance mode starts when the A band crosses the Fermi level (in the OD region). Moreover, the result that in the UD region $`\lambda `$ is similar for the both bands is understandable, since the Fermi velocities and therefore the density of states and the odd susceptibilities $`\chi _{AB}`$ and $`\chi _{BA}`$ should be comparable.

This work was financially supported in part by the Schweizerische Nationalfonds zur Foerderung der Wissenschaftlichen Forschung.


text:
# An Estimate of the Vibrational Frequencies of Spherical Virus Particles

## I Introduction

Virus particles (virions) come in a variety of sizes and shapes. However, approximately spherical shapes with diameters in the range between $`50nm`$ and $`100nm`$ are especially common. Many nearly spherical viruses are revealed by X-ray crystallography to have icosahedral symmetry. A typical virus particle contains genetic material, RNA or DNA, surrounded by a protein coat (capsid). Such an object should have reasonably distinct vibrational frequencies, the study of which may be of interest. Excitation of these vibrations could have applications in either the diagnosis or treatment of viral diseases. To this author’s knowledge, the sole discussion of these vibrational modes in the literature is that of Babincová et al Bab . These authors discuss the conjecture that ultrasound in the GHz range could be resonantly absorbed by HIV virus particles, leading to their destruction. Cooper et al Cooper  have recently reported the detection of viruses by acoustic oscillations. However, the process of “rupture event scanning”, which these authors report, involves the separation of a virus particle from antibodies by ultrasound. This is distinct from the excitation of the vibrational modes of the virus particle itself, and occurs at much lower frequencies.

There have also been some experimental studies of ultrasonic absorption by empty viral capsids Cerf79 ; Michels85 . These experiments reveal an enhanced absorption in the MHz range as proteins reassemble into a capsid, but do not find a resonant peak in this frequency range. Witz and Brown WB01  have emphasized that these and other results show that viral capsids are flexible and change size or shape in response to vibrations or to changes in temperature or pH.

The purpose of the present paper is to provide some estimates of the lowest vibrational frequencies of a spherical virus particle. The simplest estimate is to take this frequency to be of the order of a characteristic speed of sound divided by the size of the virus particle. This is the estimate used in Ref. Bab . For the purpose of giving a more accurate estimate, we will examine two models, which treat the particle (1) as a liquid drop and (2) as a uniform elastic sphere. Similar models have been used by Bulatov et al Bulatov  to estimate the vibrational frequencies of nanoclusters.

## II A Liquid Drop Model

Consider a sphere of radius $`a`$ filled with a nonviscous liquid with surface tension $`\gamma `$ and mass density $`\rho `$. The lowest vibrational mode of this sphere will be a quadrupole mode with frequency Chandra

$$\nu =\frac{1}{\pi }\sqrt{\frac{2\gamma }{\rho a^3}},$$
(1)
which can be written as

$$\nu =3.4\times 10^8Hz\left(\frac{50nm}{a}\right)^{\frac{3}{2}}\left(\frac{\rho _W}{\rho }\right)^{\frac{1}{2}}\left(\frac{\gamma }{\gamma _W}\right)^{\frac{1}{2}},$$
(2)
where $`\rho _W=10^3kg/m^3`$ and $`\gamma _W=0.073Nm`$ are the mass density and surface tension for water, respectively. The surface tension and mass density, along with the lowest vibrational frequency derived from Eq. (2) for $`a=50nm`$, are given in Table 1 for several liquids.

Recall that Eq. (2) assumes a nonviscous liquid. In fact, the viscosity of many of the liquids in Table 1 cannot be neglected for such a small drop. When viscosity is sufficiently large, the drop will not oscillate, but rather undergo overdamped motion Chandra . The main lesson from the liquid drop model is that a nonviscous liquid drop of $`a=50nm`$ with a typical surface tension and mass density would have a lowest vibrational frequency of the order of a few times $`10^8Hz`$.

## III An Elastic Sphere Model

A better model for a virus particle is to treat it as a uniform elastic sphere. The three independent parameters which characterize such a sphere can be taken to be the radius $`a`$, the speed of pressure waves, $`c_P`$, and the speed of shear waves, $`c_S`$. The oscillations of an elastic sphere are treated in detail by Pao and Mow PM . Here we quote their results in the notation of Ref. Bulatov . The eigenfrequencies of the normal modes are given by the vanishing of the determinant

$$\left|\begin{array}{cc}S_{rp}(n,x,y)\hfill & S_{rs}(n,y)\hfill \\ S_{tp}(n,x)\hfill & S_{ts}(n,y)\hfill \end{array}\right|=0,$$
(3)
where

$`S_{rp}(n,x,y)=(n^2n{\displaystyle \frac{1}{2}}y^2)^2j_n(x)+2xj_{n+1}(x),`$ (4)
$`S_{rs}(n,y)=n(n+1)[(n1)j_n(y)xj_{n+1}(y)],`$ (5)
$`S_{tp}(n,x)=(n1)j_n(x)xj_{n+1}(x),`$ (6)
$`S_{ts}(n,y)=(n^2n{\displaystyle \frac{1}{2}}x^2)^2j_n(x)xj_{n+1}(x).`$ (7)
Here $`x=\omega a/c_P`$, $`y=\omega a/c_S`$, and the $`j_n`$ are spherical Bessel functions. Let

$$b=\frac{c_S}{c_P}$$
(8)
be the ratio of the speed of the shear wave to that of the pressure wave, so that $`y=bx`$. Given $`b`$ and $`n`$, we can solve Eq. (3) for $`x`$ and hence for the frequency of the normal mode, $`\nu =xc_P/(2\pi a)`$. The smallest roots for $`x`$ for the $`n=0`$ and $`n=2`$ modes are plotted in Fig. 1. The corresponding frequencies of oscillation can be expressed as

$$\nu =4.8\times 10^9Hz\left(\frac{50nm}{a}\right)\left(\frac{c_P}{1500m/s}\right)x.$$
(9)
The frequencies obtained from this equation for the $`n=0`$ mode are given in Table II for various materials.

We can see that for a wide range of materials, the lowest mode of vibration is a purely radial mode with a frequency of the order of a few times $`10^9Hz`$ for a sphere of radius $`a=50nm`$.

## IV Conclusions

In the previous sections, we have examined two models for a spherical virus particle, a liquid drop model and an elastic sphere model. It is of interest that the two models yield estimates for the lowest vibrational frequency which differ by only about one order of magnitude. Of the two models, the elastic sphere model is probably the better description of a virus particle. An even better model might be one in which the particle has a liquid core (DNA or RNA) surrounded by an elastic outer shell (the capsid). Such a model would probably yield vibrational frequencies intermediate between those predicted by the two models discussed in this paper. In any case, we obtain an estimate for the lowest vibrational frequency of the same order of magnitude as that given in Ref. Bab , in the range of a few GHz for particles with a size of about $`100nm`$. Of course, the existence of a resonance requires that damping be below the critical value above which overdamped motion occurs. Even if this condition is fulfilled, it is difficult to predict the width of the resonance. This remains a question for experimental investigation. The existence of well defined resonances could prove valuable both for basic science and for medicine. Thus this is a potentially fruitful area for further research.

###### Acknowledgements.

I would like to thank V. Ford and D. weaver for helpful comments. This work was supported in part by the National Science Foundation under Grant PHY-9800965.


text:
# Long-term Evolution of Stellar Self-Gravitating System away from the Thermal Equilibrium : connection with non-extensive statistics

## Abstract

With particular attention to the recently postulated introduction of a non-extensive generalization of Boltzmann-Gibbs statistics, we study the long-term stellar dynamical evolution of self-gravitating systems on timescales much longer than the two-body relaxation time. In a self-gravitating $`N`$-body system confined in an adiabatic wall, we show that the quasi-equilibrium sequence arising from the Tsallis entropy, so-called stellar polytropes, plays an important role in characterizing the transient states away from the Boltzmann-Gibbs equilibrium state.

Introduction.— The long-term evolution of self-gravitating many-body system is an old problem with a rich history in astrophysics. The problem, in nature, involves the long-range nature of attractive gravity and is fundamentally connected with statistical mechanics and thermodynamics. Historically, the important consequence from the thermodynamical arguments has arisen in 1960s, known as the gravothermal catastrophe, i.e., thermodynamic instability due to the negative specific heatBinney & Tremaine (1987). Originally, the gravothermal catastrophe has been investigated in a very idealized situation, i.e., a stellar system confined in a spherical cavityAntonov (1962). Owing to the maximum entropy principle, the existence of an unstable thermal state has been found from the standard analysis of statistical mechanics with a particular attention to the Boltzmann-Gibbs(BG) entropy:

$$S_{\mathrm{BG}}=d^3𝒙d^3𝒗f(𝒙,𝒗)\mathrm{ln}f(𝒙,𝒗),$$
(1)
where $`f(𝒙,𝒗)`$ denotes the one-particle distribution function defined in phase-space $`(𝒙,𝒗)`$.

Since 1960s, the standard approach using the BG entropy has dramatically improved our view of the late-time phase of the globular cluster as a real astronomical systemMeylan & Heggie (1997), however, from a thermodynamic point of view, peculiarities of the thermal property in self-gravitating systems such as negative specific heat, as well as the non-equilibrium properties away from the BG state have not yet been understood completely.

In this Letter, aiming at a better understanding of the (non-equilibrium) thermodynamic properties of stellar self-gravitating systems, we present a set of long-term $`N`$-body simulations, the timescale of which is much longer than the relaxation time. With a particular emphasis to the recent application of the non-extensive generalization of BG statistics, we focus on the stellar dynamical evolution in an isolated star cluster before self-similar core-collapseCohn (1980). We show that the quasi-equilibrium sequence arising from the Tsallis entropyTsallis (1988) plays an important role in characterizing the non-equilibrium evolution of a self-gravitating system.

$`N`$-body simulations.— The $`N`$-body experiment considered here is the same situation as investigated in classic papers (Antonov (1962), see also Ref.Endo, Fukushige & Makino (1997)). That is, we confine the $`N`$ particles interacted via Newton gravity in a spherical adiabatic wall, which reverses the radial components of the velocity if the particle reaches the wall. Without loss of generality, we set the units as $`G=M=r_e=1`$, where $`G`$ is gravitational constant, $`M`$ and $`r_e`$ is the total mass of the system and the radius of the adiabatic wall, respectively. Note that the typical timescales appearing in this system are the free-fall time, $`T_{\mathrm{ff}}=(G\rho )^{1/2}`$, and the global relaxation time driven by the two-body encounter, $`T_{\mathrm{relax}}=(0.1N/\mathrm{ln}N)T_{\mathrm{ff}}`$ Binney & Tremaine (1987), which are basically scaled as $`T_{\mathrm{ff}}1`$ and $`T_{\mathrm{relax}}0.1N/\mathrm{ln}N`$ in our units. To perform an expensive N-body calculation, we used a special-purpose hardware, GRAPE-6, which is especially designed to accelerate the gravitational force calculations for collisional $`N`$-body systemsMakino (2002). With this implementation, the $`4`$th-order Hermite integrator with individual time stepMakino & Aarseth (1992) can work efficiently, which is suited for probing the relaxation process in denser core regions with an appropriate accuracy. We adopt the Plummer softened potential, $`\varphi =1/\sqrt{r^2+ϵ^2}`$ with a softening length $`ϵ`$ of $`1/512`$ and $`1/2048`$.

In the present numerical simulation, the choice of the initial condition is an important step. Here, we set the initial distribution to the stationary state of Poisson-Vlassov equation, i.e., dynamical equilibrium for a spherical system with isotropic velocity distribution. According to the Jeans theoremBinney & Tremaine (1987), the one-particle distribution function $`f(𝒙,𝒗)`$ can be expressed as a function of specific energy, $`\epsilon =v^2/2+\mathrm{\Phi }(r)`$ with $`r`$ and $`\mathrm{\Phi }`$ being the radius and the gravitational potential. Then keeping the energy and the mass constant, the thermal equilibrium of ordinary extensive statistics derived from the maximum entropy principle of the BG entropy (1) reduces to the exponential distribution, so-called isothermal distribution given by $`f(\epsilon )e^{\beta \epsilon }`$, which effectively satisfies the equation of state, $`P(r)\rho (r)`$, where $`P(r)`$ is pressure and $`\rho (r)`$ is mass densityBinney & Tremaine (1987); Antonov (1962).

On the other hand, as another possibility, one considers the extremum state of Tsallis’ non-extensive entropy Tsallis (1988):

$$S_q=d^3𝒙d^3𝒗[\{f(𝒙,𝒗)\}^qf(𝒙,𝒗)]/(1q),$$
(2)
which might be of particular importance in describing the quasi-equilibrium state away from the BG stateTsallis (1999). In this case, the maximum entropy principle leads to the power-law distribution, $`f(\epsilon )(\mathrm{\Phi }_0\epsilon )^{1/(q1)}`$ Plastino & Plastino (1993); Taruya & Sakagami (2002a, b, c), referred to as the stellar polytropeBinney & Tremaine (1987). It satisfies the polytropic equation of state, $`P(r)\rho (r)^{1+1/n}`$, and the polytrope index $`n`$ is related to the $`q`$-parameter as $`n=1/(q1)+3/2`$notice\_1 . Provided the polytrope index $`n`$, the equilibrium structure can be determined by solving the Láne-Emden equation Chandrasekhar (1939) and using this solution, the relationship between the dimensionless energy $`\lambda r_eE/(GM^2)`$ and the density contrast $`D\rho _c/\rho _e`$, the core density divided by the edge density, can be drawn (see Fig.1, notice\_2 ). Note that the limit $`n\mathrm{}`$ (or $`q1`$) corresponds to the isothermal distribution derived from the BG entropy (1).

Table 1 summarizes the list of the five simulation runs. A more systematic study of the systems with several initial conditions is now in progress and the details of the results will be reported elsewhere. In Table 1, in addition to the stellar polytropic initial state, we also consider the Hernquist modelHernquist (1990), which was originally introduced to account for the empirical law of observed elliptical galaxiesBinney & Tremaine (1987).

Results.— It has been recently discussed in Taruya & Sakagami (2002a, b, c) that the thermodynamic structure of the stellar polytropic distribution can be consistently characterized by the non-extensive framework of the thermostatistics. According to their results, the stellar polytrope confined in an adiabatic wall is shown to be thermodynamically stable when the polytrope index $`n<5`$. In other words, if $`n>5`$, a stable equilibrium state ceases to exist for a sufficiently large density contrast $`D>D_{\mathrm{crit}}`$, where the critical value $`D_{\mathrm{crit}}`$ given by a function of $`n`$ is determined from the second variation of entropy around the extremum state of Tsallis entropy, $`\delta ^2S_q=0`$Taruya & Sakagami (2002a, c). The dotted line in Fig.1 represents the critical value $`D_{\mathrm{crit}}`$ for each polytrope index, which indicates that the stellar polytrope at low density contrast $`D<D_{\mathrm{crit}}`$ is expected to remain stable. Apart from the BG state, one might expect that the stellar polytrope acts as a thermal equilibrium.

Of course, this naive expectation is not correct at all. Indeed, the numerical simulations reveal that the stellar polytropic distribution gradually changes with time, on the timescale of two-body relaxation. Further, it seems that the gravothermal instability appears at relatively smaller values of $`D`$ than the predicted one, $`D_{\mathrm{crit}}`$, which might be partially ascribed to the non-stationarity of the background stellar polytrope. Physically, the core-collapse phenomenon due to the gravothermal catastrophe follows from the decoupling of the relaxation timescales between the central and the outer parts, whose behavior sensitively depends on the physical property of heat transportMakino & Hut (1991). In a rigorous sense, the thermodynamic prediction might lose the physical relevance, however, focusing on the evolutionary sequence, we found that the transient state starting from the initial stellar polytrope can be remarkably characterized by a sequence of stellar polytropes (run A, B1 and B2). This is even true in the case starting from the Hernquist model (run C1).

Let us show the representative results taken from run A (Fig.2). Fig.2(a) plots the snapshots of the density profile $`\rho (r)`$, while Fig.2(b) represents the distribution function $`f(\epsilon )`$ as a function of the specific energy $`\epsilon `$. Note that just for illustrative purpose, each output result is artificially shifting to the two-digits below. Only the final output with $`T=400`$ represents the correct scales. In each figure, solid lines mean the initial stellar polytrope with $`n=3`$ and the other lines indicate the fitting results to the stellar polytrope by varying the polytrope index $`n`$notice\_3 . Note that the number of fitting parameters just reduces to one, i.e., the polytrope index, since the total energy is well-conserved in the present situation. Fig.2 shows that while the system gradually deviates from the initial polytropic state, the transient state still follows a sequence of stellar polytropes. The fitting results are remarkably good until the time exceeds $`T400`$, corresponding to $`15T_{\mathrm{relax}}`$. Afterwards, the system enters the gravothermally unstable regime and finally undergoes the core-collapse.

Now, focus on the evolutionary track in each simulation run summarized in the energy-density contrast plane (Fig.1), where the filled circle represents the initial stellar polytrope. Interestingly, the density contrast of the transient state in run A initially decreases, but it eventually turns to increase. The turning-point roughly corresponds to the stellar polytrope with index $`n56`$. Note, however, that the time evolution of polytrope index itself is a monotonically increasing function of time as shown in Fig.2(c), apart from the tiny fluctuations. This is indeed true for the other cases, indicating the Boltzmann $`H`$-theorem that any of the self-gravitating systems tends to approach the BG state. A typical example is the run C2, which finally reaches the stable BG state. However, as already shown in run A, all the systems cannot reach the BG state. Fig.1 indicates that no BG state is possible for a fixed value $`\lambda >0.335`$Antonov (1962), which can be derived from the peak value of the trajectory. Further, stable stellar polytropes cease to exist at high density contrast $`D>D_{\mathrm{crit}}`$. In fact, our simulations starting from the stellar polytropes finally underwent core-collapse (runs A, B1 and B2). Though it might not be rigorously correct, the predicted value $`D_{\mathrm{crit}}`$ provides a crude approximation to the boundary between the stability and the instability.

Fig.3 plots the snapshots of the distribution function taken from the other runs. The initial density contrast in run B1 (Fig.3(a)) is relatively low ($`D=110`$) and thereby the system slowly evolves following a sequence of stellar polytropes. After $`T=200074T_{\mathrm{relax}}`$, the system begins to deviate from the stable equilibrium sequence, leading to the core-collapse. Another noticeable case is the run C1 (Fig.3(b)). The Hernquist model as initial distribution of run C has cuspy density profile, $`\rho (r)1/r/(r+a)^3`$, which behaves as $`\rho r^1`$ at the inner part Hernquist (1990). The resultant distribution function $`f(\epsilon )`$ shows a singular behavior at the negative energy region, which cannot be described by the power-law distribution. Soon after a while, however, the gravothermal expansionEndo, Fukushige & Makino (1997) takes place and the flatter core is eventually formed. Then the system settles into a sequence of stellar polytropes and can be approximately described by the stellar polytrope with index $`n=20`$ for a long time. Thus, as long as the stellar system confined in an adiabatic wall is concerned, the stellar polytrope can be regarded as a quasi-attractor and a quasi-equilibrium state.

Of course, these remarkable features could be an outcome in a very idealized situation and one suspects that quasi-equilibrium state of stellar polytrope cannot hold if we remove the boundary of the adiabatic wall. As a demonstration, Fig.3(c) plots the results removing the boundary, in which the initial state is the same distribution as in run A. As is expected, the high-energy particles freely escape outwards from the central region and the resultant distribution function $`f(\epsilon )`$ sharply falls off at the energy region $`\epsilon 0`$, indicating that the density contrast $`D`$ becomes effectively large. Thus, compared to the system confined in the wall, the removal of the boundary makes the stellar system unstable and the core-collapse takes place earlier. Nevertheless, focusing on the inner part of the denser region, the evolution of the core is not significantly affected by the escape particles at the outer part and can be fitted by a sequence of stellar polytropes (see also the dashed line in Fig.2(c)). The successful fit to the density profile $`\rho (r)`$ almost remains the same. Hence, the stellar polytrope characterized by the Tsallis entropy can be even realized in a realistic situation removing the boundary of the adiabatic wall.

Summary & Discussions.— We have performed a set of numerical simulation of long-term stellar dynamical evolution away from the BG state and found that the transient state of the system confined in an adiabatic wall can be remarkably fitted by a sequence of stellar polytropes. This is even true in the case removing the outer boundary. Therefore, the stellar polytropic distribution can be a quasi-attractor and a quasi-equilibrium state of a self-gravitating system.

Alternative characterization of the transients away from the BG state might be possible besides the $`q`$-exponential distribution of stellar polytropes. For an empirical characterization of observed structure, the one-parameter family of truncated exponential distributions, so-called King model has been used in the literatureBinney & Tremaine (1987); Meylan & Heggie (1997); King (1966). Also, the sequence of King model has been found to characterize the evolutionary sequence of density profile for isolated stellar systems without boundaryCohn (1980). We have also tried to fit the simulation data to the King model. Similarly to the stellar polytrope, the King model accurately describes the simulated density profile $`\rho (r)`$ confined in an adiabatic wall, however, it fails to match the simulated distribution function $`f(\epsilon )`$, especially at the cutoff energy scales. Therefore, from the quantitative description of the entire phase-space structure, the power-law distribution of the stellar polytropes can be a better characterization of the quasi-equilibrium state and this could yield an interesting explanation of the origin of the empirical King model.

We are grateful to T. Fukushige for providing us the GRAPE-6 code and for his constant supports and helpful comments. We also thank to J. Makino for fruitful discussion especially on the applicability of equilibrium sequences of stellar polytrope. Numerical computations were carried out at ADAC (the Astronomical Data Analysis Center) of the National Astronomical Observatory of Japan. This research was supported in part by the grant-in-aid from Japan Society of Promotion of Science (No.1470157).


text:
# Collective dynamics of fermion clouds in cigar-shaped traps

## 1 Introduction

Current experiments on the collisional properties of a Fermi gas in two spin-polarized states under confinement have proven that this setup provides an important tool for investigating the dynamics of dilute quantum gases. For such a system in harmonic trapping the existence of zero and first sound excitations has been demonstrated experimentally  and numerically  by measuring the damping of the oscillations of a non-equilibrium density profile as a function of the collision rate. These studies have followed the evolution of the dynamical behaviour of the quantum gas from a collisional (first sound) regime to a collisionless (zero sound) regime. However, no direct measurement of the speed of sound propagation has been carried out so far in confined Fermi gases, contrary to the case of Bose-Einstein condensates in elongated traps  where the propagation of density distortions has been followed experimentally.

The main purpose of the present study is to evaluate by both theoretical and numerical means the propagation of density distortions in the zero-sound regime for a degenerate Fermi gas contained inside a cigar-shaped harmonic trap, in situations that are accessible to experiment. It may be useful to recall at this point the earlier studies of these collective excitations in homogeneous Fermi fluids. The experiments of Abel et al.  on acoustic waves propagating through liquid 3He demonstrated that the transition from collisional to collisionless dynamics is signalled by a drastic change in the frequency dependence of sound-wave attenuation and showed that the difference in the speeds of first and zero sound is quite small. This is as expected for a strongly coupled Fermi liquid from the Landau theory . The evolution of the zero-sound dispersion relation in liquid 3He into the microscopic regime was determined by Sköld et al.  in inelastic neutron scattering experiments and accounted for theoretically by Aldrich et al.  within a polarization potential approach extending Landau’s Fermi-liquid theory to larger momenta and energies. Landau’s theory has also been used by Yip and Ho  to evaluate the collisionless excitation modes corresponding to density fluctuations and to coherent interconversion of various spin species in a homogeneous Fermi gas.

At variance from the strongly coupled 3He liquid, in a dilute Fermi gas the interactions can be expected to introduce only slight modifications to the bare sound velocities, which in the homogeneous noninteracting limit are $`v_F`$ and $`v_F/\sqrt{3}`$ for zero and first sound, respectively. Here, the Fermi velocity $`v_F`$ is determined by the Fermi wave number and hence by the particle number density. These well-known results  immediately suggest how the propagation of collective density fluctuations in a confined Fermi gas may be treated within a local-density approximation.

In this work we first present an analytical study of the behaviour of the Fermi velocity and of the speed of zero sound in a two-component inhomogeneous Fermi gas as functions of the mean-field interactions in the collisionless regime. We support our analytical estimate by numerical simulations solving the Vlasov-Landau equations (VLE) for the fermionic Wigner distributions at $`T=0.2T_F`$, $`T_F`$ being the Fermi temperature . We evaluate the propagation of zero sound by simulating an experiment where perturbation in the density profile is introduced by a laser beam permanently focused on the centre of the trap . We also simulate an experiment carried on by Dutton et al.  on a Bose-Einstein condensate, where a defect with a large gradient in the density is generated by cutting off a slice of the profile by means of an ultra-compressed slow pulse of light. In their experiment the generation of quantum shock waves in the condensate was observed.

The paper is organized as follows. In Sec. 2 we recall the equations of motion for a two-component Fermi gas and evaluate in a semiclassical local-density approach the Fermi velocity and the zero sound velocity for such a system. In Sec. 3 we report illustrations of our numerical simulations for a one-component and a two-component Fermi gas, and compare the results with our theoretical estimates. Conclusions and perspectives are given in Sec. 4.

## 2 The model

The two spin-polarized fermionic components of the gas are described by the Wigner distribution functions $`f^{(j)}(𝐫,𝐩,t)`$, with $`j=1`$ or 2. These obey the kinetic equations

$$\frac{f^{(j)}}{t}+\frac{𝐩}{m_j}f^{(j)}U^{(j)}_𝐩f^{(j)}=C_{12}[f^{(j)}].$$
(1)
At high dilution the interactions between atoms of the same spin are suppressed by the Pauli principle, so that the mean-field Hartree-Fock (HF) effective potentials are $`U^{(j)}(𝐫,t)=V_{ext}^{(j)}(𝐫,t)+gn^{(\overline{j})}(𝐫,t)`$ with $`\overline{j}`$ indicating the species different from $`j`$ and $`V_{ext}^{(j)}(𝐫,t)`$ being the sum of external confining and driving potentials. Here we have set $`\mathrm{}=1`$ and $`g=2\pi a/m_r`$, with $`a`$ the $`s`$-wave scattering length ($`a>0`$ for the present case of repulsive interactions) and $`m_r`$ the reduced mass. The total spatial densities $`n^{(j)}(𝐫,t)=n_0^{(j)}(𝐫)+\delta n^{(j)}(𝐫,t)`$ are the sum of the equilibrium densities $`n_0^{(j)}(𝐫)`$ and the time-dependent perturbations $`\delta n^{(j)}(𝐫,t)`$, and are obtained by integration of $`f^{(j)}(𝐫,𝐩,t)`$ in momentum space. Finally $`C_{12}[f^{(j)}]`$ is the collision integral, which will be set to zero in the collisionless regime. Equation (1) in this regime will be solved numerically in Sec. 3, using the algorithm developed by Toschi et al.  and taking the equilibrium density profiles as given in the HF approach  by

$$n_0^{(j)}(𝐫)=\frac{d^3p}{(2\pi )^3}\left\{\mathrm{exp}\left[\beta \left(\frac{p^2}{2m_j}+U^{(j)}(𝐫)\mu ^{(j)}\right)\right]+1\right\}^1,$$
(2)
where $`\beta =1/k_BT`$ and $`\mu ^{(j)}`$ is the chemical potential of species $`j`$.

In the rest of this section we present instead a theoretical discussion of the long-wavelength excitations of the gas. The perturbation $`\delta n^{(j)}(𝐫,t)`$ moves with a group velocity whose natural unit is the Fermi velocity $`v_F^{(j)}`$. In an inhomogeneneous gas the Fermi velocities depend on position and will be related below to the local chemical potential of each species through a local density approximation. Assuming for the moment that they are known, the dynamical response of the two-component Fermi gas can be calculated within a random phase approximation  from the set of eigenvalue equations

$$\delta n^{(j)}(q,\mathrm{\Omega })=\chi _j(q,\mathrm{\Omega })[\delta V^{(j)}(q,\mathrm{\Omega })+g\delta n^{(\overline{j})}(q,\mathrm{\Omega })],$$
(3)
where

$$\chi _j=𝒩_j\left[1\frac{\mathrm{\Omega }}{2v_F^{(j)}q}\mathrm{ln}\left(\frac{\mathrm{\Omega }+v_F^{(j)}q}{\mathrm{\Omega }v_F^{(j)}q}\right)\right]$$
(4)
with $`𝒩_j=v_F^{(j)}m_j^2/2\pi ^2`$ being the density of states of each component per unit volume at the Fermi level. In writing these equations we have specifically assumed that we are treating the case of long-wavelength excitations propagating with wave number $`q`$ along the axis of an elongated cylindrical trap, with the Fourier transform $`\delta V^{(j)}(q,\mathrm{\Omega })`$ of the driving potentials being independent of the radial coordinate and the Fourier transform $`\delta n^{(j)}(q,\mathrm{\Omega })`$ of the density distortions being taken at the value $`r=0`$ of the radial coordinate. Following  and  it is easily seen that the collective modes of the Fermi gas described by Eq. (3) are determined by the solutions of the equation

$$1g^2\chi _1(q,\mathrm{\Omega })\chi _2(q,\mathrm{\Omega })=0,$$
(5)
yielding in general two eigenfrequencies $`\mathrm{\Omega }_j(q)`$.

Let us turn next to the calculation of the Fermi velocities. Each of them is determined by the corresponding single-particle dispersion relation as $`v_F=E(k)/k|_{k=k_F}`$, yielding in a local density approximation

$$v_F^{(j)}=(2\mu _{\mathrm{loc}}^{(j)}/m_j)^{1/2}$$
(6)
where $`\mu _{\mathrm{loc}}^{(j)}(𝐫)=\mu ^{(j)}U^{(j)}(𝐫)`$ is the local chemical potential. In cylindrical harmonic confinement and in the absence of interactions $`v_F^{(j)}`$ can be written as

$$v_F^{(j)}(r,z)=\left[2\left(\mu ^{(j)}V^{(j)}(r,z)\right)/m_j\right]^{1/2}$$
(7)
where $`V^{(j)}(r,z)=m_j\omega _j^2(r^2+\epsilon _j^2z^2)/2`$, with radial trap frequency $`\omega _j`$ and anisotropy $`\epsilon _j`$. Under driving the radial and axial components of $`v_F^{(j)}`$ then obey harmonic equations of motion,

$$v_F^{(j)}(t)|_{z=0}=\left(2\mu ^{(j)}/m_j\right)^{1/2}\mathrm{sin}(\omega _jt+\varphi _r)$$
(8)
and

$$v_F^{(j)}(t)|_{r=0}=\left(2\mu ^{(j)}/m_j\right)^{1/2}\mathrm{sin}(\epsilon _j\omega _jt+\varphi _z),$$
(9)
where the phases $`\varphi _r`$ and $`\varphi _z`$ are determined by the initial value of the perturbation. In the presence of mean-field interactions, on the other hand, Eq. (7) is replaced by

$$v_F^{(j)}(r,z)=\left[2\left(\mu ^{(j)}V^{(j)}(r,z)gn_0^{(\overline{j})}(r,z)\right)/m_j\right]^{1/2}.$$
(10)
In this case harmonicity under driving will be lost, except when the mean field term can be treated perturbatively and for small axial displacements from to the center of the trap. In this regime it is easy to demonstrate that the mean-field interactions lower the Fermi velocity. More generally, their net effect is to decrease the chemical potential by inducing a decrease of the fermion density.

With the above as background, we conclude this section by presenting our estimates for the speed $`c_0`$ of collisionless sound waves propagating along the trap axis at $`r=0`$ in two special cases, that will be tested in the next section by comparison with the evolution of density profiles in time from the numerical solution of Eq. (1):

(i) A single-component Fermi gas. Setting $`g=0`$, Eqs. (3) and (4) immediately yield

$$c_o=v_F(0,z).$$
(11)
This situation will arise not only in the obvious case where the interspecies mean-field interactions vanish, but also in the case where they become so strong in a two-component Fermi gas that the regime of phase separation predicted by the theory  is being approached. In this latter case the residual mean-field interactions with the second species merely shift the chemical potential of the species that resides in the inner part of the trap and hence will shift its Fermi velocity.

(ii) A symmetric two-component Fermi gas. Taking equal masses and numbers of particles for the two fermionic species inside identical traps, we can write $`c_0=\lambda v_F`$ for the local speed of sound wave propagation and determine $`\lambda `$ from Eq. (5) using the Landau Fermi-liquid theory to describe the mean-field interactions through the appropriate Landau parameters. Only the spin up-down scattering is allowed in the s-wave channel for a two-component mixture of spin polarized fermions, so that the only non-zero interaction parameters are $`F_0^s`$ and $`F_0^a`$ with the condition $`F_0^s=F_0^a`$ ($`=F_0`$, say). Equation (5) can be written in the form

$$\frac{\lambda }{2}\mathrm{ln}\left(\frac{\lambda +1}{\lambda 1}\right)=1+\frac{1}{F_0}$$
(12)
with

$$F_0=g𝒩=2ak_F/\pi .$$
(13)
In Figure 1 we have plotted the behaviour of the solution of Eq. (12) as a function of the Fermi velocity for various values of the scattering length. We also show in Figure 2 the values of $`v_F`$ and $`c_0`$ as functions of the scattering length, for the choice of the other system parameters that will be adopted in the numerical calculations reported in the next section. It is seen from these figures that, in order to obtain a very significant increase of the speed of zero sound over the Fermi velocity from repulsive mean-field interactions, one would need to work on a high-density system having huge values of the scattering length. The value of $`a`$ can to some extent be tuned by exploiting Feshbach resonances , but strong repulsions will also favour demixing of the two-component Fermi gas.

## 3 Numerical results

### 3.1 Perturbation introduced by a permanent laser beam

As a first case we consider a one-component gas of 1000 spin-polarized 40K atoms inside an axially symmetric trap with $`\omega =100`$ s-1 and $`\epsilon =0.1`$. We use 200 computational particles to represent each fermion, in order to obtain a low noise-to-signal ratio . The equilibrium density profile is perturbed by an effective potential $`U_0\mathrm{exp}(z^2/2\sigma ^2)`$ with amplitude $`U_0`$ and width $`\sigma `$, which simulates a laser beam acting continuously on the centre of the trap. The fermions interacting with the laser are pushed out of the trap centre and the density distortion moves towards the ends of the trap with a velocity which corresponds to the local Fermi velocity given by Eq. (7). Figure 3 shows the evolution in time of the perturbation along the $`z`$ axis for $`r=0`$. For such a system the maximum value of the axial zero-sound velocity, which is equal to the axial Fermi velocity, can be evaluated using Eq. (7) as $`c_0^{max}|_{r=0}=v_F^{max}|_{r=0}=0.383a_{ho}/ms`$, where $`a_{ho}=(\mathrm{}/m\omega )^{1/2}`$ is the radial harmonic oscillator length. Figure 3 also shows that the positions of the density peaks as calculated numerically at various time steps are in good agreement with those estimated from Eq. (7).

We consider next a 50-50 mixture of 2000 40K atoms which are polarized in the $`7/2`$ and $`9/2`$ spin states, under the same external potentials as in the previous case. The value of the scattering length of 40K ($`a=80`$ Bohr radii $`a_0`$) is so small that $`v_F`$ is lowered by the interactions by only 0.1$`\%`$. From Fig. 2 it is evident that the speed of zero sound is practically unchanged. To see the interaction effects we study the case $`a=2\times 10^4a_0`$. For this value of the scattering length the Fermi velocity is obtained from the numerical calculation of the chemical potential and of the effective potential as $`v_F^{max}|_{r=0}=0.343a_{ho}/ms`$. The corresponding values of $`\lambda `$ and $`c_0`$ are 1.01 and $`0.347a_{ho}/ms`$ respectively (see Fig. 2). The time evolution of the perturbation in the two components of the interacting Fermi gas is shown in Fig. 4, in comparison with that of the noninteracting gas.

At still higher couplings spatial phase separation starts in the two-component Fermi gas. The case $`a=8\times 10^4a_0`$ is illustrated in Figure 5. The two components are amost completely demixed and only one of them feels the external perturbation. In this case the interactions with the second component merely shift the chemical potential entering the evaluation of the Fermi velocity of the first component, which equals the speed of propagation of zero sound.

### 3.2 Perturbation introduced by a compressed slow light pulse

In the work of Dutton et al.  ultra-compressed slow light-pulses were used to pump a small slice of a Bose-Einstein condensate into an untrapped state, leaving behind a sharp hole in the density profile. The density defect in the condensate evolves through the formation of solitons, this process being analogous to the creation of shock waves in a classical fluid.

We generate the same type of initial profile in a two-component Fermi gas for various shapes of the density vacancy (V-shape or U-shape of various sizes) and examine the effects of the density gradient and of the size of the defect. No coherent behaviour is observed in the time evolution of the density profiles. In Fig. 6 we plot the time evolution of one of the two components of the Fermi gas for various values of the scattering length ($`a=0`$, $`800a_0`$ and $`2\times 10^4a_0`$) in the case of a small U-shape vacancy. In the two former cases the particles in filling the vacancy generate some density waves, which reach a maximum amplitude after about 150 $`ms`$. This is the time lag nedeed for the perturbation to arrive at the ends of the trap and to be reflected there, as can be roughly estimated by dividing the axial half-length $`L_z`$ of the trap by $`c_0^{max}`$ giving $`L_z/c_0^{max}140`$ ms. Because of the absence of thermal dissipation in our model, the density waves continue to travel along the profile and are again reflected at the ends of the trap with a periodicity of about $`2L_z/c_0^{max}280`$ ms. The revivals of smaller amplitude which are seen to occur at about 300 ms correspond to the situation in which the perturbation is passing through the centre of the trap. On increasing the coupling strength up to $`a=2\times 10^4a_0`$ (right panel in Fig. 6), we find that the restoring force due to the other fermionic component plays a role in depressing the intensity of the density waves relative to those seen in the weak-coupling regime.

## 4 Conclusions

In summary, we have studied the propagation of collisionless sound waves in a trapped two-component Fermi gas, assuming a mean-field model for the repulsive interactions between the two fermionic species and neglecting dissipative collisions as described by a Boltzmann quantum collision integral. The inclusion of such collisions is needed to examine the transition from the zero-sound to the first-sound regime , that we expect will be characterized by a marked change in the speed of sound propagation. We also expect that zero-sound waves will be Landau-damped in the case of attractive interactions, as is the case for a homogeneous Fermi liquid .

We have found that, both at weak coupling and in the strong-coupling regime where the two Fermi clouds are demixing from their mutual repulsions, the speed of zero-sound waves is very close to the Fermi velocity as determined by the local particle density through the local chemical potential. These results have been obtained from comparing semi-analytical estimates with numerical studies of the evolution of peaks in the density profiles under driving from an external potential simulating a blue-detuned laser that keeps expelling atoms from the central region of the trap.

We have also examined numerically the time evolution of impulsively created density defects which are characterized by sharp density gradients in a two-component Fermi gas with repulsive interactions. In the collisionless regime we have found that the filling of the vacancy at the trap centre generates density waves which travel to the ends of the trap and are back-reflected there, so that they appear as periodic revivals of distortions in the density profile. These revivals become attenuated as the strength of the coupling is increased, presumably as a consequence of the increasing efficiency of momentum-redistribution processes through non-dissipative collisions at strong coupling.

We acknowledge support from INFM through PRA2001-Photonmatter. ZA also acknowledges support from TUBITAK and from the Research Fund of the University of Istanbul under Project Number BYP-110/12122002.


text:
# Scaling Properties of Long-Range Correlated Noisy Signals

## Abstract

The Hurst coefficient $`H`$ of a stochastic fractal signal is estimated using the function $`\sigma _{MA}^2=\frac{1}{N_{max}n}_{i=n}^{N_{max}}[y(i)\stackrel{~}{y}_n(i)]^2`$, where $`\stackrel{~}{y}_n(i)`$ is defined as $`1/n_{k=0}^{n1}y(ik)`$, $`n`$ is the dimension of moving average box and $`N_{max}`$ is the dimension of the stochastic series. The ability to capture scaling properties by $`\sigma _{MA}^2`$ can be understood by observing that the function $`C_n(i)=y(i)\stackrel{~}{y}_n(i)`$ generates a sequence of random clusters having power-law probability distribution of the amplitude and of the lifetime, with exponents equal to the fractal dimension $`D`$ of the stochastic series.

Power-law distributions are important manifestations of scale invariance as observed in fractal, percolating structures and self-organized dynamical systems. The scaling exponents of power law statistics are indeed related to the universality class of the system and thus are helpful to understand the fundamental processes ruling the dynamics of complex systems. The scaling parameters have been practically deployed in fields as different as biology, geophysics, solar physics, social science. When random time records should be analyzed the Hurst exponent $`H`$, which is related to the fractal dimension of the time series by $`D=2H`$, is usually evaluated. For example, healthy and sick heart beat rate can be distinguished on the basis of the value of $`H`$. Financial series with degree of persistence higher than that of the price series have been encountered. A number of frequency, time and, even, integrated domain approaches have been thus developed to gain as accurate as possible estimate of these exponents. Such procedures generally consist in calculating appropriate statistical functions from the whole signal: in the time domain, Detrended Fluctuation Analysis (DFA) and Rescaled Range Analysis (R/S) are the most popular scaling methods to extract power-law correlation exponents from the random signals $`y(i)`$ Feder ; Mantegna ; Bak ; Sethna ; Vandewalle ; KunHu ; Rangarajan ; Heneghan ; Alessio ; Lux ; Ivanov .

A generalized variance $`\sigma _{MA}^2`$ of the noisy signal $`y(i)`$ with respect to $`\stackrel{~}{y}_n(i)`$, defined by:

$$\sigma _{MA}^2=\frac{1}{N_{max}n}\underset{i=n}{\overset{N_{max}}{}}C_n(i)^2,$$
(1)
with:

$$C_n(i)=y(i)\stackrel{~}{y}_n(i),$$
(2)
has been introduced in Alessio . In Eq.(1), $`N_{max}`$ is the size of the series, $`\stackrel{~}{y}_n(i)`$ is the moving average defined by $`\stackrel{~}{y}_n(i)=\frac{1}{n}_{k=0}^{n1}y(ik)`$ and $`n`$ is the moving average window.

The function $`\sigma _{MA}`$ corresponding to each $`\stackrel{~}{y}_n(i)`$ has been calculated for random series $`y(i)`$ with different Hurst exponent $`H`$. The values of $`\sigma _{MA}`$ are plotted as a function of $`n`$ on log-log axes shown in Fig., for a series with $`N_{max}=2^{19}`$ and $`10n10^4`$. The most remarkable property of the curves plotted in Fig. is their power-law dependence on $`n`$, i.e.:

$$\sigma _{MA}n^H.$$
(3)
The structure of the algorithm based on the $`\sigma _{MA}`$ function is analogous to the DFA technique, but a higher accuracy has been observed. This may be due to the better smogthing and detrending action of the moving average with respect to the linear or to the polynomial fit. Furthermore, the moving average filter dynamically detrends the series. Every time the discrete index $`i`$ increases by a unity, the box window $`n`$ switches its position accordingly, allowing also for online applications of the present technique. In order to gain a deeper insight of the function $`\sigma _{MA}`$, we will show that $`C_n(i)`$ generates, for each $`\stackrel{~}{y}_n(i)`$, a sequence of clusters with amplitude and lifetime distributed as power-laws. Furthermore the exponent of these power-laws is equal to tte fractal dimension $`D=2H`$ of the time series. These properties allow to understand the ability of $`\sigma _{MA}`$ to estimate the Hurst exponent $`H`$ of the stochastic series.

Consider thus the sub-set of $`y(i)`$ corresponding to the region delimited by two consecutive intersections between $`y(i)`$ and $`\stackrel{~}{y}_n(i)`$. Let us refer to such regions as cluster. In Fig., one of these clusters is shown. For each cluster $`j`$, the size $`l_j`$ :

$$l_j=\underset{i=i_c(j)}{\overset{i_c(j+1)}{}}y(i)$$
(4)
and the duration $`\tau _j`$:

$$\tau _j=i_c(j+1)i_c(j).$$
(5)
can be defined. In the previous relationships, $`i_c(j)`$ and $`i_c(j+1)`$ are the values taken by the index $`i`$ in correspondence of two subsequent intersections between $`\stackrel{~}{y}_n(i)`$ and $`y(i)`$. The size $`l_j`$ and the duration $`\tau _j`$ of the clusters are shown in Fig.(2) for a series with $`N_{max}=2^{19}`$, $`H=5.3`$ and $`n=30`$.

The probability distribution $`P(z)`$ of the cluster sizes has been determined by counting the number of clusters with length $`l`$. Fig. shows the log-log plot of the distribution $`P(l)`$ for series having different $`H`$ (respectively $`H=0.2`$, $`H=0.5`$ and $`H=0.8`$). The clusters are obtained by the intersection of the series with a moving average with a window amplitude $`n=60`$. The curves are consistent with a straight line, indicating a power-law behavior:

$$P(l)l^\alpha $$
(6)
The linear behavior over two decades indicates a scaling distribution of clusters. The exponent $`\alpha `$ has been plotted against $`H`$ can the inset of Fig..

The probability distribution $`P(\tau )`$ of the cluster lifetime has been determined by counting the number of clusters with lifetime $`\tau `$. Fig. shows the log-log plot of the distribution $`P(\tau )`$ for series having different $`H`$ (respectively $`H=0.2`$, $`H=0.4`$ and $`H=0.8`$). The clusters are obtained by the intersection of the series with a moving average with a window amplitude $`n=60`$. This leads to another line indicating a distribution of lifetimes of the form:

$$P(\tau )\tau ^\beta .$$
(7)
The values of $`\alpha `$ and $`\beta `$, for $`H`$ varying from $`0.10`$ to $`0.90`$, are reported in Table (1). The distributions $`P(l)`$ and $`P(\tau )`$ have been calculated for a wide range of values of $`N_{max}`$ and $`n`$ ( $`2^{18}<N_{max}<2^{21}`$ and $`60<n<1000`$), the exponents $`\alpha `$ and $`\beta `$ have been found to be independent of $`N_{max}`$ and $`n`$. It can be concluded that $`\alpha `$ and $`\beta `$ vary as $`2H`$, i.e. as the fractal dimension $`D`$ of the noisy signal.

These results can be understood keeping in mind the box-counting method to estimate the fractal dimension of a random signal is recalled (Ch.10 of Feder ). The stochastic signal is covered with boxes of width $`b\tau `$ (in time) and of length $`ba`$ (in amplitude), with $`\tau a`$ the minimum box. The box dimension $`D_B`$ is related to the number $`N(b;a,\tau )`$ of the boxes needed to cover the record by the relationship: $`N(b;a,\tau )b^{D_B}`$, where $`D_B=2H`$ for self-affine records.

We have demonstrated that the Hurst exponent $`H`$ of long-range correlated stochastic series $`y(i)`$ can be calculated using the function $`\sigma _{MA}`$ defined by the Eq.(1). We have found indeed the remarkable result that the function $`\sigma _{MA}`$ varies as a power-law of the amplitude $`n`$ of the moving average box. This work shed more light on the results of the papers Vandewalle ; Alessio  and reinforce the general idea that a deeper physical insight may add new perspectives to the practical methods of financial analysis.

These results, appeared in very preliminary form in Alessio , have been here validated in terms of the statistical properties of the random clusters tracked by $`\stackrel{~}{y}_n(i)`$ over the random time series. Furthermore, the power-law behavior of the distributions of cluster size and lifetime, just as expected from general arguments about dynamical system, clarifies the intrinsic capability of the function $`\sigma _{MA}`$ to capture scaling exponents of the series. Furthermore, the exponents $`\alpha `$ and $`\beta `$ of Eqs.(6,7) are insensitive to the parameters characterizing $`\stackrel{~}{y}_n(i)`$ and $`y(i)`$. This robustness is essential to the assessment of the proposed scaling technique. Last but not least, the $`\sigma _{MA}`$ algorithm yields in higher accuracy, speed of execution and possibility of online application.