text:
# Quantum Kalman Filtering and the Heisenberg Limit in Atomic Magnetometry

## Abstract

The shotnoise detection limit in current high-precision magnetometry Kominis et al. (2003) is a manifestation of quantum fluctuations that scale as $`1/\sqrt{N}`$ in an ensemble of $`N`$ atoms. Here, we develop a procedure that combines continuous measurement and quantum Kalman filtering Belavkin (1999) to surpass this conventional limit by exploiting conditional spin-squeezing to achieve $`1/N`$ field sensitivity. Our analysis demonstrates the importance of optimal estimation for high bandwidth precision magnetometry at the Heisenberg limit and also identifies an approximate estimator based on linear regression.

Magnetometry is fundamentally a parameter estimation process because, like all fields, magnetism cannot be directly observed. Rather, the strength of a magnetic field must be inferred from its influence on a probe such as an atomic spin ensemble Dupont-Roc et al. (1969). In a canonical atomic magnetometer, such an ensemble would be prepared into a coherent spin state with its bulk magnetization polarized along the $`x`$-axis, $`\widehat{𝐉}(0)=(J,0,0)`$ (such as by optical pumping). Then, a magnetic field along the $`y`$-axis with magnitude, $`B`$, would induce the atomic Bloch vector, $`\widehat{𝐉}(t)`$, to precess in the $`xz`$-plane with frequency, $`\omega _\mathrm{L}=\gamma B`$. Thus, the magnetic field could be estimated from the free induction decay of the atomic magnetization by monitoring the $`z`$-component of the Bloch vector, $`\widehat{J}_\mathrm{z}(t)=J\mathrm{exp}(t/T_2)\mathrm{sin}(\omega _\mathrm{L}t)`$, where $`T_2`$ is the transverse spin coherence time.

In practice, current atomic magnetometers operate by continuously pumping the atomic sample while a $`\widehat{J}_\mathrm{z}`$-dependent optical property of the ensemble is monitored Kominis et al. (2003); Budker et al. (2002); Kim and Lee (1998). Due to pumping, the atoms are constantly re-polarized as they Larmor precess. For small fields (the relevant case when discussing detection limits), the ensemble rapidly achieves an equilibrium that is nearly polarized along the $`x`$-axis, but with a steady-state offset, $`\widehat{J}_\mathrm{z}_{\mathrm{ss}}\gamma BJ`$. The uncertainty in measuring $`\widehat{J}_\mathrm{z}`$, is due to projection noise Itano et al. (1993), $`\mathrm{\Delta }\widehat{J}_\mathrm{z}^2\widehat{J}_\mathrm{z}^2\widehat{J}_\mathrm{z}^2`$, which has a value of $`J/2`$ for a coherent spin state. Averaging a sequence of independent measurements of $`\widehat{J}_\mathrm{z}_{\mathrm{ss}}`$ with this variance leads to the conventional shotnoise detection limit for a total measurement time of $`t_{\mathrm{tot}}`$ Kominis et al. (2003); Budker et al. (2002),

$$\delta B\frac{1}{\gamma \sqrt{JT_2t_{\mathrm{tot}}}}.$$
(1)
Since $`\mathrm{\Delta }\widehat{J}_\mathrm{z}^2`$ sets an intrinsic limit on the field sensitivity, reducing the projection noise below its standard quantum limit would improve the precision. This naturally leads one to consider spin-squeezed states Kitagawa and Ueda (1993) where uncertainty in $`\widehat{J}_\mathrm{z}`$ is reduced by redistributing it into the orthogonal spin component so that $`\mathrm{\Delta }\widehat{J}_\mathrm{y}^2>J/2`$. Since $`\widehat{J}_\mathrm{y}`$ does not directly affect the field estimation, spin-squeezing should enable one to surpass the conventional shotnoise magnetometry limit.

An improved magnetometry protocol would ideally be implemented by utilizing the conditional spin-squeezing that is automatically generated by continuous observation of an atomic sample Takahashi et al. (1999); Kuzmich et al. (2000); Thomsen et al. (2002). This dynamically generated squeezing does not occur in steady-state (narrow-band) magnetometers because of the continuous optical pumping. However, it should be possible to enable sub-shotnoise magnetometry by turning off the optical pumping once a coherent spin state has been prepared followed by continuous observation of the atoms.

But the nature of conditional spin-squeezing gives rise to potential complications that make it initially unclear how to exploit the reduced uncertainty for improved magnetometry. Fig. 1 shows simulated data (generated according to a quantum trajectory model described below) of a spin ensemble under continuous measurement with no external field, $`B=0`$. As $`\mathrm{\Delta }\widehat{J}_\mathrm{z}^2`$ decreases \[shaded region in Fig. 1(B)\] with the onset of spin-squeezing, there is no apparent change in the noise of the associated $`\widehat{J}_\mathrm{z}(t)`$ measurement, $`y(t)`$ \[Fig. 1(A)\], which is due to constant optical shotnoise.

The dynamical generation of spin-squeezing starting from an initial coherent state involves a stochastic transient at early times. As suggested by the error-ellipse diagrams of Fig. 1(B), conditional evolution gradually localizes the quantum spin state around a constant, but random, value of $`\widehat{J}_\mathrm{z}_\mathrm{c}`$. In an ensemble of continuous measurement trajectories, this constant value would be distributed with a variance of $`J/2`$ corresponding to $`\mathrm{\Delta }\widehat{J}_\mathrm{z}^2`$ of the initial coherent state. Therefore, the mean value of $`\widehat{J}_\mathrm{z}_\mathrm{c}`$ assumes a non-zero value even in the absence of an applied magnetic field, producing a stochastic offset in the photocurrent that must be distinguished from Larmor precession in a magnetometry experiment.

Fortunately, with appropriate filtering, Larmor precession of the spin state can be distinguished from the projection noise in such a way that the field estimation benefits from spin-squeezing. In this Letter, we demonstrate that quantum trajectory theory Carmichael (1993); Wiseman and Milburn (1994) allows one to construct a Kalman filter Belavkin (1999); Mabuchi (1996); Verstraete et al. (2001); Stockton et al. (2003) that optimally estimates the field magnitude from continuously observed conditional atomic dynamics. This filtering procedure enables Heisenberg limited magnetometry despite the optical shotnoise and the transient effects of spin state estimation. Furthermore, we show that for time-invariant fields, our optimal strategy approximately reduces to the simple and intuitive data analysis procedure of linear regression which is a potentially simpler experimental approach to sub-shotnoise magnetometry.

We propose a magnetometer in which the atomic ensemble undergoes a continuous quantum non-demolition (QND) observation of $`\widehat{J}_\mathrm{z}`$. It has been shown that such a measurement can be implemented by detecting $`\widehat{J}_\mathrm{z}`$-dependent changes in the phase of an off-resonant cavity mode coupled to the atomic ensemble Thomsen et al. (2002) or by the Faraday rotation of a far-detuned travelling mode Silberfarb and Deutsch (2003); Smith et al. (2003) that passes through the ensemble. In both cases the magnetometer photocurrent is given by,

$$y(t)dt=2\eta \sqrt{M}\widehat{J}_\mathrm{z}(t)_\mathrm{c}dt+\sqrt{\eta }dW(t)$$
(2)
where $`\widehat{J}_\mathrm{z}(t)_c`$ is the conditional expectation value of $`\widehat{J}_\mathrm{z}`$, $`\eta `$ is the detector efficiency and $`M`$ (in units of frequency) is an implementation-dependent constant referred to as the measurement strength. The optical shotnoise is reflected by stochastic increments, $`dW(t)`$, that obey Gaussian white-noise statistics, $`E[dW]=0`$ and $`dW^2=dt`$.

Conditional evolution of the atomic ensemble subjected to a magnetic field along the $`y`$-axis and a QND measurement of $`\widehat{J}_\mathrm{z}`$ is described by the stochastic master equation,

$`d\widehat{\rho }_c(t)`$ $`=`$ $`i\gamma B[\widehat{J}_\mathrm{z},\widehat{\rho }_c]dt+M𝒟[\widehat{J}_\mathrm{z}]\widehat{\rho }_cdt`$
$`+\sqrt{M\eta }[\widehat{J}_\mathrm{z}]\widehat{\rho }_cdW(t)`$
where $`\widehat{\rho }_c(t)`$ is the reduced atomic density operator conditioned on the measurement record Wiseman and Milburn (1994). The super-operators, $`𝒟`$ and $``$, are given by, $`𝒟[r]\widehat{\rho }=r\widehat{\rho }r^{}(r^{}r\widehat{\rho }+\widehat{\rho }r^{}r)/2`$ and $`[r]\widehat{\rho }=r\widehat{\rho }+\widehat{\rho }r^{}\mathrm{tr}[(r+r^{})\widehat{\rho }]\widehat{\rho }`$, and the initial condition is an optically-pumped coherent spin state along the $`x`$-axis, $`\widehat{\rho }(0)=|J_xJ|`$.

Each term in Eq. (Quantum Kalman Filtering and the Heisenberg Limit in Atomic Magnetometry) has a physical implication for magnetometry. First, the Hamiltonian, $`H(B)=\gamma B\widehat{J}_y`$, generates the desired Larmor precession signal used to detect the magnetic field. The second term reflects measurement induced atomic decoherence that results from coupling the ensemble to the optical shotnoise on the probe laser. As a result the length of the Bloch vector decays over time, $`J(t)=J\mathrm{exp}(Mt/2)`$, and $`M`$ can be related to a bound on the transverse spin relaxation, $`T_22M^1`$.

The significance of the third term in Eq. (Quantum Kalman Filtering and the Heisenberg Limit in Atomic Magnetometry) is best seen after employing an approximation that holds for a large net magnetization and small field ($`\omega _Lt1`$). For an ensemble polarized along the $`x`$-axis, quantum fluctuations in $`\widehat{J}_x`$ are at least second order and the operator, $`\widehat{J}_x`$, is well-approximated by the length of the Bloch vector, $`\widehat{J}_xJ`$. Physically, this assumption capitalizes on the large value of $`J`$ to treat the Bloch sphere as a locally flat phase space. This approximation is extremely good for both coherent and squeezed states with $`J1`$. In a Gaussian approximation, the first and second moments of $`\widehat{J}_\mathrm{z}`$ are sufficient to completely characterize the atomic state. Therefore, equations of motion for the mean and variance,

$`d\widehat{J}_\mathrm{z}_c`$ $`=`$ $`\gamma BJe^{Mt/2}dt+2\sqrt{M\eta }\mathrm{\Delta }\widehat{J}_\mathrm{z}^2dW(t)`$ (4)
$`d\mathrm{\Delta }\widehat{J}_\mathrm{z}^2`$ $`=`$ $`4M\eta \mathrm{\Delta }\widehat{J}_\mathrm{z}^2^2dt`$ (5)
provide a closed representation of the magnetometer’s conditional quantum dynamics (in the $`\omega _Lt1`$ and $`tM^1`$ limits). The physical significance of Eq. (4) is that the atomic Bloch vector experiences two types of motion: deterministic Larmor precession and stochastic diffusion. Equation (5) reflects the deterministic reduction of $`\mathrm{\Delta }\widehat{J}_\mathrm{z}^2`$ as the atomic state is localized by the observation process, i.e., conditional spin-squeezing.

Equations (4) and (5) can be used to implement an optimal estimation procedure that capitalizes on squeezing without mistaking measurement-induced Bloch vector rotations for true Larmor precession. Since the atomic dynamics are stochastic, the estimator must be described probabilistically—we desire a conditional probability distribution, $`p\left(B|\mathrm{\Xi }_{[0,t]}\right)`$, which measures the likelihood that the field has magnitude $`B`$ given the measurement record, $`\mathrm{\Xi }`$, defined in terms of the photocurrent, $`d\mathrm{\Xi }_ty(t)dt/(2\eta \sqrt{M})`$. The estimated magnitude, $`\stackrel{~}{B}`$, and its uncertainty, $`\mathrm{\Delta }\stackrel{~}{B}^2`$ are obtained from the moments,

$`\stackrel{~}{B}`$ $`=`$ $`{\displaystyle Bp\left(B|\mathrm{\Xi }_{[0,t]}\right)𝑑B}`$ (6)
$`\mathrm{\Delta }\stackrel{~}{B}^2`$ $`=`$ $`{\displaystyle (B\stackrel{~}{B})^2p\left(B|\mathrm{\Xi }_{[0,t]}\right)𝑑B}`$ (7)
of the conditional distribution, $`p\left(B|\mathrm{\Xi }_{[0,t]}\right)`$.

Constructing a maximum-likelihood estimator is accomplished by defining an update rule that iteratively improves $`p\left(B|\mathrm{\Xi }_{[0,t]}\right)`$ as the measurement record is acquired. Prior knowledge of the distribution of magnetic field values is encoded in $`p\left(B|\mathrm{\Xi }_0\right)`$, which may be assigned infinite variance in order to assure an unbiased estimate. Optimality requires that the conditional probability must be updated according to a Bayes’ rule,

$$dp\left(B|d\mathrm{\Xi }_t,\mathrm{\Xi }_{[0,t)}\right)=dq\left(d\mathrm{\Xi }_t|B,\mathrm{\Xi }_{[0,t)}\right)p\left(B|\mathrm{\Xi }_{[0,t)}\right)$$
(8)
where $`dq\left(d\mathrm{\Xi }_t|B,\mathrm{\Xi }_{[0,t)}\right)`$ is an infinitesimal conditional probability that describes the likelihood of the evolving measurement record, $`d\mathrm{\Xi }_t`$, given a field with magnitude $`B`$ and past history, $`\mathrm{\Xi }_{[0,t)}`$. The utility of Bayes’ rule is that $`q\left(d\mathrm{\Xi }_t|B,\mathrm{\Xi }_{[0,t)}\right)`$ can be computed using quantum trajectory theory, Eqs. (4) and (5).

Implementing this parameter estimator is best accomplished by a (recursive) Kalman filter Belavkin (1999); Mabuchi (1996); Verstraete et al. (2001). It can be shown that the filtering equations,

$$d\stackrel{~}{x}=𝐀\stackrel{~}{x}dt+D^2(𝐁+\mathrm{𝐕𝐂}^T)(d\mathrm{\Xi }𝐂\stackrel{~}{x}dt)$$
(9)
with $`\stackrel{~}{x}\left(\begin{array}{cc}\stackrel{~}{J}_\mathrm{z}& \stackrel{~}{B}\end{array}\right)^T`$, ($`\stackrel{~}{J}_\mathrm{z}`$ is the estimate of $`\widehat{J}_\mathrm{z}_c`$),

$$𝐀=\gamma Je^{\frac{Mt}{2}}\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),𝐁=\left(\begin{array}{c}\mathrm{\Delta }\widehat{J}_\mathrm{z}^2\\ 0\end{array}\right),𝐂=\left(\begin{array}{cc}1& 0\end{array}\right),$$

$`D=1/(2\sqrt{M\eta })`$, and $`\stackrel{~}{x}(0)=0`$ implement Eq. (8). We note that it is possible to extend the Kalman filter to account for time-varying or stochastic fields Stockton et al. (2003) as well as to implement quantum feedback control Stockton et al. (2003); Thomsen et al. (2002).

The conditional quantum dynamics, particularly spin-squeezing and the exponential decay of the Bloch vector, enter the estimation process via the covariance matrix,

$$𝐕(t)=\left(\begin{array}{cc}\mathrm{\Delta }\stackrel{~}{J}_\mathrm{z}^2(t)& \mathrm{\Delta }(\stackrel{~}{J}_\mathrm{z}\stackrel{~}{B})(t)\\ \mathrm{\Delta }(\stackrel{~}{B}\stackrel{~}{J}_\mathrm{z})(t)& \mathrm{\Delta }\stackrel{~}{B}^2(t)\end{array}\right)$$
(10)
which describes the uncertainty in the parameter estimations of $`\widehat{J}_\mathrm{z}`$ and $`B`$. $`𝐕(t)`$ evolves deterministically according to the matrix Riccati equation,

$`\dot{𝐕}`$ $`=`$ $`\left(𝐀D^2\mathrm{𝐁𝐂}\right)𝐕+𝐕\left(𝐀D^2\mathrm{𝐁𝐂}\right)^T`$
$`+D^2𝐕\left(𝐂^T𝐂\right)𝐕`$
subject to the initial conditions $`\mathrm{\Delta }\stackrel{~}{J}_\mathrm{z}^2(0)=0`$ and $`\mathrm{\Delta }(\stackrel{~}{J}_\mathrm{z}\stackrel{~}{B})(0)=0`$, with $`\mathrm{\Delta }\stackrel{~}{B}^2(0)`$ chosen to reflect prior knowledge on the distribution of magnetic field values. Lacking any such knowledge, one can set $`\mathrm{\Delta }\stackrel{~}{B}^2(0)\mathrm{}`$.

The smallest detectable magnetic field as a function of $`J`$ and the measurement duration, $`t`$, is determined by the estimator variance, $`\mathrm{\Delta }\stackrel{~}{B}^2`$. Solving the the matrix Riccati equation \[which is analytically soluble for the Kalman filter in Eq. (9)\], provides the time-dependent magnetic field detection threshold, $`\delta \stackrel{~}{B}\sqrt{\mathrm{\Delta }\stackrel{~}{B}^2(t)}`$,

$$\delta \stackrel{~}{B}(t)=\frac{M}{4\gamma J}\sqrt{\frac{(1+2\eta JMt)}{ae^{Mt}+4e^{Mt/2}(4\eta J+1)+b}}$$
(12)
with $`a`$ and $`b`$ given by,

$`a`$ $`=`$ $`(2\eta J(Mt+4)+1)`$
$`b`$ $`=`$ $`Mt+2\eta J(Mt4)3.`$
Expanding Eq. (12) to leading order in $`t`$ provides an expression for the detection threshold,

$$\delta \stackrel{~}{B}(t)\frac{1}{\gamma J}\sqrt{\frac{3}{M\eta t^3}},t(JM)^1.$$
(13)
that is directly comparable to Eq. (1) when the measurement strength is chosen to be $`MT_2^1`$ such that maximal spin-squeezing is achieved at time $`t=M^1`$. Such a choice for $`M`$ permits a superior $`1/J`$ (equivalently $`1/N`$) scaling that is characteristic of the Heisenberg squeezing limit Kitagawa and Ueda (1993); Thomsen et al. (2002). The optical shotnoise \[of order unity in this model, see Eq. (2)\] enters implicitly through the signal to noise ratio, SNR=$`J\sqrt{M}`$, which highlights the utility of the Kalman filter as a whitening filter— it extracts the non-stationary spin-squeezing and Larmor precession dynamics despite the presence of Gaussian noise.

Figure 2 shows numerical results that demonstrate the performance of our quantum Kalman filter (QKF). The simulations were performed for an atomic ensemble with $`J=4\times 10^6`$, $`\gamma =1`$ kHz/mG, $`\eta =1`$, and $`M=100`$ kHz in a background magnetic field of $`B=1`$ $`\mu `$G. These values nearly correspond to a magnetometer constructed from $`N=10^6`$ ground state Cs atoms coupled to a high-finesse optical cavity with single photon Rabi frequency $`g_0`$ = 10 MHz and decay rate, $`\kappa `$ = 1 MHz. The QND measurement corresponds to a phase-quadrature homodyne detection of the transmitted cavity light with a cavity mode ($`P=100`$ $`\mu `$W) that is blue detuned by $`\mathrm{\Delta }=1`$ GHz from the Cs transition at $`\lambda 852`$ nm. The initial estimator variance was chosen to be $`\mathrm{\Delta }\stackrel{~}{B}^2(0)=100`$ $`\mu `$G2 which is the initial value that one would select given prior knowledge that the magnetic field could be treated as a Gaussian random variable with this variance. For the parameters we selected, Larmor precession and spin projection noise have comparable magnitudes on timescales of order $`(JM)^1`$.

The Kalman estimation error (crosses in Fig. 2) was computed from the ensemble average, $`E[(\stackrel{~}{B}_iB)^2]`$, for $`10^5`$ trajectories, and the solid line shows the estimation uncertainty $`\mathrm{\Delta }\stackrel{~}{B}(t)`$ obtained by integrating Eq. (Quantum Kalman Filtering and the Heisenberg Limit in Atomic Magnetometry). Since our simulations were performed with $`B=1`$ $`\mu `$G, the empirical performance of the QKF closely matches the solid line. The dotted line in Fig. 2 shows the analytic Riccati solution, given by Eq. (12), for $`\mathrm{\Delta }\stackrel{~}{B}^2(0)\mathrm{}`$, which would be the expected QKF performance in a scenario with no prior knowledge of the magnetic field. Fig. 2 also shows the estimation error for simple linear regression of the measurement record (open circles). Assuming that the Bloch vector has not decayed significantly, $`\stackrel{~}{B}`$ is proportional to the slope of a line fit to the (filtered) photocurrent, $`y(t)`$. The estimation error for linear regression was obtained by computing $`E[(\stackrel{~}{B}_iB)^2]`$ for $`10^5`$ trajectories. Although the QKF is clearly superior for short times, the two estimation procedures converge for $`(JM)^1t<M^1`$ and provide a quantum parameter estimation with $`\delta \stackrel{~}{B}`$ 0.01 nG in $`t1`$ ms. The inset of Fig. 2 highlights the $`\delta \stackrel{~}{B}J^1`$ scaling that distinguishes both the QKF (crosses) and regression (open circle) estimators from the conventional shot-noise limit, Eq. 1. For sufficiently large times both estimators achieve the detection threshold in Eq. 13 (solid lines).

The QKF and linear regression differ mainly in how they treat the initial diffusive transient of $`\widehat{J}_\mathrm{z}`$ \[Fig. 1(B)\]. Since the QKF is derived from a quantum trajectory model, it is aware of the short-time diffusion and strategically under-weights the photocurrent at early times \[via the Kalman gain, $`𝐆D^2(𝐁+\mathrm{𝐕𝐂}^T)`$\]. At late times the regression analysis manages to absorb the initial diffusive transient into the $`y`$-intercept of the linear fit. Although $`\widehat{J}_\mathrm{z}(0)=0`$, the $`\widehat{J}_\mathrm{z}`$ localization process gradually determines an effective offset in the photocurrent during the interval $`t(JM)^1`$. Without explicit knowledge of the conditional dynamics the linear regression equally weights the photocurrent for $`t<(JM)^1`$. This decreases the quality of the fit, but the resulting error becomes insignificant for $`t(JM)^1`$.

Our analysis suggests that estimation procedures based on conditional quantum dynamics can play a crucial role in optimizing both the sensitivity and the bandwidth in atomic magnetometry. While conventional steady-state magnetometers can only improve their detection capabilities by increasing the number of atoms or the averaging time, the quantum estimator can achieve greater precision for the same value of $`t`$ and $`N`$ by improving the measurement strength. The significance of the Kalman filter is the optimality that is guaranteed by its derivation from a Bayes’ rule. Our finding that linear regression closely approximates the optimal procedure indicates a potentially simpler experimental procedure for sub-shotnoise magnetometry. Although Heisenberg limited spin-squeezing should be possible using current techniques in cavity quantum electrodynamics (a discussion is provided in Thomsen et al. (2002)), the experimental difficulty of achieving this limit makes it desirable to have an optimal estimator such as the QKF to fully exploit even a small amount of squeezing, to treat fluctuating fields and to achieve estimator robustness Stockton et al. (2003). In either case, estimation procedures that allow and account for conditional quantum dynamics— whether explicitly as in the QKF or implicitly as in linear regression— offer substantial improvement over steady-state procedures.

This work was supported by the NSF (PHY-9987541, EIA-0086038), the ONR (N00014-00-1-0479), and the Caltech MURI Center for Quantum Networks (DAAD19-00- 1-0374). JKS acknowledges a Hertz fellowship. Please visit http://minty.caltech.edu/Ensemble for simulation source code and detailed notes on the Kalman filter derivation.


format: text

text:
# Quantum jumps and entropy production

## I Introduction

Dynamical processes in dissipative open systems involve, in general, irreversible transformations of non-equilibrium states GROOT ; KEIZER . An important physical quantity of the theory of non-equilibrium processes is the entropy production, that is the rate at which entropy is produced as a result of irreversible processes NICOLIS . In the classical theory the entropy production provides an important Lyapunov functional of non-equilibrium stationary states and leads to the principle of minimal entropy production. It further enables the determination of non-equilibrium transport properties and Onsager coefficients (for a recent discussion of the classical theory see, e. g. DAEMS ; BAG ).

In the quantum theory of open systems TheWork  the entropy production is usually related to the negative derivative of the relative entropy with respect to an invariant state. On the ground of the weak coupling master equation for the open system’s density matrix the corresponding expressions can again be shown to have useful convexity properties and to satisfy a minimal entropy production principle SPOHN , as well as a kind of quantum Onsager relations LENDI . The concept of relative entropy is also useful in quantum information theory NIELSEN , e. g. in the quantification of entanglement VEDRAL ; WERNER  and of the success of quantum teleportation BOSE . A general account of the theory of quantum entropy may be found in WERL .

During the past decade there has been considerable interest in the application of stochastic wave function methods to the analysis and simulation of dissipative open quantum systems. By contrast to the conventional description in terms of the reduced density matrix, in these methods one considers an ensemble of pure states whose time-evolution is governed by a stochastic process in the underlying Hilbert space $``$ of the open system. In a particularly interesting method of this kind, developed by several authors DALIBARD ; DUM ; CARMICHAEL , the wave function represents a piecewise deterministic process in Hilbert space. This means that the stochastic evolution consists of smooth, deterministic parts which are interrupted by instantaneous quantum jumps. Recently, it has been shown that the stochastic wave function method can be generalized to enable the treatment of non-Markovian quantum processes BKP  and of Bosonic CARUSOTTO  and of Fermionic CHOMAZ  many-body systems by means of an appropriate stochastic representation in the doubled Hilbert space $``$.

The aim of the present paper is to connect the stochastic state vector dynamics to the production of quantum entropy. To this end, we apply the entropy balance equation of non-equilibrium thermodynamics to the von Neumann entropy of the open system. A quantum jump of the stochastic dynamics is linked to the exchange of an entropy quantum between system and environment. This enables one to express the entropy flux in terms of the random processes which count the numbers of quantum jumps belonging to the various decay modes of the open system. In this way, a stochastic representation of the production rate of quantum entropy is derived from the piecewise deterministic state vector evolution. It will be demonstrated that the obtained expression for the entropy production provides a non-negative and convex functional on the state space of the open system. Furthermore, the stochastic representation is shown to lead to a Monte Carlo simulation method for the determination of the quantum entropy production in open systems.

The paper is organized as follows. In Sec. II we briefly recall the dynamical description of open quantum systems in terms of quantum master equations and of stochastic differential equations. Sec. III contains the construction of the stochastic representation of the entropy production rate as well as the proof of its positivity and convexity. Finally, a discussion of the results and examples of Monte Carlo simulations of the entropy production are given in Sec. IV.

## II Quantum master equations and piecewise deterministic processes

We consider an open quantum system with reduced density matrix $`\rho =\rho (t)`$. The open system is coupled to a heat bath, that is to an external reservoir in a thermal equilibrium state at temperature $`T`$. The free evolution of the system is described by some Hamiltonian $`H_0`$. Furthermore, the system may be subjected to an external time-dependent perturbation, e. g. an external radiation field represented by a Hamiltonian $`H_p(t)`$. Invoking the Markov approximation, the motion of the reduced system is assumed to be representable in terms of a quantum master equation with a generator $`_t`$ in Lindblad form DAVIES ; GORINI ; LINDBLAD76  (choosing units such that $`\mathrm{}=1`$),

$$\frac{d}{dt}\rho =_t(\rho )=i[H_0+H_p(t),\rho ]+𝒟(\rho ).$$
(1)
Since we allow the perturbation to be explicitly time-dependent the dynamics does not, in general, provide a dynamical semigroup. However, since $`_t`$ is in Lindblad form for each fixed $`t`$, the time-evolution defined by Eq. (1) yields a two-parameter family of quantum dynamical maps ALICKI .

The commutator in Eq. (1) represents the unitary part of the motion of the reduced system, while $`𝒟(\rho )`$ is a super-operator known as the dissipator of the master equation,

$`𝒟(\rho )`$ $`=`$ $`{\displaystyle \underset{i}{}}\gamma _i^{}\left(A_i^{}\rho A_i^+{\displaystyle \frac{1}{2}}\{A_i^+A_i^{},\rho \}\right)`$ (2)
$`+{\displaystyle \underset{i}{}}\gamma _i^+\left(A_i^+\rho A_i^{}{\displaystyle \frac{1}{2}}\{A_i^{}A_i^+,\rho \}\right).`$
The Lindblad operators $`A_i^\pm `$ provide the coupling to the various decay channels of the open system with corresponding damping rates $`\gamma _i^\pm `$. We shall assume in the following that the master equation (1) has been derived from an underlying microscopic theory through the weak-coupling approximation. This implies that the Lindblad operators are obtained as eigen-operators of the unperturbed system Hamiltonian,

$$[H_0,A_i^\pm ]=\pm \omega _iA_i^\pm ,$$
(3)
and satisfy $`A_i^+=(A_i^{})^{}`$. The relation (3) also implies that

$$[H_0,A_i^+A_i^{}]=[H_0,A_i^{}A_i^+]=0.$$
(4)
The damping rates $`\gamma _i^\pm `$ of the decay modes are determined through certain reservoir correlation functions satisfying the Kubo-Martin-Schwinger condition which leads to the relation

$$\gamma _i^+=\mathrm{exp}\left(\omega _i/T\right)\gamma _i^{}.$$
(5)
As mentioned in the introduction the dynamics provided by the master equation (1) can be represented in terms of a piecewise deterministic process $`\psi (t)`$ (PDP) in the Hilbert space of the open system in the sense that its covariance matrix

$$\rho (t)\mathrm{E}\left[|\psi (t)\psi (t)|\right]=D\psi D\psi ^{}P[\psi ,t]|\psi \psi |$$
(6)
satisfies the master equation (1). Here, $`\mathrm{E}`$ denotes the expectation value of the stochastic process $`\psi (t)`$ which is defined in terms of the corresponding time-dependent probability density functional $`P[\psi ,t]`$.

An appropriate stochastic differential equation for the state vector dynamics is given by

$`d\psi (t)`$ $`=`$ $`iG_t(\psi (t))dt`$ (7)
$`+{\displaystyle \underset{i}{}}\left({\displaystyle \frac{A_i^{}\psi (t)}{A_i^{}\psi (t)}}\psi (t)\right)dN_i^{}(t)`$
$`+{\displaystyle \underset{i}{}}\left({\displaystyle \frac{A_i^+\psi (t)}{A_i^+\psi (t)}}\psi (t)\right)dN_i^+(t).`$
The first term on the right-hand side provides the deterministic evolution periods of the process $`\psi (t)`$. It corresponds to the deterministic non-linear Schrödinger-type equation $`\dot{\psi }=iG_t(\psi )`$, where

$$G_t(\psi )=\widehat{H}\psi +\frac{i}{2}\underset{i}{}\left(\gamma _i^{}A_i^{}\psi ^2+\gamma _i^+A_i^+\psi ^2\right)\psi $$
(8)
is a non-linear operator with a non-Hermitian Hamiltonian $`\widehat{H}=\widehat{H}(t)`$ given by

$$\widehat{H}=H_0+H_p(t)\frac{i}{2}\underset{i}{}\left(\gamma _i^{}A_i^+A_i^{}+\gamma _i^+A_i^{}A_i^+\right).$$
(9)
The quantities $`dN_i^\pm (t)`$ in Eq. (7) are the increments of integer-valued random processes $`N_i^\pm (t)`$ satisfying

$`dN_i^k(t)dN_j^l(t)`$ $`=`$ $`\delta _{ij}\delta _{kl}dN_j^l(t),`$ (10)
$`\mathrm{E}\left[dN_i^\pm (t)\right]`$ $`=`$ $`\gamma _i^\pm A_i^\pm \psi (t)^2dt,`$ (11)
where $`k,l=\pm `$. According to relation (10) these increments take on the values $`0`$ or $`1`$. Moreover, if $`dN_i^k(t)`$ is equal to $`1`$ for a particular pair $`(i,k)`$, all other increments are equal to zero. Thus, in view of the stochastic differential equation (7) the state vector then undergoes an instantaneous quantum jump of the form

$$\psi (t)\frac{A_i^k\psi (t)}{A_i^k\psi (t)}.$$
(12)
The relation (11) states that the expectation values of the increments $`dN_i^k(t)`$ increase linearly with the time increment $`dt`$ and are directly proportional to the damping rates $`\gamma _i^k`$ and to the square of the norm of $`A_i^k\psi (t)`$. Since the latter depends on time through the deterministic motion, the processes $`N_i^k(t)`$ are, in general, inhomogeneous in time. They count the number of quantum jumps of type $`(i,k)`$, i.e. the number of the jumps (12) with Lindblad operator $`A_i^k`$.

## III Definition of the entropy production rate

### III.1 Entropy flux and numbers of quantum jumps

An appropriate starting point for the definition of the entropy production rate $`\sigma `$ in an open quantum system is provided by the entropy balance equation,

$$\sigma =\frac{dS}{dt}+J.$$
(13)
Here, $`S`$ is the von Neumann entropy of the open system (the Boltzmann constant is set equal to $`1`$),

$$S=\mathrm{tr}\left\{\rho \mathrm{ln}\rho \right\}.$$
(14)
The quantity $`J`$ represents the entropy flux, that is the entropy which flows per unit of time from the open system to the reservoir, in which case we have $`J>0`$, or from the reservoir to the open system, such that $`J<0`$. Thus, according to Eq. (13) the quantity $`\sigma `$ is the rate at which entropy is produced owing to irreversible processes.

In order to define the entropy production rate $`\sigma `$ we therefore need an expression for the entropy flux $`J`$. To this end, we recall that the external reservoir causes the wave function to perform instantaneous quantum jumps given by (12). According to Eq. (3) the operator $`A_i^\pm `$ induces a transition of the system state in which the energy changes by the amount $`\pm \omega _i`$. Thus, $`A_i^+`$ describes upward and $`A_i^{}`$ downward transitions of the system. Any quantum jump with jump operator $`A_i^+`$ ($`A_i^{}`$) corresponds to the absorption (emission) of a quantum $`\omega _i`$ of energy from (into) the reservoir. Because the reservoir is assumed to be in a thermal equilibrium state at temperature $`T`$, the entropy exchange $`\delta S`$ associated with a quantum jump can be defined by

$$\delta S=\frac{\omega _i}{T}.$$
(15)
To justify this expression one first observes that within the Markov approximation of the reduced system dynamics the reservoir correlation functions decay on a time scale which is - by assumption - short compared to the relaxation time of the open system. It is therefore reasonable to assume that the energy exchange in a quantum jump with operator $`A_i^\pm `$ appears as a change $`\omega _i`$ of the heat energy of the reservoir. We further note that the reversible character of the energy exchange through a quantum jump derives from the fact that for each jump with operator $`A_i^{}`$ (emission) there also exists the reverse jump with operator $`A_i^+`$ (absorption).

Thus we see that in each quantum jump with operator $`A_i^\pm `$ an entropy quantum (15) is absorbed from or emitted into the reservoir. The integer-valued process $`N_i^\pm (t)`$ represents the number of quantum jumps with operator $`A_i^\pm `$ within the time interval $`[0,t]`$. It follows that in the time interval $`[0,t]`$, on the average, the net amount

$$\underset{i}{}\frac{\omega _i}{T}\mathrm{E}\left[N_i^{}(t)N_i^+(t)\right]$$
(16)
of entropy flows into the reservoir. Consequently, the mean entropy flux is obtained from the expression

$$J=\underset{i}{}\frac{\omega _i}{T}\frac{d}{dt}\mathrm{E}\left[N_i^{}(t)N_i^+(t)\right].$$
(17)
This relation finally leads to the following stochastic expression of the entropy production rate of the open system,

$$\sigma =\frac{dS}{dt}+\underset{i}{}\frac{\omega _i}{T}\frac{d}{dt}\mathrm{E}\left[N_i^{}(t)N_i^+(t)\right].$$
(18)
By use of this expression the entropy production rate $`\sigma (t)`$ can easily be determined by a Monte Carlo simulation of the stochastic process (7). To this end, one evaluates the von Neumann entropy by means of Eq. (6) and records the numbers $`N_i^\pm (t)`$ of the various quantum jumps carried out during the time interval $`[0,t]`$.

### III.2 Positivity of the entropy production

We note the following important properties of the entropy production rate $`\sigma `$ defined by expression (18). First, it is immediately clear that $`\sigma `$ is an extensive quantity: When combining statistically independent systems which are coupled to the same reservoir, the additivity of the von Neumann entropy and the additivity of the numbers of quantum jumps is obvious. Second, Eq. (18) implies that $`\sigma `$ does not depend on the specific representation of the state space of the reduced system, since neither the von Neumann entropy nor the number of quantum jumps of a given type changes under unitary transformations. In particular, the expression (18) is invariant under a transformation to the interaction picture with respect to the system Hamiltonian $`H_0`$. Third, the entropy production $`\sigma `$ as defined in this way is always non-negative,

$$\sigma 0,$$
(19)
which corresponds to the second law of thermodynamics.

To prove the inequality (19) we first write the entropy production rate as a functional $`\sigma =\sigma [\rho ]`$ of the reduced density matrix $`\rho `$. Invoking the master equation (1) one finds the time-derivative of the von Neumann entropy,

$$\frac{dS}{dt}=\mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}\rho \right\}.$$
(20)
To obtain an appropriate expression for the entropy flux we use Eqs. (17) and (11), as well as the representation (6) of the density matrix to get

$$J=\underset{i}{}\frac{\omega _i}{T}\mathrm{tr}\left\{\gamma _i^{}A_i^+A_i^{}\rho \gamma _i^+A_i^{}A_i^+\rho \right\}.$$
(21)
On the other hand, employing definition (2) of the dissipator and the commutation relations (3) and (4) one finds

$`\mathrm{tr}\left\{H_0𝒟(\rho )\right\}`$ (22)
$`=`$ $`{\displaystyle \underset{i}{}}\gamma _i^{}\mathrm{tr}\left\{A_i^+[H_0,A_i^{}]\rho +{\displaystyle \frac{1}{2}}[A_i^+A_i^{},H_0]\rho \right\}`$
$`+{\displaystyle \underset{i}{}}\gamma _i^+\mathrm{tr}\left\{A_i^{}[H_0,A_i^+]\rho +{\displaystyle \frac{1}{2}}[A_i^{}A_i^+,H_0]\rho \right\}`$
$`=`$ $`{\displaystyle \underset{i}{}}\omega _i\mathrm{tr}\left\{\gamma _i^{}A_i^+A_i^{}\rho \gamma _i^+A_i^{}A_i^+\rho \right\}`$
$`=`$ $`TJ.`$
The commutation relations (3) also imply that the Gibbs state $`\rho _{\mathrm{th}}=\mathrm{exp}(H_0/T)/Z`$ with partition function $`Z`$ represents a zero-mode of the dissipator of the master equation, that is

$$𝒟(\rho _{\mathrm{th}})=0.$$
(23)
Physically, this means that the Gibbs state is a stationary solution of the master equation without perturbation ($`H_p=0`$). We express the Hamiltonian $`H_0`$ in terms of the Gibbs state,

$$H_0/T=\mathrm{ln}\rho _{\mathrm{th}}+\mathrm{ln}Z,$$
(24)
which, together with Eq. (22), enables us to write the entropy flux as follows,

$$J=\mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}\rho _{\mathrm{th}}\right\}.$$
(25)
Combining this with Eq. (20) we are finally led to the following expression for the entropy production rate in the state $`\rho `$,

$$\sigma [\rho ]=\mathrm{tr}\left\{𝒟(\rho )\left(\mathrm{ln}\rho \mathrm{ln}\rho _{\mathrm{th}}\right)\right\}.$$
(26)
The inequality (19) can now be shown in two alternative ways. One way is to apply a theorem by Lindblad LINDBLAD75  to the dynamical semigroup $`\mathrm{\Lambda }_t=\mathrm{exp}(𝒟t)`$ whose generator is identical to the dissipator $`𝒟(\rho )`$ of the master equation, and to relate the functional (26) to the time-derivative of the relative entropy with respect to the Gibbs state $`\rho _{\mathrm{th}}`$. Another possibility is to use Lieb’s theorem LIEB  to conclude that the map $`\rho \mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}\rho \right\}`$ is a convex functional. Since $`\mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}\rho _{\mathrm{th}}\right\}`$ is linear in $`\rho `$, it follows that the entropy production rate $`\sigma =\sigma [\rho ]`$ represents a convex functional of the density matrix. With the help of this property it is easy to demonstrate that the entropy production is non-negative SPOHN . In fact, applying the definition of convexity to the functional $`\rho \mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}\rho \right\}`$ one is led to the following inequality which holds for all $`\lambda [0,1]`$,

$`\mathrm{tr}\left\{𝒟(\lambda \rho +[1\lambda ]\rho _{\mathrm{th}})\mathrm{ln}(\lambda \rho +[1\lambda ]\rho _{\mathrm{th}})\right\}`$
$`\lambda \mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}(\rho )\right\}[1\lambda ]\mathrm{tr}\left\{𝒟(\rho _{\mathrm{th}})\mathrm{ln}(\rho _{\mathrm{th}})\right\}.`$
Employing the linearity of $`𝒟(\rho )`$ and Eq. (23) we get

$$\lambda \mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}(\lambda \rho +[1\lambda ]\rho _{\mathrm{th}})\right\}\lambda \mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}(\rho )\right\}.$$

We divide by $`\lambda `$ and perform the limit $`\lambda 0`$ to arrive at

$$\mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}(\rho _{\mathrm{th}})\right\}\mathrm{tr}\left\{𝒟(\rho )\mathrm{ln}(\rho )\right\},$$

which, by virtue of (26), is equivalent to inequality (19).

## IV Discussion and examples

### IV.1 General properties of the entropy production

It has been shown in the preceding section that the piecewise deterministic process (7) for the stochastic state vector of an open quantum system yields a natural definition of the entropy production as a measure of the degree of the irreversibility of the dynamics. In this definition the entropy flux has been expressed in terms of the random numbers of the various types of quantum jumps corresponding to the decay channels of the open system.

The given proof of the positivity of the entropy production rate demonstrates not only that the entropy production $`\sigma `$ is non-negative (see Eq. (19)) but also that it is a convex functional and that $`\sigma (\rho _{\mathrm{th}})=0`$ is an absolute minimum of the entropy production (see Eq. (26)). Hence, in the absence of external fields ($`H_p=0`$) the entropy production vanishes when the equilibrium state is reached, in accordance with the second law of thermodynamics.

We remark that for vanishing external fields $`\sigma `$ coincides with the negative time derivative of the relative entropy with respect to the Gibbs state $`\rho _{\mathrm{th}}`$. This is no longer true when external fields are present. Contrary to the definition of entropy production in terms of the relative entropy, the above definition does not assume the existence of a stationary density matrix.

In the case of a non-vanishing external field ($`H_p0`$) one expects, in general, that the systems stays away from equilibrium and that, therefore, entropy is continuously produced. An important physical situation arises (see the examples below) if there exists a stationary solution $`\rho ^\mathrm{s}`$ of the master equation in the interaction picture. Since the expression for the entropy production is not affected by the transformation to the interaction picture one finds that the entropy production in the stationary state is given by the constant value

$$\sigma ^\mathrm{s}=\underset{i}{}\frac{\omega _i}{T}\frac{d}{dt}\mathrm{E}\left[N_i^{}(t)N_i^+(t)\right],$$
(27)
where $`\mathrm{E}`$ denotes the expectation value in the stationary state. In such a stationary non-equilibrium state entropy is thus produced at a constant rate $`\sigma ^\mathrm{s}`$. What happens physically is that the energy supplied to the system by the external field is ultimately transferred as heat energy to the reservoir leading to a constant entropy production rate.

It is well-known that the stochastic representation of the master equation (1) through a stochastic process of the form (7) is not unique: In general, one can choose, in an infinite number of different ways, the Lindblad jump operators and the Hamiltonian operator without changing the form of the Lindblad generator $`_t`$ of the master equation (1). However, different such choices do change, in general, the stochastic dynamics (7). This fact is well understood in the framework of continuous measurement theory according to which different choices for the jump operators correspond to different measurement schemes used to monitor the open system WISEMAN1 ; WISEMAN2 . Our construction of expression (18) for the entropy production rate presupposes, however, that the jump operators $`A_i^\pm `$ are eigen-operators of the system Hamiltonian, satisfying the commutation relation (3). It is this property which enables one to establish the connection between the stochastic expression (18) and the functional (26). Physically, it is the eigen-operator property which allows one to associate a certain quantum jump of type $`(i,k)`$ with the exchange of a quantum $`\omega _i`$ of energy with the thermal reservoir.

In the case of a non-degenerate frequency spectrum $`\{\omega _i\}`$ the requirement provided by the commutation relations (3) uniquely fixes the jump operators and, thus, the stochastic state vector dynamics. If a certain frequency $`\omega _i`$ is degenerate one has the freedom to replace the corresponding jump operators $`A_{i\alpha }^\pm `$ (the index $`\alpha `$ labels different Lindblad operators belonging to the same frequency $`\omega _i`$) by new operators

$$\stackrel{~}{A}_{i\alpha }^{}=\underset{\beta }{}u_{\alpha \beta }A_{i\beta }^{},\stackrel{~}{A}_{i\alpha }^+=\left(\stackrel{~}{A}_{i\alpha }^{}\right)^{},$$
(28)
where $`u_{\alpha \beta }`$ is a unitary matrix. However, in view of Eq. (21) it is immediately clear that expression (18) remains invariant under this replacement, showing the uniqueness of the stochastic expression of the entropy production rate.

### IV.2 Entropy production rate of a driven qubit

We conclude the paper by two examples which illustrate several of the features discussed above. The first example is two-state system GARDINER , a qubit consisting of an excited state $`|e`$ and a ground state $`|g`$ with corresponding energies $`E_e`$ and $`E_g`$. The qubit is coupled to a thermal reservoir and to a single-mode driving field on resonance with the transition frequency $`\omega =E_eE_g`$. The following discussion will be carried out in the interaction picture. The external field is treated in the rotating wave approximation such that the corresponding interaction Hamiltoninan $`H_p`$ becomes time-independent in the interaction picture.

The stochastic dynamics given by Eq. (7) is defined here through two jump operators, given by the raising and lowering operators

$$A^+=|eg|,A^{}=|ge|,$$
(29)
with stochastic processes $`N^\pm (t)`$ for the numbers of the corresponding jumps. The non-Hermitian Hamiltonian (9) therefore takes the form

$$\widehat{H}=\frac{\mathrm{\Omega }}{2}\left(A^++A^{}\right)\frac{i}{2}\left(\gamma ^{}A^+A^{}+\gamma ^+A^{}A^+\right),$$
(30)
where the Rabi frequency $`\mathrm{\Omega }`$ is assumed to be real. Equation (27) then yields the following entropy production rate in the stationary state,

$$\sigma ^\mathrm{s}=\frac{\omega }{T}\frac{\mathrm{\Omega }^2}{\gamma }S_3=\frac{\omega }{T}\frac{\gamma ^{}\gamma ^+}{2+(\gamma /\mathrm{\Omega })^2},$$
(31)
where $`\gamma =\gamma ^{}+\gamma ^+`$. Obviously, $`\sigma ^\mathrm{s}`$ is proportional to the inversion $`S_3=\rho _{ee}^\mathrm{s}\rho _{gg}^\mathrm{s}`$ in the stationary state.

Figure 1 shows the result of a Monte Carlo simulation of the PDP defined by Eq. (7) and compares it with the solution of the density matrix equation (1). In the stochastic simulation the entropy production rate has been determined with the help of Eq. (18), where both the von Neumann entropy $`S`$ and the expectation values of the numbers $`N^\pm (t)`$ of quantum jumps have been estimated from a sample of $`10^5`$ realizations of the PDP. The initial condition has been chosen to be the ground state $`|g`$ of the qubit. In the absence of the driving field the entropy production decreases to zero during the approach to the thermal equilibrium state. However, as can be seen from the figure, in the presence of an external driving the entropy production oscillates and converges to the stationary value given by Eq. (31).

### IV.3 $`\mathrm{\Lambda }`$-configuration with a dark state

As our second example we investigate a two-level atom with transition frequency $`\omega `$. Both levels are threefold degenerate forming manifolds with total angular momentum $`J_e=1`$ (excited level) and $`J_g=1`$ (ground level). We introduce the energy eigenstates $`|g,m_g`$ and $`|e,m_e`$ which are simultaneously eigenstates of the $`z`$-component of the atomic angular momentum operator with eigenvalues $`m_g`$ and $`m_e`$, respectively. In addition to the coupling to a thermal reservoir the atom is subjected to a resonant driving field which is linearly polarized in the $`y`$-direction. If we take some initial state in the manifold spanned by the states $`|g,m_g=\pm 1`$ we then find that the dynamics is confined to the subspace spanned by the basis states

$$|e,m_e=0,|g,m_g=+1,|g,m_g=1,$$
(32)
forming a level scheme known as $`\mathrm{\Lambda }`$-configuration SCULLY .

The dynamics of the system will again be studied in the interaction picture. It can be described by the four jump operators,

$$A_1^{}=\frac{1}{\sqrt{2}}|g,1e,0|,A_2^{}=\frac{1}{\sqrt{2}}|g,+1e,0|,$$
(33)
and $`A_1^+=(A_1^{})^{}`$, $`A_2^+=(A_2^{})^{}`$. The corresponding matrix representation of the non-Hermitian Hamiltonian (9) in the basis (32) reads

$$\widehat{H}=\frac{i}{2}\left(\begin{array}{ccc}\gamma ^{}& \mathrm{\Omega }& \mathrm{\Omega }\\ \mathrm{\Omega }& \frac{\gamma ^+}{2}& 0\\ \mathrm{\Omega }& 0& \frac{\gamma ^+}{2}\end{array}\right).$$
(34)
A numerical simulation of the dynamics is shown in Fig. 2. In the Monte Carlo simulation of the corresponding stochastic differential equation (7) the quantity $`\sigma (t)`$ has been determined with the help of Eq. (18) by recording the stochastic numbers $`N_{1,2}^\pm (t)`$ of upward and downward transitions in a sample of realizations of the process.

The stationary entropy production rate may be found from Eq. (27) with the result

$$\sigma ^\mathrm{s}=\frac{\omega }{T}\gamma ^+\rho _{23}^s=\frac{\omega }{T}\gamma ^+\frac{\gamma ^{}\gamma ^+}{2\gamma ^{}+4\gamma ^++\frac{\gamma ^+}{\mathrm{\Omega }^2}(\gamma ^{}+\frac{\gamma ^+}{2})^2}.$$
(35)
This equation shows that $`\sigma ^\mathrm{s}`$ is proportional to the matrix element $`\rho _{23}^s=g,+1|\rho ^\mathrm{s}|g,1`$ of the stationary density matrix $`\rho ^\mathrm{s}`$. We observe that, similarly to the case of the driven qubit (see Eq. (31)), $`\sigma ^\mathrm{s}`$ approaches a finite value as $`\mathrm{\Omega }\mathrm{}`$. For zero Rabi frequency, $`\mathrm{\Omega }=0`$, the stationary entropy production rate vanishes, in accordance with our general discussion. Moreover, the entropy production vanishes for infinite temperatures for both the driven qubit and the $`\mathrm{\Lambda }`$-configuration.

By contrast to the driven qubit, however, in the case of the $`\mathrm{\Lambda }`$-configuration the limit $`T0`$, corresponding to a reservoir at zero temperature, leads to a vanishing entropy production rate. In fact, in this limit we have $`\gamma ^{}\gamma _0`$, $`\rho _{23}^s1/2`$ and, therefore,

$$\underset{T0}{lim}\sigma ^\mathrm{s}=\underset{T0}{lim}\frac{\omega }{T}e^{\omega /T}\gamma _0\frac{1}{2}=0.$$
(36)
This behaviour can easily be understood if we note that the matrix element $`\rho _{23}`$ of the density matrix describes the coherence of the state

$$|\psi _\mathrm{d}=\frac{1}{\sqrt{2}}\left(|g,+1+|g,1\right).$$
(37)
This state is known as dark state since it does not couple to the driving field on account of a quantum interference effect.

According to Eq. (35) the stationary entropy production rate is thus proportional to the coherence of the dark state in the stationary state $`\rho ^\mathrm{s}`$. At zero temperature, any initial state is driven into the dark state under the dynamics. The reason of this fact is that $`|\psi _\mathrm{d}`$ represents a zero-mode of the non-Hermitian Hamiltonian (34) for $`T=0`$, that is for $`\gamma ^+=0`$. Equation (36) therefore states that the entropy production rate vanishes in the dark state. This follows immediately from the fact that $`|\psi _\mathrm{d}`$ belongs to the zero temperature ground state manifold.

Finally, we remark that this system provides an example of a system with a degenerate frequency spectrum, mentioned at the end of Sec. IV.1. Both jump operators given in (33) belong to the same transition frequency $`\omega `$. In a continuous measurement interpretation they correspond to a measurement of the angular momentum of the emitted quanta along the $`z`$-direction MOLMER . The measurement along any other direction amounts to a transformation of these jump operators which is of the form (28) and which does not change the expression for the entropy production rate. The latter is therefore seen to be independent of the measurement scheme.


text:
# 1 Cartoon picture of a phonon. Solid arrows represent excited atoms with momentum 𝐩 or -𝐩. Dashed arrows represent condensate atoms.

In a dilute, non-ideal Bose-Einstein (BE) condensate the spectrum of quasiparticle excitations is linear at low momentum: $`\stackrel{~}{E}(p)v_sp`$. Thus, long wavelength excitations are phonons — quantized sound waves — that travel at speed $`v_s`$. The standard Bogoliubov theory gives this result and it is also expected on general grounds since hydrodynamics applies to any medium at long length scales.

However, this familiar result has a puzzling aspect to it. Very low momentum atoms move very slowly while very low-momentum phonons move at a finite speed $`v_s`$. But how can a phonon in a dilute medium move so much faster than the atoms that make it up? Surprisingly, this issue seems not to be addressed in the literature. The purpose of this brief note is to offer an intuitive explanation.

It is important to stress that the puzzle arises because we are considering a dilute medium. In solids and liquids, where the atoms are “touching” (separations less than the interaction range), a tiny motion of one atom can be transmitted almost immediately to its neighbours, allowing compression and other waves to move at speeds much greater than the speed of the atoms themselves. However, in a dilute gas with only short-range interactions (and theoretically the interatomic potential can be replaced with a delta-function pseudopotential ) disturbances can only be communicated onwards by the atoms themselves travelling from one collision to the next. Thus, the disturbance should propagate no faster than the atoms’ own speeds. Indeed, for an ordinary, classical gas the speed of sound is of order the average speed (r.m.s. velocity) of the atoms for this very reason.

Of course, the above intuitive argument is inherently classical. One may immediately observe that it can fail in the BE-condensate case because the atoms are not localized entities — each atom’s wavefunction spreads out over distances much greater than the interatomic spacing. While that observation is enough to remove any logical paradox it only shows that the trick is possible; it does not really explain how the trick is done. The aim here is to offer a more tangible intuitive picture. (It should be recognized, however, that any intuitive explanation of an inherently quantum phenomenon will necessarily involve an element of “quantum mystery” that is not satisfactorily explicable in classical terms.)

In the standard theory of a uniform dilute Bose gas  the quasiparticle annihilation operators $`b_𝐤`$ are related to the atom annihilation operators $`a_𝐤`$ by a Bogoliubov transformation:

$$a_𝐤=\frac{1}{\sqrt{1\alpha _k^2}}\left(b_𝐤\alpha _kb_𝐤^{}\right),$$
(1)
with

$$\alpha _k=1+x^2x\sqrt{x^2+2},x^2k^2/(8\pi na),$$
(2)
where $`n`$ is the number density and $`a`$ is the scattering length. (We set $`\mathrm{}=1`$.) The quasiparticle spectrum is obtained as

$$\stackrel{~}{E}(p)=\left(\frac{1+\alpha _p}{1\alpha _p}\right)E(p)$$
(3)
where $`E(p)=p^2/2m`$ is the original atom spectrum. The resulting Bogoliubov spectrum is

$$\stackrel{~}{E}(p)=\frac{p}{2m}\sqrt{p^2+16\pi na}.$$
(4)
Thus, a low-momentum quasiparticle moves at a finite speed $`d\stackrel{~}{E}/dp\sqrt{4\pi na/m^2}v_s`$.

Due to interactions the ground state, even at zero temperature, does not have all the atoms in the $`𝐩=0`$ mode . The number of non-condensate atoms is

$$NN_0=\underset{𝐤0}{}a_𝐤^{}a_𝐤.$$
(5)
Substituting (1) and using the fact that the ground state is annihilated by $`b_𝐤`$ gives

$$0|a_𝐤^{}a_𝐤|0=\alpha _k^2/(1\alpha _k^2),$$
(6)
(which leads to the familiar result for the depletion $`1N_0/N=\frac{8}{3}\sqrt{\frac{na^3}{\pi }}`$.) A similar computation in a one-quasiparticle excited state $`|𝐩b_𝐩^{}|0`$ gives

$$𝐩|a_𝐤^{}a_𝐤|𝐩=\frac{1}{(1\alpha _k^2)}\left(\alpha _k^2+\delta _{𝐩,𝐤}+\alpha _p^2\delta _{𝐩,𝐤}\right),$$
(7)
which shows that an extra number $`N_{\mathrm{ex}}`$ of atoms have been excited from the condensate into the $`𝐩`$ and $`𝐩`$ modes. Specifically,

$$N_{\mathrm{ex}}=N_++N_{}=\frac{1}{1\alpha _p^2}+\frac{\alpha _p^2}{1\alpha _p^2}.$$
(8)
Note that $`N_+=N_{}+1`$, giving a net momentum $`𝐩`$ to the whole excitation.

When $`p`$ is large one has $`\alpha _p0`$, so that $`N_{\mathrm{ex}}1`$. Thus, a high-momentum quasiparticle is essentially a single moving atom. For small $`p`$, however, one has $`\alpha _p1`$, so that $`N_{\mathrm{ex}}`$ becomes very large. Thus, a low-momentum quasiparticle (a phonon) is a coherent state of many atoms with momenta $`𝐩`$ and $`𝐩`$.

Now we ask: What are all these excited atoms doing? …and how does the collective excitation of slow-moving atoms manage to move so fast? As an answer we suggest the picture in Fig. 1.

This picture shows a series of “creations” (where two condensate atoms scatter into $`𝐩`$ and $`𝐩`$ atoms) and “annihilations” (where a $`𝐩`$ and a $`𝐩`$ atom scatter back into two condensate atoms). (Each “creation” event “borrows” some energy that is repaid in the subsequent “annihilation” event.) The “creations” occur in rapid succession from left to right, resulting in an apparent motion at high speed. The situation can be seen more clearly in the spacetime picture, Fig. 2, in which each individual leg has a steep slope (the slope being the inverse of the atom speed), but the overall zig-zag path has a shallow slope (corresponding to the high sound speed).

Of course, this picture is not to be taken too literally since it assigns definite trajectories to quantum mechanical particles. \[What we have in mind is that path in Fig. 2, together with a set of similar, but more irregular zig-zag paths having the same average slope, form the class of paths that dominate the multiparticle path integral for the system.\] Note that the arrangement of “creation” and annihilation” events in quick succession requires a sort of “conspiracy” — but one that is possible given the quantum coherence of the condensate as a whole. (This “conspiracy” represents the element of “quantum mystery” referred to earlier as unavoidable in any attempt to intuitively explain an inherently quantum-mechanical phenomenon.)

Let us examine some consequences of the picture in Fig. 2 (taking it literally, despite our previous caveat). At any given instant, corresponding to the thin horizontal line, there are many excited atoms present; half left-moving and half (plus one extra) right-moving. If neighbouring creation events are separated in space by $`L`$ and separated in time by $`\tau `$ then the overall speed of the zig-zag path — the sound speed — is $`v_s=L/\tau `$. Let the time interval between a “creation” and a subsequent “annihilation” be $`T`$ (with $`T\tau `$ when $`N_{\mathrm{ex}}1`$). Each atom covers a distance $`L/2`$ in time $`T`$, giving an atom speed of $`v_a=\frac{L}{2T}`$. Hence the ratio $`\frac{v_s}{v_a}=\frac{L}{\tau }/\frac{L}{2T}=\frac{2T}{\tau }N_{\mathrm{ex}}`$. That is, the spacetime picture implies, for low momentum, the formula

$$v_s/v_a=N_{\mathrm{ex}}.$$
(9)
Since the velocities follow from differentiating the respective spectra, and since $`N_{\mathrm{ex}}`$ is given by (8), the above equation is

$$\frac{d\stackrel{~}{E}(p)}{dp}=\left(\frac{1+\alpha _p^2}{1\alpha _p^2}\right)\frac{dE}{dp}.$$
(10)
This relationship indeed follows from (3) after some algebra, employing the properties of $`\alpha _p`$ of (2). In fact Eq. (10) is true at all $`p`$, although the connection with our “cartoon picture” is valid only for small $`p`$.

The preceding considerations do not fix the absolute length scale. However, since Fig. 1 does indeed look like a density undulation with a wavelength $`L`$ we would expect $`L`$ to be $`2\pi /p`$ for momentum (wavenumber) $`p`$. We can justify this relation, in order-of-magnitude terms, as follows. Each “creation” event “borrows” an energy $`2E(p)=p^2/m`$, which is only “repaid” at the subsequent neighbouring “annihilations” a time $`T`$ later. In the usual way this violation of local energy conservation can hide under the uncertainty principle provided $`T1/(2E)`$. This in turn implies that $`L2Tv_a\frac{2m}{p^2}\frac{p}{m}`$, and so $`L`$ is of order $`1/p`$, as expected.

Experiments, notably by Ketterle’s group , have directly verified sound propagation in BE condensates. In such atom-trap experiments the longest wavelengths accessible are limited by the diameter of the cloud, hence $`p_{\mathrm{min}}\mathrm{}\pi /R`$, where $`R`$ is the Thomas-Fermi radius of the cloud in the appropriate direction. This consideration yields an upper limit

$$\frac{v_s}{v_a}0.22\frac{R}{\xi }$$
(11)
where $`\xi 1/\sqrt{8\pi na}`$ is the healing length, related to $`v_s`$ by $`v_s=\mathrm{}/(\sqrt{2}m\xi )`$. The ratio of length scales may may also be written as

$$\frac{R}{\xi }=\frac{R^2}{a_{\mathrm{osc}}^2}=\frac{a_{\mathrm{osc}}^2}{\xi ^2},$$
(12)
where $`a_{\mathrm{osc}}=\sqrt{\mathrm{}/m\omega }`$ is the nominal trap size, set by the angular frequency $`\omega `$ of the harmonic trapping potential.

For a typical $`{}_{}{}^{87}\mathrm{Rb}`$ trap with central density $`n=10^{15}\mathrm{cm}^3`$, $`\omega /2\pi =150\mathrm{Hz}`$, and $`a5.3\mathrm{nm}`$ one finds that $`v_s/v_a20`$. For the long $`{}_{}{}^{23}\mathrm{Na}`$ trap of Ref.  with $`2R450\mu \mathrm{m}`$ and $`\xi 0.2\mu \mathrm{m}`$, the maximum velocity ratio is about $`250`$. However, since the length of the pulses created was only about $`20\mu \mathrm{m}`$, the actual $`v_s/v_a`$ ratio for the phonon states studied was around $`20`$. That is still well into the region where the “puzzle” of $`v_sv_a`$ is an issue.

Acknowledgements

I thank Randy Hulet for helpful comments. This work was supported in part by the Department of Energy under Grant No. DE-FG05-97ER41031.


text:
# Star formation and Asymmetry in the spiral arms of M51 – Variable star formation caused by more than one spiral density wave

## 1 Introduction

The evolution of galaxies is influenced by how quickly gas is consumed by the formation of stars. Many empirical star formation recipes assume that the star formation rate depends on the gas surface density (e.g., Kennicutt 1998). However, observations of spiral galaxies show that the molecular gas is concentrated along spiral arms, implying that star formation is triggered by the spiral structure itself. The strongest evidence for the important role of spiral structure is based on a gravitational stability threshold, below which star formation ceases (Kennicutt, 1989). Nevertheless, above this threshold, azimuthally averaged data have shown that these empirical star formation laws are obeyed (Wong & Blitz 2002, Martin & Kennicutt 2001). In this paper we make use of recent high resolution mapping to estimate the rate and efficiency of star formation along individual spiral arms. We aim to investigate the influence of the spiral structure itself on the process of star formation.

M51 is a magnificent example of a nearby grand design spiral galaxy. Because of its large angular size and strong two-armed spiral structure, it has been the focus of studies of spiral density wave driven star formation (e.g., Tilanus & Allen 1991, Rand et al. 1992). While the stellar density contrast in the spiral arms in M51 is high compared to many spiral galaxies (Rix & Rieke, 1993), it is still far below the density contrasts seen in molecular or ionized gas. The distribution of molecular gas and associated star formation are particularly sensitive to the underlying spiral density waves propagating through the galaxy. Rix & Rieke (1993) in their Fourier analysis of the K-band light distribution measured significant deviations from bi-symmetric ($`m=2`$) symmetry in the galaxy. Because shocks occur in molecular gas, even weak perturbations in the stellar density distribution can cause large density variations in the molecular gas distribution. By measuring the shape and location of the peak densities seen in molecular gas, and comparing the shapes of the spiral arms with the star formation rates measured from the Pa$`\alpha `$ emission, we probe for differences in star formation rate and efficiency that might be caused by the spiral density waves themselves.

In this paper we take advantage of recent high angular resolution observations of M51 to measure the star formation efficiency along individual spiral arms. The BIMA Survey of Nearby Galaxies (BIMA-SONG, Helfer et al. 2003) is the first systematic interferometric survey of the disks of nearby galaxies in the CO(J=1-0) emission line. Prior to BIMA-SONG, the overwhelming majority of galaxies observed in CO had been observed with single dish telescopes, at much lower angular resolution. With a typical beam full width half max of $`6^{\prime \prime }`$, BIMA-SONG represents a major improvement in our ability to resolve the molecular gas distribution in a significant sample of nearby galaxies.

To measure the gas depletion time, which is inversely related to the star formation efficiency, we combine a measurement of the gas density with an estimate of the star formation rate. To measure the star formation rate, we use archival NICMOS data in the Pa$`\alpha `$ recombination line, at 1.87$`\mu `$m. The Pa$`\alpha `$ recombination line is better suited toward probing the possibly obscured central regions, since the near infrared is less affected by dust extinction. The high angular resolution of the NICMOS images on board HST facilitates continuum subtraction in crowded regions.

In Section 2, we describe the archival observations used for this study. To measure quantities along spiral arms, we must first determine their locations. In section 3.1 we describe our procedure for doing this and discuss our measurements that describe the shape of the arms. In section 3.2 we discuss the relative differences in gas density, star formation rate and gas depletion time measured in different spiral arm regions. In section 3.3 we discuss the effect of a possible additional spiral density wave on the gas density and on the symmetry of the spiral arms previously identified. A summary and discussion follows.

For M51 we adopt a distance of 7.7 Mpc (Tully, 1988), so $`1^{\prime \prime }`$ on the sky corresponds to a distance of 37 pc.

## 2 Observations

NICMOS narrow band (F187N and F190N) images were obtained as part of the NICMOS GTO program (Scoville et al., 2001) with a mosaic of nine fields covering an area of about $`2^{}\times 2^{}`$. See Scoville et al. (2001) for more information about these observations. The F187N and F190N filters, with effective wavelengths of 1.87 and 1.90$`\mu `$m contain Pa$`\alpha `$ line and continuum images, respectively. We scaled and subtracted the F190N continuum images from the F187N images and combined them to form one mosaic. The pixel scale of the NICMOS images is $`0.204^{\prime \prime }`$, so one pixel is 7.6 pc long in M51. The NICMOS images were calibrated using zero points listed in the HST Data Handbook. Specifically, we used $`8.723\times 10^{14}`$ erg cm-2 s-1/(DN s-1) for the continuum subtracted F187N narrow band images.

The BIMA-SONG database contains a collection of spatial velocity data cubes. For this study, we used an image made of the intensity integrated over all velocity channels. For M51, large scale fluxes were recovered in the interferometric data from single dish observations (Helfer et al., 2003).

To compare the CO image to the Pa$`\alpha `$ image we must ensure that the same areas of the sky can be measured on the same scale. Therefore, we resampled the CO image to be at the same pixel scale as the Pa$`\alpha `$ mosaic. Both images are inclination corrected and rotated so that the minor axis is at a P.A. of 170 degrees. The resulting images are shown in Figures 1 and 2.

## 3 Results

### 3.1 Logarithmic Spirals

To determine the location of the spiral arms, we model the spiral density waves in the disk with the following function, that describes a logarithmic spiral

$$\mathrm{\Sigma }\mathrm{Re}\left\{e^{i(\alpha _m\mathrm{ln}rm(\varphi \varphi _0))}\right\}$$
(1)
where $`\alpha `$ describes how tightly wound the arm is, $`\varphi _0`$ is an angular offset, and $`m`$ is an integer that corresponds to the number of spiral arms. In the case of M51, the spiral structure is predominantly two armed ($`m=2`$). One maximum is then described by

$$\alpha _m\mathrm{ln}r=m(\varphi +\varphi _0).$$
(2)
For each region we measured parameters $`\alpha _2`$ and $`\varphi _0`$ by choosing several positions of peak intensity along the spiral arms in the CO intensity map. All positions were measured from the center of the galaxy, which was determined from the centroid of an H band (F160W) image of M51 that was observed simultaneously with the Pa$`\alpha `$ data (Scoville et al., 2001). Since the nuclear region of M51 contains little molecular gas (see Figure 1), the H band centroid accurately measures the location of the stellar density peak. We defined the angle $`\varphi =0`$ to be toward the right in Figures 1 and 2. From these measurements, sets of coordinates, ($`r,\varphi `$), were measured in each arm. Linear regression was used to determine $`\alpha _2`$ and $`\varphi _0`$ in each region. Different regions in the spiral arms were fit separately, as a single spiral did not provide an adequate fit to the entire galaxy. The fits to the 8 selected spiral sub-regions are illustrated in Figure 1, and the corresponding parameters that describe them, are listed in Table 1.

In Figure 2 we compare our fit to the spiral arms traced in CO to the Pa$`\alpha `$ emission. While the Pa$`\alpha `$ is coincident with the molecular gas in some of the spiral arms, in other regions an offset is observed. Tilanus & Allen (1991) and Rand et al. (1992) noted an offset between CO and H$`\alpha `$ in the sense that H$`\alpha `$ emission was at somewhat larger radii than the CO emission. We find that the Pa$`\alpha `$ emission is offset from the CO emission at radii 1.2 to 2.8 kpc in the spiral arms 2b and 2c. However, we cannot tell if there is an offset between the Pa$`\alpha `$ emission and CO emission in spiral arm 1c because the edge of the Pa$`\alpha `$ image is near the position of the CO emission. We do note hints of Pa$`\alpha `$ emission, to the south of spiral arm 1c, lying on the very edge of the NICMOS mosaic. If there is Pa$`\alpha `$ emission to the south of the mosaic, then a significant offset between the CO and Pa$`\alpha `$ is present in this arm.

In addition to observing offsets between CO and Pa$`\alpha `$, we can look for deviations from point symmetry by plotting the reflection of the spiral arms determined in Table 1. In Figure 2 we show the fit to the spiral arms, rotated through an angle of 180 degrees. In the case of perfect $`m=2`$ symmetry, the over-plotted curve would lie on top of the peak intensities. We observe that the most inner and outer spiral arm regions appear to be symmetric, but for intermediate radii, the observed spiral arms are offset from the curves. To verify this, in Figure 3 we plot $`r`$ as a function of $`\varphi `$ for the a, b, and c spiral arms listed in Table 1, with the northern spiral arm rotated by an angle of 180 degrees. We note that the deviation from point symmetry is present over much of the inner part of the galaxy, but is most apparent for the b spiral arms, where the offset, $`\mathrm{\Delta }\varphi 20`$ degrees.

### 3.2 Star Formation Rates and Gas Depletion Times

In this section we present our measurements of the star formation rate (SFR) and CO mass, per unit length, along each segment of the spiral arms. Figures 4, 5, and 6 show the molecular gas density in each spiral arm segment and the star formation rate per unit length along the spiral arm segments. In addition, we calculated the gas depletion time from the ratio of the CO mass to the star formation rate. To measure these quantities, we used our own IDL procedure to measure the fluxes in particular regions. Each region was centered in radius on the position of a logarithmic spiral curve determined from our study described in section 3.1 (see Figures 1,2). The regions were defined by intervals in $`\varphi `$ and $`r`$, with $`\mathrm{\Delta }\varphi 0.05`$ radians and $`\mathrm{\Delta }r50`$ pixels or approximately 400 pc. The same regions were used to measure fluxes in both the CO image and the Pa$`\alpha `$ image, so that a one to one comparison could be made at each location. We ensured that $`\mathrm{\Delta }r`$ was large enough to include the Pa$`\alpha `$ emission, where offsets between Pa$`\alpha `$ and CO were observed. The measured fluxes were divided by arc length, so that integrated fluxes were computed per pc along each arm. Background sky values were determined in the Pa$`\alpha `$ image by taking the median of the regions adjacent to each box, on both sides of the spiral arms. Only spiral arm regions 1a, 1b, 2a, 2b, and 2c are compared, as regions 1c, 1d, and 2d do not lie on the NICMOS Pa$`\alpha `$ mosaic.

We next convert to physical units assuming an intrinsic flux ratio of $`F_{H\alpha }/F_{Pa\alpha }=8.46`$ for Case B recombination at a temperature of $`10^4`$K and a density $`100`$ cm-3 (line ratios are given in Osterbrock 1989). Then we converted the H$`\alpha `$ fluxes to an star formation rate, using a factor from Kennicutt, Tamblyn, & Congdon (1994),

$$SFR(M_{\mathrm{}}\mathrm{yr}^1)=\frac{L_{H_\alpha }}{(1.25\times 10^{41}\mathrm{ergs}\mathrm{s}^1)}.$$
(3)
The CO data were converted to M pc-2 using a CO to H2 conversion factor of $`2\times 10^{20}`$ H2 cm-2 (K km/s)-1 (Strong & Mattox, 1996).

In Figures 4, 5, and 6, we show the molecular gas mass, star formation rate, and gas depletion time, per unit length, for arms covered by both images. The plots are arranged such that the top plots show the molecular gas mass per pc, along the spiral arm segment. The middle plots show the star formation rate per pc, as measured from the Pa$`\alpha `$ mosaic. The bottom plots show the ratios of the points in the two plots above them. This quantity is known as the gas depletion time, the time in which all of the molecular gas would form into stars, if the star formation rate was to remain constant. In Figure 4, spiral arm 1a was placed next to its opposite spiral arm, 2a, and the same is true for arms 1b and 2b, in Figure 5. In Figure 6, spiral 2c is shown alone, because we were unable to measure Pa$`\alpha `$ emission for spiral arm 1c. Slightly negative star formation rates, the result of imperfect sky subtraction in the Pa$`\alpha `$ image, are assumed to be sites where essentially no stars are being formed.

Figures 4 and 5 can be used to look for differences in the gas density and star formation efficiency between a spiral arm and its partner on the opposite side of the galaxy. In Figure 4 we note that the gas density is twice as high on NW side (1a) than the SE side (2a). However, the star formation efficiency is about the same with a gas depletion time $`2\times 10^8`$ years. In Figure 5 we note a similarity, in the sense that the gas density is again twice as high in the western spiral arm (1b) as in the eastern spiral arm (2b). Both spiral arms 1b and 2b have a jump in gas depletion time (a decrease in star formation efficiency), at a radius of 1.5 kpc. At radii smaller than 1.5 kpc, the gas depletion time in spiral arm regions 1b and 2b is similar to that in 1a and 2a. At larger radii the gas depletion time is several times larger than within this radius. Star formation in the outer parts of spiral arms 1b and 2b is inefficient compared to that in the inner parts of the same arms. Past a radius of 2 kpc the gas depletion time drops back to $`2\times 10^8`$ years, and remains similar in spiral arm 2c. We suspect that the star formation efficiency in arm 1c would be similar to that in 2c, although we are unable to observe it.

Except for the outer parts of spiral arms 1b and 2b at a radius of about 1.5 to 2kpc, the star formation efficiency or gas depletion time is remarkably similar in all arms, and approximately $`2\times 10^8`$ yr, despite large variations in gas density and the star formation rate. In the outer parts of 1b and 2b the gas depletion time is longer, $`10^9`$ years, despite the high molecular gas densities in these arms.

### 3.3 Deviations from $`m=2`$ symmetry

In this section we consider a dynamical model that can explain the degree of deviation from $`m=2`$ symmetry that we have measured in the spiral arms of M51. By expanding galaxy images in Fourier components, Elmegreen et al. (1992) noted that many galaxies exhibit hidden three armed components and consequently suggested that multiple spiral density waves often propagate simultaneously in galaxy disks. Rix & Rieke (1993) in their Fourier expansion of the near-infrared images of M51 saw both $`m=1`$ and $`m=3`$ components as well as the dominant 2-armed structure. Rix & Rieke (1993) found that the galaxy was lopsided at all radii but also measured a weak peak in the $`m=3`$ component of the K-band images at a radius of about $`35^{\prime \prime }`$. This radius is approximately at the location where we have found the largest asymmetry in the spiral structure, although from Figure 3 we note that the asymmetry is present throughout most of the galaxy.

We now consider the possibility that more than one spiral density wave pattern is present in the disk of M51. Below we estimate the likely position of molecular gas along spiral arms resulting from more than one perturbation. We assume that the perturbations in the gravitational potential resulting from spiral modes can be described in terms of Fourier components

$$\mathrm{\Phi }(r,\varphi )=\underset{m}{}A_m\mathrm{cos}(m\varphi \alpha _m\mathrm{ln}r\beta _m)$$
(4)
where $`A_m`$ and $`\beta _m`$ are amplitude and phase shift, and $`m\varphi \alpha _m\mathrm{ln}r`$ describes an $`m`$ armed logarithmic spiral (see section 3.1). Because spiral density waves travel at pattern speeds above the sound speed in the interstellar medium, the passage of spiral modes in the disk causes shocks in the gas. The concentration of molecular clouds along spiral arms is interpreted to be a result of these galactic scale shocks. Shocks in the gas are likely to occur near maxima in the stellar density field, approximately corresponding to minima in the gravitational potential. These are found at angles where $`d\mathrm{\Phi }/d\varphi =0`$ or

$$\underset{m}{}A_mm\mathrm{sin}(m\varphi \alpha _m\mathrm{ln}r\beta _m)=0$$
(5)
This equation is also valid at potential maxima where shocks will not occur. We assume that one spiral component dominates; in the case of M51 this would be the $`m=2`$ or two-armed mode. In the absence of other modes, the potential minima (density maxima) occur near angles

$$2\varphi _{2,max}=2n\pi +\alpha _2\mathrm{ln}r+\beta _2$$
(6)
where $`n`$ is an integer. We expand Equation 5 about the angle $`\varphi _{2,max}`$, for $`\varphi =\varphi _{2,max}+\delta \varphi `$, where the offset, $`\delta \varphi `$, is the difference between the location of density maximum caused by only one $`m=2`$ wave and that resulting from two spiral density waves. To first order in $`\delta \varphi `$,

$$4A_2\delta \varphi +A_mm\left[\mathrm{sin}(m\varphi _{2,max}\alpha _m\mathrm{ln}r\beta _m)+\mathrm{cos}(m\varphi _{2,max}\alpha _m\mathrm{ln}r\beta _m)m\delta \varphi \right]=0$$
(7)
where we have assumed only one spiral density wave with $`m2`$, in addition to the dominant two-armed one. We solve for $`d\varphi `$

$$\delta \varphi =\frac{A_mm\mathrm{sin}(m\varphi _{2,max}\alpha _m\mathrm{ln}r\beta _m)}{4A_2+A_mm^2\mathrm{cos}(m\varphi _{2,max}\alpha _m\mathrm{ln}r\beta _m)}.$$
(8)
Substituting for $`\varphi _{2,max}`$ (Equation 6) and taking the limit $`A_m<A_2`$

$$\delta \varphi \frac{A_mm}{4A_2}\mathrm{cos}(nm\pi )\mathrm{sin}\left((\frac{m\alpha _2}{2}\alpha _m)\mathrm{ln}r+\frac{m\beta _2}{2}\beta _m\right).$$
(9)
We use our model to estimate the angular offset that would be measured between a spiral arm and its opposite one after a rotation of 180 degrees. For a dominant two-armed mode, each arm has a different integer $`n`$ in Equation (9); $`n=0`$ or $`1`$. Only when $`m`$ is odd will the angular offsets be asymmetric, causing reflection symmetry about the origin to be broken. The angular offset between one spiral arm and the other that has been rotated about $`\pi `$ is then given by subtracting the $`n=0`$ offset from that of the opposite one ($`n=1`$),

$$\mathrm{\Delta }\varphi \frac{A_mm}{2A_2}\mathrm{sin}\left((\frac{m\alpha _2}{2}\alpha _m)\mathrm{ln}r+\frac{m\beta _2}{2}\beta _m\right)$$
(10)
for odd $`m`$. We expect that the angular difference depends upon the amplitude of the additional spiral density wave and on its wave number. The angular offset should increase as the amplitude of the additional spiral density wave increases. In a region where the amplitudes of the two dominant spiral density waves are constant, we expect that the angular offset would be more open (less tightly wound) than the spiral density waves themselves.

We now relate the size of the angular offset to the one expected based on observed density perturbations or observed Fourier components. In the WKB approximation $`A_m2G\mathrm{\Sigma }_mr/|m\alpha _m|`$ where $`\mathrm{\Sigma }_m`$ is magnitude of the $`m`$th Fourier component of the mass surface density (Binney & Tremaine, 1987), so that

$$\frac{A_m}{A_n}\frac{\mathrm{\Sigma }_mn|\alpha _n|}{\mathrm{\Sigma }_nm|\alpha _m|}.$$
(11)
Substituting this into equation(9) and using the ratios $`\frac{\mathrm{\Sigma }_3}{\mathrm{\Sigma }_0}0.05`$, $`\frac{\mathrm{\Sigma }_2}{\mathrm{\Sigma }_0}0.15`$, and $`\frac{\mathrm{\Sigma }_1}{\mathrm{\Sigma }_0}0.15`$ at $`r=35^{\prime \prime }`$ measured by (Rix & Rieke, 1993) we estimate

$`\mathrm{\Delta }\varphi \mathrm{sin}\left(({\displaystyle \frac{\alpha _2}{2}}\alpha _1)\mathrm{ln}r+{\displaystyle \frac{\beta _2}{2}}\beta _1\right)`$ $`\mathrm{for}m=1`$
$`\mathrm{\Delta }\varphi 0.3\mathrm{sin}\left(({\displaystyle \frac{3\alpha _2}{2}}\alpha _3)\mathrm{ln}r+{\displaystyle \frac{3\beta _2}{2}}\beta _3\right)`$ $`\mathrm{for}m=3`$ (12)
where the angles are given in radians. Assuming $`\alpha _2\alpha _m`$, and for a maximum value of the angular offsets ($`\beta _2`$ with respect to $`\beta _m`$) the angular shift could be as large as 60 degrees for $`m=1`$, and 20 degrees for $`m=3`$.

We compare the estimated model angular offsets to those we measured in M51. Recalling from Figure 3, the offset measured between the 1b and 2b spiral arms is $`\mathrm{\Delta }\varphi 20`$ degrees. So the $`m=3`$ Fourier component observed by Rix & Rieke (1993) is sufficiently large to account for the angular offset we measured between these two arms. A three-armed mode, in addition to the dominant two-armed one, consistent with perturbations seen in the K-band Fourier components, could explain the asymmetry in the locations of the two spiral arms seen in molecular gas. Figure 3 shows that while the asymmetry is present over much of the galaxy, the maximum offset occurs where Rix & Rieke (1993) observe a peak in the $`m=3`$ mode. The $`m=1`$ mode, on the other hand, is strong over much of the entire galaxy and does not single out any particular radial region, though it could be responsible for asymmetry at larger and smaller radii, where the $`m=3`$ amplitude is weak. The decrease of the angular offset at small and large radii ($`r<1.2`$kpc and $`r>2.3`$kpc) can be explained by the decrease in amplitude of the $`m=3`$ wave. The nearly constant value of the angular offset between $`1`$ and $`2.2`$ kpc suggests that the wave numbers of the $`m=2`$ and $`m=3`$ modes are approximately the same size ($`\frac{3\alpha _2}{2}\alpha _3`$, see equation 12).

If there is both $`m=2`$ and $`m=3`$ waves, then at least one density peak of the $`m=3`$ wave will lie in between the dominant $`m=2`$ arms. We might expect inter-arm molecular gas and associated inter-arm star formation in the same region where we see the deviations from bi-symmetry. A more careful inspection of Figure 1 shows that there is a patch of interarm star formation within spiral arm 2b to the north east of the nucleus. As we have done for the molecular gas distribution (see Section 3.1), we crudely fit a spiral shape to the the Pa$`\alpha `$ emission in this region. The resulting fit is shown in Figure 7, along with the same arm rotated by 120 and 240 degrees. From this figure, we can determine where the $`m=3`$ and $`m=2`$ waves are likely to constructively add and where they interfere. The $`m=3`$ wave is coincident with the dominant spiral arm in the west (1a-1b), but is offset from the spiral arm in the southeast (2a-2b). We find that the $`m=3`$ wave can account for the molecular gas density enhancement in the western spiral arms, relative to the eastern spiral arms (see Figures 4, 5, and 6).

The addition of an extra wave should also introduce an asymmetry in the velocity field. If this asymmetry is sufficiently large then it could be detected in the CO velocity field. To search for this we examined the mean velocity MOM1 BIMA map of M51. We first subtracted the systemic velocity of the galaxy from the velocity field and then compared this field to the same field that has been rotated by 180 degrees and multiplied by $`1`$. We find that the line of sight velocity on the southern side of the galaxy 2kpc from the nucleus is about 20-30 km/s above that on the opposite side. The velocity asymmetry seems to be largest in the same region where we see the other asymmetries. The magnitude of the velocity asymmetry should be something that our simple model could predict. However, it depends on the pattern speed of the two waves and the details of the gas dynamics so a simulation is probably required to make a prediction sufficiently accurate to compare to the complex streaming motions seen in the CO channel maps.

## 4 Summary and Discussion

In this paper we have fit the spiral arms of M51 as traced in CO within 3kpc of the nucleus with logarithmic spirals. We have detected an asymmetry between the locations of the arms on one side of the galaxy compared to those on opposite side. The asymmetry corresponds to a maximum angular offset of about $`20`$ degrees between one arm and the opposite one rotated by 180 degrees, and is most prominent between 1 and 2.2 kpc from the nucleus. The molecular gas density in this region is also asymmetric with one spiral arm about twice as dense as that on the opposite side. In the same radial region, there is a patch of star formation (seen in the Pa$`\alpha `$ images located between the dominant spiral arms to the north east of the galaxy nucleus.

To explain the angular offset seen between one arm and the opposite one rotated by $`\pi `$, we consider a model that contains two spiral density waves; a weak $`3`$-armed wave in addition to the dominant two-armed one. When the additional wave is odd-armed, we predict that the position of the two dominant arms deviate from point symmetry in the manner we have observed in M51. We find that an additional $`m=3`$ mode, consistent with the radial location and strength of the $`m=3`$ Fourier components previously measured at K-band by Rix & Rieke (1993), is strong enough to account for the angular offset we have measured between the two dominant arms. This additional spiral density wave also is consistent with the lopsided distribution of molecular gas, and the inter-arm star formation to the north east of the galaxy nucleus. Despite molecular gas density variations of a factor of a few and corresponding star formation rate variations (measured from the Pa$`\alpha `$ emission), we see little variation in the computed gas depletion times. This suggests that the star formation efficiency is not largely affected by the large gas density variations caused by the presence of an additional spiral density wave.

Between 1 and 2.2 kpc in M51 we see large scale deviation from pure $`m=2`$ or bi-symmetric structure: an angular offset between the spiral arms, differences in density on one side compared to the other side, and inter-arm star formation. All of these three asymmetries are most pronounced in the same radial region suggesting that they are related. We have suggested here that this phenomena can be explained with one simple unifying model that includes an additional weak 3-armed spiral density wave. The additional wave does not strongly influence the star formation efficiency, however because of the large gas density variations caused by it, the wave does influence the star formation rate. In particular, large molecular gas densities and associated high star formation rates are likely to occur where and when the 3-armed and 2-armed mode both reach their density maxima simultaneously. We expect that the 3-armed wave has a different pattern speed (angular rotation rate) than the 2-armed wave. Consequently the star formation rate which is highest when the 3 and 2 armed patterns constructively add, must be highly variable. Even low amplitude waves can cause large changes in the gas density. Despite the low estimated amplitude of the 3-armed wave ($`\frac{\mathrm{\Sigma }_3}{\mathrm{\Sigma }_0}0.05`$) it is likely that it is causing large variations in the molecular gas density, and consequently in the star formation rate. The simultaneous propagation of 2 and 3 armed structure must cause associated large variations in the local star formation rate. Since our model contains more than one spiral density wave we would more strongly support scenarios that include multiple and transient spiral density waves (e.g., Toomre 1981; Fuchs 2001) rather than those that focus on a single dominant quasi-steady mode (e.g., Lowe et al. 1994).

The model containing two and three armed spiral density waves, that we have proposed to explain the asymmetries in the central 3kpc of M51, does not account for all the phenomena we have described in this paper. For example, this model does not account for the drop in star formation efficiency observed in the outer parts of spiral arms 1b and 2b between 1.5 and 2 kpc from the nucleus. It also does not account for the offsets between the CO and Pa$`\alpha `$ emission which are largest at radii greater than 2kpc. Similar offsets have previously been discussed by Rand et al. (1999); Rand (1995); Rand et al. (1992); Tilanus & Allen (1991) who have offered a variety of explanations including the triggering of star formation by spiral density waves, uneven heating and disassociation by the resulting star formation and differential clumping of the molecular gas. However, it will be tempting to search for an additional explanation for these offsets and the decrease in the star formation efficiency we see in the outer arms 1b and 2b that also involve more than one spiral density wave. Spiral structure exists in M51 over a large range of radius. In section 3.1 we attempted to fit as large as possible pieces of the spiral structure with as few as possible individual components. However, it is not clear from Figure 1,2 or 3 that the ranges of radius over which we fit individual arm pieces correspond to separate spiral density waves, traveling at individual pattern speeds. It is likely that the inner spiral arms of M51 rotate at a different pattern speed than the outer ones, and the interplay between these waves might provide an alternative explanation for the offsets observed between the CO and Pa$`\alpha `$ emission or the low star formation efficiency in the outer parts of spiral arms 1b and 2b. Future work can further test the possibility that multiple spiral density waves propagate in galaxies and explore the role that these multiple modes exert on the distribution of gas and resulting pattern of star formation. Future studies can extend this type study carried out here to more galaxies, extend the multi-wavelength coverage and exploit the velocity fields that are available in the CO and HI data.

We thank Joel Green, Judith Pipher, Dan Watson, Bruce Elmegreen, Alar Toomre, Mousumi Das and Peter Teuben for stimulating discussions. We are grateful to Tamara Helfer for providing us with access to the full BIMA SONG galaxy sample data in advance of publication.