text:
# O adsorption and incipient oxidation of the Mg(0001) surface

## I INTRODUCTION

The process of oxidation of metal surfaces is of considerable fundamental scientific interest as well of paramount technological importance.classical\_oxide\_ref  Corrosion and passivation are just two examples of either destructive or useful processes linked directly to this phenomenon, which are well-known from everyday life. Oxidation of metal surfaces begins with dissociative chemisorption of oxygen on a clean surface and is followed by the formation of a film of metal oxide. In-between these events there are many elementary processes. The study of Al and Mg oxidation is particularly important because Al and Mg belong to the group of so-called simple metals and are considered model systems for studies of oxidation of transition metals. Both metals exhibit high reactivity with oxygen and oxidize rapidly. Aluminum oxides (alumina) form several different phases where the structure of some of them was only recently identified by combined ab initio density functional theory (DFT) and experimental studies.Yourd99  Magnesium oxide is known to eventually form crystals of the rock-salt structure,Nog  but can also experience complex reconstructions at partial coverages.Else01 ; Goon02  Both Mg and Al oxides play an important role in catalyst support and find many other useful applications.

The number of experimentalGoon02 ; NamDG81 ; Ghi81 ; Hay81 ; FloMa82 ; Thi89 ; Cron89 ; MitOCS98 ; Dri99 ; Lac94 ; Esa97  studies and theoreticalElse01 ; Bung97  treatments devoted to the oxidation of magnesium is relatively limited compared to the more extensively studied oxidation of aluminum. In this paper we present a systematic study of the initial oxygen incorporation and island formation by first-principles theory calculations and we use the results to interpret high-quality x-ray photo-electron diffraction (XPD) measurements for low O2 doses.

The existing picture of the Mg(0001) surface oxidation processes dates back to the early 1980’s. Based on the extensive low-energy electron diffraction (LEED), Auger electron spectroscopy, electron energy loss spectroscopy, and work function measurements Namba et al. proposed NamDG81  a four-stage model consisting of: (1) dissociative oxygen chemisorption at random sites followed by oxygen incorporation, (2) assembly of incorporated oxygen atoms into islands and lateral growth, (3) oxide formation from surface Mg atoms and subsurface O atoms which starts at O2 exposures around 2–3 L (Langmuir) or 0.5–1 monolayer of total coverage, and (4) oxide thickening. Several other experimental studies provided arguments for a three-stage model of oxidation that was proposed at approximately the same time. In this model FloMa82 ; Thi89  atomic oxygen is directly incorporated into magnesium (below the top Mg layer) right after on-surface dissociation. In the next step there is simultaneous formation of an oxide layer and a decrease in the oxygen on the surface.Thi89  Finally, there is oxide thickening and transformation into a rocksalt structure.

This model of immediate incorporation of oxygen gained support from the x-ray photo-electron spectroscopy study of Ghijsen et al.,Ghi81  and particularly from the measured sharp decrease in the work function of the Mg(0001) surface upon initial O exposure.Hay81  The formation of an O-($`1\times 1`$) subsurface layer up to monolayer coverage and its subsequent transformation into epitaxial oxide was also reported.Hay81  The positions of the subsurface oxygen atoms were deduced to be octahedrally coordinated where the oxygen has an environment similar to that in MgO.Thi89  This is in agreement with a more direct experimental evidence for oxygen incorporation in the octahedral interstitial sites of the first two interlayer spacings of Mg during the initial stages of magnesium oxidation,Cron89  and with later inelastic ion scattering experiments on a polycrystalline Mg surface.Lac94 ; Esa97  A very recent scanning tunneling microscopy (STM) study of the oxidation of Mg(0001) suggests that at low oxygen exposures (up to 2 L) the incorporated oxygen atoms form a single layer underneath the top layer of Mg.Goon02  The idea of immediately populated subsurface sites seems also to be supported by recent measurements by Mitrovic et al.,MitOCS98  however, they show that most of the oxygen (90%) remains over the surface and only a small fraction goes below the surface. Thus, the question regarding the most favored adsorption sites is still debated and requires further analysis.

On the theoretical side, the local density-functional theory (DFT) calculations of Bungaro et al.Bung97  show that oxygen is incorporated below the surface forming a $`(1\times 1)`$ subsurface lattice. A recent DFT-based lattice gas model studyElse01  of island formation in the early Mg(0001) oxidation stages has shown that at very low filling the oxygen atoms adsorb in the top-most subsurface (tetrahedral) B-sites in an ABAB… stacking of magnesium. At higher oxygen load also the deeper sites are filled, and the adsorbed oxygen atoms form dense clusters just below the surface and in further subsurface locations.Else01 ; Bung97  For larger clusters the oxygen-oxygen interaction effects drive the oxygen atoms further into the magnesium subsurface layers and bulk.

The rest of this paper is organized as follows. In Section II systematic DFT calculations in the range of coverages between 1/16 and 3 monolayers (ML) are presented and used for the study of the initial oxygen incorporation and island formation. The previous experimental and theory results are the starting base for our first-principles calculations of the very beginning of O adsorption and oxide formation at the Mg(0001) surface. In Section III we describe our high-quality XPD measurements for low O2 doses corresponding approximately to global coverages between 0.1 and 2 ML of atomic oxygen, and in Section IV our DFT results are used to interpret the experimental XPD data.

## II THEORY

The first-principles calculations are carried out using the plane-wave density-functional dacapo codeDACAPO  with the generalized gradient approximation (GGA) for the exchange-correlation energy functionalPer92  and with ultrasoft pseudopotentialsVan90  to represent the ionic cores. The clean Mg(0001) surface is modeled by periodic slabs consisting of six to seven magnesium layers separated by 21 Å of vacuum. The plane-wave basis set with 25 Ry energy cutoff is used. The O atoms are adsorbed on one side of the slab only and the electric field arising due to the asymmetry of the system is compensated for by a consistent dipole correction.NS92 ; Ben99  A $`4\times 4\times 1`$ mesh of Monkhorst-Pack special $`k`$-points for the $`4\times 4`$ atom surface unit cell and a Fermi-surface smearing of 0.2 eV are applied to the Brillouin-zone integrations. For smaller cells the number of k-points are increased accordingly up to a $`16\times 16\times 1`$ mesh for the $`1\times 1`$ surface unit cell. The positions of atoms in the three to four topmost magnesium layers, and of all the oxygen atoms, are fully optimized until the sum of the Hellmann-Feynman forces on all unconstrained atoms converges to less than 0.05 eV/Å. The forces acting on the ions in the unit cell are derived from the converged charge densities, and the atom dynamics is determined using a preconditioned quasi-Newton method based on the Broyden-Fletcher-Goldfarb-Shanno algorithm.NumR

The study involves very low coverages of oxygen, from one O atom in a large surface cell up to coverages of three monolayers. We define the coverage $`\mathrm{\Theta }`$ as the ratio of the number of adsorbed atoms to the number of atoms in an ideal substrate layer. The binding energy per adsorbate is calculated relative to the energy of an isolated, spin-polarized O2 molecule and the clean, relaxed Mg-surface.

### II.1 The clean Mg(0001) surface

Before studying the effect of oxygen chemisorption on the Mg(0001) surface, we determined the bulk and bare-metal-surface structures. The calculated lattice constants for hcp Mg is $`a=3.19`$ Å and $`c/a=1.64`$, which agree very well with the experimental values $`3.21`$ Å and 1.623 (Ref. mglattconst, ) and other GGA calculations.WaKi01  The bulk modulus, ignoring lattice vibrations, is found to be 34.1 GPa,Ziam01  in good agreement with the measured value 35.5 GPa (obtained from the elastic constants of Ref. mgbulkmod, ) and with the GGA calculations of Ref. WaKi01, .

The relaxations of the surface interlayer spacing with respect to the bulk spacing show 1.5% expansion, 0.4% contraction, and 1.3% contraction of the first, second and third interlayer distance, respectively. The expansion of the first interlayer distance compares well with the experimentalDavHRP92  value of 1.7%, and with other GGA calculations.WaKi01  The small contraction of the second interlayer spacing disagrees with results of experiment and other calculations which predict a small expansion. Since the magnitude of this relaxation is very small this discrepancy should not have any significant effect on the results for O-adsorption. The calculated work function is 3.72 eV for the relaxed surface and agrees well with the experimental value 3.84 eV. Mich77

### II.2 Low coverages ($`\mathrm{\Theta }0.5`$) of oxygen

Having obtained reliable results for the clean, relaxed Mg(0001) surface we calculate the binding of oxygen in on-surface and subsurface sites. The locations of the possible adsorption sites of Mg(0001) are sketched with crosses in Fig. 1. In our calculations the results for single or nearest-neighbor pairs of O atoms adsorbed at the $`4\times 4`$, $`2\times 2`$, and $`1\times 2`$ surface cells were used as 1/16, 1/8, 1/4, and 1/2 ML coverage data.

In Fig. 2 we present the calculated low-coverage binding energies with respect to the energy of an isolated O2 molecule. The adsorption of oxygen atoms on the Mg(0001) surface induces further changes to the surface structure and the electronic properties. However, the positional changes of the Mg atoms are in most cases not very large. For the $`2\times 2`$ supercell with a single O-atom slightly below the surface (site B2 in Fig. 1) or in first-sublayer adsorption sites (sites A4 or A5) the Mg atoms move less than 3% vertically and 0.6% laterally from their clean-surface positions. To estimate the sensitivity of the binding energy to positional relaxations of the atoms we also calculated an even distribution of four oxygen atoms in one of the tetrahedral subsurface sites (sketched as site B2 in Fig. 1) in the $`4\times 4`$ system, and compared the binding energy to the binding energy of one site-B2 atom in the $`2\times 2`$ system. The two configurations relax to slightly different atomic positions (but both converged within the force requirement) and we found the difference in binding energy per oxygen atom to be 5 meV. We therefore estimate the relaxation-imposed inaccuracy of the binding energies to be approximately 5 meV.

On-surface adsorption. It has previously been suggested that for low coverages the O atoms are preferentially adsorbed into on-surface sites,MitOCS98  as is the case for, e.g., oxygen adsorption on Al(111),KiLu01  or that on-surface adsorption sites exist although having lower binding energies than the lower-lying adsorption sites.Bung97  For the Al(111) surface an additional O is incorporated only when 1 ML of on-surface O is complete.KiLu01  A similar situation is found for a transition metal, the Ru(0001) surface.Reuter  For this reason on-surface adsorption was closely studied. In summary, however, we find the initial oxygen adsorption to be a subsurface process.

Of the on-surface adsorption sites the natural candidates are the C1 and B2top hollows (Fig. 1). However, our calculated low-coverage binding energies, plotted in Figure 2, show that for single atom adsorption at coverages $`\mathrm{\Theta }0.5`$ (corresponding to a single O atom in the $`4\times 4`$, $`2\times 2`$, and $`1\times 2`$ surface unit cells) the only stable on-surface adsorption site is the one denoted as C1, and the binding energy of this site is significantly less than that of any of the subsurface sites. We find no single-adsorbant energy minimum near the position where we would expect the on-surface B2top-site. Instead, when relaxing dilute systems with O in the B2top sites the surrounding Mg atoms move laterally by a small amount, enough to let the O-atoms sink into the subsurface B2 (tetrahedral) site. Only if a pair of O atoms are adsorbed as nearest neighbors, or at increased O loading ($`\mathrm{\Theta }1`$), can one of the O atoms be stabilized on the surface in the B2top site, albeit at smaller binding energy than in the slightly lower-lying B2-site. Thus, for low coverages our DFT results point to the C1 site as the only possible, but energetically unfavorable, on-surface adsorption site.

These results are in agreement with the theory results of Bungaro et al.Bung97  as to subsurface adsorption sites being more favorable than on-surface ones, but contradict their prediction of a stable on-surface B2top-site for low coverages. It should be noted, however, that for coverages $`\mathrm{\Theta }<1`$ Bungaro et al.Bung97  relaxed only the O atoms but not the Mg-lattice, and that their $`k`$-point sampling was limited to the $`\mathrm{\Gamma }`$ point only. Thus our results underline the importance of relaxing the substrate.

The subsurface adsorption. By extensive searches for subsurface adsorption sites we find that within the two top Mg-layers the tetrahedral sites (sites B2, A4, and A5) have higher binding energies than the on-surface C1 and the subsurface octahedral C3-site for all coverages $`\mathrm{\Theta }0.5`$ (compare Fig. 2 and Table 1). Thus, the first O atom adsorbed in the oxidation process binds in a subsurface site of type B2, its binding energy being only 0.09 eV and 0.15 eV higher (Fig. 2) than that of O in subsurface site A4 and A5, respectively. The preference for the B2 site agrees with the results of Ref. Bung97, , but we find a binding energy about 0.6 eV larger. In addition, they find that the on-surface B2top site is much more favored than the subsurface A4 site, in clear disagreement with our results. These difference are likely due to the lack of Mg-atoms relaxation and the inferior $`k`$-points sampling used in that study.

In order to check the possibility of formation of ionic complexes (which would represent a seed for further oxide formation) consisting of the Mg atoms of the topmost layer and subsurface and on-surface oxygen, we calculated the most preferred sites for an additional O atom in the presence of an already adsorbed one. In studying possible combinations of the oxygen adsorption sites we concentrated on the ones including the B2 sites. As is evident from Table 1, even at low coverages with one O-atom isolated from other O-atoms the addition of another O atom leads to an increased O-population in the B2 sites. The two (neighboring) atoms in the B2 sites are by $``$0.03–0.04 eV more favored than one atom in the B2 site and another one in either the A4 or A5 site.

### II.3 Higher oxygen coverage ($`1\mathrm{\Theta }3`$ ML)

For 1 ML coverage with occupation in only one kind of site the energetic ordering of the most stable sites is different than for lower coverages (Fig. 2 and Table 2). The B2-type sites become less favorable than the A4 site and the actual B2 site is shifted to the B2top position. However, allowing for simultaneous occupancy of the B2 and B2top sites at 1 ML in the $`1\times 2`$ unit cell (Table 2) we find that this combined B-site structure is by 0.06 eV more favored than the single A4 occupancy shown in Fig. 2.

As more oxygen is adsorbed into the Mg surface the possible adsorbate configurations grow in number, partly due to simple combinatorics, and partly because new adsorption sites appear within an increasingly distorted Mg background lattice. This is already apparent for nearest-neighbor adsorption in the low-coverage region (Table 1) where a slight vertical distortion of the top Mg-layer in some cases leads to adsorption in the B2top site and in a ‘mixed’ B2-B2top position. In Table 2 we name the 1 ML configurations by the adsorbate position, but for higher coverages the notation is merely an approximate description of the structure, and we refer the reader to Fig. 3 for a sketch of the multitude of similar but not identical configurations at coverages 2–3 MLs.

The stable structures found in the 2 and 3 ML coverage region fall in two general classes: layered structures with oxygen and magnesium in buckled or flat layers having a hexagonal structure \[Fig. 3 (a)–(c), and (e)\], and rock-salt structures \[Fig. 3 (d)\] on top of clean magnesium, strained with respect to the clean-surface lattice constant. The buckled structures are labeled according to the approximate positions of the oxygen atoms (comparing to the low-coverage adsorption sites of Fig. 1). The buckled structures involve subsurface oxygen only. The vertical separation of O and Mg within a buckled layer is approximately 0.6 Å, and the layers can exist both on top of the clean Mg(0001) \[Fig. 3(c)\] and as subsurface layers \[Fig. 3(a) and (b)\], and at 2 or 3 ML coverage. Thus one can discriminate the structures buckled up and down, depending on the oxygen atoms staggering over or below the Mg layer, respectively. The flat structures \[Fig. 3(e)\] have vertical O-Mg separation 0.01–0.06 Å, and exist for both 2 and 3 ML coverage, but do not exist as subsurface layers. In the flat structures the top-most oxygen atom sits slightly above the top Mg layer, corresponding approximately to the position of site B2top in Fig. 1. The binding energies of the 1–3 ML structures are given in Table 2. Note that for several other configurations with 2 ML or 3 ML coverage of the tetrahedral adsorption sites no energy minimum is found. This includes any subsurface flat phase, phases of mixed flat and buckled layers, and mixed buckle-up and buckle-down layers. Further, we tested structures that include an octahedral (C3, C6, C9) site in the 2 ML and 3 ML coverage and found binding energies per atom at least 0.3–0.4 eV smaller than the energetically best structures at the same coverage, similar to the results of the low coverage results (Fig. 2).

The preference for subsurface adsorption at the Mg(0001) surface was originally suggested by a sharp work function decrease upon oxidation.Hay81  Previous local-density-approximation DFT calculationsBung97  have shown that only oxygen incorporation into B2 sites lowers the work function and for one monolayer load this lowering was found to be 0.3 eV, which is much smaller than measured experimentally. We do not find 1 ML adsorption in the B2 sites, but for the neighboring B2top sites our results presented in Table 2 confirm the lowering of the work function. However, this lowering is three times larger than that reported in Ref. Bung97,  and we note that our B2top sites are slightly above and nearly coplanar with the topmost Mg layer. This larger decrease of the work function upon oxygen adsorption agrees with recent measurements.MitOCS98

It is generally expected that an adlayer of negative O ions will increase the metal-surface dipole layer and the work function. For example, this is clearly the case for oxygen on the Al(111) surface. KiLu01  Our results suggest that 1 ML of on-surface (B2top) oxygen on the Mg(0001) surface leads to a decrease in the work function. This decrease for on-surface O is consistent with experimental observations.MitOCS98  On the other hand, for 1 ML coverage the B2top sites are not the energetically most favorable. We conjecture that the following scenario takes place. At low coverage the only on-surface sites are the fcc hollows (C1), but also the lower-lying subsurface B2 sites are exposed to O atoms arriving at the surface. Thus an O-atom that arrives to the surface will find the B2 site having a higher adsorption energy than the C1 site and will immediately diffuse to the former. The B2-sites are most favorable up to 0.5 ML coverage (Fig. 2). With growing O-coverage the energetic situation changes and the deeper tetrahedral sites become more favorable. However, there is an energy barrier to overcome in order to get into these lower-lying A4 and A5 sites so for 1 ML coverage, instead of moving further into the surface, the oxygen atoms stay in the neighborhood of the B2 sites, now shifted into the B2top position. The B2top site is energetically slightly less favored than the tetrahedral subsurface sites (by about 0.06 eV) but requires no crossing of significant energy barriers. With an increased oxygen load (around 2 ML) the energy barrier is lowered and some of the oxygen atoms can move into the A5 sites. Again, the occupation of the B2top and A5 sites (the “Flat 2ML” structure) is not the most favorable energetically (Table 2) but presumably this metastable state is most easily available for the O atoms. The work function change $`\mathrm{\Delta }\mathrm{\Phi }`$ for this configuration is similar to $`\mathrm{\Delta }\mathrm{\Phi }`$ for the energetically favorable A4+B7 subsurface structure. With a further increase of O load (3 ML) the energy barriers that separate the lower lying sites are lowered and oxygen populates the B8 sites. It seems that this scenario is supported by our XPD experiments presented in Sec. III.

Island formation. In a previous theory study based on a lattice-gas-modelElse01  it was shown that oxygen atoms that are adsorbed in the Mg(0001) surface form dense clusters in tetragonal subsurface sites. For example, at a dosing corresponding to a global coverage of 1 ML, the oxygen atoms were found to be distributed, on average, with 55% in site B7, 30% in site A4, and 15% in site A10 and with local coverages of 2–3 MLs within the clusters.

In the present, more detailed investigation we find that past the adsorption of the first, isolated O-atoms in the oxidation process the O-atoms show a tendency for island formation by pairing up and showing higher binding energies in combinations with subsurface sites. An inspection of the binding energies presented in Fig. 2 and in Table 1 shows that the higher the coverage, the larger the binding energy of the lower-lying sites (sites A4 and A5). Our results show that in general the binding energy is higher in more close-packed systems, thus showing that clusters or islands of O atoms are preferred to isolated or pairs of O atoms. In fact, for 1 ML coverage, having all O atoms in site A4 is preferred by a small amount to all other subsurface sites and to the on-surface B2top site. This agrees well with the previous lattice-gas-model study.Else01  With coverage increasing above 0.5 ML, the O-atoms from the B2 sites are pushed up into the surface B2top sites. Thus they end up, at 1 ML load, in the B2top-sites located almost coplanarly with the Mg atoms of the topmost layer, with the binding energy comparable to that of O in sites A4 and A5.

The preference for the formation of subsurface islands in tetrahedral A4 sites or lower agrees both with the lattice-gas model simulations and with the results of Bungaro et al. Bung97  but disagrees with the results of Thiry et al. Thi89  who found adsorption in octahedral sites. A new feature that was not observed by Bungaro et al., and which is probably connected with a more careful treatment of the lattice relaxation effect in our work, is the appearance of the on-surface O-islands for 1 ML coverage, with binding energy slightly lower than in subsurface islands of A4-type. Also worth noting are the relatively small binding energy differences for oxygen in different sites for $`\mathrm{\Theta }=1`$.

The tendency to form islands continues with further increase of coverage, showing some quenching in the adsorption energy at 3 ML coverage. A similar behavior was found in the lattice-gas model simulations. Else01

Layered structures. As the coverage is increased beyond a single monolayer the layered structures of Fig. 3(a)–(c) and (e) appear, along with a rocksalt surface structure, Fig. 3(d). The layered structures can be described as a Mg(0001) surface with slightly changed atomic positions to accommodate O atoms within the top layers of the Mg surface, thereby creating a layered structure. Each of the layers has a honeycomb structure with O at every second vertex and Mg at the remaining vertices. The layers are stacked in an AA stacking sequence, with O on top of Mg atoms and vice versa at a layer separation of approximately 2.3–2.7 Å.

Again, we find (Fig. 3) that subsurface layered structures are preferred. Comparing the adsorption binding energies listed in Fig. 3 it is evident that with increased O coverage, the difference in binding energies for the subsurface layered and rocksalt structures is diminishing, indicating that eventually, for still higher coverage, the rocksalt structure may become the most favorable one.

Directional bonding. Figure 4 displays the electron charge redistribution for three of the 2 ML structures. Whereas the rocksalt structure (Fig. 4(b)) has an almost isotropic charge redistribution around each atom, as is typical for ionic binding, the two layered structures show more complicated, anisotropic charge redistributions. The flat layered structure (Fig. 4(c)) shows most clearly the anisotropy, with more charge redistribution along the in-plane nearest-neighbor Mg-O lines than in other directions.

## III EXPERIMENTS

To supplement the DFT calculations and to provide a reference ground we have also carried out XPD measurements of the clean Mg(0001) surface and after low dosing of oxygen.

### III.1 Method

XPD has been chosen because of its chemical sensitivity and its sensitivity to local real space order. It is a powerful technique for surface structural investigations,Fadley  and it has been shown that full-hemispherical XPD patterns provide very direct information about the near-surface structure. At photoelectron kinetic energies above about 500 eV, the strongly anisotropic scattering by the ion cores leads to a forward focusing of the electron flux along the emitter-scatterer axis. The photoelectron angular distribution, therefore, is to a first approximation a forward-projected image of the atomic structure around the photoemitters. Analysis of the symmetry and positions of forward-focusing maxima thus permits a very straightforward structural interpretation of XPD data. Furthermore, detailed structural parameters can be determined by comparing the experimental XPD patterns to calculated ones. The relatively simple and efficient single-scattering cluster (SSC) formalismFadley  has proven adequate in most cases. The agreement between SSC calculations and experimental XPD pattern can be quantified using a reliability factor such as the R-factor $`R_{\text{MP}}`$ defined previously.RMPref ; RMPexplain

The experiments have been performed in the University of Fribourg’s VG ESCALAB Mk II spectrometer modified for motorized sequential angle-scanning data acquisition. Clean Mg(0001) surfaces have been prepared by cycles of sputtering (500 eV Ar+) and annealing (130 C). Before O2 exposure, the surface displayed a sharp $`(1\times 1)`$ LEED pattern with little background. After O2 exposure with the sample held at room temperature, the coverage was determined from the relative intensities of the O 1s and Mg 2p photoelectron peaks. Experimental O $`1s`$ XPD patterns were obtained after exposure of the Mg(0001) surface to 0.15, 0.7, 1.4 and 9.7 L of O2, corresponding to global O coverages of 0.1, 0.4, 1 and 1.7 ML, respectively.

### III.2 Oxygen adsorbed at Mg(0001)

The Mg(0001) surface with various coverages of oxygen atoms was studied with LEED and XPD. Experimental XPD patterns for O2 doses of 0.15 and 9.7 L are shown in Fig. 5. The most prominent feature of the 0.15 L pattern is a strong intensity maximum at normal emission. This observation is a clear and direct evidence for O adsorption below the topmost Mg layer: Only an atom located directly above the oxygen photoelectron emitter can give rise to a forward-focusing maximum at normal emission, i.e., in the center of the plot. A further conclusion can immediately be drawn from the evolution of the XPD pattern with O2 exposure. Apart from differences due to counting statistics, the experimental XPD patterns are strikingly similar, not only regarding the intensity maximum at normal emission but regarding all the prominent diffraction features. It must therefore be concluded that the local atomic geometry is the same over the entire exposure range.

The same conclusion is obtained from the evolution of the LEED pattern with increasing O2 exposure: The LEED pattern stays ($`1\times 1`$), and only the background rises slightly, which indicates an increasing amount of disorder with increasing coverage. Since no superstructures are observed, the presence of any other islands than islands of ($`1\times 1`$) periodicity that are larger than about 50–100 Å can therefore be excluded.

In other words, both the LEED and XPD measurements indicate that locally the geometry is the same irrespective of the global coverage (at least below 1.7 ML). Already at 0.1 ML global coverage this local structure thus contains at least 1 ML oxygen, and with increasing coverage it will grow in domain size without significant changes in geometrical structure. Therefore, determining the local geometrical structure at 1ML local coverage will determine the geometries even at the lowest global coverages accessible to the experiments.

In the further analysis we compare simulated SSC diffraction pattern to experimental XPD images of the Mg(0001) surface both with dosing 1.4 L (approx. 1 ML global coverage) and 9.7 L (approx. 1.7 ML global coverage).

## IV Comparison between experiment and theory

That the O/Mg(0001) system does not simply form a MgO(111) rocksalt structure is clearly seen by comparing a simulated diffraction pattern for MgO(111) \[Fig. 5(c)\] to the experimental 9.7 L XPD pattern shown in Fig. 5(b): The simulation does not reproduce any of the diffraction features satisfactorily. Therefore, the most relevant atomic structures, as determined by the DFT calculations, were used to simulate XPD patterns by means of SSC calculations. The resulting SSC calculations were compared to the experimental 1.4 L and 9.7 L XPD patterns, the latter of which is shown in Fig. 5(b).

None of the SSC calculations for an oxygen atom in an isolated site gave satisfactory agreement with experiment. The calculations, however, confirm the conclusion drawn from the experimental patterns that one oxygen atom is adsorbed directly below a Mg atom: Only the simulations for oxygen in the A4 and A5 sites exhibit the characteristic intensity maximum at normal emission that is seen in experiment \[Fig. 5(a),(b)\]. The fact that the experimental XPD patterns do not agree with simulated diffraction patterns of disperse, low local coverages even for low overall coverage in the experiment is consistent with the LEED measurements discussed above, indicating again that the oxygen will form island with a local coverage of 1 ML or more. Accordingly, further SSC calculations concentrated on the atomic positions obtained from DFT calculations using a $`1\times 1`$ surface cell.

Among the simulated diffraction patterns for 1 ML local coverage only the site A4 and the site A5 SSC calculations exhibit a central maximum. In contrast, the SSC calculations for each of the C3, B2top, and C1 sites (leaving all other sites empty) do not show such maximum. This means that the comparison between SSC calculations and experimental XPD patterns thus disqualifies the site B2top-only occupation at 1 ML local coverage, in agreement with the theoretical binding energy being smaller for B2top than for A4 at 1 ML coverage. Along with the binding energy obtained from the DFT calculations the comparison to the XPD patterns also rules out octahedral subsurface (site C3) occupation, in disagreement with the results of Thiry et al. Thi89

Consequently, SSC calculations were also performed for the atomic positions obtained from DFT calculations considering 2 and 3 ML coverage. In Table 2 we summarize the 1 ML, 2 ML and 3 ML local coverage SSC simulations \[one to three O atoms per $`(1\times 1)`$ cell\] by listing their reliability factors with respect to the 1.4 L and 9.7 L exposure XPD data. Best agreements between experiment and SSC calculations (lowest reliability factors) are found for the layered structures “B2mix & A4” and “B2 & A5” at 1 ML local coverage, “Flat 2ML” \[Fig. 5(d)\] and “B2 + A5” at 2 ML local coverage, and “Flat 3ML” and “B2 + A5 + B8” for 3 ML local coverage. These results are consistent with the work function changes discussed above.

The general trend is that the rocksalt structures (2 ML and 3 ML) do not agree very well with the experimental XPD patterns (high values of $`R_{\text{MP}}`$), whereas some of the layered structures agree better. In particular we notice that among the 2 ML structures the flat structure (“Flat 2ML”) has the best agreement with the 1.4 L experiment, whereas one of the buckled layer structures (“B2 + A5”) fits best to the 9.7 L experiment. This indicates that as the global coverage is increased by higher O2 dosing, more of the surface becomes covered with buckled Mg and O layers and less with the flat surface.

For the 1.4 L dosing experiment the flat surface is favored, in disagreement with the DFT results. 1.4 L dosing corresponds to $``$1 ML of global coverage, but since the distribution of oxygen is nonuniform due to the tendency of the atoms to form islands the experimentally observed surface is patched. For small surface oxide patches of 2 MLs (or more) local oxygen coverage the strain can get released over most of the island and thus, experimentally, mostly the flat structure is seen at low dosing. This gives a low $`R_{\text{MP}}`$-value for the flat structure, but still keeps the order of the buckled and rocksalt structure found in DFT calculations: buckled is also at 1.4 L more favored experimentally than rocksalt. Thus, the relevant dosing (among the measured ones) to use for comparing DFT and the XPD experimental structures is 9.7 L. This structure probably does not fully cover the surface (if clusters are at least 2 ML dense), but certainly does so to a larger extent than for dosing 1.4 L. At dosing 9.7 L the DFT calculations and the XPD experiment agree on the buckled surface.

Based on both the calculated and experimental structures discussed above we conclude that at relatively low dosage (corresponding to 2–3 ML coverage) O/Mg(0001) forms the layered oxide structure. The rocksalt structure typical of MgO starts to grow only at higher O dosing.

## V Summary and CONCLUSIONS

We have performed extended first principles calculations of oxygen adsorption and of the initial stages of Mg(0001) oxidation. A variety of configurations and a wide range of coverages were considered in order to determine the most stable structures. These were compared with x-ray photoelectron diffraction experiment and simulations. At low coverages ($`\mathrm{\Theta }0.5`$ ML) both our DFT calculations and experiment show that oxygen adsorbs in subsurface sites. Our DFT calculations show that the first O atom chemisorbed in the oxidation process binds in a subsurface tetrahedral site of the B2 type. The importance of substrate lattice relaxation in accurate determination of the most stable sites is demonstrated. At higher coverage O adsorbed in subsurface tetrahedral sites shows a tendency to form subsurface islands which results in an increased binding energy. For the 2 ML (3 ML) coverages we find some rather unanticipated surface oxide structures, consisting of two (three) mixed oxygen-magnesium layers on top of an almost undistorted Mg(0001) surface. These layered oxide structures have hexagonal symmetry and can be flat or buckled. For 2–3 ML coverage the rocksalt structure is found to be unfavorable compared to the buckled layer structures. However, our DFT results show that the rocksalt structure may become energetically competitive at an increased coverage.

###### Acknowledgements.

This work was supported in part by the Swedish Foundation for Strategic Research (SSF), the Swedish Research Council (VR), The Swedish Foundation for International Cooperation in Research and Higher Education (STINT), the Carl Tryggers Foundation, and the Polish State Committee for Scientific Research (KBN), project 5 P03B 066 21. R.F. would like to thank P. Aebi for continuous support and help with the XPD experiments. The allocation of computer time at the UNICC facility at Chalmers and Göteborg University is gratefully acknowledged.


format: text

text:
# Do many-particle neutrino interactions cause a novel coherent effect?

## 1 Introduction

The flavor composition of a neutrino system may be modified as a result of the interactions between the neutrinos and a background medium. The efficiency of this process depends, in addition to the properties of the background itself (composition, density, temperature, etc.), on the coherent (or incoherent) character of the interaction.

Here we study the case of *coherent* neutrino conversion in a background of neutrinos. The interaction of a neutrino beam with a neutrino background has the peculiar feature that, due to momentum exchange, a neutrino from the background may be scattered into the beam and vice versa. This differs from the case of scattering on electrons and/or nucleons, and implies that the effect of coherent scattering on neutrinos can not be treated in analogy with that on ordinary matter.

In the literature, the effects of the neutrino-neutrino coherent scattering on neutrino conversion have traditionally been treated following the pioneering work by Pantaleone , Sigl & Raffelt , and McKellar & Thomson . In these papers, a *single-particle* evolution equation for each neutrino is derived, which has the form of the usual oscillation equation for neutrinos in matter with an extra term designed to capture the effect of neutrino–neutrino scattering. According to this equation, the effect of the interaction with the neutrino background is significant if the neutrino-neutrino potential is at least comparable to the potential due to matter (zero-temperature and thermal terms) and/or to the vacuum oscillation terms:

$$|n_\nu n_{\overline{\nu }}|\{\begin{array}{c}|n_e^{}n_{e^+}|\\ (n_e^{}+n_{e^+})E_\nu ^2/m_w^2\\ \mathrm{\Delta }m^2/(G_FE_\nu )\end{array}$$
(1)
Here $`n_\nu `$ $`(n_{\overline{\nu }})`$ is the (anti)neutrino number density, $`n_e^{}`$ $`(n_{e^+})`$ is the electron (positron) number density, $`E_\nu `$ is the neutrino energy, $`m_w`$ is the $`W`$-boson mass and $`G_F`$ is the Fermi constant.

The condition (1) may be satisfied in the early Universe close to the Big Bang Nucleosynthesis epoch (at temperature $`TO(1)`$ MeV) or just outside a supernova core (at $`rO(100\mathrm{km})`$). Indeed, in these environments the number density of neutrinos is comparable to, or larger than, that of electrons and nucleons. At the same time, incoherent scattering is negligible and neutrinos stream freely in the medium, affected by coherent scattering (refraction) only. The evolution of the neutrino flavor composition in both cases has been extensively studied using the one-particle equation. An important result  is the strong bound on the chemical potential of non-electron relic neutrinos that follows from studying neutrino oscillations in the early Universe, with important implications on BBN. Other important applications include hypothetical oscillations between active and sterile neutrino states that may have generated a lepton asymmetry in the early Universe  and possible effects of a sterile neutrino on the $`r`$-processes in a supernova .

There is an important theoretical issue at the foundation of all of the mentioned analyses that, if unresolved, casts doubt on their validity. The issue is the *existence and the range of validity* of the description of neutrino self-refraction in terms of a set of one-particle equations. Generally speaking, the interactions between neutrinos in the ensemble may be expected to create quantum correlations, or, in other words “entangled” states of many neutrinos which are not products of individual wavefunctions. *A priori*, a system in which such entanglement exists requires a many-particle description. The importance of this point was recognized early on , but was not pursued further in the following years.

Recently, the problem of possible existence and effects of quantum correlations has been investigated using two different approaches. The first, put forth in our recent paper , is to understand the flavor evolution of a many-neutrinos system as an effect of the interference of many elementary neutrino-neutrino scattering events. Using this construction, we argued that, in the limit of infinite number of neutrinos, the many-body description factorizes into one-particle equations. These equations coincide with those given in the previous literature up to terms which change the overall phases of neutrino states. While such terms can be potentially important in non-oscillation phenomena, they do not affect the flavor evolution of the system.

As a particular example, we considered an ensemble of neutrinos which are initially in flavor eigenstates. Following our analysis, we found that in this system, somewhat counterintuitively, no coherent flavor conversion takes place. This finding is in agreement with the standard one-particle description of .

The second approach, introduced by Bell et al. , is to consider the evolution of a system of interacting neutrino plane waves in a box. The idea is to study the properties of the many-body Hamiltonian of this system with the goal of determining the rate with which statistical equilibration is achieved. The analysis is applied to a system of neutrinos which are initially in flavor eigenstates. It is concluded that the time scale $`t_{eq}`$ of the flavor evolution of a given neutrino in this system is inversely proportional to the neutrino density and to the Fermi constant: $`t_{eq}(nG_F)^1`$, which is characteristic of coherent flavor conversion. This result contrasts with the prediction of the one-particle description, and therefore has been considered as an indication of its breakdown and of the existence of new coherent effects due to entanglement. Clearly, this would have profound implications on the applications of neutrino-induced neutrino conversion to supernova physics and cosmology.

Given the importance of the subject, in this paper we present a further study of neutrino flavor conversion in a neutrino background. The aim is to clarify the origin of the contrast between Ref.  and the other literature and give a definite answer to the question of the importance of entanglement and its effects.

Since the suggested new effect comes from simultaneous interaction between many neutrinos, we study, following , the Schrödinger equation of a system of many neutrino plane waves in a box. We consider the free streaming regime and neglect vacuum oscillation terms, which can be included by a straightforward generalization of our results. We show that the problem is equivalent to that of a system of spins, for which, in the case of constant spin-spin coupling, the equilibration time can be determined completely analytically. For neutrinos initially in flavor eigenstates, in the limit of the infinite number of neutrinos, this solution shows no coherent conversion, in agreement with our earlier results.

The paper is organized as follows. In Sect. 2 we lay the foundations of our analysis, by discussing the structure of the neutrino–neutrino interaction in the flavor space and in the real space. We also perform an elementary analysis of coherence in a system of several interacting neutrinos and use it to motivate the following study. In Sect. 3 the many body approaches are discussed and the relevant results of our previous paper are reviewed. In Sect. 4 we formulate the problem of interacting neutrino plane waves in a box, present the full analytical solution and give its derivation. A discussion and conclusions follow in Sect. 5 and 6.

## 2 General considerations

### 2.1 The flavor structure of the interaction

In the range of neutrino energies that are relevant to the applications ($`10^610^7`$ eV), the interaction between neutrinos is described by the low-energy neutral current (NC) Hamiltonian

$`H_{\mathrm{NC}}={\displaystyle \frac{G_F}{\sqrt{2}}}\left({\displaystyle \underset{a}{}}\overline{\nu }_a\gamma ^\mu \nu _a\right)\left({\displaystyle \underset{b}{}}\overline{\nu }_b\gamma _\mu \nu _b\right).`$ (2)
Higher order effects, such as those coming from the expansion of the propagator or particle production (bremsstrahlung) will be neglected in our analysis.

A property of this Hamiltonian that will play a very important role later is its invariance under the rotation of the $`SU(2)`$ flavor group. Because of this property, the interaction between pairs of neutrinos – viewed in flavor space – must be equivalent to the interaction between pairs of spins.

The one-particle equations of  indeed preserve the $`SU(2)`$ structure of the problem and have a form of the interaction between spins111In , this property of the equations has been used to give an elegant physical analogy designed to explain the “synchronized” oscillations of neutrinos in a neutrino background.. It would be logically inconsistent, however, to rely on this fact in the analysis, in which these equations are being tested. Therefore, below we present an explicit proof of the equivalence, valid regardless of whether the interaction is coherent or not.

Let us consider the interaction of two neutrinos, and, in the interests of clarity, for a moment suppress the Dirac indices on the fermions and the gamma matrices. Since the NC interactions conserve flavor, in the Hamiltonian (2) the flavor space wavefunction of a given outgoing neutrino is equal to that of one of the incoming neutrinos. Let us denote the flavor wavefunctions by $`\mathrm{\Psi }_i`$ and $`\mathrm{\Phi }_i`$. There are then two possible combinations:

$`\mathrm{\Psi }_i^{}\delta _{ij}\mathrm{\Psi }_j\mathrm{\Phi }_k^{}\delta _{kl}\mathrm{\Phi }_l,`$ (3)
$`\mathrm{\Psi }_i^{}\delta _{il}\mathrm{\Phi }_l\mathrm{\Phi }_k^{}\delta _{kj}\mathrm{\Psi }_j.`$ (4)
Since $`_{i=1}^2|\mathrm{\Psi }_i|^2=_{i=1}^2|\mathrm{\Phi }_i|^2=1`$, the term in Eq. (3) equals $`1`$. The term in Eq. (4) could be transformed using the well-known property of the (4-dimensional) $`\sigma `$-matrices, $`\sigma _{\alpha \dot{\alpha }}^m\overline{\sigma }_m^{\dot{\beta }\beta }=2\delta _\alpha ^\beta \delta _{\dot{\alpha }}^{\dot{\beta }}`$ (in the notation convention of ), or

$`2\delta _{il}\delta _{jk}=\delta _{ij}\delta _{kl}+\stackrel{}{\sigma }_{ij}\stackrel{}{\sigma }_{kl}.`$ (5)
One gets

$`\mathrm{\Psi }_i^{}\delta _{il}\mathrm{\Phi }_l\mathrm{\Phi }_j^{}\delta _{jk}\mathrm{\Psi }_k={\displaystyle \frac{1}{2}}(1+\mathrm{\Psi }_i^{}\stackrel{}{\sigma }_{ij}\mathrm{\Psi }_j\mathrm{\Phi }_k^{}\stackrel{}{\sigma }_{kl}\mathrm{\Phi }_l).`$ (6)
In this form, the equivalence between a system of neutrinos and a system of spins is manifest. The complete (including both the contributions from (6) and (3)) flavor space Hamiltonian for the interaction of two neutrinos, 1 and 2, is proportional to

$`(3/2+\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2/2),`$ (7)
which coincides with the square of the operator of the *total* angular momentum, $`\widehat{L}^2(\stackrel{}{\sigma }_1/2+\stackrel{}{\sigma }_2/2)^2`$, as expected.

The strength of the interaction depends on the scattering angle in real space. Let us denote the states of the incoming neutrinos as $`\chi _\alpha `$ and $`\psi _\beta `$ and those of the outgoing states as $`\overline{\psi }_{\dot{\beta }}^{}`$ and $`\overline{\chi }_{\dot{\alpha }}^{}`$. The spatial dependence of the scattering amplitude is given by

$$\overline{\psi }_{\dot{\beta }}^{}\overline{\sigma }^{\dot{\beta }\beta m}\psi _\beta \overline{\chi }_{\dot{\alpha }}^{}\overline{\sigma }_m^{\dot{\alpha }\alpha }\chi _\alpha =1(\overline{\psi }^{}\stackrel{}{\sigma }\psi )(\overline{\chi }^{}\stackrel{}{\sigma }\chi )$$
(8)
If the neutrino wavepackets are sufficiently broad so that several wavepackets overlap, the preceding argument can be trivially generalized to show that the interactions between the “overlapping” neutrinos have the structure of the interactions between the corresponding spins. In the case of neutrino plane waves in a box, each spin interacts with all the others and the total Hamiltonian is the sum of the interactions of all pairs (see Sect. 4.3).

### 2.2 Coherence of neutrino-neutrino scattering: basic features

We next discuss general properties of coherent neutrino-neutrino scattering.

By definition, scattering is coherent when the waves scattered by different particles in a target interfere with each other. When the scatterer is distinguishable from the incident particle, the coherence condition is satisfied when the incident particle is scattered forward. For other directions, scattering is incoherent (unless the particle spacings in the target satisfy particular conditions, e.g. periodicity). In contrast, for neutrino-neutrino scattering, scattering can be coherent when the momentum of the scattered neutrino coincides with the momentum of *one* of the incident neutrinos. The corresponding values of the scattering angle $`\beta `$ are $`0`$ and $`\mathrm{\Theta }`$, where $`\mathrm{\Theta }`$ denotes the angle between the two initial-state neutrinos.

For the purpose of studying coherent effects, one may replace the full Hamiltonian (2) with a “toy” Hamiltonian which restricts scattering angles to values $`0`$ and $`\mathrm{\Theta }`$ only, for which coherent interference may occur. This simplified interaction can be viewed as a flavor exchange between neutrinos which do not change their momenta. We emphasize that, while all coherent effects in the system are fully captured in this way, most of the incoherent effects are left out. Hence, one must be careful while interpreting incoherent effects found in this framework, as will be seen later.

Following , we write the reduced Hamiltonian for a pair of interacting neutrinos with momenta $`\stackrel{}{k}`$ and $`\stackrel{}{p}`$ in the form

$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}|\nu _e(\stackrel{}{k})\nu _e(\stackrel{}{p})\\ |\nu _e(\stackrel{}{k})\nu _\mu (\stackrel{}{p})\\ |\nu _\mu (\stackrel{}{k})\nu _e(\stackrel{}{p})\\ |\nu _\mu (\stackrel{}{k})\nu _\mu (\stackrel{}{p})\end{array}\right)={\displaystyle \frac{\sqrt{2}G_F}{V}}(1\mathrm{cos}\mathrm{\Theta })\left(\begin{array}{cccc}2& 0& 0& 0\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 0& 0& 0& 2\end{array}\right)\left(\begin{array}{c}|\nu _e(\stackrel{}{k})\nu _e(\stackrel{}{p})\\ |\nu _e(\stackrel{}{k})\nu _\mu (\stackrel{}{p})\\ |\nu _\mu (\stackrel{}{k})\nu _e(\stackrel{}{p})\\ |\nu _\mu (\stackrel{}{k})\nu _\mu (\stackrel{}{p})\end{array}\right),`$ (21)
where $`V`$ is the normalization volume. The value of the angular factor, $`(1\mathrm{cos}\mathrm{\Theta })`$, follows from (8) (taking into account the anticommutativity of the spinors), while the flavor structure of the interaction is determined by Eq. (7).

The condition $`\beta =\mathrm{\Theta }`$ by itself, however, *does not guarantee* coherence. We can see that by considering the scattering of an electron neutrino on two “background” muon neutrinos. Suppose that the initial momenta of the muon neutrinos are orthogonal to that of the electron neutrino and, further, that the neutrino wave packets are sufficiently narrow so that the two scattering events are independent, with each given by Eq. (21). Then, as a result of the interaction, we get

$$|\nu _e|\nu _\mu \nu _\mu |F=(1+2ia)|\nu _e|\nu _\mu \nu _\mu +ia(|\nu _\mu |\nu _e\nu _\mu +|\nu _\mu |\nu _\mu \nu _e),$$
(22)
where $`a\sqrt{2}G_F(1\mathrm{cos}\mathrm{\Theta })t/V`$.

Both the interaction time $`t`$ and the volume $`V`$ depend on the size of the wavepackets. Taking the wavepackets for simplicity to be spheres of size $`l`$, we can write

$`t`$ $``$ $`l,Vl^3,`$ (23)
$`a`$ $``$ $`G_F(1\mathrm{cos}\mathrm{\Theta })/l^2.`$ (24)
We notice that $`|a|1`$, as long as the size of the wavepackets is much greater than (100 GeV)-1, i.e., the neutrino energy is well below the weak scale. This justifies neglecting terms of order $`a^2`$ and higher in Eq. (22).

Let us compute the probability $`P`$ that, as a result of the interaction, the incident neutrino was converted into $`\nu _\mu `$. Introduce the “$`\nu _\mu `$ number” operator $`\widehat{\mu }|\mu \mu |`$ that acts only on the state of the incident neutrino. Using the orthogonality of the states $`|\nu _e\nu _\mu `$ and $`|\nu _\mu \nu _e`$, the probability in question is

$$P=F|\widehat{\mu }|F=2a^2.$$
(25)
This shows that the process is *incoherent*: the probability we found coincides with the sum of the conversion probabilities for each scattering (a coherent process would yield $`P=(2a)^2`$). As can be easily confirmed, as the number of “background” neutrinos $`N`$ is increased, the probability of flavor conversion goes as $`Na^2`$.

It can be further shown that if the background neutrinos are in a flavor superposition state $`\nu _x=\mathrm{cos}\alpha \nu _e+\mathrm{sin}\alpha \nu _\mu `$, the probability of conversion is *not*, in general, equal to the sum of the probabilities for each scattering. Indeed, a straightforward calculation yields $`P=a^2\mathrm{sin}^2\alpha (2+2\mathrm{cos}^2\alpha )`$. Comparing this with the conversion probability for a single scattering event, $`a^2\mathrm{sin}^2\alpha `$, we see that the conversion process is *partially coherent*, and the degree of coherence increases with decreasing $`\alpha `$. A generalization to a larger number of scatterers is

$$P=a^2\mathrm{sin}^2\alpha ((N^2N)\mathrm{cos}^2\alpha +N),$$
(26)
which clearly shows the presence of partial coherence (the terms proportional to $`N^2`$).

The above example shows that whether a neutrino undergoes coherent conversion or not depends on the flavor state of the many background neutrinos it interacts with. Moreover, as can be seen from Eq. (22), each neutrino-neutrino interaction creates an entangled state. Hence, the flavor evolution of a system of several interacting neutrinos demands a many-particle description. The question is: does the wavefunction of the system somehow factorize into individual one-particle wavefunctions in the limit of large $`N`$, and, if yes, under what conditions? Below we describe two recently proposed approaches to this problem.

## 3 Single- vs. many-particle description: two approaches

### 3.1 First approach: interference of many elementary scattering amplitudes

The first method was developed in . The idea is to describe the flavor evolution in a neutrino gas by generalizing the procedure of adding elementary scattering amplitudes, as was outlined in the three-neutrino example of Sect. 2.2.

Let us consider the interactions between two orthogonal neutrino beams. Once the results for this setup are understood, they can be straightforwardly generalized to more complicated systems. Let us assume that the neutrino wavepacket size is much smaller than the particle spacing ($`ln_\nu ^{1/3}`$), so that different neutrino-neutrino interactions can be treated independently using Eq. (21). As an initial configuration, we take the first beam to be made of $`N_1`$ electron neutrinos $`\nu _e`$, and the second beam to be made of $`N_2`$ neutrinos in the flavor superposition state $`\nu _x\mathrm{cos}\alpha \nu _e+\mathrm{sin}\alpha \nu _\mu `$. It can be shown  that the result of the interaction during a small time $`\delta t`$ can be written as

$`|eee\mathrm{}|xxx\mathrm{}`$ $`\stackrel{tt+\delta t}{}`$ $`|F=|eee\mathrm{}|xxx\mathrm{}+ia|F_1,`$ (27)
$`|F_1`$ $`=`$ $`N_1N_2(1+|e|x|^2)|eee\mathrm{}|xxx\mathrm{}`$ (28)
$`+`$ $`N_2\mu |xx|e(|\mu ee\mathrm{}+|e\mu e\mathrm{}+\mathrm{})|xxx\mathrm{}`$
$`+`$ $`N_1e|xy|e|eee\mathrm{}(|yxx\mathrm{}+|xyx\mathrm{}+\mathrm{})`$
$`+`$ $`\mu |xy|e(|\mu ee\mathrm{}+|e\mu e\mathrm{}+\mathrm{})(|yxx\mathrm{}+|xyx\mathrm{}+\mathrm{}),`$
a direct generalization of Eq. (22). Here the state $`|y`$ is defined to be orthogonal to $`|x`$: $`x|y=0`$.

The evolution in Eqs. (27,28) cannot be described using only one-particle equations, since the final state in Eq. (28) is an entangled many-particle state. Nevertheless, as observed in , *the coherent part* of the evolution can be. Indeed, the last term in Eq. (28) represents an incoherent effect, since it contains the sums of mutually orthogonal terms, and could be dropped. The sum of the remaining three terms is equal, to the first order in $`a`$, to a product of rotated single-particle states.

Since for small $`\delta t`$ a state of each neutrino undergoes a small rotation and since nothing in this argument depends on the particular choice of the initial states, Eqs. (27,28) specify how the states will rotate at any point in time. This observation means that, for each neutrino, one can write a differential equation describing its evolution for any — not necessarily small — $`t`$. This equation takes on the form

$`i{\displaystyle \frac{d|\psi ^{(i)}}{dt}}`$ $`=`$ $`H_{aa}|\psi ^{(i)},`$
$`H_{aa}`$ $`=`$ $`{\displaystyle \underset{j}{}}\sqrt{2}G_Fn^{(j)}(1\mathrm{cos}\mathrm{\Theta }^{(ij)})\left[|\varphi ^{(j)}\varphi ^{(j)}|+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}|\varphi ^{(j)}|\psi ^{(i)}|^2\right],`$ (29)
where, for generality, we restored the angular factor and summed over possible angles $`\mathrm{\Theta }^{(ij)}`$. The flavor-space wavefunction of the given neutrino is denoted by $`\psi ^{(i)}`$, those of the other (“background”) neutrinos by $`\varphi ^{(j)}`$. The quantity $`n^{(j)}`$ denotes the number density of neutrinos whose momenta make an angle $`\mathrm{\Theta }^{(ij)}`$ with the given neutrino. This equation is similar to the one-particle evolution equations given in . In fact, the only difference is in the last two terms in the square brackets, which have no effect on flavor evolution.

Let us discuss a particular case when the neutrinos at $`t=0`$ are in the flavor eigenstates. We note that, according to Eq. (3.1), a coherent conversion effect in this case is absent. This result in turn could be traced to Eq. (28), where the second and the third terms vanish for $`\nu _x=\nu _\mu `$ or $`\nu _x=\nu _e`$. The only remaining source of flavor transformation is the last term in Eq. (28), which is nonzero for $`\nu _x=\nu _\mu `$. Since this term represents incoherent scattering, for a neutrino traveling through a gas of neutrinos of opposite flavor the conversion probability for small $`t`$ has the form $`P(t)=N(t)a^2`$.

At first, it is not obvious how this incoherent effect is related to the incoherent flavor exchange that occurs in a real neutrino ensemble. Indeed, as mentioned in the introduction to Eq. (21), by replacing the true interaction Eq. (2) with the forward scattering Hamiltonian, we have truncated most of the incoherent effects. Nevertheless, a connection between the two can be made. The argument, which follows below, is very instructive and will prove helpful in the subsequent Sections.

When discussing the “forward scattering probability”, we in fact mean the probability for the neutrino to be scattered in a certain (small) solid angle $`\delta \mathrm{\Omega }`$ around the direction of the incoming neutrino (the probability of neutrino scattering at any fixed angle is, of course, zero). A neutrino wave packet with the cross section area $`A`$ and wavelength $`\lambda =2\pi /E_\nu `$ is diverging due to diffraction into a solid angle $`\delta \mathrm{\Omega }_{\mathrm{diff}}=\lambda ^2/A`$. The neutrino wave which is scattered into the same solid angle $`\delta \mathrm{\Omega }`$ will interfere with the incident wave. Thus, $`\delta \mathrm{\Omega }=\delta \mathrm{\Omega }_{\mathrm{diff}}`$.

With this in mind, let us consider the result of forward scattering of a given neutrino on $`N`$ other neutrinos in the medium. According to Eq. (21), for neutrino wave packets of longitudinal size $`l`$ and cross section $`A`$, $`aG_Fl(1\mathrm{cos}\mathrm{\Theta })/Al=G_F(1\mathrm{cos}\mathrm{\Theta })/A`$. The resulting flavor conversion probability $`P`$ depends on whether the forward scattering amplitudes add up coherently or incoherently. In the first case, $`P_{\mathrm{coh}}N^2a^2`$, in the second, $`P_{\mathrm{inc}}Na^2`$. The number of interactions during time $`t`$ is $`N=Atn_\nu `$ ($`n_\nu `$ is the neutrino number density) and we find that in the case of coherent forward scattering the dependence on the size of the wave packet *drops out*,

$$P_{\mathrm{coh}}\frac{G_F^2(1\mathrm{cos}\mathrm{\Theta })^2}{A^2}A^2t^2n_\nu ^2=G_F^2(1\mathrm{cos}\mathrm{\Theta })^2t^2n_\nu ^2.$$
(30)
This is indeed seen from Eq. (3.1), which does not contain the size of the neutrino wavepacket.

On the other hand, for the case of incoherent forward scattering, the cross section of the wave packet does enter the final result,

$$P_{\mathrm{inc}}\frac{G_F^2(1\mathrm{cos}\mathrm{\Theta })^2}{A^2}Atn_\nu =\frac{G_F^2(1\mathrm{cos}\mathrm{\Theta })^2}{A}tn_\nu .$$
(31)
The dependence on the size of the wavepacket look puzzling. Nevertheless, it is easily understood once we recall that, by construction, the result in (31) represents the fraction of incoherent scattering events for which one of the two neutrinos is scattered in the forward cone $`\delta \mathrm{\Omega }\lambda ^2/A1/E_\nu ^2A`$. This means that, up to a numerical coefficient, the probability of incoherent scattering in any direction $`P_{\mathrm{inc}}^{\mathrm{TOT}}`$ is given by

$$P_{\mathrm{inc}}^{\mathrm{TOT}}P_{\mathrm{inc}}/\delta \mathrm{\Omega }G_F^2E_\nu ^2(1\mathrm{cos}\mathrm{\Theta })^2tn_\nu .$$
(32)
This has the form $`\sigma _wtn_\nu `$, where $`\sigma _w`$ is a typical weak interaction cross section, $`\sigma _wG_F^2E_\nu ^2`$, exactly as one would expect.

Thus, in summary, for the purpose of studying coherent effects, the full Hamiltonian (2) may be replaced by the forward scattering flavor exchange Hamiltonian, as done in Eq. (21), with the understanding that in the final result the terms that go like $`N^2`$ represent the *complete* coherent effect while the terms that go like $`N`$ correspond to only the *fraction* of the incoherent scattering that occurs in the forward diffraction cone.

### 3.2 Second approach: solving a many-particle equation for neutrino plane waves in a box

An alternative method of studying coherent processes in a neutrino gas was proposed in . The idea is to consider a system of interacting neutrino plane waves in a box and study its flavor evolution, again with the restriction to forward scattering only. The neutrinos are taken to be initially in flavor states and the goal is to test the existence of coherent conversion effects by looking at the rate of equilibration of the system.

The “forward scattering” Hamiltonian can be found by generalizing Eq. (21). A basis of states is formed by the initial configuration of neutrinos and all its possible permutations. Taking, e.g., $`N`$ electron and $`N`$ muon neutrinos, there are $`(2N)!/(N!)^2`$ distinct basis states. The entries of the Hamiltonian are found as follows. (i) The diagonal entries receive $`N(3N2)`$ contributions (this includes $`(2N)!/[2!(2N2)!]=N(2N1)`$ flavor-blind scattering processes with scattering angle $`\beta =0`$ and $`2N!/[2!(N2)!]=N(N1)`$ processes with scattering angle $`\beta =\mathrm{\Theta }`$). (ii) Each off-diagonal entry receives a single contribution from the exchange process that connects the two corresponding basis states; if the two basis states are not connected by a single permutation, the corresponding off-diagonal entry vanishes. (iii) Each contribution equals $`\sqrt{2}G_F(1\mathrm{cos}\mathrm{\Theta }_{ij})/V`$, where $`\mathrm{\Theta }_{ij}`$ is the angle between the momenta of the two interacting neutrinos.

For example, for the case of $`N=2`$ there are $`6`$ basis states: $`|\nu _e\nu _e\nu _\mu \nu _\mu `$, $`|\nu _\mu \nu _e\nu _e\nu _\mu `$, $`|\nu _\mu \nu _e\nu _\mu \nu _e`$, $`|\nu _e\nu _\mu \nu _e\nu _\mu `$, $`|\nu _e\nu _\mu \nu _\mu \nu _e`$, $`|\nu _\mu \nu _\mu \nu _e\nu _e`$, and the Hamiltonian is given by

$`H_{N=2}={\displaystyle \frac{\sqrt{2}G_F}{V}}\left(\begin{array}{cccccc}d_{12,34}\hfill & f_{13}\hfill & f_{14}\hfill & f_{23}\hfill & f_{24}\hfill & 0\hfill \\ f_{13}\hfill & d_{14,23}\hfill & f_{34}\hfill & f_{12}\hfill & 0\hfill & f_{24}\hfill \\ f_{14}\hfill & f_{34}\hfill & d_{13,24}\hfill & 0\hfill & f_{12}\hfill & f_{23}\hfill \\ f_{23}\hfill & f_{12}\hfill & 0\hfill & d_{13,24}\hfill & f_{34}\hfill & f_{14}\hfill \\ f_{24}\hfill & 0\hfill & f_{12}\hfill & f_{34}\hfill & d_{14,23}\hfill & f_{13}\hfill \\ 0\hfill & f_{24}\hfill & f_{23}\hfill & f_{14}\hfill & f_{13}\hfill & d_{12,34}\hfill \end{array}\right).`$ (39)
Here $`f_{ij}=1\mathrm{cos}\mathrm{\Theta }_{ij}`$ is the angular factor for the $`i,j`$ pair and $`d_{ij,kl}=f_{ij}+f_{kl}+`$, where $``$ is the sum of all distinct $`f_{ij}`$.

Let us discuss the meaning of the equilibration time $`t_{\mathrm{eq}}`$. If conversion is coherent, $`t_{\mathrm{eq}}`$ is defined as the time over which a neutrino of a given momentum is expected to change its flavor. The momentum itself stays unchanged over this time scale. If conversion is incoherent, $`t_{\mathrm{eq}}`$ is defined as the time over which the neutrino momenta will be randomized (each neutrino can no longer be identified by its momentum). In the model we are considering, however, the definition of equilibration for coherent process also applies to the incoherent process, since the interaction is chosen to preserve the neutrino momentum.

Let us determine what dependence of $`t_{\mathrm{eq}}`$ would be a sign of coherent or incoherent conversion. To do this, it is useful to repeat the exercise of Sect. 3.1 of estimating the conversion probability as a function of time. For neutrino plane waves in the box the following two modifications to the argument need to be made: first, the normalization volume $`V`$ becomes the size of the box and, second, the interaction time becomes the actual time $`t`$ elapsed from the beginning of the evolution. The scattering amplitude becomes $`aG_Ft(1\mathrm{cos}\mathrm{\Theta })/V`$ and, since each neutrino simultaneously interact with all the others, the number of interactions is $`N=Vn_\nu `$. Using these ingredients, we find

$`P_{\mathrm{coh}}`$ $``$ $`{\displaystyle \frac{G_F^2(1\mathrm{cos}\mathrm{\Theta })^2t^2}{V^2}}V^2n_\nu ^2=G_F^2(1\mathrm{cos}\mathrm{\Theta })^2t^2n_\nu ^2,`$ (40)
$`P_{\mathrm{inc}}`$ $``$ $`{\displaystyle \frac{G_F^2(1\mathrm{cos}\mathrm{\Theta })^2t^2}{V^2}}Vn_\nu ={\displaystyle \frac{G_F^2(1\mathrm{cos}\mathrm{\Theta })^2}{V}}t^2n_\nu .`$ (41)
Here, as in the case of the interacting neutrino wavepackets, we find that in $`P_{\mathrm{coh}}`$ the dependence on the size of the box cancels out (and the result is the same as Eq. (30)), while for incoherent scattering the size of the box enters the final result. Since, as was already stressed, the incoherent scattering result is, in some sense, artificial, there is no cause for alarm.

Inverting Eqs. (40,41) and Eqs. (30,31), we find that the equilibration time for coherent scattering in both cases has the dependence

$$t_{\mathrm{eq}}^{coh}=(G_Fn_\nu )^1,$$
(42)
while for incoherent scattering one finds

$$t_{\mathrm{eq}}^{inc}=(G_F^2n_\nu /A)^1$$
(43)
for the case of neutrino wavepackets, and

$$t_{\mathrm{eq}}^{inc}=(G_F\sqrt{n_\nu /V})^1$$
(44)
for plane waves in a box. In what follows, we solve the evolution of the system of neutrinos in the box completely analytically. This gives the equilibration time $`t_{eq}`$, which we compare to Eqs. (42) and (44) to determine the character of the evolution.

## 4 Plane waves in a box: analytical treatment

### 4.1 Setup

As done in , we consider a system of $`2N`$ neutrino plane waves, of which $`N`$ are initially in the $`\nu _e`$ state and the remaining $`N`$ are in the $`\nu _\mu `$ state. We write this configuration of the system as:

$$|S(0)=|\underset{N}{\underset{}{\nu _e\nu _e\mathrm{}\nu _e}}\underset{N}{\underset{}{\nu _\mu \nu _\mu \mathrm{}\nu _\mu }}.$$
(45)
Here each neutrino wave is identified by its position in the list, which corresponds to a given momentum.

The Hamiltonian for this system has the form described in Sect. 3.2. Since it is equivalent to the Hamiltonian of a system of spins (Sect. 2.1), in the following we adopt the terminology of angular momenta. As a convention, we identify the “up” (pointing along the $`z`$ axis) state of a spin with $`\nu _e`$ and the “down” state with $`\nu _\mu `$. Thus, the initial state (45) can be expressed as

$$|S(0)=|\underset{N}{\underset{}{+\frac{1}{2},+\frac{1}{2},\mathrm{}+\frac{1}{2}}},\underset{N}{\underset{}{\frac{1}{2},\frac{1}{2},\mathrm{}\frac{1}{2}}}.$$
(46)
To make an analytical treatment possible, we set

$$f_{ij}=1$$
(47)
for all neutrino pairs. The resulting Hamiltonian describes an ensemble of spins interacting with each other with equal strength222This makes it different from the familiar Ising model.. This simplifying assumption does not change the coherent (or incoherent) character of the conversion effects, as also pointed out in . With this assumption, the nonzero off-diagonal entries of the Hamiltonian become equal to $`\sqrt{2}G_F/V`$ while the diagonal entries become $`\sqrt{2}G_F/V\times N(3N2)`$.

We recall (Sect. 3.1) that for the initial configuration (45) the single-particle formalism predicts no coherent effects, as can be seen from Eq. (3.1). We investigate if this conclusion is changed once the simultaneous interactions between all neutrinos are considered.

### 4.2 Results

Let us first display our results. Their derivation and a more detailed discussion is given in Sec. 4.3.

The probability $`P_1`$ that the first spin – which is in the state $`|+1/2`$ at $`t=0`$ – is found in the same up state at time $`t`$ ($`t>0`$) can be calculated by a completely analytical procedure. It is given by the Fourier series:

$`P_1(t)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{J=0}{\overset{N1}{}}}\eta (N,J)\mathrm{cos}[(2J+2)gt],`$ (48)
where $`g\sqrt{2}G_F/V`$. The coefficients $`\eta (N,J)`$ are the squares of the Clebsch-Gordan coefficients $`j_1m_1j_2m_2|JMC(j_1,m_1,j_2,m_2,J,M)`$ for the addition of two angular momenta, $`j_1=m_1=(N1)/2`$ and $`j_2=m_2=N/2`$:

$`\eta (N,J)`$ $`=`$ $`\left[C({\displaystyle \frac{N1}{2}},{\displaystyle \frac{N1}{2}},{\displaystyle \frac{N}{2}},{\displaystyle \frac{N}{2}},J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}})\right]^2`$ (49)
$`=`$ $`{\displaystyle \frac{(1+2J)N[(N1)!]^2}{(J+N+1)!(NJ1)!}}.`$
We follow the standard notation, in which $`j_1,j_2`$ denote the values of the angular momenta being added, $`m_1,m_2`$ denote their projections along the $`\widehat{z}`$ direction, and $`J`$ and $`M`$ are the corresponding values for the total spin.

Figure 1 shows the survival probability $`P_1(t)`$, obtained by summing numerically the series (48). The horizontal axis represents time in units of the equilibration time predicted for coherent conversion, $`t_{eq}^{coh}=(2Ng)^1`$ (Eq. (42) for $`2N`$ neutrinos). The curves refer to $`N=40,160,640,2560`$. Having the same choice of units, our plot can be directly compared to fig. 2 of ref. , where, however, smaller values of $`N`$ were used333In ref.  a different form of the spin-spin (neutrino-neutrino) coupling was adopted, namely, the angular factors were randomly generated. This produces a numerical difference with respect to our results, but does not affect the conclusions about the equilibration time.. From Fig. 1, it is clear that the equilibration time $`t_{eq}`$ is larger than $`t_{eq}^{coh}`$ and the discrepancy grows with $`N`$.

In the limit $`N1`$ the series (48) is well approximated by the following function:

$`P_1(t)=1{\displaystyle \frac{\sqrt{\pi }}{2}}\sqrt{N}gt\mathrm{exp}(Ng^2t^2)\mathrm{erfi}(\sqrt{N}gt),`$ (50)
where $`\mathrm{erfi}(z)`$ is the imaginary error function

$$\mathrm{erfi}(z)i\mathrm{erf}(iz)=\frac{2}{\sqrt{\pi }}\times _0^z\mathrm{exp}(t^2)𝑑t.$$
(51)
We have checked that for the values of $`N`$ shown in Fig. 1 the approximation is very good and the curves obtained using (50) are indistinguishable from those obtained from Eq. (48).

From Eq. (50) we see that the characteristic equilibration time equals:

$$t_{eq}=(g\sqrt{N})^1.$$
(52)
This scaling of the equilibration time can be clearly seen in Fig. 1. It is precisely (up to numerical factors) the form (44) expected if the equilibration process is incoherent.

### 4.3 Derivation of results

#### 4.3.1 The eigenvalues of the Hamiltonian

From Eqs. (7) and (47) it follows that for each pair $`i`$ and $`j`$ the interaction is given by $`H_{ij}=g\widehat{L}_{ij}^2=g(\stackrel{}{\sigma }_i+\stackrel{}{\sigma }_j)^2`$ and, therefore, the Hamiltonian for the whole system is

$$H_N=g\underset{i=1}{\overset{2N1}{}}\underset{j=i+1}{\overset{2N}{}}(\stackrel{}{\sigma }_i+\stackrel{}{\sigma }_j)^2=g\left(\underset{i=1}{\overset{2N1}{}}\underset{j=i+1}{\overset{2N}{}}2\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j+\frac{3}{2}N(2N1)\right),$$
(53)
where the summation is done over all $`N(2N1)`$ pairs.

We notice that this Hamiltonian can be related to the operator of the square of the total angular momentum of the system,

$$\widehat{L}^2\left(\underset{i=1}{\overset{2N}{}}\stackrel{}{\sigma }_i\right)^2=\underset{i=1}{\overset{2N1}{}}\underset{j=i+1}{\overset{2N}{}}2\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j+\frac{3}{2}N.$$
(54)
By comparing Eqs. (53) and (54) we find

$$H_N=g[\widehat{L}^2+3N(N1)],$$
(55)
which has the eigenvalues

$$E_{N,J}=g[J(J+1)+3N(N1)],\mathrm{with}J=0,1,\mathrm{},N.$$
(56)
Since the flavor evolution will depend only on the difference of the eigenvalues (as will be shown later), in the following we will omit the second term in Eq. (56) and write

$$E_J=gJ(J+1).$$
(57)
A given value of $`E_J`$ in general corresponds to more than one state. This degeneracy can be seen by simply counting the states. For our specific setup, with $`M=0`$, there are $`(2N)!/(N!)^2`$ basis states, while there are only $`N+1`$ different values of $`J`$. In general, therefore, an orthonormal basis of eigenstates of the Hamiltonian has the form:

$$|J,k_J,$$
(58)
where the index $`k_J`$ labels eigenstates with the same value of $`J`$.

As an illustration of the present discussion, we give the Hamiltonian for the case $`N=2`$ of Sect. 3.2:

$`H_{N=2}=g\left(\begin{array}{cccccc}8\hfill & 1\hfill & 1\hfill & 1\hfill & 1\hfill & 0\hfill \\ 1\hfill & 8\hfill & 1\hfill & 1\hfill & 0\hfill & 1\hfill \\ 1\hfill & 1\hfill & 8\hfill & 0\hfill & 1\hfill & 1\hfill \\ 1\hfill & 1\hfill & 0\hfill & 8\hfill & 1\hfill & 1\hfill \\ 1\hfill & 0\hfill & 1\hfill & 1\hfill & 8\hfill & 1\hfill \\ 0\hfill & 1\hfill & 1\hfill & 1\hfill & 1\hfill & 8\hfill \end{array}\right).`$ (65)
The eigenvalues, which are degenerate as expected, are $`g(6,6,8,8,8,12)`$. They correspond to the following values of $`J`$: $`(0,0,1,1,1,2)`$, as can be seen from Eq. (56).

#### 4.3.2 The Fourier series

Let us start with proving Eq. (48).

In the interest of clarity, we first display some definitions. Let us denote as $`|+1/2_1`$ and $`|1/2_1`$ the “up” and “down” states of the first spin. Given a state $`|\psi `$ of the $`2N`$ spins of the system, it will be convenient to use the projections $`+1/2|_1|\psi `$ and $`1/2|_1|\psi `$, which are states of $`2N1`$ spins. In particular, we introduce the state $`|S^{}`$ obtained by removing the first spin from the initial configuration $`|S(0)`$ (Eq. (46)):

$`+1/2|_1|S(0)|S^{}|\underset{N1}{\underset{}{+{\displaystyle \frac{1}{2}},+{\displaystyle \frac{1}{2}},\mathrm{}+{\displaystyle \frac{1}{2}}}},\underset{N}{\underset{}{{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},\mathrm{}{\displaystyle \frac{1}{2}}}}.`$ (66)
For any state, we can define another state obtained from the first one upon a reflection along the $`z`$-axis. In this transformation, each constituent spin $`|+1/2`$ is transformed into spin $`|1/2`$ and vice versa. For instance, the application of this transformation to $`|S^{}`$ gives its “mirror” state:

$`|S^+|\underset{N1}{\underset{}{{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},\mathrm{}{\displaystyle \frac{1}{2}}}}\underset{N}{\underset{}{+{\displaystyle \frac{1}{2}},+{\displaystyle \frac{1}{2}},\mathrm{}+{\displaystyle \frac{1}{2}}}}.`$ (67)
In what follows we will exploit the properties of various states under this reflection.

We adopt the notation $`|j,m,s_n`$ to indicate the result of the following procedure: consider the state $`|s`$ and take its projection on the subspace formed by the states of $`n`$ spins that have the total angular momentum $`j`$ and projection $`m`$ along the $`\widehat{z}`$ axis; normalize the result to 1. As an example, $`|J,0,S(0)_{2N}`$ indicates the normalized projection of $`S|(0)`$ on the subspace of the states of $`2N`$ spins having total angular momentum $`J`$ and $`\widehat{z}`$-projection equal to zero.

It is also convenient to define the state:

$$|s_J^{}+1/2|_1|J,0,S(0)_{2N},$$
(68)
(which has $`2N1`$ spins) and denote as $`|s_J^+`$ its mirror state obtained upon $`\widehat{z}`$-reflection.

The proof is constituted by four parts.

1. Expansion of the initial state

The initial state can be expanded on the basis of the eigenstates of the Hamiltonian:

$$|S(0)=\underset{J,k_J}{}|J,k_JJ,k_J|S(0),$$
(69)
which then can be evolved in time as

$$|S(t)=\underset{J,k_J}{}e^{iE_Jt}|J,k_JJ,k_J|S(0).$$
(70)
Since the energy depends only on $`J`$, the summation over $`k_J`$ can be done, yielding

$`|S(0)={\displaystyle \underset{J=0}{\overset{N}{}}}A_J|J,0,S(0)_{2N},`$
$`|S(t)={\displaystyle \underset{J=0}{\overset{N}{}}}A_Je^{iE_Jt}|J,0,S(0)_{2N}.`$ (71)
The factors $`A_J`$ equal the Clebsch-Gordan coefficients:

$$A_J=C(\frac{N}{2},\frac{N}{2},\frac{N}{2},\frac{N}{2},J,0).$$
(72)
This can be understood considering that the initial state, $`|S(0)`$, equals the sum of two angular momenta with $`j_1=N/2`$, $`j_2=N/2`$ and projections $`m_1=N/2`$, $`m_2=N/2`$. The first momentum is given by the sum of all the “up” spins, while the second is the sum of all the “down” spins.

2. Decomposition into states of $`1`$ and $`2N1`$ spins

It is useful to decompose the state $`|S^{}`$ as:

$$|S^{}=\underset{J=0}{\overset{N1}{}}K_{J+1/2}|J+\frac{1}{2},\frac{1}{2},S^{}_{2N1}.$$
(73)
Similarly to Eq. (72), here we have:

$$K_{J+1/2}=C(\frac{N1}{2},\frac{N1}{2},\frac{N}{2},\frac{N}{2},J+\frac{1}{2},\frac{1}{2})=\sqrt{\eta (N,J)},$$
(74)
where $`\eta (N,J)`$ was given in Eq. (49). The argument to explain (74) is analogous to that used to justify Eq. (72): since the state $`|S^{}`$ has $`N1`$ spins “up” and N spins “down”, it is given by the sum of two angular momenta with $`j_1=(N1)/2`$, $`j_2=N/2`$ and projections $`m_1=(N1)/2`$, $`m_2=N/2`$. These momenta are the result of summing all the spins “up” and all the spins “down” separately.

Each of the states $`|J,0,S(0)_{2N}`$ can be decomposed as:

$`|J,0,S(0)_{2N}`$ $`=`$ $`a_J|+{\displaystyle \frac{1}{2}}_1|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},s_J^{}_{2N1}+b_J|+{\displaystyle \frac{1}{2}}_1|J{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},s_J^{}_{2N1}`$ (75)
$`+`$ $`c_J|{\displaystyle \frac{1}{2}}_1|J+{\displaystyle \frac{1}{2}},+{\displaystyle \frac{1}{2}},s_J^+_{2N1}+d_J|{\displaystyle \frac{1}{2}}_1|J{\displaystyle \frac{1}{2}},+{\displaystyle \frac{1}{2}},s_J^+_{2N1},`$
where the states $`|J+\frac{1}{2},+\frac{1}{2},s_J^+_{2N1}`$ and $`|J\frac{1}{2},+\frac{1}{2},s_J^+_{2N1}`$ in the second line are related to the corresponding ones in the first line by the $`z`$-reflection operation. Notice that the state $`|J,0,S(0)_{2N}`$ has parity $`(1)^{NJ}`$ under $`z`$-reflection. Requiring that the right hand side of Eq. (75) has the same parity yields

$`c_J=(1)^{NJ}a_Jd_J=(1)^{NJ}b_J.`$ (76)
We also note that

$`a_N=c_N=0b_0=d_0=0,`$ (77)
$`|a_J|^2+|b_J|^2=|c_J|^2+|d_J|^2=1/2,`$ (78)
where the last relation represents the normalization condition.

3. Properties of the decomposition (75).

We can calculate the state $`1/2|_1|S(t)`$ and demand that it vanishes at $`t=0`$. This is equivalent to requiring a null probability to find the first spin initially in the down state. Combining this condition with Eqs. (71), (75) and (76) we find

$`{\displaystyle \underset{J=0}{\overset{N1}{}}}A_Ja_J(1)^{NJ}|J+{\displaystyle \frac{1}{2}},+{\displaystyle \frac{1}{2}},s_J^+_{2N1}`$
$`+{\displaystyle \underset{J=0}{\overset{N1}{}}}A_{J+1}b_{J+1}(1)^{NJ1}|J+{\displaystyle \frac{1}{2}},+{\displaystyle \frac{1}{2}},s_{J+1}^+_{2N1}=0,`$ (79)
from which it follows that

$$A_Ja_J|J+\frac{1}{2},+\frac{1}{2},s_J^+_{2N1}=A_{J+1}b_{J+1}|J+\frac{1}{2},+\frac{1}{2},s_{J+1}^+_{2N1},$$
(80)
and, by $`z`$-reflection:

$$A_Ja_J|J+\frac{1}{2},\frac{1}{2},s_J^{}_{2N1}=A_{J+1}b_{J+1}|J+\frac{1}{2},\frac{1}{2},s_{J+1}^{}_{2N1}.$$
(81)
4. Calculation of the survival probability $`P_1`$.

Let us calculate the probability $`P_1(t)`$ that at the time $`t`$ the first spin is found in its initial state, $`|+1/2`$:

$$P_1(t)=S(t)|+1/2_1+1/2|_1|S(t).$$
(82)
The combination of Eqs. (71), (75) and (77), gives the expression of $`+1/2|_1|S(t)`$:

$`+1/2|_1|S(t)=`$ $`{\displaystyle \underset{J=0}{\overset{N1}{}}}e^{iE_Jt}A_Ja_J|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},s_J^{}_{2N1}`$ (83)
$`+{\displaystyle \underset{J=0}{\overset{N1}{}}}e^{iE_{J+1}t}A_{J+1}b_{J+1}|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},s_{J+1}^{}_{2N1},`$
which, taking into account the result (81), simplifies to:

$`+1/2|_1|S(t)={\displaystyle \underset{J=0}{\overset{N1}{}}}\left(e^{iE_Jt}+e^{iE_{J+1}t}\right)A_Ja_J|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},s_J^{}_{2N1}.`$ (84)
The final step to conclude the proof is the observation that for $`t=0`$, Eq. (84) should coincide with the expansion of the state $`|S^{}`$, given in Eq. (73). From comparison, using the equality (74), it follows that:

$`2A_Ja_J|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},s_J^{}_{2N1}`$ $`=`$ $`K_{J+1/2}|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},S^{}_{2N1}`$ (85)
$`=`$ $`\sqrt{\eta (N,J)}|J+{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},S^{}_{2N1}.`$
Combining Eqs. (82), (84), (85) and (57), we find

$$P_1(t)=\underset{J=0}{\overset{N1}{}}\frac{\eta (N,J)^2}{4}|e^{iE_Jt}+e^{iE_{J+1}t}|^2,$$
(86)
from which Eq. (48) straightforwardly follows.

#### 4.3.3 The physical (large N) limit: resummed function

To derive Eq. (50) we first elaborate $`\eta (N,J)`$, Eq. (49), in the limit $`N1`$. We approximate the factorials using the Stirling formula: $`z!\sqrt{2\pi }e^{(z1/2)\mathrm{ln}zz}`$. As a result, an exponential factor appears in the expression of $`\eta (N,J)`$. We expand this exponential in powers of $`J/N`$ and keep only the leading terms. This is justified because the various factors in Eq. (49) are significantly different from zero only for $`JN`$. From the procedure described, we obtain:

$$\eta (N,J)2N^1(1+J)\mathrm{exp}\left[(J+1)^2/N\right].$$
(87)
This function has a peak at $`J=\sqrt{N/2}1`$ and has a characteristic width $`2\sigma \sqrt{N}`$. This fact provides an a posteriori justification for taking $`JN`$ for large $`N`$.

As a second step, we notice that, for $`N1`$ the Fourier sum in Eq. (48) is well approximated by the integral:

$`{\displaystyle \underset{J=0}{\overset{N1}{}}}\eta (N,J)\mathrm{cos}[(2J+2)gt]`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑J{\displaystyle \frac{2}{N}}Je^{J^2/N}\mathrm{cos}(2Jgt)`$ (88)
$`=`$ $`1\sqrt{\pi N}gte^{Ng^2t^2}\mathrm{erfi}(\sqrt{\mathrm{N}}\mathrm{gt}),`$
where the expression (87) for $`\eta `$ has been used. The function $`\mathrm{erfi}(z)`$ is given in Eq. (51). From Eq. (88) and (48) the final result (50) follows immediately.

## 5 Discussion

1. Let us compare the main features of our solution, Eqs. (48) and (50), with the corresponding results for a gas of small neutrino wavepackets. In the latter case the survival probability exhibits a different dependence on time, characterized by a simple exponential form with the argument linear in $`t`$. To see this, consider a thought experiment in which a beam of $`2N`$ electron neutrinos enters the medium containing an equal mixture of electron and muon neutrinos. It is straightforward to see that, for incoherent forward scattering, the flavor composition of the beam will change according to

$`dN_e/dt`$ $`=`$ $`\lambda (N_eN_\mu ),`$ (89)
$`dN_\mu /dt`$ $`=`$ $`\lambda (N_\mu N_e),`$ (90)
where $`\lambda `$ denotes the frequency of the flavor changing collisions for a given neutrino. The solution is

$`N_e(t)`$ $`=`$ $`N(1+\mathrm{exp}(\lambda t)),`$ (91)
a simple exponential. The origin of the difference between Eq. (50) and (91) lies in the different physics of the interaction. In the case of the wavepackets, the time elapsed from $`t=0`$ tells us how many background neutrinos a given neutrino interacted with; the interaction time between neutrino pairs is very short. On the other hand, in the case of the interacting plane waves, any two waves in the box continuously interact with each other starting from $`t=0`$. The scattering amplitude (and the effective interaction cross section) grows with time. If the conversion process were coherent, the survival probability would not depend on the particular model adopted (see Sect. 3.2).

2. It is interesting to note that the Fourier series in Eq. (48) describes a periodic function with period $`T=\pi /g`$. This periodicity is lost when a transition to the integral is made. Indeed, in Eq. (88) the lowest frequency over which the integration is performed is zero ($`J=0`$), corresponding to infinite period. It is clear that the period $`T`$ is much larger than the equilibration time $`t_{eq}`$ (see Eq. (52)): $`T=\pi \sqrt{N}t_{eq}`$. Moreover, one can see that typically $`T`$ exceeds the age of the universe: for a volume $`V=1\mathrm{cm}^3`$ we find $`T10^{22}`$ s. Therefore the integral form (50) is an accurate approximation of the exact result over time scales (at least) of the order of $`t_{eq}`$, and applies to all the realistic physical situations.

3. Let us comment on the range of validity of our results and their possible extensions.

(a) Due to the $`SU(2)`$ invariance of the problem, the results are the same if we replace the flavor eigenstates $`\nu _e,\nu _\mu `$ with any pair of states in the flavor space, $`\nu _x,\nu _y`$, related to the flavor basis by an $`SU(2)`$ transformation. Explicitly, if we take neutrinos initially in the states $`\nu _x`$ and $`\nu _y`$, the probability of conversion $`\nu _x\nu _y`$ does not exhibit coherent effects in the large $`N`$ limit. This tells us that no coherent conversion is realized if coherence is initially absent in the system, i.e. if, for a given neutrino $`\nu _i`$ being in the flavor state $`\nu _x`$, the other neutrinos $`\nu _j`$ in the ensemble are *not* coherent superpositions of $`\nu _x`$ and $`\nu _y`$: $`\nu _i|\nu _j=0`$ or $`|\nu _i|\nu _j|=1`$.

(b) Several generalizations of our calculation are possible. For instance, one could investigate how the results depend on the initial configuration of the system. In  we have studied the case in which coherence is initially present in the system ($`|\nu _i|\nu _j|<1`$) and/or the initial state is entangled. The findings were in agreement with the one-particle description, Eq. (3.1).

(c) The assumption of constant spin-spin coupling, $`f_{ij}=1`$, could be relaxed to describe the case of, e.g., randomly distributed couplings or interactions which depend on the spin-spin distance or on the specific geometry of the system. These cases may be difficult or impossible to describe analytically; however we expect them to be characterized by the same incoherent character of the conversion.

(d) The generalization to include vacuum oscillations and interaction with ordinary matter can done by adding the appropriate (well known) terms to the Hamiltonian.

## 6 Conclusions

We have found that a system of many ($`N1`$) (forward-)interacting neutrino plane waves in a box reaches flavor equilibration after a time $`tt_{eq}=(\sqrt{2}G_F\sqrt{N}/V)^1`$ if the neutrinos are initially in flavor states. This result agrees with what is expected from incoherent scattering for this system, and therefore indicates the absence of coherent conversion effects.

Although our results were obtained under specific simplifying conditions (plane-waves in a box, constant couplings, forward scattering only, etc.), they can be used to draw important general conclusions on neutrino conversion in a neutrino background. The absence of coherent conversion has general character, since it does not depend on the size of the neutrino wavepackets (see Sect. 3.1). It contrasts with ref.  and agrees with the prediction of the traditional – and widely applied – single-particle description. This tells us that there is no reason to believe in a breakdown of this description and that, therefore, its applications have to be considered valid. The significance of this statement is clear in consideration of the important physical problems, in cosmology and in stellar physics, in which the neutrino-neutrino coherent scattering plays a role.

Finally, we remark that, in principle, the analytical solution we have found for a system of interacting spins may be interesting on its own and may find applications beyond the neutrino evolution problem considered here.

###### Acknowledgments.

We are grateful to G. Raffelt for valuable comments. A. F. was supported by the Department of Energy, under contract W-7405-ENG-36. A. F. is a Feynman Fellow at LANL. C. L. was supported by the Keck foundation and by the National Science Foundation grant PHY-0070928.