format: text

text:
# Confusing non-standard neutrino interactions with oscillations at a neutrino factory

## I Introduction

From the long-standing solar solar  and atmospheric atmos  neutrino anomalies we now have compelling evidence that an extension of the Standard Model (SM) of particle physics is necessary in the lepton sector. The simplest and most generic explanation of these experiments is provided by neutrino oscillations induced by neutrino masses and mixing. As is well known and accepted, the indication of the LSND experiment LSND  for oscillations at a large mass-squared difference can not be reconciled with solar and atmospheric data within a 3-neutrino framework. For recent four-neutrino Peltoniemi:1992ss  analyses see Ref. Maltoni:2001bc  and references therein. For this reason we choose not to consider the LSND data and focus therefore on the simplest 3-neutrino scheme with the two mass-squared differences $`\mathrm{\Delta }m_{\mathrm{sol}}^210^4`$ eV2 and $`\mathrm{\Delta }m_{\mathrm{atm}}^23\times 10^3`$ eV2 Gonzalez-Garcia:2001sq . The corresponding lepton mixing matrix is parameterized by the three angles $`\theta _{12},\theta _{23},\theta _{13}`$ and one complex CP-violating phase $`\delta `$ (see later for exact definitions) relevant in lepton-number-conserving neutrino oscillations theory ; theory1 . It is known from atmospheric neutrino data that $`\theta _{23}`$ has to be nearly maximal. On the other hand, data from present solar neutrino experiments favor a large value for the angle $`\theta _{12}`$ Gonzalez-Garcia:2000aj ; solaranalysis . An improved determination is expected from solar neutrino data and/or the results of the KamLAND experiment kamland . The value of the third angle $`\theta _{13}`$ is not known at present, there is only the bound

$$\mathrm{sin}^22\theta _{13}0.1\text{at}90\%\text{CL}$$
(1)
implied by combining the results of the reactor experiments reactor  CHOOZ and Palo Verde with atmospheric data. Together with the requirement of solar neutrino oscillations this bound implies that $`\mathrm{sin}^2\theta _{13}`$ has to be small. The value of the phase $`\delta `$ is completely unknown.

Currently a new generation of long-baseline neutrino oscillation experiments using a neutrino beam originating from the decay of muons in a storage ring is being discussed nuFacPapers . These so-called neutrino factories are considered as the ideal tool to enhance our knowledge about neutrino mixing parameters. Besides the possibility to explore CP violation in the lepton sector an important aim of a neutrino factory will be a precise determination of $`\theta _{13}`$; it is claimed that a measurement of $`\theta _{13}`$ down to values of a few$`\times 10^4`$ will be possible Freund:2001ui .

In a large class of models beyond the Standard Model non-standard interactions (NSI) of neutrinos with matter arise. The simplest NSI do not require new interactions beyond those mediated by the Standard Model electroweak gauge bosons: it is simply nature of the leptonic charged and neutral current interactions which is non-standard due to the complexity of neutrino mixing theory . On the other hand NSI can also be mediated by the exchange of new particles with mass at the weak scale such as in some super-symmetric models with R-parity violating Ross:1985yg ; Hall:1984id  interactions. Such non-standard flavor-violating physics can arise even in the absence of neutrino mass NSImodels2 ; NSImodels3  and can lead to non-universal (NU) or to flavor-changing (FC) neutrino interactions.

Non-standard interactions of neutrinos affect their propagation in matter and the magnitude of the effect depends on the interplay between conventional mass-induced neutrino oscillation features in matter MSW  and those genuinely engendered by the NSI, which do not require neutrino mass Valle:1987gv . Correspondingly their implications have been explored in a variety of contexts involving solar neutrinos MSW ; Valle:1987gv ; NSIrecent ; Bergmann:2000gp ; Guzzo:2001mi , atmospheric neutrinos Guzzo:2001mi ; NSIatm ; Val ; BergmannAtm , other astrophysical sources NSIastro ; NSIastro1  and the LSND experiment BergmannLSND ; Bergmann:2000gn . The impact of non-standard interactions of neutrinos has also been considered from the point of view of future experiments involving solar neutrinos Berezhiani:2001rt  as well as the upcoming neutrino factories Huber:2001zw ; Huber:2001de . Various aspects of NSI for a neutrino factory experiment have been considered in Refs. Huber:2001zw ; Huber:2001de ; Gago:2001xg ; Gonzalez-Garcia:2001mp ; Datta:2000ci ; johnson ; ota ; Bueno:2001jy .

In this paper we will consider from a phenomenological point of view the impact of NSI on the determination of neutrino mixing parameters at a neutrino factory. In particular we will focus on the $`\stackrel{()}{\nu }_e\stackrel{()}{\nu }_\mu `$ channels, which are supposed to be the “golden channels” for the measurement of $`\theta _{13}`$. We extend our previous work Huber:2001de  by taking simultaneously into account neutrino oscillations and the effect of NSI in neutrino source, propagation and detection ota . We will show that a certain combination of FC interactions in source and propagation can give exactly the same signal as oscillations arising due to $`\theta _{13}`$. This implies that information about $`\theta _{13}`$ can only be obtained if bounds on NSI are available. In view of the existing bounds on FC interactions, this leads to a drastic loss in sensitivity in $`\theta _{13}`$, at least two orders of magnitude.

All our considerations also apply to the determination of $`\theta _{13}`$ in long baseline experiments using upgraded conventional beams superbeams ; ota . However, due to the different production processes involved, the NSI parameters relevant in the source of a conventional beam experiment differ from the ones at a neutrino factory. The same methods discussed here can be adapted to cover also that case. A detailed numerical consideration of conventional beam experiments goes beyond the scope of this work.

The outline of the paper is as follows. In section II we briefly sketch the theoretical motivation for NSI in the context of gauge theories of neutrino mass. In section III we discuss examples of low energy four-fermion Hamiltonians, which lead to NSI in neutrino source, propagation and detection. In section IV we review some bounds on NSI parameters obtained in the literature. In section V we present the framework of our numerical calculations and discuss the appearance rate in the presence of NSI and oscillations. In section VI we derive analytical expressions for this rate and formulate the oscillation–NSI confusion theorem. In section VII we describe the simulation of a neutrino factory and our statistical method to investigate the possibilities of such an experiment to determine NSI and oscillation parameters. In section VIII we define sensitivity limits for $`\mathrm{sin}^22\theta _{13}`$ and show our numerical results for the three baselines 700 km, 3000 km and 7000 km as a function of bounds on the relevant NSI parameters. We also discuss the sensitivity limits if information from two different baselines is combined. Finally we conclude in section IX.

## II Theoretical motivation

More often than not, models of neutrino mass are accompanied by NSI, leading generically to both oscillations and neutrino NSI in matter. The simplest are those NSI which arise from neutrino-mixing. The most straightforward example of this case is when neutrino masses follow from the admixture of isosinglet neutral heavy leptons as, for example, in seesaw schemes GRS . These contain $`SU(2)U(1)`$ singlets with a gauge invariant Majorana mass term of the type $`M_{R}^{}{}_{ij}{}^{}\nu _i^c\nu _j^c`$ which breaks total lepton number symmetry, perhaps at a large $`SO(10)`$ or left-right breaking scale. The masses of the light neutrinos are obtained by diagonalizing the mass matrix

$$\left[\begin{array}{cc}M_L& D\\ D^T& M_R\end{array}\right]$$
(2)
in the basis $`\nu ,\nu ^c`$, where $`D`$ is the standard $`SU(2)U(1)`$ breaking Dirac mass term, and $`M_R=M_R^T`$ is the large isosinglet Majorana mass and the $`M_L\nu \nu `$ term is an isotriplet theory . In left-right models the latter is generally suppressed as $`M_L1/M_R`$.

The structure of the associated effective weak $`VA`$ currents is rather complex theory . The first point to notice is that the heavy isosinglets will mix with the ordinary isodoublet neutrinos in the charged current weak interaction. As a result, the mixing matrix describing the charged leptonic weak interaction is a rectangular matrix $`K`$ theory  which may be decomposed as

$$K=(K_L,K_H)$$
(3)
where $`K_L`$ and $`K_H`$ are $`3\times 3`$ matrices. The corresponding neutral weak interactions are described by a non-trivial matrix theory

$$P=K^{}K.$$
(4)
In such models non-standard interactions of neutrinos with matter arise from the non-trivial structures of the charged and neutral weak currents. Note, however, that the smallness of neutrino mass, which follows due to the seesaw mechanism $`M_{\nu \mathrm{eff}}=M_LDM_R^1D^T`$ and the condition

$$M_LM_R,$$
(5)
implies that the magnitude of neutrino NSI is expected to be negligible. However this need not be so in general. For example, since the number $`m`$ of $`SU(2)U(1)`$ singlets is arbitrary, one may consider models with Majorana neutrinos based on any value of $`m`$. One can therefore extend the lepton sector of the $`SU(2)U(1)`$ theory by adding a set of two 2-component isosinglet neutral fermions, denoted $`\nu _{}^{c}{}_{i}{}^{}`$ and $`S_i`$, in each generation NSImodels2 . In such $`m=6`$ models one can consider the $`9\times 9`$ mass matrix Gonzalez-Garcia:1989rw

$$\left[\begin{array}{ccc}0& D& 0\\ D^T& 0& M\\ 0& M^T& \mu \end{array}\right]$$
(6)
(in the basis $`\nu ,\nu ^c,S`$). The Majorana masses for the neutrinos are determined from

$$M_L=DM^1\mu M_{}^{T}{}_{}{}^{1}D^T.$$
(7)
In the limit $`\mu 0`$ the exact lepton number symmetry is recovered and neutrinos become massless NSImodels2 . This provides an elegant way to generate neutrino masses without a super-heavy scale and automatically allows one to enhance the magnitude of neutrino NSI strengths by avoiding the constraints which arise from the smallness of neutrino masses presently indicated by the oscillation interpretation of solar and atmospheric neutrino data.

The propagation of the light neutrinos is effectively described by a truncated mixing matrix $`K_L`$ which is not unitary. This may lead to oscillation effects in matter, even if neutrinos were massless Valle:1987gv . They maybe be resonant and therefore important in supernovae matter Valle:1987gv ; NSIastro . The strength of NSI is hence unrestricted by the magnitude of neutrino masses, only by universality limits, and may be large, at the few per cent level. The phenomenological implications of these models have been widely investigated Bernabeu:1987gr ; Gonzalez-Garcia:1990fb ; Rius:1990gk ; Valle:1991pk ; Gonzalez-Garcia:1992be .

An alternative way to induce neutrino NSI is in the context of the most general low-energy super-symmetry model, without R-parity conservation Ross:1985yg . In addition to bilinear violation Hall:1984id ; Hirsch:2000ef  one may also have trilinear $`L`$ violating couplings in the super-potential

$$\lambda _{ijk}L_iL_jE_k^c$$
(8)
$$\lambda _{ijk}^{}L_iQ_jD_k^c$$
(9)
where $`L,Q,E^c`$ and $`D^c`$ are (chiral) super-fields which contain the usual lepton and quark $`SU(2)`$ doublets and singlets, respectively, and $`i,j,k`$ are generation indices. The couplings in Eq. (8) give rise at low energy to the following four-fermion effective Lagrangian for neutrinos interactions with $`d`$-quark including

$$_{\mathrm{eff}}=2\sqrt{2}G_F\underset{\alpha ,\beta }{}\xi _{\alpha \beta }\overline{\nu }_{L\alpha }\gamma ^\mu \nu _{L\beta }\overline{d}_R\gamma ^\mu d_R\alpha ,\beta =e,\mu ,\tau ,$$
(10)
where the parameters $`\xi _{\alpha \beta }`$ represent the strength of the effective interactions normalized to the Fermi constant $`G_F`$. One can identify explicitly, for example, the following non-standard flavor-conserving NSI couplings

$`\xi _{\mu \mu }`$ $`={\displaystyle \underset{j}{}}{\displaystyle \frac{|\lambda _{2j1}^{}|^2}{4\sqrt{2}G_Fm_{\stackrel{~}{q}_{jL}}^2}},`$ (11)
$`\xi _{\tau \tau }`$ $`={\displaystyle \underset{j}{}}{\displaystyle \frac{|\lambda _{3j1}^{}|^2}{4\sqrt{2}G_Fm_{\stackrel{~}{q}_{jL}}^2}},`$ (12)
and the FC coupling

$$\xi _{\mu \tau }=\underset{j}{}\frac{\lambda _{3j1}^{}\lambda _{2j1}^{}}{4\sqrt{2}G_Fm_{\stackrel{~}{q}_{jL}}^2}$$
(13)
where $`m_{\stackrel{~}{q}_{jL}}`$ are the masses of the exchanged squarks and $`j=1,2,3`$ denotes $`\stackrel{~}{d}_L,\stackrel{~}{s}_L,\stackrel{~}{b}_L`$, respectively. The existence of effective neutral current interactions contributing to the neutrino scattering off $`d`$-quarks in matter, provides new flavor-conserving as well as flavor-changing terms for the matter potentials of neutrinos. Such NSI are directly relevant for solar and atmospheric neutrino propagation Guzzo:2001mi .

Clearly, such neutrino NSI are accompanied by non-zero neutrino masses. In fact one has a hybrid model for neutrino masses in which the atmospheric scale arises at the tree level while the solar neutrino scale is induced by loops which involve directly the NSI coefficients in (8). This way one obtains the co-existence of oscillations as well as NSI of neutrinos. The relative importance of NSI and oscillation features is parameter-dependent.

An alternative variant of the above scheme is provided by some radiative models of neutrino masses such as the one discussed in Fukugita:1988qe . In all such models NSI may arise from scalar interactions.

Finally, we mention that unification provides an alternative and elegant way to induce neutrino NSI. For example unified super-symmetric models lead to NSI as a result of super-symmetric scalar lepton non-diagonal vertices induced by renormalization group evolution NSImodels3 ; Barbieri:1995tw . In the special case of $`SU(5)`$ the NSI may exist without neutrino mass. In $`SO(10)`$ neutrino masses co-exist with neutrino NSI.

In what follows we shall investigate the interplay of NSI-induced and neutrino-mass-induced (oscillation-induced) conversion of neutrinos at a neutrino factory and on how it can vitiate the, otherwise very precise, determination of neutrino oscillation parameters.

## III Effective four-fermion Hamiltonians describing NSI

In this section we consider in some detail the simplest examples of effective low energy Hamiltonians, which lead to NSI in the source ($`S`$), propagation ($`P`$) and detection ($`D`$) of neutrinos in a neutrino factory experiment.

In such an experiment neutrinos are produced by the decay of the stored muons $`\mu ^+e^++\nu _e+\overline{\nu }_\mu `$ and the charge conjugated process. In the SM these processes are described by the effective four-fermion Hamiltonian

$$_{\mathrm{SM}}^S=\frac{G_F}{\sqrt{2}}\left[\overline{\nu }_\mu (1\gamma _5)\gamma _\lambda \mu \right]\left[\overline{e}(1\gamma _5)\gamma ^\lambda \nu _e\right]+\mathrm{h}.\mathrm{c}.$$
(14)
In addition to this SM term, resulting from the exchange of the W-boson, we consider now new processes $`\mu ^+e^++\nu _\alpha +\overline{\nu }_\mu `$ with any flavor $`\alpha =e,\mu ,\tau `$ for the neutrino related to the positron.111In this work we will consider the $`e\mu `$ appearance channel and hence, we are interested only in the neutrino produced together with the positron (or the electron, for the charge conjugated processes). More generally, also processes $`\mu ^+e^++\nu _\alpha +\overline{\nu }_\beta `$ with an arbitrary flavor combination $`(\alpha ,\beta )=(e,\mu ,\tau )`$ may occur. In this case the final state in the source can be different from the one in the SM. Such additional processes must be added incoherently to obtain the transition rate defined later in Eq. (26). For simplicity we will not consider this possibility further. We parametrize the corresponding NSI effective Hamiltonian by the coefficients $`ϵ_{e\alpha }^S`$:

$$_{\mathrm{NSI}}^S=\frac{G_F}{\sqrt{2}}\left[\overline{\nu }_\mu (1\gamma _5)\gamma _\lambda \mu \right]\underset{\alpha }{}ϵ_{e\alpha }^S\left[\overline{e}(1\gamma _5)\gamma ^\lambda \nu _\alpha \right]+\mathrm{h}.\mathrm{c}.$$
(15)
Note that the value of the SM Fermi constant $`G_F`$ is determined experimentally from the decay-width of the muon pdg , without measuring the flavor of the neutrinos. Therefore, we have the relation

$$G_F^{\mathrm{exp}}=G_F\left(|1+ϵ_{ee}^S|^2+\underset{\alpha =\mu ,\tau }{}|ϵ_{e\alpha }^S|^2\right)^{1/2}.$$
(16)
The term $`ϵ_{ee}^S`$ leads to exactly the same final state as the SM process and must be added coherently, whereas $`ϵ_{e\mu }^S`$ and $`ϵ_{e\tau }^S`$ lead to different final states and contribute incoherently to the decay. From Eq. (16) one learns that the high precision measurement of $`G_F^{\mathrm{exp}}`$ on its own (within an accuracy of $`9\times 10^6`$ pdg ) does not constrain any of the parameters in the Hamiltonian $`G_F`$ and $`ϵ_{e\alpha }^S`$ directly Langacker:1988cm ; only the combination shown in Eq. (16) is constrained within the accuracy of the experimental measurement.

The standard muon detectors under discussion for a neutrino factory experiment make use of charged current processes like $`\nu _\mu +d\mu ^{}+u`$. The relevant effective Hamiltonian in the SM is given by

$$_{\mathrm{SM}}^D=\frac{G_F}{\sqrt{2}}\left[\overline{d}(1\gamma _5)\gamma _\lambda u\right]\left[\overline{\nu }_\mu (1\gamma _5)\gamma ^\lambda \mu \right]+\mathrm{h}.\mathrm{c}.$$
(17)
Here $`d(u)`$ symbolizes any down-(up-)type quark. Similar to Eq. (15) we consider the following NSI four-fermion Hamiltonian:

$$_{\mathrm{NSI}}^D=\frac{G_F}{\sqrt{2}}\left[\overline{d}(1\gamma _5)\gamma _\lambda u\right]\underset{\alpha }{}ϵ_{\alpha \mu }^D\left[\overline{\nu }_\alpha (1\gamma _5)\gamma ^\lambda \mu \right]+\mathrm{h}.\mathrm{c}.$$
(18)
The coefficients $`ϵ_{\alpha \mu }^D`$ describe NU ($`\alpha =\mu `$) or FC ($`\alpha =e,\tau `$) NSI in the detector, e.g. a non-zero $`ϵ_{\tau \mu }^D`$ leads to the process $`\nu _\tau +d\mu ^{}+u`$.

In a long-baseline neutrino experiment a significant part of the neutrino path will cross the earth and hence neutrino NSI with earth matter must be taken into account. Let us consider the effective Hamiltonian describing the SM neutral current processes of neutrinos with a fermion $`f`$ due to the exchange of the Z-boson $`\nu _\alpha +f\nu _\alpha +f`$:

$$_{\mathrm{SM}}^P=\frac{G_F}{\sqrt{2}}\left[\overline{f}(g_V^fg_A^f\gamma _5)\gamma _\lambda f\right]\underset{\alpha }{}\left[\overline{\nu }_\alpha (1\gamma _5)\gamma ^\lambda \nu _\alpha \right],$$
(19)
where $`g_V^f`$ and $`g_A^f`$ are the SM vector and axial couplings of the fermion $`f`$, see e.g. Ref. pdg  Sec. 10. In the SM this interaction is the same for all flavors and hence, has no effect on the propagation of the neutrino state – in contrast to the charged current interaction of $`\stackrel{()}{\nu }_e`$ with electrons. However, if NSI are present, we must also take into account processes $`\nu _\alpha +f\nu _\beta +f`$ with arbitrary flavor combinations ($`\alpha \beta `$) described by the Hamiltonian

$$_{\mathrm{NSI}}^P=\frac{G_F}{\sqrt{2}}\left[\overline{f}(g_V^fg_A^f\gamma _5)\gamma _\lambda f\right]\underset{\alpha \beta }{}ϵ_{\alpha \beta }^f\left[\overline{\nu }_\alpha (1\gamma _5)\gamma ^\lambda \nu _\beta \right].$$
(20)
In Eqs. (15), (18), (20) we have assumed for simplicity, that the NSI have the same $`VA`$ Lorentz structure as the SM interactions. This needs not to be the case in the most general extension of the SM involving, say, left-right symmetry, where many new NSI parameters can appear (see e.g. Ref. johnson ; ota ). However, the effects of NSI with $`V+A`$ Lorentz structure are strongly suppressed since the left-right breaking scale should be rather high in order to account for the smallness of the neutrino masses indicated by solar and atmospheric experiments. Moreover, one expects that only certain combinations of parameters will be relevant for the experimental configuration we are considering here, and for any given theory our results can be mapped to the corresponding combination of parameters.

Although different processes are relevant for source, propagation and detection, in a given model relations between the coefficients $`ϵ_{\alpha \beta }^X`$ may exist. However, such relations highly depend on the underlying model. In order to be model-independent we will treat all $`ϵ_{\alpha \beta }^X`$ as independent parameters.

## IV Bounds on NSI parameters

In this section we review existing bounds on neutrino NSI obtained in the literature. The most direct upper bounds on the strength of NSI interactions arise from negative searches for neutrino oscillations johnson ; Grossman:1995wx ; Gninenko:2001id . For example the bounds on the transition probabilities $`P_{\nu _\mu \nu _\tau }3.4\times 10^4`$ and $`P_{\nu _e\nu _\tau }2.6\times 10^2`$ obtained by CHORUS chorus  yield the bounds

$$|ϵ_{\mu \tau }^{\mathrm{CHORUS}}|1.8\times 10^2,|ϵ_{e\tau }^{\mathrm{CHORUS}}|0.16.$$
(21)
With the superscript we indicate that the constrained quantity actually is a certain combination of NSI coefficients relevant in the neutrino source and detection processes for a given experiment (see later Eq. (33)), which in general is different from the NSI coefficients relevant for neutrino factory experiments.

Recently, in Ref. Val  the strong evidence for oscillations of atmospheric neutrinos has been used to set upper bounds on NSI of neutrinos with the down quarks in earth matter:

$$|ϵ_{\mu \tau }^d|3\times 10^2,|ϵ_{\tau \tau }^d|6\times 10^2.$$
(22)
Similar bounds on the magnitude of neutrino NSI with electrons and up-type quarks may be derived Guzzo:2001mi .

Besides these direct bounds on neutrino NSI there is a lot of data constraining non-standard effects in charged lepton processes. However, it is very non-trivial in general to use these data to derive model-independent bounds on neutrino NSI coefficients. For a recent discussion see, for example Refs. BergmannAtm ; BergmannLSND ; Gonzalez-Garcia:2001mp . To obtain such bounds for neutrino interactions one must convert from the bounds on charged lepton processes making some assumption about $`SU(2)_L`$ symmetry. In Ref. Bergmann:2000gp  the corresponding bounds for the charged leptons are multiplied by a factor of $`6.8`$, in order to take into account $`SU(2)_L`$ breaking effects. In this way the following bounds on neutrino NSI are derived from pure leptonic processes Bergmann:2000gp :

$$|ϵ_{e\mu }^{\mathrm{}}|7\times 10^6,|ϵ_{e\tau }^{\mathrm{}}|3\times 10^2.$$
(23)
As the neutrino production in a neutrino factory is also a pure leptonic process we take the bounds (23) as order of magnitude estimates of the NSI at the neutrino source. From bounds on $`\mu e`$ conversion in muon scattering off nuclei and from those on flavour-violating hadronic tau decays the following bounds on neutrino NSI with quarks are derived Bergmann:2000gp :

$$|ϵ_{e\mu }^q|7\times 10^5,|ϵ_{e\tau }^q|7\times 10^2.$$
(24)
We take this as an order of magnitude estimate for the NSI in propagation and detection at a neutrino factory experiment, since there also processes with quarks are involved. For the $`\mu \tau `$ channel the bounds are of order BergmannAtm

$$|ϵ_{\mu \tau }|5\times 10^2$$
(25)
and for the NU coefficients upper bounds of order 0.1 are derived.

However, in Ref. rossi  it has been stressed that in general no model-independent relation exists between NSI coefficients for charged leptons and neutrinos and only much weaker bounds of order $`50\%`$ are derived using data from $`e^+e^{}`$ colliders. This more conservative viewpoint has been exploited to show how FC neutrino interactions provide an excellent description of the solar neutrino data, while consistent with the oscillation description of the atmospheric data Guzzo:2001mi .

## V The appearance rate in a neutrino factory experiment

Let us consider the impact of NSI in source, propagation and detection on the $`e\mu `$ channel222For the sake of clarity we keep the oscillation-inspired terminology “$`e\mu `$ appearance channel” (“$`\mu \mu `$ disappearance channel”). We actually mean the production of a wrong-sign (like-sign) muon in the detector, respectively. at a neutrino factory. Starting from the decay of a $`\mu ^+`$, we make the following ansatz for the rate at which a neutrino produced together with the positron leads to the production of a $`\mu ^{}`$ in the detector ota :

$$_{e\mu }=\left|\underset{\alpha \beta }{}𝒜_{e\alpha }^S𝒜_{\alpha \beta }^P𝒜_{\beta \mu }^D\right|^2,$$
(26)
and similar for the charge conjugated process. Here we define the amplitudes describing the neutrino source and detection process as

$$𝒜_{\alpha \beta }^X\delta _{\alpha \beta }+ϵ_{\alpha \beta }^X\text{for}X=S,D,$$
(27)
and the amplitude $`𝒜_{\beta \gamma }^P`$ describes the propagation of the neutrino state from the production point to the detector. This amplitude is obtained from the solution of a Schrödinger equation with the Hamiltonian

$$H_\nu =\frac{1}{2E_\nu }U\text{diag}(0,\mathrm{\Delta }m_{\mathrm{sol}}^2,\mathrm{\Delta }m_{\mathrm{atm}}^2)U^{}+\text{diag}(V,0,0)+\underset{f}{}V_fϵ^f,$$
(28)
which takes into account neutrino oscillations and SM interactions as well as NSI with the matter crossed by the neutrino beam. Here $`E_\nu `$ is the neutrino energy and $`V=\sqrt{2}G_FN_e`$ is the matter potential due to the SM charged current interaction MSW , where $`N_e`$ is the electron number density. The last term in Eq. (28) describes the NSI with earth matter. The sum is over all fermions $`f`$ present in matter, and $`V_fϵ_{\alpha \beta }^f`$ is the coherent forward scattering amplitude of the process $`\nu _\alpha +f\nu _\beta +f`$, where $`V_f=\sqrt{2}G_FN_f`$, with the number density of the fermion $`f`$ along the neutrino path given by $`N_f`$. We define an effective NSI coefficient for the propagation by normalizing all contributions to the down-quark potential $`V_d`$:

$$ϵ_{\alpha \beta }^P\underset{f}{}\frac{V_f}{V_d}ϵ_{\alpha \beta }^f.$$
(29)
Adopting a basis where the charged lepton mass matrix is diagonal, the unitary matrix $`U`$ in Eq. (28) relates the neutrino fields in the basis where the neutrino mass matrix is diagonal to the neutrino fields in the basis where the interaction with the SM W-boson is diagonal Grossman:1995wx . We parameterize this matrix in the following way theory :

$$U=U_{23}U_{13}U_{12}=\left(\begin{array}{ccc}1& 0& 0\\ 0& c_{23}& s_{23}\\ 0& s_{23}& c_{23}\end{array}\right)\left(\begin{array}{ccc}c_{13}& 0& s_{13}\\ 0& 1& 0\\ s_{13}& 0& c_{13}\end{array}\right)\left(\begin{array}{ccc}c_{12}& e^{i\delta }s_{12}& 0\\ e^{i\delta }s_{12}& c_{12}& 0\\ 0& 0& 1\end{array}\right),$$
(30)
where $`s_{ij}=\mathrm{sin}\theta _{ij}`$ and $`c_{ij}=\mathrm{cos}\theta _{ij}`$.

In this paper we consider the following simplified scenario. First, we take all $`ϵ_{\alpha \beta }^X`$ real and we assume that they are the same for neutrinos and anti-neutrinos as in Val . Second, Eqs.(23) and (24) suggest that constraints on FC interactions in the $`e\mu `$ channel are about 3 orders of magnitude stronger than in the other channels. This motivates the approximation333Note that in Ref. Gonzalez-Garcia:2001mp  much weaker sensitivities on FCI in the $`e\mu `$ channel for the source at a neutrino factory are derived, of the order of $`10^3`$. A non-vanishing $`ϵ_{e\mu }^S`$ would imply yet an additional source for the confusion of oscillations and NSI.

$$ϵ_{e\mu }^Xϵ_{\mu e}^X0\text{for}X=S,P,D.$$
(31)
Third, we neglect the solar mass-squared difference, which implies also that the angle $`\theta _{12}`$ and the phase $`\delta `$ disappear theory . Then we are left with the following neutrino propagation Hamiltonian

$$H_\nu =U_{23}U_{13}\text{diag}(0,0,\mathrm{\Delta })U_{13}^{}U_{23}^{}+\text{diag}(V,0,0)+Vrϵ^P,$$
(32)
where we have defined $`\mathrm{\Delta }\mathrm{\Delta }m_{\mathrm{atm}}^2/2E_\nu `$ and $`rV_d/V=N_d/N_e`$ with $`r3`$ in earth matter. In the Hamiltonian (32) a sign change of $`\mathrm{\Delta }`$ is equivalent to a sign change of $`V`$, which interchanges the evolution of neutrinos and anti-neutrinos. Therefore, it is sufficient to consider only the case $`\mathrm{\Delta }>0`$, assuming that the neutrino factory is run in both polarities.

A detector close to the front end of a neutrino factory – a so-called near detector – can be a very powerful tool to constrain NSI Datta:2000ci ; Bueno:2001jy ; Mangano:2001mj . Such a detector has to be situated at a short distance (a few 100 m) from the production region of the neutrinos, such that no oscillations with $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ or $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ can develop and matter effects are negligible. In our formalism this means that $`𝒜_{\alpha \beta }^P=\delta _{\alpha \beta }`$, and the transition rate relevant for a near detector is simply given by

$$_{\alpha \beta }^{\mathrm{ND}}=\left|\underset{\gamma }{}𝒜_{\alpha \gamma }^S𝒜_{\gamma \beta }^D\right|^2.$$
(33)
It is clear that a near detector cannot provide any model-independent information on $`ϵ_{\alpha \beta }^P`$, and only a combination of $`ϵ_{\alpha \beta }^S`$ and $`ϵ_{\alpha \beta }^D`$ is constrained.

Some remarks are in order. Although the general motivation for NSI is that these accompany models of neutrino mass generation, in our following phenomenological studies we will restrict our attention only to total lepton number conserving NSI. While this will suffice to make our point, it will on the other hand greatly simplify our analysis. This happens because in this particular case it is possible to distinguish the neutrino produced together with the electron from the one produced together with the muon, if the detector can determine the charge of the muon. We insist, however, that in a generic theory for NSI and neutrino masses also lepton number violating processes Bergmann:2000gn ; Bueno:2001jy  are expected, due to the Majorana nature of neutrinos (see Sec. II). Such effects would be an additional source for the confusion of NSI and oscillations.

The off-diagonal elements $`ϵ_{\alpha \beta }^X(X=S,P,D)`$ with $`\alpha \beta `$ describe FC, whereas the diagonal elements with $`\alpha =\beta `$ lead to NU. In Eqs. (26) and (33) we consider only processes with the same final states (in source and detector) as in the SM case. If additional processes are present, with different final states, the corresponding amplitudes have to be added incoherently to the rate ota .

We want to stress that our numerical results are not restricted only to the NSI resulting from the four-fermion operators discussed in Sec. III. The results apply to all kinds of non-standard physics in source, propagation and detection in a neutrino factory experiment, which can be parametrized like in Eqs. (26), (27), (32). In general the NSI parameter combinations involved in the quantum mechanical evolution of the neutrino system are model-dependent functions of the parameters appearing in the Lagrangian of a given theory, as discussed in Sec. II.

## VI The oscillation–NSI confusion theorem

In this section we present analytical approximations for the transition rate and we show that within our simplified scenario NSI can lead to exactly the same signal at a neutrino factory as expected from genuine neutrino oscillations due to $`\theta _{13}`$. Taking into account additional parameters like $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$ or CP violating phases in the lepton mixing matrix or in NSI can only bring more serious complications for the determination of $`\theta _{13}`$ Freund:2001ui .

For the understanding of the physics relevant for the numerical results which we will present in the following sections it is useful to derive an analytic expression for the appearance rate Eq. (26). To this aim we assume a constant matter potential $`V`$ and consider the terms containing the NSI parameters $`ϵ_{\alpha \beta }^P`$ and $`s_{13}`$ as a small perturbation of the Hamiltonian and calculate eigenvalues and eigenvectors of Eq. (32) up to first order in these small quantities. Then the appearance rate Eq. (26) is of second order in $`s_{13}`$ and $`ϵ_{\alpha \beta }^X(X=S,P,D)`$ and we make the interesting observation that only the three parameters ($`s_{13},ϵ_{e\tau }^S,ϵ_{e\tau }^P`$) appear.

This is a special feature of the $`e\mu `$ channel under the approximation (31) and can be understood from Fig. 1 (a), where we show schematically the various contributions to this channel. The thick lines indicate the SM processes in source and detection and the dominating $`\mu \tau `$ oscillations due to atmospheric oscillation parameters. The thin lines show the leading contributions of small quantities, which involve only the parameters $`s_{13},ϵ_{e\tau }^S`$ and $`ϵ_{e\tau }^P`$. With very thin lines we show some channels which involve more than one small quantity, and hence do not appear up to second order in the appearance rate: we consider only the detection of muons, therefore in leading order no effects of NSI in the detector show up because of Eq. (31). Similarly no FC effects in the $`\tau \mu `$ channel show up, since transitions from $`e`$ to $`\tau `$ flavor already involve a small quantity, either $`s_{13}`$ or $`ϵ_{e\tau }^{S,P}`$. We also note that no NU coefficient $`ϵ_{\alpha \alpha }^X`$ appears in leading order Gago:2001xg ; Gonzalez-Garcia:2001mp , explaining the lack of sensitivity of neutrino factory experiments to NU parameters Huber:2001zw .

Let us introduce the abbreviations

$$ϵ_Sϵ_{e\tau }^S,ϵ_Pϵ_{e\tau }^P.$$
(34)
Then the expression for the appearance rate is a general quadratic form in the variables $`s_{13},ϵ_S`$ and $`ϵ_P`$:

$$_{e\mu }As_{13}^2+Bs_{13}ϵ_P+Cϵ_P^2+Dϵ_Pϵ_S+Eϵ_S^2+Fs_{13}ϵ_S$$
(35)
with the coefficients

$$\begin{array}{cc}\hfill \text{}A=& 4s_{23}^2\left(\frac{\mathrm{\Delta }}{\mathrm{\Delta }V}\right)^2\mathrm{sin}^2\frac{(\mathrm{\Delta }V)L}{2},\hfill \\ \hfill \text{}B=& 4s_{23}^2c_{23}r\frac{\mathrm{\Delta }}{\mathrm{\Delta }V}\left[\frac{\mathrm{\Delta }+V}{\mathrm{\Delta }V}\mathrm{sin}^2\frac{(\mathrm{\Delta }V)L}{2}+\mathrm{sin}^2\frac{VL}{2}\mathrm{sin}^2\frac{\mathrm{\Delta }L}{2}\right],\hfill \\ \hfill \text{}C=& 4s_{23}^2c_{23}^2r^2\frac{\mathrm{\Delta }}{\mathrm{\Delta }V}\left[\frac{V}{\mathrm{\Delta }V}\mathrm{sin}^2\frac{(\mathrm{\Delta }V)L}{2}+\mathrm{sin}^2\frac{VL}{2}\frac{V}{\mathrm{\Delta }}\mathrm{sin}^2\frac{\mathrm{\Delta }L}{2}\right],\hfill \\ \hfill \text{}D=& 4s_{23}^2c_{23}^2r\frac{\mathrm{\Delta }}{\mathrm{\Delta }V}\left[\mathrm{sin}^2\frac{(\mathrm{\Delta }V)L}{2}\mathrm{sin}^2\frac{VL}{2}\left(12\frac{V}{\mathrm{\Delta }}\right)\mathrm{sin}^2\frac{\mathrm{\Delta }L}{2}\right],\hfill \\ \hfill \text{}E=& 4s_{23}^2c_{23}^2\mathrm{sin}^2\frac{\mathrm{\Delta }L}{2},\hfill \\ \hfill \text{}F=& 4s_{23}^2c_{23}\frac{\mathrm{\Delta }}{\mathrm{\Delta }V}\left[\mathrm{sin}^2\frac{(\mathrm{\Delta }V)L}{2}\mathrm{sin}^2\frac{VL}{2}+\mathrm{sin}^2\frac{\mathrm{\Delta }L}{2}\right],\hfill \end{array}$$
(36)
where $`L`$ is the distance between neutrino source and detector. The appearance rate for anti-neutrinos $`_{\overline{e}\overline{\mu }}`$ is obtained by replacing $`VV`$ in Eq. (36). These analytic expressions are in agreement with numerical calculations within a few % in the relevant parameter range. In general all coefficients ($`A,\mathrm{},F`$) are of the same order of magnitude and depend on neutrino energy (via $`\mathrm{\Delta }`$), on the baseline and on the sign of $`V`$ (neutrinos or anti-neutrinos) in a nontrivial way.

Performing a similar consideration for the $`\mu `$ disappearance channel one finds that the rate $`_{\mu \mu }`$, defined in a way similar to Eq. (26), contains terms of all powers of the small quantities. This is illustrated in Fig. 1 (b). The zeroth order contribution corresponds to $`\mu \tau `$ oscillations with atmospheric neutrino oscillation parameters: in contrast to Fig. 1 (a), in the case of Fig. 1 (b) there is a channel involving only thick lines. However, we find that up to second order in small quantities only the parameters ($`ϵ_{\mu \tau }^S,ϵ_{\mu \mu }^S,ϵ_{\mu \tau }^P,ϵ_{\mu \mu }^P,ϵ_{\tau \tau }^P,ϵ_{\tau \mu }^D,ϵ_{\mu \mu }^D`$) appear. The important observation is that none of the three parameters ($`s_{13},ϵ_{e\tau }^S,ϵ_{e\tau }^P`$) relevant for the $`e\mu `$ channel appears. Therefore, no additional information on these parameters can be obtained by considering the disappearance channel. An analysis of this channel including all the parameters listed above goes beyond the scope of this paper and we will not consider it any further here.

Let us compare the transition rate for pure oscillations ($`ϵ_S=ϵ_P=0`$)

$$_{e\mu }^{\mathrm{osc}}(s_{13})=As_{13}^2$$
(37)
with the transition rate without oscillations ($`s_{13}=0`$) but non-zero NSI coefficients

$$_{e\mu }^{\mathrm{NSI}}(ϵ_S,ϵ_P)=Cϵ_P^2+Dϵ_Pϵ_S+Eϵ_S^2.$$
(38)
With the expressions for $`A,C,D,E`$ given in Eq. (36) it is easy to check that if the relation

$$ϵ_S=rϵ_P$$
(39)
holds, oscillations are indistinguishable from NSI. More precisely, we obtain

$$_{e\mu }^{\mathrm{NSI}}(rϵ_P,ϵ_P)=_{e\mu }^{\mathrm{osc}}\left(s_{13}\right)$$
(40)
with

$$s_{13}^2=r^2ϵ_P^2\frac{1+\mathrm{cos}2\theta _{23}}{2}.$$
(41)
This means that for each value of $`s_{13}`$ there is a pair of NSI parameters ($`ϵ_S,ϵ_P`$) determined by Eqs. (39) and (41) which in our approximation leads to exactly the same signal as oscillations due to $`s_{13}`$. This includes both energy and baseline dependence, for both neutrinos and anti-neutrinos. We call this the “oscillation–NSI confusion theorem”.

Of course, relation (39) represents a fine-tuning of the parameters $`ϵ_P`$ and $`ϵ_S`$. However, as long as this relation cannot be excluded one has to consider this possibility. Moreover, in a realistic experiment with finite errors and statistical uncertainties there will be a region around the point in the ($`ϵ_S,ϵ_P`$) plane corresponding to Eqs. (39) and (41) in which oscillations cannot be distinguished from NSI.

## VII The simulation of a neutrino factory experiment

In our numerical calculations we assume a neutrino factory with an energy of 50 GeV for the stored muons and $`2\times 10^{20}`$ useful muon decays of each polarity per year for a period of 5 years. We consider a magnetized iron calorimeter with a mass of $`40\mathrm{kt}`$. The neutrino detection threshold is set to 4 GeV, the energy resolution of the detector is approximated by a Gaussian resolution function with $`\mathrm{\Delta }E_\nu /E_\nu =10\%`$ and we use 20 bins in neutrino energy. We do not include any backgrounds, efficiencies and errors in the particle identification. The amplitude $`𝒜_{\alpha \beta }^P`$ describing the neutrino propagation is obtained by numerically solving the neutrino evolution equation with the Hamiltonian (32), using the average matter density along each baseline. In Ref. Freund:1999vc  it was shown that this is an excellent approximation as long as the baseline is shorter than approximately $`\mathrm{10\hspace{0.17em}000}\mathrm{km}`$, i.e. as long as it does not cross the core. Then the transition rate Eq. (26) is folded with neutrino flux, cross section and energy resolution function to in order obtain the expected event rates in the detector. For further details see Refs. Freund:2001ui ; FHL .

Our “observables” are the event rates for the appearance channel $`n_\nu ^i(n_{\overline{\nu }}^i)`$ for neutrinos (anti-neutrinos) in each energy bin $`i`$. We fix the atmospheric oscillation parameters444They will be determined at the neutrino factory with high accuracy from the disappearance channel. at their best fit values given in Gonzalez-Garcia:2001sq ; Val  $`\mathrm{\Delta }m_{\mathrm{atm}}^2=3\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta _{23}=1`$. Hence, at a given baseline, the event rates depend on the three parameters $`S_{13}`$, $`ϵ_S`$ and $`ϵ_P`$, where we have introduced the abbreviation $`S_{13}\mathrm{sin}^22\theta _{13}`$.

In order to evaluate the impact of $`ϵ_S`$ and $`ϵ_P`$ on the capability of a neutrino factory to measure $`S_{13}`$ we proceed as follows. To test a given point $`(S_{13}^0,ϵ_S^0,ϵ_P^0)`$ in the parameter space we calculate the event rates $`(n_x^i)^0n_x^i(S_{13}^0,ϵ_S^0,ϵ_P^0)`$ with $`x=\nu ,\overline{\nu }`$ and construct a $`\chi ^2`$ appropriate for a Poisson distribution555We are dealing with a counting experiment with eventually very low counts.:

$`\chi ^2(S_{13},ϵ_S,ϵ_P;`$ $`S_{13}^0,ϵ_S^0,ϵ_P^0)=`$ (42)
$`2{\displaystyle \underset{x=\nu ,\overline{\nu }}{}}{\displaystyle \underset{i}{}}\left[n_x^i(S_{13},ϵ_S,ϵ_P)(n_x^i)^0+(n_x^i)^0\mathrm{ln}{\displaystyle \frac{(n_x^i)^0}{n_x^i(S_{13},ϵ_S,ϵ_P)}}\right].`$
Thus we obtain an allowed region in the three dimensional space $`𝒫`$ spanned by $`(S_{13},ϵ_S,ϵ_P)`$ in the usual way. Note that the minimum of the $`\chi ^2`$ defined in Eq. (42) is zero and occurs at $`(S_{13},ϵ_S,ϵ_P)=(S_{13}^0,ϵ_S^0,ϵ_P^0)`$. Therefore, the allowed region at the CL $`\alpha `$ is given by the set of all points in $`𝒫`$ which fulfill

$$\chi ^2(S_{13},ϵ_S,ϵ_P;S_{13}^0,ϵ_S^0,ϵ_P^0)\mathrm{\Delta }\chi _\alpha ^2,$$
(43)
where $`\mathrm{\Delta }\chi _\alpha ^2`$ is determined by a $`\chi ^2`$-distribution with 3 degrees of freedom. For simplicity we consider only starting values with $`ϵ_S^0=ϵ_P^0=0`$, i.e. pure oscillations.

## VIII Sensitivity limits for $`\mathrm{sin}^22\theta _{13}`$

In Fig. 2 we show the allowed regions in the $`S_{13}0`$ plane for a baseline of $`\mathrm{3\hspace{0.17em}000}\mathrm{km}`$ and three different starting values for $`S_{13}^0`$. In gray we show the lines with the same total event rates (solid for neutrinos, dashed for anti-neutrinos) as in the starting point. The general shape of these lines can immediately be understood from Eq. (35). The small regions delimited by the dark solid lines arise from a global fit procedure including also the information on the spectra with the expected energy resolution as described above. One notices that these confidence regions follow closely the intersection of the lines of constant rates for neutrinos and anti-neutrinos. This means that most information is obtained from simultaneously taking into account neutrino and anti-neutrino rates Huber:2001de . This follows from the fact that the allowed regions extend as long as the lines of constant neutrino and anti-neutrino rate are close to each other, i.e. both rates are similar to the ones in the test point. We conclude that it is important to run the neutrino factory in both polarities. On the other hand we learn that most of the information is contained in the total rates; the spectral information is not very important: our results are rather insensitive to variations of the number of energy bins and of the energy resolution assumed.

Solutions in the $`S_{13}0`$ plane of the type as shown in Fig. 2 always exist, irrespective of the starting value $`S_{13}^0`$. This is a consequence of the confusion theorem we have presented in Sec. VI. However, the magnitude of the required NSI parameters $`ϵ_S`$ and $`ϵ_P`$ strongly depends on the size of $`S_{13}^0`$ as can be seen from Eq. (41) or Fig. 2. Thus, it is only possible to derive a limit on $`S_{13}`$ if there is a limit on $`ϵ_S`$ and/or $`ϵ_P`$; a neutrino factory can only test a certain value of $`S_{13}`$ if the values of $`ϵ_S`$ and/or $`ϵ_P`$, which lead to the same signal, are ruled out by some other measurement. In Fig. 3 we show the attainable sensitivity limit on $`S_{13}`$ as a function of different limits on the NSI parameters (the present estimated NSI limits are indicated by the solid horizontal lines). We define this sensitivity on $`S_{13}`$ in the following way. For a dense grid of starting values $`S_{13}^0`$ in the range $`10^610^1`$ we calculate the 90% CL allowed region in the $`S_{13}0`$ plane as shown in Fig. 2. For each value of $`S_{13}^0`$ we show the minimum value of $`ϵ_P^2`$ (left panel), $`ϵ_S^2`$ (middle panel) or $`ϵ_P^2+ϵ_S^2`$ (right panel) inside this allowed region. The neutrino factory is sensitive to this value of $`S_{13}`$ only if there is a bound on this NSI parameter (combination of parameters), which is smaller than this minimum value.

Any experiment (e.g. like a near detector at a neutrino factory) trying to measure $`ϵ_S`$ and $`ϵ_P`$ will only restrict a certain combination of the NSI parameters in the source and the detector used in this particular experiment (see Eq. (33)). In general it will be very difficult to translate such a result into a bound on $`ϵ_S`$, and even more difficult on $`ϵ_P`$, in a model-independent way. It is however to be expected that the constrained combination depends on the square of the NSI parameters since the transition rate $`ϵ^2`$. Therefore we show the results for the three simple functions $`ϵ_S^2,ϵ_P^2`$ and $`ϵ_S^2+ϵ_P^2`$ as mentioned above. The left hand panel of Fig. 3 shows the sensitivity limit if there is a bound on $`ϵ_P^2`$ and all values of $`ϵ_S`$ are allowed. It seems very difficult to obtain such a bound in a model-independent way, because it is not possible to probe directly the NSI parameters relevant in neutrino propagation. As recently stressed in Ref. rossi  in general it is not possible to use bounds on similar processes involving charged leptons. Moreover, many different processes may contribute to $`ϵ_P`$ (see Eq. (29)). In the middle panel we show the sensitivity limits for a bound on $`ϵ_S^2`$ and all values allowed for $`ϵ_P`$, which is probably the most realistic case because it should be possible to constrain $`ϵ_S`$ with a near detector setup. In the right hand panel we display the optimal situation, if a bound on the combination $`ϵ_S^2+ϵ_P^2`$ is available.

The diagonal solid line in Fig. 3 shows the theoretical bound implied by the oscillation–NSI confusion theorem Eqs. (39) and (41). This bound corresponds to the best possible situation, which can be achieved only for large values of $`S_{13}`$ due to large event numbers. For smaller values of $`S_{13}`$ the realistic bound gets worse because of statistical limitations due to small event numbers. The numerical differences between the three plots in Fig. 3 is due to the different symmetries of the used function of $`ϵ_S`$ and $`ϵ_P`$ with respect to the symmetry of the allowed regions in the $`S_{13}0`$ plane. For small values of $`ϵ_S^2+ϵ_P^2`$ (right panel) the bounds converge to the sensitivity limits obtained without taking into account NSI Freund:2001ui . The range of these limits for the three different baselines is shown as the gray vertical band.

We can understand the behavior of the sensitivity limits in Fig. 3 by considering the allowed regions in the $`S_{13}0`$ plane for different baselines, as shown in Fig. 4. For small baselines the allowed region is roughly a circle. Therefore a bound on an individual $`ϵ_P`$ or $`ϵ_S`$ is useless for a sensitive determination of $`S_{13}`$ (see left and middle panel of Fig. 3). We conclude that if only a baseline of 700 km is available it is mandatory to establish solid bounds on both NSI parameters. However, a bound on $`ϵ_P^2+ϵ_S^2`$ is most suitable for small baselines and in this case $`L=700`$ km can do even better than longer baselines (right panel of Fig. 3). For longer baselines the allowed regions in the $`S_{13}0`$ plane become smaller (see Fig. 4) and hence, also a bound on an individual NSI parameter is useful.

With the solid horizontal lines in Fig. 3 we illustrate the order of magnitude of current bounds on NSI parameters. For $`ϵ_P`$ we show the bound given in Eq. (24), while for the $`ϵ_S`$ case we use the bound given in Eq. (23). It is clearly visible that using even these rather optimistic bounds on the relevant NSI parameters the sensitivity of a neutrino factory is at its best $`S_{13}10^3`$ – compared to $`S_{13}10^5`$ \- in the absence of any NSI. The sensitivity is deteriorated by two orders of magnitude for all baselines. Let us stress again that we are not using the bounds from Eqs. (24) and (23) in our analysis because they are derived under some model-dependent assumptions from non-neutrino processes. The horizontal lines in Fig. 3 should merely “guide-the-eye” in reading the plots, they simply give a rough idea of the order of magnitude of existing bounds.

The geometry of the currently discussed muon storage rings offers the striking possibility to illuminate two detectors at different baselines with neutrinos from one neutrino factory. With this in mind, let us investigate to which extent the sensitivities for $`S_{13}`$ can be improved by combining the information of two baselines. From Fig. 4 we find that although the shape of the allowed regions in the $`S_{13}0`$ plane is very different for different baselines, they all have a common intersection. Moreover, considering the lines in the $`S_{13}0`$ plane which have the same event rates as the starting value $`S_{13}^0`$ (shown as gray solid lines for neutrinos and gray dashed lines for anti-neutrinos) we observe that all iso-rate lines meet in a single point. This is again due to our oscillation–NSI confusion theorem: Eqs. (39) and (41) do not depend on the baseline. Thus even the combination of baselines cannot lift this degeneracy.

In Fig. 5 we show sensitivity limits for $`S_{13}`$ obtained as before, but now we use the sum of the $`\chi ^2`$-functions Eq. (42) for two different baselines. In comparison with Fig. 3 we observe that the theoretical bound can now be achieved somewhat easier. However, in contrast to the case considered in our previous work Huber:2001de , where only NSI with the earth matter are taken into account, in our present more realistic situation including also effects in the neutrino source even a combination of baselines does not resolve the confusion problem. With the current bounds on the relevant NSI parameters the sensitivity is $`S_{13}10^3`$, which coincides with the sensitivity obtained at a single baseline. Again this sensitivity is two orders of magnitude worse than without NSI.

## IX Conclusions

In this paper we have considered the impact of non-standard neutrino interactions on the determination of neutrino mixing parameters at a neutrino factory. In particular we have focused on the so-called “golden channels” for the measurement of $`\theta _{13}`$, namely the $`\stackrel{()}{\nu }_e\stackrel{()}{\nu }_\mu `$ channels. We have extended our previous work Huber:2001de  by taking into account both the effects of neutrino oscillations as well as the effect of NSI at the neutrino source, propagation and detection ota . Within a very good approximation we have explicitly demonstrated how a certain combination of FC interactions in source and propagation can produce exactly the same signal as would be expected from oscillations arising due to $`\theta _{13}`$. This implies that information about $`\theta _{13}`$ can only be obtained if bounds on NSI parameters are available and that all one can achieve at a neutrino factory is a correlated oscillation–NSI study. In view of the current estimates of the bounds on FC interactions, this leads to a drastic loss in sensitivity in $`\theta _{13}`$, at least two orders of magnitude.

In order to improve the situation it is mandatory to obtain better bounds on $`ϵ_S`$ and $`ϵ_P`$ at the $`ϵ10^410^3`$ level, which is several orders of magnitude more stringent than current limits. This unexpected complication should be taken into account in the design of a neutrino factory. On the other hand, a neutrino factory may also offer the possibility to obtain these very stringent limits on the NSI parameters. Using a small near detector ($`L100\mathrm{m}`$) with very good particle identification for taus it would be possible to restrict the $`_{e\tau }`$ down to $`10^810^6`$. Since $`_{e\tau }ϵ^2`$ this translates into a bound for $`ϵ10^410^3`$ Mangano:2001mj . Thus the near site physics program of a neutrino factory is a necessary and very important part of the long baseline program.

## Acknowledgments

This work was supported by Spanish DGICYT under grant PB98-0693, by the European Commission grants HPRN-CT-2000-00148 and HPMT-2000-00124, and by the European Science Foundation *Neutrino Astrophysics Network* grant N. 86.


text:
# 1 Introduction and summary

## 1 Introduction and summary

In recent years a large number of articles have been written about exact sampling (also called perfect sampling) using Markov chains. See  for an overview. The rough idea is as follows. One wishes to sample from the unique stationary distribution $`\pi `$ of an observed irreducible Markov chain. At each transition of the chain, a decision is made whether to continue observing the chain or to stop. When the observation is stopped, a value $`S`$ is output and it is desired that, for all states $`i`$, $`P(S=i|\text{one ever stops observing the chain})=\pi _i`$. The decision about whether to stop at a particular time is made on the basis of the evolution of the chain up through that time, possibly together with some additional randomness independent of the chain.

The goal of our research leading to this paper was to determine whether or not it is possible to carry out interruptible exact sampling for finite-state chains in what Propp and Wilson  call the passive setting. (We will explain in Section 3.2 what is meant by “interruptible” and “the passive setting”.) Our central result is the following:

*Interruptible exact sampling is not possible when one observes only a single trajectory.*

This result remains true even if we assume that the chain is aperiodic and reversible. \[See Remark 6.2(b).\] However, interruptible exact sampling *is* possible for an $`N`$-state chain when one is able to observe, simultaneously, $`N`$ trajectories. Here is a guide to our specific results.

1. (positive:) We provide an algorithm (Algorithm 4.3) which, given an irreducible Markov chain on $`N`$ states as input, produces in (random) finite time an exact sample from the tree distribution, and hence also an exact sample from $`\pi `$. (The tree distribution is defined in Section 3.1.) The algorithm is interruptible, but requires $`N`$ independent synchronized trajectories from the chain. (See Theorem 4.4.)
2. (negative:) There is no algorithm in the passive setting for obtaining an observation from the stationary distribution of an irreducible aperiodic Markov chain on $`N`$ states which uses fewer than $`N`$ independent trajectories from the chain and which is both interruptible and exact. (See Theorem 5.1.)
3. (negative:) There is no algorithm in the passive setting for obtaining an observation from the common stationary distribution of any finite number of independent irreducible aperiodic Markov chains on $`N`$ states (with possibly different transition matrices) which is both interruptible and exact. (See Theorem 6.1.) This remains true even if we assume that all of the chains are reversible. \[See Remark 6.2(a).\]

## 2 Background

In 1992, Asmussen, Glynn, and Thorisson  demonstrated that exact sampling from a Markov chain is possible under certain circumstances. They also proved that it is not possible to obtain an exact sample from an arbitrary Markov chain without some prior knowledge about the chain; in particular, the size of the state space must be known. Although their paper does provide a method for generating exact samples from an $`N`$-state Markov chain when $`N`$ is known, the paper is primarily of a theoretical nature, and the method is complicated and inefficient.

In 1995, Lovász and Winkler  provided a simpler and more efficient algorithm for obtaining an exact sample from an irreducible $`N`$-state Markov chain. Although not mentioned explicitly in their paper, the method described in Section 3 of Lovász and Winkler can in fact be used to obtain an exact sample from the tree distribution of the Markov chain (as defined in Section 3.1). Aldous , Broder , and Propp and Wilson  also describe algorithms for sampling from the tree distribution. Propp and Wilson  discuss and compare these and other methods of sampling from the tree distribution, and from the stationary distribution. Their discussion includes consideration of such issues as whether or not the sampling is exact or interruptible. To our knowledge, the question of whether interruptible exact sampling is possible in the passive case (as described in Section 3.2) has not previously been considered.

## 3 Preliminaries

### 3.1 The tree distribution

Throughout this paper we consider only finite-state irreducible Markov chains. We assume that the number of states, call it $`N`$, is known; in fact, it turns out that we may as well assume (and so we do) that the state space is known to be $`[N]:=\{1,\mathrm{},N\}`$. We denote the transition matrix of such a chain generically by $`𝐏=(p_{ij})`$.

An irreducible Markov chain on $`[N]`$ can be viewed equivalently as a random walk on a connected weighted directed graph $`G`$. The vertex set of $`G`$ is $`[N]`$, and there is an edge from $`i`$ to $`j`$, with weight $`p_{ij}`$, if and only if $`p_{ij}>0`$.

For the moment, let us consider an *undirected* graph $`G`$ with vertex set $`[N]`$. Then a subgraph $`T`$ of $`G`$ is called a *spanning tree* if it contains all $`N`$ vertices and is connected and acyclic. From any spanning tree, we obtain a *directed spanning tree* by assigning a direction to each edge. A directed spanning tree is called an *arborescence* rooted at a given vertex $`r`$ if all edges are directed towards $`r`$.

We define the *weight* $`w(T)`$ of an arborescence $`T`$ with edges $`\{e_l\}`$ as $`w(T):=_{l=1}^{N1}p(e_l)`$, where $`p(e_l):=p_{ij}`$ if $`e_l`$ is directed from $`i`$ to $`j`$. For the remainder of this paper, when we say “tree” we mean an arborescence $`T`$ with $`w(T)>0`$. The *tree distribution* of the Markov chain is the probability distribution on trees obtained by normalizing the weights $`w(T)`$ so as to sum to unity.

The Markov chain tree theorem is the well-known result (see, for example,  or ) that the stationary distribution $`\pi `$ of the chain can be expressed simply in terms of the tree distribution:

$$\pi _i=w_i/w,i[N],$$

where, writing $`𝒯_i`$ for the set of trees rooted at $`i`$ and $`𝒯`$ for $`_{i[N]}𝒯_i`$,

$$w_i:=\underset{T𝒯_i}{}w(T),w:=\underset{i[N]}{}w_i=\underset{T𝒯}{}w(T).$$

In particular, any algorithm for sampling from the tree distribution provides a means of sampling from $`\pi `$: simply output the root of the tree.

### 3.2 The passive case; interruptible exact sampling

Propp and Wilson  distinguish between the active setting and the passive setting for sampling using a Markov chain. In the active setting, an algorithm is assumed to have access at all times to a transition generator, that is, to a routine which, given any input state $`i`$, generates an observation $`j`$ from the probability distribution $`(p_{ij}:j[N])`$, independent of all previously generated observations. In particular, a user can generate a trajectory from $`𝐏`$ with any desired initial state. In the passive setting, the algorithm has no control over the initial state and can only watch passively as the chain transitions from one state to the next.

We now explain what is meant by an (on-line, Markov-chain-based) interruptible exact sampling algorithm in the passive case; for simplicity, we will do this explicitly only in the case that a single trajectory from the chain is available and the desired output is an observation from the stationary distribution $`\pi `$ (rather than one from the tree distribution). Informally, an exact sampling algorithm must take as input a trajectory from the given Markov chain; possibly using external randomization to make its decisions, it watches the chain only until some finite time and then returns an observation distributed according to $`\pi `$. (Important note: The state returned is not necessarily the state of the chain at the stopping time.) More formally, we can define an exact sampling algorithm as a collection of functions $`\varphi _{k,i}:[N]^{k+1}[0,1]`$ \[with $`\varphi _{k,i}(x_0,\mathrm{},x_k)`$ to be interpreted informally as the conditional probability that the algorithm stops by time $`k`$ and outputs $`i`$, given that it sees the trajectory $`(x_0,\mathrm{},x_k)`$ through time $`k`$\] having the following properties, where (iii) and (iv) must hold for all $`\pi `$, for all $`\rho `$, and for all irreducible transition matrices $`𝐏=(p_{ij})`$ on $`[N]`$ with stationary distribution $`\pi `$:

1. $`k0(x_0,\mathrm{},x_k)[N]^{k+1}:_j\varphi _{k,j}(x_0,\mathrm{},x_k)1`$;
2. $`i[N]k0(x_0,x_1,\mathrm{})[N]^{\mathrm{}}:\varphi _{k,i}(x_0,\mathrm{},x_k)\text{ as }k`$;
3. $`lim_k\mathrm{}_{j[N]}_{x_0,x_1,\mathrm{},x_k}\rho _{x_0}p_{x_0,x_1}\mathrm{}p_{x_{k1},x_k}\varphi _{k,j}(x_0,x_1,\mathrm{},x_k)>0`$;
4. $`i[N]:lim_k\mathrm{}\frac{_{x_0,x_1,\mathrm{},x_k}\rho _{x_0}p_{x_0,x_1}\mathrm{}p_{x_{k1},x_k}\varphi _{k,i}(x_0,x_1,\mathrm{},x_k)}{_{j[N]}_{x_0,x_1,\mathrm{},x_k}\rho _{x_0}p_{x_0,x_1}\mathrm{}p_{x_{k1},x_k}\varphi _{k,j}(x_0,x_1,\mathrm{},x_k)}=\pi _i`$.

In terms of the chain $`X`$ observed and the stopping time $`\tau `$ and output state $`S`$ for the algorithm, the properties can be interpreted informally as (i) $`P(\tau k|X)1`$; (ii) $`P(\tau k|X)`$ as $`k`$; (iii) $`P(\tau <\mathrm{})>0`$; and (iv) $`P(S=i|\tau <\mathrm{})\pi _i`$. When the strengthening

1. $`lim_k\mathrm{}_{j[N]}_{x_0,x_1,\mathrm{},x_k}\rho _{x_0}p_{x_0,x_1}\mathrm{}p_{x_{k1},x_k}\varphi _{k,j}(x_0,x_1,\mathrm{},x_k)=1`$

\[interpreted as $`P(\tau <\mathrm{})=1`$\] of (iii) holds, we will call the algorithm *terminating*. When (iv) can be strengthened to

1. $`i[N]k0`$:

$`{\displaystyle \underset{x_0,x_1,\mathrm{},x_k}{}}\rho _{x_0}p_{x_0,x_1}\mathrm{}p_{x_{k1},x_k}\varphi _{k,i}(x_0,x_1,\mathrm{},x_k)`$
$`=`$ $`\pi _i\times {\displaystyle \underset{j[N]}{}}{\displaystyle \underset{x_0,x_1,\mathrm{},x_k}{}}\rho _{x_0}p_{x_0,x_1}\mathrm{}p_{x_{k1},x_k}\varphi _{k,j}(x_0,x_1,\mathrm{},x_k)`$

\[interpreted as the independence $`P(\tau k,S=i)P(\tau k)\pi _i`$ of $`\tau `$ and $`S\pi `$\], we say that the algorithm is *interruptible*. An interruptible algorithm can be aborted without biasing output; see the discussion in . For active-case algorithms, the leading example of a non-interruptible algorithm is coupling from the past , while interruptible algorithms include cycle popping , Fill’s rejection-based algorithm  , and the Randomness Recycler . The results of this paper, both positive and negative, are for interruptible algorithms.

## 4 A terminating algorithm for interruptible exact sampling in the passive case

In this section we present a terminating algorithm for interruptible exact stationary sampling in the passive case, assuming that one can watch $`N`$ synchronized copies $`X_i=(X_i(t):t=0,1,\mathrm{})`$, $`i[N]`$, of a Markov chain with state space $`[N]`$ and irreducible transition matrix $`𝐏`$. We allow arbitrary initial distribution $`\rho `$ for the $`N`$-variate chain $`𝐗:=(X_1,\mathrm{},X_N)`$, but we assume that $`X_1,\mathrm{},X_N`$ are conditionally independent given the initial state $`(X_1(0),\mathrm{},X_N(0))`$. The algorithm will produce an observation from the tree distribution corresponding to $`𝐏`$ (recall Section 3.1).

### 4.1 The algorithm in a restricted setting

In this subsection we present a terminating algorithm for interruptible exact tree-sampling in the passive case that works under the following additional restriction on $`𝐏`$:

$$\text{Assumption A: }p_{1j}>0\text{ for all }j[N].$$

While this assumption may seem unreasonably restrictive, we will show in Section 4.2 how a simple modification of the algorithm can handle the more general case.

To describe the algorithm we first define the following events for even $`t2`$:

$`A_t`$ $`:=`$ $`_{l[N]}\{X_l(t2)=1\},`$
$`B_t`$ $`:=`$ $`\{X_1(t1)=1\},`$
$`C_t`$ $`:=`$ $`\{\{X_1(t),X_2(t1),\mathrm{},X_N(t1)\}=[N]\},`$
$`D_t(T)`$ $`:=`$ $`\{\text{the graph with directed edges from }X_l(t1)\text{ to }X_l(t)\text{}2lN\text{,}`$
$`\text{is the arborescence }T\},`$
$`D_t`$ $`:=`$ $`_{T𝒯}D_t(T),`$
$`E_t(T)`$ $`:=`$ $`A_tB_tC_tD_t(T),`$
$`E_t`$ $`:=`$ $`_{T𝒯}E_t(T)=A_tB_tC_tD_t.`$
###### Algorithm 4.1 (Terminating interruptible tree-sampling, under Assumption A).

For even $`t2`$, let $`E_t(T)`$ and $`E_t`$ be defined as above, and let $`E_0:=\mathrm{}`$. The algorithm is:

* $`t0`$
* repeat
* $`tt+2`$
* until $`E_t`$ holds
* $`ST`$, for the unique $`T𝒯`$ such that $`E_t(T)`$ holds
* return $`S`$

###### Theorem 4.2.

When Asumption A holds, Algorithm 4.1 is a terminating algorithm for interruptible exact tree-sampling.

###### Proof.

Let $`\tau `$ denote the supremum of the values of the variable $`t`$ during the operation of Algorithm 4.1. Now fix a candidate value $`t`$ of $`\tau `$. Let $`T𝒯`$ be an arborescence, say with edges $`e_l`$ directed from $`i_l`$ to $`j_l`$, which we choose to index (in some arbitrary but fixed order) by $`l\{2,\mathrm{},N\}`$. The event $`D_t(T)`$ is a disjoint union of $`(N1)!`$ subevents, with each subevent corresponding to a way of mapping the $`N1`$ transitions $`(X_l(t1),X_l(t))`$ to the $`N1`$ edges $`e_l`$. These subevents will all enter symmetrically into the calculation below of $`P(\tau =t,S=T)`$. One such subevent is

$$D_t^{}(T):=_{l=2}^N\{(X_l(t1),X_l(t))=(i_l,j_l)\}.$$

Let $`\{i_1\}`$ denote the singleton $`[N]\{i_2,\mathrm{},i_N\}`$.

Define

$$A_t^{}:=\{\tau t2\}A_t,\text{even }t2.$$

Then, using the Markov property and independence of the trajectories,

$`P(\tau =t,S=T)`$
$`=`$ $`P(\{\tau t2\}E_t(T))=P(A_t^{}B_tC_tD_t(T))`$
$`=`$ $`(N1)!P(A_t^{}B_tC_tD_t^{}(T))`$
$`=`$ $`(N1)!P(A_t^{})P(B_t|X_1(t2)=1)\left[{\displaystyle \underset{l=2}{\overset{N}{}}}P(X_l(t1)=i_l|X_l(t2)=1)\right]`$
$`\times P(X_1(t)=i_1|B_t)\left({\displaystyle \underset{l=2}{\overset{N}{}}}p_{i_l,j_l}\right)`$
$`=`$ $`(N1)!P(A_t^{})p_{11}\left[{\displaystyle \underset{l=2}{\overset{N}{}}}p_{1,i_l}\right]p_{1,i_1}w(T)=(N1)!P(A_t^{})p_{11}\left({\displaystyle \underset{l=1}{\overset{N}{}}}p_{1l}\right)w(T).`$
Summing over $`T𝒯`$ we find

$$P(\tau =t)=(N1)!P(A_t^{})p_{11}\left(\underset{l=1}{\overset{N}{}}p_{1l}\right)w$$

and therefore

$$P(\tau =t,S=T)=P(\tau =t)\frac{w(T)}{w},$$

which shows that Algorithm 4.1 is an interruptible exact tree-sampling algorithm. Using the fact that $`𝐗`$ visits $`(1,\mathrm{},1)`$ at even times infinitely often (a.s.) together with the strong Markov property of $`𝐗`$, it is clear that termination occurs at the first success in an almost surely infinite sequence of Bernoulli trials with success probability $`p_{11}\left(_{l=1}^Np_{1l}\right)w>0`$ (note that this is where Assumption A is used). Thus $`P(\tau <\mathrm{})=1`$, that is, Algorithm 4.1 is terminating. ∎

### 4.2 The algorithm in the general setting

To avoid needing Assumption A, we can use the averaging technique of Lovász and Winkler . Let $`𝐏^k`$ be the $`k`$-step transition matrix of the chain $`𝐗`$. Then $`\overline{𝐏}:=\frac{1}{N}_{k=1}^N𝐏^k`$ is an irreducible transition matrix with all entries positive. Moreover, we can effectively use the original chain to sample from this “averaged” chain. The resulting more general algorithm (Algorithm 4.3) obtains, interruptibly, an exact sample $`T`$ from the tree distribution of $`𝐏`$.

To describe Algorithm 4.3, which works in the general setting described at the outset of Section 4, for $`t2N`$ we define the following events to be used in the context of the algorithm:

$`A_t`$ $`:=`$ $`_{l[N]}\{X_l(t2N)=1\},`$
$`B_t`$ $`:=`$ $`\{X_1(t2N+U_0)=1\},`$
$`C_t`$ $`:=`$ $`\{\{X_1(t2N+U_0+U_1),X_2(t2N+U_2),\mathrm{},X_N(t2N+U_N)\}=[N]\},`$
$`D_t(T)`$ $`:=`$ $`\{\text{the graph with directed edges from }X_l(t2N+U_l)`$
$`\text{to }X_l(t2N+U_l+1)\text{}2lN\text{, is the arborescence }T\},`$
$`D_t`$ $`:=`$ $`_{T𝒯}D_t(T),`$
$`E_t(T)`$ $`:=`$ $`A_tB_tC_tD_t(T),`$
$`E_t`$ $`:=`$ $`_{T𝒯}E_t(T)=A_tB_tC_tD_t.`$
In the following algorithm, successive calls to Random()are assumed to generate independent random numbers, each uniformly distributed over $`[N]`$.

###### Algorithm 4.3 (Terminating interruptible stationary sampling).

For $`t2N`$, let $`E_t(T)`$ and $`E_t`$ be defined as directly above, and let $`E_0:=\mathrm{}`$. The algorithm is:

* $`t0`$
* repeat
* $`tt+2N`$
* for $`i0\mathrm{𝚝𝚘}N`$
* $`U_i\text{Random()}`$
* until $`E_t`$ holds
* $`ST`$, for the unique $`T𝒯`$ such that $`E_t(T)`$ holds
* return $`S`$

By modifying slightly the proof of Theorem 4.2, we obtain the following result.

###### Theorem 4.4.

Algorithm 4.3 is a terminating algorithm for interruptible exact tree-sampling. ∎

###### Remark 4.5.

Our interest in providing Algorithm 4.3 is more of a theoretical nature (to establish the possibility of terminating interruptible exact sampling, given enough copies of a chain) than of a practical nature (to provide an efficient algorithm). Thus we have not fine-tuned Algorithm 4.3 to improve its performance, and we will not analyze its running time here.

###### Remark 4.6.

If we make no assumption regarding the independence of the trajectories, then interruptible sampling becomes impossible for $`N2`$ states, no matter how many trajectories are available. Indeed, it is then possible that we are in the extreme case that all the trajectories are identical, i.e., that there is “really” only one trajectory, in which case Theorem 5.1 applies.

## 5 Impossibility of interruptible exact sampling (I)

Algorithm 4.3 requires $`N`$ independent synchronized Markov chain trajectories. This may seem excessive, especially since for interesting chains $`N`$ is often enormously large. But our next main result, Theorem 5.1, shows that this is best possible. Note that to prove Theorem 5.1, we need only show that interruptible exact sampling is impossible using $`N1`$ independent trajectories. Indeed, if interruptible exact sampling is possible with $`m`$ independent trajectories, then for any $`m^{}m`$ it is possible with $`m^{}`$ independent trajectories, since extra trajectories can always be ignored.

###### Theorem 5.1.

There is no algorithm in the passive setting for obtaining an observation from the stationary distribution of an irreducible aperiodic Markov chain on $`N`$ states which uses fewer than $`N`$ independent trajectories from the chain and which is both interruptible and exact.

###### Proof.

We first establish an equation \[(5.3)\] that must hold if there exists an interruptible exact sampling algorithm for $`N`$-state chains (for given $`N2`$) that uses only a *single* trajectory; in that case the discussion of Section 3.2 applies verbatim. A similar equation, namely (5.4), must hold if interruptible exact sampling is possible using $`N1`$ trajectories. But (5.4) will lead to a contradiction via a transition-balancing argument.

So we begin with the case of a single trajectory. Suppose that functions $`\varphi _{k,i}`$ satisfying (i)–(iii) and (iv) of Section 3.2 exist. We remind the reader that (iii) and (iv) were required to hold for all initial distributions $`\rho `$; throughout the present proof it will suffice to consider trajectories starting deterministically at $`1`$. Taking $`\rho `$ to be unit mass $`\delta _1`$ at $`1`$ and $`p_{ij}`$ to be identically $`1/N`$, we find from (iii) that, for some $`0k<\mathrm{}`$,

$$\underset{j[N]}{}\underset{x_1,\mathrm{},x_k}{}\varphi _{k,j}(1,x_1,\mathrm{},x_k)>0$$
(5.1)
Let $`k_0`$ be the minimum such $`k`$, and define $`\varphi _j(𝐱)\varphi _j(x_1,\mathrm{},x_{k_0}):=\varphi _{k_0,j}(1,x_1,\mathrm{},x_{k_0})`$ and $`𝒳_j:=\{𝐱=(x_1,\mathrm{},x_{k_0}):\varphi _j(𝐱)>0\}`$ for $`j[N]`$. Again taking $`\rho `$ to be $`\delta _1`$ and $`p_{ij}`$ to be identically $`1/N`$, we find from (iv) and (5.1) that $`𝒳_i\mathrm{}`$ for $`i[N]`$. Using (iv) again, we find that for any transition matrix $`𝐏`$ with positive entries and stationary distribution $`\pi `$,

$$i[N]:\underset{𝐱𝒳_i}{}\varphi _i(𝐱)p_{1,x_1}\mathrm{}p_{x_{k_01},x_{k_0}}=\pi _i\times \underset{j[N]}{}\underset{𝐱𝒳_j}{}\varphi _j(𝐱)p_{1,x_1}\mathrm{}p_{x_{k_01},x_{k_0}},$$
(5.2)
and all terms on both sides of (5.2) are positive. Recalling the notation of Section 3.1, it now follows in particular that

$$w_2\underset{𝐱𝒳_1}{}\varphi _1(𝐱)\underset{i,j}{}p_{ij}^{n_{ij}(𝐱)}=w_1\underset{𝐱𝒳_2}{}\varphi _2(𝐱)\underset{i,j}{}p_{ij}^{n_{ij}(𝐱)},$$
(5.3)
where we write $`n_{ij}(𝐱)`$ for the number of $`ij`$ transitions in the trajectory $`(1,x_1,\mathrm{},x_{k_0})`$ and again all terms on both sides of the equation are positive.

By the same reasoning, if there exists an interruptible exact sampling algorithm for $`N`$-state chains that uses $`N1`$ independent trajectories, then there exist integer $`k0`$ and nonempty sets $`𝒳_1`$ and $`𝒳_2`$ of $`(N1)`$-tuples

$$𝐱=(x_1(1),\mathrm{},x_1(k);x_2(1),\mathrm{},x_2(k);\mathrm{};x_{N1}(1),\mathrm{},x_{N1}(k))$$

of $`k`$-tuples from $`[N]`$ such that, for any transition matrix $`𝐏`$ with positive entries,

$$w_2\underset{𝐱𝒳_1}{}\varphi _1(𝐱)\underset{i,j}{}p_{ij}^{n_{ij}(𝐱)}=w_1\underset{𝐱𝒳_2}{}\varphi _2(𝐱)\underset{i,j}{}p_{ij}^{n_{ij}(𝐱)},$$
(5.4)
where, for $`l=1,2`$ and and every $`𝐱𝒳_l`$, we have $`\varphi _l(𝐱)>0`$, and where $`n_{ij}(𝐱)`$ is the sum over $`1mN1`$ of the numbers of $`ij`$ transitions within the trajectories $`(1,x_m(1),\mathrm{},x_m(k))`$. To complete the proof, we will show that (5.4) cannot possibly hold. We will make key use of the observation that, for any $`𝐱𝒳_1𝒳_2`$,

$$0n_{1+}(𝐱)n_{+1}(𝐱)N1,$$
(5.5)
where we have introduced the notation

$$n_{1+}(𝐱):=\underset{j=2}{\overset{N}{}}n_{1j}(𝐱),n_{+1}(𝐱):=\underset{i=2}{\overset{N}{}}n_{i1}(𝐱)$$
(5.6)
for the total numbers of transitions out of and into state $`1`$, respectively. Indeed, since each trajectory $`(1,x_m(1),\mathrm{},x_m(k))`$ starts in state $`1`$, the number of transitions out of state $`1`$ within such a trajectory either equals or exceeds by one the number of transitions into state $`1`$.

To obtain the desired contradiction, we begin by observing that (5.4) can be written in the form (eliminating the diagonal variables $`p_{ii}`$) that

$$w_2f_1=w_1f_2$$
(5.7)
for all $`(p_{ij}>0:1ijN)`$ such that $`_{j:ji}p_{ij}<1`$ for every $`i[N]`$, where

$$f_l:=\underset{𝐱𝒳_l}{}\varphi _l(𝐱)\left[\underset{i,j:ij}{}p_{ij}^{n_{ij}(𝐱)}\right]\left[\underset{i}{}\left(1\underset{j:ji}{}p_{ij}\right)^{n_{ii}(𝐱)}\right],l=1,2.$$
(5.8)
Using continuity it follows that (5.7) holds for all $`(p_{ij}0:1ijN)`$ such that $`_{j:ji}p_{ij}1`$ for every $`i[N]`$.

For $`l=1,2`$, note that $`f_l`$ and $`w_l`$ are both polynomial expressions in the variables $`p_{ij}`$, $`1ijN`$ (we will denote this entire collection of $`N(N1)`$ variables by $`𝐩`$); in fact, $`w_1`$ is a polynomial expression in the $`(N1)^2`$ variables $`p_{ij}`$ with $`i,j[N]`$ and $`i\{1,j\}`$ (with a similar reduction in number of variables possible for $`w_2`$). Applying Proposition A.1 (see the Appendix) to $`F:=w_2f_1w_1f_2`$, we conclude that (5.7) holds as an equality in the ring of polynomials in the variables $`𝐩`$ over the complex field. Henceforth we shall write $`G_1G_2`$ to indicate such an identity of polynomials $`G_1,G_2`$.

According to Lemma A.2 in the Appendix, the polynomial $`w_1`$ (again, over the complex field) is irreducible; likewise, so is $`w_2`$. From the polynomial identity $`w_2f_1w_1f_2`$ at (5.7) it then follows that we can write

$$f_lw_lf,l=1,2,$$
(5.9)
for some polynomial $`f`$ in $`𝐩`$. Of course, the polynomial identities (5.9) remain true as we now reduce the number of variables to three by setting $`p_{1j}`$ to $`\alpha `$ for $`j1`$, $`p_{i1}`$ to $`\beta `$ for $`i1`$, and $`p_{ij}`$ to $`\gamma `$ if $`i1`$, $`j1`$, and $`ij`$. Observe that now

$`f_l(\alpha ,\beta ,\gamma )`$ $``$ $`{\displaystyle \underset{𝐱𝒳_l}{}}\varphi _l(𝐱)\alpha ^{n_{1+}(𝐱)}\beta ^{n_{+1}(𝐱)}\gamma ^{n_{++}(𝐱)}`$ (5.10)
$`\times [1(N1)\alpha ]^{n_{11}(𝐱)}[1\beta (N2)\gamma ]^{H(𝐱)},l=1,2,`$
recalling (5.6) and defining

$$n_{++}(𝐱):=\underset{i,j\{2,\mathrm{},N\}:ij}{}n_{ij}(𝐱),H(𝐱):=\underset{i=2}{\overset{N}{}}n_{ii}(𝐱).$$
(5.11)
Also now, by a simple generalization of the bijection argument (, Section 2.3.4.4, p. 390) showing that the number of arborescences rooted at $`1`$ is $`N^{N2}`$,

$$w_1(\beta ,\gamma )\beta [\beta +(N1)\gamma ]^{N2};$$
(5.12)
and

$$w_2(\alpha ,\beta ,\gamma )\alpha v_2(\beta ,\gamma )$$
(5.13)
for some polynomial $`v_2(\beta ,\gamma )`$ which is not divisible by $`\beta `$ \[the explanation for (5.13) being that any $`T𝒯_2`$ has precisely one directed edge leaving vertex $`1`$ and that there exists $`T𝒯_2`$ for which $`1`$ is a leaf. In fact, it can be shown that $`v_2(\beta ,\gamma )[\beta +(N1)\gamma ]^{N2}`$, but we won’t need this.\]

The idea for the remainder of the proof is to derive from the identities (5.9)–(5.10) a polynomial identity in the single variable $`\beta `$, namely (5.14), and then show that (5.14) leads to a contradiction. We will produce (5.14) by eliminating (using suitable divisibility arguments) first $`\alpha `$ and then $`\gamma `$. These arguments are carried out in the next two lemmas. ∎

###### Lemma 5.2.

Suppose that there exists an interruptible exact algorithm in the passive setting for sampling from the stationary distribution of an irreducible aperiodic Markov chain on $`N`$ states which uses fewer than $`N`$ independent trajectories from the chain. Then there exist nonempty sets $`𝒳_1^{}`$ and $`𝒳_1^{\prime \prime }`$ and a polynomial $`r`$ such that

$$\underset{𝐱𝒳_1^{\prime \prime }}{}\varphi _1(𝐱)\beta ^{n_{+1}(𝐱)m_1(\beta )}(1\beta )^{H(𝐱)}\beta ^{N2}r(\beta ),$$
(5.14)
where

$$m_1(\beta )=\underset{𝐱𝒳_1^{}}{\mathrm{min}}n_{+1}(𝐱).$$

###### Proof.

Let $`m_l(\alpha )`$ denote the highest power of $`\alpha `$ that divides $`f_l`$ at (5.10) and define $`\stackrel{~}{m}_l(\alpha ):=\mathrm{min}_{𝐱𝒳_l}n_{1+}(𝐱)`$. We claim that $`m_l(\alpha )=\stackrel{~}{m}_l(\alpha )`$, and note that this sort of highest-power observation will be used frequently—and without accompanying proof—in the sequel. \[Indeed, $`m_l(\alpha )\stackrel{~}{m}_l(\alpha )`$ is clear. To see the reverse inequality, divide $`f_l`$ by $`\alpha ^{\stackrel{~}{m}_l(\alpha )}`$ and set $`\alpha `$ to $`0`$ to obtain the expression

$$\underset{𝐱𝒳_l:n_{1+}(𝐱)=\stackrel{~}{m}_l(\alpha )}{}\varphi _l(𝐱)\beta ^{n_{+1}(𝐱)}\gamma ^{n_{++}(𝐱)}[1\beta (N2)\gamma ]^{H(𝐱)}=:g_l(\beta ,\gamma ),$$
(5.15)
which is not the zero polynomial since it has a positive value when $`\beta =1/N=\gamma `$.\]

By (5.9), (5.12), and (5.13),

$`m_2(\alpha )=m_1(\alpha )+1\text{ and}`$ (5.16)
$`g_lv_lg,l=1,2,`$ (5.17)
where $`g_l`$ is the polynomial defined at (5.15) \[recalling $`\stackrel{~}{m}_l(\alpha )=m_l(\alpha )`$\], $`v_1:=w_1`$, $`v_2`$ is defined at (5.13), and $`g`$ is obtained from $`f`$ by dividing by $`\alpha ^{m_1(\alpha )}`$ and then setting $`\alpha =0`$.

Define $`𝒳_l^{}:=\{𝐱𝒳_l:n_{1+}(𝐱)=m_l(\alpha )\}\mathrm{}`$ for $`l=1,2`$. Then, similarly, the highest power $`m_l(\beta )`$ of $`\beta `$ dividing $`g_l`$ is $`\mathrm{min}_{𝐱𝒳_l^{}}n_{+1}(𝐱)`$;

$$m_2(\beta )=m_1(\beta )1;$$
(5.18)
and, with

$$h_1(\beta ,\gamma ):=\underset{𝐱𝒳_1^{}}{}\varphi _1(𝐱)\beta ^{n_{+1}(𝐱)m_1(\beta )}\gamma ^{n_{++}(𝐱)}[1\beta (N2)\gamma ]^{H(𝐱)},$$

we have

$$h_1(\beta ,\gamma )[\beta +(N1)\gamma ]^{N2}h(\beta ,\gamma )$$
(5.19)
for some polynomial $`h`$.

The highest power $`m_1(\gamma )`$ of $`\gamma `$ dividing $`h_1`$ is $`\mathrm{min}_{𝐱𝒳_1^{}}n_{++}(𝐱)`$. Divide both sides of (5.19) by $`\gamma ^{m_1(\gamma )}`$ and set $`\gamma `$ to $`0`$ to find that (5.14) holds for some polynomial $`r`$, where $`𝒳_1^{\prime \prime }:=\{𝐱𝒳_1^{}:n_{++}(𝐱)=m_1(\gamma )\}\mathrm{}`$. ∎

###### Lemma 5.3.

The identity (5.14) cannot hold.

###### Proof.

It follows from (5.14) that $`n_{+1}(𝐱)m_1(\beta )+N2`$ for all $`𝐱𝒳_1^{\prime \prime }`$. But then, for any such $`𝐱`$ and some $`𝐱^{}𝒳_2^{}`$,

$`n_{+1}(𝐱)`$ $``$ $`m_1(\beta )+N2`$
$`=`$ $`m_2(\beta )+N1\text{ by (}\text{5.18}\text{)}`$
$`=`$ $`n_{+1}(𝐱^{})+N1`$
$``$ $`n_{1+}(𝐱^{})\text{ by the second inequality in (}\text{5.5}\text{)}`$
$`=`$ $`m_2(\alpha )`$
$`=`$ $`m_1(\alpha )+1\text{ by (}\text{5.16}\text{)}`$
$`=`$ $`n_{1+}(𝐱)+1,`$
contradicting the first inequality in (5.5). ∎

## 6 Impossibility of interruptible exact sampling (II)

Algorithm 4.3 succeeds in using $`N`$ independent synchronized Markov chain trajectories to carry out interruptible exact sampling. But the algorithm assumes that each of the trajectories has not only (i) the same stationary distribution, but also (ii) the same transition matrix. In this section we show (Theorem 6.1) that interruptible exact sampling becomes impossible when assumption (ii) is dropped, no matter how (finitely) many trajectories are available.

###### Theorem 6.1.

There is no interruptible algorithm in the passive setting for obtaining an observation exactly from the common stationary distribution of any finite number of independent irreducible aperiodic Markov chains on $`N`$ states.

###### Proof.

Let $`M`$ denote the number of trajectories available. We first prove the impossibility of interruptible exact sampling when $`M=N=2`$, then more generally when $`N=2`$ (regardless of $`M`$), and finally for general $`N`$.

For $`M=N=2`$, we note that if

$$0<p_{12},p_{21}<1\text{ and }0<\rho <1/\mathrm{max}\{p_{12},p_{21}\},$$
(6.1)
then

$$𝐏:=\left(\begin{array}{cc}1p_{12}& p_{12}\\ p_{21}& 1p_{21}\end{array}\right)\text{ and }𝐐:=\left(\begin{array}{cc}1\rho p_{12}& \rho p_{12}\\ \rho p_{21}& 1\rho p_{21}\end{array}\right)$$

are irreducible aperiodic transition matrices with common stationary distribution

$$\pi =(\frac{p_{21}}{p_{21}+p_{12}},\frac{p_{12}}{p_{21}+p_{12}}).$$

Arguing as in the proof of Theorem 5.1, if there exists an interruptible exact sampling algorithm in the present setting, then there exist integer $`k0`$, nonempty sets $`𝒵_1`$ and $`𝒵_2`$ of pairs

$$𝐳=(𝐱,𝐲)=(x_1,\mathrm{},x_k;y_1,\mathrm{},y_k)$$

of $`k`$-tuples from $`\{1,2\}`$, and positive numbers $`\psi _l(𝐳)`$ ($`𝐳𝒵_l`$, $`l=1,2`$) such that, whenever (6.1) holds,

$$p_{12}f_1=p_{21}f_2$$
(6.2)
where, using transition-count notation $`n_{ij}`$ like that in the proof of Theorem 5.1,

$`f_l`$ $`=`$ $`{\displaystyle \underset{𝐳𝒵_l}{}}\psi _l(𝐳)p_{12}^{n_{12}(𝐳)}p_{21}^{n_{21}(𝐳)}\rho ^{n_{12}(𝐲)+n_{21}(𝐲)}(1p_{12})^{n_{11}(𝐱)}(1p_{21})^{n_{22}(𝐱)}`$ (6.3)
$`\times (1\rho p_{12})^{n_{11}(𝐲)}(1\rho p_{21})^{n_{22}(𝐲)},l=1,2.`$
Using induction on the $`\rho `$-degree of the polynomial $`p_{12}f_1p_{21}f_2`$ and Proposition A.1, it is easy to show that (6.2) holds as an equality in the ring of polynomials in the variables $`p_{12},p_{21},\rho `$ over the complex field.

For $`l=1,2`$, let

$$m_l=\mathrm{min}\{[n_{12}(𝐲)+n_{21}(𝐲)]:𝐳=(𝐱,𝐲)𝒵_l\text{ for some }𝐱\}$$

denote the highest power of $`\rho `$ that divides $`f_l`$. Then, by (6.2), $`m_2=m_1`$. Divide both sides of (6.2) by $`\rho ^{m_1}`$ and then set $`\rho `$ to $`0`$ to obtain

$$p_{12}g_1p_{21}g_2,$$
(6.4)
where

$$g_l:=\underset{𝐳𝒵_l^{}}{}\psi _l(𝐳)p_{12}^{n_{12}(𝐳)}p_{21}^{n_{21}(𝐳)}(1p_{12})^{n_{11}(𝐱)}(1p_{21})^{n_{22}(𝐱)},l=1,2,$$

with $`𝒵_l^{}:=\{𝐳=(𝐱,𝐲)𝒵_l:n_{12}(𝐲)+n_{21}(𝐲)=m_1\}\mathrm{}`$. But \[cf. (5.5) with $`N=2`$\], if $`𝐳=(𝐱,𝐲)𝒵_1𝒵_2`$, then $`n_{12}(𝐲)=m_1/2`$ and $`n_{21}(𝐲)=m_1/2`$. Dividing both sides of (6.4) by $`p_{12}^{m_1/2}p_{21}^{m_1/2}`$ we obtain the polynomial identity

$$p_{12}h_1p_{21}h_2,$$
(6.5)
where

$$h_l:=\underset{𝐱𝒳_l}{}\varphi _l(𝐱)p_{12}^{n_{12}(𝐱)}p_{21}^{n_{21}(𝐱)}(1p_{12})^{n_{11}(𝐱)}(1p_{21})^{n_{22}(𝐱)},l=1,2,$$
(6.6)
with

$$𝒳_l:=\{𝐱:\text{there exists }𝐲\text{ such that }(𝐱,𝐲)𝒵_l^{}\}\mathrm{},l=1,2$$

and, for $`𝐱𝒳_l`$,

$$\varphi _l(𝐱):=\underset{𝐲:(𝐱,𝐲)𝒵_l^{}}{}\psi _l(𝐱,𝐲)>0.$$

But (6.5) is the case $`N=2`$ of (5.4), which, as shown in the proof of Theorem 5.1, cannot hold. This contradiction establishes the theorem in the case $`M=N=2`$.

We leave to the reader the routine extension of the above proof to the case of arbitrary $`M`$ and $`N=2`$. A sketch is that now there are $`M1`$ parameters $`\rho _j`$, but by using the same sort of argument for each $`\rho _j`$ in succession that we used above for $`\rho `$, one again obtains a contradiction of the form (6.5) \[with $`\varphi _l(𝐱)>0`$ for all $`𝐱𝒳_l\mathrm{}`$, $`l=1,2`$\].

We complete the proof of the theorem by showing that an algorithm for interruptible exact sampling using $`M`$ independent trajectories from chains with $`N3`$ states could be converted into one for two-state chains.

Indeed, while watching independent trajectories of $`M`$ generic irreducible aperiodic two-state chains $`X_1,\mathrm{}`$, $`X_M`$ with common (unknown) stationary distribution $`\pi =(\pi _1,\pi _2)`$, contemporaneously construct $`M`$ independent irreducible aperiodic $`N`$-state chains $`Y_1,\mathrm{},Y_M`$ by letting $`Y_i(t)=1`$ whenever $`X_i(t)=1`$ and selecting an independent uniform random value from $`\{2,\mathrm{},N\}`$ as the value of $`Y_i(t)`$ at each time $`t`$ such that $`X_i(t)=2`$. The stationary distribution for each $`Y_i`$ is $`(\pi _1,\pi _2/(N1),\mathrm{},\pi _2/(N1))`$. Applying the size-$`N`$ algorithm to $`Y_1,\mathrm{},Y_M`$, suppose the output state is $`S^{}`$. To finish the construction of the two-state algorithm, output $`S:=\mathrm{min}\{S^{},2\}`$. ∎

###### Remark 6.2.

(a) Any two-state chain is reversible, as are the chains $`Y_i`$ constructed in the preceding paragraph. Thus Theorem 6.1 remains true even if we assume that the chains are all reversible.

(b) Similarly, as mentioned in Section 1, interruptible exact sampling from the stationary distribution is not possible when one observes only a single trajectory from an irreducible aperiodic reversible finite-state chain.

(c) For $`N3`$ we do not know whether Theorem 5.1 remains true if one assumes that the chain is reversible.

Acknowledgment. We thank Dan Naiman for helpful discussions related to the Appendix.

## A Appendix: Polynomials

In this Appendix we establish two basic facts about polynomials; these were used in the proof of Theorem 5.1. Throughout the Appendix, we write $`FG`$ to indicate that $`F`$ and $`G`$ are the same element in the ring of polynomials (in some specified finite collection of variables) over the complex field.

The first fact is quite simple. For completeness, we include an elementary proof.

###### Proposition A.1.

Let

$$𝐱=(x_{ij}:1in,1jk_i)$$

be a double array of variables, where $`n0`$ and $`k_i1`$ for $`1in`$. If $`F(𝐱)`$ is a polynomial expression that vanishes whenever $`x_{ij}0`$ for all $`i,j`$ and $`_{j=1}^{k_i}x_{ij}1`$ for all $`i`$, then $`F(𝐱)0`$.

###### Proof.

Let $`K:=_{i=1}^nk_i`$. The proof is by (strong) induction on $`\kappa :=K+\mathrm{deg}F`$, for which (if $`F`$ is not the zero polynomial) the smallest possible value is $`_{1i0}1+0=0`$. The base case $`\kappa =0`$ of the induction is trivial.

For the induction step we may assume $`n1`$ and $`k_n1`$. Dividing the polynomial $`F(𝐱)`$ by $`x_{n,k_n}`$, we can write

$$F(𝐱)x_{n,k_n}F_1(𝐱)+F_2(𝐱^{})$$
(A.1)
for polynomials $`F_1`$ and $`F_2`$, where the variables collection $`𝐱^{}`$ excludes the single variable $`x_{n,k_n}`$. Setting $`x_{n,k_n}`$ to $`0`$ in (A.1), we see that $`F_2(𝐱^{})`$ is a polynomial satisfying the hypothesis of the proposition; and (in obvious notation) $`K_2=K1`$ and $`\mathrm{deg}F_2\mathrm{deg}F`$, so that $`\kappa _2<\kappa `$. By induction, $`F_2(𝐱^{})0`$, and so from (A.1) we now have $`F(𝐱)x_{n,k_n}F_1(𝐱)`$. But now $`K_1=K`$ and $`\mathrm{deg}F_1=\mathrm{deg}F1`$, so that $`\kappa _1=\kappa 1`$, and one sees that $`F_1(𝐱)`$ satisfies the hypothesis of the proposition. By induction, $`F_1(𝐱)0`$; we conclude that $`F(𝐱)0`$, as desired. ∎

As is well known (e.g., , Chapter 4), for any $`n1`$ the ring $`𝐂[x_1,\mathrm{},x_n]`$ of polynomials in the variables $`x_1,\mathrm{},x_n`$ over the complex field $`𝐂`$ is a unique factorization domain. This means that every nonzero polynomial in $`𝐂[x_1,\mathrm{},x_n]`$ can be written uniquely (up to complex scalar multiples) as a (possibly empty) finite product of irreducible polynomials. (A polynomial is said to be irreducible if it cannot be factored as the product of two nonconstant polynomials.)

###### Lemma A.2.

The polynomial $`w_1`$ \[i.e., the polynomial in the $`(N1)^2`$ variables $`p_{ij}`$ with $`i,j[N]`$ and $`i\{1,j\}`$ defined in Section 3.1\] is irreducible over the complex field.

###### Proof.

The proof is by induction on $`N`$. For $`N=1`$, the polynomial $`w_11`$ (in no variables) is certainly irreducible. For $`N=2`$, the polynomial $`w_1p_{21}`$ in the single variable $`p_{21}`$ is irreducible. To carry out the induction step for $`N3`$, we will use another induction, on $`l`$, to prove the following claim.

Claim. For $`3lN+1`$, let $`y_l`$ denote the polynomial in $`(N1)^2(N+1l)`$ variables obtained from $`w_1`$ by setting $`p_{m1}`$ to $`0`$ for $`lmN`$. Then $`y_l`$ is irreducible for $`4lN+1`$.

To prove the claim, we begin by noting that $`y_3`$ has the factorization

$$y_3p_{21}\omega _2,$$
(A.2)
where the polynomial

$$\omega _2=\omega _2((p_{ij}:2i,jN\text{ and }i\{2,j\}))$$

is obtained from the polynomial $`w_1`$ for the state space $`[N1]`$ by changing each variable name from $`p_{ij}`$ to $`p_{i+1,j+1}`$. By the induction hypothesis for our $`N`$-induction, $`\omega _2`$ is irreducible. Since $`p_{21}`$ is clearly irreducible, we conclude that (A.2) is a prime factorization of $`y_3`$.

We now treat the base case $`l=4`$ of our $`l`$-induction. Observe that $`y_40`$ (consider, e. g., the tree $`NN1\mathrm{}21`$) and that $`y_4`$ is linear in $`p_{31}`$. If $`y_4`$ is reducible, then we can write

$$y_4(g_1p_{31}+g_2)g_3,$$
(A.3)
where $`g_i`$ is a polynomial free of the variable $`p_{31}`$ ($`i=1,2,3`$) and $`g_3`$ is nonconstant. If we now set $`p_{31}`$ to $`0`$ in (A.3), the result is $`y_3g_2g_3`$. From the prime factorization (A.2) we conclude that either $`p_{21}`$ or $`\omega _2`$ divides $`g_3`$. But this is wrong: (i) $`p_{21}`$ does not divide $`g_3`$ because it clearly does not divide $`y_4`$ (consider, e. g., the tree $`NN1\mathrm{}4231`$), and (ii) $`\omega _2`$ does not divide $`g_3`$ because (we claim) it, too, fails to divide $`y_4`$. (Indeed, setting $`p_{m2}`$ to $`0`$ for $`3mN`$ causes $`\omega _2`$—but clearly not $`y_4`$—to vanish.) From this contradiction we conclude that $`y_4`$ is irreducible, establishing the $`l`$-induction base case.

For the $`l`$-induction step, let $`l5`$. If $`y_l`$ is reducible, then we can write

$$y_l(h_1p_{l1,1}+h_2)h_3,$$
(A.4)
where $`h_i`$ is a polynomial free of the variable $`p_{l1,1}`$ ($`i=1,2,3`$) and $`h_3`$ is nonconstant. If we now set $`p_{l1,1}`$ to $`0`$ in (A.4), the result is $`y_{l1}h_2h_3`$. By the $`l`$-induction hypothesis, it must be that $`h_3`$ is a nonzero complex scalar multiple of $`y_{l1}`$; from (A.4) we then deduce that $`y_{l1}`$ divides $`y_l`$. But this is wrong, because setting $`p_{m1}`$ to $`0`$ for $`2ml2`$ causes $`y_{l1}`$—but clearly not $`y_l`$—to vanish. From this contradiction we conclude that $`y_l`$ is irreducible, completing the $`l`$-induction.

Finally, set $`l`$ to $`N+1`$ in the claim to find that $`w_1y_{N+1}`$ is irreducible, completing the $`N`$-induction and the proof of the lemma. ∎