text:
# A comprehensive study of neutrino spin-flavour conversion in supernovae and the neutrino mass hierarchy

## 1 Introduction

A core-collapse supernova explosion is one of the most spectacular events in astrophysics, and it attracts a great deal of attention from many physicists and astronomers. It also produces a number of neutrinos and 99% of its gravitational binding energy is transformed to neutrinos. Therefore, neutrinos play an essential role in supernovae, and their detection by ground-based large water Čerenkov detectors, such as Super-Kamiokande (SK) and Sudbury Neutrino Observatory (SNO), would provide valuable information on the nature of neutrinos as well as supernova physics. What we can learn from the next galactic supernova has been considered in many articles (for a review, see ). For example, we can constrain the properties of neutrino oscillations, such as the mixing angle between the first and third mass eigenstates ($`\theta _{13}`$), and the mass hierarchy \[normal ($`m_1m_3`$) or inverted ($`m_1m_3`$)\] .

In addition to the nonzero neutrino masses and mixing angles, the nonzero magnetic moment is another nature of neutrinos beyond the standard model of particle physics. If neutrinos have a nonzero magnetic moment, it leads to precession between left- and right-handed neutrinos in sufficiently strong magnetic fields . In general, nondiagonal elements of the magnetic moment matrix are possible and neutrinos can be changed into different flavours and chiralities . Furthermore, with the additional effect of coherent forward scattering by matter, neutrinos can be resonantly converted into those with different chiralities  by a mechanism similar to the well-known Mikheyev–Smirnov–Wolfenstein (MSW) effect . This resonant spin-flavour (RSF) conversion induced by the neutrino magnetic moment in strong magnetic fields was first introduced to solve the solar neutrino problem, and actually gave the best fit solution before the KamLAND result . However, the recent KamLAND experiment  has shown that the large mixing angle (LMA) MSW solution is the most favourable one; the RSF mechanism is suppressed at the subdominant level. From the KamLAND negative results for the solar antineutrino search, an upper bound on the neutrino magnetic moment is obtained, $`\mu _\nu 1\times 10^{12}\mu _\mathrm{B}`$, where $`\mu _\mathrm{B}`$ is the Bohr magneton . This upper bound is comparable to the most stringent limit from the stellar cooling argument, $`\mu _\nu `$(1–4)$`\times 10^{12}\mu _\mathrm{B}`$ .

Although the RSF mechanism does not work at a dominant level in the Sun, it may occur efficiently in a denser environment with stronger magnetic field, which is actually expected in the case of core-collapse supernovae. The RSF conversion mechanism in supernovae has been investigated by many authors . Among them, Ando and Sato  have studied the RSF effect using a three-flavour formulation with the latest oscillation parameters, and pointed out that the combination of the MSW and RSF effects makes the crossing scheme very interesting to investigate. Since the RSF conversions are very sensitive to the value of $`Y_\mathrm{e}`$, which is the electron number fraction per nucleon, they have also investigated the dependence of the RSF effect on presupernova models with solar and zero metallicities . It is concluded that the efficient (either complete or incomplete, depending on presupernova models) RSF conversions occur when the supernova magnetic field is sufficiently strong, $`\mu _\nu B_0(10^{12}\mu _\mathrm{B})(10^{10}\mathrm{G})`$, where $`B_0`$ is the strength of the magnetic field at the surface of the iron core.

However, all the past studies of three-flavour RSF effect were based on the assumption that the neutrino mass hierarchy is normal ($`m_1m_3`$), although the inverted mass hierarchy ($`m_1m_3`$) has not been excluded at all. For the pure MSW effect, it is well-known that the supernova neutrino signal with the case of inverted hierarchy would be very different from that with normal hierarchy . From the analogy of the conversion mechanisms between the MSW and RSF effects, it is easily inferred that the RSF conversions will be also very sensitive to the mass hierarchy. Therefore, in this paper, we study three-flavour RSF conversions with the inverted mass hierarchy using the latest neutrino mixing parameters, and give a comprehensive discussion concerning the dependence on the mass hierarchy. Since the RSF conversions are extremely sensitive to the presupernova models as noted above, we first adopt the model of 15$`M_{}`$ with solar metallicity by Woosley and Weaver  as a reference model; the results with the model by Woosley et al with both solar and zero metallicity are addressed later. In particular, we show that the RSF conversion in the case of the inverted hierarchy with the large $`\theta _{13}`$ causes very different neutrino signal from the other models, i.e., the appearance of sharp peak of the neutronization burst in the $`\overline{\nu }_{\mathrm{}}`$ time profile. If this case were realized actually, it would be clearly confirmed by not only the neutrino spectrum but also the luminosity curve.

After the completion of our calculation, a paper appeared in which three-flavour RSF effect was studied with both normal and inverted mass hierarchy . The authors studied rather qualitatively, but they did not obtain the result that the neutronization peak in the $`\overline{\nu }_{\mathrm{}}`$ signal would be detected, when the efficient RSF conversion takes place in the case of the inverted mass hierarchy and large $`\theta _{13}`$; our numerical approach as well as the qualitative discussion clearly indicate that result. Thus, we stress that the point itself is first discovered in the present paper as well as that our study is the first numerical one which comprehensively tackles three-flavour RSF effect with both the normal and inverted mass hierarchy.

This paper is organized as following. In section 2, we give the formulation used in our calculation, which includes all three-flavour neutrinos and antineutrinos. In section 3, models of supernova neutrinos, presupernova structure, and particle properties of the neutrino are illustrated. In section 4, we give qualitative discussions concerning neutrino conversions in supernova matter both for the normal and inverted mass hierarchy, and show results of numerical calculations in section 5. Finally, in section 6, a simple discussion how to obtain information on the neutrino properties from the supernova neutrino signal is presented, and dependence on presupernova models are discussed.

## 2 Formulation

### 2.1 Interaction with matter and magnetic fields

The interaction of the magnetic moment of neutrinos and magnetic fields is described by

$$(\nu _\mathrm{i})_\mathrm{R}|H_{\mathrm{int}}|(\nu _\mathrm{j})_\mathrm{L}=\mu _{\mathrm{ij}}B_{},$$
(1)
where $`\mu _{\mathrm{ij}}`$ is the component of the neutrino magnetic moment matrix, $`B_{}`$ is the magnetic field transverse to the direction of propagation, and $`(\nu )_\mathrm{R}`$ and $`(\nu )_\mathrm{L}`$ are the right- and left-handed neutrinos, respectively. If neutrinos are Dirac particles, right-handed neutrinos and left-handed antineutrinos are undetectable (sterile neutrinos), since they do not interact with matter. On the other hand, if neutrinos are Majorana particles, $`\nu _\mathrm{R}`$ are identical to antiparticles of $`\nu _\mathrm{L}`$ and interact with matter. In this paper, we assume that neutrinos are Majorana particles. The diagonal magnetic moments are forbidden for Majorana neutrinos, and therefore only conversion between different flavours is possible, e.g., $`(\overline{\nu }_{\mathrm{}})_\mathrm{R}(\nu _{\mu ,\tau })_\mathrm{L}`$ or $`(\nu _{\mathrm{}})_L(\overline{\nu }_{\mu ,\tau })_\mathrm{R}`$.

Coherent forward scattering with matter induces an effective potential for neutrinos, which is calculated using weak interaction theory. The effective potential due to scattering with electrons is given by

$$V_{\pm \pm }=\pm \sqrt{2}G_\mathrm{F}\left(\pm \frac{1}{2}+2\mathrm{sin}^2\theta _\mathrm{W}\right)n_{\mathrm{}},$$
(2)
where $`n_{\mathrm{}}`$ is the electron number density, $`G_\mathrm{F}`$ is the Fermi coupling constant, and $`\theta _\mathrm{W}`$ is the Weinberg angle. The $`\pm `$ sign in front refers to $`\nu (+)`$ and $`\overline{\nu }()`$ and that in the parentheses to $`\nu _{\mathrm{}}(+)`$ and $`\nu _{\mu ,\tau }()`$. The difference between $`e`$ and $`\mu ,\tau `$ neutrinos comes from the existence of charged-current interaction. The subscript $`\pm \pm `$ of $`V`$ refers to the first and the second $`\pm `$ sign. The ordinary MSW effect between $`\nu _{\mathrm{}}`$ and $`\nu _{\mu ,\tau }`$ is caused by the potential difference $`V_{\mathrm{}}V_{\mu ,\tau }=V_{++}V_+=\sqrt{2}G_\mathrm{F}n_{\mathrm{}}`$, while that between $`\overline{\nu }_{\mathrm{}}`$ and $`\overline{\nu }_{\mu ,\tau }`$ by $`V_\overline{\mathrm{}}V_{\overline{\mu },\overline{\tau }}=V_+V_{}=\sqrt{2}G_\mathrm{F}n_{\mathrm{}}`$. To include the RSF effect, which causes conversion between neutrinos and antineutrinos, we should take into account the neutral-current scattering by nucleons:

$$V=\sqrt{2}G_\mathrm{F}\left(\frac{1}{2}2\mathrm{sin}^2\theta _\mathrm{W}\right)n_\mathrm{p}\sqrt{2}G_\mathrm{F}\frac{1}{2}n_\mathrm{n},$$
(3)
where $`n_\mathrm{p},n_\mathrm{n}`$ are the proton and neutron number density, respectively. For neutrinos we add $`+V`$ to the potential and for antineutrinos $`V`$. Therefore, the RSF conversion between $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_{\mu ,\tau }`$, which is important for the case considered in this paper, obeys the potential difference

$`\mathrm{\Delta }V`$ $``$ $`V_{\mathrm{}}V_{\overline{\mu },\overline{\tau }}`$ (4)
$`=`$ $`(V_{++}+V)(V_{}V)`$
$`=`$ $`\sqrt{2}G_\mathrm{F}{\displaystyle \frac{\rho }{m_\mathrm{N}}}(12Y_{\mathrm{}}),`$
where $`\rho `$ is the density, $`m_\mathrm{N}`$ is the nucleon mass, and $`Y_{\mathrm{}}=n_{\mathrm{}}/(n_{\mathrm{}}+n_\mathrm{n})`$ is the number of electrons per baryon. (When we obtained equation (4), we assumed charge neutrality $`n_{\mathrm{}}=n_\mathrm{p}`$.)

### 2.2 The simplest case: $`\nu _{\mathrm{}}\overline{\nu }_\tau ^{}`$ conversion

In this subsection, we give the simplest formulation between $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_\tau ^{}`$ with the mass squared difference and mixing angle between the first and third mass eigenstate $`\mathrm{\Delta }m_{13}^2,\theta _{13}`$, and consider the properties of the RSF conversion. Here, $`\overline{\nu }_\tau ^{}`$ represents mass eigenstate in matter which can be obtained by some linear combination of $`\overline{\nu }_\mu `$ and $`\overline{\nu }_\tau `$. Although this treatment is obviously insufficient for a realistic discussion, it is still useful for an intuitive understanding. First, we qualitatively illustrate how the RSF conversions depend on the neutrino mass hierarchy, and then more realistic three-flavour formulation is given in the next subsection. The time evolution of the mixed state of $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_\tau ^{}`$ is described by the Schrödinger equation

$$\mathrm{}\frac{}{r}\left(\begin{array}{c}\nu _{\mathrm{}}\\ \overline{\nu }_\tau ^{}\end{array}\right)=\left(\begin{array}{cc}0& \mu _\mathrm{}\tau ^{}B_{}\\ \mu _\mathrm{}\tau ^{}B_{}& \mathrm{\Delta }H\end{array}\right)\left(\begin{array}{c}\nu _{\mathrm{}}\\ \overline{\nu }_\tau ^{}\end{array}\right),$$
(5)
where $`r`$ is the radius from the center of the star, $`\mu _\mathrm{}\tau ^{}`$ is the transition magnetic moment, and $`\mathrm{\Delta }H`$ is defined by

$$\mathrm{\Delta }H\frac{\mathrm{\Delta }m_{13}^2}{2E_\nu }\mathrm{cos}2\theta _{13}\mathrm{\Delta }V.$$
(6)
The resonance occurs when $`\mathrm{\Delta }H=0`$; it does not occur if $`\mathrm{\Delta }m_{13}^2>0`$ (normal hierarchy), because $`Y_{\mathrm{}}<0.5`$ (therefore, $`\mathrm{\Delta }V<0`$) is satisfied in stellar envelope. Instead, the conversion between $`\overline{\nu }_{\mathrm{}}`$ and $`\nu _\tau `$ is affected by the resonance since the potential difference between them is $`V_\overline{\mathrm{}}V_{\mu ,\tau }=\mathrm{\Delta }V(>0)`$. On the other hand, in the case of the inverted hierarchy, the RSF conversion takes place between $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_\tau ^{}`$, but not between $`\overline{\nu }_{\mathrm{}}`$ and $`\nu _\tau ^{}`$. This situation is very similar to the case of pure MSW effect; higher MSW resonance occurs in antineutrino (neutrino) sector if the mass hierarchy is inverted (normal) .

When the resonance is adiabatic (or when the density is slowly changing at that point and the magnetic field is strong enough), significant conversion occurs. The adiabaticity of the RSF conversion is given by

$$\gamma _{\mathrm{}\overline{\tau }^{}}^{\mathrm{RSF}}=\frac{(2\mu _\mathrm{}\tau ^{}B_{})^2}{|\mathrm{\Delta }V/r|}|_{\mathrm{res}},$$
(7)
where the subscript “res” means that the value is evaluated at the resonance point. Therefore, if the magnetic field is sufficiently strong at the resonance point, the $`\nu _{\mathrm{}}\overline{\nu }_\tau ^{}`$ conversion occurs completely.

### 2.3 Three-flavour formulation

Here, we present the three-flavour (six-component) formulation of neutrino mixing, on which our discussions depend:

$$\mathrm{}\frac{}{r}\left(\begin{array}{c}\nu \\ \overline{\nu }\end{array}\right)=\left(\begin{array}{cc}H_0& B_{}M\\ B_{}M& \overline{H_0}\end{array}\right)\left(\begin{array}{c}\nu \\ \overline{\nu }\end{array}\right),$$
(8)
where

$$\nu =\left(\begin{array}{c}\nu _{\mathrm{}}\\ \nu _\mu \\ \nu _\tau \end{array}\right),\overline{\nu }=\left(\begin{array}{c}\overline{\nu }_{\mathrm{}}\\ \overline{\nu }_\mu \\ \overline{\nu }_\tau \end{array}\right),$$
(9)
$`H_0={\displaystyle \frac{1}{2E_\nu }}U\left(\begin{array}{ccc}0& 0& 0\\ 0& \mathrm{\Delta }m_{12}^2& 0\\ 0& 0& \mathrm{\Delta }m_{13}^2\end{array}\right)U^{}+\left(\begin{array}{ccc}V_{++}+V& 0& 0\\ 0& V_++V& 0\\ 0& 0& V_++V\end{array}\right),`$ (16)
$`\overline{H_0}={\displaystyle \frac{1}{2E_\nu }}U\left(\begin{array}{ccc}0& 0& 0\\ 0& \mathrm{\Delta }m_{12}^2& 0\\ 0& 0& \mathrm{\Delta }m_{13}^2\end{array}\right)U^{}+\left(\begin{array}{ccc}V_+V& 0& 0\\ 0& V_{}V& 0\\ 0& 0& V_{}V\end{array}\right),`$ (23)
$`U=\left(\begin{array}{ccc}U_\mathrm{}1& U_\mathrm{}2& U_\mathrm{}3\\ U_{\mu 1}& U_{\mu 2}& U_{\mu 3}\\ U_{\tau 1}& U_{\tau 2}& U_{\tau 3}\end{array}\right)=\left(\begin{array}{ccc}c_{12}c_{13}& s_{12}c_{13}& s_{13}\\ s_{12}c_{23}c_{12}s_{23}s_{13}& c_{12}c_{23}s_{12}s_{23}s_{13}& s_{23}c_{13}\\ s_{12}s_{23}c_{12}c_{23}s_{13}& c_{12}s_{23}s_{12}c_{23}s_{13}& c_{23}c_{13}\end{array}\right),`$ (30)
$$M=\left(\begin{array}{ccc}0& \mu _\mathrm{}\mu & \mu _\mathrm{}\tau \\ \mu _\mathrm{}\mu & 0& \mu _{\mu \tau }\\ \mu _\mathrm{}\tau & \mu _{\mu \tau }& 0\end{array}\right),$$
(31)
and $`c_{\mathrm{ij}}=\mathrm{cos}\theta _{\mathrm{ij}},s_{\mathrm{ij}}=\mathrm{sin}\theta _{\mathrm{ij}}`$. (We assume the CP phase $`\delta =0`$ in equation (30) for simplicity.)

The resonant flavour conversion basically occurs when the two diagonal elements in the matrix in equation (8) have the same value. There are four relevant resonance points, but they depends on the neutrino mass hierarchy as discussed in the previous subsection. In the case of the normal mass hierarchy, they are for $`\nu _{\mathrm{}}\nu _\mu ^{}`$ (MSW-L), $`\nu _{\mathrm{}}\nu _\tau ^{}`$ (MSW-H), $`\overline{\nu }_{\mathrm{}}\nu _\mu ^{}`$ (RSF-L), and $`\overline{\nu }_{\mathrm{}}\nu _\tau ^{}`$ (RSF-H). Here, the quantities such as $`\nu _{\mu ,\tau }^{}(\overline{\nu }_{\mu ,\tau }^{})`$ represent mass eigenstates in matter which can be obtained by linear combination of $`\nu _\mu (\overline{\nu }_\mu )`$ and $`\nu _\tau (\overline{\nu }_\tau )`$. This is because the $`\nu _{\mathrm{}}(\overline{\nu }_{\mathrm{}})`$ state coincides with mass eigenstate in matter owing to large matter potential, and the other mass eigenstates are obtained by the rotation of $`\nu _\mu (\overline{\nu }_\mu )`$ and $`\nu _\tau (\overline{\nu }_\tau )`$ basis. Suffixes ‘-L’ and ‘-H’ attached to ‘MSW’ and ‘RSF’ indicate whether the density at the resonance points are lower or higher. On the other hand, in the case of the inverted hierarchy, RSF-H occurs for $`\nu _{\mathrm{}}\overline{\nu }_\tau ^{}`$ as well as MSW-H for $`\overline{\nu }_{\mathrm{}}\overline{\nu }_\tau ^{}`$, whereas the other two resonances occurs for the same conversions. This situation is summarized in table 1.

## 3 Models

### 3.1 Original neutrino emission

We adopt, as original neutrino spectrum as well as luminosity curve, the result of numerical simulation by Thompson et al; we use the model calculated for the $`15M_{}`$ progenitor star. Their calculation has particularly focused on shock breakout and followed the dynamical evolution of the cores through collapse until the first 250 ms after bounce. They have incorporated all the relevant neutrino processes such as neutrino–nucleon scatterings with nucleon recoil as well as nucleon bremsstrahlung; these reactions have recently been recognized to give non-negligible contribution to the spectral formation. In figures 1 and 2, we show the original luminosity curve and number spectrum of neutrinos, respectively. In these figures, $`\nu _\mathrm{x}`$ represents non-electron neutrinos and antineutrinos.

Neutrino luminosity curve is quite characteristic among different flavours. In particular, there is a very sharp peak of $`\nu _{\mathrm{}}`$ called ‘neutronization burst’, whose duration is typically $`10`$ ms and peak luminosity is $`10^{53}`$ erg/s. This strong peak is illustrated as follows. As a supernova shock moves outward, it dissociates nuclei into free nucleons, which triggers deleptonization process $`\mathrm{}^{}\mathrm{p}\nu _{\mathrm{}}\mathrm{n}`$; these $`\nu _{\mathrm{}}`$ build up a sea because they are trapped and advected with matter. When the shock crosses the $`\nu _{\mathrm{}}`$ neutrinosphere, within which the created $`\nu _{\mathrm{}}`$ are trapped, they are abruptly emitted. For the other flavours $`\overline{\nu }_{\mathrm{}}`$ and $`\nu _\mathrm{x}`$, there is no such a sudden burst; both luminosities glow rather gradually and they are similar to each other. Usually, for both the pure MSW effect and the RSF effect with normal hierarchy, the most easily detected flavour $`\overline{\nu }_{\mathrm{}}`$ are transformed from original $`\overline{\nu }_{\mathrm{}}`$ and $`\nu _\mathrm{x}`$ . Thus, luminosity curve does not provide any useful information on the flavour conversion mechanism. On the other hand, for the RSF effect with inverted mass hierarchy, the conversion $`\nu _{\mathrm{}}\overline{\nu }_\tau ^{}\overline{\nu }_{\mathrm{}}`$ is considered to occur via RSF-H and MSW-H (see table 1); we expect that this case can be distinguished from the luminosity curve.

The other characteristic that provides information on flavour conversion mechanism is hierarchy of the average energy $`E_\nu _{\mathrm{}}<E_{\overline{\nu }_{\mathrm{}}}<E_{\nu _\mathrm{x}}`$ as clearly seen from figure 2; flavour conversions also cause the spectral exchange. This energy hierarchy is explained as follows. Since $`\nu _\mathrm{x}`$ interact with matter only through the neutral-current interactions in supernovae, they are weakly coupled with matter compared to $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_{\mathrm{}}`$. Thus the neutrinosphere of $`\nu _\mathrm{x}`$ locates at deeper in the core than that of $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_{\mathrm{}}`$, which leads to higher temperatures for $`\nu _\mathrm{x}`$. The difference between $`\nu _{\mathrm{}}`$ and $`\overline{\nu }_{\mathrm{}}`$ comes from the fact that the core is neutron-rich and $`\nu _{\mathrm{}}`$ couple with matter more strongly, through $`\nu _{\mathrm{}}\mathrm{n}\mathrm{}^{}\mathrm{p}`$ reaction.

### 3.2 Presupernova profiles

We use the precollapse model of massive stars of Woosley and Weaver . The model is $`15M_{}`$ progenitor star with solar metallicity and it is labeled as W95S. The density and $`Y_{\mathrm{}}`$ profiles are quite important for the flavour conversions because they determine the resonance regions as well as whether it is adiabatic or not. We show in figure 3 the $`|\rho (12Y_{\mathrm{}})|`$ (responsible for RSF) and $`\rho Y_{\mathrm{}}`$ profiles (responsible for MSW) of the W95S model. We also show $`\mathrm{\Delta }_{12}m_\mathrm{N}\mathrm{\Delta }m_{12}^2\mathrm{cos}2\theta _{12}/2\sqrt{2}G_\mathrm{F}E_\nu `$ and $`\mathrm{\Delta }_{13}m_\mathrm{N}|\mathrm{\Delta }m_{13}^2|\mathrm{cos}2\theta _{13}/2\sqrt{2}G_\mathrm{F}E_\nu `$ as two horizontal bands (the bandwidth comes from the energy range 5–70 MeV). At intersections between $`\mathrm{\Delta }_{12},\mathrm{\Delta }_{13}`$ and $`|\rho (12Y_{\mathrm{}})|,\rho Y_{\mathrm{}}`$, the RSF and MSW conversions take place.

In the previous publication, we have investigated the dependence on presupernova models ; in the paper, we have also used recent $`15M_{}`$ presupernova models by Woosley et al with both solar and zero metallicity (W02S and W02Z, where ‘S’ and ‘Z’ denotes solar and zero metallicity, respectively). This is because the RSF conversion is very sensitive to the deviation of $`Y_{\mathrm{}}`$ from 0.5 (see equation (4)), which strongly depends on the metallicities as well as on the weak interaction rates adopted in the simulation of stellar evolution. In fact, for the W02 models, in which the authors adopted more updated weak rates, the $`|\rho (12Y_{\mathrm{}})|`$ profiles suddenly drop at RSF region, yielding rather nonadiabatic conversions; these properties are shown in figure 4. Nevertheless, we use the W95S model as our reference model because of the reasons listed below. First, since the conversion probabilities in RSF-H is highly adiabatic when the magnetic field is sufficiently strong (section 5), it is easier for us to investigate the dependence on mass hierarchy. Second, although the W02 models were calculated using a recent shell model and resulted in substantial revisions to the older data sets in the W95S model, there are still many uncertainties concerning nuclear physics; thus, these models cannot be considered to be decisive one. Finally, as discussed in the previous paper , because the lifetime of massive stars which end their life by gravitational collapse is much shorter than that of the Sun, the progenitors of the galactic supernovae are considered to be younger than the Sun, and therefore, more metal rich. Because the deviation of $`Y_{\mathrm{}}`$ from 0.5 is determined by rarely existent nuclei, the suppression of $`(12Y_{\mathrm{}})`$ will be weaker and the RSF conversion will incline to be more adiabatic, even if W02S model is the correct one of the progenitor with solar metallicity.

Although we use static progenitor models in calculating the flavour transition, in fact the density profile changes drastically during a neutrino burst ($`10`$ s) owing to shock wave propagation. However, our particular interest is within first 0.25 s after core bounce because the calculation of original neutrino emission  ends at that time (3.1). In the previous paper, we have already shown that until 0.25 s, using the static presupernova and magnetic field models is considered to be a good approximation ; this is based on the numerical calculation of the Lawrence Livermore group, which is the only group succeeding the shock propagation into the outer envelope.

### 3.3 Magnetic fields

We assume that the global structure of the magnetic field is a dipole moment and the field strength is normalized at the surface of the iron core with the values $`10^{10}`$ G (nearly complete RSF conversion) as well as 0 G (pure MSW conversion). The reason for this normalization is as follows. The magnetic fields should be normalized by fields that are static and exist before the core collapse, because those of a nascent neutron star can hardly affect the far outer region, where the RSF conversions take place, within the short time scale of a neutrino burst. As discussed in the previous subsection, since the shock wave does not affect the resonance region at $`0.25`$ s after bounce, it is also expected that the magnetic field structure and strength at the resonance points are not seriously changed at that time. The strength of such magnetic fields above the surface of the iron core may be inferred from observations of the surface of white dwarfs, since both are sustained against gravitational collapse by the degenerate pressure of electrons. Observations of the magnetic fields in white dwarfs show that the strength spreads in a wide range of $`10^7`$$`10^9`$ G . Considering the possibility of the decay of magnetic fields in white dwarfs, it is not unnatural to consider magnetic fields up to $`10^{10}`$ G at the surface of the iron core. Then, in equation (8), $`B_{}=B_0(r_0/r)^3\mathrm{sin}\mathrm{\Theta }`$, where $`B_0`$ is the strength of the magnetic field at the equator on the iron core surface, $`r_0`$ the radius of the iron core, and $`\mathrm{\Theta }`$ the angle between the pole of the magnetic dipole and the direction of neutrino propagation. Hereafter, we assume $`\mathrm{sin}\mathrm{\Theta }=1`$.

### 3.4 Neutrino parameters

We adopt the realistic neutrino mixing parameters inferred from the recent experimental results: for the atmospheric neutrino parameters, $`|\mathrm{\Delta }m_{13}^2|=2.8\times 10^3\mathrm{eV}^2,\mathrm{sin}^22\theta _{23}=1.0`$, and for the solar neutrino parameters, $`\mathrm{\Delta }m_{12}^2=5.0\times 10^5\mathrm{eV}^2,\mathrm{tan}^2\theta _{12}=0.42`$.

As for still uncertain neutrino properties, we must set some assumptions, e.g., whether the mass hierarchy is normal or inverted, as well as $`\theta _{13}`$ is large, enough for the MSW-H conversion to be adiabatic, or not. There is also uncertainty concerning the neutrino magnetic moment tensor $`\mu _{\mathrm{ij}}`$; however, since the only relevant parameter is $`\mu _{\mathrm{ij}}B`$, this uncertainty is already included in that of the magnetic field strength ($`\mu _{\mathrm{ij}}=10^{12}\mu _B`$ is assumed). Thus, there are 8 parameter sets due to: whether the magnetic field is zero (label by ‘MSW’) or sufficiently strong $`B_0=10^{10}`$ G (‘RSF’); the mass hierarchy is normal (‘NOR’) or inverted (‘INV’); and $`\mathrm{sin}^22\theta _{13}=10^6`$ (‘S’) or $`\mathrm{sin}^22\theta _{13}=0.04`$ (‘L’). We label one model by connecting these labels using hyphen, e.g., MSW-NOR-S. It should be particularly noticed that although we label the model with strong magnetic field simply by RSF-, it does not mean that the pure RSF effect occurs; every model labeled by RSF- is subject to both the MSW and RSF conversions. We summarize these models in table 2.

## 4 Qualitative conversion schemes

We qualitatively illustrate the conversion scheme for each model. Since $`\overline{\nu }_{\mathrm{}}`$ is the most easily detected flavour, we focus on the conversion channel that gives $`\overline{\nu }_{\mathrm{}}`$ appearance at the detector. Figure 5 schematically shows crossings among different mass eigenstates in matter, which is helpful for the reader to understand the qualitative discussions in this section.

### 4.1 Normal mass hierarchy

First, the case of the normal mass hierarchy is addressed (figure 5(a)); this case is minutely investigated in , but we again repeat their discussion. When $`B_0=0`$, the produced $`\overline{\nu }_{\mathrm{}}`$ propagate the supernova envelope without experiencing the MSW resonance, and reach the stellar surface as the lightest mass eigenstate $`\overline{\nu }_1`$. The other mass eigenstates $`\overline{\nu }_2,\overline{\nu }_3`$ are originated from $`\overline{\nu }_{\mu ,\tau }`$ whose flux is considered to be the same at the leading order. Thus, the $`\overline{\nu }_{\mathrm{}}`$ flux at the detector is described by

$`F_{\overline{\nu }_{\mathrm{}}}`$ $`=`$ $`|U_\mathrm{}1|^2F_{\overline{\nu }_1}+|U_\mathrm{}2|^2F_{\overline{\nu }_2}+|U_\mathrm{}3|^2F_{\overline{\nu }_3}`$ (32)
$`=`$ $`|U_\mathrm{}1|^2F_{\overline{\nu }_{\mathrm{}}}^0+(1|U_\mathrm{}1|^2)F_{\nu _\mathrm{x}}^0,`$
where $`F^0`$ and $`F`$ are the flux at production and detection, respectively. Since the parameter $`\theta _{13}`$ is essential for the MSW-H conversion, the results are not sensitive to the value of $`\theta _{13}`$.

On the other hand, when $`B_0=10^{10}`$ G, RSF-H becomes almost completely adiabatic, whereas RSF-L is highly nonadiabatic. At the RSF-H point, the original $`\nu _\tau ^{}`$ are converted into $`\overline{\nu }_{\mathrm{}}`$, which propagate as the mass eigenstate owing to the large matter potential. These $`\overline{\nu }_{\mathrm{}}`$, then, cross the nonadiabatic RSF-L region (which gives no essential effects) and escape from the star as $`\overline{\nu }_1`$. The other mass eigenstates are also originated from $`\nu _\mathrm{x}`$, thus yielding

$$F_{\overline{\nu }_{\mathrm{}}}=F_{\nu _\mathrm{x}}^0.$$
(33)
The value of $`\theta _{13}`$ does not matter also in this case. In consequence, the models MSW-NOR-S,L are characterized by equation (32), whereas, RSF-NOR-S,L are characterized by equation (33).

### 4.2 Inverted mass hierarchy

The situation changes dramatically in the case of the inverted mass hierarchy (figure 5(b)). The pure MSW effect in this case has been studied in detail in the literatures ; the value of $`\theta _{13}`$ is critical for completeness of the conversions. The $`\overline{\nu }_{\mathrm{}}`$, which are produced as the lightest mass eigenstate in three antineutrino states, then cross the MSW-H resonance region. If this resonance is nonadiabatic (MSW-INV-S), the flavour conversion does not take place at the MSW-H point, i.e., $`\overline{\nu }_{\mathrm{}}\overline{\nu }_1`$. This case yields the $`\overline{\nu }_{\mathrm{}}`$ flux characterized by equation (32), which has already appeared. On the other hand for the MSW-INV-L model, MSW-H is completely adiabatic, and therefore the conversions $`\overline{\nu }_{\mathrm{}}\overline{\nu }_\tau ^{}\overline{\nu }_3`$ as well as $`\nu _\mathrm{x}\overline{\nu }_{1,2}`$ occur, which results in

$$F_{\overline{\nu }_{\mathrm{}}}=(1|U_\mathrm{}3|^2)F_{\nu _\mathrm{x}}^0+|U_\mathrm{}3|^2F_{\overline{\nu }_{\mathrm{}}}^0.$$
(34)
Since the value of $`\theta _{13}`$ (or $`|U_\mathrm{}3|^2`$) is strongly constrained to be very small by the reactor experiment , this expression essentially the same as the previous one (33).

When $`B_0=10^{10}`$ G (RSF-INV-S,L), the $`\nu _{\mathrm{}}\overline{\nu }_\tau ^{}`$ conversion takes place at the RSF-H region. This $`\overline{\nu }_\tau ^{}`$ further enter the MSW-H region. If this MSW resonance is nonadiabatic (-S model), these $`\overline{\nu }_\tau ^{}`$ becomes $`\overline{\nu }_3`$ at the stellar surface. The other relevant conversions are the same as the pure MSW model $`\overline{\nu }_{\mathrm{}}\overline{\nu }_1,\overline{\nu }_\mu ^{}\overline{\nu }_2`$. Thus, the $`\overline{\nu }_{\mathrm{}}`$ flux is given by

$$F_{\overline{\nu }_{\mathrm{}}}=|U_\mathrm{}1|^2F_{\overline{\nu }_{\mathrm{}}}^0+|U_\mathrm{}2|^2F_{\nu _\mathrm{x}}^0+|U_\mathrm{}3|^2F_\nu _{\mathrm{}}^0,$$
(35)
which is almost the same as equation (32). On the other hand, for -L model, the mass eigenstates $`\overline{\nu }_\tau ^{}`$ converted from $`\nu _{\mathrm{}}`$ reach the supernova surface as $`\overline{\nu }_1`$, because the adiabatic MSW-H conversion $`\overline{\nu }_\tau ^{}\overline{\nu }_{\mathrm{}}`$ takes place. With the remaining channels $`\overline{\nu }_\mu ^{}\overline{\nu }_2,\overline{\nu }_{\mathrm{}}\overline{\nu }_3`$, the observed flux is

$$F_{\overline{\nu }_{\mathrm{}}}=|U_\mathrm{}1|^2F_\nu _{\mathrm{}}^0+|U_\mathrm{}2|^2F_{\nu _\mathrm{x}}+|U_\mathrm{}3|^2F_{\overline{\nu }_{\mathrm{}}}^0.$$
(36)
Consequently, the models MSW-INV-S and RSF-INV-S are characterized by equation (32), MSW-INV-L is by (33) and RSF-INV-L is by (36).

In table 3, we summarize the relevant conversions stated above for each model. Although there are eight models and this rather large number is expected to complicate the discussion, each model is consequently categorized into one of three groups A, B or C, as indicated in table 3. Fortunately, this greatly reduces the complexity, but within one group, we cannot specify each model by the SK observation.

## 5 Results of numerical calculations

### 5.1 Conversion probability

We calculated equation (8) numerically with the models given in section 3, and obtained the conversion probabilities for each flavour. Figure 6(a) shows the conversion probabilities of original $`\nu _\mu `$ as a function of radius for the RSF-NOR-L model; the neutrino energy is 25 MeV. Because the $`\nu _\mu `$ state is not mass eigenstate and maximally mixing with $`\nu _\tau `$, the probability $`P(\nu _\mu \nu _\mu )`$ should oscillate abruptly. In our calculation, however, we have averaged out this local behaviour, since we are only interested in the global change of the probabilities which is concerned with resonance. Therefore, the $`\nu _\mu `$ survival probability starts from 0.5 rather than 1, because of the maximal mixing with $`\nu _\tau `$. A significant amount of $`\nu _\mu `$ change into $`\overline{\nu }_{\mathrm{}}`$, owing mainly to the RSF-H conversion which occurs around $`0.01R_{}`$, and the converted $`\overline{\nu }_{\mathrm{}}`$ propagates as the mass eigenstate due to matter effect, not being disturbed by further resonances, to $`0.1R_{}`$. Around this radius, $`\overline{\nu }_{\mathrm{}}`$ start to mix with other flavour antineutrinos, reducing the probability $`P(\nu _\mu \overline{\nu }_{\mathrm{}})`$. In figure 6(b), we show the conversion probabilities of original $`\nu _{\mathrm{}}`$ for the RSF-INV-L model. As already discussed in 4.2, the original $`\nu _{\mathrm{}}`$ are transformed into $`\overline{\nu }_{\mu ,\tau }`$ at RSF-H, and they further change to the favoured flavour $`\overline{\nu }_{\mathrm{}}`$ in the MSW-H resonance point.

We show in figure 7 the same probabilities as those shown in figure 6, but the assumed value of $`\theta _{13}`$ is small ((a) RSF-NOR-S; (b) RSF-INV-S). Figure 7(a) indicates that the conversion probabilities from the original $`\nu _\mu `$ are almost the same as those for RSF-NOR-L model; as we have already noted, the value of $`\theta _{13}`$ does not matter in the case of normal mass hierarchy. On the other hand, the behaviours in figure 7(b) is substantially different from those in figure 6(b). In particular, the most easily detected flavour $`\overline{\nu }_{\mathrm{}}`$ is not produced from the original $`\nu _{\mathrm{}}`$. All these characteristics are consistent with the simple discussions given in section 4.

### 5.2 Neutrino signals at the Super-Kamiokande detector

With the conversion probabilities given in the previous subsection and the original neutrino spectrum by Thompson et al, we calculated the flux of each flavour neutrinos on the Earth. (From this point on, we assume the galactic supernova neutrino burst and take 10 kpc as a distance to the supernova.) Using this flux and the cross section of the relevant neutrino interaction at SK as well as the sensitivity of the detector, we can calculate the expected event numbers from a future galactic supernova neutrino burst. In this paper, we adopt the most dominant reaction $`\overline{\nu }_{\mathrm{}}\mathrm{p}\mathrm{}^+\mathrm{n}`$ alone; a cross section of the reaction has been calculated in detail .

Figure 8(a) shows time evolution of the energy-integrated event as a function of time for each group A, B and C, and we show the same in figure 8(b) but using equally spaced bins. From the time evolution of neutrino events, we cannot discern the groups A and B, because they show almost the same time profile. For group C, however, since the original $`\nu _{\mathrm{}}`$ are converted into $`\overline{\nu }_{\mathrm{}}`$, the time profile shows steep neutronization peak, and the event number contained in this peak is expected to be statistically significant as clearly seen in figure 8(b); the event number included in the most prominent three bins is $`180`$. If the neutronization peak were actually obtained, it strongly indicates that the model group C would be favoured; since the group C contains only one model, RSF-INV-L, a great number of problems concerning the neutrino properties would be solved at the same time. In that case, the neutrino would have the nonzero magnetic moment, the mass hierarchy would be inverted, and the value of $`\theta _{13}`$ would be large enough for MSW-H to be adiabatic.

The expected event number per unit energy range, which is integrated during the first 0.25 s after core bounce, is shown in figures 8(c) and 8(d), in units of counts/MeV and counts/bin, respectively. From these figures, the model group C gives the softest spectrum, while B the hardest and A an intermediate one. In addition to the time evolution of the neutrino events, the number spectrum would provide useful information on the flavour conversion mechanism. Although the available data are restricted in order to avoid uncertainties concerning shock wave propagation, the obtained data would be statistically significant. Using the spectrum, degeneracy between the group A and B is expected to be broken.

## 6 Discussion

### 6.1 How far can we probe the neutrino properties from the supernova neutrino observation?

The expected neutrino signal at the SK detector has been investigated thus far. However, the mechanism of supernova explosions is quite unclear, since all the reliable numerical simulations have not succeeded in pushing the shock wave to penetrate the entire core. There may be several unknown processes which we have omitted so far, and the original neutrino spectrum as well as its luminosity curve are still controversial. Therefore, we cannot trust characteristics of the calculation by Thompson et al in detail. Instead, we use rather simple quantities in order to discuss the conversion mechanisms from the neutrino signals; this approach is expected to considerably reduce the dependence on supernova models.

We adopt the following quantities:

$`R_{\mathrm{SK}}^\mathrm{E}`$ $`=`$ $`{\displaystyle \frac{\text{Number of events for }E_{\mathrm{}}>25\mathrm{MeV}}{\text{Number of events for }E_{\mathrm{}}<15\mathrm{MeV}}},`$ (37)
$`R_{\mathrm{SK}}^\mathrm{T}`$ $`=`$ $`{\displaystyle \frac{\text{Number of events for }0<t/\mathrm{ms}<75}{\text{Number of events for }75<t/\mathrm{ms}<150}},`$ (38)
in order to represent the spectral hardness and the peak sharpness of neutronization burst, respectively. Since the time of the core collapse would never known with the neutrino signal alone, we take time origin $`t=0`$ when the first neutrino signal is detected. The place of the model groups A, B and C on the ($`R_{\mathrm{SK}}^\mathrm{E},R_{\mathrm{SK}}^\mathrm{T}`$) plane are shown in figure 9. These groups are well separated from each other, and we expect that this particular remains unchanged even if the adopted supernova model is different.

Although the degeneracy problem within each group cannot be solved by the SK observation, which mainly detect $`\overline{\nu }_{\mathrm{}}`$, it may be possible if $`\nu _{\mathrm{}}`$ could be detected efficiently. SNO is such a detector currently data taking with 1 000 tons of heavy water. The supernova $`\nu _{\mathrm{}}`$ can be detected via $`\nu _{\mathrm{}}e^{}\mathrm{pp}`$ reaction. Although the $`\overline{\nu }_{\mathrm{}}`$ are also detected through a similar reaction, $`\overline{\nu }_{\mathrm{}}e^+\mathrm{nn}`$, these events could be discriminated using delayed coincidence technique. Thus at SNO, we expect that only $`\nu _{\mathrm{}}`$ signal can be extracted. Figure 10 is the same as figure 9 for the SNO detector, but is plotted for the models in (a) group A and (b) B. As shown in this figure, the degeneracy problem could be solved with the $`\nu _{\mathrm{}}`$ data in principle, but the current smallness of the SNO detector prevents a statistically significant discussion.

### 6.2 Dependence on presupernova models

Until this point, we have adopted the presupernova model W95S as our reference model. However, as already stated in 3.2 or investigated in , the RSF conversions are highly dependent on the presupernova profiles. In fact, figures 3 and 4 show that the relevant profile for the RSF effect, $`\rho (12Y_{\mathrm{}})`$, is quite different among these models as well as their metallicities. On the other hand, the MSW conversions are insensitive to the presupernova models because their relevant profile is $`\rho Y_{\mathrm{}}`$ that is not subject to the deviation of $`Y_{\mathrm{}}`$ from 0.5. In particular for the W02 models, the value of $`\rho (1Y_{\mathrm{}})`$ suddenly drops at the RSF-H region, yielding rather nonadiabatic conversions. Since RSF becomes incomplete but MSW does not change for these models, the model group B, which contains RSF-NOR-S, RSF-NOR-L and MSW-INV-L models, should be further divided into two subgroups; we define group B I containing two RSF models, while MSW-INV-L is named group B II. Consequently, we have four characteristic model groups A, B I, B II and C; in particular, the groups B II and C contain only one specific model each. Within each group, we cannot identify each model from the SK observation.

Figures 11(a) and (b) show time evolution of neutrino events at SK, calculated with the W02S and W02Z models. Comparing these figures with figure 8(a), the neutronization peak is not as sharp as that for the W95S model. This is because the RSF-H conversion is not efficient owing to sudden drop of the $`\rho (12Y_{\mathrm{}})`$ profiles. However, since that peak is still prominent and statistically significant, the model group C would be discriminated from the other groups using the time profile detected at SK.

Figures 11(c) and (d) show the number spectrum of the neutrino signal, which is evaluated for the W02S and W02Z models, respectively. The spectra for groups A and B II are unchanged from the previous calculation with the W95S model (figure 8(c)), because they are essentially described by the pure MSW effect, which is insensitive to presupernova models. Therefore, we focus on the other two groups B I and C. For the group B I, the hardness of spectra locates between that of group A and B II, and it strongly depends on the adopted metallicities. More detailed discussions concerning this group have been already given in our previous paper , and we refer the reader to the literature. For the group C, the low energy peak disappears from the spectrum calculated with the W02S model, reflecting the $`\rho (12Y_{\mathrm{}})`$ profile of the W02S model (figure 4(a)), which significantly reduces the efficiency of the RSF-H conversions. On the other hand, this low energy peak remains for the W02Z model, also reflecting the $`\rho (1Y_{\mathrm{}})`$ profile. As clearly shown in figure 4(b), it mildly changes in the RSF-H region for low energy neutrinos. However, if the energy is beyond some critical value, RSF-H becomes completely nonadiabatic, because the $`\rho (12Y_{\mathrm{}})`$ profile abruptly drops to zero. Thus, the $`\nu _{\mathrm{}}\overline{\nu }_{\mathrm{}}`$ conversion is highly efficient for low energy neutrinos but is highly inefficient for high energy ones; for the high energy region the $`\nu _\mathrm{x}\overline{\nu }_{\mathrm{}}`$ conversion is relevant, which leads to double peak profile of the spectrum.

## 7 Conclusion

In this paper, we investigated the RSF conversions in supernovae for both the normal and inverted mass hierarchy. As the case for the pure MSW effect, we found that the RSF transitions are strongly dependent on the neutrino mass hierarchy, and also on the value of $`\theta _{13}`$. We first gave qualitative discussion on the neutrino conversions including both the RSF and MSW effect, for eight parameter sets summarized in table 2. From that consideration, it was found that these models are categorized into only three groups, each of which is expected to show characteristic neutrino signal at the SK detector; we named these groups A, B and C, and which parameter set is included in each group is summarized in table 3.

We, then, presented results of numerical calculations of flavour conversions in supernova envelope. The density and $`Y_{\mathrm{}}`$ profiles of the presupernova star calculated by Woosley and Weaver  (W95S) was adopted for the calculations. As the magnetic field structure, we assumed dipole-type and normalized its strength at the surface of the iron core. Using the conversion probabilities calculated by such a procedure and the original supernova neutrino spectrum as well as luminosity curve given by Thompson et al, the expected neutrino signal at the SK detector was estimated. As the result, it was found that there are clear difference between the model groups both for the spectral shape and the time evolution of the neutrino events. In particular, the model group C, which include the RSF-INV-L model alone, shows a sharp neutronization peak. Therefore, if this peak were detected in reality from the future galactic supernova neutrino burst, it would strongly support the RSF-INV-L model. It would indicate that many problems concerning the neutrino property should be solved at the same time, i.e., the neutrino have the nonzero magnetic moment, the mass hierarchy is inverted, and the value of $`\theta _{13}`$ is large enough to induce the adiabatic MSW-H resonance. From the spectral shape, it would be possible to discern the groups A and B, however, within each model group we need other information such as the signal of $`\nu _{\mathrm{}}`$, to discriminate one model from the others. Although the SNO detector could detect $`\nu _{\mathrm{}}`$ efficiently, current smallness of the detector would prevent the statistically significant discussions.

We also studied the dependence on presupernova models. As already investigated in our previous paper , the RSF conversion is strongly dependent on the deviation of $`Y_{\mathrm{}}`$ from 0.5, which is quite sensitive to the metallicities as well as the weak interaction rates adopted in the simulation of stellar evolution. The presupernova models by Woosley et al (W02S and W02Z) shows the considerably different $`\rho (12Y_{\mathrm{}})`$ profile from that of the W95S model. For both models, there exists sudden drop of the $`\rho (12Y_{\mathrm{}})`$ profile at the RSF-H region, which pushes the RSF conversion rather nonadiabatic, and the expected neutrino signal was found to be different from that estimated with the W95S model. Still, since the neutronization peak for the RSF-INV-L model exists at the statistically significant level, it would keep to be useful method also in these cases.

SA’s work is supported by Grant-in-Aid for JSPS Fellows. KS’s work is supported in part by Grant-in-Aid for Scientific Research provided by the Ministry of Education, Culture, Sports, Science and Technology of Japan through Research Grant No S14102004, Grant-in-Aid for Scientific Research on priority Areas No 14079202.

## References


format: text

text:
# How to apply importance-sampling techniques to simulations of optical systems

C. J. McKinstrie and P. J. Winzer

Bell Laboratories, Lucent Technologies, Holmdel, New Jersey 07733

Abstract

This report contains a tutorial introduction to the method of importance sampling. The use of this method is illustrated for simulations of the noise-induced energy jitter of return-to-zero pulses in optical communication systems.

When one rolls a die, the probability of it stopping with a particular number face-up is 1/6. To measure such a probability experimentally and accurately, one must roll the die many more than 6 times. With this fact in mind, consider the design of optical communication systems. A common way to evaluate system prototypes is to measure (in computer simulations or laboratory experiments) the distance over which they can transmit information with bit-error-ratios of order $`10^9`$ (without forward error correction). Even with current computers, it is impractical to simulate pulse transmission many more than $`10^9`$ times.

One can use importance-sampling techniques to circumvent this difficulty. Put simply, bit errors are caused by large system perturbations (deviations from ideal behavior) that occur infrequently. When one makes importance-sampled simulations of system performance one biases them in such a (controlled) way that these large perturbations occur more often than they should. Because they occur often, one can measure their statistics accurately. Subsequently, one adjusts the simulation results to remove their artificial bias. In this way one obtains results that are both unbiased and accurate.

Pioneering importance-sampled simulations of phase noise in an optical system were made by Foschini and Vannuci . Recently, importance-sampled simulations of polarization-mode dispersion were made by Biondini et al. . Importance-sampling methods were reviewed comprehensively by Smith et al.  and ways to combine data from different importance-sampled simulations were discussed by Veach . In this report we make no attempt to duplicate their discussions. Rather, we show by example how to apply importance-sampling techniques to optical systems.

It is instructive to consider the die example quantitatively. Each roll of an unbiased die is a trial, in which the probability of a successful outcome $`p=1/6`$ and the value associated with a successful outcome (the weight with which a successful outcome is counted) $`w=1`$. For each trial the expected value of the outcome is $`pw`$ and the variance is $`w^2p(1p)`$. Suppose that an unbiased die is rolled $`n_1`$ times and let $`f_1`$ be the measured 1-frequency (the number of 1s divided by the number of rolls). Then the expected value of the 1-frequency $`E(f_1)=1/6`$. By adding $`n_1`$ terms of magnitude $`5/36`$ and dividing the result by $`n_1^2`$ one finds that the variance $`V(f_1)=5/36n_1`$: To measure the 1-frequency with an accuracy of 1% (to conduct an experiment for which the standard deviation of the 1-frequency is 1% of the expected value) one would have to roll the die $`5\times 10^4`$ times. This would be a time-consuming task.

If the 5- and 6-frequencies were of no interest one could bias the die by marking the 5- and 6-faces with 1s. Since this change would increase the probability of rolling a 1 ($`p=3/6=1/2`$), one should count each successful roll with reduced weight ($`w=1/3`$). Suppose that this biased die is rolled $`n_2`$ times and let $`f_2`$ be the measured (weighted) 1-frequency. Then the expected value $`E(f_2)=1/6`$, as required, and the variance $`V(f_2)=1/36n_2`$: By causing the desired event to occur more often than it would naturally, one is able to measure its probability more accurately, or with fewer rolls. (For rare events that occur in applications, the performance improvements are much larger than the factor of 5 associated with this example.) Of course, the price one pays for this increase in accuracy is the loss of information about the 5- and 6-frequencies.

If these frequencies were of limited, but finite, interest, one could combine the measurements made with each die separately. If one were to weight the two measurements equally, by defining the combined 1-frequency $`f=(f_1+f_2)/2`$, one would find that $`E(f)=1/6`$ and $`V(f)=5/144n_1+1/144n_2`$: The accuracy of the combined measurement is limited by the less-accurate individual measurement (the one with the larger variance coefficient or the one made with fewer rolls). One should weight the individual measurements in proportion to the numbers of rolls used to make them, and in inverse proportion to their variance coefficients. Let $`f=c_1f_1+c_2f_2`$, where the condition $`c_1+c_2=1`$ ensures that $`E(f)=1/6`$. Then a short calculation shows that the optimal values of $`c_1`$ and $`c_2`$ are $`n_1/(n_1+5n_2)`$ and $`5n_2/(n_1+5n_2)`$, respectively, in which case the minimal variance $`V(f)=5/36(n_1+5n_2)`$. For the case in which $`n_1=n_2=n/2`$, where $`n`$ is the total number of rolls, $`V(f)=5/108n`$.

In this dice example simple formulas for the individual variances exist, which allow one to determine the optimal weight coefficients precisely. However, in applications such formulas might be complicated or unknown. Consequently, a different method is required. An effective method, which is called the balance heuristic , is to weight each successful outcome (roll of either die) equally. In this method, if the unbiased die is rolled $`n_1`$ times and the biased die is rolled $`n_2`$ times, the combined probability of rolling a 1 is $`(n_1+3n_2)/6(n_1+n_2)`$. To ensure that the expected value of the combined measurement is 1/6, each successful roll (of either die) should be counted with weight $`w=(n_1+n_2)/(n_1+3n_2)`$. By adding $`n_1`$ contributions of magnitude $`w^2p_1(1p_1)`$ and $`n_2`$ contributions of magnitude $`w^2p_2(1p_2)`$, and dividing the result by $`(n_1+n_2)^2`$, one finds that the combined variance is $`(5n_1+9n_2)/36(n_1+3n_2)^2`$. For the case in which $`n_1=n_2=n/2`$ the combined variance is $`7/144n`$, which is only 5% larger than the minimal variance of the preceding paragraph. Thus, one should weight each outcome in inverse proportion to its combined probability.

Now consider the noise-induced energy (amplitude) jitter of a return-to-zero (RZ) pulse. Provided that the pulse energy $`e`$ is greater than the noise energy in the surrounding bit slot, it evolves according to the stochastic ordinary differential equation (ODE)

$$d_ze=(g\alpha )e+r,$$
(1)
where $`g(z)`$ is the amplifier gain rate, $`\alpha `$ is the fiber loss rate and $`r(z)`$ is the rate at which the energy is changed by amplifier noise . This random rate of change is quantified by the equations $`r(z)=0`$ and $`r(z)r(z^{})=(2n_{\mathrm{sp}}h\nu ge)\delta (zz^{})`$, where $``$ denotes an ensemble average, $`n_{\mathrm{sp}}`$ is the spontaneous-emission factor (1.1–1.3) and $`h\nu `$ is the photon energy. Equation (1) is valid for any isolated pulse and an arbitrary combination of distributed and lumped amplification.

For definiteness, consider a 10 Gb/s system with uniformly-distributed amplification ($`g=\alpha `$), in which $`\alpha =0.21`$ dB/Km, $`\beta =0.30`$ ps2/Km (D = 0.38 ps/Km-nm) and $`\gamma =`$ 1.7/Km-W. Then a soliton with a full-width at half-maximum of 30 ps has an energy of 21 fJ (time-averaged power of 0.21 mW). If the system length $`l=10`$ Mm the output noise power in both polarizations, in a frequency bandwidth of 12 GHz (wavelength bandwidth of 0.1 nm), is 1.7 $`\mu `$W: The (optical) signal-to-noise ratio is 21 dB. Systems with nonuniformly-distributed or lumped amplification produce the same noise power in shorter distances.

For uniformly-distributed amplification Eq. (1) can be rewritten in the canonical form

$$dx=x^{1/2}dy,$$
(2)
where $`x=e/e_0`$ is the energy, normalized to the equilibrium energy (in the absence of noise), and $`y`$ is a Wiener process (Gaussian random variable) with $`y=0`$ and $`y^2=\sigma _s^2z`$, where the normalized source-strength $`\sigma _s^2=2n_{\mathrm{sp}}h\nu g/e_0`$. In the linear regime the multiplicative factor $`x^{1/2}1`$, from which it follows that $`x1+y`$: The probability-density function (PDF) of the output energies is Gaussian, with mean 1 and variance $`\sigma _s^2z`$. (From a logical standpoint the PDF of the non-negative quantity $`x`$ cannot be exactly Gaussian, because, if it were, the probability of $`x<0`$ would be finite for all $`z>0`$. From a practical standpoint this inconsistency is tolerable if the probability of $`x<0`$ is exponentially small for system lengths of interest.) For the aforementioned system the output variance $`\sigma _s^2l=6.6\times 10^3`$ (which is of order $`10^2`$) and the output deviation is $`8.1\times 10^2`$ (which is of order $`10^1`$). In the nonlinear regime the factor $`x^{1/2}`$ modifies the tails of the PDF significantly. For reference, the analytical solution of Eq. (2) has the PDF

$$P(x)\frac{\mathrm{cosh}(mx^{1/2}/v)\mathrm{exp}[(m^2+x)/2v]}{(2\pi xv)^{1/2}},$$
(3)
where $`m=1\sigma _s^2z/8`$ and $`v=\sigma _s^2z/4`$ .

Equation (2) and solution (3) model energy jitter in a (continuous) system with uniformly-distributed amplification. We simulated a (discrete) system with $`n_i=100`$ lumped amplifiers. Between the amplifiers the energy $`x`$ did not change. At the $`i`$th amplifier the energy was changed (kicked) by the random amount $`\delta y_i`$, where the properties $`\delta y_i=0`$ and $`\delta y_i\delta y_j=10^4\delta _{ij}`$ ensured that $`(_{i=1}^{n_i}\delta y_i)^2=10^2`$: The discrete system had the same characteristics as the continuous system. The output energies were assigned to energy bins of (common) width 0.02 and each bin probability $`p_j`$ was defined to be the number of pulses whose energies fell within the bin boundaries (bin count) divided by the total number of pulses. (Since probabilities cannot be measured by finite numbers of trials, these quantities should be called the relative frequencies associated with the bins. We use the term probabilities as an abbreviation for the correct term.) To facilitate comparisons to the analytical PDF (3), the simulation probabilities were defined to be the bin probabilities divided by the bin width. In these (direct) simulations the occurrence of each output energy was counted with unit weight.

The PDF associated with an ensemble of $`10^6`$ pulses is displayed in Fig. 1. Although the simulation results agree well with Eq. (3) near the peak of the PDF, they do not even begin to sample the tail of the PDF. On a 1-GHz PC these simulations, which are based on Eq. (2), take a few minutes. Simulating the transmission of many more than $`10^9`$ pulses would take many days. Realistic simulations, which are based on the nonlinear Schroedinger equation, would take even longer.

To probe the tails of the PDF one must make large energy perturbations occur more often than they would naturally. One way to achieve this goal is to increase the standard deviation of the energy kicks. Let $`q_0`$ denote the (common) unbiased kick distribution, with deviation $`\sigma _0=10^2`$, and $`q`$ denote the (common) biased kick distribution, with deviation $`\sigma >\sigma _0`$. Then, at the $`i`$th amplifier the probability that the energy is kicked by the amount $`\delta y_i`$ is increased by the factor $`q(\delta y_i)/q_0(\delta y_i)`$. Since the kicks at all the amplifiers are biased, the output energy occurs with a probability that is larger than its natural probability by the total factor $`f_t=\mathrm{\Pi }_{i=1}^{n_i}q(\delta y_i)/q_0(\delta y_i)`$, which depends on the full kick sequence $`(\delta y_1,\delta y_2,\mathrm{},\delta y_{n_i})`$. One can remove this bias by counting the output energy with reduced weight: One increments the appropriate bin probability by $`1/f_t`$, rather than 1. All other aspects of data counting remain the same. For reference, the probability factor $`1/f_t`$ is called the likelihood ratio.

The PDF associated with an ensemble of $`10^6`$ pulses is displayed in Fig. 2 for the case in which $`\sigma =1.2\sigma _0`$. The results of these (importance-sampled) simulations differ from the previous results in two ways. First, they do probe the tail of the PDF. Although the simulations reproduce the shape of the analytical PDF, the simulation probabilities are not accurate because the number of data points that sample the tail of the PDF is still small. Second, by causing large kicks to happen more often, one causes small kicks to happen less often. Consequently, the body (peak) of the PDF is not reproduced accurately. This deficiency prevents one from increasing $`\sigma `$ until the tail of the PDF is sampled accurately.

Another way to achieve the stated goal is to change the (common) mean of the kick distributions. If the mean kick $`\mu `$ is positive (negative) the energy drifts toward larger (smaller) values. The distributions associated with 3 ensembles of $`3\times 10^5`$ pulses are displayed in Fig. 3 for cases in which $`\sigma =\sigma _0`$, and $`\mu =3\times 10^3`$, 0 and $`3\times 10^3`$. These values produce data sets with mean energies of $`0.3`$, 0.0 and 0.3, respectively. Although the simulation distributions are inaccurate for energies that are far from their mean energies (because the bin counts are low), they are accurate near their mean energies. It only remains to combine the individual distributions to produce a composite distribution that is accurate for the entire domain of interest.

One way to combine the individual distributions is to weight their bin probabilities equally (which is equivalent to combining the data sets before sorting the output energies and the associated probability factors into bins). Let $`p_{jk}`$ be the $`j`$th bin probability associated with the $`k`$th data set (which was produced by kick distributions with mean $`\mu _k`$). For each $`j`$, if the individual bin probabilities $`p_{jk}`$ were all zero the combined bin probability $`p_j`$ was defined to be zero and if some of the $`p_{jk}`$ were nonzero $`p_j`$ was defined to be their average. The results of this procedure are shown in Fig. 4. Although the composite distribution does cover the domain of interest, it is inaccurate near the boundaries of the sample spaces. Combining the individual distributions with equal weight allows the bodies of the distributions (which have high bin counts) to be polluted by the tails of neighboring distributions (which have low bin counts).

The dice example suggests that it is better to weight the bin probabilities according to the associated bin counts. For each $`j`$, if the individual bin counts $`b_{jk}`$ were all zero the combined bin probability $`p_j`$ was defined to be zero and if some of the bin counts were nonzero $`p_j`$ was defined to be $`_{k=1}^{n_k}b_{jk}p_{jk}/_{k=1}^{n_k}b_{jk}`$. The results of this procedure are shown in Fig. 5. The composite distribution covers, and is accurate throughout, the domain of interest. This ad-hoc procedure allows one to combine bin probabilities generated at different times, without recourse to the data sets on which they were based.

Although the ad-hoc method works, it is not the balance heuristic. Suppose that the kick sequence $`(\delta y_1,\delta y_2,\mathrm{},\delta y_{n_i})`$ occurs during the first simulation, which is made using the biased distribution $`q_1`$ (with mean $`\mu _1`$). Then the associated probability factor $`f_{t1}=\mathrm{\Pi }_{i=1}^{n_i}q_1(\delta y_i)/q_0(\delta y_i)`$. Were the same sequence to occur during the $`k`$th simulation, the associated probability factor would be $`f_{tk}`$, which depends on the biased distribution $`q_k`$. If the simulations involve the same number of pulses, the combined probability factor $`f_t=_{k=1}^{n_k}f_{tk}/n_k`$. (It is easy to generalize this formula.) When we made the individual simulations and sorted the output energies and the individual probability factors into bins, we also sorted the combined probability factors into a separate set of bins, which was common to all the simulations. The results of this procedure are shown in Fig. 6. The balance-heuristic method works well.

In summary, we showed by example how to apply importance sampling techniques to simulations of optical communication systems. These techniques are easy to apply and increase significantly the accuracy with which rare events can be simulated.

We acknowledge useful discussions with D. Chizhik, G. Foschini, R. Moore and J. Salz.

Postscript: Independent simulations of energy jitter were made recently by Moore et al. \[Opt. Lett. 28, 105 (2003)\], who showed that the predictions of the energy equation (1) are consistent with the results of simulations based on the nonlinear Schroedinger equation.


text:
# Group velocity study in hot 87Rb vapor with buffer gas

## I Introduction.

Electromagnetically induced transparency (EIT) has been the subject of numerous theoretical and experimental studies Brandt et al., (1997); Vanier et al., (1998); Erhard and Helm, (2001); Lukin et al., (1997); Akulshin et al., (1998, 1991); Javan et al., (2002); Taichenachev et al., (2000). Materials displaying EIT have many interesting properties, such as narrow resonance width and steep dispersion resulting in ultra-low light group velocity Hau et al., (1999); Kash et al., (1999); Budker et al., 1999a . These properties make EIT interesting for several practical applications, such as frequency standards Knappe et al., (2001), precision magnetometry Novikova et al., (2001), enhanced nonlinear optics Hakuta et al., (1991); Hemmer et al., (1995); Jain et al., (1996); Zibrov et al., (1999), and quantum information storage Phillips et al., (2001); Liu et al., (2001); Zibrov et al., (2002).

The group velocity of light in media can be expressed as Kocharovskaya et al., (2001)

$$v_g=\frac{d\omega }{dk}=\frac{c\omega \frac{n(\omega ,k)}{k}}{n(\omega ,k)+\omega \frac{n(\omega ,k)}{\omega }}=\stackrel{~}{v}_gv_s$$
(1)
where $`c`$ is the speed of light in vacuum, $`\omega `$ is the frequency of the field, $`n`$ is the index of refraction in the medium, $`k`$ is the vacuum wave number, and we have defined

$$\stackrel{~}{v}_g=\frac{c}{n(\omega ,k)+\omega \frac{n(\omega ,k)}{\omega }}$$
(2)
and

$$v_s=\frac{\omega \frac{n(\omega ,k)}{k}}{n(\omega ,k)+\omega \frac{n(\omega ,k)}{\omega }}.$$
(3)
Here and below we use the convention that terms with a tilde ( $`\stackrel{~}{}`$ ) denote values in the moving reference frame and terms with no tilde denote values in the laboratory frame.

The first term, $`\stackrel{~}{v}_g`$ , is due to frequency dispersion and the second term, $`v_s`$ , is due to spatial dispersion. Because of spatial dispersion, the group velocity is different for atoms with different speeds $`v_a`$ . This is easily seen because in the moving frame $`\stackrel{~}{\omega }=\omega kv_a`$ (we take $`v_a`$ positive for an atom moving in the same direction as the light propagation), from which we see that $`v_g=\frac{d\omega }{dk}=\stackrel{~}{v}_g+v_a`$. This is just just the Galilean transformation from the moving frame of the atoms to the laboratory frame. Thus, a mono-velocity atomic beam moving in the opposite direction from the light propagation direction slows the group velocity.

In a dense Doppler broadened medium, it is possible to obtain slower or even zero group velocity Kocharovskaya et al., (2001); Rostovtsev et al., (2002) for the probe field in a $`\mathrm{\Lambda }`$ configuration of strong drive and weak probe fields when both are red detuned but maintain the two-photon resonance condition. The following conditions on the power of the drive field must be satisfied

$$\mathrm{\Omega }\sqrt{\gamma \gamma _{cb}}$$
(4)
$$\mathrm{\Omega }k_dv_T\sqrt{\frac{\gamma _{cb}}{\gamma }}$$
(5)
where $`\mathrm{\Omega }`$ is the Rabi frequency of the drive laser and $`k_d`$ is its wave vector, $`v_T`$ is the average thermal velocity, $`\gamma `$ is the decay rate of the upper level, and $`\gamma _{cb}`$ is the decay rate of coherence between the lower (ground) levels.

Equation (4) is just the usual condition for EIT for individual atoms. Equation (5) is applicable only for the case where the drive field is weak enough that EIT occurs only for a narrow spread of atomic velocities. In this case, the intensity of the drive field is not large enough to pump all atomic velocity subgroups into the dark state. This means that the optical pumping rate $`|\mathrm{\Omega }|^2\gamma /\mathrm{\Delta }^2`$, for atoms having one-photon detuning $`\mathrm{\Delta }`$, is less than the relaxation rate $`\gamma _{bc}`$, between levels $`b`$ and $`c`$. Therefore $`|\mathrm{\Omega }|^2\gamma /\mathrm{\Delta }^2<\gamma _{bc}`$ implies that EIT does not occur for all moving atoms and the light interacts with a quasi atomic beam. When these conditions are satisfied, one can choose a velocity sub-group of atoms with a particular average velocity (in the direction of light propagation) and narrow velocity spread $`\mathrm{\Delta }v\frac{\mathrm{\Omega }}{k_d}\sqrt{\gamma _{cb}/\gamma }`$. This is accomplished by changing the one-photon detuning of the drive laser field while maintaining two-photon resonance. In this case, the center of the quasi-beam of moving atoms is determined by the simple Doppler relation

$$v_a=c\frac{\mathrm{\Delta }}{\omega _d}$$
(6)
where $`\omega _d`$ is the drive laser frequency. Naturally we would expect atoms moving with the atomic quasi-beam to increase the group velocity, and atoms moving in the opposite direction to “drag” light with them, or decrease the total group velocity.

The intuitive picture described above is quantified rigorously in Ref. Kocharovskaya et al., (2001). In the present paper, we experimentally study these theoretical predictions by measuring the dependence of the group velocity on probe-field one-photon detuning for different experimental conditions.

A common method of increasing the dispersion in an EIT medium is to lengthen the ground-state coherence lifetime, thereby decreasing the linewidth of the EIT resonance. The coherence lifetime is often limited by the interaction time of the atoms with the lasers. Common methods for increasing this lifetime are by introducing a buffer gas to confine the atoms Brandt et al., (1997); Wynands and Nagel, (1998); Erhard et al., (2000); Erhard and Helm, (2001), using wall coatings in the cell so that coherence is preserved between successive interactions of the atoms and the lasers Kanorsky et al., (1995); Budker et al., (1998); Budker et al., 1999b ; Robinson et al., (1958); Bouchiat and Brossel, (1966); Alexandrov et al., (2002), and by cooling and trapping the atoms Lipsich et al., (2001); Chen et al., (1998).

We have found that for EIT conditions in a sample with buffer gas (with linewidth on the order of several kHz and group velocity on the order of a few tens to hundreds of meters/sec) the probe field has a slower group velocity when it is blue detuned with respect to resonance and higher group velocity for red detuning. This result is opposite to the intuitive picture described above, and to that of Kocharovskaya et al. in Kocharovskaya et al., (2001).

## II Experiment

### II.1 Setup

A schematic of the experimental setup is shown in Figure 1. One external cavity diode laser (ECDL) is used as the source of a strong driving field and another as a weak probe. The drive and probe lasers are combined and pass through a cell containing isotopically enhanced 87Rb vapor. The 87Rb density is varied by changing the temperature of the cell. A three layer magnetic shield (MS) screens out the laboratory magnetic field.

The lasers are phase locked to each other with a frequency offset that is tunable about the ground level hyperfine splitting of 87Rb (6835 MHz). The lasers are spatially mode-matched with a single-mode optical fiber. The drive laser is tuned to the $`5S_{1/2}(F=2)5P_{1/2}(F^{}=2)`$ transition of 87Rb and the probe laser is tuned to the $`5S_{1/2}(F=1)5P_{1/2}(F^{}=2)`$ transition as shown in Figure 2.

The configuration of drive and probe lasers shown in Figure 2 is called a $`\mathrm{\Lambda }`$ configuration. In this case, the lasers optically pump all atoms in the desired velocity subgroup into a dark state superposition of the ground levels, giving rise to strong coherence between these lower levels. Scattering of the drive field on this coherence results in the generation of a fairly strong Stokes component (new field) in the medium (shown in Figure 2). The generation of this new field near two-photon resonance of the drive and probe field is described in Ref. Lukin et al., (1997); Kash et al., (1999); Mikhailov et al., (2002).

The transmitted drive and probe fields and the generated Stokes field are detected separately by heterodyne detection on a fast photo detector. This is done by splitting off part of the drive field from the main beam before entering the cell. This component is shifted up in frequency by a small amount (60 MHz), and combined with the beams exiting the cell before the photo-detector. The photo-current thus contains beat signals at various RF frequencies, separated by 60 MHz, in the vicinity of the 6835 MHz separation of the drive and probe. By separately analyzing these components with a spectrum analyzer (SA), we extract the transmitted probe and generated new fields independently. This technique is described in Ref. Mikhailov et al., (2002); Kash et al., (1999).

### II.2 Time delay and group velocity measurement procedure.

We extract the group velocity in the medium by modulating the intensity of the probe field before the cell, and observing the time delay before this modulation is observed in the transmitted field. Experiments have been conducted with gaussian shaped (temporal) pulses and with sinusoidal modulation. We find that the delay time is independent of the probe field modulation technique as long as the bandwidth of the modulation does not exceed the transmission linewidth of the EIT resonance.

The modulation is generated by use of a frequency generator (FG) which drives an acousto-optical (AO) modulator in the probe laser. The deflected beam from the AO is blocked, so the AO serves as a probe power attenuator (PA in Figure 1). We also measure the time delay of the generated new field.

With a sinusoidal modulation of the probe, lock-in detection of the transmitted probe field provides a sensitive measure of the time delay due to the medium. When using lock-in detection in this way, we obtain the time delay from the phase shift of the transmitted probe field intensity relative to the probe intensity before the cell. For a sine wave of frequency $`f`$, this phase shift is given by

$$\psi =2\pi \tau _df\psi _0,$$
(7)
where $`\tau _d`$ is the time delay introduced by the atoms and $`\psi _0`$ is a phase shift introduced by electronics. To eliminate the unknown $`\psi _0`$ we measure the phase shift for several different modulation frequencies. The phase shift increases linearly with frequency $`f`$ and the slope of this line is $`2\pi \tau _d`$. We extract $`\tau _d`$ with a least squares fit. We then find the group velocity $`v_g`$ by setting $`v_g=L/\tau _d`$ where $`L`$ is the length of the cell. ($`L1\mathrm{cm}`$ for our experiment.)

Another experimental technique to simplify the group velocity measurement is to replace the second laser and phase-lock circuitry with an electro-optic modulator (EOM) in the drive laser beam. By applying a narrow tunable microwave signal to the EOM, upper and lower sidebands are generated at the microwave frequency, which we choose to match the 6.835 GHz ground state hyperfine splitting. We tune the laser to the drive transition, so that the upper sideband drives the probe transition and the lower sideband is off resonance. We choose the microwave amplitude to generate a sideband with power that is 1/10 of the drive power. A careful comparison of the two methods (two phased locked lasers versus one modulated laser) shows no difference in the group velocity measurements. We note that a similar technique is to modulate the laser current, which directly creates sidebands on the laser, has also been successfully employed Affolderbach et al., (1999).

We measure the dependence of $`v_g`$ on one-photon laser detuning $`\mathrm{\Delta }=\omega _d\omega _{22}`$ , where $`\omega _d`$ is the frequency of the drive laser and $`\omega _{22}`$ is the frequency of the $`F=2F^{}=2`$ transition. During each such measurement, the frequency difference of the probe and drive lasers is kept constant and equal to ground level splitting (6.835 GHz). The drive laser power is $`300\mu W`$, and probe power is $`3\mu W`$. These measurements may then be repeated for different 87Rb densities (different temperatures of the cell).

### II.3 Experimental results

We first consider the case where no buffer gas is used, and the ground-state coherence lifetime is limited by the free-flight transit time of the thermal rubidium atoms through the laser beam. In our experiment, the EIT transmission linewidth is 30 kHz and the resulting group velocity on the order of 10 km/s. Figure 3 shows the dependence of the group velocity as a function of the drive laser frequency in the vicinity of the drive resonance. The group velocity is too large to observe the spatial dispersion effect described in the introduction. Another way to see this is to note that the group velocity is much higher than the mean thermal speed of the atoms. In Figure 3 the drive detuning spans the full range of the upper-state hyperfine splitting, and the small feature on the left of the spectrum is the result of generating EIT on the upper $`F^{}=1`$ level (see Figure 2).

Next we narrow the transmission linewidth with a buffer gas. This increases the dispersion considerably, resulting in reduced group velocity. Figure 4 shows the group velocity as a function of drive laser detuning for similar conditions as Figure 3 but with the addition of 3 torr of $`\mathrm{N}_2`$ buffer gas. The EIT transmission linewidth is only a few kHz and we see that the group velocity has fallen to below 100 m/s. We also see that increasing the laser beam diameter increases the dispersion and reduces the group velocity as expected.

Figure 5 shows how the group velocity depends on drive laser detuning for increasing atomic density. We observe a lowering of the group velocity for higher density, but in no case do we observe a lowering of the group velocity as the drive laser is detuned red of resonance, which we would expect based on the prediction of dragging slow light by atoms moving opposite to the laser propagation direction Kocharovskaya et al., (2001). On the contrary, we see the group velocity increase for negative drive laser detuning, and a minimum group for one-photon detuning about 100 MHz blue of one-photon resonance (see Figure 5).

The explanation of this effect is still not clear, but it is plausible that this counter-intuitive behavior is caused by velocity changing collisions in the presence of the buffer gas. We note that when the beam diameter decreases, the group velocity dependence versus one-photon detuning becomes more shifted to the side of negative detuning (see Figure 4). In other words, the behavior becomes less counter-intuitive. We explain this by noting that for a very narrow beam, fewer velocity changing collisions take place before the atom leaves the laser beam. In any case, it is clear that the discussion of Ref. Kocharovskaya et al., (2001) about a quasi mono-velocity beam is not applicable in the case with a buffer gas, since all velocity groups are mixed by velocity changing collisions.

When the density of 87Rb atoms is increased ($`10^{12}\mathrm{cm}^3`$) highly nonlinear interaction of the drive and probe fields leads to very efficient generation of a Stokes component, or new field Lukin et al., (1997); Kash et al., (1999). We can measure the intensity of the generated field as a function of two-photon detuning and find the width to be greater than the transmission width of the EIT resonance of the probe field. (Under the conditions in our experiment it is roughly a factor of two wider Mikhailov et al., (2002).) Correspondingly, we also measure the time delay between modulation of the probe field before the cell and the resulting modulation of the new field after the cell. For the rubidium cell with buffer gas, this delay time is smaller than for the probe field, meaning that the group velocity of the new field ($`v_n`$) is greater than the group velocity of the probe field ($`v_p`$). This is no great surprise since the new field is propagating far from one-photon resonance.

As discussed above, as the rubidium density is increased the probe field group velocity decreases. Similar behavior occurs for the generated new field for low density. However, for large atomic density the group velocity for the new field starts to increase with density. These results are shown in Figure 6. This dependence is connected with a propagation effect. The effective generation of new field occurs in the part of the cell where the group velocity is small. As the light propagates through the cell, the drive field intensity decreases until the new field is decoupled from the probe. After this point, the new field propagates at nearly the vacuum speed of light. Thus the observed average speed of new field increases with atomic density.

We can also measure the new-field group velocity as a function of drive laser detuning. The results are shown in Figure 7. We find that the the group velocity is smaller for negative one-photon detuning than for positive detuning in vicinity of resonance. This behavior follows the intuitive predictions of Ref. Kocharovskaya et al., (2001).

## III Summary.

We observe counter-intuitive dependence of the probe field group velocity versus drive field one-photon detuning for different densities (temperatures) of 87Rb (see Figure 5). The group velocity decreases slightly for blue-detuned drive fields and increases slightly for red-detuned fields. We conclude that the the predictions of Ref. Kocharovskaya et al., (2001) cannot be applied to the case where the EIT linewidth is reduced with a buffer gas, since all velocities are constantly mixing via velocity changing collisions. Unfortunately, without a buffer gas we cannot achieve narrow enough EIT resonances to reach the amount of dispersion needed to get group velocities low enough to observe the effect of dragging the light by atoms.

We find that the group velocity is higher for the generated Stokes field, and that the behavior as a function of detuning is opposite that of the probe field.

## IV Acknowledgement

We thank A. B. Matsko, I. Novikova, V. Sautenkov, and M. O. Scully for useful and stimulating discussions. This work was supported by the Office of Naval Research.


text:
# Magnetic phase diagram of the antiferromagnetic pyrochlore Gd2Ti2O7

## Abstract

Gd2Ti2O7 is a highly frustrated antiferromagnet on a pyrochlore lattice, where apart from the Heisenberg exchange the spins also interact via dipole-dipole forces. We report on low-temperature specific heat measurements performed on single crystals of Gd2Ti2O7 for three different directions of an applied magnetic field. The measurements reveal the strongly anisotropic behaviour of Gd2Ti2O7 in a magnetic field despite the apparent absence of a significant single-ion anisotropy for Gd3+. The $`HT`$ phase diagrams are constructed for $`H[111]`$, $`H[110]`$ and $`H[112]`$. The results indicate that further theoretical work beyond a simple mean-field model is required.

Amongst the members of the titanium pyrochlore oxide family with the general formula R2Ti2O7, the magnetic properties of which are all presently the subject of intense investigations, gadolinium titanium oxide (GTO) holds a unique position. In contrast to the other members of this series of compounds, where magnetic anisotropy plays a crucial role and results in the appearance of some most unusual ground state configurations, including the “spin-ice” state Bramwell\_Science01 , GTO is the only compound where the magnetic rare earth ions have a negligible single-ion anisotropy. This is because in Gd2Ti2O7, the Gd3+ ions are in a state with $`S=7/2`$ and $`L=0`$. The magnetic ions are antiferromagnetically coupled. Therefore, apart from the ever present dipole-dipole interactions, GTO may be regarded as a realisation of an ideal Heisenberg antiferromagnet on a frustrated pyrochlore lattice.

The considerable theoretical interest in the 3$`d`$ pyrochlore model and its 2$`d`$ analogue (a checkerboard lattice), is mostly due to the fact that these models, both in classical and quantum descriptions, do not develop conventional magnetic order down to $`T=0`$ K Pyro\_General . Their ground state is instead described as a degenerate spin-liquid Spin\_Liquid . Whether the weaker perturbations such as a single ion anisotropy, dipolar interactions Dipolar ; Dipolar\_Enjalran , small lattice distortions Distortions , further-neighbor interactions NNN  or a magnetic field Field  break the degeneracy and force the system to adopt a particular ground state configuration is the subject of current investigations. In this context it is extremely important to identify different magnetic phases in real magnetic materials and to establish the phase diagrams.

Another unusual aspect of the thermodynamics of frustrated magnets, which may have direct consequences for GTO, has been pointed out recently by Zhitomirsky MZH\_PRB . Because a significant enhancement of the magnetocaloric effect is expected for frustrated magnets near the saturation field, GTO may have potential as an adiabatic demagnetisation refrigerator. An enhancement of the magnetocaloric effect in model systems is attributed to the presence of a significant number of so-called “soft modes” MZH\_PRB . The existance of such soft (and zero) modes in a magnetic material such as Gd2Ti2O7, still needs to be established. To date, the only direct experimental observation of soft modes has been made with inelastic neutron scattering in another frustrated compound ZnCr2O4 ZnCrO .

In this letter, we present the results of specific heat measurements on Gd2Ti2O7 in an applied magnetic field. By using single crystal samples we have been able to prove that there is an unusually high degree of anisotropy built into the ground state and to construct the $`HT`$ phase diagrams for three different directions of magnetic field, $`H[111]`$, $`H[110]`$ and $`H[112]`$.

The single crystals of Gd2Ti2O7 were grown by the floating zone technique, using an infrared image furnace Balakrishnan . The principal axes of the samples were determined using X-ray diffraction Laue photographs; the crystals were aligned to within an accuracy of 2-4. The specific heat measurements were performed in the temperature range 0.39 to 5.0 K in a field of up to 9 T using a Quantum Design PPMS calorimeter equipped with a 3He option. In order to estimate the lattice contribution to the specific heat we have measured the specific heat of a single crystal of Y2Ti2O7, a nonmagnetic compound isostructural to GTO. For all the low-temperature measurements reported in this letter ($`T<5`$ K) the lattice contribution was found to be insignificant.

We have measured the specific heat both as a function of temperature in a constant magnetic field and as a function of applied field at constant temperature. Because of the relatively high values of the specific heat of GTO and its poor thermal conductivity at low temperatures, the measurements were restricted to small ($`<1`$ mg) plate-like samples with a typical thickness of 0.1 mm. For some orientations of applied field, the demagnetizing factors were significant and had to be taken into account demag . This has been done by measuring the magnetization of the same or similarly shaped samples with either a vibrating sample or a SQUID magnetometer.

The temperature dependence of the heat capacity divided by temperature, $`C(T)/T`$, measured in zero field is shown in Fig. 1. Two sharp peaks of nearly the same amplitude observed at $`T_{N1}`$=1.02(2) K and $`T_{N2}`$=0.74(2) K mark the temperature of two consecutive phase transitions. In Fig. 1 two data sets are shown in order to demonstrate the high degree of reproducibility in these measurements. These data should be compared with the results of previous measurements performed on powder samples Raju ; Ramirez . Raju et al. Raju  reported a broad peak centred around 2 K (absent in our data) and a sharp peak just below 1 K, which were attributed to the presence of short range correlations and the development of long range magnetic order respectively. More recently, however, Ramirez et al. Ramirez  have reported the presence of another sharp peak at around 0.7 K. While our data are similar to those presented in Ref. Ramirez , they differ significantly in the absolute values of $`C(T)`$ from the values reported in Ref. Raju . For example, at the temperature of the upper phase transition, $`T_{N1}`$=1.02 K in our data, the specific heat reaches a value of 22 $`\mathrm{Jmol}^1\mathrm{K}^1`$, while the values reported in Ref. Raju  are approximately half this quantity. Moreover, Bonville et al. Bonville  have very recently repeated the measurements of specific heat on a powder sample of GTO and reported values which are in a good agreement with our data.

In terms of the overall magnetic entropy, Raju et al. observed a significant magnetic specific heat up to about 30 K, while Ramirez et al. have recovered 92% of the expected entropy, $`2R\mathrm{ln}(2S+1)=34.6\mathrm{Jmol}^1\mathrm{K}^1`$, between 0 and 5 K. The differences in the specific heat data, as well as in the $`ac`$ magnetic susceptibility behaviour, have been suggested to arise from the extreme sensitivity of the low-temperature properties to the sample preparation technique Ramirez . The inset in Fig. 1 shows the temperature dependence of the magnetic entropy, obtained by extrapolating the $`C(T)/T`$ to 0 K and numerically integrating it versus temperature. Approximately 90% of the expected value is recovered when the temperature reaches 5 K. Additional high temperature measurements have shown that the entropy is almost fully recovered by 10 K, a temperature which corresponds to the reported Curie-Weiss temperature of GTO Raju . At that point the magnetic contribution to the specific heat becomes negligible when compared to the lattice contribution. The fact that a majority of entropy is recovered only at $`T>`$$`T_{N1}`$ signifies the importance of short-range correlations in GTO.

Powder neutron diffraction measurements performed on GTO sample by Champion et al. in zero field Champion , have revealed a complex noncollinear structure with a magnetic propagation vector of $`𝐤=(\frac{1}{2}\frac{1}{2}\frac{1}{2})`$. According to these measurements, the GTO magnetic structure (the so-called P state) consists of kagome planes with “$`q=0`$”-type order with interstitial sites carrying zero magnetic moment; GTO has been suggested to be only partially ordered even at 50 mK Champion . The observed recovery of the full magnetic entropy conflicts with this suggestion. On the other hand, mean-field theory Dipolar\_Enjalran  predicts that an ordered four sublattice structure (F state) competes with the P state and becomes more stable as the temperature is reduced. Therefore we suggest that at $`T_{N1}`$ GTO undergoes a transition from a paramagnetic to the P state, while at $`T_{N2}`$ the transition is to the F state.

Our lowest experimentally accessible temperature of 0.39 K is about 50% of $`T_{N2}`$. We cannot therefore reliably extrapolate the $`C(T)`$ dependence down to 0 K. From the data shown in Fig. 1, however, it is quite clear that between 0.39 K and $`T_{N2}`$, $`C(T)/T`$ grows approximately linearly with temperature, which implies a $`T^2`$ dependence for the low-temperature part of $`C(T)`$. A similar $`T^2`$ dependence has been observed in another magnetically frustrated system, $`\mathrm{SrCr}_{9\mathrm{p}}\mathrm{Ga}_{129\mathrm{p}}\mathrm{O}_{19}`$ (SCGO) with a kagome lattice SCGO .

The temperature and field dependencies of the specific heat, measured on the single crystals of GTO for the three different directions of applied magnetic field, $`H[111]`$, $`H[112]`$ and $`H[110]`$ are shown in Figures 2 and 3. The application of a small external field causes the convergence of the transition temperatures, $`T_{N1}`$ and $`T_{N2}`$, regardless of the field orientation. In a field of about 1.5 T, the two transitions merge coincide , as seen on the left-hand panels of Fig. 2. In a small field and low-$`T`$, $`C(T)`$ still follows a $`T^2`$ dependence; above 1.5 T, however, there is no $`T^2`$ dependence in $`C(T)`$. This is in contrast with SCGO, whose specific heat shows almost no field dependence SCGO . For fields above 1.5 T, significant differences in the positions of the peaks in the $`C(T)`$ curves for GTO are clearly visible (see right-hand panels of Fig. 2).

The field dependence of the specific heat shown in Fig. 3 has been measured in order to determine precisely the positions of the phase transitions lines at low $`T`$, where in high fields the specific heat is nearly temperature independent. The higher-field peaks observed in these $`C(H)`$ curves at around 5 to 6 T, mark a transition to a saturated (paramagnetic) phase. The low-$`T`$ peaks found at approximately half the saturation field correspond to an additional field-induced transition discussed below. The most intense peaks in $`C(H)`$ are found in the higher-$`T`$ low-field region. They reach values above 70 $`\mathrm{Jmol}^1\mathrm{K}^1`$ and correspond to a transition into an intermediate temperature phase. It is interesting to compare this $`C(H)`$ data with the results of Monte Carlo simulations for classical spins on the pyrochlore lattice near the saturation field MZH\_PRB , where the presence of a macroscopic number of soft modes has been predicted. This comparison is not straightforward, as in the simulations only the exchange interactions have been included, which precludes the system from entering a state with long range magnetic order MZH\_PRB . Fig. 3 shows that for $`H>H_{sat}`$ the measured field dependence of the specific heat does not follow a simple exponential decay expected for a nonfrustrated antiferromagnet. Instead it remains almost field-independent for $`H<7`$ T and then falls rapidly at higher fields. This observation may suggest a large influence of the soft modes on the specific heat of GTO around the saturation field.

No signs of history dependence have been detected throughout the course of our measurements, e.g. the values of the measured specific heat were identical for increasing/decreasing field and temperature. This suggests a second order nature for all the phase transitions in GTO in contrast with the observations of Bonville et al. Bonville , which imply the upper phase transition at $`T_{N1}`$ has a weak first order character .

The magnetic $`HT`$ phase diagrams for GTO obtained by plotting the temperature of various phase transitions observed in both the $`Cvs.H`$ and $`Cvs.T`$ measurements are shown in Fig. 4. By analysing these phase diagrams, the highly anisotropic nature of the magnetic behaviour of GTO becomes more apparent. To start with, when $`H[111]`$, the magnetic field favors the onset of magnetic order in GTO. For example, in a field of about 4 T the transition temperature is 1.14(2) K compared to 1.02(2) K in zero field, while the transition to a saturated phase measured at 0.57 K occurs at about 5.7 T. When $`H[112]`$, an external field inhibits the onset of long range magnetic order: at 4 T the transition is reduced to 0.82 K and the saturation field at 0.57 K does not exceed 5.2 T. Application of the field along the $`[110]`$ direction has an intermediate effect.

Even more remarkably, an additional phase transition induced by the external field at about half of the saturation field is clearly present at low temperatures for $`H[111]`$, but is absent for $`H[112]`$. For $`H[110]`$, this phase transition is present, but it becomes less and less pronounced as the temperature decreases (see bottom panel of Fig. 3). The origin of this phase transition is not clear at present, although the fact that it happens at a half of $`H_{sat}`$ strongly suggests that it is likely to be related to a collinear spin state, where three spins on each tetrahedra are aligned with the field and the fourth is pointing in the opposite direction. A stabilisation of the collinear phase by quantum and thermal fluctuations has been predicted to take place at $`H=H_{sat}/2`$ for Heisenberg antiferromagnets on pyrochlore and frustrated square lattices, as well as at $`H=H_{sat}/3`$ for kagome and garnet lattices Field . Whether or not GTO, whose magnetic interactions are more complex than a simple Heisenberg model, does indeed adopt a collinear order in a magnetic field remains to be established with more direct probes of the magnetic structure, such as neutron diffraction measurements.

When compared to the phase diagram obtained from the $`C(H,T)`$ measurements on a powder sample Ramirez , Fig. 4 demonstrates the advantages of using single crystals. While it has been claimed that there are as many as four different ordered phases in GTO Ramirez , our measurements show that for any given direction of magnetic field there are at most three ordered phases. The phase diagram reported by Ramirez et al. Ramirez  seems to consist of a superposition of the single crystal phase diagrams.

When compared to the theoretical predictions based on a simple mean-field model Ramirez , our phase diagrams reveal significant discrepancies. For example, for the field along $`[111]`$, theory predicts only one transition, while our measurements clearly show two. These observations strongly suggest the need for further development of the theoretical model, which most likely should include some degree of magnetic anisotropy.

Recently Hassan et al. ESR  have reported the presence of a surprisingly large anisotropy in GTO below a temperature of 80 K, which they observed by measuring the ESR spectra. The anisotropy reveals itself as a well pronounced shift in the position of one of the resonance lines from its maximum value for $`H[111]`$ to its minimum for $`H[111]`$. The shift amounts to 3.6 T when measured at a frequency of 54 GHz. No change in the position of the resonance lines was observed when the field was rotated in a plane perpendicular to the $`[111]`$ direction ESR . The same authors have also reported dc-torque measurements dc-torque , which have shown a gradual development of macroscopic anisotropy below 80 K.

Motivated by these results, we have measured the temperature dependence of the magnetic susceptibility of a GTO crystal for $`H[111]`$ and $`H[111]`$ using a SQUID magnetometer. No significant anisotropy in the susceptibility has been detected SQUID . At present, neither the origin of the anisotropy observed in the ESR measurements ESR , nor its presence at a relatively high temperature have a satisfactory explanation. Any explanation should obviously take into consideration more than a simple Heisenberg exchange interaction between the nearest neighbor spins.