text:
# Different routes to charge disproportionation in perovskites-type Fe oxides

## Abstract

Iron perovskites CaFeO3 and La0.33Sr0.67FeO3 show charge disproportionation, resulting in charge-ordered states with Fe3+:Fe5+ $`=1:1`$ and $`=2:1`$, respectively. We have made photoemission and unrestricted Hartree-Fock band-structure calculation of CaFeO3 and compared it with La0.33Sr0.67FeO3. With decreasing temperature, a gradual decrease of the spectral weight near the Fermi level occurred in CaFeO3 as in La0.33Sr0.67FeO3 although lattice distortion occurs only in CaFeO3. Hartree-Fock calculations have indicated that both the breathing and tilting distortions are necessary to induce the charge disproportionation in CaFeO3, while no lattice distortion is necessary for the charge disproportionation in La0.33Sr0.67FeO3.

Recently charge ordering phenomena in transition-metal oxides have been extensively studied, particularly in relation to charge stripes in the high-$`T_C`$ cuprates cuprates  and giant magnetoresistance in the manganites manganites . Charge ordering in highly covalent transition-metal oxides containing $`d^4`$ (Fe4+, Mn3+) or $`d^7`$ (Ni3+) ions often exhibit so-called charge disproportionation, in which a charge state is thought to be separated into two different charge states as $`2d^nd^{n1}+d^{n+1}`$. A clear fingerprint of charge disproportionation is the breathing-type distortion of metal-oxygen octahedra since the different charge states of the transition-metal ion take different ionic radii. Indeed, a breathing-type lattice distortion has been found for YNi3+O3 YNiO3CD , NdNi3+O3 NdNiO3CD  and CaFe4+O3 by neutron diffraction studies Woodward ; Takeda  and a thermally fluctuating charge disproportionated state has been postulated for LaMnO3 at high temperatures LaMnO3 . However, charge disproportionation in La1-xSrxFeO3 with $`x`$ 0.7, which occurs as a first-order phase transition from the paramagnetic average-valence state (Fe∼3.7+) above 200 K to the antiferromagnetic charge-disproportionated state (Fe3+:Fe5+ $`=2:1`$) below that temperature Takano , is not accompanied by an appreciable lattice distortion Battle  although an electron diffraction study has shown extra spots Li  and two charge states of Fe with different hyperfine fields have been detected by Mössbauer spectroscopy Takano . The neutron diffraction study of La0.3Sr0.7FeO3 Battle  has revealed a spin-density wave (SDW) of six-fold periodicity along the $``$111$``$ direction with two inequivalent spin states of Fe, suggesting a charge-density wave (CDW) of three-fold periodicity along the same direction.

Recent photoemission studies have revealed that the electronic structure of Fe4+ oxides is rather unique: the charge-transfer energy $`\mathrm{\Delta }`$, the energy required to transfer an electron from the oxygen $`p`$ to the Fe $`3d`$ level, is extremely negative ($``$ $``$3 eV including Hund’s coupling energy) SFO\_Bocquet , and the ground state of the formal Fe4+ (“$`d^4`$”) state is in fact dominated by the $`d^5\underset{¯}{L}`$ configuration, where $`\underset{¯}{L}`$ denotes a hole in the oxygen 2$`p`$ band. The charge disproportionation is therefore more correctly described as $`2d^5\underset{¯}{L}d^5\underset{¯}{L}^2+d^5`$ rather than $`2d^4d^3+d^5`$. To understand the interesting physical properties of Fe4+ oxides such as the helical antiferromagnetic metallic state in SrFeO3 SrFeO3 , the Co-substitution induced ferromagnetism in SrFe1-xCoxO3 Kawasaki ; Abbate , and magnetic and electric phase transitions under high pressure Takano2 ; Kawakami ; Rozenberg , therefore, one has to take into account the negative $`\mathrm{\Delta }`$, namely, the oxygen-hole character of charge carriers. Especially, the $`d^5\underset{¯}{L}`$ ground state naturally explains the fact that SrFeO3 shows no charge disproportionation yet no Jahn-Teller effect since the $`d^4`$ ion necessarily undergoes a Jahn-Teller distortion.

In a recent photoemission and unrestricted Hartree-Fock band-structure calculation study of La1-xSrxFeO3 Matsuno , we have shown that the charge disproportionation is purely electronically driven and that the ordering of oxygen holes plays an important role in the charge disproportionated state of La0.33Sr0.67FeO3. In this Letter, we present a photoemission study of CaFeO3, which unlike La0.33Sr0.67FeO3 shows a breathing-type lattice distortion in the charge disproportionated state. Using photoemission we have studied how the charge disproportionation influences the electronic structure near the Fermi level ($`E_F`$) in comparison with La0.33Sr0.67FeO3. We have also performed unrestricted Hartree-Fock band-structure calculations on CaFeO3 taking into account the realistic lattice distortion. Based on those results, different driving mechanisms are proposed for the charge disproportionation in the seemingly very similar systems CaFeO3 and La0.33Sr0.67FeO3.

A polycrystalline sample of CaFeO3 was prepared by a solid state reaction and a subsequent treatment under high-pressure oxygen Kawasaki . The electrical resistivity of the present CaFeO3 sample is shown in Fig. 1 and is compared with that of La1-xSrxFeO3 ($`x`$ = 0.67) Matsuno . It shows a gradual increase below the anomaly at 290 K due the gradual charge disproportionation 2Fe4+ $``$ Fe5+ $`+`$ Fe3+ while no anomaly due to the antiferromagnetic ordering is seen at the Néel temperature $`T_N=`$ 115 K.

Ultraviolet photoemission (UPS) measurements were made using the He I resonance line ($`h\nu =`$ 21.2 eV). The He I spectra have been corrected for the He I satellite. In order to calibrate binding energies and to estimate the instrumental resolution, gold was evaporated on the sample surface after each series of measurements. The energy resolution was 48 meV. The samples were repeatedly scraped in situ with a diamond file. We have adopted the spectra taken within 40 minutes after scraping and the reproducibility of the spectra was confirmed by repeated measurements.

Figure 2 shows the temperature dependence of the He I spectra of CaFeO3 near $`E_F`$. The spectra have been normalized to the integrated intensity of the entire valence band. A clear change in the intensity from the Fermi level to 0.4 eV below it has been observed as in the case of La0.33Sr0.67FeO3. One can see that the spectral weight near $`E_F`$ show a gradual decrease of the spectral weight with decreasing temperature. The most prominent temperature dependence occurs just below 300 K, followed by more gradual changes below $``$270 K. In order to illustrate this, we have plotted the spectral weight integrated from $``$0.20 to $`+`$0.05 eV as a function of temperature in Fig. 3(a), where the previous result on La0.33Sr0.67FeO3 is also plotted Matsuno . One can see a more gradual decrease as a function of temperature in CaFeO3 than in La0.33Sr0.67FeO3 corresponding to the more gradual change in the electrical resistivity below the charge disproportionation temperature, although the total spectral change in the wide temperature range is similar between the two systems. We also note that this temperature dependence of the spectral weight is qualitatively similar to the temperature dependence of the Fe-O bond length Takeda  as shown in Fig. 3(b), suggesting that the lattice deformation is related to the change of the photoemission spectra in CaFeO3. It is remarkable that, although La0.33Sr0.67FeO3 shows negligibly small structural changes across the transition Battle , it shows spectral changes as strong as CaFeO3. These observations imply different mechanisms for the charge disproportionation between the two systems.

In Fig. 2, the difference spectrum between 300 K and 100 K is compared with the Fermi-Dirac distribution function $`f_D(ϵ)`$ at $`T`$ = 300 K. The difference spectrum seems to agree rather well with $`f_D(ϵ)`$ near $`E_F`$. This is different from La0.33Sr0.67FeO3, which shows a pseudo-gap-like behavior at $`E_F`$. In order to clarify this, the difference spectrum divided by the Fermi-Dirac distribution function is also shown in the bottom panel. The obtained spectral DOS is nearly flat around the $`E_F`$, in contrast to the corresponding spectrum for La0.33Sr0.67FeO3, where a pseudo-gap-like DOS has been observed Matsuno . This is consistent with the transport properties, which show more conducting and metallic behavior in CaFeO3 than in La0.33Sr0.67FeO3 above the transition temperature.

In order to study the driving force for the phase transition and the origin of the band gap below the transition temperature in CaFeO3, we have carried out unrestricted Hartree-Fock calculations for the multi-band $`d`$-$`p`$ lattice model, in which the full degeneracy of the Fe 3$`d`$ and oxygen 2$`p`$ orbitals are taken into account Mizokawa1 . Parameters in the model are the charge-transfer energy $`\mathrm{\Delta }`$, the multiplet averaged $`d`$-$`d`$ Coulomb interaction $`U`$ and Slater-Koster parameters $`(pd\sigma ),(pd\pi ),(pp\sigma )`$ and $`(pp\pi )`$, which represent transfer integrals between the transition metal 3$`d`$ and oxygen 2$`p`$ orbitals. The values of $`\mathrm{\Delta },U`$ and $`(pd\sigma )`$ were chosen to be 0, 6 and $`1.8`$ eV, respectively. The ratio $`(pd\sigma )/(pd\pi )`$ was fixed at $`2.16`$ and $`(pp\sigma )`$ and $`(pp\pi )`$ at $`0.60`$ and 0.15 eV, respectively Mattheiss .

First, in analogy with LaSr2Fe3O9 (a supercell model for La0.33Sr0.67FeO3Matsuno , where charge ordering with Fe3+ : Fe5+ $`=2:1`$ causes an SDW with six-fold periodicity along the $``$111$``$ direction \[Fig. 4(a)\], we assumed an SDW with four-fold periodicity, namely $``$, for CaFeO3, where ordering with Fe3+ : Fe5+ $`=1:1`$ occurs in the same direction. This assumption is compatible with the alternating charge order of Fe3+ and Fe5+ in CaFeO3. Under this assumption we obtained an insulating solution but this was not a disproportionated one. The reason why we could not explain both charge-ordered and insulating ground state within the $``$ model is schematically shown in Fig. 4(b). The tendency that holes enter oxygen orbitals between the iron sites of the parallel spins also holds in this case, but all the iron sites become equivalent in spite of the hole ordering in this geometry, as can be seen from Fig. 4(b). Thus the origin of the disproportionation in CaFeO3 and in LaSr2Fe3O9 cannot be the same. In other words, the kinetic exchange interaction derived from the particular magnetic structure as in the case of La0.33Sr0.67FeO3 cannot drive the charge disproportionation in CaFeO3.

Then we introduced the lattice distortion in CaFeO3 explicitly because in CaFeO3, lattice distortion has been found to be substantial unlike La0.33Sr0.67FeO3. CaFeO3 has the GdFeO3-type tilting distortion as well as the breathing distortion Woodward ; Takeda ; Morimoto  and the Fe-O-Fe bond angle is determined to be $``$ 158. In order to take into account these distortions, a unit cell has been selected as containing four FeO6 octahedra. The breathing distortion, which results in two kinds of Fe-O bond lengths, has been treated by scaling the transfer integrals with respect to the bond lengths ($`d`$) following Harrison’s rule Harrison  $`(pd\sigma )d^{3.5}`$. In order to simulate the effect of both the breathing and the tilting, we calculated the magnitude of the band gap as a function of the degree of the breathing distortion $`R=d_{\mathrm{long}}/d_{\mathrm{short}}`$ and the Fe-O-Fe bond angle. As for the magnetic structure, we assumed the ferromagnetic structures considering that the spin alignment between neighboring iron sites is almost ferromagnetic in the screw antiferromagnetic structure with a long wave vector \[= 0.161($`111`$)2$`\pi /a`$\] observed in CaFeO3 Woodward . The results are shown in Table 1.

We have obtained the charge-disproportionated solution for $`R1`$ as expected. The experimental values for CaFeO3 are $`R=1.022`$ and $`\mathrm{}`$Fe-O-Fe $`=158^{}`$ for CaFeO3, however, the breathing distortion of $`R=1.022`$ alone or the tilting of $`\mathrm{}`$Fe-O-Fe$`=158^{}`$ alone does not open a gap although Hartree-Fock calculations tend to overestimate band gaps. In order to open a band gap (and to lower the total energy of the system), the two types of lattice distortions should occur simultaneously. If the bond angle is 180 as in the case of SrFeO3, one needs an unreasonably large breathing distortion to open the gap. In order to confirm the co-operative nature of the two types of distortions, first-principles total-energy calculations would be necessary in future work.

The present scenario for the charge disproportionation in CaFeO3 naturally explains why the isoelectronic SrFeO3 remains free from charge disproportionation: since the ionic radius of Sr2+ is large, no tilting distortion is possible for SrFeO3 and hence the breathing is also blocked. This picture may be analogous to the case of the distorted perovskite BaBiO3, in which freezing of a breathing phonon mode is well known BKBO\_review . Liechtenstein et al. BKBO  calculated the total energy of BaBiO3 as a function of tilting and breathing distortions and found that the instability and the gap opening occurs only when the two kinds of distortions are combined; in the presence of the tilting distortion, the nesting instability of the Fermi surface emerges, leading to the alternating breathing distortion of Bi-O octahedra. If the same nesting instability is confirmed for CaFeO3 by first-principles band-structure calculations, then the quite different mechanism of the charge disproportionation would be established between La0.33Sr0.67FeO3, where the ordering of holes at oxygen sites is purely electronically driven, and CaFeO3, where the lattice distortion and associated electron-phonon coupling is important.

In conclusion, we have shown that the driving force for the charge disproportionation is quite different between the two apparently very similar compounds La0.33Sr0.67FeO3 and CaFeO3. The present finding implies that there exists subtle interplay between electron-electron interaction, magnetic interaction, and electron-lattice interaction in realizing charge ordering (or more generally charge inhomogeneity) and charge fluctuations in transition-metal oxides. This means that dominant interaction may change within the same family of compounds, and calls for a critical re-examination of charge ordering/fluctuation phenomena depending on the chemical composition (hole concentration) and crystal distortion even for an apparently similar class of materials.

Discussions with M. Seto and N. Hamada are gratefully acknowledged. This work was supported by a Grant-in-Aid for Scientific Research (A12304018) from the Ministry of Education, Culture, Sports, Science and Technology.


format: text

text:
# References

It is a pleasure and an honor to welcome participants and guests of the Third International Sakharov Conference on Physics here in the Lebedev Physical Institute. There is a number of intimate connections. Dmitry Sakharov, father of Andrei Sakharov, was a student of Petr Lebedev in the Physical Department of Moscow University. Pressure of light first observed by Lebedev was 40 years later used by Sakharov to trigger the H-bomb. There are also moral parallels: Petr Lebedev, who was more than far from politics, being a personality of high ethic principles could not help resigning from the University in a protest against Tsar Government repressions towards his colleagues. His Lab was destroyed, and although private sponsors began to build for him ”Lebedev institute” his heart did not stand the disaster, he died in March 1912 in the age 46, and supposedly because of it did not receive Nobel Prize on Physics he was nominated to for his discovery.

Actually this Conference should have taken place a year ago, being dedicated to 80 years of Sakharov who was born in 1921. The reason for this delay was serious: in May 1999 Professor Efim Fradkin died, and in June 2000 Theoretical Physics Department organized in his memory the big international conference ”Quantization, gauge fields and strings”. To organize within a year one more huge conference was impossible. There is no doubt that Andrei Dmitrievich would approve of our actions in the memory of Efim S. Fradkin, whom he highly respected as personality and as scientist. He also very much respected Professor David Kirzhnits who died recently. It was David Kirzhnits who in his memoirs about Sakharov (see in ) most transparently compared Sakharov’s methodology with the feat of Russian pilot Konstantin Artseulov during the First World War, who the first time in history of aviation deliberately dropped down his airplane into the mortal spin (thus committing a suicide in public, according to the hundreds of observers’ opinion) and safely went out of it, creating the method which saved lives of hundreds of pilots. You see: it is always quite interesting to speak about Sakharov whose ideas and actions were really very often considered as ”mortal”, but they proved to be salutary for all of us after all.

In this short Opening Word I’ll concentrate upon three ”dynamical characteristics” of Sakharov’s ”method” and his way of thinking. The first one, which may be called, ”Artseulov method” I already described above with the help of David Kirzhnits’s metaphor: after Sakharov came to certain conclusion he preferred to act most resolutely and also as a good teacher. But to come to the conclusion was not a simple dynamical process. Second feature of Sakharov’s mentality may be called ”a permanent feeling of possible personal error”. He writes in ”Memoirs”:

”My statements on general issues are often tentative, meant to provoke discussion, and subject to revision. I agree with Leszek Kolakowski when he writes: ”Inconsistency is simply a secret awareness of the contradictions of this world… a permanent feeling of possible personal error, or if not that, then of the possibility that one’s antagonist is right”. My only quarrel, - comments Sakharov after this quotation of Kulakowski, - is with the word ”inconsistency”, which I would replace with one that conveys my belief that intellectual growth and social awareness should combine dynamic self-criticism and a set of stable values… I am not a professional politician. Perhaps that is why I am always burdened by doubt about the usefulness and consequences of my actions. I inclined to the belief that a combination of moral criteria and unrestricted inquiry provides the only possible compass” (, pp. 579-580).

It is difficult to count how many times these Sakharov’s sincere words ”I am not a professional politician…” were used against him by Soviet propaganda. But we are talking now about something much more serious.

Sakharov really was a ”thinking matter”, everybody who personally knew him, I also was happy to contact him during 20 years, will confirm that his brain worked permanently and any moment he was ready to contemplate the problems of interest ”as if he had a blank sheet of paper in front of him”, - as Igor Evgenievich Tamm, his teacher and many years Chief of the Theoretical Physics Department, put it. The same was true with regards to Sakharov’s own opinions and views: that is why it was so interesting to speak to him. But then the important question arises: If there were permanent doubts and ”unrestricted inquiry” (in Sakharov’s words cited above) how he managed to come to definite conclusions, why he was not a sort of ”Chekhov’s hero” incapable of any decision-making? Sakharov definitely was not a ”Chekhov’s hero”, but his decision-makings were often rather difficult. It is not by chance that in quotations above he speaks about ”stable values” and ”moral criteria”, i.e. about invariants which may serve a certain guideline. Idea of absolute priority of individual human rights - above any social, national, religious and other so called ”great ideas” - was a product of permanent internal search of such an invariants.

But to explain why he reached so much it is also necessary to point out one more feature of his mentality \- his special talent which is very much advisable although not easy to learn from him: he was good at high mathematics, he was able, so to say, ”to count to two”. The holism (in the Oriental language) of his thinking is well known. He could consider the problem at hand all at once, in all its complexity, with all its particulars, and, this being most important, in its dynamics. A space-time picture was obtained right away, with expected answer in the end. And all this, including actions already fulfilled or planned ones, was contemplated over again and again. Why I call this ability ”to count to two”? During half of XX-th century half of Mankind admired Soviet socialism with its ideas of social justice (ONE) being incapable of noticing terror, GULAG etc. (TWO). Sakharov’s idea of convergence, although being questionable by itself, was a sigh of relief for many and many.

One more example: ”peacemaking” by western intellectual elite (there is no sense to argue against the good goal - peace is better than war, it is highest humane priority - ONE) in practice unfortunately was helping to push the Mankind closer to the thermonuclear abyss because ”elite” did not manage to notice another side of the medal - aggressiveness of militarized Soviet closed society and hence necessity for the West (just for the sake of keeping and strengthening Peace) to be sufficiently strong (TWO). It was Sakharov who had a courage to reveal the truth speaking from inside the USSR, who supported ”double solution” of NATO, which included additional deployment of rockets with nuclear warheads targeted at the Soviet cities (isn’t it the Artseulov method?). Of course because of it he himself became a target of severe revenge of Soviet authorities; then ”elite”, colleagues supported Sakharov. As a result of this dialectic ”Perestroika” began miraculously, Ronald Reagan came to the Red Square and rockets with nuclear charges began to be destroyed after all. (This time Mankind avoided the fall into the abyss, although it was real narrow escape; Sakharov was quite aware of the unacceptably high probability of the ”fall”).

Ability to see simultaneously both sides of a medal, to possess this ”mental mirrors” is really a rear talent. This Sakharov’s talent ”to count to two” now is not less demanded than before. But some of you in this Hall may ask why I cover this moral and political issues at the physical conference? The answer is two-fold:

(1) Because this is ”Sakharov” conference. Harry Lipkin from the Weizmann Institute of Science writes in : ”Andrei Dmitrievich had the remarkable ability to understand how systems function - all systems, social, political, scientific, technological, as well as interfaces between them. In a sense he was a kind of an interdisciplinary systems super-engineer”. That is why it seems reasonable even at this conference to exceed limits of physics, at least in short.

(2) Because we all know the value of speed of light, but what is the value of ”speed of darkness”? Will the Mankind manage to win the race with mortal dangers capable of ”putting off the Sun” for our civilization?

Sakharov died almost 13 years ago. Time moves on and now there are new hopes and new serious threats. World terror with its new weapon - suicide bombers, who potentially may be carriers of biological or nuclear tools of mass destruction. In this connection I draw attention to the most important and well balanced Statement by Elena Bonner (April 2002) where she outlined the global danger of this new weapon. Again plenty of people prove to be incapable ”to count to two”, got lost in two simple notions: ”struggle with terror” and ”human rights”. Those who advocate active struggle with world terror sometimes have evident difficulties with embracing demands of human rights and humanity toward peaceful population - in this way deliberately or undeliberately supporting terrorists (unceasing drama in Chechnya is a dreadful example). Others - like now in Europe - seem to advocate human rights but ignore the vital necessity of taking active measures against terrorists - in this way also actively supporting terrorists and hence most terrible violations of human rights. I feel obliged to tell that members of the Theoretical Physics Department joined the collective protest of scientists against Israel-baiting campaign in Europe. This campaign visually shows that we are still too close to the precipice, and that Sakharov’s talent to comprehend two sides of a medal simultaneously is strongly demanded, as well as his understanding that individual human rights must be considered above any national, religious and other corporate goals.

There are another dangers. General admiration of achievements in biophysics and gene engineering is accompanied by concerns about possibility of creating artificial virus capable to kill the civilization on the planet of Earth (is not our life ”a permanent narrow escape”?). Steven Hawking recently said in an interview that he considers the death of Mankind because of such viruses inevitable and the only way out is our settling in cosmos. In this connection I would like to suggest perhaps even more secure option of settling the Mankind in the extra dimensions; this may be discussed at the Conference. In physics, like in the life of human society, modern days seem to be also boiling and interesting. Of course Sakharov would be happy to participate in this Conference. And I wish it success in its work.


text:
# 1 Introduction and formulation of the model

## 1 Introduction and formulation of the model

Spatial and temporal solitons in various models with quadratic ($`\chi ^{(2)}`$) ) nonlinearity, including the second-harmonic-generating (SHG) and more general three-wave (3W) systems, have attracted a lot of attention, see reviews  and a special volume . The simplest (degenerate) SHG model involves only two waves, the fundamental-frequency (FF) one and its second harmonic (SH). In a general situation, one is dealing with a 3W system, that involves two “daughter” waves and a pump; it is known that the 3W system also gives rise to soliton solutions . A case of great practical interest is a more special version of the 3W system, known as the Type-II SHG model, in which both “daughters” represent two polarizations of the FF wave, while SH has a single polarization. The latter model corresponds to the most typical experimental conditions , and its three-component soliton solutions (also called vectorial solitons) were studied in detail .

Various modifications of the Type-II system were introduced, and the soliton solutions in them were studied \[see, e.g., Refs. \]. These modifications were focused on adding extra quadratic terms to the model’s equations. Another possibility which has not yet been considered, except for Ref.  (see below), is to introduce linear coupling (mixing) between the two FF components. The simplest physical realization of this possibility may be found in terms of the Type-II SHG process in a fiber-like birefringent waveguide, which is subjected to a twist that mixes the two FF polarizations. Note that twist-induced linear mixing between two linear polarizations is a well-known feature of optical fibers with the Kerr ($`\chi ^{(3)}`$) nonlinearity . Of course, application of the twist to a monocrystalline waveguide in which SHG takes place is problematic, but an effective twist, without applying any mechanical torque, can be created in an evident way in the practically important case when SHG itself is induced by means of periodic poling of the host medium (see a review ).

A system of equations to describe the Type-II SHG in a birefringent medium in the presence of the linear mixing can be derived as a straightforward generalization of the model developed in Ref. :

$`i\delta A_z+ib_1A_\tau +\left(1/2\right)A_{\tau \tau }\beta A+B^{}C`$ $`=`$ $`\kappa B,`$ (1)
$`i\delta ^1B_z+ib_2B_\tau +\left(1/2\right)B_{\tau \tau }\beta ^1B+A^{}C`$ $`=`$ $`\kappa A,`$ (2)
$`i\sigma C_z+\left(1/2\right)DC_{\tau \tau }\alpha C+AB`$ $`=`$ $`0,`$ (3)
where the asterisk and subscripts stand, respectively, for the complex conjugation and partial derivatives, the evolution in the temporal domain is implied, $`\tau `$ being the standard reduced-time variable , complex fields $`A(z,\tau )`$, $`B(z,\tau )`$ and $`C(z,\tau )`$ represent the two FF components and SH wave, respectively, $`\beta `$ is a phase-birefringence parameter, $`\alpha `$ measures the phase mismatch between the FF and SH waves, $`\delta >0`$ is another parameter taking into regard asymmetry between the two FF components , $`D`$ is a relative SH/FF dispersion coefficient, and $`\sigma >0`$ determines a relative propagation constant, all these parameters being real.

New (in comparison with the known model) terms introduced in Eqs. (1) and (2) are the linear-mixing ones with the real coefficient $`\kappa `$, and group-velocity-birefringence (temporal-walkoff) terms with real coefficients $`b_{1,2}`$ \[the reference frame is defined so that the group-velocity term vanishes in Eq. (3)\]. Note that in the absence of the linear mixing, $`\kappa =0`$, the walkoff terms can be easily removed from the model by means of phase transformations and change of the reference frame, therefore the model elaborated in Ref.  did not include these terms. However, these terms are irreduceable if $`\kappa 0`$.

The model may also be interpreted in the spatial domain, with $`\tau `$ is replaced by the transverse coordinate $`x`$, $`b_{1,2}`$ being spatial walkoff  parameters. In fact, the linear coupling between the FF components in the 3W model in the spatial domain (a planar optical waveguide with the $`\chi ^{(2)}`$ nonlinearity) may be induced in a simple way by means of a Bragg grating (a system of parallel scores) written on the waveguide . In the latter case, all the fields have the same polarization, the difference between the two FF components being in the direction of their Poynting vectors.

An objective of this work is to find fundamental-soliton solutions of the model and test their stability by means of precise numerical simulations. Prior to that, we will consider the model’s linear spectrum, which is necessary in order to realize what type of solitons may exist in it.

For $`\kappa 0`$ but $`b_{1,2}=0`$, stationary soliton solutions to Eqs. (1), (2), and (3) can be easily found. A dynamical test will show that the stationary solitons may be both stable and unstable in this case, depending on the relative sign of the two FF components, which is a difference from the results reported in Ref.  for the Type-II SHG model without the linear mixing ($`\kappa =0`$). With the introduction of the walkoff parameters ($`b_{1,2}0`$), numerical solution of the stationary equations becomes difficult, therefore we rely upon direct simulations in this case: propagation of an initial solitary-wave configuration, which is taken as a stable soliton for the same values of parameters but with $`b_1=b_2=0`$, generates a stable fundamental soliton, featuring intrinsic chirp, after a transient process. Additionally, we also find quasi-stable double-humped solitary waves in the latter case.

## 2 The linear spectrum

Solitons are represented by solutions to Eqs. (1) - (3) of the form

$$A(z,\tau )=e^{iKz}a(\tau ),B(z,\tau )=e^{iKz}b(\tau ),C(z,\tau )=e^{2iKz}c(\tau ),$$
(4)
where $`K`$ is a real propagation constant, and $`a(\tau )`$, $`b(\tau )`$, and $`c(\tau )`$ are even functions (complex ones, in the general case) that must vanish at $`\left|\tau \right|\mathrm{}`$. As is commonly known, (bright) soliton solutions may only exist at values of $`K`$ that do not overlap with the continuous spectrum of the FF part of the system. The spectrum is determined by the substitution of the continuous-wave solution,

$$A(z,\tau )=A_0\mathrm{exp}\left(iKzi\omega \tau \right),B(z,\tau )=B_0\mathrm{exp}\left(iKzi\omega \tau \right),$$
(5)
into linearized versions of Eqs. (1) and (2).

As for the SH component, a similar condition for the existence of solitons (the non-overlapping with the continuous spectrum) is not always necessary in it, since situations are possible when the SH equation may *not* be linearized, its quadratic term being everywhere (at all values of $`\tau `$, including $`\left|\tau \right|\mathrm{}`$) on the same order of magnitude as the linear terms. Such a situation is possible as the quadratic and linear terms in the SH equation are composed of different fields: the former one contains only the FF components, while the linear terms are expressed in terms of the SH field, see Eq. (3). Eventually, this leads to a possibility of the existence of the so-called embedded solitons, for which the SH propagation constant $`2K`$ \[see Eq. (4)\] indeed falls into the linear spectrum of the SH wave .

The combination of the linear coupling and group-velocity difference between the two FF components in Eqs. (1) and (2) makes the present system somewhat akin to $`\chi ^{(2)}`$ models in which effective dispersion and/or diffraction are induced by a Bragg grating, while the intrinsic dispersion (or diffraction) are neglected. Models of this type, which have a finite gap in their linear spectrum and, accordingly, give rise to  gap solitons supported by the quadratic nonlinearity, have attracted considerable attention, see Refs. . However, the present model, despite its similarity to the gap-soliton ones, has the linear spectrum of a different type, which contains no finite gap but, instead, the usual semi-infinite gap extending to $`K+\mathrm{}`$. Indeed, taking, in order to avoid ponderous formulas, Eqs. (1) and (2) with $`\delta =1`$, $`\beta =1`$, $`b_1=b_2b`$, and substituting the expressions (5) into the linearized equations, one can easily obtain the following expression for the linear spectrum:

$$K=\left(1+\omega ^2/2\right)\pm \sqrt{\left(b\omega \right)^2+\kappa ^2}.$$
(6)
It is straightforward to see that the combination of two branches of this spectrum does not yield any finite gap \[that would exist, in the form $`\left(K+1\right)^2<\kappa ^2`$, if the dispersion terms in Eqs. (1) and (2), and hence the term $`\omega ^2/2`$ in Eq. (6), were omitted\]. On the other hand, all sufficiently large positive values of $`K`$ do not belong to the linear spectrum, forming the above-mentioned semi-infinite gap (in the case $`\kappa ^2>b^4`$, it simply takes the form $`K>0`$).

Thus, the present system combines, to a certain extent, previously studied ordinary soliton models and those which give rise to gap solitons, which makes it interesting to search for solitons in it.

## 3 Solitons in the model without intrinsic walkoff

### 3.1 Unstable solitons

We start by seeking for soliton solutions to Eqs. (1) - (3) in the case $`b_{1,2}=0`$ (no intrinsic walkoff, while the linear mixing is present, $`\kappa 0`$). The stationary version ($`/z=0`$) of the system of Eqs. (1), (2) and (3) then becomes real, and, accordingly, stationary solutions are looked for in a real form by means of the shooting method. Note that, in the absence of the linear mixing ($`\kappa =0`$), solutions differing by a relative sign of the FF components, $`A`$ and $`B`$, are, obviously, equivalent. However, this is not the case if $`\kappa 0`$ (which is similar to the situation in a $`\chi ^{(2)}`$ model of a dual-core waveguide, where linear terms couple FF fields in two cores ). In particular, it is obvious that Eqs. (1) – (3) can be derived from a Hamiltonian, its term which accounts for the linear mixing being

$$H_{\mathrm{mix}}=\kappa _{\mathrm{}}^+\mathrm{}\left(A^{}B+AB^{}\right)𝑑\tau .$$
(7)
In the case of real solutions, and taking, by definition, $`\kappa >0`$, one may expect that the solutions with $`\mathrm{sgn}A=\mathrm{sgn}B`$ should be unstable, as they make the term (7) of the Hamiltonian positive, while solutions with $`\mathrm{sgn}A=\mathrm{sgn}B`$ may be stable, as they yield $`H_{\mathrm{mix}}<0`$ (a known principle states that a solution which makes the value of the Hamiltonian larger is likely to be unstable ).

The application of the shooting method to the stationary real equations (1) - (3) with $`b_{1,2}=0`$ indeed yields stationary solitons in a broad parametric region, provided that the relative dispersion coefficient is positive, $`D>0`$. A typical example of a solution with $`\mathrm{sgn}A=\mathrm{sgn}B`$, generated by the shooting method, is displayed in Fig. 1.

Direct simulations of the dynamical stability of the thus found stationary solitons were performed by means of the implicit Crank-Nicholson scheme. As is known , this scheme has advantages over the split-step method for very long beam-propagation simulations. To control the accuracy of the direct simulations, we made use of the fact that the underlying system (1) – (3) conserves the net energy of the three fields,

$$E=_{\mathrm{}}^+\mathrm{}\left(\delta \left|A\right|^2+\delta ^1\left|B\right|^2+\sigma \left|C\right|^2\right)𝑑\tau .$$
(8)
Results of the simulations complied with the conservation of the integral (8) to a very high accuracy.

The application of the Crank-Nicholson scheme to the soliton shown in Fig. 1 has produced a picture which is displayed, in the contour-plot form, in Fig. 2. This result, as well as many other runs of simulations for other stationary solitons, demonstrate that, in accordance with the general argument presented above, all the real solitons with $`\mathrm{sgn}A=\mathrm{sgn}B`$ are unstable. It is relevant to compare a characteristic propagation length $`z_{\mathrm{stab}}`$ before the onset of the instability, which is $`80`$ in Fig. 2, and the soliton period for the same pulse (which is defined as the propagation length necessary for the change of the internal phase of the soliton by $`\pi /2`$ ), that can be estimated for the soliton shown in Fig. 1, in terms of its width $`W`$, as $`z_{\mathrm{sol}}\pi W^215`$. This means that the unstable soliton can pass $`5`$ soliton periods as a quasi-stable object. In a typical experimental situation for spatial solitons, their diffraction length is $`1`$ mm, while the size of a sample is a few cm , hence both the unstable soliton itself and its instability may be experimentally observed in the spatial domain.

In fact, Fig. 2 displays a typical scenario of the development of the instability of solitons with $`\mathrm{sgn}A=\mathrm{sgn}B`$: after a relatively long period of a very slow “latent” growth of the instability, an abrupt explosion occurs, as a result of which the soliton sheds off some amount of radiation in the SH component, and rearranges itself into a new, completely stable, soliton (these stable solitons will be described below). Figure 3 shows the quasi-stable-propagation distance $`z_{\mathrm{stab}}`$ of the unstable solitons vs. the coupling constant $`\kappa `$, for several different values of the relative SH dispersion coefficient $`D`$. The fact that, irrespective of the value of $`D`$, the distance $`z_{\mathrm{stab}}`$ diverges as $`\kappa 0`$ is quite natural, as in this limit we get back to the usual model without the linear mixing, where the real solitons are stable irrespective of the relative sign of their two FF components .

If SH dispersion coefficient $`D`$ is very small, the character of the instability development becomes qualitatively different from that illustrated by Fig. 2: as is seen in Fig. 4, in this case only the SH component of the soliton survives. Note that, in the case $`D=0`$, Eqs. (1), (2), and (3) have an obvious solution corresponding to the eventual state observed in Fig. 4: $`A=B=0`$, $`C(z,\tau )=\mathrm{exp}\left(i\alpha \sigma ^1z\right)c(\tau )`$, where $`c(\tau )`$ is an arbitrary function.

### 3.2 Stable solitons

In the case $`b_1=b_2=0`$, another class of real stationary soliton solutions can be found (by means of the shooting method too), with opposite signs of the two FF components ($`\mathrm{sgn}A=\mathrm{sgn}B`$). A typical example of such a soliton is shown in Fig. 5. In accordance with the qualitative arguments given above \[based on the sign of the coupling term (7) in the Hamiltonian\], the solitons of this type are found to be completely stable in direct simulations of their evolution. An example is displayed in terms of the contour plots in Fig. 6 (note that the full propagation distance presented in Fig. 6 is $`100`$ soliton periods; in fact, the stability was seen, in much longer simulations, to persist indefinitely).

An accurate analysis of the stable pulses, the formation of which was observed as a result of the development of the instability of the stationary solitons with $`\mathrm{sgn}A=+\mathrm{sgn}B`$ (see Fig. 2), shows that the appearing stable pulses are identical to the stable solitons of the type considered here, with $`\mathrm{sgn}A=\mathrm{sgn}B`$. Thus, these solitons are really robust, playing a role of *attractors* in the evolution of unstable pulses (the existence of effective attractors in a conservative nonlinear-wave system is possible due to radiation losses).

## 4 Solitons in the system with the intrinsic walkoff

### 4.1 Single-humped solitons

It is quite interesting to understand how solitons are modified if the group-velocity (walkoff) terms are restored in Eqs. (1) and (2). Search for stationary soliton solutions of the full system of Eqs. (1) \- (3), which include the walkoff terms, turns out to be much harder than in the case $`b_1=b_2=0`$, as it appears quite difficult to secure convergence of results produced by the shooting method. Therefore, we adopted an approach based on direct simulations of the full equations in the following fashion: stationary solutions corresponding to stable solitons (with $`\mathrm{sgn}A=\mathrm{sgn}B`$), which were found above for the case $`b_1=b_2=0`$ (for instance, the soliton shown in Fig. 5), were used as initial conditions for simulations of the evolution equations (1) – (2) with the same values of all the parameters but $`b_1`$ and $`b_2`$. Note that the introduction of the walkoff terms must essentially rearrange the input pulses, as they are real, while stationary solutions to Eqs. (1) – (2) with $`b_{1,2}0`$ cannot be real.

A typical example of the evolution of the thus chosen input pulse is shown, by means of contour plots, in Fig. 7 \[in this figure, we display the evolution of local powers $`\left|A(\tau )\right|^2`$, $`\left|B(\tau )\right|^2`$, and $`\left|C(\tau )\right|^2`$\]. The propagation distance in Fig. 7 is extremely large ($`200`$ soliton periods), which was taken in order to make it sure that a final soliton, if any, takes a sufficiently well-established form. The result of the evolution is shown in Fig. 8.

A conclusion suggested by this and many other runs of simulations is that stable solitons with intrinsic chirp establish themselves in the presence of the intrinsic walkoff, although radiation shed off from the soliton in the course of its self-adjustment separates from it very slowly. The latter peculiarity can be understood. Indeed, the radiation tail attached to the soliton in Fig. 8 (upper panel) is all built of the SH field, which has zero group velocity in the underlying equation (3), hence it does not readily separate from the zero-velocity soliton. It seems very plausible (although detailed consideration of the issue is beyond the scope of this work) that Eqs. (1) – (3) may also generate moving (“walking”) solitons, in which case the separation of the soliton from the radiation “garbage” would probably be faster.

To further check that the (quasi-) soliton (called this way because of the radiation tail attached to it), whose formation and structure are shown in Figs. 7 and 8, has long since completed any essential evolution, in Fig. 9 we show the evolution (vs. $`z`$) of the net energy (or intensity, in the case of beams in the spatial domain) of each component of the soliton, i.e.,

$$_{\mathrm{}}^+\mathrm{}\left|A(z,\tau )\right|^2𝑑\tau ,_{\mathrm{}}^+\mathrm{}\left|B(z,\tau )\right|^2𝑑\tau ,_{\mathrm{}}^+\mathrm{}\left|C(z,\tau )\right|^2𝑑\tau $$
(9)
($`_{\mathrm{}}^+\mathrm{}`$ is realized as the integral over the whole simulation domain), and Fig. 10 displays the evolution of the total energy defined by Eq. (8). A very small initial loss of the total energy (see Fig. 10) is explained by leakage across borders of the integration domain in the process of the initial rearrangement of the soliton. Note that intensive energy exchange between the $`A`$ and $`C`$ fields (see Fig. 9) is limited to approximately the same initial stage of the evolution at which the energy loss takes place; then, any tangible evolution ceases, in terms of the integral field characteristics (9).

### 4.2 Double-humped structures

The approach described above produces stable fundamental solitons, i.e., single-humped ones. On the other hand, it is well known that $`\chi ^{(2)}`$ models readily give rise to higher-order solitons – first of all, double-humped ones  – which, however, are always unstable in standard models, including the 3W Type-II model . Search for *stable* double-humped solitons in various systems is a problem of considerable interest for physical applications, see, e.g., Refs. . In fact, the first examples of (numerically) stable one-dimensional double-humped solitons were found in 3W models combining $`\chi ^{(2)}`$ nonlinearity and linear coupling (which was induced by the Bragg grating) between two components of the FF field . Moreover, the model introduced in Ref. , that seems to be closest to the one considered in the present work, gives rise also to vast families of double- and multi-humped embedded solitons , although the stability of those solutions was not studied in detail.

We made an attempt to search for double-humped solitary-wave structures in the present model. In the absence of the linear mixing ($`b_{1,2}=0`$), they have never been found, which seems quite natural in view of the above-mentioned results obtained in allied models. However, structures of that type can indeed be found at finite (actually, quite small) values of $`\left|b_{1,2}\right|`$, and they seem to be nearly stable, although they are not completely stationary.

To this end, the stationary version of Eqs. (1) - (3), produced by the substitution of the waveforms (4), was first solved numerically with high but finite accuracy, starting from the numerically exact solution for $`b_1=b_2=0`$, such as the one shown in Fig. 5, and gradually increasing the parameter $`b_1=b_2b`$. As it was mentioned above, in the presence of the walkoff terms with $`b_{1,2}0`$ straightforward application of the shooting technique does not provide for convergence of soliton solutions to indefinitely high accuracy, this is why the accuracy was finite, as mentioned above. It was observed that if other parameters keep constant values (for instance those which are mentioned in the caption to Fig. 11), the increase of $`b`$ makes the (finite-accuracy) soliton broader, and *splitting* of the soliton’s crest into two takes place at $`b=0.00274`$. A typical example of the appearing double-humped structure is shown, for a slightly larger value $`b=0.003`$, in the upper panel of Fig. 11.

The simulated evolution of this structure over a very long propagation distance shows that this solitary wave is not a genuine steady-state solution, but it is quite close to being one. It keeps a well-pronounced double-humped shape over, at least, $`15`$ soliton periods. This implies that the double-humped structure is robust enough to be observed in an experiment.

Further increase of $`b`$ makes the double-humped pulses still less localized, and, eventually, permanent leakage of one of the FF components from the pulse starts. It is difficult to find a critical value of $`b`$ at which this pulse ceases to exist as a solitary-wave solution, as an extended “tail” of the FF field, the appearance of which signals the onset of the leakage, has a vanishingly small amplitude when it emerges.

## 5 Conclusion

We have proposed a modification of the usual three-wave second-harmonic-generation model which incorporates two features that are new to the usual model: linear mixing between two components of the fundamental-frequency wave, and a group-velocity mismatch (walkoff) between them. Although the new system is akin to gap-solitons models, its linear spectrum contains no finite gaps. In the temporal domain, the model may be interpreted as the one adding an (effective) twist of the fiber-like waveguide to the birefringence, the latter feature being typical for the Type-II $`\chi ^{(2)}`$ systems. In the spatial domain, the two FF components differ, physically, not by their polarizations, but rather by the orientation of their Poynting vectors in a planar waveguide, the linear coupling being induced by the Bragg grating.

In the absence of the intrinsic walkoff, the linear mixing induces a difference between real soliton solutions with the opposite relative signs between the two FF components, so that they are stable for one sign, and unstable for the other. The development of the instability leads to rearrangement of unstable solitons into stable ones. Adding the intrinsic-walkoff terms, we have found that the evolution leads to formation of stable chirped fundamental solitons, and, additionally, quasi-stable double-humped solitary waves were found.

## Acknowledgement

This work was supported in a part by a grant from the Research Authority of the Tel Aviv University.

Figure Captions

Fig. 1. An example of a stationary real soliton solution with $`\mathrm{sgn}A=\mathrm{sgn}B`$. The parameters are $`\delta =2`$, $`\beta =0.778`$, $`\kappa =0.2`$, $`\sigma =3`$, $`D=1`$, $`\alpha =0.156`$, $`b_{1,2}=0`$ (no walkoff between the two fundamental-frequency wave components). In this and subsequent figures, the argument $`x`$ attached to the horizontal axis replaces the variable $`\tau `$ for a case when the model is interpreted in the spatial domain, where $`x`$ is the transverse coordinate (see the text).

Fig. 2. Evolution of the soliton shown in Fig. 1.

Fig. 3. The distance necessary for the onset of the instability of the soliton with $`\mathrm{sgn}A=\mathrm{sgn}B`$ (as detected in the second-harmonic component) vs. the linear-mixing constant $`\kappa `$.

Fig. 4. Evolution of the soliton in the case $`D=0.03`$, other parameters taking the same values as in Fig. 1.

Fig. 5. An example of a stationary real soliton with $`\mathrm{sgn}A=\mathrm{sgn}B`$ for the same values of parameters as in Fig. 1.

Fig. 6. Evolution of the soliton shown in Fig. 5.

Fig. 7. Evolution of the input pulse identical to the soliton shown in Fig. 5 at the same values of parameters as in Fig. 5, except for $`b_1=b_2=0.1`$ (cf. Fig. 6, which pertains to the case $`b_{1,2}=0`$).

Fig. 8. Panels (a) and (b) display, respectively, the distribution of local powers, $`\left|A(\tau )\right|^2`$, $`\left|B(\tau )\right|^2`$, and $`\left|C(\tau )\right|^2`$, and local chirps, $`\left(\varphi _A\right)_{\tau \tau }`$, $`\left(\varphi _B\right)_{\tau \tau }`$, and $`\left(\varphi _C\right)_{\tau \tau }`$ ($`\varphi `$ stands for the phase of field), of the three waves in the (quasi)soliton generated by the evolution process displayed in Fig. 7.

Fig. 9. The net energies of the three components of the soliton, defined as per Eq. (9), vs. the propagation distance $`z`$, for the same case as in Fig. 7.

Fig. 10. The total energy (8) of all the three fields vs. $`z`$, shown for the same case as in Fig. 9. Note that the energy loss, due to some leakage through the edges of the integration domain, is very small (see numerical values on the vertical axis).

Fig. 11. An example of a quasi-stable double-humped solitary-wave structure found for $`\delta =2`$, $`\beta =0.778`$, $`\kappa =0.2`$, $`\sigma =4`$, $`D=1`$, $`\alpha =0.156`$, and $`b_1=b_2=0.003`$. The upper and lower panels show, respectively, the initial configuration at $`z=0`$, obtained as a finite-accuracy shooting solution of the stationary equations, and the final configuration obtained at $`z=4000`$. As well as in Fig. 8, the distribution of the field powers across the solitary wave is shown here.


text:
# Power of the power-laws: lessons from unification of small and large time scales for evolution

## Abstract

We develop a “unified” model that describes both “micro” and “macro” evolutions within a single theoretical framework. The eco-system is described as a dynamic network; the population dynamics at each node of this network describes the “micro”-evolution over ecological time scales (i.e., birth, ageing and natural death of individual organisms) while the appearance of new nodes, the slow changes of the links and the disappearance of existing nodes accounts for the “macro” evolution over geological time scales (i.e., the origination, evolution and extinction of species). In contrast to several earlier claims in the literature, we observe strong deviations from power law in the regime of long life times where the statistics is, usually, poor.

The recent surge in the modeling of biological evolution and extinction of species, using the concepts and techniques of statistical physics, has been stimulated partly by the claims (see drossel ; newman ; kauffman ; bak  for reviews) that the statistical distributions of several quantities associated with the extinction of species follow power laws. However, almost all of these models focus only on the “macro” evolution (i.e., the evolution of species on geological time scales). Neither the birth, ageing and, eventually, the death of the individual organisms nor the detailed population dynamics make explicit appearance in these theoretical descriptions. On the other hand, in reality, a species becomes extinct when its entire population is wiped out.

Therefore, we develop a “unified” model of an eco-system that describes both “micro” and “macro” evolutions. The eco-system is described as a dynamic network. The “micro”-evolution over ecological time scales, i.e., birth, growth (ageing) and natural death of individual organisms is described by the dynamics within each node staubook . The network itself evolves slowly with time; over sufficiently long time scales populations of some species would drop to zero, indicating their extinction, and the corresponding nodes would be deleted from the network. On the other hand, appearance of new nodes, together with their own population of individual organisms, signals origination of new species. In addition, the links of the network also change slowly to capture the adaptive evolution of the species by random mutations over geological time scales.

The dynamic network: At any arbitrary instant of time $`t`$ the model consists of $`N(t)`$ species each of which may be represented by one of the $`N`$ nodes of a dynamic network; the total number of nodes is not allowed to exceed $`N_{max}`$. Our model allows $`N(t)`$ to fluctuate with time over the range $`1N(t)N_{max}`$. The population (i.e., the total number of organisms) of a given species, say, $`i`$ at any arbitrary instant of time $`t`$ is given by $`n_i(t)`$. The limited availablity of resources, other than food, in the eco-system imposes an upper limit $`n_{max}`$ of the allowed population of each species. Thus, the total number of organisms $`n(t)`$ at time $`t`$ is given by $`n(t)=_{i=1}^{N(t)}n_i(t)`$. Both $`N_{max}`$ and $`n_{max}`$ are time-independent in the model.

The interactions: Prey-predator interactions are captured through the matrix $`𝐉`$. The influence of species $`j`$ on species $`i`$ is given by $`J_{ij}`$; in general, $`J_{ij}J_{ji}`$. The only restriction we impose initially on the elements of $`𝐉`$ is that $`J_{ii}=0`$, i.e., none of the organisms preys on any other member of the same species. Since in all practical situations, food webs specify only the sign of $`J_{ij}`$ we allow the off-diagonal elements of $`J_{ij}`$ to take only the values $`+1`$ and $`1`$. Thus, $`J_{ij}=1=J_{ji}`$ indicates that $`j`$ is the prey and $`i`$ is the predator. Similarly, the situations $`J_{ij}=1=J_{ji}`$ and $`J_{ij}=1=J_{ji}`$ correspond, respectively, to competition and cooperation between the species $`i`$ and $`j`$. We assign the values $`+1`$ or $`1`$ to the off-diagonal elemts of $`𝐉`$ randomly with equal probability in the initial state of the eco-system sole . However, our model can be easily generalized to take into account any other architecture of food webs amaral ; abramson .

We now argue that the matrix $`𝐉`$ accounts not only for the inter-species interactions but also intra-species interactions. First of all, note that if $`J_{ij}>0`$, then the species $`j`$ is a prey of the species $`i`$ if, simultaneously, $`J_{ji}<0`$ whereas the species $`j`$ cooperates with $`i`$ if, simultaneously, $`J_{ij}>0`$ and $`J_{ji}>0`$. Therefore, if $`J_{ij}>0`$, the quantity $`(J_{ij}J_{ji})/2`$ is unity if the species $`j`$ is a prey of the species $`i`$, but it vanishes if the species $`i`$ and $`j`$ mutually cooperate. Similarly, if $`J_{ij}<0`$, the quantity $`(J_{ij}J_{ji})/2`$ is unity if the species $`j`$ is a predator of $`i`$, but it vanishes if the species $`i`$ and $`j`$ compete against each other. Now, consider the two sums

$$𝒮_i^\pm =\pm \underset{j=1}{\overset{N}{}}\frac{(J_{ij}^\pm J_{ji})}{2}n_j$$
(1)
where the superscript $`\pm `$ on $`J_{ij}`$ indicates that the sum is restricted to only the positive (negative) elements $`J_{ij}`$. The sum $`𝒮_i^+`$ is a measure of the total food currently available to the $`i`$-th species whereas $`𝒮_i^{}`$ is a measure of the total population of the $`i`$-th species that would be, at the same time, consumed as food by its predators. If the food available is less than the requirement, then some organisms of the species $`i`$ will die of starvation, even if none of them is killed by any predator. This way the matrix $`𝐉`$ can account for the shortfall in the food supply and the consequent competition among the organisms of the species $`i`$ .

The collective characteristics of species: The age of an arbitrary individual organism, say, $`\alpha `$ of the species $`i`$ at time $`t`$ is denoted by the symbol $`X(i,\alpha ;t)`$. In our model each species $`i`$ is collectively characterized by stauffer ; sole :
(i) the minimum reproduction age $`X_{rep}(i)`$, (ii) the birth rate $`M(i)`$, (iii) the maximum possible age $`X_{max}(i)`$, and (iv) the elements $`J_{ij}`$ and $`J_{ji}`$ $`(j=1,2,\mathrm{},N)`$. An individual of the $`i`$-th species can reproduce only after attaining the age $`X_{rep}(i)`$. Whenever an organism of $`i`$-th species gives birth to offsprings, $`M(i)`$ of these are born simultaneously. None of the individuals of the $`i`$-th species can live longer than $`X_{max}(i)`$. Thus, even if an individual manages to escape its predators, it cannot live longer than $`X_{max}(i)`$ because of “natural death” caused by ageing.

The dynamics of the eco-system: The state of the system is updated in discrete time steps as follows:

Step I- Birth: Assuming, for the sake of simplicity, the reproductions to be asexual, each individual organism $`\alpha `$ ($`\alpha =1,\mathrm{},n_i(t)`$) of the species $`i`$ ($`i=1,2,\mathrm{}N(t)`$) is allowed to give birth to $`M(i;t)`$ offsprings at every time step $`t`$ with probability (per unit time) $`p_b(i,\alpha ;t)`$ which is non-zero only when the individual organism age $`X(i,\alpha ;t)>X_{rep}(i;t)`$.

Step II- Natural death: At any arbitrary time step $`t`$ the probability (per unit time) of “natural” death (due to ageing) of an individual organism $`\alpha `$ of species $`i`$ is $`p_d(i,\alpha ;t)`$.

Step III- Mutation: With probability $`p_{mut}`$ per unit time, all the species simultaneously re-adjust one of the incoming interactions $`J_{ij}`$ by assigning it a new value of either $`+1`$ or $`1`$ with equal probability sole .

Step IV- Starvation death and killing by prey: If $`n_i𝒮^+`$ is larger than $`𝒮^{}`$ then food shortage will be the dominant cause of premature death of a fraction of the existing population of the species $`i`$. On the other hand, if $`𝒮^{}>n_i𝒮^+`$, then a fraction of the existing population will be wiped out primarily by the predators. In order to capture these phenomena, at every time step $`t`$, in addition to the natural death due to ageing, a further reduction of the population by

$$C\mathrm{max}(𝒮^{},n_i𝒮^+)$$
(2)
is implemented where $`n_i(t)`$ is the population of the species $`i`$ that survives after the natural death step above. $`C`$ is a constant of proportionality that basically sets the time scale of the population dynamics. If implementation of these steps makes $`n_i0`$, species $`i`$ becomes extinct.

Step V- Speciation: At each time step, the niches (nodes) left empty by extinction of species are re-filled by new species, with probability $`p_{sp}`$. All the simultaneously re-filled nodes of the network originate from one common ancestor which is picked up randomly from among the surviving species. All the interactions $`J_{ij}`$ and $`J_{ji}`$ of the new species are identical to those of their common ancestor; each new species, however, either competes or cooperates with its ancestor species. The characteristic parameters $`X_{max}`$, $`X_{rep}`$, $`M_j`$ of each of the new species differ randomly by $`\pm 1`$ from the corresponding parameters for their ancestor.

Probability of birth: We chose the time-dependent probabilty of birth per unit time as

$`p_b(i,\alpha )=\left[{\displaystyle \frac{X_{max}(i)X(i,\alpha )}{X_{max}(i)X_{rep}(i)}}\right]\left(1{\displaystyle \frac{n_i}{n_{max}}}\right)`$
$`\mathrm{iff}X(i,\alpha )X_{rep}(i)\mathrm{and}X_{max}(i)>X_{rep}(i)`$ (3)
Note that in the limit of vanishingly small population, i.e., $`n_i0`$, we have $`p_b(i,\alpha )1`$ if $`X(i,\alpha )=X_{rep}(i)`$ and, then, $`p_b`$ decreases linearly austad  as the organism grows older. However, since the eco-system can support only a maximum of $`n_{max}`$ individual organisms of each species, $`p_b(i,\alpha ;t)0`$ as $`n_i(t)n_{max}`$, irrespective of the age of the individual organism $`\alpha `$ cebrat .

Probability of natural death: We assume the probability of “natural” death (due to ageing) to have the form carey

$`p_d(i,\alpha )=\left[{\displaystyle \frac{X(i,\alpha )M(i)X_{rep}(i)}{X_{max}(i)M(i)X_{rep}(i)}}\right]`$
$`\mathrm{if}X(i,\alpha )X_{rep}(i)`$ (4)
$`p_d(i,\alpha )=\left[{\displaystyle \frac{X_{rep}(i)M(i)X_{rep}(i)}{X_{max}(i)M(i)X_{rep}(i)}}\right]`$
$`\mathrm{if}X(i,\alpha )<X_{rep}(i)`$ (5)
provided $`X_{max}(i)M(i)>X_{rep}(i)`$. In all other situations, $`p_d(i,\alpha )=1`$. Note that, for a given $`X_{max}`$ and $`X_{rep}`$, the larger is the $`M`$ the higher is the $`p_d`$ for any age $`X`$. Therefore, each species have a tendency to increase $`M`$ for giving birth to larger number of offsprings whereas the higher mortality for higher $`M`$ opposes this tendency.

The longest runs in our computer simulations were continued upto a maximum of five million time steps. If each time step in our model is assumed to correspond to a real time of the order of one year, then the time scale of $`5`$ million years, over which we have monitored our model eco-system, is comparable to real geological time scales.

Since we faced difficulty in getting high quality data, with reasonably good statistics, for $`N_{max}>100`$ and $`n_{max}>1000`$, we have carried out most of our simulations with $`N_{max}=50,100`$ and $`n_{max}=100,1000`$ only. The data obtained from the different runs, each starting from a random initial condition, were averaged. Both CRAY-T3E and SUN workstations were used for the simulations.

In fig.1 we plot the total number of species, $`N(t)`$, in a particular run, starting from a single initial condition, upto half a million time steps. In fig.2 we plot the corresponding variation of the total population $`n(t)`$ over relatively short interval of $`20,000`$ time steps only. These clearly demonstrate that the evolution has periods of “stasis” during which organisms populations keep fluctuating; the stasis are interrupted by occasional bursts of rapid extinctions followed by slower recovery.

The average distributions of the lifetimes of the species are plotted in fig.3 for one set of values of the parameters. Clearly, the data are consistent with a power-law; the effective exponent, which is, approximately, $`2`$, is also consistent with the corresponding estimate quoted in the literature drossel ; newman . However, in fig.3 the power law holds only over a limited range chu . Since real eco-systems are much more complex than our model eco-system and the available fossil data are quite sparse, it is questionable whether real extinctions follow power laws and, if so, over how many orders of magnitude.

We have also observed (not shown) that the higher is the mutation probability $`p_{mut}`$ the lower is the lifetime; this is consistent with the intuitive expectation that the higher rate of mutation leads to higher levels of biological activity in the eco-system thereby leading to the extinction of larger number of species. But, $`p_{sp}`$ had weaker effect on the same data. However, if $`p_{sp}`$ is too small to maintain adequate pace of recovery of the eco-system after mass extinctions, the entire eco-system collapses.

Fig.4 shows the time-averaged age-distribution in the populations of a species as well as the distributions of $`X_{max}`$, $`X_{rep}`$ and $`M`$. We see that the minimum age of reproduction $`X_{rep}`$ is quite small, as usual in the employed ageing model stauffer . The age distribution decays stronger than a simple exponential, indicating a mortality increasing with age as it should be carey . The genetic death ages $`X_{max}72`$ are far above the upper end $`31`$ of the age distribution, as is appropriate for animals in the wild austad . Finally, fig.5 shows the distribution of $`M(i)`$; this is relatively much broader than the distributions of $`X_{max}`$ and $`X_{rep}`$.

In this letter we have presented a unified model which describes the birth, ageing and death of individuals and population dynamics on short time scales as well as the long-time evolution of species. Not only the total number of species and the inter-species interactions but also the collective characteristics, namely, $`X_{rep},X_{max}`$ and $`M`$, of each species vary following a stochastic dynamics. Thus, our model is capable of self-organization. To our knowledge, there are only a few earlier evolutionary models abramson  based on inter-species interactions which describe population dynamics of each species. The population dynamics within the framework of Lotka-Volterra equations have been considered earlier jan  for only a few species. But, these do not account for the age distributions as the entire population of each species is represented collectively by a single dynamical variable in contrast to the explicit birth, growth and death of individual organisms captured in our model.

Since we have observed strong deviations from power-law in the distributions of lifetimes of species even in the Sole-Manrubia model chowstau , in spite of the absence of detailed ”micro”-dynamics in the latter, we strongly believe that this is a generic features of evolution and extinction of species. It would be interesting to investigate the geographical effects on our model eco-system by re-formulating it on a lattice in the same spirit in which some lattice models of prey-predator systems have been formulated tome .

We thank P.M.C. de Oliveira for comments on the manuscript and the Supercomputer Center Jülich for computer time on their CRAY-T3E. This work is supported by Deutsche Forschungsgemeinschaft through a Indo-German joint research project and by Alexander von Humboldt Foundation.


text:
# 1 Introduction: radio emission from Low/Hard state Black Holes

## 1 Introduction: radio emission from Low/Hard state Black Holes

Radio emission is often observed from X-ray binaries, particularly transient systems and especially Black Hole (BH) candidates (for detailed reviews see e.g. , , , ). The Low/Hard state is one of the five ‘canonical’ X-ray states observed from BH X-ray binaries in our Galaxy. It is charaterized by a hard (spectral index $``$ 1.5) X-ray spectrum, associated with a steady, self-absorbed jet which emits synchrotron radiation; the Off state may simply be the Low state ‘turned down’ to lower accretion rates. X-ray spectra of High/Soft state BHs are dominated by thermal radiation, while the radio emission drops below detectable levels, probably corresponding to the physical disappearance of the jet. For the so called Intermediate and Very High states the connection between X-ray and radio properties is not yet completely established. Transitions between states are often associated with multiple ejections of synchrotron emitting material, possibly with high Lorentz factors.
There is an extremely strong correlation between radio and X-ray emission: the former has been directly observed to arise in outflows and to produce synchtrotron emission from a population of high energy electrons. The latter has been inferred to arise via Comptonisation by a thermal electron distribution. All the evidence points to the corona in these systems being physically related to the presence of a jet: by far the simplest interpretation therefore is that the Comptonising region is just the base of the relativistic outflow.

## 2 Jet/corona coupling: $`S_{radio}`$ versus $`S_X`$

By means of simultaneous radio/X-ray obsevations performed on the Galactic BH GX 339-4, Corbel et al. (in preparation) have found an interesting correlation between fluxes, namely the radio density flux at 8.6 GHz scales as the 3-9 keV X-ray flux raised to 0.7.
Assuming the same physics and jet/corona coupling for all sources, based on simple arguments one would predict the following scenario in the $`S_{radio}`$ vs. $`S_X`$ parameter space (see sketch in fig. 1 for clarity):

1. All sources lying on a line with the same slope if neither X-ray or radio emission were strongly beamed (i.e. low $`\beta `$).
2. Different sources lying on lines with the same slope but different normalisations if X-rays were isotropic but radio was beamed, with radio-brighter sources corresponding to higher Doppler factors.
3. As point 2. but with higher Doppler factors sources being brighter in both radio and X-ray if either were beamed.

With the purpose of verifying if and possibly which kind of relation could be shared among a wide sample of sources in the same spectral state, we collected a variety of quasi-simultaneous radio/X-ray observations of Low state BHs and plotted the radio vs. the X-ray fluxes (below 12 keV) scaled to a distance of 1 kpc (units are mJy vs. Crab). As shown in 1, there is actually clear evidence for a strong correlation, which - again - extends over more than three orders of magnitude in accretion rate, from around $`10^5\dot{m}_{Edd}`$ (calculated supposing $`M_{BH}=`$10 $`M_{}`$ ) up to close to the High/Soft state transition. It is important to stress that, in principle, the range of normalisation values can be extremely large. Assuming again the same physics for all sources, and that X-rays are isotropic while radio is beamed, the ratio $`S_{radio}/S_X`$ for two sources would equal the ratio between their Doppler factors (calculated as sum of approaching and receding jet) to the second power. For instance, with $`\beta `$=0.998 (i.e. $`\mathrm{\Gamma }`$=16) and angles ranging between 1 \- 90 , that ratio covers 12 orders of magnitude!

## 3 Special cases: V 404 Cyg and GX 339-4

Remarkably, we found that the data of V 404 Cyg and GX 339-4, i.e. the sources for which we have the largest range in luminosity at our disposal, are well fitted by a power law with the same slope, as shown in fig. 3: we obtained $`S_{radio}=k\times S_X^{+0.7\pm 0.1}`$, with $`k_{V404}=254\pm 18`$ and $`k_{GX339}=124\pm 2`$. For the same simple arguments we mentioned in the previous section, assuming isotropic X-rays and beamed radio emission, the ratio $`k_{V404}/k_{GX33}`$ should be equal to $`(\delta _{V404}/\delta _{GX339})^2`$. GX 339-4 seems to be at quite low inclination ( 15 \- 30 according with  , ) while $`i_{V404}`$=56 (, ); adopting these values the ratio $`(\delta _{V404}/\delta _{GX339})^2<1`$ for *any* value of $`\beta `$=v/c. Thus, V 404 being radio-brighter than GX 339-4 despite the higher inclination excludes both the second and the third scenario described in the previous section. The simplest possible explanation requires very low beaming for both sources, attributing the little scatter in normalisation to intrinsic errors. Obviously a great uncertainty in this argument is introduced by distance determinations, which strongly affect all the calculations. In addition we can not rule out the possibility of a much more complex physical model behind this relation.

## 4 Summary and conclusions

The results of our analysis can be summarized as follows:

* In Low/Hard state BHs the observed radio and X-ray fluxes are correlated over more than three orders of magnitude in accretion rate, implying a strong jet/corona coupling; no lower limit to the relation has been found.
* We can confidently assert that, even at accretion rates as low as $``$ $`10^5\dot{m}_{Edd}`$ a powerful jet is being formed.
* V 404 Cyg is the second source to display $`S_{radio}S_X^{+0.7}`$, from quiescence up to Soft state transition. A physical explanation for this relation is proposed by Markoff et al. ().
* Above $`10^2\dot{m}_{Edd}`$ the jet disappears within a factor of a few, probably in all sources (observed in three sources and no exceptions)
* Comparison of different sources may indicate that jets in Low/Hard have a low velocity compared to those in transient outbursts. Of course better distance and inclination determinations are strongly required to probe this conjecture.
Finally, we would like to mention that the possible detection of the Soft X-ray Transient A0620-00 at the predicted – extrapolating the relation we found – radio level would demonstrate that jet production is an ubiquitous mechanism between $`10^7`$ and $`10^2\dot{m}_{Edd}`$.


text:
# 1 Introduction

## 1 Introduction

The lack of exact radiative solutions to Einstein’s equations has promoted perturbation theory in General Relativity into a special place, as the tool for analyzing gravitational radiation emitted in physically interesting events. The use of perturbation theory in spacetimes containing black holes started with the work of Regge and Wheeler , where they addressed the stability of the Schwarzschild geometry against small deviations. This analysis was extended to include the infall of particles in a Schwarzschild black hole by Zerilli  and others (see, e.g., ).

The Kerr geometry, without spherical symmetry, proved to be more difficult to handle, but has also given some insights into relations that previously seemed a mystery . Teukolsky  was able to decouple and separate the perturbation equations for the Newman-Penrose quantities, and reduced them to a single radial equation, now known as the Teukolsky equation. However, it suffers from a couple of drawbacks, in that it has a long range potential and a source term which is in general divergent at large distances. This makes it difficult to analyze radiation emitted by test particles in generic orbits, though special cases have been handled (see for example the case of circular orbits ). Sasaki and Nakamura  improved this situation by introducing a set of transformations that bring the Teukolsky equation into the Sasaki-Nakamura equation, which has both a short range potential and an everywhere well behaved source term. Moreover, their equation reduces to the Regge-Wheeler equation when the rotation parameter is set equal to zero. Using their formalism, Sasaki, Nakamura and co-workers  have computed the gravitational radiation for generic orbits of particles falling, from rest at infinity, into a Kerr black hole. These studies have recently been revisited  not only to study the effects of radiation reaction, but also to produce accurate waveforms, to serve as templates to the various gravitational wave detectors already at work. An extension of the Sasaki-Nakamura formalism to perturbations other than gravitational is given in .

It is important to stress that this approach is only justified as long as $`\mu `$, the mass of the infalling particle, is much smaller than $`M`$, the mass of the black hole, so that it can serve as perturbation parameter. By the late 70’s however, when the first full numerical simulations of black hole collisions were being done , it became clear that taking the limit $`\mu M`$ gave unexpectedly good results in perturbation theory. For recent improvements see .

The high velocity collision of two black holes, a problem interesting in itself, has gained renewed interest recently, with the possibility of black hole formation at TeV scales . This process was studied extensively some time ago by D’Eath and Payne , by performing a perturbation expansion around the Aichelburg-Sexl metric, which describes a Schwarzschild black hole moving at the speed of light. The question arises as to whether it is possible to treat this process in the old Regge-Wheeler-Zerilli approach. Is it possible to study the collision of two black holes moving at the speed of light, by considering some expansion around a Schwarzschild metric? In a previous paper  we have argued that it is possible, provided the agreement between perturbation theory and numerical relativity is more than a coincidence. In spite of this unproven universal agreement, we do have the strong conviction that it will hold in this process, but we are ultimately justified by the excellent agreement between our results  and previous studies . The essence of our previous study was simple: consider a highly relativistic particle falling into a Schwarzschild black hole, use the Zerilli approach to compute the waveform and energy, Lorentz-boost the black hole to high velocities in the direction of the infalling particle, and end up with a collision at high velocity between a black hole and a small test particle. If one employs the assumption that this still works if $`\mu M`$ the result follows. Our results  showed an excellent agreement with results by D’Eath and Payne  and with results by Smarr , which leads us to believe that once again perturbation theory gives very good results throughout all the values of the perturbation parameter.

One would now like to extend these results to the Kerr family of rotating black holes, a less studied geometry. D’Eath and Payne’s results do not apply here for example, although with probably a great effort one could use their methods to study the collision of Aichelburg-Sexl-Kerr particles, which solution has been found in . It is true that we do not expect to see black holes colliding at the speed of light in any astrophysical scenario, but we do expect to see this at Planck energies, or, perhaps at the LHC, should the TeV scenario  be correct. In any case, the holes will most probably be rotating. Furthermore, the work done so far in the Kerr geometry dealt only with infall which starts from rest, so it will be interesting to see the outcome when the particle has a non zero velocity at infinity, and compare it with results in the Schwarzschild geometry.

Thus, the situation we consider in this paper is the following: a highly relativistic point particle impinges radially into a Kerr black hole, along its symmetry axis. We will use the Sasaki-Nakamura formalism to find the energy radiated, which will be the main result of this paper. We will then Lorentz-boost the Kerr hole to high velocity in the direction of the infalling object, which basically amounts to put $`M\gamma M`$, where $`\gamma `$ is the Lorentz factor. This will then describe the high velocity collision between a Kerr hole and a small particle. Finally, we put $`\mu M`$, which means that we are dealing with the collision at nearly the speed of light of two equal mass black holes, one rotating, the other non-rotating. The justification for doing this last step (not allowed on formal basis) comes from the excellent results obtained so far by perturbation theory .

## 2 A Green’s function solution to the Sasaki-Nakamura equation

After some manipulations, the Teukolsky equation  may be brought to the Sasaki-Nakamura  form (details about the Teukolsky formalism may be found in the original literature , and also in . For a good account of the Sasaki-Nakamura formalism we refer the reader to ,, and ):

$$\frac{d^2}{dr_{}^2}X(\omega ,r)\frac{d}{dr_{}}X(\omega ,r)𝒰X(\omega ,r)=.$$
(1)
The tortoise $`r_{}`$ coordinate is defined by $`dr_{}/dr=(r^2+a^2)/\mathrm{\Delta }`$, and ranges from $`\mathrm{}`$ at the horizon to $`+\mathrm{}`$ at spatial infinity. The functions $``$ and $`𝒰`$ can be found in the original literature . We are considering the radial infall of a highly relativistic particle into a Kerr black hole along the symmetry axis, so the situation is axisymmetric. This means that the azimuthal quantum number $`m`$ appearing in $``$ and $`𝒰`$ , may be set to zero. Also, in this simple situation it is possible to find an exact form for the function $``$ (compare this with the source term in Schwarzschild ):

$$=\frac{\mu Cϵ_0\gamma _0\mathrm{\Delta }}{2\omega ^2r^2(r^2+a^2)^{3/2}}e^{i\omega r_{}}.$$
(2)
Here, $`C=\left[\frac{8Z_l^{a\omega }}{\mathrm{sin}^2\theta }\right]_{\theta =0}`$ and $`Z_l^{a\omega }(\theta )`$ is a spin-weighted spheroidal harmonic . The function $`\gamma _0=\gamma _0(r)`$ can be found in . The parameter $`ϵ_0`$ is the energy per unit rest mass of the infalling particle, which we take to be a very large quantity ($`ϵ_0\mathrm{}`$), since we are interested in highly relativistic particles. The Sasaki-Nakamura equation (1) is to be solved under the “only outgoing radiation at infinity” boundary condition, meaning

$$X(\omega ,r)=X^{\mathrm{out}}e^{i\omega r_{}},r_{}\mathrm{}.$$
(3)
Once $`X^{\mathrm{out}}`$ is known, Teukolsky’s radial function $`R`$ can be found, when $`r_{}\mathrm{}`$ (the region of interest here) as

$$R=\frac{4\omega ^2X^{\mathrm{out}}}{\lambda (\lambda +2)12i\omega 12a^2\omega ^2}r^3e^{i\omega r_{}}=R^{\mathrm{out}}r^3e^{i\omega r_{}}.$$
(4)
Following Nakamura and Sasaki  we define the multipolar structure so as to have

$$\mathrm{\Delta }E=\frac{\pi }{4}_0^{\mathrm{}}\underset{lm}{}|h^{lm}(\omega )|^2+|h^{lm}(\omega )|^2.$$
(5)
Once $`X`$ is known, we can get the Teukolsky wavefunction near infinity from (4), and the energy radiated away from (5).

We would now like to find $`X(\omega ,r)`$ from the Sasaki-Nakamura differential equation (1). This is accomplished by a Green’s function technique, constructed so as to satisfy the usual boundary conditions, i.e., only ingoing waves at the horizon ($`Xe^{i\omega r_{}},r_{}\mathrm{}`$) and outgoing waves at infinity($`Xe^{i\omega r_{}},r_{}\mathrm{}`$). We get that, near infinity (we are interested in knowing the wavefunction in this region),

$$X^{\mathrm{out}}=\frac{\mu ϵ_0c_0C}{4i\omega ^3B}\frac{e^{i\omega r_{}}X^H}{r^2(r^2+a^2)^{1/2}}𝑑r.$$
(6)
Here $`c_0\gamma _0(r=\mathrm{})`$, and $`X^H`$ is an homogeneous solution of (1) which asymptotically behaves as

$`X^HA(\omega )e^{i\omega r_{}}+B(\omega )e^{i\omega r_{}},r_{}\mathrm{}`$ (7)
$`X^He^{i\omega r_{}},r_{}\mathrm{}.`$ (8)
We now discuss the numerical results.

## 3 Numerical Results and Conclusions

The main result of this paper is shown in Fig. 1, the energy spectra as a function of the angular quantum number $`l`$, for an “extreme” black hole with $`a=0.999M`$. The results for other rotation parameters $`a`$ are not shown, both because there are no qualitatives changes and mostly because they follow from the discussion below. The most interesting and important features of these results are:

(i) Fig. 1 clearly shows that the zero frequency limit (ZFL) of the energy spectra, $`\frac{dE}{d\omega }_{\omega 0}`$ is non-vanishing. We know  that the existence of a ZFL is closely linked to a non-zero velocity at infinity, so this comes as no surprise and was already observed  in the collision of a Schwarzchild black hole with an highly relativistic particle. The peculiar property of the ZFL arises when we look more closely into the numerical results: the numerical data shows that the ZFL is exactly (up to the numerical error) equal to the ZFL for the infall of an highly relativistic particle into a Schwarzschild black hole. This means that the ZFL seems to be independent of the spin of the colliding objects, and given by Smarr’s  expression.

(ii) The spectra is almost flat up to a critical ($`l`$-dependent) frequency, a fact also evident from Fig. 1. Again, this was also true for a non-rotating black hole . The spectra is flat up to $`\omega \omega _{QN}`$, where $`\omega _{QN}`$ is the lowest quasinormal frequency for the Kerr black hole (for work on quasinormal modes on the Kerr geometry see , for example). For $`\omega >\omega _{QN}`$, the spectra decays exponentially, according to the empirical law $`\frac{dE}{d\omega }e^{bl}`$, ($`\omega >\omega _{QN}`$), first discovered by Davis et al , and further discussed by Sasaki and Nakamura .

Figure 1. The energy spectra for a point particle moving at nearly the speed of light and colliding, along the symmetry axis, with an extreme (a=0.999M) Kerr black hole. Notice that the spectra is almost flat, the ZFL is non-vanishing and that the quadrupole carries less than half of the total radiated energy.

(iii) The total energy radiated in each multipole goes as $`\mathrm{\Delta }E_l\frac{1}{l^2}`$, which (needless to say) is exactly what happens for Schwarzschild case . This power-law dependence seems to be universal for highly relativistic collisions, and was first observed for the non-rotating case . We note that for collisions beginning at rest, the behaviour is strikingly different, for the energy goes as $`\mathrm{\Delta }E_l(\mathrm{rest}\mathrm{at}\mathrm{})e^{bl}`$ .

(iv) This power-law dependence implies that the radiation is not quadrupole ($`l=2`$) in nature, as it was for infall starting from rest. In fact, less than $`50\%`$ of the total energy is carried in the quadrupole part.

These characteristics plus the discussion in  allow one to infer that the total energy radiated will depend mainly on the behaviour of the quasinormal frequencies for the Kerr geometry. Now, for $`m=0`$, we know  that the quasinormal frequencies for this geometry are almost the same (but slightly larger) as for the Schwarzschild geometry, so we expect the total energy to be similar in both cases (slightly larger in the Kerr case). In fact, for this axisymmetric collision we find, on summing over all $`l`$’s that

$$\mathrm{\Delta }E=0.31\frac{\mu ^2ϵ_0^2}{M}a=0.999M,$$
(9)
with a $`5\%`$ error. We recall that for a Schwarzschild black hole  $`\mathrm{\Delta }E=0.26\frac{\mu ^2ϵ_0^2}{M}`$ so the infall along the symmetry axis of a Kerr hole does not enhance very much the total radiated energy (in comparison with the non-rotating case). This was also observed  for collisions along the symmetry axis, but starting from rest. As we lower the rotation parameter, the total energy decreases, and approaches the non-rotating value given above. What can we say about the collision at nearly the speed of light between a Schwarzschild and a Kerr black hole, along its symmetry axis? Supposing (with all due precautions mentioned in the Introduction) that (9) holds for $`\mu ϵ_0M`$, we should have an efficiency of $`15,5\%`$ for that process. This seems a reasonable result, but only a full numerical scheme can tell how accurate it is. For the moment, one can only say that it is a “good” result: it is higher than for the collision of two Schwarzschild black holes, but still within the upper limit imposed by the Area Theorem $`\mathrm{\Delta }E0.66M`$ and an eficiency less than $`33\%`$ . These numbers arise noting that the situation is axisymmetric, the radiation carries no angular momentum, and thus one has the same $`a=0.999M`$ for the final black hole.

## Acknowledgements

This work was partially funded by Fundação para a Ciência e Tecnologia (FCT) through project PESO/PRO/2000/4014. V.C. also acknowledges finantial support from FCT through PRAXIS XXI programme. J. P. S. L. thanks Observatório Nacional do Rio de Janeiro for hospitality.