text:
# 1 Introduction

## 1 Introduction

Relativistic jets have been observed in active galactic nuclei (AGNs) and microquasars in our Galaxy, and it is believed that they originate from the regions near accreting black holes or neutron stars \[Meier, Koide, & Uchida\]. To investigate the dynamics of accretion disks and the associated jet formation, we use a full 3-D GRMHD code. One of the most promising models for jet formation is the magnetic-acceleration model \[Blandford & Payne\]. The magnetic-acceleration mechanism has been proposed not only for AGN jets, but also for protostellar jets (see \[Meier, Koide, & Uchida\]).

Recently, Koide, Shibata, & Kudoh (1999) \[Koide, Shibata, & Kudoh\] have investigated the dynamics of an accretion disk initially threaded by a uniform poloidal magnetic field in a non-rotating corona (either in a state of steady fall or in hydrostatic equilibrium) around a non-rotating black hole. The numerical results show that the disk loses angular momentum by magnetic braking, then falls into the black hole. The infalling motion of the disk, which is faster than in the non-relativistic case because of the general-relativistic effect below $`3r_\mathrm{S}`$ ($`r_\mathrm{S}`$ is the Schwarzschild radius), is strongly decelerated at the shock formed by the centrifugal force around $`r=2r_\mathrm{S}`$ by the rotation of the disk. Plasmas near the shock are accelerated by the $`𝐉\times 𝐁`$ force, which forms bipolar relativistic jets. Inside this magnetically driven jet, the gradient of gas pressure also generates a jet over the shock region (gas-pressure-driven jet). This two-layered jet structure is formed both in a hydrostatic corona and in a steady-state falling corona. Koide et al. (2000) \[Koide et al 2000\] have also developed a new GRMHD code in Kerr geometry and have found that, with a rapidly rotating ($`a=0.95`$) black-hole magnetosphere, the maximum velocity of the jet is 0.9 c and its terminal velocity 0.85 c. All of the previous 2-D GRMHD simulations described here were made assuming axisymmetry with respect to the $`z`$-axis and mirror symmetry with respect to the plane $`z=0`$; the axisymmetric assumption suppressed the azimuthal instabilities.

## 2 3-D GRMHD Simulations: Equations and Numerical Techniques

Our basic equations are those of Maxwell for the fields and a set of general-relativistic equations representing the plasma, namely the equations of conservation of mass, momentum, and energy for a single-component conducting fluid \[Weinberg 1972, Thorne, Price, and Macdonald 1986, Koide, Shibata, & Kudoh\].

Our simulations employ the equations of general-relativistic MHD together with the Simplified Total-Variation-Diminishing (STVD) method. The STVD method was developed by Davis (1984) \[Davis 1984\] for studying violent phenomena such as shocks \[Nishikawa et al. 1997, Koide, Shibata, & Kudoh\]. This method is as Lax-Wendoroff’s method with additional diffusion term. We have checked that this method respects the energy-conservation law and its propagation properties \[Koide, Shibata, & Kudoh\].

## 3 Initial 3-D GRMHD simulations with a Schwarzschild black hole

In order to investigate how accretion disks near black holes evolve under the influence of accretion instabilities such as the magnetorotational instability \[Balbus & Hawley 1998\] and accretion-ejection instability (AEI) \[Tagger & Pellat, 1999\], the use of a fully 3-D GRMHD is essential.

### 3.1 Initial and boundary conditions

In the assumed initial state, the simulation region is divided into two parts: a background corona around a black hole, and an accretion disk (Fig. 1a). The coronal plasma is set in a state of transonic free-fall flow, as in the case of the transonic flows with $`\mathrm{\Gamma }=5/3`$ and $`H=1.3`$; here the sonic point is located at $`r=1.6r_\mathrm{S}`$. The Keplerian disk in the corona is set in the following way. The disk region is located at $`r>r_\mathrm{D}3r_\mathrm{S},|\mathrm{cos}\theta |<\delta =1/8`$. Here the density is 100 times that of the background corona (Fig. 1a), while the orbital velocity is relativistic and purely azimuthal: $`v_\varphi =v_\mathrm{K}c/[2(r/r_\mathrm{S}1)]^{1/2}`$. (Note that this equation reduces to the Newtonian Keplerian velocity $`v_\varphi =\sqrt[]{GM/r}`$ in the non-relativistic limit $`r_\mathrm{S}/r1`$). The pressure of both the corona and the disk are assumed equal to that of the transonic solution. The initial conditions for the entire plasma around the black hole are: $`\rho =\rho _{\mathrm{ffc}}+\rho _{\mathrm{dis}}`$

$$\rho _{\mathrm{dis}}=\{\begin{array}{cc}100\rho _{\mathrm{ffc}}& \mathtt{}(r>r_\mathrm{D}\mathtt{}\mathrm{and}\mathtt{}|\mathrm{cot}\theta |<\delta )\\ 0& \mathtt{}(rr_\mathrm{D}\mathtt{}\mathrm{or}\mathtt{}|\mathrm{cot}\theta |\delta )\end{array}$$
(1)
$$(v_r,v_\theta ,v_\varphi )=\{\begin{array}{cc}(0,0,v_\mathrm{K})& \mathtt{}(r>r_\mathrm{D}\mathtt{}\mathrm{and}\mathtt{}|\mathrm{cot}\theta |<\delta )\\ (v_{\mathrm{ffc}},0,0)& \mathtt{}(rr_\mathrm{D}\mathtt{}\mathrm{or}\mathtt{}|\mathrm{cot}\theta |\delta )\end{array}$$
(2)
where we set $`\delta =0.125`$; the smoothing length is $`0.3r_\mathrm{S}`$.

In addition, there is a magnetic field crossing the accretion disk perpendicularly. We set it to the Wald solution \[Wald 1974\], which represents the uniform magnetic field around a Kerr black hole: $`B_r=B_0\mathrm{cos}\theta `$, $`B_\theta =\alpha B_0\mathrm{sin}\theta `$ (where $`\alpha `$ is the lapse function, $`\alpha =(1r_\mathrm{S}/r)^{1/2}`$). At the inner edge of the accretion disk, the proper Alfvén velocity is $`v_\mathrm{A}=0.015c`$ in a typical case with $`B_0=0.3\sqrt[]{\rho _0c^2}`$, where the Alfvén velocity in the fiducial observer, $`v_\mathrm{A}B/\sqrt[]{\rho +[\mathrm{\Gamma }p/(\mathrm{\Gamma }1)+B^2]/c^2}`$. The plasma beta of the corona at $`r=3r_\mathrm{S}`$ is $`\beta p/B^2=1.40`$. The simulation is performed in the region $`1.1r_\mathrm{S}r20r_\mathrm{S}`$, $`0\theta \pi `$, $`0\varphi 2\pi `$ with $`100\times 120\times 60`$ meshes. The effective linear mesh widths at $`r=1.1r_\mathrm{S}`$ and at $`r=20r_\mathrm{S}`$ are $`5.38\times 10^3r_\mathrm{S}`$ and $`0.97r_\mathrm{S}`$, respectively, while the angular spacings along the polar and azimuthal directions are $`5.2\times 10^2`$ rad. A radiative boundary condition is imposed at $`r=1.1r_\mathrm{S}`$ and at $`r=20r_\mathrm{S}`$. The computations were made on an ORIGIN 2000 computer with 0.898 GB internal memory, and they used about 47 hours of CPU time for 10000 time steps with $`100\times 120\times 60`$ meshes.

### 3.2 Simulation results

Figure 1 shows the evolution of 3-D simulation performed in the region $`1.1r_\mathrm{S}r20r_\mathrm{S}`$, $`0\theta \pi `$, and $`0\varphi 2\pi `$ with $`100\times 120\times 60`$ meshes. The parameters used in this simulation are the same as those of the axisymmetric simulations shown in Fig. 6 of \[Koide, Shibata, & Kudoh\]. In this figure, the colored shading shows the proper mass density on logarithmic scales; the vector plots show the magnetic field.

The black circle represents the black hole. Figures 1a presents the initial conditions, which are the same as in the 2-D simulations \[Koide, Shibata, & Kudoh\]. At $`t=39.2\tau _\mathrm{S}`$ (Fig. 1b) comparing with Fig. 6c ($`t=40.0\tau _\mathrm{S}`$) in \[Koide, Shibata, & Kudoh\] the jet is less generated. At $`t=60.0\tau _\mathrm{S}`$ as shown in Fig. 1c, the jet is generated as 2-D simulation at the earlier time ($`t=40.0\tau _\mathrm{S}`$). At $`t=73.9\tau _\mathrm{S}`$, the jet is clearly created around $`r=4.5r_\mathrm{S}`$, which is shown by the enhanced density (Fig. 1d). As in the 2-D simulations the jet is generated in a hollowed cylindrical form.

The delay of jet formation seems to be due to reduction of shock formation at $`r=2r_\mathrm{S}`$ caused by the additional freedom in the azimuthal direction. Further investigation will be reported elsewhere.

## 4 Discussion

Recently, review articles on magnetohydrodynamic production of relativistic jets have reported unsolved questions related to jet formation and its mechanisms which determine the velocity of jets and time variations of jet flux \[Meier, Koide, & Uchida, Blandford 2001, Meier 2002\].

This simulation result is initial and we will perform more simulations and investigate effects of the third dimension. This simulation study will be extended to understand the different (high/soft and low/hard) states \[Fender(2001)\] and different combination of accreting rate $`\dot{m}\dot{M}/\dot{M}_{Edd}`$ and the angular momentum $`jJ/J_{\mathrm{max}}`$ \[Meier 2002\]. Further results will be reported elsewhere.

## Acknowledgments

K.N. is partially supported by NSF ATM 9730230, ATM-9870072, ATM-0100997, and INT-9981508. The simulations have been performed on ORIGIN 2000 at NCSA which is supported by NSF.


format: text

text:
# 1 Introduction

## 1 Introduction

From olfaction to vision and audition, a growing number of experiments  are examining the responses of sensory neurons to natural stimuli. Observing the full dynamic range of neural responses may require using stimulus ensembles which approximate those  occurring in nature, and it is an attractive hypothesis that the neural representation of these natural signals may be optimized in some way . Many neurons exhibit strongly nonlinear and adaptive responses that are unlikely to be predicted from a combination of responses to simple stimuli; in particular neurons have been shown to adapt to the distribution of sensory inputs, so that any characterization of these responses will depend on context . Finally, the variability of neural response decreases substantially when complex dynamical, rather than static, stimuli are used . All of these arguments point to the need for general tools to analyze the neural responses to complex, naturalistic inputs.

The stimuli analyzed by sensory neurons are intrinsically high dimensional, with dimensions $`D10^210^3`$. For example, in the case of visual neurons, the input is specified as light intensity on a grid of at least $`10\times 10`$ pixels. Each of the presented stimuli can be described as a vector $`𝐬`$ in this high dimensional stimulus space. It is important that stimuli need not be pictured as being drawn as isolated points from this space. Thus, if stimuli are varying continuously in time we can think of the stimulus $`𝐬`$ as describing a recent window of the stimulus history (e. g., the past $`K`$ frames of the movie, with dimensionality $`K`$ times larger than for the description of a single frame) and then the distribution of stimuli $`P(𝐬)`$ is sampled along some meandering trajectory in this space; we will assume this process is ergodic, so that we can exchange averages over time with averages over the true distribution as needed.

Even though direct exploration of a $`D10^210^3`$ dimensional stimulus space is beyond the constraints of experimental data collection, progress can be made provided we make certain assumptions about how the response has been generated. In the simplest model, the probability of response can be described by one receptive field (RF) or linear filter . The receptive field can be thought of as a template or special direction $`𝐯`$ in the stimulus space such that the neuron’s response depends only on a projection of a given stimulus $`𝐬`$ onto $`𝐯`$, although the dependence of the response on this projection can be strongly nonlinear. In this simple model, the reverse correlation method  can be used to recover the vector $`𝐯`$ by analyzing the neuron’s responses to Gaussian white noise. In a more general case, the probability of the response depends on projections $`s_i=\widehat{e}_i𝐬`$ of the stimulus $`𝐬`$ on a set of vectors $`\{\widehat{e}_1,\widehat{e}_2,\mathrm{},\widehat{e}_n\}`$:

$$P(\mathrm{spike}|𝐬)=P(\mathrm{spike})f(s_1,s_2,\mathrm{},s_n),$$
(1)
where $`P(\mathrm{spike}|𝐬)`$ is the probability of a spike given a stimulus $`𝐬`$ and $`P(\mathrm{spike})`$ is the average firing rate. Even though the ideas developed below can be used to analyze input–output functions $`f`$ with respect to different neural responses, such as patterns of spikes in time , we choose a single spike as the response of interest. The vectors $`\{\widehat{e}_i\}`$ may also describe how the time dependence of stimulus $`𝐬`$ affects the probability of a spike. We will call the subspace spanned by the set of vectors $`\{\widehat{e}_i\}`$ the relevant subspace (RS).

Equation (1) in itself is not yet a simplification if the dimensionality $`n`$ of the RS is equal to the dimensionality $`D`$ of the stimulus space. In this paper we will use the idea of dimensionality reduction  and assume that $`nD`$. The input–output function $`f`$ in Eq. (1) can be strongly nonlinear, but it is presumed to depend only on a small number of projections. This assumption appears to be less stringent than that of approximate linearity which one makes when characterizing neuron’s response in terms of Wiener kernels (see, for example, the discussion in Section 2.1.3 of Ref. ). The most difficult part in reconstructing the input–output function is to find the RS. Note that for $`n>1`$, a description in terms of any linear combination of vectors $`\{\widehat{e}_i\}`$ is just as valid, since we did not make any assumptions as to a particular form of nonlinear function $`f`$.

Once the relevant subspace is known, the probability $`P(\mathrm{spike}|𝐬)`$ becomes a function of only few parameters, and it becomes feasible to map this function experimentally, inverting the probability distributions according to Bayes’ rule:

$$f(\{s_i\})=\frac{P(\{s_i\}|\mathrm{spike})}{P(\{s_i\})}.$$
(2)
If stimuli are chosen from a correlated Gaussian noise ensemble, then the neural response can be characterized by the spike–triggered covariance method . It can be shown that the dimensionality of the RS is equal to the number of nonzero eigenvalues of a matrix given by a difference between covariance matrices of all presented stimuli and stimuli conditional on a spike. Moreover, the RS is spanned by the eigenvectors associated with the nonzero eigenvalues multiplied by the inverse of the a priori covariance matrix. Compared to the reverse correlation method, we are no longer limited to finding only one of the relevant directions $`\widehat{e}_i`$. Both the reverse correlation and the spike–triggered covariance method, however, give rigorously interpretable results only for Gaussian distributions of inputs.

In this paper we investigate whether it is possible to lift the requirement for stimuli to be Gaussian. When using natural stimuli, which certainly are non–Gaussian, the RS cannot be found by the spike–triggered covariance method. Similarly, the reverse correlation method does not give the correct RF, even in the simplest case where the input–output function in Eq. (1) depends only on one projection. However, vectors that span the RS clearly are special directions in the stimulus space independent of assumptions about $`P(𝐬)`$. This notion can be quantified by Shannon information, and an optimization problem can be formulated to find the RS. We illustrate how the optimization scheme works with natural stimuli for model orientation sensitive cells with one and two relevant directions, much like simple and complex cells found in primary visual cortex. It also is possible to estimate average errors in the reconstruction. The advantage of this optimization scheme is that it does not rely on any specific statistical properties of the stimulus ensemble, and can thus be used with natural stimuli.

## 2 Information as an objective function

When analyzing neural responses, we compare the a priori probability distribution of all presented stimuli with the probability distribution of stimuli which lead to a spike . For Gaussian signals, the probability distribution can be characterized by its second moment, the covariance matrix. However, an ensemble of natural stimuli is not Gaussian, so that neither second nor any other finite number of moments is sufficient to describe the probability distribution. In this situation, Shannon information provides the rigorous way of comparing two probability distributions. The average information carried by the arrival time of one spike is given by

$$I_{\mathrm{spike}}=d^D𝐬P(𝐬|\mathrm{spike})\mathrm{log}_2\left[\frac{P(𝐬|\mathrm{spike})}{P(𝐬)}\right].$$
(3)
The information per spike as written in (3) is difficult to estimate experimentally, since it requires either sampling of the high–dimensional probability distribution $`P(𝐬|\mathrm{spike})`$ or a model of how spikes were generated, i.e. the knowledge of low–dimensional RS. However it is possible to calculate $`I_{\mathrm{spike}}`$ in a model–independent way, if stimuli are presented multiple times to estimate the probability distribution $`P(\mathrm{spike}|𝐬)`$. Then,

$$I_{\mathrm{spike}}=\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}\mathrm{log}_2\left[\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}\right]_𝐬,$$
(4)
where the average is taken over all presented stimuli. As discussed in , this is useful in practice because we can replace the ensemble average $`_𝐬`$ with a time average, and $`P(\mathrm{spike}|𝐬)`$ with the time dependent spike rate $`r(t)`$. Note that for a finite dataset of $`N`$ trials, the obtained value $`I_{\mathrm{spike}}(N)`$ will be on average larger than $`I_{\mathrm{spike}}(\mathrm{})`$, with difference $`N_{\mathrm{stimuli}}/(N_{\mathrm{spike}}\mathrm{\hspace{0.17em}2}\mathrm{ln}2)`$, where $`N_{\mathrm{stimuli}}`$ is the number of different stimuli, and $`N_{\mathrm{spike}}`$ is the number of elicited spikes . The true value $`I_{\mathrm{spike}}`$ can also be found by extrapolating to $`N\mathrm{}`$ . Measurement of $`I_{\mathrm{spike}}`$ in this way provides a model independent benchmark against which we can compare any description of the neuron’s input–output relation.

Having in mind a model in which spikes are generated according to projection onto a low dimensional subspace, we start by projecting all of the presented stimuli on a particular direction $`𝐯`$ in the stimulus space, and form probability distributions

$`P_𝐯(x|\mathrm{spike})`$ $`=`$ $`\delta (x𝐬𝐯)|\mathrm{spike}_𝐬,`$ (5)
$`P_𝐯(x)`$ $`=`$ $`\delta (x𝐬𝐯)_𝐬,`$ (6)
where $`\mathrm{}|\mathrm{spike}`$ denotes an expectation value conditional on the occurrence of a spike. The information

$$I(𝐯)=𝑑xP_𝐯(x|\mathrm{spike})\mathrm{log}_2\left[\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}\right]$$
(7)
provides an invariant measure of how much the occurrence of a spike is determined by projection on the direction $`𝐯`$. It is a function only of direction in the stimulus space and does not change when vector $`𝐯`$ is multiplied by a constant. This can be seen by noting that for any probability distribution and any constant $`c`$, $`P_{c𝐯}(x)=c^1P_𝐯(x/c)`$. When evaluated along any vector, $`I(𝐯)I_{\mathrm{spike}}`$. The total information $`I_{\mathrm{spike}}`$ can be recovered along one particular direction only if $`𝐯=\widehat{e}_1`$, and the RS is one dimensional.

By analogy with (7), one could also calculate information $`I(𝐯_1,\mathrm{},𝐯_n)`$ along a set of several directions $`\{𝐯_1,\mathrm{},𝐯_n\}`$ based on the multi-point probability distributions:

$`P_{𝐯_1,\mathrm{},𝐯_n}(\{x_i\}|\mathrm{spike})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\delta (x_i𝐬𝐯_i)|\mathrm{spike}_𝐬,`$ (8)
$`P_{𝐯_1,\mathrm{},𝐯_n}(\{x_i\})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\delta (x_i𝐬𝐯_i)_𝐬.`$ (9)
If we are successful in finding all of the $`n`$ directions $`\widehat{e}_i`$ in the input–output relation (1), then the information evaluated in this subspace will be equal to the total information $`I_{\mathrm{spike}}`$. When we calculate information along a set of $`n`$ vectors that are slightly off from the RS, the answer is, of course, smaller than $`I_{\mathrm{spike}}`$ and is initially quadratic in deviations $`\delta 𝐯_i`$. One can therefore hope to find the RS by maximizing information with respect to $`n`$ vectors simultaneously. The information does not increase if more vectors outside the RS are included. For uncorrelated stimuli, any vector or a set of vectors that maximizes $`I(𝐯)`$ belongs to the RS. On the other hand, the result of optimization with respect to a number of vectors $`k<n`$ may deviate from the RS if stimuli are correlated. To find the RS, we first maximize $`I(𝐯)`$, and compare this maximum with $`I_{\mathrm{spike}}`$, which is estimated according to (4). If the difference exceeds that expected from finite sampling corrections, we increment the number of directions with respect to which information is simultaneously maximized.

The information $`I(𝐯)`$ as defined by (7) is a continuous function, whose gradient can be computed. We find

$$_𝐯I=𝑑xP_𝐯(x)\left[𝐬|x,\mathrm{spike}𝐬|x\right]\left[\frac{d}{dx}\frac{P_𝐯(x|\mathrm{spike})}{P_𝐯(x)}\right],$$
(10)
where

$$𝐬|x,\mathrm{spike}=\frac{1}{P(x|\mathrm{spike})}d^D𝐬𝐬\delta (x𝐬𝐯)P(𝐬|\mathrm{spike}),$$
(11)
and similarly for $`𝐬|x`$. Since information does not change with the length of the vector, $`𝐯_𝐯I=0`$ (which can also be seen from (10) directly).

As an optimization algorithm, we have used a combination of gradient ascent and simulated annealing algorithms: successive line maximizations were done along the direction of the gradient. During line maximizations, a point with a smaller value of information was accepted according to Boltzmann statistics, with probability $`\mathrm{exp}[(I(𝐯_{i+1})I(𝐯_i))/T]`$. The effective temperature $`T`$ is reduced upon completion of each line maximization.

## 3 Results

We tested the scheme of looking for the most informative directions on model neurons that respond to stimuli derived from natural scenes. As stimuli we used patches of black and white photos digitized to 8 bits, in which no corrections were made for camera’s light intensity transformation function. Our goal is to demonstrate that even though the correlations present in natural scenes are non–Gaussian, they can be successfully removed from the estimate of vectors defining the RS.

### 3.1 A model simple cell

Our first example is based on the properties of simple cells found in the primary visual cortex. A model phase and orientation sensitive cell has a single relevant direction $`\widehat{e}_1`$ shown in Fig. 1(a). A given frame $`𝐬`$ leads to a spike if the projection $`s_1=𝐬\widehat{e}_1`$ reaches a threshold value $`s_t`$ in the presence of noise:

$$\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}f(s_1)=\theta (s_1s_t+\xi ),$$
(12)
where a Gaussian random variable $`\xi `$ of variance $`\sigma ^2`$ models additive noise, and the function $`\theta (x)=1`$ for $`x>0`$, and zero otherwise. Together with the RF $`\widehat{e}_1`$, the parameters $`s_t`$ for threshold and the noise variance $`\sigma ^2`$ determine the input–output function.

The spike–triggered average (STA), or reverse correlation function , shown in Fig. 1(b), is broadened because of spatial correlations present in the stimuli. In a model, the effect of noise on our estimate of the STA can be eliminated by averaging the presented stimuli weighted with the exact firing rate, as opposed to using a histogram of responses to estimate $`P(\mathrm{spike}|𝐬)`$ from a finite set of trials. We have used this “exact” STA,

$$𝐯_{\mathrm{sta}}=d^D𝐬𝐬P(𝐬|\mathrm{spike})=\frac{1}{P(\mathrm{spike})}d^D𝐬P(𝐬)𝐬P(\mathrm{spike}|𝐬),$$
(13)
in calculations presented in Fig. 1(bc). If stimuli were drawn from a Gaussian probability distribution, they could be decorrelated by multiplying $`𝐯_{\mathrm{sta}}`$ by the inverse of the a priori covariance matrix, according to the reverse correlation method, $`\widehat{v}_{Gaussianest}C_{apriori}^1𝐯_{\mathrm{sta}}`$. The procedure is not valid for non–Gaussian stimuli and nonlinear input–output functions (1). The result of such a decorrelation is shown in Fig. 1(c). It clearly is missing some of the structure in the model filter, with projection $`\widehat{e}_1\widehat{v}_{Gaussianest}0.14`$. The discrepancy is not due to neural noise or finite sampling, since the “exact” STA was decorrelated; the absence of noise in the exact STA also means that there would be no justification for smoothing the results of the decorrelation. The discrepancy between the true receptive field and the decorrelated STA increases with the strength of nonlinearity in the input–output function.

In contrast, it is possible to obtain a good estimate of the relevant direction $`\widehat{e}_1`$ by maximizing information directly, see panel (d). A typical progress of the simulated annealing algorithm with decreasing temperature $`T`$ is shown in Fig. 1(e). There we plot both the information along the vector, and its projection on $`\widehat{e}_1`$. The final value of projection depends on the size of the data set, see below. In the example shown in Fig. 1 there were $`50,000`$ spikes with average probability of spike $`0.05`$ per frame, and the reconstructed vector has projection $`\widehat{v}_{max}\widehat{e}_10.9`$. Having estimated the RF, one can proceed to sample the nonlinear input-output function. This is done by constructing histograms for $`P(𝐬\widehat{v}_{\mathrm{max}})`$ and $`P(𝐬\widehat{v}_{\mathrm{max}}|\mathrm{spike})`$ of projections onto vector $`\widehat{v}_{\mathrm{max}}`$ found by maximizing information, and taking their ratio, as in Eq. (2). In Fig. 1(f) we compare $`P(\mathrm{spike}|𝐬\widehat{v}_{\mathrm{max}})`$ (crosses) with the probability $`P(\mathrm{spike}|s_1)`$ used in the model (solid line).

### 3.2 Estimated deviation from the optimal direction

When information is calculated from a finite data set, the vector $`𝐯`$ which maximizes $`I`$ will deviate from the true RF $`\widehat{e}_1`$. The deviation $`\delta 𝐯=𝐯\widehat{e}_1`$ arises because the probability distributions are estimated from experimental histograms and differ from the distributions found in the limit on infinite data size. For a simple cell, the quality of reconstruction can be characterized by the projection $`𝐯\widehat{e}_1=1\frac{1}{2}\delta 𝐯^2`$, where both $`𝐯`$ and $`\widehat{e}_1`$ are normalized, and $`\delta 𝐯`$ is by definition orthogonal to $`\widehat{e}_1`$. The deviation $`\delta 𝐯A^1I`$, where $`A`$ is the Hessian of information. Its structure is similar to that of a covariance matrix:

$$A_{ij}=\frac{1}{\mathrm{ln}2}𝑑xP(x|\mathrm{spike})\left(\frac{d}{dx}\mathrm{ln}\frac{P(x|\mathrm{spike})}{P(x)}\right)^2(s_is_j|xs_i|xs_j|x).$$
(14)
When averaged over possible outcomes of $`N`$ trials, the gradient of information is zero for the optimal direction. Here in order to evaluate $`\delta 𝐯^2=\mathrm{Tr}[A^1II^TA^1]`$, we need to know the variance of the gradient of $`I`$. By discretizing both the space of stimuli and possible projections $`x`$, and assuming that the probability of generating a spike is independent for different bins, we estimate $`I_iI_jA_{ij}/(N_{\mathrm{spike}}\mathrm{ln}2)`$. Therefore an expected error in the reconstruction of the optimal filter is inversely proportional to the number of spikes and is given by:

$$1𝐯\widehat{e}_1\frac{1}{2}\delta 𝐯^2=\frac{\mathrm{Tr}[A^1]}{2N_{\mathrm{spike}}\mathrm{ln}2}$$
(15)
In Fig. 2 we plot the average projection of the normalized reconstructed vector $`𝐯`$ on the RF $`\widehat{e}_1`$, and show that it scales correctly with the number of spikes.

### 3.3 A model complex cell

A sequence of spikes from a model cell with two relevant directions was simulated by projecting each of the stimuli on vectors that differ by $`\pi /2`$ in their spatial phase, taken to mimic properties of complex cells, as in Fig. 3. A particular frame leads to a spike according to a logical OR, that is if either $`s_1=𝐬\widehat{e}_1`$, $`s_1`$, $`s_2=𝐬\widehat{e}_2`$, or $`s_2`$ exceeds a threshold value $`s_t`$ in the presence of noise. Similarly to (12),

$$\frac{P(\mathrm{spike}|𝐬)}{P(\mathrm{spike})}=f(s_1,s_2)=\theta (|s_1|s_t\xi _1)\theta (|s_2|s_t\xi _2),$$
(16)
where $`\xi _1`$ and $`\xi _2`$ are independent Gaussian variables. The sampling of this input–output function by our particular set of natural stimuli is shown in Fig. 3(c). Some, especially large, combinations of values of $`s_1`$ and $`s_2`$ are not present in the ensemble. As is well known, reverse correlation fails in this case because the spike–triggered average stimulus is zero, although with Gaussian stimuli the spike–triggered covariance method would recover the relevant dimensions. Here we show that searching for maximally informative dimensions allows us to recover the relevant subspace even under more natural stimulus conditions.

We start by maximizing information with respect to one direction. Contrary to the result Fig. 1(e) for a simple cell, one optimal direction recovers only about 60% of the total information per spike \[Eq. (4)\]. Perhaps surprisingly, because of the strong correlations in natural scenes, even projection onto a random vector in the $`D10^3`$ dimensional stimulus space has a high probability of explaining 60% of total information per spike. We therefore go on to maximize information with respect to two directions. An example of the reconstruction of input–output function of a complex cell is given in Fig. 3. Vectors $`𝐯_1`$ and $`𝐯_2`$ that maximize $`I(𝐯_1,𝐯_2)`$ are not orthogonal, and are also rotated with respect to $`\widehat{e}_1`$ and $`\widehat{e}_2`$. However, the quality of reconstruction is independent of a particular choice of basis with the RS. The appropriate measure of similarity between the two planes is the dot product of their normals. In the example of Fig. 3, $`\widehat{n}_{(\widehat{e}_1,\widehat{e}_2)}\widehat{n}_{(𝐯_1,𝐯_2)}0.8`$.

Maximizing information with respect to two directions requires a significantly slower cooling rate, and consequently longer computational times. However, the expected error in the reconstruction, $`1\widehat{n}_{(\widehat{e}_1,\widehat{e}_2)}\widehat{n}_{(𝐯_1,𝐯_2)}`$, follows a $`N_{\mathrm{spike}}^1`$ behavior, similarly to (15), and is roughly twice that for a simple cell given the same number of spikes.

## 4 Remarks

In conclusion, features of the stimulus that are most relevant for generating the response of a neuron can be found by maximizing information between the sequence of responses and the projection of stimuli on trial vectors within the stimulus space. Calculated in this manner, information becomes a function of direction in a stimulus space. Those directions that maximize the information and account for the total information per response of interest span the relevant subspace. This analysis allows the reconstruction of the relevant subspace without assuming a particular form of the input–output function. It can be strongly nonlinear within the relevant subspace, and is to be estimated from experimental histograms. Most importantly, this method can be used with any stimulus ensemble, even those that are strongly non–Gaussian as in the case of natural images.

## Acknowledgments

We thank K. D. Miller for many helpful discussions. Work at UCSF was supported in part by the Sloan and Swartz Foundations and by a training grant from the NIH. Our collaboration began at the Marine Biological Laboratory in a course supported by grants from NIMH and the Howard Hughes Medical Institute.


text:
# Stokes operator squeezed continuous variable polarization states

## I introduction

The quantum properties of the polarization of continuous wave light are of increasing interest since they offer new opportunities for communicating quantum information with light and for transferring quantum information from atoms to photons and vice versa. In the single photon regime the quantum polarization states have been vigorously studied, theoretically and experimentally, with investigations of fundamental problems of quantum mechanics, such as Bell’s inequality Bell-theo ; Bell-exp , and of potential applications such as quantum cryptography QCr-theo ; QCr-exp . In comparison, continuous variable quantum polarization states have received little attention. Recently however, due to their apparent usefulness to quantum communication schemes, interest in them has been growing and a number of theoretical papers have been published APu89 KLLRS02 .

Continuous variable quantum polarization states can be carried by a bright laser beam, providing high bandwidth capabilities and therefore faster signal transfer rates than single photon systems. In addition, several proposals have been made for quantum networks that consist of spatially separated nodes of atoms whose spin states enable the storage and processing of information, connected by optical quantum communication channels DiVin95 KuzPolzik00 . Mapping of quantum states from photonic to atomic media is a crucial element in these networks. For continuous variable polarization states this mapping has been experimentally demonstrated HSSP99 . Even one of the fundamental phenomena of quantum physics, entanglement, has been realized for the macroscopic spin states of two gas samples JKP01 . Very recently entanglement was experimentally demonstrated for optical continuous variable polarization states for the first time BTSL02 .

Several methods for generating continuous variable polarization squeezed states have been proposed, most using nonlinearity provided by Kerr-like media and optical solitons KCh96 ; AAC98 ; KLLRS02 . The two experimental demonstrations previous to our work reported here and in BSBL02 , however, were achieved by combining a dim quadrature squeezed beam with a bright coherent beam on a polarizing beam splitter HSSP99 ; GSYL87 . In both cases only the properties of the state relevant to the experimental outcome were characterized. The full characterization of a continuous variable polarization state requires measurements of the fluctuations in both the orientation, and the length of the Stokes vector on a Poincaré sphere.

In this paper we present the complete experimental characterization of the Stokes vector fluctuations for four different quantum polarization states. We make use of ideas recently published by Korolkova et al. KLLRS02 . Their concept of squeezing more than one Stokes operator of a laser beam and a simple scheme to measure the Stokes operator variances are realized. Our results given in BSBL02  are extended and discussed in more detail. Experimental data from polarization squeezed states generated from a single quadrature squeezed beam and from two quadrature squeezed beams are compared.

The outline of this paper is as follows. We present a description of the theory involved in our experiments. Since polarization states can be decomposed into two mode quadrature states a general link between Stokes operator variances and quadrature variances is given. In the experimental section we characterize the polarization fluctuations of a single amplitude squeezed beam from an optical parametric amplifier (OPA). It can be seen that only the fluctuations of the Stokes vector length are below that of a coherent beam (ie. squeezed). Grangier et al. GSYL87  and Hald et al. HSSP99  converted this to squeezing of the Stokes vector orientation by combining the quadrature squeezed beam with a much brighter coherent beam on a polarizing beam splitter. We experimentally generate this situation and indeed show that the Stokes vector orientation is squeezed. This result is compared with measurements on polarization states generated from two quadrature squeezed beams. Two bright amplitude or phase squeezed beams from two independent OPAs are overlapped on a polarizing beam splitter KLLRS02 ; BSBL02  demonstrating “pancake-like” and “cigar-like” uncertainty volumes on the Poincaré sphere for phase and amplitude squeezed input beams, respectively. Both the orientation and the length of the Stokes vector were squeezed for the “cigar-like” uncertainty volume. In the final section several schemes for encoding information on continuous variable polarization states of light are discussed. The conventional fiber-optic communication protocol is compared with optimized coherent beam and squeezed beam protocols. We show that the channel capacity of the “cigar-like” polarization squeezed states exceeds the channel capacity of all the other states.

## II Theoretical background

The polarization state of a light beam in classical optics can be visualized as a Stokes vector on a Poincaré sphere (Fig. 1) and is determined by the four Stokes parameters Sto52 : $`S_0`$ represents the average beam intensity whereas $`S_1`$, $`S_2`$, and $`S_3`$ characterize its polarization and form a Cartesian axes system. If the Stokes vector points in the direction of $`S_1`$, $`S_2`$, or $`S_3`$ the polarized part of the beam is horizontally, linearly at 45, or right-circularly polarized, respectively. Two beams are said to be opposite in polarization and do not interfere if their Stokes vectors point in opposite directions. The quantity $`S=(S_1^2+S_2^2+S_3^2)^{1/2}`$ is the radius of the classical Poincaré sphere and describes the average intensity of the polarized part of the radiation. The fraction $`S/S_0`$ ($`\mathrm{\hspace{0.17em}0}<S/S_0<1`$) is called the degree of polarization. For quasi-monochromatic laser light which is almost completely polarized $`S_0`$ is a redundant parameter, completely determined by the other three parameters ($`S_0=S`$ in classical optics). All four Stokes parameters are accessible from the simple experiments shown in Fig. 2.

An equivalent representation of polarization states of light is given by the 4 elements of the coherence matrix (Jones matrix). The relations between these elements and the Stokes parameters can be found in WolfMandel . In contrast to the coherence matrix elements the Stokes parameters are observables and therefore can be associated with Hermitian operators. Following JauchRohrlich76  and Robson74  we define the quantum-mechanical analogue of the classical Stokes parameters for pure states in the commonly used notation:

$`\widehat{S}_0`$ $`=`$ $`\widehat{a}_H^{}\widehat{a}_H+\widehat{a}_V^{}\widehat{a}_V,`$
$`\widehat{S}_1`$ $`=`$ $`\widehat{a}_H^{}\widehat{a}_H\widehat{a}_V^{}\widehat{a}_V,`$ (1)
$`\widehat{S}_2`$ $`=`$ $`\widehat{a}_H^{}\widehat{a}_Ve^{i\theta }+\widehat{a}_V^{}\widehat{a}_He^{i\theta },`$
$`\widehat{S}_3`$ $`=`$ $`i\widehat{a}_V^{}\widehat{a}_He^{i\theta }i\widehat{a}_H^{}\widehat{a}_Ve^{i\theta },`$
where the subscripts $`H`$ and $`V`$ label the horizontal and vertical polarization modes respectively, $`\theta `$ is the phase shift between these modes, and the $`\widehat{a}_{H,V}`$ and $`\widehat{a}_{H,V}^{}`$ are annihilation and creation operators for the electro-magnetic field in frequency space footnote1 .

The commutation relations of the annihilation and creation operators

$$[\widehat{a}_k,\widehat{a}_l^{}]=\delta _{kl},\text{with}k,l\{H,V\},$$
(2)
directly result in Stokes operator commutation relations,

$$[\widehat{S}_1,\widehat{S}_2]=2i\widehat{S}_3,[\widehat{S}_2,\widehat{S}_3]=2i\widehat{S}_1,[\widehat{S}_3,\widehat{S}_1]=2i\widehat{S}_2.$$
(3)
Apart from a normalization factor, these relations are identical to the commutation relations of the Pauli spin matrices. In fact the three Stokes parameters in Eq. (3) and the three Pauli spin matrices both generate the special unitary group of symmetry transformations SU(2) of Lie algebra Kaku93 . Since this group obeys the same algebra as the three-dimensional rotation group, distances in three dimensions are invariant. Accordingly the operator $`\widehat{S}_0`$ is also invariant and commutes with the other three Stokes operators ($`[\widehat{S}_0,\widehat{S}_j]=0,\text{with}j=1,2,3`$). The non-commutability of the Stokes operators $`\widehat{S}_1`$, $`\widehat{S}_2`$ and $`\widehat{S}_3`$ precludes the simultaneous exact measurement of their physical quantities. As a direct consequence of Eq.(3) the Stokes operator mean values $`\widehat{S}_j`$ and their variances $`V_j=\widehat{S}_j^2\widehat{S}_j^2`$ are restricted by the uncertainty relations JauchRohrlich76

$$V_1V_2|\widehat{S}_3|^2,V_2V_3|\widehat{S}_1|^2,V_3V_1|\widehat{S}_2|^2.$$
(4)
In general this results in non-zero variances in the individual Stokes parameters as well as in the radius of the Poincaré sphere (see Fig. 1b)). The quantum noise in the Stokes parameters even effects the definitions of the degree of polarization APu89 ; AAC98  and the Poincaré sphere radius. It can be shown from Eqs. (II) and (2) that the quantum Poincaré sphere radius is different from its classical analogue, $`\widehat{S}=\widehat{S}_0^2+2\widehat{S}_0^{1/2}`$.

Recently it has been shown that the Stokes operator variances may be obtained from the frequency spectrum of the electrical output currents of the setups shown in Fig. 2 KLLRS02 . To calculate the Stokes operator variances we use the linearized formalism here. The creation and annihilation operators are expressed as sums of real classical amplitudes $`\alpha _{H,V}`$ and quantum noise operators $`\delta \widehat{a}_{H,V}`$ WallsMilburn

$$\widehat{a}_{H,V}=\alpha _{H,V}+\delta \widehat{a}_{H,V}.$$
(5)
The operators in Eq.(5) are non-hermitian and therefore non-physical. To express the Stokes operators of Eq. (II) in terms of hermitian operators we define the generalized quadrature quantum noise operators $`\delta \widehat{X}_{H,V}(\xi )`$

$`\delta \widehat{X}_{H,V}(\xi )=`$ $`\delta \widehat{a}_{H,V}^{}e^{i\xi }+\delta \widehat{a}_{H,V}e^{i\xi },`$ (6)
$`\delta \widehat{X}_{H,V}(\xi =0)=`$ $`\delta \widehat{X}_{H,V}^+=\delta \widehat{a}_{H,V}^{}+\delta \widehat{a}_{H,V},`$ (7)
$`\delta \widehat{X}_{H,V}(\xi =\pi /2)=`$ $`\delta \widehat{X}_{H,V}^{}=i(\delta \widehat{a}_{H,V}^{}\delta \widehat{a}_{H,V}).`$ (8)
$`\xi `$ is the phase of the quantum mechanical oscillator and $`\delta \widehat{X}_{H,V}^+`$ and $`\delta \widehat{X}_{H,V}^{}`$ are the amplitude quadrature noise operator and the phase quadrature noise operator respectively.

If the variances of the noise operators are much smaller than the coherent amplitudes then a first order approximation of the noise operators is appropriate. This yields the Stokes operator mean values

$`\widehat{S}_0`$ $`=`$ $`\alpha _H^2+\alpha _V^2=\widehat{n},`$
$`\widehat{S}_1`$ $`=`$ $`\alpha _H^2\alpha _V^2,`$ (9)
$`\widehat{S}_2`$ $`=`$ $`2\alpha _H\alpha _V\text{cos}\theta ,`$
$`\widehat{S}_3`$ $`=`$ $`2\alpha _H\alpha _V\text{sin}\theta ,`$
These expressions are identical to the Stokes parameters in classical optics. Here $`\widehat{n}`$ is the expectation value of the photon number operator. For a coherent beam the expectation value and variance of $`\widehat{n}`$ have the same magnitude, this magnitude equals the conventional shot-noise level. The variances of the Stokes parameters are given by

$`V_0=`$ $`\alpha _H^2(\delta \widehat{X}_H^+)^2+\alpha _V^2(\delta \widehat{X}_V^+)^2+2\alpha _H\alpha _V\delta \widehat{X}_H^+\delta \widehat{X}_V^+,`$
$`V_1=`$ $`\alpha _H^2(\delta \widehat{X}_H^+)^2+\alpha _V^2(\delta \widehat{X}_V^+)^22\alpha _H\alpha _V\delta \widehat{X}_H^+\delta \widehat{X}_V^+,`$
$`V_2(\theta )`$ $`=\alpha _H^2(\delta \widehat{X}_V(\theta ))^2+\alpha _V^2(\delta \widehat{X}_H(\theta ))^2,`$
$`+2\alpha _H\alpha _V\delta \widehat{X}_V(\theta )\delta \widehat{X}_H(\theta ),`$
$`V_3(\theta )`$ $`=V_2(\theta {\displaystyle \frac{\pi }{2}}).`$
It can be seen from Eqs.(II) that the variances of Stokes operators can be expressed in terms of the variances of quadrature operators of two modes. The polarization squeezed state can then be defined in a straight forward manner. The variances of the noise operators in the above equation are normalized to one for a coherent beam. Therefore the variances of the Stokes parameters of a coherent beam are all equal to the shot-noise of the beam. For this reason a Stokes parameter is said to be squeezed if its variance falls below the shot-noise of an equal power coherent beam. Although the decomposition to the $`H,V`$-polarization axis of Eqs. (II) is independent of the actual procedure of generating a polarization squeezed beam, it becomes clear that two overlapped quadrature squeezed beams can produce a single polarization squeezed beam. If two beams in the horizontal and vertical polarization mode having uncorrelated quantum noise are used Eqs. (II) can be rewritten as

$`V_0=`$ $`V_1=\alpha _H^2(\delta \widehat{X}_H^+)^2+\alpha _V^2(\delta \widehat{X}_V^+)^2,`$
$`V_2(\theta )=`$ $`\text{cos}^2\theta \left(\alpha _V^2(\delta \widehat{X}_H^+)^2+\alpha _H^2(\delta \widehat{X}_V^+)^2\right)`$ (11)
$`+\text{sin}^2\theta \left(\alpha _V^2(\delta \widehat{X}_H^{})^2+\alpha _H^2(\delta \widehat{X}_V^{})^2\right),`$
$`V_3(\theta )=`$ $`V_2(\theta {\displaystyle \frac{\pi }{2}}).`$
Here we choose the amplitude and the phase quadrature noise operators to express the variances. This corresponds to our actual experimental setup where either amplitude or the phase quadratures were squeezed. It can be seen from Eqs.(II) that in a polarization squeezed beam generated from two amplitude squeezed beams $`\widehat{S}_0`$ and two additional Stokes parameters can in theory be perfectly squeezed while the fourth is anti-squeezed if specific angles of $`\theta =0`$ or $`\theta =\pi /2`$ are used. Utilizing only one squeezed beam it is not possible to simultaneously squeeze any two of $`\widehat{S}_1`$, $`\widehat{S}_2`$, and $`\widehat{S}_3`$ to quieter than 3 dB below shot-noise ($`V_i+V_j\widehat{n}`$, with $`i,j\{1,2,3;ij\}`$).

## III Experiment

Prior to our work presented here and in BSBL02 , polarization squeezed states were generated by combining a strong coherent beam with a single weak amplitude squeezed beam GSYL87 ; HSSP99 . In both of those experiments the variance of only one Stokes parameter was determined, and therefore the polarization state was not fully characterized. In this paper we experimentally characterize the mean and variance of all four Stokes operators for these states. We extend the work to polarization squeezed states produced from two amplitude/phase squeezed beams. Fig. 3 shows our experimental setup.

### III.1 Generation of quadrature squeezed light

We produced the two quadrature squeezed beams required for this experiment in a pair of OPAs. Each OPA was an optical resonator consisting of a hemilithic MgO:LiNbO3 crystal and an output coupler. The reflectivities of the output coupler were 96% and 6% for the fundamental (1064 nm) and the second harmonic (532 nm) laser modes, respectively. Each OPA was pumped with single-mode 532 nm light generated by a 1.5 W Nd:YAG non-planar ring laser and frequency doubled in a second harmonic generator (SHG). The SHG was of identical structure to the OPAs but with 92% reflectivity at 1064 nm. The OPAs were seeded with 1064 nm light after spectral filtering in a modecleaner. The refractive indices of the MgO:LiNbO3 crystals in each resonator was modulated with an RF field, this provided error signals on the reflected seed power that were used to lock their lengths. The modulation also resulted in a phase modulation on the output beams from the SHG and each OPA. The coherent amplitude of each OPAs output was a deamplified/amplified version of the seed coherent amplitude; the level of amplification was dependent on the phase difference between pump and seed ($`\varphi _{\mathrm{SH}}`$). Therefore the second harmonic pump phase modulation resulted in a modulation of the amplification of the OPAs. Error signals could be extracted from this effect, enabling the relative phase between pump and seed to be locked. Locking to deamplification or amplification provided an amplitude or phase squeezed output, respectively. Typical measured variance spectra of the two locked quadrature squeezed beams are shown in Fig. 4. Since the squeezed states were carried by bright laser beams of approximately 1 mW, the noise reduction was degraded at lower frequencies due to the laser relaxation oscillation. At higher frequencies the squeezed spectrum was limited by the bandwidth of the OPAs.

Past experiments requiring two quadrature squeezed beams commonly used a single ring resonator with two outputs Ou ; with two independent OPAs the necessary intracavity pump power is halved, this reduces the degradation of squeezing due to green-induced infrared absorption (GRIIRA) Furukawa .

### III.2 Measuring the Stokes operators

Instantaneous values of the Stokes operators of all polarization states analyzed in this paper were obtained with the apparatus shown in Fig. 2. The uncertainty relations of Eqs. 4 dictate that $`\widehat{S}_1`$, $`\widehat{S}_2`$, and $`\widehat{S}_3`$ cannot be measured simultaneously. The beam under interrogation was split on a polarizing beam splitter and the two outputs were detected on a pair of high quantum efficiency photodiodes with 30 MHz bandwidth; the resulting photocurrents were added and subtracted to yield photocurrents containing the instantaneous values of $`\widehat{S}_0`$ and $`\widehat{S}_1`$. To measure $`\widehat{S}_2`$ the polarization of the beam was rotated by 45 with a half-wave plate before the polarizing beam splitter and the detected photocurrents were subtracted. To measure $`\widehat{S}_3`$ the polarization of the beam was again rotated by 45 with a half-wave plate and a quarter-wave plate was introduced before the polarizing beam splitter such that a horizontally polarized input beam became right-circularly polarized. Again the detected photocurrents were subtracted. The expectation value of each Stokes operator was equal to the DC output of the detection device and the variance was obtained by passing the output photocurrent into a Hewlett-Packard E4405B spectrum analyzer. Every polarization state interrogated in this work had a total power $`\widehat{S}_0`$ of roughly 2 mW.

An accurate shot-noise level was required to determine whether any given Stokes operator was squeezed. This was measured by operating a single OPA without the second harmonic pump. The seed power was adjusted so that the output power was equal to that of the beam being interrogated. In this configuration the detection setup for $`\widehat{S}_2`$ (see Fig. 2c)) functions exactly as a homodyne detector measuring vacuum noise scaled by the OPA output power, the variance of which is the shot-noise. Throughout each experimental run the power was monitored and was always within 2% of the power of the coherent calibration beam. This led to a conservative error in our frequency spectra of $`\pm `$0.05 dB.

The Stokes operator variances reported in this paper were taken over the range from 3 to 10 MHz. The darknoise of the detection apparatus was always more than 4 dB below the measured traces and was taken into account. Each displayed trace is the average of three measurement results normalized to the shot-noise and smoothed over the resolution bandwidth of the spectrum analyzer which was set to 300 kHz. The video bandwidth of the spectrum analyzer was set to 300 Hz.

As is the case for all continuous variable quantum optical experiments, the efficiency of the Stokes operator measurements was critical. The overall detection efficiency of the interrogated beams was 76%. The loss came primarily from three sources: loss in escape from the OPAs (14%), detector inefficiency (7%), and loss in optics (5%). In the experiment where a squeezed beam was overlapped with a coherent beam additional loss was incurred due to poor mode-matching between the beams and the detection efficiency was 71%. Depolarizing effects are thought to be another significant source of loss for some polarization squeezing proposals APu89 . In our scheme the non-linear processes (OPAs) are divorced from the polarization manipulation (wave plates and polarizing beam splitters), and depolarizing effects are insignificant.

### III.3 Quantum polarization states from a single squeezed beam

We first characterize the polarization state of a single bright amplitude squeezed beam provided by one of our OPAs, as shown in Fig. 5. The squeezed beam was horizontally polarized, resulting in Stokes operator expectation values of $`\widehat{S}_0=\widehat{S}_1=|\alpha _H|^2`$ and $`\widehat{S}_2=\widehat{S}_3=0`$. The variance spectra of the operators were measured and are displayed in Fig. 6. $`\widehat{S}_0`$ and $`\widehat{S}_1`$ were squeezed since the horizontally polarized amplitude squeezed beam hit only one detector in this detector setup. For the measurements of $`\widehat{S}_2`$ and $`\widehat{S}_3`$ the beam intensity was divided equally between the two detectors. The electronic subtraction yielded vacuum noise scaled by the beam intensity, thus both variance measurements were at the shot noise level. It is apparent from these measurements that only the length of the Stokes vector is well determined; the orientation is just as uncertain as it would be for a coherent state.

To obtain squeezing of the orientation of the Stokes vector Grangier et al. GSYL87  and Hald et al. HSSP99  overlapped a dim quadrature squeezed beam with a bright orthogonally polarized coherent beam. We consider this situation next (as shown in Fig. 8). Since two beams are now involved, the relative phase $`\theta `$ becomes important. A DC and an RF error-signal, both dependant on $`\theta `$, were extracted from the Stokes operator measurement device. Together, these error signals allowed us to lock $`\theta `$ to either 0 or $`\pi /2`$ rads in all of the following experiments. We mixed a bright horizontally polarized coherent beam with a dim vertically polarized amplitude squeezed beam. Since the horizontally polarized beam was much more intense than the vertically polarized beam, the Stokes operator expectation values became $`\widehat{S}_0\widehat{S}_1|\alpha _H|^2`$ and $`\widehat{S}_2\widehat{S}_3=0`$. The Stokes operator variances obtained for this polarization state are shown in Fig. 8, here $`\widehat{S}_2`$ is anti-squeezed and $`\widehat{S}_3`$ is squeezed. The variances of $`\widehat{S}_0`$ and $`\widehat{S}_1`$ were slightly above the shot-noise level because of residual noise from our laser resonant relaxation oscillation. The experiment carried out with $`\theta `$ locked to 0 rads is not shown, in this case the measured variances of $`\widehat{S}_2`$ and $`\widehat{S}_3`$ were swapped. In fact the Stokes vector was still pointing along $`\widehat{S}_1`$ but the quantum noise was rotated on the Poincaré sphere (see Fig. 13b).

### III.4 Quantum polarization states from two quadrature squeezed beams

The two experiments described in Section III.3 demonstrated how it is possible to squeeze the length and orientation of the Stokes vector. In this section we demonstrate that it is possible to do both simultaneously. The two quadrature squeezed beams produced in our OPAs were combined with orthogonal polarization on a polarizing beam splitter KLLRS02  as shown in Fig. 9. This produced an output beam with Stokes parameter variances as given by Eqs. (II). Both input beams had equal power ($`\alpha _H=\alpha _V=\alpha /\sqrt{2}`$) and both were squeezed in the same quadrature. The Stokes parameters and their variances were again determined as shown in Fig. 2. The relative phase between the quadrature squeezed input beams $`\theta `$ was locked to $`\pi /2`$ rads producing a right-circularly polarized beam with Stokes parameter means of $`\widehat{S}_1=\widehat{S}_2=0`$ and $`\widehat{S}_0=\widehat{S}_3=|\alpha |^2`$.

First both OPA pump beams were phase locked to amplification, this produced two phase squeezed beams. Fig. 10 shows the measurement results obtained; $`\widehat{S}_0`$, $`\widehat{S}_1`$ and $`\widehat{S}_3`$ were all anti-squeezed and $`\widehat{S}_2`$ was squeezed throughout the range of the measurement. The optimum noise reduction of $`\widehat{S}_2`$ was 2.8 dB below shot-noise and was observed at 4.8 MHz. Our OPAs are particularly sensitive to phase noise coupling in from the MgO:LiNbO3 crystals. We attribute the structure in the frequency spectra of $`\widehat{S}_2`$ and the poorer squeezed observed here, to this. Apart from this structure, these results are very similar to those produced by a single squeezed beam and a coherent beam; the orientation of the Stokes vector is squeezed. However, here the uncertainty in the length of the Stokes vector is greater than for a coherent state so the polarization state, although produced from two quadrature squeezed beams, is actually less certain.

Fig. 11 shows the measurement results obtained with the OPAs locked to deamplification. Therefore both OPAs provided amplitude squeezed beams. Again we interrogated the combined beams and found $`\widehat{S}_0`$, $`\widehat{S}_1`$ and $`\widehat{S}_3`$ all to be squeezed from 4.5 MHz to the limit of our measurement, 10 MHz. $`\widehat{S}_2`$ was anti-squeezed throughout the range of the measurement. Between 7.2 MHz and 9.6 MHz $`\widehat{S}_0`$, $`\widehat{S}_1`$ and $`\widehat{S}_3`$ were all more than 3 dB below shot-noise. The squeezing of $`\widehat{S}_0`$ and $`\widehat{S}_3`$ was degraded at low frequency due our lasers resonant relaxation oscillation. Since this noise was correlated it canceled in the variance of $`\widehat{S}_1`$. The maximum squeezing of $`\widehat{S}_0`$ and $`\widehat{S}_2`$ was 3.8 dB and 3.5 dB respectively and was observed at 9.3 MHz. The maximum squeezing of $`\widehat{S}_1`$ was 4.3 dB at 5.7 MHz. The repetitive structure at 4, 5, 6, 7, 8 and 9 MHz was caused by electrical pick-up in our SHG resonator emitted from a separate experiment operating in the laboratory. In this case both the orientation and the length of the Stokes vector were squeezed.

Finally we point out that the variances of $`\widehat{S}_1`$ in Figs. 8, 10 and 11 were all squeezed at frequencies down to 3 MHz and even below, whereas Fig. 4 shows a clear degradation below 4 MHz. This improved performance is due to electrical noise cancellation of correlated laser relaxation oscillation noise. This noise is effectively reduced by taking the difference of the two photo currents in our detector setup to measure $`V_1`$.

## IV Visualization of quantum correlations in continuous variable polarization states

In this section measured quantum correlations in polarization states at 8.5 MHz are visualized. Based on the theoretical formalism in Section II continuous variable polarization states can be characterized by the measurement of Stokes operator expectation values and variances using the setup shown in Fig. 2. Our noise measurement results at 8.5 MHz on five different states are summarized in Fig. 12. The noise characteristics of the Stokes parameters are mapped onto the coordinate system of the Poincaré sphere, assuming Gaussian noise statistics. Given this assumption, the standard deviation contour-surfaces shown here provide an accurate representation of the states three-dimensional noise distribution. The quantum polarization noise of a coherent state forms a sphere of noise as portrayed in Fig. 12 a). The noise volumes b) to e) visualizes the measurements on a single bright amplitude squeezed beam, on the combination of a vacuum amplitude squeezed beam and a bright coherent beam, on two locked phase squeezed input beams and on two locked amplitude squeezed beams, respectively. In all cases the Stokes operator noise volume describes the end position of the Stokes vector pointing upwards. In b) and c) the Stokes vectors are parallel to the direction of $`S_1`$, in d) and e) parallel to the direction of $`S_3`$, since we used horizontally and right-circularly polarized light, respectively. However, there was no fundamental bias in the orientation of the quantum Stokes vector in our experiment. By varying the angle of an additional half-wave plate in the polarization squeezed beam or by varying $`\theta `$ any orientation may be achieved. In fact, as mentioned earlier our experiments were also carried out with $`\theta `$ locked to 0 rads. This had the effect of rotating the Stokes vector and its quantum noise by $`\pi /2`$ around $`\widehat{S}_1`$. Nearly identical results were obtained but on alternative Stokes parameters. Fig. 13 a) shows Poincaré sphere representations of this rotation for the polarization states produced by two amplitude squeezed beams. In Fig. 13 b) the combination of a amplitude squeezed vacuum and a bright coherent beam exemplifies that different orientations of the noise volume can be generated using appropriate combination of waveplates.

## V Channel Capacity of Polarization Squeezed Beams

The reduced level of fluctuations in polarization squeezed light can be used to improve the channel capacity of communication protocols. Let us consider information encoded on the sidebands of a bandwidth limited single spatial mode laser beam. We assume that only direct detection is employed, or in other words, that phase sensitive techniques such as homodyne measurement are not available. This is not an artificial constraint since phase sensitive techniques are technically difficult to implement and are rarely utilized in conventional optical communications systems.

An upper bound to the amount of information that can be carried by a bandwidth limited additive white Gaussian noise channel is given by the Shannon capacity C shan  in bits per dimension

$$C=\frac{1}{2}\mathrm{log}_2(1+R).$$
(12)
$`R`$ is the signal to noise ratio of the channel and is given by the ratio of the spectral variance of the signal modulation $`V_s`$ and the noise spectral variance $`V_n`$

$$R=\frac{V_s}{V_n}.$$
(13)
We wish to compare the channel capacities achieveable with pure coherent and squeezed states for a given average photon number in the sidebands $`\overline{n}`$ where $`\overline{n}=\delta \widehat{a}_H^{}\delta \widehat{a}_H+\delta \widehat{a}_V^{}\delta \widehat{a}_V`$. Note that $`\overline{n}`$ takes into account both, the signal modulation and the squeezing. An overview of quantum noise limited channel capacities may be found in YHa86  and CDr94 .

First let us consider strategies which might be employed with a coherent light beam. In conventional optical communication systems the polarization degrees of freedom are ignored completely and information is encoded only on $`\widehat{S}_0`$ as intensity fluctuations. Taking $`\alpha _V=0`$ the variance of the Stokes operator $`\widehat{S}_0`$ is given by $`V_0=\alpha _H^2(V_s+1)`$ in accordance with Eq. (II). For this one-dimensional coherent channel $`V_n=|\delta X_H^+|^2=1`$ and therefore $`R=V_s`$. For this arrangement it can be shown that the average photon number per bandwidth per second is $`\overline{n}=\frac{1}{4}V_s`$ providing a photon resource limited Shannon capacity of

$$C_{\mathrm{coh}}^i=\frac{1}{2}\mathrm{log}_2(1+4\overline{n}),$$
(14)
as a function of $`\overline{n}`$. This is a non-optimal strategy however. Examining Eqs. (3) and (II) we see that it is possible to choose an arrangement for which two of the Stokes operators commute and so can be measured simultaneously. Indeed it is easy to show that such simultaneous measurements can be made using only linear optics and direct detection. In particular let us assume that $`\alpha _H=\alpha _V`$ such that $`\widehat{S}_2`$ and $`\widehat{S}_3`$ commute, and use $`\widehat{S}_2`$ and $`\widehat{S}_3`$ as two independent information channels. Then $`\widehat{S}_2=\widehat{S}_0`$ and the information in both dimensions can be simultaneously extracted by subtracting and adding the photocurrents of the same pair of detectors. Assuming equal signal to noise ratios $`R_2=R_3`$ we find that $`V_{s2}=V_{s3}=2\overline{n}`$, and the channel capacity may be written

$`C_{\mathrm{coh}}^{ii}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{log}_2(1+R_2)+{\displaystyle \frac{1}{2}}\mathrm{log}_2(1+R_3)`$ (15)
$`=`$ $`\mathrm{log}_2(1+2\overline{n}).`$
This channel capacity is always greater than that of Eq. (14) and for large $`\overline{n}`$ is 100 % greater.

For sufficiently high $`\overline{n}`$ a further improvement in channel capacity can be achieved. Consider placing signals on all three Stokes operators. Because of the non-commutation of $`\widehat{S}_1`$ with $`\widehat{S}_2`$ and $`\widehat{S}_3`$ it is not possible to read out all three signals without a measurement penalty. Suppose the receiver adopts the following strategy: divide the beam on a beamsplitter with transmitivity $`ϵ`$ and then measure $`\widehat{S}_1`$ on the reflected output and $`\widehat{S}_2`$ and $`\widehat{S}_3`$ on the other output. Division of the beam will reduce the measured signal to noise ratios due to the injection of quantum noise at the beamsplitter such as $`R_1=(1ϵ)V_{s1}`$, $`R_2=ϵV_{s2}`$ and $`R_3=ϵV_{s3}`$. We find that for large $`\overline{n}`$ an optimum is reached with $`ϵ=2/3`$ and the signal photon number in each Stokes parameter being $`\overline{n}_1=\overline{n}_2=\overline{n}_3=\overline{n}/3`$. Hence the channel capacity is

$`C_{\mathrm{coh}}^{iii}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{log}_2(1+R_1)+{\displaystyle \frac{1}{2}}\mathrm{log}_2(1+R_2)+{\displaystyle \frac{1}{2}}\mathrm{log}_2(1+R_3)`$ (16)
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{log}_2(1+{\displaystyle \frac{4}{9}}\overline{n})+\mathrm{log}_2(1+{\displaystyle \frac{8}{9}}\overline{n}).`$
This capacity beats that of Eq. (15) for $`\overline{n}>7.56`$. In summery the optimum coherent channel capacity is given by Eq. (16) for average photon numbers $`\overline{n}>7.56`$ and by Eq. (15) for lower values, see Fig. 14 traces c) and d).

Now let us examine the effect of polarization squeezing on the channel capacity. Consider first the simple case of intensity modulation on a single squeezed beam, which is equivalent to using a squeezed beam in the first coherent case considered (case i). The channel capacity can be maximized by optimizing the fraction of photons which are introduced by squeezing the quantum noise and the residual fraction of photons which actually carries the signal. For large photon numbers we find a proportioning of 0.5. For an average photon number of 1 just 1/3 of that photon should be used to reduce the quantum noise. The maximum channel capacity for a one-dimensional squeezed channel is found to be

$`C_{1\mathrm{s}\mathrm{q}\mathrm{z}}^i`$ $`=`$ $`\mathrm{log}_2(1+2\overline{n}),`$ (17)
which was previously given in YHa86 . This capacity beats only the corresponding coherent state. It is as efficient as the two-dimensional coherent channel, but less efficient than the three-dimensional coherent one for large photon numbers.

Consider now a polarization squeezed beam which is produced from a minimum uncertainty squeezed beam and a coherent beam, as in Section III.3. Suppose, as in case ii, that $`\widehat{S}_2`$ and $`\widehat{S}_3`$ commute and arrange that $`\widehat{S}_2`$ has fluctuations at the quantum noise level while $`\widehat{S}_3`$ is optimally squeezed. Again signals are encoded on $`\widehat{S}_2`$ and $`\widehat{S}_3`$. The channel capacity can be maximized by adjusting the relative signal sizes on the two Stokes operators for fixed average photon number as a function of the squeezing. Here 1/3 of the photons is used to squeeze and 2/3 is split equally for the two dimensions. The resultant maximum channel capacity is

$`C_{1\mathrm{s}\mathrm{q}\mathrm{z}}^{ii}`$ $`=`$ $`{\displaystyle \frac{3}{2}}\mathrm{log}_2(1+{\displaystyle \frac{4}{3}}\overline{n}).`$ (18)
This always beats all three coherent state cases considered here, but in the limit of large $`\overline{n}`$ the advantage is minimal since the scaling with photon number is the same as that of $`C_{\mathrm{coh}}^{iii}`$ in Eq. (16).

If the polarization squeezed beam is produced from two amplitude squeezed beams as in Section III.4 the enhancement becomes more significant. Suppose again that $`\widehat{S}_2`$ and $`\widehat{S}_3`$ commute but now that both are optimally squeezed. Again encoding on $`\widehat{S}_2`$ and $`\widehat{S}_3`$, and varying the signal strength as a function of squeezing to maximize the channel capacity for a given $`\overline{n}`$. The maximum is reached when the photons are used to squeeze the noise and transport information in equal shares. The channel capacity for this arrangement is given by

$`C_{2\mathrm{s}\mathrm{q}\mathrm{z}}^{ii}`$ $`=`$ $`2\mathrm{log}_2(1+\overline{n}),`$ (19)
which for large $`\overline{n}`$ is 33% greater than both the optimum coherent scheme and the scheme using a single quadrature squeezed beam. No further improvement of the channel capacity can be obtained by encoding the information on three Stokes parameters, as in case iii. Optimization of the beam splitter reflectivity results in the already considered two-dimensional arrangement. This is not a surprising result since the third Stokes parameter is anti-squeezed. Fig. 14 summarizes our results.

Finally we assess the channel capacities that could in principle be achieved using the polarization squeezed state generated in our experiment from two amplitude squeezed beams. The polarization squeezing achieved in Fig. 11 implies that, in the frequency range of 8 MHz to 10 MHz, 0.17 side-band photons per bandwidth and second were present in each of the two dimensions. This is an optimum quantum resource to transmit 0.68 side-band photons. Signals sufficiently high above detector dark noise would achieve a channel capacity that is around 21% greater than the ideal channel capacity achievable from a coherent beam with the same average side-band photon number (see crosses in Fig. 14).

## VI Conclusion

The field of quantum communication and computation is receiving much attention. The continuous variable polarization states investigated here, are one of the most promising candidates for carrying the information in a quantum network. In this paper we have characterized the non classical properties of these states on the basis of the Stokes operators and their variances. Different classes of polarization squeezed states have been generated and experimentally characterized. We compared the coherent polarization state in Fig. 12a) with squeezed polarization states generated from a single amplitude squeezed beam Fig. 12c) and from two amplitude squeezed beams Fig. 12e), and proved that squeezing of better than 3 dB of three Stokes parameters ($`\widehat{S}_0`$, $`\widehat{S}_1`$, and $`\widehat{S}_3`$) simultaneously is possible only in the latter case. We have theoretically analyzed the channel capacity for several communication protocols using continuous variable polarization states. For a given average photon number $`\overline{n}`$, we found the polarization state produced from two quadrature squeezed states can provide a 33% greater channel capacity than both the optimum coherent scheme and the scheme using a single quadrature squeezed beam.

We acknowledge the Alexander von Humboldt foundation for support of R. Schnabel; the Australian Research Council for financial support. This work is a part of EU QIPC Project, No. IST-1999-13071 (QUICOV).


text:
# Self-aligned fabrication process for silicon quantum computer devices

## Acknowledgement

This work was supported by the Australian Research Council, the Australian Government and by the U.S. National Security Agency (NSA), Advanced Research and Development Activity (ARDA) and the Army Research Office (ARO) under contract number DAAD19-01-1-0653.