text:
# ULTRA HIGH ENERGY COSMIC RAYS: present status and future prospects

## ULTRA HIGH ENERGY COSMIC RAYS: present status and future prospects

#### Department of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK

### 0.1 Introduction

For the purposes of this review I define ultra high-energy cosmic rays (UHECRs) as those cosmic rays having energies above $`10^{19}`$ eV. There is currently great interest in them, partly because we have little idea as to how Nature creates particles or photons of these energies. Also we know enough about their energy spectrum and arrival direction distribution to believe that we have an additional problem: their sources must be reasonably nearby (within 100 Mpc) but there is no evidence of the anisotropies anticipated if the galactic and inter-galactic magnetic fields are as weak as astronomers tell us.

The distance limit comes from a combination of well-understood particle physics and the universality of the 2.7 K radiation. Interactions of protons and heavier nuclei with this, and other, radiation fields degrade the energy of particles rather rapidly. In the case of protons, the reaction is photopion production, while heavier nuclei are photodisintegrated by the 2.7 K radiation and the diffuse infrared background. These effects were first recognised by Greisen, and by Zatsepin and Kuzmin, and lead to the expectation that the energy spectrum of cosmic rays should terminate rather sharply above $`4\times 10^{19}`$ eV (the GZK cut-off). Above $`4\times 10^{19}`$ eV about 50% of particles must come from within 130 Mpc, while at $`10^{20}`$ eV the corresponding distance is 20 Mpc.

The most recent data suggest that particles do exist with energies beyond the GZK cut-off and that the arrival direction distribution is isotropic. The mass of the cosmic rays above $`10^{19}`$ eV is not known, although there are experimental limits on the fraction of photons that constrain one of the models proposed to resolve the enigma.

### 0.2 Measurement of UHECR

The properties of UHECRs are obtained by studying the cascades, or extensive air showers (EAS), they create in the atmosphere. Many methods of observing these cascades have been explored. Currently two approaches seem to be most effective. In one, the density pattern of particles striking an array of detectors laid out on the ground is used to infer the primary energy. At $`10^{19}`$ eV the footprint of the EAS on the ground is several square kilometres so detectors can be spaced many hundreds of metres apart. Alternatively, on clear moonless nights, the fluorescence light emitted when shower particles excite nitrogen molecules in the atmosphere can be observed by massive photomultiplier cameras. This technique, uniquely, allows the rise and fall of the cascade in the atmosphere to be inferred.

The primary energy of the initiating particle or photon is deduced in different ways. For the detector arrays, Monte Carlo calculations have shown that the particle density at distances from 400 – 1200 m is closely proportional to the primary energy. Such a density can be measured accurately (usually to around 20%) and the primary energy inferred from parameters found by calculation. The estimate of the energy depends on the realism of the representation of features of particle interactions within the Monte Carlo model, at energies well above accelerator energies. The currently favoured model (QGSJET) is based on QCD and is matched to accelerator measurements. Although this model appears to describe a variety of data from TeV energies up to $`10^{20}`$ eV , one cannot be certain of the systematic error in the energy estimates.

For the fluorescence detectors, the primary energy is found by integrating the number of electrons in the cascade curve and assuming that their rate of energy loss is close to that at the minimum of the dE/dx curve for electrons, $`2.2`$ MeV per g cm-2 in the case of air. A small model-dependent correction must be made to account for the energy carried by muons and neutrinos into the ground. Ideally, one wants to compare estimates of the primary energy made in the same shower by the two techniques operating simultaneously, but this has yet to be done at these energies. So far all that has been possible is to compare estimates of the fluxes at nominally the same energy.

### 0.3  The Energy Spectrum, Arrival Direction Distribution and Mass of UHECRs

#### 0.3.1 Energy Spectrum

At the time of the Blois meeting (June 2001) data on the energy spectrum from a number of experiments had seemed in good accord . In particular, the rates of events at $`10^{19}`$ eV reported by different experiments were in agreement at the 10 – 15% level. In addition, preliminary data from the Utah-based HiRes group, reported at the International Cosmic Ray Conference in 1999, showed 7 events above $`10^{20}`$ eV, in good agreement with the number anticipated from the flux reported by the Japanese AGASA array.

The situation has now changed dramatically. At the international meeting in Hamburg (August 2001), the AGASA group reported additional data, quite consistent with their earlier work, and described 17 events above $`10^{20}`$ eV. The HiRes group reported on monocular data obtained with one of their cameras from an exposure slightly greater than that of the AGASA group. Assuming a spectrum similar to that reported earlier by the AGASA group, the HiRes team had expected to see about 20 events above $`10^{20}`$ eV, but observed only 2. This unexpected discrepancy is not yet understood.

The data from Fly’s Eye (the earliest fluorescence experiment), Haverah Park (a ground array that used water-Cherenkov detectors), HiRes and AGASA are shown in figure 1. There are several points to note. The Haverah Park energy estimates have been re-assessed  using the QGSJET model. In the range $`3\times 10^{17}`$ to $`3\times 10^{18}`$ eV there is very good agreement between the Fly’s Eye, Haverah Park and HiRes results. A recent Haverah Park analysis  suggests that protons and iron are in the ratio 35:65 in this energy range. With this mixture the spectrum agreement is even better.

This implies that the QGSJET model provides an adequate description of important features of showers up to $`10^{18}`$ eV. However, the AGASA energies have also been estimated with the QGSJET model under the assumption that the primaries are protons at energies above $`3\times 10^{18}`$ eV, the lowest AGASA energy plotted. There is no evidence as to what mass species is dominant at the highest energies but the methods used would lead to an estimate lower by only about 20% if iron nuclei were assumed. This change would be insufficient to reconcile the AGASA-HiRes differences, particularly with regard to the point at which the spectrum slope flattens above $`10^{18}`$ eV. However a combination of a change in the QGSJET model and iron primaries (for which there is no evidence) might go some way to aligning the different results at the highest energies.

There are also unanswered questions about the HiRes data. The ‘disappearance’ of the events reported as being above $`10^{20}`$ eV in 1999 is attributed to a better understanding of the atmosphere which is now claimed to be clearer than had previously been supposed. The Hamburg results were prepared using an ‘average atmosphere’ so presumably subsequently some events will be assigned larger energies and some smaller ones. Two further issues need resolving. Firstly, an accelerator-based calibration of the fluorescence yield  led to the claim ‘that the fluorescence yield of air between 300 and 400 nm is proportional to the electron dE/dx.’ This claim is not consistent with information tabulated in the paper, where it is shown that the yield from 50 keV electrons is similar to that from 1.4 MeV electrons. Also the dE/dx curve plotted there, normalised to the 1.4 MeV measurements, does not fit the accelerator data for 300, 650 and 1000 MeV electrons. The latter discrepancy is about 15 – 20% and in such a direction as would increase the HiRes energies. Secondly, Nagano et al.  has described a new measurement of the yield in air from 1.4 MeV electrons. In what seems to be a very careful study, they find that the earlier results  gave a higher yield at 356.3 nm and 391.9 nm than is found now. Nagano attributes the absence of background corrections as being responsible for at least some of the discrepancies . The longer wavelengths become increasingly important when showers are observed at the large distances common at the highest energies because of Rayleigh scattering. The magnitude of the adjustments that need to be made to the HiRes data are presently unclear and further fluorescence yield measurements are certainly required.

At the Hamburg meeting, the HiRes group also reported data from their stereo system. With 20% of the monocular exposure, they found 1 event with an energy estimated as being close to $`3\times 10^{20}`$ eV, the energy of the largest event found with the Fly’s Eye detector . My opinion is that the spectra from AGASA and HiRes will come together as further understanding is gained of the models and of the atmosphere. Knowledge of the mass composition will also help considerably. For now it seems that trans-GZK events do exist but that the flux of them is less certain than appeared a few years ago.

#### 0.3.2 Arrival Direction Distribution

The angular resolution of shower arrays and of fluorescence detectors is typically 2 – 3. The arrival direction of the 59 events with energy above $`4\times 10^{19}`$ eV registered by the AGASA group is shown in figure 2 . The distribution is isotropic and there is no preference for events to come from close to the galactic or the super-galactic planes. The AGASA group draw attention to a number of clusters, where a cluster is defined as a grouping of 2 or more events within $`2.5^{}`$. It is claimed that the number of doublets (5) and triplets (1) could have arisen by chance, with probabilities of 0.1% and 1%. The implications of such clusters would be profound but the case for them is not proven. The angular bin was not defined a priori and the data set used to make the initial claim for clusters is also being used in the ‘hypothesis testing’ phase. Furthermore, I note that the directions of the 7 most energetic events observed by Fly’s Eye, Haverah Park, Yakutsk and Volcano Ranch do not line up with any of the 6 cluster directions.

It is hard to understand the isotropy if the local extragalactic magnetic field is really just $`10^9`$ gauss. A proton of $`10^{20}`$ eV would be deflected by only about $`2^{}`$ over a distance of 20 Mpc if the field has a 1 Mpc correlation length . If the fields were much higher, as has been suggested , then the lack of anisotropy might be understood, but more energy is then stored in the magnetic field and this might create other difficulties. Similarly, if the charge of the particles initiating the showers was much higher than Z=1, the isotropy could be explained.

#### 0.3.3 Mass Composition

Interpretation of the data on UHECRs is hampered by our lack of knowledge of the mass of the incoming particles. Data from several experiments can be interpreted as indicating a change from a dominantly iron beam near $`3\times 10^{17}`$ eV to a dominantly proton beam at $`10^{19}`$ eV. But the situation is unclear and quite open at higher energies. The data are just too limited and the interpretations are ambiguous as both the fluorescence detectors and ground arrays rely on shower models to deduce composition information.

It is unlikely that the majority of the events claimed to be near $`10^{20}`$ eV have photons as parents as some of the showers have the normal numbers of muons (the tracers of primaries that are nuclei) and the profile of the most energetic fluorescence event is inconsistent with that of a photon primary . Furthermore, there is now evidence that less than 40% of the events at $`10^{19}`$ eV are photon-initiated. This limit has been set in two ways. Taking the energy spectrum as measured by Fly s Eye as being independent of the mass of the incoming particles, the rate of showers coming at large angles to the vertical can be calculated. Using Haverah Park data, it has been found that the observed rate of inclined showers is much higher than would be expected if the primary particles were mainly photons . A more traditional attack on the problem by the AGASA group, searching for showers which have significantly fewer muons than normal, has given the same upper limit . It is unlikely that many events are created by neutrinos as the distribution of zenith angles would be different from that observed. Indeed, in all aspects so far measured, events of $`10^{20}`$ eV look like events of $`10^{19}`$ eV, but ten times larger, and this can be reiterated as we go to lower and lower energies were nuclei seem certain to be the progenitors of showers.

### 0.4 Theoretical Interpretations

The UHECR enigma is attracting significant theoretical attention. Some ideas suppose a form of electromagnetic acceleration while others invoke new physics.

Currently it is popularly believed that cosmic rays with energies up to about $`10^{15}`$ eV are energised by a process known as ‘diffusive shock acceleration’. Supernovae explosions are identified as the likely sites, although so far there is no direct evidence for acceleration of nuclei by supernova remnants at any energy. The diffusive shock process, which has its roots in some early ideas of Fermi, has been extensively studied since its conception in the late 1970s. In  it is shown that the maximum energy attainable is given by E = kZeBR$`\beta `$c, where B is the magnetic field in the region of the shock, R is the size of the shock region and k is a constant less than 1. The same result has been obtained by a number of people, e.g. , and most authors agree upon it. However, some claim that the diffusive shock acceleration process can be modified to give much higher energies than indicated by the equation and that radio galaxy lobes, in particular, are probable acceleration sites. It is difficult to see how an energy of $`3\times 10^{20}`$ eV can be accounted for if the size of the shock region is 10 kpc and the magnetic field is 10 $`\mu `$G (values thought typical of lobes of radio galaxies), as even the optimum estimate of the energy is lower by a factor of 3 than the observational upper limit. It could be that the magnetic fields are stronger than is usually supposed, a line of argument that also comes from the arrival direction work mentioned above.

Proposals have been made which dispense with the need for electromagnetic acceleration. Attention has usually been focused on the highest energy events ($`>10^{20}`$ eV). However, it is my view that proposers of some of the more exotic mechanisms often overlook one or more important points. Any mechanism able to explain the highest energy events must also explain those above about $`3\times 10^{18}`$ eV, where the galactic component possibly disappears. The spectrum above this point is possibly too smooth to imagine that there are two or more radically different components — although this might be seen as an almost philosophical argument, particularly in the light of figure 1! In addition, the solutions proposed must produce particles at the top of the atmosphere that can generate showers of the type we see and now understand rather well. Finally, source energetics cannot be ignored: there seems little point in inventing a mechanism to ‘solve’ the GZK cut-off problem that requires a source region that is unrealistically energetic.

An overview of the various non-electromagnetic processes proposed can be found in  and I will only discuss one of these here. It has been suggested that UHECR arise from the decay of super-heavy relic particles. In this picture, the cold dark matter is supposed to contain a small admixture of long-lived super-heavy particles with a mass $`>10^{12}`$ GeV and a lifetime greater than the age of the Universe . It is argued that such particles can be produced during reheating following inflation or through the decay of hybrid topological defects such as monopoles connected by strings. I find it hard to judge how realistic these ideas are but the decay cascade from a particular candidate  has been studied in some detail  and . A feature of the decay cascade is that an accompanying flux of photons and neutrinos is predicted which may be detectable with a large enough installation. In particular photons are expected to be between 2 and 10 times as numerous as protons above $`10^{19}`$ eV. The anisotropy question has been examined and specific predictions have been made for the anisotropy that would be seen by a Southern Hemisphere observatory. Observation of the predicted anisotropy, plus the identification of appropriate numbers of neutrinos and photons, would be suggestive of a super-heavy relic origin. However, the experimental results on the photon/proton ratio at $`10^{19}`$ eV described above clearly do not support it.

### 0.5 Detectors of the Future

The Pierre Auger Observatory was conceived to measure the properties of the highest energy cosmic rays with unprecedented statistical precision. When completed, it will consist of two instruments, in the Northern and Southern Hemispheres, each covering 3000 km2. The design calls for a hybrid system with 1600 particle detectors and three or four fluorescence detectors at each of the sites. The particle detectors will be deep water-Cherenkov tanks arranged on a 1.5 km hexagonal grid. These detectors have been selected because water acts as a very effective absorber of the multitude of low energy electrons and photons found at distances of about 1 km from the shower axis.

At the Southern site fluorescence detectors will be set up at four locations, one near the centre of the particle array with the others on small promontories at the array edge: the site is close to the town of Malargüe in Mendoza Province, Argentina. During clear moonless nights, signals will be recorded in both the fluorescence detectors and the particle detectors, while for roughly 90% of the time only particle detector data will be available. Construction of an engineering array containing 40 water tanks and a section of a fluorescence detector has been completed (September 2001) and all of the sub-systems of the Observatory have been demonstrated. The first ‘hybrid’ events were recorded in December 2001 and there is great confidence that the Observatory will work as designed. When the Auger Observatory at Malargüe has operated for 10 years, we would expect to have recorded over 300 events above $`10^{20}`$ eV.

Achieving an exposure greater than that targeted by the Auger Observatory is a formidable challenge. A promising line is the development of an idea due to Linsley. The concept is to observe fluorescence light produced by showers from space with satellite-borne equipment. It is proposed to monitor $`10^5`$ km2 sr (after allowing for an estimated 8% on-time). Preliminary design studies have been carried out in Italy and the USA. An Italian-led collaboration has proposed a design that is under study for flight in the International Space Station. This is known as EUSO (the Extreme Universe Space Observatory), and has the potential to detect neutrinos in large numbers as well as UHECRs. Observations are scheduled to start in 2008. This type of project requires considerable technological development but may be the only way to push to energies beyond whatever energy limits are found with the Auger instruments.
Acknowledgements: I am grateful to the organisers for inviting me to the Blois meeting. On-going support of PPARC to work on ultra high-energy cosmic rays at the University of Leeds is gratefully acknowledged. I also thank my many colleagues in the Pierre Auger project for helping to make a 10-year-old dream become a reality.


format: text

text:
# A Search for Intra-group Planetary Nebulae in the M81 Group

## References

Durrell, P. R., Ciardullo, R., Feldmeier, J. J., Jacoby, G. H., & Sigurdsson, S. 2001, ApJ, submitted

Feldmeier, J. J. 2000, Ph.D. Thesis, Penn State University

Feldmeier, J. J., Ciardullo, R., Krelove, K., Jacoby, G., & Durrell, P. R. 2002, ApJ, in preparation

Garcia, A. M. 1993, A&AS, 100, 47

Jacoby, G. H., Ciardullo, R., Ford, H. C., & Booth, J. 1989, ApJ, 344, 70

Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996, Nature, 379, 613

Richstone D. O., & Malumuth, E.M. 1983, ApJ, 268, 30

Theuns, T., & Warren, S. J. 1997, MNRAS, 284, L11

Yun, M. S., Ho, P. T. P., & Lo, K. Y. 1994, Nature, 372, 530


text:
# 1 Introduction

## 1 Introduction

It is commonly known that complex cubic Ginzburg-Landau (CGL) equations constitute a class of universal models for the description of pattern formation in various nonlinear dissipative media . Equations of the CGL type are also frequently used in nonlinear fiber optics, to describe the formation, stability, and interactions of solitary pulses (SPs). CGL equations with constant coefficients apply to a long nonlinear optic-fiber link if the pulses in it are broad enough, so that the corresponding dispersion length is essentially larger than the amplification spacing. In this case, the periodically placed amplifiers and filters (which limit the gain to a relatively narrow spectral band) may be considered in the uniformly-distributed approximation, neglecting their discreteness .

The single-component CGL equation with the cubic nonlinearity possesses a well-known exact solitary-pulse solution , which includes an internal chirp (phase curvature across the pulse). However, a fundamental drawback of this solution is that it is *unstable*, as the zero solution to the cubic CGL equation, i.e., a background on top of which the pulse is built, is unstable due to the presence of linear gain in the equation. Development of physically realistic models in which solitary pulses are fully stable is a problem of an obvious interest in its own right, and it also has profound importance for fiber-optic communications (see the book  and, for instance, a recent work ), as well as for the design of soliton-generating fiber-loop lasers . In the context of optical telecommunications, an issue of fundamental significance is the development of adequate models for multi-component systems, corresponding to a wavelength-division-multiplexed (WDM) multi-channel scheme implemented in the optical fiber. An objective is then to design a system supporting pulses in all the channels, which must be stable against small perturbations and mutual collisions (see, e.g., books ,  and a recent paper ).

A single-channel system, which may suppress the instability of the zero solution, simultaneously allowing for the existence of stationary pulses and thus opening way for them to be stable, was proposed and studied by means of analytical perturbative methods in Ref. , and then tested by direct simulations in Ref. . In this system, the CGL equation is linearly coupled to an additional dissipative equation, which is a linear one in the most fundamental and physically relevant version of the model . In the context of optical fibers, the system may be realized as a *dual-core* fiber, in which an active core carries the linear gain, filtering, temporal dispersion, and Kerr nonlinearity, while the parallel-coupled core is lossy, its intrinsic nonlinearity, dispersion, and filtering being negligible. It has recently been shown  that this model may describe transmission of fully stable optical solitary pulses with an internal *chirp* (intrinsic phase structure) in an indefinitely long fiber-optic communication link.

In this paper, we study the stability and collisions of chirped solitary pulses in a system of nonlinearly coupled CGL equations stabilized by means of the aforementioned scheme, i.e., by linearly coupling each CGL equation to its own linear dissipative counterpart. In its simplest versions that are considered in this work, the model includes two or three nonlinearly coupled stabilized subsystems, which is a prototype of a *stabilized WDM system*. The most interesting issues, on which we focus in this work, are the stability of SPs in this system and collisions between them (including a simultaneous collision between three pulses belonging to three channels). Besides the obvious relevance to optical telecommunications, the obtained results are of interest in their own right, demonstrating a new type of stable traveling pulses and collisions between them in a generalized (multicomponent) CGL system.

The paper is structured as follows. In section 2, we give a detailed formulation of the model with two channels, and exact solutions for SPs in each channel. In section 3, we develop an analytical approach to the collision problem, based on the perturbation theory. In particular, a prediction for position shifts of the pulses in the case of a quasi-elastic collision is obtained in a fully analytical form. In section 4, results of systematic direct simulations of the collisions are displayed for both inelastic and quasi-elastic cases; in the latter case, the analytical predictions are found to be in good agreement with the numerical results (in a parametric region where the perturbation theory applies). In the same section, a region in the model’s parametric space is identified, in which the pulses are *fully* stable, i.e., against both small perturbations and mutual collisions. A generalization for a three-channel model is briefly considered in section 5, with a conclusion that the pulses are also stable against simultaneous collision of three of them. The paper is concluded by section 6.

## 2 The two-channel model and exact solutions for the pulses

### 2.1 The model

The simplest version of the model describes the propagation of two waves, $`u`$ and $`v`$, carried by two different wavelengths in the active core of a dual-core optical fiber. The waves interact through the cross-phase modulation induced by the Kerr effect in the active core. The fields $`u`$ and $`v`$ are assumed to be linearly coupled to two other fields, $`\varphi `$ and $`\psi `$ respectively, which propagate in the passive core placed parallel to the active one. In fact, it is not necessary to assume that all the long fiber-optic link has a dual-core structure; instead, it is sufficient to have short segments of the dual-core fiber periodically installed into the link. Then, in the same uniformly-distributed approximation which was mentioned above in relation to amplification and filtering, we may consider an *effectively homogeneous* dual-core fiber link. In fact, this approximation was already applied to the single-channel dual-core model in Ref. .

Thus, the model is based on the following system of normalized equations governing the propagation of the four above-mentioned electromagnetic waves in the two linearly coupled fiber cores:

$`i\left(u_z+cu_t\right)+\left({\displaystyle \frac{1}{2}}Di\right)u_{tt}iu+\left(|u|^2+\sigma |v|^2\right)u`$ $`=`$ $`K\varphi ,`$ (1)
$`i\left(v_zcv_t\right)+\left({\displaystyle \frac{1}{2}}Di\right)v_{tt}iv+\left(|v|^2+\sigma |u|^2\right)v`$ $`=`$ $`K\psi ,`$ (2)
$`i\left(\varphi _z+\delta \varphi _t\right)+i\mathrm{\Gamma }\varphi `$ $`=`$ $`Ku,`$ (3)
$`i\left(\psi _z\delta \psi _t\right)+i\mathrm{\Gamma }\psi `$ $`=`$ $`Kv,`$ (4)
where the complex electric field envelopes $`u`$, $`v`$, $`\varphi `$, and $`\psi `$ are functions of the propagation distance $`z`$ and retarded time $`t`$ which are defined in the usual way . In the active core, which is equipped with the amplification and filtering, the fields $`u`$ and $`v`$ obey Eqs. (1) and (2), that include the nonlinear self-phase modulation, linear gain, and effective filtering (the latter term is formally tantamount to diffusion in the $`t`$-space). Coefficients in front of the terms in Eqs. (1) and (2) accounting for these three basic effects are all normalized to be equal to 1. On the other hand, the dispersion coefficient $`D`$ is explicitly present in the equations, $`D>0`$ and $`D<0`$ corresponding, respectively, to anomalous and normal dispersion. The field envelopes $`u`$ and $`v`$ in the active core are nonlinearly coupled to each other, interacting through the cross-phase modulation (XPM) induced by the Kerr effect, which gives rise to the XPM coupling coefficient $`\sigma =2`$ in Eqs. (1) and (2). The fields $`u`$ and $`v`$ are linearly coupled, through the corresponding coefficient $`K`$, to their counterparts $`\varphi `$ and $`\psi `$ in the linear dissipative core, which is characterized by a loss coefficient $`\mathrm{\Gamma }`$. Finally, the parameters $`c`$ and $`\delta `$ in Eqs. (1)-(4) stand for the inverse-group-velocity differences between the co-propagating waves in the active and passive (lossy) cores.

A numerical value of $`c`$ (assuming that $`\delta =c`$) will play an important role below. To estimate the value that is of practical interest, we note that the difference in the inverse group velocity is simply related to the frequency separation $`\mathrm{\Delta }\omega `$ between the channels, $`c=\beta \mathrm{\Delta }\omega `$, where $`\beta =k^{\prime \prime }`$ is the standard physical dispersion coefficient . The frequency difference can be expressed in terms of the wavelength separation $`\mathrm{\Delta }\lambda `$, $`\mathrm{\Delta }\omega =\left(2\pi nc_0/\lambda ^2\right)\mathrm{\Delta }\lambda `$, where $`\lambda `$ is the wavelength proper, $`c_0`$ is the light velocity in vacuum, and $`n`$ is the refractive index. For applications, the case of interest is the one with the pulse’s temporal width $`10`$ ps  and the wavelength separation between $`1`$ nm and $`0.1`$ nm . Using these values, and typical values of the dimensional parameters for which the present model provides for the best stability of the pulses (see Eqs. (36) and (37) below), an estimate similar to that presented, e.g., in Ref.  shows that relevant dimensionless values of $`c`$ belong to an interval

$$c20200.$$
(5)
Numerical results will be presented for this region of the values of $`c`$ (see Figs. 6 and 7 below).

### 2.2 The linear spectrum

Before proceeding to the analysis of the full system, it is relevant to consider its linear spectrum in the dissipationless limit, i.e., in the case when the gain, filtering, and loss terms are dropped in Eqs. (1)-(4). In the linear limit, the two subsystems ($`u,\varphi `$) and ($`v,\psi `$) are decoupled, and, looking for the solution to the linearized equations in the ordinary form $`\mathrm{exp}\left(ikzi\omega t\right)`$, one arrives at the following dispersion relations between the propagation distance $`k`$ and frequency $`\omega `$:

$$q=\frac{1}{4D}\left\{\left(c\delta \right)^2D^2\xi ^2\pm \sqrt{\left[\left(c\delta \right)^2D^2\xi ^2\right]^2+16D^2K^2}\right\},$$
(6)
where

$$qk\delta \omega ,\xi \omega \pm \left(c\delta \right)/D,$$
(7)
In the definitions (7) of the shifted propagation constant and frequency, the upper and lower signs pertain, respectively, to the ($`u,\varphi `$) and ($`v,\psi `$) subsystems, while in the dispersion relation (6) the two different signs yield two different branches of the dispersion curve, see Fig. 1.

It follows from Eq. (6) that the spectrum shown in Fig. 1 always has a *gap*,

$$0<\mathrm{\hspace{0.17em}4}Dq<\sqrt{\left(c\delta \right)^4+16D^2K^2}\left(c\delta \right)^2,$$
(8)
inside which, following the general principles , one may expect the existence of a family of gap solitons (note that the full dissipationless version of the present model, including the nonlinear terms in Eqs. (1) and (2), is definitely nonintegrable, therefore “solitons” are meant here simply as solitary waves). However, the objective of this work is not to study that possible family, but to focus on the search for stable pulses in the full model, including the gain, filtering and loss, which is much more relevant to applications.

### 2.3 Solitary-pulse solutions

If the field is launched into one channel only, the system (1)-(4) reduces to a subsystem of Eqs. (1) and (3), or (2) and (4), each including a CGL equation linearly coupled to its linear dissipative counterpart. If, additionally, there is no difference in the group-velocity mismatch in the active and passive core (i.e., $`c=\delta `$), then, in the reference frame moving with the common group velocity, each subsystem (1), (3) or (2), (4), decoupled from the other one, possesses its own pair of exact analytical solutions for chirped SPs, which were actually found in Ref. . In the reference frame moving with the common inverse velocity $`c=\delta `$, the exact solution takes the form

$`u`$ $`=`$ $`A\mathrm{exp}(ikz)\left[\mathrm{sech}(\eta t)\right]^{1+i\mu },`$ (9)
$`\varphi `$ $`=`$ $`AK(i\mathrm{\Gamma }k)^1\mathrm{exp}(ikz)\left[\mathrm{sech}(\eta t)\right]^{1+i\mu },`$ (10)
where the definition of the retarded time $`t`$ is adjusted to the above-mentioned moving reference frame, and the chirp coefficient is

$$\mu =\frac{3}{4}D+\frac{1}{4}\sqrt{32+9D^2}.$$
(11)
The SP’s inverse width $`\eta `$ and squared amplitude $`A^2`$ (which is the peak power, in the application to optics) are expressed in terms of the wavenumber $`k`$,

$`\eta ^2`$ $`=`$ $`{\displaystyle \frac{2k\left(\mathrm{\Gamma }1\right)}{k\left(D\mu 2\right)+\mathrm{\Gamma }\left[\left(1\mu ^2\right)D+4\mu \right]}},`$
$`A^2`$ $`=`$ $`\left(3+{\displaystyle \frac{1}{4}}D^2\right)\mu \eta ^2,`$ (12)
and, finally, the wavenumber itself is determined by a cubic equation,

$$\left(\mu D2\right)(k^21)k+\left(\left(1\mu ^2\right)D+4\mu \right)\left(k^2\mathrm{\Gamma }\left(K^2\mathrm{\Gamma }^2\right)\right)=0.$$
(13)
Clearly, physical solutions of Eq. (13) are those which yield a real wavenumber $`k`$ and $`\eta ^2>0`$. Physical solutions usually exist in pairs, and only the one with a larger value of the peak power may be stable. As it was demonstrated in Refs. , the SP solution with the larger amplitude is indeed stable in a fairly vast region in the parameter space $`(\mathrm{\Gamma },K,D)`$. These pulses definitely remain stable in the framework of the full system (1) - (4). Indeed, because the two above-mentioned subsystems (1), (3) and (2), (4) are coupled solely by the nonlinear XPM terms, the only additional stability condition for a pulse belonging to either subsystem is the linear stability of the zero solution in the mate subsystem, which is always the first condition imposed on parameters of the eligible model.

Note that the exact SP solution displayed above can be extended to the full system of the four equations (1)-(4) if the group-velocity differences vanish, i.e., $`c=\delta =0`$ : the expressions (9)-(13) then yield a solution to the system of the four equations after the transformation $`u,v(u,v)/\sqrt{3}`$, $`\varphi ,\psi (\varphi ,\psi )/\sqrt{3}`$, and setting $`u=v`$ and $`\varphi =\psi `$.

In the most general case, $`c\delta `$, no exact solution for SPs is available, but pulses can be found numerically, see below. In any case, pulses generated by the decoupled subsystems (1), (3) and (2), (4) move at different velocities, hence they may collide. The strong XPM-induced nonlinear coupling between the channels, together with the dissipative character of the system (1)-(4), may give rise to complex dynamical behavior as a result of the collisions.

Our objective in this work is to study in detail collisions between SPs in the system (1)-(4) and their stability. Note, in particular, that in the case when the inverse-group-velocity differences in the active and passive cores are large and nearly equal, i.e., $`c\delta 1`$, the subsystems (1), (3) and (2), (4) nearly decouple, therefore quasi-elastic collisions are expected in this case, while at smaller values of the group-velocity mismatch collisions may be strongly inelastic. These expectations are corroborated by numerical simulations which are displayed below.

## 3 An analytical approach to collisions between solitary pulses

### 3.1 The perturbation theory

In the cases of practical interest to fiber-optic telecommunications, the model is far from any exactly integrable limit, therefore only direct numerical simulations of collisions between pulses (and of their stability), results of which will be summarized in the next section, are really relevant. Nevertheless, some qualitative insight into the collision problem can be gained from an analytical approach, assuming that pulses may be approximated as quasi-solitons. Within the framework of such an approach in its most general possible form, each pulse, in its own reference frame (in which an exact solution is given by Eqs. (9) - (13), assuming that $`\delta =c`$), is taken as

$`u`$ $`=`$ $`Af(\eta (tT))\mathrm{exp}(ikzi\omega t),`$ (14)
$`\varphi `$ $`=`$ $`\mathrm{\Phi }g(\eta (tT))\mathrm{exp}(ikzi\omega t),`$ (15)
$`v`$ $`=`$ $`Af(\eta (t+T))\mathrm{exp}(ikz+i\omega t),`$ (16)
$`\psi `$ $`=`$ $`\mathrm{\Phi }g(\eta (t+T))\mathrm{exp}(ikz+i\omega t).`$ (17)
Here $`f(\eta t)`$ and $`g(\eta T)`$ are (generally speaking, complex) functions accounting for a particular shape of the unperturbed pulses, $`\eta `$ being their inverse temporal width, $`A`$ and $`\mathrm{\Phi }`$ are amplitudes of their two components, and $`\pm \omega `$ and $`\pm T`$ are shifts of the pulses’ central frequencies and temporal positions due to the interaction between them (we consider the interaction between identical pulses, hence the symmetry between the expressions (14), (15) and (16), (17)). Each component of the pulse has its own effective mass, for instance,

$$M_u=A^2\eta ^1_{\mathrm{}}^+\mathrm{}\left|f(x)\right|^2𝑑x,M_\varphi =\mathrm{\Phi }^2\eta ^1_{\mathrm{}}^+\mathrm{}\left|g(x)\right|^2𝑑x$$
(18)
(note that, in the absence of losses and gain, $`M_u+M_\varphi `$, as well as $`M_v+M_\psi `$, are the conserved optical energies in each subsystem). For the pulse given by the solution (9) - (13), one finds

$$M_u=2\eta ^1A^2,M_\varphi =M_uK^2/\left(k^2+\mathrm{\Gamma }^2\right).$$
(19)
The XPM-induced coupling between the two subsystems gives rise to a potential force of attraction between the pulses, which can be calculated by means of well-known methods (see, e.g., Refs. ), provided that XPM may be treated as a small perturbation (conditions for applicability of this assumption will be considered below). It is also known that the filtering term in Eqs. (1) and (2) gives rise to an effective friction force, which, in the most general case, can be evaluated and combined with the potential force by means of the balance equation for the pulse’s momentum (as it was done, for instance, in Refs. ). As a result, one arrives at evolution equations for the soliton’s position and frequency shifts in the following general form,

$`{\displaystyle \frac{dT}{dz}}`$ $`=`$ $`D\omega ,`$ (20)
$`{\displaystyle \frac{d\omega }{dz}}`$ $`=`$ $`\kappa \eta ^2{\displaystyle \frac{M_u}{M_u+M_\varphi }}\omega {\displaystyle \frac{A^4}{\left(M_u+M_\varphi \right)}}U^{}(\eta \left(T+cz\right)),`$ (21)
where the prime stands for the derivative, and the friction coefficient and interaction potential are

$`\kappa `$ $`=`$ $`{\displaystyle \frac{_{\mathrm{}}^+\mathrm{}\left|df(x)/dx\right|^2𝑑x}{_{\mathrm{}}^+\mathrm{}\left|f(x)\right|^2𝑑x}},`$ (22)
$`U(y)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\left|f(xy)\right|^2\left|f(x+y)\right|^2𝑑x`$ (23)
(the XPM coefficient $`\sigma `$ was set equal to its physical value $`2`$). The additional term $`cz`$ in the argument of the potential in Eq. (21) is generated by the group-velocity difference between the two channels, and the ratio of the masses in the friction term on the right-hand side of Eq. (21) appears since the friction force acts only on the $`u`$-component of the pulse, but not on its $`\varphi `$-component. Note that these general equations are also valid in the case of dispersion management, when $`D`$ is not a constant, but a function of $`z`$ .

For the pulses with the shape given by Eqs. (9) - (13), one can find, from the expressions (22) and (23), that

$`\kappa `$ $`=`$ $`\left(4/3\right)\left(1+\mu ^2\right),{\displaystyle \frac{M_u}{M_u+M_\varphi }}={\displaystyle \frac{1}{1+{\displaystyle \frac{K^2}{k^2+\mathrm{\Gamma }^2}}}},`$ (24)
$`U(y)`$ $`=`$ $`4{\displaystyle \frac{2y\mathrm{cosh}(2y)\mathrm{sinh}(2y)}{\mathrm{sinh}^3(2y)}}.`$ (25)
It is worthy to mention that the expression (25) contains no singularity at $`y0`$.

To predict results of the collisions in the general case, the nonlinear non-autonomous ($`z`$-dependent) ODEs (ordinary differential equations) (20) and (21) with the effective potential (25) must be solved numerically. In view of the complexity of this ODE system and its approximate character, it makes sense to focus, instead, on direct simulations of the underlying PDEs (partial differential equations) (1) - (4, which will be done below. Nevertheless, some results can be obtained directly from ODEs (20) and (21). In particular, the most essential effect observed in direct simulations of the underlying PDEs is an inelastic outcome of the collision (merger or complete decay of the pulses), provided that the group-velocity difference $`c`$ is below a certain critical (threshold) value $`c_{\mathrm{cr}}`$. This value may be, very roughly, estimated as that at which the friction and potential forces in Eq. (21) are comparable, that yields

$$c_{\mathrm{cr}}\frac{|D|A^2}{\eta \left(1+\mu ^2\right)}.$$
(26)
To obtain this estimate, it was set that $`M_u/\left(M_u+M_\varphi \right)1`$, which is true in the cases considered below, the expression (24) for the friction coefficient was used, and it was naturally assumed that, for a nontrivial collision, the maximum value of the frequency shift $`\omega `$ is on the order of $`Dc`$, see Eq. (20). Below, it will be seen that the crude estimate (26) helps to understand the fact that $`c_{\mathrm{cr}}`$ is much smaller for the case of normal dispersion than for pulses propagating under anomalous dispersion.

### 3.2 Collision-induced position shifts of the pulses

The ODEs (20) and (21) can be used to obtain *quantitative* results in the limiting case of large $`c`$, so that

$$\left|\frac{dT}{dz}\right|c.$$
(27)
In fact, this is the case when the interaction of the pulses due to XPM may be treated as a small perturbation, and all the above approach is strictly valid. In this case, the term $`\eta T`$ in the argument of the potential in Eq. (21) may be omitted, hence the equation immediately becomes linear. Upon Substitution of $`\omega =D^1dT/dz`$ from Eq. (20) into Eq. (21) and integrating once, it reduces to the following first-order linear equation:

$$\frac{d\left(\eta T\right)}{d\left(c\eta z\right)}+\kappa \frac{\eta }{c}\left(\eta T\right)=\frac{A^2D}{2c^2}U(c\eta z)$$
(28)
(it is more natural to consider, as final dynamical variables, the renormalized temporal shift $`\eta T`$ and propagation distance $`c\eta z`$). To obtain Eq. (28), it was again assumed that $`M_u/\left(M_u+M_\varphi \right)1`$, which will be the confirmed below, and it was substituted $`M_u=2A^2/\eta `$, as per Eq. (19).

Equation (28) can be further simplified if, in addition to the condition (27), the group-velocity difference between the channels is large enough in comparison with an effective friction force, so that

$$\kappa \eta c.$$
(29)
In fact, this condition turns out to be less restrictive than the one (27), see below. Neglecting the friction term, Eq. (28) takes the form

$$\frac{d\left(\eta T\right)}{d\left(c\eta z\right)}=\frac{A^2D}{2c^2}U(c\eta z)$$
(30)
Equations (28) and (30) show that the collision is elastic in the present case, as the inverse-group-velocity shift $`dT/dz`$ is zero at $`z=\pm \mathrm{}`$, i.e., both before and after the collision. Nevertheless, the result of the collision is not trivial. Indeed, Eq. (30) can be used to evaluate an important characteristic of the elastic collision, viz., a residual temporal (position) shift of the pulse,

$$\eta \mathrm{\Delta }T\eta \left[T(z=+\mathrm{})T(z=\mathrm{})\right]=\frac{A^2D}{2c^2}_{\mathrm{}}^+\mathrm{}U(x)𝑑x.$$
(31)
This shift is important as it gives rise to the collision-induced temporal *jitter* of the pulses, see, e.g., Refs.  and . In particular, for the potential (25) one has $`_{\mathrm{}}^+\mathrm{}U(x)𝑑x=2`$, hence

$$\eta \mathrm{\Delta }T=A^2D/c^2.$$
(32)
This analytical prediction will be compared below with results of direct simulations.

To conclude the analytical consideration, we note that the general condition (27) takes a simple form in terms of $`\mathrm{\Delta }T`$. Indeed, in the present case the characteristic collision distance $`\mathrm{\Delta }z`$ is determined by the pulse’s temporal width $`1/\eta `$, so that $`\mathrm{\Delta }z1/\left(\eta c\right)`$, and $`dT/dz`$ may then be estimated as $`\mathrm{\Delta }T/\mathrm{\Delta }z\eta c\mathrm{\Delta }T`$. Inserting this into Eq. (27), one arrives at a simple result:

$$\eta \left|\mathrm{\Delta }T\right|1,$$
(33)
which means that the linear equation (28) applies to the description of collisions between pulses if the resulting normalized temporal shift of the pulse is small.

## 4 Numerical analysis of collisions and stability of pulses

### 4.1 The approach to the problem

We have employed the split-step Fourier algorithm to solve Eqs. (1)-(4) numerically, using, as initial conditions, a superposition of separated waveforms (9) and (10), which yield exact SP solutions for the two decoupled subsystems (1), (3) and (2), (4). Thus, the following initial configurations are used:

$`u(0,t)`$ $`=`$ $`A\left[\mathrm{sech}\left(\eta \left(tT\right)\right)\right]^{1+i\mu },v(0,t)=A\left[\mathrm{sech}\left(\eta \left(t+T\right)\right)\right]^{1+i\mu },`$ (34)
$`\varphi (0,t)`$ $`=`$ $`{\displaystyle \frac{AK}{i\mathrm{\Gamma }k}}\left[\mathrm{sech}\left(\eta \left(tT\right)\right)\right]^{1+i\mu },\psi (0,t)={\displaystyle \frac{AK}{i\mathrm{\Gamma }k}}\left[\mathrm{sech}\left(\eta \left(t+T\right)\right)\right]^{1+i\mu },`$ (35)
which incorporate an initial temporal separation $`2T`$ between the pulses. We will then be able to study both incomplete and complete collisions, corresponding to $`T=0`$ and $`T0`$ respectively (i.e., collisions between initially overlapped and separated SPs, see, e.g., Refs.  and  for the discussion of relative importance of both types of the collisions).

As for the choice of parameters, in most cases we have used the values $`\mathrm{\Gamma }=5`$ and $`K=4`$, which are located almost in the center of the stability domain of the exact SP solution to the decoupled subsystems (1), (3) and (2), (4) . For the dispersion parameter $`D`$ we have chosen the values $`D=\pm 18`$, in the anomalous- and normal-dispersion regimes, respectively. These two values of $`D`$ actually correspond to the carrier wavelength near the zero-dispersion point in a dispersion-shifted fiber  (recall that Eqs. (1) and (2) are normalized so that the effective filtering coefficients in them are set equal to $`1`$). As was shown in Ref. , for a typical physically relevant value of the filtering, the corresponding values $`|D|`$ are indeed close to $`18`$. Also, this value of $`D`$ gives rise to the best stability characteristics for SPs in the single channel model. Using these values for $`\mathrm{\Gamma }`$, $`K`$, and $`D`$, the other parameters of the exact-SP solution can be found from Eqs. (9)-(13): in the case of anomalous dispersion ($`D=+18`$),

$$\mu =0.074,\eta =1.57,k=23.35,A^2=44.56,$$
(36)
and for the normal dispersion ($`D=18`$),

$$\mu =27.074,\eta =0.058,k=23.35,A^2=22.37.$$
(37)
In the anomalous-dispersion regime, the pulses are much narrower, and (quite naturally) have a much smaller chirp than their counterparts existing in the case of the normal dispersion. Note that, for these values of the parameters, $`M_\varphi 0.028M_u`$ according to Eqs. (19), i.e., the mass of the passive-core component of the pulse is negligible in comparison with its active-core component’s mass.

Our first objective is to study in detail all possible outcomes of collisions of stable moving pulses. In the numerical simulations we have found that, depending on the value of the inverse group-velocity mismatch $`c`$ in the active core, three different outcomes of the collisions occur: (a) both SPs perish (“decay”); (b) only one SP survives the collision, while the other one is destroyed (this outcome may also be considered as a merger of two pulses into one); (c) the pulses undergo a quasi-elastic collision, so that both reappear unscathed after the collision.

Formation of a true stable bound state of two solitons as a result of the collision has never been observed in the simulations. However, it will be shown below that, in some special cases (see Fig. 5(a)), a metastable bound state is observed, which exists over a very long propagation distance, but finally collapses into a single pulse.

As far as the above-mentioned outcomes (a) and (b) are concerned, it is important to mention that, in most cases (apart from the exception corresponding to the formation of the metastable bound state, which will be specially considered below), we have found that the corresponding collision distances are quite short, $`z_{\mathrm{coll}}60`$. Thus, in these cases, after passing the short collision distance, there remains, at most, one pulse. Obviously, the only outcome acceptable for applications to the optical telecommunications is a quasi-elastic collision, when both SPs restore their shapes after the interaction.

Detailed results obtained for the incomplete and complete collisions in both normal- and anomalous-dispersion regimes, as well as results for the stability of isolated moving pulses, are summarized below. In all the simulations, the XPM coupling coefficient in Eqs. (1) and (2) was set equal to its physical value, $`\sigma =2`$.

### 4.2 Incomplete Collisions

The results for incomplete collisions ($`T=0`$ in Eqs. (34) and (35)) in the case $`c=\delta `$ are summarized in Table 1, where the three above-mentioned possible outcomes, namely “Decay”, “Merger” and “Elastic”, are indicated. As is shown, for both the normal- and anomalous-dispersion regimes, there exists a critical (i.e., minimum) value $`c_{\mathrm{cr}}`$ of the velocity $`c`$, above which the collision is always elastic. More importantly, the value of $`c_{\mathrm{cr}}`$ in the case of the normal dispersion is much smaller, by a factor $`5`$, than that for the anomalous-dispersion regime. As the smaller critical velocity difference between the channels makes it possible to have a *denser* WDM system, this result shows that the normal-dispersion regime may have an advantage over the more traditional, from the viewpoint of the soliton transmission , anomalous-dispersion regime. On the other hand, an advantage of the latter regime is that, inside a given channel, the pulses forming a data-carrying stream may be packed with a higher density, as their width is much smaller according to Eqs. (36) and (37). In fact, the best approach to the enhancement of the bit-rate of the fiber-optic telecommunication link is to use the channels in *both* the anomalous- and normal-dispersion bands.

The fact that the critical value $`c_{\mathrm{cr}}`$ is much smaller in the normal-dispersion regime can be explained by the crude estimate (26) obtained above on the basis of the analytical consideration. Indeed, the ratio of the values which the expression (26) takes for the parameters (36) and (37) corresponding to the anomalous and normal regimes is very small (however, the estimate is too crude for a detailed quantitative comparison with the numerical results).

Typical examples of the three different outcomes of the incomplete collisions are shown in Figs. 2(a)-(c) (for $`D=18`$) and Figs. 3(a)-2(c) (for $`D=+18`$) in the form of contour plots. As is readily observed, SPs in the normal-dispersion regime are indeed much broader than those in the case of the anomalous dispersion, in accordance with what is predicted by the analytical solution, see Eqs. (36)-(37). In order to test the sensitivity of $`c_{\mathrm{cr}}`$ to a possible mismatch of the interchannel inverse-group-velocity differences $`\delta `$ and $`c`$ between the two cores, we have also simulated incomplete collisions in the case $`\delta c`$, viz., $`\delta =1.1c`$ and $`\delta =2c`$. The results are summarized in Tables 2 and 3, respectively. It is found that, in both normal- and anomalous-dispersion regimes, a mismatch in the values of the inverse-group-velocity differences in the active and passive cores results in an *increase* of the critical value $`c_{\mathrm{cr}}`$, which was again found to be much larger in the case of the anomalous dispersion. As was mentioned above, for the applications it is necessary to make the critical value $`c_{\mathrm{cr}}`$ as small as possible. The results presented here clearly show that the optimum will be attained when $`\delta =c`$, i.e., when the group-velocity differences between the channels are the same in the active and passive cores, which is not difficult to understand in qualitative terms. Indeed, a group-velocity mismatch between the two cores makes it necessary for the main component (in the active core) to “drag” its counterpart in the passive core, which inevitably generates additional losses through the filtering term, thus enhancing inelasticity of collisions between the pulses, cf. the perturbative treatment of the collision in the previous section.

### 4.3 Complete Collisions

#### 4.3.1 Inelastic and elastic collisions

The results of the simulations for complete collisions (i.e., for a case of a sufficiently large initial separation between the colliding pulses) in the case $`\delta =c`$ are summarized in Table 4. In the simulations, the initial temporal separation between the solitons is taken to be equal, approximately, to five pulse widths, i.e., $`\eta T=2.34`$ (see Eqs. (9) and (10)). According to Eqs. (36) and (37), this choice implies $`T=40`$ and $`T=1.5`$, in the regions of normal and anomalous dispersion, respectively.

The results shown in Table 4 suggest that, contrary to the case of incomplete collisions, in the normal-dispersion regime the three outcomes, “decay”, “merger” and “elastic”, alternate with the increase of $`c`$, up to the threshold value $`c_{\mathrm{cr}}=13.92`$, past which only elastic collisions take place. Notice, in particular, a small interval, $`7.85c<8.09`$, where the pulses undergo elastic collisions, which is found *between* the regions where at least one pulse disappears after the collisions. Typical examples of the “merger”, “decay” and “elastic” outcomes for $`D=18`$ are shown in Figs. 4(a) through 4(c). Notice that the collision distance in all the cases does not exceed $`z5`$.

On the other hand, in the anomalous dispersion regime, outcomes of complete collisions resemble what was observed in the case of the incomplete collisions: there are only three, relatively broad intervals of $`c`$, where the outcomes are “merger”, “decay” and “elastic”, with a value of $`c_{\mathrm{cr}}`$ slightly smaller than that in the normal-dispersion regime. However, in the case of $`D=+18`$ it is important to note that, in the “merger” interval ($`0c<3.33`$), the pulses demonstrate a behavior which is remarkably different from that observed in the case of the “decay” or “merger” outcomes of incomplete or complete collisions in all the other cases. In this interval, the pulses stick together and propagate in such a quasi-bound state over a long distance (up to $`z3000`$), after which one of the pulses eventually decays. This behavior is demonstrated in Fig. 5(a) for $`c=2`$, where the colliding SPs form the locked configuration immediately (see the inset in Fig. 5(a)), and then they propagate, keeping this configuration, up to $`z3000`$, where the merger eventually takes place through the destruction of one of the pulses. Examples of other outcomes of the collision, namely, decay of both pulses and their elastic collision, are shown in Figs. 5(b) and 5(c), respectively, with the collision distance being very small, $`z0.1`$.

#### 4.3.2 Position shifts of the pulses in the case of elastic collisions

A significant post-collision effect, in the case of the complete elastic collisions between the SPs, is a temporal shift $`\mathrm{\Delta }T`$. The shift is apparent, for instance, in Figs. 4(c) and 5(c). We have performed simulations to obtain the normalized temporal shift $`\eta \mathrm{\Delta }T`$ as a function of the inverse-group-velocity difference $`c=\delta `$, for different initial separations between SPs. The results are shown in Fig. 6 (for $`D=18`$) and Fig. 7 (for $`D=+18`$); note that these figures display the region of the values of $`c`$ which is relevant to the applications, according to the estimate (5). Each curve starts from the maximum value of the temporal shift corresponding to the critical value $`c=c_{\mathrm{cr}}`$. As is seen, the smallest value of the temporal shift is attained in the case of the normal dispersion. Also, it is observed that a larger initial separation $`\eta T`$ leads to smaller values of the temporal shift in the case of the normal dispersion, while the opposite holds in the anomalous-dispersion case.

It is quite relevant to compare the numerical results displayed in Figs. 6 and 7 with the analytical prediction (32), which should be relevant under the conditions (33) and (29). First of all, the substitution of the expression (32) into the condition (33) leads to an inequality $`c^2A^2|D|`$. The values of the parameters being taken as per Eqs. (36) and (37), we conclude that this inequality holds in the parts of Figs. 6 and 7 where irregular oscillations go over into a systematic decay of $`\eta \mathrm{\Delta }T`$ with the increase of $`c`$.

Next, using the expression (24) for the effective friction coefficient which appears in the second condition (29), and the same values of the parameters from Eqs. (36) and (37), one sees that this condition readily holds in the case of the anomalous dispersion, shown in Fig. 7, and in the case of normal dispersion (Fig. 6) it holds for large values of $`c`$, where the above-mentioned systematic decay of $`\eta \mathrm{\Delta }T`$ takes place.

Thus, one may expect that the analytical result (32) may be correct for sufficiently large values of $`c`$ in both cases. The inspection of Figs. 6 and 7 corroborates this expectation: despite a considerable scatter of the values of the normalized shift, depending on the initial value of the temporal separation between the pulses, the analytically predicted dependence (32) not only qualitatively agrees with the numerical results for large $`c`$, but also, as one can readily check, numerical values of the shift, as predicted analytically and found from the simulations, are fairly close.

### 4.4 The full-stability region for the solitary pulses

The results presented above were restricted to the fixed values of the loss and coupling coefficients, $`\mathrm{\Gamma }=5`$ and $`K=4`$. It is also important to investigate pulse collisions in the full parameter space ($`\mathrm{\Gamma },K,c`$) , and, in particular, to identify parametric domains where solely elastic complete collisions occur. These domains actually represent regions of the *full stability* of SPs in the system of Eqs. (1)-(4), since collisions are natural finite perturbations in this model, against which the pulses must be stable, as well as against infinitesimal perturbations (in this extended definition of the stability, we take into regard only complete collisions, which are inherent perturbations, and disregard incomplete collisions, that strongly depend upon particular initial conditions).

The thus defined full-stability domains, found by means of systematic simulations, are displayed in Figs. 8 and 9 for the cases of the normal and anomalous dispersion, respectively, as gray regions in the ($`\mathrm{\Gamma },K`$) parametric plane for different values of the group-velocity parameter within $`0c10`$. Notice that the domains have, roughly speaking, a boomerang-like shape, resembling the corresponding domain found earlier in the single-channel stabilized CGL model, described by the decoupled subsystems (1), (3) or (2), (4) . More importantly, in both cases $`D=\pm 18`$, there exists a minimum value of $`c`$ necessary for the collisions to be elastic, which is significantly lower in the case of the anomalous dispersion, namely $`c_{\mathrm{min}}=2`$ for $`D=+18`$, and $`c_{\mathrm{min}}=6`$ for $`D=18`$. It should also be noted that the unshaded triangular region shown in Fig. 9 corresponds to the case where the two pulses co-propagate undistorted without actual interaction over a very long distance ($`z4000`$).

## 5 The three-channel model

The model considered in this work can be extended to explicitly include a larger number of WDM channels, i.e., a larger number of the CGL equations, each being coupled to its linear dissipative counterpart. The simplest generalization contains three channels, which are described by the following system of the coupled CGL equations:

$`i\left(u_z+2cu_t\right)+\left({\displaystyle \frac{1}{2}}Di\right)u_{tt}iu+\left(|u|^2+2|v|^2+2|w|^2\right)u`$ $`=`$ $`K\varphi ,`$ (38)
$`iw_z+\left({\displaystyle \frac{1}{2}}Di\right)w_{tt}iw+\left(|w|^2+2|u|^2+2|v|^2\right)w`$ $`=`$ $`K\chi ,`$ (39)
$`i\left(v_z2cv_t\right)+\left({\displaystyle \frac{1}{2}}Di\right)v_{tt}iv+\left(|v|^2+2|u|^2+2|w|^2\right)v`$ $`=`$ $`K\psi ,`$ (40)
$`i\left(\varphi _z+2c\varphi _t\right)+i\mathrm{\Gamma }\varphi `$ $`=`$ $`Ku,`$ (41)
$`i\chi _z+i\mathrm{\Gamma }\chi `$ $`=`$ $`Kw,`$ (42)
$`i\left(\psi _z2c\psi _t\right)+i\mathrm{\Gamma }\psi `$ $`=`$ $`Kv,`$ (43)
where we have assumed that each channel has the same group velocity in the active and passive cores, while the relative velocity between adjacent channels is $`2c`$, as in the above model (1)-(4). Although detailed study of the extended model is beyond the scope of this work, we give here an example of the existence of *fully stable* pulses in the system (38)-(43), which undergo elastic complete collisions with each other. As is shown in Fig. 10 (for $`\mathrm{\Gamma }=5`$, $`K=4`$, $`c=20`$, and $`\eta T=2.34`$), the three pulses after traveling a distance of $`z1`$ collide and then restore their shapes and propagate undistorted. Thus, we may conjecture that the proposed stabilized scheme may be generalized to include a larger number of channels in which case stable pulses may still experience elastic collisions.

## 6 Conclusion

In this paper we have studied in detail collisions of chirped solitary pulses in nonlinearly coupled cubic Ginzburg-Landau (CGL) equations, each being linked to a stabilizing dissipative linear equation. Primarily, the two-channel model was considered. The model may be realized as wavelength-separated data-transmission channels co-existing in a nonlinear dual-core optical fiber, which contains an active core with gain, and a passive core, where the propagation is governed by the linear dissipative equations. Each channel has its components in the active and passive cores, with a linear coupling between them. Nonlinear interaction between different channels is induced by the cross-phase modulation (XPM), which acts in the active core only. Thus, the model describes a WDM multi-channel fiber-optic transmission system, provided that the dispersion length of the pulses is essentially larger than the amplification and filtering spacing, so that the system may be considered in the approximation which assumes a uniformly distributed bandwidth-limited gain in the active core (and a continuous passive core, which, in reality, may consist of short segments periodically inserted into the long fiber-optic link, together with amplifiers and filters).

If the CGL subsystems are decoupled, they possess stable chirped pulse solutions, which can be found in an exact analytical form, provided that the group-velocity parameter is identical in the cubic and linear equations. The XPM-induced nonlinear coupling between the subsystems gives rise to interactions when the pulses belonging to the different subsystems collide. By means of direct simulations, we have studied incomplete and complete collisions in detail. Three different possible outcomes of the collision have been found, in the cases when the dispersion in the active core is anomalous or normal: destruction of both pulses, destruction of one of them, and a quasi-elastic collision. In the latter case, both pulses reappear unscathed after the collision (with some positional shifts), provided that the group-velocity difference between the channels exceeds a critical value. As a result, regions in the model’s parameter space have been identified where the pulses are stable against both small perturbations and mutual (complete) collisions. An analytical perturbation theory was developed to predict the positional shifts. The analytical results agree well with the numerical ones in the cases when applicability conditions for the perturbation theory hold.

For WDM applications, it is important not only to guarantee the quasi-elastic character of the collisions between pulses belonging to different channels, but also to have the critical group-velocity difference between adjacent channels, necessary for the elasticity of the collisions, as small as possible, so that the wavelength separation between the channels may be minimized. To this end, we have found that, in the case of incomplete collisions, the normal-dispersion regime provides an essentially smaller critical velocity, whereas in the case of complete collisions, the critical velocities are almost equal for both signs of the dispersion. However, the region in the parameter space where complete collisions are always elastic is essentially larger in the case of anomalous dispersion, and another advantage of the latter case is that the temporal width of the pulses is much smaller. On the other hand, the residual effect of elastic collisions, viz., the temporal shift of the pulses, which contributes to the soliton jitter in optical communications, is weakest in the case of normal dispersion. Actually, the best solution may be to use the channels in *both* normal-dispersion and anomalous-dispersion bands in the fiber. Finally, we have shown that the model may be extended to include more than two WDM channels, giving rise, in an appropriate region of the corresponding parameter space, to completely stable, three-pulse collisions being quasi-elastic.

Figure Captions

Figure 1. A typical form of the dispersion curves $`q(\chi )`$ (the shifted propagation constant and frequency are defined in Eqs. (7)) for the linearized system in the dissipationless approximation. The curves are shown for $`D=+18`$ (the case of the anomalous dispersion) and $`c=\delta `$. For normal dispersion ($`D=18`$), the dispersion curves are mirror images of those shown in this figure.

Figure 2. Incomplete collisions of two pulses in the normal-dispersion regime ($`D=18`$). (a) Decay of both pulses, with $`c=\delta =6`$. (b) Merger of the pulses, with $`c=\delta =9`$. (c) An elastic collision, with $`c=\delta =14`$.

Figure 3. Incomplete collisions of two pulses in the anomalous-dispersion regime ($`D=+18`$). (a) Decay of both pulses, with $`c=\delta =40`$. (b) Merger of the pulses, with $`c=\delta =2`$. (c) An elastic collision , with $`c=\delta =48`$.

Figure 4. Complete collisions of two pulses in the normal-dispersion regime ($`D=18`$). (a) Merger of the pulses, with $`c=\delta =10`$. (b) Decay of both pulses, with $`c=\delta =13`$. (c) An elastic collision, with $`c=\delta =14`$.

Figure 5. Complete collisions of two pulses in the anomalous-dispersion regime ($`D=+18`$). (a) The merger, with $`c=\delta =2`$. The two pulses get stuck almost immediately (see the inset showing the initial stage of the collision in detail), and then they propagate, keeping this shape up to $`z3000`$, where the merger (in fact, destruction of one of the pulses) eventually takes place. (b) Decay of both pulses, with $`c=\delta =5`$. (c) An elastic collision, with $`c=\delta =48`$.

Figure 6. Relative temporal shift vs. the inverse-group-velocity mismatch $`c=\delta `$ in the case of normal dispersion ($`D=18`$).

Figure 7. Relative temporal shift vs. the inverse-group-velocity mismatch $`c=\delta `$ in the case of anomalous dispersion ($`D=+18`$).

Figure 8. The full-stability region (implying the stability of the solitary pulses against both arbitrary infinitesimal perturbations, and against collisions with a pulse moving in the other channel) in the $`(K,\mathrm{\Gamma })`$ parametric plane, at different fixed values of the inverse-group-velocity difference between the channels, in the case of normal dispersion ($`D=18`$).

Figure 9. The same as in Fig. 8 for the case of anomalous dispersion ($`D=+18`$).

Figure 10. An elastic complete collision between *three* solitary pulses in the three-channel model, in the case of normal dispersion ($`D=18`$), with $`c=\delta =20`$ and $`\eta T=2.34`$.

Table Captions

Table 1. Outcomes of incomplete collisions for $`\delta =c`$.

Table 2. Outcomes of incomplete collisions for $`\delta =1.1c`$.

Table 3. Outcomes of incomplete collisions for $`\delta =2c`$.

Table 4. Outcomes of complete collisions for $`\delta =c`$ and $`\eta T=4.7`$.

| Incomplete Collisions ($`\delta =c`$) | | | |
| --- | --- | --- | --- |
| Normal Dispersion ($`D=18`$) | | Anomalous Dispersion ($`D=+18`$) | |
| Velocity | Outcome | Velocity | Outcome |
| $`c<8.35`$ | Decay | $`c<3.52`$ | Merger |
| $`8.35c<9.53`$ | Merger | $`3.52c<46.98`$ | Decay |
| $`c>9.53`$ | Elastic | $`c46.98`$ | Elastic |

Table 1. Nistazakis, Physical Review E

| Incomplete Collisions ($`\delta =1.1c`$) | | | |
| --- | --- | --- | --- |
| Normal Dispersion ($`D=18`$) | | Anomalous Dispersion ($`D=+18`$) | |
| Velocity | Outcome | Velocity | Outcome |
| $`c<9.88`$ | Decay | $`c<2.42`$ | Merger |
| $`9.88c<10.87`$ | Merger | $`2.42c<59.29`$ | Decay |
| $`10.87c<14.17`$ | Decay | $`c59.29`$ | Elastic |
| $`c14.17`$ | Elastic |  | |

Table 2. Nistazakis, Physical Review E

| Incomplete Collisions ($`\delta =2c`$) | | | |
| --- | --- | --- | --- |
| Normal Dispersion ($`D=18`$) | | Anomalous Dispersion ($`D=+18`$) | |
| Velocity | Outcome | Velocity | Outcome |
| $`c<14.26`$ | Decay | $`c<1.73`$ | Merger |
| $`c14.26`$ | Elastic | $`1.73c<122.17`$ | Decay |
|  | | $`c122.17`$ | Elastic |

Table 3. Nistazakis, Physical Review E

| Complete Collisions ($`\delta =cand\eta T=4.7`$) | | | |
| --- | --- | --- | --- |
| Normal Dispersion ($`D=18`$) | | Anomalous Dispersion ($`D=+18`$) | |
| Velocity | Outcome | Velocity | Outcome |
| $`c<1.88`$ | Decay | $`c<3.33`$ | Merger |
| $`1.88c<2.56`$ | Merger | $`3.33c<12.82`$ | Decay |
| $`2.56c<7.85`$ | Decay | $`c12.82`$ | Elastic |
| $`7.85c<8.09`$ | Elastic |  | |
| $`8.09c<11.71`$ | Merger |  | |
| $`11.71c<13.92`$ | Decay |  | |
| $`c13.92`$ | Elastic |  | |

Table 4. Nistazakis, Physical Review E


text:
# Naked Singularities in Spherically Symmetric, Self-Similar Spacetimes

## Acknowledgements

KSG thanks the University of Natal and the National Research Foundation for ongoing support. He also thanks CIRI for their kind hospitality during the course of this work.


text:
# Parton and dipole approaches in QCD

Here, we discuss QCD predictions on multiplicities in parton and dipole approaches. The most general treatment is based on the notion of the generating functions . The generating function $`G`$ is defined as

$$G(u,y)=\underset{n}{}u^nP_n(y),$$
(1)
where $`P_n`$ is the probability of the $`n`$-particle production at energy denoted by $`y`$, $`u`$ is an auxiliary variable. The mean multiplicity and higher moments of the multiplicity distribution $`P_n`$ are given by the $`u`$-derivatives of $`G`$ at $`u=1`$.

The equations for the generating functions for gluon and quark jets have been proven in QCD up to the next-to-leading order (NLO) of the modified perturbation theory . Their general structure is symbolically represented as

$$G^{}\alpha _SK[GGG]𝑑\mathrm{\Omega }.$$
(2)
It shows that the evolution of $`G`$ indicated by its derivative $`G^{}`$ over the evolution parameter (the transverse momentum or the virtuality) is determined by the cascade process of the production of two partons by a highly virtual time-like parton (the term $`GG`$) which provides new partons in the phase space volume $`d\mathrm{\Omega }`$ and by the escape of a single parton ($`G`$) from a given phase space region.

Therefore this equation contains terms corresponding to inflow and outflow of partons. In fact, it can be interpreted as the kinetic equation with the collision integral in the right hand side. The weight factors are determined by the coupling strength $`\alpha _S`$ and the splitting function $`K`$ which is defined by the interaction Lagrangian. The integral runs over all internal variables, and the symbol $``$ shows that the two created partons share the momentum of their parent. The initial condition for equation (2) is defined by the requirement for the jet to be created by a single initial parton, i.e., by

$$P_n=\delta _{n1};G_0=u.$$
(3)
It is clear from this formula that we have to deal with the non-linear integro-differential probabilistic equation with shifted arguments in the $`GG`$ term under the integral sign.

For the traditional parton approach , the energy conservation at cascade vertices is properly accounted by the shift of corresponding arguments of the generating functions in the integral whereas the transverse momentum ($`p_t`$) limitations are considered approximately as a combined effect of energy conservation and angular ordering.

In the dipole QCD evolution advocated by the Lund group , the triangle phase space in the energy - transverse momentum plane is considered. Limitations on the available phase space (both on energy and $`p_t`$) due to recoil effects, emission of additional gluons, high-$`p_t`$ processes, color reconnection etc can be explicitly implemented .

Both approaches agree up to NLO of the perturbative QCD (pQCD) . However, the corresponding equations for the generating functions differ at higher orders .

Some other modifications of these equations have been earlier proposed . Thus, the problem of formulation of general equations for generating functions can not be still considered as solved. On the way to its solution, we tried to understand how strong is the difference between the parton and dipole approaches in this respect. The partial answer is presented in this talk.

First, let us describe the present situation up to NLO where both approaches coincide. The analytical solutions of the equations successfully predict the energy behavior of mean multiplicities . However, the theoretical values of the ratio $`r`$ of multiplicities between gluon and quark jets are larger than experimental data (at Z0, by about 50$`\%`$ for the leading order (LO) which gives asymptotic values and by 30$`\%`$ for NLO). The description of higher moments is also not perfect . However, in each case NLO corrections improve the agreement compared with LO results. Moreover, pQCD has predicted  in NLO the new unexpected feature of the behavior of cumulant moments which become negative at higher ranks. Their ratio to factorial moments ($`H_q`$) as a function of $`q`$ acquires the minimum at $`q5`$. It has been confirmed by experiment. Many other features of multiparticle production have also been explained in NLO of pQCD .

To go beyond NLO, one of the strategies is to consider the parton evolution equation as a kinetic equation and find its perturbative solution. The systematic method of the Taylor series expansion  leads to the perturbative series for the ratio $`r`$ and the anomalous dimension $`\gamma `$ which determines the energy behavior of mean multiplicity. In such a way, the results up to 3NLO have been obtained . High order corrections almost do not influence conclusions about the energy dependence of mean multiplicity. However, they improve the agreement with experiment on the ratio $`r`$ (up to 15$`\%`$). High order terms completely determine the energy dependence of $`r`$ (its slope)  because the main (NLO) dependence is the same for gluon and quark jets and it cancels in their ratio. Moreover, the oscillations of $`H_q`$-moments with $`q`$ have been predicted both for running  and fixed  coupling regimes at high orders and found in experiment (first in  and then in ). However, the direct perturbative calculations of moments showed that the high order results become unreliable. Formally, this is due to the fact that the expansion parameter $`\gamma _0\alpha _S^{1/2}`$, where $`\alpha _S`$ is the QCD coupling strength, becomes multiplied by the rank $`q`$, and this product is larger than 1. In practice, this means that soft low-$`p_t`$ partons become important. While purely perturbative methods fail here, the exact solution for fixed coupling  and numerical computer solutions  show much better results. In more detail, the parton approach is reviewed in .

To cure these problems, one can try another strategy and modify the parton cascade equations so that NLO results for new and old versions coincide. Such a modification inspired by the dipole model has been proposed . It noticeably improves the agreement in 2NLO order with experimental data on $`r`$ at Z0. The agreement results from different (from parton approach) boundary condition imposed with account of $`p_t`$-conservation. A detailed comparison of dipole and parton formalism  reveals that sensitivity to cascade choice, and hence infrared cutoff quantity representing hadronization details, give uncertainties to 2NLO which reduce the value of higher order refinements on $`r`$.

However, previous experience with parton equations has taught us that higher moments are very sensitive to subtle modifications due to phase space limitations and hadronization effects. Our recent results support this statement. To be able to treat the dipole equations analytically, we use the perturbative expansion. Matching the terms of the same order on both sides of the equations, we get the moments up to (in principle, any) predefined order. The dipole evolution expanded to high orders does not converge towards oscillating $`H_q`$-moments. Instead, higher rank moments diverge severely, and the only hint of oscillations lies in the sign change of divergence for each new order added in the expansion as seen in Fig. 1 for $`H_q`$ at Z0.

The similar sign-changing curves one obtains if the Taylor series of any oscillating function (e.g., sine ) is cut off at some high order terms. Let us note that in previous work on parton equations , which revealed explicit oscillations, some kind of Pade approximation was used where terms up to a definite order were kept both in the numerator (cumulants) and denominator (factorial moments) of $`H_q`$. This is plausible in view of large expansion parameter at hand.

¿From this discussion we conclude that the slope of the ratio $`r`$ and high rank moments of multiplicity distributions are most sensitive to generalization of master evolution equations to higher orders of the perturbative expansion.

For this improvement to be successful, one must solve several problems. First, the calculations apply to parton multiplicities which then are identified with hadron data using the assumption of the local parton-hadron duality (LPHD) . We must also note that pure perturbative QCD predictions can be obtained only for infrared safe quantities. For multiplicities, one has to deal with some cut-off parameter specifying the boundary condition which is influenced by the non-perturbative region of soft partons and hadronization scheme. In general, the cut-off is determined locally in phase space, related to a matching of colors and anti-colors into small singlet sub-systems, often referred to as the ’preconfinement’ assumption. At low scales, the non-perturbative boundary conditions are essential and the energy scale where an asymptotic expansion of moments becomes useful depends on assumptions about hadronization. In the simplest approximation, ”extreme LPHD” if you like, each parton is traded for one hadron at some cut-off scale, usually around 2 to 5 $`\mathrm{\Lambda }_{QCD}`$. How such a boundary condition at a low scale influences the approach to the asymptotic result can be investigated with numerical integration on computer. One such result, $`K_2`$ in gluon dynamics, is presented in Fig. 2. We note that the differences implied by varying the boundary scale from 2 to 4 $`\mathrm{\Lambda }_{QCD}`$ prevails beyond Z0 energies.

One could of course imagine less extreme hadronization schemes, e.g., letting each parton represent a small distribution of hadrons, with some low average. In the figure is also shown the result assuming a poissonian with unit average. The difference is drastic, but it should be admitted that the poissonian has a rather unrealistic tail to large multiplicities. Note however, that factorial moments obtained with a poissonian boundary condition are identical to ’normal’ moments starting with delta-distribution (extreme LPHD). Thus the figure could serve as inspiration to a simultaneous experimental study of both factorial and normal moments, which may give additional insight to multiplicity boundary conditions.

In conclusion, analytical treatment of the iterative cascade has revealed energy conservation and the ”preconfinement” hadronization assumption as most significant for results on $`r`$ at available energies. Whether the analytical approach can be of similar use to understand in more detail higher multiplicity moments remains to be seen. Apparent is that multiplicity distributions are very powerful probes of the perturbative - nonperturbative transition in QCD. We hope that our contribution will give somebody a hint to further work on generalization of QCD master evolution equations.


text:
# Energetics and Vibrational States for Hydrogen on Pt(111)

## Abstract

We present a combination of theoretical calculations and experiments for the low-lying vibrational excitations of H and D atoms adsorbed on the Pt(111) surface. The vibrational band states are calculated based on the full three-dimensional adiabatic potential energy surface obtained from first-principles calculations. For coverages less than three quarters of a monolayer, the observed experimental high-resolution electron peaks at 31 and 68 meV are in excellent agreement with the theoretical transitions between selected bands. Our results convincingly demonstrate the need to go beyond the local harmonic oscillator picture to understand the dynamics of this system.

Hydrogen on metal surfaces, and in particular on Pt(111) has received considerable experimental  and theoretical  attention. Together with Pd, Pt is the most important material for heterogenous catalysis of hydrogenation reactions. On the fundamental side, hydrogen on metal surfaces provides a unique opportunity to observe the crossover from classical to quantum dynamics at relatively elevated temperatures .

For the vibrational excitations of H/Pt(111), Baró et al.  applied high-resolution electron energy-loss spectroscopy (HREELS) and observed for H (D) at coverages up to 0.7 monolayer (ML) two peaks at $`68(50)`$ meV (assigned to the symmetric stretch mode SS) and $`153(112)`$ meV (assigned to the asymmetric stretch mode AS). On the other hand, Richter and Ho  found three peaks at $`67(51)`$ meV, $`112(84)`$ meV and $`153(108)`$ meV in their HREELS measurements for a full ML of H (D). They assigned these to the AS mode, the SS mode, and the unresolved overtone and combination losses, respectively. To date, the disagreement between these observations is not resolved. In addition, recent measurement of the diffusion barrier for this system  yielded a very low value of about $`70`$ meV, considerably less than a previous value of about $`200`$ meV. All these discrepancies emphasize our lack of understanding of this important system.

On the theoretical side, an earlier first-principles (FP) study  at a full ML coverage has shown that H adsorbs in the three-fold hollow fcc site. More recent FP calculations  at lower coverages indicate that unlike for H/Ni(111) and H/Pd(111) , the top site for H/Pt(111) is also a local minimum only slightly higher in energy than the fcc site. In addition, the energy barrier for the fcc-hcp-fcc path is only about $`70`$ meV. In these calculations the main focus has been on adiabatic potential energy surface (APES). The curvatures determined at the potential minima  could not account quantitatively for the observed vibrational excitation energies. Given the recent results of the shallow barrier between the fcc and hcp sites, it is clear that quantum effects are essential in understanding the dynamics of this system, and the vibrational spectroscopy need to be interpreted in terms of excitations between vibrational band states as proposed earlier .

In this Letter we investigate the energetics and vibrational properties of H and D adatoms on the Pt(111) surface. The experiments consist of HREELS investigations of H and D adsorption on Pt(111) and its dependence on coverage. The experiments were performed in an ultra-high vacuum (UHV) apparatus with a base pressure of $`3\times 10^{11}`$ mbar as described elsewhere including sample preparation . The HREEL spectrometer was of the latest commercial Ibach design (Delta 0.5, SPECS, Germany). The H2 exposure was performed by backfilling the preparation chamber while the sample was cooled to $`85`$ K with liquid N2. The exposure was varied between $`0.1`$ and $`500`$ L ($`1\mathrm{L}=1.33\times 10^6`$ mbar$`\times `$s). The coverage $`\theta `$ was determined from thermal desorption (TD) assuming $`\theta =1`$ at saturation. The TD spectra of Christmann et al.  were well reproduced. We focus on the spectrum for the low coverage interval $`\theta 0.75`$ ML, within which the HREEL spectra are qualitatively the same .

The spectra for H and D for $`\theta 0.75`$ for specular geometry are shown in Fig.1. It is well known that the loss intensity for vibrations of atomic H on transition metal surfaces is quite low since it is mainly caused by impact scattering. In spite of these intensity problems (note the magnifications in Fig. 1) the spectra clearly indicate the two H derived losses at $`31`$ and $`68`$ meV. Note that these new observations are quite different from previous results obtained at a full ML . In particular, the previously reported peaks at $`112`$ meV and $`153`$ meV are not observed at these lower coverages. In Fig. 2, we also show the results for off-specular geometry where the $`31`$ meV peak disappears while the $`68`$ meV remains, indicating that the former is mostly of dipole character. For deuterium, there is a clear loss peak at $`23.5`$ meV, corresponding to the $`31`$ meV peak for hydrogen with an isotope downshift. There is a very broad higher mode between $`40`$ and $`60`$ meV not very well resolved from the noise background.

We now turn to the theoretical calculations for the vibrational excitations. The FP calculations for 3D APES were done self-consistently with a parallel plane-wave code  based on density functional theory with the GGA . The Troullier-Martins pseudopotential  was used for Pt and the Vanderbilt nonlocal ultrasoft pseudopotential  for H. Electronic wave functions were expanded in a plane-wave basis with a kinetic energy cutoff at $`35`$ Ry. For the irreducible Brillouin-zone integration eight special $`\stackrel{}{k}`$ points  were used, together with a Fermi smearing of $`0.2`$ eV. The calculated lattice constant for bulk Pt was $`a=7.568`$ a.u. (experimentally $`a=7.415`$ a.u. ). The supercell ($`2\times 2\times 4`$) has 16 Pt atoms and one H atom, corresponding to a quarter monolayer coverage. The $`z`$ axis is oriented towards the vacuum which was $`14`$ Å thick. The coordinates of the Pt atoms were relaxed before adding the H. Only the vertical coordinate of the H atom was relaxed during the calculations for each $`(x,y)`$ position of the APES. We also checked the effects of the relaxation of Pt surface atoms. We found that the potential energy differences between the fcc, bridge and hcp sites changes less than 5 meV .

We first present results for the 3D APES as shown in Fig. 3. The 2D minimum APES is shown in the inset. The top site is metastable and only $`25`$ meV above the lowest energy adsorption fcc site and is surrounded by large energy barriers. Along the fcc-hcp-fcc path, there is a shallow barrier of $`78.9`$ meV. At all sites the H (D) is chemisorbed with an electronic charge transfer from the substrate toward the adatom, similar to that reported for H/Ni(111) . These results agrees well with other recent calculations . In particular, the approach of Källén and Wahnström  is closest to the present work and the resultant APES are very similar. However, there are quantitative differences particularly concerning the curvature and anharmonicity at the fcc sites. These will lead to substantial differences for the vibrational band excitations as detailed below.

To study the vibrational modes of the H and D adatoms, we calculated the band states for the single particle Hamiltonian with the FP 3D APES. The Hamiltonian is diagonalized in a basis set $`\{|\stackrel{}{G},n\}`$ given by a direct product between $`595`$ in-plane plane waves $`\stackrel{}{r}|\stackrel{}{G}`$ and $`14`$ harmonic-oscillator states $`z|\mathrm{\Psi }_n`$ for the $`z`$ coordinate. The lower branches $`n=116`$ for H are shown schematically in Fig. 4, along with the corresponding centers $`\stackrel{}{r}_0`$ of the probability distributions for the $`\stackrel{}{k}=0`$ states, the band centers $`\epsilon _n`$, and the widths $`\mathrm{\Delta }\epsilon _n`$. In addition, many modes have mixed in-plane and vertical character and cannot be described in conventional AS or SS terms. Also, the intervals between the band centers do not match with the local harmonic oscillator model based on the curvatures at the potential minimum. Finally, the higher bands have considerable bandwidth, reflecting the delocalized nature of the adsorbed H through tunnelling processes .

The first band $`n=1`$ localized at the fcc site is very narrow ($`\mathrm{\Delta }\epsilon _1<0.1`$ meV). The next band $`n=2`$ is localized at the hcp site. Its center is at $`\epsilon _2=23.6`$ meV above the center of the lowest band and it has a bandwidth of $`\mathrm{\Delta }\epsilon _2=1.3`$ meV. The difference $`\epsilon _2\epsilon _1`$ confirms that even with the zero point energies taken into account, adsorption at low temperatures is confined to the fcc sites. To analyze the vibrational excitation energies arising from transitions starting from the lowest band $`n=1`$, we will focus only on those excited band states with finite amplitudes at the fcc site. The first excitations are toward the composite band $`34`$, which has a total bandwidth of $`4.1`$ meV. The two subbands are degenerate at $`\stackrel{}{k}=0`$, with $`E`$ symmetry. The excitation to these bands, calculated as the difference between the bandcenters $`\epsilon _{34}\epsilon _1`$, occurs at $`29.1`$ meV. Taking into account the bandwidths, this is in excellent agreement with the new mode at $`31`$ meV observed in the present experiment. Moreover, this mode has mixed in-plane and vertical character with sizable dipole elements in the normal direction. This is confirmed by the fading of the $`31`$ meV peak in the off-specular direction (see Fig.2).

The next excitation is to the band $`n=5`$, and to the $`78`$ pair of bands. The former is of $`A`$ symmetry has a bandwith of $`4.7`$ meV, while the latter has two subbands degenerate at $`\stackrel{}{k}=0`$ (with $`E`$ symmetry) and a total bandwidth of $`11.2`$ meV. The excitation energies are $`42`$ meV and $`60`$ meV, respectively. These bands have mainly an in-plane character. We suggest that the observed peak centered at $`68`$ meV could be a composite of these unresolved transitions. Another point to note about the higher bands $`78`$ is that they have considerable weight at both the fcc and hcp sites. The states in these excited bands have high mobility throughout the surface. The excitation energy $`60`$ meV needed to reach these bands is close to the activation energy observed in the quasi-elastic Helium Scattering experiment by Graham et al. . Thus, our band state results support the picture that the diffusion of H on Pt(111) happens not through classical hopping but rather via the activated tunnelling mechanism as recently proposed  to account for the similar behavior observed for H/Ni(111). The results give a direct confirmation of the protonic band picture introduced in .

When D is substituted for H in the calculations, the lower bands have a one to one correspondence with the H case with a downshift of frequency approximately by a factor of $`1/\sqrt{2}`$ as expected. For deuterium, the first excitation energy is at $`23.2`$ meV toward a composite band of width $`0.2`$ meV, in excellent agreement with the experimental observation for D as shown in Fig.1. There are transitions to higher bands at excitation energies between $`30`$ and $`60`$ meV: an excitation to an $`A`$ band at $`30.8`$ meV (bandwidth $`0.8`$ meV), one to an $`E`$ composite band at $`46`$ meV (bandwidth $`3.4`$ meV), and another to an $`E`$ composite band at $`59.8`$ meV (bandwidth $`2.9`$ meV). These could account for the broad peak centered at $`49`$ meV observed in the experiment.

There are many more excitations to the higher bands that we have not addressed here. Experimentally, they are not observed at low coverages ($`\theta 0.75`$ ML). Qualitatively, one can understand this by noting that the higher bands are much more extended than the ground states and thus the matrix elements should be much smaller. As such, the resultant weak and broad peaks are hard to observe. The picture changes at full ML because then the tunnelling of the H atom would be suppressed by on-site exclusion and effectively the H atoms are localized again. Interestingly, at a full ML, the $`31`$ meV peak disappears from the observations and two additional peaks appear at $`110`$ meV and $`153`$ meV  in the present experiment. These loss peaks agree nicely with the earlier results of Richter and Ho .

In summary, we have presented a combination of theoretical calculations and experiments on the energetics and vibrational properties of H and D adatoms on Pt(111). The experimental data for the loss peaks at low coverages ($`\theta 0.75`$ ML) differ significantly from previous studies at full H coverage. There is good agreement between the observed HREELS loss peaks and the protonic band structure calculations derived from a full 3D APES. This demonstrates that the vibrational spectroscopy of H adsorbed on metals and in particular for H/Pt(111) cannot be understood simply in terms of localized AS and SS oscillation modes at the adsorption sites, as often assumed in the literature. This picture is also crucial for a proper understanding of the recent diffusion data for H/Pt(111) . Finally, we emphasize that the dynamical behavior of this system at full ML is quite different from that at low coverages. There is also a possibility of change of adsorption sites. Further studies for the full ML case are in progress.

Acknowledgements: This work has been in part supported by the Academy of Finland through its Center of Excellence program and in part by a Galkin Fund Fellowship. We acknowledge generous computing resources from CSC – Scientific Computing Ltd., Espoo, Finland. Discussions with K. Honkala, K. E. Laasonen, M. J. Puska, R. M. Nieminen, and O. S. Trushin are also greatly appreciated. We thank J. P. Toennies for bringing our groups in contact.