text:
# Scanning Lidar Based Atmospheric Monitoring for Fluorescence Detectors of Cosmic Showers

## 1 Introduction

Contemporary fluorescence experiments (Fly’s Eye , HiRes , P. Auger ) studying cosmic rays with energies near $`10^{20}e`$V detect fluorescence light produced along the air-shower volume. As a cosmic-ray induced particle cascade develops in the atmosphere it dissipates much of its energy by exciting and ionizing air molecules. The cascade particles along the shower axis are limited to a narrow lateral distribution where excited nitrogen molecules are fluorescing in the near-UV spectral band. The efficiency of the process is rather small. However, the vast number of emitting particles in high energy shower makes this source of radiation highly significant. Ultimately, this electromagnetic (EM) cascade dissipates much of the primary particle’s energy. Fluorescence light is emitted isotropically with an intensity proportional to the number of charged particles in the shower. EM component and hence the total number of low energy EM particles is in turn fairly accurately proportional to the energy of the primary particle. Thus, the calorimetric measure of the total EM shower energy  is proportional to the integral of EM particle density $`N_{\text{em}}`$ along the shower direction $`x`$,

$$E_{\text{em}}=KN_{\text{em}}(x)x$$
(1)
with $`K2.2`$M$`e`$V cm2/g, where $`x`$ is measured in units of longitudinal air density (g/cm2). $`E_{\text{em}}`$ is a lower bound for the energy of the primary cosmic ray. The lower portion of shower development is usually obscured by the ground so that EM cascade reaching below ground is included by fitting a functional form to the observed longitudinal profile and integrating the function past surface depth. The number of photons $`N_{\text{ph}}`$ reaching fluorescence detector (FD) is proportional to EM particle density $`N_{\text{em}}(x)`$ at the point of production $`x`$, so that in turn

$$N_{\text{em}}(x)\frac{N_{\text{ph}}R^2(x)}{T(x)},$$
(2)
with $`R(x)`$ being distance between shower point $`x`$ and FD. Light originating within the shower is certainly affected by the absorption and scattering on molecules and aerosols in the atmosphere. The number of detected photons is thus reduced due to non-ideal atmospheric transmission $`T(x)<1`$, where

$$T(x)=\mathrm{exp}\left[_0^x\alpha (r)r\right]=\mathrm{}^{\tau (x)},$$
(3)
with $`\alpha (r)`$ volume extinction coefficient along the line-of-sight, and $`\tau (x)`$ the resulting atmospheric optical depth (OD) to the shower point $`x`$.

In this sense, the atmosphere can be treated as an elementary-particle detector. However, weather conditions change the atmospheric transmission properties dramatically resulting in strongly time-dependent detection efficiency. Therefore, an absolute calibration system for fluorescence light absorption is an essential part of FD .

Eq. (3) is a basis for the fluorescence-detector energy calibration. Apart from the Pierre Auger Observatory, all existing fluorescence-based experiments have suffered from a lack of a sufficient atmosphere monitoring system. Weather conditions in a desert-like atmosphere were expected to be stable enough, so that the standard attenuation-length profile should suffice to reconstruct the total EM shower energy, Eq. (1), with a controllable precision. However, this has turned out not to be the case, especially for rare events with energies above $`10^{19}e`$V. In addition, the energy reconstruction is obscured by Čerenkov radiation in the lower part of the air-shower which can not be separated from fluorescence. For more than a half of the highest-energy events measured so far the atmosphere properties are not known well enough to accurately reconstruct the primary energy. To be able to provide adequate calibration, one has to measure attenuation at the time of the event in the whole region of the air-shower.

There is also a systematic discrepancy when comparing the cosmic-ray spectra of fluorescence experiments with ground arrays. Their compatibility can be established only with a shift in energy and flux of one or the other. As in the case of fluorescence detectors, ground arrays have their own problems with the energy determination and are much more dependent on air-shower simulations. At present, it is not known whether the discrepancy is due to fluorescence-detector or ground-array method, or both. Therefore, it is of utmost importance to have in situ atmosphere monitoring system which is working coherently with a fluorescence detector.

To lower primary cosmic ray energy uncertainties, the volume extinction coefficient $`\alpha (r)`$ thus has to be well estimated over almost whole detection volume of FD. In the case of the Pierre Auger Observatory, the detection volume corresponds to ground area of 3000 km2 and height of $`15`$ km.

The paper is organized as follows. The first two sections are devoted to introductory material on lidar measurement technique and a description of our specific experimental setup. Then the atmospheric model for simulation of lidar signals is presented and the signals are evaluated by two well established inversion methods. Results of the inversions are compared to the input model and conclusions on their applicability are drawn. Next, improved approaches to FD calibration based on scanning lidar system are proposed and evaluated on real data obtained with our experimental setup.

## 2 Lidar system

One of the most suitable calibration setups for FD is the backscattering lidar system, where a short laser light pulse is transmitted from FD position in the direction of interest. With a mirror and a photomultiplier tube, backscattered light is collected and recorded as a function of time, i.e. as a function of backscatter distance. Note that light from the lidar source traverses both directions, so that in case of matching laser and fluorescence light wavelength, OD for lidar light sums to twice the OD for fluorescence. The lidar equation  describes the received laser power $`P(r)`$ from distance $`r`$ as a function of volume extinction coefficient $`\alpha (r)`$ and backscattering coefficient $`\beta (r)`$,

$$P(r)=P_0\frac{ct_0}{2}\beta (r)\frac{A}{r^2}\mathrm{}^{2\tau (r)}.$$
(4)
$`P_0`$ is the transmitted laser power and $`A`$ is an effective receiving area of the detector, proportional to the area of the mirror and proportional to an overlap between its field of view with the laser beam. $`t_0`$ is laser pulse duration. As seen from Eq. (2), measurement precision of $`\alpha `$ and corresponding $`\tau `$ directly influences the precision of primary particle energy estimation.

Simple as it may look, the lidar equation (4) is nevertheless difficult to solve for two unknown variables, $`\alpha (r)`$ and $`\beta (r)`$. All existing analysis algorithms (Klett , Fernald , and their respective variations) reviewed in one of the following sections are based on an experimental setup with static beam direction. This leads to ambiguity in determination of $`\alpha (r)`$ and $`\beta (r)`$ which can not be resolved without additional assumptions about atmospheric properties. At the FD experimental sites the atmosphere can be assumed to be almost horizontally invariant. In this case, there is an additional constraint when comparing signals coming from different directions, which solves the lidar equation for $`\alpha (r)`$ and $`\beta (r)`$ unambiguously. Even at this point, it can be stated that the need for steerable (scanning) lidar setup is unavoidable for proper solution of lidar equation.

## 3 Experimental setup

The lidar system used for verification of the analysis method is based on the Continuum MiniLite-1 frequency tripled Nd:YaG laser, which is able to transmit up to 15 shots per second, each with energy of 6 mJ and 4 ns duration (1.2 m). The emitted wavelength of 355 nm is in the $`300400`$ nm range of fluorescence spectrum. The receiver was constructed using $`80`$ cm diameter parabolic mirror with focal length of 41 cm. The mirror is made of aluminum coated pyrex and protected with SiO2.

The backscattered light is detected by a Hammamatsu R7400 photomultiplier with operating voltage up to 1000 V and gain of $`10^5`$ to $`10^6`$. To suppress background, a broadband UG-1 filter with 60% trasmitance at 353 nm and FWHM of 50 nm is used. The distance between laser beam and the mirror center is fixed to 1 m, and the system is fully steerable with $`0.1^{}`$ angular resolution.

The signal is digitized using a three-channel LICEL transient recorder TR40-160 with 12 bit resolution at 40 MHz sampling rate with 16k trace length combined with 250 MHz photon counting system. Maximum detection distance of the hardware is thus, with this sampling rate and trace length, set to 60 km. However, in real measurements, atmospheric features up to 30 km only are observed. LICEL is operated using a PC-Linux system through a National Instruments digital input-output card (PCI-DIO-32HS) with Comedi drivers  and a ROOT interface .

## 4 Lidar simulation with specific atmospheric model

In a low opacity atmosphere the attenuation and the backscattering coefficient can be written as a sum of contributions from two independent components,

$`\alpha (h)`$ $`=\alpha _\text{m}(h)+\alpha _\text{a}(h),`$ (5a)
$`\beta (h)`$ $`=P_\text{m}(180^{})\alpha _\text{m}(h)+P_\text{a}(180^{})\alpha _\text{a}(h).`$ (5b)
where $`\alpha _\text{m}`$ and $`\alpha _\text{a}`$ correspond to molecular and aerosol attenuation, respectively. The aerosol phase function $`P_\text{a}(180^{})`$ for backscattering has, apart from the wavelength, also a strong dependence on the optical and geometrical properties of the aerosol particles. Nevertheless, at wavelength of 355 nm, values in the range 0.025 and up to 0.05 sr-1 can be assumed  for aerosol phase function $`P_\text{a}(180^{})`$. The angular dependence of molecular phase function is defined by the Rayleigh scattering theory, where $`P_\text{m}(180^{})=3/8\pi \text{sr}^1`$.

For simulation purposes, the elevation dependence of the extinction coefficients is modelled as following,

$`\alpha _\text{m}(h)`$ $`={\displaystyle \frac{1}{L_\text{m}}}\mathrm{}^{h/h_\text{m}^0},`$ (6a)
$`\alpha _\text{a}(h)`$ $`={\displaystyle \frac{1}{L_\text{a}}}\{\begin{array}{cc}1,\hfill & h<h_\text{x}\hfill \\ \mathrm{}^{(hh_\text{x})/h_\text{a}^0},\hfill & hh_\text{x}.\hfill \end{array}`$ (6b)
where $`L_\text{m}`$ and $`L_\text{a}`$ are the molecular and aerosol attenuation lengths at ground level, and $`h_\text{m}^0`$ and $`h_\text{a}^0`$ are the molecular and aerosol scale height, respectively. An additional mixing height $`h_\text{x}`$ is set up for aerosols, assuming uniform concentration near the ground level and continuous transition into exponential vanishing for $`h>h_\text{x}`$. The following values of the parameters are used: $`L_\text{m}=15`$km, $`h_\text{m}^0=17.5`$km, $`L_\text{a}=2`$km, $`h_\text{x}=0.8`$km, and $`h_\text{a}^0=1.4`$km.

The atmospheric model, Eq. (6), serves as a testing ground for two widely used reconstruction methods presented in the next section. A comparison with reconstruction of the real atmosphere yields insight into the common problems of the lidar field. The Poissonian statistics of photon counting and multiplying, background noise, and effects of digitalization have been taken into account in the generation of the simulated lidar signals and under inspection match those observed in the real lidar power returns.

Model in Eq. (6) is a valid approximation to the atmospheric conditions found in real experiments. Although the vertical variation of aerosol and molecular densities is quite simple, the model still produces a non-trivial relation between total attenuation $`\alpha `$ and total backscattering coefficient $`\beta `$. Therefore, the dependence of $`\beta `$ as a function of $`\alpha `$, shown in Fig. 2, cannot be well approximated by some simple functional form.

## 5 Reconstruction of a 1D atmosphere

Concentrating on a single shot lidar measurement, the optical properties obviously have to be reconstructed in a 1D subspace of the atmosphere. Rewriting the lidar equation (4),

$$P(r)=B\frac{\beta (r)}{r^2}\mathrm{}^{2\tau (r)}$$
(7)
where the effective aperture of the system is gathered in the constant $`B`$, an auxiliary $`S`$-function can be introduced,

$$S(r)=\mathrm{ln}\frac{P(r)r^2}{P(r_0)r_0^2}=\mathrm{ln}\left[\beta (r)/\beta _0\right]2\tau (r;r_0).$$
(8)
Note that $`\tau (r;r_0)=_{r_0}^r\alpha (r^{})r^{}`$ corresponds to atmospheric OD between $`r_0`$ and $`r`$.

### 5.1 Klett inversion

Apart from the experimentally measured lidar power return $`P(r)`$, in Eq. (7) there are two unknown quantities, $`\beta `$ and $`\alpha `$ (or equivalently $`\tau `$), preventing the unique solution of the lidar equation. Nevertheless, a simple, and sometimes physically meaningful, assumption of proportionality between backscattering and extinction,

$$\beta (r)\alpha ^k(r),$$
(9)
allows for the transformation of the integral Eq. (8) to the corresponding Bernoulli’s differential equation with an existing analytical solution. Direct application of the solution (forward inversion) is numerically unstable, in some cases singular, and highly sensitive to the signal noise . Klett’s reformulation  of the solution (backward inversion) avoids these problems. The lidar backward inversion algorithm proceeds from the far point of the measured signal $`r_\text{f}`$ to the near end,

$$\alpha (r;\alpha _\text{f})=\frac{\mathrm{}^{S(r)/k}}{\mathrm{}^{S_\text{f}/k}/\alpha _\text{f}+\frac{2}{k}_r^{r_\text{f}}\mathrm{}^{S(r^{})/k}r^{}},$$
(10)
where $`S_\text{f}=S(r_\text{f})`$, and $`\alpha _\text{f}=\alpha (r_\text{f})`$ is an estimate for the attenuation at the far end of the data set. The reconstructed attenuation $`\alpha (r;\alpha _\text{f})`$ is still a one-parameter function of the unknown boundary attenuation value $`\alpha _\text{f}`$, so that independent measurement, or suitable approximation, is needed at the reference distance $`r_\text{f}`$. OD can be expressed directly from Eq. (10),

$$\tau (r;r_0,\alpha _\text{f})=\frac{k}{2}\mathrm{ln}\left[\frac{k\mathrm{}^{S_\text{f}/k}+2\alpha _\text{f}_{r_0}^{r_\text{f}}\mathrm{}^{S(r^{})/k}r^{}}{k\mathrm{}^{S_\text{f}/k}+2\alpha _\text{f}_r^{r_\text{f}}\mathrm{}^{S(r^{})/k}r^{}}\right].$$
(11)
Klett’s inversion method depends rather strongly on the assumed power law proportionality in Eq. (9). In Fig. 3, a failure of this approximation is demonstrated for the specific atmospheric model used for our simulations. The local value of the exponent,

$$k=\frac{\alpha }{\beta }\frac{\beta }{r}\left[\frac{\alpha }{r}\right]^1,$$
(12)
is shown to possess substantial range dependence. The main reason for failure of the power law proportionality stems from the inequality of the molecular and aerosol phase functions, $`P_\text{m}(180^{})`$ and $`P_\text{a}(180^{})`$, rendering the $`\alpha `$ and $`\beta `$ relationship dependent on the particular magnitude of both quantities, and consequently range dependence (see Fig. 2). Therefore, the best value of $`k`$ must be chosen using some ad hoc criterion.

Analyzing results in Fig. 4, presenting Klett inversion of simulated lidar signals, it seems that the closest reconstruction of the model profile is achieved with $`k0.5`$. From Fig. 3, showing local exponent $`k`$ obtained with use of Eq. (12), it can be seen that $`k0.5`$ is observed only in small interval around $`4\text{km}`$ whereas at other places it is substantially larger. For $`r>8\text{km}`$, dominated by molecular scattering it, slowly approaches the value of 1, most commonly adopted in the literature. Nevertheless, as can be seen in Fig. 5, reconstruction of OD with $`k=1`$ totally fails to reproduce the correct answer. Surprisingly, in case of this specific atmospheric model the most authentic result is obtained with $`k0.5`$.

Another drawback of the Klett’s method is estimation of the extinction $`\alpha _\text{f}`$ at the far end of the lidar return. In the case that $`r_\text{f}`$ corresponds to a highly elevated point, approximation $`\alpha _\text{f}\alpha _\text{m}(r_\text{f})`$, i.e. the extinction at that point is dominated by the molecular scattering, yields quite reasonable results  with qualitative convergence to the correct $`\alpha `$-profile. In general, for optically dense atmosphere (e.g. presence of moderate haze) convergence of the Klett’s method is far more rapid as in clear, optically thin case. However, sites for FD are usually chosen at locations with clear and cloudless atmosphere. For horizontal lidar measurements (zenith angle $`\varphi =90^{}`$) in a horizontally invariant atmosphere, $`\alpha _\text{f}`$ can be estimated as the one that minimizes extinction deviations from a constant value , i.e. minimizes the functional $`_{r_0}^{r_\text{f}}[\alpha (r^{})\alpha _\text{f}]^2r^{}`$.

### 5.2 Fernald inversion

Since the concentration of the molecules depends solely on the thermodynamic parameters of the atmosphere, the Rayleigh scattering on molecules is modeled separately on a basis of the meteorological data. $`\alpha _\text{m}(r)`$ acquired in that way is inserted in Eq. (5). With an estimate for the molecular and aerosol backscattering phase fraction, $`F=P_\text{m}(180^{})/P_\text{a}(180^{})`$, and modified $`S`$-function

$$\stackrel{~}{S}(r)=S(r)+2(F1)_r^{r_\text{f}}\alpha _\text{m}(r^{})r^{},$$
(13)
the lidar equation can be solved for aerosol part $`\alpha _\text{a}(r)`$ following the same steps as in Klett’s version,

$$\alpha _\text{a}(r)=F\alpha _\text{m}(r)+\frac{\mathrm{}^{\stackrel{~}{S}(r)}}{\mathrm{}^{\stackrel{~}{S}_\text{f}}/\stackrel{~}{\alpha }_\text{f}+2_r^{r_\text{f}}\mathrm{}^{\stackrel{~}{S}(r^{})}r^{}},$$
(14)
with $`\stackrel{~}{\alpha }_\text{f}=F\alpha _\text{m}(r_\text{f})+\alpha _\text{a}(r_\text{f})`$ and $`\stackrel{~}{S}_\text{f}=\stackrel{~}{S}(r_\text{f})=S(r_\text{f})`$. In the same way OD is expressed as

$`\tau (r;r_0,\stackrel{~}{\alpha }_\text{f})=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}\left[{\displaystyle \frac{\mathrm{}^{\stackrel{~}{S}_\text{f}}+2\stackrel{~}{\alpha }_\text{f}_{r_0}^{r_\text{m}}\mathrm{}^{\stackrel{~}{S}(r^{})}r^{}}{\mathrm{}^{\stackrel{~}{S}_\text{f}}+2\stackrel{~}{\alpha }_\text{f}_r^{r_\text{m}}\mathrm{}^{\stackrel{~}{S}(r^{})}r^{}}}\right]+`$
$`+(1F){\displaystyle _{r_0}^r}\alpha _\text{m}(r^{})r^{}.`$ (15)
Note that the Fernald procedure relies on three independently supplied parameters: (i) an accurate estimate of molecular part of the scattering $`\alpha _\text{m}(r)`$ along the whole range of interest, (ii) total extinction at the far end $`\stackrel{~}{\alpha }_\text{f}`$, and (iii) proper approximation for phase fraction $`F`$. As predicted by the Mie theory, it is quite difficult to obtain reasonable values for the latter. As for $`\stackrel{~}{\alpha }_\text{f}`$, conclusions are similar to those of Klett’s $`\alpha _\text{f}`$.

In Fig. 6 Fernald’s inversion of simulated lidar return is shown for different input values of $`\alpha _\text{a}(r_\text{f})`$ that enter total extinction $`\stackrel{~}{\alpha }_\text{f}`$. For upward pointing lidar measurements vanishing aerosol concentration can be assumed at the far end of atmosphere, i.e. $`\alpha _\text{a}(r_\text{f})=0`$. To test the sensitivity of the reconstructed OD on this assumption, data sets with $`\alpha _\text{a}(r_\text{f})=\pm \alpha _\text{m}(r_\text{f})/2`$ and therefore $`\stackrel{~}{\alpha }_\text{f}=(F\pm 1/2)\alpha _\text{m}(r_\text{f})`$, are also plotted. $`P_\text{a}(180^{})=0.025`$ sr-1 is used in the expression for phase fraction $`F`$. Comparing to the Klett’s method, which does not separate aerosol and molecular scattering, it is not surprising that the variation of Fernald’s results on boundary parameters is somewhat weaker. Pinning the molecular part of scattering undoubtedly stabilizes OD profiles obtained. Nevertheless, Fernald’s inversion still relies heavily on additional external parameters that are usually difficult, if not impossible, to measure.

## 6 Horizontally invariant atmosphere

Fluorescence detectors for cosmic showers are usually placed at locations with specific atmospheric conditions. In case of the Pierre Auger Observatory, the FD cameras are covering the lower part of the atmosphere over an almost perfect $`3000\text{km}^2`$ plane $`1500\text{m}`$ above the sea level with remarkable fraction of cloudless days. Due to the high elevation and dry inland climate, an optically thin atmosphere is expected. But, as noted before, in this case convergence of Klett’s method is slower and can lead to erroneous estimates of OD. Based on that, and other peculiar problems of the well established lidar inversion methods, a new approach with fewer a priori or hard-to-estimate input parameters is needed. Since the lidar equation is not uniquely solvable, a minimal set of assumptions needed for inversion has to be reconsidered. For a typical FD site it is quite reasonable to assume weak horizontal variation of the atmospheric optical properties. That is even more true for the huge plane mentioned above, with hardly any changes in elevation and vegetation coverage. Since the FD is exclusively operating at night, only atmospheric conditions at that time have to be considered. The mean night wind speeds do not exceed $`12\text{km}/\text{h}`$ , so that particularly thin layer of aerosols close to the ground is expected. At night, it is also expected that there will be a low probability for formation of convective types of atmospheric instabilities.

### 6.1 Two-angle reconstruction

Under the moderate assumptions presented above, optical parameters of atmosphere that enter the lidar equation (7) can be assumed to possess only vertical variations, while being uniform and invariant in the horizontal plane.

Thus, it makes sense to rewrite the range dependent $`S`$-function in Eq. (8) in terms of height $`h`$ and geometric factor $`\xi =1/\mathrm{cos}\varphi =\mathrm{sec}\varphi `$, when lidar shots with zenith angle $`\varphi `$ are considered. The $`S`$-function becomes

$$S(h,\xi )=\mathrm{ln}\left[\beta (h)/\beta _0\right]2\xi \tau (h;h_0)$$
(16)
with “vertical” OD $`\tau (h;h_0)=_{h_0}^h\alpha (h^{})h^{}`$ and $`\beta _0=\beta (h_0)`$. After measuring two $`S`$-functions at different zenith angles $`\xi _1=1/\mathrm{cos}\varphi _1`$ and $`\xi _2=1/\mathrm{cos}\varphi _2`$ and height $`h`$, Eq. (16) can be solved for the vertical OD,

$$\tau (h)=\frac{1}{2}\frac{S(h,\xi _1)S(h,\xi _2)}{\xi _1\xi _2},$$
(17)
and the backscatter coefficient ratio,

$$\frac{\beta (h)}{\beta _0}=\mathrm{exp}\left[\frac{\xi _2S(h,\xi _1)\xi _1S(h,\xi _2)}{\xi _1\xi _2}\right].$$
(18)
Both quantities are directly proportional to the difference of two $`S`$-functions at the same height and different angles. Therefore, choosing a small separation between zenith angles, $`\xi _1=\xi `$ and $`\xi _2=\xi +\xi `$, a differential form of Eq. (17) can be written,

$$\tau (h)=\frac{1}{2}\frac{S}{\xi }|_h.$$
(19)
Equivalently, the differential form of Eq. (18) can be obtained,

$$\frac{\beta (h)}{\beta _0}=\mathrm{exp}\left[S(h,\varphi )\xi \frac{S}{\xi }|_h\right].$$
(20)
Note that the OD is in that way determined up to the additive constant, and the backscatter coefficient up to the multiplicative factor. Nevertheless, both values should satisfy $`S(h_0)=0`$ and $`\tau (h_0)=0`$.

Taking into account the Poissonian statistics of collected photons, and neglecting all other sources of measurement uncertainties, a relative error of the obtained OD at some height depends on the lidar system parameters,

$$\frac{\sigma _\tau }{\tau }=\frac{h/h_0}{2\tau \sqrt{N_0\stackrel{~}{\beta }}}\frac{1}{|\xi _1\xi _2|}\sqrt{\mathrm{}^{2\xi _1\tau }+\mathrm{}^{2\xi _2\tau }},$$
(21)
as well as relative error of backscatter coefficient

$$\frac{\sigma _\beta }{\beta }=\frac{h/h_0}{\sqrt{N_0\stackrel{~}{\beta }}}\frac{\sqrt{\xi _2^2\mathrm{}^{2\xi _1\tau }+\xi _1^2\mathrm{}^{2\xi _2\tau }}}{|\xi _1\xi _2|},$$
(22)
where $`N_0`$ is number of detected photons in the time interval corresponding to the power return from height $`h_0`$, and $`\stackrel{~}{\beta }=\beta /\beta _0`$.

In Fig. 7, an example of $`S`$-functions and their zenith angle variation is presented. All results are obtained from real lidar measurements performed during few November nights in a typical urban atmosphere (GPS location: 4604’35” N, 01429’05” E, 312 m above sea level). For fixed primary azimuth angle $`\varphi _1=0^{}`$ and three selected secondary angles $`\varphi _2=38^{}`$, $`42^{}`$, and $`47^{}`$ results for OD (Fig. 8), backscatter coefficient (Fig. 9), and $`\alpha \beta `$ diagram (Fig. 10) are obtained from corresponding $`S`$-functions in Fig. 7. Due to presence of a thin layer of optically thick haze at $`h3\text{km}`$, a drastic change in both OD and backscattering at that height is observed. Since OD is well determined only up to an additive constant, note that the variation of results for different $`\varphi _2`$ is easily produced by the inadequate determination of $`S_0`$, in other terms, by variation of atmospheric optical properties at $`h_0`$. Compatible with a scale height of $`18\text{km}`$, the variation of backscattering in Fig. 9 is slower as found in our model, generating a gradual but still comparable $`\alpha \beta `$ diagram in Fig. 10.

In Fig. 11, a logarithmic plot of the relative error in OD is presented for typical lidar system parameters. First, the angle is fixed to $`\varphi _1=0^{}`$ while the second one, $`\varphi _2`$, is varied from a vertical to an almost horizontal shot. It is hard to avoid the fact that minimum error is produced with evaluation of two quite considerably separated lidar shots, $`\varphi _270^{}`$. Even at moderate elevations $`h`$ this can amount to large spatial separations of the two points of lidar return, and thus the requirement of horizontal invariance easily broken. In the case of an atmosphere, that is slowly horizontally modulated, a more “local” approach to the OD problem is needed.

### 6.2 Multi-angle reconstruction

For the ideal atmosphere, with true horizontal invariance, the $`\xi `$ dependence of $`S`$-function is particularly simple,

$$S(h,\xi )=\mathrm{ln}[\beta (h)/\beta _0]2\xi \tau (h;h_0),$$
(23)
with the backscatter coefficient $`\mathrm{ln}[\beta /\beta _0]`$ as offset, and OD $`\tau `$ as the slope of the resulting linear function in $`\xi `$. Therefore, the optical properties of the atmosphere can be alternatively obtained from the analysis of the $`S`$-function behavior for scanning lidar measurements. Furthermore, disagreement of the measured $`S(\xi )`$ profiles from the linear form is a suitable criterion for detection of deviations from the assumed horizontal invariance of the atmosphere.

A generalization of the two-angle equations (19) and (20) to their differential counterparts strongly suggested this way of reconstruction of optical properties, the two-angle method being a mere two-point approximation of the linear function in Eq. (23). Taking into account quite substantial uncertainties in $`S(\xi )`$ for single angle, the linear fit trough many data points seems to yield superior results and the reconstruction is no longer limited to two lidar shots, well-separated in angle. The preferred horizontal invariance is not required to take place across huge atmospheric volumes (as in case of $`\varphi _1=0^{}`$ and $`\varphi _2=60^{}`$ shots), but has to be met only in relatively small arc of interest where the continuous lidar scan is performed.

In the opposite case, when slow variation of atmospheric properties in horizontal plane is allowed, Eq. (19) is similar enough to the renown 1D “slope method”, based on assumption of small variation of $`\beta (r)`$, or equivalently $`\beta /r0`$. Bear in mind that in method presented here the variation of $`\beta `$ with height can be of any magnitude, as long as there are only modest variations in the horizontal direction.

Opposite to Fig. 7, in Fig. 12 $`S`$-function profiles with respect to zenith $`\xi `$ are drawn for fixed heights, starting with $`h=3.2`$km and up to $`7`$km with $`633`$m step. Approximate linear behavior is observed in few arc intervals, with narrow bands of minute atmospheric shifts at $`\xi =1.15`$ and 1.38. Since these shifts in profiles disappear when lifting $`h_0`$ from $`3`$km to $`3.5`$km, they are obviously due to the distortions of atmosphere in the latter interval, feature already observed in Fig. 8.

In Fig. 13 the results of fitting and extraction of OD are similar to the ones in Fig. 8. Note that in both cases OD is obtained relative to the $`h_0`$ point, so that the results may differ up to some additive constant. Therefore, comparing both figures, it is more accurate to concentrate on the same span of OD within the $`3.5`$km to $`9`$km interval. Nevertheless, the range of OD results with acceptable error bars is with multi-angle method increased up to $`12`$km.

The relative error of OD in Fig. 14 is needed for correct estimation of shower energy uncertainty. It is kept below 6% even for the OD from the far points of the range, and below 3% for modest values of OD. Fig. 15, with values for $`\beta /\beta _0`$ should be compared to Fig. 9.

## 7 Conclusions

Inversion attempts of simulated lidar returns for atmosphere, modeled by Eqs. (6), show numerous drawbacks of established numerical methods. For instance, Klett’s and Fernald’s method of section 5.1 and 5.2 do not satisfy the specific requirements of FD calibration. While they may be useful for qualitative reconstruction of atmospheric properties (spatial haze/cloud distribution, cloud base etc.), they are not applicable for absolute assessment of atmospheric transmission properties. There are many reasons for this failure. One of them is certainly strong dependence of obtained inversions on presumed extinction/backscatter functional relation, Eq. (9), in case of Klett’s method, and assumed spatial dependence of Rayleigh scattering on molecules in Fernald’s case. Another issue is the extraordinarily difficult measurement of far-side extinction rate $`\alpha _\text{f}`$, needed in Eq. (10), and phase fraction $`F`$, Eq. (13). We are therefore forced to find better solutions, even at the expense of adding scanning capabilities to an otherwise rigid lidar setup.

In contrast to that, based on the sole assumption of a horizontally invariant (or at least horizontally slowly varying) atmosphere, the two-, and especially the multi-angle, method presented in section 6, while simple in structure, nevertheless produce reliable quantitative answers with small uncertainties (e.g., see Figs. 8 and 13) to FD calibration questions. As found by our investigation of first-run measurements, the relative error of OD for distances up to 12 km stay well below 6%. This number can be reduced even further by slow angular scanning and fast multiple-shot averaging of lidar returns. Nevertheless, in that case increased interaction between FD and lidar laser source, especially FD blind time, has to be taken into account. Furthermore, concerning the specific form of the atmospheric transmission entering Eq. (2), they offer suitable starting ground for development of methods that can considerably reduce systematic errors of shower energy $`E_{\text{em}}`$ estimation with fluorescence detectors.

In the case of strict horizontal invariance, both methods deliver exact solutions of the lidar Eq. (4) with accuracy of the results limited only by the quality of the measurement. In that way, they offer reliable framework for study of the notorious lidar ratio problem (i.e., extinction to backscatter codependency), widely discussed in the pure lidar community . Since, for example in the case of the Pierre Auger Observatory, where optical properties have to be known over large volumes of atmosphere, and a scanning lidar is therefore a necessity, both mentioned methods represent natural first choice of data analysis.

Authors would like to express gratitude to O. Ullaland for the support and encouragement during our work. Authors also wish to thank G. Navarra for assistance with EAS-TOP telescopes. This work has been supported by the Slovenian Ministry of Education, Science, and Sport under program No. P0-0501-1540.


format: text

text:
# Radiation Induced Landau-Lifshitz-Gilbert Damping in Ferromagnets

## Abstract

The Landau-Lifshitz-Gilbert damping coefficient employed in the analysis of spin wave ferromagnetic resonance is related to the electrical conductivity of the sample. The changing magnetization (with time) radiates electromagnetic fields. The electromagnetic energy is then absorbed by the sample and the resulting heating effect describes magnetic dissipative damping. The ferromagnetic resonance relaxation rate theoretically depends on the geometry (shape and size) of the sample as well as temperature in agreement with experiment.

preprint: QED Ferromagnetism

The description of ferromagnetic resonance damping in single domain samples is conventionally described by a transport coefficient $`𝖱`$ in the Landau-Lifshitz-Gilbert dynamical equation of motion. In detail, if $`𝒰[𝐌,S]`$ denotes the energy of the domain as a functional of magnetization and entropy, then the equation of motion is given by

$$\frac{𝐌}{t}=\gamma (𝐌\times 𝐇_{eff}),$$
(1)
where the effective magnetic intensity

$$𝐇_{eff}=\left(𝐇\frac{\delta 𝒰}{\delta 𝐌}\right)+\tau \left(\frac{𝐌}{t}\right)$$
(2)
and $`\gamma =(eg/2mc)=(1+\kappa )(e/mc)`$ is the electronic gyromagnetic ratio. Experimental data are usually expressed in terms of the dimensionless tensor transport coefficient $`\alpha `$; i.e.

$$\alpha =|\gamma 𝐌|\tau .$$
(3)
The transport (tensor) coefficient $`\tau `$ enters into the Landau-Lifshitz-Gilbert equation simply as an experimental parameter. In spite of considerable effort, presently there is no generally accepted theoretical picture for the physical source of the irreversibility in ferromagnetic resonance. In simple magnon decay models, the ferromagnetic relaxation takes place in two stages. A magnon with virtually zero wave vector decays into one or more higher wave vector magnons. The resulting burst of magnetic energy is later distributed to phonon lattice vibrations. Phonon modulations of the dipole-dipole and spin-orbit magnetic anisotropy were among the earliest of the relaxation mechanisms. However, phonon mechanisms in clean crystals yielded magnetization decay times which were much too long when compared with experiment. The situation was partially remedied by the notion that the magnons could decay via the final density of states as determined by lattice imperfections. But the only somewhat shorter imperfect lattice relaxation times no longer had a convincingly correct experimental temperature dependence. Also, the overall relaxation always remained longer than experiment. Furthermore, none of the above efforts in understanding $`\tau `$ take note of the experimental frequency dependence of damping widths in ferromagnetic resonance. There exists (in fact) an “impedance” $`𝖹(\zeta )`$ as an analytic function of complex frequency (for $`\mathrm{}m\zeta >0`$) whose real part

$$\tau (\omega )=\mathrm{}e\left\{𝖹(\omega +i0^+)\right\}$$
(4)
describes dissipation.

The purpose of this work is to provide a simple formula for the Landau-Lifshitz-Gilbert impedance $`𝖹(\zeta )`$ in terms of (frequency dependent) electrical conductivity $`\sigma (\zeta )`$ of the ferromagnetic sample. The final result is that

$$𝖹(\zeta )=\left(\frac{\sigma (\zeta )}{c^2}\right)\mathrm{𝟣}|𝐋|^2\mathrm{𝐋𝐋},$$
(5)
where the spatial average over the sample volume $`\mathrm{\Omega }`$ is defined as $`<\mathrm{}>=\mathrm{\Omega }^1_\mathrm{\Omega }(\mathrm{})d^3𝐫`$ and the length scale $`𝐋(𝐫)`$ is defined as

$$𝐋(𝐫)=\mathrm{𝐠𝐫𝐚𝐝}_\mathrm{\Omega }\frac{d^3𝐫^{}}{|𝐑|}\mathrm{where}𝐑=𝐫𝐫^{}.$$
(6)
The derivation of our central Eqs.(5) and (6) will be given in what follows.

The physical basis of our theory of the Landau-Lifshitz-Gilbert impedance is as follows: (i) When a spin wave decays into excitations which heat the sample, the current $`𝐉=ccurl𝐌`$ radiates an electromagnetic field. (ii) The power per unit volume absorbed from the radiation by the sample via (say) simple frequency independent Ohmic heating would then be

$$P=\sigma |𝐄|^2.$$
(7)
(iii) The eddy currents (produced by the electric field $`𝐄`$ via conductivity) describe the heating mechanism for ferromagnetic resonance relaxation no matter what the excitation products. The specific excitations (for example phonon excitations with lattice impurities) are all made manifest via the conductivity $`\sigma `$.

In detail, consider the magnetic field $`𝐁`$ produced by the magnetization $`𝐌`$ in the magnetostatic limit; i.e. with

$$𝐁=curl_\mathrm{\Omega }\left(\frac{𝐌\times 𝐑}{|𝐑|^3}\right)d^3𝐫^{}.$$
(8)
If one now applies Faraday’s law,

$$ccurl𝐄=\dot{𝐁},$$
(9)
to situations in which the magnetization varies slowly in time, then one finds a radiated electric field given by Eqs.(8) and (9) as

$$c𝐄=_\mathrm{\Omega }\left(\frac{𝐑\times \dot{𝐌}}{|𝐑|^3}\right)d^3𝐫^{}.$$
(10)
If the magnetization is uniform in space within the ferromagnetic sample, then Eqs.(6) and (10) imply

$$c𝐄(𝐫,t)=𝐋(𝐫)\times \dot{𝐌}(t).$$
(11)
We note in passing that the internal electric field Eq.(11) is closely connected to the demagnetization field intensity equation

$$𝐇_d(𝐫,t)=4\pi 𝖭(𝐫)𝐌(t)$$
(12)
for which the tensor

$$4\pi 𝖭=\frac{d^3𝐫^{}}{|𝐑|}=𝐋.$$
(13)
It follows that $`Tr𝖭=div𝐋/(4\pi )=1`$. For a uniformly magnetized sample, Eq.(11) plays the same role for the internal electric field that Eq.(12) plays for the internal magnetic intensity. The power dissipation per unit volume in Eq.(7) may be evaluated using Eq.(11) yielding

$$P=\frac{1}{c^2}\sigma \left|𝐋\times \dot{𝐌}\right|^2=\dot{𝐌}\tau \dot{𝐌}$$
(14)
where $`\tau `$ is the Landau-Lifshitz-Gilbert transport coefficient tensor entering into Eqs.(1) and (2). Eq.(14) implies that

$$\tau =\frac{\sigma }{c^2}\mathrm{𝟣}|𝐋|^2\mathrm{𝐋𝐋}.$$
(15)
The above Eq.(15) holds for the case in which $`\sigma `$ and $`\tau `$ do not depend on frequency. If the conductivity depends on frequency $`\sigma (\zeta )`$, then Eq.(15) is easily generalized to Eqs.(4) and (5) which are the central results of this work. Note that the tensor nature of $`𝖹(\zeta )`$ implies that ferromagnetic damping depends on the shape and size of the sample. This theoretical consequence of our theory is of experimental importance.

For the case of a thin film, the vector $`𝐋`$ is in the direction of the normal unit vector $`𝐧`$ to the film; i.e. $`𝐋=4\pi (𝐧𝐫)𝐧`$. The Landau-Lifshitz-Gilbert transport coefficient tensor for a thin film of thickness $`d`$ is then

$$\tau =\left(\frac{4\pi ^2d^2\sigma }{3c^2}\right)(\mathrm{𝟣}\mathrm{𝐧𝐧})(\mathrm{for}\mathrm{thin}\mathrm{fims}).$$
(16)
In most experiments on ferromagnetic films $`𝐧𝐌=0`$. In the plane of the film, $`\tau `$ can then be described by an isotropic scalar $`\alpha `$; i.e. $`\alpha \mathrm{𝟣}_{}=|\gamma 𝐌|\tau _{}`$ where

$$\alpha =\left(\frac{4\pi ^2d^2\sigma }{3c^2}\right)|\gamma 𝐌|.$$
(17)
In terms vacuum impedance $`R_{vac}=(4\pi /c)`$ and the film resistance “per square” $`R=1/(\sigma d)`$, one obtains the simple result

$$\alpha =\frac{\pi }{3}\left(\frac{R_{vac}}{R}\right)\left(\frac{|\gamma 𝐌|d}{c}\right).$$
(18)
Since $`R`$ increases with increasing temperature, one expects that $`\alpha `$ should fall with increasing temperature. This is (in fact) an observed result.

We have presented above, a physical picture of how the Landau-Lifshitz-Gilbert damping of ferromagnetic resonance takes place. Just as a static magnetization $`𝐌`$ produces a demagnetizing magnetic field intensity $`𝐇_d`$, a changing magnetization $`\dot{𝐌}`$ produces an electric field $`𝐄`$. The resulting Ohms law conduction current $`𝐣=\sigma 𝐄`$ produces the Landau-Lifshitz-Gilbert damping magnetic field intensity $`𝐡=\tau \dot{𝐌}`$. Thus the damping parameter $`\alpha `$ is directly and simply related to the conductivity $`\sigma `$ of the sample. For thin films, this relationship is given in Eq.(18). Previously puzzling temperature dependences of $`\alpha `$ appear now as self evident.


text:
# Null-induced mode changes in PSR B0809+74

## 1 Introduction

Lately, it has been quiet around bright PSR B0809+74, once the canonical example of regular subpulse drifting and the centre of lively discussions.

In the thirty years since their discovery by Drake & Craft (1968), drifting subpulses feature in many discussions about the nature of pulsars and the pulsar’s emission mechanism. The subpulse-drifting phenomenon in itself is simple: when comparing adjacent individual pulses, the subpulses that comprise them are seen to shift regularly through the average pulse window. The left panel of Fig. 2, for example, shows a recent observation in which these subpulses form their driftbands.

In 1968, the interpretation of Drake & Craft that the subpulse drifting is the pulsation of the neutron star, is ground zero for much debate. Over the following years, the ever increasing quality and volume of observations of subpulses in different pulsars continue to limit possible models and require their refinement. It is in this process that observations of 0809+74 take the lead.

Alexeev et al. (1969) and Vitkevich & Shitov (1970) are the first to detect the driftbands of 0809+74. Cole (1970) then notes that the driftrate of the subpulses occasionally changes, without being able yet to identify the nulls as the trigger. One year later, Taylor & Huguenin (1971) do find that this occasional cessation of the pulsar’s emission precedes a changing driftrate, lasting a few driftbands.

In his paper solely devoted to the drifting in 0809+74, Page (1973) finds that the subpulse separation changes through the profile. He also introduces the ‘subpulse phase’, the difference between the actual position of the subpulse and its predicted position. While this makes it easier to plot long series of subpulse positions horizontally, it also introduces the notion that the subpulse position dramatically jumps in phase over a null. From the data however, it is clear that the subpulse position over the null does not change much at all – on the contrary, it is the preservation of phase that is the astonishing phenomenon.

Unwin et al. (1978) are the first to note this. They find that the position of the subpulses is identical before and after the null and conclude that, even though there is no emission from the pulsar, the emitting structures survive the null. Five years later, Lyne & Ashworth (1983) use a method devised by Ritchings (1976) to find the nulls of 0809+74 and investigate the time scales associated with the change in drift pattern around a null. The first time scales they look into are those of the decay and rise of emission around a null. A finite decay and rise time would cause the pulses neighbouring a null to be less bright than the average. A relative dimness of the last pulse before the null is found indeed, pointing at a decay time shorter than 5% of the pulse period. Contrary to the expectations, however, the first pulses after the nulls are found to outshine the average ones.

The next time scales considered are those of the slowdown and speedup of the subpulse drift around the null. Lyne & Ashworth notice that, contrary to the findings of Unwin et al., there is a some change in subpulse positions over the null. They suggest that this change originates in subpulse drift that is caused by a gradual speedup of the drifting during the null. Being several tens of pulse periods, this speedup time scale is much larger than the time scales found for the emission decay and rise, and the subpulse drift slowdown.

In 1984, simultaneous observations at 102 and 1412 MHz by Davies et al. confirm the curvature of the driftbands found by Page (1973) and show that the subpulse width varies across the pulse profile.

Although work on 0809+74 quiets down after this intense period, significant observational progress for the understanding of subpulse drifting in general is made on 0943+10. Deshpande & Rankin (1999, 2001) are able to detect periodicities in the subpulse strengths that lead to a determination of the recurrence time of individual subpulses.

On the theoretical front, the initial hypothesis that the subpulse drifting originates in pulsations of the neutron star (Drake & Craft 1968) is soon abandoned. The subpulses are now thought to be formed in the magnetosphere of the pulsar, and the topic of debate changes to the exact location of their formation. The observed periodicity is first thought to be caused by quasi-periodic fluctuations in the plasma near the co-rotation radius of the pulsar. The polarisation and beaming of the subpulses, however, later tie the emission mechanism to the surface of the pulsar.

In 1975, Ruderman & Sutherland build on the Goldreich & Julian (1969) and Sturrock (1971) models to propose a theory that incorporates the drifting subpulses and the nulling. They suggest that the pulsar emission is formed in discrete locations around the magnetic pole. These locations would rotate around the pole like a carousel. In one of many revisions of the Ruderman & Sutherland model, Filippenko & Radhakrishnan (1982) propose a model that explains the survival of the subbeam structure over the null. The work of Deshpande & Rankin (1999, 2001) visualises these models by mapping the observed subpulses onto their original positions on the carousel.

In this paper we investigate the behaviour of the pulsar emission around nulls. We will use observations of two extraordinary long sequences of post-null behaviour to explain the features of general subpulse drifting around nulls in a new framework.

## 2 Observations

### 2.1 Data

Over a period of about 18 months we have observed 0809+74 with PuMa, the Pulsar Machine, at the Westerbork Synthesis Radio Telescope (WSRT) in the Netherlands. Most observations lasted several times the scintillation time scale of this pulsar. We have only used the data in which the pulsar was bright, amounting to a total of 13 hours. With a pulse period of 1.29 seconds, this comes down to $`3.6\times 10^4`$ pulses.

The collecting area of WSRT and the low system noise made it possible for PuMa to record the data with a high timing resolution and a high signal-to-noise ratio (SNR). The WSRT consists of fourteen 25-meter dishes arranged in east-west direction. For our observations, signals from all the fourteen dishes were added after compensating for the relative geometrical delay between them, so that the array could be used like a single dish with an equivalent diameter of 94 meters. The pulsar was observed with PuMa (Kouwenhoven 2000; Voûte 2001; Voûte et al. 2002) at centre frequencies between 328 MHz and 382 MHz, with bandwidths of 10 MHz. Each 10-MHz band was split into 64 frequency channels, and the data was recorded, after being digitised to represent each sample by four bits (sixteen levels). After inspection for possible electrical interference, we corrected for the interstellar dispersion during our off-line analysis.

### 2.2 Reduction

#### 2.2.1 Finding nulls

Normally, the energy of individual pulses varies from about 0.5 to 2 times the average. When the pulsar is in the null state however, which it is for about 1.4% of the time, it is does not emit. In this case, very little energy is found where the pulse is expected. We have used this difference in on-pulse energies to identify the nulls. To do this, we compared the on-pulse energies with the energy found in off-pulse region, where no emission is expected.

To calculate the on-pulse and off-pulse energies, we split each pulse period $`P_1`$ (360o in longitude) into three sections. The first section, $`0.15P_1`$ wide, was centred around the peak of the average profile. The second section, equally wide, consisted of a part of the off-pulse region. A third section, $`0.5P_1`$ wide, contained a different, independent part of the off-pulse data and was used to estimate the noise and baseline of the signal. This baseline was subtracted from both the on-pulse and off-pulse data.

To calculate the energy in each pulse, we determined the central range in longitude that on average contained 90% of the power the pulsar emitted in the on-pulse section. For this region we summed the amplitudes of the signal to get the energy content. We did the same for an equally large piece of the off-pulse section. To correct for long term variations in the pulse energy due to scintillation, we subsequently scaled both energies with the average energy of the 100 surrounding pulses.

We defined that, in the null state, there is no observable radiation from the pulsar. This means that noise is all one observes at the time a new pulse is expected. The energy distribution of the null-state is then identical to that of the off-pulse region (high distribution in the background of Fig. 1). For our data this meant that all pulses with energies less than the highest energy found in the off-pulse distribution were nulls (small peak around zero energy in Fig. 1). Furthermore, the off-pulse distribution gives a better estimate of the expected spread in null energies than the null distribution itself, due to the larger number of pulses included.

This method, originally devised by Ritchings (1976), was first used on 0809+74 by Lyne & Ashworth (1983). In their case, there was considerable overlap between the on-pulse and off-pulse distributions, hindering their identification of the nulls. We wanted to be certain that the set of nulls we found was genuine and complete; therefore we only used the 60% of our observations in which the on-pulse and off-pulse distribution were separated by more than $`0.1E/<E>`$.

With a 3-sigma point at 15 mJy, the low noise makes it possible to fully distinguish between pulses in the normal state on the one side, and those in the null state on the other, giving unprecedented insight into the nature of the null state.

#### 2.2.2 Fitting subpulses

The major step in our data reduction was to extract the four basic elements of the individual pulses: the number of subpulses and their respective positions, widths and heights. By doing so we greatly increased the speed of subsequent operations, but at the same time kept a handle on the physics by retaining the parameters most important for visualising the data.

The high SNR of our observations showed a wealth of microstructure in the individual subpulses, which on average were Gaussian in shape. Plainly fitting Gaussians to the intricate, many featured subpulses did not immediately result in locating them reliably. Due to the changing amplitude of the signal throughout the window, the microstructure of the brightest subpulses was often more important to the goodness-of-fit than complete, weaker subpulses near the edge of the pulse window. In order to still find these weak subpulses, we devised a two-step approach of first finding the strong subpulses and then using their driftpath to suggest the position of weaker subpulses. This approach worked very well.

For our fitting, we assumed the subpulses to be Gaussian in shape, non overlapping, and to have a full width at half maximum (FWHM) less than 15o. We used a Levenberg-Marquadt method (Press et al. 1992) with multiple starting configurations to produce goodness-of-fit values for different numbers of subpulses. By comparing these $`\chi ^2`$ values, we located the subpulses that had high significance levels.

Upon finding these normal to strong subpulses, we set out to detect the weaker ones over the noise. For this purpose we composed driftbands out of the individual subpulses already identified. On a one-by-one basis, a subpulse was either added to a path that had predicted its position to within $`P_2/4`$ ($`P_2`$ being the average separation between two subpulses within a single pulse, about $`11^o`$ longitude), or taken to be the start of a new path. Paths ended if no subpulse fitted the path for more than 5 pulses, or upon reaching a null. On the first run, only paths longer than ten pulses were allowed to survive, so as to eliminate the interference of short run-away paths. This allowed a steady pattern of long paths to grow, incorporating about 90% of the subpulses found. On the second run, paths longer than 5 pulses were allowed to form: these shorter paths occur only around nulls, where the normal, long paths are interrupted.

The drifting pattern thus formed predicted the locations of the weaker subpulses. Around each predicted position, we isolated a section of the data. This section was taken as large as possible without interfering with other, previously fitted subpulses. We checked the significance of fitting a single Gaussian to indicate the presence of a subpulse.

In the same way as described above, the final driftpattern was then identified using the extended set of subpulses.

To check the method, we compared the observed average profile with the one recreated from the Gaussians fits to the subpulses. From their match, we concluded that the fitting procedure was effective.

## 3 Results

### 3.1 Driftband fitting

The driftbands are not straight, but show a systematic curvature (Davies et al. 1984). This is seen most strikingly when we plot the residuals to straight line fits: Figure 3 shows the longitudinal offset between the actual and predicted positions of the subpulses. Near the edges of the profile the subpulses arrive later than is to be expected based on a linear driftband, in the middle they arrive sooner.

This probably explains the results of Popov & Smirnova (1982). After fitting straight lines to all driftbands, they find that the driftrate of the last driftband before the null is 20% higher than that of a normal driftband. As seen in the curvature of the driftbands, however, a normal driftband consists of a fast drifting first half and a slow drifting second half. If a driftband is cut by a null, the part before the null consists only of the fast drifting first half, resulting in a higher average driftrate over this shorter driftband.

To make our results independent of the longitude, the non-linearity of the driftbands is taken into account in all the following driftrate calculations. In these cases, the positions of individual subpulses are corrected by subtracting the appropriate residual value.

### 3.2 Nulling

#### 3.2.1 Shortest nulls

The main criterion used to decide whether or not to include certain datasets was the complete separation of the on-pulse and off-pulse energy distributions. Therefore, the set of null identified is genuine and complete. This allows us to investigate the underlying statistics of null occurrence and length, and estimate the influence of short nulls.

The chance of finding 0809+74 in the null state is on average about 1.4%. If the null lengths were distributed in a Poissonian fashion, we would find a preponderance of one-pulse nulls and very few longer ones. However, in the null-length histogram (Fig. 4) we find a peak at nulls of length 2, a significant decrease towards shorter lengths and a considerable number of long nulls, showing that the occurrence of nulls is not governed by pure chance.

Comparing this histogram with the one previously found by Lyne & Ashworth (1983), we see that we identify about twice as many one-pulse nulls in the data. Still, many nulls shorter than one period must pass unnoticed, as they occur when the pulsar faces away from us.

We know that nulls, especially long ones, have a distinct impact on the drifting pattern in their vicinity. Although the impact of shorter nulls is less, a large number of unnoticeable, short nulls ($`<P_1`$) might seriously influence the drifting pattern we are trying to understand.

Using the null-length histogram we can estimate the number of these short nulls. It peaks at two-pulse nulls, and the distribution decreases towards shorter lengths. Assuming that the underlying distribution of null-lengths is continuous, this tendency of decreasing occurrence towards shorter nulls implies that there is a small number of nulls shorter than one pulse period. Extrapolating the decrease leads to an estimate of about 15 unnoticeable short nulls in the null-length interval from 0 to 0.5 pulses.

The low frequency of their occurrence (0.04%) indicates the influence of short nulls on the drifting pattern is negligible.

#### 3.2.2 Null versus burst length

The next question we address involves the interval between adjacent nulls (the so-called burst) and the duration of the nulls. Does waiting longer for a null mean it will last longer, too? We have checked these relationships and have found that the lengths of neighbouring nulls and bursts are independent.

#### 3.2.3 Position jump over nulls

If the nulling mechanism is independent of the position of the subpulses, we expect that the distribution of subpulse positions is the same for the normal pulses and the pulses that immediately precede a null. We find no proof of differences in these distributions and conclude that there is no preference for a null to start at a certain subpulse position.

The positions of subpulses change over a null. We derived this shift of the subpulse pattern for each null in our sample. As the positions of the individual subpulses were already identified, this shift was simply extracted. Each subpulse before the null matched a subpulse after the null, if the latter fell within $`3P_2/4`$ to $`P_2/4`$ of the former. This range is symmetric around the average expected jump over a null, so as to minimise the number of ambiguous cases. In the following analysis, we used the average of all individual subpulse shifts within one pulse.

Figure 5 shows this shift in the subpulse position over a null, corrected for the non-linear behaviour of the driftband. If the motion of the subpulses were independent of the emission, the subpulses would continue to drift invisibly throughout the null state, and reappear at a very different position. The associated shift in positions would then be spread around the diagonal dashed line in Fig. 5.

If, on the other hand, the subpulse drift would cease abruptly and completely during the null, the jump in position would be between zero and two times the average shift in longitude between normal pulses. The uncertainty in this estimate arises from the fact that neither the start nor end of the null are known more accurately than to one pulse period. The mean jump in position over the null would then follow the horizontal dash-dotted line in Fig. 5.

After discarding the three ambiguous cases (points near the top edge of Fig. 5), we computed averages for each null length. These averages follow neither of the two cases outlined above. There is too much change in position over the null to be accounted for by just an abrupt stop of the drift, and there is no evidence for steady (albeit slower than regular) drift during the null.

We do see that for nulls longer than one, the jump is constant, $`1.47\pm 0.16^o`$ above the offset value that we would expect in the case of no drift. This independence of null length and subpulse jump over the null is shown as the dotted horizontal line in Fig. 5.

#### 3.2.4 Driftrate around the null

We have computed the driftband slope around nulls, correcting for their general curvature.

For the driftrate before the null, we fitted straight lines to the last six subpulses of each driftband the ended with the null. We find that this driftrate just before a null does not deviate from the normal driftrate.

The driftrate after nulls is different from the normal driftrate, though. Although there is some spread, all the driftrates we find after longer nulls are lower than the normal average value (Fig. 6).

#### 3.2.5 Average pulse profile around null

To investigate whether the change from normal emission to the null state is sudden or more gradual, Lyne & Ashworth (1983) compared the energies of the pulses near a null. The finite chance that the emission from the pulsars drops or rises within the pulse window would influence the brightness of the pulses around the null. The last pulse before the null was indeed found to be less bright than a normal pulse. The first pulse after the null, however, was considerably more bright than a normal pulse.

We have not just looked at the energies of these neighbouring pulses, but also at their profiles. To this end, we have averaged all the last pulses before nulls longer than 1 period, and all the first pulses after these nulls. The results, shown in Fig. 7 and condensed in Table 1, are surprising. Not only do we find the expected differences in brightness, we also see a significant offset in the pulse position for the first pulses after nulls.

### 3.3 Slow drifting mode

Although we had set out to quantify the normally very regular drifting behaviour of 0809+74, we unexpectedly found two occasions in which the pulsar clearly deviates from its normal drifting mode. We will refer to these sequences by the year in which they were observed, 1999 and 2000. In Fig. 8a we have plotted the derived longitude (Lyne & Ashworth 1983) of subsequent pulses around the mode changes. The derived longitude of the subpulses effectively converts the different short driftbands into one long one. Figure 9 shows a grayscale plot of the two slow drifting series.

The most striking difference, as the drifting pattern is concerned, is clearly the decrease of the driftrate by about 50%.

We see that in both the observations, the drifting mode changes during or immediately following a long null. After these nulls of length 10 and 5 respectively, almost all drifting properties reappear at new values, as laid out in Table 2. As some of these properties normally already show change over time, we compare the slow drifting sequences (labeled ‘slow’) to the 600 pulses that surround them (labeled ‘normal’).

In the first three columns of Table 2, we investigate the parameters of the individual subpulses in both modes. Immediately we see a very interesting change in the relative average position of the subpulses. Right after the null, there is a definite shift in position towards earlier arrival, that can already be seen in a plot like Fig. 9. Interestingly, this shift of the pulse window is similar to the average jump in subpulse position we see over a normal null.

In the slow drifting mode, the subpulses are slightly wider, but their heights remain the same.

Next, we explore the driftband characteristics in the last three columns of Table 2. The average longitude separation of two adjacent driftbands, $`P_2`$, decreases significantly by about 15%. Hence, in the slow drifting mode the subpulses within each pulse are spaced closer together than normal.

When we compare the driftrates of the subpulses between the two modes, we find that the driftrate in the slow drifting mode is almost halved. This puts the new driftrate right in the range of driftrates we usually find after a longer null. We have indicated these driftrate values with black stars in Fig. 6.

Being dependent on $`P_2`$ and the driftrate, the fractional change in $`P_3`$ (the recurrence time of a driftband) is comparably large.

In both cases the new drifting mode is stable for about 120 pulses and then changes back to normal. To illustrate this change, we have plotted the deviation of the derived longitude from normal drifting in Fig. 8b. For the 1999 observation, we see a normal driftrate up to the null at pulse 200. After this null, the driftrate is smaller, up to pulse 330. In about 20 pulses, the drift then speeds up back to its normal value.

For reference, we have also plotted the drifting behaviour after the longest null in our sample. Again, we see a speedup to the normal driftrate in about 20 pulses at pulse 330, very similar to the slow drifting speedup time scale.

The slow drifting interval in the 2000 observation (right panel of Fig. 9) is followed and stopped by a series of frequent, longer than average nulls, out of which the pulsar emerges in its normal drifting pattern.

These surprising changes in subpulse positions, widths and separations must influence the average pulse profile during the slow drifting mode. Noting the many similarities between slow drifting pulses and the pulses that follow nulls (halved driftrate, similar offsets, similar speedup) we are interested in the possible similarities between the average profiles of these pulses, especially since the average profile of the pulses after the null is singularly bright and shifted in longitude.

In Fig. 10 and Table 1 we compare the pulsar’s normal average profile with the average profile in the slow drifting mode. In both observations, the slow drifting mode average profile is much brighter than the normal profile and offset towards earlier arrival – exactly the two peculiarities we found in the average pulse profile after a null as well.

## 4 Discussion

### 4.1 Driftband fitting

The non-linearity of the driftbands was already noted by Page (1973), who integrated several hundreds of pulses to compare the shape of the driftbands in different observations. He found considerable differences between these driftband shapes. We compare the average drift path over several thousands of pulses however, and find that the curvature of the driftbands remains the same.

The curvature of the driftband can be a direct consequence of the curved geometry of the emission region (Krishnamohan 1980). In that case, the shape of the driftband is expected to be point-symmetric: the curve will then be unchanged after mirroring it in both axes. This point-symmetry would then return in the residuals to straight line fits (Fig. 3). The shape of the driftbands we find, however, is not point-symmetric but axisymmetric: the curve is unchanged after mirroring in the y-axis.

We could still tie the curvature of the driftbands to the emission-region geometry, by assuming that part of the pulse profile of this pulsar is missing. In that case, the driftband shape we find would represent only part of the total expected shape.

This suggestion that part of the profile of 0809+74 is missing fits in with observations of this pulsar’s pulse shape at different frequencies (Bartel et al. 1981; Kuzmin et al. 1998). These imply that at lower frequencies the profile is partially “absorbed”.

### 4.2 Shortest nulls

We find that the number of invisible nulls is low. That means we can study of the drifting properties of 0809+74 straightforwardly.

### 4.3 Slow drifting mode and nulling

With our investigation of 0809+74 still ongoing, we discuss here only the phenomena themselves and defer their interpretation to a subsequent paper.

The similarity between this pulsar’s behaviour in the slow drifting mode and around nulls is striking:

* The slow drifting sequences start at nulls.
* The driftrate of the slow drifting mode is a lower limit to the driftrates found after all nulls.
* The speedup from slow drifting to normal is identical to the speedup after a null.
* The average profiles of both the post-null and slow drifting pulses are brighter than normal.
* These average profiles are both displaced to earlier arrival.
* This displacement of the average profile is caused by offsets in the individual subpulses.
* These offsets are identical to the jump of the subpulses over the null.

The natural conclusion is that the behaviour after a null and the slow drifting mode are the same, quasi-stable phenomenon. Normally the pulsar reappears from the null in the slow drifting mode. After a variable time, it quickly evolves to the normal mode. Therefore, right after a normal null we see either a short sequence of slow drifting, or the transition back to normal drifting. In the case of the long slow drifting mode sequences, the metamorphosis back to the normal mode is delayed.

This would explain all the similarities we found above. After a normal null, the pulsar is in the slow drifting mode or changing back to the normal mode. All the characteristics of the slow drifting mode can then be identified in the post-null behaviour, although they will be less pronounced; the transition to normal drifting may already be taking place.

When the driftrate increases at the end of a slow drifting sequence, this speedup is quick and identical to the speedup seen after a normal null (see Fig. 8b).

The post-null driftrate values are found in between the slow drifting and normal driftrate value, depending on how soon the transition back to normal takes place (see Fig. 6). Although the return from the low to the normal driftrate will be quick, the time the driftrate is low may vary for different nulls of the same length. The average of these different sequences will then resemble the slow exponential decay found by Lyne & Ashworth (1983).

Comparing average pulse profiles, the increased brightness and the pulse offset of the post-null average profile are attenuated versions of similar deviations seen in the slow drifting mode profile. (see Table 1, Figs. 7 and 10)

The average profile offset we see in the slow drifting mode is caused by a shift of the window in which the subpulses appear (Tables 2 and 1).

The change in position of the post-null average pulse profile must then be caused by this shift of the pulse window as well. The magnitude of this shift is identical to the subpulse-longitude jump over normal nulls.

This means that the subpulse-position jump over a null is caused by a displacement of the pulse window as a whole, like in the slow drifting mode.

Previously, this jump was thought to be the effect of the subpulse-drift speedup during the null. For this speedup to produce the observed jump in subpulse position, the time scale involved had to be long, contrasting the short time scales found for the slowdown of the subpulse drift and the rise and fall of the emission around a null.

With the displacement of the subpulses over the null accounted for, the estimated speedup time of the subpulse drifting is negligible: the preservation of the position of the subpulses over the null now only allows for a quick speedup of the subpulse drift, putting all time scales of emission and drift rise and decay around the nulls in the same range.

## 5 Conclusions

After many or all nulls, 0809+74 emits in a mode different from the normal one. This mode is quasi-stable, normally changing back to normal in several pulses. This is seen as the normal behaviour after a null, of which the reduced driftrate is the most striking characteristic. Occasionally, the quasi-stable slow drifting configuration persists for over a hundred pulses before changing back to the normal mode.

The pulses in the slow drifting mode and consequently all pulses after the null are brighter than normal pulses. In the slow drifting mode, the subpulses are closer together, they drift more slowly through the profile and the window in which they appear is offset towards earlier arrival.

This offset of the pulse window accounts for the displacement of subpulses over the null. When taking this shift of the window into account, we find that the longitude of the subpulses is perfectly conserved over a null. This indicates that the speedup time for the subpulse drift is short.


text:
# Toward a Possible Solution to the Cosmic Coincidence Problem

## Introduction

The discovery that the Universe is accelerating , probably due to dominance of vacuum energy, has pushed two uncomfortable fine-tuning problems to the forefront of the physics community. The first problem is the that scale of this vacuum energy density, (.002 eV)4 is vastly different from the MeV to Planck scale particle physics vacuum energies that must sum to give this number. In the past, one hoped that some unknown symmetry principle would require these energies to ultimately sum to zero, but the existence of a non-zero sum naively dashes these hopes. The second fine-tuning problem occurs because traditional vacuum energy density does not change with time while matter and radiation density change rapidly as the Universe expands. At early times the Universe was radiation and then matter dominated, and only recently, at $`z1.7`$, became vacuum dominated. The importance of this “cosmic coincidence” is that in the normal scenario, once vacuum energy becomes dominate, it stays dominate, and if this had occurred at any earlier epoch, the evolution of the Universe would have been completely different and most likely we wouldn’t be here to discuss it. Currently, there is no convincing fundamental physics idea for why vacuum dominance happened only recently, and this has led many workers  to conclude that some sort of anthropic principle must be at work. The anthropic idea is that there is an ensemble of universes with different values of the vacuum energy, most of which do not allow life to develop. Therefore the cosmic coincidence is “explained” by saying the existence of intelligent life selects only those values of vacuum energy density near that which has been measured in the SN Ia observations. Many non-anthropic ideas have been proposed, which reduce the coincidence, usually by having the vacuum energy “track” the matter density in some way so that the ratio is not so large, but all of these have been criticized as in fact involving fine-tuning in some way . At this point, some workers have concluded that only anthropic ideas have a chance of explaining the cosmic coincidence .

The goal of this paper is to propose a new class of solutions to the cosmic coincidence problem which is not anthropic and which may allow removal of fine-tuning. The idea presented here is similar to, but somewhat more general than, that of the tracking oscillating potential model of Dodelson, Kaplinghat, and Stewart . This class of models also makes it possible to change the current cosmological constant problem from why the vacuum energies sum to a small number, back to the older problem of why the vacuum energies sum to zero. We do not have anything to say about the basic cosmological constant problem of why this sum should be zero in first place. To understand our class of potential solutions an analogy may be helpful.

## An Analogy

Suppose we demand a non-anthropic solution to the cosmic coincidence problem. What would such a solution look like? Physics offers several examples of potential cosmic coincidences that we do not consider as such. For example, U238 is a common radioactive substance with a half life of 4.5 Gyr, almost exactly the age of the solar system. The reason that this is not a cosmic coincidence is clear in this case. There are thousands of nuclides with half lives ranging over a enormous range of time scales, from microseconds (e.g. U222, $`\tau =1\mu `$s), to seconds (e.g. U226, $`\tau =0.5`$s), to days (e.g. U231, $`\tau =4.2`$d), to millennia (e.g. U233, $`\tau =1.6\times 10^5`$ yr) to the age of the Universe (e.g. U238), to $`10^{20}`$ years (e.g. Se82, $`\tau =1.4\times 10^{20}`$ yr), and up.

Thus there is no surprise that for any given time scale, such as the current age of the Solar System, there is some nuclide which is decaying on just this time scale. Note that if U238 were the only nuclide in existence and everything were made of it, then there would be a cosmic coincidence problem quite similar to the vacuum energy cosmic coincidence problem. Also, if the range of radioactive decay constants were not exponentially distributed over such a wide range of time scales then again even with hundreds of nuclides, finding one with $`\tau =4.6`$ Gyr would be unlikely. In this example, it is the exponential sensitivity of radioactive decay to the nuclear wavefunction that allows a small change, such as the addition of one neutron, to make a large change in the half life.

Thus we are inspired to seek a possible class of models to solve the vacuum energy cosmic coincidence by positing not just one scalar field whose vacuum energy is making its appearance today, but an ensemble of fields, whose vacuum energies span an exponentially large range of energy densities, some of which have dominated the Universe for short periods of times at many time scales in the past. We note that we have no a priori reason for the existence of such an ensemble of vacuum fields, but are using the cosmic coincidence problem as a clue that such an ensemble might exist.

In this class of models, the specific field (or fields) responsible for the current acceleration is not special, but just happens to be the one dominant at this time. We note that in this context, the purported inflaton field responsible for cosmic inflation in the very early Universe could be just another one (or more than one) of these fields.

This type of solution to the cosmic coincidence problem makes some interesting predictions:

* The vacuum energy is not a cosmological “constant”. The current phase of acceleration is temporary and will eventually finish; $`w=p/\rho `$ is not $`w=1`$ exactly.
* There probably were several periods of acceleration and vacuum dominance in the past, followed by radiation and/or matter domination, and finally vacuum dominance again today.
* There will likely be additional periods of acceleration in the future. Thus predicting the ultimate fate of the Universe will not be possible without an understanding of the origin of all these fields and their vacuum energies.
* The sum of all these changing vacuum energies may well eventually be zero; that is, the minimum of the potential of all these fields may be a zero that we are evolving toward. Just as even very long lived nuclides will eventually decay, it is thus possible that the actual cosmological “constant” is zero, and we are just part way there. We offer no suggestion here as to why the sum should be exactly zero, just note that in this scheme it is possible. However, just as some nuclides are stable, it is also possible that the final cosmological constant is not zero.

We note that the tracking oscillating model of Dodelson, Kaplinghat, and Stewart  is similar to our class of models. It uses a single scalar field with a potential with many wiggles and also predicts many periods of acceleration in the past. It addresses the cosmic coincidence problem in a way similar to ours. One difference is that the single potential must be quite complicated to allow such behavior, and it seems interesting and somewhat more generic and flexible to consider an ensemble of scalar fields with simple potentials.

## Constraints

There are several observations that constrain vacuum energy density in the Universe and which any model of an accelerating Universe must satisfy. Our model contains several fields and so differs somewhat from single field quintessence/dark energy models . However, given that we do not have a specific model, we will not attempt an accurate tallying of the constraints, but only mention a few constraints that should be considered.

First, if the density of scalar field energy, $`\rho _{vac}`$, is more than a few percent of the radiation density, $`\rho _{rad}`$, at the time of big bang nucleosynthesis (BBN), the universal expansion rate will increase enough to make a significant difference in the predicted abundances of helium and deuterium. Since these abundances are fairly well measured, we demand $`\rho _{vac}/\rho _{rad}<.02`$ when $`10^8<z<10^{10}`$.

Second, the equation of state $`w`$ must be sufficiently negative for the scalar field energy solutions to match the current type Ia supernova measurements and be consistent with the cosmic microwave background (CMB) and large scale structure measurements. We require $`w<0.79`$ (95% CL) . In addition we can directly calculate the absolute magnitude-redshift relation and compare to the supernovae reported in Perlmutter, et al. , demanding that, within errors, it be as good a fit as a cosmological constant model. Also, since the growth of large scale structure can be reduced by periods of vacuum dominance, or near vacuum dominance, there are several observational constraints on the power spectrum that are possible. For simplicity we will merely calculate the linear growth factor from the time of matter-radiation equality until the present, and compare it to that expected from a cosmological constant model.

Next, additional vacuum energy at the time of photon decoupling can shift the well measured acoustic peaks in the CMB, providing the constraint $`\rho _{vac}/\rho _{other}<0.64`$, at $`z1100`$  where $`\rho _{other}`$ includes dark matter, radiation, and baryons. To preserve the peak positions, we can also demand that the angular diameter distance to surface of last scattering at $`z1100`$ not be too different from that implied by a cosmological constant model.

We easily find models in which these constraints are satisfied, basically by chosing the field content and parameters such that none of the additional scalar fields is very important during either the decoupling or nucleosynthesis epochs, and such that the field that dominates today has $`w`$ in an appropriate range. From one point of view it thus requires fine tuning of many parameters to satisfy these constraints and so these restrictions are a weakness of our idea. However, from another point of view, these constraints and the data from which they derive are actually just measurements of the initial scalar field content and parameters. That is, had these values been different, then the present the values of our current cosmological parameters and the the positions of the CMB peaks would now be different.

## Example model

Since this class of solutions is motivated purely by solving the cosmic coincidence problem, there are no restrictions on the types of fields or forms of potentials that may be used. As a simple example consider an ensemble of $`N`$ scalar fields $`\varphi _i`$, $`i=1,\mathrm{},N`$ that do not interact with ordinary matter or each other, and which have potentials of the form

$$V_i(\varphi _i)=\lambda _i\varphi _{i}^{}{}_{}{}^{\alpha _i}.$$
(1)
The total scalar field potential is then just $`V=_iV_i(\varphi _i)`$. To first approximation (ignoring gravity, finite temperature effects, parametric resonance, etc.) each field is governed by the standard equations: $`\ddot{\varphi _i}+3H\dot{\varphi _i}+V_i^{}(\varphi _i)=0`$, $`H^2=\frac{1}{3M_{\mathrm{Pl}}^2}[\rho _{\mathrm{other}}+\frac{1}{2}_i\dot{\varphi _i}^2+_iV_i(\varphi _i)]`$, where $`H(z)\dot{a}/a`$, and the Planck mass is $`M_{\mathrm{Pl}}=(8\pi G_N)^{1/2}=2.44\times 10^{18}\mathrm{GeV}`$, and where $`a`$ is the scale factor of the Universe, $`\rho _{\mathrm{other}}`$ is the total energy density of the other contributing fields, the dot represents differentiation with respect to time, and the prime indicates differentiation with respect to $`\varphi _i`$. Note that with $`\alpha =4`$, the $`\lambda _i\varphi _{i}^{}{}_{}{}^{\alpha _i}`$ form for the potential is that used in chaotic inflation, and if $`\varphi _i`$ starts at a non-zero value, it will eventually approach zero under the influence of its potential as long as $`\alpha _i`$ is an even integer. The $`\varphi _i`$ field may come to dominate the energy density of the Universe depending upon the magnitude of its initial value and the speed at which it goes toward zero, which is determined by $`\lambda _i`$ and $`\alpha _i`$. Thus in this model, the period of vacuum dominance will not last forever, since eventually $`\varphi _i`$ reaches zero and then oscillates around zero. During the oscillation phase, $`\varphi _i`$ effectively behaves as material with $`\overline{w}_i=(\alpha _i2)/(\alpha _i+2)`$, and density scaling as $`\rho _ia^{3(1+\overline{w}_i)}`$ , unless it couples and decays into ordinary particles. Thus any potential with $`\alpha _i>4`$ will eventually fade away faster than radiation or matter unless particles are produced.

Naively, with more than one scalar field this dominance of the “false vacuum” can happen several times, with first $`\varphi _1`$ coming to dominate and then fade away, and then $`\varphi _2`$ coming to dominance and so on. Since the values of $`\varphi _i(init)`$, $`\lambda _i`$, and $`\alpha _i`$ determine when and how long each $`\varphi _i`$ dominates, an appropriate ensemble of such values could give a series of periods of vacuum domination and universal acceleration, followed by periods of radiation and/or matter domination, depending upon how each $`\varphi _i`$ decays and when the next field rises to dominance.

Very roughly the time of vacuum dominance for any $`\varphi _i`$ occurs when its energy density $`\rho _i=V_i(\varphi _i)+\frac{1}{2}\dot{\varphi _i}^2`$ equals the radiation (or matter) energy density $`\rho _{rad}8.6\times 10^5a^4\rho _{crit0}`$. That is, supposing for simplicity that the scalar field kinetic energy is not large compared with the potential energy, we have $`a_i(vacuumdom)\left(\lambda _i\varphi _{i}^{}{}_{}{}^{\alpha _i}(init)/8.6\times 10^5\rho _{crit0}\right)^{1/4}`$. Thus if the values of $`V_i(init)`$ are distributed over a wide range of values so that the $`a_i(vacuumdom)`$ are distributed over epochs from the Planck time at $`a=10^{27}`$ to today at $`a=1`$, then our suggested scenario might take place. Of course there are many ways that the $`\lambda _i`$ and the $`\varphi _i(init)`$ might be distributed to make this work, and we will not speculate on the precise form of the distribution since we have no understanding of it, except for possibly at two points: the inflaton in the very early Universe and the current vacuum dominated epoch.

For illustration purposes we consider two simple two-field cases, both of which satisfy the above constraints on $`w`$, fit the SN1A data as well as a cosmological constant model with $`\mathrm{\Omega }_V=\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, have nearly the same distance to CMB last scattering, and satisfy the constraints on the fraction of vacuum energy at decoupling and BBN. The first example, shown in Figure 1a, has $`\alpha _1=6`$, $`\lambda _1=10^{175}`$, $`V_1(init)=5\times 10^4`$, $`\alpha _2=4`$, $`\lambda _2=10^{125}`$, $`V_2(init)=10^8`$, where $`V_i(init)`$ is given in units of $`m_{Pl}^2\mathrm{Mpc}^2`$ to help with numerics and $`\varphi _i`$ is found from equation (1). In these units $`\rho _{crit0}=1.6\times 10^8`$, $`M_{\mathrm{Pl}}=1.94\times 10^{28}`$, $`\varphi `$ has units of ($`m_{Pl}^{1/2}\mathrm{Mpc}^{1/2}`$), and $`\lambda `$ has units $`(m_{Pl}^{1/2}\mathrm{Mpc}^{1/2})^{4\alpha }`$. This example has has the current accelerating expansion coming from $`\varphi _2`$, but earlier, at $`z1000`$, $`\varphi _1`$ started to become important, reaching more than 30% of the matter density at $`z300`$, and then fading away. Today this model gives $`\mathrm{\Omega }_V=.67`$ and $`w=.97`$. This example might be ruled out since the linear growth factor calculated from matter-radiation equality is about 30% smaller than for the $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ model (about the same as for an $`\mathrm{\Omega }_\mathrm{\Lambda }=0.78`$ model). Of course the value of $`V_1(init)`$ could be reduced to lessen the size of these effects.

The second example shows a typical period of complete vacuum dominance in the early Universe (figure 1b). Here the parameters used are $`\alpha _1=10`$, $`\lambda _1=10^{275}`$, $`V_1(init)=10^{16}`$, $`\alpha _2=4`$, $`\lambda _2=10^{124}`$, $`V_2(init)=10^8`$. Here the false vacuum dominates between redshifts of $`10^5`$ and $`10^7`$ and also near $`z=0`$, and is again a small contributor at the time of BBN and photon decoupling. Today this model gives $`\mathrm{\Omega }_V=.65`$ and $`w=.91`$. The linear growth factor is nearly the same as for a cosmological constant model. Note in all these plots we took the initial values of $`\dot{\varphi }_i(init)`$ to be zero at $`a_{init}=10^{10}`$, but very similar results obtain when we use equipartition of kinetic and potential energy by setting $`\dot{\varphi _i}(init)=\pm (2V_i(init))^{1/2}`$.

## Discussion

This paper suggests a class of models that might explain the cosmic coincidence problem without invoking the anthropic principle. The basic idea is that there is an ensemble of scalar (or pseudoscalar) fields with exponentially distributed parameters that cause them to dominate the universal expansion at random times throughout the history of the Universe. This paper does not attempt to solve the basic cosmological constant problem of why the energies of all the scalar fields sum to near zero today, but does allow for consistency between the present accelerating expansion and a zero sum. The simple polynomial models used here do require fine-tuning of their initial values and coupling constants. The initial values determine the time of domination and the values of the $`\lambda `$’s determine whether or not vacuum domination takes place and for how long. Note that in models with attractor potentials there is no need to fine-tune the initial values, and so these models have received the bulk of the attention of the community. In our scenario, we could probably remove this fine-tuning by considering an ensemble of hybrid potentials which have tracking behavior at early times and then asymptote to polynomials with $`\alpha >4`$ after dominance, but for simplicity’s sake we did not pursue this option in this paper. The basic point is that if one wants a non-anthropic solution to the cosmic coincidence problem, one probably wants an ensemble of fields with properties such as we discussed. This clue may motivate field or string theorists to find a way of naturally generating such an unusual set of initial conditions and coupling constants, or to find an ensemble of hybrid potentials that switch from tracking to decay in the right way.

In summary, we have not investigated any models in detail and have no first principle reason for why such an ensemble of fields should exist or why their parameters should be properly distributed, but are using the idea of a non-anthropic solution to the cosmic coincidence problem as the main motivation. However, we do note that the Higgs field in the standard model of particle physics is a scalar field that contributes to the universal vacuum energy and whose contribution must be cancelled by a large negative contribution from either a cosmological constant or from the vacuum energy of another field (or fields). Including the inflaton and the field giving rise to the current universal acceleration we thus probably have at least three important scalar fields contributing to standard cosmology, making our suggestion of a large ensemble of such fields more palatable.

This idea has several advantages, including a “unification” of the inflaton field and the field that currently causes universal acceleration, as well as perhaps detectable periods of vacuum dominance in the past, and predictions that the current era of vacuum dominance will end and that $`w`$ is not precisely unity. This last prediction might be testable in proposed experiments to measure the value and time derivative of $`w`$.

There are many possibilities and open questions that should be addressed. A general question is what kind of potentials and initial values can give rise to realistic implementations of this idea, and what kind of theories could give rise to such an ensemble of fields and initial data? Other important questions are probably difficult to address without a more specific field theoretic framework; for example, one should consider how each of these small periods of vacuum dominance end; are particles created? do substantial adiabatic or iso-curvature fluctuations result? is the Universe reheated? Is the power spectrum affected? Even in a model dependent way it would be interesting to explore how much and what kind of late time vacuum domination, or near domination, is allowed by current observations. It may even be that periods of vacuum dominance, or near dominance, could help the fit between theory and observation, or that the left over oscillating fields could make up some of the dark matter. In general one could ask what sets of fields and initial data could give rise to our current Universe and what observable effects would remain today?

## Acknowledgements

We thank Andy Albrecht, Neal Dalal, Saul Perlmutter, Ewan Stewart, Martin White, and Art Wolfe for discussion and valuable suggestions. This work was supported in part by the U.S. Department of Energy, under grant DOE-FG03-97ER40546.