text:
# Supernova Observation Via Neutrino-Nucleus Elastic Scattering in the CLEAN Detector

## I Introduction

Rich information on neutrino properties, oscillations, the supernova mechanism, and very dense matter is contained in the neutrinos from core-collapse supernovae. Existing detectors such as Super-Kamiokande should accurately measure the $`\overline{\nu }_e`$ component of the supernova signal. However the very interesting $`\nu _\mu `$, $`\nu _\tau `$, $`\overline{\nu }_\mu `$, and $`\overline{\nu }_\tau `$ (collectively $`\nu _x`$) components may be detected without direct energy information and or in the presence of significant backgrounds from other neutrino induced reactions. Therefore, additional $`\nu _x`$ detectors could be very useful.

Perhaps the “ultimate” supernova detector involves neutrino-nucleus elastic scattering . The count rate in such a detector could be very high because the coherent elastic cross section is large and all six neutrino components ($`\nu _e`$,$`\overline{\nu }_e`$, and the four $`\nu _x`$) contribute to the signal. In particular, the detector is sensitive to $`\nu _x`$, which are expected to have a high energy and large cross section. Elastic scattering detectors can have yields of a few or more $`\nu _x`$ events per ton for a supernova at 10 kpc ($`3.1\times 10^{20}`$ m). This is an increase by a factor of 20 or more over existing light-water detector yields of hundreds of $`\overline{\nu }_e`$ and tens of $`\nu _x`$ events per kiloton.

Furthermore, the energy of nuclear recoils provides direct information on the $`\nu _x`$ spectrum. Existing detectors measure $`\nu _x`$ via neutral-current inelastic reactions on oxygen, deuterium , or carbon . Here the observed energy deposition does not depend on the neutrino energy as long as it is above threshold. Perhaps neutrino-proton elastic scattering can be detected in KamLAND. This is similar to neutrino-nucleus elastic scattering but has a smaller cross section.

Alternatively, it may be possible to detect $`\nu _x`$ using inelastic excitations of Pb. Proposals include using lead perchlorate, as suggested by Elliott, OMNIS  and LAND . Here some information on $`\nu _x`$ energies may be obtained by measuring the ratio of single- to two-neutron knockout. However, the inelastic Pb cross sections are somewhat uncertain. In contrast, neutrino-nucleus elastic cross sections can be calculated accurately with very little theoretical uncertainty.

The $`\nu _x`$ spectrum depends on how neutrinos thermalize with matter in a supernova, and is somewhat uncertain. Keil, Raffelt, and Janka have studied the effects of NN bremsstrahlung, pair annihilation, and nucleon recoil on the $`\nu _x`$ spectrum . These effects can be measured with an elastic-scattering detector.

Obtaining direct information on $`\nu _x`$ energies may be very important because the difference in energies for $`\nu _x`$ compared to $`\nu _e`$ or $`\overline{\nu }_e`$ is the primary lever arm for observing neutrino oscillations. For example, $`\nu _x\nu _e`$ oscillations could lead to high-energy $`\nu _e`$. However, deducing the oscillation probability may depend crucially on knowing how hot the $`\nu _x`$ were to begin with. Neutrino-nucleus elastic scattering itself is “flavor blind”. Therefore, the signal should be independent of neutrino oscillations (among active species). Thus elastic scattering may provide a baseline with which to characterize the supernova source. Comparing this information to other flavor-dependent information and theoretical simulations may provide the best evidence of oscillations.

The kinetic energy of the recoiling nuclei is low, typically below 100 keV. It is difficult to detect such low-energy events in the presence of radioactive backgrounds. Furthermore, scintillation signals from the nuclear recoils may be reduced by quenching because of the very high ionization density. However, recent progress in designing detectors for low-energy solar neutrinos suggests that detection may be feasible. In general, backgrounds for solar neutrinos with a low count-rate signal may be more severe than those for a supernova, where all of the events are concentrated in a few-second interval.

Cryogenic Low Energy Astrophysics with Noble gases (CLEAN) is a proposed detector for low-energy solar neutrinos based on scintillation in an ultrapure cryogenic liquid. This will detect electrons from neutrino-electron scattering at energies comparable to the recoil energy of nuclei from supernova neutrinos. In this paper, we discuss the utility of CLEAN for supernovae detection via neutrino-nucleus elastic scattering.

There is considerable interest in the direct detection of weakly interacting massive particles (WIMPS). These are expected to produce recoiling nuclei with a spectrum somewhat similar to that of supernova neutrino-nucleus elastic scattering. Again, backgrounds for WIMP detection may be larger than for supernovae because of the low WIMP count rates. Present WIMP detectors, for example CDMS, have small target masses. However, future detectors may be larger.

It is important to search for neutrinoless double-beta decay, as this can distinguish Dirac from Majorana neutrinos. Existing 76Ge experiments use multikilogram masses . The next-generation experiments such as Majorana  or Genius  may employ up to a ton of Ge. The need for good energy resolution, to tell neutrinoless from two-neutrino double-beta decay, often aids in the detection of low-energy recoils. We calculate that the largest double-beta decay experiments may soon be sensitive to galactic supernovae via elastic scattering.

Finally, micropattern gas detectors  may have a threshold low enough to detect nuclear recoils. This may allow the observation of neutrino-nucleus elastic scattering using reactor antineutrinos.

Thus the technical requirements for detecting low-energy solar neutrinos, WIMPS, double-beta decay, and supernovae via nuclear-elastic scattering may be similar. One detector or approach for low-threshold, low-background, large mass measurements may have applications in multiple areas including supernova detection via neutrino-nucleus scattering.

In section II, we fold elastic scattering cross sections with a model supernova neutrino spectrum to produce recoil spectra and yields. We consider a range of noble-gas targets from 4He to 132Xe along with 12C, 28Si, 76Ge, 114Cd, 130Te, and 208Pb. We also discuss extrapolating yields to nearby isotopes. Section III focuses on the liquid-Ne-based CLEAN detector, which appears to be very promising. A Monte Carlo simulation of backgrounds in CLEAN is presented and compared to the expected supernova signal. We discuss the large signal-to-background ratio, choice of fiducial volume, and possible detector thresholds. We conclude in section IV.

## II Supernova Signals in Various Detectors

The neutrino-nucleus elastic-scattering cross section $`d\sigma /d\mathrm{\Omega }`$ is ,

$$\frac{d\sigma }{d\mathrm{\Omega }}=\frac{G^2}{4\pi ^2}k^2(1+\mathrm{cos}\theta )\frac{Q_w^2}{4}F(Q^2)^2,$$
(1)
for a neutrino of energy $`k`$ scattering at angle $`\theta `$. The Fermi constant is $`G`$. This coherent cross section depends on the square of the weak charge $`Q_w`$

$$Q_w=N(14\mathrm{s}\mathrm{i}\mathrm{n}^2\mathrm{\Theta }_W)Z$$
(2)
of a nucleus with $`N`$ neutrons and $`Z`$ protons. The weak mixing angle is sin$`{}_{}{}^{2}\mathrm{\Theta }_{W}^{}0.231`$ . We assume a spin-zero target. Finally, the ground-state elastic form factor $`F(Q^2)`$ at momentum transfer $`Q`$,

$$Q^2=2k^2(1\mathrm{cos}\theta ),$$
(3)
is

$$F(Q^2)=\frac{1}{Q_w}d^3r\frac{\mathrm{sin}(Qr)}{Qr}[\rho _n(r)(14\mathrm{s}\mathrm{i}\mathrm{n}^2\mathrm{\Theta }_W)\rho _p(r)].$$
(4)
Here $`\rho _n(r)`$ is the neutron density and $`\rho _p(r)`$ is the proton density. The form factor is normalized $`F(Q^2=0)=1`$. We neglect a small correction from the single-nucleon form factors.

The inclusion of $`F(Q^2)`$ is crucial for heavier targets. However, we evaluate it at relatively small $`Q^2`$ so the exact form of the densities is not important. The proton density is often well constrained by measured charge densities. For simplicity we use theoretical densities from simple relativistic-mean-field calculations using the successful NL3 effective interaction . These calculations assume spherical ground states and do not include pairing corrections. The use of other densities is not expected to change our results significantly.

We now consider a simple “standard model” for the supernova-neutrino spectra, see for example . This model is close to what others have used. A total energy of $`3\times 10^{53}`$ ergs (1 erg $`=10^7`$ J) is assumed to be radiated in neutrinos from a supernova at a distance $`d`$ of 10 kpc ($`3.1\times 10^{20}`$ m). For simplicity we use Boltzmann spectra at temperatures of $`k_BT=3.5`$, 5 and 8 MeV for the $`\nu _e`$, $`\overline{\nu }_e`$ and $`\nu _x`$ components respectively. Here $`k_B`$ is the Boltzmann constant, which we set to one in the rest of the paper. The use of Fermi Dirac spectra (at zero chemical potential) should give similar results. However, neutrino-nucleus elastic scattering is sensitive to the high-energy tails in the spectra. Therefore non thermal spectra could modify our results somewhat and should be investigated in future work.

We assume equal partition in energy among the $`\nu _e`$ $`\overline{\nu }_e`$ and the four $`\nu _x`$ components. Therefore this standard supernova radiates a total of $`N_{\nu _e}=3.0\times 10^{57}`$, $`N_{\overline{\nu }_e}=2.1\times 10^{57}`$ and $`N_{\nu _x}=5.2\times 10^{57}`$ neutrinos. The time integral of the neutrino flux at Earth $`\varphi _i(k)`$ for a neutrino of energy $`k`$ is

$$\varphi _i(k)=\frac{1}{4\pi d^2}\frac{N_i}{2T_i^3}k^2\mathrm{e}^{k/T_i},$$
(5)
for $`i=\nu _e,\overline{\nu }_e`$, or $`\nu _x`$.

Microscopic simulations of supernovae suggest that equal partition of energy may be good only to $``$ 25 %. Furthermore, there is important uncertainty in the $`\nu _x`$ spectrum, with estimates of $`T_{\nu _x}`$ ranging from $`6`$ to 8 MeV. It is an important goal of elastic-scattering detectors to measure $`T_{\nu _x}`$. Therefore, the predictions of our supernova spectrum have significant uncertainties. Nevertheless, this simple model should provide order-of-magnitude estimates and may allow easy comparisons to other calculations.

The yield of recoiling nuclei with energy $`E`$ and mass $`M`$ is

$$Y(E)=2\pi N_t\mathrm{\Sigma }_{i=\nu _e,\overline{\nu }_e,\nu _x}_0^{\mathrm{}}𝑑k\varphi _i(k)_1^1𝑑\mathrm{cos}\theta \delta (E\frac{Q^2}{2M})\frac{d\sigma }{d\mathrm{\Omega }},$$
(6)
where $`N_t`$ is the total number of target atoms. For our Boltzmann spectra this integral is simple:

$$Y(E)=\frac{G^2}{\pi }\frac{Q_w^2}{4}MF^2(2ME)\left(\frac{N_t}{4\pi d^2}\right)\mathrm{\Sigma }_{i=\nu _e,\overline{\nu }_e,\nu _x}N_i(t_i+1)\mathrm{e}^{t_i},$$
(7)
with $`t_i=(ME/(2T_i^2))^{1/2}`$. For large recoil energy $`E`$, $`Y(E)`$ is proportional to

$$Y(E)F^2(2ME)\mathrm{e}^{(\frac{M}{2T_{\nu _x}^2})^{1/2}E^{1/2}}.$$
(8)
For light nuclei the high-energy tail continues to hundreds of keV and is produced by the scattering of very-high-energy $`\nu _x`$. However, for heavier nuclei the tail is sharply reduced by the nuclear form factor.

We consider first noble-liquid detectors from 4He to 132Xe and then a range of other detectors in order of increasing mass number $`A`$ from 12C to 208Pb. The yield $`Y(E)`$ from Eq. (7) is shown in Figure 1 for detectors made of 4He, 20Ne, 40Ar, 84Kr, and 132Xe. We don’t mean to imply that detectors would be feasible with all of these liquids. However, we show these nuclei to illustrate how the yield depends on $`A`$ for a broad range of $`A`$. The spectra in Figure 1 are peaked at low recoil energy $`E`$. Increasing $`A`$ raises the cross section because coherent scattering is proportional to $`N^2`$. Thus at low energies $`Y(E)`$ increases significantly with $`A`$. However, as $`A`$ increases the spectrum is strongly shifted to lower energies by the form factor and the large target mass. The energy integral of $`Y(E)`$, or total yield, is given in Table I in events per ton of detector. Also listed are events above a threshold of 5, 10, 25, or 50 keV. Finally, the average recoil energy of the nuclei is given. This average is influenced by a small number of events at high energies, while the spectrum is peaked at low energies.

The optimal choice of $`A`$ may involve a trade-off between cross section, which favors high $`A`$, and recoil energy, which favors small $`A`$. This choice may depend on the attainable threshold. Furthermore, the choice of target material depends on a host of other practical considerations including the presence of possible backgrounds from radioactive isotopes. Although 4He has a relatively small cross section it has high recoil energies. We find a total yield in Table I of 0.85 events per ton. Helium-based solar-neutrino detectors include HERON  and TPC  (or HELLAZ ).

Perhaps nuclei near $`A=40`$ give a reasonable balance between cross section and recoil energy. However, Ar may have backgrounds from radioactive 39Ar and 42Ar, while Kr may have backgrounds from 85Kr. Xenon is being used in several dark-matter, double-beta decay and solar-neutrino detectors including XMASS , XENON  and ZEPLIN. The total yield is very large, 31 events per ton. However, because 132Xe is heavy, there is a strong premium on obtaining a low threshold. Background from the $`2\nu `$ double-beta decay of 136Xe should not be a problem for a supernova detector.

Above a threshold of 25 keV the yield is a relatively slow function of $`A`$. Therefore one may have considerable freedom in the choice of target material. We consider the CLEAN liquid-Ne solar-neutrino detector  at some length. The total yield of 3.99 events per ton is dominated by 3.08 $`\nu _x`$ events with only 0.38 $`\nu _e`$ and 0.53 $`\overline{\nu }_e`$ events. One is very sensitive to the $`\nu _x`$ spectrum. For example, if the $`\nu _x`$ temperature is not the expected 8 MeV but instead is near that for $`\overline{\nu }_e`$ $`T_x=5`$ MeV, the spectrum in Figure 1 is greatly changed. This verifies that the recoil spectrum contains direct information on the $`\nu _x`$ energies. A Monte Carlo simulation of backgrounds in CLEAN is presented in section III.

We now consider a number of other targets. Detectors based on organic scintillator such as Borexino and KamLAND  have a yield of 2.50 events per ton of 12C, see Table I. With any carbon-based detector there will be backgrounds from 14C. Borexino should have a ratio 14C/ 12C of the order $`10^{18}`$ . At a concentration of $`10^{18}`$ there will be about 2 14C decays per ton during the 10 seconds of a supernova neutrino burst. This is comparable to the number of elastic recoils. However, the 14C background should have a different spectrum and can be well measured at other times. Therefore, this background may not prevent the use of carbon as a supernova detector, even if it does prevent the observation of pp solar neutrinos. Quenching may be a serious problem for organic scintillator. The amount of light produced from recoiling C ions may be much less than that for recoiling electrons .

An organic scintillator will also have events from neutrino-proton elastic scattering . We estimate a yield of about 0.33 events per ton of CH2. Because of the light proton mass these events will have a larger recoil energy. The proton elastic-scattering cross section has a theoretical uncertainty of 10 to 20 % from possible strange quark contributions to the nucleon’s axial current and spin. It would be very useful to have better laboratory measurements of neutrino-proton elastic scattering. In contrast, strange quarks are not expected to make significant contributions for supernova neutrino-nucleus elastic scattering. Indeed, there is almost no theoretical uncertainty in the neutrino-nucleus elastic cross sections.

Silicon detectors such as those described in  have a yield of 5.5 events per ton and the recoil spectrum is shown in Figure 2, while 76Ge detectors have a yield of 18.6 events per ton and an average recoil energy of 9.5 keV. The double-beta decay experiment Majorana  is proposed to have a 500 kg mass, while Genius  has a proposed mass of one ton. These detectors should have very low backgrounds and low thresholds. Therefore they should be sensitive to a supernova at 10 kpc. It is remarkable that such small target masses can yield statistics for our own galaxy comparable to those for the historic IMB and Kamiokande signals from SN1987A in the Large Magellanic Cloud.

Finally, 114Cd, 130Te and 208Pb yields are listed in Table I. The yield and spectrum for 130Te is very close to 132Xe (see Figures 1 and 2) since they both have 78 neutrons. The heavy nucleus 208Pb has a very large yield of 47.5 events per ton. However, the average recoil energy is only 2.6 keV. Backgrounds and the need for a very low threshold may make a Pb detector very difficult to build.

The nuclei in Figures 1 and 2 display a range of recoil spectra. For point nuclei there would be a single universal spectral shape, with the recoil energy decreasing and the yield increasing with increasing $`A`$. However, the different nuclear form factors modify the spectra for heavy nuclei.

Finally, we provide a simple formula to extrapolate the yields in Table I to nearby isotopes. If one ignores small changes in the form factors of nearby nuclei, then the yield will be approximately proportional to the square of an effective weak charge, $`Q_{eff}`$,

$$Q_{eff}^2=\delta _{A,odd}3g_a^2+Q_W^2,$$
(9)
with $`Q_W`$ from Eq. (2), and $`\delta _{A,odd}=1`$ for odd $`A`$ nuclei and $`\delta _{A,odd}=0`$ for even $`A`$ nuclei. This factor takes into account the axial current of the last nucleon with $`g_a=1.26`$. Because this term adds in quadrature with $`Q_W^2`$, it makes a very small contribution, except for very light systems such as the protonWe ignore the slightly different angular distribution of this term.. For example, using Eq. (9) we find that 21Ne and 22Ne have cross sections respectively 1.29 and 1.48 times that of 20Ne. Natural Ne is 0.3 % 21Ne and 8.8 % 22Ne so this will lead to a slight increase in yield over that for pure 20Ne.

## III The CLEAN detector

CLEAN, a detector concept based on liquid Ne, was originally proposed for the detection of low-energy solar neutrinos. It will also have high sensitivity to weakly interacting massive particles (WIMPs). Liquid neon has a high scintillation yield, has no long-lived radioactive isotopes, and can be easily purified using cold traps. In addition, neon is inexpensive, dense, and transparent to its own scintillation light, making it practical for use in a large self-shielding apparatus. Here we consider a CLEAN detector in which a stainless steel tank holds 200 tons of liquid neon, half of which is exposed to a wavelength shifter to convert the ultraviolet light to the visible. Inside the tank and suspended in the liquid neon are several thousand photomultipliers. A diagram of the proposed CLEAN design is shown in Figure 3.

CLEAN will also be sensitive to supernova neutrinos, detected through neutrino-nuclear scattering. In this case the entire active neon mass inside the wavelength shifter can be used, with a possible modest fiducial volume cut to reduce radioactive backgrounds.

Of prime importance for determining the sensitivity of CLEAN to neutrino-nuclear scattering is the determination of the light yield of liquid neon for nuclear recoils. Because the density of excitation in scintillator is typically higher for nuclear recoils than for electron recoils, the chemically excited species are more likely to interact, increasing the likelihood that energy will be lost through mechanisms that do not produce light. This quality is often expressed as a “quenching factor”, the ratio of light emitted for a nuclear recoil to the light emitted from an electron recoil, per unit deposition energy. While the quenching factor for liquid neon has not yet been measured, we expect it to be similar in magnitude to the quenching factor measured for liquid xenon. Recent (and widely disagreeing) measurements of the liquid xenon quenching factor are 22 % and 43 %. The amount of quenching for liquid neon should be less than that for liquid xenon, as the density (1.2) of liquid neon is less than that of liquid xenon, while the scintillation mechanism in liquid neon is qualitatively similar. For the following simulations, we assume that the quenching factor for liquid neon is 0.25. Clearly the scintillation quenching factor in liquid neon would have to be accurately measured in order to properly interpret any supernova data.

As for all neutrino detectors, a prime design consideration in CLEAN is the reduction of radioactive backgrounds. We expect that any radioactive species suspended in the liquid neon will be removed by passing the neon through charcoal or similar adsorbant; however there will remain a high rate of gamma rays entering the liquid neon after being emitted by the surrounding photomultipliers, photomultiplier support structure, wavelength shifter, and stainless steel tank containing the liquid neon.

Figure 4 shows the expected supernova neutrino recoil spectrum for a CLEAN detector with 100 tons of active liquid neon, assuming 3750 scintillation photons per MeV, 100 % efficiency for the wavelength shifter, photomultiplier coverage of 75 %, and a photomultiplier quantum efficiency of 15 %. Also shown is the expected radioactive background for 10 seconds of observing time, assuming that the combined gamma and x-ray emission is dominated by the photomultiplier glass and the wavelength shifter substrate. The simulation assumes photomultipliers 20 cm in diameter, each with mass 650 g, with 30 ng per g of U and Th and 60 $`\mu `$g per g of K in the glass. The wavelength shifter is assumed to be evaporated on quartz wafers of 1 mm thickness, with a U and Th concentration of 1 ng per g. From $`10^4`$ seconds of simulated data, we find a background of 62 $`\pm `$ 8 events in 10 seconds within the energy range of 0 to 200 keV. Here the uncertainty corresponds to the $`\pm `$ 1-sigma interval for a particular 10-second observation. Since liquid neon has no long-lived radioactive isotopes that would limit its practical threshold, the CLEAN detector could conceivably trigger on as few as two photoelectrons. The accidental coincidence rate in CLEAN will be low, as the photomultiplier dark count rate will be suppressed at liquid-neon temperature (27 K). Thus we expect that the full 100 tons of liquid neon could be viewed with a threshold of 2 photoelectrons, equivalent to a recoil energy of about 5 keV. The expected 62 background events is smaller than the expected supernova signal of 330 Ne elastic events above 5 keV in 100 tons.

The radioactive backgrounds, while small in comparison to the supernova signal, can be lowered further through position resolution, since most gamma rays will deposit their energy in the outer edges of the liquid neon sphere. Recently, Monte Carlo simulations have shown that the location of ionizing radiation events in CLEAN can be determined by analyzing the pattern of photomultiplier hits. These simulations are described in detail in an upcoming publication. Here we show only some results relevant to supernova neutrino detection. We find that a mild fiducial radius cut, leaving a mass of 70 tons, virtually eliminates gamma-ray background for the purposes of supernova neutrino detection (only 2.6 $`\pm `$ 1.6 events in the fiducial volume in 10 seconds). These results are shown in Figure 5. In the case of a fiducial volume cut, detecting a few more photoelectrons will improve the convergence of the position resolution algorithm. Currently we can analyze events that produce as few as 8 photoelectrons (recoil energy of 21 keV) if position cuts are applied. In general, as the position cut is increased or the number of photoelectrons is reduced, algorithm convergence should be carefully checked. The expected 2.6 background events are much smaller than the expected supernova signal of 140 Ne events above 21 keV in 70 tons.

The signal-to-noise ratio reported here, though already quite large, might be significantly improved by the development of photomultipliers with lower inherent radioactivity. Such photomultipliers are under investigation by the CLEAN and XENON collaborations. In addition, plans are being made to measure the quenching factor for nuclear recoils in liquid neon, as this measurement is important for the determination of the sensitivity of CLEAN to WIMP particles, as well as for supernova neutrinos.

## IV Conclusions

Detectors with reduced radioactive backgrounds may be able to study supernovae via neutrino-nucleus elastic scattering. This could provide important information on the flux and spectrum of $`\nu _x`$ ($`\nu _\mu `$ and $`\nu _\tau `$). Elastic-scattering detectors could see a few or more events per ton for a supernova at 10 kpc. This is 20 or more times the number of events per ton of existing water detectors. The CLEAN experiment, based on the detection of scintillation in liquid Ne, is a prime example of a detector that would be sensitive to the elastic scattering of supernova neutrinos. Its active mass of 100 tons may yield almost 400 $`\nu _x`$ events. In addition, many other detectors, including the largest dark-matter and double-beta decay experiments, may also be sensitive to supernova neutrinos via elastic scattering.

Observation of neutrino-nuclear elastic scattering will complement supernova signals from other detectors. Water detectors such as SNO or Super-K will detect $`\nu _x`$ without $`\nu _x`$ spectral information. This energy information could be important for neutrino oscillations. KamLAND may be able to measure $`\nu p`$ elastic scattering. This contains energy information just like $`\nu A`$ scattering. The small cross section for $`\nu p`$ scattering may be compensated by a large detector mass. This may yield only slightly smaller statistics than CLEAN. We strongly encourage development of detectors based on both $`\nu p`$ and $`\nu A`$ elastic scattering since they may have different backgrounds, thresholds, and systematic errors. Furthermore, the very large $`\nu A`$ elastic cross sections may allow even larger statistics in future detectors.

###### Acknowledgements.

We thank John Beacom, Steve Elliott, and Harry Miley for useful discussions and acknowledge financial support from DOE grant DE-FG02-87ER40365 and NSF grant 0226142. We acknowledge the hospitality of the Aspen Center for Physics, where some of this work was done.
| Table I: Yield. | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Target | Y | Y $`>`$ | Y $`>`$ | Y $`>`$ | Y $`>`$ | $`<E>`$ |
|  |  | 5 keV | 10 keV | 25 keV | 50 keV | (keV) |
| 4He | 0.85 | 0.82 | 0.79 | 0.72 | 0.62 | 240 |
| 12C | 2.5 | 2.2 | 2.0 | 1.6 | 1.1 | 83 |
| 20Ne | 4.0 | 3.3 | 2.9 | 2.0 | 1.2 | 46 |
| 28Si | 5.5 | 4.2 | 3.4 | 2.1 | 1.1 | 31 |
| 40Ar | 9.4 | 6.6 | 5.0 | 2.5 | 0.99 | 21 |
| 76Ge | 18.6 | 9.6 | 5.8 | 1.7 | 0.30 | 9.5 |
| 84Kr | 19.8 | 9.5 | 5.5 | 1.4 | 0.20 | 8.4 |
| 114Cd | 26.3 | 9.7 | 4.6 | 0.70 | 0.041 | 5.7 |
| 130Te | 31.8 | 10.1 | 4.3 | 0.47 | 0.014 | 4.8 |
| 132Xe | 31.1 | 9.8 | 4.1 | 0.43 | 0.012 | 4.8 |
| 208Pb | 47.5 | 7.3 | 1.7 | 0.022 | 0.001 | 2.6 |

Yield in events per ton for a supernova at 10 kpc assuming different target materials. Also listed is the number of events above thresholds of 5, 10, 25 or 50 keV. Finally the average recoil energy $`<E>`$ is given.


format: text

text:
# Statistics for transition of a plasma turbulence with multiple characteristic scales

## I Introduction

Transition phenomena with sudden changes of states are observed in turbulent plasmas. Since these transition phenomena like the L-H transition play crucial roles for magnetic confinement of fusion plasmas, the transition associated with formation of transport barriers is one of the main subjects of high-temperature plasma physics.

Traditional theories for transition of high-temperature plasmas are formulated in terms of averaged physical quantities Itoh and Itoh (1988); Shaing and E. C. Crume (1989). Fluctuations around the averages are ignored and transition phenomena are described deterministically.

However, high-temperature plasmas are strongly non-linear systems with huge number of degrees of freedom and hence their behavior should be chaotic and unpredictable in a deterministic way. In fact, broad distribution of critical values of parameters where transitions occur and intermittent transport called “avalanche phenomena” are observed in recent experiments ITER H-Mode Database Working Group (1994); Politzer (2000); Yoshizawa et al. (2002). Occurrence of these behaviors, which cannot be described only by averaged quantities, are considered due to strong turbulent fluctuations. It is inevitable to describe turbulent plasmas in terms of probability or ensembles like numerical forecast of weather, since the magnitudes of fluctuations around averages are of the same order as that of the averages in turbulent states.

Therefore, transition takes place as a statistical process in the presence of stochastic noise sources induced by turbulence interactions. As a generic feature, transition is expected to occur with a finite probability when a control parameter approaches the critical value.

Statistical theories for plasma turbulence have been developed and the framework to calculate the probability density function (PDF), the transition probability etc. has been made Bowman et al. (1993); Krommes (1996, 1999, 2002); Itoh and Itoh (1999, 2000); Itoh et al. (2002a); Yoshizawa et al. (2002). In the statistical theories, the time-development of the system is described by a set of differential equations with random forces, called the “Langevin equations”. All the information on the statistical properties of the system is obtained by solving the Langevin equation.

The framework has been applied to cases where only one turbulent mode is excited and the turbulence is characterized by one space-time scale Bowman et al. (1993); Krommes (1996, 1999, 2002); Itoh and Itoh (2000); Itoh et al. (2002b); Kawasaki et al. (2002a, b).

However, it is well known that there are many kinds of turbulent fluctuations in high-temperature plasmas and that different characteristic length scales coexist. The importance of interactions between modes with different scale lengths has recently been recognized. For instance, the dynamics of the meso-scale structure of the radial electric field Biglari et al. (1990); Itoh et al. (1999); Terry (2000) is known to cause variation in the dynamics of microscopic fluctuations like in the electric filed domain interface Itoh et al. (1991); Diamond et al. (1997), zonal flow Smolyakov and Diamond (1999) and streamer Drake et al. (1995). Coexistence of multiple scale turbulence has also been investigated by use of the direct numerical simulations Yagi et al. (2002); Kishimoto et al. (2002). Statistical theory on zonal flow dynamics Krommes and Kim (2000) and that of the L-H transition theory has been developed Itoh et al. (2002b).

In the present paper, we apply the statistical theoretic algorithm to a model of high-temperature plasma where two characteristic scales coexist. One is the current diffusive interchange mode (CDIM) micro turbulence Itoh et al. (1999), whose characteristic length scale is of the order of the collisionless skin depth $`\delta =c/\omega _p`$. The other is the ion-temperature gradient (ITG) mode turbulence, whose characteristic wave length is of the order of the ion gyroradius $`\rho _i`$, as an example of the drift wave fluctuations considered to dictate a considerable part of the turbulent transport Horton (1999). Hereafter, we call these two modes “the high wave number mode” and “the low wave number mode” respectively for its simplicity. We assume that the condition $`\rho _i\delta `$ holds. Both turbulences are considered to cause the anomalous transport and hence the coexistence of the high wave number mode and the low wave number mode turbulences and their interplay should be taken into account.

It is known that the subcritical transition occurs in this system, when the pressure-gradient and the radial electric field are changed Itoh and Itoh (2001). However, the turbulent fluctuations are ignored in the analysis. In the present paper, with the statistical theory, we analyze the effect of the turbulent fluctuations on stochastic properties of the transition and show that the fluctuation changes the phase structure of the system completely. More precisely, we show that the transition is smeared out by noises, i.e., the physical quantities changes gradually without clear transition.

The present paper is organized as follows; the statistical theory and the model are formulated in Sec. II. The results of deterministic analyses including occurrence of the subcritical bifurcation are summarized in Sec. III. The Section IV, i.e., the main part of this paper, presents the statistical properties of the system. Summary and discussions are given in Sec. V.

## II Theoretical framework and the model

In this section, we briefly review the theoretical framework and the model of turbulent plasmas where two different characteristic scales coexist Itoh and Itoh (2001).

The starting point of the theory is the reduced MHD for the three fields: the electro-static potential, the current and the pressure. The Langevin equation that the statistical theory is based upon is derived as one of model equations which reproduce the two-time correlation functions and the response functions obtained by the renormalization perturbation theory (the direct-interaction approximation) for the reduced MHD Krommes (2002).

The Langevin equation describes time-development of two variables characterizing the system. One is the spectrum amplitude of the electro-static potential for the characteristic wave number $`k_h`$ of the high wave number mode, $`x_hk_h^2\varphi _{k_h}^{}\varphi _{k_h}/D_h^2`$, and the other is that for the low wave number mode, $`x_lk_l^2\varphi _{k_l}^{}\varphi _{k_l}/D_l^2`$, whose characteristic wave number is $`k_l`$. Here, $`D_h`$ and $`D_l`$ denote the renormalized transport coefficients for the high wave number mode and the low wave number mode respectively when there is no interactions between two modes. The characteristic time constants for the two modes are defined as

$$\omega _lk_l^2D_l,\omega _hk_h^2D_h.$$
(1)
See Itoh and Itoh (2001) for the details of the notation.

The Langevin equation gives the time-development of these two variables as

$`{\displaystyle \frac{dx_l}{d(\omega _lt)}}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\sqrt{x_h}}{r}}+\sqrt{{\displaystyle \frac{x_h}{r^2}}+4x_l}2\right)x_l+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\sqrt{x_h}}{r}}+\sqrt{{\displaystyle \frac{x_h}{r^2}}+4x_l}\right)x_lw_l=0,`$ (2)
$`{\displaystyle \frac{dx_h}{d(\omega _ht)}}+\left[\sqrt{x_h}\sqrt{1+{\displaystyle \frac{\sqrt{x_l}}{\sqrt{x_h}/r+\sqrt{x_h/r^2+4x_l}}}}/(1+pr^2x_l)\right]x_h+x_h^{3/2}w_h+\sqrt{ϵ}w_t=0.`$ (3)
Nonlinear drags and magnitudes of the noises have been evaluated in Itoh and Itoh (2001).

In these equations, the nonlinear terms of the reduced MHD are divided into two parts: One part is coherent with the test field and is renormalized into the deterministic terms, i.e., the second terms. The other is incoherent and is modeled by random forces $`w_h,w_l`$. Another random force $`w_t`$ denotes the thermal noise and the magnitude of the noise, $`ϵ`$, is a small quantity compared to the magnitudes of other noises. Technically, introduction of $`ϵ`$ is needed in order to exclude singularity at $`x_h=0`$ of the PDF.

Since $`w_h,w_l`$ and $`w_t`$ are forces which fluctuates randomly in time, the Langevin equation describes the stochastic time-development of the turbulent fluctuation of the system. For simplicity, we assume that the random forces are Gaussian and white;

$$w(t)w(t^{})=\delta (tt^{}).$$
(4)
This system is sustained by the space inhomogeneities, the curvature of the magnetic field $`\mathrm{\Omega }^{}`$, the pressure gradient $`dP_0/dx`$ and the gradient of the radial electric field $`dE_r/dr`$. Here, the shear of the magnetic field is given as the slab configuration: $`𝐁=(0,B_0sx,B_0)`$ where $`B_0(x)=\text{const}\times (1+\mathrm{\Omega }^{}x+\mathrm{})`$. The pressure is assumed to change in $`x`$direction. These driving forces are characterized by the parameter $`r`$ and $`p`$ as $`rD_l/D_h(\mathrm{\Omega }^{}dP_0/dx)^{1/2},pD_h^2/I_{eff}^{hl}(dE_r/dr)^6\mathrm{\Omega }^{}dP_0/dx`$. Here, $`I_{eff}^{hl}`$ denotes the critical strength of the nonlinear interactions between the low and high wave number modes. The pressure gradient controls the growth rate of the low wave number mode and $`\mathrm{\Omega }^{}dP_0/dx`$ excites both the high and low wave number mode turbulences. The gradient of the radial electric field suppresses turbulences Biglari et al. (1990); Itoh et al. (1999); Terry (2000).

We assume that the relation $`\rho _i\delta `$ holds and hence the characteristic length-scales for the two modes are widely separated as

$$k_hk_l.$$
(5)
The mutual interactions between the low and high wave number modes are asymmetric, since the spatial structure of the low wave number mode is a large-scale inhomogeneity for the high wave number mode. The assertion Eq. (5) also means that the time-scales are widely separated since the time scales are given by Eq. (1);

$$\omega _h\omega _l.$$
(6)
By analyzing the Langevin equation, Eqs. (2, 3), a number of statistical properties of turbulent plasmas can be derived. For example, the analytical formulae of the rate of change of states of plasmas, the transition rates, were derived. Furthermore, since the renormalized transport coefficients and the random forces have the same origin, i.e., nonlinear interactions in MHD turbulence, relations between the fluctuation levels of turbulence and the transport coefficients like the viscosity and the diffusivity were derived.

## III Bifurcation without random forces

With the theoretical framework briefly described in the previous section, we analyze the model of the inhomogeneous plasma with two characteristic scales.

At first, we show the steady state solutions when random forces are ignored, in order to compare to the results with the random forces obtained later. The steady state solutions are obtained by solving the set of nonlinear equations when the random forces are turned off:

$`\left({\displaystyle \frac{\sqrt{x_h}}{r}}+\sqrt{{\displaystyle \frac{x_h}{r^2}}+4x_l}2\right)x_l=0,`$ (7)
$`\left[\sqrt{x_h}\sqrt{1+{\displaystyle \frac{\sqrt{x_l}}{\sqrt{x_h}/r+\sqrt{x_h/r^2+4x_l}}}}/(1+pr^2x_l)\right]x_h=0.`$ (8)
Figure 1 shows the $`r`$-dependence when $`p=30`$. When $`r0.49\mathrm{}`$, the low wave number mode turbulence is suppressed. As $`r`$ is increased, the system experiences the subcritical transition to the state where the low wave number mode turbulence is excited. When $`0.49\mathrm{}<r<1`$, there are two stable solutions and it means that the system is bi-stable. From the deterministic point of view, the transition is expected to occur only at the ridge point and the bifurcation point. The qualitative behavior does not depend on the value of $`p`$ as far as $`p>1.9\mathrm{}`$.

The phase diagram is shown in Fig. 2. The subcritical bifurcation is observed when the value of the parameter $`p`$ is larger than $`1.9\mathrm{}`$. On the other hand, when $`p1.9\mathrm{}`$, the bifurcation is supercritical.

## IV The stochastic properties

In the rest of the present paper, we analyze the stochastic properties of the model, Eqs. (2, 3), to investigate the effect of the turbulent fluctuations, i.e., the random forces.

### IV.1 The adiabatic approximation and the Fokker-Planck equation

At first, we approximate the Langevin equation Eqs. (2, 3) with making use of the time scale separation between the low and high wave number mode turbulences.

The scale separation, Eq. (6), means that the high wave number mode $`x_h`$ quickly relaxes to the steady state determined by the value of the low wave number mode variable $`x_l`$ which is fixed at the value $`x_l=x_l(t)`$ at the time.

We analyze the steady state of the high wave number mode $`x_h`$ when $`x_l`$ is fixed. A state of a stochastic system is described by the probability that a state variable takes a certain value. The time-development of the probability density function (PDF) $`P(x_h,t)`$ of $`x_h`$ is determined by Kramers-Moyal expansion applied to the Langevin equation of $`x_h`$, Eq. (3). The Kramers-Moyal expansion is given as

$$\frac{P(x_h,t)}{t}=\frac{}{x_h}C_1(x_h)P(x_h,t)+\frac{1}{2}\frac{^2}{x_h^2}C_2(x_h)P(x_h,t).$$
(9)
The expansion can be truncated at the second order, since the random forces in the Langevin equation are assumed to be Gaussian. Here, the coefficient $`C_n(x_h)`$ is given by

$$C_n(x_h)\underset{\mathrm{}t0}{lim}\frac{1}{\mathrm{}t}[x_h(t+\mathrm{}t)x_h(t)]^n,$$
(10)
where $``$ denotes the average over the all realizations of the random forces and the average is taken under the condition $`x_h(t)=x_h`$. The resulting equation of motion of the probability density function is called the Fokker-Planck equation and is written as

$`{\displaystyle \frac{P(x_h,t)}{(\omega _ht)}}={\displaystyle \frac{}{x_h}}\left[\sqrt{x_h}\sqrt{1+{\displaystyle \frac{\sqrt{x_l}}{\sqrt{x_h}/r+\sqrt{x_h/r^2+4x_l}}}}/(1+pr^2x_l)\right]`$
$`\times x_hP(x_h,t)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{x_h^2}}(x_h^3+ϵ)P(x_h,t).`$ (11)
The steady state solution of Eq. (11) when $`x_l`$ is fixed is shown in Fig. 3. It is seen that the peak of the PDF is relatively narrow and it means that the high wave number mode turbulence spends most of time at the peak value $`x_h^{}`$. Hence, we can say that it is good approximation to replace $`x_h`$ in the Langevin equation of the low wave number mode with the peak value $`x_h^{}`$.

The equation to determine the location of the peak $`x_h^{}`$ is given by the condition $`dP(x_h^{})/dx_h^{}=0`$ as

$$\sqrt{1+\frac{\sqrt{x_l}}{\sqrt{x_h^{}}/r+\sqrt{x_h^{}/r^2+4x_l}}}/(1+pr^2x_l)\frac{3}{2}x_h^{}\sqrt{x_h^{}}=0.$$
(12)
It is important to note that Eq. (12) is essentially different from Eq. (8) in existence of the second term of Eq. (12), which comes from the random force $`x_h^{3/2}w_h`$ of Eq. (3). It implies that the random forces change the steady state of the high wave number mode.

Consequently, the adiabatically approximated Langevin equation for the low wave number mode is given by Eq. (2) where $`x_h`$ is replaced with $`x_h^{}`$ determined by Eq. (12). The reduced Langevin equation is written as

$$\frac{dx_l}{d(\omega _lt)}+\frac{1}{2}\left(\frac{\sqrt{x_h^{}}}{r}+\sqrt{\frac{x_h^{}}{r^2}+4x_l}2\right)x_l+\frac{1}{2}\left(\frac{\sqrt{x_h^{}}}{r}+\sqrt{\frac{x_h^{}}{r^2}+4x_l}\right)x_lw_l=0.$$
(13)
The corresponding Fokker-Planck equation, which determines the time-development of the PDF of $`x_l`$, $`P(x_l,t)`$, is given by

$`{\displaystyle \frac{P(x_l,t)}{(\omega _lt)}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{x_l}}\left({\displaystyle \frac{\sqrt{x_h^{}}}{r}}+\sqrt{{\displaystyle \frac{x_h^{}}{r^2}}+4x_l}2\right)x_lP(x_l,t)`$ (14)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{x_l^2}}\left[{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\sqrt{x_h^{}}}{r}}+\sqrt{{\displaystyle \frac{x_h^{}}{r^2}}+4x_l}\right)^2x_l^2\right]P(x_l,t).`$
### IV.2 The probability density functions and the effect of the random forces

In this subsection, we investigate properties of the steady state PDF of the low wave number mode, $`P(x_l)`$, by numerically solving Eq. (14). Although we will show figures when the parameter $`p`$ is fixed at $`p=30`$, qualitative behavior is the same as that for other $`p`$ values larger than $`1.9\mathrm{}`$.

At first, we show in Fig. 4 the steady state PDF in the small $`r`$ region.

In this region, we have seen that the low wave number mode is suppressed when the turbulent fluctuations are ignored. Figure 4 shows that the probability that the low wave number mode is quiet is large and hence the analysis ignoring random forces is a good approximation in this region.

However, as the value of the parameter $`r`$ is increased, the characteristics change completely. The steady state PDF when $`r>0.49\mathrm{}`$ is shown in Fig. 5.

In the region $`0.49\mathrm{}<r<1`$, we have seen that the subcritical transition occurs and the system is bi-stable if the random forces are ignored. Although there are two peaks in the PDF which are compatible with the previous result without random forces, the valley between the peaks is too shallow to identify the two states. In other words, even if one observes the time-series of $`x_l(t)`$, the value of $`x_l(t)`$ strolls around the two peaks without a sudden change. It means that the bi-stability of the system, i.e., the subcritical transition, is smeared out by the turbulent fluctuations.

The steady state PDF for $`r>1`$ is also shown in Fig. 5. It is seen that the PDF has a single peak. The single peak is compatible with the result of the deterministic analysis, where only one state exists. However, the peak is wide and hence $`x_l`$ fluctuates widely around the location of the peak. The variance of the fluctuation is as large as the peak value and it means that the fluctuation cannot be ignored. Furthermore, the peak value $`0.29\mathrm{}`$ obtained when $`r=1.5,p=30`$ is different from the root of Eqs. (7, 8), which is given as $`x_l0.98`$.

In order to see these consequences explicitly, in Fig. 6, we show the comparison of the contour plot of the PDF and the bifurcation diagram obtained by the analysis ignoring the random forces.

It is seen that the location of the region which has large probability and mainly denoted with the thin solid lines in Fig. 6 is shifted gradually as the value of the parameter $`r`$ is changed. Furthermore, the location of the peak of the PDF is far from the result when the random forces are ignored when $`r0.49\mathrm{}`$.

Next, in order to investigate the meaning of the average in such a turbulent system, we analyze the tail of the steady state PDF $`P(x_l)`$ shown in Fig. 7.

It is shown that the tail of the PDF is well-approximated by the power-law and hence the probability is distributed broadly over large values of the fluctuation. The exponent is about $`0.34`$ for the case shown in Fig. 7. It means that the average, i.e., “the center of mass” of the PDF, is shifted in the large $`x_l`$ direction with the long tail and hence the average is not equal to the value that $`x_l`$ takes with large probability. The values that $`x_l`$ takes with large probability are characterized by the peaks of the PDF and are called the most probable values. The difference between the average and the most probable value is shown in Fig. 8.

Furthermore, hysteresis cannot be captured by observing the average, since the quantity is single-valued from its definition. On the other hand, the most probable values depict the hysteresis as shown in Fig. 9. Although relatively steep variation of the average is observed, the hysteresis can be seen only for the most probable values.

## V Summary and discussion

Finally, we summarize our results and consider their implications. We applied the statistical theory to the model of the inhomogeneous turbulent plasma where two turbulences well-separated in those space-time scales coexist; the “high wave number mode” turbulence (the CDIM micro turbulence) and the “low wave number mode” turbulence (the ITG mode semi-micro turbulence). We derived the development equations of the PDFs of the spectrum amplitudes of the electro-static potential for the characteristic wave numbers. By numerically solving the adiabatically approximated Fokker-Planck equation, the steady state PDFs for the low wave number mode turbulence were obtained. Although the subcritical bifurcation is observed when the turbulent fluctuations are ignored, the shape of the PDF shows that the transition is smeared out by the fluctuations. It means that the approximation ignoring the turbulent fluctuation like the traditional transition theories could overestimate the range of cusp catastrophe.

We also compared the average values with the most probable values and showed that these two characteristic values of stochastic nature are different due to the long power-law tail of the PDF. It means that the average does not mean the value which is expected to realize most probably and hence description of the state of the turbulent systems needs not only the average but also other statistical quantities like the most probable values and the variances.

Consequently, these our results warn that the deterministic description of high-temperature plasmas cannot capture important information of the turbulent systems and the statistical analyses with the Langevin equations and the PDFs are inevitably needed.

###### Acknowledgements.

Nice discussions and critical reading of the manuscript by Prof. Yagi are acknowledged. This work was supported by the Grant-in-Aid for Scientific Research of Ministry of Education, Culture, Sports, Science and Technology, the collaboration programmes of RIAM of Kyushu University and the collaboration programmes of NIFS.


text:
# Physical Review Letters Version of Dynamics and Wrinkling of Radially Propagating Fronts Inferred from Scaling Laws in Channel Geometries

## Abstract

Flame Propagation is used as a prototypical example of expanding fronts that wrinkle without limit in radial geometries but reach a simple shape in channel geometry. We show that the relevant scaling laws that govern the radial growth can be inferred once the simpler channel geometry is understood in detail. In radial geometries (in contrast to channel geometries) the effect of external noise is crucial in accelerating and wrinkling the fronts. Nevertheless, once the interrelations between system size, velocity of propagation and noise level are understood in channel geometry, the scaling laws for radial growth follow.

The main idea of this Letter is that in order to derive scaling laws for unstable front propagation in radial geometry, it is useful to study noisy propagation in channel geometries, in which the noiseless dynamics results usually in simple shapes of the advancing fronts . Famous examples of such situations are Laplacian growth patterns and unstable flame fronts. Laplacian growth in a channel geometry results in a single finger whose width is uniquely determined by the channel width $`L`$, the velocity of the tip $`v`$ and the surface tension. Laplacian growth in radial geometries results in convoluted and complex structures whose full characterization still eludes repeated theoretical attacks. Flame propagation offers a similar situation; in channel geometries the flame front looks like a giant cusp whose velocity $`v`$ is constant. In radial geometries the flame front accelerates with time, and while propagating it wrinkles by the addition of hierarchies of cusps of all sizes.

The understanding of radial geometries requires control of the effect of noise on the unstable dynamics of propagation. It is particularly difficult to achieve such a control in radial geometries due to the vagueness of the distinction between external noise and noisy initial conditions. Channel geometries are simpler when they exhibit a stable solution for growth in the noiseless limit. One can then study the effects of external noise in such geometries without any ambiguity. If one finds rules to translate the resulting understanding of the effects of noise in channel growth to radial geometries, one can derive the scaling laws in the later situation in a satisfactory manner. We will exemplify the details of such a translation in the context of premixed flames that exist as self sustaining fronts of exothermic chemical reactions in gaseous combustion. But our contention is that similar ideas should be fruitful also in other contexts of unstable front propagation. Needless to say, there are aspects of the front dynamics and statistics in the radial geometry that cannot be explained from observations of fronts in a channel geometry; examples of such aspects are discussed at the end of this Letter.

Mathematically our example is described  by an equation of motion for the angle- dependent modulus of the radius vector of the flame front, $`R(\theta ,t)`$:

$`{\displaystyle \frac{R}{t}}`$ $`=`$ $`{\displaystyle \frac{U_b}{2R_{0}^{}{}_{}{}^{2}(t)}}\left({\displaystyle \frac{R}{\theta }}\right)^2+{\displaystyle \frac{D_M}{R_{0}^{}{}_{}{}^{2}(t)}}{\displaystyle \frac{^2R}{\theta ^2}}`$ (1)
$`+`$ $`{\displaystyle \frac{\gamma U_b}{2R_0(t)}}I(R)+U_b.`$ (2)
Here 0$`<\theta <2\pi `$ is an angle and the constants $`U_b,D_M`$ and $`\gamma `$ are the front velocity for an ideal cylindrical front, the Markstein diffusivity and the thermal expansion coefficient respectively. $`R_0(t)`$ is the mean radius of the propagating flame:

$$R_0(t)=\frac{1}{2\pi }_0^{2\pi }R(\theta ,t)𝑑\theta .$$
(3)
The functional $`I(R)`$ is best represented in terms of its Fourier decomposition. Its Fourier component is $`|k|R_k`$ where $`R_k`$ is the Fourier component of $`R`$. Simulations of this equation, as well as experiments in the parameter regime for which this equation is purportedly relevant, indicate that for large times $`R_0`$ grows as a power in time

$$R_0(t)=(const+t)^\beta ,$$
(4)
with $`\beta >1`$, and that the width of the interface $`W`$ grows with $`R_0`$ as

$$W(t)R_0(t)^\chi ,$$
(5)
with $`\chi <1`$. Note that this exponent may differ significantly from the exponent which characterizes the scaling of correlation functions . This phenomenon is common to cases in which large features dominate the width of the graph. In a previous publication  we showed that $`\beta `$ and $`\chi `$ satisfy the scaling relation $`\beta =1/\chi `$. In this Letter we determine the numerical value of these exponents from the analysis of noisy channel growth.

In channels there is a natural lengthscale, the width $`\stackrel{~}{L}`$ of the channel. The translation of channel results to radial geometry will be based on the identification in the latter context of the time dependent scale $`(t)`$ that plays the role of $`\stackrel{~}{L}`$ in the former. To do this we need first to review briefly the main pertinent results for noisy channel growth. In channel geometry the equation of motion is written in terms of the position $`h(x,t)`$ of the flame front above the $`x`$-axis. After appropriate rescalings  it reads:

$$\frac{h(x,t)}{t}=\frac{1}{2}\left[\frac{h(x,t)}{x}\right]^2+\nu \frac{^2h(x,t)}{x^2}+I\{h(x,t)\}+1.$$
(6)
It is convenient to rescale the domain size further to $`0<\theta <2\pi `$, and to change variables to $`u(\theta ,t)h(\theta ,t)/\theta `$. In terms of this function we find

$$\frac{u(\theta ,t)}{t}=\frac{u(\theta ,t)}{L^2}\frac{u(\theta ,t)}{\theta }+\frac{\nu }{L^2}\frac{^2u(\theta ,t)}{\theta ^2}+\frac{1}{L}I\{u(\theta ,t)\}$$
(7)
where $`L=\stackrel{~}{L}/2\pi `$. In noiseless conditions this equation admits exact solutions that are represented in terms of $`N`$ poles whose position $`z_j(t)x_j(t)+iy_j(t)`$ in the complex plane is time dependent:

$$u(\theta ,t)=\nu \underset{j=1}{\overset{N}{}}\mathrm{cot}\left[\frac{\theta z_j(t)}{2}\right]+c.c.,$$
(8)
The steady state for channel propagation is unique and linearly stable; it consists of $`N(L)`$ poles which are aligned on one line parallel to the imaginary axis. The geometric appearance of the flame front is a giant cusp, analogous to the single finger in the case of Laplacian growth in a channel. The height of the cusp is proportional to $`L`$, and the propagation velocity is a constant of the motion. The number of poles in the giant cusp is linear in $`L`$,

$$N(L)=\left[\frac{1}{2}\left(\frac{L}{2\pi }+1\right)\right].$$
(9)
The introduction of additive random noise to the dynamics changes the picture qualitatively. It is convenient to add noise to the equation of motion in Fourier representation by adding a white noise $`\eta _k`$ for every $`k`$ mode. The noise correlation function satisfies the relation $`<\eta _k(t)\eta _k^{}(t^{})>=\delta (k+k^{})\delta (tt^{})2f/L`$. The noise in our simulations is taken from a flat distribution in the interval $`[\sqrt{2f/L},\sqrt{2f/L}]`$; this guarantees that when the system size changes, the typical noise per unit length of the flame front remains constant. It was shown  that for moderate but fixed noise levels the average velocity $`v`$ of the front increases with $`L`$ as a power law. In our present simulations we found

$$vL^\mu ,\mu 0.35\pm 0.03.$$
(10)
For a fixed system size $`L`$ the velocity has also a power law dependence on the level of the noise, but with a much smaller exponent: $`vf^\xi ,\xi 0.02`$. These results were understood theoretically by analysing the noisy creation of new poles that interact with the poles defining the giant cusp. It was shown that the system is always linearly stable or marginally stable for the introduction of new poles, but with increasing $`L`$ all the eigenvalues (which are negative for a positive measure of $`L`$ values) all decrease in absolute magnitude like $`1/L^2`$, making the system more and more susceptible to noisy perturbations . The picture used remains valid as long as the poles that are introduced by the noisy perturbation do not destroy the identity of the giant cusp. Indeed, the numerical simulations show that in the presence of moderate noise the additional poles appear as smaller cusps that are constantly running towards the giant cusp. Our point here, is not to predict the numerical values of the scaling exponents in the channel (this was done in ref., but to use them to predict the scaling exponents characterizing the acceleration and the geometry of the flame front in radial geometry.

Superficially it seems that in radial geometry the growth pattern is qualitatively different. In fact, close observation of the growth patterns (see Fig.1) shows that most of the time there exist some big cusps that attract other smaller cusps, but that every now and then “new” big cusps form and begin to act as local absorbers of small cusps that appear randomly. The understanding of this phenomenon gives the clue how to translate results from channels to radial growth.

In the radial geometry we non-dimensionalize Equation (2) by rescaling $`Rr,R_0r_0`$ and $`t\tau `$, and finally rewriting the equation in terms of $`ur/\theta `$:

$$\frac{u}{\tau }=\frac{u}{r_0^2}\frac{u}{\theta }+\frac{1}{r_0^2}\frac{^2u}{\theta ^2}+\frac{\gamma }{2r_0}I\{u\}.$$
(11)
To complete this equation we need a second one for $`r_0(t)`$. This equation is obtained by averaging (2) over the angles (i.e. operating on the equation with $`\frac{1}{2\pi }_0^{2\pi }𝑑\theta `$):

$$\frac{dr_0}{d\tau }=\frac{1}{2r_0^2}\frac{1}{2\pi }_0^{2\pi }u^2𝑑\theta +1.$$
(12)
These equations again admit exact solutions in terms of poles, of the form of Eq.(7). It is easy write down the equations of motion of the poles and check that the poles are attractive along the real direction (which means physically that they are attracted along the angular coordinate) but they are repulsive along the imaginary direction, which is associated with the radial coordinate. If it were not for the stretching that is caused by the increase of the radius (and with it the perimeter), all the poles would have coalesced into one giant cusp. Thus we have a competition between pole attraction and stretching. Since the attraction decreases with the distance between the poles in the angular directions, there is always an initial critical length scale above which poles cannot coalesce their real coordinates when time progresses.

Suppose now that noise adds new poles to the system. The poles do not necessarily merge their real positions with existing cusps. If we have a large cusp made from the merging of the real coordinates $`x_c`$ of $`N_c`$ poles, we want to know whether a nearby pole with real coordinate $`x_1`$ will merge with this large cusp. The answer will depend of course on the distance $`Dr_0|x_cx_1|`$. A direct calculation , using the equation of motions for the poles shows that there exists a critical length $`(r_0)`$ such that if $`D>(r_0)`$ the single pole never merges with the giant cusp. The result of the calculation is that

$$r_0^{1/\beta }.$$
(13)
Note that a failure of a single pole to be attracted to a large cusp means that tip-splitting has occurred. This is the exact analog of tip-splitting in Laplacian growth.

It is now time to relate the channel and radial geometries. We identify the typical scale in the radial geometry as $`Wr_0^\chi `$. On the one hand this leads to the scaling relation $`\beta =1/\chi `$. On the other hand we use the result established in a channel, (10), with this identification of a scale, and find $`\dot{r}_0=r_0^{\chi \mu }`$. Comparing with (4) we find:

$$\beta =\frac{1}{(1\chi \mu )}.$$
(14)
This result leads us to expect two dynamical regimes for our problem. Starting from smooth initial conditions, in relatively short times the roughness exponent remains close to unity. This is mainly since the typical scale $``$ is not relevant yet, and most of the poles that are generated by noise merge into a few larger cusps. In later times the roughening exponent settles at its asymptotic value, and all the asymptotic scaling relations used above become valid. We thus expect $`\beta `$ to decrease from $`1/(1\mu )`$ to an asymptotic value determined by $`\chi =1/\beta `$ in (14):

$$\beta =1+\mu 1.35\pm 0.03.$$
(15)
The expected value of $`\chi `$ is thus $`\chi =0.74\pm 0.03`$.

We tested these predictions in numerical simulations. We integrated Eq.(7), and in Fig. 2 we display the results for the growth velocity as a function of time. After a limited domain of exponential growth we observe a continuous reduction of the time dependent exponent. In the initial region we get $`\beta =1.65\pm 0.1`$ while in the final decade of the temporal range we find $`\beta =1.35\pm 0.1`$. We consider this a good agreement with (15). A second important test is provided by measuring the width of the system as a function of the radius, see Fig.2b. Again we observe a cross-over related to the initial dynamics; In the last temporal decade the exponent settles at $`\chi =0.75\pm .1`$. We conclude that at times large enough to observe the asymptotics our predictions are verified.

Finally, we stress some differences between radial and channel geometries. Fronts in a channel exhibit mainly one giant cusp which is only marginally disturbed by the small cusps that are introduced by noise. In the radial geometry, as can be concluded from the discussion above, there exist at any time cusps of all sizes from the smallest to the largest. This broad distribution of cusps (and scales) must influence correlation function in ways that differ qualitatively from correlation functions computed in channel geometries.

To make the point clear we exhibit in Fig.3a the structure function

$$F(y)\sqrt{|R(x+y)R(x)|^2}$$
(16)
computed for a typical radial front, with $`x=R\theta `$. To stress the scaling region we exhibit the second derivative of this function in Fig.3b. The low end of the graph can be fitted well by a power law $`y^\alpha `$ with $`\alpha 0.6`$. This indicates that $`F(y)Ay+By^{2\alpha }`$. In a channel geometry we get entirely different structure functions that do not exhibit such scaling functions at all. The way to understand this behavior in the radial geometry is to consider a distribution of cusps that remain distinct from each other but whose scales are distributed according to some distribution $`P(\mathrm{})`$. For each of these cusps there is a contribution to the correlation function of the form $`f(y,\mathrm{})\mathrm{}g(y/\mathrm{})`$ where $`g(x)`$ is a scaling function, $`g(x)x`$ for $`x<1`$ and $`g(x)`$ constant for $`x>1`$. The total correlation can be estimated (when the poles are distinct) as

$$F(y)\underset{\mathrm{}}{}P(\mathrm{})\mathrm{}g(y/\mathrm{}).$$
(17)
The first derivative leaves us with $`_{\mathrm{}}P(\mathrm{})g^{}(y/\mathrm{})`$, and using the fact that $`g^{}`$ vanishes for $`x>1`$ we estimate $`F^{}(y)=_{\mathrm{}=y}^WP(\mathrm{})`$. The second derivative yields $`F^{\prime \prime }(y)P(y)`$. Thus the structure function is determined by the scale distribution of cusps, and if the latter is a power law, this should be seen in the second derivative of $`F(y)`$ as demonstrated in Fig.3. The conclusion of this analysis is that the radial case exhibits a scaling function that characterizes the distributions of cusps, $`P(\mathrm{})\mathrm{}^\alpha `$.

In summary, we demonstrated that it is possible to use information about noisy channel dynamics to predict nontrivial features of the radial evolution, such as the acceleration and roughening exponents. It would be worthwhile to examine similar ideas in the context of Laplacian growth patterns.

###### Acknowledgements.

This work has been supported in part by the Basic Research Fund administered by the Israel Academy for Science and Humanities. IP thanks Mitchell Feigenbaum for the hospitality at the Rockefeller University.


text:
# References

A System of Differential Equations for the Airy Process

Craig A. Tracy
Department of Mathematics
University of California, Davis, CA 95616
email address: tracy@math.ucdavis.edu

Harold Widom
Department of Mathematics
University of California, Santa Cruz, CA 95064
email address: widom@math.ucsc.edu

Abstract

The Airy process $`\tau A_\tau `$ is characterized by its finite-dimensional distribution functions

$$\text{Pr}\left(A_{\tau _1}<\xi _1,\mathrm{},A_{\tau _m}<\xi _m\right).$$

For $`m=1`$ it is known that $`\text{Pr}\left(A_\tau <\xi \right)`$ is expressible in terms of a solution to Painlevé II. We show that each finite-dimensional distribution function is expressible in terms of a solution to a system of differential equations.

I. Introduction

The Airy process $`\tau A_\tau `$, introduced by Prähofer and Spohn , is the limiting stationary process for a certain $`1+1`$-dimensional local random growth model called the polynuclear growth model (PNG). It is conjectured that the Airy process is, in fact, the limiting process for a wide class of random growth models. (This class is called the $`1+1`$-dimensional KPZ universality class in the physics literature .) The PNG model is closely related to the length of the longest increasing subsequence in a random permutation . This fact together with the result of Baik, Deift and Johansson  on the limiting distribution of the length of the longest increasing subsequence in a random permutation shows that the distribution function $`\text{Pr}\left(A_\tau <\xi \right)`$ equals the limiting distribution function, $`F_2(\xi )`$, of the largest eigenvalue in the Gaussian Unitary Ensemble . $`F_2`$ is expressible either as a Fredholm determinant of a certain trace-class operator (the Airy kernel) or in terms of a solution to a nonlinear differential equation (Painlevé II). The finite-dimensional distribution functions

$$\text{Pr}\left(A_{\tau _1}<\xi _1,\mathrm{},A_{\tau _m}<\xi _m\right)$$

are expressible as a Fredholm determinant of a trace-class operator (the extended Airy kernel) . It is natural to conjecture  that these distribution functions are also expressible in terms of a solution to a system of differential equations. It is this last conjecture which we prove. For $`m=2`$ this conjecture was also proved, in a different form, by Adler and van Moerbeke .

II. Statement

The Airy process is characterized by the probabilities

$$\mathrm{Pr}\left(A_{\tau _1}<\xi _1,\mathrm{},A_{\tau _m}<\xi _m\right)=det(IK),$$

where $`K`$ is the operator with $`m\times m`$ matrix kernel having entries

$$K_{ij}(x,y)=L_{ij}(x,y)\chi _{(\xi _j,\mathrm{})}(y)$$

and

$$L_{ij}(x,y)=\{\begin{array}{cc}_0^{\mathrm{}}e^{z(\tau _i\tau _j)}\mathrm{Ai}(x+z)\mathrm{Ai}(y+z)𝑑z\hfill & \mathrm{if}ij,\hfill \\ & \\ _{\mathrm{}}^0e^{z(\tau _i\tau _j)}\mathrm{Ai}(x+z)\mathrm{Ai}(y+z)𝑑z\hfill & \mathrm{if}i<j.\hfill \end{array}$$

We assume throughout that $`\tau _1<\mathrm{}<\tau _m`$, and think of $`K`$ as acting on the $`m`$-fold direct sum of $`L^2(\alpha ,\mathrm{})`$ where $`\alpha <\mathrm{min}\xi _j`$.

To state the result we let $`R=K(IK)^1`$ and let $`A(x)`$ denote the $`m\times m`$ diagonal matrix $`\mathrm{diag}(\mathrm{Ai}(x))`$ and $`\chi (x)`$ the diagonal matrix $`\mathrm{diag}(\chi _j(x))`$, where $`\chi _j=\chi _{(\xi _j,\mathrm{})}`$. Then we define the matrix functions $`Q(x)`$ and $`\stackrel{~}{Q}(x)`$ by

$$Q=(IK)^1A,\stackrel{~}{Q}=A\chi (IK)^1$$

(where for $`\stackrel{~}{Q}`$ the operators act on the right). These and $`R(x,y)`$ are functions of the $`\xi _j`$ as well as $`x`$ and $`y`$. We define the matrix functions $`q,\stackrel{~}{q}`$ and $`r`$ of the $`\xi _j`$ only by

$$q_{ij}=Q_{ij}(\xi _i),\stackrel{~}{q}_{ij}=\stackrel{~}{Q}_{ij}(\xi _j),r_{ij}=R_{ij}(\xi _i,\xi _j).$$

Finally we let $`\tau `$ denote the diagonal matrix $`\mathrm{diag}(\tau _j)`$.

Our differential operator is $`𝒟=_j_j`$, where $`_j=/\xi _j`$, and the system of equations is

$`𝒟^2q`$ $`=`$ $`\xi q+2q\mathrm{\Theta }\stackrel{~}{q}q2[\tau ,r]q,`$ (1)
$`𝒟^2\stackrel{~}{q}`$ $`=`$ $`\stackrel{~}{q}\xi +2\stackrel{~}{q}q\mathrm{\Theta }\stackrel{~}{q}2\stackrel{~}{q}[\tau ,r],`$ (2)
$`𝒟r`$ $`=`$ $`q\mathrm{\Theta }\stackrel{~}{q}+[\tau ,r].`$ (3)
Here the brackets denote commutator, $`\xi `$ denotes the diagonal matrix $`\mathrm{diag}(\xi _j)`$ and $`\mathrm{\Theta }`$ the matrix with all entries equal to one.

This can be interpreted as a system of ordinary differential equations if we replace the variables $`\xi _1,\mathrm{},\xi _m`$ by $`\xi _1+\xi ,\mathrm{},\xi _m+\xi `$, where $`\xi _1,\mathrm{},\xi _m`$ are fixed and $`\xi `$ variable. Then $`𝒟=d/d\xi `$, and the $`\xi _j`$ are regarded as parameters.

To get a representation for $`det(IK)`$ observe that

$$_jK=L\delta _j,$$
(4)
where the last factor denotes multiplication by the diagonal matrix with all entries zero except for the $`j^{\mathrm{th}}`$, which equals $`\delta (x\xi _j)`$. We deduce that

$$_j\mathrm{log}det(IK)=\mathrm{Tr}(IK)^1_jK=R_{jj}(\xi _j,\xi _j).$$

Hence $`𝒟\mathrm{log}det(IK)=\mathrm{Tr}r`$, and so it follows from (3) that

$$𝒟^2\mathrm{log}det(IK)=\mathrm{Tr}q\mathrm{\Theta }\stackrel{~}{q}$$

since the trace of $`[\tau ,r]`$ equals zero. This gives the representation

$$det(IK)=\mathrm{exp}\left\{_0^{\mathrm{}}\eta \mathrm{Tr}q(\xi +\eta )\mathrm{\Theta }\stackrel{~}{q}(\xi +\eta )𝑑\eta \right\}.$$

Here the determinant is evaluated at $`(\xi _1,\mathrm{},\xi _m)`$ and in the integral $`\xi +\eta `$ is shorthand for $`(\xi _1+\eta ,\mathrm{},\xi _m+\eta ).`$

If $`m=1`$ the commutators drop out, $`q=\stackrel{~}{q}`$, equations (1) and (2) are Painlevé II and these are the previously known results.

III. Proof

The proof will follow along the lines of the derivation in  for the case $`m=1`$. There the kernel was “integrable” in the sense that its commutator with $`M`$, the operator of multiplication by $`x`$, was of finite rank. The same was then true of the resolvent kernel, which was useful. But now our kernel is not integrable, so there will necessarily be some differences.

With $`D=d/dx`$ we compute that

$$[D,K]_{ij}=\mathrm{Ai}(x)\mathrm{Ai}(y)\chi _j(y)+L_{ij}(x,\xi _j)\delta (y\xi _j)+(\tau _i\tau _j)K_{ij}(x,y).$$

Equivalently,

$$[D,K]=A(x)\mathrm{\Theta }A(y)\chi (y)+L\delta +[\tau ,K],$$

where $`\delta =_j\delta _j`$, multiplication by the matrix $`\mathrm{diag}(\delta (x\xi _j))`$, and $`L`$ is the operator with kernel $`L_{ij}(x,y)`$. (For clarity we sometimes write the kernel of an operator in place of the operator itself.) To obtain $`[D,R]`$ we replace $`K`$ by $`KI`$ in the commutators and left- and right-multiply by $`\rho =(IK)^1`$. The result is

$$[D,R]=Q(x)\mathrm{\Theta }\stackrel{~}{Q}(y)+R\delta \rho +[\tau ,\rho ].$$
(5)
We have already defined the matrix functions $`Q`$ and $`\stackrel{~}{Q}`$ and we define

$$P=(IK)^1A^{},u=(\stackrel{~}{Q},\mathrm{Ai})=\stackrel{~}{Q}(x)\mathrm{Ai}(x)𝑑x.$$

It follows from (5) and the fact that $`\tau `$ and $`A`$ commute that

$$Q^{}=PQ\mathrm{\Theta }u+R\delta Q+[\tau ,Q].$$
(6)
Next, it follows from (4) that

$$_jR=R\delta _j\rho ,$$
(7)
and it follows from this that $`_jQ=R\delta _jQ.`$ Summing over $`j`$, adding to (6) and evaluating at $`\xi _k`$ give

$$𝒟Q(\xi _k)=P(\xi _k)Q(\xi _k)\mathrm{\Theta }u+[\tau ,Q(\xi _k)].$$

If we define $`p_{ij}=P_{ij}(\xi _i)`$ then we obtain

$$𝒟q=pq\mathrm{\Theta }u+[\tau ,q].$$
(8)
Next we use the facts that $`D^2M`$ commutes with $`L`$ and that $`M`$ commutes with $`\chi `$. It follows that

$$[D^2M,K]=[D^2M,L\chi ]=L[D^2M,\chi ]=L[D^2,\chi ]=L(\delta D+D\delta ).$$

It follows from this that

$$[D^2M,\rho ]=\rho L\delta D\rho +\rho LD\delta \rho .$$

Applying both sides to $`A`$ and using the fact that $`(D^2M)A=0`$ we obtain

$$Q^{\prime \prime }(x)xQ(x)=\rho L\delta Q^{}+\rho LD\delta Q.$$
(9)
The first term on the right equals $`R\delta Q^{}`$. For the second term observe that

$$\rho LD\chi =\rho L\chi D+\rho L[D,\chi ]=RD+\rho L\delta ,$$

so we can interpret that term as $`R_y\delta Q`$ (the subscript denotes partial derivative) where $`R_y(x,y)`$ is interpreted as not containing the delta-function summand which arises from the jumps of $`R`$. With this interpretation of $`R_y`$ we can write the second term on the right as $`R_y\delta Q`$. Thus,

$$Q^{\prime \prime }(x)xQ(x)=R\delta Q^{}R_y\delta Q.$$

Using this we obtain from (6)

$$P^{}=xQ(x)+R\delta Q^{}R_y\delta Q+Q^{}\mathrm{\Theta }uR_x\delta Q[\tau ,Q^{}],$$

and then from (6) once more

$$P^{}=xQ(x)+R\delta (PQ\mathrm{\Theta }u+R\delta Q+[\tau ,Q])R_y\delta Q$$

$$+(PQ\mathrm{\Theta }u+R\delta Q+[\tau ,Q])\mathrm{\Theta }uR_x\delta Q[\tau ,PQ\mathrm{\Theta }u+R\delta Q+[\tau ,Q]].$$

It follows from (5) that

$$R_x+R_y=Q(x)\mathrm{\Theta }\stackrel{~}{Q}(y)+R\delta R+[\tau ,\rho ].$$

(We replaced $`R\delta \rho `$ by $`R\delta R`$ since, recall, $`R_y`$ does not contain delta-function summands.) We use this and also the identity $`R\delta [\tau ,Q][\tau ,R\delta Q]=[\tau ,R\delta ]Q`$, and the fact that $`\delta `$ and $`\tau `$ commute. The result is that

$$P^{}=xQ(x)+R\delta P+Q(x)\mathrm{\Theta }\stackrel{~}{Q}\delta Q+(PQ\mathrm{\Theta }u+[\tau ,Q])\mathrm{\Theta }u$$

$$2[\tau ,R]\delta Q[\tau ,PQ\mathrm{\Theta }u+[\tau ,Q]].$$

It follows from (7) that $`_jP=R\delta _jP`$. Summing over $`j`$, adding to the above and evaluating at $`\xi _k`$ give

$$𝒟P(\xi _k)=\xi _kQ(\xi _k)+Q(\xi _k)\mathrm{\Theta }\stackrel{~}{Q}\delta Q+(P(\xi _k)Q(\xi _k)\mathrm{\Theta }u+[\tau ,Q(\xi _k)])\mathrm{\Theta }u$$

$$2[\tau ,R(\xi _k,)]\delta Q[\tau ,P(\xi _k)Q(\xi _k)\mathrm{\Theta }u+[\tau ,Q(\xi _k)]].$$

Hence $`𝒟p`$ is equal to

$$\xi q+q\mathrm{\Theta }\stackrel{~}{q}q+(pq\mathrm{\Theta }u+[\tau ,q])\mathrm{\Theta }u2[\tau ,r]q[\tau ,pq\mathrm{\Theta }u+[\tau ,q]].$$

Equivalently, in view of (8),

$$𝒟p=\xi q+q\mathrm{\Theta }\stackrel{~}{q}q+𝒟q\mathrm{\Theta }u2[\tau ,r]q[\tau ,𝒟q].$$
(10)
Let us compute $`𝒟u`$. We have

$$u_{ij}=\mathrm{Ai}(x)\chi _i(x)\rho _{ij}(x,y)\mathrm{Ai}(y)𝑑x𝑑y,$$

and so

$$_ku_{ij}=\delta _{ik}\mathrm{Ai}(\xi _k)\rho _{kj}(\xi _k,y)\mathrm{Ai}(y)𝑑y$$

$$\mathrm{Ai}(x)\chi _i(x)[R_{ik}(x,\xi _k)\rho _{kj}(\xi _k,y)]\mathrm{Ai}(y)𝑑x𝑑y,$$

where we use (7) again. This is equal to

$$\delta _{ik}\mathrm{Ai}(\xi _k)Q_{kj}(\xi _k)\left(\stackrel{~}{Q}_{ik}(\xi _k)\delta _{ik}\mathrm{Ai}(\xi _k)\right)Q_{kj}(\xi _k),$$

and so

$$_ku_{ij}=\stackrel{~}{Q}_{ik}(\xi _k)Q_{kj}(\xi _k).$$
(11)
This gives

$$𝒟u=\stackrel{~}{q}q.$$
(12)
Next, we find from (7) and (5) that

$$𝒟R(\xi _j,\xi _k)=Q(\xi _j)\mathrm{\Theta }\stackrel{~}{Q}(\xi _k)+[\tau ,R(\xi _j,\xi _k)].$$

This gives $`𝒟r=q\mathrm{\Theta }\stackrel{~}{q}+[\tau ,r]`$, which is equation (3).

To get equation (1) we apply $`𝒟`$ to (8) and use (10) and (12). We find that

$$𝒟^2q=\xi q+q\mathrm{\Theta }\stackrel{~}{q}q+𝒟q\mathrm{\Theta }u2[\tau ,r]q[\tau ,𝒟q]𝒟q\mathrm{\Theta }u+q\mathrm{\Theta }\stackrel{~}{q}q+[\tau ,𝒟q]$$

$$=\xi q+2q\mathrm{\Theta }\stackrel{~}{q}q2[\tau ,r]q,$$

which is (1).

Finally, to get equation (2) we use the fact that $`\chi _j(y)\rho _{jk}(y,x)`$ is equal to $`\chi _k(x)`$ times $`\rho _{kj}^{}(x,y)`$, where $`\rho ^{}`$ is the resolvent kernel for the matrix kernel with $`i,j`$ entry $`L_{ji}(x,y)\chi _j(y)`$. Hence $`\stackrel{~}{Q}_{jk}(x)`$ is equal to $`\chi _k(x)`$ times the $`Q_{kj}(x)`$ associated with $`L_{ji}`$. Consequently for all the differentiation formulas we have for the $`Q_{kj}(\xi _k)`$, etc., there are analogous formulas for the $`\stackrel{~}{Q}_{jk}(\xi _k)`$, etc.. The difference is that we have to reverse subscripts and replace $`r`$ by $`r^t`$ and $`\tau `$ by $`\tau `$. The upshot is that, by computations analogous to those used to derive (1), we derive another equation which can be obtained from (1) by making the replacements $`q\stackrel{~}{q}^t`$, $`\stackrel{~}{q}q^t`$, $`rr^t`$, $`\tau \tau `$ and then taking transposes. The result is equation (2).

Acknowledgment

This work was supported by the National Science Foundation through grants DMS-9802122 and DMS-9732687.


text:
# Energy Level Statistics of a Critical Random Matrix Ensemble

## 1 Introduction

Random matrix ensembles (RME) provide a natural frame work for the statistical description of quantized energy levels of complex systems. This is probably best explained by paraphrasing Dyson : in ordinary statistical mechanics a renunciation of knowledge about the system is made i.e., by assuming that all states of a large ensemble are equally probable. Therefore only the overall behavior of the system is accessible. Thus orthodox statistical mechanics is clearly inadequate for the discussion of energy level statistics of complex systems. In energy level statistics one wishes to make statements about the fine detail of the energy level structure, and such statements cannot be made in terms of an ensemble of states. Random matrix theory (RMT) turns out to be a suitable tool in making such statements about the eigen-spectrum of complex systems e.g., energy levels of highly excited nuclei . The underlying assumption in RMT is that one renounces exact knowledge not of the state but of the nature of the system itself. This is the fundamental difference between ordinary statistical mechanics and energy level statistics.

There are three RME that emerge from an analysis of real eigen-spectra of random Hermitian matrices ; the Gaussian orthogonal ensemble (GOE) for a system with time-reversal without spin-$`\frac{1}{2}`$ interactions, the Gaussian symplectic ensemble (GSE) for a system with time-reversal and spin-$`\frac{1}{2}`$ interaction, and the Gaussian unitary ensemble (GUE) for a system without time-reversal. These ensembles (generally called Wigner-Dyson (WD) statistics) exhibit strong energy-level correlations leading to the phenomenon of level repulsion  in their respective eigen-spectra. One could add the Poisson statistics (PS) as the fourth ensemble which does not exhibit level repulsion in it’s corresponding eigen-spectra.

In general, the WD statistics describe a system of delocalized states and the PS describe a system of localized states. The transition between these two cases controlled by the strength of disorder (i.e the Anderson model ) is a quantum critical phenomenon. At a critical disorder the states acquire a special property of multifractality  which marks a qualitative difference as compared to truly delocalized states in a metal and truly localized states in an insulator. These critical states correspond to the critical level statistics (CLS) .

There are numerous suggestions in literature  describing CLS by ensembles of random matrices. Initially, there were two main requirements for the critical RME: the first was that the level statistics be system-size independent and secondly, the ensemble should be able to interpolate between WD and PS . Later on it was realized  that the key property of the critical states relevant for the CLS is multifractality. The power-law banded RME suggested in  emerged as the best candidate for the critical RME as its level statistics is known to possess these two properties and its eigenfunction statistics is multifractal. Thus it has been suggested  to be the critical random matrix ensemble (CrRME). In the limit of large bandwidth this model turns out to be the best deformation of the WD-RME that generates correlations in the eigenstates but retains all the basic properties of WD theory, including the factorizability of higher spectral correlation functions and mapping onto the Calogero-Sutherland model of fictitious interacting fermions .

In the absence of an exact solution of CrRME various mappings on and approximate correspondences with exactly solvable models have been used in order to study the properties and CrRME. Here we mention the mapping onto the nonlinear sigma-model , and the approximate equivalence  between the CrRME and the exactly solvable RME suggested in Refs.. Making use of the exact solution  it was possible  to compute the level spacing distribution function (LSDF) for this model and compare it with the exact diagonalization of the 3d Anderson model for the critical value of disorder. The coincidence of results (for certain values of the parameter $`b`$ which were found from the fitting of the far tails of LSDF) was amazing : in all three Dyson symmetry classes the deviation was not larger than between the Wigner surmise and the exact WD distribution function.

In this paper we study the level number variance in the CrRMT suggested . The ultimate goal here is to compare the numerical simulations of this critical ensemble and the analytical results for the exactly solvable ensemble in Ref.  (which will be referred to as MNS) for all values of the control parameter $`b`$. For large $`b1`$ (weak multifractality) these two ensembles are known to be equivalent . The first correction to the Poisson spectral statistics at small $`b1`$ (strong multifractality) is also the same . We would like to check these predictions numerically for the unitary CrRME and to see how large are the deviations between these two RME for $`b1`$.

The paper is organized as follows: in Section 2 we give a general overview of the level-number variance ($`\mathrm{\Sigma }_2`$-statistics) and the calculation of the level compressibility, in Section 3 we give some numerical details and finally a summary and outlook of this work is given in Section 4.

## 2 $`\mathrm{\Sigma }_2`$-Statistics and Level Compressibility

Our model is the power-law random banded matrix defined as an ensemble of random Hermitian $`N\times N`$ matrices belonging to GUE. The real and imaginary parts of the matrix elements are independently distributed Gaussian random variables with zero mean $`H_{ij}=0`$ and the variance of the matrix elements is given by;

$$\left|H_{ij}\right|^2=\left[1+\frac{1}{b^2}\frac{\mathrm{sin}^2(\pi \left|ij\right|/N)}{(\pi /N)^2}\right]^1\times \{\begin{array}{cc}\frac{1}{\beta }\hfill & \text{if }i=j\hfill \\ \frac{1}{2}\hfill & \text{if }ij\hfill \end{array}$$
(1)
where $`b`$ is a parameter that controls the fractal dimensions of the critical eigenfunctions and $`\beta `$ is the symmetry parameter i.e., $`\beta =1,2`$ and $`4`$ for the GOE, GUE and GSE respectively. This model remains critical for any arbitrary value of $`b`$.

From a physical point of view the GUE can be realized by putting a system in a strong external magnetic field i.e., provided the splitting of levels by the magnetic field is of the same order of magnitude as the average level spacing in the absence of the field.

In this paper we study the level compressibility $`\chi (N,b)`$ via the level-number variance $`\mathrm{\Sigma }_2\left(n\right)`$ of the above model. The level-number variance is a statistical quantity that provides a quantitative measure of the long-range ”rigidity” of the energy spectrum ;

$$\mathrm{\Sigma }_2\left(n\right)=n^2n^2$$
(2)
The behavior of $`\mathrm{\Sigma }_2`$ is known (as a function of the mean number of levels $`n1`$) in the delocalized and localized phases respectively. However, in the critical regime, the level-number variance has been conjectured to be Poisson-like ;

$$\mathrm{\Sigma }_2\{\begin{array}{cc}\frac{2}{\pi ^2\beta }\mathrm{ln}\left(n\right)\hfill & \text{delocalized,}\hfill \\ \chi n\hfill & \text{critical,}\hfill \\ n\hfill & \text{localized,}\hfill \end{array}$$
(3)
where the level compressibility $`\chi `$ is another important parameter to characterize the localization-delocalization transition (LDT) and takes values $`0\chi 1`$, being zero in the delocalized state and unity in the localized state. Formally, one computes $`\chi `$ as;

$$\chi \underset{n\mathrm{}}{lim}\underset{N\mathrm{}}{lim}\frac{d\mathrm{\Sigma }_2\left(n\right)}{dn}$$
(4)
The limits in Eqn.(4) do not commute and their non-commutativity is attributed to the violation of the normalization sum rule . It is remarkable (see  for details) that $`\chi `$ can be expressed in terms of the fractal or the correlation dimension $`D_2`$ and the spatial dimension $`d`$. In Ref. an explicit and very simple relationship between the critical spectral and eigenfunction statistics has been suggested:

$$\chi =\frac{1}{2}\left(1\frac{D_2}{d}\right)$$
(5)
There is vast analytical and numerical evidence that this formula is correct at weak multifractality. However, its status for the case of strong multifractality (small $`D_2`$) raises much doubts . We later see that whilst our numerical data are in agreement with Eqn.(5) in the limit $`b^10`$ they are in disagreement  for the opposite limit. At the same time the overall agreement with the analytical result Eqn.(10) for the MNS model is unexpectedly good for all values of $`b`$ and the moderately large matrix sizes up to $`N=400`$. We were not able to detect any systematic dependence on the size of the matrix $`N`$ in this range.

In the next section we turn to the numerical details of how to compute the level compressibility.

## 3 Results and Discussion

We performed numerical simulations of the RME described by Eqn.(1) for values of $`b^1[0.4,40]`$. The matrices are not sparse and as such the use of fast and efficient algorithms i.e., the Lanczos algorithm is not possible. However, we calculated the spectra of $`N\times N`$ matrices for $`N=50`$ with $`10000`$ realizations to $`N=400`$ with $`2000`$ realizations by exact numerical diagonalization and hence accumulating $`1\times 10^5`$ eigenvalues for each value of $`N`$. The results we discuss here are for the GUE only.

The average density of states as a function of energy $`E`$ is well described by;

$$\rho \left(E\right)=\frac{N}{\pi ^2b}_0^+\mathrm{}\frac{dt}{1+\left(1/z\right)\mathrm{exp}\left[\left(\beta /2\right)\left(E^2+t^2\right)\right]}$$
(6)
where $`z=\mathrm{exp}\left(\pi \beta b\right)1`$. The two limiting cases of Eqn.(6) are;

$$\rho \left(E\right)=\{\begin{array}{cc}\frac{N}{\pi ^2b}\sqrt{2\pi bE^2}\hfill & \text{Semi-circle if }\pi \beta b1\hfill \\ N\sqrt{\frac{\beta }{2\pi }}\mathrm{exp}\left(\beta E^2/2\right)\hfill & \text{Gaussian otherwise}\hfill \end{array}$$
(7)
All our calculations are done at the band center (containing $`n1`$) where the density of states is nearly constant. The calculation of $`\chi `$ from $`\mathrm{\Sigma }_2`$ is non-trivial as there are finite-size effects $`n^2/N`$ in $`\mathrm{\Sigma }_2`$ originating from the long tail $`R(s)1/N`$ of the two-level correlation function. This is circumvented by fitting $`\mathrm{\Sigma }_2`$ data with a polynomial in $`n`$ of degree $`m`$:

$$\mathrm{\Sigma }_2\left(n\right)\chi _0+\chi n+\underset{k=2}{\overset{m}{}}\chi _kn^k$$
(8)
Then we identify the coefficient of the linear term as $`\chi `$ and $`\chi _0`$ is the constant term. Both these terms are expected to have only weak $`N`$-dependence. This is in contrast to that of the coefficients $`\chi _k1/N^{k1}`$ with $`k>1`$ that have a strong $`N`$-dependence which describe finite-size effects in $`\mathrm{\Sigma }_2`$. By making the polynomial fit Eqn.(8) we separate the finite-size effects in $`\mathrm{\Sigma }_2`$ and the weak $`N`$-dependence in $`\chi `$ and in $`\chi _0`$. This weak $`N`$-dependence in $`\chi `$ and in $`\chi _0`$ could be of fundamental interest and needs further investigation.

In Fig.1 we plot $`\mathrm{\Sigma }_2`$ for $`N=400`$ and $`\mathrm{\Sigma }_2\left(n\right)`$ is linear as a function of $`n`$ as is expected for a critical ensemble. The lines through the data points are in agreement with the Poisson-like behavior i.e., $`\mathrm{\Sigma }_2\chi n+\chi _0`$. For large $`b1/2\pi `$ the value of $`\chi `$ is in good agreement with Eqn.(5) and Eqn.(9) that follows from Eqn.(5) and the correlation dimension $`D_2`$ found in Refs. :

$$\chi =1/4\pi b$$
(9)
However, our results are in disagreement with Eqn.(5) for small $`b1/2\pi `$ . The analytic formula for $`\chi `$ (MNS in Fig.2) that follows from  for all values of $`b`$ reads:

$$\chi =\frac{d\mathrm{ln}\left(\right)}{d\mathrm{ln}\left(z\right)}=\frac{\mathrm{PolyLog}[1/2,1\mathrm{exp}\left(2\pi \mathrm{b}\right)]}{\mathrm{PolyLog}[+1/2,1\mathrm{exp}\left(2\pi \mathrm{b}\right)]}$$
(10)
where

$$=_0^{\mathrm{}}\frac{dp}{\left(1+z^1\mathrm{exp}\left(p^2\right)\right)}$$

and $`\mathrm{PolyLog}[\mathrm{n},\mathrm{z}]`$ is an analytic continuation of $`_{k=1}^{\mathrm{}}z^k/k^n`$. This function can easily be computed by use of mathematica.

The main result of this work is that the values of $`\chi `$ (for modestly large $`N<400`$) of the critical ensemble due to  agree unexpectedly well with the theoretical formula due to , Eqn.(10). This agreement is, however, not perfect since the data points lie below Eqn.(10) in the range $`3b^110`$, where the function $`\chi \left(b^1\right)`$ in Eqn.(10) has a peculiar curved shape.

## 4 Conclusion

We have been able to compute $`\chi `$ as a function of $`b^1`$ for the GUE critical ensemble. Our results are in agreement with the theoretical formula Eqn.(10) except for the range $`3b^110`$ where the numerical data are below Eqn.(10).

Furthermore, $`\chi _0`$ and $`\chi `$ are almost $`N`$-independent which is a signature of criticality. It is important to note that our simulations exclude the strong $`N`$-dependence of $`\chi `$ but they do not exclude the weak $`\mathrm{ln}(N)`$-dependence of $`\chi `$ which is not trivial to realize numerically. Recently, it has been shown analytically  that in the CrRME the corrections to $`\chi `$ of order $`b^2`$ have a logarithmic $`N`$-term i.e $`b^2\mathrm{ln}(N)`$ in the crossover between the GOE and the GUE. In both the GUE and the GOE this logarithmic factor is absent up to order $`b^2`$. It’s presence is not excluded in higher $`b`$ corrections, e.g as $`b^3\mathrm{ln}(N)`$. Such a logarithmic dependence is of fundamental interest. This is because the presence of $`\mathrm{ln}(N)`$-corrections would signal the incomplete criticality of the model. The logarithmic terms may drive the system to a different, stable, set of critical points where the $`\chi `$ may be drastically different from the MNS result. Unfortunately, the accuracy of our simulations is not enough to detect such a weak logarithmic dependence. This point needs further investigation.

We note also that the constant term $`\chi _0`$ carries important information about CLS since it has a finite limit as $`N\mathrm{}`$. As is seen in figure 3, the value of $`\chi _0`$ as a function of $`b^1`$ decays to zero with increasing $`b^1`$ which is (so far unnoticed) an important characteristic feature of CLS.

In summary we have been able to compare via numerical simulations of the GUE critical ensemble  and the theoretical formula Eqn.(10). The two results agree but in the region $`3b^110`$, the simulations are below the theoretical calculations. However, we emphasize that our results are preliminary and further work is in progress.

## 5 Acknowledgments

MLN would like to thank R. A. Römer for encouraging discussions, funding from the Deutschen Forschungsgemeinschaft through the Sonderforschungsbereich 393 and the ICTP where this work was started.

## References