text:
# Quantum State Tomography of Complex Multimode Fields using Array Detectors

## I Introduction

The subject of quantum state tomography (QST) has been of great interest in the recent years after the experimental reconstruction of the complete wave function of the vacuum and squeezed states by balanced homodyne (BHD) method. Much progress has been made by way of theoretical exploration of alternate schemes of QST like the self-homodyne tomography and their experimental verification. In essence, the BHD method combines a strong local oscillator (LO) with the signal in a 50:50 beam-splitter (BS). The two outputs of the BS are measured by two photodetectors and electronically subtracted. The resultant quantity is directly proportional to the rotated quadrature $`[\widehat{a}\mathrm{exp}(i\varphi )+\widehat{a}^{}\mathrm{exp}(i\varphi )]/\sqrt{2}`$ of the signal field and the angle $`\varphi `$ is fixed by the LO. By performing the experiment many times for a given $`\varphi `$ and repeating it for various values of $`\varphi `$ in the range of $`0\varphi \pi `$, the probability density $`p(x,\varphi )`$ for the quadrature rotated through $`\varphi `$ can be measured. From the measured probability density the Wigner function for the signal can be constructed and hence the elements of the density matrix. The serious problems that may arise in the numerical reconstruction of density matrix elements and the methods to avoid such pitfalls are discussed in the review article by Welsch et al  and in Ref. . The BHD can be extended  to multimode signals and to measure other distributions like the positive P-distribution. The BHD is for optical fields and methods have been developed for reconstructing the vibrational state of molecules and spin systems. For measruing the quantum states of the modes of a cavity, atoms with properly chosen energy levels can be used as probes .

A serious drawback of the BHD using single detectors is the requirement for mode matching the signal and LO mode functions. A quantitative measure of mode matching is the overlap between the signal and LO mode functions and the efficiency of the BHD scheme is directly proportional to it,

$$\mathrm{efficiency}U_{\mathrm{signal}}^{}(x)U_{\mathrm{LO}}(x)𝑑x.$$
(1)
Here $`U_{\mathrm{signal}}`$ and $`U_{\mathrm{LO}}`$ are the respective mode functions of the signal and LO. If there is perfect mode matching, i.e., the overlap integral is unity, the BHD is very efficient. In the case of signal and the LO modes being orthogonal to each other, the BHD simply fails to determine the quantum state of the signal.

An ingenious way of circumventing the problem of mode matching was suggested by Beck. The key idea is to replace the single detectors at the BS output ports by array detectors. Beck considers array detectors with pixels small enough such that the mode functions of the LO and signal over any pixel is constant. The output from a given pixel labeled $`(j,j^{})`$ of one of the array detectors is given by

$`N_1(j,j^{},\varphi )={\displaystyle \frac{\delta x\delta ycT}{2L}}[{\displaystyle \frac{|\beta |^2}{D_xD_y}}+{\displaystyle \underset{n,m}{}}\widehat{a}_n^{}\widehat{a}_mU_n^{}(j,j^{})U_m(j,j^{})`$ (2)
$`+{\displaystyle \frac{\beta }{\sqrt{D_xD_y}}}{\displaystyle \underset{m}{}}[\widehat{a}_m\mathrm{exp}(i\varphi )U_m(j,j^{})+H.c]].`$ (3)
Here the signal is taken to be multimode and the LO is assumed to be a single mode coherent state $`|\beta \mathrm{exp}^{i\varphi }`$ ($`\beta `$ real). The summation is over the various modes of the signal field. The function $`U_m(j,j^{})`$ is the value of the $`m`$th mode function of the signal over the pixel $`(j,j^{})`$. The operators $`\widehat{a}_m^{}`$ and $`\widehat{a}_m`$ are the creation and annihilation operators for the $`m`$th mode of the signal. The mode function for the LO is taken to be the constant $`\frac{1}{\sqrt{D_xD_y}}`$ over the entire detector surface. The dimension of each pixel is $`\delta x\times \delta y`$ and that of the detector is $`D_x\times D_y`$ (=$`N\delta x\times N\delta y`$). The constant $`c`$ is the speed of light vacuum, $`T`$ is the duration of counting and $`L`$ is the longitudinal quantization length. Similarly, the output of the corresponding pixel at the other output port of BS is given by

$`N_2(j,j^{},\varphi )`$ $`=`$ $`{\displaystyle \frac{\delta x\delta ycT}{2L}}[{\displaystyle \frac{|\beta |^2}{D_xD_y}}+{\displaystyle \underset{n,m}{}}\widehat{a}_n^{}\widehat{a}_mU_n^{}(j,j^{})U_m(j,j^{})`$ (5)
$`{\displaystyle \frac{\beta }{\sqrt{D_xD_y}}}{\displaystyle \underset{m}{}}[\widehat{a}_m\mathrm{exp}(i\varphi )U_m(j,j^{})+H.c]].`$
On taking the difference of $`N_1(j,j^{},\varphi )`$ and $`N_2(j,j^{},\varphi )`$, it is seen that

$$N_d(j,j^{},\varphi )=\frac{\delta x\delta y}{\sqrt{D_xD_y}}\beta \underset{m}{}[\widehat{a}_m\mathrm{exp}(i\varphi )U_m(j,j^{})+H.c],$$
(6)
where we have set $`L=cT`$. If one of the mode functions, say $`U_l`$, is real, then the above expression yields, on summing over all the pixel indices and using orthonormal relationship among the modes,

$$\underset{j,j^{}}{}N_d(j,j^{},\varphi )U_l(j,j^{})=\frac{1}{\sqrt{D_xD_y}}\beta [\widehat{a}_l\mathrm{exp}(i\varphi )+\widehat{a}_l^{}\mathrm{exp}(i\varphi )].$$
(7)
It is to be noted that the mode matching factor does not enter in the expression relating the measured photocurrent difference and the rotated quadrature. This is of great practical utility as mode matching is extremely difficult in an actual experiment. The above expression for quadrature works for any mode function that is real and hence the quantum state of all such modes can be determined. However, their joint state cannot be estimated by the method as proposed by Beck. But correlated quantum states are of prime importance for conceptual understanding of quantum theory as well as applications. Hence, it is all the more essential that the array detector scheme be extended to correlated multimode fields whose mode functions need not be real.

In the present work we extend the Beck’s method to two-mode correlated states by interposing a linear device which introduces a variable relative phase between the two modes. The stringent condition to have real mode functions is relaxed by making measurements with the LO phase varying from zero to $`3\frac{\pi }{2}`$ and combining the measurements to get the value of the quadratures roatated through zero to $`\pi `$ as required for state reconstruction. Further, we establish that one of the outputs of BS as measured by an array detector is sufficient to measure the quantum state. This paper is organized as follows. In Section II the extension of Beck’s work to two-mode correlated state is given. Section III of the paper contains the results of Monte Carlo simulations to demonstrate the power of the method to construct joint probability densities with complicated structures. In Section IV of the paper we explicitly show that for mode functions satisfying certain conditions, a single array detector is sufficient to measure the quantum state. We have provided an appendix to make transparent the arguments given Section II for realxing the requirement for real mode functions.

## II Joint quantum state of two-mode field

In the present section a method of estimating the quantum state of a correlated two mode signal is given. The basic requirement is to construct the joint probability density $`P(x_1,\varphi _1,x_2,\varphi _2)`$ for the two quadratures rotated through $`\varphi _1`$ and $`\varphi _2`$ respectively. It is, therefore, essential to rotate the two quadratures independently. If the signal modes are spatially seperable, they can be mixed with two independent LOs (usually made from a single source by using a BS) and measurements are carried out using two BHD arrangements. In the other case, the signal is passed through a medium which mixes the modes with a relative phase. The mixed modes are taken to be spatially seperable and one of the seperated modes is used in the input port of BS. This requires only one LO and one BHD arrangment. In both the cases, however, spatial seperation of modes is required at some stage and mode matching is essential. In the present section we discuss how to use the BHD with array detectors in the latter method. The mixing of modes is achieved by using a linear device which introduces a relative phase between the two modes. For instance, the linear device could be a medium in which the two modes interact via the Hamiltonian given by

$$H_{int}\widehat{a}_1^{}\widehat{a}_2\mathrm{exp}(i\theta )+H.c$$
(8)
The field operators $`(\widehat{a}_1^{},\widehat{a}_2^{})`$ at the output of the linear device can be written as

$`\widehat{a}_1^{}`$ $`=`$ $`\mathrm{cos}\nu \widehat{a}_1i\mathrm{sin}\nu \mathrm{exp}(i\theta )\widehat{a}_2,`$ (9)
$`\widehat{a}_2^{}`$ $`=`$ $`\mathrm{cos}\nu \widehat{a}_2i\mathrm{sin}\nu \mathrm{exp}(i\theta )\widehat{a}_1,`$ (10)
where $`\nu `$ is a constant determined by the length of the linear device, strength of interaction, etc. The device mixes the two modes of the signal. For later use we define the rotated quadratures for the transformed field:

$$\widehat{X}_k^{}(\varphi ,\theta ,\nu )=\frac{\widehat{a}_k^{}\mathrm{exp}(i\varphi )+\widehat{a}_k^{}\mathrm{exp}(i\varphi )}{\sqrt{2}},k=1,2$$
(11)
The positive frequency part of the electric field operator entering the signal port of BS is

$$\widehat{E}^{(+)^{}}(x,y)=\left[\frac{2\pi \mathrm{}\stackrel{~}{\omega }}{L}\right]^{1/2}\left[U_1(x,y)\widehat{a}_1^{}+U_2(x,y)\widehat{a}_2^{}\right].$$
(12)
The electric field operator depends on $`\varphi `$,$`\theta `$ and $`\nu `$ in addition to the spatial coordinates $`(x,y)`$. However, in the sequel we include the dependence explicitly in the expressions for difference counts and other derived quantities. With the field operators given by Eq.9 entering the signal port of BS and the LO phase fixed at $`\varphi `$, the difference count $`N_d(j,j^{},\varphi ,\theta ,\nu )`$ is

$`N_d(j,j^{},\varphi ,\theta ,\nu )`$ $`=`$ $`\beta {\displaystyle \frac{\delta x\delta y}{\sqrt{D_xD_y}}}[\mathrm{exp}(i\varphi )[U_1(j,j^{})\widehat{a}_1^{}+U_2(j,j^{})\widehat{a}_2^{}]`$ (14)
$`+\mathrm{exp}(i\varphi )[U_1^{}(j,j^{})\widehat{a}_1^{}+U_2^{}(j,j^{})\widehat{a}_2^{}]].`$
and the difference count with the LO phase rotated further by $`\pi /2`$ is

$`N_d(j,j^{},\varphi +{\displaystyle \frac{\pi }{2}},\theta ,\nu )`$ $`=`$ $`i\beta {\displaystyle \frac{\delta x\delta y}{\sqrt{D_xD_y}}}[\mathrm{exp}(i\varphi )[U_1^{}(j,j^{})\widehat{a}_1^{}+U_2^{}(j,j^{})\widehat{a}_2^{}]`$ (16)
$`\mathrm{exp}(i\varphi )[U_1(j,j^{})\widehat{a}_1^{}+U_2(j,j^{})\widehat{a}_2^{}]].`$
The quantities $`N_d(j,j^{},\varphi ,\theta ,\nu )`$ and $`N_d(j,j^{},\varphi +\frac{\pi }{2},\theta ,\nu )`$ are measured in the experiment. From these measured quantities we construct

$`R(j,j^{},\varphi ,\theta ,\nu )`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}\sqrt{{\displaystyle \frac{D_xD_y}{2}}}\{N_d(j,j^{},\varphi ,\theta ,\nu )`$ (18)
$`iN_d(j,j^{},\varphi +{\displaystyle \frac{\pi }{2}},\theta ,\nu )\},`$
$`=`$ $`{\displaystyle \frac{\delta x\delta y}{\sqrt{2}}}\mathrm{exp}(i\varphi )[\widehat{a}_1^{}U_1^{}(j,j^{})+\widehat{a}_2^{}U_2^{}(j,j^{})].`$ (19)
Multiplying $`R(j,j^{},\varphi ,\theta ,\nu )`$ by $`U_k(j,j^{})`$ $`(k=1,2)`$ and summing over $`(j,j^{})`$ yields

$`{\displaystyle \underset{j,j^{}}{}}R(j,j^{},\varphi ,\theta ,\nu )U_k(j,j^{})`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\mathrm{exp}(i\varphi )\widehat{a}_k^{}.`$ (20)
The quadrature $`\widehat{X}{}_{}{}^{}{}_{k}{}^{}(\varphi )`$ can now be written in terms $`R(j,j^{},\varphi ,\theta ,\nu )`$ and its conjugate as

$`X{}_{}{}^{}{}_{k}{}^{}(\varphi )={\displaystyle \underset{j,j^{}}{}}[R(j,j{}_{}{}^{},\varphi ,\theta ,\nu )U_k(j,j^{})+C.c].`$ (21)
The equations determine the mixed quadratures in terms of measured difference counts. Using Eqs. 9-10, the quadratures for the mixed mode can be written in terms of the signal quadratures as

$`\widehat{X}{}_{}{}^{}{}_{1}{}^{}(\varphi ,\theta ,\nu )`$ $`=`$ $`\mathrm{cos}\nu \widehat{X}_1(\varphi )+\mathrm{sin}\nu \widehat{X}_2(\varphi +\theta +{\displaystyle \frac{\pi }{2}}),`$ (22)
$`\widehat{X}{}_{}{}^{}{}_{2}{}^{}(\varphi ,\theta ,\nu )`$ $`=`$ $`\mathrm{cos}\nu \widehat{X}_2(\varphi )+\mathrm{sin}\nu \widehat{X}_1(\varphi \theta +{\displaystyle \frac{\pi }{2}}).`$ (23)
If measurements are carried out for two values of $`\nu `$, say $`\nu _1`$ and $`\nu _2`$, then the signal quadratures can be evaluated from Eqs.(22)-(23) and the resulting expressions are:

$`\widehat{X}_2(\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}(\nu _2\nu _1)}}[\mathrm{sin}\nu _1\widehat{X}{}_{}{}^{}{}_{2}{}^{}(\theta ,\varphi ,\nu _2)\mathrm{sin}\nu _2\widehat{X}{}_{}{}^{}{}_{2}{}^{}(\theta ,\varphi ,\nu _1)],`$ (24)
$`\widehat{X}_1(\varphi \theta +{\displaystyle \frac{\pi }{2}})`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}(\nu _1\nu _2)}}[\mathrm{cos}\nu _1\widehat{X}{}_{}{}^{}{}_{2}{}^{}(\theta ,\varphi ,\nu _2)\mathrm{cos}\nu _2\widehat{X}{}_{}{}^{}{}_{2}{}^{}(\theta ,\varphi ,\nu _1)],`$ (25)
$`\widehat{X}_1(\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}(\nu _2\nu _1)}}[\mathrm{sin}\nu _1\widehat{X}{}_{}{}^{}{}_{1}{}^{}(\theta ,\varphi ,\nu _2)\mathrm{sin}\nu _2\widehat{X}{}_{}{}^{}{}_{1}{}^{}(\theta ,\varphi ,\nu _1)],`$ (26)
$`\widehat{X}_2(\varphi +\theta +{\displaystyle \frac{\pi }{2}})`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}(\nu _1\nu _2)}}[\mathrm{cos}\nu _1\widehat{X}{}_{}{}^{}{}_{1}{}^{}(\theta ,\varphi ,\nu _2)\mathrm{cos}\nu _2\widehat{X}{}_{}{}^{}{}_{1}{}^{}(\theta ,\varphi ,\nu _1)].`$ (27)
The first two expressions yield the values of the quadratures $`\widehat{X_2}`$ and $`\widehat{X_1}`$ rotated through angles $`\varphi `$ and $`\varphi \theta +\frac{\pi }{2}`$ respectively. The last two yield the value of two quadratures rotated through $`\varphi _1=\varphi `$ and $`\varphi _2=\varphi +\theta +\frac{\pi }{2}`$. By changing the LO phase $`\varphi `$ and the phase $`\theta `$ in the interaction Hamiltonian, the quadratures can be measured to construct the joint probability distribution $`P(x_1,\varphi _1,x_2,\varphi _2)`$. Note that there is no need to assume that the mode functions are real. Of course, the difference count measurements are to be carried out over a range of $`03\pi /2`$ for the LO phase; in the case of real mode functions it is sufficient to measure the difference count when the LO phase is varied from $`0\pi `$. From the measured joint probability density, the Wigner function for the two-mode state can be constructed which, in turn, can be used to determine the density matrix. The measured joint probability can also be used to construct the elements of the density matrix in two-mode number state basis.

The scheme suggested here is similar to what is described as the method of generalized rotation in phase space. However, in our scheme there is no need to spatially separate the signal modes at the output of the linear device. This is possible as mode matching is not essential in the present scheme. In fact, not separating them is useful to reduce the number of experimental runs. In one run of the experiment, the pairs $`\{\widehat{X}_1(\varphi ),\widehat{X}_2(\varphi +\theta +\frac{\pi }{2})\}`$ and $`\{\widehat{X}_1(\varphi \theta +\frac{\pi }{2}),\widehat{X}_2(\varphi )\}`$ are simultaneously estimated.

## III Monte Carlo Simulations

In this section we study the experimental feasibility of the above method by Monte Carlo (MC) simulations with the signal in the two-mode Perelomov state. This state is generated from the two-mode vacuum $`|0,0`$ as follows:

$$|\zeta =\mathrm{exp}(\zeta \widehat{a}_1^{}\widehat{a}_2^{}\zeta ^{}\widehat{a}_1\widehat{a}_2)|0,0.$$
(28)
The number state expansion for the above state is given by

$$|\zeta =\frac{1}{\mathrm{cosh}r}\underset{n=0}{\overset{\mathrm{}}{}}[\mathrm{exp}(i\gamma )\mathrm{tanh}r]^n|n,n,$$
(29)
and the joint probability density is

$$p(x_1,\varphi _1,x_2,\varphi _2)=\frac{2}{\pi \mathrm{AB}}\mathrm{exp}\left[\frac{(x_1+x_2)^2}{\mathrm{A}}\frac{(x_1x_2)^2}{\mathrm{B}}\right].$$
(30)
The constants $`r`$, $`\mathrm{A}`$ and $`\mathrm{B}`$ are related to the squeeze parameter $`\zeta `$:

$`\zeta `$ $`=`$ $`r\mathrm{exp}(i\gamma ),`$ (31)
$`A`$ $`=`$ $`{\displaystyle \frac{|1+\mathrm{tanh}r\mathrm{exp}[i(\varphi _1+\varphi _2+\gamma )]|^2}{1\mathrm{tanh}^2r}}`$ (32)
$`B`$ $`=`$ $`{\displaystyle \frac{|1\mathrm{tanh}r\mathrm{exp}[i(\varphi _1+\varphi _2+\gamma )]|^2}{1\mathrm{tanh}^2r}};`$ (33)
Random numbers distributed according to the distribution given in Eq. 30 were generated from Gaussian distributed random numbers by von Neumann’s rejection method. These numbers were then used to generate the output of an actual experiment. In Fig. 1 the result of the simulation to reconstruct the joint probability density is given for $`\varphi _1=\pi /4`$, $`\varphi _2=\pi /2`$, $`\gamma =\pi /4`$ and $`r=1.0`$ along with its contour plot . For comparison, the theoretical distribution is also given. We find that to reconstruct the probability density with reasonable accuracy, as many as $`1.6\times 10^5`$ experimental runs would be required for each set of values of $`\varphi _1`$ and $`\varphi _2`$. In Fig. 2 we present the simulated results for identical values for the parameters except $`\varphi _2`$ which is set equal to $`\pi /4`$. In Fig. 1 the distribution is narrower compared to what is shown in Fig, 2. As expected, from the figures it is clear that the reconstruction is better, for a given number of experiments, if the probability density is narrower.

In order to see whether an experiment can capture the complicated spatial structures of probability density, the MC simulations were done for the signal in truncated Perelomov state. These are defined as

$$|c_1,c_2=c_1|0,0+c_2\mathrm{exp}(i\delta )|1,1.$$
(34)
The superposition coefficients $`c_1`$ and $`c_2`$ are taken to be real and they satisfy the normalization condition $`c_1^2+c_2^2=1`$. The joint probability density for the quadratures is

$`p(x_1,\varphi _1,x_2,\varphi _2)={\displaystyle \frac{\mathrm{exp}(x_1{}_{}{}^{2}x_2{}_{}{}^{2})}{\pi }}[c_1^2+4c_2^2x_1{}_{}{}^{2}x_{2}^{}{}_{}{}^{2}+4x_1x_2c_1c_2\mathrm{cos}(\varphi _1+\varphi _2+\delta )].`$ (35)
The result of the MC simulation for the truncated Perelomov state is depicted in Fig. 3. We have used $`c_1=c_2=1/\sqrt{2}`$, $`\delta =\pi /8`$, $`\varphi _1=\pi /4`$ and $`\varphi _2=\pi /4`$. It requires $`2\times 10^5`$ points to generate the pattern shown.

## IV  QST using a single array detector

In this section we describe how to measure the joint probability density of a two-mode state using an array detector in one of the output ports of BS. No measurements are made at the other output. In the discussions to follow the pixels are labeled by a single index instead of two indices. We first introduce some notations to simplify the expressions:

$`U_1(k)`$ $`:`$ $`\text{value of the mode function }U_1\text{ on the pixel labeled }k,`$ (36)
$`U_2(k)`$ $`:`$ $`\text{value of the mode function }U_2\text{ on the pixel labeled }k,`$ (37)
$`n_d(k)`$ $`:`$ difference count between the pixels labeled $`1`$ and $`k`$$`(k=2,3\mathrm{}9)`$ (38)
$`V^T`$ $`:`$ transpose of vector $`V`$ whose elements are (41)
$`\{\widehat{a}_1^{}\widehat{a}_1^{},\widehat{a}_2^{}\widehat{a}_2^{},\widehat{a}_1^{}\widehat{a}_2^{},\widehat{a}_1^{}\widehat{a}_2^{},\mathrm{exp}(i\varphi )\widehat{a}_1^{},\mathrm{exp}(i\varphi )\widehat{a}_1^{},`$
$`\mathrm{exp}(i\varphi )\widehat{a}_2^{},\mathrm{exp}(i\varphi )\widehat{a}_2^{}\}`$
$`M`$ $`:`$ an $`8\times 8`$ matrix whose $`(k1)`$th $`(k=2,\mathrm{..9})`$ row is given by (46)
$`\{U_1^2(1)U_1^2(k),U_2^2(1)U_2^2(k),U_1^{}(1)U_2(1)U_1^{}(k)U_2(k),`$
$`U_1(1)U_2^{}(1)U_1(k)U_2^{}(k),S[U_1(1)U_1(k)],`$
$`S[U_1^{}(1)U_1^{}(k)],S[U_2(1)U_2(k)],S[U_2^{}(1)U_1^{}(k)]\},`$
$`\text{where}S={\displaystyle \frac{\beta }{\sqrt{D_xD_y}}}`$
The expression for $`n_d(k)(k=2,3\mathrm{}9)`$ is

$`n_d(k)`$ $`=`$ Counts in pixel 1 - Counts in pixel $`k`$ (47)
$`=`$ $`{\displaystyle \frac{\delta x\delta y}{2}}[{\displaystyle \underset{n,m=1}{\overset{2}{}}}\widehat{a}_n^{}\widehat{a}_m^{}[U_n^{}(1)U_m(1)U_n^{}(k)U_m(k)]`$ (49)
$`+{\displaystyle \frac{\beta }{\sqrt{D_xD_y}}}{\displaystyle \underset{m}{}}[\mathrm{exp}(i\varphi )\widehat{a}_m^{}[U_m(1)U_m(k)]+H.c]]`$
$`=`$ $`{\displaystyle \frac{\delta x\delta y}{2}}{\displaystyle \underset{j=1}{\overset{8}{}}}M_{kj}V_j.`$ (50)
The elements of the matrix are various combinations of the mode functions $`U_1`$ and $`U_2`$ over a selected set of nine pixels of the array detector. If the elements of $`M`$ are such that the matrix is invertible, the elements of the vector $`V`$ can be determined from the equation

$$V_j=\underset{k=2}{\overset{9}{}}(M^1)_{j,k1}n_d(k)j=1,2\mathrm{}8.$$
(51)
Note that the vector $`V`$ has as its elements the creation and annihilation operators of the two modes and their quadratic forms. The above equation then implies that the elements of $`V`$ are determined from the measured quantities $`n_d(k)`$. It is evident that the mixed quadratures can be determined as

$`\widehat{X}{}_{}{}^{}{}_{1}{}^{}(\varphi ,\theta ,\nu )`$ $`=`$ $`{\displaystyle \frac{V_5+V_6}{\sqrt{2}}},`$ (52)
$`\widehat{X}{}_{}{}^{}{}_{2}{}^{}(\varphi ,\theta ,\nu )`$ $`=`$ $`{\displaystyle \frac{V_7+V_8}{\sqrt{2}}}.`$ (53)
Once the mixed quadratures are estimated, the signal mode qudratures are obtained using Eqs.22-23. It is clear that if we can choose nine pixels so that the matrix $`M`$ is invertible, then the joint probability density of two-mode states can be determined. Despite the fact that only one of the output ports of BS is used, the LO fluctuations are eliminated as in the case of BHD. Another interesting feature is that there is no need to assume that the mode functions are real and the LO phase need not be changed beyond the usual range $`0\pi `$. The method can be extended to more than two modes. But the requirements on the mode functions will become more stringent.

## V Summary

The use of array detectors in BHD eliminates the need for difficult-to-achieve mode-matching. For two-mode signals, a linear device to introduce a relative phase between two modes can be used in the signal port of the BS. By varying the LO phase from zero to $`\frac{3\pi }{2}`$ and the relative phase between the modes, the joint probability density for the quadratures of the two modes can be experimentally determined. The QST measurements can carried out with one array detector instead of detectors at both the output ports of BS. The method retains the advantage of BHD in the sense that the LO fluctuations do not affect the measurement. However, the success of the method depends of the shape of the mode functions.

## A QST of a single, complex mode signal

In this appendix we show how to measure the quantum state of a single, complex mode signal. This appendix is provided to make the arguments of Section II more transparent. In Beck’s method it is required to have the mode function real for determining its quantum state. If the mode function is not real, then it cannot be factored out in the RHS of Eq. (6) and hence the measured difference count cannot be related to the quadrature. Note that this problem does not arise in the case of the conventional BHD using single detectors. We set $`N_d(j,j^{},\varphi )`$ as $`N_d(\varphi )`$ and $`U(j,j^{})`$ as $`U`$ to avoid lengthy expressions. Hence, Eq. (6) is rewritten as

$$N_d(\varphi )=\frac{\delta x\delta y}{\sqrt{D_xD_y}}\beta [U\mathrm{exp}(i\varphi )\widehat{a}+U{}_{}{}^{}\mathrm{exp}(i\varphi )\widehat{a}^{}].$$
(A1)
With the LO phase rotated through $`\varphi +\frac{\pi }{2}`$, the difference count is

$$N_d(\varphi +\frac{\pi }{2})=i\frac{\delta x\delta y}{\sqrt{D_xD_y}}\beta [\widehat{a}^{}\mathrm{exp}(i\varphi )U{}_{}{}^{}U\mathrm{exp}(i\varphi )\widehat{a}].$$
(A2)
Comibining the expressions for $`N_d(\varphi )`$ and $`N_d(\varphi +\frac{\pi }{2})`$ we get

$$N_d(\varphi )iN_d(\varphi +\frac{\pi }{2})=U{}_{}{}^{}\mathrm{exp}(i\varphi )\widehat{a}^{}$$
(A3)
and

$$N_d(\varphi )+iN_d(\varphi +\frac{\pi }{2})=U\mathrm{exp}(i\varphi )\widehat{a}.$$
(A4)
Having measured the operators $`\mathrm{exp}(i\varphi )\widehat{a}`$ and its adjoint, the quadrature can be estimated at once as it is a linear combination of the two operators. To achieve this we multiply the last two equations by $`U`$ and $`U^{}`$ respectively and use the normalization relation $`\delta x\delta y_{j.j^{}}U^{}U=1`$ to yield

$$\widehat{X}(\varphi )=\frac{\sqrt{D_xD_y}}{2\sqrt{2}\beta }\underset{j,j^{}}{}[N_d(\varphi )(U+U^{})+iN_d(\varphi +\frac{\pi }{2})(U^{}U)].$$
(A5)
As asserted in the beginning of the appendix, there is no need to assume the mode functions to be real in the above expression. When $`U`$ is real the second term on the RHS vanishes and we recover, as expected, Eq.7.


format: text

text:
# The Far-Infrared/Radio Correlation in Nearby Abell Clusters

## 1 Introduction

The effect of the cluster environment on member galaxies is a well-established field of study. Butcher & Oemler (1978) studied the fraction of blue galaxies in clusters as a function of redshift, and thereby demonstrated evolution in the cluster environment. In the more local universe, studies have compared cluster environments to their field counterparts. The HI content of cluster spirals has been shown to be deficient by about an order of magnitude (Haynes et al., 1984; Cayette et al., 1990), and several studies have investigated star formation rates in clusters (Kennicutt & Kent, 1983; Kennicutt et al., 1984; Gavazzi et al., 1991; Moss & Whittle, 2000).

A natural extension of the studies on the effect of the cluster environment is that of the far-infrared (FIR)/radio correlation. Though identified much earlier (van der Kruit, 1973), the FIR/radio correlation was most solidly demonstrated with data from IRAS. Using the IRAS 60$`\mu `$m and 100$`\mu `$m flux densities to estimate the total emission between about 40$`\mu `$m and 120$`\mu `$m, the FIR flux of normal star forming galaxies is strongly correlated with their radio flux (most frequently measured at 20cm). The postulated cause of this correlation is the same underlying stellar population giving rise to two very different emission mechanisms (for a review see Condon, 1992). It is believed that the FIR emission is thermal and originates in dusty HII regions heated by massive stars. The radio emission is mainly non-thermal (synchrotron) but presumably arises as the same massive stars accelerate cosmic ray electrons in supernovae. If the cluster environment were to affect these two emission mechanisms in differing ways, it would be noticeable in a deviation from the standard FIR/radio correlation.

Several studies have investigated this or related effects. Some of the earlier work was done by Gavazzi & Jaffe (1986), who investigated the 1.4GHz (20 cm) radio emission of spiral galaxies in Abell 1367 and Abell 1656. They found that the radio emission of the cluster spirals was up to a factor of ten times greater than that observed in comparable field spirals. Noting that some of the galaxies had high H$`\alpha `$ equivalent widths, they suggested that interaction of cluster spirals with the intracluster medium (ICM) enhanced star formation. Both H$`\alpha `$ and FIR studies have largely refuted that cluster spirals exhibit enhanced star formation (e.g., Kennicutt et al., 1984; Bicay & Giovanelli, 1987), requiring some other explanation for the increased radio emission relative to field spirals. Gavazzi et al. (1991) suggested that while the star formation rates of the cluster galaxies may not be greatly effected, compression of the galaxies’ star forming gas could result from ram pressure with the ICM. Assuming pressure balance, the interstellar media of the cluster galaxies would reach $`p_{ram}=\rho _{ICM}v^2`$. This would increase the magnetic field strength and consequently the emitted radio power for the same star formation rate. Thus, tests of the ram pressure model focus on position in the clusters (as a proxy for ICM density) and velocity of the galaxies. In addition to position, cluster richness can serve as an indicator of density as richer cluster presumably have more significant ICMs.

Andersen & Owen (1995) compared the FIR/radio correlation in field, rich cluster, and poor cluster environments. As the ram pressure model would predict, the rich clusters had radio over-luminous spirals (i.e., lower FIR/radio ratios) while there was no statistical difference between poor clusters and the field. Another piece of evidence for the ram pressure compression model was noted in Gavazzi & Boselli (1999), which again argued that the spirals in richer nearby clusters (Abell 1367 and Abell 1656) were more radio luminous, by about a factor of five. Furthermore, they claimed that the radio over-luminous galaxies had higher velocities relative to the cluster systemic velocities, which would constitute further evidence for the ram pressure compression model.

In this paper we investigate the FIR/radio correlation in nearby Abell clusters using the recently completed radio galaxy catalog of Miller & Owen (2001a). This catalog provides a number of advantages for studies of this nature. Its large angular extent within each cluster allows radial trends in the FIR/radio correlation to be explored. If enhancement of radio emission is a consequence of interaction with the ICM, it should be more prevalent in centrally-located cluster galaxies. Furthermore, the large extent within each cluster provides a built-in field sample for comparison. This removes any luminosity bias which could result in studies with field samples located nearer than the cluster samples (e.g., Andersen & Owen, 1995). The velocity data in the catalog also minimize the effect of non-cluster galaxies seen in projection on the cluster, and allow testing of the radio properties as a function of line-of-sight velocity relative to the cluster systemic velocity. And finally, since the catalog includes 329 radio galaxies spread over 18 nearby clusters, it provides a large database from which more statistically-significant conclusions can be drawn.

In addition to the advantages of the Miller & Owen (2001a) catalog, additional observational data help solidify this project. We have reasonably high signal-to-noise optical spectra of roughly half of the radio galaxies, including most of the anamolous objects (i.e., those with unusually high radio emission relative to their FIR emission). These spectra are used to identify AGN in the radio galaxy sample. The fact that ellipticals tend to reside in more crowded environments (an observation known as the morphology-density relation; Dressler, 1980) and that such galaxies are more likely to host radio AGN means such identification is critical. Inclusion of these galaxies would produce a strong signal that centrally-located galaxies are radio over-luminous, but this would clearly be unrelated to any effect of environment on star formation. Since the physical scales probed by the optical spectroscopy could include an AGN as well as circumnuclear star formation, we have performed high resolution radio imaging of a subset of the sample. These radio data help constrain the net effect an AGN could have on the total radio emission from the galaxies.

The paper is organized as follows. In Section 2, we outline the data used. This includes a brief summary of the Miller & Owen (2001a) catalog and parameters for testing, plus optical spectroscopy and high resolution radio imaging methods. The statistical tests and their results are presented in Section 3, for the complete sample and defined subsamples. Section 4 discusses the results, including the implications for the ram pressure compression model and the potential sources of bias in this investigation. Summary conclusions are provided in Section 5. For all calculations which require them, we use $`H_o`$ = 75 km s<sup>-1</sup> and $`q_o=0.1`$.

## 2 Data

### 2.1 The Sample and Test Parameters

The sample is drawn from the nearby cluster radio galaxy catalog of Miller & Owen (2001a). This catalog was created using the NRAO VLA Sky Survey (NVSS; Condon et al., 1998) and covers eighteen nearby Abell clusters (ACO; Abell, Corwin, & Olwin, 1989). Advantages of the catalog are its large physical extent within each cluster – from the cluster cores to projected radial distances of 3$`h_{75}^1`$Mpc – and spectroscopic redshift confirmation of cluster members. The catalog also includes FIR flux densities from IRAS, taken either from the IRAS catalogs or determined through cross-scan addition. These flux densities, coupled with the NVSS 1.4GHz radio fluxes, are used to evaluate $`q`$, a parametrization of the FIR-radio correlation defined by Helou et al. (1985) as

$$q\mathrm{log}\left(\frac{FIR}{3.75\times 10^{12}\text{W}\text{m}^2}\right)\mathrm{log}\left(\frac{S_{1.4GHz}}{\text{W}\text{m}^2\text{Hz}^1}\right)$$
(1)
where FIR is defined as

$$\left(\frac{\text{FIR}}{\text{W}\text{m}^2}\right)1.26\times 10^{14}\left(\frac{2.58S_{60\mu \text{m}}+S_{100\mu \text{m}}}{\text{Jy}}\right).$$
(2)
Star forming galaxies - both normal and starburst - are tightly grouped around $`q=2.3`$, having $`\sigma 0.2`$.

To assess the importance of velocity in the radio emission of the cluster galaxies, we have used the velocity data in Miller & Owen (2001a) and Struble & Rood (1999) to calculate the line-of-sight deviation of each galaxy’s velocity from its cluster velocity. This is simply expressed as

$$\delta v=\frac{c|z_{galaxy}z_{cluster}|}{\sigma _{cluster}}.$$
(3)
### 2.2 Optical Spectroscopy: Star Formation or AGN?

One possible explanation for cluster galaxies with low $`q`$ values is that there is a contribution to the net radio flux from an AGN. While this would be an interesting result in its own right, we wish to test whether interaction with the ICM has increased the cluster spirals’ radio emission via some alteration in star formation environment. Consequently, we need to identify AGN and separate them from the star forming galaxies for subsequent analysis. This was done with optical spectroscopy whenever possible, using the line diagnostics of Veilleux & Osterbrock (1987). Specifically, we adopt the basic definitions of Ho (1996): star forming (\[OI\] $`\lambda 6300<0.08`$ H$`\alpha `$, \[NII\] $`\lambda 6584<0.6`$ H$`\alpha `$, \[SII\] $`\lambda \lambda 6717+6731<0.4`$ H$`\alpha `$) and AGN (\[OI\] $`\lambda 63000.08`$ H$`\alpha `$, \[NII\] $`\lambda 65840.6`$ H$`\alpha `$, \[SII\] $`\lambda \lambda 6717+67310.4`$ H$`\alpha `$). Note that our star forming galaxies correspond to Ho’s HII nuclei, and our AGN correspond to both Seyferts and LINERs in Ho’s classification. We have also included galaxies with the absorption-line spectrum indicative of an old stellar population in our AGN class.

The spectroscopy was carried out with medium resolution gratings at both Kitt Peak and Apache Point. Basic details of the spectroscopy and two-dimensional reductions may be found in Miller & Owen (2001a). Classification of galaxies as star forming or AGN are presented in Miller & Owen (2001b), although a brief summary is provided here. For each galaxy, a 2″ nuclear aperture and a 15″ aperture were extracted. These enabled classification of each galaxy based on its nucleus and an off-nuclear region (the difference of the two apertures). Corrections for stellar Balmer absorption were made based on the strength of the 4000Å break ($`D_{4000}`$, the ratio of average $`F_\nu `$ evaluated in the ranges 4050Å – 4250Å to 3750Å – 3950Å; Bruzual, 1983) using the model results of Poggianti & Barbaro (1997). Line ratio diagnostics using the corrected values are then analyzed, with the location and significance of each ratio relative to its corresponding dividing line between star forming and AGN determined. Star forming galaxies are those for which no single test produces an AGN result. This scheme is thereby conservative in that the galaxies identified as star forming are unambiguously classified as such.

The target galaxies were chosen to satisfy several science drivers, and therefore cannot be easily categorized as some subset of the complete sample. In general, the radio galaxies were observed if they either lacked public velocity measurements or had $`q`$ values which did not allow for simple classification as star forming or AGN (in practice, $`1.0q2.1`$). This has resulted in observation of all but fourteen galaxies for which $`1.0<q<2.0`$, eight of which are classified as E or S0 in NED. Galaxies with values of $`q`$ less than 1.0 are generally the powerful radio galaxies residing in elliptical hosts. For all observations, the nucleus of the target galaxy is placed in the slit.

This spectroscopic selection introduces a potential source of bias. It is expected that optical spectroscopy of all the sample radio galaxies would greatly increase the number of objects in the confirmed star forming sample, and slightly increase the number in the AGN sample. This is because in general, galaxies with normal values of $`q`$ were not observed and these are usually normal star forming galaxies. The net effect of this on the investigation of the FIR/radio correlation in cluster star forming galaxies is somewhat uncertain. Inclusion of the those star forming galaxies for which long-slit spectra were not obtained could improve the significance of the results through stronger number statistics. However, these unobserved galaxies will tend to have $`q2.3`$ independent of location within the clusters and their inclusion might in effect mask the postulated centrally-located low $`q`$ population.

The spectroscopic selection also may effect the results for the AGN. Since galaxies with $`q<1.00`$ were assumed to be powerful radio galaxies and few optical spectra were collected for these objects, the subset corresponding to AGN is not complete. In effect, this enhances the apparent correlation of FIR and radio fluxes for the AGN sample by removing those AGN with very high radio fluxes relative to their FIR fluxes. However, the AGN sample should include almost all of the Seyferts and LINERs in the complete sample.

### 2.3 High Resolution Radio Imaging

While the optical spectroscopy can identify AGN components of the cluster radio galaxies, it suffers two drawbacks for the current study. First, the $`2\mathrm{}`$ extraction aperture corresponds to a linear resolution of about one kpc for the average cluster galaxy. The scale of a circumnuclear starburst can be significantly less than this figure, so if the emission line contribution due to an AGN were relatively weak compared to a concurrent nuclear starburst the AGN would, in effect, be missed. As such AGN could explain the enhancement of radio luminosity implied in the low $`q`$ galaxies, another test is necessary. Second, should AGN be present in low $`q`$ radio galaxies it would be useful to place limits on their net contribution to the radio luminosities.

To achieve these goals, we have conducted a program of high resolution radio imaging on a sample of radio galaxies in Abell 1367 and Abell 1656. Using the VLA in its A configuration at a frequency of 8.46GHz, we obtain a resolution of roughly $`0.25\mathrm{}`$ or approximately 100 pc. In general, the observational goals and experimental design were the same as those used by Condon et al. (1991) in their high-resolution radio survey of ultraluminous IR galaxies. The resolution of the observations implies that the radio contribution from a compact AGN would be unresolved. Condon (1992) finds that the most luminous starbursts have compact radio emission with $`D200`$ pc, so in general these radio observations would resolve radio emission from star formation in all but the most compact nuclear starbursts.

The observations were performed during a 12-hour track of time in June 1999. Even at frequencies as low as 8.46GHz, temporal changes in atmospheric water vapor can lead to significant phase variations on the longer baselines. This is particularly true during summer days. Consequently, the observations were performed using the technique of fast switching (see discussion in Carilli, Holdaway, & Sowinski, 1996) which relies on ‘nodding’ the VLA dishes between the target field and a nearby phase calibrator on time scales short enough to calibrate out the atmospheric phase variations. We adopted a net cycle time of 100 seconds, with 70 seconds devoted to the target field and 30 seconds devoted to the phase calibrator, which was typically within five degrees of the target field. Note that these times include move time for the dishes, so the net time on source is around 20 seconds less per cycle. Flux calibration was performed via observation of 3C286.

A total of 20 radio galaxies – 10 from each cluster – were observed. The primary science targets were identified radio galaxies with low $`q`$ values, though additional comparison sources were also observed. These comparison sources included several normal star forming galaxies with $`q`$ values near the canonical value of 2.3, and several known AGN with low $`q`$ values. The results will be presented in the next section.

## 3 Analysis

### 3.1 Statistical Results

As can be seen in Figure 1a, the FIR/radio correlation holds very well in the nearby Abell clusters of the sample. The $`q`$ values are mainly clustered around the expected value of 2.3 (see Condon, 1992). Perhaps the most remarkable aspect of this is that the full sample contains a large fraction of AGN and powerful radio galaxies. The results specific to various categories of galaxies will be discussed below.

The FIR/radio correlation, as parametrized by $`q`$, was tested for linear correlation with both radial separation from the cluster core and relative velocity. The purpose of these tests was to assess whether $`q`$ values of galaxies were a function of location or velocity within clusters. Because the area surveyed in each cluster reaches from the core out to large radii more consistent with the field, a significant correlation between $`q`$ and $`r`$ would suggest that the cluster environment did effect the FIR/radio correlation. Conversely, no correlation would suggest that any deviations in $`q`$ were normal fluctuations not caused by environment. In the case of velocity, we wish to test whether higher velocity objects are more likely to have enhanced radio emission, as would be expected in a simple model wherein ram pressure compressed the galaxies’ magnetic fields.

The statistical test used was Kendall’s Tau, with censoring of the $`q`$ values when they represented upper limits. This test examines the null hypothesis that no correlation is present between the two variables being tested. For this and all subsequent tests the IRAF ‘statistics’ package, which implements the “Astronomical Survival Analysis” programs (see Feigelson & Nelson, 1985; Lavalley, Isobe, & Feigelson, 1992), was used. Upper limits on $`q`$ values resulted from cases where the galaxy was not detected at 60$`\mu `$m, 100$`\mu `$m, or both and a value of three times the rms noise was substituted for the undetected flux(es). In a few cases, lower limits for $`q`$ were applied. These resulted when the low resolution of the NVSS blended the radio emission ascribed to one galaxy with another. The tests for linear correlations were performed on the sample as a whole and subsamples designed to isolate the star forming galaxies. Since we only possess optical spectroscopy for about half of the sample, we also created a subsample on the basis of radio powers alone. Based on the local RLF of Condon (1988), star forming galaxies are more common than AGN for $`\mathrm{log}(P_{1.4GHz})22.7`$, so we have used this figure as a rough separator of star forming galaxies and AGN. Note that this subsample does contain AGN, many of them identified in the optical spectroscopy (for more discussion on this point, see Miller & Owen, 2001b). However, it provides the benefits of a simple, objective classification for galaxies and a large sample size for statistical testing. Table 1 presents the results for each of the four samples (the complete sample, the subsample for $`\mathrm{log}(P_{1.4GHz})22.7`$, the spectroscopically-confirmed star forming galaxies, and the spectroscopically-confirmed AGN). Graphical representations of the data are provided in Figure 1 and Figure 2 for $`q`$ vs. $`r`$ and $`q`$ vs. $`\delta v`$, respectively. It should also be noted that the FIR and radio fluxes for the AGN are also well correlated. As can be seen in Figure 1, the scatter in this correlation is larger than it is for the star forming galaxies and the radio fluxes are in excess of those predicted by the FIR/radio correlation for normal star forming galaxies. This will be discussed further in Section 4.

These linear correlation tests produced evidence that radio over-luminous galaxies were centrally located in clusters, but not that they were higher velocity objects. For the complete sample (Figure 1a), the correlation of $`q`$ and $`r`$ was highly significant. This is easily explained as the complete sample consists of all radio galaxies in the clusters, including the powerful radio galaxies. These are almost always ellipticals, so the strength of correlation is simply a consequence of the morphology-density relation. The correlation is also quite strong for the subsample defined by radio power (Figure 1b). Most interestingly, the spectroscopically-confirmed star forming galaxies show a possible correlation of $`q`$ with $`r`$ ($`91.2\%`$; see Figure 1c). This is the result to be expected if the ICM affects these galaxies, and is explored in greater detail in subsequent tests. The equivalent plots for $`q`$ and velocity are shown in Figure 2a–2d. None of the tested samples showed significant evidence for a correlation of $`q`$ with velocity (see Table 1).

Given the large radial extent of the survey, a linear correlation of $`q`$ and $`r`$ could easily be washed out. This is compounded by projection effects, as we can reasonably expect the galaxies at low apparent radial separation to include a fraction of galaxies with much greater true cluster-centric separations. This population could easily overwhelm the contribution of the real cluster core members. In order to test whether the $`q`$ values of centrally-located cluster galaxies differed from the more field-like cluster galaxies at the periphery of the clusters, we grouped the data radially and compared the $`q`$ value distributions. Thus, within each sample we compared the distributions of $`q`$ values for galaxies within one Mpc of the cluster cores (in projection) with those located between one and two Mpc, and with those between two and three Mpc. The statistical tests used were the Kolmogorov-Smirnoff (KS) test and several variants on the Wilcoxon test. The KS test evaluates the maximum difference between the cumulative distributions of two samples, and tests the null hypothesis that the samples are drawn from the same parent distribution. The Wilcoxon test ranks the data of two samples, and tests whether the samples are drawn from two separate distributions having different means. Application of the KS test used no censoring of the data (i.e., information on upper and lower limits for $`q`$ was ignored) while the Wilcoxon tests applied several weighting schemes for censored data (see Feigelson & Nelson, 1985). Table 2 presents the results for these tests. For the Wilcoxon tests, three types of censoring are presented – Gehan permutation variance, logrank, and Peto & Prentice – though the results are nearly identical. Figure 3 has the histograms of the $`q`$ distributions for each subsample and radial bin.

The results of the distributional tests are consistent with the tests of linear correlation. Each of the three bins differed significantly from the others in the analysis of the complete sample, with the strongest difference being that between the inner and outer Mpc (at a significance level much greater than $`99.9\%`$; see Figure 3a). By the outer regions of the cluster the populations are more similar, with the middle and outer Mpc bins differing at only $`90\%`$ significance. As noted before, this is largely a result of the centrally-concentrated radio powerful ellipticals. The galaxies of the spectroscopically-confirmed star forming subsample also show a significant difference between the inner and outer Mpc bins (Figure 3c). The Wilcoxon tests place the significance level at $`97\%`$, with the mean $`q`$ value of the inner bin being less than that of the outer bin. Furthermore, the KS test of these two groups suggests that they differ at greater than $`99\%`$ significance. The distributions of the AGN also seem to show some importance for the cluster environment. The Wilcoxon tests show that the middle and outer bins are most disparate (see Figure 3d), the result of the middle bin having the lowest mean $`q`$ value and the outer bin having the highest. As can be seen from Figure 1d, by the outer Mpc the AGN have $`q`$ values much closer to those expected for normal star forming galaxies.

The possible dependence of $`q`$ on velocity was also tested further. Since only the line-of-sight component of the velocity can be determined from the redshift data (analagous to the problems encountered when examining radial separations), distributional tests were performed. To test whether the low $`q`$ star forming galaxies have higher velocities than the normal $`q`$ objects, we arbitrarily divided the spectroscopically-confirmed star forming sample into those for which $`q<2.00`$ and those for which $`q2.00`$; a dividing line where radio emission is a factor of two greater than the canonical value for the FIR/radio correlation. These two subsamples were then compared via the KS and Wilcoxon tests, and can be seen graphically in Figure 4. There is at best marginal evidence for the low $`q`$ galaxies to have higher relative velocities - the KS test indicates a $`64\%`$ probability that the two samples are drawn from different velocity distributions whereas the Wilcoxon tests are significant at about the $`90\%`$ level.

### 3.2 Results of High Resolution Radio Imaging

The 8.46GHz observations of radio galaxies in Abell 1367 and Abell 1656 are summarized in Table 2. Only five of the twenty targets were detected, and these correspond to five AGN designated as a comparison sample. Consistent with the theory that their radio emission originates predominately from a compact active nucleus, four of these five galaxies were unresolved and the fifth consisted of a strong unresolved component within a weaker, resolved emission complex. The implied brightness temperatures of these sources also point weakly to their being AGN ($`T_b10^4K`$). The five star forming galaxies for which $`q2.3`$ were not detected, nor were any of the ten program galaxies.

The sample selection implies that all of the objects will have radio emission at 8.46GHz, so the lack of detections are the result of the beam size and rms noise for each observation. Two important details can be obtained from this knowledge. First, an estimate of the maximum contribution to the 1.4GHz emission caused by an AGN can be determined. This is simply obtained from the rms noise per beam in the 8.46GHz map. Second, we can determine a minimum physical size over which the radio emission originates by assuming it is spread uniformly. In total, these results argue that the radio emission of the low $`q`$ galaxies does not arise from either an AGN or a compact nuclear starburst.

The non-detections demonstrate that the low $`q`$ values are not the result of radio contributions from active nuclei. Adopting a spectral index of 0.8 (i.e., $`S_\nu \nu ^{0.8}`$) we can calculate a 3$`\sigma `$ upper limit to the 1.4GHz contribution of an AGN. Based on the rms noise results for each observation (see Table 3), this ranges from 0.5 to 1.0 mJy. Removing this from the NVSS flux for each galaxy has a minimal effect on $`q`$, and even for sources at the NVSS detection limit (about 2.5 mJy) can not explain the factor of 2–3 amplification in radio fluxes implied by the lower $`q`$ values. For the worst case scenario of a galaxy just detected by the NVSS at 2.5 mJy, the difference in $`q`$ caused by 1.0 mJy of the total flux arising from a compact core is only $`0.2`$. For the strong star forming galaxies such as 114349+195811, the difference is very nearly zero ($`0.01`$). Note that the assumption of a spectral index of 0.8 is probably conservative in this regard. Should the spectral index be flatter, the net contribution of the AGN to the overall 1.4GHz flux would be even less.

Interestingly, two of the test objects exhibit AGN spectra but were not detected by the 8.46GHz observations. The optical spectrum of each of these galaxies is that of a weak LINER, dominated by an old stellar population but with weak emission of \[NII\] and \[SII\]. Both are resolved by the NVSS, suggesting that the radio emission arises over significantly larger scales than those probed by the 8.46GHz observations. The radio emission could arise from a low surface brightness jet, and is also consistent with the theory in which shock ionization over fairly large scales is the emission mechanism in LINERs (Veilleux et al., 1995).

The non-detections also place interesting limits on the size scale over which star formation is occuring in the galaxies. If the 8.46GHz emission were relatively concentrated, it would be detectable. In fact, this is the case for the one comparison sample object with a small resolved component. Assuming the 8.46GHz emission is spread uniformly, we can calculate the size at which we would just fail to detect the emission. Setting our detection threshold at 3$`\sigma `$, this corresponds to regions of about 400 pc diameter for the weaker sources and nearly 1.5 kpc for the strongest sources. Given the size of a compact starburst ($``$200 pc) and the propagation distance of the cosmic ray electrons presumed to be accelerated in supernovae within the star forming regions ($``$1 kpc), we can largely rule out that the star formation occurs in a compact circumnuclear starburst (Condon, 1992).

## 4 Discussion

### 4.1 Statistical Results and Implications for Models

In all, 24 of the galaxies in the spectroscopically-confirmed star forming sample had $`q<2.00`$. In three of these galaxies, the $`q`$ value is a lower limit where the net radio flux arises from multiple blended sources. Of the remaining low $`q`$ galaxies, about half had $`1.90q<2.00`$ (10 galaxies, though two represent upper limits). Thus, the population of galaxies with the enhanced radio emission implied by large deviations from the FIR/radio correlation is fairly rare. It should also be noted that the majority of the low $`q`$ galaxies come from richer clusters. Although seven of the low $`q`$ galaxies are members of Abell Richness Class 0 clusters (Abell 262, Abell 347, Abell 397, and Abell 569), the velocity dispersions of these clusters are generally around 600 km s<sup>-1</sup>, suggesting that they are fairly massive. The cluster richest in galaxies for which $`q<2.00`$ is Abell 1656 (Coma), with seven such galaxies. While this hints at the importance of a richer ICM we note that Abell 426 (Perseus), a Richness Class 2 cluster like Coma, has no low $`q`$ galaxies. This cluster is thought to be dynamically relaxed (although see Dupke & Bregman, 2000) and possesses a very low fraction of spiral galaxies.

The AGN prove interesting on two counts. First, as noted previously their FIR and radio fluxes are also strongly correlated, as has also been found in other studies (e.g., studies of Seyferts by Baum et al., 1993; Roy et al., 1998). This seems to be evidence for star formation in such galaxies. In fact, Baum et al. (1993) were able to recover the standard FIR-radio correlation by subtracting off the radio emission from the central kpc of Seyfert galaxies, suggesting normal star formation in their disks. Second, there is evidence for a change in the FIR-radio correlation of AGN within clusters. This may be understood if some aspect of the cluster environment affects the radio emission associated with star formation, producing a radial effect in correlations for both purely star forming galaxies and for AGN with concurrent star formation.

Our spectroscopic classification of galaxies underscores a potential pitfall in studies of the radio emission from cluster spirals. Figure 3d and Table 2 demonstrate that AGN also are affected by the cluster environment, with the $`q`$ values of AGN in the inner regions of the clusters being lower than those in the outer regions of the clusters. Removing the AGN which have spectra dominated by an old stellar population (i.e., creating a subsample of only emission-line AGN including both LINERs and Seyferts) does not alter this conclusion. Performing the Wilcoxon tests on the emission-line AGN within 2 Mpc (about one Abell radius) vs. those outside 2 Mpc, we find the $`q`$ value distributions differ at greater than $`98\%`$ significance. Studies which rely on imaging to select star forming galaxies are thereby prone to overestimate the effect of the cluster environment on such types. These studies frequently use wide-band filter images to select spiral galaxies, and follow up narrow-band H$`\alpha `$ imaging to determine which of the spirals are forming stars. The H$`\alpha `$ images can not differentiate between H$`\alpha `$ and \[NII\], and it is likely that some of the supposed star forming galaxies are actually AGN. Inclusion of such galaxies will produce a false signal of the effect of the cluster environment on the radio emission of cluster star forming galaxies.

The correlation and distributional tests consistently suggest the importance of location in a cluster, but not velocity. This is perhaps best exemplified in the distributional tests for the star forming galaxies (see Table 2 and Figure 3c). Galaxies seen in projection within the inner Mpc of the clusters differ from those over 2 Mpc from the cluster center at a very high significance level. This is strong evidence that local cluster environment has an effect on star forming galaxies, especially when one considers that the effect of projection is to lessen the magnitude of such an effect. But comparing the velocities of these low $`q`$ galaxies with their normal counterparts does not produce significant results. This mixed evidence for the ram pressure model will be discussed below.

The total external pressure on a galaxy due to interaction with its environment is the sum of pressures due to specific effects. For simplicity, we will include only two such pressures: ram pressure and the thermal pressure due to the surrounding ICM. The sum of these pressures may then be compared to a representative ISM pressure of the galaxies to determine the relative importance for external pressure and estimate how much compression could occur. Thus,

$$P_{external}=\rho v^2+\frac{\rho kT}{\mu m_p}$$
(4)
We can simply evaluate where compression of the ISM could be occuring by making some assumptions. The above equation is a function of only ICM density, temperature, and the relative velocity of a galaxy to the ICM. Each of these can be estimated for every galaxy in the sample. To evaluate the local ICM density, we use the results of models based on X-ray data. Briel et al. (1992) investigated the Coma cluster using ROSAT data, and found a reasonable fit to the cluster gas distribution by a King profile (King, 1966):

$$\rho (r)=\frac{\rho _o}{\left[1+(r/r_c)^2\right]^{3/2}}.$$
(5)
Using this formula, we estimate an upper limit to the local density simply by the projected radial separation from the cluster core. We adopt $`\rho _o=6.5\times 10^{27}g/cm^3`$ and $`r_c=0.2Mpc`$ (see Briel et al., 1992). In addition, the temperature is assumed to be 8 keV (about the temperature of Coma; see David et al., 1993). Given these parameters and assuming $`\mu =1`$ (i.e., the ICM is purely hydrogen), we can solve Equation 4 for the velocity at which the external pressure equals set increments of the representative ISM pressure for the galaxies. These increased pressures would presumably compress the galaxies’ ISM, thereby increasing the magnetic field strength and emitted synchrotron power.

The adopted representative ISM magnetic field strength was 5 $`\mu G`$. This value is typical of the regular field strength in the disks of spiral galaxies (Beck et al., 1996). Although extremal values in spiral arms may reach several times this value, the non-thermal radio emission associated with star formation occurs as cosmic rays travel several kpc from their progenitor supernovae (Condon, 1992). Thus, choice of the regular field strength is more appropriate. For comparison, the same models with a 10 $`\mu G`$ field strength are also considered. Hummel & Beck (1995) quote this value for the regular field strength of NGC 2276 and note that it is likely higher due to interaction with the local ICM.

In Figure 5 we plot the projected radii and line-of-sight velocities for galaxies from the spectroscopically-confirmed star forming sample. To differentiate the anomalous galaxies, we code the different $`q`$ values by symbols. Representations of different values of external pressure are provided for several cases. The top curves assume the parameters given above for an ICM representative of the Coma cluster. The solid curve designates where the external pressure is equivalent to the magnetic pressure provided by a $`5\mu G`$ field, with the dashed curve designating the same for a $`10\mu G`$ field. This means that a galaxy lying along a given curve will experience twice the pressure it would experience in the absence of the cluster effects (assuming it ordinarily had a pressure consistent with the specified magnetic field strength). This would compress the galaxy’s ISM, increasing the magnetic field strength, and thereby increasing the emitted synchrotron power by a factor of two. Similarly, one could interpret the 10 $`\mu G`$ curve as the solution to the case where a galaxy which would have a 5 $`\mu G`$ field in the absence of cluster effects experienced five times its non-cluster pressure. The bottom two curves adopt the same representative magnetic field strengths but assume an ICM more consistent with a poor cluster. The central density of this gas is an order of magnitude lower and the assumed temperature is 2.5 keV.

The first thing to notice from Figure 5 is that the curves are generally very flat. This provides one possible explanation for the lack of an observed correlation of $`q`$ with relative velocity: in general, the thermal pressure (which depends only on position within the cluster) is greater than the ram pressure. In fact, since both $`P_{ram}`$ and $`P_{thermal}`$ depend directly on $`\rho `$, it is easy to solve for the velocity at which they are equal. For the example of the richer cluster given above this velocity is 880 km s<sup>-1</sup>, and for the poorer cluster it is 490 km s<sup>-1</sup>. Only for galaxies with relative velocities above these values will ram pressure dominate over the thermal pressure. There are very few, if any, galaxies on the plot located in regions which require a significant contribution from ram pressure to the total external pressure.

This discussion is obviously highly dependent upon projection effects. To obtain an understanding of the potential magnitude of such effects, we simply apply Equation 5 to the galaxy distribution in the clusters. Integrations of the equation provide the total number of galaxies within a given radius as well as the projected number of galaxies within a given radius. Inside a radius equal to the core radius, half of all galaxies seen in projection actually lie outside the core radius. Similarly, roughly a third of all galaxies seen in projection at $`r2r_c`$ are at higher true separation, about a quarter for $`r3r_c`$, one-fifth for $`r4r_c`$, etc. The ratio flattens out for larger radii, meaning it is relatively constant at around $`15\%`$ for radii equal to or greater than $`10r_c`$.

The actual distribution of the star forming galaxies is somewhat uncertain. The best-fit core radius for all galaxies in the clusters is fairly small, with Adami et al. (1998) citing $`r_c0.2`$ Mpc. However, this fit is likely dominated by the strong concentration of early-type galaxies in the cluster cores. The star forming galaxies have a much larger best-fit core radius ($`r_c0.8`$ Mpc), potentially including a ‘hole’ in their centers (e.g., see Morrison & Owen, 2001; Miller & Owen, 2001b). This argues that the statistical significance of the radial trend in $`q`$ values is likely much stronger than it appears.

Only in the scenario of a richer ICM do many galaxies lie in the region where external pressure is equal to or greater than the representative pressure of the galaxies’ ISM. While there is a tendency for these to have lower $`q`$ values, a number of the low $`q`$ galaxies are at too large radial separations and too small relative velocities to be explained by the above model. Of course, in addition to the effect of projection noted above the true velocity of the galaxies could be higher than their line-of-sight velocities. The effect of these two biases is that galaxies will move to the right in the diagram and possibly up, which does not greatly alter the qualitative conclusions drawn.

The observation that richer clusters have more low $`q`$ galaxies than poorer clusters is actually more consistent with thermal pressure being the source of compression in the galaxies than ram pressure. As can be seen in Figure 5, the thermal pressure provided by a poor cluster ICM is at most comparable to the ambient pressure of the galaxies. This is demonstrated by the bottom curve. The thermal pressure due to the poor cluster ICM is at most only $`40\%`$ of the pressure provided by a 10 $`\mu G`$ field, so only in regions where ram pressure is greater than thermal pressure can such a galaxy be compressed by the factor of two implied by the curve. This is the cause of the slight down turn in this curve on Figure 5.

Ram pressure most likely does affect some star forming galaxies, as seen from HI images of cluster spirals (Cayette et al., 1990) and morphological studies of specific galaxies (e.g., NGC 4522 and UGC6697; Kenney & Koopman, 1998; Nulsen, 1982). The simple treatment above allows a few possibilities by which ram pressure could be the important factor in explaining the low $`q`$ galaxies. Numerical models consistently predict that most of the galaxies infalling into clusters will be on radial orbits (e.g., Ghigna et al., 1998). This suggests that while few galaxies are currently seen in regions where ram pressure should be important, many galaxies will be affected by ram pressure at some point in their lives as they pass through the cluster cores. A cluster galaxy moving at 1000 km s<sup>-1</sup> travels roughly 1 Mpc in about $`10^9`$ years. These representative values for low $`q`$ galaxies indicate that compression due to ram pressure could occur during core passage, followed by a return to equilibrium slower than about a Gyr.

Another possible reason for the lack of correlation of $`q`$ with relative velocity is that the treatment above assumed a static ICM through which galaxies moved. Bulk flows of the ICM are likely, and are frequently invoked to explain the degree of bending of extended radio sources. In fact, Dupke & Bregman (2000) detect significant flows in Perseus through X-ray spectroscopy of iron lines performed with ASCA. Using one pointing centered on the cluster and seven additional pointings surrounding this central one, they find convincing evidence for a net rotation of the ICM. Two opposite pointings of the outer ring registered significantly different velocities for the measured line, consistent with rotation with a circular velocity of $`4100_{3100}^{+2200}`$ km s<sup>-1</sup> (90$`\%`$ confidence interval). Flows of this magnitude in a cluster previously thought to be dynamically relaxed could easily produce the lack of correlation seen for $`q`$ with the projected relative velocity of galaxies.

## 5 Conclusions

In this paper, we have used the cluster galaxy database of Miller & Owen (2001a) to investigate the FIR/radio correlation in nearby Abell clusters. We find that in general the FIR/radio correlation holds remarkably well in nearby clusters of galaxies. For the star forming galaxies in the clusters we find that:

* Star forming galaxies with greatly enhanced radio emission are fairly rare objects. In seventeen studied clusters, only four galaxies had a radio excess of a factor of five or greater, the degree of amplification suggested by other studies. In fact, most of the anomalous objects suggested an enhancement of radio emission relative to FIR emission of only a factor of two to three.
* There is some evidence that cluster star forming galaxies with enhanced radio emission (relative to that expected from the standard FIR/radio correlation) are centrally-located. The evidence for this effect is significant largely because of the sample size used in this study.
* High resolution radio images of a subset of the anomalous galaxies demonstrate that the enhanced radio emission can not be explained by a contribution from a compact AGN. In addition, the population of radio over-luminous spiral galaxies in clusters does not appear to be related to compact starbursts.

In addition to these results, there were similar findings for cluster AGN. Although the radio powers of AGN are greater than those predicted by the canonical FIR/radio correlation value obtained for star forming galaxies, their FIR emission is also well correlated with their radio emission. As for the star forming galaxies, there is evidence that AGN located more toward the centers of the clusters have slightly greater radio emission. This argues that a portion of the net radio emission from such objects is associated with star formation.

These results have a number of implications. First, it is critical when performing studies of this nature to identify and remove AGN from the sample. Since AGN usually have stronger radio emission relative to their FIR emission, their inclusion would overstate the magnitude of the effect of the cluster environment on star formation as evidenced by departures from the FIR/radio correlation. This is likely the cause of the discrepancy between the factor of amplification of radio emission in cluster spirals found in this study (2–3) vs. those suggested by prior studies (5–10). Second, no significant correlation of increased radio emission with galaxy velocity was found. This, in conjunction with the radial distribution of the galaxies, suggests that the ram pressure model alone has difficulty explaining all of the galaxies which lie off of the normal FIR/radio correlation. While a contribution from ram pressure is possible, simple thermal pressure due to the ICM is a more likely candidate for compression of the cluster galaxies’ magnetic fields.

It is hoped that additional studies will shed more light on the nature of the unusual galaxies. We are currently in the process of analyzing their spectra in detail, in hopes of learning more about their star formation histories through relative strengths of Balmer emission and absorption features. It would be interesting if the radio over-luminous star forming galaxies located farther from the cluster cores had stronger post-starburst features, which might underscore the importance of evolutionary changes as galaxies undergo their first core crossing. There is also a need for a larger sample of objects to observe at high resolution in the radio to solidify our conclusions about potential radio enhancement via AGN.

We thank Min Yun for valuable discussions and comments, and an anonymous referee for identifying useful clarifications to the manuscript. N.A.M. would like to thank the NRAO predoctoral program for support of this research.


text:
# Axino as a sterile neutrino and 𝑅 parity violation

## Abstract

We point out that axino can be a natural candidate for sterile neutrino which would accommodate the LSND data with atmospheric and solar neutrino oscillations. It is shown that the so called $`3+1`$ scheme can be easily realized when supersymmetry breaking is mediated by gauge interactions and also $`R`$-parity is properly broken. Among the currently possible solutions to the solar neutrino problem, only the large angle MSW oscillation is allowed in this scheme. The weak scale value of the Higgs $`\mu `$ parameter and the required size of $`R`$-parity violation can be understood by means of spontaneously broken Peccei-Quinn symmetry.

preprint: KAIST-TH 01/01, SNUTP 01-002, hep-ph/0101026

Current data from the atmospheric and solar neutrino experiments are beautifully explained by oscillations among three active neutrino species glb . Another data in favor of neutrino oscillation has been obtained in the LSND experiment lsnd . Reconciliation of these experimental results requires three distinct mass-squared differences, implying the existence of a sterile neutrino $`\nu _s`$ moha . In the four-neutrino oscillation framework, there are two possible scenarios: the $`2+2`$ scheme in which two pairs of close mass eigenstates are separated by the LSND mass gap $`1`$ eV and the $`3+1`$ scheme in which one mass is isolated from the other three by the LSND mass gap. It has been claimed that the LSND results can be compatible with various short-baseline experiments only in the context of the 2+2 scheme 2by2 . However, according to the new LSND results lsnd , the allowed parameter regions are shifted to smaller mixing angle, accepting the $`3+1`$ scheme smirnov ; barger ; giunti . Although it can be realized in a rather limited parameter space, the 3+1 scheme is attractive since the fourth (sterile) neutrino can be added without changing the most favorable picture that the atmospheric and solar neutrino data are explained by the predominant $`\nu _\mu \nu _\tau `$ and $`\nu _e\nu _\mu ,\nu _\tau `$ oscillations, respectively. In particular, the $`3+1`$ scheme with the heaviest $`\nu _s`$ would be an interesting explanation of all existing neutrino data. On the other hand, it is rather difficult to find a well-motivated particle physics model which would yield the desired four-neutrino masses and mixing in a consistent way yanagida .

In this paper, we show that the 3+1 scheme can be easily realized in supersymmetric model with $`U(1)`$ Peccei-Quinn (PQ) symmetry when supersymmetry (SUSY) breaking is mediated by gauge interactions. In this model, axino can be as light as 1 eV, so plays the role of sterile neutrino chun . A proper axino-neutrino mixing can be induced by $`R`$-parity violating couplings which appear as a consequence of spontaneously broken $`U(1)_{PQ}`$. It turns out that only the large angle MSW solution to the solar neutrino problem is allowed in our model. The weak scale value of the Higgs $`\mu `$-parameter and the required size of $`R`$-parity violation can be understood by means of the Frogatt-Nielsen mechanism nielson  of spontaneously broken $`U(1)_{PQ}`$.

The model under consideration includes three sectors: the observable sector, the SUSY-breaking sector, and the PQ sector. The first contains the usual quarks, leptons, and two Higgs superfields, i.e. the MSSM superfields. The second contains a gauge-singlet Goldstino superfield $`X`$ and the gauge-charged messenger superfields $`Y,Y^c`$. Finally the third contains gauge-singlet superfields $`S_k`$ which break $`U(1)_{PQ}`$ by their vacuum expectation values (VEV), as well as gauge-charged superfields $`T,T^c`$ which have the Yukawa coupling with some of $`S_k`$. The superpotential of the model is given by

$$W=hS_3(S_1S_2f_{PQ}^2)+\kappa S_1TT^c+\lambda XYY^c+W_{\mathrm{MSSM}}+W_{\mathrm{SB}}$$
(1)
where $`W_{\mathrm{MSSM}}`$ involves the MSSM fields, and $`W_{\mathrm{SB}}`$ describes SUSY breaking dynamics enforcing $`X`$ develop a SUSY breaking VEV: $`\lambda X=M_X+\theta ^2F_X`$. This VEV generates soft masses of the MSSM fields gmsb , $`m_{\mathrm{soft}}=\frac{\alpha }{2\pi }F_X/M_X`$, which are assumed to be of order the weak scale. One can easily arrange the symmetries of the model, $`U(1)_{PQ}`$ and an additional discrete symmetry to get

$`W_{\mathrm{MSSM}}`$ $`=y_{ij}^EH_1L_iE_j^c+y_{ij}^DH_1Q_iD_j^c+y_{ij}^UH_2Q_iU_j^c+{\displaystyle \frac{y_0}{M_{}}}S_1^2H_1H_2`$ (2)
$`+{\displaystyle \frac{y_i^{}}{M_{}^2}}S_1^3L_iH_2+{\displaystyle \frac{\gamma _{ijk}}{M_{}}}S_1L_iL_jE_k^c+{\displaystyle \frac{\gamma _{ijk}^{}}{M_{}}}S_1L_iQ_jD_k^c,`$
where the Higgs, quark and lepton superfields are in obvious notations. Here $`M_{}`$ denotes the fundamental scale of the model which will be taken to be the grand unification scale, $`M_{GUT}10^{16}`$ GeV. Full details of the model will be presented elsewhere cch .

Let $`A=(\varphi +ia)+\theta \stackrel{~}{a}+\theta ^2F_A`$ denote the axion superfield where $`a`$, $`\varphi `$ and $`\stackrel{~}{a}`$ are the axion, saxion and axino, respectively. It is then convenient to parameterize $`S_1`$ and $`S_2`$ as $`S_1=Se^{A/f_{PQ}}`$ and $`S_2=Se^{A/f_{PQ}}`$. Note that the VEV of $`S`$ is uniquely determined to be $`S=f_{PQ}`$ which would correspond to the axion decay constant if the saxion is stabilized at $`e^{\varphi /f_{PQ}}1`$. After integrating out the SUSY-breaking sector as well as the heavy fields in the PQ sector, the low energy effective action includes the Kähler potential and the superpotential of $`A`$,

$`K_{\mathrm{eff}}`$ $`=`$ $`f_{PQ}^2\{e^{(A+A^{})/f_{PQ}}+e^{(A+A^{})/f_{PQ}}\}+\mathrm{\Delta }K_{\mathrm{eff}}`$
$`W_{\mathrm{eff}}`$ $`=`$ $`\mu _0e^{2A/f_{PQ}}H_1H_2+\mu _i^{}e^{3A/f_{PQ}}L_iH_2`$ (3)
$`+e^{A/f_{PQ}}(\lambda _{ijk}L_iL_jE_k^c+\lambda _{ijk}^{}L_iQ_jD_k^c),`$
where $`\mathrm{\Delta }K_{\mathrm{eff}}`$ is $`A`$-dependent loop corrections involving the SUSY-breaking effects and

$`\mu _0=y_0f_{PQ}^2/M_{},\mu _i^{}=y_i^{}f_{PQ}^3/M_{}^2,`$
$`\lambda _{ijk}=\gamma _{ijk}f_{PQ}/M_{},\lambda _{ijk}^{}=\gamma _{ijk}^{}f_{PQ}/M_{}.`$ (4)
This shows that in our framework the weak-scale value of $`\mu _0`$ and also the smallness of $`R`$-parity violating couplings can be understood by means of the Frogatt-Nielsen mechanism of $`U(1)_{PQ}`$ with $`f_{PQ}M_{}`$ nielson . Although not written explicitly, the coefficients of $`B`$-violating operators $`U_i^cD_j^cD_k^c`$ can be easily arranged to be small enough to avoid too rapid proton decay, including the decays into axino or gravitino pdecay . The best lower bound on $`f_{PQ}`$ is from astrophysical arguments implying $`f_{PQ}10^9`$ GeV pqbound . To accommodate the LSND data, we need the axino-neutrino mixing mass of order 0.1 eV. It turns out that this value is difficult to be obtained for $`f_{PQ}>10^{10}`$ GeV. We thus assume $`f_{PQ}=10^910^{10}`$ GeV with $`M_{}=M_{GUT}`$ for which $`\mu _0`$ takes an weak scale value (with appropriate value of $`y_0`$) and all $`R`$-parity violating couplings are appropriately suppressed.

An important issue here is the saxion stabilization. One dominant contribution to the saxion effective potential comes from $`\mathrm{\Delta }K_{\mathrm{eff}}`$ which is induced mainly by the threshold effects of $`T,T^c`$ having the $`A`$-dependent mass $`M_T=\kappa f_{PQ}e^{A/f_{PQ}}`$. If $`M_T\lambda X`$, one finds

$$\mathrm{\Delta }K_{\mathrm{eff}}\frac{N_T}{16\pi ^2}\frac{M_TM_T^{}}{𝒵_T𝒵_{T^c}}\mathrm{ln}\left(\frac{\mathrm{\Lambda }^2𝒵_T𝒵_{T^c}}{M_TM_T^{}}\right),$$
(5)
where $`N_T`$ is the number of superfields in $`T`$, $`𝒵_T`$ is the Kähler metric of $`T`$, and $`\mathrm{\Lambda }`$ is a cutoff scale which is of order $`M_X`$. If $`M_T\lambda X`$, $`\mathrm{\Delta }K_{\mathrm{eff}}`$ would be enhanced by $`|\lambda X/M_T|^2`$, however then the resulting saxion effective potential can not stabilize $`\varphi `$ at the desired VEV cch . We thus assume $`M_T\lambda X`$. With $`𝒵_T|_{\theta ^2\overline{\theta }^2}m_{\mathrm{soft}}^2`$, the above $`\mathrm{\Delta }K_{\mathrm{eff}}`$ gives a negative-definite saxion potential

$$V_\varphi ^{(1)}\frac{N_T}{16\pi ^2}m_{\mathrm{soft}}^2|\kappa f_{PQ}|^2e^{2\varphi /f_{PQ}}.$$
(6)
There is another (positive-definite) potential from the $`A`$-dependent $`\mu `$-parameter:

$$V_\varphi ^{(2)}e^{4\varphi /f_{PQ}}|\mu _0|^2(|H_1|^2+|H_2|^2).$$
(7)
If $`\kappa f_{PQ}`$ is of order few TeVs, $`\varphi `$ is stabilized by $`V_\varphi ^{(1)}+V_\varphi ^{(2)}`$ at the desired value $`e^{\varphi /f_{PQ}}1`$. This requires a rather small Yukawa coupling $`\kappa 10^6`$ which may be a consequence of some flavor symmetry. For $`e^{\varphi /f_{PQ}}1`$, the saxion and axino masses are estimated to be

$`m_\varphi ^2=(1010^2\mathrm{keV})^2+\mathrm{\Delta }m_\varphi ^2,`$
$`m_{\stackrel{~}{a}}=(10^410^2\mathrm{eV})+\mathrm{\Delta }m_{\stackrel{~}{a}},`$ (8)
where the numbers represent the gauge-mediated contributions for $`f_{PQ}=10^910^{10}`$ GeV, and $`\mathrm{\Delta }m_\varphi `$ and $`\mathrm{\Delta }m_{\stackrel{~}{a}}`$ are the supergravity-mediated contributions which are generically of order $`m_{3/2}`$ lukas . In gauge-mediated SUSY breaking models gmsb , the precise value of $`m_{3/2}`$ depends on the details of SUSY breaking sector. However most of models give $`m_{3/2}1`$ eV, implying $`m_{\stackrel{~}{a}}`$ is dominated by supergravity contribution. In this paper, we assume

$$m_{\stackrel{~}{a}}\mathrm{\Delta }m_{\stackrel{~}{a}}1\mathrm{eV}$$
(9)
which would allow the axino to be a sterile neutrino for the LSND data. We note that although it is generically of order $`m_{3/2}`$, $`\mathrm{\Delta }m_{\stackrel{~}{a}}`$ can be significantly smaller than $`m_{3/2}`$ when the supergravity Kähler potential takes a particular form, e.g. the no-scale form lukas .

Having defined our supersymmetric axion model, we discuss the full $`4\times 4`$ axino-neutrino mass matrix:

$$\frac{1}{2}m_{AB}\nu _A\nu _B$$
(10)
where $`A,B=s,e,\mu ,\tau `$ and $`\nu _s\stackrel{~}{a}`$ with $`m_{ss}=m_{\stackrel{~}{a}}`$. We will work in the field basis in which $`\mu _i^{}L_iH_2`$ $`(i=e,\mu ,\tau )`$ in $`W_{\mathrm{eff}}`$ are rotated away by an appropriate unitary rotation of $`H_1`$ and $`L_i`$. It is straightforward to see that the axino-neutrino mixing mass is given by

$$m_{is}=\frac{ϵ_i\mu _0H_2}{f_{PQ}}0.1\left(\frac{ϵ_i}{10^5}\right)\left(\frac{\mu _0}{600\mathrm{GeV}}\right)\left(\frac{10^{10}\mathrm{GeV}}{f_{PQ}}\right)\mathrm{eV},$$
(11)
where $`ϵ_i=\mu _i^{}/\mu _0`$ for $`\mu _i^{}`$ in $`W_{\mathrm{eff}}`$ before $`\mu _i^{}L_iH_2`$ are rotated away. Note that this axino-neutrino mixing survives under the unitary rotation eliminating $`\mu _i^{}L_iH_2`$ since $`L_i`$ and $`H_1`$ have different $`U(1)_{PQ}`$-charges. This charge difference is also responsible for the suppression of $`R`$-parity violating couplings as well as the weak-scale value of $`\mu _0`$. The $`3\times 3`$ mass matrix of active neutrinos is induced by $`R`$-parity violating parameters hall . At tree-level,

$$m_{ij}\frac{g_a^2\stackrel{~}{\nu }_i^{}\stackrel{~}{\nu }_j^{}}{M_a},$$
(12)
where $`M_a`$ denote the gaugino masses. The sneutrino VEV’s $`\stackrel{~}{\nu }_i`$ are determined by the bilinear $`R`$-parity violations in the SUSY-breaking scalar potential: $`m_{L_iH_1}^2L_iH_1^{}+B_i^{}L_iH_2`$. In our model, nonzero values of $`m_{L_iH_1}^2`$ and $`B_i^{}`$ at the weak scale arise through renormalization group evolution (RGE), mainly by the coupling $`\lambda _{i33}^{}y_b`$ where $`y_b`$ is the $`b`$-quark Yukawa coupling hwang . Moreover, $`BH_1H_2`$ arises also through RGE which predicts a large $`\mathrm{tan}\beta 4060`$ sarid . We then find

$$m_{ij}10^2t^4\left(\frac{\lambda _{i33}^{}y_b}{10^6}\right)\left(\frac{\lambda _{j33}^{}y_b}{10^6}\right)\mathrm{eV}$$
(13)
where $`t=\mathrm{ln}(M_X/m_{\stackrel{~}{l}})/\mathrm{ln}(10^3)`$ for the slepton mass $`m_{\stackrel{~}{l}}`$. Here we have taken $`m_{\stackrel{~}{l}}300`$ GeV and $`\mu _02m_{\stackrel{~}{l}}`$ which has been suggested to be the best parameter range for correct electroweak symmetry breaking sarid . There are also bunch of loop corrections to $`m_{ij}`$ from $`R`$-parity violating couplings kang , however in our case they are too small to be relevant.

Let us now see how nicely all the neutrino masses and mixing parameters are fitted in our framework. The analysis of Ref. giunti  leads to the four parameter regions, R1–R4 of Table 1, accommodating the LSND with short baseline results. In our model, Eqs. (9) and (11) can easily produce these LSND values of the largest mass eigenvalue $`m_4m_{ss}=m_{\stackrel{~}{a}}1`$ eV and the mixing matrix elements $`U_{i4}m_{is}/m_{ss}0.1`$. The effective $`3\times 3`$ mass matrix of active neutrinos is then given by

$$m_{ij}^{\mathrm{eff}}=m_{ij}m_{is}m_{js}/m_{ss}.$$
(14)
Upon ignoring the small loop corrections, this mass matrix has rank two, and can be written as

$$m_{ij}^{\mathrm{eff}}=m_\mathrm{x}\widehat{\mathrm{x}}_i\widehat{\mathrm{x}}_j+m_\mathrm{y}\widehat{\mathrm{y}}_i\widehat{\mathrm{y}}_j$$
(15)
where $`\widehat{\mathrm{x}}_i`$ and $`\widehat{\mathrm{y}}_j`$ are the unit vectors in the direction of $`m_{is}`$ and $`\stackrel{~}{\nu }_j`$, respectively. Remarkably, the mass scale $`m_\mathrm{x}(m_{is}/m_{ss})^2m_{ss}10^2`$ eV gives the right range of the atmospheric neutrino mass. Eq. (13) shows that $`m_\mathrm{y}`$ is also in the range of $`10^2`$ eV, so $`m^{\mathrm{eff}}`$ would be able to provide the right solar neutrino mass unless $`\mathrm{\Delta }m_{sol}^210^4`$ eV2. Note from Eq. (Axino as a sterile neutrino and $`R`$ parity violation) that the typical size of $`ϵ_i,\lambda _{ijk},\lambda _{ijk}^{}`$ is around $`10^6`$ for $`f_{PQ}10^{10}`$ GeV and $`M_{}10^{16}`$ GeV. To discuss the consequeces of the mass matrix (15) on the atmospheric and solar neutrino oscillations of three active neutrinos, we diagonalize it to find the mass eigenvalues

$$m_{2,3}=\frac{1}{2}\left(m_\mathrm{x}+m_\mathrm{y}\pm \sqrt{(m_\mathrm{x}+m_\mathrm{y}\mathrm{cos}^22\xi )^2+m_\mathrm{y}^2\mathrm{sin}^22\xi }\right)$$
(16)
and the $`3\times 3`$ mixing matrix

$$U=(\widehat{\mathrm{z}}^T,\widehat{\mathrm{w}}^Tc_\theta \widehat{\mathrm{x}}^Ts_\theta ,\widehat{\mathrm{w}}^Ts_\theta +\widehat{\mathrm{x}}^Tc_\theta ).$$
(17)
where $`\widehat{\mathrm{z}}\widehat{\mathrm{x}}\times \widehat{\mathrm{y}}/|\widehat{\mathrm{x}}\times \widehat{\mathrm{y}}|`$, $`\widehat{\mathrm{w}}\widehat{\mathrm{x}}\times \widehat{\mathrm{z}}/|\widehat{\mathrm{x}}\times \widehat{\mathrm{z}}|`$, and $`\mathrm{tan}2\theta =m_\mathrm{y}\mathrm{sin}2\xi /(m_\mathrm{x}+m_\mathrm{y}\mathrm{cos}2\xi )`$ for $`\mathrm{cos}\xi =\widehat{\mathrm{x}}\widehat{\mathrm{y}}`$. The Super-Kamiokande data sk-atm  combined with the CHOOZ result chooz  imply that $`U_{\mu 3}^2U_{\tau 3}^21/2`$ and $`U_{e3}^21`$. The solutions to the solar neutrino problem can have either a large mixing angle (LA): $`U_{e1}^2U_{e2}^21/2`$, or a small mixing angle (SA): $`U_{e1}^21`$. This specify the first column $`\widehat{\mathrm{z}}^T`$ of $`U`$ as $`\widehat{\mathrm{z}}(1/\sqrt{2},1/2,1/2)`$ for LA and $`\widehat{\mathrm{z}}(1,0,0)`$ for SA up to sign ambiguities. Since $`\widehat{\mathrm{x}}\widehat{\mathrm{z}}=0`$, the pattern $`\widehat{\mathrm{z}}(1,0,0)`$ implies $`\widehat{\mathrm{x}}_e0`$. This leads to a too small $`U_{e4}m_{es}/m_{ss}10^2`$, so the SA solution is not allowed within our model. Among various LA solutions to the solar neutrino problem, only the large-angle MSW solution with $`\mathrm{\Delta }m_{sol}^210^4\mathrm{eV}^2`$ can be naturally fitted since $`m_\mathrm{x}m_\mathrm{y}10^2`$ eV in our scheme. It is remarkable that $`f_{PQ}10^{10}`$ GeV and $`M_{}M_{GUT}`$ lead to the right size of $`R`$-parity violation yielding the desired values of $`m_{is}`$ and $`m_{ij}`$ also for the atmospheric and solar neutrino masses.

To see the feasibility of our whole scheme, we scanned our parameter space to reproduce the allowed LSND islands R1–R4 of Table 1 together with the following range of atmospheric and solar neutrino parameters: $`\mathrm{\Delta }m_{31}^2=(18)\times 10^3\mathrm{eV}^2,\mathrm{\Delta }m_{21}^2=(0.18)\times 10^4\mathrm{eV}^2`$ and $`\mathrm{tan}^2\theta _{23}=0.333.8`$, $`\mathrm{tan}^2\theta _{12}=0.23.0`$, $`\mathrm{tan}^2\theta _{13}0.055`$ glb . Our parameter space consists of $`m_{ss},m_{is},\lambda _{i33}^{}y_b`$ whose values are centered around 1 eV, $`0.1`$ eV, $`10^6`$, respectively. For R1 and R4, we could find some limited parameter space, however they need a strong alignment between $`\widehat{\mathrm{x}}`$ and $`\widehat{\mathrm{y}}`$ and also a large cancellation between $`m_\mathrm{x}`$ and $`m_\mathrm{y}`$. On the other hand, R2 and R3 do not require any severe fine tuning of parameters, although the bi-maximal mixing is obtained by some accident. We provide in Fig. 1 the scatter plots on the planes of $`(m_\mathrm{y},m_\mathrm{x})`$ and $`(\mathrm{cos}\xi ,m_\mathrm{x})`$ for the LSND island R2. The scatter plots for R3 have similar shape. The full details of the analysis will be presented elsewhere cch .

To conclude, we have shown that the 3+1 scheme of four-neutrino oscillation can be nicely obtained in supersymmetric model endowed with the PQ solution to the strong CP problem and supersymmetry breaking mediated by gauge interactions.

Acknowledgement: This work is supported by BK21 project of the Ministry of Education, KOSEF through the CHEP of KNU and KRF Grant No. 2000-015-DP0080.