text:
# Unitary transformation approach for the trapped ion dynamics

## I Introduction

Trapped ions interacting with laser beams  have become an extremelly interesting system for the investigation of fundamental physics, e.g. the generation of states of the harmonic oscillator  as well as for potential applications such as precision spectroscopy  and quantum computation . The theoretical treatment of the interaction of a trapped ion with one or several laser beams constitutes a complicated problem, being the full Hamiltonian highly nonlinear. Therefore approximations are normally required, such as, for instance, the Lamb-Dicke approximation, in which the ion is considered to be confined within a region much smaller than the laser wavelength. This makes possible to obtain Hamiltonians of the Jaynes-Cummings type, in which the center-of-mass of the trapped ion plays the role of the field mode in cavity QED. Recently , it has been suggested a new approach to this problem, based on the application of a unitary transformation which linearizes the total ion-laser Hamiltonian. Moreover, in that approach  it is obtained a Jaynes-Cummings type Hamiltonian including counter-rotating terms. A rotating wave approximation (RWA) may be performed provided we are in the special “resonant” regime, or $`\mathrm{\Omega }\nu `$, where $`\mathrm{\Omega }`$ is basically related to the laser intensity and $`\nu `$ is the frequency of the ion in the magneto-optical trap. The RWA performed this way allows the exact diagonalization of the transformed Hamiltonian, but on the other hand poses an upper limit to the possible values of the Lamb-Dicke parameter. In this paper we propose a method of diagonalization of the ion-laser Hamiltonian entirely based on unitary transformations. As a first step we linearize the Hamiltonian, as it is done in . The resulting Hamiltonian is further transformed, and cast in a form that it is suitable for diagonalization. Although the diagonalization procedure is an approximate one, based on recursion relations, it is not based in an expansion on any of the parameters of the problem. Therefore, the found solution does not depend on any specific regime, such as the Lamb-Dicke regime, for instance.

This paper is organized as follows: In section 2 we obtain the linearized form of the full Hamiltonian. In section 3 we show how the resulting Hamiltonian may be further transformed in order to be diagonalized. In section 4 we summarize our conclusions.

## II Linearization of the Hamiltonian

We consider a single trapped ion interacting with two laser (classical) plane waves (frequencies $`\omega _1`$ and $`\omega _2`$), in a Raman-type configuration ($`\omega _L=\omega _1\omega _2`$). The lasers effectively drive the electric-dipole forbidden transition $`|g|e`$ (frequency $`\omega _0`$), with $`\delta =\omega _0\omega _L`$. We end up with an effective two-level system for the internal degrees of freedom of the atom and the vibrational motion, after adiabatically eliminating the third level in the Raman configuration. This situation is described, in the atomic basis, by the following effective Hamiltonian

$$\widehat{H}=\mathrm{}\left(\begin{array}{cc}\nu \widehat{n}+\frac{\delta }{2}& \mathrm{\Omega }e^{i\eta \widehat{X}}\\ \mathrm{\Omega }e^{i\eta \widehat{X}}& \nu \widehat{n}\frac{\delta }{2}\end{array}\right),$$
(1)
being $`\widehat{X}=\widehat{a}+\widehat{a}^{}`$, $`\widehat{n}=\widehat{a}^{}\widehat{a}`$, $`\eta `$ the Lamb-Dicke parameter, and $`\widehat{a}^{}`$ ($`\widehat{a}`$) the ion’s vibrational creation (annihilation) operator.

By applying the unitary transformation

$$\widehat{T}=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}\widehat{D}^{}(\beta )& \widehat{D}^{}(\beta )\\ \widehat{D}(\beta )& \widehat{D}(\beta )\end{array}\right),$$
(2)
to the Hamiltonian in equation (1), where $`\widehat{D}(\beta )=\mathrm{exp}(\beta \widehat{a}^{}\beta ^{}\widehat{a})`$ is Glauber’s displacement operator, with $`\beta =i\eta /2`$, we obtain the following transformed Hamiltonian

$$\widehat{}\widehat{T}^{}\widehat{H}\widehat{T}=\mathrm{}\left(\begin{array}{cc}\nu \widehat{n}+\nu \frac{\eta ^2}{4}+\mathrm{\Omega }& i\lambda \widehat{Y}\frac{\delta }{\eta \nu }\\ i\lambda \widehat{Y}\frac{\delta }{\eta \nu }& \nu \widehat{n}+\nu \frac{\eta ^2}{4}\mathrm{\Omega }\end{array}\right),$$
(3)
where $`\widehat{Y}=\widehat{a}\widehat{a}^{}`$, and $`\lambda =\frac{1}{2}\eta \nu `$. This result holds for any value of the Lamb-Dicke parameter $`\eta `$. We have thus transformed our original Hamiltonian in equation (1) into a Jaynes-Cummings-type Hamiltonian, and therefore its exact diagonalization is allowed provided we perform the rotating wave approximation (RWA). However this imposes limitations on the values of the Lamb-Dicke parameter. Because the “effective coupling constant” of the transformed Hamiltonian is given by $`\lambda =\frac{1}{2}\eta \nu `$, we shall have $`\frac{1}{2}\eta \nu \nu `$ if we want to neglect the counter rotating terms (RWA). Moreover, we still have to be “in resonance”, or $`\nu 2\mathrm{\Omega }`$. In this regime the trapped ion-laser system suitably prepared exhibits long-time-scale revivals (superrevivals), as discussed in . Here we are not going to perform the RWA. We propose instead a way of further transforming the Hamiltonian in equation (3) in order to allow its diagonalization in an approximate, although nonperturbative way. In what follows we are going to present a novel method of transforming the Hamiltonian, showing how it may be diagonalized.

## III Diagonalization of the Hamiltonian

We have succeeded in transforming the ion-laser Hamiltonian in Equation (2) into a more tractable form (3). Nevertheless, its exact diagonalization has nor been achieved yet. The problem is either treated exactly after performing the RWA, or by means of approximate methods such as perturbative expansions or even numerically . Here we propose a different approach to that problem, based on unitary transformations. We are going to restrict ourselves to the resonant case $`\delta =0`$. After discarding the constant term $`\frac{1}{4}\mathrm{}\nu \eta ^2`$, which just represents an overall phase factor, we obtain

$$\widehat{}=\mathrm{}\left(\begin{array}{cc}\nu \widehat{n}+\mathrm{\Omega }& i\lambda \widehat{Y}\\ i\lambda \widehat{Y}& \nu \widehat{n}\mathrm{\Omega }\end{array}\right),$$
(4)
Now we define the following unitary transformations

$$\widehat{T}_1=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),$$
(5)
and

$$\widehat{T}_2=\left(\begin{array}{cc}(1)^{\widehat{n}}& 0\\ 0& 1\end{array}\right).$$
(6)
These transformations have interesting effects when applied to the Jaynes-Cummings type Hamiltonian in (4). First we apply $`\widehat{T}_1`$, which yields

$$\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\left(\begin{array}{cc}\mathrm{\Omega }& \lambda \widehat{Y}\\ \lambda \widehat{Y}& \mathrm{\Omega }\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=\left(\begin{array}{cc}i\lambda \widehat{Y}& \mathrm{\Omega }\\ \mathrm{\Omega }& i\lambda \widehat{Y}\end{array}\right),$$
(7)
i.e., the operators are transposed to the diagonal. Then we apply $`\widehat{T}_2`$ to the matrix which resulted form the operation above, such that

$$\left(\begin{array}{cc}(1)^{\widehat{n}}& 0\\ 0& 1\end{array}\right)\left(\begin{array}{cc}i\lambda \widehat{Y}& \mathrm{\Omega }\\ \mathrm{\Omega }& i\lambda \widehat{Y}\end{array}\right)\left(\begin{array}{cc}(1)^{\widehat{n}}& 0\\ 0& 1\end{array}\right)=\left(\begin{array}{cc}i\lambda \widehat{Y}& \mathrm{\Omega }(1)^{\widehat{n}}\\ \mathrm{\Omega }(1)^{\widehat{n}}& i\lambda \widehat{Y}\end{array}\right),$$
(8)
where it has been used the fact that $`(1)^{\widehat{n}}\widehat{a}(1)^{\widehat{n}}=\widehat{a}`$ (the same for $`\widehat{a}^{}`$). We have now that $`\widehat{Y}`$ is multiplied by the identity matrix. At this stage we apply $`\widehat{T}_1`$ again, rearranging the terms in such a way that we obtain a Hamiltonian diagonal in the atomic state basis, or

$$\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\left(\begin{array}{cc}i\lambda \widehat{Y}& \mathrm{\Omega }(1)^{\widehat{n}}\\ \mathrm{\Omega }(1)^{\widehat{n}}& i\lambda \widehat{Y}\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=i\lambda \widehat{Y}+\left(\begin{array}{cc}\mathrm{\Omega }(1)^{\widehat{n}}& 0\\ 0& \mathrm{\Omega }(1)^{\widehat{n}}\end{array}\right).$$
(9)
Diagonalization of the total Hamiltonian becomes an easy task now, because the atomic part is already in a diagonal form. The total transformed Hamiltonian may be conveniently expressed in an operator form as

$$\widehat{}^{}=\widehat{T}_1^{}\widehat{T}_2^{}\widehat{T}_1^{}\widehat{T}^{}\widehat{H}\widehat{T}\widehat{T}_1\widehat{T}_2\widehat{T}_1=\mathrm{}\left(\nu \widehat{n}+i\lambda \widehat{Y}\mathrm{\Omega }\sigma _z(1)^{\widehat{n}}\right).$$
(10)
A general expression for the eigenstates of the Hamiltonian in expression (10) is

$$|\mathrm{\Psi }_l=|\phi _l^g|g+|\phi _l^e|e,$$
(11)
with $`|\phi _l^e=_nC_{n,l}^e|n`$ and $`|\phi ^g=_nC_{n,l}^g|n`$. From the eigenvalues equation $`\widehat{}^{}|\mathrm{\Psi }_l=\mathrm{\Lambda }_l|\mathrm{\Psi }_l`$ we obtain

$$\left(\nu \widehat{n}+i\lambda \widehat{Y}\mathrm{\Omega }(1)^{\widehat{n}}\right)\underset{n}{}C_{n,l}^e|n=\mathrm{\Lambda }_l\underset{n}{}C_{n,l}^e|n,$$
(12)
and

$$\left(\nu \widehat{n}+i\lambda \widehat{Y}+\mathrm{\Omega }(1)^{\widehat{n}}\right)\underset{n}{}C_{n,l}^g|n=\mathrm{\Lambda }_l\underset{n}{}C_{n,l}^g|n.$$
(13)
The expansion coefficients $`(C^{}s)`$ may be obtained by means of recursion relations. For instance, for the coefficient $`C_{l,n}^e`$. From equation (11) we may write

$$\underset{n}{}\left(nC_{n,l}^e|n+i\lambda C_{n,l}^e\sqrt{n}|n1i\lambda C_{n,l}^e\sqrt{n+1}|n+1\mathrm{\Omega }(1)^nC_{n,l}^e|n\right)=\mathrm{\Lambda }_l\underset{n}{}C_{n,l}^e|n.$$
(14)
After rearranging some of the terms, we have the following relation between the coefficients

$$C_{n+2,l}^e=\frac{\sqrt{n+1}}{\sqrt{n+2}}C_{n,l}^e\frac{i}{\lambda }\frac{\left[\mathrm{\Lambda }+\mathrm{\Omega }(1)^{n+1}\nu (n+1)\right]}{\sqrt{n+2}}C_{n+1,l}^e.$$
(15)
By multiplying expression (14) by $`0|`$, we also obtain the first term

$$C_{1,l}^e=i\frac{\mathrm{\Omega }+\mathrm{\Lambda }}{\lambda }C_{0,l}^e.$$
(16)
Similar relations may be found for the coefficients $`C_{n,l}^g`$. Of course the normalization condition $`_n|C_{n,l}^g|^2+_n|C_{n,l}^e|^2=1`$ should be satisfied. Having diagonalized the transformed Hamiltonian, the evolution of the state vector becomes a trivial task. For that we have to express a generic state $`|\psi `$ in terms of the basis states, or

$$|\psi =\underset{l}{}A_l|\mathrm{\Psi }_l.$$
(17)
The choice of a specific initial state determines the set of coefficients $`A_l`$. For instance, if we initially prepare the trapped ion in the state $`|\mathrm{\Psi }(0)=\frac{1}{\sqrt{2}}|\beta \left(|g|e\right)`$, the transformed state will read

$$|\psi (0)=\widehat{T}_1^{}\widehat{T}_2^{}\widehat{T}_1^{}\widehat{T}^{}|\mathrm{\Psi }(0)=|0|e=\underset{l}{}A_l|\mathrm{\Psi }_l.$$
(18)
The coefficients $`A_l`$ may be determined by multiplying the left hand side of Equation (18) by $`n|e|`$ and $`n|g|`$, so that we obtain a set of coupled equations for the coefficients $`A_l`$. It is now easy to calculate the time-evolution of the state vector, namely

$$|\psi (t)=\mathrm{exp}(i\widehat{}^{}t/\mathrm{})\underset{l}{}A_l|\mathrm{\Psi }_l=\underset{l}{}A_l\mathrm{exp}(i\mathrm{\Lambda }_lt/\mathrm{})|\mathrm{\Psi }_l.$$
(19)
The next step is to apply on $`|\psi (t)`$ the sequence of transformations $`\widehat{T}_1^{}\widehat{T}_2^{}\widehat{T}_1^{}\widehat{T}^{}`$ backwards in order to recover the wanted solution for the state vector which describes the trapped ion system interacting with the laser beams.

## IV Conclusions

We have presented a novel way of treating the problem of the interaction of trapped ions with laser beams entirely based on unitary transformations. The system Hamiltonian is successively modified and cast in a more tractable form. We propose a method of diagonalization of the transformed Hamiltonian by means of the construction of recursive relations for expansion coefficients in the Fock state basis. Despite of not being an exact diagonalization, there is no need of performing any approximation, such as the Lamb-Dicke limit, for instance. Interesting possibilities for the investigation of such a system in different regimes are opened up, given that there are no restrictions on the relevant parameters and no approximations have been made so far.

###### Acknowledgements.

This work was partially supported by Consejo Nacional de Ciencia y Tecnología (CONACyT), México, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).


format: text

text:
# The Detection of Multimodal Oscillations on 𝛼 UMa

## 1 Introduction

Over the past few decades, our understanding of the interior of the Sun – its thermodynamic structure, internal rotation, and dynamics – has been revolutionized by the technique of helioseismology, the study of the frequencies and amplitudes of seismic waves that penetrate deep into the solar interior (Leibacher et al. 1985; Duvall et al. 1988; Schou et al. 1998). The high quality of the modeling in the solar case is made possible by the large number ($`10^7`$) of modes visible in the Sun. Unfortunately, the lack of spatial resolution inherent in stellar observations limits the number of detectable modes in stars to only a few (those with low degree $`l`$). Nonetheless, successful detection of even a few modes has the potential to provide greatly improved values for fundamental stellar parameters such as mass, abundance, and age (Gough 1987; Brown et al. 1994)

While oscillations have been successfully detected on roAp stars, $`\delta `$ Scuti stars, and white dwarfs, these stars all show oscillation amplitudes several orders of magnitude larger than is expected from solar-like stars. More recently, however, a number of authors (Hatzes & Cochran 1994a, 1994b; Edmonds & Gilliland 1996) have reported the detection of periodic variability at levels of a few hundred $`\mathrm{m}\mathrm{s}^1`$ and/or several millimagnitudes in a number of K giants. However, each detection is of only a single mode, and multimodal oscillations have yet to be unambiguously detected in any cool star other than the Sun.

In March 1999, NASA launched the Wide-Field Infrared Explorer satellite, with the intent of carrying out an infrared sky survey to better understand galaxy evolution. Unfortunately, within days of launch the primary science instrument on WIRE failed due to loss of coolant. However, the satellite itself continues to function nearly perfectly, and in May we began a program of asteroseismology using the WIRE’s onboard 52 mm aperture star camera. Below we report WIRE’s probable detection of multimodal oscillations on a cool star, which is the first of its kind.

## 2 Instrument Description and Data Reduction Technique

The WIRE satellite star camera, a Ball Aerospace model CT-601, consists of a $`512\times 512`$ SITe CCD with $`27\mu `$m pixels (1 arc minute on the sky) and gain of $`15\mathrm{e}^{}/\mathrm{ADU}`$, fed by a 52 mm, f/1.75 refractive optical system. The read noise of the system is 30 electrons. The pixel data is digitized with a 16 bit ADC, and up to five $`8\times 8`$ pixel fields can be digitized and transmitted to the ground. For this work, we used only one field, which permitted us to read out the CCD at a rate of 10 Hz. The stellar image is somewhat defocused, but essentially all of the light falls on the central $`2\times 2`$ pixel spot. The spectral response of the system is governed entirely by the response of the CCD plus the optical system, and is approximately equivalent to the V+R bandpass.

The WIRE satellite is in a sun-synchronous orbit which, when combined with constraints imposed by the solar panels, limits pointing to two strips, each approximately $`\pm 30^{}`$ wide, located perpendicular to the Earth-Sun line. In addition, continuous observing is not possible. Early in the program, we were able to obtain only 7 or 8 minutes worth of data during every 96-minute orbit. Later, after viewing constraints were relaxed (which involved scheduling software and onboard data-table modifications), observing efficiency rose to as much as 40 minutes per target per orbit – up to two targets are possible during any orbit. Thus, with integrations every 0.1 s, we acquired continuous data segments of up to 24,000 observations.

Bias correction was performed on board the satellite, and further data reduction was accomplished using software developed at IPAC. Each $`8\times 8`$ pixel field was extracted by summing the flux in the central $`4\times 4`$ pixel region. Although scattered light in the field is limited by the one-meter sun-shield mounted on the star camera, we performed a background subtraction using the flux from a 20-pixel octagonal annulus surrounding the central region of the image. Finally, we converted our fluxes to an instrumental magnitude. After removal of thermal effects (see below), the rms noise in the final reduced time series was typically comparable to the 1.8 mmag noise expected from pure photon statistics, although non-Poisson noise is certainly present as well. The lack of a good flat field for the instrument was a concern, which we dealt with by rejecting those frames in which the mean centroid of the stellar image lay more than $`4\sigma `$ from the mean position, where $`\sigma `$ is the mean standard deviation of the image centroid. This criterion applied to approximately 3.8% of the observations, and the vast majority of these observations were at the start or end of an orbital segment. Overall, the satellite displayed excellent attitude stability during our run, with $`\sigma `$ measured to be typically 0.7 arc sec or less (Laher et al. 2000).

## 3 Observations and Data Analysis

$`\alpha `$ UMa was the primary target for WIRE from 18 May through 23 June 1999. It is a K0 III star (Taylor 1999) with an angular diameter of $`6.79\mathrm{mas}`$ (Hall 1996; Bell 1993). At the Hipparcos distance of 38 pc, this corresponds to a stellar radius of $`28\mathrm{R}_{\mathrm{}}`$, in substantial agreement with the earlier value of $`25\mathrm{R}_{\mathrm{}}`$ derived by Bell (1993). The effective temperature of the star is variously reported as $`T_e=3970\mathrm{K}`$ to $`T_e=4660\mathrm{K}`$ (Cayrel de Stroebel 1992), with the latter value being the most recent (Taylor 1999).

$`\alpha `$ UMa is a member of a binary system, with a total system mass of $`5.94\mathrm{M}_{\mathrm{}}`$ (Soderhjelm 1999). The spectral type of $`\alpha `$ UMa B is somewhat unclear, with F7 V most often cited in the literature (see, e.g., the SAO catalog). However, on the basis of IUE observations, Kondo, Morgan, & Modisette (1977) estimate the secondary to be “late A” (see also Ayres, Marstad, & Linsky 1981). In either case, however, the secondary makes only a small contribution ($`6\%`$) to the total system luminosity, and should show oscillation frequencies quite different from those of the primary.

During the observation period, WIRE made a total of $`4,036,448`$ observations of $`\alpha `$ UMa, after removal of those observations with poor pointing characteristics. The data were first phased at the period of the spacecraft, to determine if any obvious instrumental periodicities existed. This phasing showed the existence of a strong sinusoidal variation, with amplitude $`8`$ mmag. A thermistor is mounted on the star camera, and examination of data from this thermistor showed that these variations were in fact correlated with temperature variations of a few tenths of a degree. Although a thermoelectric cooler (TEC) is mounted on the CCD, the star tracker thermal design lowered the CCD temperature below the default setpoint, so the TEC never actually turned on. In order to reduce the impact of this signal on our data analysis, we prewhitened the data by fitting and subtracting from the entire time series a sinusoid constrained to have the satellite orbital period (the phase was allowed to vary). The best-fit sinusoid has an amplitude of 8.63 mmag, and its subtraction results in an rms residual of 2.05 mmag, compared to the 1.84 mmag expected from photon statistics. We also explored other means of fitting and removing the thermal signature, including high-order polynomial fits to the phased data. Such approaches, while more complex than sine fitting, did not lead to any appreciable improvement in the fit, and were therefore discarded. The use of any of these thermal fitting procedures did not affect the peaks in the amplitude spectrum described below. In fact, these peaks were visible even before the application of any fit to the thermal variation, though removal of the thermal signature, by removing the largest single peak, did help to render them more easily visible.

Data from the first portion of the run (18 May - 7 June) were all in short ($`t<8\mathrm{minutes}`$) segments, which obviously extended over a short range of orbital phase ($`0.08`$). We were concerned that including these data in our sine fit would bias the result, so we decided to exclude them from our analysis. In addition, NASA responded to our expressed concern about the TEC by lowering the set point on 18 June, and the CCD behavior subsequent to this had not yet reached an equilibrium point before the end of our run. Thus, we excluded data taken after 18 June from the analysis as well, leaving only the 8 June – 18 June window of approximately $`9.2\times 10^5\mathrm{s}`$. Shortening the time series clearly adversely affects our frequency resolution, but we felt that confidence in the data quality was the overriding issue. The data that were subjected to analysis are shown in Figure 1, after removal of the thermal variation.

Data were searched for periodicities using Discrete Fourier Transform (DFT; Foster 1996), Lomb-Scargle periodogram (Scargle 1982; Horne & Baliunas 1986), and epoch-folding techniques (Davies 1990) which are essentially equivalent to phase dispersion modulation (PDM; see Stellingwerf 1978; Schwarzenberg-Czerny 1998). The Scargle periodogram analysis was conducted as described in Scargle (1982; see also Hatzes & Cochrane 1994ab); the data were windowed using a Parzen function and the resulting periodogram was oversampled by a factor of 8 in frequency space. The DFT analysis was similarly oversampled, and implemented the CLEAN algorithm of Roberts, Lehar, & Dreher (1987) to remove alias peaks. When using the DFT, the data were windowed using a Parzen function, and 100 iterations of CLEAN were performed. Due to its relative slowness, epoch-folding analysis was conducted only within the frequency range of interest, as identified by the Scargle and DFT analysis, and was used only to aid in interpretation of the results from the other two algorithms. In the discussion that follows, we will concentrate on the Lomb-Scargle periodogram results, though, in general, the three techniques gave similar results.

Figure 2 shows the window function for the time series. The upper frequency limit for the figure has been arbitrarily set at $`5\mathrm{mHz}`$ to enhance the visibility of the amplitude spectrum at frequencies near 1 mHz, while the lowest frequencies are shown in the inset. No significant features are present in the window function above 5 mHz. The large evenly spaced peaks correspond to the satellite orbital frequency and its aliases.

The Lomb-Scargle periodogram for the time series is shown in Figure 3 on the same frequency scale as the window function. The low-frequency inset shows ten significant peaks, and the frequencies and amplitudes of these peaks are given in Table 1 (derived from Lorentzian fits), along with conservative formal error estimates derived from the half-width of the periodogram peaks. We note that periodograms of portions of the time series (halves and thirds) give results similar to that of the whole, with decreased frequency resolution, and that sine fits to these frequencies show coherent phasing in these different portions.

Unfortunately, the implementation of the on-board data collection on WIRE means that we lack simultaneous observations of a comparison star, and are thus essentially performing absolute photometry with an instrument not designed for that purpose. However, we have observed stars other than $`\alpha `$ UMa, and the sun-synchronous orbit of WIRE implies that most instrumental effects should be similar for all sources. The dashed line in the Figure 3 inset shows the periodogram from a time series of $`\alpha \mathrm{Leo}`$, a B7 V star not expected to show significant low-frequency oscillations. The $`\alpha `$ Leo data set, which was obtained from 23 May through 3 June 1999, consists of segments similar in length to those of $`\alpha `$ UMa, has similar rms noise to the $`\alpha `$ UMa time series, and was reduced in exactly the same manner. Our object here is not to perform analysis of the $`\alpha `$ Leo data (which might well benefit from a different approach than we have used for $`\alpha `$ UMa), but rather to show that the particular low frequency peaks in the periodogram of $`\alpha `$ UMa do not arise from either instrumental effects or the data reduction procedure itself. The larger peaks in the $`\alpha `$ Leo periodogram may arise from an imperfect removal of both long-term and orbital variations; none show coherent phasing across different segments of the time series. The dissimilarity of the two periodograms increases our confidence that the peaks visible in the $`\alpha `$ UMa periodogram are due to the star itself, although it is of course possible that the instrumental behavior changes significantly for different targets.

As is apparent from Figures 2 and 3, the family of peaks visible at low frequencies is repeated at higher frequencies, which leads to the difficulty of determining which set of peaks is the correct one. We can easily eliminate peaks above $`200\mu \mathrm{Hz}`$ by examining the summed power spectrum of the individual orbital segments, which shows no significant power at these frequencies. We are therefore left with the problem of selecting between the set of peaks in the $`250\mu \mathrm{Hz}`$ range and the similar set around $`200\mu \mathrm{Hz}`$. We believe that the low-frequency peaks are the physical solution and the higher-frequency set an alias for the following two reasons:

1. The low-frequency peaks are always of larger amplitude than the alias peaks. While a resonance of stochastic noise with an alias frequency can enhance an individual alias peak such that it is larger than the corresponding true peak, this is unlikely to occur simultaneously for multiple peaks. We have performed simulations which confirm this reasoning.
2. Hatzes (1999, private communication) has searched for oscillations in $`\alpha `$ UMa using ground-based spectroscopic methods. He reports finding frequencies of 1.36 and 6.0 $`\mu \mathrm{Hz}`$ (the latter less convincingly), although his observing run was too short to lend much confidence to the exact values. While not identical to the frequencies we report here, these frequencies are certainly comparable to our lower-frequency peaks rather than to the alias peaks.

Though neither of these factors is conclusive on its own, we believe that together they indicate that we are on solid ground interpreting the observed peaks in the amplitude spectrum as stellar oscillations. Of course, it remains possible that the observed variations are due to instrumental effects or to non-oscillatory stellar phenomena such as granulation.

## 4 Interpretation

Below we discuss the astrophysical implications of our results. A more detailed discussion, in the context of a complete stellar interiors model for $`\alpha `$ UMa A, can be found in Guenther et al. (1999).

### 4.1 Mode Frequencies and Spacings

The frequency of the fundamental mode is determined primarily by structure in the envelope and so its determination requires a complete stellar model. However, we can easily determine the range in which it should lie. The fundamental period $`P_0`$ is given by (Christy 1966; 1968)

$$Q_0=P_0\sqrt{\frac{\rho }{\rho _{\mathrm{}}}}$$
(1)
where $`Q_0`$ should lie between the value of 0.038 for a polytrope with $`\gamma =4/3`$ and 0.116 (for $`\gamma =5/3`$). Using the values of $`R=28R_{\mathrm{}}`$ from interferometry and the mass $`M4M_{\mathrm{}}`$ appropriate to a K0 giant yields a fundamental mode of between 2.8 and 8.6 days. The lowest frequency that we see corresponds to a period of 6.35 days, so we identify it as the fundamental mode for $`\alpha `$ UMa.

As noted above, the average mode spacing for the first 8 modes is $`2.94\mu \mathrm{Hz}`$. The last two modes have much larger spacings, which we interpret as signifying that we are not detecting all of the possible oscillation modes for the star, presumably because they are not excited to large amplitudes. The large separation $`\mathrm{\Delta }\nu _0`$ is related to the mean stellar density, as shown by Cox (1980):

$$\mathrm{\Delta }\nu _0=135\left(\frac{\rho }{\rho _{\mathrm{}}}\right)^{1/2}\mu \mathrm{Hz}$$
(2)
Using the values appropriate to a typical K0 giant yields a predicted spacing of $`1.82\mu \mathrm{Hz}`$, about half the observed value. Once again, this discrepancy can be accounted for by assuming that all modes are not excited. In particular, if only even or odd-valued radial $`n`$ modes are excited, we would expect the observed large separation to be twice the predicted value. The simplest explanation is then that only the $`l=0`$ modes are excited, and the frequencies we observe correspond to radial oscillations of the star.

### 4.2 Mode Amplitudes

The Kjeldsen & Bedding (1995) scaling law

$$\delta L/L=\frac{L/L_{\mathrm{}}}{(\lambda /550\mathrm{nm})(T_e/5777\mathrm{K})^2(M/M_{\mathrm{}})}\times 5.1\mu \mathrm{mag}$$
(3)
predicts oscillation amplitudes of $`500\mu \mathrm{mag}`$ for $`\alpha `$ UMa, which is essentially in agreement with our results. It should be noted that the WIRE data are obtained in white light and, consequently, phase differences in the oscillation amplitudes as a function of wavelength would tend to combine to reduce the observed amplitude. Consequently, it is not surprising that the observed amplitudes are somewhat smaller than those predicted by theory. In addition, of course, extending a relationship derived for lower main sequence stars to giants is a risky enterprise! Nonetheless, the near-agreement between theory and observation may imply that the excitation mechanism for oscillations in $`\alpha `$ UMa is fundamentally similar to the solar mechanism (presumably convection; see, e.g. Bogdan et al 1993), unlike oscillations observed in other K giants, which show amplitudes an order of magnitude greater than those we have detected.

We gratefully acknowledge the support of Dr. Harley Thronson and Dr. Phillipe Crane at NASA Headquarters for making this unusual use of WIRE possible. J.C., T.C., R.L., T.G., and D.S. would like to thank Drs. Carol Lonsdale and Perry Hacking for the opportunity to work with them on the WIRE project, and their support of the WIRE asteroseismology effort. The hard work of many people, including the WIRE operations and spacecraft teams at GSFC, and the timeline generation team team at IPAC, was essential to making this project a reality. We would also like to acknowledge the contributions of the anonymous referee, whose criticisms helped to greatly improve the presentation of our results.


text:
# Velocity Autocorrelation and Harmonic Motion in Supercooled Nondiffusing Monatomic Liquids

## 1 Introduction

Recent work by Wallace and Clements  has uncovered several important properties of the many-body potential underlying the motion of liquid sodium systems. Specifically, it has been shown that (a) the potential surface consists of a large number of intersecting nearly harmonic valleys, (b) these valleys can be classified as symmetric (crystalline, microcrystalline, or retaining some nearest-neighbor remnants of crystal symmetry) or random, with the random valleys vastly outnumbering the symmetric ones, (c) the frequency spectra of different random valleys are nearly identical (while those of the symmetric valleys vary widely), and (d) below 35 K the system remains in a single valley throughout the longest molecular dynamics (MD) runs that were performed. Results (a) through (c) verify predictions made by Wallace in his theory of liquid dynamics , which has been successfully applied to account for the high-temperature specific heats of monatomic liquids  and a study of the velocity autocorrelation function . These four results together suggest that below 35 K the motion of the atoms in liquid sodium is purely harmonic to a high degree of approximation, again as predicted by Wallace in , and we would like to test this hypothesis further. One check is to compare the mean square displacement from MD with the prediction from purely harmonic motion, which is done in Fig. 12 of , where the two are found to agree closely. However, it would be more convincing if the theory could be shown to reproduce an entire scalar function calculated from MD (instead of just a single number), such as the normalized velocity autocorrelation function $`\widehat{Z}(t)`$. That is the aim of this paper. We will show that purely harmonic motion of the atoms in a potential valley produces a $`\widehat{Z}(t)`$ which matches that of MD calculations to within the calculations’ accuracy; thus we will conclude that the motion of atoms in a nondiffusing supercooled liquid state is very nearly entirely harmonic. For completeness, in Sec. 2 we briefly review the calculation of $`\widehat{Z}(t)`$ assuming harmonic motion, and in Sec. 3 we compare this result with MD. Finally, in Sec. 4 we make contact with work by others in this field, as well as comparing these results to Wallace’s earlier effort  mentioned above.

## 2 Harmonic Theory

If an $`N`$-body system is moving in a potential valley, the potential can be expanded about the valley minimum with the resulting Hamiltonian

$$H=\underset{Ki}{}^{}\frac{p_{Ki}^2}{2M}+\underset{Ki,Lj}{}^{}\mathrm{\Phi }_{Ki,Lj}u_{Ki}u_{Lj}+\mathrm{\Phi }_A$$
(1)
where $`u_{Ki}`$ is the $`i`$th component of the $`K`$th particle’s displacement from equilibrium, $`p_{Ki}`$ is the corresponding momentum, and the anharmonic term $`\mathrm{\Phi }_A`$ contains all of the higher order parts of the expansion. The primed sum indicates that the sum is performed under the constraint that the center of mass of the system is stationary. (As a result, the system has only $`3N3`$ independent degrees of freedom.) The matrix $`\mathrm{\Phi }_{Ki,Lj}`$ is called the dynamical matrix of the system. If the valley is approximately harmonic, we can neglect $`\mathrm{\Phi }_A`$. If coordinates $`q_\lambda `$ are defined by the relation

$$u_{Ki}=\underset{\lambda }{}w_{Ki,\lambda }q_\lambda $$
(2)
where the $`w_{Ki,\lambda }`$ form a $`3N\times 3N`$ orthogonal matrix, satisfying

$$\underset{Ki}{}w_{Ki,\lambda }w_{Ki,\lambda ^{}}=\delta _{\lambda \lambda ^{}},$$
(3)
then the Hamiltonian in these new coordinates is

$$H=\underset{\lambda }{}\frac{p_\lambda ^2}{2M}+\underset{Ki,Lj}{}^{}\underset{\lambda \lambda ^{}}{}w_{Ki,\lambda }\mathrm{\Phi }_{Ki,Lj}w_{Lj,\lambda ^{}}q_\lambda q_\lambda ^{}$$
(4)
where the $`p_\lambda `$ are the momenta conjugate to the $`q_\lambda `$. Now one can always choose the $`w_{Ki,\lambda }`$ to diagonalize $`\mathrm{\Phi }_{Ki,Lj}`$, so that

$$\underset{Ki,Lj}{}w_{Ki,\lambda }\mathrm{\Phi }_{Ki,Lj}w_{Lj,\lambda ^{}}=M\omega _\lambda ^2\delta _{\lambda \lambda ^{}}.$$
(5)
(This equation defines the frequencies $`\omega _\lambda `$ in terms of the eigenvalues of $`\mathrm{\Phi }_{Ki,Lj}`$.) With this choice, the Hamiltonian becomes

$$H=\underset{\lambda }{}\left(\frac{p_\lambda ^2}{2M}+\frac{1}{2}M\omega _\lambda ^2q_\lambda ^2\right).$$
(6)
Three of the $`\omega _\lambda `$ are zero; these modes correspond to uniform motion of the center of mass. Since we have restricted the center of mass position and velocity to zero, these modes are not excited. The classical equations of motion for the remaining modes are solved by

$$q_\lambda (t)=a_\lambda \mathrm{sin}(\omega _\lambda t+\alpha _\lambda ),$$
(7)
or, returning to the original coordinates,

$$u_{Ki}(t)=\underset{\lambda }{}w_{Ki,\lambda }a_\lambda \mathrm{sin}(\omega _\lambda t+\alpha _\lambda ),$$
(8)
with the understanding that the sum on $`\lambda `$ ranges from $`1`$ to $`3N3`$. The velocities of the particles are

$$v_{Ki}(t)=\underset{\lambda }{}w_{Ki,\lambda }\omega _\lambda a_\lambda \mathrm{cos}(\omega _\lambda t+\alpha _\lambda ).$$
(9)
We compute the $`𝒗(t)𝒗(0)`$ in $`Z(t)`$ by calculating $`𝒗_K(t)𝒗_K(0)`$, summing over $`K`$ and dividing by $`N1`$ (remember that only $`3N3`$ coordinates are independent), and averaging over the amplitudes $`a_\lambda `$ and phases $`\alpha _\lambda `$. Thus

$`Z(t)`$ $`=`$ $`{\displaystyle \frac{1}{3}}𝒗(t)𝒗(0)`$ (10)
$`=`$ $`{\displaystyle \frac{1}{3N3}}{\displaystyle \underset{Ki}{}}{\displaystyle \underset{\lambda \lambda ^{}}{}}w_{Ki,\lambda }w_{Ki,\lambda ^{}}\omega _\lambda \omega _\lambda ^{}a_\lambda a_\lambda ^{}\mathrm{cos}(\omega _\lambda t+\alpha _\lambda )\mathrm{cos}(\alpha _\lambda ^{})`$
$`=`$ $`{\displaystyle \frac{1}{3N3}}{\displaystyle \underset{\lambda \lambda ^{}}{}}\delta _{\lambda \lambda ^{}}\omega _\lambda \omega _\lambda ^{}a_\lambda a_\lambda ^{}\mathrm{cos}(\omega _\lambda t+\alpha _\lambda )\mathrm{cos}(\alpha _\lambda ^{})`$
$`=`$ $`{\displaystyle \frac{1}{3N3}}{\displaystyle \underset{\lambda }{}}\omega _\lambda ^2a_\lambda ^2\mathrm{cos}(\omega _\lambda t+\alpha _\lambda )\mathrm{cos}(\alpha _\lambda )`$
$`=`$ $`{\displaystyle \frac{1}{6N6}}{\displaystyle \underset{\lambda }{}}\omega _\lambda ^2a_\lambda ^2\mathrm{cos}(\omega _\lambda t).`$
By the equipartition theorem,

$$\frac{1}{2}M\omega _\lambda ^2q_\lambda ^2=\frac{1}{2}kT$$
(11)
for any nonzero $`\omega _\lambda `$, from which it follows that

$$a_\lambda ^2=\frac{2kT}{M\omega _\lambda ^2},$$
(12)
so

$$Z(t)=\frac{1}{3N3}\frac{kT}{M}\underset{\lambda }{}\mathrm{cos}(\omega _\lambda t).$$
(13)
Notice that $`Z(0)=kT/M`$, so $`\widehat{Z}(t)`$ defined by

$$Z(t)=Z(0)\widehat{Z}(t)$$
(14)
is given in this theory by

$$\widehat{Z}(t)=\frac{1}{3N3}\underset{\lambda }{}\mathrm{cos}(\omega _\lambda t).$$
(15)
This is the result we wish to compare with MD.

To do so, we need the frequencies $`\omega _\lambda `$, which are related to the eigenvalues of the dynamical matrix $`\mathrm{\Phi }_{Ki,Lj}`$ as indicated in Eq. (5). These were evaluated for five separate random valleys in  by quenching all the way down to a valley minimum and diagonalizing $`\mathrm{\Phi }_{Ki,Lj}`$ there; as pointed out in Sec. 1, these eigenvalues were found to be independent of the specific random valley chosen. All five sets of eigenvalues are shown in Fig. 7 of , and we picked one set at random to use in performing the sum in Eq. (15); the other sets produce identical graphs of $`\widehat{Z}(t)`$.

We can also use the set of eigenvalues to reconstruct the density of frequencies $`g(\omega )`$; the results are shown in Fig. 1. Note that we do not actually integrate over this $`g(\omega )`$ to

evaluate $`\widehat{Z}(t)`$ below; we directly sum over the given set of frequencies as indicated in Eq. (15). The Figure is provided only to convey a sense of the shape of the frequency distribution. Also note that this $`g(\omega )`$ is determined from fully mechanical considerations; as a result, it is not temperature-dependent as are the frequency spectra used in Instantaneous Normal Mode (INM) studies . We will expand on this point in the Conclusion.

## 3 Comparison with MD

The MD setup used to calculate $`\widehat{Z}(t)`$ to compare with Eq. (15) is essentially that described in : $`N`$ particles interact through a potential that is known to reproduce accurately a wide variety of experimental properties of metallic sodium (see discussion in  for details). The two significant changes are that we used $`N=500`$ for all runs and that the MD timestep was reduced to $`\delta t=0.2t^{}`$, where $`t^{}=7.00\times 10^{15}`$ s is the natural timescale defined in . (The system’s mean vibrational period $`\tau =2\pi /\omega _{\mathrm{rms}}`$, where the rms frequency $`\omega _{\mathrm{rms}}`$ is calculated in , is approximately $`300\delta t`$.) We cooled the sodium sample to 22.3 K and 6.69 K, and then we ran each at equilibrium to collect velocities $`𝒗_K(t)`$ to be used to calculate $`Z(t)`$ by the formula

$$Z(t)=\frac{1}{3N}\underset{K}{}\frac{1}{n+1}\underset{t^{}=0}{\overset{n}{}}𝒗_K(t+t^{})𝒗_K(t^{}).$$
(16)
We then divided by $`Z(0)`$ to obtain $`\widehat{Z}(t)`$. The number $`n`$ was chosen as large as possible without running beyond the data calculated in the MD run. We know that at these temperatures the sodium is nondiffusing for two reasons: Both temperatures are below the 35 K threshold , and $`\widehat{Z}(t)`$ from either MD run (shown below in Fig. 2) integrates to zero, yielding zero diffusion coefficient.

The formula above may fail to produce reliable values of $`\widehat{Z}(t)`$ for three reasons. First, the number of data points in the time average may be too small; if the MD simulation is run out to time $`t_{\mathrm{max}}`$, then for a given value of $`t`$ in Eq. (16), the maximum possible value of the upper limit $`n`$ is $`t_{\mathrm{max}}t`$. Thus we require $`t_{\mathrm{max}}>>t`$; we have chosen $`t_{\mathrm{max}}=50,000`$ timesteps and we have calculated $`\widehat{Z}(t)`$ only to $`t=1000`$. To ensure that this value of $`t_{\mathrm{max}}`$ is large enough, we also performed MD runs out to 200,000 timesteps and calculated $`\widehat{Z}(t)`$ from them; the differences from the 50,000 timestep result were of order $`10^3`$. Hence we are confident that 50,000 timesteps is enough if we calculate $`\widehat{Z}(t)`$ to only 1000 timesteps. Second, it is possible that reducing the timestep (thus increasing the accuracy of the simulation) might improve the accuracy of $`\widehat{Z}(t)`$. To test this, we performed another MD run with $`\delta t`$ reduced to $`0.05t^{}`$, keeping the “real” time of the run the same; this also produced differences in $`\widehat{Z}(t)`$ of order $`10^3`$. Thus we are sure that our timestep is small enough. Finally, there is the possibility of finite size effects. Since the MD system has periodic boundary conditions, an acoustic wave sent out from the system at $`t=0`$ could propagate across the simulation region and return to its point of origin in a finite time, producing spurious correlations that would show up in $`\widehat{Z}(t)`$ but would not be present in a large-$`N`$ system. To see if this effect is relevant, we estimated the time it would take for an acoustic wave to cross the region, using the numbers from . The speed of sound in sodium at its melting temperature is $`2.5\times 10^5`$ cm/s, and the volume of the region occupied by one atom is $`278a_0^3`$, so from the fact that there are 500 atoms one finds that the time required for an acoustic wave to cross the region is $`783\delta t`$, or about 800 timesteps. (The speed of sound in sodium at our lower temperatures varies from that at the melting point by roughly 5%, so this result is valid to the same accuracy.) In the Figure below, the MD result for $`\widehat{Z}(t)`$ begins to show small oscillatory revivals at about this time; we conclude that this is a finite size effect, but it does not affect the data before that time.

In Fig. 2, Eq. (15) is plotted on top of the MD data for $`\widehat{Z}(t)`$ of sodium at the two temperatures. Although both temperatures compare exceptionally well to the harmonic theory, the match is visibly poorer for the lower temperature. However, repeated MD runs at the lower temperature revealed that overall variations in $`\widehat{Z}(t)`$ amount to $`10^2`$ on average, which is the same order as the differences between theory and MD in this Figure. By 500 timesteps the theory is slightly out of phase with the MD data, and this small difference persists out to more timesteps.

## 4 Conclusions

These results show that the motion of a liquid in a single potential valley is harmonic to an extremely high approximation; the harmonic prediction for the function $`\widehat{Z}(t)`$ matches the calculation from MD very closely. Any contributions due to anharmonicity (which are certainly present) are at most of the same order as the accuracy of the MD calculations.

Some form of harmonic approximation, such as the one used here, has been taken up by many workers attempting to understand the dynamics of liquids, and it is helpful to compare their models with our approach. One of the most popular is the theory of Instantaneous Normal Modes (INM), introduced by Rahman, Mandell, and McTague  and LaViolette and Stillinger  and developed extensively by Stratt (for example, ). Stratt expands the many-body potential in the neighborhood of an arbitrary point to second order in displacements from that point, and he expresses the potential as a quadratic sum of normal modes, in which the frequencies may be either real or imaginary. He then replaces the frequencies by their thermal averages over the potential surface, resulting in a temperature-dependent density of frequencies. From this point he calculates the system’s motion and considers various time correlation functions, including $`Z(t)`$. He observes that his results are accurate to order $`t^4`$ for short times, but his predictions also diverge from MD results very rapidly, in a time shorter than half of one vibrational period. The agreement with MD at long times can be improved by omitting the imaginary frequencies from the calculation of $`Z(t)`$, but of course this makes the short time behavior inexact. (The work of Vallauri and Bermejo  follows Stratt’s procedure.) Efforts to improve the long time behavior of the correlation functions calculated using INM have been made by Madan, Keyes, and Seeley , who have attempted to extract from the imaginary part of the INM spectrum a damping factor for $`Z(t)`$ of the general type suggested by Zwanzig . Also taking their cue from Zwanzig, Cao and Voth  have followed a slightly different path, replacing the actual potential by a set of temperature-dependent effective normal modes which, as they emphasize, bears little resemblance to the mechanical normal modes of a single many-particle valley. In fact, they state quite explicitly that a theory based on purely mechanical normal modes will have little success in accounting for equilibrium or dynamical properties of liquids.

An obvious difference between our theory and INM is the nature of our approximation. In INM, one approximates the potential quadratically at an arbitrary point, with the result that the motion so predicted is accurate only for very small times; in our theory, we expand the potential only at very special points where we know the predicted motion will be valid for very long times. Both theories then face the problem of extending their validity beyond the initial approximation, of course, and we will briefly mention our extension in the final paragraph below, but there is one particular reason why we strongly prefer the approach taken here: The other models all replace the true potential by a temperature-dependent potential determined by one or another thermal averaging process. A temperature-dependent potential does not provide a true Hamiltonian, and therefore it cannot be used to calculate the quantum or classical motion, i.e., it cannot be used in the Schrödinger equation or Newton’s law. (On the dynamical level, temperature is not even a well-defined concept.) Further, the Hamiltonian resulting from a temperature-dependent potential cannot be used to do statistical mechanics, except through uncontrolled self-consistent procedures. We prefer to build our theory in terms of the actual potential, hence in terms of its true Hamiltonian, and to find at least approximate solutions for the Hamiltonian motion, so we can apply the standard procedures of equilibrium and nonequilibrium statistical mechanics.

Further, we would argue that Cao and Voth’s skepticism regarding purely dynamical approaches is unfounded, given the results here. It is difficult to compare our $`\widehat{Z}(t)`$ results to those of others, because their MD-simulated states are not always characterized as diffusing or nondiffusing. We are fairly confident that Vallauri and Bermejo’s Fig. 2b is a comparable state (glassy Cs at 20 K), and we believe our fit to MD is slightly better. Madan, Keyes, and Seeley’s Fig. 3b is an ambiguous case (it is likely that a glass transition has occurred), but there also we are confident that our match with MD is better. Hence we would claim that this method shows as much promise as the others currently available, and with the physical potential as opposed to a thermal average potential.

It is also instructive to compare the results of this paper with a model for $`\widehat{Z}(t)`$ previously proposed by Wallace  in which a single particle oscillates in a three-dimensional harmonic valley, and at each turning point it may with probability $`\mu `$ “transit” to an adjacent valley. To apply that model to a nondiffusing case, we set $`\mu =0`$ (indicating no transits), yielding $`\widehat{Z}(t)=\mathrm{cos}(\omega t)`$. Clearly this would not fit the MD data for any $`\omega `$, and it is easy to see why: Wallace included only one frequency in his earlier model, whereas our Eq. (15) contains contributions from many frequencies, all of which are necessary to raise the first minimum in $`\widehat{Z}(t)`$ above $`1`$ and then damp $`\widehat{Z}(t)`$ out by dephasing. This suggests an alternate path to understanding diffusing states: Begin with a mean atom trajectory model that by construction reproduces the correct result for $`\widehat{Z}(t)`$ in the nondiffusing regime (Eq. (15)), and then incorporate Wallace’s notion of transits into this model. Our work in this direction, with comparison to MD data for higher-temperature diffusing states of liquid sodium, will be described in a subsequent paper .


text:
# Spectral Analysis of the Ly𝛼 Forest Using Wavelets

## 1 Introduction

Measurements of QSO spectra show that the Intergalactic Medium (IGM) is composed of highly inhomogeneous structures. Ever since their identification by Lynds (1971) and the pioneering survey of Sargent et al. (1980), these inhomogeneities have been described as discrete absorption systems, the Ly$`\alpha `$forest. With the view that the systems arise from individual intervening gas clouds, the Ly$`\alpha `$forest has been characterized using traditional absorption line statistics, most notably the line equivalent widths and, as the spectra improved in resolution and signal–to–noise ratio, the Doppler widths and $`\mathrm{I}`$column densities through Voigt profile line fitting to the features.

In the past few years, numerical simulations have successfully modelled many of the measured properties of the forest, showing that the absorption systems may arise as a consequence of cosmological structure formation (Cen et al. 1994; Zhang, Anninos & Norman 1995; Hernquist et al. 1996; Bond & Wadsley 1997; Zhang et al. 1997; Theuns, Leonard & Efstathiou 1998). The simulations have shown, contrary to the picture in which the systems are isolated intergalactic gas clouds, that most of the systems originate in an interconnected web of sheets and filaments of gas and dark matter (Cen et al. 1994; Bond & Wadsley 1997; Zhang et al. 1998). Alternative statistical methods were subsequently introduced for describing the forest using the more direct measurements of the induced light fluctuations. These include the 1-point distribution of the fluctuations (Miralda-Escudé et al. 1996; Zhang et al. 1997), and a quantity related to the 2-point distribution based on a weighted difference of the light fluctuations in neighbouring wavelength pixels (Miralda-Escudé et al. ). A direct estimate of the 2-point transmission correlation function was made by Zuo & Bond (1994).

While the newer methods for analysing the Ly$`\alpha `$forest avoid the identification of absorption lines and the fitting of Voigt profiles, they are not necessarily fundamentally different in their description of the spectra. For instance, Zhang et al. (1998) find that the distribution of optical depth per pixel in their simulation may be recovered by modelling the spectra entirely by discrete absorption lines with Voigt profiles. Rather the more direct methods circumvent a difficulty that has long plagued attempts to characterize the absorbers in terms of Voigt profiles: the sensitivity of the resulting line statistics to noise and to the fitting procedure. Absorption line fitting of necessity requires arbitrary decisions to be made regarding the setting of the continuum level, the deblending of features, and a decision on the acceptability of a fit. Different observational groups report different distributions for the line parameters. Most discrepant has been the inferred distribution of line widths. Even with the highest quality data gathered to date using the Keck HIRES, agreement is still lacking, with Hu et al. (1995) finding a narrower Doppler parameter distribution with a significantly higher mean than found by Kirkman & Tytler (1997). The differences are important, as cosmological simulations predict comparable differences for a range of plausible cosmological models (Machacek et al. 2000; Meiksin et al. 2000).

The purpose in this paper is to develop a method that provides an alternative objective description of the statistics of the Ly$`\alpha `$forest. Ultimately the goal is to employ the same method for analysing both observational data and data derived from numerical simulations in order to compare the two on a fair basis. Because of the large number of synthetic spectra generated from a simulation necessary to provide a correct average description of the forest, two principal requirements of the procedure are that it be fast and easily automated. Although automated or semi–automated Voigt profile fitting procedures exist (AutoVP, Davé et al. 1997; VPFIT, developed by Carswell and collaborators), these procedures still require arbitrary decisions to be made to obtain acceptable fits. The complexity of the codes makes it difficult to assess the statistical significance of differences between the measured distributions of the absorption line parameters and those predicted. The codes also are computationally expensive, making very costly their application to the large number of simulated spectra required to obtain a statistically valid average of the line parameters. For these reasons, a faster less complex method would be desirable. The Voigt profile fitting codes yield important parameters, like the linewidths, which contain physical information (eg, gas temperature and turbulent velocities), that the direct-analysis methods do not. It would thus be desirable for an alternative method to retain some of this information. The method presented here utilizes wavelets to characterize the absorption statistics of the Ly$`\alpha `$forest. It is not intended to be a replacement for Voigt profile fitting, but a fast alternative that allows a ready comparison between the predictions of numerical models and measured spectra and a clear statistical analysis of the results.

The outline of the paper is as follows: in §2 it is shown how the statistics of the Ly$`\alpha `$forest may be characterized using wavelets. In §3 the method is applied to the measured spectrum of a high redshift QSO. The results are summarized in §4.

## 2 Analysing the Ly$`\alpha `$Forest with Wavelets

### 2.1 Terminology

Although wavelets have been used in signal processing, image analysis, and the study of fluid dynamics for a decade, they are only beginning to enter the vernacular of astronomers. Accessible introductions are provided in Press et al. (1992), and in Slezak, Bijaoui & Mars (1990) and Pando & Fang (1996), who apply wavelets to study the clustering of galaxies and Ly$`\alpha `$absorbers, respectively. More complete accounts of wavelet methodology are Chui (1992), Daubechies (1992), and Meyer (1993). The description here is confined to those elements necessary to introduce the notation and terminology that will be used below.

Wavelets are defined variously in the literature. The definition of most use here, somewhat restrictive but appropriate to a multiresolution analysis using the Discrete Wavelet Transform (DWT), is (Meyer):

> A wavelet is a square–integrable function $`\psi (x)`$ defined in real space such that $`\psi _{jk}2^{j/2}\psi (2^jxk)`$, where $`j`$ and $`k`$ are integers, is an orthonormal basis for the set of square–integrable functions.

The wavelet $`\psi (x)`$ satisfies $`_{\mathrm{}}^{\mathrm{}}𝑑x\psi (x)=0`$, and is generally chosen to be concentrated near $`x=k2^j`$. Its defining properties permit it to perform two operations governed by the values of $`j`$ and $`k`$. Smaller values of $`j`$ correspond to coarser variations in $`f(x)`$, while differing values of $`k`$ correspond to shifting the centre of the transform.

The wavelet coefficients of a function $`f(x)`$ are defined by

$$w_{jk}𝑑xf(x)\psi _{jk}(x).$$
(1)
The set of coefficients $`\{w_{jk}\}`$ comprises the wavelet transform of the function $`f(x)`$. The function may then be recovered through the inverse transform

$$f(x)=\underset{j,k}{}w_{jk}\psi _{jk}(x),$$
(2)
since the set of functions $`\psi _{jk}`$ forms a complete orthonormal basis. The wavelet coefficients at a level $`j`$ express the changes between the smoothed representations of $`f(x)`$ at the resolution scales $`j+1`$ and $`j`$.

Several functions may serve as wavelets. A set that has proven particularly useful was developed by Daubechies (Daubechies 1992). These functions are constructed to have vanishing moments up to some value $`p`$, and the functions themselves vanish outside the range $`0<x<2p+1`$. The wavelet coefficients decrease rapidly with $`p`$ for smooth functions. Accordingly, the higher order Daubechies wavelets are the most suitable for analyzing smooth data. The DWT is computed using the pyramidal algorithm as implemented in Numerical Recipes (Press et al. ). The Daubechies wavelet of order 20 is chosen throughout.

### 2.2 Monte Carlo simulations

The properties of the wavelet transform of the Ly$`\alpha `$forest are examined by performing Monte Carlo realizations of spectra. The spectra are constructed from discrete lines with Voigt profiles using the $`\mathrm{I}`$column density and Doppler parameter distributions found by Kirkman & Tytler. Specifically, the $`\mathrm{I}`$column densities $`N_{\mathrm{HI}}`$are drawn from a power law distribution of slope 1.5 between $`12.5<\mathrm{log}_{10}N_{\mathrm{HI}}<16`$ and the Doppler parameters $`b`$ from a gaussian with mean 23 $`\mathrm{km}\mathrm{s}^1`$and standard deviation 14 $`\mathrm{km}\mathrm{s}^1`$. A cut–off in $`b`$ is imposed according to $`b>14+4(\mathrm{log}_{10}N_{\mathrm{HI}}12.5)`$ $`\mathrm{km}\mathrm{s}^1`$. The resulting average Doppler parameter is 31 $`\mathrm{km}\mathrm{s}^1`$. The number density of lines per unit redshift matches that of Kirkman & Tytler at $`z=3`$. The resolution is set at $`\lambda /d\lambda =5\times 10^4`$, and gaussian noise is added according to a specified continuum signal–to–noise ratio per pixel. This is the fiducial model used in all the simulations unless stated otherwise. Segments 128 pixels wide were found adequate for extracting the statistical properties of the wavelet coefficients.

A representative spectrum and its discrete wavelet transform are shown in Figure 1. A block at resolution $`j`$ is $`128/2^j`$ pixels wide and $`2^j`$ pixels long for $`j=1`$ to 6. The resolution becomes finer as $`j`$ increases from 1 to 6 (downwards). The uppermost level ($`j=0`$) corresponds to smoothed averages of the spectrum. The wavelet coefficients tend to increase in magnitude with decreasing resolution (decreasing $`j`$). The low values indicate that only small changes occur in the spectrum when smoothed at one resolution level to the next higher. The small values are desirable, as they signify the dominant absorption features in the spectra are adequately resolved.

Because the wavelet functions form a complete set of basis functions, the full set of wavelet coefficients completely describes the spectrum: the spectrum may be reconstructed identically from the inverse transform. For noisy spectra, however, it will generally be unnecessary to retain the full set of coefficients. Indeed, this is the motivation for multi–resolution data compression. By employing a judicious set of basis functions, a signal may be compressed into only a small fraction of its original size. The method of chosing the optimal basis set such that the compressed signal matches the original as closely as possible in a least squares sense with the least number of retained basis elements is known as Proper Orthogonal Decomposition or the Karhunen–Loéve procedure (see Berkooz, Holmes & Lumley 1993 for a review). The basis set, however, will in general differ from signal to signal if its components are highly variable, as in the case of the Ly$`\alpha `$forest. Although not optimal in the least squares sense, the wavelet basis nonetheless achieves a large amount of data compression and has the advantage of generality. Next is described how wavelets may be applied to assessing the amount of useful information in a spectrum.

Two measures of the information content of a noisy spectrum are considered, one based on $`\chi ^2`$ and the second on entropy. If $`s(x_i)`$ is the original spectrum defined at $`N`$ points $`x_i`$ (eg, wavelength or velocity), and $`s_n(x_i)`$ is the spectrum reconstructed from the $`n`$ largest (in magnitude) wavelet coefficients, then

$$\chi ^2=\underset{i=1}{\overset{N}{}}\left[\frac{s(x_i)s_n(x_i)}{\sigma _i}\right]^2$$
(3)
where $`\sigma _i`$ is the measurement error associated with pixel $`i`$. For gaussian distributed measurements, the expectation value of $`\chi ^2`$ is the number of degrees–of–freedom. If $`n`$ wavelet coefficients are retained, the number of degrees–of–freedom is $`Nn`$. (Hence, for example, $`\chi ^2=0`$ is expected for $`n=N`$.) The reduced $`\chi _{\mathrm{red}}^2=\chi ^2/(Nn)`$ then defines the optimal value of $`n`$ for truncating the wavelet coefficients.

The information content may also be expressed in terms of the wavelet coefficients directly as an “entropy”111Meyer (1993) defines the entropy to be the exponential of $`S`$.

$$S=\underset{jk}{}\alpha _{jk}^2\mathrm{log}\alpha _{jk}^2,$$
(4)
where the $`\alpha _{jk}`$ are the normalized coefficients

$$\alpha _{jk}=\frac{w_{jk}}{\left(_{jk}w_{jk}^2\right)^{1/2}}.$$
(5)
This quantity behaves like a physical entropy in the sense that it is maximum when the signal is completely random so that the full set of coefficients $`\{w_{jk}\}`$ is required to describe it, while it vanishes when the signal may be entirely described by a single coefficient.

The reduced $`\chi ^2`$ for an ensemble of Monte Carlo realizations is shown in Fig. 2 as a function of the fraction $`(Nn)/N`$ of the wavelet coefficients discarded. As the signal–to–noise ratio per pixel increases, the value of $`\chi _{\mathrm{red}}^2`$ for a given $`n`$ increases. In all cases, however, there is some $`n<N`$ for which $`\chi _{\mathrm{red}}^2=1`$. This suggests that an acceptable fit to a noisy spectrum may be provided by only a fraction $`n/N`$ of the full set of coefficients, with the fraction required increasing as the noise level decreases.

The entropy $`S`$ is shown in Figure 3. The entropy stays nearly constant out to $`\chi _{\mathrm{red}}^2=1`$, indicating that little information has been lost by discarding the small coefficients. As $`\chi _{\mathrm{red}}^2`$ increases, eventually the entropy decreases as information is lost. Due to the greater information content of the less noisy spectra, as the noise level is decreased, the entropy remains constant to increasingly higher values of $`\chi _{\mathrm{red}}^2`$ before declining.

### 2.3 Statistics of the wavelet coefficients

It was shown above how wavelets may be used to characterize the noise properties of a spectrum. The wavelet coefficients, however, may also be used to characterize the statistics of the Ly$`\alpha `$forest itself.

The distributions of the coefficients (in absolute value) for the several resolution levels are shown in Fig. 4 for a set of simulated spectra with $`S/N=50`$, typical of the Keck HIRES spectra. The number of coefficients at a level $`j`$ is $`2^j`$, with $`j=1`$ corresponding to the coarsest resolution, and $`j=6`$ to the finest for the $`2^7=128`$ pixels used in a spectrum. The finest resolution ($`j=6`$) curve is the steepest. As the resolution becomes increasingly coarse, the amplitude of the coefficients increases, as was found in Fig. 1. This indicates that most of the information in the spectrum is carried by the coarser levels (as well as by the two course scale averages, not shown). The finest level has resolved the spectral structures, with little difference between the smoothed representations of the spectrum at resolution levels $`j=5`$ and 6. Applying a cut–off in the coefficients corresponding to $`\chi _{\mathrm{red}}^2=1`$ yields for the average number retained of the initially $`2^j`$ coefficients for $`j=1\mathrm{}6`$ the respective values 1.9, 3.8, 7.4, 11.6, 6.9, and 3.4. While almost all of the coefficients for $`j3`$ are needed, a decreasing fraction is required to describe the spectra at higher resolution.

The distributions are insensitive to the signal–to–noise ratio, as shown in Fig. 5. Except for the lowest ratio of 10, the curves coincide, showing that they may be measured accurately even for a varying signal–to–noise ratio in a spectrum, provided it is not too low.

To demonstrate that the wavelet coefficient distributions may be used to discriminate between different predictions for the statistical properties of the Ly$`\alpha `$forest, a second set of Monte Carlo realizations with alternative column density and Doppler parameter distributions is generated. The parameters adopted are those reported by Hu et al. They found that the forest statistics are consistent with an $`\mathrm{I}`$column density distribution with a slope of $`1.5`$ for clouds with $`12.3<\mathrm{log}N_{\mathrm{HI}}<14.5`$ and a Gaussian Doppler parameter distribution with mean 28 $`\mathrm{km}\mathrm{s}^1`$, standard deviation 10 $`\mathrm{km}\mathrm{s}^1`$, and a sharp cut–off below 20 $`\mathrm{km}\mathrm{s}^1`$. The resulting average Doppler parameter is 37 $`\mathrm{km}\mathrm{s}^1`$. The simulation is normalized to the same line density per unit redshift at $`z=3`$ as found by Hu et al. , but the column density distribution is extended to $`\mathrm{log}N_{\mathrm{HI}}=16`$ to be consistent with the previous set of simulations. The resulting wavelet coefficient distributions are compared with those from the previous simulation for $`j=3,4,`$ and 5 in Fig. 6. The wavelet coefficients are able to distinguish between the two models.

The Kolmogorov–Smirnov test may be used to assess the probability that the wavelet coefficients of a measured spectrum match a given distribution for each resolution level $`j`$. The most stringent test, however, is given by combining the probabilities for all the distributions. Because any given absorption feature may be expected to affect the coefficients at more than a single resolution level $`j`$, it is possible that the coefficients corresponding to a given set of nested blocks for different $`j`$ (see Fig. 1) may be correlated. In this case, the probabilities of matching the various distributions may not be combined as if they were independent. To determine the degree to which the distributions may be treated as independent, the correlations are measured for coefficients between the various levels $`j`$ corresponding to the same hierarchy of blocks, and then averaged over all the hierarchies, for a set of Monte Carlo realizations using the fiducial forest model. The results are shown in Table 1. (The level $`j=0`$ refers to the correlations with the pair of coefficients corresponding to the course scale averages.) A signal–to–noise ratio of 50 is assumed, and a cut-off in the coefficients is applied to ensure $`\chi _{\mathrm{red}}^2=1`$. The error on the correlations is $`0.1`$%. Although the correlations are small, they are not absent. They are sufficiently small, however, that treating the probabilities for the different distributions as independent should be an adequate approximation for model testing.

### 2.4 Data compression

One of the key features of wavelets is their ability to compress data. Figs. 2 and 3 show that it is possible to fit a spectrum using only a subset of the wavelets used in its DWT at a statistically acceptable level ($`\chi _{\mathrm{red}}^2=1`$), without significantly degrading the information content of the spectrum as measured by the wavelet entropy. This suggests that filtering the spectrum in this way may provide a usable spectrum that is relatively noiseless and suitable for absorption line fitting.

This is illustrated by performing Voigt profile fitting to Monte Carlo realisations of the fiducial line model, with an assumed signal–to–noise ratio of 50. A wavelet filtered representation of each realised spectrum is generated with coefficients truncated to give a reduced $`\chi _{\mathrm{red}}^2=1`$ for the difference between the original and wavelet filtered spectra. This corresponds on average to retaining only 30% of the full set of coefficients. A representative spectrum is shown in Fig. 7.

Absorption lines are then identified in the filtered spectrum and fit using AutoVP. The results of $`10^4`$ realisations are shown in Figs. 8 and 9. Also shown are the distributions obtained from AutoVP using the original spectra with no wavelet filtering applied. The distributions are nearly identical. A negligible loss is incurred in the recovery of the line parameters despite the exclusion of 70% of the information in the original spectra.

## 3 Application to Q1937–1009

In this section, the Discrete Wavelet Transform is used to analyse the Ly$`\alpha `$forest as measured in the $`z=3.806`$ QSO Q1937–1009. The spectrum was taken with the Keck HIRES at a resolution of $`8.5`$ $`\mathrm{km}\mathrm{s}^1`$(Burles & Tytler 1997). The signal–to–noise ratio per pixel was $`50`$. The spectrum covers the range between Ly$`\alpha `$and Ly$`\beta `$ in the QSO restframe. (The region analysed is restricted to the redshift interval $`3.055<z<3.726`$ to avoid any possible influence by the QSO.)

The distribution of wavelet coefficients is shown in Fig. 10 for $`j=2`$, 3, 4, and 5. As in Fig. 4, the high frequency coefficients are generally smaller than at lower frequencies, indicating that the fluctuations that dominate the spectra have been resolved.

The cumulative distributions of the coefficients are compared with the predicted distributions for the fiducial model in Fig. 11. The predicted distributions were generated by simulating spectra with the same pixelization, resolution, signal–to–noise ratio and wavelength coverage as for the measured spectrum of Q1937–1009. An increase in line density per unit redshift proportional to $`(1+z)^{2.6}`$ (Kirkman & Tytler) was included to match to the redshift range of Q1937–1009. While the distributions generally agree well, a large variation is found for $`j=4`$, corresponding to fluctuations on the scale of $`1734`$ $`\mathrm{km}\mathrm{s}^1`$, suggesting some differences from the line model of Kirkman & Tytler. Effects neglected in the simulations that could produce a difference are the presence of metal systems and redshift correlations between the Ly$`\alpha `$absorption systems. The changes that would be produced, however, are most likely small: the number of metal systems is small, and the correlations appear weak or absent (Meiksin & Bouchet 1995; Kim et al. 1997). Still, the sensitivity of the wavelet coefficient distributions to these effects may be worth more careful consideration.

## 4 Summary

Wavelets may be usefully employed to provide a statistical characterizaton of the absorption properties of the Ly$`\alpha `$forest. An approach is presented that performs a multiresolution analysis of the forest using the Discrete Wavelet Transform of the QSO spectrum. The transform decomposes the local frequency dependence of the light fluctuations into an orthogonal hierarchy of basis functions, the wavelets. It is found that in spectra of better than 10 $`\mathrm{km}\mathrm{s}^1`$resolution, most of the information of the spectrum is carried by the lower frequency wavelets. For a signal–to–noise ratio typical of even the highest quality spectra ($`S/N=10100`$), only 10–30% of the wavelets are required to provide a statistically acceptable description of the spectrum, corresponding to a data compression factor of 3–10. It is shown that a Voigt profile line analysis performed on the wavelet filtered spectra yields nearly identical line parameter distributions as obtained from the original unfiltered spectra.

The distributions of the wavelet coefficients offer an alternative statistical description of the Ly$`\alpha `$forest while retaining information on the line widths. It is demonstrated that the correlations of coefficients between different levels in the wavelet hierarchy are weak (a few percent or smaller). Consequently, each of the distributions may be treated as statistically independent to good approximation.

The method is applied to a Keck HIRES spectrum of Q1937–1009. The wavelet coefficient distributions behave qualitatively similarly to those found in Monte Carlo simulations based on the line parameter distributions reported by Kirkman & Tytler. The measured distributions, however, show some differences on the scale $`1734`$ $`\mathrm{km}\mathrm{s}^1`$.

The results demonstrate that Multiresolution Analysis using the Discrete Wavelet Transform provides an alternative objective, easily automated procedure for analysing the Ly$`\alpha `$forest suitable for basing a comparison between the measured properties of the Ly$`\alpha `$forest and the predictions of numerical models.

The author thanks S. Burles and D. Tytler for kindly providing the spectrum of Q1937–1009, and R. Davé for permission to use AutoVP.


text:
# Non perturbative renormalization group potentials and quintessence

\[

## Abstract

New solutions to the non preturbative renormalization group equation for the effective action of a scalar field theory in the Local Potential Approximation having the exponential form $`e^{\pm \varphi }`$ are found. This result could be relevant for those quintessence phenomenological models where this kind of potentials are already used, giving them a solid field theoretical derivation. Other non perturbative solutions, that could also be considered for the quintessence scenario, are also found. Apart from this particular cosmological application, these results could be relevant for other models where scalar fields are involved, in particular for the scalar sector of the standard model.

pacs 98.80.Cq, 11.10.Hi , 11.10.Gh

\] One of the most challenging problems in fundamental physics has been the search for a theoretical argument, hopefully a symmetry, that could explain the vanishing of the cosmological constant. Recent observations from high redshift supernovae, combined with the data on the fluctuation of the cosmic microwave background, have changed our perspective. Apparently $`\mathrm{\Omega }_M`$, the ratio of the barionic and cold dark matter density to the critical density, is about $`\frac{1}{3}`$. This means that either the universe is open or the missing energy is provided by some new form of matter. The simplest candidate is a cosmological constant term. Alternative models where the missing energy is given by a scalar field slowly rolling down its potential, the “quintessence”, have recently attracted lot of attention.

For a truly constant vacuum energy term the “old” problem of explaining the vanishing of the cosmological constant is replaced by the equally difficult one of explaining why it has the small observed value of about $`(3\times 10^3ev)^4`$. In the quintessence scenario the equivalent problem is the so called “coincidence problem”. The matter and the scalar fields evolve differently but we observe today an order of magnitude “coincidence” between the matter energy density $`\mathrm{\Omega }_M`$ and the quintessence energy density $`\mathrm{\Omega }_\varphi `$ that requires a fantastic fine tuning of the initial conditions. The notion of “tracker field” , a quintessence scalar field that evolves to an attractor solution during its rolling down, has been introduced to circumvent this problem. For a wide range of initial conditions the attractor is stable. An explanation for the coincidence can be obtained this way and this has been advocated as an argument in favour of the quintessence scenario. Recently another interesting argument has been given in the framework of the brane-world picture.

Even restricting ourselves to tracker fields there is still a large amount of arbitrariness in the possible form of this potential. Different proposals  have been essentially based on their capability to reproduce the observational data. Also interesting attempts have been done  to derive its form from particle physics models, but they still leave the door open to several possibilities and it is not possible to discriminate between the different proposals.

In this Letter I show that the renormalization group equation for the effective action of a single component scalar field theory in the Local Potential Approximation (LPA) possesses non perturbative solutions in addition to the well known perturbative one. These are of the form of exponential potentials that are among the favoured candidates in the quintessence scenario. I think that this result can give a strong motivation to these potentials. I will comment more on this point later.

The exact renormalization group equation for the effective action has been found in Ref.. By considering only the potential term, an approximation to this equation (the LPA) is obtained, and for a single component scalar field theory in $`d=4`$ dimensions it reads:

$$k\frac{}{k}U_k(\varphi )=\frac{k^4}{16\pi ^2}ln\left(\frac{k^2+U_k^{^{\prime \prime }}(\varphi )}{k^2+m^2}\right).$$
(1)
Here $`k`$ is the current scale, $`m^2`$ is a constant with dimension $`(\mathrm{mass})^2`$, $`U_k`$ is the potential at the scale $`k`$ and the  means derivative with respect to the field $`\varphi `$.

Eq.(1) is a non perturbative evolution equation for $`U_k`$. It is immediate to see that the k-independent solution to Eq.(1), i.e. the fixed point potential, is the trivial gaussian one:

$$U_f(\varphi )=\frac{1}{2}m^2\varphi ^2+\alpha \varphi +\beta .$$
(2)
One simple way to show that Eq.(1) admits the well known perturbative solution is as follows. Let’s consider first a small deviation of the potential $`U_k(\varphi )`$ around $`U_f(\varphi )`$,

$$U_k(\varphi )=U_f(\varphi )+\delta U_k(\varphi ).$$
(3)
We develop now the logarithm in Eq.(1) in powers of $`\delta U_k`$ and expand $`\delta U_k`$ in powers of the field (for the sake of simplicity we consider a potential with the $`Z(2)`$ symmetry $`\varphi \varphi `$)

$$\delta U_k(\varphi )=\frac{\lambda _2(k)}{2}\varphi ^2+\frac{\lambda _4(k)}{4!}\varphi ^4+\mathrm{}.$$
(4)
At any order in $`\delta U_k`$ we obtain this way an infinite system of equations for the coupling constants. By truncating this system to the first two equations for $`\lambda _2(k)`$ and $`\lambda _4(k)`$ and solving iteratively up to the second order in $`\delta U_k`$  we get the known perturbative one loop RG flow for the coupling constants. Extrapolating down to $`k=0`$, this flow identifies the gaussian fixed point as infrared stable, i.e. we have the standard result about the triviality of the theory.

We ask now the question about the possibility of having non perturbative solutions. The only k-independent solution to Eq.(1), i.e. the only fixed point potential, is the gaussian one found in Eq.(2).

We want to linearize now Eq.(1) around $`U_f`$ and look for a small but non perturbative $`\delta U_k`$. We have (for notational simplicity I write $`U_k`$ rather than $`\delta U_k`$):

$$k\frac{}{k}U_k(\varphi )=\frac{k^4}{16\pi ^2}\frac{1}{k^2+m^2}U_k^{^{\prime \prime }}(\varphi ).$$
(5)
There is a class of non perturbative solutions to Eq.(5) that is very easy to find. Let’s seek for solutions of the form:

$$U_k(\varphi )=f(k)g(\varphi ).$$
(6)
Once inserted in Eq.(5), the ansatz (6) gives:

$$\frac{d^2g(\varphi )}{d\varphi ^2}\frac{1}{\mu ^2}g(\varphi )=0$$
(7)
$$\frac{df(k)}{dk}+\frac{1}{16\pi ^2\mu ^2}\frac{k^3}{k^2+m^2}f(k)=0,$$
(8)
where $`\mu ^2`$ is any constant with dimension $`(\mathrm{mass})^2`$ that allows to separate Eq.(5) in the two Eqs. (7) and (8).

Solving now these equations is a simple exercise. For positive values of the constant $`\mu ^2`$ the solutions to Eq.(5) have the form

$$U_k(\varphi )=M_1^4e^{\frac{1}{32\pi ^2\mu _1^2}\left(k^2m^2\mathrm{ln}\frac{k^2+m^2}{m^2}\right)}e^{\frac{\varphi }{\mu _1}}$$
(9)
$$U_k(\varphi )=M_2^4e^{\frac{1}{32\pi ^2\mu _2^2}\left(k^2m^2ln\frac{k^2+m^2}{m^2}\right)}e^{+\frac{\varphi }{\mu _2}},$$
(10)
for arbitrary values of the mass dimension constants $`M_i`$ and $`\mu _i`$.

For negative values of $`\mu ^2`$ (calling $`\overline{\mu }^2=\mu ^2>0`$), the solutions to Eq.(5) have the form:

$$U_k(\varphi )=M_3^4e^{+\frac{1}{32\pi ^2\overline{\mu }_3^2}\left(k^2m^2ln\frac{k^2+m^2}{m^2}\right)}\mathrm{cos}\left(\frac{\varphi }{\overline{\mu }_3}\right)$$
(11)
$$U_k(\varphi )=M_4^4e^{+\frac{1}{32\pi ^2\overline{\mu }_4^2}\left(k^2m^2ln\frac{k^2+m^2}{m^2}\right)}\mathrm{sin}\left(\frac{\varphi }{\overline{\mu }_4}\right).$$
(12)
For solutions of the kind (9) and (10) the gaussian potential $`U_f(\varphi )`$ is an ultraviolet fixed point. From that we immediately understand the completely different nature of these solutions with respect to the perturbative one.

As already mentioned these exponential potentials $`e^{\pm \varphi }`$, as well as linear combinations of them, are among the favourite candidates for the quintessence scenario.

The attempts that have been done to derive different forms of quintessence potentials all started from some sort of “fundamental” higher energy model. For example, inverse power-law potentials have been motivated from Supersymmetric QCD . For another derivation of the same kind of potentials see. Exponential potentials of the form found above arise naturally in several higher energy/higher dimensional theories. As we simply don’t know the theory that describes our world at very high energies, the use of phenomenologically motivated potentials is certainly well justified. In that respect indications coming directly from the effective theory of the quintessence field should be considered as very welcome. Actually, whatever the structure of the higher energy/higher dimensional theory, i.e. whatever the fundamental origin of the scalar quintessence field is, the “low energy” effective theory for this field should be very well described by Eq.(1). This happens because the higher energy degrees of freedom decouple from the quintessence field. There could still be the problem of the coupling of this field to ordinary matter, i.e. to ordinary standard model fields. As the long range forces that these couplings would generate are not observed, we can suppose that they are suppressed through some mechanism as for instance the one proposed in Ref.. In this case the renormalization group Eq.(1) gives the flow equation for the quintessence field irrespectively of the nature of the higher energy theory and the above results (9) and (10) give a solid motivation from the “low energy side” to the phenomenological exponential potentials $`e^{\pm \varphi }`$.

Concerning the solutions of the kind (11) and (12), we see that the gaussian potential $`U_f(\varphi )`$ is neither an infrared nor an ultraviolet fixed point for them, even though we can multiply these solutions times a small dimensionless constant so that they still make sense as linearized solutions of the Eq.(1). I want to mention here that cosine potentials have also been considered as possible quintessence candidates.

We want to seek now for other non perturbative solutions. Actually there is at least another class of such solutions that can be easily found. To see that, we switch first to the dimensionless form of Eq.(1). If we define the dimensionless field $`\phi `$, the dimensionless scale parameter $`t`$ and the dimensionless potential $`v(\phi ,t)`$ from

$$\phi =\frac{1}{4\pi }\frac{\varphi }{k},t=\mathrm{ln}\frac{\mathrm{\Lambda }}{k},U_k(\varphi )=\frac{k^4}{16\pi ^2}v(\phi ,t),$$
(13)
where $`\mathrm{\Lambda }`$ is a boundary value for $`k`$, Eq.(1) becomes:

$$\frac{v}{t}+\phi \frac{v}{\phi }4v=ln\left(\frac{1+\frac{^2v}{\phi ^2}}{1+\frac{m^2}{\mathrm{\Lambda }^2}e^{2t}}\right).$$
(14)
The dimensionless potential $`v_f(\phi ,t)`$ that corresponds to the gaussian potential $`U_f(\varphi )`$ of Eq.(2) is:

$$v_f(\phi ,t)=\frac{1}{2}\frac{m^2}{\mathrm{\Lambda }^2}e^{2t}\phi ^2+\frac{4\pi \alpha }{\mathrm{\Lambda }^3}e^{3t}\phi +\frac{16\pi ^2\beta }{\mathrm{\Lambda }^4}e^{4t},$$
(15)
and solves the equation:

$$\frac{v}{t}+\phi \frac{v}{\phi }4v=0.$$
(16)
Actually Eq.(14) is more often written as:

$$\frac{v}{t}+\phi \frac{v}{\phi }4v=ln\left(1+\frac{^2v}{\phi ^2}\right),$$
(17)
i.e. by setting $`m^2=0`$. This simply corresponds to choose a massless fixed point potential and from now on we also restrict ourselves to this case. We consider now a small fluctuation around the fixed point potential,

$$v(\phi ,t)=v_f(\phi ,t)+\delta v(\phi ,t),$$
(18)
and linearize Eq.(17) around $`v_f`$ (again we write $`v(\phi ,t)`$ rather than $`\delta v(\phi ,t)`$) to get:

$$\frac{v}{t}+\phi \frac{v}{\phi }4v=\frac{^2v}{\phi ^2}.$$
(19)
By following the same strategy as before we look for solutions to Eq.(19) with factorized $`t`$ and $`\phi `$ dependence:

$$v(\phi ,t)=f(t)g(\phi ).$$
(20)
Inserting the ansatz (20) in Eq.(19) we have:

$$\frac{d^2g(\phi )}{d\phi ^2}\phi \frac{dg(\phi )}{d\phi }+4g(\phi )=\alpha g(\phi )$$
(21)
$$\frac{df(t)}{dt}=\alpha f(t),$$
(22)
where $`\alpha `$ is an arbitrary dimensionless constant that allows to separate Eq.(19) in the two Eqs. (21) and (22).

Equation (22) is trivially solved and gives:

$$f(t)=Ae^{\alpha t},$$
(23)
where $`A`$ is the integration constant.

The solution to Eq.(21) can be found by series. Writing

$$g(\phi )=\underset{n=0}{\overset{\mathrm{}}{}}c_n\phi ^n,$$
(24)
inserting (24) in Eq.(21) and exploiting the recurrence relations between the coefficients $`c_n`$, we get the two linearly independent solutions ($`a=\alpha 4`$):

$$g_1(\phi )=c_0\left(1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{\Pi }_{j=1}^n(a+2j2)}{(2n)!}\phi ^{2n}\right)$$
(25)
$$g_2(\phi )=c_1\left(\phi +\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{\Pi }_{j=1}^n(a+2j1)}{(2n+1)!}\phi ^{2n+1}\right)$$
(26)
where I have explicitly kept $`c_0`$ and $`c_1`$, the first two coefficients in Eq.(24), as integration constants.

After some trivial algebra we see that Eqs.(25) and (26) can be written in terms of confluent hypergeometric functions and the general solution to Eq.(21) takes the compact form

$$g(\phi )=c_0M(\frac{a}{2},\frac{1}{2},\frac{\phi ^2}{2})+c_1\phi M(\frac{a+1}{2},\frac{3}{2},\frac{\phi ^2}{2}).$$
(27)
We recall here that the confluent hypergeometric function $`M(a,b,x)`$ is defined as

$$M(a,b,x)=1+\frac{a}{b}x+\frac{1}{2!}\frac{a(a+1)}{b(b+1)}x^2+\mathrm{}.$$
(28)
Of course once dimensionless solutions to Eq.(19) are known, it is a trivial exercise to reconstruct dimensionfull potentials from Eq.(13). Some of the potentials obtained this way could also be considered as quintessence candidates.

It is worth to compare the general solution (27) to Eq.(21) with a non perturbative result that has been obtained some years ago by Halpern and Huang following a different but essentially equivalent approach. Searching for alternatives to the trivial $`\varphi ^4`$ theory, the authors expanded the potential in even powers of the field and derived an infinite system of differential equations for the coupling constants. Looking for new eigendirections in the parameter space, they actually ended with the solution $`g_1(\phi )`$ to Eq.(21). More precisely, after making the trivial changes to match the two different notations and restricting ourselves to consider the $`N=1`$ and d=4 case as in the present paper, we can immediately check that eq.(49) of Ref. coincides with the $`g_1(\phi )`$ solution above. As these authors considered potentials containing only even powers of the field, obviously they could only get the solution $`g_1(\phi )`$.

At a first sight it could seem that the solution $`g_2(\phi )`$ should be discarded as it contains odd powers of the field and as such it is unbounded from below. We would conclude in this case that the only physically acceptable general solution to Eq.(19) is $`g_1(\phi )`$, i.e. the Halpern-Huang result. But this is not always true. It is immediate to see from Eq.(27) that for all the positive or negative integer odd values of $`\alpha `$ such that $`\alpha <4`$, $`g_2(\phi )`$ is a polynomial in $`\phi `$ and it can be combined with $`g_1(\phi )`$ to give a bounded from below potential. It is also easy to give examples where even when the hypergeometric function in $`g_2(\phi )`$ keeps all the infinite terms, still a linear combination of $`g_1(\phi )`$ and $`g_2(\phi )`$ gives a bounded from below potential. Take for instance the case $`\alpha =5`$ for which both $`g_1`$ and $`g_2`$ are not polynomials. From the asymptotic behaviour of the hypergeometric function we easily see that $`c_0`$ and $`c_1`$ can be chosen in such a way that the resulting potential is bounded from below. The class of physically acceptable potentials, that are solutions of the Eq.(21), is larger than that spanned by $`g_1`$ only.

I should mention at this point that an attempt to solve Eq.(19) has been recently made in Ref. where the solution $`v(\phi ,t)=e^{5t}e^{\frac{\phi ^2}{2}}`$ is presented. The author says that this solution is asymptotically similar to those of Ref. but that the connection between the two results is unclear to him. From Refs. and from Eq.(27) above (setting $`c_1=0`$) we immediately see that this solution is just the Halpern-Huang result for $`\alpha =5`$. It is also a trivial exercise to verify that the other solutions presented in (see Eq.(21) of that paper) are particular cases of the general solution, Eq.(27), obtained for the integer values $`\alpha =6,7,8,\mathrm{}`$.

To summarize in this letter I have presented new solutions to the LPA of the exact renormalization group equation for the effective action of a single component scalar field theory. These potentials have the exponential form $`e^{\pm \varphi }`$ and have been recently used in phenomenological quintessence models. As the effective theory for the quintessence field should be governed by Eq.(1) whatever its higher energy origin, I argue that the fact that they arise as solutions to the renormalization group equation gives a derivation of these potentials that is alternative and complementary to those based on higher energy theories and should allow to discriminate between different proposals. I have also presented other solutions to the LPA renormalization group equation some of which already partially known.

Apart from the application to the particular cosmological problem suggested in the present paper, I think that these results could also be relevant in other frameworks where scalar fields play a role. For the non perturbative solutions to Eq.(19) of the kind $`e^{\pm \varphi }`$, the gaussian potential is an UV stable fixed point. The same is true for those solutions of the kind (27) when $`\alpha >0`$. This result is opposite to the perturbative one and its implications in particle physics models are certainly worth to explore. The scalar sector of the standard model is itself an open problem. It is not a priori clear whether the existence of these solutions could have some relevance for the theory. I hope to come back to this issue in a future paper. Work is in progress in this direction.