text:
# T violation in Λ_𝑏→Λℓ⁺ℓ⁻ decays

## Abstract

We examine the T-odd transverse lepton and $`\mathrm{\Lambda }`$ polarizations in the baryonic decays of $`\mathrm{\Lambda }_b\mathrm{\Lambda }l^+l^{}(l=e,\mu )`$ to study the T violating effects. We show that the lepton polarizations are suppressed but the $`\mathrm{\Lambda }`$ ones can be as large as $`50\%`$ in CP violating theories beyond the standard model such as the SUSY models, which can be tested in various future hadron colliders.

Recently, time-reversal violation (TV) has been measured experimentally in the $`K^0`$ system , and thus complements the information on CP violation (CPV) that has been steadily accumulating for the past thirty seven years. The data is in accordance with the CPT theorem which is fundamental to local quantum field theories with Lorentz invariance and the usual spin-statistics connection. However, the origin of the violation remains unclear. In the standard model, CPV or TV arises from a unique physical phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix . This paradigm also predicts CPV effect in the b-quark system. To test the accuracy of this paradigm and to search for other sources of CPV one needs to look for new processes. Indeed this is an important quest of the B-factories. In addition we deem it particularly interesting if the time reversal symmetry violation can be directly detected in the b-system rather than inferring it as a consequence of CPT invariance. In this talk, we will present our recent studies on T violation in the baryonic decays of $`\mathrm{\Lambda }_b\mathrm{\Lambda }l^+l^{}(l=e,\mu )`$ .

It is known that for a general three-body decay of a baryon the triple spin-momentum correlations, such as $`\stackrel{}{s}_k\left(\stackrel{}{p}_i\times \stackrel{}{p}_j\right)`$, are T-odd observables, where $`\stackrel{}{s}_k`$ and $`\stackrel{}{p}_{i,j}`$ are the spin and momentum vectors of the final particles, respectively. There are a number of different sources that might give rise to these T-odd observables. The most important ones being the weak CPV such as the CKM phase of the SM. However, final state interactions such as QCD for non-leptonic decays or the electromagnetic (EM) interaction among the final state particles can also make contributions. These are usually less interesting and they could even hide the signals from the weak CPV. We note that the T-odd triple correlations do not need non-zero strong phases unlike some of CP violating observables, such as the rate asymmetry between a particle and its antiparticle. In $`\mathrm{\Lambda }_b\mathrm{\Lambda }l^+l^{}`$ decays, we can use either the polarization of the lepton ($`s_l`$) or the $`\mathrm{\Lambda }`$ baryon ($`s_\mathrm{\Lambda }`$) to study the T-odd correlations. It is interesting to note that the decays can be used to analyze the helicity structure of the interaction, which is impossible in the B meson system.

We start with the effective Hamiltonian for $`bsl^+l^{}`$ by including the right-handed couplings in the hadronic sector are given by

$`(bsl^+l^{})={\displaystyle \frac{G_F\alpha _{em}}{\sqrt{2}\pi }}V_{tb}V_{ts}^{}(H_{1\mu }L_V^\mu +H_{2\mu }L_A^\mu )`$
$`H_{1\mu }`$ $`=`$ $`\overline{s}\gamma ^\mu (C_9^LP_L+C_9^RP_R)b`$
$`{\displaystyle \frac{2m_b}{q^2}}\overline{s}i\sigma _{\mu \nu }q^\nu (C_7^LP_R+C_7^RP_L)b,`$
$`H_{2\mu }`$ $`=`$ $`\overline{s}\gamma ^\mu (C_{10}^LP_L+C_{10}^RP_R)b,`$
$`L_V^\mu `$ $`=`$ $`\overline{l}\gamma ^\mu l,L_A^\mu =\overline{l}\gamma ^\mu \gamma _5l,`$ (1)
where $`C_i^L`$ and $`C_i^R\left(i=7,9,10\right)`$ denote the effective Wilson coefficients of left- and right-handed couplings, respectively. In the standard model,

$`C_9^L=C_9^{eff},`$ $`C_9^R=0,`$
$`C_{10}^L=C_{10},`$ $`C_{10}^R=0,`$
$`C_7^L=C_7^{eff},`$ $`C_7^R={\displaystyle \frac{m_s}{m_b}}C_7^{eff},`$ (2)
where $`C_9^{eff}`$, $`C_{10}`$, and $`C_7^{eff}`$ are the standard Wilson coefficients .

To study the exclusive decays of $`\mathrm{\Lambda }_b\mathrm{\Lambda }l^+l^{}`$, one needs to know the form factors in the transition of $`\mathrm{\Lambda }_b(p_{\mathrm{\Lambda }_b})\mathrm{\Lambda }(p_\mathrm{\Lambda })`$, parametrized generally as follows:

$`\mathrm{\Lambda }\left|\overline{s}\mathrm{\Gamma }_\mu b\right|\mathrm{\Lambda }_b=f_1^{(T)}\overline{u}_\mathrm{\Lambda }\gamma _\mu u_{\mathrm{\Lambda }_b}`$ (3)
$`+`$ $`f_2^{(T)}\overline{u}_\mathrm{\Lambda }i\sigma _{\mu \nu }q^\nu u_{\mathrm{\Lambda }_b}+f_3^{(T)}q_\mu \overline{u}_\mathrm{\Lambda }u_{\mathrm{\Lambda }_b},`$
$`\mathrm{\Lambda }\left|\overline{s}\mathrm{\Gamma }_\mu \gamma _5b\right|\mathrm{\Lambda }_b=g_1^{(T)}\overline{u}_\mathrm{\Lambda }\gamma _\mu \gamma _5u_{\mathrm{\Lambda }_b}`$
$`+`$ $`g_2^{(T)}\overline{u}_\mathrm{\Lambda }i\sigma _{\mu \nu }q^\nu \gamma _5u_{\mathrm{\Lambda }_b}+g_3^{(T)}q_\mu \overline{u}_\mathrm{\Lambda }\gamma _5u_{\mathrm{\Lambda }_b},`$
where $`\mathrm{\Gamma }_\mu =\gamma _\mu (i\sigma _{\mu \nu })`$, and $`f_i^{(T)}`$ and $`g_i^{(T)}`$ are the form factors of vector (tensor) and axial-vector (axial-tensor) currents, respectively. In the heavy quark effective theory (HQET) , the form factors in Eq. (3) can be simplified by using

$`\mathrm{\Lambda }(p_\mathrm{\Lambda })\left|\overline{s}\mathrm{\Gamma }b\right|\mathrm{\Lambda }_b(p_{\mathrm{\Lambda }_b})`$ (4)
$`=`$ $`\overline{u}_\mathrm{\Lambda }(F_1(q^2)+/vF_2(q^2))\mathrm{\Gamma }u_{\mathrm{\Lambda }_b},`$
where $`\mathrm{\Gamma }`$ denotes the Dirac matrix, $`v=p_{\mathrm{\Lambda }_b}/M_{\mathrm{\Lambda }_b}`$ is the four-velocity of $`\mathrm{\Lambda }_b`$, and $`q=p_{\mathrm{\Lambda }_b}p_\mathrm{\Lambda }`$ is the momentum transfer, and the relations among the form factors can be found in Ref. . Explicitly, under the HEQT, we have

$`f_1`$ $`=`$ $`g_2=f_2^T=g_2^T=F_1+\sqrt{r}F_2,`$
$`\rho `$ $``$ $`M_{\mathrm{\Lambda }_b}\left({\displaystyle \frac{f_2+g_2}{f_1+g_1}}\right)={\displaystyle \frac{M_{\mathrm{\Lambda }_b}}{q^2}}\left({\displaystyle \frac{f_1^T+g_1^T}{f_1+g_1}}\right)`$ (5)
$`={\displaystyle \frac{F_2}{F_1+\sqrt{r}F_2}}.`$
From experiments  and theories , one finds that $`|F_1/F_2|>1`$ and $`|\rho |<1`$.

To explore the T violating effects using the lepton or $`\mathrm{\Lambda }`$ spin polarization, we write $`i`$ ($`i=l`$ or $`\mathrm{\Lambda }`$) four-spin vector in terms of a unit vector, $`\widehat{\xi }^i`$, along the spin in its rest frame, as

$`s_{i0}={\displaystyle \frac{\stackrel{}{p}_i\widehat{\xi }_i}{M_i}},\stackrel{}{s}_i=\widehat{\xi }_i+{\displaystyle \frac{s_{i0}}{E_i+M_i}}\stackrel{}{p}_i,`$ (6)
and choose the unit vectors along the longitudinal, normal, transverse components of the $`i`$ polarization, to be

$`\widehat{e}_L^i`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{p}_i}{|\stackrel{}{p}_i|}},\widehat{e}_N^i={\displaystyle \frac{\stackrel{}{p}_i\times \left(\stackrel{}{p}_l^{}\times \stackrel{}{p}_\mathrm{\Lambda }\right)}{\left|\stackrel{}{p}_i\times \left(\stackrel{}{p}_l^{}\times \stackrel{}{p}_\mathrm{\Lambda }\right)\right|}},`$
$`\widehat{e}_T^i`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{p}_l^{}\times \stackrel{}{p}_\mathrm{\Lambda }}{\left|\stackrel{}{p}_l^{}\times \stackrel{}{p}_\mathrm{\Lambda }\right|}},`$ (7)
respectively. Hence, the differential decay rates with polarized $`i`$ is given by

$`d\mathrm{\Gamma }^i(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}d\mathrm{\Gamma }^0(t)\left[1+\stackrel{}{P}^i\widehat{\xi }^i\right],`$ (8)
where $`t=E_\mathrm{\Lambda }/M_{\mathrm{\Lambda }_b}`$, $`\stackrel{}{P}^i`$ is the $`i=l`$ or $`\mathrm{\Lambda }`$ polarization vector, defined by

$`\stackrel{}{P}^i=P_L^i\widehat{e}_L^i+P_N^i\widehat{e}_N^i+P_T^i\widehat{e}_T^i,`$ (9)
and $`d\mathrm{\Gamma }^0(t)`$ is given in Ref. . With the T odd transverse $`i`$ polarization defined by

$`P_T^i`$ $`=`$ $`{\displaystyle \frac{d\mathrm{\Gamma }^i\left(\widehat{\xi }_i\widehat{e}_T^i=1\right)d\mathrm{\Gamma }^i\left(\widehat{\xi }_i\widehat{e}_T^i=1\right)}{d\mathrm{\Gamma }^i\left(\widehat{\xi }_i\widehat{e}_T^i=1\right)+d\mathrm{\Gamma }^i\left(\widehat{\xi }_i\widehat{e}_T^i=1\right)}},`$
we obtain

$`P_T^l`$ $``$ $`{\displaystyle \frac{m_l}{M_{\mathrm{\Lambda }_b}}}Im\left(C_9^LC_{10}^L+C_9^RC_{10}^R\right),`$ (10)
$`{\displaystyle \frac{m_l}{M_{\mathrm{\Lambda }_b}}}Im\left(C_7^LC_{10}^L+C_7^RC_{10}^R\right)`$
$`P_T^\mathrm{\Lambda }`$ $``$ $`Im\left(C_9^RC_{10}^LC_9^LC_{10}^R\right),`$ (11)
$`Im\left(C_7^LC_{10}^RC_7^RC_{10}^L\right)`$
It is easy to see that from Eq. (10) $`P_T^l`$ is suppressed since it is always associated with the small lepton mass and to have a non-zero value of $`P_T^\mathrm{\Lambda }`$, from Eq. (11), it is necessary to have conditions of (a) the existence of $`C_k^R`$ and (b) a phase of $`C_i^LC_j^R(ij)`$. We remark that these conditions for $`P_T^\mathrm{\Lambda }`$ are clearly different from those of the T odd transverse lepton polarizations in both inclusive decay of $`bsl^+l^{}`$  and exclusive ones, such as $`BK^{()}l^+l^{}`$  and $`\mathrm{\Lambda }_b\mathrm{\Lambda }l^+l^{}`$ . Such property distinguishes the T odd transverse $`\mathrm{\Lambda }`$ polarizations from various other T odd observables.

In the standard model, since there are no $`C_9^R`$ and $`C_{10}^R`$ as seen from Eq. (2), we have

$`P_T^l`$ $``$ $`{\displaystyle \frac{m_l}{m_b}}Im\left(C_7^{eff}C_{10}\right)`$
$`P_T^\mathrm{\Lambda }`$ $``$ $`{\displaystyle \frac{m_s}{m_b}}Im\left(C_7^{eff}C_{10}\right)`$ (12)
which are all suppressed. We note that the contribution to $`P_T^i`$ from the EM final state interaction are $`<O(10^3)`$. Moreover, the long-distance (LD) effects in the one-loop matrix elements of $`O_{1,2}`$  and $`\mathrm{\Lambda }_b\mathrm{\Lambda }J/\mathrm{\Psi }`$ with $`J/\mathrm{\Psi }l^+l^{}`$ are absorbed to $`C_9^L`$. On the other hand, it is clear that a large value of $`P_T^\mathrm{\Lambda }`$, according to Eq. (11), can be obtained if a theory contains $`C_9^R`$ or $`C_{10}^R`$ or a large $`C_7^R`$, with a non-zero phase. Many theories beyond the standard model could give rise to $`C_9^R`$ or $`C_{10}^R`$; examples are the left-right symmetric and supersymmetric models.

To illustrate our result, we use SUSY models  both with and without R-parity as shown in Ref. . In Figure 1,

we give the differential branching ratios (BRs) of $`\mathrm{\Lambda }_b\mathrm{\Lambda }\mu ^+\mu ^{}`$ with respect to $`E_\mathrm{\Lambda }`$ with and without including resonant states of $`\mathrm{\Psi }`$ and $`\mathrm{\Psi }^{}`$, and we find that the integrated BRs for the latter are $`(2.10,\mathrm{\hspace{0.17em}1.66},\mathrm{\hspace{0.17em}4.41})\times 10^6`$ for the standard and SUSY with and without R-parity models, respectively. In Figures 2 and 3,

we show $`P_T^\mathrm{\Lambda }(\mathrm{\Lambda }_b\mathrm{\Lambda }\mu ^+\mu ^{})`$ as a function of $`E_\mathrm{\Lambda }/M_{\mathrm{\Lambda }_b}`$ for the two types of models, with and without the LD contributions, respectively. As seen from the figures, even though the derivations of BRs to the standard model result are insignificant, the transverse $`\mathrm{\Lambda }`$ polarization asymmetries can be over $`50\%`$ in both SUSY models with and without R-parity. Similar results are also expected for the decay of $`\mathrm{\Lambda }_b\mathrm{\Lambda }e^+e^{}`$. We remark that measuring a large $`P_T^\mathrm{\Lambda }`$ in $`\mathrm{\Lambda }_b\mathrm{\Lambda }l^+l^{}`$ is a clean indication of T violation as well as new CP violation mechanism beyond the standard model. To measure $`P_T^\mathrm{\Lambda }(\mathrm{\Lambda }_b\mathrm{\Lambda }\mu ^+\mu ^{})10\%`$ at $`3\sigma `$ level, at least $`4.5\times 10^7`$ $`\mathrm{\Lambda }_b`$ decays are required if we use $`BR(\mathrm{\Lambda }_b\mathrm{\Lambda }\mu ^+\mu ^{})2\times 10^6`$. Clearly, the measurement could be done in the second generation of B-physics experiments, such as LHCb, ATLAS, and CMS at the LHC, and BTeV at the Tevatron, which produce $`10^{12}b\overline{b}`$ pairs per year . This is certainly within reach of a super B factory under discussion now . Finally, we note that similar results are expected for using polarized $`\mathrm{\Lambda }_b`$ instead of $`\mathrm{\Lambda }`$.

Acknowledgments This work was supported in part by the National Center for Theoretical Science, National Science Council of the Republic of China under Contract Nos. NSC-90-2112-M-001-069 and NSC-90-2112-M-007-040, and National Science and Engineering Research Council of Canada.


format: text

text:
# References

hep-ph/0210197

RAL-TR-2002-023

30 Sep 2002

$`\mu \tau `$ Reflection Symmetry in Lepton Mixing and

Neutrino Oscillations

P. F. Harrison

Physics Department, Queen Mary University of London

Mile End Rd. London E1 4NS. UK 111E-mail:p.f.harrison@qmul.ac.uk

and

W. G. Scott

Rutherford Appleton Laboratory

Chilton, Didcot, Oxon OX11 0QX. UK 222E-mail:w.g.scott@rl.ac.uk

## Abstract

We examine the possibility that the lepton mixing matrix may exhibit a strong form of $`\mu \tau `$ universality, in which corresponding elements of the $`\mu `$\- and $`\tau `$-neutrino flavour eigenstates have equal moduli, and find that it is consistent with present data on neutrino oscillations. We point out that in the Dirac case, this is equivalent to symmetry under $`\mu \tau `$ reflection, ie. the combined operation of $`\mu \tau `$ flavour exchange in the MNS matrix with a $`CP`$ transformation on the whole leptonic sector. We give the most general form for such a mixing matrix, examine the lepton mass matrices under such a symmetry, and explore the observable manifestations of the symmetry in neutrino oscillations.

(To be published in Physics Letters B)

1 Strong $`\mu \tau `$ Universality in Lepton Mixing The data on atmospheric neutrinos   indicate that the muon neutrino couples to the heaviest neutrino mass eigenstate, $`\nu _3`$, with a strength $`|U_{\mu 3}|1/\sqrt{2}`$, where $`U`$ is the conventional MNS lepton mixing matrix . Independently, the data from the CHOOZ  and PALO VERDE  reactor neutrino experiments indicate that the electron neutrino has at most a small coupling to the heaviest neutrino, $`|U_{e3}|^2\stackrel{<}{}0.03`$. These two facts taken together with the normalisation condition on the $`\nu _3`$ eigenstate imply that $`|U_{\tau 3}|1/\sqrt{2}`$ and hence, in particular, that $`|U_{\mu 3}||U_{\tau 3}|`$. Since $`U_{e3}`$ is small and $`|U_{\mu 3}||U_{\tau 3}|`$, it follows immediately from the orthogonality of the $`\nu _3`$ eigenstate with the $`\nu _1`$ and $`\nu _2`$ eigenstates respectively that $`|U_{\mu i}||U_{\tau i}|`$ for $`i=1,2`$ also. Thus each of the moduli of the three $`\mu `$-flavour mixing elements is at least approximately equal to that of the corresponding $`\tau `$-flavour mixing element.

In this paper, we examine the possibility that

$$|U_{\mu i}|=|U_{\tau i}|,\mathrm{exactly},\mathrm{for}\mathrm{all}i=1,2,3$$
(1)
and find that it is consistent with all the present data on neutrino oscillations (excluding the LSND results , which are awaiting independent confirmation). This strong form of $`\mu \tau `$ universality (as distinct from the familiar lepton flavour universality implicit in the Standard Model, in which $`|U_{\alpha 1}|^2+|U_{\alpha 2}|^2+|U_{\alpha 3}|^2=1,\alpha `$) is interesting, as it removes two of the otherwise free parameters of the lepton mixing matrix. Several previously suggested hypotheses for the lepton mixing matrix satisfy Eq. (1), in particular, the trimaximal , bimaximal  and Altarelli-Ferruglio  models, (as well as the more recently proposed tri-bimaximal  and tri-$`\chi `$maximal  schemes). The generalisation of these hypotheses to a unitary mixing matrix simply satisfying Eq. (1) was already identified as interesting in  and such a mixing matrix was implicit in . In this paper, we examine this class of mixing matrices and the associated phenomenology in more detail.

2 Mixing Matrix Parameterisation and $`\mu \mathbf{}\tau `$ Reflection Symmetry We develop here a simple parameterisation of mixing matrices satisfying Eq. (1) in terms of only moduli of the mixing matrix elements, as distinct from the standard parameterisation in terms of mixing angles and a phase . A unitary matrix cannot contain two identical rows. Hence, in a unitary mixing matrix satisfying Eq. (1), the elements of the two rows corresponding to the $`𝝂_𝝁`$ and $`𝝂_𝝉`$ respectively must contain relative phases. For Dirac neutrinos, we are free to rephase independently the rows and columns of the mixing matrix with no observable consequence, and we can exploit this freedom (for columns) to arrange that $`𝑼_{𝝁𝒊}\mathbf{=}𝑼_{𝝉𝒊}^{\mathbf{}}`$ for each such pair of elements. We can then write:

$`𝑼\mathbf{=}\mathbf{\left(}\begin{array}{c}\underset{\mathbf{¯}}{𝒖}\\ \underset{\mathbf{¯}}{𝒗}\\ \underset{\mathbf{¯}}{𝒗}^{\mathbf{}}\end{array}\mathbf{\right)}\mathbf{=}\mathbf{\left(}\begin{array}{ccc}𝒖_\mathrm{𝟏}& 𝒖_\mathrm{𝟐}& 𝒖_\mathrm{𝟑}\\ 𝒗_\mathrm{𝟏}& 𝒗_\mathrm{𝟐}& 𝒗_\mathrm{𝟑}\\ 𝒗_\mathrm{𝟏}^{\mathbf{}}& 𝒗_\mathrm{𝟐}^{\mathbf{}}& 𝒗_\mathrm{𝟑}^{\mathbf{}}\end{array}\mathbf{\right)}`$ (2)
which is sufficient to ensure Eq. (1). The normalisation of any column of $`𝑼`$ leads to

$$\mathbf{|}𝒗_𝒊\mathbf{|}\mathbf{=}\mathbf{\left(}\frac{\mathrm{𝟏}\mathbf{}\mathbf{|}𝒖_𝒊\mathbf{|}^\mathrm{𝟐}}{\mathrm{𝟐}}\mathbf{\right)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}\mathbf{,}\mathbf{}𝒊$$
(3)
while the orthogonality of any pair of columns gives

$$𝒖_𝒊𝒖_𝒋^{\mathbf{}}\mathbf{=}\mathbf{}\mathrm{𝟐}\mathrm{𝐑𝐞}\mathbf{(}𝒗_𝒊𝒗_𝒋^{\mathbf{}}\mathbf{)}\mathbf{}𝒊\mathbf{}𝒋\mathbf{.}$$
(4)
Eq. (4) can be satisfied only if all three elements of $`\underset{\mathbf{¯}}{𝒖}`$ have a common phase and we use the freedom to rephase the electron neutrino flavour eigenstate to set $`\underset{\mathbf{¯}}{𝒖}`$ real and positive. Using Eq. (3), we may now re-write Eq. (4) as

$$\mathrm{𝐜𝐨𝐬}𝜶_{𝒊𝒋}\mathbf{=}\frac{\mathbf{}𝒖_𝒊𝒖_𝒋}{\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒊^\mathrm{𝟐}\mathbf{)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒋^\mathrm{𝟐}\mathbf{)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}}\mathbf{}𝒊\mathbf{}𝒋$$
(5)
where $`𝜶_{𝒊𝒋}`$ is the relative phase between $`𝒗_𝒊`$ and $`𝒗_𝒋`$111Unitarity triangles formed from the orthogonality condition of $`\underset{¯}{u}`$ with $`\underset{¯}{v}`$ or $`\underset{¯}{v}^{}`$ are always acute-angled triangles. Those formed from pairs of columns are isosceles triangles.. The row $`\underset{\mathbf{¯}}{𝒗}`$ (and hence also the row $`\underset{\mathbf{¯}}{𝒗}^{\mathbf{}}`$) is now completely determined up to an overall phase by the magnitudes of the elements of $`\underset{\mathbf{¯}}{𝒖}`$ (only two of which are independent, $`\underset{\mathbf{¯}}{𝒖}`$ being normalised to unity).

If the MNS matrix is a matrix of the form, Eq. (2), with $`\underset{\mathbf{¯}}{𝒖}`$ real, the standard action is invariant under the simultaneous interchange of the $`𝝂_𝝁`$ and $`𝝂_𝝉`$ flavour eigenstates ($`\underset{\mathbf{¯}}{𝒗}`$ and $`\underset{\mathbf{¯}}{𝒗}^{\mathbf{}}`$ respectively) and complex conjugation of the MNS matrix. However, we know  that complex conjugation of the MNS matrix is formally equivalent to making a $`𝑪𝑷`$ transformation on the whole leptonic sector of the Standard Model (modified to include Dirac neutrino masses and standard MNS mixings). We will refer to the combined operation of $`𝝁\mathbf{}𝝉`$ flavour exchange in the MNS matrix and $`𝑪𝑷`$ transformation on the leptonic sector as “$`𝝁\mathbf{}𝝉`$ reflection”. Symmetry under $`𝝁\mathbf{}𝝉`$ reflection implies an MNS matrix which may be written with the $`𝝂_𝒆`$ flavour eigenstate real and the $`𝝂_𝝁`$ and $`𝝂_𝝉`$ flavour eigenstates complex conjugates of each other, Eq. (2). This is sufficient to ensure the strong $`𝝁\mathbf{}𝝉`$ universality expressed by Eq. (1), the two definitions being equivalent (for Dirac neutrinos222For Majorana neutrinos, $`\mu \tau `$ reflection symmetry is more constraining on the observable parameters of the mixing matrix than the strong $`\mu \tau `$ universality of Eq. (1). In the following sections of this paper however, there is no distinction between the Majorana and Dirac cases, and the universality and the symmetry may be taken to be equivalent.). As we show at the end of this section, $`𝝁\mathbf{}𝝉`$ reflection symmetry is consistent with all confirmed neutrino oscillation measurements currently available.

For a mixing matrix which is $`𝝁\mathbf{}𝝉`$ reflection symmetric, as in Eq. (2), Jarlskog’s $`𝑪𝑷`$-violating invariant  is given by

$`𝓙`$ $`\mathbf{=}`$ $`\mathrm{𝐈𝐦}\mathbf{(}𝒖_𝒊𝒖_𝒋^{\mathbf{}}𝒗_𝒊^{\mathbf{}}𝒗_𝒋\mathbf{)}`$ (6)
$`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}}}𝒖_𝒊𝒖_𝒋\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒊^\mathrm{𝟐}\mathbf{)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒋^\mathrm{𝟐}\mathbf{)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}\mathrm{𝐬𝐢𝐧}𝜶_{𝒊𝒋}`$ (7)
$`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}}}𝒖_𝒊𝒖_𝒋\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒊^\mathrm{𝟐}\mathbf{}𝒖_𝒋^\mathrm{𝟐}\mathbf{)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}`$ (8)
$`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathrm{𝟐}}}𝒖_\mathrm{𝟏}𝒖_\mathrm{𝟐}𝒖_\mathrm{𝟑}\mathbf{,}`$ (9)
a strikingly simple result (we recall the fact that the $`𝒖_𝒊`$ have been taken to be real). We note that quite generally, for a given, fixed pair of magnitudes of elements of the electron flavour eigenstate, $`\underset{\mathbf{¯}}{𝒖}`$, requiring $`𝛍\mathbf{}𝛕`$ reflection symmetry, is equivalent to requiring maximal $`𝐂𝐏`$ violation in the leptonic sector (the proof is not given here but follows readily by differentiation of the complete expression for $`𝓙`$ keeping $`\mathbf{|}𝒖_\mathrm{𝟐}\mathbf{|}`$ and $`\mathbf{|}𝒖_\mathrm{𝟑}\mathbf{|}`$ fixed; it may also been seen in the standard parameterisation , where Eq. (1) corresponds to fixing $`𝜽_{\mathrm{𝟐𝟑}}\mathbf{=}𝝅\mathbf{/}\mathrm{𝟒}`$ and $`𝜹\mathbf{=}\mathbf{\pm }𝝅\mathbf{/}\mathrm{𝟐}`$).

In order to make mixing matrices of the form Eq. (1) $`𝝁\mathbf{}𝝉`$ reflection symmetric, we have defined the phase convention of the first row and all the columns, and maintaining these conventions restricts somewhat our freedom to rephase the remaining rows of the matrix. However, we are still free to make a phase transformation on the matrix $`𝑼\mathbf{}𝚽𝑼`$ where $`𝚽\mathbf{=}\mathrm{𝐞𝐱𝐩}\mathbf{\{}𝒊\mathrm{𝐝𝐢𝐚𝐠}\mathbf{(}\mathrm{𝟎}\mathbf{,}\mathit{\varphi }\mathbf{,}\mathbf{}\mathit{\varphi }\mathbf{)}\mathbf{\}}`$, without any observable consequence, and without disturbing the $`𝝁\mathbf{}𝝉`$ reflection symmetry. We may therefore choose $`\mathit{\varphi }`$ such that any chosen neutrino mass eigenstate (any column) is wholly real. With such a phase convention, the matrix $`𝑼`$ has only two free parameters, which can be taken to be any pair of the (real) elements of $`\underset{\mathbf{¯}}{𝒖}`$.

It is perhaps of the most practical convenience to choose $`\mathit{\varphi }`$ such that the $`𝝂_\mathrm{𝟑}`$ eigenstate is real (so that the three parameters presently constrained directly by experiment, $`𝑼_{𝒆\mathrm{𝟐}}`$, $`𝑼_{𝒆\mathrm{𝟑}}`$ and $`𝑼_{𝝁\mathrm{𝟑}}`$ are all real) and to take as the two free parameters, $`𝑼_{𝒆\mathrm{𝟐}}\mathbf{=}𝒖_\mathrm{𝟐}`$ and $`𝑼_{𝒆\mathrm{𝟑}}\mathbf{=}𝒖_\mathrm{𝟑}`$. We know from the HOMESTAKE , SUPER-K  and SNO  solar neutrino meaurements that $`\mathbf{|}𝑼_{𝒆\mathrm{𝟐}}\mathbf{|}\mathbf{}\mathrm{𝟏}\mathbf{/}\sqrt{\mathrm{𝟑}}`$ and from the CHOOZ  and PALO VERDE  reactor experiments that $`\mathbf{|}𝑼_{𝒆\mathrm{𝟑}}\mathbf{|}\stackrel{\mathbf{<}}{\mathbf{}}\mathbf{0.17}`$. Hence, the hypothesis of $`𝝁\mathbf{}𝝉`$ reflection symmetry accommodates all atmospheric, solar, and reactor neutrino data currently available (each of these results is confirmed by more than one experiment). The LSND result , cannot be accommodated simultaneously with all these confirmed data using a $`\mathrm{𝟑}\mathbf{\times }\mathrm{𝟑}`$ MNS lepton mixing matrix, and would require a fourth family of neutrinos . We have not considered the question of $`𝝁\mathbf{}𝝉`$ reflection symmetry in such an extended picture.

3 $`\mu \mathbf{}\tau `$ Reflection Symmetry in the Flavour Basis Any symmetry of the lepton mixing matrix should have its origin in a corresponding symmetry of the lepton mass matrices expressed in a weak basis, namely a basis in which the weak interaction is flavour diagonal and universal. Such a basis is not uniquely defined however, and we choose to consider the special case333Another special case, the “circulant basis”, was already considered in , and a third alternative, in which the neutrino mass matrix is diagonal, is examined in Appendix A. in which the charged lepton mass matrix is diagonal (the flavour basis), which is defined up to some unobservable phase transformation. We will furthermore consider the Hermitian squares, $`𝑴^\mathrm{𝟐}\mathbf{}𝑴𝑴^{\mathbf{}}`$, of the mass matrices, $`𝑴`$, in this basis, in order to avoid an arbitrariness in the determination of the mass matrices themselves due to the non-observation of right-handed charged currents.

In the case of $`𝝁\mathbf{}𝝉`$ reflection symmetry (ie. assuming $`𝑼`$ has the form Eq. (2), with $`\underset{\mathbf{¯}}{𝒖}`$ real), the neutrino mass matrix (squared) is given in the flavour basis by

$$𝑴_𝝂^\mathrm{𝟐}\mathbf{=}𝑼𝑫_𝝂^\mathrm{𝟐}𝑼^{\mathbf{}}$$
(10)
where $`𝑫_𝝂^\mathrm{𝟐}\mathbf{=}\mathrm{𝐝𝐢𝐚𝐠}\mathbf{(}𝒎_\mathrm{𝟏}^\mathrm{𝟐}\mathbf{,}𝒎_\mathrm{𝟐}^\mathrm{𝟐}\mathbf{,}𝒎_\mathrm{𝟑}^\mathrm{𝟐}\mathbf{)}`$, is the diagonal neutrino mass matrix squared. Then

$`𝑴_𝝂^\mathrm{𝟐}\mathbf{=}\mathbf{\left(}\begin{array}{ccc}𝒛& 𝒘& 𝒘^{\mathbf{}}\\ 𝒘^{\mathbf{}}& 𝒙& 𝒚\\ 𝒘& 𝒚^{\mathbf{}}& 𝒙\end{array}\mathbf{\right)}\mathbf{,}`$ (11)
a form  which again explicitly respects $`𝝁\mathbf{}𝝉`$ reflection symmetry (in this basis, the rows and columns of $`𝑴_𝝂^\mathrm{𝟐}`$ are labelled by the flavour indices $`\mathbf{(}𝒆\mathbf{,}𝝁\mathbf{,}𝝉\mathbf{)}`$):

$$𝑴_𝝂^\mathrm{𝟐}\mathbf{}\mathbf{(}𝑴_𝝂^\mathrm{𝟐}\mathbf{)}^{\mathbf{}}\mathbf{}𝑬_{𝝁𝝉}𝑴_𝝂^\mathrm{𝟐}𝑬_{𝝁𝝉}^{\mathbf{}}\mathbf{=}\mathbf{(}𝑴_𝝂^\mathrm{𝟐}\mathbf{)}^{\mathbf{}}\mathbf{,}$$
(12)
where $`𝑬_{𝝁𝝉}`$ is the $`𝝁\mathbf{}𝝉`$ flavour exchange operator:

$`𝑬_{𝝁𝝉}\mathbf{=}\mathbf{\left(}\begin{array}{ccc}\mathrm{𝟏}& \mathrm{𝟎}& \mathrm{𝟎}\\ \mathrm{𝟎}& \mathrm{𝟎}& \mathrm{𝟏}\\ \mathrm{𝟎}& \mathrm{𝟏}& \mathrm{𝟎}\end{array}\mathbf{\right)}\mathbf{.}`$ (13)
$`𝑴_𝝂^\mathrm{𝟐}`$ has 5 observable degrees of freedom: $`𝒛\mathbf{,}𝒙\mathbf{,}\mathbf{|}𝒘\mathbf{|}\mathbf{,}\mathbf{|}𝒚\mathbf{|}`$ and the phase of the combination $`𝒘^\mathrm{𝟐}𝒚`$, which together correspond to the three neutrino masses, and the two mixing degrees of freedom in the matrix $`𝑼`$. In general, $`𝑪𝑷`$ violation enters when the combination $`𝒘^\mathrm{𝟐}𝒚`$ has a non-trivial phase: $`\mathrm{𝐈𝐦}\mathbf{(}𝒘^\mathrm{𝟐}𝒚\mathbf{)}\mathbf{=}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}𝚫𝒎_{\mathrm{𝟑𝟐}}^\mathrm{𝟐}𝓙`$ (where $`𝚫𝒎_{𝒊𝒋}^\mathrm{𝟐}\mathbf{=}𝒎_𝒊^\mathrm{𝟐}\mathbf{}𝒎_𝒋^\mathrm{𝟐}`$). The elements of $`𝑴_𝝂^\mathrm{𝟐}`$ are given by

$`𝒛`$ $`\mathbf{=}`$ $`𝒎_\mathrm{𝟏}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}𝒖_\mathrm{𝟐}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}𝒖_\mathrm{𝟑}^\mathrm{𝟐}`$ (14)
$`𝒙`$ $`\mathbf{=}`$ $`𝒎_\mathrm{𝟏}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}\mathbf{|}𝒗_\mathrm{𝟐}\mathbf{|}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}\mathbf{|}𝒗_\mathrm{𝟑}\mathbf{|}^\mathrm{𝟐}`$ (15)
$`𝒘`$ $`\mathbf{=}`$ $`𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}𝒖_\mathrm{𝟐}𝒗_\mathrm{𝟐}^{\mathbf{}}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}𝒖_\mathrm{𝟑}𝒗_\mathrm{𝟑}^{\mathbf{}}`$ (16)
$`𝒚`$ $`\mathbf{=}`$ $`𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}𝒗_\mathrm{𝟐}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}𝒗_\mathrm{𝟑}^\mathrm{𝟐}\mathbf{.}`$ (17)
These equations may be combined with Eqs. (3) and (5) to find closed expressions in terms of any desired pair of elements of the vector $`\underset{\mathbf{¯}}{𝒖}`$. We have left the unphysical phase of $`\underset{\mathbf{¯}}{𝒗}`$ (and $`\underset{\mathbf{¯}}{𝒗}^{\mathbf{}}`$) undefined, hence the overall phases of $`𝒘`$ and $`𝒚`$ are similarly unphysical and undefined, and may be fixed by the choice of the phase $`\mathit{\varphi }`$ discussed in Section 2.

Tri-$`𝝌`$maximal mixing  is a special case of $`𝝁\mathbf{}𝝉`$ reflection symmetry defined by $`𝒖_\mathrm{𝟐}\mathbf{=}\mathrm{𝟏}\mathbf{/}\sqrt{\mathrm{𝟑}}`$, $`𝒖_\mathrm{𝟑}\mathbf{=}\sqrt{\frac{\mathrm{𝟐}}{\mathrm{𝟑}}}\mathrm{𝐬𝐢𝐧}𝝌`$, and (in a particular phase convention) corresponds to

$`𝒛`$ $`\mathbf{=}`$ $`𝒎_\mathrm{𝟏}^\mathrm{𝟐}\mathbf{+}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟑}\mathbf{)}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}\mathbf{+}\mathbf{(}\mathrm{𝟐}\mathbf{/}\mathrm{𝟑}\mathbf{)}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}𝝌`$ (18)
$`𝒙`$ $`\mathbf{=}`$ $`𝒎_\mathrm{𝟏}^\mathrm{𝟐}\mathbf{+}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟑}\mathbf{)}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}\mathbf{[}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathbf{}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟑}\mathbf{)}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}𝝌\mathbf{]}`$ (19)
$`𝒘`$ $`\mathbf{=}`$ $`\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟑}\mathbf{)}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}\mathbf{[}\mathbf{}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟑}\mathbf{)}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}𝝌\mathbf{}\mathrm{𝟑}𝒊𝓙\mathbf{]}`$ (20)
$`𝒚`$ $`\mathbf{=}`$ $`\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟑}\mathbf{)}𝚫𝒎_{\mathrm{𝟐𝟏}}^\mathrm{𝟐}\mathbf{+}𝚫𝒎_{\mathrm{𝟑𝟏}}^\mathrm{𝟐}\mathbf{[}\mathbf{}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathrm{𝐜𝐨𝐬}^\mathrm{𝟐}𝝌\mathbf{+}\mathbf{(}\mathrm{𝟏}\mathbf{/}\mathrm{𝟔}\mathbf{)}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}𝝌\mathbf{}\mathrm{𝟑}𝒊𝓙\mathbf{]}`$ (21)
where $`𝓙\mathbf{=}\mathrm{𝐬𝐢𝐧}\mathrm{𝟐}𝝌\mathbf{/}\mathbf{(}\mathrm{𝟔}\sqrt{\mathrm{𝟑}}\mathbf{)}`$. Tri-bimaximal mixing  is obtained by setting $`𝝌\mathbf{=}\mathrm{𝟎}`$.

4 $`\mu \mathbf{}\tau `$ Reflection Symmetry in Neutrino Oscillations The amplitude $`𝑨_{𝜶𝜷}`$ for a neutrino of flavour $`𝜶`$ to be detected as a neutrino of flavour $`𝜷`$ is given as a function of propagation distance $`𝑳`$ by the (matrix) equation:

$$𝑨\mathbf{=}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}𝒊𝑯𝑳\mathbf{)}$$
(22)
where $`𝑯`$ is the effective neutrino Hamiltonian in the flavour basis (in the rest of this paper, all calculations are assumed to be in this flavour basis). We may take:

$$𝑯\mathbf{=}𝑴_𝝂^\mathrm{𝟐}\mathbf{/}\mathrm{𝟐}𝑬$$
(23)
where $`𝑴_𝝂^\mathrm{𝟐}`$ is the Hermitian square of the neutrino mass matrix with $`𝝁\mathbf{}𝝉`$ reflection symmetry given in Eq. (11), and $`𝑬`$ is the neutrino energy. Obviously, since $`𝑴_𝝂^\mathrm{𝟐}`$ is $`𝝁\mathbf{}𝝉`$ reflection symmetric, then so is $`𝑯`$, and we note further that any power or polynomial (with real coefficients) of $`𝑯`$ shares this symmetry.

Eq. (22) may be re-written :

$$𝑨\mathbf{=}\underset{𝒊}{\mathbf{}}𝑿^𝒊\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}𝒊𝝀_𝒊𝑳\mathbf{)}$$
(24)
where the $`𝝀_𝒊`$ are the eigenvalues of $`𝑯`$ and the $`𝑿^𝒊`$ are hermitian (projection) operators given by second-order polynomials (with real coefficients) of $`𝑯`$:

$$𝑿^𝒊\mathbf{=}\frac{\mathbf{}_{𝒋\mathbf{}𝒊}\mathbf{(}𝑯\mathbf{}𝝀_𝒋\mathbf{)}}{\mathbf{}_{𝒋\mathbf{}𝒊}\mathbf{(}𝝀_𝒊\mathbf{}𝝀_𝒋\mathbf{)}}\mathbf{.}$$
(25)
Clearly, the $`𝑿^𝒊`$ share the $`𝝁\mathbf{}𝝉`$ reflection symmetry of Eq. (11). Comparison with the alternative approach in which the Hamiltonian is diagonalised before exponentiation allows us to identify the elements of the $`𝑿^𝒊`$ with the familiar combinations of the lepton mixing matrix elements :

$$𝑿_{𝜶𝜷}^𝒊\mathbf{=}𝑼_{𝜶𝒊}𝑼_{𝜷𝒊}^{\mathbf{}}$$
(26)
(no summation over $`𝒊`$ is implied).

The squared amplitude $`\mathbf{|}𝑨_{𝜶𝜷}\mathbf{|}^\mathrm{𝟐}`$ for $`𝜶\mathbf{}𝜷`$ ($`𝜶\mathbf{=}𝜷`$) gives the appearance (survival) probabilities as a function of neutrino energy and propagation distance and may be decomposed as a sum of two matrices:

$$𝓟_{𝜶𝜷}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}\mathbf{|}𝑨_{𝜶𝜷}\mathbf{|}^\mathrm{𝟐}\mathbf{=}𝓢_{𝜶𝜷}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{+}𝓐_{𝜶𝜷}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}$$
(27)
where the $`𝑪𝑷`$ (and $`𝑻`$) symmetric and anti-symmetric terms are respectively:

$`𝓢_{𝜶𝜷}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`𝜹_{𝜶𝜷}\mathbf{}\mathrm{𝟒}{\displaystyle \underset{𝒊\mathbf{<}𝒋}{\mathbf{}}}𝑲_{𝜶𝜷}^{𝒊𝒋}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}\mathbf{(}𝚫_{𝒊𝒋}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}`$ (28)
$`\mathrm{𝐚𝐧𝐝}𝓐_{𝜶𝜷}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`\mathrm{𝟖}𝑱_{𝜶𝜷}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝚫_{\mathrm{𝟏𝟐}}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝚫_{\mathrm{𝟐𝟑}}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝚫_{\mathrm{𝟑𝟏}}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}`$ (29)
(the dependence on neutrino energy is implicit in the factors $`𝚫_{𝒊𝒋}\mathbf{=}\mathbf{(}𝒎_𝒊^\mathrm{𝟐}\mathbf{}𝒎_𝒋^\mathrm{𝟐}\mathbf{)}\mathbf{/}\mathrm{𝟐}𝑬`$). The respective $`𝑪𝑷`$-even coefficients are:

$$𝑲_{𝜶𝜷}^{𝒊𝒋}\mathbf{=}\mathrm{𝐑𝐞}\mathbf{(}𝑼_{𝜶𝒊}𝑼_{𝜷𝒊}^{\mathbf{}}𝑼_{𝜶𝒋}^{\mathbf{}}𝑼_{𝜷𝒋}\mathbf{)}\mathbf{=}\mathrm{𝐑𝐞}\mathbf{(}𝑿_{𝜶𝜷}^𝒊𝑿_{𝜶𝜷}^𝒋\mathbf{}\mathbf{)}\mathbf{,}$$
(30)
and the $`𝑪𝑷`$-odd one is:

$$𝑱_{𝜶𝜷}\mathbf{=}\mathrm{𝐈𝐦}\mathbf{(}𝑼_{𝜶𝒊}𝑼_{𝜷𝒊}^{\mathbf{}}𝑼_{𝜶𝒋}^{\mathbf{}}𝑼_{𝜷𝒋}\mathbf{)}\mathbf{=}\mathrm{𝐈𝐦}\mathbf{(}𝑿_{𝜶𝜷}^𝒊𝑿_{𝜶𝜷}^𝒋\mathbf{}\mathbf{)}\mathbf{=}\mathit{ϵ}_{𝜶𝜷}𝓙$$
(31)
where the flavour-antisymmetric matrix $`\mathit{ϵ}`$ is

$`\mathit{ϵ}\mathbf{=}\mathbf{\left(}\begin{array}{ccc}\mathrm{𝟎}& \mathrm{𝟏}& \mathbf{}\mathrm{𝟏}\\ \mathbf{}\mathrm{𝟏}& \mathrm{𝟎}& \mathrm{𝟏}\\ \mathrm{𝟏}& \mathbf{}\mathrm{𝟏}& \mathrm{𝟎}\end{array}\mathbf{\right)}`$ (32)
and $`𝓙`$ is Jarlskog’s $`𝑪𝑷`$\- and $`𝑻`$-odd invariant. For anti-neutrinos, $`𝓟_{\overline{𝜶}\overline{𝜷}}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ is given by the same formulation modified by $`𝓙\mathbf{}\mathbf{}𝓙`$.

As already noted, the projection operators $`𝑿^𝒊`$ are symmetric under $`𝝁\mathbf{}𝝉`$ reflection:

$$𝑿^𝒊\mathbf{}\mathbf{(}𝑿^𝒊\mathbf{)}^{\mathbf{}}\mathbf{}𝑬_{𝝁𝝉}𝑿^𝒊𝑬_{𝝁𝝉}^{\mathbf{}}\mathbf{=}\mathbf{(}𝑿^𝒊\mathbf{)}^{\mathbf{}}$$
(33)
where $`𝑬_{𝝁𝝉}`$ was defined in Eq. (13). Hence, the real matrices $`𝑲^{𝒊𝒋}`$ and $`𝑱`$ derived from them are respectively symmetric (anti-symmetric) under $`𝝁\mathbf{}𝝉`$ flavour exchange (as well as under the $`𝑪𝑷`$-transformation), and are both therefore symmetric under full $`𝝁\mathbf{}𝝉`$ reflection 444The complex matrices (in flavour-space) of plaquettes  $`\mathrm{\Pi }^{ij}=K^{ij}+iJ(i<j)`$ are symmetric under $`\mu \tau `$ reflection, themselves having the form Eq. (11).. Hence, the $`𝑪𝑷`$-even and $`𝑪𝑷`$-odd components of $`𝓟\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ defined in Eq. (29) are separately $`𝝁\mathbf{}𝝉`$ reflection symmetric, $`𝓢\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ having the analogous form to that of Eq. (11) with all parameters real, while $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$, is proportional to $`\mathit{ϵ}𝓙`$ whose factors flip sign under $`𝝁\mathbf{}𝝉`$ flavour exchange and under $`𝑪𝑷`$ respectively.

By unitarity (applied separately for $`𝝂`$ and $`\overline{𝝂}`$), the rows and columns of $`𝓢\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ each sum to unity (at each value of $`𝑳\mathbf{/}𝑬`$), while those of $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ each sum to zero. By $`𝑪𝑷𝑻`$ symmetry, $`𝓢\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ is a symmetric matrix ($`𝓢\mathbf{=}𝓢^𝑻`$; for brevity, we drop the dependence on $`𝑳\mathbf{/}𝑬`$), while $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ is anti-symmetric $`\mathbf{(}𝓐\mathbf{=}\mathbf{}𝓐^𝑻\mathbf{)}`$. In the general vacuum case, these constraints reduce the number of independent functions of $`𝑳\mathbf{/}𝑬`$ which are needed to specify $`𝓟\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ to three $`𝑪𝑷`$-even (eg. $`𝓢_{𝜶𝜷}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{,}𝜶\mathbf{<}𝜷`$) and one $`𝑪𝑷`$-odd (eg. $`𝓐_{𝒆𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$). Hence survival probabilities can be written entirely in terms of the $`𝑪𝑷`$-conserving parts of appearance probabilities ($`𝓟_{𝜶𝜶}\mathbf{=}\mathrm{𝟏}\mathbf{}𝓢_{𝜶𝜷}\mathbf{}𝓢_{𝜶𝜸}\mathbf{,}𝜶\mathbf{}𝜷\mathbf{}𝜸`$) and vice versa ($`𝓢_{𝜶𝜷}\mathbf{=}\frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝓟_{𝜶𝜶}\mathbf{}𝓟_{𝜷𝜷}\mathbf{+}𝓟_{𝜸𝜸}\mathbf{)}`$).

Adding the constraints of $`𝝁\mathbf{}𝝉`$ reflection symmetry further reduces the number of independent $`𝑪𝑷`$-even functions of $`𝑳\mathbf{/}𝑬`$ to two and we may therefore write:

$`𝓟\mathbf{=}\mathbf{\left(}\begin{array}{ccc}\mathrm{𝟏}\mathbf{}\mathrm{𝟐}𝓢_{𝒆𝝁}& 𝓢_{𝒆𝝁}& 𝓢_{𝒆𝝁}\\ 𝓢_{𝒆𝝁}& \mathrm{𝟏}\mathbf{}𝓢_{𝒆𝝁}\mathbf{}𝓢_{𝝁𝝉}& 𝓢_{𝝁𝝉}\\ 𝓢_{𝒆𝝁}& 𝓢_{𝝁𝝉}& \mathrm{𝟏}\mathbf{}𝓢_{𝒆𝝁}\mathbf{}𝓢_{𝝁𝝉}\end{array}\mathbf{\right)}\mathbf{+}𝓐`$ (34)
where $`𝓐\mathbf{=}\mathit{ϵ}𝓐_{𝒆𝝁}`$ and $`\mathit{ϵ}`$ was given in Eq. (32), or equivalently, in terms of survival probabilities ($`𝓢_{𝜶𝜶}\mathbf{=}𝓟_{𝜶𝜶}`$):

$`𝓟\mathbf{=}\mathbf{\left(}\begin{array}{ccc}𝓢_{𝒆𝒆}& \frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝓢_{𝒆𝒆}\mathbf{)}& \frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝓢_{𝒆𝒆}\mathbf{)}\\ \frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝓢_{𝒆𝒆}\mathbf{)}& 𝓢_{𝝁𝝁}& \frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{+}𝓢_{𝒆𝒆}\mathbf{}\mathrm{𝟐}𝓢_{𝝁𝝁}\mathbf{)}\\ \frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝓢_{𝒆𝒆}\mathbf{)}& \frac{\mathrm{𝟏}}{\mathrm{𝟐}}\mathbf{(}\mathrm{𝟏}\mathbf{+}𝓢_{𝒆𝒆}\mathbf{}\mathrm{𝟐}𝓢_{𝝁𝝁}\mathbf{)}& 𝓢_{𝝁𝝁}\end{array}\mathbf{\right)}\mathbf{+}𝓐\mathbf{.}`$ (35)
Hence there are only three functions to evaluate in order to completely specify the matrix $`𝓟\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$.

For completeness, we provide here the expressions for the various components of $`𝓟`$, (noting from above that $`𝓢_{𝒆𝝁}`$ and $`𝓢_{𝝁𝝉}`$ are not independent of $`𝓢_{𝒆𝒆}`$ and $`𝓢_{𝝁𝝁}`$):

$`𝓢_{𝒆𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`\mathrm{𝟏}\mathbf{}\mathrm{𝟒}{\displaystyle \underset{𝒊\mathbf{<}𝒋}{\mathbf{}}}𝒖_𝒊^\mathrm{𝟐}𝒖_𝒋^\mathrm{𝟐}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}\mathbf{(}𝚫_{𝒊𝒋}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}`$ (36)
$`𝓢_{𝝁𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`\mathrm{𝟏}\mathbf{}{\displaystyle \underset{𝒊\mathbf{<}𝒋}{\mathbf{}}}\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒊^\mathrm{𝟐}\mathbf{)}\mathbf{(}\mathrm{𝟏}\mathbf{}𝒖_𝒋^\mathrm{𝟐}\mathbf{)}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}\mathbf{(}𝚫_{𝒊𝒋}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}`$ (37)
$`𝓢_{𝒆𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`\mathrm{𝟐}{\displaystyle \underset{𝒊\mathbf{<}𝒋}{\mathbf{}}}𝒖_𝒊^\mathrm{𝟐}𝒖_𝒋^\mathrm{𝟐}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}\mathbf{(}𝚫_{𝒊𝒋}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}`$ (38)
$`𝓢_{𝝁𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`{\displaystyle \underset{𝒊\mathbf{}𝒋\mathbf{}𝒌}{\mathbf{}}}\mathbf{(}𝒖_𝒌^\mathrm{𝟐}\mathbf{}𝒖_𝒊^\mathrm{𝟐}𝒖_𝒋^\mathrm{𝟐}\mathbf{)}\mathrm{𝐬𝐢𝐧}^\mathrm{𝟐}\mathbf{(}𝚫_{𝒊𝒋}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}`$ (39)
$`𝓐_{𝒆𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`\mathrm{𝟒}\mathbf{(}𝒖_\mathrm{𝟏}𝒖_\mathrm{𝟐}𝒖_\mathrm{𝟑}\mathbf{)}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝚫_{\mathrm{𝟏𝟐}}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝚫_{\mathrm{𝟐𝟑}}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathrm{𝐬𝐢𝐧}\mathbf{(}𝚫_{\mathrm{𝟑𝟏}}𝑳\mathbf{/}\mathrm{𝟐}\mathbf{)}\mathbf{.}`$ (40)
Of course, only two elements of $`\underset{\mathbf{¯}}{𝒖}`$ are independent, and the third may be simply determined from a given pair using the normalisation of $`\underset{\mathbf{¯}}{𝒖}`$.

5 Matter Effects Matter effects are of considerable phenomenological importance. In this section, we discuss how they affect the predictions of $`𝝁\mathbf{}𝝉`$ reflection symmetry under quite general circumstances. When (anti-)neutrinos propagate in matter, the Hamiltonian in the flavour basis is modified: $`𝑯\mathbf{}\stackrel{\mathbf{~}}{𝑯}`$, with

$$\stackrel{\mathbf{~}}{𝑯}\mathbf{=}𝑯\mathbf{\pm }\mathrm{𝐝𝐢𝐚𝐠}\mathbf{(}𝒂\mathbf{(}𝑳\mathbf{)}\mathbf{,}\mathrm{𝟎}\mathbf{,}\mathrm{𝟎}\mathbf{)}$$
(41)
where $`𝒂\mathbf{(}𝑳\mathbf{)}\mathbf{=}\sqrt{\mathrm{𝟐}}𝑮_𝑭𝑵_𝒆\mathbf{(}𝑳\mathbf{)}`$, $`𝑵_𝒆\mathbf{(}𝑳\mathbf{)}`$ is the (in general position-dependent) electron number density of the matter, and the “$`\mathbf{}`$” sign applies in the case of anti-neutrinos. This modifies both the eigenvalues and eigenstates of the Hamiltonian in a non-trivial way, in general as a function of $`𝑳\mathbf{/}𝑬`$. However, in the case of a $`𝝁\mathbf{}𝝉`$ reflection symmetric Hamiltonian, Eq. (11), its symmetry is unaffected, as only the parameter $`𝒛`$ is changed ($`𝒛\mathbf{}\stackrel{\mathbf{~}}{𝒛}\mathbf{=}𝒛\mathbf{\pm }\mathrm{𝟐}𝑬𝒂\mathbf{(}𝑳\mathbf{)}`$). Hence, matter effects respect $`𝛍\mathbf{}𝛕`$ reflection symmetry and any result which is a consequence of the symmetry, is unaltered when neutrinos pass through matter.

We first consider the most general case of arbitrary matter profile, and no $`𝝁\mathbf{}𝝉`$ reflection symmetry. It is still convenient to decompose the matrix $`𝓟\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ of survival and appearance probabilities into $`𝑪𝑷`$-even and $`𝑪𝑷`$-odd parts, $`𝓢\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ and $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ respectively. Without solving the Schroedinger Equation for the particular $`𝑳`$-dependent Hamiltonian of Eq. (41) the functional dependences of $`𝓢`$ and $`𝓐`$ on $`𝑳\mathbf{/}𝑬`$ are unknown. Moreover, in arbitrary matter, several of the familiar constraints of the vacuum case are lost, in particular because $`𝑪𝑷𝑻`$ symmetry does not in general hold. For example, the matrix $`𝓢\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ is no longer symmetric, and $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ is no longer anti-symmetric; the survival probabilities $`𝓟_{𝝁𝝁}`$ and $`𝓟_{𝝉𝝉}`$ acquire $`𝑪𝑷`$-violating contributions . However, we can still learn a considerable amount from general principles: appealing to unitarity, as in the previous section, we see that in general, there are now four independent $`𝑪𝑷`$-even functions of $`𝑳\mathbf{/}𝑬`$, eg. we may write:

$`𝓢\mathbf{=}\mathbf{\left(}\begin{array}{ccc}\mathrm{𝟏}\mathbf{}𝓢_{𝒆𝝁}\mathbf{}𝓢_{𝒆𝝉}& 𝓢_{𝒆𝝁}& 𝓢_{𝒆𝝉}\\ 𝓢_{𝝁𝒆}& \mathrm{𝟏}\mathbf{}𝓢_{𝝁𝒆}\mathbf{}𝓢_{𝝁𝝉}& 𝓢_{𝝁𝝉}\\ 𝓢_{𝒆𝝉}\mathbf{}\mathbf{(}𝓢_{𝝁𝒆}\mathbf{}𝓢_{𝒆𝝁}\mathbf{)}& 𝓢_{𝝁𝝉}\mathbf{+}\mathbf{(}𝓢_{𝝁𝒆}\mathbf{}𝓢_{𝒆𝝁}\mathbf{)}& \mathrm{𝟏}\mathbf{}𝓢_{𝒆𝝁}\mathbf{}𝓢_{𝝁𝝉}\end{array}\mathbf{\right)}\mathbf{.}`$ (42)
Similarly, there are three $`𝑪𝑷`$-odd ones (we use the fact that $`𝓐_{𝒆𝒆}\mathbf{=}\mathrm{𝟎}`$ ) eg.:

$`𝓐\mathbf{=}\mathbf{\left(}\begin{array}{ccc}\mathrm{𝟎}& 𝓐_{𝒆𝝁}& \mathbf{}𝓐_{𝒆𝝁}\\ 𝓐_{𝝁𝒆}& 𝓐_{𝝁𝝁}& \mathbf{}𝓐_{𝝁𝒆}\mathbf{}𝓐_{𝝁𝝁}\\ \mathbf{}𝓐_{𝝁𝒆}& \mathbf{}𝓐_{𝒆𝝁}\mathbf{}𝓐_{𝝁𝝁}& 𝓐_{𝒆𝝁}\mathbf{+}𝓐_{𝝁𝒆}\mathbf{+}𝓐_{𝝁𝝁}\end{array}\mathbf{\right)}\mathbf{.}`$ (43)
The $`𝝁\mathbf{}𝝉`$ reflection symmetry makes the following predictions555It is interesting to note that in the case of the element $`𝒫_{\mu \tau }=𝒮_{\mu \tau }+𝒜_{\mu \tau }`$, the $`\mu \tau `$ exchange operator, $`E_{\mu \tau }`$, has the same effect as the time-reversal operator, $`T`$, and that the resulting constraints on these matrix elements due to $`\mu \tau `$ reflection symmetry are identical to those of $`CPT`$ symmetry. This accidental mimicking of $`CPT`$ symmetry is, however, restricted to the $`𝒫_{\mu \tau }`$ and $`𝒫_{\tau \mu }`$ elements, and $`\mu \tau `$ reflection symmetry generally leads to constraints distinct from those of $`CPT`$ symmetry.:

$`𝓢_{𝒆𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}𝓢_{𝒆𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ ; $`𝓢_{𝝁𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}𝓢_{𝝉𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ (44)
$`𝓢_{𝝁𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}𝓢_{𝝉𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ ; $`𝓢_{𝝁𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}𝓢_{𝝉𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ (45)
$`𝓐_{𝒆𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}\mathbf{}𝓐_{𝒆𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ ; $`𝓐_{𝝁𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}\mathbf{}𝓐_{𝝉𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ (46)
$`𝓐_{𝝁𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}\mathbf{}𝓐_{𝝉𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ ; $`𝓐_{𝝁𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}\mathbf{}𝓐_{𝝉𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ (47)
(clearly, some though not all of them are redundant relative to the constraints of unitarity already invoked, eg. Eq. (46)). Although following simply from $`𝝁\mathbf{}𝝉`$ reflection symmetry, we have verified these predictions in detail using the elegant formalism  for arbitrary matter profiles.

It is easy to see that the constraints from $`𝝁\mathbf{}𝝉`$ reflection symmetry given in Eqs. (44) and (45) reduce $`𝓢\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ to the form already given in Eq. (34), so that now it has only two independent components, as before. Similarly, constraints (46) and (47) simplify the form of $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ so that its complete specification now requires two functions to be determined, eg. $`𝓐_{𝝁𝝉}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ and $`𝓐_{𝝁𝝁}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ (or whichever are most convenient). Now, (using unitarity) we may write:

$`𝓐\mathbf{=}\mathbf{\left(}\begin{array}{ccc}\mathrm{𝟎}& 𝓐_{𝝁𝝉}\mathbf{}𝓐_{𝝁𝝁}& 𝓐_{𝝁𝝁}\mathbf{}𝓐_{𝝁𝝉}\\ \mathbf{}𝓐_{𝝁𝝉}\mathbf{}𝓐_{𝝁𝝁}& 𝓐_{𝝁𝝁}& 𝓐_{𝝁𝝉}\\ 𝓐_{𝝁𝝉}\mathbf{+}𝓐_{𝝁𝝁}& \mathbf{}𝓐_{𝝁𝝉}& \mathbf{}𝓐_{𝝁𝝁}\end{array}\mathbf{\right)}\mathbf{.}`$ (48)
For $`𝑻`$-symmetric matter profiles, $`𝓐_{𝝁𝝁}\mathbf{=}\mathrm{𝟎}`$ , and the symmetry of $`𝓐\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ reduces to that of the vacuum case, Eq. (32).

6 Application to Atmospheric Neutrinos The observation from atmospheric neutrinos  that $`\mathbf{|}𝑼_{𝝁\mathrm{𝟑}}\mathbf{|}\mathbf{}\mathrm{𝟏}\mathbf{/}\sqrt{\mathrm{𝟐}}`$, was important to the argument in Section 1 that Nature respects strong $`𝝁\mathbf{}𝝉`$ universality, at least approximately. However, $`\mathbf{|}𝑼_{𝝁\mathrm{𝟑}}\mathbf{|}\mathbf{=}\mathrm{𝟏}\mathbf{/}\sqrt{\mathrm{𝟐}}`$ is not a prediction of $`𝝁\mathbf{}𝝉`$ reflection symmetry in general. On the contrary, from Eq. (3), if the symmetry is good, we have that $`\mathbf{|}𝑼_{𝝁\mathrm{𝟑}}\mathbf{|}\mathbf{=}\mathbf{\left(}\frac{\mathrm{𝟏}\mathbf{}\mathbf{|}𝑼_{𝒆\mathrm{𝟑}}\mathbf{|}^\mathrm{𝟐}}{\mathrm{𝟐}}\mathbf{\right)}^{\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}}`$, which gives $`\mathbf{|}𝑼_{𝝁\mathrm{𝟑}}\mathbf{|}\mathbf{=}\mathrm{𝟏}\mathbf{/}\sqrt{\mathrm{𝟐}}`$ only if $`𝑼_{𝒆\mathrm{𝟑}}\mathbf{=}\mathrm{𝟎}`$, a possibility which is not excluded, but which is by no means essential. In fact, if $`\mathrm{𝟎}\mathbf{}\mathbf{|}𝑼_{𝒆\mathrm{𝟑}}\mathbf{|}\stackrel{\mathbf{<}}{\mathbf{}}\mathbf{0.17}`$ as is currently thought to be the case , then $`𝝁\mathbf{}𝝉`$ reflection symmetry predicts

$$\mathbf{0.696}\stackrel{\mathbf{<}}{\mathbf{}}\mathbf{|}𝑼_{𝝁\mathrm{𝟑}}\mathbf{|}\mathbf{=}\mathbf{|}𝑼_{𝝉\mathrm{𝟑}}\mathbf{|}\mathbf{}\mathrm{𝟏}\mathbf{/}\sqrt{\mathrm{𝟐}}\mathbf{,}$$
(49)
a rather precise prediction, and certainly in agreement with the data.

One piece of experimental data which we have not used in the present paper, is the non-observation of electron neutrino disappearance (or appearance) in the atmospheric neutrino oscillation data . It is tempting to think that the closeness of the measured atmospheric electron neutrino rate to the expected value, independent of zenith angle, is evidence for the fact that $`\mathbf{|}𝑼_{𝒆\mathrm{𝟑}}\mathbf{|}\mathbf{<<}\mathrm{𝟏}`$. We will show here that under $`𝝁\mathbf{}𝝉`$ reflection symmetry, and in the (realistic) approximation of an initial flux ratio of muon neutrinos to electron neutrinos, $`\mathit{\varphi }_{𝝂_𝝁}\mathbf{/}\mathit{\varphi }_{𝝂_𝒆}\mathbf{=}\mathrm{𝟐}`$, the measured electron neutrino survival probability is close to unity, independent of the value of $`\mathbf{|}𝐔_{𝐞\mathrm{𝟑}}\mathbf{|}`$, and hence that there is essentially no sensitivity to $`\mathbf{|}𝑼_{𝒆\mathrm{𝟑}}\mathbf{|}`$ in the rate of atmospheric electron neutrino events.

From Eq. (34), it can be seen that a prediction of $`𝝁\mathbf{}𝝉`$ reflection symmetry is that the $`𝑪𝑷`$-conserving part of the appearance probability for electron neutrinos in a muon neutrino beam, $`𝓢_{𝝁𝒆}`$, is exactly half of the disappearance probability from an electron neutrino beam, $`\mathrm{𝟏}\mathbf{}𝓟_{𝒆𝒆}`$, for all $`𝑳\mathbf{/}𝑬`$ (it is even true for arbitrary matter profiles, Eq. (44)). In the limit $`\mathit{\varphi }_{𝝂_𝝁}\mathbf{/}\mathit{\varphi }_{𝝂_𝒆}\mathbf{=}\mathrm{𝟐}`$, these two terms compensate each other exactly in the measured rate of atmospheric electron neutrinos. Defining the cross-section-weighted neutrino flavour fractions $`𝒓\mathbf{=}𝝈_{𝝂_𝒆}\mathit{\varphi }_{𝝂_𝒆}\mathbf{/}\mathbf{(}𝝈_{𝝂_𝒆}\mathit{\varphi }_{𝝂_𝒆}\mathbf{+}𝝈_{\overline{𝝂}_𝒆}\mathit{\varphi }_{\overline{𝝂}_𝒆}\mathbf{)}`$, $`𝒔\mathbf{=}𝝈_{𝝂_𝒆}\mathit{\varphi }_{𝝂_𝝁}\mathbf{/}\mathbf{(}𝝈_{𝝂_𝒆}\mathit{\varphi }_{𝝂_𝒆}\mathbf{+}𝝈_{\overline{𝝂}_𝒆}\mathit{\varphi }_{\overline{𝝂}_𝒆}\mathbf{)}`$ and $`\overline{𝒔}\mathbf{=}𝝈_{\overline{𝝂}_𝒆}\mathit{\varphi }_{\overline{𝝂}_𝝁}\mathbf{/}\mathbf{(}𝝈_{𝝂_𝒆}\mathit{\varphi }_{𝝂_𝒆}\mathbf{+}𝝈_{\overline{𝝂}_𝒆}\mathit{\varphi }_{\overline{𝝂}_𝒆}\mathbf{)}`$, the ratio of the observed electron neutrino rate to that expected with no oscillations is given by:

$`𝑹_𝒆\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}`$ $`\mathbf{=}`$ $`𝒓\mathbf{[}𝓟_{𝒆𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{+}\mathbf{(}\mathit{\varphi }_{𝝂_𝝁}\mathbf{/}\mathit{\varphi }_{𝝂_𝒆}\mathbf{)}𝓟_{𝝁𝒆}\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{]}\mathbf{+}\mathbf{(}\mathrm{𝟏}\mathbf{}𝒓\mathbf{)}\mathbf{[}𝝂\mathbf{}\overline{𝝂}\mathbf{]}`$ (50)
$`\mathbf{=}`$ $`𝒓\mathbf{+}\mathbf{(}\mathrm{𝟏}\mathbf{}𝒓\mathbf{)}\mathbf{+}𝒔𝓐_{𝝁𝒆}\mathbf{+}\overline{𝒔}𝓐_{\overline{𝝁}\overline{𝒆}}\mathrm{𝐟𝐨𝐫}\mathit{\varphi }_{𝝂_𝝁}\mathbf{/}\mathit{\varphi }_{𝝂_𝒆}\mathbf{=}\mathit{\varphi }_{\overline{𝝂}_𝝁}\mathbf{/}\mathit{\varphi }_{\overline{𝝂}_𝒆}\mathbf{=}\mathrm{𝟐}`$ (51)
$`\mathbf{}`$ $`\mathrm{𝟏}`$ (52)
where the $`𝑪𝑷`$-violating effects are suppressed by the smallness of $`𝑼_{𝒆\mathrm{𝟑}}`$ and $`𝚫𝒎_{\mathrm{𝟏𝟐}}^\mathrm{𝟐}`$, which would explain why present data are consistent with unity. While such compensation theorems have been noted previously by several authors for specific models , we emphasise here the perfect compensation of the $`𝑪𝑷`$-even terms for the generality of models with $`𝝁\mathbf{}𝝉`$ reflection symmetry, and for arbitrary matter profiles. We note furthermore, that although naively there is a partial cancellation in Eq. (52) between the terms involving $`𝓐_{𝝁𝒆}`$ and $`𝓐_{\overline{𝝁}\overline{𝒆}}`$ (owing to a relative minus sign in their dependences on $`𝓙`$), matter effects have different influences on their respective $`𝑳\mathbf{/}𝑬`$ dependences, and the cancellation may be incomplete, even if $`𝒔\mathbf{=}\overline{𝒔}`$.

Any observed deviation from $`𝑹_𝒆\mathbf{(}𝑳\mathbf{/}𝑬\mathbf{)}\mathbf{=}\mathrm{𝟏}`$ is therefore evidence for one of three things: i) a deviation from $`\mathbf{(}\mathit{\varphi }_{𝝂_𝝁}\mathbf{/}\mathit{\varphi }_{𝝂_𝒆}\mathbf{)}\mathbf{=}\mathrm{𝟐}`$ and/or from $`\mathbf{(}\mathit{\varphi }_{\overline{𝝂}_𝝁}\mathbf{/}\mathit{\varphi }_{\overline{𝝂}_𝒆}\mathbf{)}\mathbf{=}\mathrm{𝟐}`$, either of which would appear in the leading oscillations, given sufficient $`𝑳\mathbf{/}𝑬`$ resolution; ii) a deviation from $`𝝁\mathbf{}𝝉`$ reflection symmetry, which would also appear in the leading oscillation; iii) $`𝑪𝑷`$ violation, which is potentially distinguishable from i) and ii) as it would be manifest only on the longer distance scale of $`𝚫𝒎_{\mathrm{𝟏𝟐}}^\mathrm{𝟐}𝑳\mathbf{/}\mathrm{𝟒}𝑬\mathbf{\ll ̸}\mathrm{𝟏}`$.

7 Conclusions Current data from atmospheric, solar and reactor neutrino oscillation experiments are all consistent with the hypothesis of $`𝝁\mathbf{}𝝉`$ reflection symmetry for the leptonic sector of the Standard Model (augmented with neutrino masses and standard MNS mixings). The $`𝝁\mathbf{}𝝉`$ reflection transformation is the product of the $`𝝁\mathbf{}𝝉`$ exchange operator, $`𝑬_{𝝁𝝉}`$, applied to the neutrino flavour eigenstates, and the $`𝑪𝑷`$ operator. $`𝑪𝑷𝑬_{𝝁𝝉}`$ is clearly a unitary transformation and is also Hermitian with $`\mathbf{(}𝑪𝑷𝑬_{𝝁𝝉}\mathbf{)}^\mathrm{𝟐}\mathbf{=}𝑰`$, so that its eigenvalues are $`\mathbf{\pm }\mathrm{𝟏}`$. If $`𝝁\mathbf{}𝝉`$ reflection symmetry is a good symmetry of the neutrino oscillation Hamiltonian, one would expect the $`𝑪𝑷𝑬_{𝝁𝝉}`$ quantum number to be conserved in such processes. One might call this quantum number “$`𝝁\mathbf{}𝝉`$ reflectivity” or simply “mutativity”. The eigenstates of mutativity are clearly:

$`\mathrm{𝐦𝐮𝐭𝐚𝐭𝐢𝐯𝐢𝐭𝐲}\mathbf{+}\mathrm{𝟏}`$ $`\mathbf{:}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\sqrt{\mathrm{𝟐}}}}\mathbf{(}𝝂_𝝁\mathbf{+}\overline{𝝂}_𝝉\mathbf{)}\mathrm{𝐚𝐧𝐝}{\displaystyle \frac{\mathrm{𝟏}}{\sqrt{\mathrm{𝟐}}}}\mathbf{(}\overline{𝝂}_𝝁\mathbf{+}𝝂_𝝉\mathbf{)}`$ (53)
$`\mathrm{𝐦𝐮𝐭𝐚𝐭𝐢𝐯𝐢𝐭𝐲}\mathbf{}\mathrm{𝟏}`$ $`\mathbf{:}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\sqrt{\mathrm{𝟐}}}}\mathbf{(}𝝂_𝝁\mathbf{}\overline{𝝂}_𝝉\mathbf{)}\mathrm{𝐚𝐧𝐝}{\displaystyle \frac{\mathrm{𝟏}}{\sqrt{\mathrm{𝟐}}}}\mathbf{(}\overline{𝝂}_𝝁\mathbf{}𝝂_𝝉\mathbf{)}`$ (54)
although it is not easy to see how one might prepare a beam of neutrinos in such an eigenstate to test directly the conservation law in a non-trivial way.

As we have explored in the main body of this paper, in addition to the above postulated conservation law, the hypothesis of $`𝝁\mathbf{}𝝉`$ reflection symmetry makes a set of specific predictions about the $`𝑪𝑷`$-even and $`𝑪𝑷`$-odd parts of neutrino appearance and survival probabilities. Further detailed measurements are required to falsify this hypothesis, by finding deviations from the predictions, Eqs. (44)-(47), valid even in the presence of terrestrial matter, and Eqs. (1), (3), (5), (9), (49) and (52).

Acknowledgements It is a pleasure to thank Tom Weiler for useful suggestions. This work was supported by the UK Particle Physics and Astronomy Research Council (PPARC).

Appendix A: $`\mu \mathbf{}\tau `$ Reflection Symmetry in the Neutrino Mass Basis It is also interesting to consider the manifestations of $`𝝁\mathbf{}𝝉`$ reflection symmetry in another special weak basis, namely that in which the neutrino mass matrix is diagonal (the neutrino mass basis). In this basis, the lepton mixing degrees of freedom are all contained in the charged lepton mass matrix-squared (the weak interaction, and the neutrino mass matrix are both diagonal by definition in this basis), which is then given by

$$𝑴_{\mathbf{}}^\mathrm{𝟐}\mathbf{=}𝑼^{\mathbf{}}𝑫_{\mathbf{}}^\mathrm{𝟐}𝑼$$
(55)
where $`𝑫_{\mathbf{}}^\mathrm{𝟐}\mathbf{=}\mathrm{𝐝𝐢𝐚𝐠}\mathbf{(}𝒎_𝒆^\mathrm{𝟐}\mathbf{,}𝒎_𝝁^\mathrm{𝟐}\mathbf{,}𝒎_𝝉^\mathrm{𝟐}\mathbf{)}`$, from which we obtain

$$\mathbf{(}𝑴_{\mathbf{}}^\mathrm{𝟐}\mathbf{)}_{𝒊𝒋}\mathbf{=}\mathbf{\left(}\frac{𝒎_𝝁^\mathrm{𝟐}\mathbf{+}𝒎_𝝉^\mathrm{𝟐}}{\mathrm{𝟐}}\mathbf{\right)}𝜹_{𝒊𝒋}\mathbf{+}\mathbf{\left[}𝒎_𝒆^\mathrm{𝟐}\mathbf{}\mathbf{\left(}\frac{𝒎_𝝁^\mathrm{𝟐}\mathbf{+}𝒎_𝝉^\mathrm{𝟐}}{\mathrm{𝟐}}\mathbf{\right)}\mathbf{\right]}𝒖_𝒊𝒖_𝒋\mathbf{+}𝒊\mathit{ϵ}_{𝒊𝒋𝒌}\mathbf{\left(}\frac{𝒎_𝝉^\mathrm{𝟐}\mathbf{}𝒎_𝝁^\mathrm{𝟐}}{\mathrm{𝟐}}\mathbf{\right)}𝒖_𝒌\mathbf{,}$$
(56)
whose real part is symmetric under the interchange $`𝝁\mathbf{}𝝉`$ and whose imaginary part is anti-symmetric (in this basis, the rows and columns of $`𝑴_{\mathbf{}}^\mathrm{𝟐}`$ are labelled by the neutrino mass eigenstate indices, the flavour indices labelling only the charged lepton mass eigenvalues). This matrix explicitly respects $`𝝁\mathbf{}𝝉`$ reflection symmetry, as the effects of $`𝝁\mathbf{}𝝉`$ flavour exchange are clearly reversed by a $`𝑪𝑷`$ transformation, which is again equivalent to complex conjugation. The fact that the symmetry is manifest also in this basis demonstrates that the symmetry is not simply an accident in the flavour basis considered in the main text. $`𝑴_{\mathbf{}}^\mathrm{𝟐}`$ encodes 5 degrees of freedom - the three charged lepton masses, and any pair of elements of the $`𝝂_𝒆`$ flavour eigenstate, $`\underset{\mathbf{¯}}{𝒖}`$.


text:
# Electron energy loss and induced photon emission in photonic crystals

## I Introduction

Near-field spectroscopy can be performed using spatially-resolved electron energy loss spectroscopy in scanning transmission electron microscopes, where the electron beam provides an external evanescent field with which to probe the sample within a spatial range on the nanometer scale. In particular, the relatively intense, low-energy part ($`<50`$ eV) of the loss spectrum can be used with minimum sample damage to provide some insight on plasmons and other collective excitations.  For the electron velocities typically employed in microscopes (above half the speed of light) and for samples that are homogeneous across distances of a few nanometers, frequency-dependent dielectric functions are sufficient to describe the materials that are involved and the loss spectra reflect the geometry of the sample interfaces.

The geometry becomes particularly important in photonic crystals, where the periodic spatial modulation of the dielectric function affects the propagation of radiation to the point of forbidding it within band-gap energy regions. As a consequence, photonic crystals are known to inhibit the spontaneous emission of light within the band gap.  They can also be used to make omnidirectional dielectric mirrors that reflect light from all directions without absorption  and wave guides able to deflect light around sharp corners on the scale of the wavelength.  These are applications of photonic crystals that involve free external radiation, but equally remarkable effects are expected to accompany evanescent fields like those of external electrons.

In this work, we examine the energy loss spectra of electrons moving near or inside photonic crystals over a wide range of lattice parameters $`a`$. When $`a`$ is much smaller than the wavelength of the radiation associated to a given energy loss, $`\lambda `$, it is shown that the crystal can be described to some extend by an effective dielectric function.  For larger lattice constants, the photonic band structure becomes more complicated and this is reflected in the loss spectra. The relation between the band structure, the reflectance of photonic crystal slabs, and the electron energy loss spectra is discussed in detail in what follows.

In a previous development, Pendry and Martín-Moreno  calculated the energy loss probability for electrons moving near crystals made of either metallic spheres or metallic cylinders in the $`a<<\lambda `$ limit using the transfer matrix approach  to solve Maxwell’s equations. Their results exhibit a complicated energy-loss structure even for relatively dilute crystals , which has not been reproduced in the present work. Therefore, in order to test the convergence of our method, we have calculated the energy loss near the surface of a crystal of dilute spheres and have found very good agreement between the results derived from the theory described in this work and those obtained from the analytical expression for isolated spheres.

In addition to producing energy loss, the interaction between the electron and the crystal gives rise to the emission of the so-called Smith-Purcell radiation.  This effect has already been discussed for one-dimensional  and three-dimensional  crystals. Radiation emission is one of the contributions to the total energy loss, and in frequency regions where the crystal is transparent the probability of these two must coincide. Examples of this are offered below.

Here, the energy loss and photon emission probability are calculated from the induced electric field, which is derived using the reflection coefficients of the crystal,  as explained in Sec. II. Results for the case of small lattice constants ($`a<<\lambda `$) are given in Sec. III and for larger ones ($`a\lambda `$) in Secs. IV and V. The main conclusions are summarized in Sec. VI. Gaussian atomic units (a.u., that is, $`\mathrm{}=m=e=1`$) will be used from now on, unless otherwise specified.

## II Theory

### A The reflection coefficients of the crystal

We shall consider crystals that are composed of a number of layers perpendicular to the $`z`$ direction and extending from $`z=0`$ towards $`z<0`$. The host region outside the crystal will be assumed to be described by a frequency-dependent dielectric function $`ϵ_h(\omega )`$ and a magnetic permeability $`\mu _h(\omega )`$. Each crystal layer consists of the repetition of a given object with certain two-dimensional translational symmetry that is shared by all layers. The crystal will be characterized here by its reflection coefficients, that is, the amplitudes of the reflected plane wave components for a given incident wave.

In general, an external plane wave of frequency $`\omega `$ that propagates near the crystal can be represented as

$`𝐄_𝐐^\pm \mathrm{exp}(\mathrm{i}𝐊_𝐐^\pm 𝐫),`$ (1)
where $`𝐐=(Q_x,Q_y)`$,

$`𝐊_𝐐^\pm =(𝐐,\pm \mathrm{i}\mathrm{\Gamma }_Q),`$ (2)
$`\mathrm{\Gamma }_Q^2=Q^2k_h^2,`$ (3)
and $`k_h=(\omega /c)(ϵ_h\mu _h)^{1/2}`$ is the momentum of the plane wave. The real part of $`\mathrm{\Gamma }_Q`$ is chosen to be positive, and the $`+`$ ($``$) sign in these expressions stands for a wave moving towards $`z>0`$ ($`z<0`$).

The momentum $`𝐐`$ parallel to the surface has been singled out to make explicit use of crystal symmetry: invoking momentum conservation, a wave incident from the $`z>0`$ region with momentum $`𝐊_𝐐^{}`$ will only produce a discrete set of reflected waves of momentum $`𝐊_{𝐐+𝐆}^+`$, where $`𝐆`$ runs over reciprocal surface lattice vectors.

The transversal character of these waves (i.e., the fact that $`𝐊_𝐐^\pm 𝐄_𝐐^\pm =0`$) can be exploited to express the electric field in terms of s and p components according to

$`𝐄_𝐐^\pm =E_{𝐐,s}^\pm \widehat{𝐞}_{𝐐,s}^\pm +E_{𝐐,p}^\pm \widehat{𝐞}_{𝐐,p}^\pm ,`$ (4)
where the vectors

$`\widehat{𝐞}_{𝐐,s}^\pm ={\displaystyle \frac{1}{Q}}(Q_y,Q_x,0)`$ (5)
and

$`\widehat{𝐞}_{𝐐,p}^\pm ={\displaystyle \frac{\mathrm{i}}{k_hQ}}(\pm \mathrm{\Gamma }_QQ_x,\pm \mathrm{\Gamma }_QQ_y,\mathrm{i}Q^2)`$ (6)
satisfy the identities $`\widehat{𝐞}_{𝐐,s}^\pm \widehat{𝐞}_{𝐐,s}^\pm =\widehat{𝐞}_{𝐐,p}^\pm \widehat{𝐞}_{𝐐,p}^\pm =1`$ and $`\widehat{𝐞}_{𝐐,s}^\pm \widehat{𝐞}_{𝐐,p}^\pm =\widehat{𝐞}_{𝐐,s}^\pm 𝐊_𝐐^\pm =\widehat{𝐞}_{𝐐,p}^\pm 𝐊_𝐐^\pm =0`$. Also, $`\widehat{𝐞}_{𝐐,s}^\pm `$ is perpendicular to the plane defined by $`𝐐`$ and the surface normal.

The amplitudes of the reflected waves depend linearly on the amplitudes of s and p components of the incident wave, and the coefficients of the linear relation between them are the reflection coefficients $`R_{\mathrm{𝐐𝐆}}^{\sigma \sigma ^{}}`$, implicitly defined by

$`[E_{𝐐+𝐆,\sigma }^+]^r={\displaystyle \underset{\sigma ^{}}{}}R_{\mathrm{𝐐𝐆}}^{\sigma \sigma ^{}}[E_{𝐐,\sigma ^{}}^{}]^i,`$ (7)
where $`\sigma `$ and $`\sigma ^{}`$ run over polarization directions s and p, and the super-indices $`r`$ and $`i`$ stand for reflected and incident components, respectively.

We have used the layer KKR method to calculate the reflection coefficients $`R_{\mathrm{𝐐𝐆}}^{\sigma \sigma ^{}}`$ both for $`a<<\lambda `$ and for $`a\lambda `$. In the layer KKR method, the transmission and reflection coefficients are calculated exactly for each single layer using multiple scattering in a basis set of multipoles centered around each object of the layer. Scattering among layers is then expressed in terms of those coefficients. The maximum multipole order $`l_{\mathrm{max}}`$ and the number of reflected and transmitted beams (i.e., the number of $`𝐆`$ vectors) are the only convergence parameters, which have been tested in all calculated results shown below. Stefanou et al.  developed this method for spheres and we have extended it to be able to deal with non-spherical objects and arbitrary values of $`l_{\mathrm{max}}`$. Further details of the method will be given elsewhere.

### B The field of the electron

We shall consider an electron moving with constant velocity $`v`$ along a trajectory parallel to the crystal surface and described by $`𝐫_t=(vt,y_0,z_0)`$, with $`z_0>0`$, so that the electron moves in the host medium described by $`ϵ_h`$ and $`\mu _h`$. Neglecting the crystal for the moment, the electron electric field can be expressed in frequency space $`\omega `$ in terms of vector and scalar potentials as

$`𝐄_0={\displaystyle \frac{\mathrm{i}\omega }{c}}𝐀_0\varphi _0.`$ (8)
In the Lorentz gauge, Maxwell’s equations can be recast as

$`(^2+k_h^2)\varphi _0={\displaystyle \frac{4\pi }{ϵ_h}}\rho `$ (9)
and

$`(^2+k_h^2)𝐀_0={\displaystyle \frac{4\pi \mu _h}{c}}𝐣,`$ (10)
where $`\rho (𝐫,\omega )`$ is the electron charge density, and $`𝐣(𝐫,\omega )=(v\rho ,0,0)`$ is its charge current.

Using the relation

$`(^2+k_h^2){\displaystyle \frac{d𝐩}{2\pi ^2}\frac{\mathrm{e}^{\mathrm{i}𝐩𝐫}}{p^2k_h^2\mathrm{i0}^+}}=4\pi \delta (𝐫)`$ (11)
and noticing that $`\rho (𝐫,t)=\delta [𝐫(vt,y_0,z_0)]`$, the electric field is found to be

$`𝐄_0(𝐫,\omega )`$ $`=`$ $`[{\displaystyle \frac{}{ϵ_h}}{\displaystyle \frac{\mathrm{i}\omega v\mu _h}{c^2}}\widehat{𝐱}]`$ (13)
$`\times {\displaystyle }{\displaystyle \frac{d𝐩}{2\pi ^2}}{\displaystyle }dt\mathrm{e}^{\mathrm{i}\omega t}{\displaystyle \frac{\mathrm{e}^{\mathrm{i}𝐩[𝐫(vt,y_0,z_0)]}}{p^2k_h^2\mathrm{i0}^+}},`$
where the time integral represents the inverse Fourier transform that permits obtaining $`\rho (𝐫,\omega )`$ in terms of $`\rho (𝐫,t)`$. The above integral can be reduced to

$`𝐄_0(𝐫,\omega )={\displaystyle 𝑑Q_y\mathrm{e}^{\mathrm{i}𝐊_𝐐^\pm [𝐫(0,y_0,z_0)])}𝐄_𝐐^\pm },`$ (14)
where $`𝐐=(\omega /v,Q_y)`$ and

$`𝐄_𝐐^\pm ={\displaystyle \frac{\mathrm{i}}{\mathrm{\Gamma }_Q}}({\displaystyle \frac{𝐊_𝐐^\pm }{vϵ_h}}{\displaystyle \frac{\omega \mu _h}{c^2}}\widehat{𝐱}).`$ (15)
The $`+`$ ($``$) sign must be used in these expressions when $`z>z_0`$ ($`z<z_0`$), so that the integrand of Eq. (14) represents a plane wave that propagates towards positive (negative) $`z`$’s. When the electron is moving in vacuum, $`\mathrm{\Gamma }_Q`$ is real and the waves in the integrand of Eq. (14) are evanescent.  However, when $`k_h`$ is real and larger than $`\omega /v`$, some of those waves describe Cherenkov radiation that propagates without attenuation; this is the case of electrons that travel faster than light in the medium.

### C The field induced by interaction of the electron and the crystal

The decomposition of these plane waves into s and p components is readily found to be

$`[E_{𝐐,s}^\pm ]^i={\displaystyle \frac{\mathrm{i}Q_y\omega \mu _h}{Q\mathrm{\Gamma }_Qc^2}}\mathrm{e}^{\mathrm{i}𝐊_𝐐^\pm (0,y_0,z_0)}`$ (16)
and

$`[E_{𝐐,p}^\pm ]^i=\pm {\displaystyle \frac{k_h}{vQϵ_h}}\mathrm{e}^{\mathrm{i}𝐊_𝐐^\pm (0,y_0,z_0)}.`$ (17)
Each of the incident plane waves \[i.e., each value of $`Q_y`$ in Eq. (14)\] gives rise to a set of reflected waves whose amplitudes are obtained from Eq. (7). Therefore, the electric field in the region near the ion can be constructed as the sum of $`𝐄_0`$ and the reflected field,

$`𝐄(𝐫,\omega )=𝐄_0(𝐫,\omega )+𝐄_r(𝐫,\omega ),`$ (18)
where

$`𝐄_r(𝐫,\omega )=`$ (19)
$`{\displaystyle \underset{𝐆,\sigma \sigma ^{}}{}}{\displaystyle 𝑑Q_y\mathrm{e}^{\mathrm{i}𝐊_{𝐐+𝐆}^+𝐫}R_{\mathrm{𝐐𝐆}}^{\sigma \sigma ^{}}[E_{𝐐,\sigma ^{}}^{}]^i\widehat{𝐞}_{𝐐+𝐆,\sigma }^+},`$ (20)
and the integral over $`Q_y`$ has been copied directly from Eq. (14) in virtue of the linearity of Maxwell’s equations.

### D Electron energy loss in front of the crystal

The electron energy loss can be calculated from the retarded force exerted by the induced part of the electric field $`𝐄^{\mathrm{ind}}`$ acting back on the electron. Integrating this force along the trajectory and dividing by the total path length $`L`$, one finds

$`{\displaystyle \frac{\mathrm{\Delta }E}{\mathrm{\Delta }x}}={\displaystyle \frac{v}{L}}{\displaystyle 𝑑tE_x^{\mathrm{ind}}(𝐫_t,t)}={\displaystyle _0^{\mathrm{}}}\omega 𝑑\omega P(\omega ),`$ (21)
where

$`P(\omega )={\displaystyle \frac{v}{\pi \omega L}}{\displaystyle 𝑑t\mathrm{Re}\{\mathrm{e}^{\mathrm{i}\omega t}E_x^{\mathrm{ind}}(𝐫_t,\omega )\}}`$ (22)
is the loss probability per unit of path length.

For an electron moving parallel to a crystal surface, Eq. (18) permits separating the loss probability as

$`P=P_0+P_r,`$ (23)
where

$`P_0(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{\pi v^2}}\mathrm{Im}\{({\displaystyle \frac{v^2}{c^2}}\mu _h{\displaystyle \frac{1}{ϵ_h}})\mathrm{log}[{\displaystyle \frac{q_c^2k_h^2}{(\omega /v)^2k_h^2}}]\}`$ (24)
corresponds to the contribution of $`𝐄_0`$ in the absence of the crystal. Here, $`q_c`$ is a momentum cut-off related to energy conservation. This contribution vanishes in vacuum. However, it gives rise to Cherenkov losses when $`ϵ_h`$ and $`\mu _h`$ are real, in which case  $`P_0=|\mu _h|[1/c^21/(v^2ϵ_h\mu _h)]`$, subject to the Cherenkov condition $`v^2ϵ_h\mu _h>c^2`$.

The remaining contribution $`P_r`$ is due to the retarding force exerted by the reflected field $`𝐄_r`$. Inserting Eq. (20) into Eq. (22), the time $`t`$ and the lateral impact parameter $`y_0`$ appear only through the factor $`\mathrm{exp}\{\mathrm{i}(G_xvt+G_yy_0)\}`$, so that time integration eliminates all vectors $`𝐆`$ with $`G_x0`$ from the sum in Eq. (20). Furthermore, averaging over $`y_0`$ leaves only the $`G_y=0`$ term and one obtains

$`P_r(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{\pi v^2}}{\displaystyle \frac{dQ_y}{Q^2}}`$ (25)
$`\times `$ $`\mathrm{Im}\{[R_{𝐐0}^{pp}+ϵ_h\mu _h({\displaystyle \frac{Q_yv}{\mathrm{\Gamma }_Qc}})^2R_{𝐐0}^{ss}`$ (27)
$`+{\displaystyle \frac{\mathrm{i}Q_yk_hv}{\omega \mathrm{\Gamma }_Q}}(R_{𝐐0}^{sp}R_{𝐐0}^{ps})]{\displaystyle \frac{\mathrm{\Gamma }_Q\mathrm{e}^{2\mathrm{\Gamma }_Qz_0}}{ϵ_h}}\}.`$
The last term inside the square bracket of this expression gives no contribution when the electron trajectory is contained in a plane of specular symmetry of the crystal surface.

As a test one can apply Eq. (25) to an electron moving in vacuum ($`ϵ_h=\mu _h=1`$) at a distance $`z_0`$ from the surface of a non-magnetic medium described by $`ϵ`$, in which case $`R_{𝐐0}^{sp}=R_{𝐐0}^{ps}=0`$. Then, Eq. (25) reduces to

$`P(\omega )`$ $`=`$ $`{\displaystyle \frac{2}{\pi v^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Gamma }_QdQ_y}{Q^2}}\mathrm{e}^{2\mathrm{\Gamma }_Qz_0}`$ (28)
$`\times `$ $`\mathrm{Im}\{{\displaystyle \frac{ϵ\mathrm{\Gamma }_Q\mathrm{\Gamma }_Q^{}}{ϵ\mathrm{\Gamma }_Q+\mathrm{\Gamma }_Q^{}}}+({\displaystyle \frac{Q_yv}{\mathrm{\Gamma }_Qc}})^2{\displaystyle \frac{\mathrm{\Gamma }_Q\mathrm{\Gamma }_Q^{}}{\mathrm{\Gamma }_Q+\mathrm{\Gamma }_Q^{}}}\},`$ (29)
where $`\mathrm{\Gamma }_Q^{}=\sqrt{Q^2ϵ\omega ^2/c^2}`$, $`\mathrm{\Gamma }_Q=\sqrt{Q^2\omega ^2/c^2}`$, and the reflection coefficients of Eq. (25) have been taken from Fresnel’s equations.  In the non-relativistic limit, $`\mathrm{\Gamma }_Q=\mathrm{\Gamma }_Q^{}=Q`$, so that the contribution of the reflection of s waves vanishes and the loss probability becomes proportional to $`\mathrm{Im}\{1/(ϵ+1)\}`$.

## III Energy loss for small lattice constants

When the crystal lattice constant is much smaller than the wavelength corresponding to a given frequency component $`\omega `$, the details of the crystal lattice cannot be resolved by the electron, so that the medium behaves like a uniform material, characterized by an effective dielectric constant and magnetic permeability. For instance, a dilute system of spheres can be regarded as a set of interacting dipoles, which leads to the well-known Maxwell-Garnett formula.  For more compact systems the details of the microscopic structure becomes relevant via important multipolar interactions.

These are the cases considered in Fig. 1 for a crystal composed of small aluminum spheres that are disposed in a simple cubic lattice of constant $`a=5`$ nm. The solid curves of Fig. 1(a) show the loss probability for an electron moving with velocity $`v=0.4c`$ parallel to the surface of a slab made up of six $`(100)`$ layers of such crystal. The lost probability has been calculated by means of Eq. (25) for two different filling fractions $`f`$ of the aluminum ($`6.5\%`$ and $`30\%`$, respectively). The dielectric function of aluminum has been approximated by a Drude expression $`ϵ(\omega )=1\omega _p^2/\omega (\omega +\mathrm{i}\eta )`$ with $`\omega _p=15`$ eV and $`\eta =1`$ eV. The wavelengths considered in the figure lie in the range $`\lambda =2\pi c/\omega =88.6310`$ nm, so that we are in the $`\lambda >>a`$ limit. The maximum orbital quantum number used to achieve convergence on the figure is $`l_{\mathrm{max}}=6`$.

Crystalline effects come from the interaction among spheres, and a method to determine their relative role consists in comparing these results with a calculation in which that interaction is suppressed. This is what the broken curves stand for. It is very clear that the interaction among spheres is almost negligible for $`f=6.5\%`$, whereas it becomes very important at $`f=30\%`$. In the latter case, the sphere surfaces are closer to each other and their mutual electromagnetic coupling becomes relevant. The non-interacting case shows a peak at around $`\omega _1=8.7`$ eV that corresponds to the dipole Mie resonance in the small-sphere limit, given by the expression $`ϵ(\omega _1)=2`$ (i.e., $`\omega _1=\omega _p/\sqrt{3}`$). The main effect of the sphere-sphere interaction (solid curve for $`f=30\%`$) consists in splitting this peak, in a similar way as splitting of degenerate levels (the Mie resonances here) occurs in atomic bonding.

This is actually observed in the projected photonic band structure of this crystal, represented in Fig. 2 for $`f=30\%`$. The structure is dominated by nearly flat bands corresponding to localized excitations near the Mie resonances of the small isolated spheres, $`\omega _l=\omega _p\sqrt{l/(l+1)}`$, which lie in the 8.7-10.6 eV range (see the labels on the right hand side of the figure for the dipole and quadrupole Mie modes). The interaction between spheres gives rise to a complex structure that encompasses regions of relative transparency outside that energy range. The evanescent plane wave components of the external electron field are subject to the condition $`\omega =𝐐𝐯`$, which defines the straight lines shown in the figure. The loss spectrum can then be understood as originating from both absorption and direct coupling to propagating modes of the crystal. This last effect is clearly seen as a bump in the loss probability near 6.7 eV, connected to the crossing of the noted straight lines with a low-energy region of allowed propagating modes.

The dashed curves of Fig. 1(a) (uncoupled spheres) have been obtained in two different ways. The first one consists in calculating the reflection coefficients that appear in Eq. (25) using the layer KKR method but neglecting the interaction among spheres within each layer and also the interaction among different layers. A second procedure consists in making use of the analytical expression for the energy loss probability of an electron near an isolated sphere (see Refs. \[\] and \[\] for non-relativistic and fully-relativistic formulas, respectively) and summing over all impact parameters of the electron trajectory with respect to the spheres of the crystal. The results coming out of these two very different procedures cannot be distinguished on the scale of the figure, and this is a strong indication of the degree of convergence of our numerical calculations with respect to the number of plane waves used in the layer KKR method.

These results do not support, however, previous calculations by Pendry and Martín-Moreno  for this exact system, where the number of features and their energy position for the filling fractions under consideration differ from ours. A possible lack of convergence of the transfer matrix method used by the authors for this three-dimensional system of metallic scatterers might be the reason of this discrepancy. In particular, some of their low-energy features around 6.5 eV could originate in the modes of the wedge associated to their space discretization procedure (resonances near that energy have also been found for the $`90^{}`$ aluminum wedge ).

The example offered in Fig. 1 illustrates what happens in the $`a<<\lambda `$ limit, where it is reasonable to define an effective dielectric function $`ϵ_{\mathrm{eff}}`$ for the infinite crystal and to compare the results of the detailed, exact calculation with those obtained for an electron moving parallel to a homogeneous surface of a material described by such a dielectric function. Several recipes for defining $`ϵ_{\mathrm{eff}}`$ for granular materials have been given in the literature, ranging from simple effective medium theories like Maxwell-Garnett’s  to more elaborate ones that take into account the actual shape of the constituents, both for disordered composites  and for crystals.  Among the latter, one finds extensions of the Maxwell-Garnett theory that go beyond dipolar interactions,  spectral representations,  or direct derivation of the light dispersion relation.  Here, we have used a different method that consists in finding the dielectric function $`ϵ_{\mathrm{eff}}`$ of the equivalent homogeneous medium that has the same reflectance as the crystal surface;  the reflectance calculated for the crystal is well reproduced by Fresnel equations with a single value of $`ϵ_{\mathrm{eff}}(\omega )`$ for each frequency $`\omega `$ within a $`2\%`$ under the present conditions. The results are represented by the solid curves of Fig. 1(b), where they are compared with Maxwell-Garnett’s theory (broken curves). For $`f=6.5\%`$ these two models are nearly identical, as Maxwell-Garnett formula gives the dipole of the small isolated sphere correctly. However, for $`f=30\%`$ the deviation between both models is significant, since higher momenta are involved in the interaction between neighboring spheres . The magnitude actually represented in the figure is the surface loss function, $`\mathrm{Im}\{1/(1+ϵ_{\mathrm{eff}})\}`$, which is directly comparable to the loss probability of Fig. 1(a). The agreement in the position of the peaks between the detailed energy loss calculation and our effective medium theory is reasonably good, indicating that an effective homogeneous medium describes the solid appropriately within this energy range at $`v=0.4c`$; however, the relative weight of the features for $`f=30\%`$ changes completely at lower velocities (dotted curve, $`v=0.06c=8.2`$ a.u.).

Another example of a crystal that can be represented by an effective dielectric constant is given in Fig. 3, consisting of a simple cubic lattice with the same parameters as in Fig. 1, except that finite cylinders have been used instead of spheres. As a result, the medium is strongly anisotropic and characterized by different bulk plasmon modes when the electric field is directed parallel or perpendicular to the cylinders \[see Fig. 3(c)-(d)\]. The loss probability \[Fig. 3(a)\] is shown to share most of the features of the surface loss function for anisotropic media \[Fig. 3(b)\], except for the 7.8 eV peak in the loss spectrum with $`f=30\%`$, which might be connected to the proximity effect when the electron starts sensing the non-uniform character of the surface via evanescent waves.

## IV Energy loss for lattice constants comparable to the wavelength

For lattice constants comparable to the wavelength associated to the energy losses under consideration, one can no longer define the effective dielectric constant of an equivalent homogeneous medium. Then, it is useful to relate the loss spectra directly to the photonic band structure. This has been done in Fig. 4 for a crystal of aluminum spheres immersed in a dielectric with $`ϵ=3`$. The features of the loss spectra \[Fig. 4(c)\] are strongly correlated with the band structure \[Fig. 4(a)-(b)\]. In particular, the pseudo-gaps near 3.2 eV are translated into a dip in the loss probability.

Using electrons to analyze the crystal permits exploring regions of the band structure that are not accessible to external light, but that can be reached via the evanescent waves contained in the perturbing field of the electrons. This is the case of Fig. 4 for $`v/c=0.5`$ \[solid curve in Fig. 4(c)\], well below the Cherenkov threshold $`v/c=1/\sqrt{3}`$, that permits to observe how the peak in the loss spectrum near 2.7 eV results from the coupling with propagating modes of the crystal as a result of the intersection of the external field components \[straight solid lines in Fig. 4(a)\] and the allowed regions of propagation (shaded areas).

At larger velocities \[$`v/c=0.75`$, dashed curve in Fig. 4(c)\], Cherenkov radiation would be produced in the absence of the crystal, that would result in a loss probability independent of $`\omega `$ (see little arrow near the horizontal axis of the figure). This is clearly seen in the big overlap of the external field with the regions of propagation within the crystal inside the low-energy region results in an enhanced loss probability, as compared to $`v/c=0.5`$. However, the presence of the crystal modulates the loss spectrum.

Part of the energy lost by the electron must be converted into the so-called Smith-Purcell radiation.  In particular, when transparent materials are used to build the crystals, the light emission probability must coincide with the energy-loss probability. For a crystal of finite thickness, like the one considered in Fig. 5, consisting of 8 (111)-layers of an inverted Si opal ($`ϵ=11.9`$) with a filling fraction of $`67\%`$, part of this emission occurs towards the side of the crystal opposite to the electron trajectory (here, the electron is taken to be moving parallel to the crystal surface). The intensity of the emitted light \[Fig. 5(c)\] has been calculated from the integral of the Poynting vector far from the crystal (see Ref. \[\] for more details), and it presents a strong dip near $`a/\lambda 0.75`$, where $`a=1220`$ nm is the lattice constant. The region of emission depletion is actually contained within a full band-gap of the infinite crystal \[see Fig. 5(a)\]; this is also seen in the transmission of light incident on the crystal both normal to the surface \[dashed curve in Fig. 5(b)\] or with an angle corresponding to the Cherenkov radiation produced by the electron in the plane defined by the trajectory and the surface normal (solid curve).

## V Energy loss in the bulk of a crystal

So far we have considered an electron moving in front of the surface of a photonic crystal. When the electron is moving in the bulk of an infinite crystal, the electric field can also be written in terms of reflectance coefficients if the trajectory is contained in a low-index plane that does not intersect any of the crystal objects. The reflectance in question is that of the lower and upper semi-infinite crystals into which the entire crystal is divided by the noted plane. The corresponding reflectance matrices will be denoted $`R_1`$ and $`R_2`$, respectively. They contain the amplitudes of reflected beams, labeled by vectors $`𝐆`$ of the 2D reciprocal lattice of the plane under consideration. In particular, for a fixed value of $`𝐐`$, $`R_{1,\mathrm{𝐆𝐆}^{}}^{\sigma \sigma ^{}}`$ is given by $`R_{𝐐+𝐆^{},𝐆𝐆^{}}^{\sigma \sigma ^{}}`$, as defined by Eq. (7). Similarly, $`R_{2,\mathrm{𝐆𝐆}^{}}^{\sigma \sigma ^{}}`$ connects $`𝐊_{𝐐+𝐆}^+`$ wave components with waves reflected from the upper surface with momentum $`𝐊_{𝐐+𝐆^{}}^{}`$.

This is represented schematically in Fig. 6, where the reflectance matrices of the lower and upper semi-infinite crystals are $`R_1`$ and $`R_2`$, respectively, and upwards (downwards) arrows represent waves of momentum $`𝐊_{𝐐+𝐆}^+`$ ($`𝐊_{𝐐+𝐆}^{}`$). The upper and lower surfaces have been separated in the figure for the sake of clarity, but they will be considered to lie on the same plane (a plane that contains the electron trajectory) in what follows.

Using Eq. (22) for the loss probability and a straightforward extension of Eq. (20) for the reflected electric field including all terms of Fig. 6, the loss probability per unit of path length averaged over all impact parameters parallel to the crystal surfaces is found to be

$`P_r(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{\pi v^2}}{\displaystyle \frac{dQ_y}{Q^2}}`$ (30)
$`\times `$ $`\mathrm{Im}\{[A_{00}^{pp}+B_{00}^{pp}C_{00}^{pp}D_{00}^{pp}`$ (34)
$`+ϵ_h\mu _h({\displaystyle \frac{Q_yv}{\mathrm{\Gamma }_Qc}})^2(A_{00}^{ss}+B_{00}^{ss}+C_{00}^{ss}+D_{00}^{ss})`$
$`+{\displaystyle \frac{\mathrm{i}Q_yk_hv}{\omega \mathrm{\Gamma }_Q}}(A_{00}^{sp}B_{00}^{sp}+C_{00}^{sp}D_{00}^{sp}`$
$`A_{00}^{ps}+B_{00}^{ps}+C_{00}^{ps}D_{00}^{ps})]{\displaystyle \frac{\mathrm{\Gamma }_Q}{ϵ_h}}\},`$
where the matrices $`A_{\mathrm{𝐆𝐆}^{}}^{\sigma \sigma ^{}}`$, $`B_{\mathrm{𝐆𝐆}^{}}^{\sigma \sigma ^{}}`$, $`C_{\mathrm{𝐆𝐆}^{}}^{\sigma \sigma ^{}}`$, and $`D_{\mathrm{𝐆𝐆}^{}}^{\sigma \sigma ^{}}`$ are defined in Fig. 6, and only the elements $`𝐆=𝐆^{}=0`$ enters this expression.

The above formula has been applied to calculate the energy loss probability for an electron moving in between two parallel crystal slabs of 32 layers each surrounded by Si, as shown in Fig. 7. The crystal parameters are the same as in Fig. 5(c), so that an absolute band gap shows up near $`a/\lambda =0.75`$. The electron energy (100 keV) is large enough to produce Cherenkov light, which must be confined in between the two crystals for wavelengths lying in the gap region. This is actually observed as a pronounced dip, connected to the fact that no electromagnetic modes can be created that escape through the crystals, except for a small transmission due to their finite thickness. However, some radiation can escape though modes that are trapped in the slab formed by the two crystals, so that the loss probability does not actually reach a zero value in the dip.

## VI Concluding remarks

The energy loss of fast electrons interacting with photonic crystals has been calculated for dielectric, metallic, and metalo-dielectric systems.

When the wavelength of the radiation associated to the energy loss is much larger than the lattice constant, the crystal can be regarded as a continuous medium, characterized by an effective dielectric function, as it has been shown for crystals of small aluminum spheres in Fig. 1 and small finite cylinders in Fig. 3. The effective medium is highly anisotropic in the latter case. The effective dielectric function has been calculated in both cases from the reflectance coefficient of the crystal, and this has been shown to contain most of the information needed to understand the calculated energy-loss spectra.

When the lattice constant is comparable to the wavelength, the idea of an equivalent effective continuous medium is no longer valid, and one has to rely on the detailed band structure of the crystal to understand the loss spectra \[Figs. 4, 5, and 7\].

Finally, the interaction of the electron with the crystal produces Smith-Purcell radiation, which contributes to the energy loss. If the crystal is composed of non-absorbing materials, the energy loss probability and the light emission probability must coincide, as shown in Fig. 5 for an inverted Si opal, in which case the light emitted after transmission through the crystal shows dips that are directly connected to the presence of photonic band gaps. For electrons moving in a region surrounded by photonic crystals a dip is also observed in the loss probability within the gap energies (Fig. 7).

Our hope is that the present work can provide a stimulus to use fast electrons in the analysis of photonic crystals as a way to bring a source of evanescent radiation (the electromagnetic field of the electron in vacuum) into close contact with the crystal and also to probe regions of the crystal that would not be easily accessible to other sources of external electromagnetic radiation.

###### Acknowledgements.

The authors gratefully acknowledge help and support from the Basque Departamento de Educación, Universidades e Investigación, the University of the Basque Country UPV/EHU (contract No. 00206.215-13639/2001), and the Spanish Ministerio de Ciencia y Tecnología (contract No. MAT2001-0946).