text:
# Fractons and Luttinger liquids

## Abstract

We consider the concept of fractons as particles or quasiparticles which obey a specific fractal statistics in connection with a one-dimensional Luttinger liquid theory. We obtain a dual statistics parameter $`\stackrel{~}{\nu }=\nu +1`$ which is identified with the controlling parameter $`e^{2\phi }`$ of the Luttinger model. In this way, a bosonic system characterized by a fractal index $`i_f[h]=i_f[2]=1`$ is considered as a conformal field theory with central charge $`c[\nu =0]=1=i_f[2]`$ with a compactified radius $`R=\frac{1}{\sqrt{\stackrel{~}{\nu }}}=1`$. Thus, we have a mapping of a bosonic theory to a fermionic one and vice-versa, i.e. the duality symmetry $`\stackrel{~}{h}=3h`$ of the universal class $`h`$ of fractons defined in the interval $`1`$$`<`$$`h`$$`<`$$`2`$ is satisfied.

We have introduced in the concept of universal classes of particles or quasiparticles, labeled by a fractal parameter ( or Hausdorff dimension ) $`h`$, defined in the interval $`1`$$`<`$$`h`$$`<`$$`2`$. These particles termed fractons This word was used by Alexander and Orbach in another context. obey fractal statistics ( the distribution functions are fractal functions )

$`n[h]={\displaystyle \frac{1}{𝒴[\xi ]h}},`$ (1)
where the function $`𝒴[\xi ]`$ satisfies the equation

$`\xi =\left\{𝒴[\xi ]1\right\}^{h1}\left\{𝒴[\xi ]2\right\}^{2h},`$ (2)
and $`\xi =\mathrm{exp}\left\{(ϵ\mu )/KT\right\}`$, has the usual definition. The particles are collected in each class taking into account the fractal spectrum

$`h1=1\nu ,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<\nu <1;h1=\nu 1,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}<\nu <2;`$ (3)
$`etc.`$ (4)
and the spin-statistics relation $`\nu =2s`$.

Fractons are charge-flux systems which live in two dimensional multiply connected space and carry rational or irrational quantum numbers, as charge and spin. The fractal parameter $`h`$ is associated to the path of the quantum-mechanical particle and so, we have a quantum-geometrical description of the statistical laws of Nature. This means that the fractal characteristic of the quantum path, which reflects the Heisenberg uncertainty principle, is embodied in the fractal statistics expression. Thus, our formulation generalizes in a natural way the bosonic and fermionic distributions. Besides this, we verify that the classes $`h`$ of particles have a duality symmetry defined by $`\stackrel{~}{h}=3h`$, such that fermions ($`h=1`$) and bosons ($`h=2`$) are dual objects. As a result, we extract a fractal supersymmetry explicited by the theorem: If the particle with spin $`s`$ is within the class $`h`$, then its dual $`s+\frac{1}{2}`$ is into the class $`\stackrel{~}{h}`$.

Now in we have obtained a relation between fractons and conformal field theory(CFT)-quasiparticles ( edge excitations ), through the concept of fractal index associated with the classes $`h`$

$$i_f[h]=\frac{6}{\pi ^2}_{\mathrm{}(T=0)}^{1(T=\mathrm{})}\frac{d\xi }{\xi }\mathrm{ln}\left\{\mathrm{\Theta }[𝒴(\xi )]\right\}$$
(5)
where

$`\mathrm{\Theta }[𝒴]={\displaystyle \frac{𝒴[\xi ]2}{𝒴[\xi ]1}}`$ (6)
is the single-particle partition function of the universal class $`h`$ and hence

$`n[h]`$ $`=`$ $`\xi {\displaystyle \frac{}{\xi }}\mathrm{ln}\mathrm{\Theta }[𝒴].`$ (7)
We have verified a connection between the central charge $`c[\nu ]`$(a dimensionless number which characterizes conformal field theories in two dimensions) and the universal classes $`h`$ of particles. For that we considered the particle-hole duality $`\nu \frac{1}{\nu }`$ for integer-value $`\nu `$ of the statistical parameter which some systems of the condensed matter satisfy, for example, the Calogero-Sutherland model and others. The central charge for $`\nu `$ even is defined by

$`c[\nu ]=i_f[h,\nu ]i_f[h,{\displaystyle \frac{1}{\nu }}]`$ (8)
and for $`\nu `$ odd is defined by

$`c[\nu ]=2\times i_f[h,\nu ]i_f[h,{\displaystyle \frac{1}{\nu }}],`$ (9)
where $`i_f[h,\nu ]`$ means the fractal index of the universal class $`h`$ which contains the particles with distinct spin values which obey a specific fractal statistics.

In this short note we show the connection between fractons and Luttinger liquids, taking into account the parameters which characterize these models. They are one-dimensional interacting many-body systems, which at low-energy or low-temperature are characterized by a controlling parameter, $`e^{2\phi }`$.

We also have obtained a Fermi velocity given by $`v_F=\frac{v_N}{\stackrel{~}{\nu }}`$, where $`\stackrel{~}{\nu }=\nu +1`$ ($`\nu 0`$) is a dual statistical parameter and $`v_N=\frac{2\mathrm{}}{mL}c[\nu ]`$ is an excitation related to the change in the particle number which depends on the central charge. The quantities $`\mathrm{}`$ is the Planck constant, $`m`$ is the mass of particle and $`L`$ is the length of the edge. Another excitation related to the current is defined as $`v_J=\frac{v_F}{\stackrel{~}{\nu }}`$. Haldane has just considered a relation between these two types of excitation for characterize the Luttinger liquid $`v_F=\sqrt{v_Nv_J}`$. Thus, our definition satisfies this condition and on one way, we identify the controlling parameter with the dual statistical parameter $`\stackrel{~}{\nu }=e^{2\phi }`$, such that $`v_N=v_Fe^{2\phi }`$ and $`v_J=v_Fe^{2\phi }`$.

Now, we observe that $`\stackrel{~}{\nu }=\nu +1`$($`\nu 0`$), satisfies the duality symmetry, i.e. we can map a fermionic theory to a bosonic one and vice-versa. We also know that the fractal index $`i_f[h]`$ coincides with the central charge $`c[\nu ]`$, for $`0\nu 1`$ and we have distinguished two concepts of central charge, one related to the universal classes $`h`$ and the other related to the particles which belong to these classes. Thus, the fractal index for $`\nu =0`$( boson, $`\nu `$ even ) is $`i_f[h,\nu ]=i_f[2,0]=1`$ and for the dual statistical parameter $`\stackrel{~}{\nu }=1`$ ( fermion, $`\stackrel{~}{\nu }`$ odd ) is $`i_f[\stackrel{~}{h},\stackrel{~}{\nu }]=i_f[1,1]=\frac{1}{2}`$. In this way, we define for a bosonic theory $`i_f[2]=1`$ identified with a conformal field theory with central charge $`c[\nu =0]=1`$, a compactified radius $`R`$ determined by the dual statistical parameter $`R^2=\frac{1}{\stackrel{~}{\nu }}=e^{2\phi }`$ and $`\phi =\mathrm{ln}R`$. The radius $`R`$ expresses an angular invariance of the theory.

Finally, we would like to emphasize that our results were obtained independently of any bosonization idea as discussed in. In contrast, we have used a symmetry principle, i.e. the concept of duality symmetry between the universal classes $`h`$ of particles. On the other hand, our definition of fractal statistics differs of the notion of exclusion statistics, given that in our approach we consider ab initio the spin-statistics relation $`\nu =2s`$ for particles with rational or irrational spin quantum number $`s`$. This claim is supported by the fractal spectrum which is a mirror symmetry in our formulation. Another consequence is the fractal supersymmetry which is realized, for example, in the context of fractional quantum Hall effect.


format: text

text:
# Energy separation of single-particle and continuum states in a 𝑆=1/2 weakly-coupled chains antiferromagnet.

\[

## Abstract

Inelastic neutron scattering is used to study transverse-polarized magnetic excitations in the quasi-one-dimensional $`S=1/2`$ antiferromagnet BaCu2Si2O7, where the saturation value for the Néel order parameter is $`m_0=0.12`$ $`\mu _\mathrm{B}`$ per spin. At low energies the spectrum is totally dominated by resolution-limited spin wave-like excitations. An excitation continuum sets in above a well-defined threshold frequency. Experimental results are discussed in the context of current theories for weakly-interacting quantum half-integer spin chains.

\]

One of the outstanding problems in quantum magnetism today is that of dimensional crossover in systems composed of weakly-coupled spin chains. In particular, for $`S=1/2`$ Heisenberg antiferromagnets (AFs), the interplay between classical dynamics in Néel-ordered three-dimensional materials (single-particle order-parameter excitations, or spin waves) and critical dynamics of the quantum-disordered $`S=1/2`$ one-dimensional (1D) model (2-spinon excitation continuum ) is not fully understood. In quasi-1D $`S=1/2`$ systems, for arbitrary weak inter-chain interactions, long-range order is restored. However, according to recent quantum-mechanical calculations, as long as the 3D order parameter remains small, the excitation spectrum has a unique dual nature. At low frequencies the diffuse continuum is “cleaned up” and replaced by long-lived single-particle “magnon” excitations. What remains of the critical dynamics in isolated chains, is now seen at higher energies as a series of multi-magnon states. The lowest-energy contribution is from two-magnon processes. At the 1D AF zone-center the continuum is thus expected to set in at twice the characteristic magnon energy, and should be separated from the single-particle part of the spectrum by a “second gap”.

To date, this separation of single-particle and continuum dynamics has not been clearly observed experimentally in any quasi-1D material. Among other technical obstacles, is the limited choice of suitable model compounds. Most of what is experimentally known of coupled $`S=1/2`$ chains comes from neutron scattering studies of KCuF3 . Here the ordered moment is rather large, making the spin waves very intense at low energies, and the continuum contribution difficult (though not impossible) to isolate. Materials like Sr2CuO3  or SrCuO2, on the other hand, are very good 1D systems, and it is the spin waves that are hard to identify. Below we report inelastic neutron scattering studies of spin dynamics in the $`S=1/2`$ quasi-1D antiferromagnet BaCu2Si2O7. Using neutron setups with complimentary resolution characteristics, we clearly observe the frequency-separation of single-particle and continuum contributions. In doing so, we pay special attention to the polarization of magnetic excitations, and compare the results to existing theoretical models.

The recently characterized BaCu2Si2O7 (orthorhombic space group $`Pnma`$, $`a=6.862`$ Å, $`b=13.178`$ Å, $`c=6.897`$ Å) appears to be an ideal model material for the present study. Compared to KCuF3, it has a stronger 1D character: with an in-chain (along the $`c`$-axis) Heisenberg exchange constant $`J=279`$ K, the Néel temperature is only $`T_\mathrm{N}=9.2`$ K. The structure of the magnetically ordered phase was previously guessed from bulk data and an analysis of a single magnetic Bragg reflection . As part of the present study, we performed a comprehensive neutron diffraction study of the magnetic structure at $`T=1.5`$ K using the D23 single-crystal diffractometer at Institut Laue-Langevin (ILL) . Overall, 35 inequivalent magnetic Bragg intensities were measured. The spin arrangement was found to be exactly as proposed in Ref. . The obtained value for the ordered moment $`m_0=0.12`$ $`\mu _\mathrm{B}`$ is much smaller then in KCuF3 ($`m_00.5`$ $`\mu _\mathrm{B}`$), yet large enough to expect a measurable spin-wave contribution to inelastic scattering.

Low-energy magnetic excitations in BaCu2Si2O7 were studied in high-resolution inelastic neutron scattering experiments on a $`5\times 5\times 50`$ mm3 BaCu2Si2O7 single-crystal sample at the IN14 cold-neutron spectrometer at ILL . Measurements were performed in the $`(0,k,l)`$ scattering plane with a fixed final neutron energy of 3 meV. Typical constant-$`Q`$ and constant-$`E`$ scans collected near the 1D AF zone-center $`l=1`$ are shown in Fig. 1. Figure 2 shows a series of constant-$`Q`$ scans for different momentum transfers perpendicular to the chain-axis. The polarization dependence of the magnetic neutron scattering cross section ensures that within this range, to a good approximation, only fluctuations of spin components perpendicular to the chains (and thus the ordered moment, which is parallel to $`c`$) are observed. Very similar data sets (not shown) were collected using a comparable setup on the TASP cold-neutron 3-axis spectrometer at Paul Scherrer Institut, Switzerland, for momentum transfers along the $`(h,0,1)`$ and $`(h,h,1)`$ reciprocal-space rods .

The main feature seen in all constant-$`Q`$ scans is a sharp asymmetric peak that we attribute to single-magnon excitations. The intensity onset on the low-energy side is very steep, and suggests the peaks are resolution-limited. The extended high-energy “tail” is due to a steep dispersion along the chain axis and a non-zero wave vector resolution in that direction. The data were analyzed using a single-mode approximation (SMA) cross section, derived from the quantum chain-Mean Field (chain-MF) model . The dispersion measured along $`(0,k,1)`$, $`(h,0,1)`$ and $`(h,h,1)`$, we found that at least three inter-chain exchange constants are needed, between nearest-neighbor spins along the $`(1,0,0)`$, $`(0,1,0)`$ and $`(1,1,0)`$ real-space directions, that we denote as $`J_x`$, $`J_y`$ and $`J_3`$, respectively. The SMA cross section was convoluted with the calculated experimental resolution function and used to fit all measured scans to determine $`J_x`$, $`J_y`$ and $`J_3`$. The in-chain coupling constant $`J`$ was fixed at $`J=24.1`$ meV, as previously deduced from bulk susceptibility curves . An excellent global fit to all scans is obtained with $`J_x=0.463(2)`$ meV, $`J_y=0.161(1)`$ meV and $`2J_3=0.145(1)`$ meV. Simulations based on these values are shown in solid lines in Figs. 1,2. The resulting dispersion relation along the $`(0,k,1)`$ direction is shown in a solid line in the $`(\mathrm{}\omega )`$-$`k`$ plane in Fig. 2. Additional analysis also gives an upper estimate for the intrinsic excitation width: $`0.03`$ meV. These results conclusively demonstrate that up to about 4 meV energy transfer the magnetic excitation spectrum is fully accounted for by the single-particle picture, in the entire range of momentum transfers perpendicular to the chain-axis.

Extending the cold-neutron study to higher energy transfers was hindered by certain geometrical constrains imposed by the design of IN14 spectrometer. Instead, the intermediate-energy spectrum was studied using the IN22 thermal-neutron spectrometer at ILL using 14.7 meV final-energy neutrons. Most of the data were collected in constant-$`E`$ scans along the chain-axis, in the vicinity of the $`(0,0,3)`$ 1D AF zone-center. The background originating from the empty Al sample-holder was measured separately for each scan. Typical background-subtracted data sets are presented in Fig. 3. None of these show the distinct 2-peak structure as at 3 meV energy transfer (Fig. 1c). To verify that this is not a resolution effect, we performed a least-squares fit of the SMA cross section, convoluted with the calculated resolution, to each measured profile. Very poor fits are obtained, as shown in dashed lines in Fig. 3. It is clear that if the observed intensity was due to a single mode, two separate peaks would have been easily resolved. For lack of a convenient analytical result for interacting chains, we fitted the measured scans to the Müller-ansatz continuum form , known to work well for isolated chains. The agreement in this case is much better (solid lines in Fig. 3). For contrast, the dashed line in Fig. 1c shows the best Müller-ansatz fit to the 3 meV high-resolution const-$`E`$ scan. From this comparative analysis we conclude that, unlike at low energies, a large contribution to the observed scattering at high energies must come from continuum excitations.

Additional evidence was obtained in analyzing the measured intensities. In the chain-MF model the intensity of magnon scattering has the same $`1/\omega `$ dependence as for classical antiferromagnetic spin waves, in full agreement with the bulk of our cold-neutron data. Having confirmed the single-mode picture at low energies, we exploited this scaling relation to estimate spin wave intensities at higher energy transfers. A normalization factor was obtained by fitting the SMA cross section to the intensity measured using the thermal setup at $`𝒒=(0,0,3)`$ between 2 and 3 meV transfer, where, according to our cold-neutron results, the SMA picture is still valid. The thus calculated single-mode contributions are represented by shaded areas in Fig. 3.

Clearly, magnon excitations account for only a minor fraction of the dynamic susceptibility at high frequencies. To quantify this behavior, in Fig. 4a we plot the measured $`q_z`$-integrated intensity $`I_z`$, scaled by energy transfer $`\mathrm{}\omega `$, as a function of the latter. The advantage of this integration procedure that the actual wave-vector dependence of continuum scattering becomes unimportant. The expected spin wave contribution is practically constant above 3 meV in this plot (Fig. 4a, dashed line). This follows from that in a fixed-final energy configuration one directly measures $`S(𝒒,\omega )`$ without any additional energy-dependent scaling factors ($`\lambda /2`$ effects in the monitor are expected to be less than 5%, thanks to the use of a thermal guide). In contrast, experimentally, $`I_z\omega `$ increases substantially with frequency.

The continuum fraction can be isolated by subtracting the extrapolated SMA part from the experimental data. The result is shown in Fig. 4b. The energy resolution being about 1.5 meV, we conclude that the continuum sets in between 3.5 and 5.5 meV energy transfer. We have verified that this result is quite robust, and remains valid even if we assume that there is a 50% systematic error in the estimate of the single-mode contribution. From Fig. 4b we see that above the threshold the $`q_z`$-integrated continuum intensity is only slightly energy-independent. Looking back at the cold-neutron measurements, we find that they actually do contain direct evidence of continuum scattering as well. Indeed, the increase of intensity seen around 4.5 meV in Fig. 1 is not accounted for by the SMA picture, and thus represents the onset of the continuum, in full agreement with the thermal-neutron results.

To understand the observed long-lived and continuum excitations, we shall recall the simple yet profound physical picture provided by the quantum chain-MF theory . For an isolated chain in a staggered field, the magnons have an energy gap $`\mathrm{\Delta }`$ that scales as $`H_\pi ^{2/3}`$ . In the coupled-chains case, magnon dispersion perpendicular to the chain axis leads to a softening of the gap at the 3D zone-center for transverse-polarized excitations. These gapless modes correspond to conventional spin waves in the classical spin wave theory (SWT), and are the long-lived excitations that we see in BaCu2Si2O7 at low energies. Excitations at exactly $`\mathrm{\Delta }`$ can be observed in points of reciprocal space where inter-chain-interactions cancel out at the MF level. In BaCu2Si2O7 this occurs at $`𝒒=(0.5,0.5,1)`$. Using the measured exchange constants we can estimate $`\mathrm{\Delta }=2.37`$ meV. It is straightforward to verify that the experimental values for $`\mathrm{\Delta }`$, $`J`$, $`m_0`$ and $`T_\mathrm{N}`$ satisfy the relations predicted by the chain-MF model  remarkably well.

A unique quantum-mechanical feature of a half-integer Heisenberg AF chain in a staggered field, totally absent in conventional SWT, is a longitudinal mode, i.e., a magnon excitation polarized parallel to the direction of static staggered magnetization. The existence of such an excitation was recently confirmed experimentally in KCuF3 . Focusing on transverse spin correlations in the present study, we do not directly observe this branch in BaCu2Si2O7. It does, nonetheless, play an important role in the observed transverse continuum. In the chain-MF and field-theoretical models for weakly-interacting chains, low-energy continuum excitations are seen as multi-magnon states . The existence of a longitudinal mode allows two-particle excitations in the transverse channel, composed of one longitudinal and one transverse magnon. The threshold for such excitations at the 1D AF zone-center is at twice the characteristic magnon energy: $`\mathrm{\Delta }_c=2\mathrm{\Delta }`$. This prediction is consistent with our experimental results for BaCu2Si2O7, where $`\mathrm{\Delta }_c`$ is between 3.5 and 5.5 meV, and $`2\mathrm{\Delta }=4.8`$ meV.

An alternative model that predicts multi-magnon continua is SWT with kinematic corrections, which explicitly takes into account non-linearities in the spin wave Hamiltonian. Continuum excitations observed in KCuF3 are remarkably well described by this model . In BaCu2Si2O7, however, the ordered moment is much smaller, and the ground state resembles the Néel phase, i.e., the starting point of SWT calculations, considerably less. For this material SWT completely fails at the quantitative level. For example, the spin wave correction to sublattice magnetization in BaCu2Si2O7 that we calculated using the known dispersion bandwidths is well over 100%. For the transverse-polarized excitation continuum, the inadequacy of SWT becomes even more apparent. Indeed, the longitudinal mode being absent in this model, two-magnon transverse excitations are prohibited by selection rules. The lowest-energy contribution is from three-magnon states. The lower bound of the continuum should thus be at roughly $`3\mathrm{\Delta }`$. The 3-magnon SWT contribution can be exactly evaluated, provided the spin wave dispersion is known. Following the recipe given in Ref. , the spin wave bandwidths measured in BaCu2Si2O7 were used to calculate the 3-magnon profile at $`𝒒=(0,0,3)`$. The result (arbitrary scaling) is shown in Fig. 4b in a dash-dot line. The calculated continuum onset is clearly at a larger energy, by 2 to 3 meV, than observed experimentally.

In summary, we find that the transverse excitation spectrum in BaCu2Si2O7 is divided into two well-defined regions. Below 4 meV the weight is consolidated into long-lived single-particle spin wave-like excitations, and no sign of multi-particle continuum is found. These excitations carry only a small fraction of the total spectral weight at higher energies, where an excitation continuum becomes the dominant contribution. While further studies of the lower continuum bound are needed, the present results suggest a threshold between 3.5 and 5.5 meV. For this strongly 1D material the conventional spin wave- theoretical description becomes inadequate, even if kinematic corrections are taken into account. In contrast, the chain-MF model, based on quantum-mechanical properties of individual chains, is in very good agreement with experiment.

We would like to thank Dr. P. Böni (PSI Villigen) for allowing us to refer to the yet unpublished TASP data, Dr. L.-P. Regnault (CEA Grenoble) for his assistance with experiments at ILL, Prof. A. Tsvelik, Prof. R. A. Cowley(Oxford University) and Dr. I. Zaliznyak (BNL) for illuminating discussions, and Mr. R. Rothe (BNL) for technical support. This work is supported in part by the U.S. -Japan Cooperative Program on Neutron Scattering, Grant-in-Aid for COE Research “SCP coupled system” from the Ministry of Education, Science, Sports, and Culture. Work at Brookhaven National Laboratory was carried out under Contract No. DE-AC02-98CH10886, Division of Material Science, U.S. Department of Energy. One of the authors (M. K.) is supported by a TMR-fellowship from the Swiss National Science Foundation under contract no. 83EU-053223.


text:
# Untitled Document

CASIMIR ENERGY IN NON-COVARIANT GAUGES

Giampiero Esposito1,2, Alexander Yu. Kamenshchik3,4 and K. Kirsten5 1Istituto Nazionale di Fisica Nucleare, Sezione di Napoli, Complesso Universitario di Monte S. Angelo, Via Cintia, Edificio N’, 80126 Napoli, Italy 2Università degli Studi di Napoli Federico II, Dipartimento di Scienze Fisiche, Complesso Universitario di Monte S. Angelo, Via Cintia, Edificio N’, 80126 Napoli, Italy 3L.D. Landau Institute for Theoretical Physics, Russian Academy of Sciences, Kosygina Str. 2, Moscow 117334, Russia 4Landau Network-Centro Volta, Villa Olmo, Via Cantoni 1, 22100 Como, Italy 5The University of Manchester, Department of Physics and Astronomy, Theory Group, Schuster Laboratory, Oxford Road, Manchester M13 9PL, England

Abstract. The zero-point energy of a conducting spherical shell is studied by imposing the axial gauge via path-integral methods, with boundary conditions on the electromagnetic potential and ghost fields. The coupled modes are then found to be the temporal and longitudinal modes for the Maxwell field. The resulting system can be decoupled by studying a fourth-order differential equation with boundary conditions on longitudinal modes and their second derivatives. Complete agreement is found with a previous path-integral analysis in the Lorenz gauge, and with Boyer’s value. This investigation leads to a better understanding of how gauge independence is achieved in quantum field theory on backgrounds with boundary.

The quantization programme of fundamental interactions in non-covariant gauges has been always of great importance, as is well known from the work in Refs. . On the other hand, recent progress in the measurement of Casimir energies \[3–6\] motivates new theoretical efforts, including the comparison of covariant and non-covariant gauges for the evaluation of Casimir energies.

We here rely entirely on the path-integral analysis of Ref. , and avoid review for length reasons. Recall therefore that, in the case of a spherical shell of radius $`R`$, the axial gauge forces normal modes of the potential to vanish everywhere, jointly with ghost modes. One is then left with transverse modes $`T_l(r)=\sqrt{\pi /2}\sqrt{r}I_{l+1/2}(Mr)`$ which obey homogeneous Dirichlet conditions at $`r=R`$, while temporal and longitudinal modes (denoted by $`a_l`$ and $`c_l`$, respectively), obey the coupled system

$$\left[\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\frac{l(l+1)}{r^2}\right]a_l\frac{Ml(l+1)}{r^2}c_l=0,$$
$`(1)`$
$$\left[\frac{d^2}{dr^2}M^2\right]c_lMa_l=0.$$
$`(2)`$
First, we now point out that this system can be decoupled if $`a_l`$ is expressed from Eq. (2) and the result is inserted into Eq. (1). On setting $`yMr`$, this leads to the fourth-order equation

$$\left[\frac{d^4}{dy^4}+\frac{2}{y}\frac{d^3}{dy^3}\left(1+\frac{l(l+1)}{y^2}\right)\frac{d^2}{dy^2}\frac{2}{y}\frac{d}{dy}\right]c_l(y)=0.$$
$`(3)`$
Moreover, we know from Ref.  that, at $`r=R`$, both $`a_l`$ and $`c_l`$ should vanish. By virtue of Eq. (2) this implies that, at $`r=R`$, both $`c_l`$ and its second derivative should vanish. The origin $`r=0`$ is a regular singular point of our boundary-value problem, and hence we require that both $`a_l`$ and $`c_l`$ should vanish therein. By virtue of Eq. (2), this implies that also the second derivative of $`c_l`$ with respect to $`y`$ vanishes at $`y=0`$. Now we set $`F_l(y)c_l^{}(y)`$, and cast Eq. (3) in the form

$$\frac{1}{y^2}\frac{d}{dy}\left[y^2\left(\frac{d^2}{dy^2}\left(1+\frac{l(l+1)}{y^2}\right)\right)F_l(y)\right]=0.$$
$`(4)`$
On denoting by $`C`$ a constant, the desired $`F_l(y)`$ solves therefore the inhomogeneous equation

$$\left[\frac{d^2}{dy^2}\left(1+\frac{l(l+1)}{y^2}\right)\right]F_l(y)=\frac{C}{y^2},$$
$`(5)`$
whose regular solution reads ($`d`$ being a constant)

$$F_l(y)=CF_{sp}(y)+d\sqrt{\pi /2}\sqrt{y}I_{l+1/2}(y),$$
$`(6)`$
where $`F_{sp}`$ is chosen to be a particular solution of Eq. (5) with $`C=1`$ whose first derivative vanishes at $`0`$ and at $`bMR`$. The vanishing of $`c_l`$ at the boundary fixes the relation between $`C`$ and $`d`$, i.e.

$$C=d\frac{_0^b𝑑z\sqrt{\pi /2}\sqrt{z}I_{l+1/2}(z)}{_0^b𝑑zF_{sp}(z)}.$$
$`(7)`$
Eventually, the vanishing of $`c_l^{\prime \prime }`$ at the boundary leads to

$$I_{l+1/2}^{}+\frac{1}{2b}I_{l+1/2}(b)=0,$$
$`(8)`$
because, by construction, $`CF_{sp}`$ is a particular solution of Eq. (5) whose first derivative vanishes at the boundary. Equation (8), jointly with the vanishing at the boundary of transverse modes, yields the same set of eigenvalue conditions, with the same degeneracies, found for the interior problem in the Lorenz gauge in Ref. . Thus, complete agreement  with Boyer’s value  for the Casimir energy is recovered.

Although non-covariant gauges break relativistic covariance, they make it possible to decouple Faddeev–Popov ghosts from the gauge field. Thus, ghost diagrams do not contribute to cross-sections and need not be evaluated, and this property has been regarded as the main advantage of non-covariant gauges . In the case of Casimir energies, the ghost field is forced to vanish everywhere by virtue of the boundary conditions appropriate for the axial gauge , and the analysis of temporal and longitudinal modes, although rather involved, has been proved to lead to the same Casimir energy  for the interior problem in the Lorenz gauge . The particular solution of the inhomogeneous equation (5) plays a non-trivial role in ensuring that, despite some technical points, the resulting Casimir energy is the same as in the Lorenz gauge, hence proving explicitly the equivalence of a covariant and a non-covariant gauge for a conducting spherical shell. A better theoretical understanding of gauge independence in quantum field theory has been therefore gained, after the encouraging experimental progress of recent years in the measurement of Casimir forces \[3–6\]. Acknowledgments. Klaus Kirsten thanks Stuart Dowker for interesting discussions. The work of Klaus Kirsten has been supported by EPSRC, grant no GR/M08714.

References

 G. Leibbrandt, Rev. Mod. Phys. 59, 1067 (1987).

 G. Leibbrandt, Noncovariant Gauges (World Scientific, Singapore, 1994).

 S. K. Lamoreaux, Phys. Rev. Lett. 78, 5 (1997).

 U. Mohideen and A. Roy, Phys. Rev. Lett. 81, 4549 (1998).

 A. Roy and U. Mohideen, Phys. Rev. Lett. 82, 4380 (1999).

 A. Roy, C. Y. Liu and U. Mohideen, Phys. Rev. D 60, 111101 (1999).

 G. Esposito, A. Yu. Kamenshchik and K. Kirsten, Int. J. Mod. Phys. A 14, 281 (1999).

 G. Esposito, A. Yu. Kamenshchik and K. Kirsten, “Casimir Energy in the Axial Gauge” (DSF preprint 2000/21, HEP-TH 0006220).

 T. H. Boyer, Phys. Rev. 174, 1764 (1968).


text:
# HIGH ENERGY GAMMA-RAY OBSERVATIONS OF THE CRAB NEBULA AND PULSAR WITH THE SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT

## 1 INTRODUCTION

The Crab Nebula is the most studied object in the gamma-ray sky. Numerous ground-based experiments have detected emission from the Crab at energies from $`300\mathrm{GeV}`$ through $`50\mathrm{TeV}`$. In addition, the EGRET detector on the Compton Gamma Ray Observatory has measured the spectrum of unpulsed gamma rays from the Crab up to an energy of $`10\mathrm{GeV}`$ (Nolan et al., 1993). High energy gamma rays are believed to result from inverse Compton scattering of relativistic electrons on synchrotron photons, thermal dust photons, and cosmic microwave photons. Such inverse Compton models have been described by de Jager & Harding (1992), Atoyan & Aharonian (1996), and de Jager et al. (1996). Under this scenario, high energy electrons are accelerated at the pulsar wind termination shock. They then propagate through the nebula, producing both synchrotron and inverse Compton emission.

In addition to unpulsed gamma rays, EGRET has also seen pulsed gamma-ray emission from the Crab (Nolan et al., 1993; Ramanamurthy et al., 1995; Fierro et al., 1998), which is assumed to be produced by the pulsar. Pulsed emission has not been seen at higher energies accessible to ground-based experiments (Lessard et al., 1999). The measurements require that pulsed emission cut off somewhere between $`10\mathrm{GeV}`$ and $`250\mathrm{GeV}`$. Polar Cap models of pulsar emission (Daugherty & Harding, 1982) predict sharp cutoffs near $`10\mathrm{GeV}`$, whereas Outer Gap models (Cheng, Ho, & Ruderman, 1986; Romani, 1996) predict pulsed emission as high as $`50\mathrm{GeV}`$.

No previous gamma-ray detector has had sensitivity to gamma rays at energies between $`10\mathrm{GeV}`$ and $`300\mathrm{GeV}`$. Nonetheless, there are strong motivations to observe the Crab in this energy range, not only to test inverse Compton models at previously unexplored energies, but also to search for pulsed gamma-ray emission that could distinguish between pulsar emission models. Satellite experiments have been limited to energies below $`10\mathrm{GeV}`$ by low statistics because of their small apertures and the rapidly falling flux of gamma rays. Ground-based experiments, which detect gamma rays by the Cherenkov light produced in extensive air showers initiated by the gamma ray, typically have energy thresholds of $`300\mathrm{GeV}`$. The energy threshold is limited by the ability of the instrument to distinguish dim flashes of Cherenkov light amidst the fluctuations in night sky background light. From signal-to-noise considerations, the energy threshold of an atmospheric Cherenkov telescope can be shown to behave like

$$E_{th}\sqrt{\frac{\mathrm{\Phi }\mathrm{\Omega }\tau }{ϵA}}.$$
(1)
Here, $`\mathrm{\Phi }`$ is the flux of night sky background light, $`\mathrm{\Omega }`$ is the field of view of the instrument, $`\tau `$ is the length of the electronic trigger coincidence window, $`ϵ`$ is the efficiency for detecting Cherenkov photons, and $`A`$ is the total mirror collection area of the instrument (Weekes, 1988). In particular, the number of Cherenkov photons collected from a shower will scale with the mirror area $`A`$, while fluctuations in the night sky background light will increase with the square root of the light collection, and hence with $`\sqrt{A}`$. Thus, the energy threshold of an atmospheric Cherenkov detector decreases approximately in proportion to the inverse square root of the total mirror area.

The Solar Tower Atmospheric Cherenkov Effect Experiment (STACEE) is a new atmospheric Cherenkov detector for gamma-ray astronomy which uses the large mirror area of a solar power research facility to achieve a low energy threshold. STACEE uses an array of steerable mirrors, called heliostats, to collect Cherenkov light produced in extensive air showers. Heliostats focus the light onto secondary mirrors on a central tower, which image the light onto photomultiplier tube cameras (see Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT). Previous tests demonstrate that heliostat arrays can be used to detect Cherenkov light from extensive air showers (Ong et al., 1996; Chantell et al., 1998). Here we report on observations of the Crab Nebula and Pulsar taken with a prototype version of the STACEE instrument. This prototype, called STACEE-32, used 32 heliostats and a preliminary electronic trigger system. The full instrument incorporates 64 heliostats, and is presently under construction. It will start regular observations in the fall of 2000.

## 2 DETECTOR

STACEE is located at the National Solar Thermal Test Facility (NSTTF), at Sandia National Laboratories in Albuquerque, NM, USA (34.962$`^{}`$N, 106.509$`^{}`$W, 1705$`\mathrm{m}`$ above sea level). The NSTTF is a solar power research facility of the U.S. Department of Energy. It includes 212 heliostat mirrors, each $`37\mathrm{m}^2`$ in area, arrayed across a $`300\mathrm{m}\times 120\mathrm{m}`$ area (see Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT). For solar power research, heliostats track the Sun and concentrate sunlight onto a central receiver tower. At night, STACEE-32 used a subset of the heliostats to track a gamma-ray source. These heliostats reflected Cherenkov light towards two $`1.9\mathrm{m}`$ diameter secondary mirrors near the top of the tower. The secondary mirrors concentrated light from the heliostats onto phototube cameras. Each photomultiplier tube (PMT) received light from a single heliostat in the array. The heliostats were tilted inwards to track “shower maximum,” the point of maximum air shower development (about $`11\mathrm{km}`$ above sea level for vertical showers). Details of the site, heliostats, secondary mirrors, and cameras are given elsewhere (Chantell et al., 1998; Oser, 2000).

Each camera contained 16 photomultiplier tubes and received light from a separate portion of the heliostat field (see Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT). Solid acrylic light collectors based on the Dielectric Total Internal Reflecting Concentrator design (Ning, Winston, & O’Gallagher, 1987) were placed in front of the PMTs. The collectors concentrated incident light onto $`5\mathrm{cm}`$ diameter photomultiplier tubes and restricted the field of view of each PMT on the secondary mirror. The field of view of a heliostat was determined by the diameter viewed by its PMT on the secondary mirror, divided by the distance between the secondary mirror and the heliostat. The fields of view were adjusted so that each heliostat saw a solid angle of approximately $`1.2\times 10^4`$$`\mathrm{sr}`$ on the sky.

The photomultiplier tubes (Photonis XP2282B model) had bi-alkali photocathodes with peak sensitivity near $`400\mathrm{nm}`$ and excellent time response ($`1\mathrm{ns}`$ resolution). The single photoelectron rate from night sky noise on each tube was $`1.5\mathrm{GHz}`$, and the typical operating gain was $`1.1\times 10^5`$. Signals from the PMTs were AC-coupled and amplified ($`\times 100`$) and then routed to electronics racks inside the tower.

In order to achieve a low energy threshold, the PMT signals must be combined in a short time coincidence (see eq. ). The heliostats are at varying distances from the tower, and light from an air shower strikes different heliostats at different times, depending on the orientation of the shower. Therefore, signals from different heliostats must be delayed by varying amounts in order to put them into time coincidence. STACEE-32 used a two-level digital trigger to form this coincidence. The 32 heliostats were divided into four compact subclusters of eight heliostats each, as indicated in Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT. PMT signals were discriminated, and the discriminator outputs were individually delayed to put the eight signals within a subcluster into time coincidence with one another. If five of eight tubes fired within a 12.5$`\mathrm{ns}`$ window, then that subcluster triggered. The four subcluster triggers were themselves delayed and put into time coincidence. Three of four subclusters firing within a 15$`\mathrm{ns}`$ window resulted in an experiment trigger, which initiated readout. Thus, a total of at least fifteen heliostats had to fire to trigger the detector. The discriminator threshold was chosen so that the accidental trigger rate due to night sky background fluctuations was negligible ($`<0.002`$$`\mathrm{Hz}`$). The typical discriminator level was $`170\mathrm{mV}`$, which was approximately 5.5 times the mean single photoelectron amplitude. The mean discriminator counting rate in each channel was $`2\mathrm{MHz}`$.

Multi-hit time-to-digital converters (Lecroy 3377 TDC) measured the arrival times of pulses from each discriminator channel. Charge-integrating analog-to-digital converters (Lecroy 2249SG ADC) on 24 channels measured the pulse charge over a 37$`\mathrm{ns}`$ integration gate. The long delay times in the trigger formation ($`1\mu \mathrm{s}`$) required that the ADC inputs be first routed through $`1200\mathrm{ns}`$ of high quality RG213 coaxial cable. These cables delayed the arrival of the pulses at the ADC until after the trigger had formed. Scaler counters recorded the PMT and subcluster trigger rates, and a GPS clock recorded the Universal Time for each event trigger. The PMT currents were also read out for every event. An experiment trigger asserted a veto, which inhibited further triggers until the event was read out and the veto was cleared by the data acquisition (DAQ) system. Counters measured the livetime of the experiment (the fraction of time for which triggers were enabled). The experiment was also triggered with an external pulser every two seconds to collect a sample of non-Cherenkov events for calibration purposes. All readout and control was done via an interface to a Silicon Graphics workstation. Details of the electronics are given in Oser (2000).

Calibration systems for STACEE-32 included a nitrogen laser and an LED flasher circuit. Optical fibers carried pulses from the laser to small wavelength shifter plates mounted at the center of each secondary mirror. The laser light pulsed all phototubes in each camera simultaneously. A filter wheel was used to attenuate the laser light level, and the slewing responses of the PMTs (time versus pulse amplitude) were so determined. On an occasional basis, the LED flasher was placed in front of each PMT to measure its gain in situ. Phototubes, delay modules, and other electronics systems were also calibrated extensively in the lab before installation. In addition, the orientations and status of all heliostats were logged during tracking, and any malfunctions were noted offline and the affected data were identified.

## 3 EVENT RECONSTRUCTION AND BACKGROUND REJECTION

Although the energy threshold of an atmospheric Cherenkov detector is ultimately limited by night sky background light, the sensitivity of the detector is limited by a different kind of background, namely, cosmic ray air showers. STACEE-32’s high-multiplicity trigger coincidence effectively eliminated all accidental triggers due to fluctuations in night sky light. The high multiplicity of the trigger also rejected approximately 98% of all Cherenkov showers initiated by cosmic rays (see § 3.2). The remaining cosmic ray triggers form the background from which gamma-ray events must be distinguished.

Traditional atmospheric Cherenkov telescopes use a single mirror with a pixellated phototube camera. Each element in the camera is mapped to a different angular region on the sky, and thus one records an image of the Cherenkov shower in the sky. However, STACEE is not an imaging detector. Instead, the experiment is a lateral array which samples the Cherenkov wavefront at many locations within the light pool. In this respect, event reconstruction must be handled quite differently than for an imaging Cherenkov telescope. STACEE records the arrival time and photon density of the Cherenkov wavefront at each heliostat for every event. From these quantities the shape and lateral density profile of the Cherenkov light pool are determined. For STACEE-32, the charge resolution on the ADC measurements was limited by the high night sky light levels and relatively long integration gates. The charge resolution had a typical value of 8 photoelectrons. This resolution was not adequate for reconstructing the lateral density profile for most events, and so the ADCs were used only for slewing corrections and diagnostics. The full STACEE detector will use 1$`\mathrm{GHz}`$ waveform digitizers, which will provide greatly improved resolution and allow reconstruction of the lateral density profile. For STACEE-32, event reconstruction was based solely upon timing information.

### 3.1 Timing Reconstruction

Multi-hit TDCs measured the arrival times of PMT pulses. For each channel, the signal propagation times through the electronics were calibrated, and thus the expected arrival time of each Cherenkov pulse could be calculated, relative to the trigger time. We defined a time window with 12$`\mathrm{ns}`$ width centered around the expected hit location. A heliostat which had a TDC hit within its time window is said to have an “in-time hit.” Some channels did not have in-time hits, if those heliostats were not hit by the air shower, or if the Cherenkov pulse failed to exceed the discriminator threshold. For every in-time hit, we reconstructed the arrival time of the Cherenkov pulse at the heliostat by correcting for the transit time of light between the heliostat and the PMT, and the calibrated transit time of the PMT pulse through the electronics. ADC measurements were used to apply a slewing correction to those channels which had ADCs. The values of the slewing correction were determined from laser calibrations.

From the shape and orientation of the reconstructed timing wavefront at the heliostats, we determined the incident direction of the primary which initiated the air shower. At high energies ($`E>1\mathrm{TeV}`$), the Cherenkov-emitting core of the air shower extends through the atmosphere like a line source, and the resulting Cherenkov wavefront has a conical shape. At lower energies, however, the shower does not penetrate very far into the atmosphere, and most of the Cherenkov light is produced near the location of shower maximum. For this situation, the air shower approximates a point source of Cherenkov light, and the resulting wavefront has a more spherical shape. A ray drawn between the center of the light pool on the ground and the center of this sphere will point back towards the gamma-ray source.

The timing wavefront of each triggering event was fit to a spherical shape. Only those heliostats which had an in-time hit and a slewing-corrected time were used in the fit. The timing resolution for each channel was estimated from its discriminator rate on an event-by-event basis. Calibration runs taken with cosmic ray triggers showed that the timing resolution for channel $`i`$ was related to its rate $`R_i`$ by an empirical relation:

$$\sigma _{t,i}=\sigma _{0,i}\sqrt{1+\left(\frac{R_i}{R_{0,i}}\right)^2}.$$
(2)
The constants $`\sigma _{0,i}`$ and $`R_{0,i}`$ were determined for each channel from the calibration data. Using the estimated timing resolution, the corrected times were fit to a sphere by a least squares fitting procedure. The location of the sphere’s center was allowed to vary, but the radius of the sphere was constrained to a distance corresponding to a slant depth of 271$`\mathrm{g}\mathrm{cm}^2`$ from the top of the atmosphere. This radius corresponds to the average location of shower maximum for a 100$`\mathrm{GeV}`$ gamma ray. The fit was done three times in an iterative manner, first excluding points lying 4 standard deviations from the original fit, and then excluding any remaining points lying more than 4.5 standard deviations from the second fit.

For STACEE-32, the angular resolution was dominated by uncertainties in the determination of the center, or core, of the Cherenkov light pool. For gamma-ray air showers, the lateral density profile is quite smooth, and it is difficult to determine its center. For this analysis, we assumed the core location for all events to be at the geometric center of the array. Monte Carlo simulations indicate that STACEE-32 had an angular resolution of $`0.25^{}`$ for 200$`\mathrm{GeV}`$ gamma rays. The field of view of the instrument was only $`0.35^{}`$ (half-angle), and therefore selecting events based on the reconstructed direction would not result in any significant improvement in the significance of a detection. However, the mean direction of air showers reconstructed by STACEE-32 was within 0.1$`^{}`$of the expected pointing direction, which verified the absolute pointing of the instrument to that accuracy.

### 3.2 Hadronic Rejection

Traditional imaging atmospheric Cherenkov telescopes for gamma-ray astronomy reject cosmic ray background events based on the orientation and width of the shower image in the camera. As a non-imaging Cherenkov telescope, STACEE uses very different means for statistically identifying and removing cosmic ray events.

Cosmic ray showers differ from gamma-ray showers in the total amount of Cherenkov light they produce. A gamma-ray shower produces, on average, significantly more Cherenkov photons than a cosmic ray of the same energy (see Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT). This effect happens because the nuclear cascade of a cosmic ray air shower contains fewer Cherenkov-producing particles than the electromagnetic cascade of a gamma-ray air shower, and because the energy threshold for Cherenkov light production is higher for nucleons than for electrons. Furthermore, the difference in the Cherenkov yield between gamma-ray showers and cosmic ray showers increases as the primary energy decreases. Gamma-ray air showers are brighter and trigger the detector with much greater efficiency than cosmic ray showers.

In addition, there are large differences between the lateral density profiles of the Cherenkov light in gamma-ray and cosmic ray events. Air shower simulations show that gamma-ray showers produce smooth and uniform light pools on the ground. In contrast, cosmic ray showers produce non-uniform light distributions. In these showers, a large proportion of the light is often concentrated in a small area on the ground. STACEE-32’s trigger required 15 heliostats to be hit simultaneously across the array, and it therefore imposed a uniformity requirement. The high multiplicity requirement of the trigger thus selected the smooth, gamma-like events, and rejected most of the irregular, hadron-like events. Measured cosmic ray rates and Monte Carlo simulations suggest a trigger rejection factor of $`50`$ for cosmic rays, compared to gamma rays with the same energy spectrum. Nonetheless, there were additional hadronic cuts that were applied offline to further suppress the hadronic background.

Two hadronic cuts were developed for STACEE-32. The first cut was a trigger re-imposition cut. After in-time hits were found for each event (as described in § 3.1), the trigger multiplicity requirement was re-imposed. The in-time hit window (12$`\mathrm{ns}`$ wide) was smaller than the trigger coincidence width, and thus re-imposing the trigger in software removed some events that passed the hardware trigger. Hadronic air showers have wider, more irregular timing profiles than gamma-ray showers, and so re-imposing the trigger with a narrower time window modestly increases the signal-to-noise for gamma rays. This cut also reduces the impact of any random PMT hits that occur in coincidence with actual air showers.

The second event cut used the shape of the timing wavefront to distinguish between gamma-ray showers and cosmic ray showers. As described above, gamma-ray showers have smooth, approximately spherical, timing profiles. Cosmic ray showers tend to have more irregular, less spherical wavefronts (see Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT). The goodness-of-fit of a spherical shape to the wavefront was used to distinguish between gamma rays and cosmic rays. We selected events based on the value of the chi-squared per degree of freedom ($`\chi ^2`$/ndf) of the spherical fit used for angular reconstruction. This selection should further increase the signal-to-noise of a gamma-ray signal.

## 4 OBSERVATIONS

The STACEE-32 instrument observed the Crab extensively between 1998 November 15 and 1999 February 18. Data were taken on clear and moonless nights. We required the Crab to be within 45$`^{}`$of zenith. We took calibration data during times when Crab was low in the sky. Weather conditions and the presence of clouds were monitored regularly.

STACEE uses an ON-OFF observing strategy. Off-source runs are used to estimate the cosmic ray background level. A signal for gamma-rays shows up as an excess of ON events, compared to OFF events. On-source runs of 28 minutes length were taken with the heliostats tracking the gamma-ray source. Off-source runs were taken with the instrument tracking a point in the sky displaced by one half hour in right ascension from the gamma-ray source, but at the same declination. For each on-source run we took an off-source run of the same length. Together these two runs form an ON-OFF pair. The ON half of the pair preceded or followed the OFF half by exactly one half hour, and so both halves of the pair tracked the same trajectory in local coordinates on the sky. During these observations, STACEE-32 acquired a total of 141 ON-OFF pairs, corresponding to $`65`$ hours of on-source observing time.

### 4.1 Run Cuts

In order to do a valid background subtraction, the ON and OFF halves of each pair must be closely matched in terms of detector properties and weather conditions. A series of run cuts were imposed to remove pairs affected by changing conditions.

Run cuts were designed to remove data with poor or changing weather, detector malfunctions, or high or fluctuating PMT rates. Run pairs affected by weather conditions were identified from observing logs. Given the relative proximity of the city of Albuquerque, clouds reflect artificial lighting, and change the apparent sky brightness. Observing logs identified six ON-OFF pairs taken under partially cloudy conditions. These were removed from the data set. Another pair was removed because atmospheric haze was noted. Yet another pair coincided with the arrival of a major high-pressure front, which changed the atmospheric conditions on a short timescale. Finally, on a single night, frost developed on the heliostats, and four run pairs from that night were discarded. In all, twelve ON-OFF pairs were removed because of weather conditions.

Major detector malfunctions affected a small number of runs. Two pairs had corrupted data because of failures in the DAQ system. Runs in which multiple heliostats malfunctioned in a single subcluster were also discarded. In these runs, the affected subcluster had a relatively low efficiency for participating in the trigger because its effective coincidence level was tightened from 5 of 8 tubes to 5 of 6 tubes. Five pairs were removed on this basis. Runs in which a single heliostat malfunctioned, or in which two heliostats in different subclusters malfunctioned, were left in the data set, although the affected heliostats were removed from event reconstruction for both halves of the ON-OFF pair, as described in § 4.2 below. Heliostat malfunctions were almost always repaired the following day.

High PMT rates degrade data in several ways. As equation  shows, the timing resolution worsens with increasing rate. The rate of accidental coincidences, completely negligible under normal operating conditions, increases rapidly with rising PMT rates. Finally, high rates increase deadtime in the subcluster trigger delay electronics. Runs with high rates were identified from their subcluster trigger rates. Typical runs have subcluster rates of $`1\mathrm{kHz}`$. All runs in which the mean trigger rate for any subcluster exceeded 20$`\mathrm{kHz}`$ were removed. If the RMS variation of a subcluster’s rate within the run exceeded 3$`\mathrm{kHz}`$, the run was also removed. These stringent cuts ensure that the accidental trigger rate remained $`<0.01\mathrm{Hz}`$, and that the deadtime in the trigger delay of each subcluster was $`<1\%`$. Twenty ON-OFF pairs were removed by these cuts.

Finally, the data were scanned for runs with abnormally low event trigger rates. A single anomalous run was identified with a trigger rate about eight standard deviations below that of similar runs. The low rate is indicative of a detector malfunction, the cause of which is under investigation. The pair containing this anomalous run was removed from the final data set.

Run cuts were applied in a blind fashion, before the significance of any possible gamma-ray excess was determined, so as not to bias the result. After application of run cuts, 101 ON-OFF pairs remain in the data set. Table 1 summarizes the various run cuts.

### 4.2 PMT Cuts

It was not uncommon for a single heliostat or electronics channel to malfunction. Removing the entire run in this circumstance is rather draconian, and unnecessarily reduces the data set. Instead, bad channels were identified on a run-by-run basis. Offending channels were then removed from the offline analysis in both the ON and OFF halves of the run pair, and the data were analyzed as if the PMT were simply turned off. The most common problems were heliostat tracking errors, which could be identified from the heliostat log files. Removing the bad channels from both halves of a pair ensures that the ON and OFF pairs are balanced in terms of detector response.

## 5 RESULTS

After run cuts, the final STACEE-32 data set consists of 101 ON-OFF pairs, with a total on-source observing time of 155,335$`\mathrm{s}`$ ($`43`$ hours). Each run was processed separately, and event reconstruction proceeded as described in § 3. For each event, the GPS event time, fitted shower direction, and $`\chi ^2`$/ndf of the spherical fit were found. The livetime fraction for each run was calculated, and the event totals for each run were corrected for deadtime by dividing by the livetime fraction (typically between 88% and 92%).

A gamma-ray signal shows up as an excess of events in the ON runs, compared to the OFF runs. As described in § 3.2, two kinds of event cuts were used to enhance any potential gamma-ray signal. The first cut was to re-impose the trigger multiplicity in software. The second cut was based on the $`\chi ^2`$/ndf of a spherical fit to the wavefront’s shape. Monte Carlo simulations suggest that a cut of $`\chi ^2`$/ndf$`<1`$ will maximize the signal-to-noise.

### 5.1 Unpulsed Emission

First we performed a search for unpulsed emission, presumed to originate in the nebula of the Crab. The number of events in the on-source and off-source runs were tallied, with and without event cuts. Table 2 summarizes the results. The raw data show a strong excess of 4860 events, before event cuts, against a background of $`420,000`$ events. The Li-Ma statistical significance (Li & Ma, 1983) of this excess is $`+5.27\sigma `$. Of the 101 ON-OFF pairs, 65 pairs show an excess in their raw rates, while just 36 show a deficit.

Next, the trigger multiplicity was re-imposed. The total number of excess events after trigger re-imposition was 4551, with a background of $`350,000`$, for a total significance of $`+5.44\sigma `$. Seventy of the 101 pairs showed an excess of events from the source.

Finally, a cut on $`\chi ^2`$/ndf $`<1`$ was applied on top of the trigger re-imposition cut. The remaining number of excess events was 4062. The background was reduced by a factor of approximately two, to $`180,000`$ events, and the total significance of the signal increased markedly to $`+6.75\sigma `$. Of the 101 pairs, 77 now show an excess of events from the source region, while just 24 show a deficit.

The observed excess is statistically strong, being present at greater than the five standard deviation level in the raw trigger rates alone. Because STACEE is a new detector using a novel technique, there is the question of whether this excess is actually due to a gamma-ray signal, or whether it could be due to some unforeseen systematic effect that is present in spite of our efforts to closely match ON and OFF halves of each pair. This question will ultimately be decided by confirmation of the result by other experiments. Detection of other sources with this technique should provide further confirmation of the method. Nonetheless, there are consistency checks which can be performed to strengthen the conclusion that STACEE-32 sees gamma rays from the Crab.

The first and most powerful check is that the significance of the excess increases as selection cuts are applied to the data. Both the trigger re-imposition cut, and especially the cut on the $`\chi ^2`$/ndf, increase the significance of the signal by amounts in agreement with expectations from simulation. In short, the observed excess behaves just as a gamma-ray signal should.

Secondly, one can look at the distribution of pairwise significances. The pairwise significance is defined as the observed excess or deficit, in standard deviations, for each ON-OFF pair. For a steady source with a constant flux, the distribution of pairwise significances should be normally distributed, with unit width and a shifted mean. As is seen in Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT, the distribution of pairwise significances for our data does have the expected form. This fact further supports the interpretation of the excess as a gamma-ray signal.

In Table 3, the data have been broken down by month. Since fewer runs were taken in January and February of 1999, these months have been combined. An excess with a significance of $`4\sigma `$ is seen in each era. Also noteworthy is that the background rate (OFF events per unit time) decreased by 20% between the November data and the December-February data. We interpret this decrease as an increase in the energy threshold of the experiment. Closer examination of the data shows that the rates were steady within the month of November, and within December as well. Whatever change happened to the detector must therefore have occurred during the intervening full Moon period during which STACEE did not operate. Although the reason for this increase in threshold is not known with certainty, we suspect that the angular pointing of the heliostats drifted slightly out of alignment, decreasing the optical throughput.

### 5.2 Pulsed Emission Search

Having established the presence of a gamma-ray signal in our data, we then carried out a search for pulsed emission in phase with the Crab Pulsar’s emission. A GPS clock recorded the arrival time of every event with an accuracy of $`1\mu \mathrm{s}`$. The arrival time of each event at Earth was corrected to barycentric dynamical time (TDB) using the JPL Planetary and Lunar Ephemerides DEC-200 package (Standish, 1997). The corresponding radio phases were obtained from the Jodrell Bank ephemeris (Lyne, Pritchard, & Roberts, 1999). The phase values were interpolated between monthly epochs by expanding the phase in a Taylor series, and by requiring continuity of the phase, period, and time derivatives of the period at the midway point between the two months. The barycentering and ephemeris calculations were checked by applying them to optical pulsar data recorded by the Whipple gamma-ray collaboration. For these data, we verified that the optical pulse was extracted properly.

The arrival times of on-source events were then folded with the pulsar phase to produce a phase histogram for the STACEE-32 data. Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT shows the phase histogram for on-source events passing all event cuts. No obviously significant structure is seen. We have applied the H-test to test the uniformity of the phase histogram (de Jager, Swanepoel, & Raubenheimer, 1989; de Jager, 1994). The test yields a value for the $`h`$ parameter of 5.37, which has a chance probability of 0.12. We therefore conclude that STACEE-32’s Crab phase histogram is consistent with being uniform in phase.

Concluding that our data set contains no strong evidence of pulsed emission, we then set an upper limit on the pulsed fraction of the total observed excess. Because the pulse profile of the pulsed emission at STACEE’s energies is unknown, we make an assumption about the phase intervals in which the gamma-ray emission would occur. We assume that the pulsed emission at STACEE energies would occur in phase with the pulsed emission seen by EGRET (Nolan et al., 1993; Ramanamurthy et al., 1995; Fierro et al., 1998). We defined an “on-pulse” region for the phase histogram, which includes an interval centered about the main radio pulse at 0.94-0.04 of the phase histogram, and a second interval about the intrapulse from 0.32-0.43. Using the method of Helene (1983), we derive an upper limit on the pulsed fraction of $`<5.5\%`$ of the observed excess, at the 90% confidence level. Table 4 contains the event totals for the on-pulse and off-pulse regions.

### 5.3 Energy and Integral Flux Result

The energy threshold and integral flux for these data must be determined from careful calibration and simulations of all parts of the detector. STACEE-32 used three simulation packages to model the response of the instrument to gamma rays. The well-established MOCCA code was used to model extensive air showers and Cherenkov light production (Hillas, 1985). A complete ray-tracing program followed the paths of photons through the optics. Finally, a detailed electronics simulation which made use of a library of digitized PMT pulses was used to model the performance of the PMTs and electronics. All parts of the simulations have been verified against calibration data.

For the purpose of determining an energy threshold and integral flux, only data taken in November 1998 were used. As explained in § 5.1, these data had the lowest energy threshold, as measured by the trigger rates. Also, the detector was carefully aligned during this month. During the subsequent months, it is believed that small degradations in the alignment increased the energy threshold, but we do not have enough information to track the effect reliably. Using only data from November 1998 therefore provides the lowest possible energy threshold and minimizes possible systematics associated with optical alignment. The increased statistical uncertainty from using a subset of the data as opposed to the entire data set is still smaller than the systematic uncertainty on the flux. Likewise, we have used the total rates before event cuts in determining the energy threshold and integral flux, since the efficiencies of these cuts depend on the simulation, and so would introduce additional systematic uncertainties if included.

Monte Carlo simulations were used to determine the sensitivity of STACEE-32 to gamma rays. The sensitivity was calculated for gamma rays of various energies coming from multiple positions along the Crab’s trajectory across the sky. The results were expressed as an effective area curve, which is the gamma-ray collection area of STACEE-32 as a function of the gamma-ray energy. The average effective area curve for the November 1998 data was calculated by weighting each position on the sky by its exposure in the data set. Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT shows the average effective area as a function of energy. The sensitivity of STACEE-32 starts at an energy below 100$`\mathrm{GeV}`$, and quickly rises with energy. The effective area approaches a plateau value of $`28,000\mathrm{m}^2`$ above 1$`\mathrm{TeV}`$.

To determine the energy threshold and integral flux, we assume that the Crab’s energy spectrum follows a differential power law of the form

$$\frac{dN}{dE}=CE^{2.4}.$$
(3)
The power law index of $`2.4`$ is estimated from measured Crab spectra, interpolated to STACEE-32’s energy range (Hillas et al., 1998). Convolving this spectrum with the effective area curve yields the differential trigger rate per unit energy. The peak rate occurs at an energy of 190$`\mathrm{GeV}`$. This is our estimate for the so-called “spectral threshold energy.” The energy threshold changes by less than $`\pm 10\mathrm{GeV}`$ if the assumed spectral index is varied by $`\pm 0.2`$.

Estimates of the uncertainty in each element of the simulations (e.g. uncertainties in measured mirror reflectivities, PMT gains, etc.) yield a total systematic uncertainty on the energy of 32% ($`\pm 60\mathrm{GeV}`$). No single factor dominates the systematic uncertainty.

For the November data, a total of 183,501 on-source events were seen in the raw data (without event cuts), versus 181,349 off-source events, in 56,056$`\mathrm{s}`$ of observing time. The excess rate is therefore $`0.038\pm 0.011\mathrm{Hz}`$. By integrating the assumed differential flux with the effective area curve, and equating this to the observed rate, we determine the integral flux. We find that the integral flux of gamma rays from the Crab Nebula above our energy threshold of $`E_{th}=190\pm 60\mathrm{GeV}`$ is

$$I(E>E_{th})=(2.2\pm 0.6\pm 0.2)\times 10^{10}\mathrm{photons}\mathrm{cm}^2\mathrm{s}^1.$$
(4)
Here the first error is statistical, and the second error is the systematic error on the flux itself, not including the effects of uncertainty in the energy threshold. The dominant systematic error is the uncertainty in the energy threshold $`E_{th}`$, which does not change the flux value itself, but which does change the energy at which that flux is reported.

## 6 DISCUSSION

Figure HIGH ENERGY GAMMA-RAY OBSERVATIONS  OF THE CRAB NEBULA AND PULSAR WITH THE  SOLAR TOWER ATMOSPHERIC CHERENKOV EFFECT EXPERIMENT shows STACEE-32’s integral flux value. Also shown for comparison are measurements from the CAT (Barrau et al., 1997) and Whipple (Hillas et al., 1998) experiments. Our measured flux is consistent with an extension to lower energies of the power law spectrum seen by Whipple. Inverse Compton models generally predict a hardening of the spectral slope as the energy decreases. Given the systematic uncertainties, our measurement is consistent with these expectations.

No evidence of pulsed emission from the Crab Pulsar is seen, and an upper limit on the pulsed fraction is derived at $`<5.5\%`$ (90% CL) of the observed signal. At present, the energy threshold of this limit, although lower than any previous limit, is not low enough to differentiate between Polar Cap and Outer Gap models. As the energy threshold of STACEE is lowered further, however, future limits should constrain the theoretical models.

This detection is the first for the STACEE instrument. The CELESTE collaboration, which is operating its own solar heliostat experiment at the Thémis site in France, has reported a preliminary detection of gamma rays from the Crab using a similar technique. Their data analysis is in progress (Smith et al., 1998). Although STACEE-32 was a prototype instrument, it achieved an unprecedentedly low energy threshold. The full STACEE detector will be completed in the year 2000, and will feature twice as many heliostats, improved trigger electronics, and 1$`\mathrm{GHz}`$ sampling waveform digitizers. Based upon the STACEE-32 results, we expect that the full STACEE instrument will obtain its design goal of an energy threshold of $`50\mathrm{GeV}`$.

In the future, solar heliostat experiments such as STACEE, CELESTE (Giebels et al., 1998), and Solar Two (Zweerink et al., 1999) will provide complete spectral coverage at energies between $`50`$ and $`500\mathrm{GeV}`$ for the Northern Hemisphere. These sorts of low threshold experiments will provide continuity between satellite measurements and ground-based detectors at TeV energies, and will naturally complement future satellite experiments such as GLAST.

We are grateful to the staff at the NSTTF for their excellent support. We thank the Physics Division of Los Alamos National Laboratory for loans of electronic equipment. We also thank the SNO collaboration for providing us with acrylic for light concentrators, and the machine shop staffs at Chicago and McGill for their assistance. The Whipple collaboration kindly provided us with optical pulsar data for testing our barycentering analysis. Many thanks to Tom Brennan, Katie Burns, Jaci Conrad, Jim Hinton, William Loh, Anthony Miceli, Gora Mohanty, Alex Montgomery, Heather Ueunten, and François Vincent. This work was supported in part by the National Science Foundation, the Natural Sciences and Engineering Research Council, FCAR (Fonds pour la Formation de Chercheurs et l’Aide à la Recherche), the Research Corporation, and the California Space Institute. CEC is a Cottrell Scholar of Research Corporation.


text:
# Memory in aged granular media

## Abstract

Stimulated by recent experimental results, we simulate “temperature”-cycling experiments in a model for the compaction of granular media. We report on the existence of two types of memory effects: short-term dependence on the history of the sample, and long-term memory for highly compact (aged) systems. A natural interpretation of these results is provided by the analysis of the density heterogeneities.

The study of the rheology of granular media represents nowadays an important chapter in the more general context of the so-called soft condensed matter . The general questions one is interested in concern the response behaviour of a generic granular medium submitted to an external perturbation. In this spirit it has been recognized that the response of such systems depends in a non-trivial way on the mechanical properties of the grains composing the medium, on the boundary conditions, on the driving procedure (the way to inject energy to perturb the system) and, last but not the least, on the past history of the system, e.g. on the procedures the system has undergone until the moment we study it.

In this context the similarities between granular media and the phenomenology of glassy systems have been exploited. One of the main results has been to recognize that in a compaction procedure, where a certain system increases its density if shaken or tapped, a granular medium displays ageing  and its properties, e.g. the two-times correlation function, depend on the time where one performs the measurement, i.e. on the age of the system.

The analogy with glassy systems has mainly allowed to export tools and techniques from one field to the other. Typical experiments on glassy systems consists in monitoring the effect of changes imposed to the external temperature on the macroscopic observables . Translated in the framework of granular matter, i.e. of a non-thermal system, this means to analyse the effect of changes in the driving procedure, e.g. in the tapping amplitude $`\mathrm{\Gamma }`$, which, in a typical compaction experiment, represents the control-parameter.

In a previous work , we have presented a detailed study, in the framework of a microscopic model for granular media, focusing on the response to small perturbations in the driving rate. We have pointed out the importance of the interplay between the response properties and the spatial structures that spontaneously emerge as a consequence of the dynamics imposed to the system. In particular we have studied the relevance of the densities heterogeneities in the process of granular compaction, and shown that the density is not the only relevant information about the system. Configurations with the same density but reached with different dynamical procedure show completely different rheological properties. This was nicely confirmed by a set of recent experiments: three systems prepared at the same density but in different ways display different behaviours if the same tapping acceleration is applied to them. This feature concerns the coding of the system history in the microscopic configurations.

In the same set of experiments, similar in spirit to the temperature-cycling experiments in spin-glasses the authors of have put in evidence the existence of memory effects in the compaction of granular media. In particular they performed the study of the effect of relatively large and abrupt changes in the tapping acceleration. This is an efficient way to obtain reliable experimental data since the effect of small changes in the tapping acceleration is difficult to detect in real experiments. What was shown in  is that, during the compaction, the response to an abrupt change in the tapping acceleration $`\mathrm{\Gamma }`$ is opposite (at least at short times) to what could be expected from the long-time behaviour of the compaction: while at constant $`\mathrm{\Gamma }`$ a higher compaction rate is obtained for larger $`\mathrm{\Gamma }`$, a sudden increase in $`\mathrm{\Gamma }`$ leads to a decompaction, and a sudden decrease gives rise to a temporary increase in the compaction rate. This behaviour is however transient, and the usual compaction rate is recovered after a while: the authors therefore call this effect “short-term memory”. The same qualitative behaviour has been obtained also in experiments of compaction under shear .

In order to push forward this kind of studies we have extended our work  to study the response to large changes in the driving. In particular we explore the behaviour of the Tetris Model  for abrupt changes in the driving rates in various time regimes, that were not considered in the preliminary, time-limited, experiments . We also explore the return to the initial tapping acceleration and the presence or absence of memory, in a way similar to the temperature-cycling experiments in spin-glasses. We perform these temperature-cycling experiments in non-equilibrium (aging) as well as in quasi-stationary (time-translation invariant) situations. We hope in this way to stimulate new and more precise experiments along the same lines.

The model we consider, the so-called Tetris Model, consists of particles on a tilted (two-dimensional) square lattice, which has lateral periodic boundary conditions 111the use of closed lateral boundary conditions do not change the results and a closed boundary at the bottom. With respect to its original definition we consider a random version where each particle can be schematized in general as a cross with $`4`$ arms (in general the number of arms is equals to the coordination number of the lattice) of different lengths, denoted by $`l_{NE}`$, $`l_{NW}`$, $`l_{SE}`$, $`l_{SW}`$, chosen in a random way (see  for details).

Gravity and shaking are implemented as follows: the particles diffuse on the lattice, with geometrical constraints, and probability $`p_{up}`$ to move upwards (with $`0<p_{up}<1`$), and a probability $`p_{down}=1p_{up}`$ to move downwards. The quantity $`x=p_{up}/p_{down}`$ can be related to the adimensional acceleration $`\mathrm{\Gamma }`$ used in compaction experiments through the relation $`\mathrm{\Gamma }1/log(1/\sqrt{x})`$. The sizes (width $`\times `$ height) used ranged from $`200\times 60`$ to $`20,40,60\times 200`$, in order to test various aspect ratios.

The procedure used in the experiments and in our simulations is inspired by classic experiments in spin-glasses. The system is initialized by random addition of particles, and let evolve with a constant external forcing $`x_1`$ up to a certain time $`t_w`$. Two copies of the system are then made: one copy is kept evolving with $`x_1`$ at all times for reference, while the second copy evolves with a different forcing $`x_2`$ for a time $`\mathrm{\Delta }t`$, after which the forcing is set back to $`x_1`$. We monitor the mean height, the bulk density (defined as the average density in the lower $`50`$ of the system), and the density profiles of both copies during the whole evolution for various values of the difference $`x_1x_2`$, $`t_w`$ and $`\mathrm{\Delta }t`$. We recall that the density profiles are defined averaging the density along the horizontal direction, giving thus indications on the heterogeneities along the vertical direction (see figure 1), and allowing to roughly distinguish two regions, the bulk and the interface. Before analyzing the different regimes, let us recall a few general features of the compaction at constant forcing. At short times, the density is an increasing function of the forcing while at long times, there exists an optimal value of $`x`$ that allows to obtain the maximal density in the fraction of the system considered for the measurement. This phenomenology has been shown to be linked to the presence of heterogeneities in the system .

First case: $`𝐱_\mathrm{𝟏}>𝐱_\mathrm{𝟐}`$. When, at $`t_w`$, $`x`$ is suddenly lowered, the compaction rate of the perturbed system first increases. This behaviour, shown in figure 2 and 3, is in agreement with the experimental results, and it is opposed to the results one would expect looking at the compaction data at constant forcing. This shows that the system has some memory of its history at $`t_w`$. After a transient, however, this memory is lost: compaction slows down, and the rate of compaction crosses over to the one observed at constant forcing: the curves of the reference and perturbed systems therefore cross.

While the transient, corresponding to the short-term memory, is observed for all values of the parameters investigated, an important difference appears for different $`t_w`$ at long times: if $`t_w`$ is small, the system is not very compact, and is able to evolve a lot during $`\mathrm{\Delta }t`$. The system does not display memory for times larger than $`t_w+\mathrm{\Delta }t`$. On the other hand, if the system is already quite compact at $`t_w`$ (i.e. if it is aged enough), the system evolves very slowly during $`\mathrm{\Delta }t`$, and for $`tt_w+\mathrm{\Delta }t`$ the compaction curves can be translated and superimposed to the reference one, as shown in the inset of figure 3. During $`\mathrm{\Delta }t`$ the system is able to keep memory of its state at $`t_w`$. This shows that a sufficiently aged system displays a long term memory, in addition to the short-term memory existing at all times.

The interpretation of these results is quite straightforward using the results of  and looking at the density profiles along the vertical direction: when $`x`$ is abruptly lowered, the first effect is that the particles tend to go down, and the interface becomes steeper and more compact. Therefore the density first increases with respect to the unperturbed case. At larger times however, the evolution is slowed down by the creation of a dense layer at the interface, which blocks the bulk rearrangements needed for the compaction . After $`t_w+\mathrm{\Delta }t`$, the increase in the forcing allows to suppress the dense layer, and the compaction can become again fast. Moreover, if $`t_w`$ is large enough (typically $`10^3`$ or $`10^4`$ MC steps), the bulk of the system is already quite compact, and therefore the smaller value of the forcing during $`\mathrm{\Delta }t`$ leads to a compaction of the interface but the bulk almost does not evolve. At $`t_w+\mathrm{\Delta }t`$, the forcing is again increased: the relaxation of the interface being fast, this leads the system back to its state at $`t_w`$. The inset of Figure 3 illustrates this behaviour: the curves of compaction after $`\mathrm{\Delta }t`$, shifted backwards by $`\mathrm{\Delta }t`$, superimpose with the reference curves. This is not the case for lower values of $`t_w`$, for which the bulk is not yet highly compact at $`t_w`$, and evolves a lot even with a small driving $`x_2`$.

This memory effect, similar to what can be obtained in spin-glasses when the temperature is first lowered then increased back to its previous value, was unfortunately not investigated in : it can indeed occur only at long times. It would certainly be very interesting to have experimental checks on this aspect.

Second Case: $`𝐱_\mathrm{𝟏}<𝐱_\mathrm{𝟐}`$. When $`x_1<x_2`$ we observe a short-term memory effect, as in the previous case. First, as the forcing is increased, one observes a decompaction; later on, the fact that $`x`$ is larger prevails and the compaction proceeds faster, at the normal rate for constant $`x=x_2`$. The memory of the history up to $`t_w`$ is lost after a transient (see figure 4). These features are again similar to the experimental results. Moreover, at long times no memory is kept, since the whole system evolves a lot during $`\mathrm{\Delta }t`$.

Once again, the study of the density profile allows for an interpretation of these results. As the tapping intensity is increased, the first effect is a decompaction, especially at the interface. The fact that the interface is less compact then allows for a much better compaction of the bulk. Note that this behaviour is closely related to the existence of a response function that is positive at short times and negative at longer times . At $`t=t_w+\mathrm{\Delta }t`$, the bulk has been deeply modified, so the system cannot have any memory of its configuration at $`t_w`$.

As a first partial conclusion we see that the interplay between the heterogeneities of the system and the influence of the forcing on different zones of the sample  leads to two quite different memory effects. On one hand the short-time memory is always present and it occurs in a way which is counterintuitive with respect to the effect expected on the basis of the behaviour at constant $`x`$. On the other hand the long-time memory only exists for a cycle when $`x_2<x_1`$ and for long times, i.e. for aged systems.

Perturbation of a stationary state The numerical experiments mentioned above have been performed during the compaction process, i.e. while the system ages. Let us now focus instead on the behaviour obtained after a slow decrease of the vibration intensity, starting from large forcing: this is the usual way to obtain a very compact system  which, under weak driving $`x_1`$, is now in a quasi-stationary state, i.e. a state where time-translation invariance holds. This state can now be perturbed by applying a stronger forcing $`x_2`$ for a time interval $`\mathrm{\Delta }t`$, and the relaxation can be studied when the forcing is set back to the original value $`x_1`$. In , such experiments were performed for small $`\mathrm{\Delta }t`$, showing non-exponential relaxations, slower the larger is $`\mathrm{\Delta }t`$. We are able here to explore different regimes, varying $`\mathrm{\Delta }t`$ over four orders of magnitude.

One observes a relaxation occurring in two steps (see figure 5). First one has a rapid decrease (during typically 10 time steps) corresponding to the relaxation of the topmost layers of the interface; this can be checked by looking carefully at the density profiles during the relaxation. Later on a much slower, almost logarithmic, decay occurs corresponding then to the re-compaction of layers under the interface. As $`\mathrm{\Delta }t`$ grows, deeper and deeper layers get involved during the decompaction process; as $`x`$ is lowered again, the first relaxation, which compactifies the interface, makes it more difficult for the lower layers to get as compact as before: the relaxation therefore slows down.

Although the very dense medium, if shaken at constant $`x`$, is in a stationary state, the perturbation at larger $`x`$ leads to same mechanisms that are also at work during the aging. Therefore, it is not so surprising to find a slow re-compaction process after the perturbation.

The precise form of the relaxation can be expected to be model dependent, at least at short times: for example, in the case of “tapping”, the first, rapid part of the relaxation does not exist. Therefore, we do not investigate it in details, and underline only the most relevant features: even though the system seems in “equilibrium” (as shown in , the dense system is stationary and does not display aging), the relaxation is non-exponential and very slow.

In summary, we have shown that the Random Tetris Model reproduces the complex phenomenology appearing in real experiments for large and abrupt changes of the driving: at small times after the change, the reaction of the system is opposite to what could be expected looking at the constant $`x`$ compaction curves: if $`x`$ is decreased, the system compactifies better, while it gets decompactified for an increase in $`x`$. At larger times, the compaction rate crosses over to the one at constant $`x`$: this effect was therefore called ’short-term memory’.

Moreover, we have investigated the return to the previous value of the driving, in various time regimes, and shown that an increase in $`x`$ erases quickly the memory of the previous configurations while, at long times, a decrease in $`x`$ can “quench” the system until the previous value of $`x`$ is restored and the compaction starts again from a configuration statistically equivalent to the one perturbed at $`t=t_w`$.

It is noticeable that such a behaviour shares some similarities with the phenomena of rejuvenation and memory observed in the context of spin glasses . This feature pushes forward the similarities between (non-thermal) granular media and (thermal) glasses and spin-glasses. For example, note that a trap model , widely exploited for the understanding of the phase space of spin glasses, has recently been used to describe granular media , as well as a random-walk model . These approaches also reproduce qualitatively the experimental results of ; however they describe the behaviour of a system in phase space. It seems to us that the study of the heterogeneities (and therefore of real space models) could be crucial for the understanding of the basic mechanisms underlying the dynamics of these systems. In this respect granular media seem more tractable than spin glasses and detailed experiments about the local rearrangements during compaction could provide valuable information for the understanding of the basic mechanisms at work in dense granular media.