Reproduce: Johansson et al., PRA 87 043804 (2013)

In [1]:
from sympy import *
init_printing()
In [2]:
# Note: This notebook requires and experimental branch of sympy available at: https://github.com/jrjohansson/sympy
from sympy.physics.quantum import *

Introduction

In this notebook I reproduce some results from Johansson et al., PRA 87, 043804 (2013), which is one of my own papers. Normally this would not qualify as a reproduced results, but in this case I will use a different method to obtain the same result. Here I will use SymPy and its quantum physics module to demonstrate how an analytical result can be obtained using symbolic operator manipulation, instead of calculations with pen and paper.

Here we are interested in evaluating the nonclassicality indicator [see, e.g., Miranowicz et al. Phys. Rev. A 82, 013824 (2010)]

$$ \langle: f^\dagger f :\rangle $$

which is positive for classical states, and where $f$ is a function of the field operator $b$. We will use the following choise of $f$, which is ideally suited for detecting the presence of nonclassial two-mode squeezing (nonclassical quadrature correlations between photons at difference frequencies):

$$ f = e^{i\theta} b_- + e^{-i\theta} b_-^\dagger + i(e^{i\theta} b_+ - e^{-i\theta} b_+^\dagger ) $$

where $b_\pm$ are the field operators at frequencies $\omega_\pm = \omega_d/2 \pm \delta\omega$. Note that the structure of $f$ is $f = x_- - p_+$ where $x_\pm$ and $p_\pm$ are the $x$ and $p$ quadrature operators of modes with frequency $\omega_\pm$, which suggest that this form is suitable for two-mode quadrature squeezing where exactly this operator combination is squeezed below its vacuum level.

Below we will evaluate $\langle: f^\dagger f :\rangle$ for the state of the field $b$ that is produced by the dynamical Casimir effect when the boundary condition of the field is varied in time. The relation between the output field $b$ and the input field $a$ (that impinges on the time-varying boundary condition) is, in the perturbative regime, given as

$$ b_{\pm} = -a_{\pm} - i\frac{\delta\!L_{\rm eff}}{v}\sqrt{\omega_-\omega_+} a_{\mp}^\dagger $$

Here $\frac{\delta\!L_{\rm eff}}{v}\sqrt{\omega_-\omega_+}$ is a small parameter in the perturbation calculation, which we for brievity will denote $x$ in the following (see PRA 87 043804 for details).

We start by defining symbols and operators:

In [3]:
theta = Symbol("theta", real=True)
x = Symbol("x", real=True)
In [4]:
dL, v, wp, wm = symbols("{\delta\!L}, v, omega_-, omega_+")
In [5]:
# input field operators
am = BosonOperator("a_-")
ap = BosonOperator("a_+")
In [6]:
# output field operators
bm = BosonOperator("b_-")
bp = BosonOperator("b_+")
In [7]:
# relation between input field operators and output field operators
bm_sub = -am - 1j * x * Dagger(ap)
bp_sub = -ap - 1j * x * Dagger(am)

bmd_sub = Dagger(bm_sub)
bpd_sub = Dagger(bp_sub)

$$ f = e^{i\theta} b_- + e^{-i\theta} b_-^\dagger + i(e^{i\theta} b_+ - e^{-i\theta} b_+^\dagger ) $$

In [8]:
f = exp(1j*theta) * bm + exp(-1j*theta) * Dagger(bm) + 1j * (exp(1j*theta) * bp - exp(-1j*theta) * Dagger(bp))

f
Out[8]:
$$e^{- 1.0 i \theta} {{b_-}^\dagger} + e^{1.0 i \theta} {b_-} + 1.0 i \left(- e^{- 1.0 i \theta} {{b_+}^\dagger} + e^{1.0 i \theta} {b_+}\right)$$

Now we construct $f^\dagger f$

In [9]:
fdf = Dagger(f) * f
In [10]:
fdf.expand()
Out[10]:
$$- 1.0 i e^{- 2.0 i \theta} {{b_+}^\dagger} {{b_-}^\dagger} - 1.0 e^{- 2.0 i \theta} \left({{b_+}^\dagger}\right)^{2} - 1.0 i e^{- 2.0 i \theta} {{b_-}^\dagger} {{b_+}^\dagger} + e^{- 2.0 i \theta} \left({{b_-}^\dagger}\right)^{2} + 1.0 i e^{2.0 i \theta} {b_+} {b_-} - 1.0 e^{2.0 i \theta} \left({b_+}\right)^{2} + 1.0 i e^{2.0 i \theta} {b_-} {b_+} + e^{2.0 i \theta} \left({b_-}\right)^{2} + 1.0 {{b_+}^\dagger} {b_+} - 1.0 i {{b_+}^\dagger} {b_-} + 1.0 {b_+} {{b_+}^\dagger} + 1.0 i {b_+} {{b_-}^\dagger} + 1.0 i {{b_-}^\dagger} {b_+} + {{b_-}^\dagger} {b_-} - 1.0 i {b_-} {{b_+}^\dagger} + {b_-} {{b_-}^\dagger}$$

Now we calculate the normal order expression $:f^\dagger f:$, with respect to the operators in the output field ($b_\pm$)

In [11]:
fdf_normal = normal_order(fdf.expand())
In [12]:
fdf_normal
Out[12]:
$$- 1.0 i e^{- 2.0 i \theta} {{b_+}^\dagger} {{b_-}^\dagger} - 1.0 e^{- 2.0 i \theta} \left({{b_+}^\dagger}\right)^{2} - 1.0 i e^{- 2.0 i \theta} {{b_-}^\dagger} {{b_+}^\dagger} + e^{- 2.0 i \theta} \left({{b_-}^\dagger}\right)^{2} + 1.0 i e^{2.0 i \theta} {b_+} {b_-} - 1.0 e^{2.0 i \theta} \left({b_+}\right)^{2} + 1.0 i e^{2.0 i \theta} {b_-} {b_+} + e^{2.0 i \theta} \left({b_-}\right)^{2} + 2.0 {{b_+}^\dagger} {b_+} - 2.0 i {{b_+}^\dagger} {b_-} + 2.0 i {{b_-}^\dagger} {b_+} + 2 {{b_-}^\dagger} {b_-}$$

Now substitute the output field operators with the input field operators:

In [13]:
fdf_1 = fdf_normal.expand().subs({bm: bm_sub, Dagger(bm): bmd_sub, bp: bp_sub, Dagger(bp): bpd_sub}).expand()

fdf_1
Out[13]:
$$1.0 i x^{2} e^{- 2.0 i \theta} {a_+} {a_-} - 1.0 x^{2} e^{- 2.0 i \theta} \left({a_+}\right)^{2} + 1.0 i x^{2} e^{- 2.0 i \theta} {a_-} {a_+} + 1.0 x^{2} e^{- 2.0 i \theta} \left({a_-}\right)^{2} - 1.0 i x^{2} e^{2.0 i \theta} {{a_+}^\dagger} {{a_-}^\dagger} - 1.0 x^{2} e^{2.0 i \theta} \left({{a_+}^\dagger}\right)^{2} - 1.0 i x^{2} e^{2.0 i \theta} {{a_-}^\dagger} {{a_+}^\dagger} + 1.0 x^{2} e^{2.0 i \theta} \left({{a_-}^\dagger}\right)^{2} + 2.0 x^{2} {a_+} {{a_+}^\dagger} + 2.0 i x^{2} {a_+} {{a_-}^\dagger} - 2.0 i x^{2} {a_-} {{a_+}^\dagger} + 2.0 x^{2} {a_-} {{a_-}^\dagger} - 1.0 x e^{- 2.0 i \theta} {{a_+}^\dagger} {a_+} + 1.0 i x e^{- 2.0 i \theta} {{a_+}^\dagger} {a_-} - 1.0 x e^{- 2.0 i \theta} {a_+} {{a_+}^\dagger} - 1.0 i x e^{- 2.0 i \theta} {a_+} {{a_-}^\dagger} - 1.0 i x e^{- 2.0 i \theta} {{a_-}^\dagger} {a_+} - 1.0 x e^{- 2.0 i \theta} {{a_-}^\dagger} {a_-} + 1.0 i x e^{- 2.0 i \theta} {a_-} {{a_+}^\dagger} - 1.0 x e^{- 2.0 i \theta} {a_-} {{a_-}^\dagger} - 1.0 x e^{2.0 i \theta} {{a_+}^\dagger} {a_+} + 1.0 i x e^{2.0 i \theta} {{a_+}^\dagger} {a_-} - 1.0 x e^{2.0 i \theta} {a_+} {{a_+}^\dagger} - 1.0 i x e^{2.0 i \theta} {a_+} {{a_-}^\dagger} - 1.0 i x e^{2.0 i \theta} {{a_-}^\dagger} {a_+} - 1.0 x e^{2.0 i \theta} {{a_-}^\dagger} {a_-} + 1.0 i x e^{2.0 i \theta} {a_-} {{a_+}^\dagger} - 1.0 x e^{2.0 i \theta} {a_-} {{a_-}^\dagger} + 2.0 i x {{a_+}^\dagger} {{a_-}^\dagger} + 2.0 x \left({{a_+}^\dagger}\right)^{2} - 2.0 i x {a_+} {a_-} + 2.0 x \left({a_+}\right)^{2} + 2.0 i x {{a_-}^\dagger} {{a_+}^\dagger} - 2.0 x \left({{a_-}^\dagger}\right)^{2} - 2.0 i x {a_-} {a_+} - 2.0 x \left({a_-}\right)^{2} - 1.0 i e^{- 2.0 i \theta} {{a_+}^\dagger} {{a_-}^\dagger} - 1.0 e^{- 2.0 i \theta} \left({{a_+}^\dagger}\right)^{2} - 1.0 i e^{- 2.0 i \theta} {{a_-}^\dagger} {{a_+}^\dagger} + e^{- 2.0 i \theta} \left({{a_-}^\dagger}\right)^{2} + 1.0 i e^{2.0 i \theta} {a_+} {a_-} - 1.0 e^{2.0 i \theta} \left({a_+}\right)^{2} + 1.0 i e^{2.0 i \theta} {a_-} {a_+} + e^{2.0 i \theta} \left({a_-}\right)^{2} + 2.0 {{a_+}^\dagger} {a_+} - 2.0 i {{a_+}^\dagger} {a_-} + 2.0 i {{a_-}^\dagger} {a_+} + 2 {{a_-}^\dagger} {a_-}$$

Write all input operators in normal ordered form:

In [14]:
fdf_2 = normal_ordered_form(fdf_1.expand(), independent=True)

fdf_2
Out[14]:
$$1.0 i x^{2} e^{- 2.0 i \theta} {a_+} {a_-} - 1.0 x^{2} e^{- 2.0 i \theta} \left({a_+}\right)^{2} + 1.0 i x^{2} e^{- 2.0 i \theta} {a_-} {a_+} + 1.0 x^{2} e^{- 2.0 i \theta} \left({a_-}\right)^{2} - 1.0 i x^{2} e^{2.0 i \theta} {{a_+}^\dagger} {{a_-}^\dagger} - 1.0 x^{2} e^{2.0 i \theta} \left({{a_+}^\dagger}\right)^{2} - 1.0 i x^{2} e^{2.0 i \theta} {{a_-}^\dagger} {{a_+}^\dagger} + 1.0 x^{2} e^{2.0 i \theta} \left({{a_-}^\dagger}\right)^{2} + 4.0 x^{2} + 2.0 x^{2} {{a_+}^\dagger} {a_+} - 2.0 i x^{2} {{a_+}^\dagger} {a_-} + 2.0 i x^{2} {{a_-}^\dagger} {a_+} + 2.0 x^{2} {{a_-}^\dagger} {a_-} - 2.0 x e^{- 2.0 i \theta} - 2.0 x e^{- 2.0 i \theta} {{a_+}^\dagger} {a_+} + 2.0 i x e^{- 2.0 i \theta} {{a_+}^\dagger} {a_-} - 2.0 i x e^{- 2.0 i \theta} {{a_-}^\dagger} {a_+} - 2.0 x e^{- 2.0 i \theta} {{a_-}^\dagger} {a_-} - 2.0 x e^{2.0 i \theta} - 2.0 x e^{2.0 i \theta} {{a_+}^\dagger} {a_+} + 2.0 i x e^{2.0 i \theta} {{a_+}^\dagger} {a_-} - 2.0 i x e^{2.0 i \theta} {{a_-}^\dagger} {a_+} - 2.0 x e^{2.0 i \theta} {{a_-}^\dagger} {a_-} + 2.0 i x {{a_+}^\dagger} {{a_-}^\dagger} + 2.0 x \left({{a_+}^\dagger}\right)^{2} - 2.0 i x {a_+} {a_-} + 2.0 x \left({a_+}\right)^{2} + 2.0 i x {{a_-}^\dagger} {{a_+}^\dagger} - 2.0 x \left({{a_-}^\dagger}\right)^{2} - 2.0 i x {a_-} {a_+} - 2.0 x \left({a_-}\right)^{2} - 1.0 i e^{- 2.0 i \theta} {{a_+}^\dagger} {{a_-}^\dagger} - 1.0 e^{- 2.0 i \theta} \left({{a_+}^\dagger}\right)^{2} - 1.0 i e^{- 2.0 i \theta} {{a_-}^\dagger} {{a_+}^\dagger} + e^{- 2.0 i \theta} \left({{a_-}^\dagger}\right)^{2} + 1.0 i e^{2.0 i \theta} {a_+} {a_-} - 1.0 e^{2.0 i \theta} \left({a_+}\right)^{2} + 1.0 i e^{2.0 i \theta} {a_-} {a_+} + e^{2.0 i \theta} \left({a_-}\right)^{2} + 2.0 {{a_+}^\dagger} {a_+} - 2.0 i {{a_+}^\dagger} {a_-} + 2.0 i {{a_-}^\dagger} {a_+} + 2 {{a_-}^\dagger} {a_-}$$

Assume thermal input field

Elimitate operators that do not contribute to expectation values assuming thermal input field. The only operators that contribute are $a_-^\dagger a_-$ and $a_+^\dagger a_+$, and expressions such as for example $\langle a_-a_-\rangle = 0$.

In [15]:
# assume thermal input: first remove nonthermal operator combinations
fdf_3 = fdf_2.subs({am**2: 0, ap**2: 0,
                    ap*am: 0, am*ap: 0, 
                    Dagger(ap) * am: 0, Dagger(am) * ap: 0, 
                    Dagger(ap) * Dagger(am): 0, Dagger(am) * Dagger(ap): 0, 
                    Dagger(ap)**2: 0, 
                    Dagger(am)**2: 0})

fdf_3
Out[15]:
$$4.0 x^{2} + 2.0 x^{2} {{a_+}^\dagger} {a_+} + 2.0 x^{2} {{a_-}^\dagger} {a_-} - 2.0 x e^{- 2.0 i \theta} - 2.0 x e^{- 2.0 i \theta} {{a_+}^\dagger} {a_+} - 2.0 x e^{- 2.0 i \theta} {{a_-}^\dagger} {a_-} - 2.0 x e^{2.0 i \theta} - 2.0 x e^{2.0 i \theta} {{a_+}^\dagger} {a_+} - 2.0 x e^{2.0 i \theta} {{a_-}^\dagger} {a_-} + 2.0 {{a_+}^\dagger} {a_+} + 2 {{a_-}^\dagger} {a_-}$$

Now replace the operators that do contibute with their expectation values:

$$ \langle a_-^\dagger a_-\rangle = n_{\rm th}^- $$

$$ \langle a_+^\dagger a_+\rangle = n_{\rm th}^+ $$

In [16]:
n_th_m, n_th_p = symbols("{n^-_{th}}, {n^+_{th}}", positive=True)
In [17]:
# thermal operator combinations 
fdf_4 = fdf_3.subs({Dagger(am) * am: n_th_m, Dagger(ap) * ap: n_th_p})

fdf_4
Out[17]:
$$2.0 x^{2} {n^{+}_{{th}}} + 2.0 x^{2} {n^{-}_{{th}}} + 4.0 x^{2} - 2.0 x {n^{+}_{{th}}} e^{- 2.0 i \theta} - 2.0 x {n^{+}_{{th}}} e^{2.0 i \theta} - 2.0 x {n^{-}_{{th}}} e^{- 2.0 i \theta} - 2.0 x {n^{-}_{{th}}} e^{2.0 i \theta} - 2.0 x e^{- 2.0 i \theta} - 2.0 x e^{2.0 i \theta} + 2.0 {n^{+}_{{th}}} + 2 {n^{-}_{{th}}}$$

Replace $e^{\pm i\theta}$ with trigonometric expressions, and simplify the expression:

In [18]:
fdf_5 = fdf_4.collect(x).replace(exp(2j*theta), cos(2*theta) + 1j * sin(2*theta)).\
                         replace(exp(-2j*theta), cos(2*theta) - 1j * sin(2*theta))
    
fdf_6 = fdf_5.simplify().collect(cos(2*theta)).collect(x).simplify()

fdf_6
Out[18]:
$$x^{2} \left(2.0 {n^{+}_{{th}}} + 2.0 {n^{-}_{{th}}} + 4.0\right) - 4.0 x \left({n^{+}_{{th}}} + {n^{-}_{{th}}} + 1\right) \cos{\left (2 \theta \right )} + 2.0 {n^{+}_{{th}}} + 2.0 {n^{-}_{{th}}}$$

Drop the $x^2$ term because our relation between the output field and input field is only valid to first order in $x$:

In [19]:
fdf_7 = fdf_6.subs(x**2, 0)

fdf_7
Out[19]:
$$- 4.0 x \left({n^{+}_{{th}}} + {n^{-}_{{th}}} + 1\right) \cos{\left (2 \theta \right )} + 2.0 {n^{+}_{{th}}} + 2.0 {n^{-}_{{th}}}$$

Substitute the small parameter $x$ with the physical parameters it represents:

In [20]:
fdf_8 = fdf_7.subs(x**2, 0).subs(x, dL/v * sqrt(wp * wm))

fdf_8
Out[20]:
$$2.0 {n^{+}_{{th}}} + 2.0 {n^{-}_{{th}}} - \frac{4.0 {\delta\!L}}{v} \sqrt{\omega_{+} \omega_{-}} \left({n^{+}_{{th}}} + {n^{-}_{{th}}} + 1\right) \cos{\left (2 \theta \right )}$$

This is equation 14 in Phys. Rev. A 87, 043804 (2013).

The interesting regime is when this expression is negative, because it indicates that there are nonclassical correlations between the quadrature at the different frequencies $\omega_\pm$. It is clear that $\theta = 0$ minimizes the expression, and that for finite temperature the expression is initially positive (the prefactor to the other term is small). We can find the transition point from positive to negative by solving $\langle:f^\dagger f:\rangle = 0$ for $x = \frac{\delta\!L_{\rm eff}}{v}\sqrt{\omega_-\omega_+}$:

In [21]:
# switch back to using x as small paramter, and set theta = 0
fdf_9 = fdf_8.subs(dL/v * sqrt(wp*wm), x).subs(theta, 0)

fdf_9
Out[21]:
$$- 4.0 x \left({n^{+}_{{th}}} + {n^{-}_{{th}}} + 1\right) + 2.0 {n^{+}_{{th}}} + 2.0 {n^{-}_{{th}}}$$
In [22]:
fdf_neg = solve(fdf_9, x)[0]

fdf_neg
Out[22]:
$$\frac{0.5 {n^{+}_{{th}}} + 0.5 {n^{-}_{{th}}}}{{n^{+}_{{th}}} + {n^{-}_{{th}}} + 1.0}$$

Which in the low tempeature case ($n_{\rm th}^+ + n_{\rm th}^- \ll 1$) simply gives

$$ \frac{\delta\!L_{\rm eff}}{v}\sqrt{\omega_-\omega_+} > \sim (n_{\rm th}^+ + n_{\rm th}^-)/2 $$

as a condition for when the output field becomes nonclassical.

Versions

In [23]:
#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%reload_ext version_information

%version_information sympy
Out[23]:
SoftwareVersion
Python3.3.2+ (default, Oct 9 2013, 14:50:09) [GCC 4.8.1]
IPython1.1.0
OSposix [linux]
sympy0.7.4.1-git
Wed Dec 18 14:05:32 2013 KST