format: text

text:
# Superconductivity and Quantum Criticality in CeCoIn5

## Abstract

Electrical resistivity measurements on a single crystal of the heavy-fermion superconductor CeCoIn5 at pressures to 4.2 GPa reveal a strong crossover in transport properties near $`P^{}1.6`$ GPa, where $`T_c`$ is a maximum. The temperature-pressure phase diagram constructed from these data provides a natural connection to cuprate physics, including the possible existence of a pseudogap.

The relationship between unconventional superconductivity and quantum criticality is emerging as an important issue in strongly correlated materials, including cuprates, organics and heavy-fermion intermetallics orenstein00 ; chubukov02 . In each of these, either Nature or an externally imposed parameter, such as pressure or chemical substitutions, tunes some magnetic order to a $`T=0`$ transition where quantum fluctuations introduce new excitations that control thermodynamic and transport properties over a broad range of temperature-parameter phase space. A frequent, but controversial anderson02 , assumption is that these same excitations may mediate Cooper pairing in the anisotropic superconductivity that appears in the vicinity of a quantum-critical point (QCP). Indeed, the effective potential generated by proximity to a quantum-critical spin-density wave coleman01a  also can lead montoux91  to anisotropic superconductivity.

Though the issue of quantum criticality and superconductivity came to prominence in the context of the cuprates (see, eg varma01 ), there have been several recent examples of this interplay in heavy-fermion compounds mathur98 . Like the cuprates, unconventional superconductivity in these Ce-based heavy-fermion materials develops out of a distinctly non-Fermi-liquid normal state that evolves in proximity to a continuous $`T=0`$ antiferromagnetic transition; but unlike the cuprates, this state can be accessed cleanly by applied hydrostatic pressure and without introducing extrinsic disorder associated with chemical substitutions. In addition to chemical inhomogeneity in the cuprates, a further ’complication’ is the existence of a pseudogap state above $`T_c`$ for dopings less than optimal timusk99 . Whether evolution of the pseudogap with doping produces another QCP near the optimal $`T_c`$ in the cuprates is an open question as is the origin of the pseudogap varma99 . Thus far, a pseudogap has not been detected in canonical quantum-critical heavy-fermion systems, such as CeIn3 or CePd2Si2 under pressure, but this may not be too surprising. In these cases, superconductivity appears only at very low temperatures, less than 0.4 K, and over a narrow window of pressures $`(0.8\mathrm{GPa})`$ centered around the pressure where $`T_N`$ extrapolates to zero. In analogy to the cuprates, we would expect a pseudogap in these Ce compounds to exist only over a narrow temperature window, $`0.10.2`$ K above $`T_c`$, and in a fraction of the narrow pressure window, which certainly would make detection of a pseudogap difficult.

CeCoIn5 offers the possibility of making a connection between other Ce-based heavy-fermion materials and the cuprates. Bulk, unconventional superconductivity is present in CeCoIn5 at atmospheric pressure petrovic01b ; movshovich01 ; izawa01 ; kohori01  and develops out of a heavy-fermion normal state in which the resistivity is approximately linear in temperature. Application of a magnetic field sufficient to quench superconductivity reveals a specific heat diverging as $`T\mathrm{ln}T`$ and associated low-temperature entropy consistent with its huge zero-field specific heat jump at $`T_c`$ $`(\mathrm{\Delta }C/\gamma T_c=4.5)`$. Together, these properties suggest that CeCoIn5 may be near an antiferromagnetic quantum-critical point at $`P=0`$. Its isostructural, antiferromagnetic relative CeRhIn5 becomes a pressure-induced unconventional superconductor near 1.6 GPa, with a $`T_c=2.1`$ K close to that of CeCoIn5 at $`P=0`$ hegger00 ; fisher01 . This and other similarities between CeCoIn5 at $`P=0`$ and CeRhIn5 at $`P=1.6`$ GPa reinforce speculation nicklas01 ; shishido02  that the nearby antiferromagnetic QCP in CeCoIn5 may be at an inaccessible slightly negative pressure. The 1.7% smaller cell volume of CeCoIn5 compared to that of CeRhIn5 is consistent with this view. We have studied the effect of pressure on the electrical resistivity and superconductivity of CeCoIn5 and find striking correlations between them that are reminiscent of cuprate behaviors.

Four-probe AC resistivity measurements, with current flowing in the tetragonal basal plane, were made on a single crystal of CeCoIn5 grown from excess In flux. Hydrostatic pressures to 5 GPa were generated in a toroidal anvil cell khvostantsev98  in which a boron-epoxy gasket surrounds a teflon capsule filled with a glycerol-water mixture (3:2 volume ratio) that served as the pressure transmitting fluid. Pressure inside the capsule was determined at room temperature and at low temperatures from the pressure-dependent electrical resistivity and $`T_c`$ of Pb, respectively eiling81 . The width of superconducting transition of Pb did not exceed 15 mK, indicating good hydrostatic conditions and providing an estimate of the pressure-measurement uncertainty, $`\pm 0.04`$ GPa.

The response of $`\rho (T)`$ to pressure over a broad temperature scale nicklas01 ; shishido02  is typical of many Ce-based heavy-fermion compounds and reflects a systematic increase in the energy scale of spin fluctuations. Most interesting is the low-temperature response (Fig. 1) where $`T_c`$ and the resistivity just above $`T_c`$ $`(\rho (T_c^+))`$ decrease rapidly with pressures $`P>1.35`$ GPa. We have analyzed the low temperature resistivity, plotted over a broader temperature scale in the inset of Fig. 1, by either fitting to a form $`\rho (T)=\rho _0+AT^n`$, with $`\rho _0`$, $`A`$ and $`n`$ fitting parameters, or from plots of $`(\rho (T)\rho _0)/T^n`$. A drawback to the former approach is that parameters can be sensitive to the temperature range chosen for fitting. The more sensitive and revealing approach is the latter in which $`\rho _0`$ and $`n`$ are adjusted to give the best horizontal line. This procedure requires one less parameter, $`A`$, which is provided directly by the magnitude of $`(\rho (T)\rho _0)/T^n`$. Representative examples of this approach are plotted in Fig. 2. Parameters obtained from these plots are well defined and agree, typically to 10% or better, with parameter values extracted from straightforward fits. As indicated by arrows on the lowest pressure curve in Fig. 2, $`(\rho (T)\rho _0)/T^n`$ deviates from the horizontal trend at a temperature $`T_{pg}`$ well above the onset of superconductivity. This departure implies either a decrease in the scattering rate or increase in carrier density. Resistivity measurements at $`P=0`$ and in a magnetic field (6 T) greater than $`H_{c2}(0)`$ find movshovich02  a similar departure of $`(\rho (T)\rho _0)/T^n`$ from a $`T`$-independent curve at $`T_{pg}`$=3 K. Thus, $`T_{pg}`$ does not originate from superconducting fluctuations or from trace amounts of free indium. We also note that this behavior was present in previous measurements on other single crystals CeCoIn5 but was missed in inspections of $`\rho `$ versus $`T`$ curves nicklas01 . We will return to a discussion of $`T_{pg}`$ later.

Values of $`\rho _0`$, $`n`$ and $`A`$ obtained from these plots are summarized in Figs. 3a and b. There is an unmistakable crossover in the pressure dependence of $`\rho _0`$ and $`n`$ at $`P^{}1.6`$ GPa. Below $`P^{}`$, $`n=1.0\pm 0.1`$ for $`T_{pg}T10`$ K. This value of $`n`$ is expected for a 2-dimensional, antiferromagnetic quantum-critical system sachdev99  and is not too surprising given the layered crystal structure moshopoulou01  and anisotropy in electronic states of CeCoIn5 petrovic01b ; settai01 . At pressures greater than $`P^{}`$, $`n`$ rapidly approaches the Fermi-liquid value of 2.0, which holds from just above $`T_c`$ to $`T_{FL}`$ (see Fig. 2). Though not shown in either Fig. 2 or 3a, at these higher pressures $`n`$ assumes a value of $`1.5\pm 0.1`$, characteristic of a 3-D antiferromagnetic QCP, in a temperature interval from slightly greater than $`T_{FL}`$ to as high as $``$ 60 K. For $`PP^{}`$, $`\rho _0`$ decreases reversibly by an order of magnitude to a very small value of about 0.2 $`\mu \mathrm{\Omega }\mathrm{cm}`$. Clearly, pressure does not remove impurities from the sample; the large decrease in $`\rho _0`$ must be due to a pressure-induced change in inelastic scattering processes. Theories of electronic transport in quantum-critical systems show rosch99 ; miyake01  that impurity scattering can be strongly enhanced by quantum-critical fluctuations. This sensitivity provides a natural explanation for a decrease in scattering with increasing pressure, if there is a crossover near $`P^{}`$ from quantum-critical $`(P<P^{})`$ to a Fermi-liquid-like state for $`(P>P^{})`$. Such a crossover is suggested by the pressure variation of $`n`$.

Experimental values of the $`T^n`$-coefficient of resistivity, plotted in Fig. 3b, also decrease strongly with increasing pressure for $`PP^{}`$. For comparison, we show the variation of the electronic specific heat divided by temperature $`(C/T\gamma )`$ at 3 K for pressures less than 1.5 GPa sparn01  and, at higher pressures, we plot values of $`\gamma `$ inferred from the empirical relationship $`A=1\times 10^5\gamma ^2`$, where $`A`$ is the $`T^2`$-coefficient of resistivity, followed by several heavy-fermion compounds kadowaki86 . Though a $`T^2`$ coefficient proportional to $`\gamma ^2`$ is expected for a Landau Fermi-liquid varma01 , it is not obvious a priori why $`A(P)`$ approximately tracks the directly measured $`\gamma `$ at low pressures, but this seems to be the case. In any event, both $`A(P)`$ and $`\gamma (P)`$ exhibit qualitatively different magnitudes and functional dependencies above and below $`P^{}`$. The dramatic differences in $`\gamma (P)`$ at low and high pressures suggest that either the empirical relationship is invalid and underestimates $`\gamma `$ in the high-pressure regime by at least two orders of magnitude or there is a very rapid crossover in the density of low-energy excitations in the vicinity of $`P^{}`$. In the Fermi-liquid-like state above $`P^{}`$, we would expect the temperature scale $`T_{FL}`$ (see Fig. 2) to be proportional to $`A^{0.5}`$ instead of the approximate $`T_{FL}\sqrt{A}`$ found experimentally. We do not understand this discrepancy. Direct measurements of $`C/T`$ at pressures above $`P^{}`$ would be valuable.

Finally, we note that the resistive transition to the superconducting state is equally sensitive to the crossover at $`P^{}`$. As shown in Fig. 3c, the relative transition width $`\mathrm{\Delta }T_c/T_c`$ passes through a pronounced minimum near $`P^{}`$. To some extent this behavior is accentuated by the pressure variation of $`T_c`$, but $`\mathrm{\Delta }T_c`$ itself is a minimum near $`P^{}`$. It is as if the optimally homogeneous state is singular in the vicinity of $`P^{}`$.

From the data presented above, we construct a temperature-pressure phase diagram for CeCoIn5 given in Fig. 4a. This phase diagram is more like that of the cuprates around optimal doping than the canonical behavior of heavy-fermion systems exemplified, for example, by CeIn3 or CePd2Si2 mathur98 . There is no obvious long-range ordered state in CeCoIn5 at atmospheric or higher pressures. Rather, in view of its relationship to CeRhIn5, a reasonable speculation is that the antiferromagnetic QCP of CeCoIn5 is at a slightly negative pressure. The resistivity exponent $`n`$ and specific heat (ref. sparn01 ) are those of a non-Fermi-liquid up to $`P^{}`$ where $`T_c`$ is a maximum. Beyond $`P^{}`$, $`T_c`$ drops rapidly and parameters characterizing electronic transport are qualitatively different, Fermi-liquid-like at temperatures just above $`T_c`$ and those of a 3-D quantum-critical antiferromagnetic state at higher temperatures.

The generic $`TP`$ phase diagram in Fig. 4b provides a broader perspective of our experimental observations and their relationship to CeRhIn5. In the spirit of this diagram, the non-Fermi-liquid transport and thermodynamic properties of CeCoIn5 below $`P^{}`$ are controlled by an antiferromagnetic QCP at the inaccessible negative pressure denoted by $`P_c`$ in Fig. 4. We have chosen $`P_c`$ to correspond approximately to the critical pressure at which superconductivity appears in CeRhIn5. Spin-lattice relaxation measurements on CeRhIn5 identify kawasaki01  a pseudogap whose signature first appears at pressures slightly less than $`P_c`$ and at a temperature of $``$5 K and then decreases toward $`T_c`$ with increasing pressure until its signature is lost at $`PP_c`$. If we associate the temperature $`T_{pg}(P)`$ in Fig. 4a with a resistive signature for a pseudogap and extrapolate $`T_{pg}(P)`$ in CeCoIn5 to slightly negative pressures, there is a natural connection between our observations and properties of CeRhIn5. Neither the $`1/T_1T`$ measurements on CeRhIn5 nor our transport measurements on CeCoIn5 permits definitive statements about the origin or nature of the pseudogap pseudogap . We note, however, that 4-fold anisotropy of the thermal conductivity in the basal plane of CeCoIn5 persists izawa01  to $`3.2\mathrm{K}T_{pg}`$ at atmospheric pressure. If this in-plane modulation of the normal-state thermal conductivity is due to the presence of a pseudogap, it implies d-wave symmetry.

The schematic phase diagram, Fig. 4b, reflects the relationship among various phases of CeCoIn5 and CeRhIn5 and is similar to the $`T`$-doping phase diagram of the cuprates. This diagram and experimental data from which it is inferred indicate that the physics of heavy-fermion systems may be more closely related to that of the cuprates than previously appreciated. Aside from the structural layering in CeCoIn5 and CeRhIn5 and their relatively high $`T_c`$’s, they are similar in many respects to other Ce-based heavy-fermion materials in which superconductivity develops near a QCP. Further study of CeCoIn5, in parallel with the cuprates and other pressure-induced heavy-fermion superconductors, holds promise for bridging our understanding of the inter-relationship between unconventional superconductivity and quantum criticality across these interesting classes of correlated matter.

We thank R. Movshovich for helpful discussions. VS acknowledges S. M. Stishov for his support during the initial stage of this research. Work at Los Alamos was performed under the auspices of the US Department of Energy.


text:
# Non-invasive single-bunch matching and emittance monitor

## I Introduction and background

A quadrupole pick-up is a non-invasive device that measures the quadrupole moment

$$\kappa =\sigma _x^2\sigma _y^2+\overline{x}^2\overline{y}^2.$$
(1)
of the transverse beam distribution. Here, $`\sigma _x`$ and $`\sigma _y`$ are the r.m.s. beam dimensions in the horizontal and vertical directions, while $`\overline{x}`$ and $`\overline{y}`$ denote the beam position.

The practical use of quadrupole pick-ups was pioneered at SLACMiller et al. (1983), where six such pick-ups, distributed along the linac, were used. The emittance and Twiss parameters of a passing bunch were obtained from the pick-up measurements by solving a matrix equation, derived from the known transfer matrices between pick-ups.

In rings, the use of quadrupole pick-ups has largely focused on the frequency content of the raw signal. Beam width oscillations produce sidebands to the revolution frequency harmonics in the quadrupole signal, at a distance of twice the betatron frequency, and this can be used to detect injection mismatch. This was done at the CERN Antiproton Accumulator, where the phase and amplitude of the detected sidebands were also used to find a proper correction, using an empirical response matrixChohan et al. (1990). However, this measurement was complicated by the fact that the same sidebands can be produced by position oscillations, which demanded that position oscillations were kept very small.

In this paper, the idea behind the SLAC method is applied and further developed for use in rings. The quadrupole pick-ups used for the measurements presented here were specially developed for the CERN PS and optimized to measure the quadrupole momentJansson and Williams (2002).

They consist of four induction loops oriented to be sensitive to the magnetic flux in the radial direction (see FIG. 1). Since the field from a centered round beam has a flux only in the azimuthal direction, only deviations from roundness or position induce a signal in the loops. Therefore each loop is directly sensitive to the quadrupole moment, unlike previous instruments where the quadrupole moment was extracted by detecting tiny differences between four large electrode signals.

Two pick-ups have been installed in consecutive straight sections of the machineJansson et al. (2001). The optical parameters at their locations are given in TABLE 1. As shown later, it is crucial that the pick-ups be installed at locations with different ratios between horizontal and vertical beta value. The phase advance between pick-ups is also an important input parameter in the data analysis. In order to minimize the dependence of this phase advance on the programmed machine tunes and the beam intensity (space charge detuning), the pick-ups were installed as close as possible to each other.

The PS pick-ups provide both beam position and quadrupole moment information, with bunch-by-bunch resolution, over several hundred turns. Since the beam position is also measured, its contribution to the quadrupole moment can be subtracted, leaving only the beam-size related part, $`\sigma _x^2\sigma _y^2`$. Throughout the rest of this paper, when referring to $`\kappa `$, it will be assumed that this ‘artificial centering’ has been performed, unless stated otherwise.

An example of a position-corrected measurement is shown in FIG. 2, where the usefulness of the correction is clear. The initial beam size oscillation due to injection mismatch is clearly visible in the corrected signal. Note that beam-size oscillations are sensitive to the direct space charge, which means that they have a larger tune-spread, and therefore decohere much faster than beam position oscillations. The difference in decoherence time between beam size and position oscillations is therefore a rather direct measure of the incoherent tune shift. The detuning of the quadrupole signal frequencies can also be used to measure the incoherent tune shift, as has been done in the Low Energy Antiproton Ring (LEAR) at CERNChanel (1996).

## II Signal Acquisition and treatment

At the input to the data acquisition system, located in a building next to the machine, the analogue signals from the pick-up have the form

$`\mathrm{\Xi }(t)`$ $`=`$ $`Z_\mathrm{\Xi }\kappa (tt_\mathrm{\Xi })i(tt_\mathrm{\Xi })\text{quad. moment}`$ (2)
$`\mathrm{\Delta }_x(t)`$ $`=`$ $`Z_\mathrm{\Delta }\overline{x}(tt_{\mathrm{\Delta }_x})i(tt_{\mathrm{\Delta }_x})\text{hor. position}`$ (3)
$`\mathrm{\Delta }_y(t)`$ $`=`$ $`Z_\mathrm{\Delta }\overline{y}(tt_{\mathrm{\Delta }_y})i(tt_{\mathrm{\Delta }_y})\text{ver. position}`$ (4)
where the $`Z`$s are the transfer impedances, $`\kappa `$, $`\overline{x}`$ and $`\overline{y}`$ are defined as before, and $`i`$ is the beam current. In the PS, the beam current is not measured by the pick-up itself, so a separate beam current reference signal

$$\mathrm{\Sigma }(t)=Z_\mathrm{\Sigma }i(tt_\mathrm{\Sigma })\text{sum signal}$$
(5)
is taken from a nearby wall-current monitor. These analogue signals are sampled by digital oscilloscopes. The digitized signals are then re-sampled111Eventually, the digital re-sampling will be replaced by analogue delay lines to improve the noise performance. to correct for the signal timing differences $`t_\mathrm{\Xi }`$, $`t_{\mathrm{\Delta }_x}`$, $`t_{\mathrm{\Delta }_x}`$ and $`t_\mathrm{\Sigma }`$. These are mainly due to cable length differences, and have been measured both with a synthetic signal and using the beam.

The analysis of the data is made in a LabView program. In order to resolve single bunches, the data is treated in the time domain, considering each bunch passage separately. The first step in the analysis is to rid the signal of its intensity dependence, by normalizing to the measured beam current. The analysis is performed in two different ways, depending on whether the position and quadrupole moment are expected to be constant or varying along the bunch.

### II.1 Position and size constant along bunch

If there is no variation in position and size along the bunch, and one assumes that the quadrupole pick-up and the wall current monitor have the same frequency response, then the shape a given pulse must be exactly the same in all signals (apart from a baseline offset and noise effects). The normalization problem then consists in determining the scaling factor between a pulse in the beam current signal and the corresponding pulse on the pick-up outputs.

To do this, time slices of about one RF period centered on the bunch are selected. Each selected slice is a vector of $`N`$ samples and, under the above assumption, corresponding slices are proportional to each other. The quadrupole moment can therefore be found as the least squares solution to an overdetermined matrix equation, which in the case of the quadrupole signal has the form

$$\left(\begin{array}{cc}\mathrm{\Sigma }_1& 1\\ \mathrm{\Sigma }_2& 1\\ \mathrm{}& \mathrm{}\\ \mathrm{\Sigma }_N& 1\end{array}\right)\left(\begin{array}{c}\kappa \\ c\end{array}\right)=\frac{Z_\mathrm{\Sigma }}{Z_\kappa }\left(\begin{array}{c}\mathrm{\Xi }_1\\ \mathrm{\Xi }_2\\ \mathrm{}\\ \mathrm{\Xi }_N\end{array}\right).$$
(6)
The constant $`c`$ depends on the base line difference and is not used. The same calculation is performed for the position signals, and the position contribution to the quadrupole moment is then subtracted.

An attractive feature of this method, apart from noise suppression, is that the base line is automatically, and unambiguously, corrected for. Differences in frequency response of the two instruments could be corrected by filtering the signals, if these responses are known. However, such sophisticated corrections would enhance noise, and are not necessary in the PS.

### II.2 Position or size varying along bunch

Sometimes, there can be a variation in oscillation amplitude and phase along the bunch. At injection into the PS, there are two main causes for this

* The injection kicker pulse is not perfectly flat, which causes a variation of initial position along the bunch. The result is a fast position oscillation in those parts of the bunch that did not receive the correct kick.
* If the injected beam is longitudinally mismatched, the mismatched bunch will rotate in the bucket with the synchrotron frequency, causing the bunch length to oscillate. When the bunch is tilted in longitudinal phase space there is a correlation between energy and time, apparent as a variation of the mean energy along the bunch. The degree of correlation varies as the bunch rotates, and at a position with non-zero dispersion this gives rise to a slow head-tail oscillation at twice the synchrotron frequency. Both PS pick-ups are installed in dispersive regions, and are therefore sensitive to this effect.
* There can also be a variation of the beam dimensions along the bunch, as discussed toward the end of this paper.

In these cases, the basic assumption behind the algorithm described in the previous section is no longer valid. In fact, if the position varies along the bunch, any algorithm that calculates the average position and quadrupole moment of the bunch will give an erroneous result. Since

$$<x^2><x>^2$$
(7)
one can not simply use the average bunch position in Eq. (1) when correcting for the position. The correction must be done point-by-point along the bunch. For this purpose, a second normalization algorithm is used, which first establishes and subtracts the base line, and then calculates the position

$$x(t)=\frac{Z_\mathrm{\Sigma }}{Z_{\mathrm{\Delta }_x}}\frac{\mathrm{\Delta }_x(t)}{\mathrm{\Sigma }(t)}$$
(8)
as well as the quadrupole moment in each point. After this correction, an average beam quadrupole moment can be calculated, but it is also possible to study variations of the beam size along the bunch.

## III Beam-based calibration

### III.1 Internal signal consistency

One can take advantage of the position dependence of the quadrupole moment to make a consistency check between the position and quadrupole moment measurement of the pick-up, using data with large beam position oscillations but stable beam size. Since the beam size oscillations damp away much faster than beam position oscillations, such data can easily be obtained at injection by an appropriate trigger delay. A plot of expected versus measured variation of the quadrupole moment with beam position is shown in FIG. 3, showing a good agreement. This test can easily be automated, and is a good indicator of whether the beam position correction works well.

### III.2 Comparison with wire-scanners

The standard method for emittance measurement on a circulating beam in the PS is the fast wire-scanner. In order to test the calibration of the pick-ups, measurements were done on several different stable beams, approximately 15 ms after injection. The quadrupole pick-up signal was acquired over 200 machine turns, at the same time as the wire traversed the beam. The comparative measurement was performed on all the operational beams available in the machine, with the exception of the very high intensity beams that saturate the pick-up amplifiers. Thus there was a significant difference in both beam and machine parameters between the different measurements. This was done in an attempt to randomize any systematic errors. The beam parameters are given in TABLE 2, where the different beams have been tagged with their operational names.

The r.m.s. variation in the measured quadrupole moments from turn to turn was of the order of 0.2-0.5 mm2, depending on the beam intensity. Assuming that the beam size was perfectly stable, this gives an estimate of the single-turn resolution of the pick-up measurement. Also the wire-scanner measurements were stable, although for some beams there was a systematic disagreement between the two wire-scanners measuring in the same plane.

To compare the two instruments, the emittances measured with the wire-scanners were used to calculate the expected quadrupole moment at the locations of the pick-ups. The momentum spread required for both the wire scanner measurement and the subsequent calculation was obtained by a tomographic analysis of the bunch shapeHancock et al. (1998). The propagated systematic error was estimated on the assumption that the wire-scanner accuracy is 5% in emittance, the beta function at the pick-ups is known to 5%, the dispersion to 10% and the momentum spread to 3% accuracy. These estimates are rather optimistic, but give considerable propagated errors for certain measurement points. For simplicity, possible correlations between errors (e.g. beta function errors at different locations in the machine) were ignored, and all different error sources were added in quadrature. To accentuate the cases with wire-scanner disagreement, each of the four different ways of combining the two horizontal and two vertical wire-scanners was calculated separately and displayed as separate points. The result is shown in FIG. 4.

Overall, the measured data seem to indicate that the offsets are slightly smaller than measured in the lab, which could be explained by the fact that the pick-ups were dismantled in the lab to be moved to the machine. However, the effect is within the error-bar, and no strong conclusion can therefore be made. Moreover, the pick-ups have been dismantled and rebuilt in the lab, without effect on the measured offsets.

The point corresponding to the EASTC beam appears to disagree somewhat in both planes, although the effect is just about within the error-bar. There are a number of possible explanations for this:

* The PS is operated in a time-sharing mode, where a so-called super-cycle containing a certain number (usually 12) of beam cycles is repeated over and over again. At the time of the measurement, the super-cycle contained several instances of the EASTC beam, and it is known from experience that the position within the super-cycle can affect the beam characteristics. For this particular measurement, it is not guaranteed that the measurements with the two instruments were done on the same instance of the beam, whereas for all other measurements there was either only one instance of the beam in the super-cycle, or the acquisition was locked to a certain instance. Some fluctuations of the measured value were also observed.
* The EASTC beam has a large momentum spread and a horizontal tune close to an integer resonance. Theory indicates that the correction quadrupoles used to obtain this working point can perturb the dispersion function by more than 15%Carli , which would affect both the accuracy of the wire-scanner measurement, and the subsequent calculation of the expected quadrupole moment. Studies of this effect are planned for the 2002 run.

The general conclusion from the measurement series is that the wire-scanner and quadrupole pick-up agree within the error bar. The systematic errors due to optics parameters make it impossible to detect with certainty any difference in pick-up behavior between the laboratory measurements with a simulated beam, and the measurements on real beam in the machine. In order to calibrate the pick-ups accurately using the beam, the wire-scanners and the pick-up should be situated in the same straight section, which is excluded in the PS due to space limitations.

### III.3 Comparison with turn-by-turn profile measurement

Comparative measurements of injection matching have been done using a SEM grid with a fast acquisition systemBenedikt et al. (2001), that can measure beam profiles turn-by-turn for a single bunch. This is a destructive device and can only be used in rare dedicated machine development sessions. It is also limited both in bandwidth and maximum beam intensity, and therefore it has not been possible to make a full systematic study on beams with different characteristics. Instead, a special beam was prepared, with low intensity to spare the grid, and long bunches due to the bandwidth limitations.

The SEM grid data was used to calculate the expected value of the quadrupole moment at the pick-up locations, using the beta values, dispersion, and relative phase advance in Table 1. The results are shown in FIG. 5, and show a rather good agreement with what was actually measured with the pick-ups. The small differences can be accounted for by systematic error sources, i.e. the optical parameters used in the comparison.

## IV Emittance measurement

When the circulating beam is stable, the quadrupole moments of a given bunch, as measured by the two pick-ups, are constant and given by

$`\kappa _1`$ $`=`$ $`ϵ_x\overline{\beta }_{x\mathrm{\hspace{0.17em}1}}ϵ_y\overline{\beta }_{y\mathrm{\hspace{0.17em}1}}+\overline{D}_{x\mathrm{\hspace{0.17em}1}}^2\sigma _\mathrm{p}^2`$
$`\kappa _2`$ $`=`$ $`ϵ_x\overline{\beta }_{x\mathrm{\hspace{0.17em}2}}ϵ_y\overline{\beta }_{y\mathrm{\hspace{0.17em}2}}+\overline{D}_{x\mathrm{\hspace{0.17em}2}}^2\sigma _\mathrm{p}^2`$ (9)
where $`ϵ`$ denotes the emittance, $`\beta `$ the beta value, $`D`$ the dispersion, and $`\sigma _\mathrm{p}`$ the relative momentum spread. The bar over certain parameters indicate that these are properties of the lattice, to be distinguished from the corresponding beam properties (typeset without the bar).

When the momentum spread is known, the system of equations can be solved for the emittances if

$$\frac{\overline{\beta }_{x\mathrm{\hspace{0.17em}1}}}{\overline{\beta }_{y\mathrm{\hspace{0.17em}1}}}\frac{\overline{\beta }_{x\mathrm{\hspace{0.17em}2}}}{\overline{\beta }_{y\mathrm{\hspace{0.17em}2}}}$$
(10)
which explains the earlier statement about the requirement on the beta functions at the pick-up locations. If the ratio between horizontal and vertical beta function is significantly different at the two locations, the equations are numerically stable. Thus measuring the emittance of a stable circulating beam with quadrupole pick-ups is in fact rather straightforward.

Statistical errors due to random fluctuations in the measurement of $`\kappa `$ can, although they are usually small, be reduced by averaging over many consecutive beam passages. The dominant errors are therefore systematic, coming from offsets in the pick-ups and errors in the beta functions, lattice dispersion and momentum spread. The pick-up offsets are, however, known from test bench measurements. Furthermore, by comparing the amplitude of position oscillations as measured by the two pick-ups, the ratios $`\overline{\beta }_{x\mathrm{\hspace{0.17em}1}}/\overline{\beta }_{x\mathrm{\hspace{0.17em}2}}`$ and $`\overline{\beta }_{y\mathrm{\hspace{0.17em}1}}/\overline{\beta }_{y\mathrm{\hspace{0.17em}2}}`$ can be determined.

The main uncertainty is thus the absolute value of the beta function, as for almost any other emittance measurement (e.g. wire-scanner). The accuracy can therefore be expected to be comparable to that of a wire-scanner. An emittance measurement using the pick-up system is shown in FIG. 6, and compares well with wire-scanner results.

Note that with three pick-ups, suitably located, the momentum spread could also be measured.

## V Matching measurement

Even though quadrupole pick-ups can be used to measure filamented emittance, the main reason for installing such instruments in the machine is to be able to measure betatron and dispersion matching at injection, as no other instrument (apart from the destructive SEM-grid) is able to do this. One would like not only to detect mismatch, but also to quantify the injection error in order to be able to correct it.

### V.1 Matrix inversion method

To determine the parameters of the injected beam, the SLAC methodMiller et al. (1983) based on matrix inversion could be directly applied, since the quadrupole moment is measured on a single-pass basis. An advantage when performing this measurement in a ring, as compared to measuring in a linac, is that each pick-up can be used several times on the same bunch. Therefore it is enough to use two pick-ups instead of six, which reduces both the hardware cost and the systematic error sources. It is also straightforward to improve on statistics by increasing the number of measured turns, thereby reducing noise. Another advantage in a ring is that the periodic boundary conditions reduce the number of parameters needed to calculate the matrix. Many of these parameters (tunes, phase advance between pick-ups, ratios between beta function values) can also be easily measured, which means that the matrix can be experimentally verified.

However, the matrix method was developed for a linac, and does not take full advantage of the properties of a ring. Also, it does not include dispersion effects, and it is necessary to make assumptions on the space charge detuning when calculating the matrix.

### V.2 Parametric fit method

In a ring, the turn by turn evolution of the beam envelope, and therefore the quadrupole moment, can be expressed in a rather simple analytical formula. Expanded in terms of the optical parameters, the quadrupole moment of a beam is given by

$$\kappa =\sigma _x^2\sigma _y^2=\epsilon _x\beta _x\epsilon _y\beta _y+\sigma _p^2D_x^2\sigma _p^2D_y^2$$
(11)
assuming linear optics with no coupling between planes (note that there are no bars, i.e. the optical parameters here refer to the beam). If the beam is initially mismatched in terms of Twiss functions or dispersion, the value of $`\kappa `$ will vary with the number of revolutions $`n`$ performed asJansson (2001)

$$\begin{array}{c}\kappa _n=\overline{\beta }_x(\epsilon _x+\mathrm{\Delta }\epsilon _x)\overline{\beta }_y(\epsilon _y+\mathrm{\Delta }\epsilon _y)+\overline{D}_x^2\sigma _\mathrm{p}^2\hfill \\ \hfill +\overline{\beta }_x\epsilon _x\delta _{\beta _x}\mathrm{cos}(2\nu _xn\varphi _{\beta _x})+\overline{\beta }_x\sigma _\mathrm{p}^2\delta _{D_x}^2\mathrm{cos}(2\nu _xn2\varphi _{D_x})\\ \hfill \overline{\beta }_y\epsilon _y\delta _{\beta _y}\mathrm{cos}(2\nu _yn\varphi _{\beta _y})\overline{\beta }_y\sigma _\mathrm{p}^2\delta _{D_y}^2\mathrm{cos}(2\nu _yn2\varphi _{D_y})\\ \hfill +\sqrt{\overline{\beta }_x}\sigma _\mathrm{p}^2\overline{D}_x\delta _{D_x}\mathrm{cos}(\nu _xn\varphi _{D_x})\end{array}$$
(12)
Here, $`\nu _x`$ and $`\nu _y`$ are the fractional tunes expressed in radians, and $`\mathrm{\Delta }ϵ`$ denotes the emittance increase caused by the mismatch. The first line contains constant terms, and also gives the steady state value that will be reached when the oscillating components have damped away.

The two middle lines of Eq. (12) are signal components at twice the horizontal and vertical betatron frequencies. They arise from both dispersion and betatron mismatch. The betatron mismatch is parametrized by

$$\stackrel{}{\delta }_{\beta _x}=\left(\begin{array}{c}\frac{\beta _x}{\overline{\beta }_x}\frac{\overline{\beta }_x\gamma _x+\overline{\gamma }_x\beta _x2\overline{\alpha }_x\alpha _x}{2}\\ \frac{\overline{\alpha }_x\beta _x\alpha _x\overline{\beta }_x}{\overline{\beta }_x}\end{array}\right)\left(\begin{array}{c}\frac{\mathrm{\Delta }\beta _x}{\overline{\beta }_x}\\ \overline{\alpha }_x\frac{\mathrm{\Delta }\beta _x}{\overline{\beta }_x}\mathrm{\Delta }\alpha _x\end{array}\right)$$
(13)
where the last approximation is valid for small mismatch. Here, the shorthand notation $`\mathrm{\Delta }\beta =\beta \overline{\beta }`$ and $`\mathrm{\Delta }\alpha =\alpha \overline{\alpha }`$ is used for the difference between lattice and beam value.

The fourth line of Eq. (12) is a signal at the horizontal betatron frequency, which is due to dispersion matching. This mismatch is parametrized by the vector

$$\stackrel{}{\delta }_{D_x}=\left(\begin{array}{c}\frac{\mathrm{\Delta }D_x}{\sqrt{\overline{\beta }_x}}\\ \frac{\overline{\beta }_x\mathrm{\Delta }D_x^{}+\overline{\alpha }_x\mathrm{\Delta }D_x}{\sqrt{\overline{\beta }_x}}\end{array}\right)$$
(14)
where, again, shorthand notation ($`\mathrm{\Delta }D=D\overline{D}`$ and $`\mathrm{\Delta }D^{}=D^{}\overline{D}^{}`$) is used. There is no corresponding signal at the vertical betatron frequency due to the absence of vertical lattice dispersion. Therefore, it is not possible to distinguish vertical dispersion mismatch from vertical betatron mismatch by studying the quadrupole signal. However, one does not usually expect a large vertical dispersion mismatch.

The steady state (filamented) emittance is given by

$$\begin{array}{c}\epsilon _x+\mathrm{\Delta }\epsilon _x=\epsilon _x\frac{1}{2}\left(\overline{\beta }_x\gamma _x+\overline{\gamma }_x\beta _x2\overline{\alpha }_x\alpha _x\right)+\hfill \\ \hfill +\sigma _\mathrm{p}^2\frac{(\mathrm{\Delta }D_x)^2+(\overline{\beta }_x\mathrm{\Delta }D_x^{}+\overline{\alpha }_x\mathrm{\Delta }D_x)^2}{\overline{\beta }_x}\\ \hfill \epsilon _x+\epsilon _x\frac{|\stackrel{}{\delta }_{\beta _x}|^2}{2}+\sigma _\mathrm{p}^2\frac{|\stackrel{}{\delta }_{D_x}|^2}{2}\end{array}$$
(15)
where, again, the last approximation is valid for small betatron mismatch222There is also a contribution to the emittance increase due to injection miss-steering that is not included here, since normally coherent dipole oscillations filament much slower than quadrupole oscillations, and the beam position contribution is subtracted from the signal..

By fitting the above function to the data, the injected emittances, the betatron mismatches in both planes, and the horizontal dispersion mismatch are directly obtained. The tunes can also be free parameters in the fit, which automatically estimates and corrects for space charge detuning. An example of a fit to measured data is shown in FIG. 7. A requirement for a good fit convergence is, as when measuring filamented emittance, that the ratio between beta functions should be different at the pick-up locations. Also, the tunes must be such that enough independent data points are obtained. In other words, if the quadrupole signal is repetitive, it must have a period larger than the minimum number of turns required for the fit. In the PS, this means that the working point $`Q_h=Q_v=6.25`$, which is close to the bare tune, should be avoided. The fit result is also less stable in the vicinity of this working point, and when the tune in only one of the planes is close to 6.25. With two pick-ups, at least five machine turns (10 data points) are required for the fit, if the tunes are also free parameters. Some more turns can be used to check the error, but the maximum number of turns is limited by decoherence, as discussed below.

Note that since the beam size oscillations due to dispersion mismatch are also detuned by space charge, measuring the dispersion component separately (by changing the energy of the beam and measuring the coherent response) would result in an accumulated phase error in the dispersion term.

### V.3 The effect of decoherence

The fit function above does not include the effect of decoherence (damping) of the beam width oscillations. Fortunately, due to the physics of the decoherence process, the decay of the oscillation amplitudes is not exponential as for many other damping phenomena. If the beam is approximated by an ensemble of harmonic oscillators with a tune distribution $`\rho (\mathrm{\Delta }Q)`$ and an average tune Q, its coherent response to an initial displacement is

$$x(s)=e^{i2\pi Qs}\underset{A(s)}{\underset{}{A_0_{\mathrm{}}^{\mathrm{}}e^{i2\pi \mathrm{\Delta }Qs}\rho (\mathrm{\Delta }Q)d(\mathrm{\Delta }Q)}}$$
(16)
and the derivative of the amplitude function

$$\frac{A}{s}_{\mathrm{}}^{\mathrm{}}e^{i2\pi \mathrm{\Delta }Qs}\mathrm{\Delta }Q\rho (\mathrm{\Delta }Q)d(\mathrm{\Delta }Q)$$
(17)
is zero at $`s=0`$, i.e. initially the amplitude is unchanged by the decoherence process. A plot of the amplitude versus time for some tune distributions is given in FIG. 8, showing that the initial behavior is also largely independent of the distribution.

In reality, the tune of each individual particle is changing with time (e.g. due to synchrotron motion), and therefore the decoherence pattern is more complicated. However, synchrotron motion is negligible for the first few turns. Therefore, data analysis is greatly simplified and accuracy is improved, if one limits the number of turns to a rather small value. This also demonstrates an advantage of the fit method over an FFT analysis of the signals, since an FFT needs many points to achieve good frequency resolution.

### V.4 Measurement results

To test the injection matching measurement, a series of measurements was done with different settings of some focusing elements of the PS injection line. An example of such a measurement is shown in FIG. 9, where a quadrupole was changed in steps of $`10`$A, and the resulting variation of the fit parameters recorded. The variation of the different error vectors expected from beam optics theory is also shown, and there is a rather good agreement, both in direction and magnitude of the changes. The injected emittances are unchanged, as expected.

By using the theoretical response matrix for dispersion and betatron matching, a proper correction to the measured error can be calculatedGiovannozzi et al. (1998). So far, actual corrections of the measured mismatches have not been made, since the dominant error (the dispersion mismatch) can not be corrected without a complete change of optics of the entire line. Studies for a new dispersion matched optics are underway.

While the dispersion mismatch is large for all beams, due to the transfer line design, the level of betatron mismatch varies between different operational beams. Most high intensity beams measured were observed to be fairly well matched, whereas some lower intensity beams had a significant mismatch. This might be explained by the fact that mismatch is likely to cause losses for aperture limited beams, and therefore the process of intensity optimization leads to well-matched beams, although the mismatch is never directly measured. This indirect matching mechanism is absent for the future bright LHC beam, and it can therefore be expected to develop a relatively large mismatch if not continuously monitored and corrected.

## VI Measurement within the bunch

As mentioned earlier, the transverse mean position can sometimes vary along the bunch. However, in some cases, also the beam size itself varies along the bunch. This is notably the case for high intensity beams that are highly non-Gaussian. For a Gaussian beam distribution, the transverse bunch width is constant along the bunch. This is because the multi-dimensional Gaussian is just a product of one-dimensional Gaussians. However, for a parabolic beam this is no longer true, as may be easily verified analytically. With the pick-ups, it is possible to measure the quadrupole moment as a function of position within the bunch. The measurement is good over most of the bunch, but naturally gets very noisy and prone to systematic errors in the head and tail, since these regions are sparsely populated. A measurement made on a stable beam is shown in FIG. 10. The plot also shows the same measurement with the dispersion contribution subtracted333The momentum spread as a function of position within the bunch was obtained from a tomographic analysisHancock et al. (1998) of the bunch shape data., indicating that the variation of beam size along the bunch is mainly due to variations in momentum spread. This fits with the fact that the longitudinal bunch distribution is usually non-Gaussian. Applying the methods discussed earlier on the dispersion corrected data, it is also possible to calculate the emittance variations along the bunch.

## VII Conclusions

The quadrupole pick-ups recently built and installed in the PS machine have been evaluated in a series of measurements. These pick-ups measure both injection matching and emittance for a single, selected bunch in the machine. The measurement can be made parasitically, without perturbing the beam, because the devices are non-intercepting.

Comparison with other instruments in the machine show good agreement. All observed deviations are within the estimated systematic error bars. The systematic errors come mainly from the imperfect knowledge of beta value and dispersion needed to evaluate the data. Systematic errors are indeed expected to dominate the total error in the quadrupole pick-up measurement, as is the case for most emittance measurement devices.

For matching applications, the pick-ups can be used to determine phase and amplitude of horizontal and vertical betatron mismatch, as well as horizontal dispersion mismatch. This analysis can be done individually on each injected bunch. Since the mismatch is detected as an oscillation, the effect of systematic errors (e.g. pick-up offsets) is not very important.

As emittance measurement devices, the pick-ups have some interesting properties. The single turn resolution makes it possible to measure and follow the evolution of the emittance over many turns (limited only by acquisition memory). When measuring filamented emittance, it possible to reduce the effect of noise by averaging over many turns, and also to check that the beam is stable during the measurement, something that is assumed but not actually verified during a wire scanner measurement. More important, the pick-ups have no moving parts that wear out, as is the case for a wire-scanner. This makes it possible to create a watchdog application to monitor the evolution of the emittances pulse by pulse over a long period. In such an application, systematic errors are again of lesser importance, since variations rather than absolute values are sought.

The pick-ups can also be used to study variations of the emittance along the bunch, although this may be mainly of academic interest.

###### Acknowledgements.

The author would like to thank D.J. Williams and L. Søby for their support and important contributions to the pick-up hardware; H. Koziol for reading and commenting on an early draft to this paper; U. Raich and C. Carli for contributing to the turn-by-turn SEM-grid data acquisition and analysis.


text:
# Generic Sandpile Models Have Directed Percolation Exponents

## Abstract

We study sandpile models with stochastic toppling rules and having sticky grains so that with a non-zero probability no toppling occurs, even if the local height of pile exceeds the threshold value. Dissipation is introduced by adding a small probability of particle loss at each toppling. Generically, for models with a preferred direction, the avalanche exponents are those of critical directed percolation clusters. For undirected models, avalanche exponents are those of directed percolation clusters in one higher dimension.

In recent years, there has been a lot of interest in the study sandpile models as models of real granular matter , and also as paradigms of self-organized critical systems in general. Following the well-known work of Bak, Tang and Wiesenfeld, many different types of sandpile models with different toppling rules have been studied  : deterministic and stochastic, with or without preferred direction, different instability criteria , or particle distribution rules , with fixed energy  etc.. However, a clear picture of the factors that determine different universality classes of critical behavior is yet to emerge .

A different paradigm for non-equilibrium critical phenomena has been directed percolation (DP) which is believed to describe active to absorbing state transition in a wide class of reaction-diffusion systems . The activity of avalanches in sandpile models can grow, diffuse or die, and any stable configuration is an absorbing state. Thus, this should be in the universality class of DP with many absorbing states . Several models of self-organized criticality which show critical exponents related to DP have been studied earlier. However, these models do not involve any conserved field. The critical exponents of known models with conservation of sand are very different from those of DP, and this is presumably due to role of the local conservation of sand in the model. In , a model with conservation of particles showing DP exponents was studied, but this study was mainly numerical.

In this Letter we study several sandpile models with stochastic toppling rules, both directed and undirected. The grains are ‘sticky’ in our models in the sense that there is a nonzero probability that any grains arriving at a site during the avalanche process just get stuck there. We find that the distribution function of avalanches in these models has the same power law distribution as that of the critical DP clusters. Our theoretical arguments, supported by numerical simulations, show that generically these models belong to the DP universality class. The previously studied deterministic models , and the stochastic toppling model introduced by Manna  are unstable to perturbations and flow to the DP fixed point.

The relation of sandpile models to DP was attempted earlier in  using sticky grains. However, in that paper, no bulk dissipation was present, and sand was dissipated only at the boundaries. The boundaries break the translational invariance of the steady state. The density becomes space dependent, and the density profile affects the statistics of avalanches. The critical exponents of avalanches are not those of critical DP clusters, but expressible in terms of DP exponents . Introducing bulk dissipation in this paper, we are able to get rid of these problems, and the relationship between the sandpile and DP problems becomes more transparent.

For simplicity, we start by defining the directed model on a $`(1+1)`$-dimensional square lattice. Various generalizations to higher dimensions, undirected case, and other toppling rules are straight forward, and briefly discussed at the end. The sites on an $`L\times M`$ torus are labelled by euclidean coordinates $`(i,j)`$ with $`(i+j)`$ even and $`j`$ increasing downwards. At each site $`(i,j)`$, there is a non-negative integer $`h_{i,j}`$ to be called the height of the pile at that site. Initially all $`h_{i,j}`$ are zero. The system is driven by choosing a site at random and increasing the height at that site by one. If one or more particles are added to a site at time $`t`$ (from outside or from other sites), and its height becomes greater than $`1`$, then it is said to become unstable at time $`t`$.

Any site $`(i,j)`$ made unstable at time $`t`$ relaxes at the time $`(t+1)`$ stochastically: With probability $`(1p)`$, it becomes stable without losing any grains. We say that the added particle(s) sticks to the existing grains. Otherwise (with probability $`p`$), the relaxation occurs by toppling in which the height at the site decreases by two, and the site becomes stable. We introduce bulk dissipation, by assuming that in each toppling, there is a small probability $`\delta `$ that both grains from the toppling are lost, other wise (with probability $`1\delta `$), the two grains are transferred to the two downwards neighbors $`(i\pm 1,j+1)`$.

Note that there is a nonzero probability that a stable site can have arbitrarily large heights. We relax all the unstable sites by parallel dynamics. A site made unstable at time $`t`$ is relaxed in one step at time $`t+1`$, independent of whether it received one or more grains at the previous time step. Once a site has relaxed, it remains stable until perturbed again by new grains coming to the site. This relaxation process is repeated until all sites become stable, and then a new grain is added.

The model is specified by two parameters $`p`$ and $`\delta `$. The case $`p=1,\delta =0`$ is exactly soluble, and its critical exponents are known in all dimensions . In this case, one has to introduce open boundary conditions to ensure the existence of a steady state. In two dimensions, the probability that adding a particle will cause an avalanche of s topplings varies as $`s^{\tau _s}`$ for large $`s`$ with $`\tau _s=4/3`$. The probability that the duration of avalanche equals $`T`$ varies as $`T^{\tau _t}`$ with $`\tau _t=3/2`$. The case $`\delta =0,p`$ arbitrary was studied earlier in  discussed in the introduction.

Fig. 1 The line $`p=p^{}(\delta )`$ separating the steady state and no-steady-state regions in the $`p`$-$`\delta `$ plane. The inset shows a plot of $`[p_cp^{}(\delta )]`$ versus $`\delta `$. The straight line shows the theoretical slope $`1/\gamma `$.

The probability distribution of different configurations in the steady state for general values of $`p`$ and $`\delta `$ have an interesting structure. Define a variable $`g_{i,j}=h_{i,j}(mod2)`$. We group together different stable configurations $`\{h_{i,j}\}`$ corresponding to the same values of $`\{g_{i,j}\}`$. There are $`2^{LM}`$ such equivalence classes. On addition of a particle to a site, the $`g`$ variable at that site flips. A toppling leads to flipping of $`g`$’s at the two downward sites. From detailed balance, it follows that in the steady state, each of the $`2^{LM}`$ equivalent classes occurs with equal probability .

If $`p=1`$, with $`\delta `$ arbitrary, the only allowed height values are $`0`$ and $`1`$. In this case, we get a full characterization of the steady state. The $`n`$-point correlation functions satisfy linear equations, and exact solution of  can be generalized to arbitrary $`\delta `$. We omit the details here. The distribution of avalanche-sizes has an exponential decay for non-zero $`\delta `$.

As $`p`$ is decreased below $`1`$, heights greater than $`1`$ appear with non-zero probability in the steady state. The mean height $`h_{i,j}`$ increases as $`p`$ is decreased, and there exists a critical value $`p^{}(\delta )`$ such that for $`pp^{}(\delta )`$, the height of the pile increases without bound, and there is no steady state (Fig. 1).

The absence of a steady state is obvious along the the line $`p=0`$, with $`\delta `$ arbitrary, as in this case, different sites cannot influence each other, and each added particle just sticks to the existing pile. A similar decoupling occurs along the line $`\delta =1`$, with $`p`$ arbitrary. In this case, the average particle loss per added grain is $`2(1c_0)p`$, where $`c_0`$ is the density of sites with height $`0`$. In the steady state this must equal $`1`$. As $`c_00`$, a steady state can exist only for $`p>1/2`$. Thus $`p^{}(\delta =1)=1/2`$.

We now derive the equation for the boundary line $`p=p^{}(\delta )`$. At the phase boundary, clearly $`c_0=0`$. In the growth of an avalanche in the steady state of the system, any site which receives at least one grain from its upward neighbors sends grains to downward neighbors with a probability $`p(1\delta )`$. Thus the probability $`Prob(s)`$ that an added grain will cause an avalanche in which at least $`s`$ sites transfer particles to downward neighbors is same as the probability $`Prob_{DP}(s|\stackrel{~}{p})`$ of a cluster of at least $`s`$ sites in a directed site percolation process with concentration of active sites$`=\stackrel{~}{p}`$. A

$$Prob(s)=Prob_{DP}(s|\stackrel{~}{p}=p(1\delta )).$$
(1)
Note that Eq.(1) also holds in the entire ‘no steady state’ phase, where mean height continues to increase, but $`Prob(s)`$ tends to a limiting distribution for large times.

By particle conservation, in the steady state, the average number of topplings in an avalanche must be equal to $`1/(2\delta )`$. Let $`n_{_{DP}}(\stackrel{~}{p})`$ is the sum of the average number of occupied and perimeter sites in the cluster corresponding to a randomly picked site in the DP problem. Now, the equation for the phase boundary can be expressed in terms of the function $`n_{_{DP}}(\stackrel{~}{p})`$ as

$$n_{_{DP}}((1\delta )p^{}(\delta ))=1/(2\delta ).$$
(2)
For small $`\delta `$, the average size of clusters is large, and $`\stackrel{~}{p}`$ is near the critical probability $`p_c`$ for the directed site percolation on this lattice, the function $`n_{DP}`$ is known to vary as $`(p_c\stackrel{~}{p})^\gamma `$. Substituting this in Eq.(2) we get

$$p^{}(\delta )=p_cA\delta ^{1/\gamma }+\mathrm{terms}\mathrm{of}\mathrm{higher}\mathrm{order}\mathrm{in}\delta ,$$
(3)
where $`A`$ is some constant, and $`\gamma `$ is the susceptibility exponent of the DP problem. In particular, we have $`p^{}(\delta =0)=p_c`$. The inset in Fig. 1 shows a log-log plot of the numerically determined values of $`p_cp^{}(\delta )`$ plotted versus $`\delta `$. We get good agreement with Eq.(3) using the existing estimates $`p_c=0.70548515`$ and $`\gamma =2.277730`$ .

Since the average number of topplings in an avalanche in the steady state is $`1/(2\delta )`$, we will see self- organized criticality with long-ranged correlations only in the limit $`\delta 0^+`$, and $`p_cp1`$, (marked with a heavy line in Fig. 1). For $`\delta =0^+`$ and $`1p>p_c`$, in the steady state $`c_0>0`$. Each site is characterized by two probabilities, $`p_1`$ and $`p_2`$, which correspond to the probability that a site topples when one or two particles are added to it respectively. Clearly, $`p_1=(1c_0)p`$ and $`p_2=p`$. The correlations between heights at different sites are small, and if we ignore them , the evolution of avalanches (Fig. 2) in the steady state is as in the Domany-Kinzel model of DP with parameters $`(p_1,p_2)`$ .

Fig. 2 Picture of a typical avalanche for the 2-d directed model (A) and the time-evolution of the 1-d undirected model (B) for $`p=.873>p_c`$ and $`\delta =.0001`$ are compared with the clusters at the critical line of the Domany-Kinzel model of DP with $`p_2=p`$.

Even if some short-range correlations are present, they should not modify critical behavior, which is expected to be same as in DP. In Fig. $`3`$, we have compared the the probability distribution of avalanches with that of DP clusters. For $`s1`$, the latter is expected to satisfy the scaling form

$$Prob(s)\frac{A}{s^{\tau 1}}f[s/s^{}],$$
(4)
where $`A`$ is some amplitudes, $`\tau `$ is a critical exponent for cluster size distribution in DP and the function $`f(x)`$ tends to a finite constant as $`x`$ tends to zero, and decreases exponentially with $`x`$ for large $`x`$. We assume that $`s^{}`$ varies as $`\delta ^\varphi `$. Using the constraint $`s\delta ^1`$, we get $`\varphi =1/(2\tau )`$. Using the known numerical estimate $`\tau =1.108`$ in $`d=1+1`$ , we find a very good collapse when $`s^{\tau 1}Prob(s)`$ is plotted versus $`s\delta ^\varphi `$ for two different values $`p=p_c`$ and $`p=0.873`$, and two values of $`\delta =10^3`$ and $`\delta =10^4`$. In the inset, the scaling function $`f`$ is compared with that for DP clusters. We get an excellent collapse, a strong evidence that the two functions are the same, and the correlations in heights in the steady state are irrelevant.

It is straightforward to extend the previous discussion to higher dimensions. Thus the avalanche exponents in the $`(d+1)`$-dimensional model are the same as the exponents of cluster size distribution in the $`(d+1)`$-dimensional DP at critical point. The upper critical dimension is $`d=4`$.

Consider now the undirected version of the problem on a $`d`$-dimensional hypercubical lattice. The rules are the same as before, except that on toppling, a particle is transferred to each of the $`2d`$ neighbors of the toppling site. Clearly, in this case also, there is no steady state for small $`p`$. For $`p=p^{}(\delta )`$, the mean height per site diverges. Then, any site which has one or more particles added to it, sends particles to its neighbors with probability $`p(1\delta )`$. If we look at the space-time history of the evolution of the avalanche, we get a directed site-percolation cluster on a $`(d+1)`$-dimensional body centered hypercubical (bch) lattice. The phase boundary in this undirected model is also the same as that of the directed model on the $`(d+1)`$-dimensional bch lattice.

Fig. 3 Scaling collapse of $`s^{\tau 1}\mathrm{Prob}(s)`$ versus $`s\delta ^\varphi `$ for four different combinations of parameters $`p,\delta `$. The inset compares the scaling function for the directed and undirected models with $`(p,\delta )=(0.873,10^4)`$ with that for DP at $`p=0.7`$.

For $`p>p^{}(\delta )`$, the undirected model differs from the directed model in that the height of the pile at a site does not change between two topplings. This may give rise to possible long term memory effects. However, in our model this effect is rather small, as the probability of toppling depends on height only if the height is zero. Along the line $`p=p^{}(\delta )`$, the memory affect is strictly absent, as the density of sites with zero height goes to zero. We find that even for $`p`$ as large as $`.873`$, the avalanches are qualitatively similar to DP (Fig. 2), and the distribution function is also indistinguishable from that of near-critical DP clusters (see inset of Fig. 3). There is a crossover from deterministic limit ($`p=1,\delta =0`$)  to DP for $`p1`$.

Why does the conservation of particles not change the critical behavior away from DP in our problem? Consider a simple DP process of particles of type X, on a lattice with $`n_X(\stackrel{}{r})`$ particles of X at site $`\stackrel{}{r}`$. We now attach a register $`n_Y(\stackrel{}{r})`$ at each site $`\stackrel{}{r}`$, which decreases (increases) by 1 each time a particle X is created ( destroyed) at $`\stackrel{}{r}`$. Then clearly, we have a local conservation of $`n_X+n_Y`$. Clearly, if the dynamics of the X particle is not affected by the book keeping, the process still belongs to the DP universality class. In our model, $`n_Y`$ is the height of the pile. It fluctuates about its mean value, but there is an influence of $`n_Y`$ on the dynamics of X particles as birth of X particles is not allowed if $`n_Y`$ is zero.

The phenomenological evolutions for the coarse-grained density fields $`n_X(\stackrel{}{r})`$ and $`n_Y(\stackrel{}{r})`$ in the conservative limit $`\delta =0^+`$ may be written as

$`_tn_X=^2n_X+an_X`$ $`bn_X^2+cn_Xg(n_Y)+\eta (\stackrel{}{r},t)`$ (5)
$`_t(n_X+n_Y)=^2n_X`$ (6)
where $`\eta (\stackrel{}{r},t)`$ is a noise term, and $`a,b`$ 6and $`c`$ are phenomenological constants. The crucial term here is the coupling term $`cn_Xg(n_Y)`$. Vespignani et al  chose $`g(n_Y)`$ proportional to $`n_Y`$. In our case, $`n_Y`$ has a threshold and the effect on $`n_X`$ saturates to a finite value even as $`n_Y`$ increases to infinity. The simplest choice of $`g(n_Y)`$ to model this is to choose $`g(n_Y)=\theta (n_Yh_c)`$, where $`\theta `$ is the step function and $`h_c`$ is the threshold value. In naive power counting, this term has the same scaling as the linear term in $`n_X`$. We are not able to treat this analytically. But the results of our simulations strongly indicate that this perturbation does not change the critical exponents .

Fig. 4 A schematic flow diagram of renormalization group flows between different fixed points of sandpile models.

The DP fixed point is expected to be rather robust against perturbations. We have tested several variations of toppling rules in simulations. One can make the particle transfer process stochastic with each transferred particle going to a randomly chosen downward neighbor or one can allow multiple topplings at a site with each site toppling twice, thrice, etc., with decreasing probability, so long as the height is $`>1`$. With both multiple topplings and stochasticity in particle transfer, in the limit of no stickiness this becomes Manna model . For sticky grains, in all these models, we get the DP behavior. The schematic renormalization group flows are shown in Fig. 4.

To summarize, we have studied several sandpile models which show DP-like behaviour. A feature, which is common in all these models is ‘stickiness’, $`i.e.`$ with a nonzero probability a site can remain stable even with a height greater than the threshold. This behavior seems to be robust against perturbations, and is the generic behavior of sandpile models, both directed, and undirected.

We thank M. Barma, R. Dickman and A. Vespignani for their comments on the manuscript.


text:
# Linear and nonlinear marginal stability for fronts of hyperbolic reaction diffusion equations

## I Introduction

Reaction diffusion equations are used to model transport phenomena in a variety of contexts such as population dynamics, transmission lines, flame propagation among others. The prototype of such equations, for which a thorough understanding KPP37 ; AW78  of its properties has been achieved is the parabolic reaction diffusion equation

$$u_t=u_{xx}+f(u),\text{with}f(0)=f(1)=0,$$
(1)
where subscripts denote derivatives. For positive reaction terms, sufficiently localized initial conditions evolve into a monotonic decaying travelling front joining the stable $`u=1`$ to the unstable $`u=0`$ equilibrium points. For bistable reaction terms, which satisfy $`f<0`$ in $`(0,a)`$ and $`f>0`$ in $`(a,1)`$, with $`_0^1f>0`$, it is possible to find initial conditions for which the system will evolve into a monotonic decaying front joining the two stable equilibrium points $`u=1`$ to $`u=0`$. In the first case there is a continuum of speeds for which there exist monotonic fronts, the system evolves into the front of minimal speed. In the second case there is a unique speed.

The use of Eq. (1) to model a physical process involves assumptions on the stochastic process that describes the motion of the individuals. More specifically, Brownian motion is assumed mackean . If instead a more realistic process is considered Dunbar , then, in one spatial dimension, the differential equation that describes the motion is a hyperbolic reaction diffusion equation of the form

$$u_{tt}+\varphi (u)u_t=u_{xx}+f(u).$$
(2)
It has been shown that, as in the usual reaction diffusion equation, for positive $`\varphi (u)`$, the hyperbolic equation admits monotonic decaying fronts of speed $`c_h<1`$ and that the speed of the front is determined from a related parabolic equation Hadeler2 .

Equations of this form have been studied for particular cases of $`\varphi (u)`$ obtained assuming specific functions for memory effects in the diffusion term Mendez ; Manne ; Sancho ; Abramson . Recent work has dealt with the stability of the traveling front Gallay . For constant $`\varphi (u)`$ and positive, concave, $`f`$ the traveling fronts are stable to local perturbations. The asymptotic behavior of the front of minimal speed has also been found Gallay . Numerical solutions for positive but non concave $`f`$ indicate that the transition from linear to nonlinear marginal stability occurs at the same parameter values than for the parabolic equation Eq. (1)Sancho . Detailed analysis for a piecewise linear reaction term has given insight into the nature of possible travelling waves Manne . In systems with exponential transport memory $`\varphi (u)`$ is of the form

$$\varphi (u)=sf^{}(u).$$

For this case a variational principle for the front of minimal speed was constructed for positive reaction terms $`f(u)`$; it was proved that for concave reaction terms, as in the parabolic equation, linear marginal stability holds Mendez .

In this work we construct a variational principle for any positive function $`\varphi (u)`$ and general reaction terms for which the parabolic equation admits a front. From this principle upper and lower bounds are constructed which permit the a priori determination of the transition from linear to nonlinear marginal stability VS88 ; VS89 . We recover all known results as particular examples, and construct others that show different types of behavior depending on the explicit functions $`f(u)`$ and $`\varphi (u)`$. Our main result is the following, consider Eq. (2), with $`\varphi (u)>0`$, and $`f(u)`$ a reaction term for which fronts of the parabolic equation exist, as described above. The minimal (or unique, for the bistable case) speed of the fronts joining the stable to the unstable points $`u=1`$ to $`u=0`$ is given by

$$\frac{c^2}{1c^2}=\underset{g}{sup}2\alpha \frac{_0^1g(K(u))f(u)𝑑u}{_0^1g^2(u)/g^{}(u)du}$$
(3)
where

$$\frac{1}{\alpha }=_0^1\varphi (u)𝑑u$$
(4)
and

$$K(u)=\alpha _0^u\varphi (u)𝑑u.$$
(5)
The supremum is taken over all positive decaying functions $`g(u)`$ for which the integrals exist. From this result it follows that, for positive reaction terms which satisfy $`f^{}(0)>0`$, the following bound holds,

$$4\frac{f^{}(0)}{\varphi ^2(0)}\frac{c^2}{1c^2}4\underset{u}{sup}\frac{\alpha f(u)}{\varphi (u)K(u)}$$

which enables one to characterize functions for which marginal stability holds. The lower bound, the marginal stability value, is that imposed by linear considerations alone. Improved lower bounds valid for all reaction terms are obtained by direct use of the variational principle.

## II Speed of fronts

Consider the hyperbolic reaction diffusion equation

$$u_{tt}+\varphi (u)u_t=u_{xx}+f(u),$$
(6)
with

$$\varphi (u)>0\text{and}f(0)=f(1)=0,$$

where $`f(u)`$ and $`\varphi (u)`$ $`C^1[0,1]`$. We shall assume that $`f`$ belongs to the class for which monotonic fronts joining the equilibrium points $`u=1`$ to $`u=0`$ exist. The precise conditions have been spelled above in the introduction, and in more precisely elsewhere AW78 ; frentes .

We wish to find the minimal or unique speed for which there is a monotonic decaying travelling wave solution $`u=u(xc^{}t)`$ of Eq.(6). The speed satisfies $`(1c^2)u_{zz}+c\varphi (u)u_z+f(u)=0`$, $`limu_z\mathrm{}=1`$, $`limu_z\mathrm{}=0`$, $`u_z<0`$, where $`z=xct`$. It is known that when $`f>0`$ there is a continuum of fronts for a range of speeds $`c^{}<c<1`$, where $`c^{}`$ is the minimal speed of a related parabolic equation Hadeler2 . There is a unique speed in the bistable case. We may, as usual, consider the trajectory in phase space by defining $`p=u_z(u)`$. The monotonic decaying front obeys

$`(1c^2)p{\displaystyle \frac{dp}{du}}c\varphi (u)p+f(u)=0,`$ (7a)
$`p(0)=p(1)=0,p>0\text{in}(0,1).`$ (7b)
Before going any further we recall the constraints posed by linearization around the fixed points. Linearizing around $`u=0`$ we find that $`p`$ approaches zero as $`mu`$, where $`m`$ is the positive root of

$$(1c^2)m^2c\varphi (0)m+f^{}(0)=0,$$

that is,

$$m=\frac{1}{2}\frac{c\varphi (0)}{1c^2}+\frac{1}{2(1c^2)}\sqrt{c^2\varphi ^2(0)\frac{4f^{}(0)}{1c^2}}.$$
(8)
When $`f^{}(0)>0`$, since it is known that monotonic fronts exist with $`c^2<1`$, $`m`$ is real if the term in square brackets is positive, that is,

$$\frac{c^2}{1c^2}\frac{4f^{}(0)}{\varphi ^2(0)}\frac{c_L^2}{1c_L^2},$$
(9)
where we have called $`c_L`$ the minimal speed derived from linear theory or, as it is known, linear marginal stability value.

When $`f^{}(0)0`$, linear theory imposes no constraint on the speed. Analysis near the equilibrium point $`u=1`$ imposes no constarint on the speed as we assume that $`f`$ is positive or of the bistable type.

The simplest, but not unique, method to obtain a variational characterization of the speed is to introduce a stretching of coordinates which reduces the equation to the standard parabolic reaction diffusion equation Hadeler2 . Since $`\varphi (u)>0`$, we may introduce a new independent coordinate defined by the following transformation,

$$y=K(u)$$
(10)
with $`K(u)`$ defined in Eq.(5). The transformation is invertible, and $`y`$ varies between 0 and 1. In the new coordinates, Eq.(7) reads

$$p(y)\frac{dp}{dy}\frac{c}{\alpha (1c^2)}p(y)+\frac{f(K^1(y))}{\alpha (1c^2)\varphi (K^1(y))}=0,$$
(11)
$$p(0)=p(1)=0,p>0\text{in}(0,1).$$

This equation is the equation for fronts of the parabolic equation, of speed

$$\widehat{c}=\frac{c}{\alpha (1c^2)}$$

and a reaction term

$$F(y)=\frac{f(K^1(y))}{\alpha (1c^2)\varphi (K^1(y))}.$$

We now check that $`F(y)`$ satisfies all the requirements of the existence of fronts, if $`f`$ does. Effectively, $`F(0)=F(1)=0`$; for $`\varphi >0`$ and $`1c^2>0`$, sign(F) = sign(f), and

$$_0^1F(y)𝑑y=\frac{1}{1c^2}_0^1f(u)𝑑u.$$

So that if $`f`$ is of the bistable type, for which fronts exist, the same holds for $`F`$ if $`1c^2>0`$ which we know to hold. Moreover

$$F^{}(y=0)=\frac{1}{\alpha ^2(1c^2)}f^{}(0).$$

We may then apply directly the variational principle for fronts of the parabolic reaction diffusion equation frentes , that is,

$$\widehat{c}^2=\underset{g}{\mathrm{max}}\frac{2_0^1g(y)F(y)𝑑y}{_0^1(g^2(y)/g^{}(y))𝑑y}$$
(12)
where the maximum is taken over all positive decaying functions $`g`$ in (0,1) for which the integrals exist. The maximum is attained for

$$g=\widehat{g}=\mathrm{exp}\left(_{y_0}^y\frac{\widehat{c}}{p(s)}𝑑s\right),0<y_0<1$$
(13)
provided that $`F^{}(0)0`$, or, if $`F^{}(0)>0`$, when $`\widehat{m}/\widehat{c}>1/2`$. This is,

$$\frac{\widehat{m}}{\widehat{c}}>\frac{1}{2}[1+\sqrt{\widehat{c}^24F^{}(0)}].$$

Written in terms of the original quantities, $`c,f,\varphi `$, this condition reads

$$\frac{c^2}{1c^2}>\frac{4f^{}(0)}{\varphi ^2(0)}.$$

As in the usual parabolic case the maximum is attained except at the linear marginal stability value. Here too one can show that the linear marginal stability value is obtained by taking the supremum instead of the maximum. In terms of the original functions $`f`$ and $`\varphi `$ , Eq.(12) reads,

$$\frac{c^2}{1c^2}=\underset{g}{\mathrm{max}}\frac{2\alpha _0^1(g(y)f(K^1(y))/\varphi (K^1(y)))𝑑y}{_0^1(g^2(y)/g^{}(y))𝑑y}.$$
(14)
Notice that the simpler form Eq.(3) is obtained by changing to the original independent variable $`u`$.

## III Criterion for Linear and Nonlinear Marginal Stability

For the parabolic reaction diffusion equation Eq(1), for $`f>0`$, we know that the minimal speed, $`c_{PF}`$ of the propagating front is bounded below and above. When these two bounds coincide, then one can determine the speed unambiguosly. A first estimate comes from the bound KPP37 ; AW78

$$2\sqrt{f^{}(0)}c_{PF}2\underset{u}{sup}\sqrt{\frac{f(u)}{u}}.$$

When these two bounds coincide linear marginal stability holds. It may still hold if these two bounds do not coincide. This can be decided making use of the integral variational principle which improves the lower bound frentes , and of a minimax variational principle which improves the upper bound HR75 . Combined use of the two permits the exact determination of the speed. For the bistable case, we can obtain the speed making use of the integral variational principle. In this section the analogue of the above results is obtained for the hyperbolic fronts.

### III.1 Upper bound

We present here a simple derivation of the upper bound, valid both for positive and bistable reaction terms. The following evident inequality holds,

$$\frac{c^2}{1c^2}=\underset{g}{\mathrm{max}}\frac{2\alpha _0^1\frac{g(y)f(K^1(y))}{\varphi (K^1(y))}𝑑y}{_0^1(g^2(y)/g^{}(y))𝑑y}$$

$$\underset{g}{\mathrm{max}}\underset{y}{sup}\left[\frac{2\alpha f(K^1(y))}{y\varphi (K^1(y))}\right]\frac{_0^1yg(y)𝑑y}{_0^1g^2(y)/g^{}(y)dy}.$$

But, as shown in the Appendix,

$$_0^1yg(y)𝑑y2_0^1g^2(y)/g^{}(y)dy,$$

from where we obtain

$$\frac{c^2}{1c^2}\underset{y}{sup}\left[\frac{4\alpha f(K^1(y))}{y\varphi (K^1(y))}\right].$$
(15)
Going back to the original independent variable, this is

$$\frac{c^2}{1c^2}\underset{u}{sup}\left[\frac{4\alpha f(u)}{K(u)\varphi (u)}\right].$$
(16)
Observe that this upper bound coincides with the linear value when the supremum occurs at $`u=0`$.

### III.2 Lower Bounds

The lower bounds have already been obtained, one is the linear bound Eq.(9), and improved lower bounds can be obtained by choosing particular trial functions $`g`$. Equality, in the expression below, holds for a certain $`g`$.

$$\frac{c^2}{1c^2}\mathrm{max}(\frac{c_L^2}{1c_L^2},\frac{2\alpha _0^1\frac{g(y)f(K^1(y))}{\varphi (K^1(y))}𝑑y}{_0^1(g^2(y)/g^{}(y))𝑑y}).$$
(17)
## IV Applications

### IV.1 Constant $`\varphi (u)`$

As a first example consider the case $`\varphi (u)=\varphi _0`$, a constant. As explained in the introduction the stability of the traveling front for this equation has been proven recently Gallay . Numerical investigations have shown that the transition from linear to nonlinear marginal stability occur at the same parameter values as for the parabolic equation with the same reaction term Sancho . This last fact follows directly from the bounds above. For a constant $`\varphi (u)`$, we obtain $`\alpha =1/\varphi _0`$ and $`y=u`$. Then, the bounds read

$$4f^{}(0)\frac{\varphi _0^2c^2}{1c^2}4\underset{u}{sup}\frac{f(u)}{u}.$$

The upper and lower bounds are identical to those for the speed of fronts of the parabolic equation parabolic equation; the speed is given by

$$\frac{\varphi _0^2c^2}{1c^2}=\underset{g}{\mathrm{max}}\frac{2_0^1g(u)f(u)𝑑u}{_0^1(g^2(u)/g^{}(u))𝑑u}=C_{PF}^2,$$
(18)
The speed of the front is determined in terms of the speed of the parabolic front with the same reaction term, therefore the transition from linear to nonlinear marginal stability occurs at the same parameter values as for the parabolic case. The speed itself, however, is lower than the speed of the front of the parabolic equation,

$$c^2=\frac{c_{PF}^2}{\varphi _0^2+c_{PF}^2}<c_{PF}^2.$$
(19)
### IV.2 Suppresion of nonlinear marginal stability

As a second example, which shows that substantially different behavior may occur, take the reaction function $`f(u)=u(1u)(1+au)`$. The front of minimal speed for the parabolic equation is given by

$`C_{PF}`$ $`=`$ $`2\text{for}0<a<2`$
$`C_{PF}`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{a}}}+\sqrt{{\displaystyle \frac{a}{2}}}\text{for}a>2.`$
If we now consider the hyperbolic equation with $`\varphi (u)=1+au`$, the upper and lower bounds coincide and the speed

$$\frac{c^2}{1c^2}=4\text{for all}a.$$

Effectively we obtain

$$4=4f^{}(0)\frac{c^2}{1c^2}4\underset{u}{sup}\frac{1u}{1+au/2}=4.$$

The minimal speed is the linear marginal stability value, no transition from linear to nonlinear marginal stability occurs. A significant slowdown of the speed takes place.

### IV.3 Time delayed diffusion

In systems with time delayed diffusion the function $`\varphi (u)`$ adopts the form

$$\varphi (u)=sf^{}(u)$$

form which arises when considering diffusion with exponential time delay. Here $`s`$ is a parameter. For this function $`\alpha =1/s`$ and $`K(u)=uf(u)/s`$. The speed of the front was studied for positive reaction terms, and it was shown that for concave reaction terms marginal stability holds Mendez . This can be proved easily from the present approach. As Eq.(16) shows, linear marginal stability holds for all positive reaction terms with $`f^{}(0)>0`$ provided that the supremum of

$$B(u)=\frac{4\alpha f(u)}{K(u)\varphi (u)}$$
(20)
occurs on $`u=0`$, since $`B(0)=4f^{}(0)/\varphi ^2(0)`$, is equal to the lower bound. To prove then that linear marginal stability holds it suffices to show that $`B^{}(u)<0`$ for $`u(0,1]`$. Taking the derivative we obtain, after some algebra,

$$B^{}(u)=4\frac{\alpha }{K(u)}\left[\frac{ff^{\prime \prime }}{\varphi ^2(u)}+\frac{s\alpha }{\varphi (u)K(u)}h(u)\right]$$

where $`h(u)=uf^{}f`$. It follows immediately that for concave reaction terms $`B^{}(u)<0`$. For concave terms $`f^{\prime \prime }<0`$, so the first term in the square brackets is negative. And $`h(u)`$ is negative for concave functions as well, $`h(0)=0`$, and $`h^{}=uf^{\prime \prime }<0`$. So we have recovered the result that for concave functions with $`f^{}(0)>0`$, linear marginal stability holds. As it occurs for the parabolic case. it may hold for noncave functions as well, fact which can be determined by use of the variational principle.

Consider again the reation term $`f(u)=u(1u)(1+au)`$ as in subsection B. In the parabolic equation transition from linear to nonlinear marginal stability occurs at $`a=2`$. Choose here $`s=2`$ which guarantees $`\varphi >0`$ for $`0<a<(5+\sqrt{21})/24.79`$ . We show in Fig. 1 that the transition from linear to nonlinear marginal stability occurs at a lower value of $`a`$ than that for the parabolic case. We chose as a trial function $`g(u)=(1u)^7/u^{1.8}`$. The dot-dash line shows the lower bound imposed by the linear marginal stability value, the solid line the lower bound from the variational principle. The transition occurs at at $`a`$ as low as 1.6 at least. The precise transition value can be found by using a better trial function. The dashed line shows the upper bound, which guarantees that linear marginal stability holds for $`0<a<1`$, values for which $`f`$ is concave.

### IV.4 A bistable reaction term

The variational principle and the bounds obtained above hold for positive and bistable reaction systems, no assumption on the sign of $`f`$ is made in the derivations. We may therefore apply them to a bistable reaction term. We take as an example a case for the speed of the parabolic front can be determined exactly,

$$f(u)=u(1u)(ua)\text{for}0<a<1/2$$

The unique value of the speed for which a front exists is

$$c_{PF}=\frac{1}{\sqrt{2}}a\sqrt{2}$$

If $`\varphi (u)=\varphi _0`$ a constant, the speed of the front is given by Eq(19), and the exact speed is known with $`c_{PF}`$ given above.

For this bistable reaction term, with $`\varphi (u)=1f^{}(u)`$, which is positive for $`0<a<1/2`$, with a simple trial function $`g(u)=1u`$, all integrals are elementary, we obtain, using Eq.(14),

$$c^2\frac{1849a+14a^2}{8849a+14a^2}.$$

Numerical determination of the speed can be made with any desired accuracy by use of improved trial functions.

We wish to point that, in recent work, a variational principle valid only for positive reaction terms has been applied to a bistable example Mendez . Bistable reaction terms, by definition, are not positive throughout the whole interval, so the hypothesis of their derivation is violated. We attribute the nearly perfect agreement between the speed of the front obtained from their (nonapplicable) variational principle to that obtained from direct numerical integration to a coincidence for the specific trial function chosen. It is not difficult to find an acceptable trial function, sharply peaked at the origin, for which the integrand of their variational expression, and so the speed, becomes complex.

## V Conclusion

We have studied the speed of travelling fronts of hyperbolic reaction diffusion equations (6) for reaction terms for which the parabolic reaction diffusion equation (1) admits fronts. A variational principle for the speed of the fronts has been obtained for a wide class of systems. Explicit upper and lower bounds, obtained from the variational expression allow one to characterize sytems for which marginal stability holds. We have seen that depending on the nonlinearities of the equation the range of validity of linear marginal stability, can be increased or decreased, compared to the comparable range for a parabolic equation with the same reaction term. The stability of the fronts has been proved for a limited class of such equations Gallay , numerical solutions have been explored for systems with memory effect. These results will be more relevant if the stability of fronts of more general systems is proved. At present numerical results indicate that the front of minimal or unique speed can be selected. The general initial value problem for hyperbolic reaction diffusion equations is an open problem.

## Appendix A

Let $`g(y)`$ be a positive decaying function, which satisfies $`g(1)=0`$, and call $`h(y)=g^{}(y)`$, so that $`h>0`$. Then

$`\left[{\displaystyle _0^1}yg(y)𝑑u\right]^2`$ $`=`$ $`\left[{\displaystyle _0^1}\left({\displaystyle \frac{g}{\sqrt{h}}}\right)(\sqrt{h}y)𝑑u\right]^2`$
$``$ $`{\displaystyle _0^1}{\displaystyle \frac{g^2(y)}{h(y)}}𝑑y{\displaystyle _0^1}y^2h(y)𝑑y,`$
where we used Schwarz’s inequality. The second integral on the right side is, integrating by parts,

$$_0^1y^2h(y)𝑑y=_0^1y^2(g^{}(y))𝑑y=2_0^1yg(y)𝑑y,$$

where we used the fact that $`g`$ diverges at $`y=0`$ slower than $`(1/y^2)`$ (this is equivalent to $`m/c>1/2`$) and that $`g(1)=0`$ which can be seen form Eq.(13). Replacing this in the inequality above we have the desired result,

$$_0^1yg(y)𝑑y2_0^1g^2(y)/h(y)𝑑y.$$


text:
# Near-Infrared Adaptive Optics Spectroscopy of Binary Brown Dwarf HD 130948B and C1footnote 11footnote 1Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.

## 1 Introduction

HD 130948 is a triple system discovered by Potter et al. (2002) using the University of Hawaii adaptive optics system (Hōkūpa‘a) in the course of an imaging survey of low-mass companions to nearby young stars using the Gemini North telescope. The HD 130948 (Gliese 564, HR 5534) system consists of a G2V primary star and a binary companion at a separation of 2$`\stackrel{}{\mathrm{.}}`$63 from the primary. The binary companions are separated by 0$`\stackrel{}{\mathrm{.}}`$134, or 2.4 AU at a distance of 17.9 pc determined by Hipparcos. Potter et al. (2002) confirmed the common motion of the members of the triple system with a baseline of 7 months to reject the chance projection of background sources. By using the AO system on the Subaru Telescope, we obtained individual spectra of each component to confirm the substellar nature of the faint companions, HD 130948B and C.

## 2 Observations and Data Reduction

The spectroscopic observation was made on UT 2001 May 3 using the Infrared Camera and Spectrograph (IRCS; Tokunaga et al. 1998; Kobayashi et al. 2000) with the 8.2 m Subaru Telescope in conjunction with its adaptive optics (AO) system (Takami et al., 1998; Gaessler et al., 2002). The Subaru AO system, which is equipped with a 36-element curvature sensor, is installed at the front end of the telescope Cassegrain port. A medium-resolution grism was used with a 0$`\stackrel{}{\mathrm{.}}`$10 slit in the 22 mas camera section of the IRCS to provide spectra with a resolving power of 800–1000 in the $`H`$ and $`K`$ bands. HD 130948A was used as the reference source for the AO system. The slit was put along the position angle of HD 130948B and C to cover them simultaneously. The spectra were recorded by nodding the telescope by 1″ along the slit to subtract the sky emission. The total on-source integration time was 1200 s in both the $`H`$ and $`K`$ bands. A nearby B8V star HD 164352 was observed as a spectroscopic standard at similar airmass. The standard star is bright enough to give nearly the same AO correction as for HD 130948A. The spectroscopic flat field was obtained at the end of the night with a halogen lamp. The seeing was 0$`\stackrel{}{\mathrm{.}}`$4–0$`\stackrel{}{\mathrm{.}}`$5 at 2.2 µm throughout the observing period.

Because of the close proximity and the large difference in luminosity between the primary and companions, the sky background in the spectrum is affected by the gradient in the wing of the point-spread function (PSF) of the bright primary. After a pair subtraction of the sky the residual background slope was removed by interpolating the background on both sides of the companion.

We obtained the one-dimensional spectra using the IRAF333IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. aperture extraction package after flat-fielding and bad-pixel correction. Special attention was paid to the adjustment of the aperture width. Figure Near-Infrared Adaptive Optics Spectroscopy of  Binary Brown Dwarf HD 130948B and C111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. shows crosscuts of the observed spectrogram in the spatial direction with an overlay of the standard star profile normalized at the peak of each object. The FWHM of a standard star is as sharp as 0$`\stackrel{}{\mathrm{.}}`$1 and 0$`\stackrel{}{\mathrm{.}}`$08 in the $`H`$ and $`K`$ bands, respectively. However, the intensity of the first diffraction ring is significant. The separation between the centroid and the peak of the first diffraction ring for an ideal diffraction-limited image of the 8.2 m circular opening should be 0$`\stackrel{}{\mathrm{.}}`$07 and 0$`\stackrel{}{\mathrm{.}}`$09 in the $`H`$ and $`K`$ bands, respectively. This is similar to the angular separation of 0$`\stackrel{}{\mathrm{.}}`$13 between HD 130948B and C. Consequently, the energy of HD 130948B in the first ring as well as the halo overlaps with HD 130948C. This may produce similar spectral signatures in both components that may be misleading.

We calculated the contamination expected with the standard star PSF. The flux contamination of C by B within the 3 pixel (0$`\stackrel{}{\mathrm{.}}`$067) aperture width is 25% and 22% in the $`H`$ and $`K`$ bands, respectively. The contamination of B by C is $`<`$20%. The extraction aperture width of the standard star is set to the same as that of the objects. To provide the correct weight to the optimal extraction method (Horne, 1986), the statistical error in the 2D spectrogram image is calculated by combining the readout noise and the shot noise from the sky emission and PSF spillover from the primary. The wavelength calibration was performed by fitting about 20 emission lines of the argon arc lamp of the Subaru Telescope calibration unit with a linear function.

The intrinsic stellar lines in the standard star spectra, which are mostly atomic hydrogen lines, were difficult to remove because of their superposition with the telluric absorption lines. Only six of them at 1.588 µm (Br10), 1.641 µm (Br8), 1.681 µm (Br7), 1.737 µm (Br6), 2.166 µm (Br$`\gamma `$), and 2.374 µm (Pf20) were fit with a Lorentzian profile and subtracted before division. The small discrepancy in the airmass between the standard and the object was corrected by rescaling the standard star spectra according to Beer’s law. Rough flux calibration relative to the standard star was made in each band by assuming the intrinsic spectra of a B8V star is represented by a Planck function of $`T_{\mathrm{eff}}=`$ 10700 K (Tokunaga, 2000).

## 3 Result and Discussion

### 3.1 Spectral Type of HD 130948B and C

To determine the spectral types of HD 130948B and C, we compared the object spectra with template spectra of low-mass stars published by Reid et al. (2001), Geballe et al. (2002), Leggett et al. (2000), and Leggett et al. (2001). We found none of the template spectra match well with the object spectra. The continua of the object spectra were too steep to be a low-mass star, though the low-mass nature is apparent in the strong absorption of H2O at both ends of the $`H`$ and $`K`$ bands. We concluded that the continuum slope was affected by the AO correction. Not only a little offset from a narrow slit, but also slight difference in seeing can result in the variation of the continuum slope because of the sensitive wavelength dependency of the encircled energy on the AO performance (M. Goto et al., in preparation). The variation of the continuum slope from exposure to exposure was sometimes as much as 10% per 0.1 $`\mu `$m. Most of the time the difference was well represented by a simple linear function.

Because of the uncertainty of the continuum slope, we had to correct the object spectra by multiplying a linear function before comparison with templates. The strength of H2O absorption increases with spectral type, and the spectral change between any two spectral types is not a linear function. In other words, we cannot convert the spectrum of an L1 star to an L8 star by just multiplying a linear function. Therefore, in principle, we can determine the spectral type of an object uniquely by matching both the H2O absorption with a template spectrum, even if we have a linear continuum slope uncertainty.

First, we assumed HD 130948B (or C) is a certain spectral type, and estimated how much linear continuum correction is required to best match the object spectrum with the template spectrum of the assumed spectral type. We calculated the cross-correlation between the object spectrum and the template spectrum as a function of the steepness of the applied linear function correction. The slope of the function is determined so that the cross-correlation between the corrected object spectrum and the template spectrum is highest. We repeated this procedure for each spectral type from L1 to L8.

Second, we investigated how well it matched each spectral type from L1 to L8. The value of the cross correlation calculated above peaks at L3–L5 in the $`H`$ band and at L1–L5 in the $`K`$ band. Considering the simultaneous match in the $`H`$ and $`K`$ bands, we conclude both HD 130498B and C agree best with the template for L4 $`\pm `$ 1. This is consistent with the spectral type of L2 $`\pm `$ 2 determined independently by Potter et al. (2002).

The final object spectra after correction are shown in Figure Near-Infrared Adaptive Optics Spectroscopy of  Binary Brown Dwarf HD 130948B and C111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. with overlay of the best-matched template spectra, 2MASSW J0036+18 (L4). The spectra were binned by 2 pixels along the dispersion direction. The spectra of HD 130948B and C turn out to be very similar to each other in detail and are almost twins, except that HD 130948B is 0.2–0.3 mag brighter than C. The deep absorption bandhead of the CO molecule at 2.3–2.4 $`\mu `$m is conspicuous. The presence of absorption lines of Na I (2.209 $`\mu `$m) and K I (1.516 $`\mu `$m) supports the idea that these are low mass stars. The absorption band of CH4 is not seen. This is in good agreement with the conclusion that HD 130948B and C are not later than mid L-type stars. A series of the “unidentified” absorption features reported by Reid et al. (2001) at 1.58, 1.613, and 1.627 $`\mu `$m can be all attributed to the FeH molecule (Wallace & Hinkle, 2001). Ca I lines noted by Jones et al. (1994) and Tinney, Mould, & Reid (1993) are marked in the figures.

### 3.2 Mass of HD 130948B and C

The range of spectral types derived above corresponds to an effective temperature of $`T_{\mathrm{eff}}=1900\pm 75`$ (Leggett et al., 2001). On the basis of moderate X-ray activity, a fast rotation period ($``$ 7.8 days), and photospheric lithium abundance, the age of the primary should be 0.6 $`\pm `$ 0.2 Gyr (Gaidos, 1998; Gaidos, Henry, & Henry, 2000). On the other hand K. Fuhrmann (2001, private communication) finds that HD 130948A may be a member of the Ursa Major Stream, which has estimated age of 0.3 to 0.65 Gyr (Giannuzzi, 1979; Soderblom & Mayor, 1993; Palouš & Hauck, 1986; Chen et al., 1997).

Figure Near-Infrared Adaptive Optics Spectroscopy of  Binary Brown Dwarf HD 130948B and C111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. is a reproduction of the evolutionary tracks of various low-mass stars of solar metallicity calculated by Baraffe et al. (1998). Regardless of the uncertainty in the age determination of the HD 130948 system, the estimated masses for HD 130948B and C are under the sustainable hydrogen burning limit, 0.075 M for solar metallicity, for both companions, and the most likely mass of the object is 0.040–0.065 M.

## 4 Summary

We have presented moderate-resolution near-infrared spectra of the close binary companions recently discovered near HD 130948. The general spectral features of the companions, HD 130948B and C, are very similar. Based on a direct comparison with template spectra of low-mass stars, we have determined the spectral type of HD 130948B and C is L4 $`\pm `$ 1. This is consistent with the result of Potter et al. (2002), who found a spectral type of L2 $`\pm `$ 2 for the binary companions. If we take the young age of the primary star (0.3–0.8 Gyr) into account, then HD 130948B and C are binary brown dwarfs.

We thank T. Geballe, S. Leggett, I. Neill Reid, and their collaborators for making the template spectra available in machine readable form for use in Figure Near-Infrared Adaptive Optics Spectroscopy of  Binary Brown Dwarf HD 130948B and C111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.. We also thank I. Baraffe and collaborators for making the tabular data available at their ftp site for use in Figure Near-Infrared Adaptive Optics Spectroscopy of  Binary Brown Dwarf HD 130948B and C111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.. We thank all the staff and crew of the Subaru Telescope and NAOJ for their valuable assistance obtaining these data and their continuous support for IRCS construction. M. Goto is supported by a Japan Society for the Promotion of Science fellowship.