format: text

text:
# Exact c-number Representation of Non-Markovian Quantum Dissipation

## Abstract

The reduced dynamics of a quantum system interacting with a linear heat bath finds an exact representation in terms of a stochastic Schrödinger equation. All memory effects of the reservoir are transformed into noise correlations and mean-field friction. The classical limit of the resulting stochastic dynamics is shown to be a generalized Langevin equation, and conventional quantum state diffusion is recovered in the Born–Markov approximation. The non-Markovian exact dynamics, valid at arbitrary temperature and damping strength, is exemplified by an application to the dissipative two-state system.

Irreversible quantum processes are important in almost every field of condensed-matter physics and chemistry. In phenomena as different as diffusion of light particles in solids and light harvesting in biological systems, the dissipation of energy, the destruction of phase coherence, and the generation of entropy play a key role. The question how such processes can be described by the stochastic propagation of quantum states has recently attracted increasing attention . For classical systems with linear dissipation, Langevin equations provide a theoretical (and numerical) tool to accurately describe the interaction of a system with a complex thermal reservoir in comparably simple terms of stochastic forces and memory friction. For open quantum systems, no such simple and exact approach in terms of an equation of motion has been established so far. The dynamical quantity of interest is the reduced density matrix, obtained by tracing out the reservoir degrees of freedom from the dynamics. Significant information about correlations between system and reservoir is lost in this operation, hence no deterministic differential equation of motion can be derived for the reduced density matrix without approximations. Traditional approaches treat these correlations perturbatively, yielding approximate equations of motion such as Redfield or Master equations. While these have been used very successfully in the fields of quantum optics and magnetic resonance, there are many condensed-matter problems where they are qualitatively wrong. Large coupling constants and long correlation timescales both need to be treated non-perturbatively. This can be accomplished in a formally exact manner by path integrals for open quantum systems, but an analytic evaluation of these functional integrals is usually restricted to special cases and approximations, e.g., the semiclassical limit. Furthermore, present exact numerical techniques such as Quantum Monte Carlo (QMC) methods  have to cope with a sign problem arising from the fact that probability amplitudes with varying phases rather than probabilities have to be added in a quantum computation.

Quantum State diffusion (QSD)  has been established as an alternative theory of quantum dissipation in the perturbative regime. Through the stochastic propagation of pure quantum states it yields an intuitively appealing picture of the behavior of individual quantum trajectories in open systems and permits effective numerical calculations. Generalizing this formalism to the case of non-perturbative, non-Markovian dynamics holds the promise to overcome most of the above-mentioned limitations of currently known techniques. This approach has recently been taken by Diósi, Strunz, and Gisin , but at the price of incorporating non-Markovian retardation effects in a memory functional. In practice this means that a general solution of the resulting equation of motion, even numerically, is almost as elusive as that of the underlying path integral.

In this Letter we derive and discuss stochastic Schrödinger equations for open systems which allow the treatment of linear dissipation of arbitrary strength and correlation time scales. Without approximation, all effects of the system-reservoir interaction are re-cast in the form of c-number noise and friction forces with a suggestive physical interpretation: The ensemble of phase-space points described by a classical Langevin equation is generalized to a stochastic ensemble of quantum states. Both descriptions are linked by a correspondence principle for open systems.

Any general theory of quantum dissipation has to start from an open system embedded in a large reservoir system, whose degrees of freedom are treated fully quantum mechanically, but can later be eliminated from the dynamics. The Hamiltonian of such a model consists of system and reservoir terms and an interaction potential,

$$H=H_0+H_R+H_I.$$
(1)
Let us consider a thermal correlation function

$$AB(t)=tr\left\{\mathrm{e}^{\beta H}A\mathrm{e}^{\frac{\mathrm{i}}{\mathrm{}}Ht}B\mathrm{e}^{\frac{\mathrm{i}}{\mathrm{}}Ht}\right\}/tr\mathrm{e}^{\beta H},$$
(2)
where $`A`$ and $`B`$ are system operators. The time evolution operators can be joined with the “imaginary-time” propagator $`\mathrm{e}^{\beta H}`$ to form a propagator with complex time argument which describes the time evolution of a quantum system along the contour $`𝒞`$ depicted in Fig. 1. This time evolution can be expressed by a path integral over the configuration space of the system

$$AB(t)𝒟[𝐪]\mathrm{e}^{\frac{\mathrm{i}}{\mathrm{}}S_0[𝐪]}I[𝐪]𝐪_{0^+}|A|𝐪_0^{}𝐪_{t^+}|B|𝐪_t^{},$$
(3)
where a normalization factor has been suppressed for simplicity. The action functional

$$S_0[𝐪]=_𝒞d\tau L(𝐪,\dot{𝐪})$$
(4)
depends in the usual way on the classical Lagrange function. The measure $`\mathrm{d}\tau `$ is the complex time differential on $`𝒞`$, and $`\dot{𝐪}=\frac{\mathrm{d}𝐪}{\mathrm{d}\tau }`$. In the generalized influence functional

$$I[𝐪]T_𝒞\mathrm{exp}\left(_𝒞d\tau \frac{\mathrm{i}}{\mathrm{}}H_I(𝐪)\right)_\beta $$
(5)
the thermal average of the contour-ordered exponential is taken with respect to the unperturbed reservoir. This exact result can be simplified  by performing a cumulant expansion of the expectation value in $`I[𝐪]`$. The most commonly discussed case is that of linear dissipation, equivalent to truncating the cumulant expansion after the second order. This procedure is exact for a model reservoir of harmonic oscillators, but also applies to any microscopic model whose collective response to a perturbation $`H_I`$ is dominated by its linear term.

In the following we restrict ourselves to linear dissipation, for which $`I[𝐪]`$ is a Gaussian functional determined by the correlation function $`H_I(𝐪(\tau ))H_I(𝐪(\tau ^{}))`$. After an expansion of $`H_I(𝐪)`$ in terms of an arbitrary set of reservoir operators $`X_j`$,

$$H_I=\underset{j}{}f_j(𝐪)X_j,$$
(6)
the correlation matrix $`L_{jk}(\tau \tau ^{})=X_j(\tau )X_k(\tau ^{})_\beta `$, characterizing the isolated reservoir, completely determines the influence functional

$`I[𝐪]`$ $`=`$ $`\mathrm{exp}({\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle \underset{j,k}{}}{\displaystyle _𝒞}\mathrm{d}\tau {\displaystyle _{\tau ^{}\tau }}\mathrm{d}\tau ^{}f_j(𝐪(\tau ))`$ (8)
$`L_{jk}(\tau \tau ^{})f_k(𝐪(\tau ^{}))).`$
Here “$``$” denotes the order relation induced by the contour orientation.

The aim of the present work is to transform the exact description of open quantum systems by the influence functional technique into an equivalent stochastic propagation of the quantum states of the open system. We first show that the influence functional $`I[𝐪]`$ can be constructed as a noise average of the form

$$I[𝐪]=\mathrm{exp}\left(\frac{\mathrm{i}}{\mathrm{}}\underset{j}{}_𝒞d\tau z_j(\tau )f_j(𝐪(\tau ))\right)_W,$$
(9)
where the average is taken with a Gaussian probability measure $`W[z_j]`$. If the stochastic covariance matrix of the complex-valued noise coordinates $`z_j(\tau )`$ matches the quantum correlation matrix of the coupling operators $`X_j`$,

$$z_j(\tau )z_k(\tau ^{})_W=\{\begin{array}{cc}\frac{1}{2}L_{jk}(\tau \tau ^{}),\hfill & \tau \tau ^{}\hfill \\ \frac{1}{2}L_{jk}(\tau ^{}\tau ),\hfill & \tau \tau ^{},\hfill \end{array}$$
(10)
Eqs. (8) and (9) constitute a formal identity between Gaussian functionals. Noting that Eq. (10) does not fully determine the noise statistics, we find that correlations of the type $`z_j(\tau )z_k^{}(\tau ^{})`$ can be chosen such that the covariance matrix of all real-valued noise components represents a positive quadratic form, i.e., Gaussian noise $`\{z_j\}`$ with the desired properties exists.

For each sample of the noise, the propagation is now governed by a time-local action functional, i.e., the path-integral dynamics can be translated into an equivalent Schrödinger equation. The propagation along the contour $`𝒞`$ is governed by the Hamiltonian $`H_0`$ and a stochastic potential term,

$$\mathrm{i}\mathrm{}|\dot{\psi }=H_0|\psi \underset{j}{}z_j(\tau )f_j(𝐪)|\psi .$$
(11)
This remarkably simple equation shows that linear dissipation can be described exactly by a linear QSD theory containing no memory terms.

In order to compare this new finding with previous results, we need to choose the Feynman-Vernon influence functional, for which the noise vanishes on the imaginary axis. For simplicity, we discuss a one-dimensional coupling $`H_I=qx`$, where $`x`$ is a reservoir operator with correlation function $`L(\tau \tau ^{})`$, associated with a single noise coordinate $`z(\tau )`$. The propagation along the two real-time segments of the contour can be re-stated in the form of a stochastic Liouville equation for the reduced density matrix,

$$\mathrm{i}\mathrm{}\frac{\mathrm{d}\rho }{\mathrm{d}t}=[H_0,\rho ]z_1q\rho +z_2^{}\rho q,$$
(12)
where $`z_1(t)=z(t)`$ and $`z_2(t)=z(t\mathrm{i}\mathrm{}\beta )`$, and where the sample $`\rho `$ of the reduced density matrix is separable, $`\rho =|\psi _1\psi _2|`$. Hence Eq. (12) is just a compact notation for two stochastic Schrödinger equations for $`|\psi _1`$ and $`|\psi _2`$.

The noise forces may be represented as the sum of statistically independent terms, $`z_{1,2}=z+v_{1,2}`$. The terms $`v_1`$ and $`v_2`$ have identical statistics, but are uncorrelated. Performing an average over $`v_{1,2}`$ takes us back to the partial decomposition of the influence functional given by Diósi, Strunz and Gisin . Using the same representation, the stochastic Schrödinger equation of Markovian QSD can be recovered. The same limits and approximations used to derive Lindblad-type Master equations and their QSD counterparts from a system-reservoir model, i.e., weak damping and time coarse-graining, can be applied to Eq. (12). Coarse-graining is done at an intermediate timescale which is short compared to the relaxation time of the damped system but long compared to its oscillation periods and the relaxation times of the reservoir. In this approximation the dynamics is only sensitive to segments of the noise spectrum of $`z`$ which are centered around the natural transition frequencies of the undamped system. In the Born-Markov limit the coarse-graining time scales to zero, and each of these segments naturally transforms into a distinct Markovian noise force coupled to a Lindblad operator associated with the respective quantum transition. Finally, the work of Stockburger and Mak  using real noise and earlier work concentrating on strictly ohmic damping can be recovered from the present decomposition of the Feynman-Vernon influence functional by applying the stochastic decomposition only to the real part of $`L(\tau \tau ^{})`$ while explicitly evaluating memory effects induced by the imaginary part.

As in the linear versions of conventional QSD, $`tr\rho `$ is not conserved by the time evolution of Eqs. (11) and (12). A transparent physical interpretation of the noise forces and stochastic samples arises in a modified nonlinear theory for which the trace of the reduced density matrix is preserved not only in the stochastic average, but for each sample.

It will be advantageous to introduce new noise variables $`\xi =(z_1+z_2^{})/2`$ and $`\mathrm{}\nu =z_1z_2^{}`$, with resulting covariances

$`\xi (t)\xi (t^{})_W`$ $`=`$ $`\mathrm{Re}L(tt^{})`$ (13)
$`\xi (t)\nu (t^{})_W`$ $`=`$ $`(2\mathrm{i}/\mathrm{})\mathrm{\Theta }(tt^{})ImL(tt^{})`$ (15)
$`=\mathrm{i}\chi _R(tt^{})`$
$`\nu (t)\nu (t^{})_W`$ $`=`$ $`0,`$ (16)
where $`\chi _R`$ is the response function of the reservoir. The noise force $`\xi (t)`$ can be chosen real, and this case will be discussed here for simplicity. For comparison with phenomenological friction models, one conventionally augments the coupling Hamiltonian by a counter-term in order to make the coupling translationally invariant. The additional term is a quadratic potential modification $`\mu q^2/2`$ with $`\mu =_0^{\mathrm{}}dt\chi _R(t)`$, which eliminates the static response of the reservoir. Eq. (12) is thus transformed into

$$\mathrm{i}\mathrm{}\dot{\rho }=[H_0,\rho ]+\frac{\mu }{2}[q^2,\rho ]\xi [q,\rho ]\frac{\mathrm{}}{2}\nu \{q,\rho \}.$$
(17)
For the normalized density matrix sample $`\widehat{\rho }=\rho /tr\rho `$, this translates into the quasilinear equation of motion

$$\mathrm{i}\mathrm{}\dot{\widehat{\rho }}=[H_0,\widehat{\rho }]+\frac{\mu }{2}[q^2,\widehat{\rho }]\xi [q,\widehat{\rho }]\frac{\mathrm{}}{2}\nu \{q\overline{q},\widehat{\rho }\}$$
(18)
with $`\overline{q}=trq\widehat{\rho }`$. When averaging the normalized density matrix $`\widehat{\rho }`$, the factor

$$tr\rho =\mathrm{exp}\left(\mathrm{i}_0^tdt^{}\overline{q}(t^{})\nu (t^{})\right)$$
(19)
needs to be incorporated into the integration measure. Because the original integration measure $`W[\xi ,\nu ,\nu ^{}]`$ is Gaussian, the new measure $`W_t[\xi ,\nu ,\nu ^{}]`$ including the factor $`tr\rho `$ can be rewritten by “completing the square” in the exponent of the Gaussian functional, making the exponent a quadratic form of shifted variables, and yielding $`W_t[\xi ,\nu ,\nu ^{}]=W[\xi _t,\nu _t,\nu _t^{}]`$ with

$`\xi _t(t^{})=\xi (t^{})+{\displaystyle _0^t}ds\chi _R(t^{}s)\overline{q}(s)`$ (20)
$`\nu _t(t^{})=\nu (t^{})`$ (21)
$`\nu _t^{}(t^{})=\nu ^{}(t^{})\mathrm{i}{\displaystyle _0^t}ds\nu ^{}(t^{})\nu (s)_W\overline{q}(s)`$ (22)
Using $`W[\xi _t,\nu _t,\nu _t^{}]`$ as integration measure allows us to interpret the shifted variables as noise to be used for stochastic integration. The Jacobian determinant of the variable change is unity: Because $`\nu _t=\nu `$, and because $`\overline{q}`$ is independent of $`\nu ^{}`$, it simplifies to

$$\left|\begin{array}{ccc}\frac{\delta \xi _t}{\delta \xi }& \frac{\delta \xi _t}{\delta \nu }& \frac{\delta \xi _t}{\delta \nu ^{}}\\ \frac{\delta \nu _t}{\delta \xi }& \frac{\delta \nu _t}{\delta \nu }& \frac{\delta \nu _t}{\delta \nu ^{}}\\ \frac{\delta \nu _t^{}}{\delta \xi }& \frac{\delta \nu _t^{}}{\delta \nu }& \frac{\delta \nu _t^{}}{\delta \nu ^{}}\end{array}\right|=\left|\begin{array}{ccc}\frac{\delta \xi _t}{\delta \xi }& \frac{\delta \xi _t}{\delta \nu }& 0\\ 0& 1& 0\\ \frac{\delta \nu _t^{}}{\delta \xi }& \frac{\delta \nu _t^{}}{\delta \nu }& 1\end{array}\right|=\left|\frac{\delta \xi _t}{\delta \xi }\right|.$$
(23)
Using the fact that both the bath response and $`\overline{q}[\xi ,\nu ]`$ are causal, one finds that $`\left|\delta \xi _t/\delta \xi \right|`$ is represented by a triangular matrix with unit diagonal elements, i. e., $`\left|\delta \xi _t/\delta \xi \right|=1`$.

After the variable change (with subscripts $`t`$ now dropped) the equation of motion is rewritten as

$`\mathrm{i}\mathrm{}\dot{\widehat{\rho }}`$ $`=`$ $`[H_0,\widehat{\rho }]\xi [q,\widehat{\rho }]+m\gamma (t)\overline{q}(0)[q,\widehat{\rho }]`$ (26)
$`+m{\displaystyle _0^t}dt^{}\gamma (tt^{})\dot{\overline{q}}(t^{})[q,\widehat{\rho }]`$
$`+{\displaystyle \frac{\mu }{2}}[(q\overline{q})^2,\widehat{\rho }]{\displaystyle \frac{\mathrm{}}{2}}\nu \{q\overline{q},\widehat{\rho }\}.`$
Here we have integrated by parts, using the relation $`m\dot{\gamma }(t)=\chi _R(t)`$ between the reservoir response function and the friction kernel $`\gamma (t)`$. Remarkably, the mere choice of a constraint for $`tr\widehat{\rho }`$ has taken us from the abstract mathematical decomposition (9) to an equation of motion with a number of terms looking quite familiar—as in classical dissipation, the dynamics is governed by c-number force terms representing noise and a mean-field version of classical friction. Additional terms ascertain that the noise-averaged $`\widehat{\rho }`$ contains the effect of all system-reservoir correlations.

Further examining this analogy, we find a “sample-by-sample” correspondence between quantum mechanical and classical stochastic dynamics. The time evolution of Eq. (26) can be described in a reduced Heisenberg picture, which is defined by interpreting expectation values of the form $`tr\{A^{}\widehat{\rho }\}`$ as a scalar product $`(A,\widehat{\rho })`$ of two operators. This allows the introduction of Heisenberg operators $`A(t)`$ through the propagating superoperator $`𝒰(t)`$ and its adjoint,

$$(A,\widehat{\rho }(t))=(A,𝒰(t)\widehat{\rho }_0)=(𝒰^{}(t)A,\widehat{\rho }_0)=(A(t),\rho _0).$$
(27)
In the classical limit, terms proportional to $`q\overline{q}`$ in Eq. (26) can be neglected, and the equations of motion governing the evolution of Heisenberg operators $`A(t)`$ become those of a classical stochastic system. For a potential model, these reduce to the generalized Langevin equation

$$m\ddot{q}(t)=V^{}(q(t))m_0^tdt^{}\gamma (tt^{})\dot{q}(t^{})+\xi (t).$$
(28)
The limit $`\mathrm{}0`$ turns our stochastic ensemble of quantum states gradually into an ensemble of phase-space points with classical thermal noise and memory friction, i.e., the quantum dynamics defined in Eq. (26) obeys a stochastic correspondence principle.

Fig. 2 shows numerical results for the symmetric spin-boson system, where $`q`$ is the Pauli matrix $`\sigma _z`$ and $`H_0=(\mathrm{}\mathrm{\Delta }/2)\sigma _x`$, compared to a QMC algorithm which directly samples a path integral . For the coherent dynamics shown in Fig. 2a, the QMC method does not fully converge even if run longer than our simulation. Due to an only marginally ergodic Metropolis random walk, the accuracy of the QMC method degrades significantly with increasing $`t`$. Problems of this kind never appear in our algorithm because our samples are statistically independent by construction. The incoherent dynamics shown in Fig. 2b, a forte of the QMC method, is treated equally well by both approaches.

We have derived exact stochastic equations of motion for the reduced density matrix in linear form (11) and trace-conserving form (26). In the latter, a close correspondence with classical Langevin equations is evident. Examples of numerical results show that the stochastic equations of motion presented here provide a promising alternative to currently available exact methods based directly on the path integral representation. Given the freedom of choice for the noise covariances left undetermined by the physics of quantum dissipation, future improved numerical methods are to be anticipated.

The authors would like to thank Walter Strunz for valuable discussions and correspondence. This work was supported by Deutsche Forschungsgemeinschaft under SPP 470.

To be published in Phys. Rev. Lett. (2002)


text:
# Stellar Contents and Globular Cluster Candidates in the Sculptor Group Galaxy NGC 300

## 1 INTRODUCTION

Sculptor group is the nearest galaxy group to the Local Group containing five major spiral galaxies (NGC 55, NGC 247, NGC 253, NGC 300, and NGC 7793) and $``$20 dwarf galaxies (Côté et al. 1997, Whiting 1999; cf. van den Bergh 1999). Due to the rather small distance to the Sculptor group, there have been many studies on the galaxies in this group.

NGC 300 (=ESO 295-20, IRAS 00523-3756, PGC 3238) is an SA(s)d spiral galaxy in the Sculptor group. Basic information of this galaxy is summarized in Table 1. It is a rather bright ($`M_B=18.6`$) and nearly face-on galaxy, similar to M33 seen in the northern hemisphere (Blair & Long 1997). There have been several studies on NGC 300, and some of them are summarized in the following. Pietrzyński et al. (2001) have found 117 OB associations in NGC 300, Davidge (1998) has studied the evolved red stellar contents from $`JHK`$ photometry, Breysacher et al. (1997) have studied Wolf-Rayet stars detected in five associations of NGC 300, and Soffner et al. (1996) have searched for and found 34 planetary nebulae and 88 H ii regions. Freedman et al. (1992) obtained the distance modulus to this galaxy of $`(mM)_0=26.66\pm 0.10`$ mag (d = $`2.1\pm 0.1`$ Mpc) using 16 known Cepheid variables. Pannuti et al. (2000) have conducted supernova remnant (SNR) search in NGC 300 using X-ray, optical, and radio observations, resulting in total 44 SNR candidates including 16 new objects. Recently Read & Pietsch (2001) have studied X-ray emission sources in NGC 300 using $`ROSAT`$ data. Freedman et al. (2001) have revised the distance modulus to NGC 300 using the Hubble Space Telescope Key Project Cepheid data to be $`(mM)_0=26.53\pm 0.07`$ mag (d = $`2.02\pm 0.07`$ Mpc), which is adopted in this study.

There have been a few surveys for globular clusters(GCs) in the Sculptor group galaxies. From excellent seeing plates and spectrophotometry, Da Costa & Graham (1982) have found three GCs in edge-on galaxy NGC 55. Liller & Alcaino (1983a) have found 51 GC candidates in NGC 55 from excellent-seeing plates taken with ESO 3.6m telescope. Liller & Alcaino (1983b) also have found 63 GC candidates in NGC 253 from excellent-seeing plates taken with ESO 3.6m telescope. In NGC 253, Blecha (1986) have found 25 GC candidates from using high-resolution-wide field McMullan electronographic camera attached to the Danish 1.5m telescope. Later, for the 57 and 58 GC candidates in NGC 55 and NGC 253, respectivey, Beasley & Sharples (2000) have performed spectroscopic GC survey and identified one probable GC and 14 GCs in NGC 55 and NGC 253, respectivey. And they have estimated total number of GCs in NGC 55 and NGC 253 to be $``$20 and $``$60, respectivey. Watson et al. (1996) have discovered four young, luminous, compact stellar clusters in the central region of the starburst spiral galaxy NGC 253, which are believed to be excellent candidates of young GCs. Recently Jerjen & Rejkuba (2001) have found one GC candidate in dwarf galaxy ESO 540-032 using Very Large Telescope UT1+FORS1. However, little is known about GCs in NGC 300. Only Bright & Olsen (2000) have indicated that they obtained 100 – 150 GC candidates in each of seven brightest galaxies in the Sculptor group. But the detailed information is not yet published. In this paper we report the first list of GC candidates in NGC 300.

This paper is composed as follows. Section II describes observations and data reduction, and Section III presents photometric diagrams and photometric analyses of 18 OB associations in NGC 300. Section IV presents the search method and the result for new GC candidates, and Section V discusses mass function of bright stars, brightest stars, some characteristics of GCs, and the nucleus of NGC 300. Finally, a summary and conclusions are given in Section VI.

## 2 OBSERVATIONS

$`UBVI`$ CCD photometry was performed on 1997 November 23 at the Siding Spring Observatory with the 40 inch (1m) telescope (f/8) and a thinned SITe 2048 $`\times `$ 2048 CCD (24$`\mu `$m pixels). The scale was 0.′′602 pixel-1, giving 20.5 on a side. Exposure times used in the observations were 60 s and 120 s in $`I`$, 300 s in $`V`$, 600 s in $`B`$, 1200 s in $`U`$. The night was photometric and the seeing was 1.′′4 in the 300 sec $`V`$ image. A finding chart for the observed area is shown in Fig. 1. In the figure we marked only stars with $`V`$ = 16.5 – 21 mag to show the spatial distribution of bright stars in NGC 300. In addition we also marked 18 rich OB associations and 17 GC candidates.

All the preprocessing, such as overscan correction, bias subtraction and flat fielding were done using the IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc./CCDRED package. Instrumental magnitudes were obtained using the IRAF version of DAOPHOT II (Stetson 1990) via point spread function (PSF) fitting. All the instrumental magnitudes were transformed to the standard Johnson-Cousins $`UBVI`$ system using SAAO observations of equatorial standards (Menzies et al. 1991) and blue and red standards by Kilkenny et al. (1998). The primary extinction coefficients were determined and used. Details involving standard transformations can be found in Sung & Bessell (2000).

The total number of stars which were measured down to $`V`$ = 22 mag is 6027. Table 2 lists the photometry of 167 stars with $`V<19.0`$ mag and contains the running number, $`\alpha `$ (J2000.0), $`\delta `$ (J2000.0), the weighted mean values of the magnitude and color, the weighted errors, and the number of independent measurements. The variable stars found by Graham (1984), the photometric sequence of Graham (1981) and Walker (1995), and GC candidates newly found in this paper are identified in the last two columns.

The mean values of the errors are listed in Table 3 (“N” in Table 3 represents the number of stars used in the statistics). The $`UBVI`$ magnitude and colors for $`V_{\mathrm{CCD}}20`$ mag are compared with previous photoelectric or CCD photometry by Graham (1981) and Walker (1995) in Fig. 2, showing differences in the sense, previous photometry minus our new photometry. Our photometry is consistent with that by Graham (1981) or by Walker (1995) in $`BV`$ and $`VI`$ ($`\mathrm{\Delta }(BV)=+0.008\pm 0.013`$ and $`\mathrm{\Delta }(VI)=0.007\pm 0.014`$), but slightly fainter in $`V`$ ($`\mathrm{\Delta }V=0.025\pm 0.024`$) and redder in $`UB`$ ($`\mathrm{\Delta }(UB)=0.042\pm 0.016`$).

## 3 PHOTOMETRIC DIAGRAMS AND OB ASSOCIATIONS

### 3.1 Photometric Diagrams

As NGC 300 is a member of Sculptor group of galaxies ($`d`$ 2.02 Mpc : Freedman et al. 2001), the stars in the galaxy are very faint. In addition many bright blue stars are members of OB associations, and therefore their photometric values contain inevitably large errors due to crowding effect even though we paid much attention to the photometry. The resulting color-magnitude diagrams (CMDs) from our photometry are presented in Fig. 3. In the figure we marked the GC candidates (open circle with dot; see Section IV for the selection criteria) and variable stars found by Graham (1984) (open square with dot). Among 34 variable stars, 22 were identified from the coordinates and charts given by Graham (1984) (The coordiantes given by him are not accurate enough to identify the variables unambiguously). Among the remaining 12 stars, 9 stars are outside our field of view, 3 stars are too faint to be measured from our images.

The ($`V,VI`$) diagram (Fig. 3 (a)) do not give much information on the photometric characteristics of the stars in NGC 300 due partly to the short exposures of $`I`$ images and partly to the characteristic nature of $`VI`$ color, which is insensitive to the temperature variation of blue stars. On the other hand, $`UB`$ is rather sensitive to the temperature variation of hot stars, but the blue stars in NGC 300 are very faint in $`U`$ pass band. Therefore the ($`V,UB`$) CMD also does not give much information on the faint stars. Among the three CMDs, the ($`V,BV`$) diagram shows well-defined blue stars more clearly. We will use only the ($`V,BV`$) diagram in the analysis of the photometric characteristics of stars in NGC 300. In this figure, a large number of faint blue stars can easily be found. They are bona fide members of NGC 300. Apparently stars brighter than $`V`$ = 17.5 mag seem to be foreground thick disk or halo stars in our Galaxy.

We presented the ($`UB,BV`$) diagram in Fig. 4. As the most bright blue stars are suspected to be blue supergiants in NGC 300, we drew the intrinsic relation for supergiants (Sung 1995) in the figure. Even though the photometric error in $`BV`$ as well as $`UB`$ is large, most blue stars are scattered around the intrinsic relation for supergiants, indicating that the internal reddening is very small. There are several yellow stars ($`BV`$ = 0.4 – 0.8 mag). A few of them are GC candidates, but the remainings are foreground stars in the Galaxy. As they are suspected as thick disk or halo stars in the Galaxy, they show evident UV excesses. The others with abnormal red colors in the diagram are faint red stars. Many of them ($`V<`$ 17.5 mag) are foreground red stars in the galaxy, but a few faint red stars may be red supergiant stars in NGC 300. Their abnormal colors are mainly due to large photometric errors. As NGC 300 is located near the south galactic pole ($`b=79.\mathrm{°}42`$), the foreground reddening due to the interstellar cloud in the Galaxy is practically zero. In addition NGC 300 is nearly face-on, and therefore the internal reddening may not be very large. For the reddening correction of stars in NGC 300 we have used the value $`E(BV)`$ = 0.013 mag obtained by Schlegel, Finkbeiner, & Davis (1998).

### 3.2 OB Associations

Recently Pietrzyński et al. (2001) published a catalog of OB associations in NGC 300. They found 117 OB associations using Path Linkage Criterion technique. The telescope they used is the ESO/MPI 2.2 m, and therefore the limiting magnitude of their photometry is about 1 mag fainter than ours. To derive the age distribution of OB associations in NGC 300, the OB associations having more than 10 stars in our photometry were selected. In the analysis, we limited our attention to the OB associations having more than 30 member stars in the catalog ($`N_{}>30`$). Among 117 OB associations, 18 OB associations are satisfied the above criterion. We present 18 CMDs with appropriate isochrones of the Padova group (Bertelli et al. 1994) in Fig. 5. Due to the large photometric errors of faint stars in NGC 300, the CMD of individual OB association shows a large scatter. The color of faint stars are more affected by crowding effect (see Sung et al. (1999) for the effect of crowding on color). The large scatter in $`BV`$ is therefore due partly to the crowding effect and partly to the contamination of faint field stars in NGC 300. In addition we cannot exclude the possibility of the contamination of foreground stars in the Galaxy or GCs in NGC 300 among yellow or red stars in the line of sight. We determined the age of each association by paying more attention to the bright blue stars in the association.

To obtain the reliable ages of OB associations or young open clusters, spectral type and luminosity class should be known because the photometric indices such as $`BV`$ or even $`UB`$ are inadequate for the accurate estimate of temperature of O- and early B-type stars. Actually no information on the spectral type of blue stars in NGC 300 is available currently. The ages of OB associations in NGC 300 have, therefore, inevitably large errors if we determine the age from photometric colors alone. We should determine the age of each OB association with such a limitation. Among the 18 OB associations, 13 associations have ages younger than 10 Myr ($`\mathrm{log}\tau <7`$). Many of them may have much younger ages (only about a few Myr) because there are a few very bright blue stars ($`M_V<8`$ mag). The associations that belong to this category are AS\_016, \_029, \_052, \_082, \_084, \_086, and \_102. For these associations the ages from the CMDs are highly uncertain. There are five associations (AS\_018, \_026, \_027, \_058, and \_096) with ages older than 10 Myr. They have apparently no O-type stars and seem to be old B associations. A few OB associations seem to contain a red supergiant candidate (AS\_002, \_009, and less luminous red supergiant candidates in AS\_101).

Spatially these rich OB associations ($`N_{}>30`$) are more concentrated in the NW part of the galaxy, but the brightest OB association in NGC 300 is AS\_102, about 4 E of the center of NGC 300. Two out of three brightest blue stars in NGC 300 are in this association. In the H$`\alpha `$ image of NGC 300, AS\_102 is the brightest in H$`\alpha `$, which indicates that AS\_102 is associated with the brightest H ii region. AS\_102 is the most active star forming region in NGC 300. There are many bright blue stars in the CMD of AS\_102. Even though the photometric error is large, there is few red supergiant candiates in the association. This fact again implies that the association is extremely young. The second and third brightest OB associations are AS\_052 and \_084. Even though AS\_052 is the richest OB association ($`N_{}=101`$) in NGC 300, the brightest star in the associaton is about 0.8 mag fainter that the brightest star in AS\_102.

## 4 GLOBULAR CLUSTER CANDIDATES

### 4.1 Search Method

Globular clusters in NGC 300 are barely resolved in our images. We have used photometric information, morphological parameters and visual inspection to survey globular cluster candidates as follows. From the ($`V,BV`$) CMD in Fig. 3, we have first taken the point sources with $`V<20.0`$ mag and $`0.3<(BV)<2.0`$ to search for globular cluster candidates. We have considered $`V=20.0`$ mag as the reliable magnitude limit and $`0.3<(BV)<2.0`$ to be the safe color range that contains all the GC candidates in NGC 300. For these selected point sources we have calculated the morphological image classification parameters ($`r_1`$ and $`r_2`$ radial moments) in three ($`V`$-band and two $`I`$-band) best seeing images. The two radial moments are effective image radii in pixel units, where $`r_1`$ is an indicator of image wing spread and $`r_2`$ indicates the image central concentration (see Harris et al. 1991; Kim et al. 2002; Lee et al. 2002). Because the images have degradations at the corners, we only considered the objects inside of a circle centered at the image center and with the diameter of the image width ($`20.^{}5`$).

Fig. 6 shows the $`r_1`$ and $`r_2`$ parameters calculated on the $`V`$-band image. The horizontal sequence of objects are mainly stars, objects with a little larger moment values should be GCs in NGC 300, and some objects with largest moment values are probably background galaxies. The box region contains the objects for which we have performed visual examination on each image and investigated radial, contour, and surface profiles. The circled objects denote the resulting GC candidates and the size of the circle is proportional to the probability to be a GC; the largest circles denote probable GC candidates (class 1), the intermediate size circles denote the possible candidates (class 2), and the smallest circles denote doubtful candidates (class 3).

Fig. 7 shows the $`V`$-band mosaic images and radial profiles of a typical star and a class 1 GC candidate discovered in this study. The FWHM of the star in Fig. 7 is $`2.4`$ pixel, while that of the GC candidate is $``$ 2.8 pixel. Though the GC candidate is fainter than the star ($`V`$=18.2 versus $`V`$=17.4), it is clear that the GC candidate is more extended than the star.

We also have performed smoothing the images using the IRAF/RMEDIAN package and subtracted the results from the originals to form the final images on which the GC search was made again. While these images clearly show the fine structures of NGC 300 especially in the nuclear region and patchy regions, they are not much better in the search for GCs than the original images.

There are several sources that we could not be able to investigate the image structure due to the crowding of two or more objects in small region and large pixel size of our images (1 px $`=`$ $`0.^{\prime \prime }602`$ = 5.89 pc assuming d = 2.02 Mpc, and so $`1^{\prime \prime }`$ $`=`$ 9.79 pc). Imaing data of better quality (e.g. Bright & Olsen 2000) will be helpful to uncover the nature of these objects.

### 4.2 Search Result

In Table 4 we present the newly discovered GC candidates in NGC 300. Column 1 gives the GC candidate ID numbers. Columns 2 and 3 are the RA and DEC coordinates (epoch J2000), columns 4-7 are $`V`$ magnitude, and $`VI`$, $`BV`$, and $`UB`$ colors, respectively, columns 8-11 are the errors in $`V`$, $`VI`$, $`BV`$ and $`UB`$, respectively, columns 12-15 are the number of independent measurements and the last column is the class of the candidates.

The positions of these GC candidates in the CMD and ($`UB,BV`$) diagram are shown in Fig. 3 and Fig. 4, respectively. The $`V`$ magnitudes of most of the candidates are between $``$18 to 20 mag. The $`BV`$ colors of the candidates are between 0.3 and 1.1 and the $`VI`$ are between 0.6 and 1.2. Some characteristics of these GC candidates will be discussed in the next Section.

## 5 DISCUSSION

### 5.1 Stellar Mass Function

The upper part of stellar mass function (MF) of NGC 300 can be derived from the photometry of bright stars in the galaxy. But as mentioned in Section III (b) the most massive part of MF is highly uncertain because we have only the photometric data for the stars. Due to the degeneracy of $`BV`$ color for the hottest stars (T$`{}_{eff}{}^{}>`$20,000K), it is nearly impossible to estimate accurately the mass of blue stars in NGC 300. To show this effect more clearly, we plotted two CMDs with theoretical evolution tracks in Fig. 8. The evolution tracks in Fig. 8 (a) are those of the Padova group (Bertelli et al. 1994), and are transformed to the observational plane using the calibration of Sung (1995). The model shows that most massive star (m = 100$`M_{}`$) does not evolve to the brightest part of the CMD. On the other hand, that of the Geneva group (Schaller et al. 1992) shows that 120 $`M_{}`$ star evolves to the brightest part of CMD. The $`\eta `$ Car, the most massive star in the Galaxy (m $`100M_{}`$) is one of the brightest stars in the Galaxy. This means that at least a few of the bright blue stars in NGC 300 have masses around $`100M_{}`$.

To derive the MF of NGC 300, we should know the mass - luminosity (ML) relation. To derive the ML relation, we used the evolution tracks by the Geneva group (Schaller et al. 1992) and the Padova group (Bertelli et al. 1994), and the calibration of Sung (1995). We simply divide the stars in the CMD into 4 groups (4 color bins - I: $`BV`$ 0.2, II: 0.2 $`<BV`$ 0.8, III: 0.8 $`<BV`$ 1.3, and IV: $`BV>`$ 1.3). And we calculated the mean $`M_V`$ for a given model star in each color bin, i.e. $`<M_V>M_V(j)𝑑\tau (j)/𝑑\tau (j)`$ during the star evolves in the given color bin. The ML relation derived from the above scheme is strongly dependent on the absolute magnitude which the star spends most of time for the given color bin, i.e. for the color bin II, III, and IV the position of helium burning phase. The theoretical ML relation for given color bin is presented in Fig. 9 (a). For the blue stars (group I), they are actually evolved blue stars, but the theoretical ML relation derived using the above scheme is that for main sequence (MS) stars, not for evolved blue (super)giants because the star of a given mass spends most of their life time in the MS stage. Therefore the theoretical ML relation for group I is inadequate for the actual case. For group I stars, we simply determined the mid-point of evolution tracks. The adopted ML relation for the group I is also presented in Fig. 9 (a).

To derive the luminosity function (LF) of each group, we calculated the number of stars for each group. And we subtracted the number of GC candidates and the number of foreground stars based on the estimates by Ratnatunga & Bahcall (1985) for the Sculptor dwarf galaxy. Ratnatunga & Bahcall (1985) presented the expected number of foreground stars in the direction of the Sculptor dwarf galaxy in $`\mathrm{\Delta }V`$ = 2 mag and three bins in $`BV`$ ($`<`$ 0.8, 0.8 – 1.3, and $`>`$ 1.3 mag). We estimated the number of stars for given color bin by keeping the trend of the variation of star number. For blue stars (group I) most of them are bona fide members of NGC 300. Even though there are apparently many bright stars in the group II (see Fig. 8), the final number of stars in group II is very few. The resulting LF of stars in NGC 300 is presented in Fig. 9 (b). In the transformation of the LF to the MF, we used the ML relations derived above using the models by the Geneva group, i.e. $`\xi (\mathrm{log}m)d\mathrm{log}mN(M_V)dM_V`$ (Salpeter 1955). The resulting MF of each group is presented in Fig. 9 (c).

To derive the initial mass function (IMF) of a heterogeneous group of stars, such as a galaxy like NGC 300, we should know the star formation history. It is practically impossible to know the history of star formation in NGC 300. To derive the IMF of NGC 300, we assumed that the star formation rate is constant with time. This assumption is valid at least for the massive stars. For this case, the transformation of MF to IMF is very simple, i.e. only applying the difference in the total life span between stars with different masses. By interpolating the MF of each group, we derived the MF of bright stars in NGC 300 by summing up the MF of each group. After that by dividing the difference in life span, we derived the IMF of NGC 300. The total MF and IMF is presented in the same figure. The slope of IMF ($`\mathrm{\Gamma }d\mathrm{log}\xi (\mathrm{log}m)/d\mathrm{log}m`$) of NGC 300 is $`2.6`$ $`\pm `$ 0.3. This value is steeper than that of the solar neighborhood ($`\mathrm{\Gamma }=1.35`$: Salpeter 1955), but less steeper than the general field of the magellanic clouds ($`\mathrm{\Gamma }`$ = $`4.1`$ $`\pm `$ 0.2 for LMC and $`3.7`$ $`\pm `$ 0.5 for SMC: Massey et al. 1995). The IMF for NGC 300 is that for the bright stars measured in this study, i.e. a mixed group of stars, and the slope of the IMF is reasonably agreed to the result obtained by Massey et al. (1995). In addition, if the star formation rate is a decreasing function of time, then the resulting slope of IMF is steeper than the actual value of the slope.

### 5.2 The Brightest Stars

The birghtest stars in galaxies have been used as distance indicator (Lyo & Lee 1997). In order to avoid the contamination of foreground stars, the mean brightness of three brightest stars in a galaxy is used for this purpose. For NGC 300, the identification of red supergiant stars is very difficult due to the possible contamination of foreground thick disk or halo stars in the Galaxy as well as the GCs in NGC 300. For blue stars in NGC 300, they are well-separated from the foreground stars, and therefore they are nearly free from the contamination.

Three brightest blue stars in NGC 300 are easily selected: \# 50, 51, and 52 with similar magnitudes. The mean magnitude and color of these stars are derived to be $`<V>`$ = 17.62 $`\pm `$ 0.01 mag, $`<BV>`$ = $`0.21`$ $`\pm `$ 0.05 mag. Assuming the distance to NGC 300 adopted in this study, the mean absolute magnitude of three brightest blue stars is derived to be $`<M_V^{BSG}(3)>`$ = $`8.95`$ mag. The brightness of these blue supergiants in NGC 300 is very similar to that of three blue supergiants in M33 (Lyo & Lee 1997). The calibration for the blue supergiants was given by Lyo & Lee (1997): $`<M_V^{BSG}(3)>=0.32(\pm 0.05)M_B^T2.69(\pm 0.94)`$. Using this calibration, the absolute magnitude of the brightest blue supergiants expected for the luminosity of NGC 300 ($`M_B^T=18.59`$) is derived to be $`<M_V^{BSG}(3)>=8.63(\pm 0.54)`$, which agrees within calibration uncertainty with the value obtained in this study.

### 5.3 Characteristics of the Globular Cluster Candidates

#### 5.3.1 Luminosity function

Fig. 10 shows the GC luminosity function (GCLF) of NGC 300 obtained for the 17 GC candidates found in Section IV relative to the GCLF of M31 GCs. Only class 1 objects are plotted in Fig. 10 (a), objects of classes of 1 and 2 are plotted in Fig. 10 (b), and objects of all three classes are plotted in Fig. 10 (c). Milky Way Galaxy (MWG) GC data (N=146) taken from Harris (1996) are shown in dotted lines and M31 GC data (N=429) taken from Barmby et al. (2000) are also shown in dashed lines for comparison.

GCLF is well fitted by Gaussin or $`t_5`$ function (Barmby, Huchra, & Brodie 2001; Secker 1992; Secker & Harris 1993). However our GCLF of NGC 300 shown in Fig. 10 (c) is far from Gaussin shape and clearly lacks symmetry. Considering the fact that the faintest bin of our GCLF is much affected by the photometric errors and assuming the turnover of NGC 300 GCLF is also at $`M_V7.4`$ as in the case of the Milky Way (Harris 2001), it can be said that we have obtained the bright half of the LF.

From the largest number GCLF of NGC 300 (Fig. 10 (c)) it is concluded that the faintest bin of GCLF of NGC 300 may possibly contain some non-GCs, and there maybe more fainter (than the turnover) GCs not discovered yet in NGC 300. Further search for the rest of GCs is needed for this galaxy.

#### 5.3.2 Color distribution

We present the $`BV`$ and $`VI`$ color distribution of our GC candidates in Fig. 11. The data of MWG GCs (Harris 1996; N=117 for $`BV`$ and N=97 for $`VI`$) and M31 GCs (Barmby et al. 2000; N=383 for $`BV`$ and N=242 for $`VI`$) are also plotted with dotted and dashed lines, respectively. The labels in the ordinates are for the GC candidates in NGC 300 and are in arbitrary scales for the GCs in M31 and MWG.

In general the color distributions of the three galaxies agree well, though there are some differences: (1) the peak is skewed to the blue side ($`(BV)_{peak}0.5`$) in the $`(BV)`$ color distribution (Fig. 11 (c)) and there are double peaks ($`(VI)_{peak}0.6`$ and 1.2) in the $`(VI)`$ color distribution (Fig. 11 (f)), and (2) there are more blue GC candidates in NGC 300 than red GC candidates compared to MWG or M31. The reason for the latter might be either ($`i`$) NGC 300 has actually larger rate of blue to red clusters than MWG or M31 or ($`ii`$) we have discarded more red GC candidates than blue ones taking them to be background galaxies or foreground stars during the search process. It is not easy to find out which is the cause at the present. It is needed to perform a new, deep search for GCs in NGC 300 to find more clusters and to distinguish which is the cause.

#### 5.3.3 Spatial distribution

In Fig. 12 we show the histogram of the projected galactocentric distances of the GC candidates from the nuclues of NGC 300 in arcmin. Most candidates distribute evenly in the range $`1^{}<R_{gc}<11^{}`$, while there are more objects at $`4^{}<R_{gc}<8^{}`$. It is not clear if this underabundance of GC candidates in the galaxy center is real or due to the low discovery rate of our search in the galactic center region.

### 5.4 Nucleus

The nucleus of NGC 300 was originally included in our GC candidates as Class 1 or 2 object, which means it could be regarded as a star cluster in photometric images. Based on little emission from ionized gas from the nuclear region of NGC 300, Soffner et al. (1996) have suggested that the nucleus of NGC 300 is very probably an unresolved compact stellar cluster, similar to the nucleus of M33 (Kormendy & McClure 1993). A detailed aperture photometry and surface photometry of the nucleus region and the surface photometry of the whole CCD area of NGC 300 (just like in the case of NGC 205: Kim & Lee 1998) will follow in a subsequent paper(Kim et al. in preparation).

## 6 SUMMARY AND CONCLUSIONS

We have presented the $`UBVI`$ CCD photometry of $`20.^{}5\times 20.^{}5`$ area of the Sculptor group galaxy NGC 300. The primary results from this study of the stellar populations and the GC candidates in this galaxy are summarized as follows.

(1) The ($`V,BV`$) CMD shows well-defined blue stars extending up to $`V17.5`$ mag.

(2) The ($`UB,BV`$) diagram shows that most bright blue stars in NGC 300 are blue supergiants.

(3) CMDs of 18 OB associations with more than 30 member stars are presented together with ages from the Padova isochrones.

(4) The initial mass function (IMF) of bright stars in NGC 300 is derived to have the slope of $`\mathrm{\Gamma }=2.6`$ $`\pm `$ 0.3, which is steeper than that of the solar neighborhood ($`\mathrm{\Gamma }=1.35`$), but less steeper than those of the general field of the magellanic clouds ($`\mathrm{\Gamma }`$ = $`4.1`$ for LMC and $`3.7`$ for SMC).

(5) Assuming the distance of NGC 300 to be $`(mM)_0=26.53\pm 0.07`$, the mean absolute magnitude of three brightest blue stars is estimated to be $`<M_V^{BSG}(3)>`$ = $`8.95`$ mag.

(6) From the search for GCs in NGC 300, we discovered 17 GC candidates and they are 4 probable (class 1), 7 possible (class 2), and 6 doubtful (class 3) candidates. Some characteristics of these newly found GC candidates are discussed.

This work was supported in part by the Korea Research Foundation Grant (KRF-2000-DP0450).


text:
# Near-optimal two-mode spin squeezing via feedback

## I Introduction

Two-mode spin squeezed (TMSS) states are important to consider for a number of reasons. First, they may be used to demonstrate entanglement experimentally using only spin measurements Julsgaard . Second, as they are entangled, they may be used as a resource for quantum information, for example, quantum teleportation Kuzmich ; Berry , quantum computing, superdense coding, etc. NC . In addition, TMSS states are analogous to two-mode squeezed states for light Caves85 , which have proven to be useful for quantum information purposes Furusawa .

There are several approaches to produce TMSS states. Recently they have been produced experimentally for samples of atoms by quantum nondemolition (QND) measurements of the spin components Julsgaard . In Ref. Kuzmich  it was proposed to produce TMSS states of light by combining the two Einstein-Podolsky-Rosen (EPR) output modes of an optical parametric oscillator with coherent light at polarizing beam splitters.

It is also possible to adapt proposals for generating single-mode spin squeezing to the two-mode case. For example, it has been proposed that single-mode spin squeezed states can be produced via an adiabatic approach Sor01 . In this approach, one starts with a Hamiltonian for which the ground state is easily achieved, and slowly varies the Hamiltonian towards one that has as its ground state the optimal spin squeezed state. In principle, this procedure can be adapted to obtain optimal TMSS states. A drawback to this method is that, in practice, the appropriate Hamiltonian may be difficult to realize.

Another technique for generating TMSS states that may be generalized from the single-mode case is feedback. Feedback has proven to be a powerful technique for state preparation Wise94 ; Laura  and measurement Wise98 ; Berry01 . It has been shown that it is possible to produce near-optimal spin squeezing in the single-mode case by continual measurement of the spin operator $`J_z`$, and using feedback to bring $`J_z`$ close to zero Laura . Here we adapt this method to the QND measurement approach discussed in Ref. Julsgaard .

The method for production of TMSS states in Ref. Julsgaard  suffers the drawback that, although the variances in the sums of the spin components are small, their means are not close to zero. As is shown in Ref. Berry , the optimal TMSS states have zero means. We would therefore wish to obtain states with zero means in order for them to be as close as possible to the optimal states. This is analogous to the problem with the production of single-mode spin squeezing as considered by Thomsen, Mancini, and Wiseman (TMW) Laura .

## II two-mode spin squeezing

In Ref. Julsgaard , TMSS states (conditioned on the measurement record) are obtained via QND measurements of the spin of the atomic samples. In order to perform this QND measurement, an off-resonant pulse is transmitted through two cesium gas samples. The interaction Hamiltonian may be expressed as Kuzmich ; KMJYEB

$$H_{\mathrm{int}}=\mathrm{}aS_z(z,t)J_z(z,t)𝑑z,$$
(1)
where $`a`$ is a coupling constant, $`S_z(z,t)`$ is the $`z`$ component of the instantaneous Stokes vector for the light, and $`J_z(z,t)`$ is the $`z`$ component of the continuous spin operator. These variables are explained further in Appendix A.

In the Schrödinger picture $`J_z(z,t)`$ is independent of $`t`$, and the total Stokes vector $`𝐒(z)`$ is independent of $`z`$. The total time evolution operator is, therefore,

$`U`$ $`=\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}{\displaystyle H_{\mathrm{int}}𝑑t}\right)`$
$`=\mathrm{exp}\left(ia{\displaystyle S_z(z,t)J_z(z)𝑑z𝑑t}\right)`$
$`=\mathrm{exp}\left(ia{\displaystyle S_zJ_z(z)𝑑z}\right)`$
$`=\mathrm{exp}(iaS_zJ_z).`$ (2)
In Ref. Julsgaard  the light is initially linearly polarized, so the light is in the maximally weighted $`S_x`$ eigenstate. In addition, the total spin $`s`$ is large, so we may apply a SU(2)$``$HW(2) contraction. That is, we make the substitutions Holstein ; Berry

$`S_x`$ $`sb^{}b,`$ (3)
$`S_y`$ $`\sqrt{{\displaystyle \frac{s}{2}}}(b+b^{}),`$ (4)
$`S_z`$ $`{\displaystyle \frac{1}{i}}\sqrt{{\displaystyle \frac{s}{2}}}(bb^{}),`$ (5)
where $`b^{}b`$ is the number operator for photons subtracted from the $`x`$ linear polarization mode (and added to the $`y`$ linear polarization mode). We will assume that the probe beam is in a coherent spin state oriented in the $`x`$ direction (i.e., $`x`$ linearly polarized), and take the limit $`s\mathrm{}`$ so that this contraction is exact.

Under this contraction, $`[S_y,S_z]is`$. This means that, under the Heisenberg picture, the operator $`S_y`$ transforms as

$`S_y^{\mathrm{out}}=S_y^{\mathrm{in}}+[iaS_zJ_z,S_y]S_y^{\mathrm{in}}+{\displaystyle \frac{an}{2}}J_z,`$ (6)
where $`n=2s`$ is the total number of photons in the pulse. Now we consider two samples, where we will use the notation $`J_k^{(1)}`$ and $`J_k^{(2)}`$, where $`k\{x,y,z\}`$, for the spin operators for samples 1 and 2, respectively, and $`J_k^{(\pm )}=J_k^{(1)}\pm J_k^{(2)}`$. We may replace $`J_z`$ with $`J_z^{(+)}`$, the total $`z`$ component of spin for the two samples, giving

$$S_y^{\mathrm{out}}S_y^{\mathrm{in}}+\frac{an}{2}J_z^{(+)}.$$
(7)
Thus a measurement of $`S_y^{\mathrm{out}}`$ gives a QND measurement of $`J_z^{(+)}`$. In considering detection, it is more convenient to consider the Stokes vector for the light integrated over time interval $`\delta t`$. The transformation of this component is

$$\delta S_y^{\mathrm{out}}\delta S_y^{\mathrm{in}}+\frac{a\delta n}{2}J_z^{(+)},$$
(8)
where $`\delta n`$ is the total number of photons in time $`\delta t`$. In Ref. Julsgaard , TMSS states are obtained by performing QND measurements on $`J_y^{(+)}`$ as well. They achieve this by applying a magnetic field in the $`x`$ direction, in order to induce Larmor precession of the $`y`$ and $`z`$ components of the spin. This means that Eq. (8) becomes

$$\delta S_y^{\mathrm{out}}(t)\delta S_y^{\mathrm{in}}(t)+\frac{a\delta n}{2}[J_z^{(+)}\mathrm{cos}(\mathrm{\Omega }t)+J_y^{(+)}\mathrm{sin}(\mathrm{\Omega }t)],$$
(9)
where $`\mathrm{\Omega }`$ is the frequency of the precession. In this case measurement of $`\delta S_y^{\mathrm{out}}`$ alternately gives QND measurements of $`J_z^{(+)}`$ and $`J_y^{(+)}`$, thus giving reduced uncertainty in both $`J_z^{(+)}`$ and $`J_y^{(+)}`$.

For the TMSS states considered in Ref. Julsgaard  there is the requirement that $`J_x^{(1)}=J_x^{(2)}1`$. Here we will instead take the convention that both $`J_x^{(1)}`$ and $`J_x^{(2)}`$ are large and close to the total spin $`j`$, and replace $`J_y^{(+)}`$ with $`J_y^{()}`$. This is equivalent to a trivial rotation of coordinates for one of the modes and does not require a physical alteration to the experiment. If we now define

$$J_z^{(\mathrm{\Omega })}(t)J_z^{(+)}\mathrm{cos}(\mathrm{\Omega }t)+J_y^{()}\mathrm{sin}(\mathrm{\Omega }t),$$
(10)
then Eq. (9) becomes

$$\delta S_y^{\mathrm{out}}(t)\delta S_y^{\mathrm{in}}(t)+\frac{a\delta n}{2}J_z^{(\mathrm{\Omega })}(t).$$
(11)
In order for this to be accurate, we require $`\delta n`$ to be large, so that we may perform the SU(2)$``$HW(2) contraction. Under this contraction, $`\delta S_y=\sqrt{\delta s/2}X`$, where $`X=b+b^{}`$. Therefore the transformation in the quadrature $`X`$ is

$$X^{\mathrm{out}}(t)=X^{\mathrm{in}}(t)+a\sqrt{\delta n}J_z^{(\mathrm{\Omega })}(t).$$
(12)
As the initial state is equivalent to the vacuum state \[under the SU(2)$``$HW(2) contraction\], it has a variance of 1 in its measured value. The same is true for $`X^{\mathrm{out}}`$, as the extra term only changes the mean.

Measurements are made on the output beam by splitting the beam into $`+45^{}`$ and $`45^{}`$ linearly polarized components at a polarizing beam splitter and measuring the intensities at the two outputs. We then define the photocurrent as

$$I_c(t)=\underset{\delta t0}{lim}\frac{\delta N_+\delta N_{}}{\sqrt{\nu }\delta t},$$
(13)
where $`\delta N_+`$ and $`\delta N_{}`$ are the photon counts at the $`+45^{}`$ and $`45^{}`$ linearly polarized outputs and $`\nu `$ is the photon flux $`\delta n/\delta t`$. We may interpret $`I_c(t)`$ as the measured value of $`X`$ divided by $`\sqrt{\delta t}`$, so the photocurrent is given by

$$I_c(t)=a\sqrt{\nu }J_z^{(\mathrm{\Omega })}_c+\xi (t),$$
(14)
where the subscript $`c`$ indicates the average for the conditioned evolution, and $`\xi (t)`$ is a real Gaussian white noise term such that $`\xi (t)\xi (t^{})=\delta (tt^{})`$.

We take the limit of small time intervals $`\delta t`$, though for finite photon flux this would mean that the SU(2)$``$HW(2) contraction is no longer valid, as $`\delta n`$ would go to zero. In order for the limit to be rigorously correct, we need to take the limit of large $`\nu `$ as well as small $`\delta t`$. Nevertheless, this limit should be an accurate approximation for large photon flux. The conditioned master equation is then

$$d\rho _c=\mathrm{\Gamma }𝒟[J_z^{(\mathrm{\Omega })}]\rho _cdt+\sqrt{\mathrm{\Gamma }}dW(t)[J_z^{(\mathrm{\Omega })}]\rho _c,$$
(15)
where $`\mathrm{\Gamma }=a^2\nu /4`$, $`dW(t)=\xi (t)dt`$ is an infinitesimal Wiener increment, $`𝒟[r]\rho =r\rho r^{}(r^{}r\rho +\rho r^{}r)/2`$, and $`[r]\rho =r\rho +\rho r^{}\mathrm{Tr}[(r+r^{})\rho ]\rho `$.

## III Feedback

The quadrature (12) is measured by detecting the photocurrent (13), thereby yielding TMSS states. The drawback is that this will yield nonzero values of $`J_z^{(+)}`$ and $`J_y^{()}`$, whereas for the optimal TMSS states Berry , the expectation values are zero. In order to obtain states closer to the optimal TMSS states, we apply feedback to bring these expectation values closer to zero. For example, when $`\mathrm{cos}(\mathrm{\Omega }t)=1`$, we may apply a Hamiltonian proportional to $`J_y^{(+)}`$ in order to bring $`J_z^{(+)}`$ towards zero, analogous to the case considered by TMW. When $`\mathrm{sin}(\mathrm{\Omega }t)=1`$, we wish to apply a Hamiltonian proportional to $`J_z^{()}`$ in order to bring $`J_y^{()}`$ towards zero. For arbitrary $`t`$, we will apply the Hamiltonian

$$H_{\mathrm{fb}}=F(t)I_c(t),$$
(16)
where

$$F(t)=\frac{\lambda (t)}{\sqrt{\mathrm{\Gamma }}}J_y^{(\mathrm{\Omega })}(t),$$
(17)
and

$$J_y^{(\mathrm{\Omega })}(t)J_y^{(+)}\mathrm{cos}(\mathrm{\Omega }t)J_z^{()}\mathrm{sin}(\mathrm{\Omega }t).$$
(18)
We choose this definition for $`J_y^{(\mathrm{\Omega })}`$ because $`[J_z^{(\mathrm{\Omega })},J_y^{(\mathrm{\Omega })}]=iJ_x^{(+)}`$ \[we will omit the $`(t)`$ from this point for brevity\]. As $`J_x^{(+)}`$ is close to $`2j`$ for the states we are considering, this means that rotation about $`J_y^{(\mathrm{\Omega })}`$ will alter the expectation value of $`J_z^{(\mathrm{\Omega })}`$. This Hamiltonian may be implemented by a radio frequency magnetic field Sang . This generates a Hamiltonian proportional to $`J_y`$, which becomes $`J_y^{(\mathrm{\Omega })}`$ due to the Larmor precession of the atomic spin.

Using this feedback, we obtain the unconditioned master equation master

$$\dot{\rho }=i[\frac{1}{2}\{c^{}F(t)+F(t)c\},\rho ]+𝒟[ciF(t)]\rho ,$$
(19)
with

$$c=\sqrt{\mathrm{\Gamma }}J_z^{(\mathrm{\Omega })}.$$
(20)
In the single-mode case TMW pointed out that the term $`c^{}F(t)+F(t)c`$ is proportional to $`J_zJ_y+J_yJ_z`$, which is the two-axis countertwisting Hamiltonian Kit . In a similar way we can derive a Hamiltonian for the two-mode case that produces TMSS states and is analogous to the two-axis countertwisting Hamiltonian. It is straightforward to show that

$`c^{}F(t)+F(t)c`$ $`=\lambda (t)\{[(J_y^{(1)})^2+(J_y^{(2)})^2(J_z^{(1)})^2`$
$`(J_z^{(2)})^2]\mathrm{sin}(2\mathrm{\Omega }t)+(J_z^{(1)}J_y^{(1)}+J_y^{(1)}J_z^{(1)}`$
$`+J_z^{(2)}J_y^{(2)}+J_y^{(2)}J_z^{(2)})\mathrm{cos}(2\mathrm{\Omega }t)`$
$`+2(J_z^{(1)}J_y^{(2)}+J_y^{(1)}J_z^{(2)})\}.`$ (21)
In the limit of large $`\mathrm{\Omega }`$, the contribution to the evolution from the terms containing $`\mathrm{sin}(2\mathrm{\Omega }t)`$ and $`\mathrm{cos}(2\mathrm{\Omega }t)`$ is negligible. Omitting these terms yields

$$c^{}F(t)+F(t)c2\lambda (t)(J_z^{(1)}J_y^{(2)}+J_y^{(1)}J_z^{(2)}),$$
(22)
which indicates that it is possible to produce two-mode spin squeezing using the Hamiltonian

$$HJ_z^{(1)}J_y^{(2)}+J_y^{(1)}J_z^{(2)}.$$
(23)
### III.1 Simple feedback

Now we will consider the appropriate feedback strength, $`\lambda (t)`$. It is easily shown that the variation in $`J_z^{(\mathrm{\Omega })}_c`$ due to the conditioned evolution (15) is

$$\mathrm{Tr}(J_z^{(\mathrm{\Omega })}d\rho _c)=2\sqrt{\mathrm{\Gamma }}dW[(J_z^{(\mathrm{\Omega })})^2_cJ_z^{(\mathrm{\Omega })}_c^2].$$
(24)
The variation in $`J_z^{(\mathrm{\Omega })}_c`$ due to the feedback is

$`i[H_{\mathrm{fb}},J_z^{(\mathrm{\Omega })}]_c`$ $`=i{\displaystyle \frac{\lambda }{\sqrt{\mathrm{\Gamma }}}}I_c[J_y^{(\mathrm{\Omega })},J_z^{(\mathrm{\Omega })}]_c`$
$`={\displaystyle \frac{\lambda }{\sqrt{\mathrm{\Gamma }}}}I_cJ_x^{(+)}_c.`$ (25)
If $`J_z^{(\mathrm{\Omega })}_c=0`$, then the appropriate feedback to keep this equal to zero is

$$\lambda (t)=\frac{2\mathrm{\Gamma }(J_z^{(\mathrm{\Omega })})^2_c}{J_x^{(+)}_c}.$$
(26)
If the unconditioned state is close to being pure, then we may also use this expression with the unconditioned averages without significant loss of accuracy:

$$\lambda (t)=\frac{2\mathrm{\Gamma }(J_z^{(\mathrm{\Omega })})^2}{J_x^{(+)}}.$$
(27)
Note that this result is comparable to that found for the case of feedback for single-mode spin squeezing by TMW:

$$\lambda (t)=\frac{2\mathrm{\Gamma }J_z^2}{J_x}.$$
(28)
This is due to the fact that the commutation relations for $`J_x`$, $`J_y`$, and $`J_z`$ are identical to those for $`J_x^{(+)}`$, $`J_y^{(\mathrm{\Omega })}`$, and $`J_z^{(\mathrm{\Omega })}`$.

### III.2 Analytic feedback

Next we will consider the possibilities for using an analytic expression for the feedback strength $`\lambda (t)`$. In Ref. Laura  the approximate relations

$`J_z^2`$ $`(4\mathrm{\Gamma }t+2/j)^1,`$ (29)
$`J_x`$ $`j\mathrm{exp}(\mathrm{\Gamma }t/2),`$ (30)
are used to derive the analytic feedback scheme

$$\lambda (t)=\frac{\mathrm{\Gamma }\mathrm{exp}(\mathrm{\Gamma }t/2)}{1+2j\mathrm{\Gamma }t}.$$
(31)
In the two-mode case, it is possible to obtain the corresponding relation for $`J_x^{(+)}`$, but not for $`(J_z^{(\mathrm{\Omega })})^2`$.

In these derivations it is convenient to define the scaled variables $`v=\mathrm{\Gamma }t`$ and $`\mathrm{\Lambda }(t)=\lambda (t)/\mathrm{\Gamma }`$. In terms of these variables the master equation (19) becomes

$`\dot{\rho }`$ $`={\displaystyle \frac{i\mathrm{\Lambda }(t)}{2}}[(J_z^{(\mathrm{\Omega })}J_y^{(\mathrm{\Omega })}+J_y^{(\mathrm{\Omega })}J_z^{(\mathrm{\Omega })}),\rho ]+𝒟[J_z^{(\mathrm{\Omega })}`$
$`i\mathrm{\Lambda }(t)J_y^{(\mathrm{\Omega })}]\rho ,`$ (32)
where the dot indicates a derivative with respect to $`v`$. Using these variables, the evolution is independent of $`\mathrm{\Gamma }`$. The value of $`\mathrm{\Gamma }`$ does not qualitatively affect the evolution and merely provides a scaling to the time and to $`\lambda (t)`$.

It is straightforward to show that

$`{\displaystyle \frac{d}{dv}}J_x^{(+)}`$ $`=\mathrm{Tr}(J_x^{(+)}\dot{\rho })`$
$`=\frac{1}{2}J_x^{(+)}+2\mathrm{\Lambda }(J_z^{(\mathrm{\Omega })})^2\frac{1}{2}\mathrm{\Lambda }^2J_x^{(+)}`$
$`=\frac{1}{2}J_x^{(+)}+{\displaystyle \frac{2(J_z^{(\mathrm{\Omega })})^2^2}{J_x^{(+)}}}.`$ (33)
Initially this derivative is equal to zero. However, the value of $`(J_z^{(\mathrm{\Omega })})^2`$ falls rapidly and that of $`J_x^{(+)}`$ only falls slowly, so the second term quickly becomes negligible. If we neglect the second term, the evolution of $`J_x^{(+)}`$ is

$$J_x^{(+)}2j\mathrm{exp}(v/2),$$
(34)
which is equivalent to that in the single-mode case. This technique fails to find the corresponding expression for $`(J_z^{(\mathrm{\Omega })})^2`$, but we have found numerically that a good approximation is given by

$$(J_z^{(\mathrm{\Omega })})^2\frac{e^{v/4}}{2v+1/j}.$$
(35)
Using these results in the expression for the feedback (27), we obtain

$$\lambda (t)=\frac{\mathrm{\Gamma }\mathrm{exp}(\mathrm{\Gamma }t/4)}{1+2j\mathrm{\Gamma }t}.$$
(36)
### III.3 Optimal feedback

The third alternative that we will consider is optimal feedback. In order to derive this, it is convenient to define the variables

$$\zeta =\frac{(J_z^{(+)})^2+(J_y^{()})^2}{2j},\chi =\frac{J_x^{(+)}}{2j}.$$
(37)
The initial state used is that with individual coherent spin states in each arm, so both $`\zeta `$ and $`\chi `$ are equal to 1. As time progresses both $`\chi `$ and $`\zeta `$ decrease. In order for the feedback to be optimal, the value of $`\zeta `$ should be as small as possible for each value of $`\chi `$. Therefore, if $`\chi `$ decreases by an amount $`\delta \chi `$, then the amount by which $`\zeta `$ decreases should be maximized. This means that the feedback that is optimum will be that which produces the maximum slope. Therefore, in order to determine the optimal feedback, we wish to solve

$$\frac{d}{d\mathrm{\Lambda }}\frac{\zeta }{\chi }=0.$$
(38)
This can be solved using

$$\frac{d}{d\mathrm{\Lambda }}\left[\frac{\frac{d}{dv}(J_z^{(+)})^2+(J_y^{()})^2}{\frac{d}{dv}J_x^{(+)}}\right]=0.$$
(39)
Using the master equation (19), it is straightforward to show that

$$\mathrm{Tr}\left[(J_x^{(+)})\dot{\rho }\right]=\frac{1}{2}[1+\mathrm{\Lambda }^2]J_x^{(+)}+2\mathrm{\Lambda }(J_z^{(\mathrm{\Omega })})^2.$$
(40)
Determining the corresponding expression for the numerator in Eq. (39) is more difficult, and if it is performed directly it leads to an extremely complicated expression. However, there is a more convenient way of determining this expression. Note that

$`(J_z^{(\mathrm{\Omega })})^2`$ $`=(J_z^{(+)})^2\mathrm{cos}^2(\mathrm{\Omega }t)+(J_y^{()})^2\mathrm{sin}^2(\mathrm{\Omega }t)+(J_y^{()}J_z^{(+)}`$
$`+J_z^{(+)}J_y^{()})\mathrm{cos}(\mathrm{\Omega }t)\mathrm{sin}(\mathrm{\Omega }t).`$ (41)
For large $`\mathrm{\Omega }`$ we may average over $`t`$, which gives

$$(J_z^{(\mathrm{\Omega })})^2\frac{1}{2}(J_z^{(+)})^2+(J_y^{()})^2.$$
(42)
Numerically this approximation was found to be very accurate. Using this result we find that

$`\frac{1}{2}\mathrm{Tr}[((J_z^{(+)})^2+(J_y^{()})^2)\dot{\rho }]\mathrm{Tr}[(J_z^{(\mathrm{\Omega })})^2\dot{\rho }]`$
$`=\mathrm{\Lambda }^2(J_x^{(+)})^2(J_z^{(\mathrm{\Omega })})^2\mathrm{\Lambda }4J_z^{(\mathrm{\Omega })}J_x^{(+)}J_z^{(\mathrm{\Omega })}+J_x^{(+)}.`$ (43)
Thus we find that Eq. (39) becomes

$$\frac{d}{d\mathrm{\Lambda }}\left\{\frac{\mathrm{\Lambda }^2d\mathrm{\Lambda }e}{f\mathrm{\Lambda }^2f+\mathrm{\Lambda }g}\right\}=0,$$
(44)
where

$`d`$ $`=(J_x^{(+)})^2(J_z^{(\mathrm{\Omega })})^2,`$
$`e`$ $`=4J_z^{(\mathrm{\Omega })}J_x^{(+)}J_z^{(\mathrm{\Omega })}+J_x^{(+)},`$
$`f`$ $`=\frac{1}{2}J_x^{(+)},`$
$`g`$ $`=2(J_z^{(\mathrm{\Omega })})^2.`$ (45)
Solving this gives

$`\mathrm{\Lambda }={\displaystyle \frac{fd\pm \sqrt{(fd)^2+ef(fedg)}}{fedg}}.`$ (46)
Numerically it is found that the correct solution that maximizes the slope is that with the positive sign.

Note that this derivation of the optimum feedback does not rely on any commutation relations other than the usual su(2) commutation relations. It therefore is also applicable to the case of feedback for single-mode spin squeezing, with $`J_x^{(+)}`$, $`J_y^{(\mathrm{\Omega })}`$, and $`J_z^{(\mathrm{\Omega })}`$ replaced by $`J_x`$, $`J_y`$, and $`J_z`$.

## IV Numerical results

### IV.1 Unconditioned master equation

The unconditioned master equation (III.2) was solved using a simple, finite step method with step sizes of $`\delta v=1/1000`$. It was found that better results were obtained as $`\mathrm{\Omega }`$ was increased. In order to estimate the best result in the limit of large $`\mathrm{\Omega }`$, $`\mathrm{\Omega }`$ was assigned the maximum value possible with this step size, $`\pi /2\delta v`$. The initial condition used was that where the individual modes were in independent coherent spin states oriented along the $`x`$ axis.

At each time step the quantities $`\zeta `$, $`\chi `$, and the purity

$$P=\mathrm{Tr}[\rho ^2]$$
(47)
were calculated. The results for the simple feedback scheme of Eq. (27) and a spin of $`j=5`$ are shown in Fig. 1. Initially the system behaves as we would expect, with the sum of the variances $`\zeta `$ decreasing, until a time of approximately $`v=3`$. Then the variances dramatically increase, and the purity drops almost to zero. After this, however, the system stabilizes. Note also that the value of $`\chi `$ drops regularly as time increases until the change at $`v=3`$.

Recall that the system must be entangled if $`(J_z^{(+)})^2+(J_y^{()})^2<J_x^{(+)}`$ Berry , which is equivalent to $`\zeta <\chi `$. Therefore a state with $`\zeta <\chi `$ can be described as TMSS. Evidently the feedback produces quite strong spin squeezing. In order to determine how close to optimal the states produced by this feedback are, the value of $`\zeta `$ was plotted against $`\chi `$ and compared with the plot for optimal TMSS states in Fig. 2. As time progresses the state travels from the upper right corner towards the lower left corner of the figure.

The state starts out very close to optimum and does not deviate significantly from optimum until $`\zeta <0.4`$. Eventually the value of $`\zeta `$ increases dramatically, and the system spirals towards its equilibrium state at $`(0.5,0.5)`$. The results for $`j=5`$ are typical, and similar results are obtained for other spins $`j>1/2`$.

For the purposes of teleportation, the important issue is how close it is possible to get to the optimal TMSS states considered for teleportation in Ref. Berry . This state is indicated by the cross in Fig. 2. As can be seen, the states obtained by feedback closely approach this optimal state. Similar results are obtained for higher spin, but for smaller spin the states obtained by feedback are further from the TMSS states considered for teleportation.

The plots of $`\zeta `$ versus $`\chi `$ for analytic feedback and optimal feedback are also shown in Fig. 2. The results for analytic feedback are very similar to those for simple feedback, with $`\zeta `$ dramatically increasing at later times. In contrast, under optimal feedback the value of $`\zeta `$ monotonically decreases to an asymptotic value for $`\chi =0`$. For small times the three feedback schemes produce very similar results, and the states obtained using optimal feedback do not pass significantly closer to the optimal state for teleportation.

The last alternative that we consider in this section is the Hamiltonian (23), which is analogous to the two-axis countertwisting Hamiltonian in the single-mode case. The results for this Hamiltonian are also shown in Fig. 2. For this case there is significant spin squeezing, but the squeezing is generally poorer than for feedback. Nevertheless, at later times the value of $`\zeta `$ does not rise dramatically (until $`\chi `$ falls below zero), as opposed to the results for the simple or analytic feedback.

### IV.2 Conditioned master equation

Next we consider the results for the conditioned master equation. There are two examples where it is necessary to consider the conditioned master equation. In the case of simple feedback, at the time that $`\zeta `$ increases dramatically the purity of the state has greatly decreased. This means that the feedback based on the unconditioned averages, Eq. (27), is a poor approximation of the feedback based on conditioned averages, Eq. (26), and will not accurately keep the average $`J_z^{(\mathrm{\Omega })}_c`$ equal to zero. For this reason we may obtain better results using the conditioned master equation and the feedback (26).

The second case that we consider here is that without any feedback. In this case there is no reduction in the variances for the unconditioned equation. If we consider the conditioned equation, however, there is a reduction in the variances, though the means $`J_z^{(+)}_c`$ and $`J_y^{()}_c`$ are nonzero. For this reason we will consider the conditioned master equation for this case, with $`\zeta `$ defined by

$$\zeta =\frac{(J_z^{(+)}J_z^{(+)}_c)^2+(J_y^{()}J_y^{()}_c)^2_c}{2j}.$$
(48)
To minimize the random differences between the two cases the same random numbers were used for each. The results for these two cases are plotted in Fig. 3.

The results for the feedback (26) are extremely close to the line for optimal states, with only temporary excursions from it. In particular, the value of $`\zeta `$ does not increase dramatically at large times. This indicates that improved results may be obtained using feedback based on the conditioned state.

On the other hand, implementing this feedback may be far more difficult. In the forms of feedback considered for the unconditioned equation, the feedback was proportional to the measurement record, with a proportionality constant $`\lambda `$ that is a function of time only. This feedback may be performed efficiently by calculating $`\lambda (t)`$ beforehand and programming it into a field programmable gate array FPGA . In contrast, the value of $`\lambda `$ given by Eq. (26) is dependent on the measurement record, and therefore would need to be calculated during the measurement, introducing a significant time delay.

The results without feedback quickly diverge from the line for optimal states, although they do not differ from it greatly. The major difference is at later times. While, with feedback, the value of $`\zeta `$ is reduced very close to zero, without feedback the value of $`\zeta `$ cannot be reduced as far. This means that, although a significant reduction in the variances can be achieved via feedback for strong QND measurements (where $`\mathrm{\Gamma }t>1`$), feedback does not give much improvement for weak measurements.

## V Comparison with single-mode case

There is clearly a great similarity between the single-mode and two-mode cases. The operators $`J_x^{(+)}`$, $`J_y^{(\mathrm{\Omega })}`$, and $`J_z^{(\mathrm{\Omega })}`$ obey exactly the same commutation relations as $`J_x`$, $`J_y`$, and $`J_z`$ in the single-mode case. This means that, for example, the optimal feedback also applies to the single-mode case (as was mentioned in Sec. III.3). As we may replace the sum of the variances $`(J_z^{(+)})^2+(J_y^{()})^2`$ with $`2(J_z^{(\mathrm{\Omega })})^2`$, it might appear that this case is identical to the single-mode case. Nevertheless, there is a subtle difference due to the fact that the operator $`J_z^{(\mathrm{\Omega })}`$ is time dependent.

For the case of spin 1/2, there is complete equivalence between the single-mode case (for spin $`j=1`$) and the two-mode case. We will take the Hilbert space for the single-mode case to be that for two spin 1/2 particles, so that it is the same as for the two-mode case. The optimal TMSS states for $`j=1/2`$ are also optimal single-mode squeezed states in terms of the total spin. That is, they minimize $`(J_z^{(+)})^2`$ for a given $`J_x^{(+)}`$. The converse is also true: the states optimized for single-mode spin squeezing are automatically optimized for two-mode spin squeezing.

In addition, both the single- and two-mode countertwisting Hamiltonians produce optimal spin squeezed states. Not only this, but the feedback as given by Eq. (28) in the single-mode case, or Eq. (27) in the two-mode case, and the optimal feedback described in Sec. III.3, give optimal spin squeezed states. The only feedback that does not give optimal states is the analytic feedback for the single- or two-mode case. These results are depicted in Fig. 4. In this figure the variables $`\zeta `$ and $`\chi `$ are defined for the single-mode case as

$$\zeta =2J_z^2/j,\chi =J_x/j.$$
(49)
To explain these results analytically, first note that, for spin 1/2,

$$J_z^{(+)}J_y^{(+)}+J_y^{(+)}J_z^{(+)}=2(J_z^{(1)}J_y^{(2)}+J_y^{(1)}J_z^{(2)}).$$
(50)
This means that the countertwisting Hamiltonians for the one- and two-mode cases are identical, and therefore produce identical states.

In addition we find that $`(J_z^{(+)})^2(J_y^{()})^2`$ commutes with these Hamiltonians. This means that for the states produced by countertwisting, the values of $`(J_z^{(+)})^2`$ and $`(J_y^{()})^2`$ will be identical. This means that, if the states produced are optimal single-mode spin squeezed states (i.e., $`(J_z^{(+)})^2`$ is minimized), then they must also be optimal TMSS states.

To show that optimal single-mode spin squeezed states are produced by the countertwisting Hamiltonian, note that for total spin 1 we have the differential equations

$`{\displaystyle \frac{d}{dv}}J_x`$ $`=2J_z^2J_y^2,`$ (51)
$`{\displaystyle \frac{d}{dv}}J_z^2`$ $`=J_x,`$ (52)
$`{\displaystyle \frac{d}{dv}}J_y^2`$ $`=J_x.`$ (53)
The solution of these equations is

$$J_x^2=4(J_z^2J_z^2^2).$$
(54)
This is the relation for optimal single-mode spin squeezed states as given by Sørensen and Mølmer Sor01 . This shows that the states produced by the one- and two-mode countertwisting Hamiltonians and the one- and two-mode optimal states are all identical.

For the case of feedback, it is sufficient to show that the feedbacks of Eqs. (28) and (27) give optimal states, as this will imply that the optimal feedback gives optimal spin squeezed states. The derivation in this case is lengthy, and is given in Appendix B. It is interesting to note that in this case the optimal feedback is given by Eq. (65), and differs dramatically from the analytic feedback used for higher spin.

The equivalence between the single- and two-mode cases does not hold for any total spin higher than 1. As shown in Fig. 5, there are significant differences between the optimal single- and two-mode spin squeezed states for total spin above 1. As can be seen, for a total spin of 2 or more the values of $`\zeta `$ are significantly less for the single-mode case than for the two-mode case. This reflects the fact that in the two-mode case we wish to minimize $`(J_y^{()})^2`$ as well as $`(J_x^{(+)})^2`$.

As there are such strong similarities between the single-mode and two-mode cases, it is reasonable to apply the squeezing parameter as considered in Refs. Laura ; Sorensen ; Wang  to the two-mode case. In the single-mode case the squeezing parameter was defined by

$$\xi _z^2=2jJ_z^2/J_x^2.$$
(55)
With the definitions of $`\zeta `$ and $`\chi `$ for the single-mode case above, this can be expressed as $`\xi _z^2=\zeta /\chi ^2`$. It would seem reasonable to use the same definition for the two-mode case.

For the two-mode case we find that the minimum value of $`\xi _z^2`$, namely $`\xi _{\mathrm{min}}^2`$, varies with spin $`j`$ as shown in Fig. 6(a). In this figure the values of $`\xi _{\mathrm{min}}^2`$ have been multiplied by $`j+1`$ in order to more easily compare the results for different spins. Very similar results are obtained for the simple feedback of Eq. (27) and the optimal feedback of Eq. (46). Using the analytic feedback of Eq. (36) gives results that are similar to, but slightly above, those for the other two feedback schemes. All of these feedback schemes give results that are significantly above the result for the optimal states. Using the countertwisting Hamiltonian (23) gives results that are higher than those for feedback.

Similar results are obtained for the single-mode case \[see Fig. 6(b)\]. The results for feedback are noticeably above those for optimal states, and the results for two-axis countertwisting are significantly above those for feedback. To summarize, the scaling constants for each of the cases are as given in Table 1. In this table we can see the similarities between the single-mode and two-mode cases: the scaling constants obtained by feedback are significantly above those for optimal states, and the scaling constants for two-axis countertwisting (or its two-mode equivalent) are significantly above those for feedback. In both cases there are only small differences between the results for the different types of feedback. The only qualitative difference between the scaling constants for the single- and two-mode cases are that in the single-mode case the analytic feedback gives slightly better results than the simple feedback of Eq. (27). In contrast, in the two-mode case the scaling constants are indistinguishable.

## VI Experimental prospects

There are some complications in applying this theory to the experiment of Ref. Julsgaard . In this experiment, $`a5\times 10^{13}`$, and $`\nu 2\times 10^{16}\mathrm{s}^1`$. This means that $`\mathrm{\Gamma }1.4\times 10^9\mathrm{s}^1`$, which implies that the time required for maximal squeezing is on the order of $`10^9\mathrm{s}`$ (over 20 yr). Therefore a far more intense beam or a stronger interaction would be required to obtain the two-mode spin squeezing described here.

It must be emphasized that this long time is required for the measurements, not the feedback. Feedback may be applied to the states obtained in Ref. Julsgaard  with a negligible increase in the time required for the experiment. Nevertheless, for the conditions of this experiment, although the feedback will bring the means of $`J_z^{(+)}`$ and $`J_y^{(+)}`$ towards zero, it will not significantly reduce the variances. It is only when the QND measurement can be performed strongly enough (i.e., with a strong enough interaction over a long enough time) to obtain spin squeezing close to maximum that an improvement in the variance is obtained by using feedback.

Another problem is spontaneous losses due to absorption of QND probe light. The loss rate due to this is $`N\gamma g^2n/4\mathrm{\Delta }^2`$, where $`N=4j`$ is the total number of atoms and $`g`$ is the one-photon Rabi frequency. As in the case of single-mode spin squeezing Laura , this loss will be very large over the time period $`1/\mathrm{\Gamma }`$ for free space. Both this problem and the problem of the long interaction time may be overcome by performing the experiment in a cavity in the strong-coupling regime.

Two other common experimental problems are inefficient detectors and time delays. As for the case of single-mode spin squeezing, inefficient detectors do not affect the scaling. The system has the potential to be far more sensitive to time delays than in the single-mode case, however. The problem is that, as the spin component that is being measured is rotating rapidly, the feedback may be correcting for measurements of a different spin component.

In order to correct the right spin component, the rotation of the spin component that is measured should have completed an integral number of rotations during the time delay. Experimentally, this means that the magnetic field should be adjusted such that the frequency $`\mathrm{\Omega }`$ is an integral multiple of $`2\pi /\tau `$, where $`\tau `$ is the time delay. Provided that this is done, time delays should not be more of a problem than in the single-mode case.

## VII Conclusions

Two-mode spin squeezed states are important states to produce for quantum teleportation. Here we have shown that it is possible to produce states very close to the optimal TMSS states by adapting the feedback for single-mode spin squeezing considered by TMW. These states are not conditioned on the measurement record, in contrast to the conditional two-mode spin squeezing discussed in Ref. Julsgaard .

Using the simple feedback scheme (27) it is possible to obtain states that are quite close to optimal TMSS states for small times, but strongly diverge from these states at later times. In particular, for spins above about 5 it is possible to obtain states very close to the TMSS states considered for teleportation in Ref. Berry . An analytic feedback scheme (36) also gives similar results. This feedback scheme is more practical experimentally, as the appropriate feedback strength to be used is easily calculated.

We have derived the optimal feedback that produces the maximum possible spin squeezing. This feedback is also applicable to the case of feedback for single-mode spin squeezing. This feedback gives states even closer to the optimal TMSS states.

We have also derived a Hamiltonian for the two-mode case that is equivalent to the two-axis countertwisting Hamiltonian introduced in Ref. Kit . This Hamiltonian produces spin squeezing, but not as strongly as the measurements with feedback.

In the case of spin 1/2 both the feedback (except for the analytic feedback) and the countertwisting Hamiltonian produce optimal TMSS states. These states are equivalent to optimal single-mode spin squeezed states if the two modes are considered as a single spin-1 system. In the single-mode case optimal spin squeezed states are produced both by feedback and by the countertwisting Hamiltonian.

###### Acknowledgements.

The authors acknowledge valuable discussions with Laura Thomsen and Howard Wiseman. We are also grateful for constructive criticism from Eugene Polzik. This project has been supported by an Australian Research Council Large Grant.

## Appendix A Relation of variables to experiment

Here we give further explanation of the variables introduced in Eq. (1). The coupling constant $`a`$ is given by

$$a=\frac{\sigma }{A(I+1/2)}\frac{\gamma }{\mathrm{\Delta }}\alpha _v,$$
(56)
where $`\sigma `$ is the resonant absorption cross section, $`A`$ is the area of the transverse cross section of the light beam, $`I`$ is the nuclear spin, $`\gamma `$ is the spontaneous emission rate of the upper atomic level, $`\mathrm{\Delta }`$ is the detuning, and $`\alpha _v`$ is the dynamic vector polarizability. We have omitted the bounds on the integral (1) for greater generality, so we may apply this expression to multiple samples. Explicit bounds are unnecessary, as there is no contribution to the integral from regions where there are no atoms.

The continuous spin operators for the sample are defined as

$$J_k(z,t)=\underset{\delta z0}{lim}\frac{1}{\delta z}\underset{\mu }{}\frac{1}{2}\sigma _k^\mu ,$$
(57)
where $`k\{x,y,z\}`$ and $`\sigma _k^\mu `$ is the Pauli operator for particle $`\mu `$. The sum is over all particles between $`z`$ and $`z+\delta z`$ over the cross section of the sample. The operator for the entire sample is obtained by integrating over $`z`$.

The field is described by the continuous-mode annihilation operators $`a_k(z,t)`$, where $`k=x`$ and $`y`$ for the $`x`$ polarized and $`y`$ polarized modes, respectively. The instantaneous Stokes parameters are

$`S_x(z,t)`$ $`=\frac{1}{2}[a_x^{}(z,t)a_x(z,t)a_y^{}(z,t)a_y(z,t)],`$
$`S_y(z,t)`$ $`=\frac{1}{2}[a_x^{}(z,t)a_y(z,t)+a_y^{}(z,t)a_x(z,t)],`$
$`S_z(z,t)`$ $`={\displaystyle \frac{i}{2}}[a_x^{}(z,t)a_y(z,t)a_y^{}(z,t)a_x(z,t)].`$ (58)
The Stokes vector for the entire pulse at position $`z`$ is given by $`𝐒(z)=𝐒(z,t)𝑑t`$.

## Appendix B Feedback for total spin 1

Here we show that the feedback of either Eqs. (28) or (27) gives optimal spin squeezed states for a total spin of 1. To see this, note first that the first term in Eq. (19) is just the same as that produced by the countertwisting Hamiltonian, and so will produce optimal states. In order to show that the feedback produces optimal states, it remains to be shown that the second term, $`𝒟[ciF]\rho `$, does not alter the evolution. Specifically, in the single-mode case

$$\mathrm{Tr}\left\{J_x𝒟[J_zi\mathrm{\Lambda }J_y]\rho \right\}=\frac{1}{2}(1+\mathrm{\Lambda }^2)J_x+\mathrm{\Lambda }J_z^2+J_y^2.$$
(59)
If the state $`\rho `$ is an optimal state, then $`J_z^2+J_y^2=1`$, so this simplifies to

$$\mathrm{Tr}\left\{J_x𝒟[J_zi\mathrm{\Lambda }J_y]\rho \right\}=\frac{1}{2}(1+\mathrm{\Lambda }^2)J_x+\mathrm{\Lambda }.$$
(60)
Similarly we can show that

$$\mathrm{Tr}\left\{J_z^2𝒟[J_zi\mathrm{\Lambda }J_y]\rho \right\}=\frac{\mathrm{\Lambda }}{2}J_x+\mathrm{\Lambda }^2J_x^2J_z^2.$$
(61)
For optimal states, $`J_x^2=1`$, so this becomes

$$\mathrm{Tr}\left\{J_z^2𝒟[J_zi\mathrm{\Lambda }J_y]\rho \right\}=\frac{\mathrm{\Lambda }}{2}J_x+\mathrm{\Lambda }^2(1J_z^2).$$
(62)
It is simple to show that both Eqs. (60) and (62) will be zero if Eq. (54) is satisfied, and the feedback is given by Eq. (28).

Therefore, if the state is already in an optimal state, and the feedback as given by Eq. (28) is used, then the state will continue to be in an optimal state. On the other hand, if some other feedback is used, then optimal states will not be obtained. For example, the analytic feedback considered by LMW does not give optimal states. In order to determine analytic feedback that will give optimal states, note that the differential equations for $`J_x`$ and $`J_z^2`$ are

$`{\displaystyle \frac{d}{dv}}J_x`$ $`=\mathrm{\Lambda }(2J_z^21),`$
$`{\displaystyle \frac{d}{dv}}J_z^2`$ $`={\displaystyle \frac{\mathrm{\Lambda }}{2}}J_x.`$ (63)
Solving these using the feedback (28) and Eq. (54) gives

$`J_x`$ $`=\sqrt{2e^ve^{2v}},`$
$`J_z^2`$ $`=\frac{1}{2}e^v.`$ (64)
This implies that the value of $`\mathrm{\Lambda }`$ should change with time as

$$\mathrm{\Lambda }(v)=\frac{1}{\sqrt{2e^v1}}.$$
(65)
This is quite different from the analytic expression applied for larger spins.

The case for feedback for two-mode spin squeezing is analogous to this, except that $`J_x`$, $`J_y`$, and $`J_z`$ are replaced with $`J_x^{(+)}`$, $`J_y^{(\mathrm{\Omega })}`$, and $`J_z^{(\mathrm{\Omega })}`$. This therefore shows that optimal states are obtained in the one- and two-mode cases using the feedback given by Eqs. (28) and (27), respectively.