text:
# Gamma-rays and cosmic-rays from a pulsar in Cygnus OB2

## 1 Introduction

The Cygnus OB2 is the largest massive galactic association at a relatively small distance of 1.7 kpc. It has the diameter of $`60`$ pc and the core radius of $`14\pm 2`$ pc. The total H2 mass of the Cyg OB2 derived from the CO survey is $`3.3\times 10^5`$ M (Butt et al. 2003), and the total stellar mass is $`(410)\times 10^4`$ M (Knödlseder 2000). It contains many massive stars of the OB type ($`2600\pm 400`$) and O type ($`120\pm 20`$) (Knödlseder 2000). The typical main sequence evolution lifetime of the massive O type stars is of the order of a few $`10^6`$ yrs and B type stars a few tens of $`10^6`$ yrs. Massive O type stars should pass through the Wolf-Rayet phase and explode as supernovae. It is expected that very fast pulsars are created in explosions of such massive stars at a rate of about one per a few $`10^4`$ yrs within the Cyg OB2. This rate can be much higher if most of these massive stars are of similar age. The pulsars, created $`10^4`$ yrs ago, should resemble the Vela type pulsars (e.g. Vela pulsar, PSR 1706-44), which are strong sources of pulsed $`\gamma `$-rays in the hundred MeV-GeV energies (EGRET observations, Kanbach et al. 1994, Thompson et al. 1996). These pulsars are also surrounded by relatively small pulsar wind nebulae (PWNe), which are sources of TeV $`\gamma `$-rays (Chadwick et al. 2000, Kifune et al. 1995, Chadwick et al. 1998).

The Cygnus region is difficult to observe in high energies due to its very complex structure and location on the galactic plane. The analysis of the EGRET data from different viewing periods revealed a few sources in a relatively small field of view containing the Cyg OB2: 2EG J2033+4112 (Mori et al. 1997), GeV J2035+4214 (Lamb & Macomb 1997), GRO J2034+4203 (Reimer, Dingus & Nolan 1997), 3EG J2033+4118 (Hartman et al. 1999). The first source is interpreted by the authors as caused by the massive binary Cyg X-3. The other sources may be due to the same object in the Cyg OB2, according to Butt et al. (2003).

Recently, the HEGRA group has reported the discovery of an unidentified, steady, TeV source at the edge of the 95% error circle of the source 3EG J2033+4118 (Aharonian et al. 2002). The spectrum of this source is rather flat (differential spectral index $`1.9\pm 0.3_{\mathrm{stat}}\pm 0.3_{\mathrm{sys}}`$), with the flux about two orders below the flux observed from 3EG J2033+4118. Moreover, there is an indication that the TeV source is extended with a radius of $`5.6^{}`$ (Aharonian et al. 2002), confirmed in the new analysis of the HEGRA data by Rowell & Horns (2002).

The discovery of TeV source has led to a more careful analysis of the X-ray observations of Cyg OB 2 from the ROSAT and Chandra telescopes (Mukherjee et al. 2003, Butt et al. 2003). None of the several X-ray sources detected by the Chandra in the region with $`11^{}`$ diameter around the TeV source was particularly prominent (Butt et al. 2003). Mukherjee et al. (2003) notes only one interesting variable source which lies $`7^{}`$ from the TeV centroid with the flux of $`2\times 10^{12}`$ ergs cm-2 s-1. It might be a quiescent low-mass X-ray binary containing a neutron star or a black hole. However, weak diffuse X-ray emission has also been observed, with the intensity of $`1.3\times 10^{12}`$ ergs cm-2 s-1 in 0.5-2.5 keV and $`3.6\times 10^{12}`$ ergs cm-2 s-1 in 2.5-10 keV (Butt et al. 2003). This is more than ten times as bright as the sum of all the point-like sources.

Recent analysis of the arrival directions of cosmic rays with energies of $`10^{18}`$ eV by the AGASA group shows the excess of particles from the direction of the Galactic Center (Hayashida et al. 1999a, confirmed by the SUGAR analysis, Bellido et al. 2001), and from the direction of the Cygnus region. The second excess has a significance of 3.9 $`\sigma `$ (3401 events with the background of 3148) in the energy range $`10^{18}10^{18.4}`$ eV (Hayashida et al. 1999b). We suggest that this excess of EHE particles can also be related to the high energy processes occurring in the massive Cyg OB2 association.

The possibility that $`\gamma `$-rays are produced by massive young stellar clusters, as a result of acceleration of particles in the winds or shocks produced by the early type stars, has been frequently discussed in the past (with the application to the Cyg OB2, see e.g. Chen, White & Bertsch 1992, Manchanda et al. 1996, Benaglia et al. 2001). Aharonian et al. (2002) discusses such general scenario, suggesting that the TeV emission may originate in collisions of hadrons with a dense cloud. These authors also propose that $`\gamma `$-ray emission may originate in the inverse Compton scattering by electrons, which are accelerated by a jet-driven termination shock expelled from a microquasar, e.g. Cyg X-3 (see Aharonian & Atoyan 1998). Butt et al. (2003) argue that the TeV emission is more likely to originate in collisions of hadrons with the matter than in the leptonic processes. Hadrons might be accelerated to PeV energies in collisions of strong stellar winds.

In this paper we interpret the variety of high energy processes connected with the Cyg OB2 association by applying the model in which the basic role is played by a young pulsar, born in this association about $`10^4`$ yrs ago. It is shown that the observed GeV emission can be due to the Vela type pulsar present inside or close to the Cyg OB2 (Sect. 3). The TeV $`\gamma `$-ray emission originates in the pulsar wind nebula as a result of interactions of hadrons and leptons (Sect. 4). The observed excess of cosmic rays from the direction of the Cygnus region at energies of $`10^{18}`$ eV is due to neutrons, which are dissolved from heavy nuclei in their interactions with the matter of the Cyg OB2 (Sect. 5). We show that if a small part of nuclei, which escaped from the PWNa, is captured in a strong magnetic field of dense regions of the Cyg OB2 association, the flux of neutrons can explain the observed cosmic ray excess from the Cygnus region.

## 2 A model for high energy processes in the Cyg OB2

Let us assume that a very energetic pulsar is formed in a core collapse of one of the massive stars in the Cyg OB2 association. The massive O type stars should evolve through the Wolf-Raeyt phase during which the star loses a significant part of its envelope, rotates fast, and has strong surface magnetic fields. The Wolf Raeyt stars are probably progenitors of neutron stars with extreme parameters (period, surface magnetic field) which allow acceleration of particles to very high energies. Following the recent works (e.g. Gallant & Arons 1994, Blasi, Epstein & Olinto 2000) we assume that such neutron stars can accelerate heavy nuclei in their wind zones (e.g. in the mechanism called magnetic slingshot, Gunn & Ostriker 1969). Arons and collaborators (see e.g. Arons 1998) argue that the iron nuclei can reach energies equal to a significant part, $`\chi `$, of the total potential drop through the pulsar polar cap region,

$`E_{\mathrm{Fe}}=\chi Ze\mathrm{\Phi }_{\mathrm{open}},`$ (1)
where, $`Ze`$ is the charge of the iron nuclei, $`\mathrm{\Phi }_{\mathrm{open}}=\sqrt{E_{\mathrm{rot}}/c}`$ is the maximum potential drop across the polar cup of the pulsar, $`E_{\mathrm{rot}}`$ is the rate of rotational energy lost by the pulsar, $`\chi 1`$, and $`c`$ is the velocity of light. These nuclei are accelerated with efficiency corresponding to the Goldreich & Julian current at the light cylinder radius (Goldreich & Julian 1969). Due to the rotational energy losses of the pulsar, the energies of nuclei injected into the expanding nebula drop with time. Based on the modelling of the Crab Nebula, Arons (1998) argues that $`\chi `$ should be not very far from one. Assuming $`\chi =0.5`$, we use the following prescription for the evolution of the Lorentz factor of injected Fe nuclei with time

$`\gamma _{\mathrm{Fe}}=E_{\mathrm{Fe}}/m_{\mathrm{Fe}}4\times 10^9B_{12}P_{\mathrm{ms}}^2,`$ (2)
where $`m_{\mathrm{Fe}}`$ is the mass of the iron nuclei, $`B=10^{12}B_{12}`$ G is the surface magnetic field of the pulsar, and $`P=10^3P_{\mathrm{ms}}`$ s is the pulsar period which evolves in time $`t`$ (in seconds) due to the rotational energy losses according to

$`P_{\mathrm{ms}}^2=P_{0,\mathrm{ms}}^2+2\times 10^9tB_{12}^2,`$ (3)
where $`P_{0,\mathrm{ms}}`$ is an initial period of the pulsar.

In order to obtain the equilibrium spectra of these nuclei inside the nebula at a specific time, we construct a simple time dependent model for the expanding nebula, basing on the early papers which discuss the evolution of a supernova envelope under the influence of a young pulsar (e.g. Ostriker & Gunn 1971). This model takes into account possible acceleration of the nebula due to the energy supply by the pulsar. This has influence on the main parameters of the nebula, i.e. an expansion velocity, a shock radius, an outer radius, and a magnetic field inside the nebula. The details of such a model are described in Sect. 2 of Bednarek & Bartosik (2003). The equilibrium spectra of nuclei are obtained by integrating their injection spectra at a specific time over the activity period of the pulsar, and taking into account different processes such as: (1) inelastic collisions of nuclei with the matter (important at an early phase); (2) the adiabatic energy losses of particles due to expansion of the nebula; and (3) the escape of nuclei from the nebula. We include all these processes as it has been already done in our previous papers (see sect. 4 in Bednarek & Protheroe 2002). However, in order not to complicate the model too much, the gravitational energy losses of the pulsar are neglected. We also include partial disintegration of iron nuclei, extracted from the neutron star surface, during their propagation through the radiation field of the pulsar outer gaps in the inner pulsar magnetosphere following Bednarek & Protheroe (1997). Even for the pulsars with the age $`2\times 10^4`$ years, the iron nuclei still lose on average one nucleon in collisions with the nonthermal radiation in the pulsar’s outer gap. Nucleons (protons, neutrons) dissolved from the iron nuclei in this process are added to the spectrum of nuclei injected into the nebula. The example equilibrium spectrum of nuclei, which are still inside the nebula at $`2\times 10^4`$ years after supernova explosion, is shown in Fig. 1. A higher energy bump in the spectrum of hadrons contains heavy nuclei of the iron group. Their energies are limited to a relatively narrow range due to the adiabatic energy losses of nuclei injected at the earlier phase and the escape of more energetic nuclei from the nebula. A lower energy bump in the spectrum of hadrons contains protons from disintegrations of heavy nuclei during their propagation in the pulsar outer gap and from collisions of nuclei with a dense matter of the nebula at an early stage after supernova explosion. All the medium mass nuclei with relatively large Lorentz factors, produced in collisions of nuclei with the matter of the nebula at an early phase after the explosion, have already managed to escape from the nebula. The equilibrium spectra of nuclei inside the nebula at a specific time allow us to calculate the spectra of $`\gamma `$-rays from decay of pions, produced by nuclei in collisions with the matter of the nebula.

According to Hoshino et al. (1992), relativistic heavy nuclei can transfer a part of their energy, $`\xi `$, to relativistic positrons in the process of resonant scattering in the region above the pulsar wind shock. As a result, positrons are accelerated with the power law spectrum with the index $`\alpha 2`$ between $`m_\mathrm{e}c^2\gamma _{\mathrm{Fe}}`$ and $`\gamma _{\mathrm{Fe}}m_{\mathrm{Fe}}c^2/Z`$, where $`m_\mathrm{e}`$ is the mass of the electron. The injection spectra of these positrons evolve in time due to the change of energies and injection rate of iron nuclei with the pulsar period. The equilibrium spectra of positrons inside the nebula at a specific time are calculated by taking into account their radiation energy losses on bremsstrahlung, synchrotron, and inverse Compton (IC) processes, and due to the adiabatic losses caused by an expansion of the nebula as it was shown in Sect. 4.2 in Bednarek & Bartosik (2003). The efficiencies of specific energy loss processes are determined by the evolving in time parameters of the nebula. The bremsstrahlung process depends on the average density of matter inside the nebula which can be easily evaluated for the known total mass of the nebula and its radius at a specific time. The synchrotron losses of leptons are determined by the magnetic field strength inside the nebula which is estimated from the known dependence of the magnetic field inside the pulsar inner magnetosphere and the pulsar wind zone (see for details Eq. 7 in sect. 2 in Bednarek & Bartosik 2003). The evolution of the magnetic field strength in the pulsar wind postshock region is shown in Fig. 2b, for the example parameters of the pulsar and the nebula. The IC losses of leptons is more difficult to estimate. For relatively old nebulae the energy density of background radiation is dominated by the microwave background radiation (MBR), and possibly infrared background and optical background. Since the optical background created by the massive stars inside the Cyg OB2 becomes important only for low energy leptons, $`100`$ GeV, we do not consider their influence on the equilibrium spectrum of leptons inside the nebula. The radiation field inside relatively young nebulae is dominated by the synchrotron photons produced by leptons inside the nebula which spectrum is determined by the conditions inside the nebula in the past. Therefore, in order to estimate the density of synchrotron photons at a specific time we apply simple model for the evolution of nebula and include the above defined prescription for the injection of particles and the evolution of the radiation processes. It is assumed that only a part of the whole nebula radiate efficiently. This might be motivated by the observed anisotropic injection of particles into the pulsar wind nebulae which manifests itself by the presence of X-ray torus in the equatorial plane of the pulsars (e.g. the X-ray torus around $`20^o`$ of the equatorial plane observed in the Crab Nebula). We estimate the density of synchrotron photons inside the nebula assuming that leptons are injected and radiate only within the part of the nebula limited by the above mentioned $`20^o`$ torus. The evolution of the synchrotron background inside the nebula at a specific time after supernova explosion is shown in Fig. 2c. It is clear that inside the nebulae older than $`3\times 10^3`$ yrs, the IC energy losses of leptons on synchrotron photons becomes negligible. Only scattering of the MBR and possibly the infrared radiation, produced by the hot gas inside the nebula, plays important role. Fig. 2a shows the equilibrium spectra of leptons inside the nebula at a specific time after supernova explosion. The lower energy part of the spectrum of leptons is formed by leptons which were injected relatively early after supernova explosion and lose energy at latter time mainly due to the adiabatic expansion of the nebula. Therefore, the low energy cut-off in the spectrum of leptons, present between 0.1 - 0.01 GeV (Fig. 2a), shifts to lower energies with the evolution of the nebula. The efficient synchrotron cooling determines the shape of the spectrum in the higher energy range when the nebula is young and the magnetic field is strong. This process is responsible for the power low spectrum with the index close to 2 (seen in Fig. 2a) at the time $`10^2`$ yrs. On the other hand, the spectrum of leptons at higher energy part at the latter phase of expansion of the nebula is formed by freshly injected leptons which are not able to cool to low energies fast enough due to the weak magnetic and photon fields inside the nebula.

Described above pulsar-supernova model is applied to a variety of high energy processes observed from the direction of the Cyg OB2 association. To describe self-consistently all the observed radiation processes, we postulate that the pulsar was born with the initial period $`P_0=2`$ ms and the surface magnetic field $`B=6\times 10^{12}`$ about $`t=2\times 10^4`$ yrs ago, within or close to the dense core of the OB2 association. If this pulsar loses energy mainly on emission of electromagnetic radiation, then its present period should be $`210`$ ms (see Eq. 3). A pulsar with such parameters resembles the Vela type pulsars, i.e. it should produce the pulsed $`\gamma `$-rays with luminosity comparable to the Vela pulsar itself or e.g. PSR 1706-44. Note that these pulsars are surrounded by the synchrotron pulsar wind nebulae from which the steady TeV $`\gamma `$-ray emission is observed. It is assumed that the supernova exploded close to the dense regions of the Cyg OB2 in the medium with the density of at least $`30`$ cm3 (Butt et al. 2003), and the initial velocity of the bulk mass of the supernova envelope is 2000 km s-1 and its mass is 4 M (as derived for the Crab Nebula, Davidson & Fesen 1985). For these parameters of the pulsar and the supernova, the velocity of the supernova front wave after $`2\times 10^4`$ yrs is $`100`$ km s-1. Such a low velocity front wave can be difficult to observe due to its distortions in the interactions with the clumpy high density medium of the Cyg OB2 association. In fact, the supernova remnant shell-like structure has been revealed near the TeV source region in the CO, HI, and radio emission maps (Butt et al. 2003). The average density of matter inside this shell ($`30`$ cm-3) is much smaller than the average density of matter in the core of the Cyg OB2 equal to $`300`$ cm-3, calculated from the total mass of the Cyg OB2 $`10^5`$ M (Butt et al. 2003) and its core radius $`15`$ pc (Knödlseder 2000).

In the forthcoming sections we argue that the GeV and TeV $`\gamma `$-ray emission from the Cygnus region can be due to the Vela type pulsar and its nebula, and the excess of $`10^{18}`$ eV cosmic rays is due to the neutrons from disintegration of very relativistic nuclei which were accelerated soon after the pulsar formation and are still present inside the high density region of the Cyg OB2 association.

## 3 3EG J2033+4118 - a Vela type pulsar

The spectrum of 3EG J2033+4118 was originally described as consistent with a simple power law with the index of $`1.96\pm 0.10`$, but a double power law or a power law with the exponential cut-off give better fits (Bertsch et al. 2000, Reimer & Bertsch 2001). This source is classified as being a non-variable (see discussion in Butt et al. 2003), in contrary to its earlier classification as a variable by Mukherjee et al. (2003). 3EG J2033+4118 is believed to be located within the Cyg OB2 association (Aharonian et al. 2002, Mukherjee et al. 2003, Butt et al. 2003).

Interestingly, the $`\gamma `$-ray emission of this source resembles the spectra of the middle aged pulsars (the break at $`1`$ GeV and the sharp cut-off at a few tens of GeV, required in the spectrum of 3EG J2033+4118 by the level of TeV emission observed by the HEGRA telescope). Therefore, we suggest that the $`\gamma `$-ray source 3EG J2033+4118 is caused by a yet undiscovered Vela type pulsar within the Cyg OB2 association. In fact, also the level of the $`\gamma `$-ray emission from 3EG J2033+4118 is very close to that observed from the Vela type pulsars. Fig. 3 shows the spectra of the Vela pulsar (V), after its re-normalization from the Vela pulsar distance $`500`$ pc to the distance to the Cyg OB2 $`1.7`$ kpc, and the PSR 1706-44 (P), which is at a similar distance as the Cyg OB2. The spectrum of the Vela pulsar is measured by the EGRET (Thompson et al. 1997), COMPTEL (Schönfelder et al. 2000), OSSE (Strickman et al. 1996), and RXTE (Harding et al. 2002). The spectrum of PSR 1706-44 is measured by the EGRET (Thompson et al. 1996) and constrained by the upper limits on the pulsed emission in the X-rays by the ROSAT and ASCA (Becker et al. 1995, Finley et al. 1998). It is clear that the $`\gamma `$-ray spectra of these sources are very similar. The present parameters of the pulsar inside the Cyg OB2 can be like those ones in the example considered in the previous section, i.e. the age of $`2\times 10^4`$ yrs, the present pulsar period of $`210`$ ms, and the surface magnetic field of $`6\times 10^{12}`$ G. We encourage the searches of such order of periodicity in the $`\gamma `$-ray data of the source 3EG J2033+4118.

## 4 TeV J2032+4130 - a pulsar wind nebula

As it has already been mentioned, the Vela type pulsars are surrounded by the non-thermal nebulae observed in the X-rays and TeV $`\gamma `$-rays. The present parameters of the pulsar, which we propose as respossible for the high energy processes in the Cyg OB2, are the closest to the pulsar PSR 1706-44. It lies at a similar distance, equal to $`1.8`$ kpc, has the present period of 102 ms, and the surface magnetic field of $`3.1\times 10^{12}`$ G. The level of non-thermal X-ray emission from the nebula around PSR 1706-44, observed by the ROSAT and ASCA (Becker, Brazier & Trümper 1995, Finley et al. 1998), is close to the level of diffuse X-ray emission observed by the Chandra from the direction of the $`\gamma `$-ray source TeV J2032+4130 in the Cyg OB2 (Butt et al. 2003).

The X-ray to TeV $`\gamma `$-ray emission from the nebula around PSR 1706-44 has recently been interpreted in terms of the hadronic, leptonic model for the high energy radiation from the pulsar wind nebulae (Bednarek & Bartosik 2003). The relatively low level of X-ray emission from this nebula (as compared to e.g. the Crab Nebula) is explained if the efficiency of positron acceleration by the nuclei injected by the pulsar is relatively low $`\xi 0.08`$. On the other hand, in order to explain the level of the TeV $`\gamma `$-ray emission, an additional target of infrared photons inside the nebula is required, apart from the low energy synchrotron photons and the MBR (see Fig. 8 in Bednarek & Bartosik 2003). The $`\gamma `$-rays, produced in collisions of hadrons with the matter of the nebula, do not contribute essentially to the TeV $`\gamma `$-ray spectrum observed from this nebula due to the relatively low density of matter around PSR 1706-44.

In contrast to the PSR 1706-44 and its nebula, which were modeled by the pulsar with the initial period of 15 ms, here we postulate that the pulsar in the Cyg OB2 was born in a recent supernova explosion of one of the massive stars with a much shorter initial period, of the order of 2 ms. This supernova exploded in a much denser surrounding, estimated at $`30`$ cm-3 by Butt et al. (2003), than that one considered for the progenitor of the PSR 1706-44. Note that the average density of the Cyg OB2 association is significantly higher than the above mentioned value. This can be explained by the sweeping effect of the expanding strong wind created by the massive, pre-supernova star. Because of the shorter initial period of the pulsar and the higher density medium, the nebulae around the PSR 1706-44 and the pulsar in the Cyg OB2 should differ in some aspects.

In terms of the hadronic-leptonic model discussed by Bednarek & Bartosik (2003), we calculate the synchrotron and the IC spectra produced by leptons, which gain energy in the resonant interactions with nuclei, and the $`\gamma `$-ray spectra from decay of pions, produced in the interactions of nuclei with the matter inside the nebula. An example equilibrium spectra of leptons and nuclei inside the nebula $`2\times 10^4`$ yrs after supernova explosion are shown in Fig. 1. In these calculations it is assumed that the surface magnetic field of the pulsar born in the Cyg OB2 association is equal to $`6\times 10^{12}`$ G, and its initial period of 2 ms. Note, that the rotational energy lost by such a pulsar contributes significantly to the total energy budget of the expanding nebula.

The nuclei with the equilibrium spectrum interact with the matter of the PWNa producing $`\gamma `$-rays via decay of pions. These $`\gamma `$-ray spectra are calculated by applying the scaling break model for hadronic interactions proposed by Wdowczyk & Wolfendale (1987). The $`\gamma `$-ray spectra emitted $`2\times 10^4`$ yrs after supernova explosion are shown in Fig. 3. Note that the level of this $`\gamma `$-ray emission is consistent with the flux reported by the HEGRA group. The situation with the synchrotron and IC spectra produced by leptons is more complicated since the low energy photon target inside the nebula is not precisely known. First, we consider the case when the soft photon field is created only by the microwave background radiation (MBR) (the energy density of synchrotron photons is already negligible for a $`2\times 10^4`$ years old nebula). The results of calculations of the synchrotron and IC spectra from the nebula with the age of $`2\times 10^4`$ yrs are shown in Fig. 3, assuming that the efficiency for lepton acceleration by the nuclei is equal to $`\xi =0.05`$. Such low efficiency is required in order not to exceed the limit on the diffuse X-ray emission from the Cyg OB2 observed by the Chandra telescope (Butt et al. 2003). Next, we increase the soft photon field inside the nebula by adding the infrared photons produced by the gas inside or around the nebula. The temperature of this radiation is 100 K and the average energy density is 5 times higher than that of the MBR. The synchrotron and IC spectra calculated with the presence of this additional infrared target are shown by the thick dashed curve in Fig. 3, using this same efficiency of lepton acceleration as before. It is evident that only the model with the additional infrared target inside the nebula can produce IC $`\gamma `$-ray flux on the level consistent with reports by the HEGRA group. However, this emission has a much steeper spectrum above 1 TeV, spectral index $`2.5`$, than that one from the decay of pions. Therefore, the observations of TeV J2032+4130 at energies below 1 TeV should decide which processes, hadronic or leptonic, play the main role in this source. The comparable contribution of both processes to the $`\gamma `$-ray flux at $`1`$ TeV requires fine tuning of the acceleration efficiency of leptons by the nuclei to the value mentioned above.

## 5 Cygnus CR excess - a remnant of very young pulsar

The pulsar with the initial parameters: $`P_\mathrm{o}=2`$ ms and $`B=6\times 10^{12}`$ G, can generate electric potentials in which in principle nuclei can be accelerated above $`10^{18}`$ eV per nucleon (see Eqs. 1 and 2). Nuclei with such energies will soon leave the nebula around the pulsar and propagate in the surrounding medium. The model discussed above for the evolution of the PWNa allows us to calculate the radius of the nebula and the magnetic field inside the nebula at an arbitrary time (see Fig. 2b and also Fig. 1 in Bednarek & Bartosik 2003). Nuclei injected inside the nebula at the time, $`t_{\mathrm{acc}}`$, escape from it at the time, $`t_{\mathrm{esc}}`$, if their diffusion distance, $`R_{\mathrm{dif}}`$, in the magnetic field of the nebula is equal to the dimension of the nebula, $`R_{\mathrm{Neb}}`$, at the time $`t_{\mathrm{esc}}`$. When estimating the escape time of nuclei from the nebula we take into account the structure of the magnetic field inside the nebula (its dependence on the radius), described by the model of Kennel & Coroniti (1984). Also the time evolution of the magnetic field at the pulsar wind shock is taken into account. The diffusion of nuclei through the nebula is determined by the diffusion coefficient

$`D_{\mathrm{dif}}=(cR_L/3)^{1/2},`$ (4)
where $`R_L=E_{Fe}/eZB(t)`$ is the Larmor radius of nuclei, $`B(t)`$ is the magnetic field inside the nebula at the time t. The distance travelled by the nuclei is then obtained by integration over the time from the moment of injection up to the escape (see for details Bednarek & Protheore 2002),

$`R_{\mathrm{dif}}={\displaystyle _{t_{\mathrm{acc}}}^{t_{\mathrm{esc}}}}\sqrt{{\displaystyle \frac{3D}{2t^{}}}}𝑑t^{}.`$ (5)
The above mentioned condition, $`R_{\mathrm{dif}}=R_{\mathrm{Neb}}`$, allows us to determine the moment when nuclei with different energies and injection time escape from the nebula. The spectra of nuclei which managed to escape from the nebula during the first $`2\times 10^4`$ yrs after explosion of supernova are shown in Fig. 4, assuming the initial parameters of the pulsar and nebula as considered above. Our calculations show that nuclei with the Lorentz factors above $`10^6`$ should diffuse out of the nebula into the surrounding medium during the first $`2\times 10^4`$ yrs.

We consider the situation in which the pulsar (and its nebula) is formed inside or close to the high density regions of the Cyg OB2 association. The nuclei with the highest energies, Lorentz factors $`10^9`$, which escaped into the galactic medium with the magnetic field of $`5\times 10^6`$ G, have the Larmor radii of $`200`$ pc, and diffuse during $`2\times 10^4`$ yrs up to $`R_{\mathrm{dif}}^{\mathrm{ISM}}600`$ pc from the Cyg OB2. It is likely that a part of these nuclei can be captured by dense regions for $`10^4`$ years. The magnetic fields in dense molecular clouds are much stronger than in the Galactic medium, and can reach the values of $`10^410^3`$ G (e.g. Crutcher 1999). In such magnetic fields, the Larmor radius and the diffusion distance of nuclei with the Lorentz factors $`10^9`$ are $`10`$ pc and $`140`$ pc during the time of $`2\times 10^4`$ years. These distances are close to the radius of the central core in the Cyg OB2 estimated on $`R_{\mathrm{OB2}}15`$ pc and to its diameter of the equal to $`60`$ pc. Therefore, since the pulsar is close to the Cyg OB2, we estimate that up to $`\xi (R_{\mathrm{OB2}}/R_{\mathrm{dif}}^{\mathrm{ISM}})^310^5`$ of all injected nuclei can be captured by the nearby Cyg OB2.

The spectra of nuclei which escaped from the PWNa are shown in Fig. 4. Nuclei with different mass numbers were produced in fragmentation of heavy nuclei (the iron group), occurring during their collisions with the matter of the expanding supernova remnant during the first several years after explosion, when the column density of matter was still large. Note that two groups of heavy nuclei can be distinguished in the total spectrum (see Fig. 4). The higher energy heavy nuclei escape from the supernova envelope just after their injection into the nebula by the pulsar since their Lorentz factors are large enough that corresponding Larmor radii are comparable to the dimension of the nebula at that time. On the other hand, the lower energy heavy nuclei escape from the nebula at a later time ($`10^3`$ years), after suffering significant adiabatic energy losses due to the expansion of the nebula.

We calculate the flux of neutrons produced in collisions of nuclei with the matter inside the Cyg OB2 association. The average density of matter in the Cyg OB2 is estimated at $`300`$ cm-3, for the core radius of the Cyg OB2 equal to 15 pc and its total mass of $`10^5`$ M. It is assumed that in a single collision, the nuclei are split into two lower mass stable nuclei with the rest of nucleons released free. The energies of these nucleons are reduced by a factor of two in comparison to energies of nucleons in primary nuclei. The flux of neutrons produced in this process is shown in Fig. 3. Interestingly, the neutron flux calculated in terms of the model discussed here, is consistent with the excess of cosmic ray events at energies $`10^{18}`$ eV, reported from the Cygnus region by the AGASA group equal to $`8\%`$ of the cosmic ray flux (Hayashida et al. 1999b).

## 6 Conclusion

We have shown that different types of high energy phenomena observed from the direction of the Cyg OB2 association can be self-consistently explained by the model of energetic pulsar created in the recent explosion of one of the massive stars. The pulsar should have the initial period close to 2 ms. For the present age of $`2\times 10^4`$ yrs, its period should be close to $`210`$ ms, if its surface magnetic field is $`6\times 10^{12}`$ G. Such a pulsar belongs to the family of the so-called Vela type pulsars, producing pulsed $`\gamma `$-rays with the intensity comparable to the EGRET source, 3EG J2033+4118, observed in the direction of the Cyg OB2. In fact, the spectral features of 3EG J2033+4118 (possible break in the spectrum, required high energy cut-off) are characteristic for the spectra of the Vela pulsar and the PSR 1706-44. Therefore, we encourage the searches of the hundred millisecond periodicity in the $`\gamma `$-ray sources reported in the direction of the Cyg OB2.

The pulsar in the Cyg OB2 should also be surrounded by the pulsar wind nebula, as observed in the case of other Vela type pulsars. We propose that the TeV $`\gamma `$-ray source reported by the HEGRA group is due to the $`\gamma `$-rays produced by leptons and/or hadrons accelerated by the PWNa inside or close to the Cyg OB2. Therefore, in spite of the close physical relation, the TeV $`\gamma `$-ray emission from the Cyg OB2 (produced in the PWNa) do not need to connect smoothly with the GeV $`\gamma `$-ray emission (produced in the pulsar magnetosphere).

We also suggest that the cosmic ray excess, observed from the Cygnus region by the AGASA experiment, is also physically related to this energetic pulsar. The highest energy nuclei, injected by the pulsar soon after its formation, escape from the PWNa. A part of these nuclei can be captured by the strong magnetic fields in the high density regions of the Cyg OB2. These nuclei collide with the matter of the Cyg OB2, injecting neutrons which create the observed excess of cosmic rays in the direction of the Cyg OB2.

The $`\gamma `$-rays from the Cyg OB2, produced in collisions of hadrons with the matter, should also be accompanied by the comparable fluxes of neutrinos. However, the predicted muon neutrino event rate in a 1 km2 detector of the IceCube type during a 1 yr observation is rather low (based on the calculations of the likelihood of detecting neutrinos by such a detector obtained by Gaisser & Grillo (1987), and the absorption coefficients of neutrinos in the Earth derived by Gandhi (2000)). We estimate that hadrons captured in the dense regions of the Cyg OB2 should give between $`0.14`$ to $`0.05`$ neutrino events per year, for the source observed close to the horizon (not absorbed by the Earth) and from the Nadir direction (absorbed by the Earth), respectively. The hadrons confined in the PWNa should produce $`0.48`$ ($`0.46`$) events per yr for the source located at the horizon (the Nadir). The spectra of neutrinos produced by hadrons in the PWNa and the Cyg OB2 are shown in Fig. 5. Neutrino spectra produced in the Cyg OB2 are above the atmospheric neutrino background expected within 1o of the source but they are below the sensitivity of the 1 km2 neutrino detector. On the other hand, the neutrino spectra produced inside the PWNa are closer to the sensitivity limit of the 1 km2 detector but they are significantly below the atmospheric neutrino background. Therefore, we do not predict any detection of neutrinos from the direction of the Cyg OB2 association.

Finally, we conclude that the presence of at least one Vela type pulsar in the Cyg OB2 association is quite likely on statistical grounds, since the supernova rate in this association should be one per a few tens of hundreds years on average (the Cyg OB2 contains about 100 O type stars with the lifetime of a few milion years).

## Acknowledgements

I would like to thank the referee, F. A. Aharonian, for useful comments. This work is supported by the Polish KBN grants No. 5P03D02521 and PBZ KBN 054/P03/2001.


format: text

text:
# 1 Introduction

## 1 Introduction

Hadroproductions of charmed particles have been measured in fixed-target experiments  and indicate the large leading/non-leading asymmetries, i.e. asymmetries between production cross sections of leading particles and those of non-leading particles defined as

$`A={\displaystyle \frac{\sigma (leading)\sigma (nonleading)}{\sigma (leading)+\sigma (nonleading)}}.`$ (1)
The leading particle shares a valence quark in common with the incident hadron, while the non-leading one does not. The leading particles are copiously produced at large Feynman $`x`$ as compared with the non-leading particles in the forward region of the incident hadron. This is called as leading particle effect. For example, the asymmetry between $`D^{}(d\overline{c})`$ and $`D^+(c\overline{d})`$ in $`\pi ^{}(d\overline{u})`$ interaction with a nucleon is defined as $`A_{\pi ^{}N}(D^{},D^+)=(d\sigma (D^{})d\sigma (D^+))/(d\sigma (D^{})+d\sigma (D^+))`$ and increases from zero to nearly one with $`x`$ in the $`\pi ^{}`$ fragmentation region. In the conventional perturbative QCD at leading order (LO), the factorization theorem predicts that $`c`$ and $`\overline{c}`$ quarks are produced symmetrically and then fragment into $`D`$ and $`\overline{D}`$ independently. As a consequence, LO perturbative QCD predicts no asymmetry in contrast with the experimental data. The next to leading order (NLO) calculations generates asymmetries but they are much smaller than the data.

To explain the leading/non-leading asymmetry, many attempts have been investigated: $`k_T`$ factorization, string fragmentation, intrinsic charm contributions, recombination process with beam remnants, recombination with participants in the hard scattering process, recombination with surrounding light quarks in projectile and target, recombination using valon concept, light quark fragmentation, meson cloud model and so on. Productions of $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ in $`\mathrm{\Sigma }^{}A`$  and $`\mathrm{\Sigma }^{}A,pA`$ and $`\pi ^{}A`$ collisions have been measured and analyzed in Refs. and . These models explain the leading/non-leading asymmetry of charmed hadrons successfully, however, there are some uncertainties in the perturbative QCD calculations as reviewed in Refs. and . Hadronization of heavy quarks is still an open problem, since it is intrinsically non-perturbative process.

We have applied the covariant quark-diquark cascade model for hadroproductions of charmed mesons as well as hadrons composed $`u,d`$ and $`s`$ quarks. In the present paper, we investigate the hadroproductions of $`\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c,D`$ and $`\overline{D}`$ and the leading/non-leading asymmetry in the framework of the constituent quark-diquark cascade picture and discuss the role of constituent diquark.

## 2 The model

### 2.1 Valence quark distributions in the incident hadron

We assume that a baryon is composed of a constituent quark and a constituent diquark, and a meson is composed of a constituent quark and a constituent anti-quark. The constituent quark (diquark) is a quark-gluon cluster consisting of the valence quark (diquark), sea quarks and gluons. Hereafter, we simply call them as quark and diquark. When the collision between the incident hadrons $`A`$ and $`B`$ occurs, it is assumed that the incident hadrons either break up into two constituents or emit gluons followed by a quark-antiquark pair creation. There are four interaction types: a) non-diffractive dissociation, b) and c) single-diffractive dissociations of $`B`$ and $`A`$, and d) double-diffractive dissociation types as shown in Fig.1. The probabilities of these types to occur are $`(1P_{gl})^2,P_{gl}(1P_{gl}),P_{gl}(1P_{gl})`$ and $`P_{gl}^2`$, respectively. Here we denote the quark-antiquark pair emitted from $`A`$ ($`B`$) via gluons as $`M_A`$ ( $`M_B`$). The probabilities of $`M_A`$ ( $`M_B`$) to be $`u\overline{u},d\overline{d},s\overline{s}`$ and $`c\overline{c}`$ are are denoted by $`P_{u\overline{u}},P_{d\overline{d}},P_{s\overline{s}}`$ and $`P_{c\overline{c}}`$, respectively.

The light-like variables of $`A`$ and $`B`$ are defined as follows:

$`x_{0\pm }^A={\displaystyle \frac{E^A\pm p_{cm}}{\sqrt{s_0}}},x_{0\pm }^B={\displaystyle \frac{E^Bp_{cm}}{\sqrt{s_0}}},`$ (2)
where $`\sqrt{s_0}=E^A+E^B`$ is the center of mass energy of the incident hadrons $`A`$ and $`B`$. The momentum fractions of $`A,B`$ and $`M_A`$ for the process shown in Fig.1b are given as follows: The momentum fraction of $`M_A`$ is fixed by the distribution function

$`H_{M_A/A}(z)=z^{\beta _{gl}1}(1z)^{\beta _{ld}1}/B(\beta _{gl},\beta _{ld}),`$ (3)
and the uniform distribution $`R`$ in the interval from zero to one as,

$`x_+^{M_A}=x_{0+}^Az,x_{}^{M_A}=x_0^AR.`$ (4)
Then the momentum fractions of the incident particles $`A`$ and $`B`$ become as follows:

$`x_+^A=x_{0+}^A(1z),x_{}^A=m_A^2/(x_+^As_0),`$
$`x_{}^B=x_0^B(x_{}^Ax_0^A(1R)),x_+^B=x_{0+}^B,`$ (5)
where the mass shell condition is considered and transverse momenta are neglected. Analogous expressions describe the momentum fractions of $`A,B`$ and $`M_B`$ for the process shown in Fig.1c , upon replacing superscripts $`AB`$ and subscripts $`+`$. The single-diffractive dissociation type mechanisms result in diffractive peaks of spectra of the same kind of the incident particles $`A`$ and $`B`$ at $`|x|1`$. The probability of the incident hadron to emit gluons followed by a quark-antiquark pair $`P_{gl}`$ and the parameters $`\beta _{gl}`$ and $`\beta _{ld}`$ in (3) are chosen to reproduce the diffractive peaks of the spectra such as $`\pi ^\pm ,K^\pm `$ and $`p`$ in $`\pi ^\pm p,K^\pm p`$ and $`pp`$ collisions at $`|x|1`$. The parameters are chosen as $`P_{gl}=0.15,\beta _{gl}=0.1`$ and $`\beta _{ld}=3.0`$. However, the main contributions to charmed hadrons come from the non-diffractive dissociation type mechanism.

The distribution functions of the constituents in the projectile $`A`$ composed of $`a`$ and $`a^{}`$ are described as

$`H_{a/A}(z)=H_{a^{}/A}(1z)={\displaystyle \frac{z^{\beta _a1}(1z)^{\beta _a^{}1}}{B(\beta _a,\beta _a^{})}}.`$ (6)
Similarly, the distribution function of $`q`$ in $`M_A`$ is given as $`H_{q/M_A}(z)=z^{\beta _q1}(1z)^{\beta _q1}/B(\beta _q,\beta _q)`$. The dynamical parameters $`\beta `$’s in (6), which determine the momentum sharing of the constituents, are related to the intercepts of the Regge trajectories as $`\beta _u=\beta _d=1\alpha _{\rho \omega }(0)0.5,\beta _s=1\alpha _\varphi (0)1.0,\beta _c=1\alpha _{J/\psi }(0)4.0`$. For flavor anti-symmetric diquark $`[qq^{}]`$ and symmetric diquark $`\{qq^{}\}`$, the dynamical parameters $`\beta _{[qq^{}]}`$ and $`\beta _{\{qq^{}\}}`$ are chosen as $`\beta _{[qq^{}]}=\gamma _{[]}(\beta _q+\beta _q^{})`$ and $`\beta _{\{qq^{}\}}=\gamma _{\{\}}(\beta _q+\beta _q^{})`$, respectively. The commonly shared constituents of $`\mathrm{\Lambda }_c`$ with the incident $`p,\mathrm{\Sigma }^{}`$ and $`\pi ^{}`$ beams are $`ud`$-diquark, $`d`$ and $`d`$-quark respectively. $`\overline{\mathrm{\Lambda }}_c`$ has a commonly shared $`\overline{u}`$ quark with $`\pi ^{}`$ and $`K^{}`$ beams and $`\mathrm{\Lambda }_c`$ is a non-leading particle in $`K^{}`$ beam. These distributions of the constituents in the incident hadrons affect the shapes of spectra of produced hadrons significantly. Hadrons are produced on the chain between a constituent valence quark($`q_A`$)(anti-quark($`\overline{q}_A`$)) from the incident hadron $`A(M_A)`$ and the constituent valence diquark($`qq_B(\overline{q}_B)`$)(quark($`q_B`$)) from the incident hadron $`B(M_B)`$ for a meson-baryon collision as illustrated in Fig.1.

The small intrinsic transverse momenta of constituents in the incident hadrons are introduced by changing the direction of the incident constituents by an angle $`\phi _{in}`$ in the rest frame of the chain. Here we choose the distributions for $`z=\mathrm{cos}\phi _{in}`$ as

$`D_{in}(z)={\displaystyle \frac{\beta _{in}+1}{2^{\beta _{in}+1}}}(1+z)^{\beta _{in}}`$ (7)
in the region $`1<z<1`$. Here we choose very small intrinsic transverse momenta and use the value $`\beta _{in}=40`$ used in Ref. in order to reproduce $`p_T^2`$ dependences of $`D`$ mesons.

We take account of the annihilation of the valence antiquark from beam and valence quark from target in a slightly different way from the one pointed out in Refs.  and . When the incident constituents in $`A`$($`M_A`$) and $`B(M_B)`$ are quark and its anti-quark, it is assumed that the annihilation process occurs with the probability $`P_{an}`$ and the non-annihilation process occurs with $`1P_{an}`$. Let us consider the case of $`\pi ^{}p`$ collision. The annihilation process

$`\overline{u}+u\overline{q}+q\mathrm{"}hadrons\mathrm{"}+\mathrm{"}hadrons\mathrm{"}`$ (8)
occurs with the probability $`P_{an}`$, when $`\overline{u}`$ in $`\pi ^{}`$ interacts with $`u`$ in $`p`$. Branching ratios of $`\overline{u}u\overline{u}u,\overline{d}d,\overline{s}s,`$ and $`\overline{c}c`$ are chosen to be equal to $`P_{u\overline{u}},P_{d\overline{d}},P_{s\overline{s}}`$ and $`P_{c\overline{c}}`$ for the channels allowed energetically. The produced $`\overline{q}`$ and $`q`$ in (8) are supposed to be non-free in terms of hadronization mechanism. Thus it is assumed that $`\overline{q}`$ has tendency to be produced in the forward direction of $`\overline{u}`$ in $`\pi ^{}`$ beam. Therefore, here we choose the distributions for $`z=\mathrm{cos}\phi _{an}`$ as

$`D_{an}(z)={\displaystyle \frac{\beta _{an}+1}{2^{\beta _{an}+1}}}(1+z)^{\beta _{an}}`$ (9)
in the region $`1<z<1`$, where $`\phi _{an}`$ is the angle between the directions of $`\overline{u}`$ and $`\overline{q}`$. To reproduce the $`A_{\pi ^{}N}(D^0,\overline{D}^0)`$, we choose the value $`P_{an}=0.2`$ and use a value of $`\beta _{an}=20`$ which roughly means $`D_{an}(z)\delta (z1)`$ .

### 2.2 Quark-diquark cascade process

Hadrons are produced on the cascade chain by processes

$`q`$ $``$ $`M(q\overline{q}^{})+q^{},`$ (10)
$`q`$ $``$ $`B(q[q^{}q^{\prime \prime }])+\overline{[q^{}q^{\prime \prime }]},B(q\{q^{}q^{\prime \prime }\})+\overline{\{q^{}q^{\prime \prime }\}},`$ (11)
$`\overline{[q^{}q^{\prime \prime }]}`$ $``$ $`\overline{B}(\overline{q}\overline{[q^{}q^{\prime \prime }]})+q,`$ (12)
$`\overline{[q^{}q^{\prime \prime }]}`$ $``$ $`M(q\overline{q}^{})+\overline{[qq^{\prime \prime }]},M(q\overline{q}^{})+\overline{\{qq^{\prime \prime }\}},`$ (13)
$`\overline{\{q^{}q^{\prime \prime }\}}`$ $``$ $`\overline{B}(\overline{q}\overline{\{q^{}q^{\prime \prime }\}})+q,`$ (14)
$`\overline{\{q^{}q^{\prime \prime }\}}`$ $``$ $`M(q\overline{q}^{})+\overline{[qq^{\prime \prime }]},M(q\overline{q}^{})+\overline{\{qq^{\prime \prime }\}},`$ (15)
where $`q`$ denotes $`u,d,s`$ and $`c`$ and $`[q^{}q^{\prime \prime }]`$ does $`[ud],[us],[uc],[ds],[dc],[sc]`$ and so on. The symbols $`[]`$ and $`\{\}`$ denote the flavor anti-symmetric and symmetric diquarks, respectively. Meson production probabilities from $`q,\overline{[q^{}q^{\prime \prime }]}`$ and $`\overline{\{q^{}q^{\prime \prime }\}}`$ are $`1ϵ,\eta _{[]}`$ and $`\eta _{\{\}}`$, respectively. The parameter $`ϵ`$ is the probability of baryon productions from a quark and is related to the baryon productions from $`\pi `$ and $`K`$ beams resulting in $`ϵ=0.07`$. The processes (13) and (15) contribute to production of baryons which share no valence quark or diquark in common with the incident hadron. Although both $`\mathrm{\Omega }`$ and $`\overline{\mathrm{\Omega }}`$ share no valence quark common with the incident proton, the spectrum of $`\mathrm{\Omega }`$ is harder than that of $`\overline{\mathrm{\Omega }}`$ in the proton fragmentation region. In our model this is described as follows: after emitting $`K`$ mesons through these processes, the incident diquark in the incident proton converts into $`\{ss\}`$ diquark and produce $`\mathrm{\Omega }`$ by the process (14). The difference between the spectra of $`\mathrm{\Omega }`$ and $`\overline{\mathrm{\Omega }}`$ in $`pp`$ collision is described by the processes (13),(15) and (14) and the parameters are chosen as $`\eta _{[]}=\eta _{\{\}}=0.25`$.

For the cascade process $`aH(a\overline{b})+b`$ in the beam (target) side, the distribution function of light-like fraction of $`b`$ is assumed as

$`F_{ba}(z)=z^{\gamma \beta _a1}(1z)^{\beta _a+\beta _b1}/B(\gamma \beta _a,\beta _a+\beta _b)`$ (16)
and the light-like fraction of $`b`$ becomes $`x_+^b=zx_+^a(x_{}^b=zx_{}^a)`$. We fix the transverse momentum of the hadron $`H`$ as a difference between two random points $`𝒑_T^H=𝒑_{T2}𝒑_{T1}`$ with their lengths fixed by the distribution function

$`G(𝒑_T^2)={\displaystyle \frac{\sqrt{m_H}}{C}}\mathrm{exp}({\displaystyle \frac{C}{\sqrt{m_H}}}𝒑_T^2)`$ (17)
in $`𝒑_T^2`$ space, where $`m_H`$ denotes the mass of $`H`$. The parameter $`C`$ is fixed from the experimental data on $`𝒑_T^2`$ distributions of pions. We use the values of $`\gamma =1.75`$ and $`C=1.8`$GeV$`^{\frac{3}{2}}`$, from previous analyses.

We now see how the cascade processes occurs between the incident constituents, for example, $`q`$ from $`A`$ $`(M_A)`$ and $`\{q^{}q^{\prime \prime }\}`$ from $`B`$. We redefine the light-like fractions of $`q`$ and $`\{q^{}q^{\prime \prime }\}`$ in the rest frame of the incident constituents $`q`$ and $`\{q^{}q^{\prime \prime }\}`$. The momentum sharing of the cascade process $`q+\{q^{}q^{\prime \prime }\}H(q\overline{b})+b+\{q^{}q^{\prime \prime }\}`$ from $`q`$ with $`x_\pm ^q`$ and $`\{q^{}q^{\prime \prime }\}`$ with $`x_\pm ^{\{q^{}q^{\prime \prime }\}}`$ takes place as follows: Using Eq.(16) with $`a=q`$, we fix the light-like fractions of $`b`$ and $`H(q\overline{b})`$ as $`x_+^b=x_+^qz`$ and $`x_+^H=x_+^qx_+^b`$, respectively and put $`x_{}^b=x_{}^q`$. The transverse momentum of $`H`$ is fixed by Eq.(17). Then $`x_{}^H`$ is fixed from the onshell condition $`x_+^Hx_{}^H=(m_H^2+𝒑_T^2)/s^{}`$, where $`\sqrt{s^{}}`$ is the subenergy of the incident $`q`$ and $`\{q^{}q^{\prime \prime }\}`$ system. The transverse momentum of $`b`$ is $`𝒑_T^b=𝒑_T^q𝒑_T^H`$. The light-like fraction of $`\{q^{}q^{\prime \prime }\}`$ is decreased to $`\stackrel{~}{x}_{}^{\{q^{}q^{\prime \prime }\}}=x_{}^{\{q^{}q^{\prime \prime }\}}x_{}^H`$. If the energy of the $`b+\{q^{}q^{\prime \prime }\}`$ system is enough to create another hadron, the cascade process such as $`b+\{q^{}q^{\prime \prime }\}b+H(c\{q^{}q^{\prime \prime }\})+\overline{c}`$ takes place in the opposite side. In the final step, we assume that the constituents recombine into one or two hadrons according to the processes: $`q+\overline{q}^{}M(q\overline{q}^{}),q+[q^{}q^{\prime \prime }]B(q[q^{}q^{\prime \prime }]),\{qq^{}\}+\overline{[q^{\prime \prime }q^{\prime \prime \prime }]}M(q\overline{q^{\prime \prime }})+M(q^{}\overline{q^{\prime \prime \prime }})`$ and so on. The momenta of the recombined hadrons are the sum of those of the final constituents and are offshell. In the cms of incident hadrons $`A`$ and $`B`$, we have energy-momentum conservation relations $`\underset{i}{}𝒑_i=0`$ and $`\underset{i}{}p_i^0=\sqrt{s_0}`$, where $`𝒑_i`$ denotes the three momentum of $`i`$-th produced particle. In order to put recombined hadrons onshell, we multiply the three momenta of all produced hadrons by a common factor $`f`$ so that the summation $`\underset{i}{}\sqrt{(f𝒑_i)^2+m_i^2}`$ would be equal to $`\sqrt{s_0}`$.

We take into account pseudoscalar $`(Ps)`$, vector $`(V)`$ and tensor $`(T)`$ mesons. The production probabilities for them are assumed to be $`P_{PS}=1/9,P_V=3/9`$ and $`P_T=5/9`$. For baryons we consider lower lying baryons: octet ($`B_8`$) and decuplet ($`B_{10}`$) baryons composed of $`u,d,`$ and $`s`$ flavors and the corresponding ones with charm flavor. Octet and decuplet baryons are described as

$`|B_8>=\mathrm{cos}\theta |q[q^{}q^{\prime \prime }]>+\mathrm{sin}\theta |q\{q^{}q^{\prime \prime }\}>,`$ (18)
$`|B_{10}>=|q\{q^{}q^{\prime \prime }\}>.`$ (19)
We assume the production probabilities for octet and decuplet baryons from symmetric diquarks to be $`P_{B8}=1/3`$ and $`P_{B10}=1P_{B8}=2/3`$, respectively. Then the emission probabilities of individual cascade processes are determined from the above probabilities. For examples, the probabilities to produce $`\pi ^+,\rho ^+,a_2^+,\mathrm{},\mathrm{\Delta }^{++},p,\mathrm{\Delta }^+,\mathrm{},\mathrm{\Xi }_{cc}^{++}`$ and $`\mathrm{\Xi }_{cc}^{++}`$ from a $`u`$ quark are as follows:

$`(1ϵ)P_{d\overline{d}}P_{PS},(1ϵ)P_{d\overline{d}}P_V,(1ϵ)P_{d\overline{d}}P_T,\mathrm{},ϵP_{u\overline{u}}^2,ϵP_{u\overline{u}}P_{d\overline{d}}\frac{1}{3}P_{B8}\mathrm{sin}^2\theta /(\frac{1}{3}P_{B8}\mathrm{sin}^2\theta +\frac{2}{3}P_{B10}),ϵP_{u\overline{u}}P_{d\overline{d}}\frac{2}{3}P_{B10}/(\frac{1}{3}P_{B8}\mathrm{sin}^2\theta +\frac{2}{3}P_{B10}),\mathrm{},ϵP_{c\overline{c}}^2\frac{2}{3}P_{B8}\mathrm{sin}^2\theta /(\frac{2}{3}P_{B8}\mathrm{sin}^2\theta +\frac{1}{3}P_{B10}),ϵP_{c\overline{c}}^2\frac{1}{3}P_{B10}/(\frac{2}{3}P_{B8}\mathrm{sin}^2\theta +\frac{1}{3}P_{B10}),`$

where the factors 1/3 and 1/2 are flavor SU(4) factors. Directly produced resonances decay into stable particles.

The mixing angle of symmetric and anti-symmetric diquarks in octet baryon $`\mathrm{cos}\theta `$ and the pair creation probabilities of $`u\overline{u},d\overline{d}`$ and $`s\overline{s}`$ are determined from the analyses of non-charmed hadron productions in $`\pi p,Kp`$ and $`pp`$ collisions. In Ref., we used $`\mathrm{cos}^2\theta =\frac{2}{3}`$ and $`P_{u\overline{u}}=0.45,P_{d\overline{d}}=0.45`$ and $`P_{s\overline{s}}=0.10`$. To reproduce the data on $`\pi ^{}pD^\pm X`$ cross section, we choose $`\mathrm{cos}^2\theta =\frac{1}{2},P_{c\overline{c}}=0.0003,P_{u\overline{u}}=0.435,P_{d\overline{d}}=0.435`$ and $`P_{s\overline{s}}=0.1297`$. The dynamical parameters $`\gamma _{[]}`$ and $`\gamma _{\{\}}`$ as well as the mixing angle of symmetric and anti-symmetric diquarks in octet baryon is related to the $`\pi ^\pm `$ spectra in proton fragmentation region and are chosen as $`\gamma _{[]}=1.5`$ and $`\gamma _{\{\}}=2.0`$.

### 2.3 Interaction with nucleus

We now consider hadron-nucleus interaction. We regard the hadron-nucleus collision as a sum of hadron-nucleon collisions and neglect the intra-nuclear cascade of produced hadrons. The projectile hadron successively collides with nucleons inside the nucleus. The probability $`P_{hA}(\nu )`$ that the incident hadron $`h`$ collides with $`\nu `$ nucleons is calculated by the Glauber-type multiple collision model. The number $`\nu `$ is determined by the distribution of nucleons in the nucleus and the cross section of the incident hadron with a nucleon $`\sigma _{hN}`$. We take the values $`\sigma _{pN}=40,\sigma _{\mathrm{\Sigma }^{}N}=34,\sigma _{\pi ^{}N}=24`$ and $`\sigma _{K^{}N}=21`$mb in the hadron-nucleus collisions. The nucleon number density of the nucleus with a mass number $`A`$ is assumed as

$`\rho (r)=\rho _0{\displaystyle \frac{A}{1+\mathrm{exp}(\frac{rr_A}{d})}},`$ (20)
where we choose $`r_A=1.19A^{\frac{1}{3}}1.61A^{\frac{1}{3}}`$ (fm) and $`d=0.54`$ (fm).

Then it is assumed that $`M_1N_1,M_2N_2,\mathrm{},`$ $`M_{\nu 1}N_{\nu 1}`$ and $`hN_\nu `$ interactions take place, where $`M_i`$ is $`q\overline{q}`$ state emitted from $`h`$ as $`M_A`$ discussed in section 2.1. At each collision, the projectile hadron $`h`$ looses its momentum. From the data on $`A`$ dependences of the spectra of the projectile hadrons, we put the degradation rate of the momentum of $`h`$ as $`P(z)=0.25z^{0.251}`$ which is used in Ref..

## 3 Comparison with experiment

In this section, we give the results of our model for charmed hadron productions in $`p,\mathrm{\Sigma }^{},\pi ^{}`$ and $`K^{}`$ beam interactions with nuclei. To see the target mass number $`(A)`$ dependence, we show the $`x`$ distributions of $`\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c,D`$ and $`\overline{D}`$ productions in $`p,\mathrm{\Sigma }^{},\pi ^{}`$ and $`K^{}`$ interactions with $`N`$ and $`Pb`$ targets at $`p_L=600`$ GeV/c in Fig.2. All spectra increase with $`A`$ in the region $`x\stackrel{<}{}0`$. The spectra of $`\mathrm{\Lambda }_c`$ become soft with $`A`$ in the baryon beam fragmentation region. This is due to the momentum degradation of the incident baryon through the multiple interactions with nucleons in the target.

The spectra of $`\overline{\mathrm{\Lambda }}_c`$ scarcely depend on $`A`$ in the baryon beam fragmentation regions. In meson beam fragmentation regions, both $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ slightly depend on $`A`$ and become soft with increasing $`A`$. In baryon beam fragmentation regions, $`A`$-dependences of the spectra of $`\overline{D}`$ mesons are smaller than those for $`D`$ mesons. The shapes of the spectra for $`D`$ mesons in baryon beam fragmentation regions are steep as compared with those of $`\overline{D}`$ mesons in baryon and $`D`$ and $`\overline{D}`$ mesons in meson beam fragmentation regions.

In Fig.3, we show the $`x`$ distributions of $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ in $`p,\mathrm{\Sigma }^{},\pi ^{}`$ and $`K^{}`$ beam interactions with $`C`$ and $`Cu`$ targets in arbitrary unit and compare the results with the experimental data E781. The shapes of $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ in $`p,\mathrm{\Sigma }^{}`$ and $`\pi ^{}`$ fragmentation regions are well explained by our model. In our model, the shapes of the spectra are strongly related to the distributions of valence constituents in the incident hadrons, which are characterized by the parameters $`\beta `$’s in (6). Diquarks in the incident baryons have a tendency to bring a large momentum fraction of the incident baryons. The valence constituents $`d`$ and $`\overline{u}`$ in $`\pi ^{}`$ beam have the same distribution i.e. energetic $`d`$ and wee $`\overline{u}`$ or vice versa. The valence constituents $`d`$ and $`\overline{u}`$ in $`\pi ^{}`$ are harder than the valence quarks in $`\mathrm{\Sigma }^{}`$ and $`p`$ beams. The $`ds`$ diquark in $`\mathrm{\Sigma }^{}`$ is harder than the $`ud`$ diquark in $`p`$. The valence $`\{ud\}`$ diquark in $`p`$ beam is the hardest among the commonly shared constituents between $`\mathrm{\Lambda }_c`$ and the incident $`p,\mathrm{\Sigma }^{},\pi ^{}`$ and $`K^{}`$ beams. Therefore $`x`$-distribution of $`\mathrm{\Lambda }_c`$ in proton beam is harder than those in other hadron beams.

The asymmetries between the spectra $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ in $`pA,\mathrm{\Sigma }^{}A,\pi ^{}A`$ and $`K^{}A`$ collisions are shown in Fig.4. The results are in good agreement with the experimental data except for $`\pi ^{}`$ beam. For $`\pi ^{}`$ beam, $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ share the common valence constituents $`d`$ and $`\overline{u}`$, respectively. Therefore the same spectra are expected for $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ productions in the region of $`\pi ^{}`$ beam fragmentation. In our model, however, the excess of $`\overline{\mathrm{\Lambda }}_c`$ over $`\mathrm{\Lambda }_c`$ is seen at $`0.3\stackrel{<}{}x`$. This is due to the fact that energetic $`\mathrm{\Lambda }_c`$ and $`\overline{\mathrm{\Lambda }}_c`$ are produced in the $`d`$-quark-diquark and the $`\overline{u}`$-antiquark-quark chains, respectively. The former chain is slightly shifted to the target fragmentation region as compared with the latter one due to the large momentum fraction of the diquark. Therefore our model gives the negative asymmetry $`A_{\pi ^{}A}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$ at $`0.3\stackrel{<}{}x`$. The valence diquark in $`\mathrm{\Sigma }^{}`$ beam contributes to $`\mathrm{\Lambda }_c`$ production via a processes such as $`dsdu+K^{}\overline{c}+\mathrm{\Lambda }_c+K^{}`$ in $`\mathrm{\Sigma }^{}`$ fragmentation region. The contributions to $`\mathrm{\Lambda }_c`$ production from the valence diquarks dominate the contributions from the valence quarks in $`\mathrm{\Sigma }^{}`$ beam.

In Fig.5, We show the results of asymmetries between charmed and anti-charmed mesons in $`p`$ and $`\mathrm{\Sigma }^{}`$ interactions with $`N,C`$ and $`Cu`$ targets. The asymmetries $`A_{pA}(D^{},D^+),A_{pA}(\overline{D}^0,D^0),A_{\mathrm{\Sigma }^{}A}(D_s^{},D_s^+)`$ and $`A_{\mathrm{\Sigma }^{}A}(D^{},D^+)`$ are small as compared with the asymmetries $`A_{pA}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$ and $`A_{\mathrm{\Sigma }^{}A}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$. Although both $`D_s^+`$ and $`D_s^{}`$ in $`p`$ and $`D^0`$ and $`\overline{D}^0`$ in $`\mathrm{\Sigma }^{}`$ beam are non-leading particles, the results of asymmetries $`A_{pN}(D_s^{},D_s^+)`$ and $`A_{\mathrm{\Sigma }^{}A}(\overline{D}^0,D^0)`$ are non-zero and have large values. This is because in our model $`\overline{D}`$ mesons are produced mainly from the chain between the valence diquark in the incident baryon and the valence quark in the target nucleon. For example, $`D_s^{}`$ in proton fragmentation region is produced via a process such as $`udsd+K^+cd+D_s^{}+K^+`$. The more energetic the valence diquark becomes, the more copiously $`\overline{D}`$ mesons are produced in the baryon fragmentation region. In our model the valence $`ds`$ diquark in $`\mathrm{\Sigma }^{}`$ beam has a larger momentum fraction than the valence diquarks in $`p`$ beam. Consequently the asymmetries $`A_{\mathrm{\Sigma }^{}A}(\overline{D}^0,D^0)`$ are large as compared with $`A_{pA}(D_s^{},D_s^+)`$. The asymmetry $`A_{\mathrm{\Sigma }^{}A}(D_s^{},D_s^+)`$ is larger than other charmed meson asymmetries in $`p`$ and $`\mathrm{\Sigma }^{}`$ beams. This is also due to the large momentum fraction of the valence $`ds`$ diquark in $`\mathrm{\Sigma }^{}`$. In the target fragmentation region, both $`D`$ and $`\overline{D}`$ mesons increase similarly with the mass number $`A`$ and the asymmetries show small mass number dependences. Since the spectra of $`D`$ mesons increase with $`A`$ and those of $`\overline{D}`$ mesons change little, there are considerable target mass number dependence in the asymmetries $`A(\overline{D},D)`$ at $`0\stackrel{<}{}x`$.

## 4 Conclusions

1. We have examined hadroproductions of charmed baryons in the frame work of the quark-diquark cascade model. Their spectra are well explained by using the values of the dynamical parameters fixed from Regge intercepts. The large leading/non-leading asymmetries $`A_{pA}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$ and $`A_{\mathrm{\Sigma }^{}A}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$ are naturally explained.

2. The valence diquark in the incident baryon, which brings a large momentum fraction of the beam, plays an important role for $`\overline{D}`$ productions as well as for $`\mathrm{\Lambda }_c`$ production. The non-leading anti-charmed meson productions in baryon beams are largely affected by the valence diquark distribution in the incident baryons through the processes (13) and/or (15).

3. In meson beam fragmentation regions, $`A`$-dependences of the spectra of non-leading charmed mesons are small as compared with those for leading charmed mesons. There are considerable $`A`$-dependences in $`\mathrm{\Lambda }_c`$ and $`D`$ meson productions, while those for $`\overline{\mathrm{\Lambda }}_c`$ and $`\overline{D}`$ productions are small in baryon beam fragmentation regions.

4. There are discrepancies between the data and calculations of asymmetries $`A_{\pi ^{}A}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$ at $`\pi ^{}`$ fragmentation region as seen above. Piskounova discussed this problem successfully by introducing the transfer of string junction in quark-gluon string model. To explain the experimental data in our model, it is necessary, for example, a energetic $`d`$ quark from $`\pi ^{}`$ beam and a wee $`uc`$ diquark from target nucleon to recombine into a $`\mathrm{\Lambda }_c`$. Here the $`uc`$ diquark is converted from the valence diquark in the target nucleon by emitting several mesons through the processes (13) and/or (15). The significant value of the asymmetry $`A_{\pi ^{}A}(\mathrm{\Lambda }_c,\overline{\mathrm{\Lambda }}_c)`$ as well as the small value of $`A_{\pi ^{}N}(D^0,\overline{D}^0)`$ in meson beam fragmentation is a challenging problem.

## Acknowledgments

This work was partly supported by grant of the Academic Frontier Project promoted by the Ministry of Education, Culture, Sports, Science and Technology, Japan.


text:
# Atomic and Electronic Structures of Unreconstructed Polar MgO(111) Thin Film on Ag(111)

## I INTRODUCTION

Polar surfaces have attracted wide attention not only for fundamental science but also for technological applications, because several interesting properties, such as novel catalytic activity, two dimensional electron system, are expected for the polar surface. Figure 1 shows the atomic geometry of non-polar (100) and polar (111) faces of rocksalt crystals schematically. The (111) surface consists of only one atomic species, either cations or anions, while the (100) surface consists of the same number of cations and anions. In a rocksalt structure crystal, each atom is surrounded by six atoms of different species. The coordination number decreases to five for atoms on a (100) surface, while it changes drastically to three for atoms on a (111) surface. Therefore, surface energy of the (111) face is much higher than that of the (100) face, which makes the (111) surface unstable. In other words, alternate stacking of cation and anion layers forms a dipole layer along the  direction (Fig. 2-a). Accumulation of the dipole layers would produce a macroscopic electric field, which causes appearance of a flat (111) surface to be quite unpreferable in the rocksalt structure compounds and the flat (111) surface does not occur in nature.

From a theoretical viewpoint, the macroscopic electric field in the polar surface can be cancelled either by surface reconstruction or by reduction of effective charge of the surface layer. According to Wolf, the polar (111) surface of rocksalt crystals would be stable under an octapole termination, which leads to formation of a 2$`\times `$2 reconstruction on the (111) surface (Fig. 2-b). On the other hand, Tsukada and Hoshino pointed out that the change in charge state of surface atoms could stabilize the (111) surface of rocksalt structure compounds. If the charge of top atoms is reduced to half of the bulk atoms, the macroscopic electric field would not appear.

Under these backgrounds, several attempts to grow the polar surface of rocksalt structure compounds have been performed so far. Metal oxides have been extensively studied to create crystals or thin films having (111) surfaces. For a NiO(111) surface, p(2$`\times `$2) reconstruction was found, and the surface was determined to be a single Ni termination with double steps by the measurement of grazing-incidence X-ray diffraction. The unrecostrcuted surface is stabilized by some adsorbates. Langell and Berrie reported that the NiO(111) surface is stabilized by OH adsorbates and that desorption of the hydroxyl changes the surface to a thermodynamically more stable NiO(100). Stabilization by impurities was found also for Pb and Si on NiO(111).

In contrast with the surface of single crystals, the electrostatic energy does not reach exorbitant values for the small thickness of a film (Fig. 2-c). When the thin film is grown on a metal substrate, further reduction of energy could come from the image charge induced in the supporting metal surface. Ventrice et al. have studied the growth of NiO on Au(111) by LEED and STM, and obtained a 6 ML thick stable p(2$`\times `$2) reconstructed NiO(111) film. For a FeO(111) thin film, interesting results were reported by Koike and Furukawa. They obtained the FeO (111) p(2$`\times `$2) surface by oxidizing Fe(110) and measured the polarization of secondary electrons emitted from the FeO (111) surface. Ferromagnetic ordering was found for the FeO (111) p(2$`\times `$2) surface, although FeO itself was an antiferromagnetic material with the Neel temperature of 198 K. They thought that the ferromagnetic ordering comes from the reconstruction of a polar FeO (111) surface to reduce the large electrostatic surface energy. On the other hand, the relaxation, that is, the decrease in interlayer spacing, was observed for the FeO (111) grown on Pt (111). In this case, the layer by layer growth was limited to a maximum of 2.5 ML.

Thus, although the (111) surface of rocksalt structure oxides has been studied by many groups, there have been few studies on unreconstructed adsorbate-free (111) surfaces of a rocksalt structure until now. Furthermore, the electronic structure of the polar surface has not been clear. In the preset study, we have examined growth of a MgO thin film on Ag(111) by supplying Mg and O, alternately. MgO has a rocksalt structure with a lattice constant of 4.21 Å, while Ag has a fcc structure with a lattice constant of 4.09 Å. As the lattice misfit of MgO to Ag is only -2.9 %, the first Mg layer is expected to become a template for the growth of the MgO film along the  direction. The commensurate bonding between Mg and Ag atoms at the interface might help alternate stacking of Mg and O layers, leading to a flat (111) surface. Actually, we have observed rather streak RHEED pattern, implying formation of a flat (111) MgO surface. We have characterized the grown surface by using EELS and ultraviolet photoemission spectroscopy (UPS) and discussed stabilization mechanism of the polar surface.

## II EXPERIMENTAL

The experiments were performed in a custom-designed ultrahigh-vacuum (UHV) system with a base pressure of 2$`\times `$10-8 Pa. A mechanically and electrochemically polished Ag(111) surface was cleaned by repeated cycles of Ar+ sputtering and annealing at 900 K. After repeated preparation cycles, a sharp RHEED pattern was observed, and no contamination was detected by AES. A MgO film was grown by alternate adsorption of Mg and O2 on Ag(111) at substrate temperature of 300 K. First, 1 ML (2 Å) Mg was deposited by evaporating high purity (99.98 %) Mg onto Ag(111). The growth rate was monitored using a quartz crystal oscillator. The Mg film was dosed with 10 L O2. Real-time observation of the crystallinity and orientation of films was done by RHEED. Surface compositions and the electronic structure of the grown films were investigated in situ by AES, EELS, and UPS. The analyses were performed with a double pass CMA (PHI 15-255G) in the analysis chamber. Adoption of a pulse counting detector reduced the probing electron current down to 1 nA, which reduced surface damage as much as possible.

## III RESULTS

### A Heteroepitaxial growth of MgO on Ag(111)

Figure 3 shows the Auger spectrum of the 10 ML thick MgO film grown on Ag(111). Disappearance of the Ag LMM Auger peak (351 eV) indicated that the grown MgO film covered the Ag(111) substrate completely because the inelastic mean free path of 351 eV electron was around 1 nm. The chemical state of Mg atoms of the MgO film was known by the energy of the Mg LMM Auger peak, because the energy of the Mg LMM Auger peak is 45 eV and 32 eV for metallic Mg and MgO, respectively. The energy of the Mg LMM Auger peak was 32 eV for the grown MgO film, showing that Mg was completely oxidized. The stoichiometry of the MgO film was also examined by the ratio of the Mg KLL Auger peak intensity to the O KLL Auger peak intensity. Considering the Auger electron emission probabilities, the ratio of the amount of Mg to that of O was 1 : 0.92$`\pm `$0.18, supporting the idea that a stoichiometric MgO film was grown on Ag(111).

Figure 4 shows the typical sequence of RHEED pattern during the growth at a substrate temperature of 300 K. The incident electron beam was parallel to the \[1$`\overline{1}`$0\] azimuth of the substrate. The result of RHEED patterns indicated that the MgO film grew heteroepitaxially on Ag(111). The epitaxial orientation of the MgO film was determined to be (111)MgO//(111)Ag and \[1$`\overline{1}`$0\]MgO//\[1$`\overline{1}`$0\]Ag. The half-order streaks did not appear during the growth, showing that the (1$`\times `$1) unreconstructed MgO(111) film was grown on Ag(111). Streaks in RHEED patterns indicated that a rather flat (111) surface could be obtained. The RHEED pattern became blurred with increasing film thickness, suggesting that a thick MgO(111) film was unstable.

The in-plane lattice constant of the MgO(111) film was calculated from the spacing between streaks in the RHEED pattern. For the 10 ML thick MgO(111)/Ag(111), the in-plane lattice constant was determined to be 3.28$`\pm `$0.03 Å, which was +10 % larger than that of bulk one (2.97 Å). The in-plane lattice constant was 3.25$`\pm `$0.03 Å for the 2 ML thick film and did not change with film thickness, indicating that the expansion was uniform throughout the epitaxial layer. The increase of the in-plane lattice constant leads to reduction of the surface electron density, which might help a decrease in the electrostatic energy of the MgO(111) film.

### B Stability of the (1$`\times `$1) MgO(111) film on Ag(111)

Having established the existence of the (1$`\times `$1) unreconstructed MgO(111) film on Ag(111), we should discuss why this structure is stable in spite of the previously mentioned arguments on the instability of a polar surface of ionic crystals. We think that the stabilization of the (1$`\times `$1) unreconstructed MgO(111) film can be explained in terms of interaction between Ag and Mg atoms at the interface and charge state of surface atoms.

In the previous study, a single crystalline MgO(100) film grew heteroepitaxially on Ag(100) by evaporating MgO congruently from an electron beam evaporator(EB). On Ag(111), however, a MgO film could not grow heteroepitaxially by EB. Heteroepitaxial growth of a single crystalline MgO film was achieved by alternate adsorption of Mg and O2 on Ag(111). These facts indicate that the strong interaction between Mg and Ag atoms plays a decisive role in stabilizing the (1$`\times `$1) MgO(111) film. The free energy of the interface ($`E_{inter}`$) is smaller for the case with coherent Mg-Ag bonds as compared with the case, in which fine particles with thermodynamically more stable {100} faces grow on Ag(111) without coherent Mg-Ag bonds. The structure of the grown film is determined to minimize the sum of the electrostatic energy of the film ($`E_{electro}`$) and $`E_{inter}`$. Here, we should notice that $`E_{electro}`$ of a (111) film does not reach exorbitant values for ultrathin films (see Fig. 2-c), because $`E_{electro}`$ is proportional to the film thickness. Therefore, the gain in interaction between Mg and Ag overcomes the disadvantage of the electrostatic energy, and the (1$`\times `$1) unreconstructed MgO(111) film is grown on Ag(111). In addition to the strong interaction with Mg atoms, the Ag substrate plays another important role for stabilization of the MgO(111) film. The MgO film was prepared by alternate adsorption of Mg and O2. Therefore, the film is unstable due to breakdown of charge neutrality (Fig. 2-d), when the film thickness is odd number (Mg top). On the metal substrate, however, the image charge is induced in a metal, helping stabilization of the film with odd layers.

Besides the strong interaction between Mg and Ag atoms at the interface, the stabilization of MgO(111) can be explained in terms of charge state. Due to the reduction of the coordination number, Madelung potential largely decreases for the surface O2- on the (111) surface. The binding energy of the 2p band, thus, decreases for the surface O2-, and the Fermi level might be located in the 2p band. In such a case, electrons flow out of the upper valence states of the surface O2-, and the charge of the surface O2- is reduced. The charge of the Mg2+ at the interface would also decrease to keep the film neutral (Fig. 2-e). The decrease in the charge reduces a macroscopic electric field, and the instability of the film should be remedied. That is the second point to stabilize the polar MgO(111) film on Ag(111). Here we should notice that the reduction in charge of surface O2- leads to reduction of Madelung potential. Therefore, the amount of the charge transfer is determined under the condition, in which both Madelung potential and the charge of surface O2- are self consistent. Tsukada and Hoshino studied the O-terminated MgO(111) surface by DV-X$`\alpha `$ calculations, and revealed that the charge of the top and bottom atoms is reduced to half of the bulk atoms. When the charge of the surface atom is reduced to half of the bulk atoms, the macroscopic electric field would not appear, irrespective of the film thickness. In the present study, a thick MgO(111) film could not be grown on Ag(111), implying that the charge of the topmost O atoms did not reduce to half of the bulk one.

Finally, we would discuss the surface reconstruction of the MgO(111) film. Since the surface energy of a {100} face is lower than that of a {111} face in rocksalt structure compounds, {100} faces are expected to appear for growth normal to the (111) face. In case of NaCl/GaAs(111)A, triangular pyramids with three exposed {100} faces grow on the substrate. As discussed in Introduction, a p(2$`\times `$2) reconstructed NiO(111) film is grown on Au(111), and the p(2$`\times `$2) structure corresponds to the smallest triangular pyramids surrounded by {100} faces. Here, we should notice that both films were grown by congruent evaporation of NiO or NaCl at the substrate temperatures higher than 420 K. Therefore, the grown film forms in a thermodynamically stable structure. In the present case, however, we have grown a metastable flat MgO(111) film by alternate adsorption of Mg and O2 at 300 K. Once a flat (111) film is formed on Ag(111), activation energy is needed to change the metastable (111) film into stable triangular pyramids surrounded by {100} faces. Therefore, the flat unreconstructed (111) MgO film remains in a metastable form on Ag(111).

### C Electronic structure of MgO(111)

As we could prepare the unreconstructed polar MgO(111) surface, the electronic structure of this surface was investigated by EELS and UPS, comparing with that of a non-polar MgO(100) surface. Absence of Ag LMM Auger peak in the 10 ML MgO(111)/Ag(111) assures that EELS and UPS probe the MgO film with influences of the substrate negligible, because their probing depth is smaller than or equal to that of AES.

Figure 5 shows UPS spectrum for the MgO(111) surface measured with a He I (21.2 eV) source. For comparison, the spectrum of the 10 ML thick MgO(100) film on Ag(100) is included in the same figure. The main features for MgO(111) were almost similar to those for MgO(100) except around the Fermi level. For MgO(100), there were no density of states (DOS) near the Fermi level, while an appreciable DOS appeared in the region from 2 eV to the Fermi level for MgO(111). These finite DOS at the Fermi level did not originate from metallic Mg, since component of metallic Mg was not observed in AES. Furthermore, we have revealed that the finite DOS at the Fermi level disappeared even for the incompletely oxidized film in the previous study. Therefore, the finite DOS at Fermi level did not originate from metallic Mg or the incompletely oxidized Mg, the surface states of MgO(111). The UPS results suggested that the polar MgO(111) surface was a metallic surface.

Figure 6 shows the EELS of the MgO(111) surface measured with primary electron energy of 60 eV. For comparison, spectrum of the MgO(100) surface is also shown in the figure. The structure above 7 eV for MgO(111) was almost similar to that of MgO(100), showing that a stoichiometric MgO film grew on Ag(111). On the other hand, difference appears in the structure below 5 eV between MgO(100) and (111). In contrast with MgO(100), the clear band gap region was not observed for MgO(111). In addition, the tail of the elastic peak became larger and a new peak (p1) appeared at 1.2 eV, indicating existence of the inelastically scattered electrons with small loss-energy. Since the AES signal indicative of Mg metal was not observed at all, the tail or the new peak did not originate from Mg aggregates in the MgO film.

Due to the low growth temperature and the significant mismatch between the substrate and the oxide lattice, the film may contain a significant concentration of defects. Therefore, the spectroscopic features could be due to emission from defect states in the MgO(111) film. Defective MgO surfaces have been studied by a variety of techniques such as UPS, EELS, theoretical calculation, etc. In EELS of defective MgO surfaces, sharp peaks are observed at 2.3 eV and 5 eV, which are attributed to Mg vacancy (V center) and O vacancy (F center), respectively. Sharp structures are also observed for MgO powders by diffuse reflectance spectra, and these structures at 5.8 and 4.6 eV have been attributed to ions with four and three-fold ligand coordination. On the other hand, such sharp peaks were not observed in the EELS of MgO(111) below 5 eV. Recently, we have revealed that metal induced gas states(MIGS) were formed at the insulator-metal interface. Since the new states of the MgO(111) film were observed independently of the film thickness, the states did not originate from MIGS. Therefore, the structure below 5 eV could be assigned to the excitation of electrons in the unfilled band derived from O 2p, or the excitations of newly appearing electronic states on the (111) surface.

The UPS and EELS results suggested that the MgO(111) surface was not an insulating surface but a semiconducting or metallic surface. This electronic structure of the MgO(111) surface can be explained in terms of Madelung potential as discussed in the previous section. Compared to the binding energy of isolated ions, the binding energy of Mg2+ 3s orbital decreases, and the binding energy of O2- 2p orbital increases by the Madelung potential in a crystal phase. Because of the decrease in the coordination number, the Madelung potential largely decreases for the (111) surface, compared with the (100) surface. Therefore, the band gap is reduced and the MgO(111) surface changes into a semiconductor or metallic surface. However, a possibility of certain new electronic state characteristic of the (111) surface cannot be excluded.

## IV CONCLUSIONS

Heteroepitaxial growth of MgO on Ag(111) has been investigated by RHEED, AES, UPS and EELS. The flat unreconstructed polar MgO(111) surface was obtained on Ag(111) by alternate adsorption of Mg and O2 at substrate temperature of 300 K. The EELS and UPS of the (111) surface were different from those of the (100) surface. For MgO(111), a clear band gap region was not observed in EELS, while finite density of states appeared in UPS, implying that the MgO(111) surface was a semiconducting or metallic surface.

###### Acknowledgements.

The authors are grateful for the financial support of the Grant-in-Aid from Ministry of Education, Culture, Sports, Science and Technology of Japan.


text:
# Magnetic bipolar transistor

## Electrical spin injection.

Our goal is to calculate how much spin polarization $`\alpha _c=s_c/N_{dc}`$ will accumulate in the collector in response to the nonequilibrium spin $`\alpha _e`$ and the equilibrium spin $`\alpha _{0b}`$. To that end we model the emitter-base junction as a forward biased magnetic p-n junction with equilibrium spin polarization $`\alpha _{0b}`$ in the $`p`$ region (base) and source spin polarization $`\alpha _c`$ in the $`n`$ region (emitter). Our theory of magnetic p-n junctions Fabian et al. (2002) determines the electron $`n_{be}`$ and spin $`s_{be}`$ densities at the $`be`$ region (see Fig. 1) at the depletion layer:

$`n_{be}=`$ $`n_{0b}e^{qV_{be}/k_BT}\left(1+\alpha _e\alpha _{0b}\right),`$ (2)
$`s_{be}=`$ $`n_{0b}e^{qV_{be}/k_BT}\left(\alpha _{0b}+\alpha _e\right).`$ (3)
The nonequilibrium electron density injected into the base depends on the product $`\alpha _e\alpha _{0b}`$, realizing the Silsbee-Johnson spin-charge coupling Silsbee (1980); Johnson and Silsbee (1985). If $`\alpha _e=0`$, Eq. (2) reduces to the standard Shockley’s equation Shockley (1950) for the nonequilibrium minority electron density in a biased p-n junction.

We next model the base-collector junction as another magnetic p-n junction. This junction is reverse biased, and has both the equilibrium spin polarization $`\alpha _{0b}`$ and the source spin density $`s_{be}`$ in the $`p`$ region (base). This is the case of a magnetic solar cell Žutić et al. (2001b), since the electron and the source spin densities in the $`p`$ region mimic the carrier and spin generation by light. For this case our theory Fabian et al. (2002) gives

$$s_c\gamma _1s_{be}=\gamma _1n_{0b}e^{qV_{be}/k_BT}\left(\alpha _{0b}+\alpha _e\right),$$
(4)
where $`\gamma _1\left(L_s/w_b\right)\mathrm{tanh}\left(w_c/L_s\right)`$. The accumulated spin polarization, which is the measure of the electrical spin injection efficiency, is $`\alpha _c=s_c/N_{dc}`$. Typically the spin diffusion length in the collector $`L_{sc}w_b`$, which means that $`\alpha _c`$ can be a considerable fraction (say, 10%) of $`\alpha _e`$ or $`\alpha _{0b}`$. What is interesting in Eq. (4) is the fact that $`\alpha _{0b}`$ plays the same role as $`\alpha _e`$ in the spin injection: the equilibrium spin can cause spin accumulation in the low injection limit, because it leads first to nonequilibrium spin $`s_{be}`$. This has no analog in magnetic diodes, where spin accumulation cannot result from the presence of just an equilibrium spin polarization.

## Spin control of current amplification.

When written in terms of $`n_{be}`$, the formulas for the currents $`j_e`$ and $`j_c`$ are the same as for the standard (nonmagnetic) bipolar transistors derived by Shockley Shockley (1950). After we write those formulas for the active forward regime, we substitute Eq. (2) for $`n_{be}`$ and obtain the dependence of the currents (and of $`\beta `$) on $`\alpha _e`$ and $`\alpha _{0b}`$.

The emitter current is

$$j_e=j_{gb}^n(n_{be}/n_{0b})+j_{ge}^p(p_e/p_{0e}),$$
(5)
where $`j_{gb}^n`$ is the electron generation current $`(qD_n/L_n)n_{0b}\mathrm{coth}(w_b/L_n)`$, the hole generation current $`j_{ge}^p`$ is $`(qD_p/L_p)p_{0e}\mathrm{coth}(w_e/L_p)`$, and the injected hole density in the emitter $`p_e=p_{0e}\mathrm{exp}(qV_{be}/k_BT)`$. The collector current comprises only electrons (Fig. 1):

$$j_c=j_{gb}^n(n_{be}/n_{0b})\mathrm{cosh}(w_b/L_n).$$
(6)
After evaluating $`j_b=j_ej_c`$ and substituting Eq. (2) for $`n_{be}`$, it is straightforward to show that in the narrow base limit ($`w_bL_n,L_s`$) the gain is

$$\beta =1/(\alpha _T^{}+\gamma ^{}),$$
(7)
where we use the standard transistor notation Tiwari (1992):

$`\alpha _T^{}`$ $`=`$ $`(w_b/L_n)^2/2,`$ (8)
$`\gamma ^{}`$ $`=`$ $`{\displaystyle \frac{N_bD_p}{N_eD_n}}{\displaystyle \frac{w_b}{L_p\mathrm{tanh}(w_e/L_p)}}{\displaystyle \frac{\sqrt{1\alpha _{0b}^2}}{\left(1+\alpha _e\alpha _{0b}\right)}}.`$ (9)
The factor $`\alpha _T^{}`$ determines how much electrons will recombine in the base, thus not reaching the collector. This factor is not affected by the presence of spin, and is the same as in the standard transistors. The factor $`\gamma ^{}`$ is related to the emitter injection efficiency, since it measures the proportion of the electron flow in the emitter current (where both electrons and holes contribute). This factor does depend on the spin. To get the maximum amplification, both $`\alpha _T^{}`$ and $`\gamma ^{}`$ need to be small. For the most efficient spin control of $`\beta `$, one needs $`\alpha _T^{}\gamma ^{}`$, the case of Si-based transistors which have slow carrier recombination. In this case

$$\beta (\alpha _e,\alpha _{0b})=\beta (\alpha _e=0,\alpha _{0b}=0)\times \frac{1+\alpha _e\alpha _{0b}}{\sqrt{1\alpha _{0b}^2}}.$$
(10)
The current amplification is affected by both $`\alpha _e`$ and $`\alpha _{0b}`$.

As an illustration we calculate $`\beta `$ as a function of $`\alpha _{0b}`$ for an MBT with $`\alpha _e=0.9`$ and with generic materials parameters, specified for a Si-like transistor. The nominal widths of the emitter, base, and collector are 2, 1.5, and 2 $`\mu `$m, respectively. The dopings are $`N_e=10^{17}`$, $`N_b=10^{16}`$, and $`N_c=10^{15}`$ cm-3. Electron (hole) diffusivities at room temperature are taken to be $`D_n=100`$ and $`D_p=10`$ cm2/s. The bias voltages are $`V_{be}=0.5`$ and $`V_{bc}=0`$ volts. The intrinsic carrier density $`n_i=10^{10}`$ cm-3 and the dielectric constant (needed to calculate the effective widths $`w`$) is 12. The carrier and spin diffusion lengths (note that Si has long recombination and spin relaxation times Lepine (1970)) are taken to be $`L_n=30`$ $`\mu `$m, $`L_p=L_s=10`$ $`\mu `$m. The calculated $`\beta `$ varies strongly with the spin, following closely the approximate $`\beta `$ given by Eq. (10). The amplification is largest (smallest) for the parallel (antiparallel) orientation of the source and equilibrium spins.

We conclude that spin can be injected through MBTs and that current amplification can be controlled by both the source and the equilibrium spin, making MBTs attractive for spintronic applications.

This work was supported by DARPA, the NSF-ECS, and the US ONR.