text:
# Detecting Very High Energy Neutrinos by the Telescope Array

## 1 Introduction

The detection of very high energy cosmic neutrinos is important from at least three points of view. First, they can carry astrophysical information about point sources which are optically thick at almost all wavelengths and relics of an early inflationary phase in the history of the universe. Second, the connection between the emission of cosmic nuclei, gamma rays, and neutrinos from astrophysical accelerators is of considerable interest with regards the problem of the origin of extragalactic extremely high energy cosmic rays. Third, due to vacuum oscillations, a flux of high energy cosmic neutrinos is expected to be almost equally distributed among the three flavors due to vacuum oscillations, provided that the neutrinos are produced by distant objects such as gamma ray bursts (GRBs) and active galactic nuclei (AGN), and hence the study of oscillation effects on high energy neutrino fluxes can be used to study neutrino mixing and distinguish between different mass schemes .

In this paper we present the potential of the Telescope Array (TA), which will be able to detect cosmic neutrinos with energies above 1016 eV, so far unobserved. They have great potential as probes of astrophysics and particle physics phenomena. They have escaped from regions of dense matter and point back to their sources and they provide a unique window into the most violent events in the universe. Hadrons are accelerated to extremely high energies, regardless of their source, accelerating mechanism, and composition (i.e. $`pp`$ and/or $`p\gamma \pi ^\pm `$), and a subsequent decay of accompanying pions is expected to result in neutrino fluxes through the decay process $`\pi \mu +\nu _\mu e+\nu _e\nu _\mu +\nu _\mu `$, with a resultant ratio $`\nu _\mu /\nu _e`$=2. Observations  of cosmic rays with energies beyond the Greisen-Zatsepin-Kuzmin (GZK) cutoff  raise interest in extremely high energy neutrinos. The GZK mechanism can lead to neutrino fluxes well above those from proposed exotic sources for the super-GZK events around the energy of 1018 eV.

A curious coincidence between the energy flow of the highest energy cosmic rays and gamma ray bursts (GRB) suggests a possible common source . Another speculation is that very high energy cosmic rays may be the result of the annihilation of topological defects left over from the early universe . The energy scales of such events are of the order of 1024 eV (GUT scale). In addition to these ideas is the proposal that the highest energy cosmic rays are evidence for novel particle physics or astrophysics . These super-GZK events provide motivation for examining and confirming predictions of modes of photo production mechanisms. The neutrino fluxes of some representative sources and the range of atmospheric neutrino background as the zenith angle changes from $`0^{}`$ and $`90^{}`$ are given in Fig.1. Note the atmospheric neutrino background  appears to have so large spectral index that it cannot largely contribute to the neutrino detection in the higher energy region above 1016 eV in any source models.

In each of these source cases, neutrino induced air-showers penetrate the atmosphere to a great depth. The imaging technique used in the TA detector has an advantage in this respect as it can observe the slant depth of a deeply penetrating air-shower.

The recent discovery of near-maximal $`\nu _e`$-$`\nu _\mu `$ and $`\nu _\mu `$-$`\nu _\tau `$ mixing  has a significant impact on the detection strategy. The proposed astrophysical sources produce predominantly $`\nu _\mu `$ and $`\nu _e`$ with very small admixtures of $`\nu _\tau `$. However, these sources are so far away that, even at the high energies of interest here, cosmic neutrino fluxes with a $`\nu _e:\nu _\mu :\nu _\tau `$ ratio at the source of $`1:2:0`$, inevitably oscillate in the three neutrino framework to a ratio of $`1:1:1`$, irrespective of the mixing angle.

By looking for upward moving showers from Earth-skimming tau neutrinos, one can test for the existence of neutrino oscillations . Once created in the Earth, a charged lepton loses energy through Bremsstrahlung, pair production, and photonuclear interactions. At the Earth’s surface, with a density of about 2.7 g/cm3, tau leptons and muons are expected to travel 11 km and 1.5 km, respectively, before losing a decade in energy; electrons however do not travel any significant distance. Detectable neutrinos therefore skim the Earth at angles approximately 1 above the horizontal and these events are expected to be predominantly showers induced by tau neutrinos.

Another important signature of very high energy tau neutrinos is the “double bang” which they produce in the atmosphere. The initial shower is produced by the original interaction which creates a tau lepton and a hadronic shower. This is followed by the decay of the tau lepton, producing the second shower bang. The two bangs are separated by a distance of $``$4.9 km (Eτ/1017 eV) where Eτ is the energy of the tau lepton converted from the tau neutrino. This decay length is measurable with the imaging telescopes in the TA detector.

## 2 The TA Detector

The TA detector (Figs.3 and 3) has been designed in order to investigate the origin of the highest energy cosmic rays . For this purpose, the detector is required to obtain statistics of a much greater number of events than the rate of one super-GZK event per year that the AGASA is able to detect. It should also provide particle identification as well as accurate determination of the direction of the primary cosmic ray to test the source models at a high confidence level. The basic measurement principle use a huge effective aperture to detect the image of fluorescence light produced by air showers and from this the longitudinal shower development can be reconstructed.

The TA detector consists of 10 observational stations installed over an area at about 30-40 km intervals as shown in Fig. 3. Each station consists of 40 telescopes with a 3 m-diameter f/1 mirror system on 2 layers of supports. 256 2-inch photo multiplier tubes (PMTs) mounted on the focal plane serve as pixels of the fluorescence sensor of each telescope. Each PMT covers a visual field of angular aperture 1.1$`{}_{}{}^{}\times `$1.0. We expect the detection rate by the TA detector for events with energies beyond the GZK-cutoff to exceed that of AGASA by about 30 times. We note that this is similar to that expected for another proposed detector, the Pierre Auger Observatory (AUGER) . The accuracy of determining the energy and arrival direction for the highest energy events are roughly 25% and 0.6 respectively.

In order for the detector stations to have good visibility, atmospheric transparency, no significant nearby sources of light pollution, and be away from significant traffic, it is preferable to install them on top of small mountains or hills, while still maintaining the arrangement in 30-40 km intervals. We have identified sites in the south-west area of Delta in southern Utah. The planned detector station array covers an area of 200 km2.

## 3 Neutrinos Deeply Penetrating the Atmosphere

Neutrinos can produce air showers through interactions with the atmosphere. The number of the shower depends on the particular interaction. We consider here both inelastic charged and neutral current interactions, which always produce hadronic showers. Note however that the case of charged current electron neutrino interactions, in addition to the emerging electron, a pure electromagnetic shower carrying a large fraction of the incoming particle energy is produced.

In our discussions of neutrino induced showers we make use of new calculations of the cross sections of charged-current and neutral current interactions of neutrinos with nuclei , according to the CTEQ4-DIS (deep inelastic scattering) parton distributions . The CTEQ4-DIS parton distributions take account of new information about parton distributions within the nucleon  using accurate and extensive DIS data from the New Muon Collaboration (NMC)  and DESY ep collider HERA , as well as new data from E665 .

For a neutrino flux $`dI_\nu /dE_\nu `$ interacting through a particular process with differential cross section $`d\sigma /dy`$, where *y* is the fraction of the incident particle energy transferred to the target, the event rate for deeply penetrating showers can be obtained by a simple convolution:

$`Rate[E_\nu >E_{th}]`$ $`=`$
$`N_A\rho _{air}{\displaystyle _{E_{th}}^{\mathrm{}}}`$ $`dE_\nu `$ $`{\displaystyle _0^1}𝑑y{\displaystyle \frac{dI_\nu }{dE_\nu }}(E_\nu ){\displaystyle \frac{d\sigma }{dy}}(E_\nu ,y)ϵ(E_{sh}),`$
where $`N_A`$ is Avogadro’s number and $`\rho _{air}`$ is the air density. The energy integral corresponds to the primary neutrino energy $`E_\nu `$ which is related to the shower energy $`E_{sh}`$, the relationship being different for each interaction. $`ϵ`$ is the detector acceptance, a function of shower energy, which corresponds to the volume and solid angle integrals for different shower positions and orientations with respect to the detector. The function is different for both showers induced by charged current electron neutrino interactions and those arising in neutral current or muon neutrino interactions. This is because hadronic and electromagnetic showers have different particle distribution functions. For $`(\nu _e+\overline{\nu }_e)N`$ charged current interactions, we take the shower energy to be the sum of hadronic and electromagnetic energies, $`E_{sh}=E_\nu `$. For $`(\nu _\mu +\overline{\nu }_\mu )N`$ charged current interactions and for neutral current interactions, we take the shower energy to be the hadronic energy, $`E_{sh}=yE_\nu `$. We take the inelasity $`y`$ to be a function of $`E_\nu `$ as in reference , while $`<y>`$0.25 does not depend strongly on the primary energy beyond $`10^{16}`$ eV in both charged and neutral current interactions.

In the case of $`\nu _\tau `$, we must also take into account the decay length of secondary tau leptons, subsequent tau lepton decay, and a further energy input into the air shower from decaying electrons, photons, and hadrons. The tau decay branching fractions and tau polarization effect are reasonably well accounted for simply by assuming that the energies of neutral mesons fully contribute to the electromagnetic showers. The adopted decay branching fractions and the momentum spectra of the decay with a polarization of $`\pm `$1 in the colinear approximation in the laboratory frame are shown in Table 1 and Fig. 4 respectively. The $`\tau `$ decay length is often long enough that two subsequent energy inputs into the air showers can be seen in the field of view.

Simulations for a given primary cosmic ray (electron neutrino, muon neutrino, tau neutrino, and proton) were performed at fixed energies. The shower energy for neutrinos depends on the generation, as described above. For each shower energy, the mean depth of the proton shower maximum was determined from simulations . For each primary particle at each energy, the mean of the interaction length $`X_1`$ was determined from the above interaction cross sections of neutrinos with nuclei or of protons, $`X_1`$=83.1$`(E/\text{GeV})^{0.052}`$ g/cm2 . The shower energy determines the shower size at a maximum, $`N_{max}`$. Given $`N_{max}`$, a slant depth at the maximum shower size $`X_{max}`$, and $`X_1`$, the complete longitudinal profile can be described by the Gaisser-Hillas function . The Nishimura-Kamata-Greisen (NKG) lateral distribution function  normalized with the Gaisser parameterization has been used to find the total number of electrons and positrons in hadronic (electromagnetic) cascade showers in order to determine the location where fluorescence is produced. In our simulations we took into account the fluctuations of the first interaction depth, impact point, and directional angles of air shower cores, with appropriate distributions, but not air shower size fluctuations.

An important factor that must be calculated is the light yield arriving at the detector site. For this the Rayleigh and Mie scattering processes were simulated, with full account taken of the spectral characteristics of the light. The isotropically emitted fluorescence light, as well as direct and scattered Cerenkov light, is propagated, and the night sky background noise was added to the signal. All processes that affect the overall optical efficiency, mirror area and reflectivity, optical filter transmission, and PMT quantum efficiency factors were included in the light spectrum to give the photoelectron yield in each PMT, due to both signal and noise.

A “fired” PMT is defined such that its instantaneous photoelectron current is greater than the 4$`\sigma `$ noise level of the night sky background. We preselected events such that at least one of the 10 stations contains at least 6 firing PMTs. To ensure track quality, we cut events for which shower maxima were not viewed by any station. Finally, to reject proton background events, we selected only events with $`X_{max}>1700`$ g/cm2. Once the selection cuts were made the detector aperture could be determined and the acceptances calculated by multiplying the appropriate interaction length of a neutrino with a nucleon. Figure 5 shows the shower slant depth distributions from Monte Carlo simulations by which neutrino events can be distinguished from proton events. The results are shown in Fig. 6 for electron, muon, and tau neutrinos as a function of the primary energy. The acceptances of the TA for high-energy neutrinos deeply penetrating air compare well with those of IceCube and AUGER  in lower and higher energy regions respectively, even if we take into account the 10 % duty factor of TA, which is a rate of the observation time at moonless night in good weather.

## 4 Earth-skimming Tau Neutrinos

Very high energy neutrinos penetrate the Earth and convert to charged leptons which then travel through the Earth. This sequence is illustrated for an event with a nadir angle $`\theta `$ in Fig. 7.

We define the critical angle $`\theta _c`$ such that the chord thickness at the nadir angle $`\theta _c`$ corresponds to the charged current interaction length $`L_{CC}^\nu (E_\nu )`$ determined by the interaction cross section for a neutrino traveling with energy $`E_\nu `$. For nadir angles smaller than $`\theta _c`$, neutrinos are shadowed by the Earth, and for larger nadir angles, they rarely interact to produce charged leptons. Table 2 shows charged current cross sections, interaction lengths, and $`90^{}\theta _c`$ at various neutrino energies. For neutrino energies above 1017 eV, $`90^{}\theta _c`$ is small and both the neutrinos and the created leptons travel essentially horizontally.

Once created, a charged lepton loses energy through Bremsstrahlung, pair production, and photonuclear interactions. Electrons lose their energy too quickly in the Earth and hence cannot be detected by fluorescence. Assuming that a lepton loses energy uniformly and continuously, its energy loss can be parametrized by: $`dE_l/dz=\beta _l\rho E_l,`$ where the Earth’s density $`\rho `$ can be simplified to be 2.65 g/cm3 uniformly for small nadir angles. The values of the constant $`\beta _l`$ are $`\beta _\mu 6.0\times 10^6cm^2/g`$ and $`\beta _\tau 0.8\times 10^6cm^2/g`$ for muons and tau leptons respectively for the energies of interest here . Hence, at the Earth’s surface muons and tau leptons can travel 1.5 km and 11 km, respectively, before losing a decade in energy. Therefore, tau neutrinos contribute dominantly to Earth-skimming events observed above the Earth’s surface. Figure 8 shows distributions on the plane of energy $`E_\tau `$ and $`90^{}\theta `$ of tau leptons that exit the Earth from primary neutrinos with energies of 1019 eV and 1017 eV.

The survival probability $`P_{svv}`$ for a tau losing energy as it moves through the Earth is described by the energy loss $`dE_\tau /dz`$ as above, and

$$\frac{dP_{svv}}{dz}=\frac{P_{svv}}{c\tau _\tau E_\tau /m_\tau },$$

where $`c`$ is the speed of light and $`m_\tau `$ and $`\tau _\tau `$ are tau lepton’s rest mass and lifetime respectively. These can be solved and the survival probability given by

$$P_{svv}=\mathrm{exp}[\frac{m_\tau }{c\tau _\tau \beta _\tau \rho }(\frac{1}{E_\nu }\frac{1}{E_\tau })].$$

For tau leptons, this factor plays a significant role. The competition between $`L_{CC}^\nu `$ and $`P_{svv}`$ makes the peak of the distribution at the survived tau lepton energy of $`10^{17}`$ eV in Fig. 8, which is independent of the primary tau neutrino energy. Due to this effect, lower energy sensitivity of an air fluorescence telescope makes a key role to detect the Earth skimming tau neutrinos.

The flux of survived taus with energy $`E_{min}<E_\tau <E_{max}`$ is

$$\mathrm{\Phi }_\tau =\frac{ln(E_{max}/E_{min})}{2R\beta _\tau \rho }\mathrm{\Phi }_\nu ,$$

where $`\mathrm{\Phi }_\nu `$ is the flux of neutrinos. As a result, the flux in any given decade in tau energy is $`\mathrm{\Phi }_\tau =8.5\times 10^4\mathrm{\Phi }_\nu `$, that is, 1 in every 1200 neutrinos that skims the Earth emerges as a tau lepton with the required energy . This remarkably simple and robust statement is quite useful because it is independent of the neutrino energy and microscopic details. We have examined our Monte Carlo simulation results by comparison with this analytic estimate.

The tau decay length is given by $`L_\tau =c\tau _\tau (E_\tau /m_\tau )`$ 4.9 km ($`E_\tau /10^{17}`$ eV) , using the lifetime $`c\tau _\tau =87.11\mu `$m the tau mass $`m_\tau =1777.03`$ MeV from PDG . The complete treatment for tau decays can be performed as described in the previous section. Effective detection apertures and acceptances have been estimated for Earth-skimming tau events changing the neutrino energies from the detailed Monte Carlo simulation. Figure 9 shows the effective detection acceptances using only one TA station. The tau decay length as a function of the tau lepton energy which is observed within the view of the telescope makes the detection acceptance slightly decrease as increasing the kinetic energy of Earth skimming tau neutrino above the primary energy of 1019 eV.

## 5 Detection Event Rates

Models that assumes shock acceleration in the Active Galactic Neuclei (AGN) cores predict relatively flat fluxes up to energies of about 1015 eV . The GeV to TeV gamma ray emissions observed in AGN corresponds to the blazar class. Most recent models for the proton blazars site the acceleration in the jets themselves. We use the prediction of reference , which illustrates that the emitted neutrinos may extend well into the EeV region. Annual event rates have been calculated for the TA detector for neutrino induced air showers with fluxes from GRB model  and Greisen neutrino  models as well as from AGN-core  and AGN-jet , assuming a duty cycle of 10%, and are shown in Table 3. The testability of the bottom-up scenario models is statistically simple even when considering only shower events deeply penetrating the air.

We also compare the detectability of Earth-skimming tau neutrino events with the downward events described above. Table 3 also lists annual event rates for neutrino induced air-showers from various models including top-down scenario sources of shower events deeply penetrating the air and Earth-skimming events. The duty factor of 10% is assumed as usual. The enhancement of the statistics by considering Earth-skimming events is significant, particularly in the ultra-high energy region. Hence it is feasible to test the top-down scenario models such as topological defects and Z-bursts in a statistically significant way by considering both downward and upward neutrino-induced air-shower events in complementary energy regions.

## 6 Conclusions

The next generation of cosmic ray detectors such as AUGER and the TA have target volumes of atmosphere competitive with or better than that of IceCube and their annual detection rates assuming AGN-jet proton acceleration models are statistically sizable . The TA using an air fluorescence technique is advanced as a method for discriminating primary neutrinos proton induced air showers. Recently a novel strategy for detecting Earth-skimming extremely high energy neutrinos has been quantitatively proposed for which significant improvements in detection can be made using the TA detector compared to down-going neutrino detection . Encouragingly, extremely high energy cosmic neutrino sources, to say nothing of AGN proton acceleration models, such as Greisen photo production , topological defects , long-lived super heavy particles, and Z-bursts , may be experimentally tested with the TA.

## 7 Acknowledgments

We are indebted to my colleagues in the TA Collaboration, especially those who have contributed to the TA Design Report, for assistance in the preparation of this paper. We would like to especially acknowledge the contributions of Prof. T. Matsuda and Prof. C. Fukunaga for useful discussions.


format: text

text:
# JOURNAL OF THE PHYSICAL SOCIETY OF JAPAN 70 (2001) 2833 Anisotropic optical conductivity of Nd2-xCexCuO4 thin films

## I Introduction

Oxide high-$`T_c`$ superconductors have been investigated intensively over the last decade.?? $`d`$-wave superconductivity is now well established for hole-doped superconductors. However, there is a class of high-$`T_c`$ superconductors doped with electrons,?? for which both $`s`$-wave? and $`d`$-wave pairing??? have been reported. Nd2-xCexCuO4 is a typical example of electron-doped materials and the symmetry of Cooper pairs has been a controversial issue. It is important to examine the symmetry of Cooper pairs in the study of high-$`T_c`$ superconductors.

Far-infrared (FIR) spectroscopy is a powerful technique to investigate the nature of the superconducting gap. The conventional FIR spectroscopy based on a Kramers-Kronig transformation, however, is rather unfavorable for studying electronic properties in the low energy regime of the FIR region.?? Since the superconducting gap $`\mathrm{\Delta }`$ in Nd2-xCexCuO4 is very small, there have been no reports on the study of the nature of the superconducting gap of Nd2-xCexCuO4 through such techniques, although there have been a number of reports on the FIR spectroscopy of Nd2-xCexCuO4.?? Recently, a new method to examine FIR spectroscopy has been developed without the need for evaluating the Kramers-Kronig transformation.? In this method, the optical conductivity is estimated from the data of reflectance spectra $`R(\omega )`$ and transmittance spectra $`T(\omega )`$ by substituting them into a set of coupled equations. The new method is free from the conventional difficulties in the FIR region since we do not need the aid of the Kramers-Kronig transformation.? This method is referred to as the R-T method in this paper.

The purpose of this paper is to investigate FIR optical properties of Nd2-xCexCuO4 obtained by the R-T method from the viewpoint of unconventional superconductors. We will show that the available data for the optical conductivity and transmittance are well explained by the $`d`$-wave pairing model in the clean limit. The value of the superconducting gap is estimated as $`2\mathrm{\Delta }5060`$ cm-1, which is consistent with the available value estimated by scanning tunneling spectroscopy.?

## II Theory and Measured Conductivity

The frequency-dependent conductivity $`\sigma (\omega )`$ has been calculated by Mattis and Bardeen,? Abrikosov et al.? and Skalski et al.? for isotropic superconductors. The original Mattis-Bardeen theory was tested for a conventional type-I $`s`$-wave superconductor, where the coherence length $`\xi `$ and magnetic penetration depth $`\lambda `$ satisfy $`\xi \lambda `$. The opposite limit $`\xi \lambda `$ (London limit) was also examined for $`s`$-wave pairing by field theoretical treatments.? For the high-$`T_c`$ compounds of type-II superconductor with small coherence length, the formula in the London limit is appropriate for optical conductivity measurements. Recently, the conductivity $`\sigma (\omega )`$ of an unconventional superconductor has been examined theoretically in the London limit.????? The current response function, from which the optical conductivity is derived, is given by?

$$K_{ij}(𝐪,\omega _m)=\frac{e^2k_F^2}{m^2c}\underset{k}{}\widehat{k}_i\widehat{k}_j\frac{1}{\beta }\underset{n}{}Tr[G(𝐤_+,ϵ_n+\omega _m)G(𝐤_{},ϵ_n)],$$
(1)
where $`𝐤_\pm =𝐤\pm 𝐪/2`$ and $`ϵ_n=(2n+1)\pi /\beta `$. The single-particle matrix Green’s function is given by

$$G(𝐤,ϵ_n)=\frac{i(ϵ_n\mathrm{\Sigma }(ϵ_n))\tau ^0+\xi _k\tau ^3+\mathrm{\Delta }_k\tau ^1}{(ϵ_n\mathrm{\Sigma }(ϵ_n))^2+\xi _k^2+\mathrm{\Delta }_k^2},$$
(2)
where $`\mathrm{\Delta }_k`$ is the anisotropic order parameter and $`\mathrm{\Sigma }(ϵ_n)`$ is the self-energy due to impurity scattering. $`\tau ^i`$ ($`i=0,1,\mathrm{}`$) denote Pauli matrices. Since we consider the case where $`\xi \lambda `$ holds, the real part of optical conductivity is well approximated by the formula in the London limit:

$$\sigma _{ij}(\omega )=\frac{c}{\omega }\underset{q0}{lim}\mathrm{Im}K_{ij}(𝐪,\omega ).$$
(3)
The expressions of $`K_{ij}(𝐪,\omega )`$ for real and imaginary parts were given in ref. 20. Our focus is the collisionless limit of the normalized conductivity to compare it with the data for Nd2-xCexCuO4 since $`\xi \mathrm{}`$ holds for the mean-free path $`\mathrm{}`$. For anisotropic superconducting order parameter $`\mathrm{\Delta }_k`$ such that the average over the Fermi surface vanishes $`\mathrm{\Delta }_k=0`$, the expression for $`\sigma _{xx}(\omega )`$ in the collisionless limit on the plane is simply given by?

$$\frac{\sigma _{xx}(\omega )}{\sigma _n}=\frac{1}{2\omega }_{\mathrm{}}^{\mathrm{}}𝑑x\mathrm{Re}\frac{|x|}{(x^2\mathrm{\Delta }_k^2)^{1/2}}\mathrm{Re}\frac{|x\omega |}{[(x\omega )^2\mathrm{\Delta }_k^2]^{1/2}}[\mathrm{tanh}(\frac{\beta x}{2})\mathrm{tanh}(\frac{\beta (x\omega )}{2})],$$
(4)
which is an angle-dependent generalization of the Mattis-Bardeen formula. For the $`d`$-wave symmetry, the average over the Fermi surface denoted by the angular brackets is defined as

$$\mathrm{Re}\frac{x}{(x^2\mathrm{\Delta }_k^2)^{1/2}}=\mathrm{Re}\frac{d\varphi }{2\pi }\frac{x}{[x^2(\mathrm{\Delta }\mathrm{cos}(2\varphi ))^2]^{1/2}},$$
(5)
where the order parameter is factorized as $`\mathrm{\Delta }_k=\mathrm{\Delta }\mathrm{cos}(2\varphi )`$. In Fig.1, we show the behaviors of $`\sigma _s(\omega )\sigma _{xx}(\omega )`$ as a function of $`\omega `$ for several values of temperature $`T`$. The FIR behaviors reflect the lines of nodes on the Fermi surface.

FIR reflection $`R(\omega )`$ and transmission $`T(\omega )`$ measurements were performed for Nd2-xCexCuO4 ($`x=0.15`$) thin films deposited by laser ablation onto (001) MgO substrates. The thickness of the Nd2-xCexCuO4 thin film was about 40 nm. $`T_c`$ was estimated to be $``$20K. The electric field of the FIR radiation was predominantly parallel to the $`a`$-$`b`$ plane. The conductivity spectra were evaluated by the R-T method from the data for $`R(\omega )`$ and $`T(\omega )`$ at $`T=4.3`$ and 30K.?

The R-T method provides us with reliable spectroscopic data in the FIR region due to which a comparison between the experimental data and that of theoretical analysis is possible. In the R-T method, both the refrectance spectra $`R(\omega )`$ and the transmittance spectra $`T(\omega )`$ are measured experimentally from which a set of coupled equations is derived describing the transmittance and reflectance of a thin film on a substrate. The coupled equations are solved numerically by the Newton method to determine the optical conductivity. This method is free from the difficulties in the FIR region which occur commonly in the conventional method employing a Kramers-Kronig transformation. In Fig.2, we show the observed data and theoretical curves at $`T=0`$ for $`2\mathrm{\Delta }=50`$, 60 and 70 cm-1. The experimental data $`\sigma _1(3.4K)/\sigma _1(30K)`$ normalized by the normal state values at $`T=30`$K are shown in Fig.2. It is obvious from the experimental results that there is no evidence of a true gap, which is suggestive of an anisotropic superconducting gap, since the spectral weight of conductivity should vanish for $`\omega 2\mathrm{\Delta }`$ at $`T=0`$ in conventional isotropic superconductors. It is also shown in Fig.2 that they are well fitted by the curve with $`2\mathrm{\Delta }=50`$ cm-1, which is consistent with the value estimated by scanning tunneling spectroscopy measurements.?

A transmission curve is also presented in Fig.3, where $`T_S/T_N`$, the ratio of the transmission in the superconducting ($`T=`$3.4K) to that in the normal state ($`T=`$30K), is the experimentally measured quantity. The following phenomenological expression for $`T_S/T_N`$ is employed to determine the transmission curve theoretically,?

$$\frac{T_S}{T_N}=\frac{1}{[T_N^{1/2}+(1T_N^{1/2})(\sigma _1/\sigma _n)]^2+[(1T_N^{1/2})(\sigma _2/\sigma _n)]^2},$$
(6)
where $`\sigma _1`$ and $`\sigma _2`$ are real and imaginary parts of the conductivity $`(c/\omega )K(𝐪,\omega )`$ for $`𝐪0`$, respectively. $`T_N`$ is determined as $`T_N0.08`$ from the expression for the ratio of the power transmitted with a film to that with no film given as

$$T_N=1/[1+\sigma _nd\frac{Z_0}{n+1}]^2.$$
(7)
Here, $`d`$ is the film thickness, $`n`$ is the index of refraction of the substrate, and $`Z_0`$ is the impedance of free space. We have assigned the following values; $`d=4\times 10^6`$cm, $`n=3.13`$, $`Z_0=377\mathrm{\Omega }`$, $`\sigma _n7\times 10^3\mathrm{\Omega }^1`$cm-1. Apparently, the $`\omega `$-dependence of measured transmittance agrees with the theoretical curve for $`2\mathrm{\Delta }=50`$60cm-1 as shown in Fig.3. The agreement between the observed quantities and theoretical curve is significant, which should be compared to the isotropic BCS prediction calculated from the Mattis-Bardeen equations.?

## III Summary

We have successfully made a comparison between experimental data and theoretical data for the optical conductivity of Nd2-xCexCuO4 in the FIR region for the first time. We have shown that there is a reasonable agreement between the optical conductivity $`\sigma _1(\omega )`$ observed by the R-T method and theoretical analysis without adjustable parameters except the superconducting gap. An estimate of 50$``$60 cm-1 for the superconducting gap is consistent from both the experimental and theoretical aspects. The FIR optical conductivity suggests that the superconducting gap of electron-doped Nd2-xCexCuO4 is an unconventional one with nodes on the Fermi surface. The anisotropic nature of electron-doped superconductors is consistent with the recent researches performed for the one-band and three-band Hubbard models.???? If the superconducting gap is anisotropic for the electron-doped superconductors, there is a possibility that both the hole-doped and electron-doped cuprates superconductors are governed by the same superconductivity mechanism.


text:
# The fractality of the relaxation modes in deterministic reaction-diffusion systems

## I Introduction

Recently, many studies have been devoted to fractal structures induced by the chaotic dynamics in the phase-space of non-equilibrium statistical systems and, in particular, to the fractal character of the hydrodynamic modes governing the exponential relaxation of the system to the thermodynamic equilibrium  . This fractal character has been related to the entropy production in the approach of equilibrium, for a diffusive multibaker model . Previous works had shown similar results for the entropy production in non-equilibrium steady states . On the other hand, the fractal dimension characterizing the hydrodynamic modes of diffusion has been related to the diffusion coefficient and to the positive Lyapunov exponent by a relation, first obtained for diffusive multibaker maps , and proved for general chaotic systems with two degrees of freedom . The construction of the hydrodynamic modes in Refs.  gives us a new approach to relate macroscopic transport coefficients and microscopic chaotic properties. It establishes a link between the irreversible transport processes and the reversible microscopic dynamics.

Two other approaches are the so-called escape rate formalism and the thermostated system approach. The first approach can be applied to open systems with absorbing boundaries . The escape of trajectories leads to the formation of a fractal repeller, consisting in the set of orbits forever trapped within the absorbing boundaries. Under such conditions, the transport properties can be related to the positive Lyapunov exponents and the Kolmogorov-Sinai entropy per unit time as well as to fractal dimensions . In the second approach, an external field is applied to the system, and in order to keep the kinetic energy constant, a special force acts on the particles as a heat pump . The system is still time-reversible but does not preserve phase-space volumes anymore. The trajectories converge to a fractal attractor. In thermostated systems, the transport properties have been related to the sum of Lyapunov exponents . In the new approach of Refs. , the advantage is that the fractal curve considered is directly related to the hydrodynamic modes of relaxation towards thermodynamic equilibrium. It does not require absorbing boundaries neither a thermostat.

In the present paper, we extend the new approach of Refs.  to reaction-diffusion processes. Reaction-diffusion processes are of fundamental importance to understand self-organization in physico-chemical systems . In chemical systems, the microscopic mechanisms of the relaxation toward the thermodynamic equilibrium are still poorly understood already for simple reactions. For the purpose of contributing to this important question, we consider here a simple reactive process of isomerization. The particle is supposed to perform a deterministic motion among static scatterers and to carry a color $`A`$ or $`B`$. When colliding on some special scatterers, playing the role of catalysts, its color changes instantaneously with a given probability $`p_0`$,

$`A+\text{catalytic scatterer}`$ $``$ $`B+\text{catalytic scatterer},\text{with probability}p_0,`$ (1)
$`A+\text{other scatterer}`$ $``$ $`A+\text{other scatterer},`$ (2)
$`B+\text{other scatterer}`$ $``$ $`B+\text{other scatterer}.`$ (3)
For this class of systems, we obtain an expression of the Hausdorff dimension of the reactive modes of relaxation in terms of the reactive dispersion relation and of a function $`Q(\alpha ,\beta )`$ generalizing the Ruelle topological pressure $`P(\beta )`$. In the long wavelength limit, we then infer a relation between the Hausdorff dimension of the modes, the reactive diffusion coefficient, the reaction rate, and two derivatives of the function $`Q(\alpha ,\beta )`$. In the limit $`p_00`$, we recover a relation previously derived for the diffusive case . This new expression relates the macroscopic transport and reaction processes to the microscopic underlying dynamics.

The paper is organized as follows. In Section II, we generalize the work of Refs.  to the reactive case considered here. In Section III, we test our relation on a reactive multibaker model. We next study the case of a two-dimensional periodic reactive Lorentz gas in Section IV. Conclusions are drawn in Section V.

## II Relaxation modes: general results

We consider a system such as a Lorentz gas in which the particles move independently of each other along trajectories of a deterministic dynamical system given by

$$\frac{d𝐗}{dt}=𝐅(𝐗)$$
(4)
The variables of this system are for instance the positions and velocities of the particle: $`𝐗=(𝐫,𝐯)`$. The vector field $`𝐅(𝐗)`$ is defined in this space. The system is spatially periodic in the positions $`𝐫`$. The particle moves in a periodic lattice of scatterers given by a certain potential of interaction.

Moreover, each particle carries a color which is either $`A`$ or $`B`$. The color changes when the particle interacts with some of the scatterers, called the catalytic scatterers according to the reaction (3). We notice that we may also consider that the particle has a spin one-half which flips upon collision on special scatterers.

We can model the reactive events by a change of color with probability $`p_0`$ when the trajectory meets a certain hypersurface $`\mathrm{\Sigma }(𝐗)=0`$ defined in the space of the variables $`𝐗`$. This hypersurface surrounds each catalytic scatterer and constitutes the locus where the reaction occurs. The introduction of the color has for consequence that the phase space is composed of two copies of the space of the variables $`𝐗`$, which are glued together along the hypersurface $`\mathrm{\Sigma }`$. This hypersurface can be viewed as the gate from copy corresponding to the color $`A`$ to the one corresponding to color $`B`$ and vice versa. These two copies now form two subspaces of the whole phase space of the reactive system. In this framework, a reactive event corresponds to the crossing of the hypersurface $`\mathrm{\Sigma }`$ from one color subspace to the other.

A statistical ensemble of particles of both colors is described by two probability densities $`q(𝐗,A,t)`$ and $`q(𝐗,B,t)`$, to find a particle at the point $`𝐗`$ with color $`c=A`$ or $`B`$, at the current time $`t`$. The time evolution of these probability densities is ruled by two coupled Liouville equations:

$$\{\begin{array}{cc}_tq(𝐗,A,t)+\mathrm{div}\left[𝐅(𝐗)q(𝐗,A,t)\right]=p_0\mathrm{\Xi }\left[\mathrm{\Sigma }(𝐗)\right]\left[q(𝐗,A,t)q(𝐗,B,t)\right],\hfill & \\ _tq(𝐗,B,t)+\mathrm{div}\left[𝐅(𝐗)q(𝐗,B,t)\right]=+p_0\mathrm{\Xi }\left[\mathrm{\Sigma }(𝐗)\right]\left[q(𝐗,A,t)q(𝐗,B,t)\right].\hfill & \end{array}$$
(5)
where $`\mathrm{\Xi }`$ is a special distribution defined as

$$\mathrm{\Xi }\left[\mathrm{\Sigma }(𝐗)\right]h(𝐗)\underset{ϵ0^+}{lim}\left|𝐅(𝐗)\frac{\mathrm{\Sigma }}{𝐗}(𝐗)\right|\delta \left[\mathrm{\Sigma }(𝐗)\right]h\left[𝐗ϵ𝐅(𝐗)\right]$$
(6)
$`\delta `$ being the Dirac distribution. The distribution (6) expresses the fact that the particles change their color upon crossing of the hypersurface $`\mathrm{\Sigma }`$. Accordingly, the probability densities just after the hypersurface $`\mathrm{\Sigma }`$ in the direction of the vector field $`𝐅`$ are the following linear combinations of the densities just before $`\mathrm{\Sigma }`$:

$$\{\begin{array}{cc}q(𝐗_+,A,t)=(1p_0)q(𝐗_{},A,t)+p_0q(𝐗_{},B,t),\hfill & \\ q(𝐗_+,B,t)=p_0q(𝐗_{},A,t)+(1p_0)q(𝐗_{},B,t),\hfill & \end{array}$$
(7)
where $`𝐗_\pm =𝐗_{}\pm ϵ𝐅(𝐗_{})`$ with $`\mathrm{\Sigma }(𝐗_{})=0`$ and an arbitrarily small $`ϵ>0`$. These linear combinations can be expressed in terms of the transmission matrix

$$𝖳=\left(\begin{array}{cc}1p_0& p_0\\ p_0& 1p_0\end{array}\right).$$
(8)
We denote by $`𝐗=𝚽^t𝐗_0`$ the trajectory of the differential equation (4) from the initial condition $`𝐗_0`$. Moreover, we gather both probability densities in the vector

$$𝐪(𝐗,t)=\left(\begin{array}{c}q(𝐗,A,t)\\ q(𝐗,B,t)\end{array}\right).$$
(9)
The solution of the coupled Liouville equations (5) from the initial probability densities $`𝐪_0`$ can then be written in the form of the following Frobenius-Perron operator:

$$𝐪(𝐗,t)=\left(\widehat{𝖯}^t𝐪\right)(𝐗)=\frac{𝖳^{N_t(𝐗)}𝐪_0\left(𝚽^t𝐗\right)}{\left|\frac{𝚽^t}{𝐗}(𝚽^t𝐗)\right|}$$
(10)
with the matrix (8) and where $`N_t(𝐗)`$ is the number of crossings through the hypersurface $`\mathrm{\Sigma }`$ performed by the segment of trajectory from $`𝐗_0=𝚽^t𝐗`$ to $`𝐗=𝚽^t𝐗_0`$. If the flow of Eq. (4) is volume preserving, we have that $`\mathrm{div}𝐅=0`$ which implies that $`\left|\frac{𝚽^t}{𝐗}\right|=1`$, which is the case for the systems considered in the following.

The system of Liouville’s equations (5) defines a random process and it is important to characterize the dynamical randomness generated by this process. For this purpose, we calculate its Kolmogorov-Sinai entropy per unit time which is the amount of information produced per unit time by the process. The Kolmogorov-Sinai entropy is calculated by partitioning the whole phase space into small cells and by computing the probability for a trajectory to visit successively different cells of the partition . If we suppose that the differential equation (4) defines a chaotic dynamical system with Lyapunov exponents $`\{\lambda _i\}`$, the dynamical instability of the deterministic motion is a first contribution to the Kolmogorov-Sinai entropy given by the sum of positive Lyapunov exponents according to Pesin’s theorem . A second contribution comes from the random change of color upon crossing of the hypersurface $`\mathrm{\Sigma }`$. At each crossing, the particle jumps from one color subspace to the other with probability $`p_0`$. If the number of crossings is equal to $`N_t`$ during the lapse of time $`t`$, the frequency of $`\mathrm{\Sigma }`$ crossings is equal to :

$$\nu \underset{t\mathrm{}}{lim}\frac{N_t}{t},$$
(11)
where $``$ denotes the average with respect to the equilibrium state. At each crossing, we have an independent dichotomous random variable with probabilities $`(p_0,1p_0)`$, which defines a Bernoulli process with the rate $`\nu `$. Accordingly, the overall Kolmogorov-Sinai entropy per unit time of the process (5) is given by

$$h_{\mathrm{KS}}=\underset{\lambda _i>0}{}\lambda _i+\nu \left[p_0\mathrm{ln}\frac{1}{p_0}+(1p_0)\mathrm{ln}\frac{1}{1p_0}\right].$$
(12)
Pesin’s expression of the Kolmogorov-Sinai entropy is recovered in the deterministic limits $`p_0=0`$ and $`p_0=1`$.

For $`0<p_0<1`$, we remark that the Kolmogorov-Sinai entropy (12) remains positive and finite as for a fully deterministic dynamical system. The way the reaction is modelled here is therefore compatible with a deterministic dynamics that would take place near each catalytic scatterer on a time scale much shorter than the free flight between each scatterer. We can imagine that the reaction follows a complicated trajectory on such a short time scale in a complicated potential, which is here simplified by capturing the complexity of the reaction into a sole transition probability $`p_0`$. Still the so-defined random process has a finite Kolmogorov-Sinai entropy contrary to stochastic processes which have an infinite Kolmogorov-Sinai entropy per unit time. Indeed, it has been shown elsewhere  that stochastic processes such as the birth-and-death processes or the linearized Boltzmann processes have a dynamical randomness on arbitrarily small scales $`ϵ`$ with an $`ϵ`$-entropy per unit time increasing as $`h(ϵ)\mathrm{ln}(1/ϵ)`$. Stochastic processes defined by Langevin stochastic differential equations have an even larger dynamical randomness with an $`ϵ`$-entropy per unit time increasing as $`h(ϵ)1/ϵ^2`$. The Kolmogorov-Sinai entropy is the supremum of the $`ϵ`$-entropy in the limit $`ϵ0`$ and is therefore infinite. In contrast with such stochastic processes, the present random process has a finite Kolmogorov-Sinai entropy. This dynamical randomness is comparable with the randomness of the classical dynamics on graphs which has been recently characterized . For the classical graphs, the motion is one-dimensional so that the only source of randomness is the random bifurcations at the vertices of the graph. Here, there is an extra source of randomness which is the dynamical instability due to the chaotic motion of the particle in the multi-dimensional phase space of Eq. (4).

The overall dynamical randomness (12) induces a relaxation of the probability densities toward the thermodynamic equilibrium. In order to study this relaxation we have to find the relaxation rates which are given by the eigenvalues of the Frobenius-Perron operator (10). It turns out that the Liouville equations (5) can be decoupled if we introduce the probability density $`f`$ to find a particle at the point $`𝐗`$ no matter which is its color, and the difference $`g`$ between both color densities:

$`f(𝐗,t)`$ $`=`$ $`q(𝐗,A,t)+q(𝐗,B,t)`$ (13)
$`g(𝐗,t)`$ $`=`$ $`q(𝐗,A,t)q(𝐗,B,t).`$ (14)
These new functions are gathered in the vector

$$𝐪_\mathrm{d}(𝐗,t)=\left(\begin{array}{c}f(𝐗,t)\\ g(𝐗,t)\end{array}\right)$$
(15)
which is related to the previous vector (9) by the relation

$$𝐪_\mathrm{d}=𝖲𝐪$$
(16)
in terms of the matrix

$$𝖲=\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).$$
(17)
This matrix diagonalizes the transition matrix (8) into

$$𝖳_\mathrm{d}=𝖲𝖳𝖲^1=\left(\begin{array}{cc}1& 0\\ 0& 12p_0\end{array}\right).$$
(18)
Accordingly, the Liouville equations (5) decouple into

$$\{\begin{array}{cc}_tf(𝐗,t)+\mathrm{div}\left[𝐅(𝐗)f(𝐗,t)\right]=0,\hfill & \\ _tg(𝐗,t)+\mathrm{div}\left[𝐅(𝐗)g(𝐗,t)\right]=2p_0\mathrm{\Xi }\left[\mathrm{\Sigma }(𝐗)\right]g(𝐗,t).\hfill & \end{array}$$
(19)
which shows that the probability density $`f`$ obeys the deterministic Liouville equation while the difference $`g`$ solely controls the reaction. Similarly, the Frobenius-Perron operator (10) splits in the two following operators:

$`f(𝐗,t)`$ $`=`$ $`{\displaystyle \frac{f_0\left(𝚽^t𝐗\right)}{\left|\frac{𝚽^t}{𝐗}(𝚽^t𝐗)\right|}}`$ (20)
$`g(𝐗,t)`$ $`=`$ $`(12p_0)^{N_t(𝐗)}{\displaystyle \frac{g_0\left(𝚽^t𝐗\right)}{\left|\frac{𝚽^t}{𝐗}(𝚽^t𝐗)\right|}}.`$ (21)
An eigenvalue problem can now be considered separately for both operators (20) and (21). If we denote them generically by $`\widehat{P}^t`$, we are looking for eigenstates given by some distributions which are quasiperiodic in the position space $`𝐫`$ with a wavenumber $`𝐤`$:

$$\mathrm{\Psi }_𝐤(𝐗)=e^{i𝐤𝐫(𝐗)}\mathrm{{\rm Y}}_𝐤(𝐗)$$
(22)
where $`\mathrm{{\rm Y}}_𝐤`$ is periodic in $`𝐫`$. The wavenumber $`𝐤`$ describes spatial inhomogeneities of wavenumber $`\mathrm{}=2\pi /𝐤`$. The distribution (22) is supposed to be a solution of the eigenvalue problem

$$\widehat{P}^t\mathrm{\Psi }_𝐤=e^{s_𝐤t}\mathrm{\Psi }_𝐤$$
(23)
associated with the operator $`\widehat{P}^t`$. The quantity $`s_𝐤`$ is the relaxation rate of the eigenstate (22).

The idea is that these eigenstates can be obtained in principle by applying the operator $`\widehat{P}^t`$ over a long time to some expression (22) and by renormalizing the result in order to compensate the relaxation according to

$$\mathrm{\Psi }_𝐤(𝐗)=\underset{t\mathrm{}}{lim}e^{s_𝐤t}\widehat{P}^t\left[e^{i𝐤𝐫(𝐗)}\mathrm{{\rm Y}}(𝐗)\right]$$
(24)
even starting from some arbitrary function $`\mathrm{{\rm Y}}`$. This renormalization-semigroup idea has been developed in detail for diffusion in Refs. . For chaotic dynamical systems (4), such eigenstates turn out to be singular distributions without a density function but with cumulative functions defined by

$$F_𝐤(\theta )=_0^\theta 𝑑\theta ^{}\mathrm{\Psi }_𝐤\left[𝐗(\theta ^{})\right]$$
(25)
where $`\theta `$ is a parameter (such as an angle) of a curve $`𝐗(\theta )`$ defined in the phase space.

The hydrodynamic modes of diffusion as well as the reactive modes can be identified as the eigenstates corresponding to the slowest relaxation rates and they can be constructed thanks to the above idea of renormalization semigroup, as we shall show below.

### A Diffusive modes

In order to construct a diffusive mode of wavenumber $`𝐤`$ we consider the Frobenius-Perron operator (20) and an initial distribution $`\mathrm{{\rm Y}}(𝐗)`$ of particles of initial positions $`𝐫_0`$. We integrate the trajectories by using the equations of motion (4) in order to obtain the position $`𝐫_t`$ after a long time $`t`$. As a consequence of the aforementioned idea of renormalization semigroup, the relaxation rate of the diffusive mode is given by the decay rate of the Van Hove intermediate incoherent scattering function  as

$$s_𝐤=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\mathrm{exp}[i𝐤(𝐫_t𝐫_0)]$$
(26)
where $``$ is an average over initial conditions distributed according to some distribution $`\mathrm{{\rm Y}}`$ . In a similar way, the cumulative function of the diffusive mode can be obtained as

$$F_𝐤(\theta )\underset{t\mathrm{}}{lim}\frac{_0^\theta 𝑑\theta ^{}\mathrm{exp}\{i𝐤[𝐫_t(\theta ^{})𝐫_0(\theta ^{})]\}}{_0^{2\pi }𝑑\theta ^{}\mathrm{exp}\{i𝐤[𝐫_t(\theta ^{})𝐫_0(\theta ^{})]\}}$$
(27)
where the initial conditions are distributed according to a distribution $`\mathrm{{\rm Y}}`$ which is here taken as a uniform distribution along a one-dimensional line parametrized by the angle $`\theta [0,2\pi [`$ .

For two-degree-of-freedom chaotic systems with one positive Lyapunov exponent $`\lambda `$, the cumulative function (27) forms a fractal curve in the complex plane and its Hausdorff dimension is given in terms of the diffusion coefficient and the Lyapunov exponent according to

$$D_\mathrm{H}(𝐤)=1+\frac{𝒟}{\lambda }k^2+𝒪(k^4)$$
(28)
as shown in Refs. . The derivation of Ref.  is here summarized because we need it for the generalization to the reactive modes.

We assume that we can use the result by Sinai, Bowen and Ruelle that averages can be performed in terms of unstable trajectories covering the phase space within a certain resolution . The probability weight given to each of these trajectories is inversely proportional to their stretching factor $`\mathrm{\Lambda }_t^{(i)}`$ which characterizes their dynamical instability. Accordingly, Eq. (26) can be transformed into the condition

$$\underset{j}{}|\mathrm{\Lambda }_t^{(j)}|^1\mathrm{exp}(s_𝐤t)\mathrm{exp}\{i𝐤[𝐫_t^{(j)}𝐫_0^{(j)}]\}_t\mathrm{}1.$$
(29)
If the sum in Eq. (29) is restricted to the trajectories issued from initial conditions in the interval $`[0,\theta ]`$, we obtain at time $`t`$ a polygonal approximation of the cumulative function (27): the integral $`_0^\theta `$ is an average over trajectories with initial conditions in $`[0,\theta ]`$ and the denominator is proportional to the factor $`\mathrm{exp}(s_𝐤t)`$. At time $`t`$, the curve $`(\mathrm{Re}F_𝐤,\mathrm{Im}F_𝐤)\text{}`$ is thus approximated by a polygon of sides given by the small complex vectors

$$\mathrm{\Delta }F^{(j)}=|\mathrm{\Lambda }_t^{(j)}|^1\mathrm{exp}(s_𝐤t)\mathrm{exp}\left\{i𝐤\left[𝐫_t^{(j)}𝐫_0^{(j)}\right]\right\}.$$
(30)
Each side has the length

$$\epsilon _j=|\mathrm{\Delta }F^{(j)}|=|\mathrm{\Lambda }_t^{(j)}|^1\mathrm{exp}(\mathrm{Re}s_𝐤t).$$
(31)
In the limit $`t\mathrm{}`$, this polygon converges to a fractal curve, characterized by a Hausdorff dimension given by $`_j\epsilon _j^{D_\mathrm{H}}1`$. Accordingly, the Hausdorff dimension of the hydrodynamic mode should satisfy the condition

$$\underset{j}{}|\mathrm{\Lambda }_t^{(j)}|^{D_\mathrm{H}}\mathrm{exp}(D_\mathrm{H}\mathrm{Re}s_𝐤t)_t\mathrm{}1.$$
(32)
On the other hand, Ruelle’s topological pressure is defined in dynamical systems theory by

$$P(\beta )\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}|\mathrm{\Lambda }_t|^{1\beta }$$
(33)
where the average $``$ is carried out over the equilibrium invariant measure . The mean positive Lyapunov exponent of the system is given by $`\lambda =dP/d\beta |_{\beta =1}`$ and in this system without escape, $`P(1)=0`$. Equation (33) can be transformed into

$$\underset{j}{}|\mathrm{\Lambda }_t^{(j)}|^\beta \mathrm{exp}\left[P(\beta )t\right]_t\mathrm{}1.$$
(34)
Comparing with Eq. (29), we obtain

$$P(D_\mathrm{H})=D_\mathrm{H}\mathrm{Re}s_𝐤.$$
(35)
If the wavenumber $`𝐤`$ vanishes, the cumulative function (27) becomes $`F_{𝐤=0}(\theta )=\theta /2\pi `$, which forms a straight line in the complex plane. In this equilibrium limit, the relaxation rate vanishes, $`s_{𝐤=0}=0`$, so that $`P(D_\mathrm{H})=0`$. It means for this system without escape that $`D_\mathrm{H}=1`$, as it should be for a straight line. For a non-vanishing but small wavenumber, the Hausdorff dimension is expected to deviate from unity. Inserting $`D_\mathrm{H}=1+\delta `$ and the dispersion relation $`s_𝐤=𝒟k^2+𝒪(k^4)`$ in Eq. (35), we can expand both sides in powers of the wavenumber, using the aforementioned properties of the topological pressure. This straightforward calculation shows that the Hausdorff dimension of the hydrodynamic mode is given by Eq. (28) .

### B Reactive modes

#### 1 General results

We are now going to extend this result to the reactive modes of the isomerization process (3). The modes of interest here are the ones describing the evolution of the difference between the two colors $`A`$ and $`B`$, a quantity which relax to zero when the system reaches equilibrium. In this case, we consider the operator (21) which differs from the diffusive Frobenius-Perron operator (20) by the extra factor $`(12p_0)^{N_t}`$. Accordingly, this factor should multiply the exponential $`\mathrm{exp}(i𝐤𝐫)`$ in the expression (26) in order to obtain the relaxation rate of the reactive modes, which is thus given by

$$s_𝐤=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}(12p_0)^{N_t}\mathrm{exp}[i𝐤(𝐫_t𝐫_0)]$$
(36)
where $`N_t`$ is the number of collisions performed on catalysts during the time $`t`$ . The cumulative function of the reactive modes are similarly given by

$$F_𝐤(\theta )\underset{t\mathrm{}}{lim}\frac{_0^\theta 𝑑\theta ^{}(12p_0)^{N_t(\theta ^{})}\mathrm{exp}\{i𝐤[𝐫_t(\theta ^{})𝐫_0(\theta ^{})]\}}{_0^{2\pi }𝑑\theta ^{}(12p_0)^{N_t(\theta ^{})}\mathrm{exp}\{i𝐤[𝐫_t(\theta ^{})𝐫_0(\theta ^{})]\}}.$$
(37)
We notice that Eqs. (36) and (37) consistently reduce to Eqs. (26) and (27) in the purely diffusive limit $`p_00`$.

As for the diffusive case, Eq. (36) can be transformed into the condition

$$\underset{j}{}|\mathrm{\Lambda }_t^{(j)}|^1(12p_0)^{N_t^{(j)}}\mathrm{exp}(s_𝐤t)\mathrm{exp}\{i𝐤[𝐫_t^{(j)}𝐫_0^{(j)}]\}_t\mathrm{}1.$$
(38)
The curve $`(\mathrm{Re}F_𝐤,\mathrm{Im}F_𝐤)\text{}`$ is approximated at time $`t`$ by a polygon of sides given by

$$\epsilon _j=|\mathrm{\Delta }F^{(j)}|=|\mathrm{\Lambda }_t^{(j)}|^1|12p_0|^{N_t^{(j)}}\mathrm{exp}(\mathrm{Re}s_𝐤t).$$
(39)
In the limit $`t\mathrm{}`$, this polygon converges to a fractal curve, characterized by a Hausdorff dimension given by $`_j\epsilon _j^{D_\mathrm{H}}1`$. Accordingly, the Hausdorff dimension of the reactive mode should satisfy the condition

$$\underset{j}{}|\mathrm{\Lambda }_t^{(j)}|^{D_\mathrm{H}}|12p_0|^{N_t^{(j)}D_\mathrm{H}}\mathrm{exp}(D_\mathrm{H}\mathrm{Re}s_𝐤t)_t\mathrm{}1.$$
(40)
By analogy with Eq. (34), we can define a function $`Q(\alpha ,\beta )`$ by

$$\underset{j}{}|\mathrm{\Lambda }_t^{(j)}|^\beta |12p_0|^{\alpha N_t^{(j)}}\mathrm{exp}\left[Q(\alpha ,\beta )t\right]_t\mathrm{}1$$
(41)
which can be rewritten as

$$Q(\alpha ,\beta )=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}(12p_0)^{\alpha N_t}|\mathrm{\Lambda }_t|^{1\beta }.$$
(42)
Comparing Eq. (41) with Eq. (38), we obtain

$$Q(D_\mathrm{H},D_\mathrm{H})=D_\mathrm{H}\mathrm{Re}s_𝐤.$$
(43)
The relation between this new function and the topological pressure is given by

$$P(\beta )=Q(0,\beta ).$$
(44)
If the wavenumber $`𝐤`$ vanishes, the cumulative function (37) is reduced to

$$F_{𝐤=0}(\theta )=\underset{t\mathrm{}}{lim}\frac{_0^\theta 𝑑\theta ^{}(12p_0)^{N_t(\theta ^{})}}{_0^{2\pi }𝑑\theta ^{}(12p_0)^{N_t(\theta ^{})}}.$$
(45)
This is not a straight line as in the diffusive case: each initial condition is here weighted by a factor $`(12p_0)^{N_t}`$. Provided that $`p_0<0.5`$, this means that the contribution of a given initial condition to $`F_0`$ will be smaller if the corresponding trajectory meets many catalysts, and vice versa. But the Hausdorff dimension of $`F_0`$ is anyway equal to unity. Indeed, for $`𝐤=0`$, the relaxation rate is given by

$$s_{𝐤=0}=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}(12p_0)^{N_t},$$
(46)
which is also equal to $`Q(1,1)`$ according to the definition (42). Using these equalities, we can deduce from (43) that $`D_\mathrm{H}=1`$ for vanishing wavenumber. Again the Hausdorff dimension will deviate from unity for a non-vanishing wavenumber. The reactive dispersion relation is given by $`s_𝐤=2\kappa 𝒟^{(\mathrm{r})}k^2+𝒪(k^4)`$. If we compare it with (46) at vanishing wavenumber, we obtain

$$Q(1,1)=2\kappa =\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}(12p_0)^{N_t}.$$
(47)
Using this equality, the dispersion relation and the expression $`D_\mathrm{H}(k)=1+ak^2+𝒪(k^4)`$, Eq. (43) can be expanded in powers of the wavenumber, giving an expression of the Hausdorff dimension of the reactive mode

$$D_\mathrm{H}(𝐤)=1\frac{𝒟^{(\mathrm{r})}k^2}{2\kappa +\frac{Q}{\alpha }(1,1)+\frac{Q}{\beta }(1,1)}+𝒪(k^4).$$
(48)
An expression for the two partial derivatives of $`Q(\alpha ,\beta )`$ can be directly obtained from Eq. (42), provided that $`p_0<0.5`$:

$$\frac{Q}{\alpha }(1,1)=\mathrm{ln}(12p_0)\underset{t\mathrm{}}{lim}\frac{N_t(12p_0)^{N_t}}{t(12p_0)^{N_t}}$$
(49)
and

$$\frac{Q}{\beta }(1,1)=\underset{t\mathrm{}}{lim}\frac{(12p_0)^{N_t}\mathrm{ln}|\mathrm{\Lambda }_t|}{t(12p_0)^{N_t}}.$$
(50)
Similarly, the reactive diffusion coefficient is given by

$$𝒟^{(\mathrm{r})}=\underset{t\mathrm{}}{lim}\frac{(12p_0)^{N_t}(𝐫_t𝐫_0)^2}{2dt(12p_0)^{N_t}}$$
(51)
where $`d`$ is the dimension of the space of positions $`𝐫`$ where diffusion occurs.

#### 2 Limit of small reaction probability $`p_01`$

Let us now consider the limit case $`p_01`$. The three terms of the denominator appearing in Eq. (48) can be expanded in powers of $`p_0`$. The reaction rate is given by

$`2\kappa `$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}[\mathrm{\hspace{0.33em}2}p_0N_t+2p_0^2N_t+{\displaystyle \frac{8}{3}}p_0^3N_t`$ (54)
$`2p_0^2\left(N_t^2N_t^2\right)4p_0^3\left(N_t^2N_t^2\right)`$
$`+4p_0^3({\displaystyle \frac{1}{3}}N_t^3+{\displaystyle \frac{2}{3}}N_t^3N_tN_t^2)+𝒪(p_0^4)]`$
$`2\kappa `$ $`=`$ $`2\nu p_0+2(\nu \sigma ^2)p_0^2+4({\displaystyle \frac{2}{3}}\nu \sigma ^2+\tau )p_0^3+𝒪(p_0^4)`$ (55)
in terms of the frequency (11) of collisions on the catalysts $`\nu `$ and of the quantities

$`\sigma ^2`$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\left(N_t^2N_t^2\right)`$ (56)
$`\tau `$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\left({\displaystyle \frac{1}{3}}N_t^3+{\displaystyle \frac{2}{3}}N_t^3N_t^2N_t\right).`$ (57)
The two further terms are given by

$`{\displaystyle \frac{Q}{\alpha }}(1,1)`$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}[2p_0N_t+2p_0^2N_t+{\displaystyle \frac{8}{3}}p_0^3N_t`$ (60)
$`4p_0^2\left(N_t^2N_t^2\right)8p_0^3\left(N_t^2N_t^2\right)`$
$`+12p_0^3({\displaystyle \frac{1}{3}}N_t^3+{\displaystyle \frac{2}{3}}N_t^3N_t^2N_t)+𝒪(p_0^4)]`$
$`{\displaystyle \frac{Q}{\alpha }}(1,1)`$ $`=`$ $`2\kappa +2\sigma ^2p_0^2+4(\sigma ^22\tau )p_0^3+𝒪(p_0^4)`$ (61)
and

$`{\displaystyle \frac{Q}{\beta }}(1,1)=\lambda `$ $`+`$ $`2p_0\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}C_{\kappa \lambda }(t)`$ (62)
$`+`$ $`2p_0^2\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\left[\left(1+2N_t\right)C_{\kappa \lambda }(t)D_{\kappa \lambda }(t)\right]`$ (63)
$`+`$ $`4p_0^3\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}[({\displaystyle \frac{2}{3}}+2N_t+2N_t^2N_t^2)C_{\kappa \lambda }(t)`$ (65)
$`(1+N_t)D_{\kappa \lambda }(t)+{\displaystyle \frac{1}{3}}E_{\kappa \lambda }(t)]+𝒪(p_0^4)`$
$`{\displaystyle \frac{Q}{\beta }}(1,1)`$ $`=`$ $`\lambda +\alpha _{\kappa \lambda }p_0+\beta _{\kappa \lambda }p_0^2+\gamma _{\kappa \lambda }p_0^3+𝒪(p_0^4)`$ (66)
where

$`\lambda `$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\mathrm{ln}|\mathrm{\Lambda }_t|`$ (67)
$`C_{\kappa \lambda }(t)`$ $`=`$ $`N_t\mathrm{ln}|\mathrm{\Lambda }_t|N_t\mathrm{ln}|\mathrm{\Lambda }_t|`$ (68)
$`D_{\kappa \lambda }(t)`$ $`=`$ $`N_t^2\mathrm{ln}|\mathrm{\Lambda }_t|N_t^2\mathrm{ln}|\mathrm{\Lambda }_t|`$ (69)
$`E_{\kappa \lambda }(t)`$ $`=`$ $`N_t^3\mathrm{ln}|\mathrm{\Lambda }_t|N_t^3\mathrm{ln}|\mathrm{\Lambda }_t|`$ (70)
$`\alpha _{\kappa \lambda }`$ $`=`$ $`2\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}C_{\kappa \lambda }(t)`$ (71)
$`\beta _{\kappa \lambda }`$ $`=`$ $`2\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\left[\left(1+2N_t\right)C_{\kappa \lambda }(t)D_{\kappa \lambda }(t)\right]`$ (72)
$`\gamma _{\kappa \lambda }`$ $`=`$ $`4\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}[({\displaystyle \frac{2}{3}}+2N_t+2N_t^2N_t^2)C_{\kappa \lambda }(t)`$ (74)
$`(1+N_t)D_{\kappa \lambda }(t)+{\displaystyle \frac{1}{3}}E_{\kappa \lambda }(t)].`$
For these coefficients, the index $`\kappa `$ refers to the reaction, related to the number of collisions on a catalyst $`N_t`$, while the index $`\lambda `$ refers to the dynamical instability, described by the stretching factor $`\mathrm{\Lambda }_t`$ and the Lyapunov exponent $`\lambda `$.

The reactive diffusion coefficient (51) can also be expanded in powers of the reaction probability $`p_0`$ as

$$𝒟^{(r)}=𝒟+\mathrm{\Delta }p_0+𝒪(p_0^2)$$
(75)
in terms of the standard diffusion coefficient

$$𝒟=\underset{t\mathrm{}}{lim}\frac{1}{2dt}(𝐫_t𝐫_0)^2$$
(76)
and of the coefficient

$$\mathrm{\Delta }=\underset{t\mathrm{}}{lim}\frac{1}{dt}\left[N_t(𝐫_t𝐫_0)^2N_t(𝐫_t𝐫_0)^2\right]$$
(77)
which characterizes the correlation between the transport by diffusion and the reaction.

We notice that, in the limit $`p_00`$, the reaction rate (55) vanishes while the reactive diffusion coefficient (51) tends to the diffusion coefficient (76), as expected.

Using (61) and (66), the Hausdorff dimension (48) becomes

$`D_\mathrm{H}(𝐤)`$ $`=`$ $`1+{\displaystyle \frac{𝒟^{(\mathrm{r})}k^2}{\lambda \left[1\frac{\alpha _{\kappa \lambda }}{\lambda }p_0\frac{\beta _{\kappa \lambda }}{\lambda }p_0^2\frac{\gamma _{\kappa \lambda }}{\lambda }p_0^32\frac{\sigma ^2}{\lambda }p_0^24\frac{\sigma ^2}{\lambda }p_0^3+8\frac{\tau }{\lambda }p_0^3+𝒪(p_0^4)\right]}}+𝒪(k^4)`$ (78)
$`=`$ $`1+{\displaystyle \frac{𝒟^{(\mathrm{r})}k^2}{\lambda }}[1+{\displaystyle \frac{\alpha _{\kappa \lambda }}{\lambda }}p_0+2{\displaystyle \frac{\sigma ^2}{\lambda }}p_0^2+{\displaystyle \frac{\beta _{\kappa \lambda }}{\lambda }}p_0^2+{\displaystyle \frac{\alpha _{\kappa \lambda }^2}{\lambda ^2}}p_0^2+4{\displaystyle \frac{\sigma ^2}{\lambda }}p_0^38{\displaystyle \frac{\tau }{\lambda }}p_0^3`$ (80)
$`+{\displaystyle \frac{\gamma _{\kappa \lambda }}{\lambda }}p_0^3+2{\displaystyle \frac{\alpha _{\kappa \lambda }\beta _{\kappa \lambda }}{\lambda ^2}}p_0^3+4{\displaystyle \frac{\alpha _{\kappa \lambda }\sigma ^2}{\lambda ^2}}p_0^3+{\displaystyle \frac{\alpha _{\kappa \lambda }^3}{\lambda ^3}}p_0^3+𝒪(p_0^4)]+𝒪(k^4).`$
or

$$D_\mathrm{H}(𝐤)=1+\left[\frac{𝒟^{(\mathrm{r})}}{\lambda }+𝒪(p_0)\right]k^2+𝒪(k^4).$$
(81)
We immediately see that for $`p_00`$, Eq. (80) reduces to the expression (28) obtained for the diffusive case.

We are now going to test this formula for two examples: a reactive triadic multibaker and a reactive two-dimensional periodic Lorentz gas.

## III The reactive triadic multibaker map

### A Description of the model

The multibaker map is a simple model used to mimic a diffusive motion . In the case of the triadic symmetric multibaker, we consider an infinite one-dimensional chain of squares,

$$𝒮_l=\{(x,y,l):\mathrm{\hspace{0.33em}0}x1,\mathrm{\hspace{0.33em}0}y1\},$$
(82)
where $`l\text{}`$ is the label of the square. A point particle performs jumps from square to square, according to the transition rule $`l(l1)`$ if $`x1/3`$, $`ll`$ if $`1/3<x2/3`$ and $`l(l+1)`$ if $`x>2/3`$. In order to introduce a reaction, the point particle carries a color $`A`$ or $`B`$, and the central part $`]\frac{1}{3},\frac{2}{3}]`$ of one square over $`L`$ is catalytic: if the point particle is in this portion of the square, it changes its color with a probability $`p_0`$ (see Fig. 1). The map of the model is thus given by

$$𝝋(x,y,l,c)=\{\begin{array}{cc}(3x,\frac{y}{3},l1,c),\hfill & 0x1/3\hfill \\ (3x1,\frac{y+1}{3},l,\epsilon (l)c),\hfill & 1/3<x2/3\hfill \\ (3x2,\frac{y+2}{3},l+1,c),\hfill & 2/3<x1\hfill \end{array}$$
(83)
where $`c`$ is a discrete variable taking the values $`+1`$ or $`1`$ when the color of the point particle is respectively $`A`$ or $`B`$, and

$$\epsilon (l)=\{\begin{array}{cc}1\hfill & \text{if }l\text{ is a catalytic square, with a probability }p_0\text{ ,}\hfill \\ +1\hfill & \text{otherwise .}\hfill \end{array}$$
(84)
We consider the probability density $`q(x,y,l,c)`$ which evolves under the Frobenius-Perron equation

$$q_{t+1}(x,y,l,c)=\{\begin{array}{cc}q_t(\frac{x}{3},3y,l+1,c),\hfill & 0y1/3\hfill \\ [1p_0(l)]q_t(\frac{x+1}{3},3y1,l,c)+p_0(l)q_t(\frac{x+1}{3},3y1,l,c),\hfill & 1/3<y2/3\hfill \\ q_t(\frac{x+2}{3},3y2,l1,c),\hfill & 2/3<y1\hfill \end{array}$$
(85)
where

$$p_0(l)=\{\begin{array}{cc}p_0\hfill & \text{if }l\text{ is a catalytic square ,}\hfill \\ 0\hfill & \text{otherwise .}\hfill \end{array}$$
(86)
The Kolmogorov-Sinai entropy per iteration (12) of the reactive triadic multibaker is given by

$$h_{\mathrm{KS}}=\mathrm{ln}3+\frac{p_0}{3L}\mathrm{ln}\frac{1}{p_0}+\frac{1p_0}{3L}\mathrm{ln}\frac{1}{1p_0}$$
(87)
because the frequency (11) is equal to the ratio of the area of a reactive rectangle to the area of $`L`$ squares: $`\nu =1/3L`$. We notice that the result (87) can also be obtained by using a generating partition which establishes the equivalence between the reactive multibaker map (83) and a Markov chain (see Ref. ).

If we define the difference between the densities of the colors $`A`$ and $`B`$,

$$g_t(x,y,l)=q_t(x,y,l,+1)q_t(x,y,l,1),$$
(88)
we obtain the following equation for its evolution

$$g_{t+1}(x,y,l)=\{\begin{array}{cc}g_t(\frac{x}{3},3y,l+1),\hfill & 0y1/3\hfill \\ [12p_0(l)]g_t(\frac{x+1}{3},3y1,l),\hfill & 1/3<y2/3\hfill \\ g_t(\frac{x+2}{3},3y2,l1),\hfill & 2/3<y1\hfill \end{array}$$
(89)
which is the form of the reactive operator (21) in the case of the present multibaker model.

### B Relaxation rates and eigenmodes of the reactive operator

We restrict our description to a unit cell of length $`L`$, $`l=0`$ being the catalytic square, and we choose quasiperiodic boundary conditions, assuming that the solution of the equation (89) is quasiperiodic on the chain with a wavenumber $`k`$. Moreover, we suppose exponentially decaying solutions, with decaying factor $`\chi (k)=\mathrm{exp}s_k`$ where $`|\chi (k)|1`$ or $`\mathrm{Re}s_k0`$

$$g_t(x,y,l)\chi ^t\mathrm{exp}(ikl).$$
(90)
The following operator $`\widehat{Q}_k`$ governs the time evolution of $`g_t(x,y,l)`$ :

$$\widehat{Q}_k\{\begin{array}{cc}g_{t+1}(x,y,0)\hfill & =\theta (\frac{1}{3}y)g_t(\frac{x}{3},3y,1)+(12p_0)\theta (y\frac{1}{3})\theta (\frac{2}{3}y)g_t(\frac{x+1}{3},3y1,0)\hfill \\ & +e^{ikL}\theta (y\frac{2}{3})g_t(\frac{x+2}{3},3y2,L1)\hfill \\ g_{t+1}(x,y,1)\hfill & =\theta (\frac{1}{3}y)g_t(\frac{x}{3},3y,2)+\theta (y\frac{1}{3})\theta (\frac{2}{3}y)g_t(\frac{x+1}{3},3y1,1)\hfill \\ & +\theta (y\frac{2}{3})g_t(\frac{x+2}{3},3y2,0)\hfill \\ & \mathrm{}\hfill \\ g_{t+1}(x,y,L2)\hfill & =\theta (\frac{1}{3}y)g_t(\frac{x}{3},3y,L1)+\theta (y\frac{1}{3})\theta (\frac{2}{3}y)g_t(\frac{x+1}{3},3y1,L2)\hfill \\ & +\theta (y\frac{2}{3})g_t(\frac{x+2}{3},3y2,L3)\hfill \\ g_{t+1}(x,y,L1)\hfill & =e^{ikL}\theta (\frac{1}{3}y)g_t(\frac{x}{3},3y,0)+\theta (y\frac{1}{3})\theta (\frac{2}{3}y)g_t(\frac{x+1}{3},3y1,L1)\hfill \\ & +\theta (y\frac{2}{3})g_t(\frac{x+2}{3},3y2,L2)\text{.}\hfill \end{array}$$
(91)
We now want to obtain the eigenvalues and the eigenstates of this operator

$$\widehat{Q}_k\{\mathrm{\Psi }(x,y,l)\}_{l=0}^{l=L1}=\chi \{\mathrm{\Psi }(x,y,l)\}_{l=0}^{l=L1}.$$
(92)
We suppose that the leading eigenstates are uniform along the unstable direction $`x`$

$$\mathrm{\Psi }(x,y,l)=𝒴(y,l).$$
(93)
If we now cumulate over $`x`$ and $`y`$, we get

$$G_t(x,y,l)=_0^x𝑑x^{}_0^y𝑑y^{}\mathrm{\Psi }(x^{},y^{},l)xC(y,l)$$
(94)
where

$$C(y,l)=_0^y𝑑y^{}𝒴(y^{},l)$$
(95)
is the reactive cumulative function (37) for the multibaker. We replace $`g_t(x,y,l)`$ by $`𝒴(y,l)`$ and $`g_{t+1}(x,y,l)`$ by $`\chi 𝒴(y,l)`$ in Eq. (91) and we integrate over the interval $`[0,y]`$. We obtain

$`C(y,0)`$ $`=`$ $`\{\begin{array}{cc}\frac{1}{3\chi }C(3y,1),\hfill & 0y1/3\hfill \\ \frac{1}{3\chi }[C(1,1)+(12p_0)C(3y1,0)],\hfill & 1/3<y2/3\hfill \\ \frac{1}{3\chi }[C(1,1)+(12p_0)C(1,0)+e^{ikL}C(3y2,L1)],\hfill & 2/3<y1\hfill \end{array}`$ (96)
$`C(y,1)`$ $`=`$ $`\{\begin{array}{cc}\frac{1}{3\chi }C(3y,2),\hfill & 0y1/3\hfill \\ \frac{1}{3\chi }[C(1,2)+C(3y1,1)],\hfill & 1/3<y2/3\hfill \\ \frac{1}{3\chi }[C(1,2)+C(1,1)+C(3y2,0)],\hfill & 2/3<y1\hfill \end{array}`$ (97)
$`\mathrm{}`$ (98)
$`C(y,L2)`$ $`=`$ $`\{\begin{array}{cc}\frac{1}{3\chi }C(3y,L1),\hfill & 0y1/3\hfill \\ \frac{1}{3\chi }[C(1,L1)+C(3y1,L2)],\hfill & 1/3<y2/3\hfill \\ \frac{1}{3\chi }[C(1,L1)+C(1,L2)+C(3y2,L3)],\hfill & 2/3<y1\hfill \end{array}`$ (99)
$`C(y,L1)`$ $`=`$ $`\{\begin{array}{cc}\frac{e^{ikL}}{3\chi }C(3y,0),\hfill & 0y1/3\hfill \\ \frac{1}{3\chi }[e^{ikL}C(1,0)+C(3y1,L1)],\hfill & 1/3<y2/3\hfill \\ \frac{1}{3\chi }[e^{ikL}C(1,0)+C(1,L1)+C(3y2,L2)],\hfill & 2/3<y1\hfill \end{array}.`$ (100)
The eigenvalues are obtained by setting $`y=1`$ in Eq. (100), which leads to the eigenvalue equation

$$\left(\begin{array}{cccccc}(12p_03\chi )& 1& 0& \mathrm{}& 0& e^{ikL}\\ 1& (13\chi )& 1& \mathrm{}& 0& 0\\ 0& 1& (13\chi )& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& (13\chi )& 1\\ e^{ikL}& 0& 0& \mathrm{}& 1& (13\chi )\end{array}\right)\left(\begin{array}{c}C(1,0)\\ C(1,1)\\ C(1,2)\\ \mathrm{}\\ C(1,L2)\\ C(1,L1)\end{array}\right)=0.$$
(101)
We thus have to solve the determinant of an Hermitian matrix to find $`\chi (k)`$. For this purpose, we assume a solution of the form

$$C(1,l)=Ae^{+i\theta l}+Be^{i\theta l}.$$
(102)
After substitution in the system (101), we find that

$$\chi =\frac{1+2\mathrm{cos}\theta }{3}$$
(103)
and

$$\mathrm{cos}kL=\mathrm{cos}\theta L+p_0\frac{\mathrm{sin}\theta L}{\mathrm{sin}\theta }.$$
(104)
For the present reactive system, the relaxation rate $`s(k)=\mathrm{ln}\chi (k)`$ depends on the wavenumber $`k`$ as $`s(k)=2\kappa 𝒟^{(\mathrm{r})}k^2+𝒪(k^4)`$, where $`\kappa `$ is the reaction rate and $`𝒟^{(\mathrm{r})}`$ is the reactive diffusion coefficient.

Solving Eq. (104) for its leading root, we can thus obtain the reaction rate and the reactive diffusion coefficient either in powers of the transition probability $`p_0`$:

$`\kappa `$ $`=`$ $`{\displaystyle \frac{p_0}{3L}}+\left({\displaystyle \frac{1}{L^2}}{\displaystyle \frac{1}{3}}\right){\displaystyle \frac{p_0^2}{6}}+𝒪(p_0^3),`$ (105)
$`𝒟^{(\mathrm{r})}`$ $`=`$ $`{\displaystyle \frac{1}{3}}+{\displaystyle \frac{2p_0}{9L}}+𝒪(p_0^2),`$ (106)
or in the asymptotic limit $`L\mathrm{}`$:

$`\kappa `$ $`=`$ $`{\displaystyle \frac{\pi ^2}{6L^2}}{\displaystyle \frac{2\pi ^2}{3p_0L^3}}+𝒪\left({\displaystyle \frac{1}{L^4}}\right),`$ (107)
$`𝒟^{(\mathrm{r})}`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{3p_0L}}+𝒪\left({\displaystyle \frac{1}{L^2}}\right).`$ (108)
In the limit $`p_00`$, Eq. (105) shows that the reaction rate vanishes linearly with the transition probability and Eq. (106) that the reactive diffusion coefficient tends to the purely diffusive one $`𝒟=1/3`$, as expected. On the other hand, Eq. (107) shows that the reaction rate decreases as $`L^2`$ when the catalytic sides are diluted in the limit $`L\mathrm{}`$, as expected for a diffusion-controled reaction in a one-dimensional system . Moreover, the vanishing of the reactive diffusion coefficient (108) as $`L^1`$ for $`L\mathrm{}`$ is consistent with the results previously obtained in Ref. . We notice the remarkable property that, in the diffusion-controled regime of the reaction, the leading term of the asymptotic expression for the reaction rate (107) does not depend on the transition probability $`p_0`$, which only appears in the next terms. Therefore, the leading term of (107) no longer expresses the property that the reaction rate vanishes in the limit $`p_00`$, as shown by Eq. (105). This is consistent with the fact that the reaction rate is largely independent of the local reaction rate at the catalysts in the diffusion-controled regime, as already noticed in Refs. .

Let us consider more explicitly the simple case $`L=1`$. In this case, Eq. (101) reduces to

$$3\chi 2\mathrm{cos}k1+2p_0=0$$
(109)
which gives the solution

$$\chi (k)=\frac{12p_0+2\mathrm{cos}k}{3}.$$
(110)
For $`k1`$, $`\chi (k)`$ is then given by

$`\chi (k)`$ $``$ $`e^{2\kappa }e^{𝒟^{(\mathrm{r})}k^2}`$ (111)
$``$ $`\chi _0(1𝒟^{(\mathrm{r})}k^2).`$ (112)
Expanding (110) in powers of $`k`$ up to the second degree, we obtain

$$\chi \frac{32p_0}{3}\left(1\frac{1}{32p_0}k^2\right)$$
(113)
which allows us to identify $`\chi _0=\frac{32p_0}{3}=e^{2\kappa }`$ and $`𝒟^{(\mathrm{r})}=\frac{1}{32p_0}`$.

The dependence of the reactive diffusion coefficient $`𝒟^{(\mathrm{r})}`$ on $`p_0`$ is shown in Fig. 2 for different values of the chain length $`L`$. The value of the diffusion coefficient $`𝒟`$ is indicated for comparison.

### C Fractal dimension of the cumulative functions of the eigenmodes

We now want to study the cumulative functions $`C(y,l)`$, solutions of Eq. (100), and more specifically to calculate their dimension. These functions are represented in the complex plane in the case $`L=1`$ in Figs. 3 and 4. Fig. 3 shows their dependence on the reaction probability $`p_0`$, the wavenumber being fixed, $`k=0.1`$. Their dependence on the wavenumber for $`p_0=0.3`$ is shown in Fig. 4. As explained in Section II, these complex cumulative functions can be approximated by a polygon at a given time $`t`$, which means at a certain iteration step $`n`$ in this discrete case. This concept will here be applied using a symbolic dynamics.

The symbolic dynamics is constructed by partitioning the phase space of the multibaker maps depicted in Fig. 1. Each of the $`L`$ squares of the chain is divided in three parts for $`x[0,\frac{1}{3}]`$, $`x]\frac{1}{3},\frac{2}{3}]`$ or $`x]\frac{2}{3},1]`$. These $`3L`$ rectangles are successively assigned the symbols $`\omega _n\{0,1,\mathrm{},3L2,3L1\}`$ from left to right along the chain. At the $`n^{\mathrm{th}}`$ step of the construction described in Subsection II B, the cumulative function is a polygon in the complex plane with its sides corresponding to a subrectangle associated with a sequence of $`n`$ symbols $`\{\omega _{n1}\omega _{n2}\mathrm{}\omega _1\omega _0\}`$. Moreover, we notice that the sequence of symbols also determines the square $`l`$ where the point is located along the chain. At the next steps of the construction, the cumulative function keeps its values already obtained at the borders of each of the aforementioned subrectangles. We can thus denote the values of the cumulative function at the lower and upper borders $`y_\pm (\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)`$ of one of the subrectangles in the square $`l`$ more shortly by

$$C_\pm (\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)=C[y_\pm (\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0),l].$$
(114)
The complex vector defining the side of the polygon at the $`n^{\mathrm{th}}`$ step of the construction is thus given by

$$\mathrm{\Delta }C(\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)=C_+(\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)C_{}(\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0).$$
(115)
The succession of two symbols $`\omega _1\omega _0`$ corresponds to a transition in the multibaker phase-space, with which a certain probability is associated . Using (100), we can construct a $`3L\times 3L`$ transition matrix $`𝖬`$ of elements $`M_{\omega _1\omega _0}`$ describing the factor connecting the portion of square corresponding to $`\omega _0`$ to the one corresponding to $`\omega _1`$. Using this matrix, we can describe the time evolution of the intervals $`\mathrm{\Delta }C(\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)`$ as

$$\mathrm{\Delta }C(\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)=M_{\omega _{n1}\omega _{n2}}M_{\omega _{n2}\omega _{n3}}\mathrm{}M_{\omega _2\omega _1}M_{\omega _1\omega _0}\mathrm{\Delta }C(\omega _0).$$
(116)
We define

$`\mathrm{\Gamma }_n^D(C)`$ $`=`$ $`{\displaystyle \underset{\omega _0\mathrm{}\omega _{n1}}{}}\left|\mathrm{\Delta }C(\omega _{n1},\omega _{n2},\mathrm{},\omega _1,\omega _0)\right|^D`$ (117)
$`=`$ $`{\displaystyle \underset{\omega _0\mathrm{}\omega _{n1}}{}}\left|M_{\omega _{n1}\omega _{n2}}\right|^D\mathrm{}\left|M_{\omega _1\omega _0}\right|^D\left|\mathrm{\Delta }C(\omega _0)\right|^D`$ (118)
$`=`$ $`{\displaystyle \underset{\omega _0\omega _{n1}}{}}(^{n1})_{\omega _{n1}\omega _0}\left|\mathrm{\Delta }C(\omega _0)\right|^D`$ (119)
where the matrix $``$ is composed of the elements:

$$_{\omega \omega ^{}}=\left|M_{\omega \omega ^{}}\right|^D.$$
(120)
Supposing that this matrix admit the spectral decomposition, $`=_\mu |\mu \mu \mu |`$, Eq. (119) can be rewritten as

$$\mathrm{\Gamma }_n^D(C)=\underset{\mu }{}𝐮|\mu \mu ^{n1}\mu |𝐯$$
(121)
where $`𝐮`$ denotes the vector of unit elements $`u_\omega =1`$ and $`𝐯`$ the vector with the elements $`v_\omega =|\mathrm{\Delta }C(\omega )|^D`$. According to the definition of the Hausdorff dimension, for $`n\mathrm{}`$, the quantity (121) must converge to a finite non-vanishing value if $`D=D_\mathrm{H}`$. This will only be the case if $`\mu =1`$. $`D_\mathrm{H}`$ is thus solution of the equation obtained by imposing $`1`$ as eigenvalue of the matrix $``$. Let us consider explicitly the case $`L=1`$. We have to calculate the $`3\times 3`$ determinant

$$\left|\begin{array}{ccc}\left|\frac{1}{3\chi }\right|^{D_\mathrm{H}}1& \left|\frac{1}{3\chi }\right|^{D_\mathrm{H}}& \left|\frac{1}{3\chi }\right|^{D_\mathrm{H}}\\ \left|\frac{12p_0}{3\chi }\right|^{D_\mathrm{H}}& \left|\frac{12p_0}{3\chi }\right|^{D_\mathrm{H}}1& \left|\frac{12p_0}{3\chi }\right|^{D_\mathrm{H}}\\ \left|\frac{1}{3\chi }\right|^{D_\mathrm{H}}& \left|\frac{1}{3\chi }\right|^{D_\mathrm{H}}& \left|\frac{1}{3\chi }\right|^{D_\mathrm{H}}1\end{array}\right|=0.$$
(122)
This gives the equation

$$\left(\frac{12p_0}{3\chi }\right)^{D_\mathrm{H}}+2\left(\frac{1}{3\chi }\right)^{D_\mathrm{H}}1=0.$$
(123)
We will solve (123) in the limit case $`p_00`$, using the expression of $`\chi `$ in the case $`L=1`$ (110). For $`p_0=0`$, we have $`D_\mathrm{H}^{(0)}=\frac{\mathrm{ln}3}{\mathrm{ln}(1+2\mathrm{cos}k)}`$. On the other hand, the Hausdorff dimension takes the unit value for $`k=0`$. Substituting the expansion $`D_\mathrm{H}=1+ak^2+𝒪(k^4)`$ in Eq. (123) and expanding in powers of the wavenumber $`k`$, we can solve for the unknown coefficient $`a`$ and obtain

$$D_\mathrm{H}=1+\frac{k^2}{(32p_0)\mathrm{ln}(32p_0)(12p_0)\mathrm{ln}(12p_0)}+𝒪(k^4).$$
(124)
Now expanding in powers of the transition probability $`p_0`$, we finally obtain

$$D_\mathrm{H}=1+\frac{1}{\mathrm{ln}3}\frac{1}{3}\left(1+\frac{2}{3}p_0+\frac{4}{9}p_0^2+\frac{8}{27}p_0^3\right)k^2+\frac{1}{(\mathrm{ln}3)^2}\left(\frac{4}{27}p_0^2+\frac{80}{243}p_0^3\right)k^2+𝒪(p_0^4k^2)+𝒪(k^4)$$
(125)
We can again compare this expression with the general result (80) of Section II. For the triadic multibaker, the Lyapunov exponent is given by $`\lambda =\mathrm{ln}3`$. Moreover, the stretching factor $`\mathrm{ln}|\mathrm{\Lambda }_t|`$ is constant in the triadic multibaker, which implies that the coefficients $`C_{\kappa \lambda }`$, $`D_{\kappa \lambda }`$, $`E_{\kappa \lambda }`$ and therefore $`\alpha _{\kappa \lambda }`$, $`\beta _{\kappa \lambda }`$, $`\gamma _{\kappa \lambda }`$ defined in Eqs. (68)-(74) are vanishing. Whereupon, Eq. (80) is here reduced to

$$D_\mathrm{H}(k)=1+\frac{𝒟^{(\mathrm{r})}k^2}{\lambda }\left[1+2\frac{\sigma ^2}{\lambda }p_0^2+4\frac{\sigma ^2}{\lambda }p_0^38\frac{\tau }{\lambda }p_0^3+𝒪(p_0^4)\right]+𝒪(k^4).$$
(126)
Using the expansion of the reactive diffusion coefficient up to the third order in $`p_0`$

$$𝒟^{(\mathrm{r})}=\frac{1}{3}\left[1+\frac{1}{3}p_0+\frac{4}{9}p_0^2+\frac{8}{27}p_0^3+𝒪(p_0^4)\right],$$
(127)
Eq. (126) becomes

$`D_\mathrm{H}(k)`$ $`=`$ $`1+{\displaystyle \frac{1}{\mathrm{ln}3}}{\displaystyle \frac{1}{3}}\left(1+{\displaystyle \frac{1}{3}}p_0+{\displaystyle \frac{4}{9}}p_0^2+{\displaystyle \frac{8}{27}}p_0^3\right)k^2+{\displaystyle \frac{1}{(\mathrm{ln}3)^2}}\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{2}{9}}p_0\right)\mathrm{\hspace{0.33em}2}\sigma ^2p_0^2k^2`$ (129)
$`+{\displaystyle \frac{1}{(\mathrm{ln}3)^2}}{\displaystyle \frac{1}{3}}\mathrm{\hspace{0.33em}4}\sigma ^2p_0^3k^2{\displaystyle \frac{1}{(\mathrm{ln}3)^2}}{\displaystyle \frac{1}{3}}\mathrm{\hspace{0.33em}8}\tau p_0^3k^2+𝒪(p_0^4k^2)+𝒪(k^4).`$
The second term of the right members of (125) and (129) are equal. They correspond to $`\frac{𝒟^{(\mathrm{r})}k^2}{\lambda }`$. Let us now compare the other corrections in $`p_0`$. The only correction in $`p_0^2`$ in (129) is $`\frac{1}{(\mathrm{ln}3)^2}\frac{1}{3}\mathrm{\hspace{0.33em}2}\sigma ^2p_0^2k^2`$. We can evaluate $`\sigma ^2=lim_t\mathrm{}\frac{1}{t}(N_t^2N_t^2)`$ in the case of the triadic reactive multibaker $`L=1`$. Indeed here, colliding on a catalyst corresponds to being in the region $`x]\frac{1}{3},\frac{2}{3}]`$. The mean number of collisions on a catalyst per unit time is then

$$\frac{1}{t}N_t=\frac{1}{3}0+\frac{1}{3}1+\frac{1}{3}0=\frac{1}{3}$$
(130)
and

$$\frac{1}{t}N_t^2=\frac{1}{3}0^2+\frac{1}{3}1^2+\frac{1}{3}0^2=\frac{1}{3}$$
(131)
which finally gives us

$$\sigma ^2=\frac{1}{3}\left(\frac{1}{3}\right)^2=\frac{2}{9}.$$
(132)
The coefficient of the correction in $`p_0^2`$ is then equal to $`\frac{1}{(\mathrm{ln}3)^2}\frac{4}{27}`$, which is exactly the coefficient obtained in (125). The correction in $`p_0^3`$ in (129) contains three different contributions: $`\frac{1}{(\mathrm{ln}3)^2}\frac{4}{9}\sigma ^2p_0^3k^2+\frac{1}{(\mathrm{ln}3)^2}\frac{4}{3}\sigma ^2p_0^3k^2\frac{1}{(\mathrm{ln}3)^2}\frac{8}{3}\tau p_0^3k^2`$. Again we can evaluate $`\tau `$

$`\tau `$ $`=`$ $`\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}\left({\displaystyle \frac{1}{3}}N_t^3+{\displaystyle \frac{2}{3}}N_t^3N_t^2N_t\right)`$ (133)
$`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{3}}+{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{27}}{\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{3}}={\displaystyle \frac{2}{81}}.`$ (134)
We obtain $`\frac{1}{(\mathrm{ln}3)^2}\frac{80}{243}`$, which is the coefficient expected for the correction in $`p_0^3`$. Accordingly, we have confirmed that Eq. (80) for the Hausdorff dimension of the reactive modes is valid for the reactive multibaker map. The validity of the approach is furthermore tested in Fig. 5, where the Hausdorff diemansion $`D_\mathrm{H}`$ is plotted versus the square of the wavenumber $`k^2`$. The symbols are obtained numerically by solving Eq. (123) and compared with the expression (125). The agreement is good up to $`p_0=0.3`$. We notice that, beyond this value, further terms would be required in the Taylor expansion (125).

## IV The 2D reactive periodic Lorentz gas

We are now going to consider as second example a two-dimensional reactive periodic Lorentz gas. This model has been introduced in Ref. . In this system, a point particle, carrying a color $`A`$ or $`B`$ undergoes elastic collisions on hard disks fixed in the plane and forming a regular triangular lattice characterized by the lattice fundamental vectors $`𝐞_1=d(1,0)`$ and $`𝐞_2=d(\frac{1}{2},\frac{\sqrt{3}}{2})`$. The radius of the disks is assumed to be equal to unity and $`d`$ is the distance between their centers. We shall work in the finite horizon regime, $`2<d<\frac{4}{\sqrt{3}}`$, for which the diffusion coefficient is finite. Some of the disks play the role of catalysts: when the point particle collides on one of them, it changes its color instantaneously with a probability $`p_0`$. The catalysts form a regular triangular superlattice over the disks lattice, with as fundamental vectors $`𝐄_1=nd(0,\sqrt{3})`$ and $`𝐄_2=nd(\frac{3}{2},\frac{\sqrt{3}}{2})`$, where $`n`$ is an integer parameter controlling the density of catalysts: in the direction of $`𝐄_1`$ and $`𝐄_2`$, one disk over $`n`$ is a catalyst. Globally, one disk over $`N=3n^2`$ is a catalyst. The configurations with $`n=1`$ and $`n=2`$ are depicted in Fig. 6.

Our goal here is to study the cumulative reactive eigenfunctions of this system. As done in Ref. , we will construct the Frobenius-Perron operator governing the time evolution of the color density describing the reactive process in this system. We will then perform a spectral analysis of this operator.

### A The reactive operator

The variables needed to describe the point particle in this system are its position, its velocity and its color, i.e. $`(x,y,v_x,v_y,c)`$, where $`c`$ is a discrete variable taking the values $`+1`$ or $`1`$ when the particle color is $`A`$ or $`B`$. However, the collisions on the disks being elastic, the energy of the point particle is conserved. The magnitude of its velocity is a constant of motion, that we suppose equal to one, $`𝐯=1`$. The velocity can therefore be described by only one coordinate, which can be for example the angle between the velocity and the $`x`$-axis. This system is periodic, which means that instead of describing the dynamics in the full space, we can restrict our description to an elementary cell of the superlattice of catalysts, containing $`N`$ disks (see Fig. 6). Moreover, the relevant information concerning the dynamics is contained in the collision dynamics: between two collisions, the point particle performs a free flight. The three-dimensional flow dynamics can therefore be reduced to a Poincaré-Birkhoff map, describing the dynamics from collision to collision. We will use the Birkhoff coordinates $`𝐱=(j,\theta ,\varpi )`$, $`1jN=3n^2`$, $`0\theta <2\pi `$, $`1\varpi 1`$, where $`j`$ is the index of the disk of the elementary cell on which the collision takes place, $`\theta `$ is an angle giving the position of the impact on this disk and $`\varpi `$ is the sine of the angle between the outgoing velocity and the normal at the impact. In these coordinates, the mapping of the reactive Lorentz gas is given by

$$\{\begin{array}{cc}𝐱_{n+1}=𝝋(𝐱_n),\hfill & \\ t_{n+1}=t_n+T(𝐱_n),\hfill & \\ 𝐥_{n+1}=𝐥_n+𝐚(𝐱_n),\hfill & \\ c_{n+1}=\epsilon (𝐱_n)c_n.\hfill & \end{array}$$
(135)
$`t_n`$ is the time of the $`n^{\mathrm{th}}`$ collision. $`T(𝐱)`$ is the first-return time function. The vector $`𝐥_n`$ gives the cell of the superlattice visited at time $`t_n`$ and $`𝐚(𝐱_n)`$ is a vector giving the jump carried out on the superlattice during the free flight from $`𝐱_n`$ to $`𝐱_{n+1}`$. The change of color is controlled by the function:

$$\epsilon (𝐱_n)=\{\begin{array}{cc}1\hfill & \text{if }j_n\text{ is a catalyst, with a probability }p_0\text{ ,}\hfill \\ +1\hfill & \text{otherwise .}\hfill \end{array}$$
(136)
We notice that the dynamical randomness of the reactive Lorentz gas is characterized by the following Kolmogorov-Sinai entropy per unit time

$$h_{\mathrm{KS}}=\lambda +\nu \left[p_0\mathrm{ln}\frac{1}{p_0}+(1p_0)\mathrm{ln}\frac{1}{1p_0}\right]$$
(137)
where $`\lambda `$ is the unique positive Lyapunov exponent of the two-dimensional Lorentz gas and $`\nu `$ is the frequency of encounters with a catalytic disk.

We now want to obtain the Frobenius-Perron operator associated with this map. We first start from a description of the particle in the full space. The current state and position of the particle are given by a suspended flow in terms of the coordinates $`(𝐱,\tau ,𝐥,c)`$, where $`0\tau <T(𝐱)`$ is the time elapsed since the last collision, $`𝐥`$ is a vector of the superlattice and $`c=\pm 1`$. The phase-space probability density depending on these variables $`q(𝐱,\tau ,𝐥,c)`$ evolves in time under the Frobenius-Perron operator defined by Eq. (10) . This continuous-time operator can be reduced to an operator describing the dynamics from collision to collision under the mapping (135). The details are presented in Ref. . The operator for the modes of wavenumber $`𝐤`$ is given by

$`\widehat{R}_{𝐤,s}\stackrel{~}{q}_{𝐤,s}(𝐱,c)`$ $`=`$ $`\mathrm{exp}[sT(𝝋^1𝐱)i𝐤𝐚(𝝋^1𝐱)]`$ (139)
$`\times \{[1p(𝝋^1𝐱)]\stackrel{~}{q}_{𝐤,s}(𝝋^1𝐱,c)+p(𝝋^1𝐱)\stackrel{~}{q}_{𝐤,s}(𝝋^1𝐱,c)\}`$
where

$$p(𝐱=\{j,\theta ,\varpi \})=\{\begin{array}{cc}p_0\hfill & \text{if }j\text{ is a catalyst,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(140)
The operator $`\widehat{R}_{k,s}`$ is similar to a Frobenius-Perron operator of the mapping $`𝝋`$ up to important factors which take into account the varying time of flight between the collisions, the spatial modulation by the wavenumber $`𝐤`$, as well as the color change on a part of the domain of definition of the Poincaré-Birkhoff mapping. If we now consider the total phase-space density $`f(𝐱)=\stackrel{~}{q}(𝐱,A)+\stackrel{~}{q}(𝐱,B)`$ and the difference of phase-space densities $`g(𝐱)=\stackrel{~}{q}(𝐱,A)\stackrel{~}{q}(𝐱,B)`$, their time evolutions are decoupled and we can write an operator governing the evolution of $`g(𝐱)`$, describing the reactive modes of wavenumber $`𝐤`$ in this system

$$\widehat{W}_{𝐤,s}g(𝐱)=[12p(𝝋^1𝐱)]\mathrm{exp}[sT(𝝋^1𝐱)i𝐤𝐚(𝝋^1𝐱)]g(𝝋^1𝐱),$$
(141)
which we call the reactive operator.

### B Spectral analysis of the reactive operator

To identify the eigenmode of the reactive operator (141) which controls the slowest relaxation of the color, we consider the eigenvalue problem for $`\widehat{W}_{𝐤,s_𝐤}`$ and its adjoint $`\widehat{W}_{𝐤,s_𝐤}^{}`$ imposing $`1`$ as eigenvalue

$`\widehat{W}_{𝐤,s_𝐤}\psi _𝐤(𝐱)`$ $`=`$ $`\psi _𝐤(𝐱),`$ (142)
$`\widehat{W}_{𝐤,s_𝐤}^{}\stackrel{~}{\psi }_𝐤(𝐱)`$ $`=`$ $`\stackrel{~}{\psi }_𝐤(𝐱),`$ (143)
assuming the normalization condition

$$\stackrel{~}{\psi }_𝐤^{}\psi _𝐤=1.$$
(144)
A formal solution for $`\psi _𝐤`$ (respectively $`\stackrel{~}{\psi }_𝐤`$) can be obtained by applying successively $`\widehat{W}_{𝐤,s}`$ (respectively $`\widehat{W}_{𝐤,s}^{}`$) to the unit function:

$`\psi _𝐤(𝐱)`$ $`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \underset{m=1}{\overset{n}{}}}[12p(𝝋^m𝐱)]\mathrm{exp}[s_𝐤T(𝝋^m𝐱)i𝐤𝐚(𝝋^m𝐱)],`$ (145)
$`\stackrel{~}{\psi }_𝐤(𝐱)`$ $`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \underset{m=0}{\overset{n1}{}}}[12p(𝝋^m𝐱)]\mathrm{exp}[s_𝐤^{}T(𝝋^m𝐱)+i𝐤𝐚(𝝋^m𝐱)].`$ (146)
The normalization condition gives an equation to obtain $`s_𝐤`$:

$$1=\underset{n\mathrm{}}{lim}\underset{m=n+1}{\overset{n}{}}[12p(𝝋^m𝐱)]\mathrm{exp}[s_𝐤T(𝝋^m𝐱)i𝐤𝐚(𝝋^m𝐱)].$$
(147)
We notice that if we consider segments of trajectories of time $`t`$ from initial conditions $`𝐫_0`$ we have that

$`t`$ $``$ $`{\displaystyle \underset{m=n+1}{\overset{n}{}}}T(𝝋^m𝐱),`$ (148)
$`𝐫_t𝐫_0`$ $``$ $`{\displaystyle \underset{m=n+1}{\overset{n}{}}}𝐚(𝝋^m𝐱),`$ (149)
$`(12p_0)^{N_t}`$ $``$ $`{\displaystyle \underset{m=n+1}{\overset{n}{}}}[12p(𝝋^m𝐱)],`$ (150)
so that Eq. (147) can be transformed into Eq. (36) for the reaction rate.

We are interested in the smallest relaxation rate $`s_𝐤`$ which dominates at long times $`t`$. For the reactive case, the dispersion relation is given by

$$s_𝐤=2\kappa 𝒟^{(\mathrm{r})}𝐤^2+𝒪(𝐤^4),$$
(151)
where $`\kappa `$ is the reaction rate and $`𝒟^{(\mathrm{r})}`$ is the reactive diffusion coefficient. Solving Eq. (147) allows us to obtain the values of these two coefficients. Fig. 7 shows the dependence of the reactive diffusion coefficient $`𝒟^{(\mathrm{r})}`$ on the reaction probability $`p_0`$, for two different configurations of catalysts, where $`n=1`$ and $`n=2`$. The value of the diffusion coefficient $`𝒟`$ is indicated for comparison. We notice the remarkable result that the behavior is similar to the one observed for the multibaker model (see Fig 2). The explanation is the following. The case $`n=1`$, with a single catalytic disk in each cell where the motion is ballistic, is similar to the multibaker map for $`L=1`$ (see Fig. 1). In this case, the reactive diffusive coefficient rapidly increases with the transition probability $`p_0`$ in both the reactive multibaker model with $`L=1`$ and the reactive Lorentz gas with $`n=1`$. In contrast, the reaction starts to be controled by the diffusion as soon as $`n2`$ in the reactive Lorentz gas  and, similarly, for the reactive multibaker with $`L2`$. In these cases, the reactive diffusion coefficient has a milder dependence on the transition probability $`p_0`$ and becomes even smaller than the diffusion coefficient for large values of $`p_0`$.

### C Fractal dimension of the cumulative eigenmodes

The cumulative function of the reactive eigenmode is obtained by integration of (145) over the angle $`\theta `$ between the initial position and the $`x`$-axis

$$F_𝐤(\theta ;j,\varpi )=\frac{1}{2\pi }_0^\theta 𝑑\theta ^{}\underset{n\mathrm{}}{lim}\underset{m=1}{\overset{n}{}}[12p(𝝋^m𝐱^{})]\mathrm{exp}[s_𝐤T(𝝋^m𝐱^{})i𝐤𝐚(𝝋^m𝐱^{})],$$
(152)
which is equivalent to the formula (37) according to the correspondences (148)-(150).

Fig. 8 shows these functions in the complex plane and as a function of the wavenumber $`𝐤`$, in the case $`n=1`$, $`d=2.3`$ and $`p_0=0.3`$. The Hausdorff dimension of these curves can be computed by a box-counting algorithm for different values of the wavenumber $`𝐤`$. The values obtained are plotted versus $`k^2`$ in Fig. 9, in the case $`n=1`$, $`d=2.3`$, $`p_0=0.1`$ and $`p_0=0.3`$. We have numerically computed the coefficients $`\alpha _{\kappa \lambda }`$, $`\beta _{\kappa \lambda }`$ and $`\sigma ^2`$, using the expressions (71), (72) and (56). We have obtained $`\alpha _{\kappa \lambda }0.23`$, $`\sigma ^20.18`$ and $`\beta _{\kappa \lambda }0.15`$. The solid lines in Fig. 9 are given by Eq. (80), considering the corrections up to the second order in $`p_0`$. We observe the remarkable agreement between the theoretical prediction and the numerical results.

## V Conclusions

In this paper, we have studied in detail simple reaction-diffusion systems in which an isomerization is induced by chaotic dynamics. We have constructed the reactive modes of relaxation toward the thermodynamic equilibrium. These modes are given by singular distributions having a cumulative function at small values of the wavenumber. The cumulative functions form fractal curves in the complex plane because of the chaotic dynamics and in spite of the randomness introduced by the transition probability $`p_0`$ of the color in the isomerization. We explain this remarkable result by the fact that the randomness of the system is of a mild type since its Kolmogorov-Sinai entropy per unit time (12) remains finite.

We have here derived a formula given by Eq. (48) for the Hausdorff dimension of the reactive relaxation modes in chaotic reaction-diffusion systems with two degrees of freedom (43). We have defined a function $`Q(\alpha ,\beta )`$ generalizing the Ruelle topological pressure $`P(\beta )`$. For small values of the wavenumber, we obtain a relationship between the Hausdorff dimension $`D_\mathrm{H}`$, the reactive diffusion coefficient $`𝒟^{(\mathrm{r})}`$, the reaction rate $`\kappa `$, and two derivatives of $`Q(\alpha ,\beta )`$, namely $`\frac{Q}{\alpha }(1,1)`$ and $`\frac{Q}{\beta }(1,1)`$. In the limit of a small reaction probability $`p_0`$, the Hausdorff dimension is related to $`𝒟^{(\mathrm{r})}`$ and to the Lyapunov exponent $`\lambda `$, with corrections in powers of $`p_0`$. For $`p_0=0`$, we recover the relationship (28) obtained in Ref.  for the diffusive hydrodynamic modes. We have shown that the reaction rate, the reactive diffusion coefficient and the Hausdorff dimension can be expressed in terms of quantities which characterize the statistical correlations between the number $`N_t`$ of encounters with the catalysts, the distance $`𝐫_t𝐫_0`$ travelled in the diffusive motion across the system, and the stretching rates $`\mathrm{ln}|\mathrm{\Lambda }_t|`$ characterizing the dynamical instability and the chaos. Our formula (48) thus establishes a new relationship between macroscopic reaction-diffusion processes and the microscopic dynamical chaos.

We have tested our results on two reaction-diffusion systems with a simple isomerization reaction: a triadic multibaker map and a two-dimensional reactive periodic Lorentz gas. In both cases, we study the limit of small wavenumbers and small reaction probabilities. The agreement with the theory is excellent.

Moreover, the study confirms and clarifies results previously obtained for similar models of reaction-diffusion . In particular, our results confirm that the reaction-diffusion process will obey on a large spatial scale and a long time scale the macroscopic reaction-diffusion equations:

$`{\displaystyle \frac{\rho _A}{t}}`$ $``$ $`𝒟_{AA}{\displaystyle \frac{^2\rho _A}{𝐥^2}}+𝒟_{AB}{\displaystyle \frac{^2\rho _B}{𝐥^2}}\kappa (\rho _A\rho _B),`$ (153)
$`{\displaystyle \frac{\rho _B}{t}}`$ $``$ $`𝒟_{BA}{\displaystyle \frac{^2\rho _A}{𝐥^2}}+𝒟_{BB}{\displaystyle \frac{^2\rho _B}{𝐥^2}}+\kappa (\rho _A\rho _B).`$ (154)
for the color densities $`\rho _A`$ and $`\rho _B`$. According to results (75) and (77), the diffusion coefficients of Eqs. (154) are given by

$`𝒟_{AA}`$ $`=`$ $`𝒟_{BB}={\displaystyle \frac{𝒟+𝒟^{(\mathrm{r})}}{2}}=𝒟+{\displaystyle \frac{\mathrm{\Delta }}{2}}p_0+𝒪(p_0^2)`$ (155)
$`𝒟_{AB}`$ $`=`$ $`𝒟_{BA}={\displaystyle \frac{𝒟𝒟^{(\mathrm{r})}}{2}}={\displaystyle \frac{\mathrm{\Delta }}{2}}p_0+𝒪(p_0^2)`$ (156)
which explains that the coefficients of cross-diffusion can be small if the transition probability $`p_0`$ is small. Eq. (77) shows that this cross-diffusion depends on the correlation between the number $`N_t`$ of catalysts encountered and the square $`(𝐫_t𝐫_0)^2`$ of the distance travelled by the particle in its diffusive motion. Moreover, Figs. 2 and 7 and Eqs. (106) and (108) for the multibaker model show that the coefficient of cross-diffusion can become significant for high transition probability $`p_0`$ and if the density of catalysts is very high, in the case $`L=1`$ for the multibaker and $`n=1`$ for the Lorentz gas. We may here say that the cross-diffusion, i.e. the coupling between the diffusions of colors $`A`$ and $`B`$, is induced by the reaction.

In conclusion, the modes of relaxation toward equilibrium turn out to display unexpected fractal properties which have their origin in the chaotic dynamics of the microscopic motion of particles, even if this motion is time-reversal symmetric and volume-preserving. This result has been previously obtained for diffusive processes and, here, we have shown that this result extends to reaction-diffusion systems where it allows us to understand the relaxation to equilibrium at the level of the microscopic motion.

Acknowledgements. We thank Professor G. Nicolis for support and encouragement in this research. The authors are financially supported by the National Fund for Scientific Research (F. N. R. S. Belgium). This research is supported, in part, by the Interuniversity Attraction Pole program of the Belgian Federal Office of Scientific, Technical and Cultural Affairs, by the Royal Academy of Belgium, and by the F. N. R. S. .