text:
# Spin-dependent Bohm trajectories for Pauli and Dirac eigenstates of hydrogen

## 1 INTRODUCTION

In Bohm’s original causal interpretation of quantum mechanics , the motion of a quantum particle is determined by its Schrödinger wave function $`\psi `$, which acts as a kind of guidance wave . If the wave function is written as

$$\psi (𝐱,t)=R(𝐱,t)e^{iS(𝐱,t)/\mathrm{}},$$
(1)
where $`R`$ and $`S`$ are real-valued, then the trajectory of the particle is determined by the guidance relation

$$𝐩=S.$$
(2)
The momentum is related to the well-known Schrödinger current j as follows,

$$𝐩=\frac{m}{\rho }𝐣$$
(3)
where $`\rho =\psi ^{}\psi =R^2`$. Comprehensive discussions of the de Broglie-Bohm causal interpretation can be found in  and .

It is quite natural to examine the hydrogen atom, one of the simplest quantum systems, in terms of the de Broglie-Bohm theory. Indeed, as originally discussed in , the Schrödinger guidance relation (2) predicts that $`𝐩=\mathrm{𝟎}`$ for all real eigenstates, including the ground and all higher $`s`$ hydrogenic states. However, as Holland  pointed out, Eq. (2) is valid only for spinless particles. For particles with spin, the condition of Lorentz covariance on the law of motion implies that the momentum of a particle with spin s must be given by

$$𝐩=S+\mathrm{log}\rho \times 𝐬,$$
(4)
where $`\rho =\psi ^{}\psi `$ . Only in this way can the theory be embedded in a relativistic formulation. Indeed, in papers , it was shown that in order for it to be consistent with Dirac theory, the Schrödinger equation must be regarded as describing an electron in a definite eigenstate of spin. In these papers, the current vector associated with (4),

$$𝐣=\frac{1}{m}\rho S+\frac{1}{m}\rho \times 𝐬,$$
(5)
was referred to as the Pauli current, the nonrelativistic limit of the Dirac current.

The guidance law (4) no longer implies that $`𝐩=\mathrm{𝟎}`$ for real eigenstates so it is natural to ask how it applies to the hydrogen atom. In , we showed that for an electron in a spin eigenstate with $`s_z=\pm \frac{1}{2}`$, the spin-dependent term in (4) will be responsible for a motion in a plane perpendicular to the $`z`$-axis and along a contour of constant $`\rho `$ value. For the case of an electron in a Schrödinger energy/angular momentum eigenstate, $`\psi _{nlm}`$, this implies circular motion about the $`z`$-axis.

In this paper, we examine de Broglie-Bohm trajectories for an electron in a hydrogen atom as described by the Pauli and Dirac equations using appropriate (spin-dependent) currents. The electron is assumed to be in a Pauli or Dirac eigenstate of energy and total angular momentum. Note that, in general, this does not imply that the electron is in a spin eigenstate of known $`s_z`$ value. We also show that under appropriate nonrelativistic limits, the angular rotation rates for Dirac trajectories become those of Pauli trajectories. In the cases that the electron is in a spin eigenstate (e.g., $`1s`$, $`2s`$, $`2p_0`$), the Pauli rotation rates agree with the Schrödinger trajectories obtained in . The result is a coherent application of de Broglie-Bohm theory to relativistic and nonrelativistic hydrogen atom eigenstates.

For both the Dirac and Pauli cases, the Schrödinger guidance formula in (4) can be generalized using the relationship (3) where j is the appropriate (Dirac/Pauli) current and $`\rho =\psi ^{}\psi `$. First, consider the Dirac equation,

$$i\mathrm{}\frac{}{t}\psi =(e\varphi +\beta E_o+𝜶(c𝐩+e𝐀))\psi .$$
(6)
Here, $`\psi =(\psi _1,\psi _2,\psi _3,\psi _4)`$ is a four-component wave function, $`\varphi `$ and $`𝐀`$ are the scalar and vector potentials, $`E_0`$ and $`𝐩`$ are the rest mass energy and momentum operators, $`e`$ is the electric charge, and $`𝜶`$ and $`\beta `$ are the Dirac operators. In this study, $`𝐀=\mathrm{𝟎}`$ and the current is given by

$$𝐣=c\psi ^{}𝜶\psi =(j_x,j_y,j_z),$$
(7)
where the $`𝜶`$ are the $`4\times 4`$ Dirac matrices,

$$𝜶=\left(\begin{array}{cc}0& 𝝈\\ 𝝈& 0\end{array}\right),$$
(8)
and the $`𝝈`$ are the $`2\times 2`$ Pauli matrices.

If the particle is in a potential such that $`e\varphi mc^2`$, then there exist stationary states for which the average velocity $`\overline{v}`$ is nonrelativistic, and $`EE_0=mc^2`$. In this case, the latter two components of the 4-component Dirac state are smaller in magnitude than the first two components by a factor of $`\overline{v}/c`$. The Dirac equation may then be reduced to the Pauli equation involving the two components $`\psi _1`$ and $`\psi _2`$:

$$i\mathrm{}\frac{\psi }{t}=\frac{1}{2m}(i\mathrm{}+e𝐀)^2\psi +\frac{e\mathrm{}}{2m}𝝈𝐁\psi eV\psi .$$
(9)
Here $`\psi =(\psi _1,\psi _2)`$ is a two-component Pauli spinor wave function. Once again, we assume that $`𝐀=\mathrm{𝟎}`$ so that the associated Pauli current is given by

$$𝐣=𝐣_A+𝐣_B=\frac{\mathrm{}}{2mi}(\psi ^{}\psi \psi \psi ^{})+\frac{\mathrm{}}{2m}\times (\psi ^{}\sigma \psi ).$$
(10)
Note that if one assumes that the system is in an eigenstate of the spin operator, then (10) and (3) together reduce to (4).

In the case of the Schrödinger equation for hydrogen, it is usually assumed that spin interactions are negligible so that the wave function can be written as a product of spatial- and spin-dependent terms, i.e.,

$$\psi =\psi (𝐫,t)\chi _s.$$
(11)
As is well known, the spatial hydrogenic energy Schrödinger eigenfunctions,

$$\psi _{n,\mathrm{},m}(r,\theta ,\varphi )=R_{n,\mathrm{}}(r)Y_{lm}(\theta ,\varphi ),$$
(12)
solutions to the (spinless) time-independent Schrödinger equation, are also simultaneous eigenstates of the orbital angular momentum operator $`L^2`$, with eigenvalues $`\mathrm{}^2\mathrm{}(\mathrm{}+1)`$, and the operator $`L_z`$, with eigenvalues $`\mathrm{}m`$.

In the Pauli equation, where spin-orbit interactions are excluded, the orbital angular momentum operator $`L^2`$ commutes with the hamiltonian. This is not the case for the Dirac equation. For both the Pauli and Dirac equations, however, each component of $`M`$, the total angular momentum operator, commutes with the Hamiltonian $`H`$, implying that $`M^2`$ commutes with $`H`$ as well. For this reason, it is conventional to choose eigenstates of $`H`$, $`M^2`$ and $`M_z`$, with eigenvalues $`E_n`$, $`\mathrm{}^2j(j+\frac{1}{2})`$ and $`\mathrm{}m`$, respectively.

There is one further subtlety: Although the orbital angular momentum does not commute with the hamiltonian in the Dirac case, it can be shown that $`\mathrm{}`$ is ‘almost’ a good quantum number (see ). That is to say, eigenstates can be found for which

$$𝐋^2\psi =\mathrm{}^2\mathrm{}(\mathrm{}+1)\psi +\mathrm{}^2w,$$
(13)
where the spinor $`w`$ is negligible. (Its large components actually vanish.) Hence for both the Dirac and Pauli cases considered in this paper, eigenstates are presented in terms of quantum numbers $`n`$, $`\mathrm{}`$, $`j`$ and $`m`$ for purposes of comparison, even though $`\mathrm{}`$ is not strictly a good quantum number in the Dirac case.

Finally, in the following discussions, the time-dependent phase factor $`e^{iE_nt/\mathrm{}}`$ that accompanies the eigenfunctions in the solution of the time-dependent Pauli and Dirac equations will be ignored since it contributes nothing to the associated currents.

## 2 PAULI EIGENSTATES

In this section we examine the Pauli current (10) for some hydrogen atom eigenstates. These eigenstates, two-component solutions to the Pauli equation, are given by

$$\begin{array}{cc}\hfill \psi _{n,\mathrm{},j=\mathrm{}+\frac{1}{2},m}& =\frac{1}{\sqrt{2\mathrm{}+1}}R_n\mathrm{}(r)\left(\begin{array}{c}\sqrt{\mathrm{}+m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m\frac{1}{2}}(\theta ,\varphi )\\ \sqrt{\mathrm{}m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m+\frac{1}{2}}(\theta ,\varphi )\end{array}\right)\hfill \\ \hfill \psi _{n,\mathrm{},j=\mathrm{}\frac{1}{2},m}& =\frac{1}{\sqrt{2\mathrm{}+1}}R_n\mathrm{}(r)\left(\begin{array}{c}\sqrt{\mathrm{}m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m\frac{1}{2}}(\theta ,\varphi )\\ \sqrt{\mathrm{}+m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m+\frac{1}{2}}(\theta ,\varphi )\end{array}\right)\hfill \end{array}$$
(14)
Here, the $`R_{nl}(r)`$ are the standard radial wave functions for the hydrogen atom and the $`Y_{l,m\pm \frac{1}{2}}(\theta ,\varphi )`$ are the spherical harmonics.111We use the following convention for the relevant functions, for consistency with :    $`Y_\mathrm{}m(\theta ,\varphi )={\displaystyle \frac{1}{\sqrt{2\pi }}}P_\mathrm{}m(\mathrm{cos}\theta )e^{im\varphi }`$    $`P_\mathrm{}m(x)=\sqrt{{\displaystyle \frac{2\mathrm{}+1}{2}}{\displaystyle \frac{(\mathrm{}m)!}{(\mathrm{}+m)!}}}{\displaystyle \frac{1}{2^{\mathrm{}}\mathrm{}!}}(1x^2)^{m/2}{\displaystyle \frac{d^{\mathrm{}+m}}{dx^(\mathrm{}+m)}}(x^21)^{\mathrm{}},m0`$    $`P_{\mathrm{},m}(x)=(1)^mP_\mathrm{}m(x),m<0.`$   We use spherical polar coordinates in which $`r`$ is the radius, $`\varphi `$ is the angle measured counterclockwise from the $`x`$-axis and $`\theta `$ is the angle measured down from the $`z`$-axis.

As mentioned earlier, the wave functions given in (14) are eigenfunctions of $`L^2`$ (the orbital angular momentum) with eigenvalue $`\mathrm{}^2\mathrm{}(\mathrm{}+1)`$, $`M^2`$ (the total angular momentum) with eigenvalue $`\mathrm{}^2j(j+1)`$ and $`M_z`$ with eigenvalue $`m\mathrm{}`$. In general, however, they are not eigenstates of $`s_z`$, the projection of spin along the $`z`$-axis.

The eigenfunctions can be classified as follows:

* For each $`n`$ value, $`\mathrm{}`$ can assume the values $`\mathrm{}=0,1,\mathrm{},n1`$.
* For each $`\mathrm{}`$ value, $`m`$ can assume the values $`m=\mathrm{}+1/2,\mathrm{}+3/2,\mathrm{},\mathrm{}1/2`$.
* For each of the above there are two possibilities, $`j=\mathrm{}+1/2`$ and $`j=\mathrm{}1/2`$.

This accounts for all eigenfunctions listed in (14).

From (10), the two contributions to the velocity are given by

$$𝐯_a=\frac{𝐣_a}{\rho }=\frac{\mathrm{}}{2m_ei}\frac{(\psi ^{}\psi \psi \psi ^{})}{\psi ^{}\psi }=\frac{\mathrm{}}{m_e}\frac{\text{Im}(\psi ^{}\psi )}{\psi ^{}\psi }$$
(15)
and

$$𝐯_b=\frac{𝐣_b}{\rho }=\frac{\mathrm{}}{2m_e}\frac{\times 𝐬}{\psi ^{}\psi },$$
(16)
where

$$𝐬=\psi ^{}\sigma \psi $$
(17)
is the ‘spin vector’ and $`m_e`$ is the mass of the electron. In order to determine trajectories for the above hydrogen eigenstates we must compute the velocities $`𝐯_a`$ and $`𝐯_b`$ for the wave functions given in (14).

We first examine the velocity $`𝐯_a`$ arising from the Schrödinger current $`𝐣_a`$. Writing

$$\psi =\frac{1}{\sqrt{2\mathrm{}+1}}R_n\mathrm{}(r)\left(\begin{array}{c}v_1(\theta ,\varphi )\\ v_2(\theta ,\varphi )\end{array}\right),$$
(18)
the term $`\text{Im}\{\psi ^{}\psi \}`$ can be shown to be

$$\begin{array}{c}\text{Im}\{\psi ^{}\psi \}=\frac{1}{2\mathrm{}+1}\text{Im}\{R_n\mathrm{}R_n\mathrm{}^{}(|v_1|^2+|v_2|^2)\widehat{e}_r+\frac{1}{r}R_n\mathrm{}^2(v_1^{}\frac{v_1}{\theta }+v_2^{}\frac{v_2}{\theta })\widehat{e}_\theta \hfill \\ \hfill +\frac{1}{r\mathrm{sin}\theta }R_n\mathrm{}^2(v_1^{}\frac{v_1}{\varphi }+v_2^{}\frac{v_2}{\varphi })\widehat{e}_\varphi \}.\end{array}$$
(19)
In the above, the only nonzero term comes from the $`\widehat{e}_\varphi `$ component:

$$v_k^{}\frac{v_k}{\theta }=i(m\pm \frac{1}{2})|v_k|^2,k=1,2,$$
(20)
so that (15) yields

$$𝐯_a=\frac{\mathrm{}}{m_er\mathrm{sin}\theta }\left(m+\frac{1}{2}(\frac{|v_2|^2|v_1|^2}{|v_1|^2+|v_2|^2})\right)\widehat{e}_\varphi .$$
(21)
It remains to compute $`𝐣_b`$ and the corresponding velocity $`𝐯_b`$, with reference to (16). To do this, we first find the three components of the spin vector $`𝐬=\psi ^{}\sigma \psi `$:

$$\psi ^{}\sigma \psi =\frac{1}{2\mathrm{}+1}R_n\mathrm{}(r)^2(2\text{Re}\{v_1^{}v_2\},2\text{Im}\{v_1^{}v_2\},|v_1|^2|v_2|^2).$$
(22)
From the form of the wave functions (14),

$$2\text{Re}(v_1^{}v_2)=2c_1c_2N_1N_2P_{\mathrm{}}^{m\frac{1}{2}}(\theta )P_{\mathrm{}}^{m+\frac{1}{2}}(\theta )\mathrm{cos}\varphi $$
(23)
and

$$2\text{Im}(v_1^{}v_2)=2c_1c_2N_1N_2P_{\mathrm{}}^{m\frac{1}{2}}(\theta )P_{\mathrm{}}^{m+\frac{1}{2}}(\theta )\mathrm{sin}\varphi ,$$
(24)
where the $`c_i`$ are given in (14) and the $`N_i`$ are the normalization constants of the relevant spherical harmonics. For simplicity of notation, define

$$a=c_1N_1P_{\mathrm{}}^{m\frac{1}{2}},b=c_2N_2P_{\mathrm{}}^{m+\frac{1}{2}}.$$
(25)
If we write

$$\psi ^{}\sigma \psi =\frac{1}{2\mathrm{}+1}R_n\mathrm{}^2(r)𝐰,$$
(26)
then the vector $`𝐰`$ can be expressed in Cartesian form as

$$(w_x,w_y,w_z)=(2ab\mathrm{cos}\varphi ,2ab\mathrm{sin}\varphi ,a^2b^2).$$
(27)
We may also express $`𝐰`$ in spherical polar form, i.e.,

$$\begin{array}{cc}\hfill w_x& =r_s\mathrm{sin}\theta _s\mathrm{cos}\varphi _s\hfill \\ \hfill w_y& =r_s\mathrm{sin}\theta _s\mathrm{sin}\varphi _s\hfill \\ \hfill w_z& =r_s\mathrm{cos}\theta _s,\hfill \end{array}$$
(28)
where the orientation of the spin vector $`𝐬=\psi ^{}\sigma \psi `$ is given by the angles $`\theta _s`$ and $`\varphi _s`$. (See also , p. 62-63 for a brief discussion of the spin vector.) A comparison of (28) with (27) suggests that we might let $`r_s=a^2+b^2`$, $`\varphi _s=\varphi `$ and then compute $`\theta _s`$ in terms of $`\theta `$ using the relations

$$\mathrm{cos}\theta _s=\frac{a^2b^2}{a^2+b^2},\mathrm{sin}\theta _s=\frac{2ab}{a^2+b^2}.$$
(29)
However, this is consistent with the definition of spherical coordinates only if $`2ab0`$ since $`\theta _s`$ is restricted to the interval $`[0,\pi ]`$. When this condition is not met, i.e., $`2ab<0`$, then the polar coordinates for $`𝐰`$ are given by $`r_s=a^2+b^2`$, $`\varphi _s=\varphi +\pi `$ and

$$\mathrm{cos}\theta _s=\frac{a^2b^2}{a^2+b^2},\mathrm{sin}\theta _s=\frac{|2ab|}{a^2+b^2}=\frac{2ab}{a^2+b^2}.$$
(30)
In either of the above cases, the spin vector $`𝐬`$ lies in a plane defined by the position vector $`𝐫`$ and the $`z`$ axis, which is in agreement with . After some manipulation, we find that

$$𝐬=s_r\widehat{e}_r+s_\theta \widehat{e}_\theta ,$$
(31)
where

$$\begin{array}{cc}\hfill s_r& =s\mathrm{cos}\theta (\frac{a^2b^2}{a^2+b^2})+s\mathrm{sin}\theta (\frac{2ab}{a^2+b^2})\hfill \\ \hfill s_\theta & =s\mathrm{sin}\theta (\frac{a^2b^2}{a^2+b^2})+s\mathrm{cos}\theta (\frac{2ab}{a^2+b^2}).\hfill \end{array}$$
(32)
(Here $`\widehat{e}_r`$, $`\widehat{e}_\theta `$ and $`\widehat{e}_\varphi `$ are the spherical polar unit vectors at the position of the electron.) Evolution of the position coordinates as the particle follows the trajectory implies that the spin vector precesses about the $`z`$-axis. This was originally described by Holland .

From the above result we find that

$$𝐯_b=\frac{\mathrm{}}{2m_e}\frac{\times 𝐬}{\psi ^{}\psi }=\frac{\mathrm{}}{2mrs}\left(s_\theta +r\frac{s_\theta }{r}\frac{s_r}{\theta }\right)\widehat{e}_\varphi .$$
(33)
In other words, as was the case for $`𝐯_a`$ in (21), the contribution to the velocity from $`𝐯_b`$ is again only in the $`\widehat{e}_\varphi `$ direction. Therefore, for all eigenstates of the form in (14), the motion of the electron is in the $`\widehat{e}_\varphi `$ direction, i.e. rotational motion about the $`z`$ axis. This is in qualitative agreement with the Schrödinger results.

The total speed in the $`\widehat{e}_\varphi `$ direction is given by $`v=v_a+v_b`$, i.e.,

$$v=\frac{\mathrm{}}{m_er\mathrm{sin}\theta }\left(m+\frac{1}{2}(\frac{|v_2|^2|v_1|^2}{|v_1|^2+|v_2|^2})\right)+\frac{\mathrm{}}{2mrs}\left(s_\theta +r\frac{s_\theta }{r}\frac{s_r}{\theta }\right).$$
(34)
In what follows it will be useful to understand the relationship between the velocities for positive and negative (corresponding) values of $`m`$. Recall that for $`j=\mathrm{}+\frac{1}{2}`$,

$$\psi _{n,\mathrm{},j=\mathrm{}+\frac{1}{2},m}=\frac{1}{\sqrt{2\mathrm{}+1}}R_n\mathrm{}(r)\left(\begin{array}{c}\sqrt{\mathrm{}+m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m\frac{1}{2}}(\theta ,\varphi )\\ \sqrt{\mathrm{}m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m+\frac{1}{2}}(\theta ,\varphi )\end{array}\right)$$
(35)
and

$$Y_\mathrm{}m=\frac{1}{\sqrt{2\pi }}P_\mathrm{}me^{im\varphi },P_{\mathrm{},m}(x)=(1)^mP_\mathrm{}m(x).$$
(36)
Also, from the derivation above, the spin vector s is proportional to

$$𝐰=(2\text{Re}\{v_1^{}v_2\},2\text{Im}\{v_1^{}v_2\},|v_1|^2|v_2|^2).$$
(37)
When $`m`$ is replaced by $`m`$, we have (denoting the new term with a superscript $`()`$ and the old with $`(+)`$)

$$v_1^{()}=\sqrt{\mathrm{}m+{\scriptscriptstyle \frac{1}{2}}}Y_{\mathrm{},m\frac{1}{2}}=\sqrt{\mathrm{}m+{\scriptscriptstyle \frac{1}{2}}}\frac{1}{\sqrt{2\pi }}(1)^{m+\frac{1}{2}}P_{\mathrm{},m+\frac{1}{2}}e^{i(m\frac{1}{2})\varphi }=v_2^{(+)}.$$
(38)
Similarly,

$`v_2^{()}`$ $`=\sqrt{\mathrm{}+m+{\scriptscriptstyle \frac{1}{2}}}{\displaystyle \frac{1}{\sqrt{2\pi }}}P_{\mathrm{},m+\frac{1}{2}}e^{i(m+\frac{1}{2})\varphi }`$ (39)
$`=\sqrt{\mathrm{}+m+{\scriptscriptstyle \frac{1}{2}}}{\displaystyle \frac{1}{\sqrt{2\pi }}}(1)^{m\frac{1}{2}}P_{\mathrm{},m\frac{1}{2}}e^{i(m\frac{1}{2})\varphi }=v_1^{(+)}.`$ (40)
Therefore,

$$(|v_1|^2|v_2|^2)^{()}=(|v_1|^2|v_2|^2)^{(+)},$$
(41)
and furthermore,

$$(v_1^{}v_2)^{()}=v_2^{(+)}(v_1^{(+)})=(v_1^{}v_2)^{(+)}.$$
(42)
All three components of $`𝐰`$ change sign when $`m`$ is replaced with $`m`$ (the other eigenvalues are left unchanged), so that the spin vector in (26) changes sign. Therefore,

$$𝐯_a^{()}=\frac{\mathrm{}}{m_er\mathrm{sin}\theta }(m+\frac{1}{2}(\frac{(|v_2|^2|v_1|^2)^{()}}{(|v_1|^2+|v_2|^2)^{()}})\widehat{e}_\varphi =𝐯_a^{(+)}$$
(43)
and

$$𝐯_b^{()}=\frac{\mathrm{}}{2m_e}\frac{\times 𝐬^{()}}{\psi ^{}\psi }=𝐯_b^{(+)}.$$
(44)
Thus, both $`𝐯_a`$ and $`𝐯_b`$ change sign when $`m`$ changes sign, so that the overall velocity simply changes direction. This simplifies the computation of the rates of revolution. A similar proof holds for the case $`j=\mathrm{}\frac{1}{2}`$.

Before concluding this section, we mention that in their treatment of the Pauli equation using Euler angles, Bohm and Schiller  (p. 80) deduced that the electron in a hydrogen atom eigenstate would execute circular motion about the principal axis with constant angular velocity. However, no angular velocities were computed in the paper. In the next section we compute the angular velocities for the first few Pauli hydrogen eigenstates.

### 2.1 ANGULAR VELOCITIES FOR $`n=1`$ and $`n=2`$ PAULI EIGENSTATES

We have computed explicitly the rates of revolution $`d\varphi /dt`$ for the first few Pauli hydrogen eigenstates following the procedure described above. In each case, one computes the velocity $`𝐯_a`$, followed by the spin vector $`s=\psi ^{}\psi `$, finding $`s_r`$ and $`s_\theta `$ from (32). Then $`𝐯_b`$ is computed to give the total velocity v. Since, for all cases, v points in the $`\widehat{e}_\varphi `$ direction, the angular velocity $`d\varphi /dt`$ is given by

$$\frac{d\varphi }{dt}=\frac{v}{r\mathrm{sin}\theta }.$$
(45)
The results of our computations are presented in Table 1.

The first three results presented in Table 1 correspond to wave functions that are also eigenstates of $`s_z`$ because of the special coupling of spin and orbital angular momentum vectors. As expected, these rates of revolution agree with those computed in  for, respectively, the $`1s`$, $`2s`$ and $`2p_0`$ Schrödinger eigenstates. However, the final two states in Table 1 are not spin eigenstates. As such, they have no analogue in the Schrödinger case so that no comparisions of rates can be made.

## 3 DIRAC EIGENSTATES

We now consider the 4-component Dirac eigenstates for hydrogen. Following , they are given as follows: For $`j=\mathrm{}+\frac{1}{2}`$,

$$\begin{array}{cc}\hfill \psi _1& =g(r)\sqrt{\frac{\mathrm{}+m+\frac{1}{2}}{2\mathrm{}+1}}Y_{\mathrm{},m\frac{1}{2}}(\theta ,\varphi )\hfill \\ \hfill \psi _2& =g(r)\sqrt{\frac{\mathrm{}m+\frac{1}{2}}{2\mathrm{}+1}}Y_{\mathrm{},m+\frac{1}{2}}(\theta ,\varphi )\hfill \\ \hfill \psi _3& =if(r)\sqrt{\frac{\mathrm{}m+\frac{3}{2}}{2\mathrm{}+3}}Y_{\mathrm{}+1,m\frac{1}{2}}(\theta ,\varphi )\hfill \\ \hfill \psi _4& =if(r)\sqrt{\frac{\mathrm{}+m+\frac{3}{2}}{2\mathrm{}+3}}Y_{\mathrm{}+1,m+\frac{1}{2}}(\theta ,\varphi )\hfill \end{array}$$
(46)
and for $`j=\mathrm{}\frac{1}{2}`$,

$$\begin{array}{cc}\hfill \psi _1& =g(r)\sqrt{\frac{\mathrm{}m+\frac{1}{2}}{2\mathrm{}+1}}Y_{\mathrm{},m\frac{1}{2}}(\theta ,\varphi )\hfill \\ \hfill \psi _2& =g(r)\sqrt{\frac{\mathrm{}+m+\frac{1}{2}}{2\mathrm{}+1}}Y_{\mathrm{},m+\frac{1}{2}}(\theta ,\varphi )\hfill \\ \hfill \psi _3& =if(r)\sqrt{\frac{\mathrm{}+m\frac{1}{2}}{2\mathrm{}1}}Y_{\mathrm{}1,m\frac{1}{2}}(\theta ,\varphi )\hfill \\ \hfill \psi _4& =if(r)\sqrt{\frac{\mathrm{}m\frac{1}{2}}{2\mathrm{}1}}Y_{\mathrm{}1,m+\frac{1}{2}}(\theta ,\varphi ).\hfill \end{array}$$
(47)
The $`Y_{\mathrm{},m}`$ are the usual spherical harmonics and $`f(r)`$ and $`g(r)`$ are the normalized radial Dirac eigenfunctions (see  p. 69). Recall that even though $`\mathrm{}`$ is not a good quantum number, the eigenstates are written in terms of $`\mathrm{}`$ because it is ‘almost’ a good quantum number, and also because these solutions to the Dirac equation are built from the corresponding Pauli eigenstates. (For a complete discussion, see .)

### 3.1 TRAJECTORIES FOR GENERIC DIRAC EIGENSTATES

In this section, we show that Bohm trajectories for Dirac hydrogen share common features. First, the components of the Dirac current in (7) may be expressed in terms of the components of the wave function as follows,

$$\begin{array}{cc}\hfill \frac{1}{c}j_x& =2\text{Re}\{\psi _1^{}\psi _4\}+2\text{Re}\{\psi _2^{}\psi _3\}\hfill \\ \hfill \frac{1}{c}j_y& =2\text{Im}\{\psi _1^{}\psi _4\}2\text{Im}\{\psi _2^{}\psi _3\}\hfill \\ \hfill \frac{1}{c}j_z& =2\text{Re}\{\psi _1^{}\psi _3\}2\text{Re}\{\psi _2^{}\psi _4\}.\hfill \end{array}$$
(48)
We now compute these components using the hydrogenic wave functions given in (46) and (47). Starting with the $`z`$ component, we find that in the $`j=\mathrm{}+\frac{1}{2}`$ case,

$$\psi _1^{}\psi _3=if(r)g(r)\sqrt{\frac{\mathrm{}+m+\frac{1}{2}}{2\mathrm{}+1}}\sqrt{\frac{\mathrm{}m+\frac{3}{2}}{2\mathrm{}+3}}Y_{\mathrm{},m\frac{1}{2}}(\theta ,\varphi )^{}Y_{\mathrm{}+1,m\frac{1}{2}}(\theta ,\varphi )$$
(49)
and

$$\psi _2^{}\psi _4=if(r)g(r)\sqrt{\frac{\mathrm{}m+\frac{1}{2}}{2\mathrm{}+1}}\sqrt{\frac{\mathrm{}+m+\frac{3}{2}}{2\mathrm{}+3}}Y_{\mathrm{},m+\frac{1}{2}}(\theta ,\varphi )^{}Y_{\mathrm{}+1,m+\frac{1}{2}}(\theta ,\varphi ).$$
(50)
Both $`\psi _1^{}\psi _3`$ and $`\psi _2^{}\psi _4`$ are imaginary since the phases of the spherical harmonics cancel, implying that $`j_z=0`$. This is also the case for $`j=\mathrm{}\frac{1}{2}`$. Therefore, in all cases, motion of the electron is constrained to planes of constant $`z`$. While this is a simple result, it applies to all hydrogen eigenstates of the forms (46) and (47) and is therefore of general interest.

We find the other components of the current in a similar fashion. For the $`j=\mathrm{}+\frac{1}{2}`$ case,

$$\frac{1}{c}j_x=2\mathrm{sin}\varphi f(r)g(r)\left(\sqrt{{\scriptscriptstyle \frac{\mathrm{}m+\frac{1}{2}}{2\mathrm{}+1}}}\sqrt{{\scriptscriptstyle \frac{\mathrm{}m\frac{1}{2}}{2\mathrm{}1}}}P_{\mathrm{},m\frac{1}{2}}P_{\mathrm{}+1,m+\frac{1}{2}}+\sqrt{{\scriptscriptstyle \frac{\mathrm{}m+\frac{1}{2}}{2\mathrm{}+1}}}\sqrt{{\scriptscriptstyle \frac{\mathrm{}m+\frac{3}{2}}{2\mathrm{}+3}}}P_{\mathrm{},m+\frac{1}{2}}P_{\mathrm{}+1,m\frac{1}{2}}\right).$$
(51)
We define $`F(\mathrm{cos}\theta )`$ to be the quantity in brackets so that

$$\frac{1}{c}j_x=2\mathrm{sin}\varphi f(r)g(r)F(\mathrm{cos}\theta ).$$
(52)
Because of its similarity in form, $`j_y`$ also has the $`F(\mathrm{cos}\theta )`$ term:

$$\frac{1}{c}j_y=2\mathrm{cos}\varphi f(r)g(r)F(\mathrm{cos}\theta ).$$
(53)
From (3) and (7), the motion of the electron in a plane of constant $`z`$ is given by the following system of DEs:

$$\begin{array}{cc}\hfill \dot{x}& =\frac{j_x}{\psi ^{}\psi }=\frac{2cf(r)g(r)F(\mathrm{cos}\theta )\mathrm{sin}\varphi }{\psi ^{}\psi }\hfill \\ \hfill \dot{y}& =\frac{j_y}{\psi ^{}\psi }=\frac{2cf(r)g(r)F(\mathrm{cos}\theta )\mathrm{cos}\varphi }{\psi ^{}\psi }.\hfill \end{array}$$
(54)
From the polar forms of $`x`$ and $`y`$, we have

$$x\dot{x}+y\dot{y}=\frac{d}{dt}(x^2+y^2)=0.$$
(55)
In other words, the motion is circular about the $`z`$-axis.

A similar proof applies to the $`j=\mathrm{}\frac{1}{2}`$ case. For the sake of the computations in the next section, the components of the current for this case are $`j_z=0`$,

$$\begin{array}{cc}\hfill \frac{1}{c}j_x& =2\mathrm{sin}\varphi f(r)g(r)\left(\sqrt{{\scriptscriptstyle \frac{\mathrm{}m+\frac{1}{2}}{2\mathrm{}+1}}}\sqrt{{\scriptscriptstyle \frac{\mathrm{}m\frac{1}{2}}{2\mathrm{}1}}}P_{\mathrm{},m\frac{1}{2}}P_{\mathrm{}1,m+\frac{1}{2}}+\sqrt{{\scriptscriptstyle \frac{\mathrm{}+m+\frac{1}{2}}{2\mathrm{}+1}}}\sqrt{{\scriptscriptstyle \frac{\mathrm{}+m\frac{1}{2}}{2\mathrm{}1}}}P_{\mathrm{},m+\frac{1}{2}}P_{\mathrm{}1,m\frac{1}{2}}\right)\hfill \\ & =2\mathrm{sin}\varphi f(r)g(r)G(\mathrm{cos}\theta )\hfill \end{array}$$
(56)
and

$$\frac{1}{c}j_y=2\mathrm{cos}\varphi f(r)g(r)G(\mathrm{cos}\theta ),$$
(57)
where $`G(\mathrm{cos}\theta )`$ is the term in brackets in (56). Once again we find that motion about the $`z`$-axis is circular.

In summary, we have shown that electron trajectories associated with Dirac hydrogen eigenstates are circular, as was the case for Pauli eigenstates. In the next section, we compute some rates of revolution for these trajectories and their nonrelativistic limits.

Finally, note that if $`m`$ changes from positive to negative, both $`F(\mathrm{cos}\theta )`$ and $`G(\mathrm{cos}\theta )`$ simply undergo an overall sign change. This means that the angular rotation simply changes direction, but maintains the same functional form for $`m=\pm \frac{1}{2}`$, as was the case for the Pauli trajectories.

### 3.2 ANGULAR VELOCITIES FOR $`n=1`$ and $`n=2`$ DIRAC EIGENSTATES

We find the angular rate of revolution in general from (54) using the relation

$$\dot{y}=(r\mathrm{sin}\theta \mathrm{cos}\varphi )\dot{\varphi }.$$
(58)
From (54), it follows that

$$\frac{d\varphi }{dt}=\frac{2cf(r)g(r)F\mathrm{cos}(\theta )}{\psi ^{}\psi r\mathrm{sin}\theta \mathrm{cos}\varphi }.$$
(59)
Although this equation is deceptively simple in appearance, the functions $`f(r)`$, $`g(r)`$ – and therefore $`\psi ^{}\psi `$ – are quite complicated in form. Since we already know the qualititative motion, explicit computations of the rates of revolution for the general case, beyond the result given in (59), are not particularly enlightening.

However, it is enlightening is to examine the nonrelativistic limits of (59) and compare the results to the values computed from the Pauli equation. If the de Broglie-Bohm picture is to give a coherent account of the hydrogen atom, these results must agree. In what follows, we examine the nonrelativistic limits of (59) for the $`n=1`$ and $`n=2`$ eigenstates. Note that we compute only the positive $`m`$ value since the angular velocity simply changes sign for negative $`m`$.

#### 3.2.1 n=1

In the ground state, we have

$$\begin{array}{cc}\hfill g(r)& =\left(\frac{2}{a}\right)^{3/2}\sqrt{\frac{1+ϵ_1}{2\mathrm{\Gamma }(2\gamma _1+1)}}e^{\rho _1/2}\rho _1^{\gamma _11}\hfill \\ \hfill f(r)& =\sqrt{\frac{1ϵ_1}{1+ϵ_1}}g=\delta g,\hfill \end{array}$$
(60)
where

$$\gamma _1=\sqrt{1\alpha ^2},\rho _1=2r/a,ϵ_1=\left(1+\frac{\alpha ^2}{\gamma _1^2}\right)^{1/2},\delta =\sqrt{\frac{1ϵ_1}{1+ϵ_1}},$$
(61)
$`\alpha `$ is the fine structure constant, and $`a`$ is the Bohr radius. The ground state wave function is given by

$$\psi _1=\frac{g}{\sqrt{4\pi }},\psi _2=0,\psi _3=\frac{1}{\sqrt{4\pi }}if\mathrm{cos}\theta ,\psi _4=\frac{1}{\sqrt{4\pi }}if\mathrm{sin}\theta e^{i\varphi }.$$
(62)
This gives

$$\rho =\psi ^{}\psi =\frac{1}{4\pi }(1+\delta ^2)g^2\text{and}F(\mathrm{cos}\theta )=\frac{1}{4\pi }\mathrm{sin}\theta .$$
(63)
Substitution into (59) yields, after cancellation,

$$\frac{d\varphi }{dt}=\left(\frac{2}{r}\right)\frac{\delta c}{1+\delta ^2}.$$
(64)
In the nonrelativistic limit, $`c\mathrm{}`$, which implies that $`\alpha =e^2/\mathrm{}c0`$ and $`\gamma _11`$. Furthermore, this implies that $`ϵ_11`$ and $`\delta 0`$. In order to determine the behaviour of $`\delta c`$, we expand $`ϵ_1`$ as

$$ϵ_11\frac{1}{2}(\frac{\alpha ^2}{\gamma _1^2})$$
(65)
so that

$$\frac{1ϵ_1}{1+ϵ_1}\frac{1}{4}\alpha ^2\text{as}c\mathrm{}.$$
(66)
Therefore $`\delta \frac{1}{2}\alpha `$. Substitution into (64) yields

$$\frac{d\varphi }{dt}\frac{1}{r}\alpha c=\frac{e^2}{r\mathrm{}},$$
(67)
which, when written in terms of the Bohr radius $`a`$, becomes

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{m_ear}.$$
(68)
This is the angular rotation rate for the ground state Dirac wave function. It is in agreement with the Schrödinger rate given by Holland , and it is also in agreement with the rate found for the Pauli equation in Table 1.

#### 3.2.2 n=2

1. $`2S_{1/2}`$ ($`n=2,\mathrm{}=0,j=\frac{1}{2},m=\frac{1}{2}`$)

As in the $`1s`$ case, we have the wave function

$$\psi _1=\frac{g}{\sqrt{4\pi }},\psi _2=0,\psi _3=\frac{1}{\sqrt{4\pi }}if\mathrm{cos}\theta ,\psi _4=\frac{1}{\sqrt{4\pi }}if\mathrm{sin}\theta e^{i\varphi },$$
(69)
where the functions $`f(r)`$ and $`g(r)`$ are suitably modified for the $`n=2`$ case. Again, their exact functional form is not relevant, but the relationship between $`f`$ and $`g`$ is important; here, rather than (61) we have

$$\rho _2=\frac{2r}{N_2a},N_2=\sqrt{2(1+\gamma _1)},ϵ_2=\left(1+(\frac{\alpha }{1+\gamma _1})^2\right)^{1/2},\delta =\sqrt{\frac{1ϵ_2}{1+ϵ_2}}A,$$
(70)
and the number $`A`$ is given by

$$A=\frac{(2\gamma _1+1)(N_2+2)(N_2+1)\rho _2}{(2\gamma _1+1)N_2(N_2+1)\rho _2}.$$
(71)
After some cancellation, substitution of the wave function into (59) gives, as in (64),

$$\frac{d\varphi }{dt}=\left(\frac{2}{r}\right)\frac{\delta c}{1+\delta ^2}.$$
(72)
Once again, we examine how the quantities in (70) behave in the nonrelativistic limit $`c\mathrm{}`$. In this case,

$$ϵ_21\frac{1}{2}(\frac{\alpha }{1+\gamma _1})^2$$
(73)
so that

$$\sqrt{\frac{1ϵ_2}{1+ϵ_2}}\frac{\alpha }{4}.$$
(74)
From the properties $`\gamma _11`$, hence $`N_22`$, the limit of $`A`$ in (70) is

$$A\frac{4\rho _2}{2\rho _2}.$$
(75)
This implies that

$$\delta c\alpha c\frac{4\rho _2}{4(2\rho _2)}$$
(76)
and (72) becomes

$$\frac{d\varphi }{dt}=\frac{e^2(4\rho _2)}{2\mathrm{}r(2\rho _2)},$$
(77)
which can be rewritten as

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{2m_ear}\left(1\frac{1}{\frac{r}{2a}1}\right).$$
(78)
This is the angular rotation rate for the $`2s`$ Schrödinger state given in  and is also in agreement with the Pauli result of Table 1.
2. $`2P_{3/2}`$ ($`n=2,\mathrm{}=1,j=\frac{3}{2},m=\frac{3}{2}`$)

In this case the wave function is given by

$$\psi _1=\sqrt{\frac{3}{8\pi }}g\mathrm{sin}\theta e^{i\varphi },\psi _2=0,\psi _3=if\sqrt{\frac{3}{8\pi }}\mathrm{cos}\theta \mathrm{sin}\theta e^{i\varphi },\psi _4=if\sqrt{\frac{3}{8\pi }}\mathrm{sin}^2\theta e^{3i\varphi }.$$
(79)
The functions $`f`$ and $`g`$, and the relationship between them, will be the same as in the above case, because only $`m`$ has changed. However, we now have $`\psi ^{}\psi =\frac{3}{8\pi }\mathrm{sin}^2\theta g^2`$ in the limit as $`\delta 0`$. Therefore, the expression corresponding to (64) is

$$\frac{d\varphi }{dt}=\frac{2c\delta }{r}.$$
(80)
Substitution of the nonrelativistic limiting expressions from the previous case yields

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{2m_ear}.$$
(81)
This is the angular rotation rate for the $`2p_1`$ Schrödinger state given in  and is also in agreement with the Pauli result in Table 1.
3. $`2P_{1/2}`$ ($`n=2,\mathrm{}=1,j=\frac{1}{2},m=\frac{1}{2}`$)

This case is similar to the previous one although the functional forms of $`f`$ and $`g`$ are different. We have

$$\psi _1=\frac{1}{\sqrt{4\pi }}g\mathrm{cos}\theta ,\psi _2=\frac{1}{\sqrt{4\pi }}g\mathrm{sin}\theta e^{i\varphi },\psi _3=if\frac{1}{\sqrt{4\pi }},\psi _4=0$$
(82)
and

$$\frac{d\varphi }{dt}=\left(\frac{2}{r}\right)\frac{\delta c}{1+\delta ^2}$$
(83)
and most of the definitions of (70) remain the same. In this case, the term $`A`$ is given by

$$A=\frac{(2\gamma _1+1)N_2(N_21)\rho _2}{(2\gamma _1+1)(N_22)(N_21)\rho _2}.$$
(84)
In the non-relativistic limit,

$$A\frac{6\rho _2}{\rho _2}$$
(85)
so that

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{m_er^2}\left(3\frac{r}{2a}\right).$$
(86)
In this case, no comparison can be made with any Schrödinger state.
4. $`2P_{3/2}`$ ($`n=2,\mathrm{}=1,j=\frac{3}{2},m=\frac{1}{2}`$)

This case is somewhat different. Here, the wave function is

$$\begin{array}{cc}\hfill \psi _1=\frac{1}{\sqrt{2\pi }}g\mathrm{cos}\theta ,& \psi _2=\frac{1}{\sqrt{8\pi }}g\mathrm{sin}\theta e^{i\varphi }\hfill \\ \hfill \psi _3=if\frac{1}{\sqrt{8\pi }}(3\mathrm{cos}^2\theta 1),& \psi _4=if\sqrt{\frac{9}{8\pi }}\mathrm{sin}\theta \mathrm{cos}\theta e^{i\varphi }\hfill \end{array}$$
(87)
and

$$\psi ^{}\psi =\frac{1}{8\pi }(4\mathrm{cos}^2\theta +\mathrm{sin}^2\theta )g^2+\frac{1}{8\pi }((3\mathrm{cos}^2\theta 1)^2+3\mathrm{sin}^2\theta \mathrm{cos}^2\theta )f^2.$$
(88)
Furthermore, the function $`F(cos\theta )`$ becomes

$$F(cos\theta )=\frac{1}{8\pi }\mathrm{sin}\theta (8\mathrm{cos}^2\theta \mathrm{sin}^2\theta ).$$
(89)
Again, $`f=\delta g`$, with

$$\delta =\sqrt{\frac{1ϵ_3}{1+ϵ_3}},ϵ_3=\left(1+\frac{\alpha ^2}{\gamma _2^2}\right)^{1/2},\gamma _2=\sqrt{4\alpha ^2}$$
(90)
and we find from (59) that

$$\frac{d\varphi }{dt}=\frac{2\delta c(8\mathrm{cos}^2\theta \mathrm{sin}^2\theta )}{r(4\mathrm{cos}^2\theta +\mathrm{sin}^2\theta )}.$$
(91)
As $`c\mathrm{}`$,

$$ϵ_31\frac{1}{2}(\frac{\alpha }{2})^2$$
(92)
so that $`\delta c\frac{\alpha }{4}`$. After substituting and rewriting $`\alpha `$, we obtain

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{2m_ear}\left(\frac{8\mathrm{cos}^2\theta \mathrm{sin}^2\theta }{4\mathrm{cos}^2\theta +\mathrm{sin}^2\theta }\right).$$
(93)
Once again, no comparison can be made with any Schrödinger state.

In each case presented above, the nonrelativistic limit of the Dirac angular velocity agrees with the corresponding Pauli result given in Table 1. We expect this, since the Pauli equation is the nonrelativistic limit of the Dirac equation. However, the results are not obvious, since the expressions in (7) and (10) for, respectively, the Dirac and Pauli currents are quite different.

## 4 CONCLUDING REMARKS

In this paper, we have determined the general features of de Broglie-Bohm trajectories for energy/total angular momentum eigenstates of the Pauli and Dirac hamiltonians for hydrogen. In all cases, the electron, assumed to be in an eigenstate of $`M_z`$, the $`z`$-component of the total angular momentum, $`M`$, is confined to a plane of constant $`z`$-value and executes circular motion about the $`z`$-axis with a constant angular velocity $`d\varphi /dt`$. As well, we have outlined a procedure to compute these angular velocities for general eigenstates and have explicitly computed them for the $`n=1`$ and $`n=2`$ Pauli and Dirac hydrogen eigenstates.

In the cases where the Pauli eigenstates are also eigenstates of the $`s_z`$ operator, our results from the Pauli equation agree with earlier computations of the trajectories of corresponding Schrödinger eigenstates . Furthermore, the nonrelativistic limits of the Dirac results agree with the Pauli results. We have therefore shown that the de Broglie-Bohm causal picture can be applied coherently to the hydrogen atom, moving from the Schrödinger to the Pauli and ultimately to the Dirac equation.

Finally, one may well wish to consider trajectories for Pauli or Dirac wave functions other than those considered in this paper. For example, it may be interesting to examine trajectories for particular linear combinations of eigenstates. (In , we examined Bohm trajectories for the familiar Schrödinger $`2p_x`$ and $`2p_y`$ orbitals used in descriptions of chemical bonding. As well, we examined trajectories associated with a time-varying linear combination of $`1s`$ and $`2p_0`$ hydrogenic wave functions that simulated an electronic transition induced by an oscillating electric field.) The method of computing Bohm trajectories outlined in Sections 2 and 3 can be extended in a straightforward manner to treat such linear combinations, although the computations may well become quite complicated.

## ACKNOWLEDGEMENTS

We gratefully acknowledge that this research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) in the form of a Postgraduate Scholarship (CC) and an Grant in Aid of Research (ERV). CC also acknowledges partial financial support from the Province of Ontario (Graduate Scholarship) as well as the Faculty of Mathematics, University of Waterloo.


format: text

text:
# Purity of Gaussian states: measurement schemes and time–evolution in noisy channels

## I Introduction

Nonclassical features of atomic and radiation systems play a relevant role in quantum information, communication and high precision measurements, as well as in many fundamental experiments to test quantum mechanics pati ; heiss . In particular, pure Gaussian states of continuous variable (CV) systems, such as coherent and squeezed–coherent states, are the key ingredients of secure optical communication yuk ; yue ; gran ; ralph  and Heisenberg limited quantum interferometry cav ; bon ; yur ; ind ; int . The characterization of several properties of Gaussian states has been the subject of intense recent work kim ; ban ; eisert ; bin ; fil ; rad ; hol , stimulated by the seminal analysis on their entanglement properties simon ; cirac .

Any attempt to exploit Gaussian states for quantum information and quantum measurement schemes must however face the obvious difficulty that pure states are unavoidably corrupted by the interaction with the environment. Therefore, CV Gaussian states that are available for experiments are usually mixed states, and it becomes crucial to establish their degree of purity (or mixedness) determined by the environmental noise. In the present paper, we study the purity of Gaussian states for single–mode continuous variable systems focusing on two aspects: its experimental characterization, and its time–evolution in noisy channels. We first show that the joint detection of two conjugate quadratures is a necessary and sufficient measurement to determine the purity of a Gaussian state with reliable experimental statistics; we then derive an evolution equation for the purity of Gaussian states in a noisy channel, considering the instances of a thermal bath and of a squeezed thermal bath, and determine the evolutions that at any given time maximize the purity.

Let us refer to $`\mu =\text{Tr}\left[\varrho ^2\right]`$ as to the purity of a generic quantum state $`\varrho `$; the conjugate quantity $`S_L=(1\mu )d/(d1)`$, where $`d`$ is the dimension of the Hilbert space of the system under investigation, is known as linear entropy or mixedness. In general, $`\mu `$ ranges from one, which is the value for a pure state, to $`\mu =1/d`$ for a completely mixed state. Throughout the paper we will consider CV systems, i.e. infinite dimensional Hilbert spaces, and therefore we will have $`0<\mu 1`$. Since $`\mu `$ is a nonlinear function of the density matrix, it cannot be connected to an observable quantity if we perform repeated measurements on single copies of the state. That is, it cannot be the expectation value of a single-system self-adjoint operator, nor it can be related to a single-system probability distribution obtained from a positive operator-valued measure (POVM). On the other hand, if collective measurements on two copies of the state are possible, then the purity may be measured directly art ; fil2 . For instance, collective measurements of overlap and fidelity have been experimentally realized for qubits encoded into polarization states of photons hey .

In general, purity can be determined by the knowledge of the quantum state of the system, which in turn can be obtained by quantum tomography review . However, in this case, the statistics is usually poor, since the measurement of a whole quorum of observables is needed, unavoidably leading to large fluctuations add . On the other hand, if we focus our attention on the class of Gaussian states, it is indeed possible to find an operative method to experimentally determine $`\mu `$ . In fact, Gaussian states are uniquely defined by their first two statistical moments, which can be measured by the joint detection of two conjugate quadratures, say position and momentum or quadrature phases of the electromagnetic field. Such a measurement corresponds to an estimate of the Husimi $`Q`$–function $`Q(\alpha )=\alpha |\varrho |\alpha `$, $`|\alpha `$ being a coherent state of the harmonic oscillator. We will show that the measurement of the $`Q`$–function is the optimal minimal measurement for the purity, in the sense that it is necessary and sufficient to determine $`\mu `$ and requires the minimum number of observables to be measured.

The joint measurement of two conjugate quadratures is possible for a single–mode radiation field as well as for a single atom kelly ; ste ; wil . Remarkably, for these systems, the class of Gaussian states include almost all the states that can be reliably produced, and employed in communication or measurement protocols.

Finally, we will show that the previous discussion allows to unravel the dynamics of purity only in terms of observable quantities. Indeed, the time–evolution of the purity of an initial Gaussian state in a noisy channel can be uniquely expressed as a function of the initial observable parameters of the input state and of the asymptotic observable parameters of the environment. This property allows then to determine and engineer optimal evolutions, i.e. evolutions that maximize the purity at any given time.

The paper is structured as follows. In Section II we show how purity is related to observable quantities for Gaussian states, and how it can be obtained either from the $`Q`$–function or by three single-quadrature detection. In Section III we present the results of a systematic numerical analysis, establishing the statistical reliability and the range of applicability of the method by means of Monte Carlo simulated experimental runs. We also show that the $`Q`$-function based determination of purity is a more reliable method than single-quadrature detection. Section IV is devoted to derive and solve an evolution equation for the purity of an initial Gaussian state in a noisy channel, both for thermal and squeezed thermal baths. We show that, even though the asymptotic value of purity is not related to the initial conditions, its behavior at finite times does depend on the initial squeezing and thermal excitations, and we determine the evolutions that maximize the purity at any finite time. We show in particular that purity is maximized for an initial coherent state evolving in a thermal bath, or for an initial squeezed state evolving in a squeezed thermal bath whose asymptotic squeezing is orthogonal to that of the input state. Finally, in Section V we present some concluding remarks.

## II Purity of Gaussian states

We begin by reviewing some fundamental properties of the Wigner phase-space representation barnett  which will be useful throughout the paper. The Wigner representation of an arbitrary operator $`O`$ is defined as follows

$$O(\alpha )=_{}\frac{\mathrm{d}^2\gamma }{\pi ^2}\mathrm{e}^{\overline{\gamma }\alpha \gamma \overline{\alpha }}\mathrm{Tr}[OD(\gamma )],$$
(1)
where $`D(\gamma )=\mathrm{exp}(\gamma a^{}\overline{\gamma }a)`$ is the displacement operator, and $`\mathrm{Tr}[OD(\gamma )]`$ is usually referred to as the characteristic function of the operator $`O`$. Let $`O_1`$ and $`O_2`$ be operators that admit regular Wigner representations $`O_1(\alpha )`$ and $`O_2(\alpha )`$. Then the trace $`\text{Tr}\left[O_1O_2\right]`$ can be computed as an integral over phase space according to

$$\text{Tr}[O_1O_2]=\pi _{}\mathrm{d}^2\alpha O_1(\alpha )O_2(\alpha ).$$
(2)
From now on, we will move to the phase–space variables $`x`$ and $`p`$, corresponding to quadrature phases $`\widehat{x}=(a+a^{})/\sqrt{2}`$ and $`\widehat{p}=i(a^{}a)/\sqrt{2}`$ of the field $`a`$, whose expectation values $`\widehat{x}x`$ and $`\widehat{p}p`$ are related to $`\alpha `$ by $`\alpha =(x+ip)/\sqrt{2}`$.

The Wigner representation $`W(\alpha )`$ of the density matrix $`\varrho `$ of a quantum state is referred to as the Wigner function of the state. The class of Gaussian states is defined as the class of states with Gaussian Wigner function, namely

$$W(x,p)=\frac{\mathrm{e}^{\frac{1}{2}X𝝈^1X^T}}{\pi \sqrt{\mathrm{Det}[𝝈]}},$$
(3)
where $`X`$ is the displaced vector $`X=(xx_0,pp_0)`$ and $`𝝈`$ is the covariance matrix

$$\sigma _{ij}=\frac{1}{2}\widehat{x}_i\widehat{x}_j+\widehat{x}_j\widehat{x}_i\widehat{x}_i\widehat{x}_j,$$

where $`\widehat{x}_1=\widehat{x},\widehat{x}_2=\widehat{p}`$. The density matrix of the most general Gaussian state can be written as adam

$`\varrho =D(\alpha _0)S(r,\phi )\nu _{\overline{n}}S^{}(r,\phi )D^{}(\alpha _0),`$ (4)
where $`\alpha _0=(x_0+ip_0)/\sqrt{2}`$, and $`\nu _{\overline{n}}`$ is a thermal state with average photon number $`\overline{n}`$:

$$\nu _{\overline{n}}=\frac{1}{1+\overline{n}}\underset{k=0}{\overset{\mathrm{}}{}}\left(\frac{\overline{n}}{1+\overline{n}}\right)^k|kk|,$$

$`D(\alpha _0)`$ denotes the displacement operator and $`S(r,\phi )=\mathrm{exp}(\frac{1}{2}r\mathrm{e}^{i2\phi }a^2\frac{1}{2}r\mathrm{e}^{i2\phi }a^2)`$ the squeezing operator. A convenient parametrization of Gaussian states can be achieved replacing the $`\sigma _{ij}`$’s by $`n`$, $`r`$, $`\phi `$, which have a more direct phenomenological interpretation. By applying the phase-space representation of squeezing barnett ; walls , the following relations are easily derived

$`\sigma _{xx}`$ $`=`$ $`{\displaystyle \frac{2\overline{n}+1}{2}}[\mathrm{cosh}(2r)\mathrm{sinh}(2r)\mathrm{cos}(2\phi )],`$
$`\sigma _{pp}`$ $`=`$ $`{\displaystyle \frac{2\overline{n}+1}{2}}[\mathrm{cosh}(2r)+\mathrm{sinh}(2r)\mathrm{cos}(2\phi )],`$
$`\sigma _{xp}`$ $`=`$ $`{\displaystyle \frac{2\overline{n}+1}{2}}\mathrm{sinh}(2r)\mathrm{sin}(2\phi ).`$ (5)
Exploiting Eq. (2), one can write

$$\mu \mathrm{Tr}[\varrho ^2]=\frac{\pi }{2}_{}_{}dxdpW^2(x,p),$$
(6)
so that, for a Gaussian state

$$\mu =\frac{1}{2\sqrt{\mathrm{Det}[𝝈]}}=\frac{1}{2\sqrt{\sigma _{xx}\sigma _{pp}\sigma _{xp}^2}}.$$
(7)
In terms of $`\overline{n}`$, $`r`$ and $`\phi `$, Eq. (7) can then be recast as marian1 ; dodonov

$$\mu =\frac{1}{2\overline{n}+1}.$$
(8)
Eq. (8) shows that the purity of a generic Gaussian state depends only on the average number of thermal photons, as one should expect since displacement and squeezing are unitary operations. Therefore, the measurement of the purity of a Gaussian state is equivalent to the measurement of its average number of thermal photons.

As the last step in connecting $`\mu `$ to observables we report the expression of the $`\sigma _{ij}`$’s in terms of the $`Q`$–function $`Q(\alpha )`$. This follows from the antinormally ordered expression of the second moments. We have, for instance

$$\widehat{x}^2=\frac{a^2+a^2+2aa^{}𝕀}{2},$$

which, in terms of phase–space variables, corresponds to $`x^2\frac{1}{2}`$. Therefore, we eventually get

$`\widehat{x}^2`$ $`=`$ $`\mathrm{Tr}[\varrho {\displaystyle \frac{(a+a^{})^2}{2}}]`$
$`=`$ $`{\displaystyle _{}}{\displaystyle _{}}dxdpQ(x,p)(x^2{\displaystyle \frac{1}{2}}),`$
where we have moved from variables $`\alpha `$ and $`\overline{\alpha }`$ to variables $`x`$ and $`p`$, previously defined. In much the same way, we obtain

$`\widehat{p}^2`$ $`=`$ $`{\displaystyle _{}}{\displaystyle _{}}dxdpQ(x,p)(p^2{\displaystyle \frac{1}{2}}),`$ (9)
$`{\displaystyle \frac{1}{2}}\widehat{x}\widehat{p}+\widehat{p}\widehat{x}`$ $`=`$ $`{\displaystyle _{}}{\displaystyle _{}}dxdpQ(x,p)xp.`$ (10)
Since first moments are naturally antinormally ordered, evaluation of first moments of quadratures is easily obtained, and the $`\sigma _{ij}`$’s can be eventually computed.

Gaussian states may be effectively characterized as well by single-quadrature measurements obtained by balanced homodyne detection parlik . Thus a question arises whether or not one really needs to resort to joint measurement of two conjugate quadratures to determine the purity. In particular, since Gaussian states are fully characterized by the first and second moments, it suffices to measure the rotated quadrature $`x_\theta =(a^{}\mathrm{e}^{i\theta }+a\mathrm{e}^{i\theta })/\sqrt{2}`$ for three different values of $`\theta `$ to have a complete characterization of the state, including the measure of its purity. This fact can be proven by reminding that the probability distribution $`p(x,\theta )`$ of a measurement of $`x_\theta `$ on a state of the form (4) is a Gaussian centered in $`x_0=\text{Re}[\alpha _0\mathrm{e}^{i\theta }]`$, with variance

$$\sigma _\theta =\frac{1}{2\mu }\left[\mathrm{e}^{2r}\mathrm{cos}^2(\theta \phi )+\mathrm{e}^{2r}\mathrm{sin}^2(\theta \phi )\right].$$
(11)
By measuring three quadratures we directly obtain the purity $`\mu `$ by comparison of variances. By choosing $`\theta =0,\pi /2,\pi /4`$ we have

$`\mu =\left[4\sigma _{\pi /4}(\sigma _0+\sigma _{\pi /2}\sigma _{\pi /4})(\sigma _0\sigma _{\pi /2})^2\right]^{\frac{1}{2}}.`$ (12)
In the next Section, we will compare the two different experimental schemes on the basis of Monte Carlo simulated experiments.

## III Monte Carlo simulated experiments

As we have seen, in order to evaluate the $`\sigma _{ij}`$’s and then the purity, we need to estimate averages over the $`Q`$–function. These estimates can be obtained if one disposes of data distributed according to the $`Q`$–function $`Q(x,p)`$ itself. Indeed, such a distribution can be experimentally reconstructed for a single–mode radiation field through heterodyne sha , eight–port homodyne ott ; leo  or six–port homodyne detectors tri , and for atoms by coupling the atom with two light fields and measuring the corresponding phase–shifts wil .

In order to test the effectiveness of the proposed scheme, we have performed a systematic numerical analysis by means of Monte Carlo simulated experiments. The simulations are needed to show the actual independence of the method on the squeezing and displacing parameters, in compliance with Eq. (8). Moreover, they provide a crucial test on the actual possibility of getting reliable (i.e. with reduced fluctuations) determinations of $`\mu `$ in realistic experimental settings and even for most unfavorable states.

The purity $`\mu `$ and its dispersion $`\mathrm{\Delta }\mu `$ have been evaluated from samples of the $`Q`$–function, varying the values of the parameters of the simulated Gaussian state. Besides $`\overline{n}`$, $`r`$, $`\phi `$, and $`\alpha _0`$, the experimental determination of the $`Q`$–function depends on the number $`N_x`$ of collected data.

We find that $`\mu `$ and $`\mathrm{\Delta }\mu `$ are essentially independent on the complex displacement parameter $`\alpha _0`$ and on the squeezing angle $`\phi `$. On the other hand, $`\mathrm{\Delta }\mu `$ does depend on $`\overline{n}`$ and $`r`$, decreasing with increasing $`\overline{n}`$ and increasing with increasing $`r`$.

In Fig. 1 we report the determination of purity for a strongly squeezed thermal state as a function of the number of data. The error on purity is of the order of a few percent for samples made of $`N_x10^5`$ data.

In order to compare the determination of $`\mu `$ by the $`Q`$–function with that coming from single–quadrature detection, we have simulated the measurement of three quadratures $`x_\theta `$, $`\theta =0,\pi /2,\pi /4`$ by balanced homodyne detection. In Fig. 2 we report the estimated purity \[using Eq. (12)\] for the same strongly squeezed thermal state of Fig. 1 as a function of the number of data. Some features are immediately evident. First of all one can see that the determination is biased: in the present case the estimated $`\mu `$ is always larger than the true value, while the opposite case occurs by inverting the phase of the squeezing. Therefore the method is very sensitive to the choice of the phase. Moreover, the relative error is not a smooth function of the number of data i.e. the method is not statistically reliable as the joint–measurement one. This is again due to the remarkable dependence of the variances on the phase of the squeezing, a dependence which is instead smoothed out in the measurement of the $`Q`$–function. Summing up, for some specific states (as the example considered here) single-quadrature detection may be asymptotically even more efficient than the heterodyne one. However, in general, the number of data needed for the relative error to be below the joint–measurement level is strongly state–dependent. We conclude that the measurement of the $`Q`$–function is statistically more reliable and thus more suited for a systematic analysis of the purity of Gaussian states.

Let us now go back to the analysis of the $`Q`$–function determination of purity. A smaller number of data is needed to obtain a given precision for states with smaller squeezing. The effect of the squeezing parameter on the determination of purity is illustrated in Fig. 3, where we report $`\mu `$ and $`\mathrm{\Delta }\mu /\mu `$ versus $`r`$ for Gaussian states with $`\alpha =0`$, $`\phi =0`$, and $`\overline{n}=0.5`$, and for a number of data $`N_x=310^4`$. Notice that in Fig. 3 the range of $`r`$ corresponds to a quite large number of mean squeezed photons $`0\mathrm{sinh}^2r15`$.

In the deep quantum regime, i.e. for small $`\overline{n}`$, fluctuations of $`\mu `$ become more relevant. This is not surprising, since $`\mu `$ is a highly nonlinear function of the second–order moments. However, simulations show that even for highly squeezed (up to $`15`$ mean squeezed photons) and slightly mixed (down to $`\overline{n}0.1`$) states, realistic experimental conditions allow a statistically reliable determination of $`\mu `$ that complies with the theoretical expectation (8), up to an error of a few percent. In Fig. 4 we plot the determination of purity for different squeezed thermal states as a function of the average number of thermal photons $`\overline{n}`$, for samples made of $`N_x=10^5`$ data.

From the above analysis we conclude that the joint measurement of two conjugate quadratures provides a statistically reliable method to determine the purity of a generic Gaussian state. This is best achieved with experimental schemes that involve data distributed according to the Husimi $`Q`$–function, such as heterodyne and multi–port homodyne detection schemes.

## IV Evolution of purity in a noisy channel

Let us consider the time evolution of an initial, pure or mixed, generic single–mode Gaussian state in presence of noise and damping (and/or pumping) toward a final squeezed thermal state. If $`\mathrm{\Gamma }^1`$ is the photon lifetime in the noisy channel, the evolution of a state is described, in the interaction picture, by the following master equation

$`\dot{\varrho }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}NL[a^{}]\varrho +{\displaystyle \frac{\mathrm{\Gamma }}{2}}(N+1)L[a]\varrho `$ (13)
$``$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\left(\overline{M}D[a]\varrho +MD[a^{}]\varrho \right),`$
where the dot stands for time–derivative and the Lindblad superoperators are defined by

$`L[O]\varrho `$ $``$ $`2O\varrho O^{}O^{}O\varrho \varrho O^{}O,`$ (14)
$`D[O]\varrho `$ $``$ $`2O\varrho OOO\varrho \varrho OO.`$ (15)
$`M`$ is the correlation function of the bath (which is usually referred to as the squeezing of the bath); it is in general a complex number $`M=M_1+iM_2`$, and $`\overline{M}`$ denotes its complex conjugate, while $`N`$ is a phenomenological parameter related, as we shall see, to the purity of the asymptotic state. Positivity of the density matrix imposes the important constraint $`|M|^2N(N+1)`$. At thermal equilibrium, i.e. for $`M=0`$, $`N`$ coincides with the average number of thermal photons in the bath. The master equation (13) can be transformed into a Fokker–Planck equation for the Wigner function $`W(x,p,t)`$. Using the differential representation of the superoperators walls ; paris  in Eq. (13), the corresponding Fokker-Planck equation reads as follows

$`\dot{W}(x,p,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}(2+x_x+p_p+{\displaystyle \frac{2N+1}{2}}(_{xx}^2+_{pp}^2)`$ (16)
$`+`$ $`M_1(_{xx}^2_{pp}^2)+2M_2_{xp})W(x,p,t).`$
For a general single–mode Gaussian state of the form (3) one has, in compact notation

$`\dot{W}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}(2X𝝈^1\left({\displaystyle \genfrac{}{}{0pt}{}{x}{p}}\right)+{\displaystyle \frac{2N+1}{2}}𝒮_{X𝝈}(𝕀)`$ (17)
$`+`$ $`M_1𝒮_{X𝝈}(𝔸)+M_2𝒮_{X𝝈}(𝔹))W,`$
where $`𝒮_{X𝝈}(𝜸)`$ denotes the seralian (or siralian) operator, a scalar function of the matrix $`𝜸`$ given by $`𝒮_{X𝝈}(𝜸)X𝝈^1𝜸𝝈^1X^T\mathrm{Tr}[𝜸𝝈^1]`$. The displaced vector $`X`$ and the covariance matrix $`𝝈`$ have been previously defined, whereas $`𝕀,𝔸,𝔹`$ form a basis in the space of $`2\times 2`$ real symmetric matrices:

$$𝕀=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),𝔸=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),𝔹=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$

For any given real matrix $`𝜸`$ and generic Gaussian states, the seralian operator shows the remarkable property:

$$_{}dx_{}dp𝒮_{X𝝈}(𝜸)W(x,p,t)=0.$$
(18)
It can be easily shown that this property assures that the last three terms of Eq. (17) \[corresponding to diffusion terms in the Fokker–Planck equation (16)\] do not enter in the time–evolution equations for the first statistical moments $`x_0`$ and $`p_0`$. Such evolution is governed by the drift terms and is described by the following equation for the vector $`X_0\left(\genfrac{}{}{0pt}{}{x_0}{p_0}\right)`$

$`\dot{X}_0={\displaystyle _{}}dx{\displaystyle _{}}dp\left({\displaystyle \genfrac{}{}{0pt}{}{x}{p}}\right){\displaystyle \frac{\dot{W}}{2}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}X_0.`$ (19)
First moments are damped through the noisy channel: this effect should be expected since it is the mathematical evidence of the absorption of the state’s coherent photons.

The evolution of the covariance matrix of the state can be described by monitoring different sets of variables. A good choice of variables is given by the $`\sigma _{ij}`$’s, in terms of which the evolution equations decouple. The relations between the variables $`\sigma _{ij}`$ and the variables $`\mu `$, $`r`$, and $`\phi `$ are given in Eqns. (5) and (7). Here we recall some further relations that will be useful in the following:

$$\mathrm{Det}[𝝈]=\sigma _{xx}\sigma _{pp}\sigma _{xp}^2=\frac{(2\overline{n}+1)^2}{4}=\frac{1}{4\mu ^2},$$
(20)
$$\mathrm{Tr}[𝝈]=\sigma _{xx}+\sigma _{pp}=(2\overline{n}+1)\mathrm{cosh}(2r)=\frac{\mathrm{cosh}(2r)}{\mu },$$
(21)
$$\sigma _{pp}\sigma _{xx}=\frac{\mathrm{sinh}(2r)\mathrm{cos}(2\phi )}{\mu }.$$
(22)
As we have seen, in the Wigner phase–space picture the expectation values can be computed as phase–space integrals. The first–order evolution equation for the covariance matrix $`𝝈`$ is thus obtained by straightforward integration, and reads:

$`\dot{𝝈}`$ $`=`$ $`\mathrm{\Gamma }\left(𝝈_{\mathrm{}}𝝈\right),`$ (23)
$`\mathrm{with}𝝈_{\mathrm{}}`$ $``$ $`\left(\begin{array}{cc}\frac{(2N+1)+2M_1}{2}& M_2\\ & \\ M_2& \frac{(2N+1)2M_1}{2}\end{array}\right).`$ (27)
The matrix $`𝝈_{\mathrm{}}`$, determined by the bath parameters alone, turns out to be the asymptotic covariance matrix. In fact, integration of Eq. (23) yields

$$𝝈(t)=𝝈_{\mathrm{}}\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)+𝝈(0)\mathrm{e}^{\mathrm{\Gamma }t}.$$
(28)
Eq. (28) shows a simple example of a Gaussian completely positive map eisplen . The Gaussian character of the evoluted Wigner function can be proven, a posteriori, by verifying that a function of the form (3), with covariance matrix given by Eq. (28), indeed solves Eq. (16). In order to be a bona fide covariance matrix, $`𝝈(t)`$ must satisfy the usual condition encoding the $`\widehat{x}\widehat{p}`$ uncertainty relations eisplen ; simsud

$$𝝈(t)+\frac{i}{2}𝐉0,\mathrm{with}𝐉=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$
(29)
It is promptly seen that such a condition is satisfied at any time by the convex combination giving $`𝝈(t)`$ in Eq. (28) iff $`𝝈_{\mathrm{}}`$ is a legitimate covariance matrix. This last requirement is assured by the necessary constraint $`N(N+1)|M|^2`$ that guarantees positivity of the density matrix.

By introducing

$$\mu _{\mathrm{}}\left[(2N+1)^24|M|^2\right]^{1/2},$$

and exploiting Eqns. (2022) we can eventually express $`\mu `$, $`r`$ and $`\phi `$ as functions of time

$`\mu (t)`$ $`=`$ $`\mu _0[{\displaystyle \frac{\mu _0^2}{\mu _{\mathrm{}}^2}}(1\mathrm{e}^{\mathrm{\Gamma }t})^2+\mathrm{e}^{2\mathrm{\Gamma }t}`$ (30)
$`+`$ $`2\mu _0({\displaystyle \frac{\sqrt{1+4\mu _{\mathrm{}}^2|M|^2}\mathrm{cosh}(2r_0)}{\mu _{\mathrm{}}}}+2\mathrm{sinh}(2r_0)(M_1\mathrm{cos}(2\phi _0)M_2\mathrm{sin}(2\phi _0)))(1\mathrm{e}^{\mathrm{\Gamma }t})\mathrm{e}^{\mathrm{\Gamma }t}]^{1/2},`$
$$\mathrm{cosh}[2r(t)]=\mu (t)\left(\frac{\sqrt{1+4\mu _{\mathrm{}}^2|M|^2}\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)}{\mu _{\mathrm{}}}+\mathrm{e}^{\mathrm{\Gamma }t}\frac{\mathrm{cosh}(2r_0)}{\mu _0}\right),$$
(31)
$$\mathrm{tan}[2\phi (t)]=\frac{M_22\mu _0\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)+\mathrm{sinh}(2r_0)\mathrm{sin}(2\phi _0)\mathrm{e}^{\mathrm{\Gamma }t}}{M_12\mu _0\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)+\mathrm{sinh}(2r_0)\mathrm{cos}(2\phi _0)\mathrm{e}^{\mathrm{\Gamma }t}},$$
(32)
where $`\mu _0`$, $`r_0`$ and $`\phi _0`$ are, respectively, the initial purity and the initial squeezing parameters.

Let us first consider the case $`M=0`$, for which the initial state is damped toward a thermal state with mean photon number $`N`$ barnett ; marian . In this case, see Eq. (32), $`\phi `$ is constant in time and does not enter in the expression of $`\mu `$. The corresponding solutions for $`\mu (t)`$ and $`r(t)`$ read then as follows

$`\mu (t)`$ $`=`$ $`\mu _0\left[{\displaystyle \frac{\mu _0^2}{\mu _{\mathrm{}}^2}}\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)^2+2{\displaystyle \frac{\mu _0}{\mu _{\mathrm{}}}}\mathrm{e}^{\mathrm{\Gamma }t}\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)\mathrm{cosh}(2r_0)+\mathrm{e}^{2\mathrm{\Gamma }t}\right]^{1/2},`$ (33)
$`\mathrm{cosh}[2r(t)]`$ $`=`$ $`\mu (t)\left({\displaystyle \frac{1\mathrm{e}^{\mathrm{\Gamma }t}}{\mu _{\mathrm{}}}}+\mathrm{e}^{\mathrm{\Gamma }t}{\displaystyle \frac{\mathrm{cosh}(2r_0)}{\mu _0}}\right).`$ (34)
The quantities $`\mu (t)`$ and $`r(t)`$ in Eqns. (33) and (34) solve the following system of coupled equations

$`\dot{\mu }`$ $`=`$ $`\mathrm{\Gamma }\left(\mu {\displaystyle \frac{\mu ^2\mathrm{cosh}(2r)}{\mu _{\mathrm{}}}}\right),`$
$`\dot{r}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}{\displaystyle \frac{\mu }{\mu _{\mathrm{}}}}\mathrm{sinh}(2r),`$ (35)
which, in turn, can be directly found working out the basic evolution equation $`\dot{\mu }=2\mathrm{T}\mathrm{r}[\dot{\varrho }\varrho ]`$ as a phase–space integral and exploiting Eqns. (2022). It is easy to see that, as $`t\mathrm{}`$, $`\mu (t)\mu _{\mathrm{}}=(2N+1)^1`$ and $`r(t)0`$, as one expects, since the channel damps (pumps) the initial state to a thermal state with mean photon number $`N`$. Therefore, the only constant solution of Eq. (35) is $`\mu =\mu _{\mathrm{}}`$, $`r=0`$, i.e. only initial non–squeezed states are left unchanged by the evolution in the noisy channel.

In fact Eq. (33) shows that $`\mu (t)`$ is a decreasing function of $`r_0`$: in a non–squeezed channel ($`M=0`$), a squeezed state decoheres more rapidly than a non–squeezed one (see Figs. 5 and 6). Let us consider, for instance, an initially pure state in a channel with $`N=1`$ (so that $`\mu _{\mathrm{}}=\frac{1}{3}`$); after a time $`t=\mathrm{\Gamma }^1`$, the ratio of the purity of a state with $`r_0=1.5`$ to the purity of a state with $`r_0=0`$ is 53.7%. This dependence could therefore be relevant for practical purposes. The optimal evolution for the purity, obtained letting $`r=0`$ in Eq. (33), reads

$$\mu (t)=\frac{\mu _0\mu _{\mathrm{}}}{\mu _0+\mathrm{e}^{\mathrm{\Gamma }t}(\mu _{\mathrm{}}\mu _0)}.$$
(36)
Obviously, $`\mu (t)`$ is not necessarily a decreasing function of time: if $`\mu _0<\mu _{\mathrm{}}`$ then the initial state will undergo a certain amount of purification, asymptotically reaching the value $`\mu _{\mathrm{}}`$ which characterizes the channel, as shown in Fig. 6. In addition, $`\mu (t)`$ is not a monotonic function for any choice of the initial conditions. Letting $`\dot{\mu }=0`$ in Eq. (35), and exploiting Eqns. (33) and (34), one finds the following condition for the appearance of a zero of $`\dot{\mu }`$ at finite positive times: $`\mathrm{cosh}(2r_0)>\mathrm{max}(\frac{\mu _0}{\mu _{\mathrm{}}},\frac{\mu _{\mathrm{}}}{\mu _0})`$. If this condition is satisfied, then $`\mu (t)`$ shows a local extremum, in fact a minimum since, differentiating the first of Eqns. (35) and letting $`\dot{\mu }=0`$, one obtains $`\ddot{\mu }>0`$. This behavior is shown in Fig. 6.

Let us now treat the more general instance $`M0`$ of a squeezed thermal bath. Recalling the definition of $`\mu _{\mathrm{}}`$ and exploiting Eqns. (3032), one easily finds the asymptotic values of the physical parameters $`\mu `$, $`r`$ and $`\phi `$

$`\mu _{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{(2N+1)^24|M|^2}}},`$ (37)
$`\mathrm{cosh}(2r_{\mathrm{}})`$ $`=`$ $`\sqrt{1+4\mu _{\mathrm{}}^2|M|^2},`$ (38)
$`\mathrm{tan}(2\phi _{\mathrm{}})`$ $`=`$ $`{\displaystyle \frac{M_2}{M_1}}.`$ (39)
These values characterize the squeezed channel. Eq. (37) shows that, if $`M0`$, then $`N`$ is not simply the mean thermal photon number $`\overline{n}`$ of the asymptotic state. One has:

$$N=\frac{\sqrt{\left(2\overline{n}+1\right)^2+4|M|^2}1}{2}.$$

In order to understand the dynamics of purity when $`M0`$, it is convenient to write again the expression (30) for $`\mu (t)`$, using Eqns. (38) and (39) to switch from the complex parameter $`M=M_1+iM_2`$ to the asymptotic values of the squeezing parameters $`r_{\mathrm{}}`$ and $`\phi _{\mathrm{}}`$; one obtains

$`\mu (t)`$ $`=`$ $`\mu _0[{\displaystyle \frac{\mu _0^2}{\mu _{\mathrm{}}^2}}(1\mathrm{e}^{\mathrm{\Gamma }t})^2+\mathrm{e}^{2\mathrm{\Gamma }t}`$ (40)
$`+`$ $`2{\displaystyle \frac{\mu _0}{\mu _{\mathrm{}}}}(\mathrm{cosh}(2r_{\mathrm{}})\mathrm{cosh}(2r_0)+\mathrm{sinh}(2r_{\mathrm{}})\mathrm{sinh}(2r_0)\left(\mathrm{cos}(2\phi _{\mathrm{}}2\phi _0)\right))(1\mathrm{e}^{\mathrm{\Gamma }t})\mathrm{e}^{\mathrm{\Gamma }t}]^{1/2}.`$
We see from Eq. (40) that $`\mu (t)`$ is a monotonically decreasing function of the factor $`\mathrm{cos}(2\phi _{\mathrm{}}2\phi _0)`$, which gives the only dependence on the initial phase $`\phi _0`$ of the squeezing. Thus, for any given $`\phi _{\mathrm{}}`$ characterizing the squeezing of the bath, $`\phi _0=\phi _{\mathrm{}}+\frac{\pi }{2}`$ is the most favorable value of the initial angle of squeezing, i.e. the one which allows the maximum purity at a given time. For such a choice, $`\mu (t)`$ reduces to

$$\mu (t)=\mu _0\left[\frac{\mu _0^2}{\mu _{\mathrm{}}^2}\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)^2+\mathrm{e}^{2\mathrm{\Gamma }t}+2\frac{\mu _0}{\mu _{\mathrm{}}}\mathrm{cosh}(2r_{\mathrm{}}2r_0)\left(1\mathrm{e}^{\mathrm{\Gamma }t}\right)\mathrm{e}^{\mathrm{\Gamma }t}\right]^{1/2}.$$
(41)
This is a decreasing function of the factor $`\mathrm{cosh}(2r_{\mathrm{}}2r_0)`$, so that the maximum value of the purity at a given time is achieved for the choice $`r_0=r_{\mathrm{}}`$, and the evolution of the purity of a squeezed state in a squeezed channel is identical to the evolution of the purity of a non–squeezed state in a non–squeezed channel expressed by Eq. (36) and illustrated in Fig. 6.

In conclusion, for the most general instance of a channel characterized by arbitrary $`\mu _{\mathrm{}}`$, $`r_{\mathrm{}}`$, $`\phi _{\mathrm{}}`$ and $`\mathrm{\Gamma }`$, the initial Gaussian state for which purity is best preserved in time must have a squeezing parameter $`r_0=r_{\mathrm{}}`$ and a squeezing angle $`\phi _0=\phi _{\mathrm{}}+\frac{\pi }{2}`$, i.e. it must be antisqueezed (orthogonally squeezed) with respect to the bath. The net effect for the evolution of the purity is that the two orthogonal squeezings of the initial state and of the bath cancel each other exactly, thus reproducing the optimal purity evolution of an initial non–squeezed coherent state in a non–squeezed thermal bath.

## V Conclusions

We have shown that the purity of Gaussian states for continuous variable systems can be operationally determined by the joint measurement of two conjugate quadratures. In order to perform such a measurement, the minimal, necessary and sufficient requirement is that the measurement apparatus records data distributed according to the Husimi quasi–probability function. We have then verified by Monte Carlo simulated experiments the statistical reliability of the associated measurement schemes, thus proving the possibility of an experimentally realizable characterization of the purity of Gaussian states. We have compared as well the scheme based on the $`Q`$–function with the one based on single–quadrature detection, and showed that the former provides a more reliable statistics. Moreover, we have derived an evolution equation for the purity of Gaussian states in noisy channels, both in the case of a thermal and of a squeezed thermal bath. Our analysis shows that the purity is maximized at any given time for an initial coherent state evolving in a thermal bath, or for an initial squeezed state evolving in a squeezed thermal bath whose squeezing is orthogonal to that of the input state. We have focused our attention on the purity of single–mode Gaussian states. The time–evolution of the purity for specific initial non Gaussian states of great physical relevance can be studied, as well as the extension to Gaussian states of multimode systems, both pure and mixed. These topics are currently being explored and will be the subject of forthcoming work.

###### Acknowledgements.

The work of MGAP has been sponsored by INFM through the project PRA-2002-CLON. FI, AS and SDS thank INFM and INFN for financial support. SDS thanks the ESF COSLAB program. FI thanks the ESF BEC2000+ program.