arXiv:1612.05144v1 [quant-ph] 15 Dec 2016

Maximizing optomechanical entanglement with
optimal control
Dionisis Stefanatos
Hellenic Army Academy, Vari, Athens 16673, Greece
E-mail: dionisis@post.harvard.edu
October 2016
Abstract.
In this article, we formulate the generation of optomechanical entanglement
between the linearly coupled cavity field and the mechanical resonator as an optimal
control problem in hyperbolic space H 3 , with control input the coupling rate of the
two oscillators. Next, we use optimal control theory to find the allowed optimal values
of the coupling which maximize the amount of generated entanglement for a fixed
duration of the interaction. Finally, we employ a numerical optimization method to
obtain the exact optimal pulse sequences for several illustrative examples. In the strong
coupling regime, where the coupling rate is comparable or larger than the frequency
of the mechanical resonator, a substantial amount of entanglement can be generated
within a fraction of a single oscillator period.

PACS numbers: 42.50.Wk, 03.67.Bg, 02.30.Yy, 02.60.Pn

Keywords: optomechanical entanglement, quantum control, optimal control
1. Introduction
Cavity optomechanics, the field of physics studying the interactions of light confined
in an optical cavity with the mechanical resonator forming the movable end of the
cavity, is a very active area of research with many important applications [1]. These
include quantum information processing, for the interconversion between storage and
communication qubits, high-precision measurements of tiny masses, displacements
and forces, and fundamental tests of quantum physics at a macroscopic scale. One
fascinating application is the creation of entanglement between the electromagnetic field
trapped in the cavity and the motion of the macroscopic mechanical resonator. Several
theoretical works have suggested the creation of this optomechanical entanglement in
the so-called steady state regime, see for example [2, 3]. In these and all the related
studies, a continuous-wave light field is applied to the optomechanical system and drives

Maximizing optomechanical entanglement with optimal control

2

it to a steady state, where the cavity filed and the mechanical motion of the oscillator
are entangled.
The major drawback of the steady state schemes is that the amount of generated
entanglement is limited by the requirement of working in a stable stationary state.
In order to overcome this limitation, an alternative pulsed approach was suggested
[4]. According to this protocol, the entanglement is not created with the application
of a continuous-wave field but with a light pulse, avoiding thus the limits imposed
by stability requirements. This pulsed scheme was implemented experimentally for
generating entanglement between mechanical motion and the microwave radiation field
[5].
In the current paper we also consider a pulsed scheme, where the coupling strength
between the cavity filed and the mechanical resonator is modulated in order to create
the optomechanical entanglement. This is obviously different than the pulsed protocol
discussed above. Another important difference is that here we concentrate in the
strong coupling regime, where the coupling strength is comparable or larger than the
frequency of the mechanical oscillator, while the pulsed protocol [4] focuses in the
week coupling regime (coupling much lower than the mechanical frequency), where the
rotating wave approximation can be employed. An important consequence of working
with large coupling is that a substantial amount of entanglement can be created within
a fraction of the period of the mechanical motion. Note that the pulsed modulation of
the coupling has been mostly used for the numerical optimization of optomechanical
cooling [6, 7, 8, 9]. Recently, it was also employed to numerically investigate the
fast entanglement creation between two oscillators (for example two mechanical modes)
coupled via their interactions with a third oscillator (optical mode) [9]. Here we also
formulate entanglement generation as an optimal control problem, but before recoursing
to numerical optimization we use optimal control theory to draw as much information
as possible for the optimal solution that maximizes entanglement.
The present article is formulated as follows. In the next section we present some
basic facts about optomechanical entanglement and identify the target state, while
in section 3 we formulate entanglement generation as an optimal control problem in
hyperbolic space H 3 . In section 4 we use optimal control theory to study the resultant
problem and obtain the allowed optimal values of the coupling strength. In section 5
we use a numerical optimization method to find the exact optimal pulse sequences for
some illustrative examples, while section 6 concludes this work.
2. Optomechanical entanglement
In a laser driven optical cavity, the light from the laser entering the cavity is reflected
between the mirrors and the enhanced electromagnetic field built inside the cavity exerts
a force to the movable mirror due to radiation pressure. The interaction between the
cavity (photon) field and the mechanical vibration of the mirror (phonon field) can be
described by the following Hamiltonian [1], sufficient for the study of many interesting

Maximizing optomechanical entanglement with optimal control

3

phenomena in cavity optomechanics,
H = −~∆â+ â + ~ωm b̂+ b̂ + ~g̃(t)(â+ + â)(b̂+ + b̂),

(1)

where ~ is Planck’s constant, ∆ is roughly the laser detuning from the cavity frequency
and ωm is the frequency of the mechanical resonator. The operators â, b̂ are the
annihilation operators for the photon and the phonon fields, respectively, while â+ , b̂+
are the corresponding creation operators. The first two terms in (1) express the energy
of each individual oscillator (optical and mechanical), while the last term the interaction
between them. The coupling rate g̃(t) has dimensions of frequency and is the available
control which can be altered with time, with the aim to entangle the mechanical motion
of the mirror with the cavity field (optomechanical entanglement).
The state of the system can be described by the density matrix ρ(t), which satisfies
the Liouville-von Neumann equation [10]
ρ̇ = −i[H/~, ρ].

(2)

Note that here we consider only the coherent evolution and ignore relaxation, which
describes the undesirable interaction of the system with its environment, since we are
interested in the fast pulsed regime, where the evolution takes place within a fraction
of a single oscillation period, and not the steady state regime. This is a legitimate
practice in minimum-time quantum control problems [8, 9], the alternative being to
include relaxation and maximize the fidelity of the final point [6]. We will concentrate
on the so-called “blue-detuned regime” ∆ = ωm , the case which is more relevant for
entanglement generation as we immediately explain. Observe that the interaction term
in (1) can be decomposed into two terms as follows
(â+ + â)(b̂+ + b̂) = (â+ b̂ + âb̂+ ) + (â+ b̂+ + âb̂).
The first term is suitable for optomechanical cooling and state swapping, while the
second one for entanglement generation [4]. For a blue-detuned laser, the second term
is resonant with the cavity, and this is why we concentrate on this case for studying the
generation of optomechanical entanglement. If we normalize time using the frequency
ωm , i.e. set tnew = ωm told , then the Liouville-von Neumann equation takes the following
form, with g(t) = g̃(t)/ωm
ρ̇ = −i[H, ρ], H = −â+ â + b̂+ b̂ + g(t)(â+ + â)(b̂+ + b̂).

(3)

We consider that at t = 0 the quantum system starts from the vacuum state
|ψ(0)i = |0ia |0ib.

(4)

The entangled state which we would like to generate at the final time t = T is closely
related to the well known two-mode squeezed vacuum state [11]
∞
X
1
+ +
tanhn r|nia |nib ,
(5)
|φi = er(â b̂ −âb̂) |0ia |0ib =
cosh
r
n=0

where r ≥ 0 is the squeezing parameter. Note that the above expansion in terms of the
number states can be obtained following the procedure described in [11]. Before turning

Maximizing optomechanical entanglement with optimal control

4

our attention to the target state that we actually consider in this article, we present
some characteristics of state (5) which we will need in the analysis that follows.
First of all note that, since (5) is a pure state, its entanglement is quantified by the
von Neumann entropy of the reduced state [12]
∞
X
S=−
pn ln pn = cosh2 r ln(cosh2 r) − sinh2 r ln(sinh2 r),
(6)
n=0

where the probabilities
2

tanhn r
pn =
cosh r

are taken from the expansion (5) and the final result is obtained by summing the series
appropriately. If we set x = sinh2 r, then
S(x) = (1 + x) ln (1 + x) − x ln x,
which is an increasing function of x > 0. But x is also an increasing function for r ≥ 0,
thus the entanglement is an increasing function of the squeezing parameter r.
Next, we specifically identify the quadratures which are squeezed. The annihilation
operators âφ , b̂φ of the two modes in state |φi, which is related to the vacuum state by
the unitary transformation given in (5), are
+ b̂+ −âb̂)

âφ = e−r(â

−r(â+ b̂+ −âb̂)

b̂φ = e

+ b̂+ −âb̂)

= â cosh r + b̂+ sinh r,

(7a)

r(â+ b̂+ −âb̂)

= â+ sinh r + b̂ cosh r,

(7b)

âer(â
b̂e

where the last terms in the above equations are obtained from the Taylor expansion of
the middle terms. Now, if we denote the quadrature observables of a mode â with
X̂a =

â + â+
√ ,
2

X̂ φ =

X̂aφ − X̂bφ
√
,
2

P̂a =

â − â+
√
i 2

(8)

and define
P̂ φ =

P̂aφ + P̂bφ
√
,
2

(9)

where X̂aφ , P̂aφ and X̂bφ , P̂bφ are the quadrature observables corresponding to âφ , b̂φ ,
respectively, then it is not hard to verify that
e−r
e−r
X̂ φ = √ (X̂a − X̂b ), P̂ φ = √ (P̂a + P̂b )
2
2
From the last relation we find
e−2r
e−2r
h∆(X̂ φ )2 i =
[∆X̂a2 + ∆X̂b2 ] =
,
2
2
e−2r
e−2r
[∆P̂a2 + ∆P̂b2 ] =
,
h∆(P̂ φ )2 i =
2
2

(10)

(11a)
(11b)

since ∆X̂a2 = ∆X̂b2 = ∆P̂a2 = ∆P̂b2 = 1/2 in the vacuum state. Thus, X̂ φ and P̂ φ are
squeezed, while note that
[X̂ φ , P̂ φ ] = 0,

(12)

Maximizing optomechanical entanglement with optimal control

5

and this is why we used the calligraphic symbol in P̂ φ .
As we explain in the next section, state (5) cannot be generated from the initial
vacuum state under evolution (3), in the blue-detuned regime. For this reason, we
consider the following modified target state at the final time t = T
∞
X
π
π
1
+
+
tanhn re−in 2 |nia |nib .
(13)
|ψ(T )i = e−i 4 (â â+b̂ b̂) |φi =
cosh
r
n=0
The annihilation operators âψ , b̂ψ of the two modes in state |ψ(T )i are
π

+ â+b̂+ b̂)

âψ = ei 4 (â

i π4 (â+ â+b̂+ b̂)

b̂ψ = e

π

+ â+b̂+ b̂)

âφ e−i 4 (â

−i π4 (â+ â+b̂+ b̂)

b̂φ e

π

= e−i 4 âφ ,

(14a)

−i π4

(14b)

=e

b̂φ ,

where the last terms in the above equations are obtained from the Taylor expansion of
the middle terms. If we define
X̂aψ + P̂bψ
X̂b + P̂a
√
X̂ ψ =
, P̂ ψ = ψ√ ψ
(15)
2
2
where X̂aψ , X̂bψ , P̂aψ , P̂bψ are given by (8) for âψ , b̂ψ , then it is not hard to verify, using
(14a), (14b) and (9), that
X̂ ψ =

X̂ φ + P̂ φ
√
,
2

P̂ ψ =

−X̂ φ + P̂ φ
√
.
2

(16)

Thus
1
e−2r
h∆(X̂ ψ )2 i = h∆(P̂ ψ )2 i = [∆(X̂ φ )2 + ∆(P̂ φ )2 ] =
2
2

(17)

[X̂ ψ , P̂ ψ ] = 0,

(18)

and

where we have used (11a), (11b) and (12). In the target state (13) the squeezed
quadratures are X̂ ψ and P̂ ψ , defined in (15).
In order to motivate why we would like to generate such a state, we close this
section by briefly reminding its role in the optomechanical teleportation protocol [4].
According to this protocol, a first light pulse is entangled with the mechanical motion of
the mirror forming state |ψ(T )i given in (13), while a second light pulse, prepared in the
state |ψ ′ i which will be teleported, interacts with the first pulse in a beam splitter. Two
homodyne detectors at the output ports of the beam splitter measure the quadratures
q̂X = P̂aψ + X̂aψ′ , q̂P = X̂aψ + P̂aψ′ . Subsequently, the mirror is displaced in position and
momentum by the outcome of these measurements. This feedback brings the mirror to
a state |f i, with position and momentum given by
√
X̂bf = X̂bψ + q̂X = X̂aψ′ + (X̂bψ + P̂aψ ) = X̂aψ′ + 2P̂ ψ ,
√
P̂bf = P̂bψ + q̂P = P̂aψ′ + (X̂aψ + P̂bψ ) = P̂aψ′ + 2X̂ ψ .
For large values of the squuzing parameter r → ∞ the terms proportional to X̂ ψ , P̂ ψ
are suppressed and X̂bf → X̂aψ′ , P̂bf → P̂aψ′ , in other words the state |ψ ′ i is teleported
to the mirror.

Maximizing optomechanical entanglement with optimal control

6

3. Entanglement creation as an optimal control problem in hyperbolic
space H 3
For the evolution described by (3), a closed set of equations can be obtained for the
second moments of the creation and annihilation operators of the two resonators, for
example â+ â, b̂+ b̂, â+ b̂, â2 etc., see [6]. In order to formulate entanglement generation as
a control problem, we do not follow the standard approach to use directly the moment
operators [6], but rather use specific linear combinations of them, a set of ten operators
which are the generators of the symplectic group Sp(4) [13], introduced by Dirac [14]
1
1
Jˆ0 = (â+ â + b̂b̂+ ), Jˆ1 = (â+ â − b̂+ b̂)
2
2
1
1
Jˆ2 = (â+ b̂ + âb̂+ ), Jˆ3 = (â+ b̂ − âb̂+ )
2
2i
i
1 + +
K̂1 = (â b̂ + âb̂), Q̂1 = (â+ b̂+ − âb̂)
2
2
i
1 + 2
K̂2 = − [(â ) + â2 − (b̂+ )2 − b̂2 ], K̂3 = [(â+ )2 − â2 + (b̂+ )2 − b̂2 ],
4
4
1
i + 2
Q̂2 = − [(â ) − â2 − (b̂+ )2 + b̂2 ], Q̂3 = − [(â+ )2 + â2 + (b̂+ )2 + b̂2 ].(19)
4
4
The matrix representation of these operators in terms of Pauli matrices can be found in
[15]. They satisfy the following commutation relations
[Jˆi , Jˆj ] = iǫijk Jˆk , [Jˆi , K̂j ] = iǫijk K̂k , [Jˆi , Q̂j ] = iǫijk Q̂k
[K̂i , Q̂j ] = iδij Jˆ0 , [K̂i , K̂j ] = [Q̂i , Q̂j ] = −iǫijk Jˆk ,
[Jˆi , Jˆ0 ] = 0, [K̂i , Jˆ0 ] = iQ̂i , [Q̂i , Jˆ0 ] = −iK̂i .

(20)

where ǫijk is the Levi-Civita symbol, which is 1 if (i, j, k) is an even permutation of
(1, 2, 3), −1 if it is an odd permutation, and 0 if any index is repeated, while δij is
Kronecker’s delta.
In terms of operators (19), the optomechanical Hamiltonian H (3) can be written
as
H = 2[−Jˆ1 + g(t)(K̂1 + Jˆ2 )].

(21)

The expectation value of an operator Ô is O = hÔi = Tr(ρÔ) and, using expression (21)
for the Hamiltonian in the Liouville-von Neumann equation (3) and the commutation
relations (20), we obtain the following systems for the expectation values of the above
operators
 



0
0 2g(t) −2g(t)
Q1
Q̇1
 



0
0
2
0
 Q̇2  
  Q2 
(22)
=


,
0
0
 Q̇3   −2g(t) −2
  Q3 
J˙0
−2g(t) 0
0
0
J0

Maximizing optomechanical entanglement with optimal control
 

0
0
2g(t)
0
0
0
K̇1
 K̇  
0
0
2
0
0
2g(t)
 2  
 K̇   −2g(t) −2
0
0
−2g(t)
0
 3  
 ˙ =
 J1  
0
0
0
0
0
2g(t)
 

˙
 J2  
0
0
−2g(t)
0
0
2
˙
J3
0
2g(t)
0
−2g(t)
−2
0

7










K1
K2
K3
J1
J2
J3







.




(23)

We next move to express the initial and target states of the above systems. Recall
that we start from the vacuum state (4), thus the corresponding initial conditions for
the expectation values of the operators defined in (19) are


 
1
(24a)
Q1 Q2 Q3 J0 = 0 0 0 2

 

(24b)
K1 K2 K 3 J 1 J 2 J 3 = 0 0 0 0 0 0 .

The target entangled state is (13) and using (14a), (14b) and (7a), (7b) we find the
following expectation values of operators (19) in this state


 
cosh 2r
sinh 2r
0 0
Q1 Q2 Q3 J0 = − 2
2

 

K1 K2 K 3 J 1 J 2 J 3 = 0 0 0 0 0 0 .

Now observe from (24b) that the initial values of K1 , K2 , K3 , J1 , J2 , J3 are zero.
System (23) is linear in these variables and homogeneous, thus these expectation values
remain zero during the whole time evolution. Consequently, when starting from the
vacuum state, the generation of the target entangled state is solely described by system
(22). If we normalize the expectation values Q1 , Q2 , Q3 , J0 with 12 , but keep the same
notation for these variables, then system (22) remains unchanged, while the starting
and target points become

 

(25a)
Q1 Q2 Q3 J0 = 0 0 0 1

 

(25b)
Q1 Q2 Q3 J0 = − sinh 2r 0 0 cosh 2r .

Using the initial condition (25a) and Eq. (22) we find that the system variables satisfy
Q21 + Q22 + Q23 − J02 = −1,

(26)

which is the equation of a hyperboloid of two sheets. But from the definition of Jˆ0 in
(19) it is obviously J0 > 0, thus system (22) actually evolves in the upper sheet, the
hyperbolic space H 3 [16].
This space is closely related to the Minkowski space. We use the usual notation for
contravariant and covariant four-vectors


xµ = Q1 Q2 Q3 J0
(27a)


xµ = Q1 Q2 Q3 −J0
(27b)

and for the scalar product between two four-vectors a, b
a · b = aµ bµ ,

(28)

Maximizing optomechanical entanglement with optimal control

8

where the Einstein summation convention over repeated indices is assumed. The
geodesic distance d between two points in the space, represented by a and b, is given by
d = cosh−1 (−a · b),

(29)

where note that for points in H 3 it can be shown that it is always
a · b ≤ −1,

(30)

thus the geodesic distance is well defined. For the initial vacuum state (25a) and the
final squeezed state (25b) we have
x(0) · x(T ) = − cosh 2r,

(31)

and from (29) we find the geodesic distance to be
d = 2r,

(32)

twice the squeezing parameter, an aesthetically pleasant expression.
Up to now we have implicitly assumed that the value of the squeezing parameter r
at the final state is given. In such a case, a meaningful question to ask is how should we
pick the time-varying coupling g(t), which is bounded due to experimental limitations
as below
− G ≤ g(t) ≤ G,

G > 0,

(33)

in order to reach the target point in minimum time T . Although in the subsequent
sections we do not actually study this minimum-time problem but a closely related
one, here we make some simple observations about it. First of all note that, since the
dynamics are governed by system (22), which imposes certain restrictions on the allowed
paths and speed, the path corresponding to the shortest geodesic distance may not be
permitted or may not be the fastest. Next, consider the quantity
ẋ2 = ẋµ ẋµ = Q̇21 + Q̇22 + Q̇23 − J˙02
n
o
= 4 g 2 (t)(Q3 − J0 )2 + 2g(t)Q1 Q2 + Q22 + Q23 .

(34)

With some algebra we can show that ẋ2 ≥ 0 for points in H 3 , and note that there is no
contradiction with (30), since ẋ belongs to the tangent space of H 3 . Since the quantity
√
ẋ2 actually corresponds to the speed in the flat (3 + 1) Minkowski space with metric
tensor diag(1, 1, 1, −1), where H 3 is embedded, it is in principle desirable to move along
paths where this quantity is maximized. But from (34) observe that ẋ2 is actually a
quadratic function of the control g(t), with the coefficient of g 2 (t) being positive
(Q3 − J0 )2 > 0,
since Q3 = J0 would imply in (26) Q21 + Q22 = −1, which is obviously not true. Thus
the instantaneous speed ẋ2 is a convex function of the coupling, which is restricted as
in (33), and is maximized when g(t) takes values at the boundaries ±G. Note of course
that, since the choice of g(t) also affects the path travelled, there might be cases where
the minimum-time path is not the one along which the speed is maximized.

Maximizing optomechanical entanglement with optimal control

9

Inspired from the above analysis and the well known fact that optomechanical state
swap can be achieved for specific constant values of the coupling rate and the duration
T [8], it is tempting to examine whether a constant coupling g(t) = g can transfer the
system to the desired target state. We show that, unfortunately, this is not the case
for squeezing. It can be verified using the system equation (22) that for constant g(t)
the quantity ẋ2 given in (34) is also constant. But using the initial and final conditions
(25a) and (25b) we obtain
ẋ2 (T ) = 4g 2 cosh2 2r ≥ 4g 2 = ẋ2 (0),
where the equality holds only for r = 0. The conclusion is that in order to reach the
target state the coupling has to change in time.
The problem that we actually study in the rest of the paper is closely related to
the minimum-time problem stated above. Specifically, we fix the duration T and ask
what is the maximum value of the squeezing parameter, corresponding to the maximum
entanglement, which can be obtained at the final time t = T . In order to precisely
formulate this problem and also reduce the dimension of the system from four to three,
we use the constant of the motion (26). Recall that J0 > 0, thus
q
(35)
J0 = 1 + Q21 + Q22 + Q23 .

Using this relation we end up with the three-dimensional system

q
2
2
2
1 + Q1 + Q2 + Q3 − Q3
Q̇1 = −2g

(36a)

Q̇2 = 2Q3

(36b)

Q̇3 = −2gQ1 − 2Q2 .

(36c)

Starting from


Q1 Q2 Q3



=



0 0 0



,

(37)

we would like to find the time-varying control g(t), bounded as in (33), which minimizes
Q1 (T ) while it drives the other two state variables to the final conditions
Q2 (T ) = Q3 (T ) = 0.

(38)

Note that, since Q1 (T ) = − sinh 2r, its minimization corresponds to maximizing the
squeezing parameter r and the entanglement.
We finally show that for two durations T, T ′ , with T ′ > T , the corresponding
minimum values of Q1 satisfy Q1 (T ′) ≤ Q1 (T ), i.e. the minimum value of Q1 is
a nonincreasing function of the final time, thus the squeezing parameter and the
entanglement are nondecreasing functions of it. Indeed, if g(t) is the optimal control for
0 ≤ t ≤ T , which drives the system at the point (Q1 (T ), 0, 0), then the control
(
g(t), 0 ≤ t ≤ T
,
(39)
g ′(t) =
0, T < t ≤ T ′
drives the system at the same point at the final time t = T ′ . Thus, for the larger duration
T ′ we can obtain at least the value Q1 (T ), corresponding to the minimum for the smaller

Maximizing optomechanical entanglement with optimal control

10

duration T . Consequently, the minimum corresponding to T ′ satisfies Q1 (T ′ ) ≤ Q1 (T ).
This monotonic behavior of the objective function indicates the dual character of our
maximum-entanglement problem and the minimum-time problem stated above. If we
solve a series of maximum-entanglement problems for increasing duration T , then the
solution of the minimum-time problem is the shortest time for which the desired level
of entanglement is obtained.
4. Analysis of the optimal solution
In this section we use optimal control theory [17] to obtain some characteristics of the
optimal solution for the maximum-entanglement problem. The control Hamiltonian for
this problem is defined as
Hc = λ1 Q̇1 + λ2 Q̇2 + λ3 Q̇3
q
= 2g(−λ1 1 + Q21 + Q22 + Q23 + λ1 Q3 − λ3 Q1 ) + 2(λ2 Q3 − λ3 Q2 ),(40)

where the costates λi , i = 1, 2, 3 satisfy the equations
∂Hc
= 2g
λ̇1 = −
∂Q1

λ1 Q1
p

1 + Q21 + Q22 + Q23

+ λ3

!

,

!
λ1 Q2
∂Hc
= 2 gp
+ λ3 ,
λ̇2 = −
∂Q2
1 + Q21 + Q22 + Q23
"
!
#
λ1 Q3
∂Hc
=2 g p
− λ1 − λ2 .
λ̇3 = −
∂Q3
1 + Q21 + Q22 + Q23

(41a)
(41b)
(41c)

Since we want to minimize Q1 (T ), the following terminal condition should be satisfied
by the corresponding costate [17]
λ1 (T ) = 1.

(42)

According to Pontryagin’s Minimum Principle, the optimal control g(t) is chosen
such that the control Hamiltonian is minimized for 0 ≤ t ≤ T . But the control appears
only linearly in Hc , while it is bounded as in (33), thus the optimal choice of g depends
on the sign of the quantity which multiplies it in (40), the so-called switching function
q

2
2
2
Φ = −λ1
1 + Q1 + Q2 + Q3 − Q3 − λ3 Q1 .
(43)

Specifically, the optimal control is


G,
Φ<0

g(t) =
singular, Φ = 0 .

 −G,
Φ>0

(44)

Observe that when the switching function is nonzero, the optimal control takes the
corresponding boundary value which minimizes Hc and we call this a bang pulse. When
the switching function is zero during a time interval of nonzero measure, the control is

Maximizing optomechanical entanglement with optimal control

11

called singular and is calculated using the derivatives of Φ, as we demonstrate in the
next paragraph.
During a time interval where the control is singular, the switching function and its
time derivatives are zero
Φ = Φ̇ = Φ̈ = . . . = 0.
For the problem at hand we find, using the system and costate differential equations
q

2
2
2
Φ=0⇒
1 + Q1 + Q2 + Q3 − Q3 λ1 + Q1 λ3 = 0
(45a)
Φ̇ = 0 ⇒ −Q2 λ1 + Q1 λ2 = 0


q
Φ̈ = 0 ⇒ −Q3 λ1 + g Q3 − 1 + Q21 + Q22 + Q23 λ2 + (Q1 − gQ2 )λ3 = 0

(45b)
(45c)

Now observe that the costates λi cannot be simultaneously zero at any moment, since in
such case the linear homogeneous system of first order differential equations (41a)-(41c)
for the costates would imply λi (t) = 0 for 0 ≤ t ≤ T , violating the terminal condition
(42). This non-triviality of the costates requires the determinant corresponding to the
linear homogeneous system (45a)-(45c) to be zero, which leads to the condition
Q1 = 0.

(46)

When the control is singular, the system evolves in the above surface, which is called
singular surface.
On this surface it is also Q̇1 = 0, and from system equation (36a)
p
p
1 + Q21 + Q22 + Q23 − Q3 ⇒ g = 0, since obviously 1 + Q21 + Q22 + Q23 −
we find g
Q3 > 0. Thus, the singular control is
g(t) = 0.

(47)

From the system equations it is obvious that a singular trajectory corresponds to a
rotation around Q1 -axis on the plane Q1 = 0.
Note that Pontryagin’s Minimum Principle provides necessary conditions for
optimality. In the case of singular control, there is an additional necessary condition, the
so-called generalized Legendre-Clebsch condition [17]. For the minimization problem at
hand it takes the form
 2 

d
∂Hc
1 ∂
≥0
(−1)
∂g dt2 ∂g


q
2
2
2
⇒ λ2 Q3 − 1 + Q1 + Q2 + Q3 − λ3 Q2 ≤ 0.
(48)

We will use this condition in the next section as an extra test of optimality.
We close this section by observing that, although we have identified the possible
values of the optimal control, finding the exact pulse sequence is a difficult problem which
depends on the duration T and the control bound G. For this reason, in the next section
we recourse to numerical optimization in order to obtain the optimal pulse sequences
for specific values of these parameters. But before doing so, here we exclude certain
pulse sequences and this will lead us to the most frequently encountered optimal control

Maximizing optomechanical entanglement with optimal control

12

profiles. First of all, recall that in the previous section we showed that a simple constant
pulse cannot drive the system to a final state of the desired form. Next, observe that
a bang-singular pulse sequence is forbidden because the initial bang pulse cannot reach
the singular surface (46), since during this pulse Q̇1 has constant sign, from (36a) with
g(t) = ±G, while Q1 (0) = 0. The sequence singular-bang is equivalent to a simple bang
pulse, since the singular control is g(t) = 0 and the system remains at the starting point
(the origin) during this part, thus it is also excluded. We examine next a bang-bang
sequence, composed by two bang pulses (±G ∓ G). Observe that, when the duration
T is fixed, there is only one free parameter, the duration of one of the pulses, while
there are two final conditions (38) to be satisfied. Thus, such a pulse sequence cannot
in general be optimal and this may happen only for certain values of T and G. In order
to satisfy the two final conditions (38), a sequence of at least three pulses is necessary,
where the durations of the first two can be taken as free parameters for fixed T . Indeed,
in the next section we present a numerical example where the optimal control is a bangbang-bang sequence. Sequences of the form bang-bang-singular, singular-bang-bang
and bang-singular-bang can be excluded using similar reasoning as above. The simplest
pulse sequence containing a singular arc has the form bang-bang-singular-bang and this
is the optimal pulse sequence that we encountered in most of our simulations.
5. Numerical examples and discussion
In this section, we use the freely available optimal control solver BOCOP [18] to obtain
numerically the optimal solutions for several illustrative examples. In the BOCOP
software package, the continuous-time optimal control problem is approximated by
a finite-dimensional optimization problem, using time discretization. The resultant
nonlinear programming problem is subsequently solved using the nonlinear solver Ipopt.
In the first example that we consider we take G = 1, a maximum coupling equal
to the frequency of the mechanical resonator, and T = π, which is half the period of
this oscillator. We use a time discretization of 20000 points. In Fig. 1(a) we plot
the numerically obtained normalized optimal control g(t)/G. Observe that it has the
form bang-bang-bang. The corresponding maximum value of the squeezing parameter
is r = 1.8990. The accuracy for the boundary constraints is of the order of 10−28 ,
while that for the dynamic constraints (the discretized system equations) of the order
of 10−15 . In Fig. 1(c) we plot the corresponding switching function (43) (blue solid
line) superimposed on the optimal control (blue dashed line). Observe that g = G when
Φ < 0 and g = −G when Φ > 0, in accordance to the Minimum Principle (44). In Fig.
1(e) we plot the optimal trajectory which consists of three segments, each corresponding
to a bang control.
The next case that we consider is with G = 2 and T = π/2. Again, we use a time
discretization of 20000 points. In Fig. 1(b) we plot the normalized optimal control,
which now has the form bang-bang-singular-bang. The maximum value of the squeezing
parameter is r = 0.8336, with accuracies for boundary and dynamic constraints of the

13

1

1

0.5

0.5
g(t)/G

g(t)/G

Maximizing optomechanical entanglement with optimal control

0

0

−0.5
−1
0

−0.5

1

2

−1
0

3

1
1.5
time
(b) Optimal control, G = 2, T = π/2

time
(a) Optimal control, G = 1, T = π

0.5

2
optimality tests

switching function

0
−5
−10

−2
−4
−6

−15
0

0

1

2

−8
0

3

1
1.5
time
(d) Optimality tests, G = 2, T = π/2

time
(c) Swithcing function, G = 1, T = π

0.5

0
2
Q1

−10
Q1

0
−20
−30
10

−2
2

0
5
−10 0
Q3
Q2
(e) Optimal trajectory, G = 1, T = π

10

−4
4

2

1
0

Q2
−2 0
Q3
(f) Optimal trajectory, G = 2, T = π/2

Figure 1. (Color online) Normalized optimal control g(t)/G for (a) G = 1, T = π
and (b) G = 2, T = π/2. (c) Switching function (blue solid line) superimposed on
the optimal control (blue dashed line), for G = 1, T = π. (d) Switching function
(blue solid line) and generalized Legendre-Clebsch condition (red dashed-dotted line)
superimposed on the optimal control (blue dashed line), for G = 2, T = π/2. Optimal
trajectories for (e) G = 1, T = π and (f) G = 2, T = π/2.

order of 10−29 and 10−16 , respectively. Observe that this accuracy is better than in the
previous example, since we use the same discretization for a shorter time interval. In
Fig. 1(d) we plot the switching function (blue solid line) and the left hand side of the

14

Maximizing optomechanical entanglement with optimal control

2.5

r

2
1.5
1
0.5

5

10
G

15

20

Figure 2. Squeezing parameter r as a function of the maximum control amplitude G
for fixed T = π/2.

generalized Legendre-Clebsch condition (48) (red dadhed-dotted line), superimposed on
the optimal control (blue dashed line). Observe that g = G for Φ < 0 and g = −G
for Φ > 0, as in the previous case, while g = 0 when Φ = 0 (singular arc). Note that
on the singular arc the generalized Legendre-Clebsch condition is also satisfied. In Fig.
1(f) we plot the optimal trajectory where four segments are clearly distinguished, each
corresponding to a constant control interval. Observe that the singular third segment
is actually a rotation around Q1 -axis on the plane Q1 = 0. Here we point out that
singular optimal controls appear in the context of quantum control in the contrast
imaging problem of Nuclear Magnetic Resonance [19], while they are quite common in
aerospace applications, see for example the famous Goddard problem of maximizing the
altitude of a vertically ascending rocket with a finite amount of fuel [20].
In Fig. 2 we plot the squeezing parameter r as a function of the control bound
G, for 1 ≤ G ≤ 20 with a step ∆G = 0.1 and for fixed duration T = π/2. For all
these optimizations we have found that the optimal control has the form bang-bangsingular-bang. The percentage of the time occupied by the singular arc is increasing
with increasing G. Finally, in Fig. 3 we plot the squeezing parameter r as a function
of T , for 1 ≤ T ≤ 3 with a step ∆T = 0.1 and for fixed control bound G = 2. Observe
that this last plot can be used to find the necessary time to obtain a desired level
of squeezing-entanglement. From the last two figures it is obvious that in the strong
coupling regime, where G is comparable or larger than the frequency of the mechanical
resonator, a substantial amount of entanglement can be created within a fraction of a
single oscillation period, which is 2π.

15

Maximizing optomechanical entanglement with optimal control

3

r

2

1

0
1

1.5

2
T

2.5

3

Figure 3. Squeezing parameter r as a function of duration T for fixed G = 2.

6. Conclusion
In this paper, we formulated the creation of optomechanical entanglement as an optimal
control problem and used optimal control theory to study the optimal solution. We
subsequently used numerical optimization to obtain the optimal pulse sequences for
several examples. Although here we considered as the initial state of both oscillators
the vacuum, it is our intention to extend the present work to the more practical situation
where the mechanical resonator is initially in a thermal state. Note that in this case
the problem is complicated by the fact that both systems (22) and (23) have nonzero
initial conditions.
References
[1] Aspelmeyer M, Kippenberg T J and Marquardt F 2014 Cavity optomechanics Rev. Mod. Phys. 86
1391
[2] Vitali D, Gigan S, Ferreira A, Böhm H R, Tombesi P, Guerreiro A, Vedral V, Zeilinger A and
Aspelmeyer M 2007 Optomechanical entanglement between a movable mirror and a cavity field
Phys. Rev. Lett. 98 030405
[3] Paternosto M, Vitali D, Gigan S, Kim M S, Brukner C, Eisert J and Aspelmeyer M 2007 Creating
and probing multipartite macroscoping entanglement with light Phys. Rev. Lett. 99 250401
[4] Hofer S G, Wieczorek W, Aspelmeyer M and Hammerer K 2011 Quantum entanglement and
teleportation in pulsed cavity optomechanics Phys. Rev. A 84 052327
[5] Palomaki T A, Teufel J D, Simmonds R W and Lehnert K W 2013 Entangling mechanical motion
with microwave fields Science 342 pp 710–3
[6] Wang X, Vinjanampathy S, Strauch F W and Jacobs K 2011 Ultraefficient cooling of resonators:
beating sideband cooling with quantum control Phys. Rev. Lett. 107 177204

Maximizing optomechanical entanglement with optimal control

16

[7] Machnes S, Cerrillo J, Aspelmeyer M, Wieczorek W, Plenio M B and Retzker A 2012 Pulsed laser
cooling for cavity optomechanical resonators Phys. Rev. Lett. 108 153601
[8] Stefanatos D 2016 Fast cavity optomechanical cooling Automatica 73 71–5
[9] Jacobs K, Wu R, Wang X, Ashhab S, Chen Q.-M and Rabitz H 2016 Fast quantum communication
in linear networks EPL 114 40007
[10] Merzbacher E 1997 Quantum Mechanics (Wiley)
[11] Lvovsky A I 2015 Squeezed light Fundamentals of Photonics and Physics (Photonics vol 1) ed D
L Andrews (West Sussex, United Kingdom: Wiley) pp 121–64
[12] Van Enk S J 1999 Discrete formulation of teleportation of continuous variables Phys. Rev. A 60
5095
[13] Başkal S, Kim Y S and Noz M E 2016 Entangled harmonic oscillators and space-time entanglement
Symmetry 8 55
[14] Dirac P A M 1963 A remarkable representation of the 3 + 2 de Sitter group J. Math. Phys. 4 901
[15] Han D and Kim Y S 1998 Squeezed states as representations of symplectic groups arXiv:
physics/9803017
[16] Bengtsson I 1998 Anti-de Sitter Space (Lecture notes, Stockholm University, Sweden)
[17] Bryson A E and Ho Y C 1975 Applied Optimal Control (New York, U.S.A.: Taylor & Francis)
[18] Bonnans F, Martinon P, Giorgi D, Grélard V, Maindrault S and Tissot O 2015 BOCOP User
Guide, INRIA-Saclay
[19] Bonnard B, Cots O, Glaser S J, Lapert M, Sugny D and Zhang Y 2012 Geometric optimal control
of the contrast imaging problem in Nuclear Magnetic Resonance IEEE Trans. Autom. Control
57 pp 1957–69
[20] Tsiotras P and Kelley H J 1992 Goddard problem with constrained time of flight J. Guid. Control
Dyn. 15 pp 289–96

