text:
# MINIMAL SEESAW MECHANISM

## I Introduction

There is now a strong evidence for neutrino oscillations, especially through the SuperKamiokande, K2K, and SNO, KamLAND experiments exp . Neutrino oscillations are naturally accounted for if neutrinos have small masses, so that leptons can mix in a similar way as quarks do mix . Moreover, small neutrino masses can be achieved by means of the seesaw mechanism ss . In this framework, the effective (Majorana) mass matrix of neutrinos $`M_L`$ is related to the Dirac neutrino mass matrix $`M_\nu `$ and the heavy Majorana neutrino mass matrix $`M_R`$ by the relation

$$M_LM_\nu M_R^1M_\nu ^T.$$
(1)
As a matter of fact, the seesaw formula (1) is valid at the high $`M_R`$ scale, and therefore one should determine both $`M_\nu `$ and $`M_L`$ at that scale, in order to find a consistent model. The effective matrix $`M_L`$ is partially described at the low scale through the analysis of several experiments fs1 . On the other hand, the Dirac matrix $`M_\nu `$ is based on theoretical hints. Both have to be renormalized to the $`M_R`$ scale. Then the problem is to find models for $`M_\nu `$ and $`M_R`$ which reproduce the phenomenological forms of $`M_L`$ according to the master relation (1). Such a problem has been addressed in many papers (see, for instance, the review bd ). In the present article we consider a structure for $`M_\nu `$ which is usually adopted for charged fermion mass matrices, and minimal models for $`M_R`$. We select minimal forms which are compatible with phenomenology.

In section II we discuss the effective neutrino mass matrix. In section III we describe the Dirac and Majorana mass matrices of our minimal framework. Then, in section IV, the seesaw formula is applied and the resulting neutrino mass matrix is compared to phenomenology. A brief discussion is finally proposed.

## II The effective neutrino mass matrix

Experimental informations on neutrino oscillations imply that the lepton mixing matrix is given by

$$U\left(\begin{array}{ccc}\frac{\sqrt{2}}{\sqrt{3}}& \frac{1}{\sqrt{3}}& ϵ\text{e}^{\text{i}\delta }\\ \frac{1}{\sqrt{6}}\frac{1}{\sqrt{3}}ϵ\text{e}^{\text{i}\delta }& \frac{1}{\sqrt{3}}\frac{1}{\sqrt{6}}ϵ\text{e}^{\text{i}\delta }& \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{6}}\frac{1}{\sqrt{3}}ϵ\text{e}^{\text{i}\delta }& \frac{1}{\sqrt{3}}\frac{1}{\sqrt{6}}ϵ\text{e}^{\text{i}\delta }& \frac{1}{\sqrt{2}}\end{array}\right)\text{diag}(\text{e}^{\text{i}\phi _1/2},\text{e}^{\text{i}\phi _2/2},1)$$
(2)
in the standard parametrization, where $`ϵ<0.16`$, $`0<\delta <2\pi `$, $`0<\phi _1,\phi _2<2\pi `$ (see, for example, Ref.tribi ). Moreover, neutrino oscillations determine the following square mass differences,

$$\mathrm{\Delta }m_{32}^2=m_3^2m_2^2310^3\text{eV}^2,$$
(3)
$$\mathrm{\Delta }m_{21}^2=m_2^2m_1^2710^5\text{eV}^2,$$
(4)
where $`m_1,m_2,m_3`$ are the effective neutrino masses. Then, in the basis where the charged lepton mass matrix is diagonal, $`M_L`$ is obtained by means of the transformation

$$M_L=U^{}D_LU^{},$$
(5)
with $`D_L=\text{diag}(m_1,m_2,m_3)`$. Neglecting $`ϵ`$ in $`U`$, except for $`U_{e3}`$, the calculation leads to

$$M_{ee}ϵ^2m_3+\frac{m_2}{3}+2\frac{m_1}{3}$$

$$M_{e\mu }ϵ\frac{m_3}{\sqrt{2}}+\frac{m_2}{3}\frac{m_1}{3}$$

$$M_{e\tau }ϵ\frac{m_3}{\sqrt{2}}\frac{m_2}{3}+\frac{m_1}{3}$$

$$M_{\mu \mu }M_{\tau \tau }\frac{m_3}{2}+\frac{m_2}{3}+\frac{m_1}{6}$$

$$M_{\mu \tau }\frac{m_3}{2}\frac{m_2}{3}\frac{m_1}{6}$$

not writing phases, which can be inserted by $`ϵϵ\text{e}^{\text{i}\delta }`$, $`m_1m_1\text{e}^{\text{i}\phi _1}`$, $`m_2m_2\text{e}^{\text{i}\phi _2}`$.

Since $`\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{32}^2`$, we have two main mass spectra for light neutrinos, the normal hierarchy $`m_1<m_2m_3`$, and the inverse hierarchy $`m_1m_2m_3`$. For the normal hierarchy the dominant elements are given by

$$M_L\left(\begin{array}{ccc}0& 0& 0\\ 0& \frac{1}{2}& \frac{1}{2}\\ 0& \frac{1}{2}& \frac{1}{2}\end{array}\right)m_3,$$
(6)
with $`m_3\sqrt{\mathrm{\Delta }m_{32}^2}`$, and for the inverse hierarchy they are given by

$$M_L\left(\begin{array}{ccc}1& 0& 0\\ 0& \frac{1}{2}& \frac{1}{2}\\ 0& \frac{1}{2}& \frac{1}{2}\end{array}\right)m_{1,2},$$
(7)
with $`m_{1,2}\sqrt{\mathrm{\Delta }m_{32}^2}`$. Both contain a democratic $`\mu \tau `$ block, due to near maximal $`U_{\mu 3}`$. The difference stands in the element $`ee`$, which is suppressed in (6) but dominant in (7). Now, according to Ref.fs2 , the general structure of $`M_L`$ is not changed by renormalization. Therefore, we can take matrices (6) and (7) as simple forms at the high scale, the zero elements meaning suppressed with respect to dominant elements. They correspond to distinct predictions for the double beta decay parameter $`M_{ee}`$, since for the normal hierarchy we get $`M_{ee}\sqrt{\mathrm{\Delta }m_{21}^2}`$, while for the inverse hierarchy we have $`M_{ee}\sqrt{\mathrm{\Delta }m_{32}^2}`$.

## III Dirac and Majorana mass matrices

In order to apply the seesaw formula, we need the expression of the Dirac and Majorana mass matrices. We take a typical form for the three mass matrices of charged fermions cf :

$$M_d\left(\begin{array}{ccc}0& \sqrt{m_dm_s}& 0\\ \sqrt{m_dm_s}& m_s& \sqrt{m_dm_b}\\ 0& \sqrt{m_dm_b}& m_b,\end{array}\right)$$
(8)
$$M_u\left(\begin{array}{ccc}0& \sqrt{m_um_c}& 0\\ \sqrt{m_um_c}& m_c& \sqrt{m_um_t}\\ 0& \sqrt{m_um_t}& m_t\end{array}\right)$$
(9)
for down and up quarks, and

$$M_e\left(\begin{array}{ccc}0& \sqrt{m_em_\mu }& 0\\ \sqrt{m_em_\mu }& m_\mu & \sqrt{m_em_\tau }\\ 0& \sqrt{m_em_\tau }& m_\tau \end{array}\right)$$
(10)
for charged leptons. Then, since $`M_eM_d`$, a natural choice is also $`M_\nu M_u`$. In fact, we have br  $`m_d/m_sm_s/m_b\lambda ^2`$, $`m_e/m_\mu m_\mu /m_\tau \lambda ^2`$, and $`m_u/m_cm_c/m_t\lambda ^4`$, where $`\lambda =0.22`$ is the Cabibbo parameter. The renormalization of quark mass matrices does not affect their expression in terms of powers of $`\lambda `$ bbo . Therefore, for the Dirac neutrino mass matrix we take

$$M_\nu \left(\begin{array}{ccc}0& a& 0\\ a& b& c\\ 0& c& 1\end{array}\right)m_t,$$
(11)
with $`abc1`$. As order in $`\lambda `$ we have $`a\lambda ^6`$, $`c\lambda ^4`$. Expressions (8) and (9) lead to small quark mixings, while lepton mixings $`U_{e2}`$ and $`U_{\mu 3}`$ are large. The matrix $`M_R`$ should produce, through the seesaw formula, large lepton mixings smir .

For this Majorana mass matrix we consider minimal forms. These include matrices with four texture zeros:

$$M_R=\left(\begin{array}{ccc}0& A& 0\\ A& 0& 0\\ 0& 0& B\end{array}\right)m_R,$$
(12)
$$M_R=\left(\begin{array}{ccc}A& 0& 0\\ 0& 0& B\\ 0& B& 0\end{array}\right)m_R,$$
(13)
$$M_R=\left(\begin{array}{ccc}0& 0& A\\ 0& B& 0\\ A& 0& 0\end{array}\right)m_R,$$
(14)
and matrices with three texture zeros, that is the diagonal form

$$M_R=\left(\begin{array}{ccc}A& 0& 0\\ 0& B& 0\\ 0& 0& C\end{array}\right)m_R,$$
(15)
and the Zee-like form stech

$$M_R=\left(\begin{array}{ccc}0& A& B\\ A& 0& C\\ B& C& 0\end{array}\right)m_R.$$
(16)
Parameters $`A,B,C`$ can take values $`1`$ and $`r<1`$. Therefore, in such matrices there is a scale factor, related to $`m_R`$, and possibly one hierarchy parameter $`r`$.

## IV Seesaw mechanism

In this section we calculate the effective neutrino mass matrix by means of the seesaw formula (1) on mass matrices discussed in the previous section. Then we look for structure of the kind (6) and (7). We exclude possible cancellations during our analysis.

### IV.1 Four texture zeros

Matrix (12) leads to

$$M_L\left(\begin{array}{ccc}0& \frac{a^2}{A}& 0\\ \frac{a^2}{A}& \frac{2ab}{A}+\frac{c^2}{B}& \frac{ac}{A}+\frac{c}{B}\\ 0& \frac{ac}{A}+\frac{c}{B}& \frac{1}{B}\end{array}\right)\frac{m_t^2}{m_R}.$$
(17)
Condition $`M_{\mu \tau }M_{\tau \tau }`$ gives $`A/Bac`$. Then $`A=rac`$ and $`B=1`$. Condition $`M_{\mu \mu }M_{\tau \tau }`$ is satisfied as a consequence. See also Ref.bo  for a discussion on this structure.

Matrix (13) leads to

$$M_L\left(\begin{array}{ccc}0& \frac{ac}{B}& \frac{a}{B}\\ \frac{ac}{B}& \frac{a^2}{A}+\frac{2bc}{B}& \frac{b}{B}+\frac{c^2}{B}\\ \frac{a}{B}& \frac{b}{B}+\frac{c^2}{B}& \frac{2c}{B}\end{array}\right)\frac{m_t^2}{m_R}.$$
(18)
Here condition $`M_{\mu \tau }M_{\tau \tau }`$ is satisfied while $`M_{\mu \mu }M_{\tau \tau }`$ requires $`A/Ba^2/c`$. Hence $`A=ra^2/c`$ and $`B=1`$.

Matrix (14) leads to

$$M_L\left(\begin{array}{ccc}\frac{a^2}{B}& \frac{ab}{B}& \frac{ac}{B}\\ \frac{ab}{B}& \frac{b^2}{B}+\frac{2ac}{A}& \frac{bc}{B}+\frac{a}{A}\\ \frac{ac}{B}& \frac{bc}{B}+\frac{a}{A}& \frac{c^2}{B}\end{array}\right)\frac{m_t^2}{m_R},$$
(19)
so that $`M_{\mu \tau }M_{\tau \tau }`$ gives $`B/Ac^2/a`$. Then $`A=1`$ and $`B=rc^2/a`$. The condition $`M_{\mu \mu }M_{\tau \tau }`$ is valid as a consequence.

The normal hierarchy is achieved in all cases. However, note the three different values for the scale $`m_R`$, that is $`m_t^2/m_3`$, $`cm_t^2/m_3`$, $`am_t^2/m_3`$, respectively.

### IV.2 Diagonal form

In this case the effective neutrino mass matrix is given by

$$M_L\left(\begin{array}{ccc}\frac{a^2}{B}& \frac{ab}{B}& \frac{ac}{B}\\ \frac{ab}{B}& \frac{a^2}{A}+\frac{b^2}{B}+\frac{c^2}{C}& \frac{bc}{B}+\frac{c}{C}\\ \frac{ac}{B}& \frac{bc}{B}+\frac{c}{C}& \frac{c^2}{B}+\frac{1}{C}\end{array}\right)\frac{m_t^2}{m_R}.$$
(20)
Here condition $`M_{\mu \tau }M_{\tau \tau }`$ gives $`B/Cc^2`$. Then $`C=1`$ and $`B=rc^2`$. Both $`A=r`$ and $`A=1`$ are consistent with $`M_{\mu \mu }M_{\tau \tau }`$. The normal hierarchy is obtained. The scale $`m_R`$ is given by $`m_t^2/m_3`$. Large lepton mixing can indeed be obtained even by means of small mixing in $`M_\nu `$ and zero mixing in $`M_R`$ (see Ref.afm ). In particular, for $`ABc^2`$ we get

$$M_L\left(\begin{array}{ccc}k^2& k& k\\ k& 1& 1\\ k& 1& 1\end{array}\right)\frac{m_t^2}{m_R},$$
(21)
with $`k=a/c`$. This form of $`M_L`$ has already been proposed several times vis . Moreover, the same form is realized in (19) for $`Bc^2/a`$, but with the scale $`m_R`$ suppressed by the factor $`a`$ with respect to (20).

### IV.3 Zee-like form

In this case, apart from an overall scale $`m_t^2/2m_R`$, we get the following approximate effective matrix

$$\left(\begin{array}{ccc}[\frac{a^2B}{AC}]& [\frac{abB}{AC}+\frac{a^2}{A}+\frac{ac}{C}]& [\frac{acB}{AC}+\frac{a}{C}]\\ & [\frac{a^2C}{AB}\frac{b^2B}{AC}\frac{c^2A}{BC}+\frac{2ab}{A}+\frac{2ac}{B}+\frac{2bc}{C}]& [\frac{bcB}{AC}\frac{cA}{BC}+\frac{ac}{A}+\frac{a}{B}+\frac{c^2}{C}+\frac{b}{C}]\\ & & [\frac{c^2B}{AC}\frac{A}{BC}+\frac{2c}{C}]\end{array}\right).$$
(22)
In order to have a useful $`\mu \tau `$ block, the leading terms must be those with $`AC`$ in the denominator. Then the normal hierarchy is achieved for $`A=C=rc`$ and $`B=1`$. Here the scale $`m_R`$ is about $`m_t^2/m_3`$.

## V Discussion

We have studied the seesaw mechanism assuming simple forms of the fermion mass matrices and in particular minimal forms for the heavy neutrino mass matrix, which contain four or three texture zeros, a scale factor and a hierarchy parameter. Our minimal framework allows only the normal hierarchy for light neutrinos and not the inverse hierarchy. Large lepton mixing is achieved by tuning the hierarchy parameter $`r`$ in the heavy neutrino mass matrix.

There is another possible mass spectrum for neutrinos, the degenerate spectrum, $`m_1m_2m_3`$, which gives the dominant elements $`M_L\text{diag}(1,1,1)m_{1,2,3}`$ or

$$M_L\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right)m_{1,2,3}.$$
(23)
One can easily check that such forms are not reproduced in our minimal framework. However, both the inverse hierarchy and the degenerate spectrum can be achieved in some nonminimal models fal . Indeed, generally it is quite hard to yield degeneracy in $`M_L`$ from hierarchy in $`M_\nu `$ by means of $`M_R`$ in the seesaw formula.

The present framework could also be embedded in a unified $`SO(10)`$ model. In fact, the relations $`M_eM_d`$ and $`M_\nu M_u`$ can be the result of a quark-lepton symmetry, and the high $`M_R`$ scale can as well be related to the unification or intermediate breaking scale of the supersymmetric or nonsupersymmetric model, respectively falc . Then, matrix (13) and especially matrix (14) correspond to the nonsupersymmetric model, while matrices (12), (15) and possibly (16) correspond to the supersymmetric model.

We thank F. Buccella for discussions.


format: text

text:
# Numerical Estimation of the Asymptotic Behaviour of Solid Partitions of an Integer

## I Introduction

Combinatorial enumeration problems arise naturally in many problems of statistical physics. The number of partitions of an integer (see andrews ; stanley  for an introduction) is one such enumeration problem with a history dating back to Euler. Examples of applications to physical problems include the $`q\mathrm{}`$ Potts model wurollet , compact lattice animals wurollet ; bhatia , crystal growth temperley1 , lattice polygons rensburg , Bose-Einstein statistics holthaus ; temperley2  and dimer coverings elser . The solution to the integer partitioning problem is known for $`1`$-dimensional and $`2`$-dimensional partitions. However, not much is known about higher dimensional partitions. Numerical estimation of the asymptotic behaviour of these higher dimensional partitions could lead to theoretical insights. In this paper, we determine numerically the leading asymptotic behaviour of $`3`$-dimensional partitions by exact enumeration and Monte Carlo techniques.

A $`1`$-dimensional or linear partition of an integer is a decomposition into a sum of positive integers in which the summands are ordered from largest to smallest. A $`2`$-dimensional or plane partition of an integer is a decomposition into a sum of smaller positive integers which are arranged on a plane. The ordering property generalises to the summands being non-increasing along both the rows and the columns. Generalisation to $`d`$-dimension is straightforward. Consider a $`d`$-dimensional hyper cubic lattice with sites labelled by $`𝐢=(i_1,i_2,\mathrm{},i_d)`$, where $`i_k=1,2,\mathrm{}`$. An integer height $`h(𝐢)`$ (corresponding to a summand) is associated with site $`𝐢`$. A $`d`$-dimensional partition of a positive integer $`n`$ is a configuration of heights such that

$`h(𝐢)`$ $``$ $`0,`$
$`h(𝐢)`$ $``$ $`\underset{1kd}{\mathrm{max}}h(i_1,i_2,\mathrm{},i_k+1,\mathrm{},i_d),`$ (1)
$`{\displaystyle \underset{𝐢}{}}h(𝐢)`$ $`=`$ $`n.`$
The second condition in Eq. (1) means that the heights $`h(𝐢)`$ are non-increasing in each of the $`d`$ lattice directions. As an illustration, the linear and plane partitions of $`4`$ are shown in Fig. 1(a) and 1(b) respectively.

Let $`p_d(n)`$ denote the number of partitions of $`n`$ in $`d`$-dimensions. The generating function $`G_d(q)`$ is then defined as

$$G_d(q)=\underset{n=0}{\overset{\mathrm{}}{}}p_d(n)q^n,$$
(2)
where $`p_d(0)1`$. The generating function for linear partitions is due to Euler and is

$$G_1(q)=\underset{k=1}{\overset{\mathrm{}}{}}(1q^k)^1,$$
(3)
and $`p_1(n)`$ for large $`n`$ varies as hardyramanujan

$$p_1(n)\frac{1}{4n\sqrt{3}}\mathrm{exp}\left(\pi \sqrt{\frac{2n}{3}}\right),n1.$$
(4)
The corresponding formulae for plane partitions are macmahon

$$G_2(q)=\underset{k=1}{\overset{\mathrm{}}{}}(1q^k)^k,$$
(5)
and

$$p_2(n)\frac{c_2}{n^{25/36}}\mathrm{exp}\left(\alpha _2n^{2/3}\right),n1,$$
(6)
where $`c_2=0.40099\mathrm{}`$ and $`\alpha _2=2.00945\mathrm{}`$ wright1 . While the generating functions for three and higher dimensional partitions are not known, it is known that $`lim_n\mathrm{}\mathrm{ln}\left(p_d(n)\right)/n^{d/(d+1)}`$ has a finite non-zero limit bhatia . We define $`\alpha _d`$ to be

$$\alpha _d=\underset{n\mathrm{}}{lim}\frac{\mathrm{ln}\left(p_d(n)\right)}{n^{d/(d+1)}}.$$
(7)
In this paper, we numerically estimate $`\alpha _3`$ for solid ($`3`$-dimensional) partitions to be

$$\alpha _3=1.79\pm 0.01.$$
(8)
Generalising the results for linear and plane partitions to higher dimensions, MacMahon suggested that the generating function for $`d`$-dimensional partitions could be macmahon

$$G_d^{(m)}(q)=\underset{k=1}{\overset{\mathrm{}}{}}(1q^k)^{\left(\genfrac{}{}{0pt}{}{k+d2}{d1}\right)}.$$
(9)
Equation (9) is usually known as the MacMahon conjecture. Clearly, $`G_d^{(m)}(q)`$ is the correct result for $`d=1,2`$. However, it is known that $`G_d^{(m)}(q)`$ is different from $`G_d(q)`$ for all $`d3`$ atkin ; wright2 . In particular, in $`3`$-dimensions

$$G_3^{(m)}(q)=\underset{k=1}{\overset{\mathrm{}}{}}(1q^k)^{k(k+1)/2},$$
(10)
gives the wrong answer for the number of solid partitions of $`6,7,8,\mathrm{}`$.

Let $`p_d^{(m)}(n)`$ denote the coefficient of $`q^n`$ in $`G_d^{(m)}(q)`$. The asymptotic behaviour of $`p_3^{(m)}(n)`$ for large $`n`$ can be determined from Eq. (10) by the method of steepest descent. In Appendix A, we present a heuristic derivation of the large $`n`$ behaviour of the coefficient of $`q^n`$ in the expansion of the infinite product

$$F(q)=\underset{k=1}{\overset{\mathrm{}}{}}(1q^k)^{a_1k^2a_2ka_3}.$$
(11)
Substituting $`a_1=1/2`$, $`a_2=1/2`$ and $`a_3=0`$ in Eq. (25), we obtain

$$p_3^{(m)}(n)\frac{c_3^{(m)}}{n^{61/96}}\mathrm{exp}\left(\alpha _3^{(m)}n^{3/4}+\beta _3^{(m)}n^{1/2}+\gamma _3^{(m)}n^{1/4}\right),$$
(12)
where $`c_3^{(m)}`$ is a constant and

$`\alpha _3^{(m)}`$ $`=`$ $`{\displaystyle \frac{2^{7/4}\pi }{15^{1/4}3}}=1.7898\mathrm{},`$ (13)
$`\beta _3^{(m)}`$ $`=`$ $`{\displaystyle \frac{\sqrt{15}\zeta (3)}{\sqrt{2}\pi ^2}}=0.3335\mathrm{},`$ (14)
$`\gamma _3^{(m)}`$ $`=`$ $`{\displaystyle \frac{15^{5/4}\zeta (3)^2}{2^{7/4}\pi ^5}}=0.0414\mathrm{}.`$ (15)
Comparing the values for $`\alpha _3^{(m)}`$ in Eq. (8) and $`\alpha _3`$ in Eq. (13), we conclude that the MacMahon conjecture, though not exact, could still give the correct leading asymptotic behaviour for solid partitions. The value of $`\alpha _3^{(m)}`$ is a function of only $`a_1`$ in Eq. (11). Thus, if we assume that the asymptotic behaviour for solid partitions is correctly captured by a product form as in Eq. (11), then it should have the form $`_k(1q^k)^{(1/2\pm 0.012)k^2}`$.

The rest of the paper is organised as follows. In Sec. II, we present the results of the exact enumeration study. In Sec. III, we describe the Monte Carlo algorithm and the simulation results for plane and solid partitions. Finally, we conclude with a summary and conclusions in Sec. IV.

## II Exact enumeration

Previous attempts at studying solid partitions on the computer have been based on exact enumeration atkin ; knuth ; bratley ; huangwu . Tables of $`p_3(n)`$ exist for $`n`$ up to $`28`$ knuth . The table is extended up to $`n=50`$ in this paper by using the standard back tracking algorithm martin . The algorithm is made faster by the following. Partitions that are related to each other by symmetry operations are counted only once and multiplied by the corresponding symmetry factor. Also, parts of the partition that are restricted to planes are generated by using the known generating functions for plane partitions. In Table 1, we list the solid partitions from $`n=29`$ to $`n=50`$. For solid partitions up to $`n=28`$, we refer to Ref. knuth .

We compare the exact enumeration results with the answer predicted by the MacMahon conjecture. In Fig. 2, we show the variation of $`\mathrm{ln}[p_3(n+1)/p_3(n)]`$ with $`n`$ for both $`p_3(n)`$ as well as $`p_3^{(m)}(n)`$. While there seems to be a good agreement, we are unable to determine the precise asymptotic behaviour of $`p_3(n)`$ from these $`50`$ numbers. This is possibly due to the presence of strong corrections to the leading asymptotic behaviour. It is difficult to further extend the table of solid partitions due to the large computing times involved. One possible method of probing larger values of $`n`$ is to use Monte Carlo simulations. These are described in Sec. III.

## III Monte Carlo simulation

### III.1 Algorithm

We use an algorithm proposed recently by Wang and Landau for measuring density of states in spin systems landau . The algorithm is described below. Consider a $`N_x\times N_y\times N_z`$ lattice with initial height $`h(𝐢)`$ assigned to each lattice point in such a way that the configuration is a valid solid partition. To each positive integer $`n`$ is associated a histogram $`H(n)`$ and the number of solid partitions $`p_3(n)`$. The histogram $`H(n)`$ keeps track of the number of times solid partitions of $`n`$ have been visited during the simulations. The algorithm is based upon the fact that if the probability of transition to a solid partition $`n`$ is proportional to $`[p_3(n)]^1`$, then a flat distribution is generated for the histogram $`H(n)`$. At the start of the program $`H(n)=0`$ and $`p_3(n)=1`$.

A site is chosen randomly and as a trial move the height $`h(𝐢)`$ is increased or decreased by one with equal probability, provided that the new state is an allowed partition. If the new state is an allowed partition, then the move is accepted with probability

$$\text{Prob}(n_{old}n_{new})=\mathrm{min}[\frac{p_3(n_{old})}{p_3(n_{new})},1],$$
(16)
where $`n_{old}`$ and $`n_{new}`$ are the sum of heights for the old and new states respectively, i.e, $`n_{new}=n_{old}+1`$, $`n_{old}1`$ or $`n_{old}`$ depending on whether the height increased by one, decreased by one or remained the same. The histogram $`H(n)`$ and $`p_3(n)`$ are updated as

$`H(n_{new})`$ $`=`$ $`H(n_{new})+1,`$ (17)
$`p_3(n_{new})`$ $`=`$ $`f_ip_3(n_{new}),`$ (18)
where $`f_i`$ is a modification factor greater than $`1`$.

These steps are repeated until a flat histogram is created; in practice, this means that $`H(n)_{min}>cH(n)_{ave}`$, where $`c`$ is a flatness criteria typically between $`0.75`$ and $`0.9`$ while $`H(n)_{min}`$ is the minimum of the $`H(n)`$’s and $`H(n)_{ave}`$ is the average of the $`H(n)`$’s. When the histogram becomes flat, the modification factor $`f_i`$ is changed to

$$f_{i+1}=f_i^a,$$
(19)
and the histogram is reset to zero. The exponent $`a`$ is less than $`1`$ and defines the smoothness of the iteration. The program runs until $`f`$ is less than a predetermined value $`f_{final}`$.

Note that the algorithm does not obey detailed balance during the simulations. However, in the limit $`f_i1`$ when $`p_3(n)`$ takes its correct value, the system does obey detailed balance with the weight of a solid partition of $`n`$ being proportional to $`[p_3(n)]^1`$.

The algorithm can be made faster by adopting certain ideas from the $`N`$-fold method binder ; bortz . In this modification, sites which cannot undergo a valid move are never chosen. To do so, we define four classes. (i) c1: sites at which the height can only increase. (ii) c2: sites at which the height can either increase or decrease. (iii) c3: sites at which the height can only decrease. (iv) c4: an auxiliary class to help implementation of detailed balance.

First, we update the histogram $`H(n)`$ and $`p_3(n)`$ by $`\mathrm{\Delta }`$ times with $`\text{Prob}(\mathrm{\Delta }=k)=p^k(1p^{})`$, where $`p^{}=(|c_4|+|c_1|/2+|c_3|/2)/(|c_1|+|c_2|+|c_3|+|c_4|)`$ and $`|c_k|`$ denotes the number of elements in the class $`c_k`$. We then choose one of the classes c1, c2, c3 with probabilities $`|c_1|/(|c_1|+2|c_2|+|c_3|)`$, $`2|c_2|/(|c_1|+2|c_2|+|c_3|)`$ and $`|c_3|/(|c_1|+2|c_2|+|c_3|)`$ respectively. A site is picked up randomly from the chosen class and a trial move is decided, for example if site $`𝐢`$ from class c1 is chosen, the trial move is $`h(𝐢)=h(𝐢)+1`$, while in the case of class c2 the height increases or decreases with equal probability. Finally, we either accept or reject the trial move according to Eq. (16) and update the histogram and $`p_3(n)`$ according to Eqs. (17) and (18). With this construction, only valid trial moves are chosen at each time step and the algorithm becomes considerably faster. The role of the class c4 is to make the algorithm obey detailed balance asymptotically, i.e. when $`f_i1`$. We define $`|c4|=C(|c_1|+|c_2|+|c_3|)`$, where $`C`$ is some large enough constant.

Further speeding up can be done by dividing the interval 1–8000 to smaller slightly overlapping intervals ($`14`$ in our case) and the simulation is done for each interval. After the simulations are over, these intervals can be joined together to produce $`p_3(n)`$. The distribution is finally normalised by setting $`p_3(1)=1`$.

The parameters we have used for the simulations are $`c=0.85`$, $`a=1.0/1.4`$, $`C=2000`$, $`f_0=2.5`$ and $`f_{final}=1.0000099`$ (corresponding to $`35`$ iterations). Changing these parameters slightly does not change the final outcome of the simulation. The lattice sizes used for plane and solid partitions were $`100\times 100\times 1`$ and $`50\times 50\times 50`$ respectively. Random numbers were generated using standard RANMAR algorithm. With this setup, a typical run producing one $`p_3(n)`$ for $`n`$ between 1–8000 takes about 12 hours with a Pentium 4 processor. For statistics we performed 20 runs for plane partitions and 24 runs for solid partitions using different random number sequences. Since $`p_3(n)`$ is typically a very large number, the quantity that we keep track of in the simulations is $`\mathrm{ln}(p_3(n))`$.

### III.2 Simulation results for plane partitions

We first test the algorithm against the known case of plane partitions. In Fig. 3, we compare the simulation results with the exact answer. The two curves are almost indistinguishable. Consider the relative error defined by

$$\delta (n)=\frac{|\mathrm{ln}(p_2^{(s)}(n))\mathrm{ln}(p_2(n))|}{\mathrm{ln}(p_2(n))},$$
(20)
where $`p_2^{(s)}(n)`$ is the value obtained from simulations. We show the variation of $`\delta (n)`$ with $`n`$ in the inset of Fig. 3. The relative error goes to zero for large $`n`$. Thus, we conclude that the algorithm does give the correct leading asymptotic behaviour. In principle, the simulations can be made arbitrarily precise, and the correction terms can be determined. However, in our case, the statistical errors are not small enough to allow a reliable determination of the correction terms to the leading asymptotic behaviour.

### III.3 Simulation results for solid partitions

For solid partitions, we calculated $`p_3(n)`$ numerically by averaging over $`24`$ different runs. In Fig. 4, we show the results from simulation while in the inset of Fig. 4, we show the relative error. We estimate the asymptotic behaviour by fitting the data to an assumed form by the method of least squares fit. We fit $`\mathrm{ln}[p_3^{(s)}(n+1)/p_3^{(s)}(n)]`$ in the range 20–8000 to the form $`0.75\alpha _3n^{1/4}+0.5\beta _3n^{1/2}+0.25\gamma _3n^{3/4}+d_3n^1`$, where $`p_3^{(s)}(n)`$ is the value for solid partitions obtained from simulations. We choose this form since the MacMahon conjecture has the same functional form. As a test for the fitting routine, we test it against the simulation results as well as against the exact results for plane partitions using the fitting form $`0.667\alpha _2n^{1/3}+0.333\beta _2n^{2/3}+d_2n^1`$. The results for $`\alpha _{(2,3)}`$ are presented in Table. 2.

For plane partitions, the value of $`\alpha _2`$ obtained from the fitting routine is in excellent agreement with the exact answer 2.00945.. wright1 . Agreement with the correct answer is also seen for $`\alpha _3^{(m)}`$ (see Eq. (13)). Hence, we conclude that $`\alpha _3=1.79\pm 0.01`$.

## IV Summary and conclusions

In summary, we studied numerically the problem of solid partitions of an integer. Using exact enumeration methods, we extended the table of solid partitions for integers up to $`50`$. However, we were unable to determine the precise asymptotic behaviour of solid partitions from these $`50`$ numbers. Solid partitions for larger values of $`n`$ were studied using Monte Carlo simulations. From these simulations, we showed that $`lim_n\mathrm{}n^{3/4}\mathrm{ln}\left(p_3(n)\right)=1.79\pm 0.01`$. This value is consistent with the MacMahon value for solid partitions. Thus, if we assume that the asymptotic behaviour for solid partitions is correctly captured by a product form as in Eq. (11), then it should have the form $`_k(1q^k)^{(1/2\pm 0.012)k^2}`$.

## Acknowledgements

We would like to thank D. Dhar, D. B. Abraham and R. Stinchcombe for useful discussions. VM was partially supported by the Academy of Finland, Research Centre for Computational Science and Engineering, project No. 44897 (Finnish Centre of Excellence Programme 2000-2005), Wihuri Foundation, Finnish Cultural Foundation and Tekniikan edistämissäätiö. RR would like to acknowledge the financial support of EPSRC, UK.

## Appendix A Asymptotics for the MacMahon conjecture

In this appendix, we present a heuristic derivation of the asymptotic behaviour of the coefficient of $`q^n`$ in the expansion of the product

$$F(q)=\underset{k=1}{\overset{\mathrm{}}{}}(1q^k)^{a_1k^2a_2ka_3}.$$
(21)
Let $`q=e^ϵ`$. Taking logarithms on both sides of Eq. (21) and converting the resulting summation into an integral by using the Euler-Maclaurin summation formula (for example, see bender ), we obtain

$`\mathrm{ln}\left(F(e^ϵ)\right)`$ $`=`$ $`{\displaystyle \frac{2a_1\zeta (4)}{ϵ^3}}+{\displaystyle \frac{a_2\zeta (3)}{ϵ^2}}+{\displaystyle \frac{a_3\zeta (2)}{ϵ}}`$ (22)
$`+{\displaystyle \frac{a_2+6a_3}{12}}\mathrm{ln}(ϵ)+O(ϵ^0),`$
where $`\zeta (n)`$ is the Riemann zeta function. Let the coefficient of $`q^n`$ in $`F(q)`$ be denoted by $`c(n)`$. Then

$$c(n)=\frac{1}{2\pi i}\frac{F(q)}{q^{n+1}}.$$
(23)
For large $`n`$, we evaluate $`c(n)`$ by the method of steepest descent. The saddle point is the maximum of $`ϵn+\mathrm{ln}(F(ϵ))`$. This occurs at $`ϵ_0`$, where

$`ϵ_0`$ $`=`$ $`{\displaystyle \frac{a_1^{1/4}\pi }{15^{1/4}}}n^{1/4}+{\displaystyle \frac{\sqrt{15}a_2\zeta (3)}{2\sqrt{a_1}\pi ^2}}n^{1/2}`$ (24)
$`+`$ $`{\displaystyle \frac{15^{1/4}(a_1a_3\pi ^645a_2^2\zeta (3)^2)}{24a_1^{5/4}\pi ^5}}n^{3/4}+O(n^1).`$
Evaluating the integral about this saddle point, we obtain

$`\mathrm{ln}[c(n)]={\displaystyle \frac{4a_1^{1/4}\pi }{15^{1/4}3}}n^{3/4}`$ (25)
$`+{\displaystyle \frac{\sqrt{15}a_2\zeta (3)}{\sqrt{a_1}\pi ^2}}n^{1/2}+{\displaystyle \frac{5^{1/4}(a_1a_3\pi ^645a_2^2\zeta (3)^2)}{2a_1^{5/4}3^{3/4}\pi ^5}}n^{1/4}`$
$`\left({\displaystyle \frac{5}{8}}+{\displaystyle \frac{a_2}{48}}+{\displaystyle \frac{a_3}{8}}\right)\mathrm{ln}(n)+O(n^0).`$


text:
# Constraining photon-axion oscillations using quasar spectra

## 1 Introduction

One of the most direct and conceptually simple ways to probe the energy densities of the universe is through the redshift-distance relation. The change in the expansion rate caused by the different energy components can be measured by comparing recession velocities (i.e., redshifts) at different epochs of the universe (i.e., different cosmological distances). Type Ia supernovae (SNe) have proved to be very useful distance indicators since, after proper empirical corrections, they have small intrinsic dispersion and are luminous enough to be observed at very large distances. This method has been used by independent groups to show the existence of a dominant energy component with negative pressure, e.g., a cosmological constant . Future supernova (SN) surveys, like the CFHTLS and Supernova/Acceleration Probe (SNAP)  will increase the statistics radically, especially at large redshifts, thus making the minimization and control of different systematic effects crucial in order to be able to take full advantage of the increased statistics.

Concerns have been raised about, e.g., evolution of SN properties, gravitational lensing, dust absorption and photon-axion oscillations.

Brightness evolution of SNe can be studied by, e.g., comparing SN lightcurves and spectra at different redshifts and for different host galaxy types. Studies indicate that SNe have quite similar properties, even when formed in very different environments. Therefore cosmic evolution of SN properties or host galaxy extinction is not believed to affect the conclusion of the universe being dominated by a negative pressure component .

Systematic effects from gravitational lensing can be diagnosed by studying the induced asymmetry in the luminosity distribution whereas intergalactic dust absorption can be disclosed and corrected for from observed SN colours and spectrophotometry .

Photons oscillating into very light axions may cause distant SNe to appear dimmer . The net effect could be confused with the existence of a dark energy component with a negative equation of state. The presence of a significant component of dark energy has recently been independently inferred from other cosmological tests, most notable from recent WMAP background radiation measurements  in combination with limits on the Hubble parameter from the HST Key Project  and large scale structure measurements from the 2dFGRS survey . Therefore the proposed photon-axion mixing does not replace the need for dark energy, but its equation of state could be different.

The broad-band luminosity distribution and wavelength dependent attenuation caused by photon-axion oscillations have not been studied in any great detail so far. In this paper we examine to what degree it is possible, despite the complicated wavelength dependence of the effect, to put limits on the magnitude of the attenuation by studying a sample of quasar stellar object (QSO) spectra.

## 2 Photon-axion oscillations

Spin 1 particles, as the photon, can mix with zero spin bosons, as the hypothetical axion, in the presence of a mixing agent that allows for the conservation of quantum numbers during the process. The mixing agent could be a magnetic field transverse to the propagation direction of the photon.

The interaction between the photon and the axion is described by the Lagrangian

$$_{\mathrm{int}}=\frac{a}{M_\mathrm{a}}\stackrel{}{E}\stackrel{}{B},$$
(1)
where $`a`$ is the axion field, $`M_\mathrm{a}`$ is a mass scale determining the strength of the coupling and $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ are the electrical and magnetic field, respectively. The mass scale is given by

$$M_\mathrm{a}\frac{\pi f_\mathrm{a}}{\alpha }$$
(2)
where $`f_\mathrm{a}`$ is the decay constant of the axion. Besides the axion decay $`a2\gamma `$, this coupling allows for the conversion $`a\gamma `$ in an external electromagnetic field as depicted in Fig. 1.

Due to the interactions with gluons inducing transitions to neutral pion, axions acquire an effective mass, $`m_\mathrm{a}`$, given by

$$m_\mathrm{a}f_\mathrm{a}m_\pi f_\pi ,$$
(3)
where $`m_\pi =135`$ MeV is the pion mass and $`f_\pi 93`$ MeV its decay constant. Note however that the axion mass depends on the quark mass ratios as well as higher-order corrections which are not well known. In this paper, we assume that there exist an axion with $`m_\mathrm{a}f_\mathrm{a}m_\pi f_\pi `$.

The most important astrophysical limits on $`f_\mathrm{a}`$ are based on the requirement that the axionic energy loss of, e.g., globular cluster stars or the core of SN1987A, is not too efficient. A lower limit on the decay constant is given by $`f_\mathrm{a}10^9`$ GeV (corresponding to $`m_\mathrm{a}10^2`$ eV), indicating that axions, if they exist, are very light and very weakly interacting .

## 3 Simulation method

In order to predict the attenuation due to photon oscillations the magnetic field strength configuration and the electron plasma along the line of sight need to be modeled. Earlier work focused on refining the configuration model. However, the tools used to calculate the attenuation were too crude to adequately predict the wavelength dependence of the effect. The most common method used involves a calculation of the attenuation over a very short distance and then these small attenuation elements are added to get results over cosmological distances. Because this method does not to treat the oscillatory nature of the interactions, it fails to characterize its differential features along the wavelength axis. A more careful calculation shows that the effect can vary quite fast with wavelength, especially in the optical and IR regime. Thus, the simple approach will introduce errors larger than that caused by the uncertainties in the magnetic field strength and electron plasma density. As demonstrated in Ref. , a full density matrix calculation is necessary to get reliable results.

### 3.1 Magnetic fields and plasma densities

Only model dependent upper limits on the strength of intergalactic magnetic fields are available. This has to do with the fact that the most important observational technique used to trace far-away intergalactic magnetic fields – Faraday rotation – requires knowledge of the intergalactic electron density, $`n_\mathrm{e}`$ as well as the field configurations . The Faraday rotation is the effect of the polarization vector rotating when light travels through a magnetized medium. The polarization vector will rotate by an angle

$$\psi =\frac{e^3}{2\pi m_\mathrm{e}^2c^4}_0^ln_\mathrm{e}(l)B_{}(l)\left[\frac{\lambda (l)}{\lambda _{\mathrm{obs}}}\right]^2𝑑l,$$
(4)
where $`\lambda (l)`$ is the wavelength at position $`l`$ along the light path. Since neither the intergalactic electron density nor the magnetic field domain size is well known, we have chosen to work with a simple model where the domain sizes correspond to typical galaxy-galaxy separations with an average comoving density, $`n_\mathrm{e}`$, with a Gaussian dispersion of 50 % between domains. In accordance with Ref. , we assume that the magnetic field is frozen into the plasma and subsequently given by $`B_0n_\mathrm{e}^{2/3}`$ with random direction. For each model, we have tested whether our simulations are consistent with current Faraday rotation measurements. Since the magnetic domains are small compared to the total travel length and the magnetic fields have random directions, the Faraday rotation induced is too small to be used to rule out any of our investigated models.

### 3.2 Density matrix formalism

The equation to solve for the evolution of the density matrix $`\rho `$ is given by

$$\mathrm{i}\delta _\mathrm{t}\rho =\frac{1}{2\omega }[M,\rho ],$$
(5)
with initial conditions

$$\rho _0=\left(\begin{array}{ccc}\frac{1}{2}& 0& 0\\ 0& \frac{1}{2}& 0\\ 0& 0& 0\end{array}\right).$$
(6)
Here the three diagonal elements refer to two different polarization intensities and the axion intensity, respectively and,

$$M=\left(\begin{array}{ccc}\mathrm{\Delta }_{}& 0& \mathrm{\Delta }_\mathrm{M}\mathrm{cos}\alpha \\ 0& \mathrm{\Delta }_{}& \mathrm{\Delta }_\mathrm{M}\mathrm{sin}\alpha \\ \mathrm{\Delta }_\mathrm{M}\mathrm{cos}\alpha & \mathrm{\Delta }_\mathrm{M}\mathrm{sin}\alpha & \mathrm{\Delta }_\mathrm{m}\end{array}\right).$$
(7)
The quantities appearing in this matrix are given by

$`\mathrm{\Delta }_{}`$ $`=`$ $`3.6\times 10^{25}\left({\displaystyle \frac{\omega }{1\mathrm{eV}}}\right)^1\left({\displaystyle \frac{n_\mathrm{e}}{10^8\mathrm{cm}^3}}\right)\mathrm{cm}^1,`$
$`\mathrm{\Delta }_{}`$ $`=`$ $`\mathrm{\Delta }_{},`$
$`\mathrm{\Delta }_\mathrm{M}`$ $`=`$ $`2\times 10^{26}\left({\displaystyle \frac{B_{0,}}{10^9\mathrm{G}}}\right)\left({\displaystyle \frac{M_\mathrm{a}}{10^{11}\mathrm{GeV}}}\right)^1\mathrm{cm}^1,`$
$`\mathrm{\Delta }_\mathrm{m}`$ $`=`$ $`2.5\times 10^{28}\left({\displaystyle \frac{m_\mathrm{a}}{10^{16}\mathrm{eV}}}\right)^2\left({\displaystyle \frac{\omega }{1\mathrm{eV}}}\right)^1\mathrm{cm}^1,`$ (8)
where $`B_{0,}`$ is the strength of the magnetic field perpendicular to the direction of the photon, $`M_\mathrm{a}`$ is the inverse coupling between the photon and the axion, $`n_\mathrm{e}`$ is the electron density, $`m_\mathrm{a}`$ is the axion mass and $`\omega `$ is the energy of the photon. The angle $`\alpha `$ is the angle between a fixed polarization vector and the perpendicular (projected) magnetic field.

We solve the system of 9 coupled (complex) differential equations numerically111We have used the lsoda package from Netlib, http://www.netlib.org, by following individual light paths through a large number of cells where the strength of the magnetic field and the electron density is determined from predefined distributions and the direction of the magnetic field is random. Through each cell the background cosmology and the wavelength of the photon are updated, as are the matrices $`\rho `$ and $`M`$. In all of our simulations, we use a $`[\mathrm{\Omega }_\mathrm{M}=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,h=0.7]`$-cosmology. However, the cosmology dependence is quite weak.

### 3.3 Parameter dependence

In order to study the qualitative behavior of the solutions without any regard to the polarization state of the photons, we rewrite $`M`$ as a $`2\times 2`$ matrix,

$$M^{2\mathrm{D}}=\left(\begin{array}{ccc}\mathrm{\Delta }& \mathrm{\Delta }_\mathrm{M}& \\ \mathrm{\Delta }_\mathrm{M}& \mathrm{\Delta }_\mathrm{m}& \end{array}\right),$$
(9)
where $`\mathrm{\Delta }=\mathrm{\Delta }_{}=\mathrm{\Delta }_{}`$ and $`\mathrm{\Delta }_\mathrm{M}`$ is the component of the magnetic field parallel to the average polarization vector of the photon beam. We solve Eq. (5) for the density matrix $`\rho ^{2\mathrm{D}}`$ with initial conditions

$$\rho _0=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right),$$
(10)
where the diagonal elements refer to the photon and the axion intensity respectively. Assuming a homogeneous magnetic field and electron density, we can solve the two-dimensional system analytically. For the $`\rho _{11}^{2\mathrm{D}}`$ component, referring to the photon intensity, we get

$`\rho _{11}^{2\mathrm{D}}`$ $`=`$ $`1\left({\displaystyle \frac{\mathrm{\Delta }_\mathrm{M}}{2\omega \mathrm{\Omega }}}\right)^2(1\mathrm{cos}\mathrm{\Omega }t),`$
$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{\sqrt{(\mathrm{\Delta }\mathrm{\Delta }_\mathrm{m})^2+2\mathrm{\Delta }_\mathrm{M}^2}}{2\omega }}.`$ (11)
For values close to the typical set of parameter-values as indicated in Eq. (3.2), we can set $`\mathrm{\Omega }\mathrm{\Delta }/(2\omega )`$ to get

$$\rho _{11}^{2\mathrm{D}}1\left(\frac{\mathrm{\Delta }_\mathrm{M}}{\mathrm{\Delta }}\right)^2(1\mathrm{cos}\frac{\mathrm{\Delta }t}{2\omega }).$$
(12)
The oscillation length is of the order $``$ Mpc, i.e., comparable to the size of the domains used in the numerical integration. Note however that this number is very sensitive to the specific input parameter values used. In general smaller magnetic domains (corresponding to smaller cell sizes) yield lower oscillation probabilities, in accordance with the results of Ref. .

For $`m_\mathrm{a}m_{\mathrm{max}}38\sqrt{n_\mathrm{e}/(10^8\mathrm{cm}^3)}\mathrm{\hspace{0.17em}10}^{16}\mathrm{eV}`$, the oscillations are suppressed as $`m_\mathrm{a}^4`$. For $`m_\mathrm{a}m_{\mathrm{max}}`$, the effect is quite insensitive to the values of the axion mass which numerical simulations show is true even for $`m_\mathrm{a}m_{\mathrm{max}}`$ . Thus, in all simulations we set $`m_\mathrm{a}=10^{16}\mathrm{eV}`$.

Already from Eq. (3.2), it is seen that only the combination $`B_0/M_\mathrm{a}`$ affects the oscillation probability. Defining $`M_\mathrm{a}^{11}=\frac{M_\mathrm{a}}{10^{11}\mathrm{GeV}}`$, for small mixing angles, the effect is roughly proportional to $`(B_0/M_\mathrm{a}^{11})^2`$ whereas the effect is rather insensitive to the exact value of $`B_0/M_\mathrm{a}^{11}`$ in cases of close to maximal mixing. Note that even if the magnetic field strength $`B_0`$ would be well constrained from independent measurements, the relevant quantity $`B_0/M_\mathrm{a}^{11}`$ will not be, due to our lack of knowledge of $`M_\mathrm{a}^{11}`$.

Generally the effect is stronger for low values of the electron density. The value of $`n_\mathrm{e}`$ also has a strong impact on the wavelength dependence of the effect in the respect that higher plasma densities generally cause a more rapid wavelength dependence. In Fig. 2, we show the median QSO spectrum obtained in Ref.  unaffected by photon-axion oscillations (upper panel), with \[$`B_0=10^9,n_\mathrm{e}10^8\mathrm{cm}^3`$\] (middle panel) and \[$`B_0=10^9,n_\mathrm{e}10^9\mathrm{cm}^3`$\] (lower panel) for a redshift of $`z=1`$. For $`B_0=10^9`$ G and $`n_\mathrm{e}10^{10}\mathrm{cm}^3`$, we have an oscillation “wavelength” of $`1000`$ Å. For $`n_\mathrm{e}10^8\mathrm{cm}^3`$, we have variations over an interval of $`10100`$ Å whereas for $`n_\mathrm{e}10^7\mathrm{cm}^3`$ the variations occur over intervals of $`1`$ Å making them hard to discern from noise.

## 4 Analysis method

One of the main difficulties in constraining photon-axion oscillations is that the effect can look dramatically different for different set of parameter values, as shown in Fig. 2.

Given a set of very similar low-noise, high redshift spectra, we should nevertheless be able to put useful constraints on this kind of effect by studying the dispersion around the mean spectrum.

### 4.1 The sample

We have used a data sample of 3814 QSO spectra from the Sloan Digital Sky Survey (SDSS) Early Data Release (EDR) . The QSOs have redshifts between $`z=0.15`$ and $`5.03`$ distributed according to Fig. 3.

We divide the QSOs into redshift bins of size $`\mathrm{\Delta }z=0.2`$ in the interval $`0.1<z<2.9`$. Each spectrum is smoothed by binning the data into wavelengths bins of size 4 Å (in the restframe) and calculating the mean flux in the data points that fall within each bin. In each redshift bin, we calculate a mean restframe spectrum by normalizing the flux so that we have a mean flux of one (arbitrary units) over the common wavelength interval and average over the smoothed spectra. The spectrum-to-spectrum differential variation depends on the wavelength and redshift but is lower than 20 % over large intervals. In Fig. 4, we show the mean spectrum (upper panel) and spectrum-to-spectrum variation (lower panel) in the redshift bin $`1.1<z<1.3`$ consisting of 383 QSOs with a mean redshift of $`z=1.203`$.

### 4.2 Simulated spectra

Using the density matrix method described in Sec. 3.2, we have made a large number of Monte Carlo simulations of the effect on QSO spectra of photons oscillating into axions over cosmological distances. Our parameters are $`B_0/M_\mathrm{a}^{11}`$ and the plasma density $`n_\mathrm{e}`$. For each set of parameter values and redshift, we add the effect of photon-axion oscillations to the mean spectra derived from the real data set (e.g., for $`1.1<z<1.3`$, upper panel in Fig. 4).

### 4.3 Statistical analysis

In each redshift bin $`z`$, we calculate the dispersion (in flux), $`\sigma _{z,i}^r`$, around the mean for each spectrum ($`i`$). This is done in two steps, where after the first pass, spectral points more than 3 sigma off are rejected not to put excess weight on outliers caused by specific intrinsic features in the spectra. We thus obtain a distribution of $`\sigma _{z,i}^r`$ in each redshift bin, i.e., one number for each spectrum included in the redshift bin. The routine is repeated for the simulated data, with a set of parameters $`\stackrel{}{\theta }=(B_0/M_\mathrm{a}^{11},n_\mathrm{e})`$ giving a corresponding set of $`\sigma (\stackrel{}{\theta })_{z,i}`$. Since photon-axion oscillations only add to the intrinsic dispersion, a very conservative approach is to rule out any scenario for which $`\sigma _{z,i}^r<\sigma (\stackrel{}{\theta })_{z,i}`$. Comparing the cumulative distributions $`S(\sigma _{z,i}^r)`$ and $`S(\sigma (\stackrel{}{\theta })_{z,i})`$ using the Kolmogorov-Smirnov (K-S) test and only taking into consideration the maximum value of the difference $`S(\sigma (\stackrel{}{\theta })_{z,i})S(\sigma _{z,i}^r)`$ we can use standard statistics (see, e.g., ) to obtain a probability $`p_z`$ that $`\sigma _{z,i}^r`$ generally is larger than $`\sigma (\stackrel{}{\theta })_{z,i}`$. By multiplying results from each redshift bin, we obtain a final probability for the configuration vector

$$P(\stackrel{}{\theta })=\underset{z}{}p_z(\stackrel{}{\theta }),$$
(13)
that $`\sigma ^r>\sigma (\stackrel{}{\theta })`$ and can thus rule out a specific set of input parameter values if, e.g., $`P(B_0/M_\mathrm{a}^{11},n_\mathrm{e})<0.05`$.

We are also interested in what observational effects each set of parameter values have when integrating over broad band filters. We have investigated this effect by integrating the attenuation over the restframe B-band for Type Ia SNe. In Fig. 5, we have combined the results from this investigation with the results from the K-S test. The grey area indicates the allowed region in the $`[B_0/M_\mathrm{a}^{11},n_\mathrm{e}]`$-plane from the K-S test using a $`P<0.05`$ cut, i.e., roughly corresponding to a 95 % confidence level. Note however that this is a very conservative estimation since we have not taken into consideration the intrinsic variation of QSO spectra.

The four contour levels labeled \[0.2,0.1,0.05,0.02\] correspond to the amount of attenuation caused by photon-axion oscillations when integrating over the restframe B-band of a set of Type Ia SNe at a redshift of $`z=0.8`$. There is a band in the $`[B_0/M_\mathrm{a}^{11},n_\mathrm{e}]`$-plane where more than 0.05 mag attenuation results beyond what may be excluded from this analysis. For low values of the plasma density, the allowed attenuation could be as large as 0.2 mag. Also indicated (dashed line) is the baryon density inferred from recent WMAP measurements , indicating the maximum allowed value of the plasma density, $`n_\mathrm{e}2.710^7\mathrm{cm}^3`$. To test the effect of the data quality on the results, the analysis was also done only including spectra with S/N$`>`$10, with a somewhat larger exclusion region as a result, emphasizing the conservative nature of our limits. The small scale features in Fig. 5 are due to low resolution in the parameter grid.

Previous results indicate that low electron densities are required to get an attenuation that increases with redshift. However, since the notion of photon-axion oscillations being solely responsible for the observed faintness of high redshift SNe seems less and less probable when independent cosmological tests support the interpretation of a dominant negative pressure component, we do not restrict our interest to such scenarios. Rather, we are trying to constrain the possible systematic effect for future high precision cosmological tests, e.g., the Supernova/Acceleration Probe (SNAP)  aiming at a total error budget of $`\mathrm{\Delta }\mathrm{mag}<0.02`$. In Fig. 6 we used the SNOC Monte-Carlo simulation package  to show the amount of attenuation caused by photon-axion oscillations when integrating over the restframe B-band of a set of Type Ia SNe as a function of redshift for three sets of parameter values passing the K-S test. Note that the result for low plasma densities are very similar to the one for a negative pressure component and that the overall normalization can be set by varying the strength of the magnetic fields and/or the photon-axion coupling strength.

A possible concern is what implications photon-axion oscillations can have for the CMB observations. As noted by previous authors, assuming that the magnetic fields form at relatively low redshifts, the background radiation will be redshifted to such long wavelengths that the effect becomes negligible. This can be seen by studying Eq. (3.3) for small photon energies. In doing this, we can also see that for infrared (IR) energies, we do have a sizable effect for low plasma densities which can be used to put constraints complementary to what we have derived using only optical spectra, see Fig. 7.

### 4.4 Polarization

Due to the mixing between the different polarization states and the axion, polarized light might lose some of the polarization when travelling cosmological distances. In contrast to Faraday rotation, this effect can be large even when having a large number of magnetic domains with uncorrelated field strength and direction. Since the amount of depolarization depends on the values of $`B_0/M_\mathrm{a}^{11}`$ and $`n_\mathrm{e}`$, observations of high redshift polarized sources can be used to put limits on these parameter values. However, the highest degree of polarization observed for high redshift sources is $`20`$ % (see, e.g., Ref.  for a short compilation of observations) and by comparing this values with the results depicted in Fig. 8 it is clear that these observations can not be used to add any additional constraints on the magnetic field strength and/or plasma density.

## 5 Discussion

The lack of oscillatory dependence in the differential attenuation of high-$`z`$ QSO optical light allows us to set severe limits on the parameter space for photon-axion oscillations. However, the present very conservative analysis does not exclude that such effects may be affecting the measured flux of $`z1`$ objects at the $``$10 % level. From Fig. 6, it is clear that assuming a low electron plasma density ($`n_\mathrm{e}10^{10}\mathrm{cm}^3`$), a low mass axion with sufficient coupling strength can imitate the effect of a negative pressure component with respect to the observed luminosity of high redshift sources. Since there are independent cosmological tests supporting the notion of a dominant negative pressure component in the Universe today, such a scenario seems unlikely. However, to ensure the use of, e.g., high-$`z`$ Type Ia supernovae for precise estimates of cosmological parameters, this issue needs to be addressed further. This could be done by extending our present study of quasar spectra to also include IR information. In the future, instruments like JWST or SNAP may be able to rule out this scenario or provide independent evidence for the existence of a very light axion.

## Acknowledgments

The authors would like to thank Georg Raffelt for helpful comments. A.G. is a Royal Swedish Academy Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation.


text:
# Flow of thin liquid films on chemically structured substrates

## 1 Introduction

In recent years substantial efforts have been invested in miniaturizing chemical processes by building microfluidic systems. The “lab on a chip concept” integrates a great variety of chemical and physical processes into a single device  in a similar way as an integrated circuit incorporates many electronic devices into a single chip. These microfluidic devices do not only allow for cheap mass production but they can operate with much smaller quantities of reactants and reaction products than standard laboratory equipments. This is particularly important for rare and expensive substances such as certain biological substances and for toxic or explosive materials. Even though most available microfluidic devices today have micron sized channels further miniaturization is leading towards the nano-scale. Besides meeting technical challenges, new theoretical concepts are needed to understand the basic physical processes underlying this new technology. Whereas the ultimate limits for the miniaturization of electronic devices are set by quantum fluctuations, in a chemical chip these limits are determined by thermal fluctuations and can be explored by methods of classical statistical mechanics.

There are two main lines of development for microfluidic systems. The first one encompasses systems with closed channels. A common technique to produce these devices is to cast poly(dimethylsiloxane) (PDMS) over a topographically structured master, to peel off the polymer after curing, and to seal the resulting topographically structured material onto a flat surface. This way rather complicated devices including valves and pumps have been fabricated . However, closed channel systems have the disadvantage that they can be easily clogged by solute particles such as colloids or large bio-polymers.

The second type of systems are open with a free liquid-vapor interface and the fluid is not confined by physical but by chemical walls. The idea is that the liquid will be guided by lyophilic stripes on an otherwise lyophobic substrate. There are two sub-types of this technique: using a single chemically patterned substrate  (see also Fig. 1) or a chemically patterned slit pore , respectively. The substrate surfaces can be structured chemically by printing or photographical techniques. All the techniques are confined to two dimensions.

Liquids in contact with chemically structured substrates are not only of technological use but also enjoy high interest in basic research. Wetting on chemically structured substrates has been studied extensively on all length scales. For example adsorption on stripe patterns shows a rich morphology on the macro-scale  as well as on the nano-scale . In this context flow of liquids has also been studied to a large extend. Most of this research has been focused on homogeneous substrates, in particular on dynamic contact angles, dewetting of thin films, and flow of thin films in many circumstances. For a review see Ref. .

In the following we first describe briefly those equilibrium features of wetting on chemically patterned substrates which are most relevant for microfluidic applications. After discussing the characteristics of thin film flow on homogeneous substrates, we state the three main conceptual differences between micro-scale and nano-scale fluid transport in chemical channels. We conclude with formulating four design principles for nanofluidic devices.

## 2 Wetting on chemically patterned substrates

Equilibrium wetting phenomena on chemically patterned substrates have been analyzed theoretically in great detail on both the macroscopic and the microscopic scale. Microscopic theories, such as the successful density functional theory, take into account the finite range of inter-molecular attractions and short-ranged repulsions explicitly . Density functional theories do not only allow to study the order of wetting transitions and the equilibrium shape of the wetting film but also the detailed microscopic structure of the liquid in the vicinity of the substrate and at the liquid-vapor interface .

In macroscopic theories, however, the inter-molecular interactions are approximated by local descriptions. For a wetting film this means that the free energy of the film is described by a bulk term proportional to the volume of the fluid, a surface tension term proportional to the area of the liquid-vapor interface, an interface term proportional to the area of the liquid-substrate interface, and a line tension term proportional to the length of the three phase contact line between liquid, vapor, and substrate. In most studies the line tension term is neglected. Macroscopic theories have been used to describe the shape of droplets on homogeneous and structured substrates . The equilibrium droplet shape in chemically patterned slit-like pores has also been studied with the same technique .

It is a great challenge to describe the intermediate scale between the microscopic and the macroscopic one. In most cases it is impossible to obtain analytical results from microscopic theories and numerical calculations are prohibitive for large systems. Current research focuses on extending the scope of macroscopic theories down to the meso-scale by incorporating microscopic effects such as the wave-length dependence of the surface tension  and detailed properties of the three phase contact line .

## 3 Flow of thin liquid films

Flow of thin liquid films on homogeneous substrates has been studied extensively, in particular the motion of the three-phase contact line, dewetting of thin films, the stability of falling liquid films, and Marangoni flow. Apart from a few molecular dynamics simulations, in most studies the liquid flow is described in terms of meso-scale hydrodynamics. This means that the hydrodynamic equations are augmented with long-ranged liquid-substrate interactions and hydrodynamic slip of the liquid at the substrate (with slip lengths on the nanometer scale). For a review see Ref. . In a phase-field description even compressibility effects and the finite width of the liquid-vapor interface have been taken into account . The finite interface width in the phase field models is the only signature of thermal fluctuations taken into account in hydrodynamic thin film models up to now.

The main analytic tool in this approach is the so-called lubrication approximation, which is a small-gradient expansion for the film thickness and which leads to a fourth order in space, first order in time parabolic partial differential equation for the film thickness. This equation is also often referred to as the thin film equation and it has been successful in describing dewetting processes quantitatively . Flow over chemical substrate inhomogeneities has been studied in the context of dewetting of unstable films  but to little extent for the actual situations relevant to microfluidics and nanofluidics .

## 4 Conceptual differences in microfluidics and nanofluidics

In this section we state three important conceptual differences between microfluidics and nanofluidics.

### 4.1 Triple line and chemical steps

One of the most critical issues in open microfluidic systems is to keep the liquid in the desired areas such as channels, reactors, and reservoirs. On a macroscopic scale the liquid will stay on the lyophilic channels for low filling and the three phase contact line will lie on the channel area or it will be pinned at the chemical step. Spill-over onto the lyophobic areas occurs once the contact angle of the liquid at the chemical step exceeds the advancing contact angle $`\alpha _a`$ on the lyophobic area (which is in general larger than the equilibrium contact angle $`\alpha _e`$ due to surface defects).

On the nano-scale the situation is quite different. First, the concepts of a contact line and contact angle have to be revised. A sharp contact line is replaced by a smooth transition from a mesoscopic wetting film to the precursor film  which is only some few molecular diameters thick and spreads ahead of the main portion of the moving liquid. Moreover even an atomically sharp boundary between lyophilic and lyophobic areas on the substrate will lead to a smooth lateral variation of the interaction potential between the liquid particles and the substrate.

Thus the macroscopic and sharp criterion for a liquid staying on a chemical channel, namely whether the triple line crosses the channel boundary or not, becomes fuzzy at the nano-scale. Since there is always a certain amount of liquid on the lyophobic part of the substrate, one has to address the issue which fraction of the liquid is outside the channels rather than whether there is liquid outside the channels.

### 4.2 Interaction between neighboring channels

The question of spillage is of course closely related to the question of interaction between neighboring channels. Macroscopically (neglecting evaporation and recondensation) two neighboring channels will interact once the two liquid films merge. One has to keep in mind that in a macroscopic description an empty and a filled channel next to each other (Fig. 3(e)) or two filled but not interacting channels (Fig. 3(a)) can be metastable states. In the first case two equally filled channels and in the second case a liquid bridge can be the equilibrium state. How such configurations are affected by flow has not been investigated yet.

On the nano-scale, however, tails of the liquid from two neighboring channels can leak onto the lyophobic area between the channels. If the channels are too close these tails will overlap and thus form a bridge (see Fig. 3(b) and (d)). Liquid can flow through such a bridge and particles immersed in the fluid can diffuse through these bridges. Thus keeping the two flows separated is a question of time-scales and not a question that can be answered definitively.

Also the presence of a monomolecular precursor film on the lyophobic area can lead to an exchange of molecules (see Fig. 3(f)). Macroscopically this is certainly negligible but since the ratio between the film thickness in the channel and the thickness of the precursor film can be of the order of $`10`$, this can be significant at the nano-scale.

Fluctuations also become more important at small scales. In the absence of a precursor film on the area between the channels, a filled channel next to an empty one can be a metastable state as in the macroscopic case discussed above. But the energy barrier which has to be overcome to connect the two channels is much smaller than for macroscopic distances and much smaller fluctuations of the film height are needed.

Not only the liquid-substrate interactions but also the liquid-liquid interaction has a range of up to $`100`$ nm. Thus the direct interaction between liquid streams in parallel channels will influence the flow and also the thermal fluctuations.

### 4.3 Flow and the atomistic structure of liquids

The conclusions in the last two subsections are mostly based on quasi-static considerations. Transport mechanisms and dynamic properties have not been discussed. Experience tells that down to length scales of about $`1`$$`10`$ nm hydrodynamic theories provide a quite good description of liquid flow. However, even at the micron scale atomistic properties of the liquid show up via the slip length at the liquid-substrate interface and via the details of the regularization of the stress singularity at the moving triple line (see Fig. 4(a)).

In nanofluidic systems there is a window of length scales within which long-ranged intermolecular forces play a role and hydrodynamics is still applicable. This window is centered around film thicknesses of about $`100`$ nm. Below this length scale the atomistic structure of the liquid comes into play. Currently only molecular dynamics simulations explore this atomic length scale region (see for example Fig. 4(b)).

## 5 Design issues for nanofluidics

For liquid flow inside chemically patterned micron-sized channels several design issues have been addressed . Due to the conceptual differences between flow on the micron scale and flow on the nano-scale the answers given for the micron scale cannot be transferred directly to the nano-scale. However, the basic issues are the same: (1) How much liquid can a chemical channel contain before considerable spillage onto the lyophobic areas occurs? (2) How wide must the channels be to support flow? (3) How small can the radius of curvature of bends in the channel be? (4) What is the minimum distance between liquid streams below which they interact?

In micron scale channels these issues can be addressed by considering only surface tensions, contact angles, and line tensions. On the nano-scale the situation is more complex. In particular the details of the interaction between liquid and substrate, the interaction among the liquid molecules, the influence of the atomistic structure of the fluid on the transport properties, and thermal fluctuations have to be taken into account in addition.


text:
# An MT-Style Optical Package for Optical Data Transmission

## 1 Introduction

Optical fibers are now replacing traditional copper wires for signal transmission in the detectors for high energy and nuclear physics research. Optical fibers are more compact than copper wires, freeing up valuable detector space, and allow signal transmission over long distances with small attenuation. More importantly, the fibers eliminate electromagnetic interference and ground loops.

Commercially available VCSEL (Vertical Cavity Surface Emitting Laser) or PIN packages for transmitting and receiving optical signals are quite bulky. Detectors for high energy and nuclear physics research using optical packages usually need to custom design their packages due to space constraints. The problem is especially severe for the vertex detector located within the vicinity of the beam colliding region. In this paper, we present a compact MT-style optical package for housing VCSEL and PIN diodes. The diodes couple to the fibers in a commercial MT connector. This design simplifies the testing and assembly of the optical components because the MT connector with the long fibers attached can be remounted with ease while preserving good light coupling efficiency. In the next section, we present the design of the package. This is followed by a section on the prototype results and then a summary.

## 2 Optical Package Design

The optical package is a two-piece design as shown in Fig. 1. The U-shaped receptor contains two holes for the guide pins. Each hole has a diameter of $`700\pm 1\mu `$m and the distance between the two holes is $`4600\pm 3\mu `$m. The tolerances for the two dimensions are from the MT connector specification , hence the same tolerances should be used in the optical package fabrication. The receptor can be made of plastic with mold injection.

The base inside the receptor is fabricated using alumina. The two pieces are glued together with epoxy. On the side of the base facing the fibers, there are two gold traces for each VCSEL or PIN diode. One of the traces is for the cathode and the other for the anode. A diode is attached to one of the traces with conducting epoxy and a wire bond connects the other trace to the pad on the top of the diode. The two traces go over the corner to the other side of the base and connect via wire bonds to the optical transmitter or receiver chip placed in close proximity. The top traces are 1.5 mm in length, hence the chip is connected to the diode with small stray capacitance and inductance, minimizing degradation of the electrical signal. This is one of the advantages of this optical package design. The base should be placed at a location such that the distance between the diode and the fiber is $`200\mu `$m for good light coupling efficiency.

To fabricate the bases, alumina sheet is ground to the proper thickness of the bases and then cut into strips for deposition of three-dimensional traces . The trace width is 150 $`\mu `$m and the gap between traces is 75 $`\mu `$m. Most of the deposited traces ($``$ 94%) have good connectivity across the corner of the strips. Strips with a large number of traces of good connectivity are then diced into individual bases.

The main technical challenge in the fabrication of the package is the stringent alignment tolerance of the VCSEL with respect to the fiber location. Since the light from the VCSEL is highly collimated, the VCSEL must be placed radially within $`10\mu `$m of the center of the fiber for good light coupling efficiency. The distance between the VCSEL and the end of the fiber is much less critical and we choose 200 $`\mu `$m for the design shown in Fig. 1.

The alignment tolerance for the PIN diode is not stringent because the light sensitive area of the PIN is large compared with the silica core of the fiber. Consequently, we have demonstrated the principle of the design by prototyping the VCSEL optical package only.

For the placement of the VCSEL and PIN diodes, a small drop of conducting epoxy is deposited on a trace at the location of interest. A diode is then placed on top of the epoxy and aligned. Since the location of the fibers in an MT connector are defined precisely with respect to the two guide pins, one can align the diode with respect to the center of the two holes in the receptor under an optical comparator. Alternatively, one can mount an MT connector with fibers attached on the optical package and shine light into the other end of the fibers. With the MT connector pulled back slightly from the optical package, one can see a bright spot on the face of the diode. The diode should then be moved so that the most sensitive area is in the center of the beam spot.

There are many design variations depending on the requirements of an experiment. If only one or two channels are needed in an optical package, one can use the much smaller mini-MT connector. A typical MT connector has either 8 or 12 fibers in a row with 250 $`\mu `$m between the centers of adjacent fibers. The number of VCSEL/PIN diodes that can be placed on a base depends on the physical dimension of the diodes. In the example shown in Fig. 1, one can place six diodes, thus only half of the fibers are used. Alternatively, one can fabricate MT connectors with wider fiber spacing to accommodate the physical size of the diodes. There are MT connectors that can connect up to 72 fibers in a 6 $`\times `$ 12 matrix . By staggering the diodes instead of arranging them in a row to circumvent the diode physical size limitation, one can place more diodes in a package. If only a small number of diodes is used in a package, one can add a wide ground trace between two signal traces to act as a ground plane for shielding to reduce potential crosstalk. One can of course use a mixture of VCSEL and PIN diodes in the same package so that the package provides two-way communication. One can also use a common cathode VCSEL or PIN array in the package instead of individual diodes. In this case, the number of traces is reduced by half and a trace is added for the common cathode connection.

## 3 Prototype Results

We fabricated several prototype VCSEL optical packages in order to test the principle of the design. As stated in the previous section, we only demonstrated the principle on the VCSEL optical package and not on the PIN package, due to the more stringent VCSEL alignment tolerances. We used VCSELs manufactured by Truelight . The physical dimension of the VCSEL was $`280\mu \mathrm{m}\times 280\mu \mathrm{m}`$, hence we could only place six VCSELs on a base for use with a 12-fiber MT connector as shown in Fig. 1. In principle, we could fabricate the U-shaped receptor using mold injection. However, this would be somewhat time consuming and not cost effective for a prototype test. We therefore cut an MT connector into three pieces and machined out a pocket for the base in each piece. This automatically guaranteed that the guide pin hole diameter and separation were within the mechanical specification. In fact, this is a viable alternative to mold injection for any project requiring only a few hundred packages. Care must be exercised in the machining or else the plastic receptor could crack.

Figure 2 shows the coupled VCSEL power for a 10 mA current through the VCSEL for two prototype packages. All channels have good coupled power, 400 $`\mu `$W or above. The power coupling efficiency is shown in Fig. 3. All channels have good efficiency, 50% or higher. We have also tested the feasibility of replacing one of the VCSELs as may be needed in a production when a VCSEL is damaged or misaligned. Figure 4 shows the coupled optical power of all VCSELs before and after the repair. It is evident that the coupled power of the VCSELs after the repair are consistent with those before the operation. This demonstrates that a VCSEL can be replaced without damaging other channels.

The optical signal from the VCSEL is of high quality as shown in Fig. 5. Both the rise and fall times of the signal are fast, $`600`$ ps.

## 4 Summary

In summary, we have demonstrated the principle of an MT-style optical package. The VCSELs in the optical package have good coupled power (efficiency). The package is quite compact and individual VCSEL or PIN diodes in the package can be replaced if necessary. This package design simplifies the testing and assembly of the optical components because the MT connector with the long fibers attached can be remounted with ease while preserving good light coupling efficiency.

## 5 Acknowledgement

The author wishes to thank K. Arms, J. Burns, R.D. Kass, S. Smith, and R. Wells for their contributions.