- Email: [email protected]
Computational procedure to determine quantum state evolution in Fock space
Computer Physics Communications 180 (2009) 226–230
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Computational procedure to determine quantum state evolution in Fock space D. Portes Jr. a , H. Rodrigues a , B. Baseia b , S.B. Duarte c,∗ a b c
Centro Federal de Educação Tecnológica do Rio de Janeiro, Av. Maracanã 229, CEP 20271-110, Rio de Janeiro, RJ, Brazil Instituto de Física, Universidade Federal de Goiás, PO Box 131, 74.001-970, Goiania, GO, Brazil Centro Brasileiro de Pesquisas Físicas /CNPq, Rua Dr. Xavier Sigaud 150, CEP 22.290-180, Rio de Janeiro, RJ, Brazil
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 13 March 2008 Received in revised form 10 September 2008 Accepted 24 September 2008 Available online 27 September 2008
A new algorithm to solve the quantum state evolution of a system described by a general quadratic Hamiltonian form in creation and the annihilation operators of Fock space is presented. The nonlinear equation for the dynamic operators are obtained in the matrix representation, and by a recursive relation the time evolution operator in the Fock basis is constructed. The method permits to obtain the evolution of entangled quantum states of interacting subsystems when the Hamiltonian of the whole system is in the above mentioned form. Numerical solution with the method is sufficiently accurate to safely analyze the important question of quantum state transfer between the interacting subsystems. A qubits transfer is discussed as an illustrative example when the method is applied to a system described by a particular quadratic Hamiltonian form. © 2008 Elsevier B.V. All rights reserved.
Keywords: Quantum state transfer Fock space Qubits
1. Introduction Many optical and interacting quantum field-matter systems are modeled by a Hamiltonian which can be reduced to a particular quadratic form of the creation and annihilation operators in the Fock space. This happens, for example, when one considers the interaction of a quantized electromagnetic field with a two-level atom inside a cavity. The field states are described by coupled field-modes [2] due to the effect of the atom-field interaction [1]. In the realm of quantum optics and atomic physics this type of model Hamiltonian is widely applied to discuss the transference of quantum state from one subsystem to another along the time evolution of the whole system. The Fock space solution for the time evolution is important to directly recognize when we have a complete‘[3] or only a partial state transfer [4]. This work presents a procedure to solve the time evolution of a state vector of a quantum system described by the most general quadratic form of the creation and annihilation operators, a and a+ [5] for two coupled subsystems, given by H h¯
=
2 j =1
ω j a+j a j +
α ∗j 2
(a+j )2 +
αj 2
(1)
with ω j real and α j , λ, γ being c-numbers. Note that this Hamiltonian involves ten free parameters, which can be even time-
*
Corresponding author. E-mail address: [email protected] (S.B. Duarte).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2008.09.012
©
2008 Elsevier B.V. All rights reserved.
H
=
h¯
1 2
W i j xi x j ,
(2)
ij
where W is a 4 × 4 real matrix and x = {q1 , p1 ; q2 , p2 }, with qi and pi being the canonical coordinate and momentum associated to each oscillator (i = 1, 2). Also, we have [xi , x j ] = i h¯ S i j , with S i j being the elements of matrix,
⎡
0 ⎢ −1 S=⎣ 0 0
1 0 0 0 0 0 0 −1
⎤
0 0⎥ ⎦. 1 0
(3)
The Heisenberg equations of motion are easily obtained and read d dt
a2j
+ ∗ ∗ + + + λa+ 1 a2 + λ a1 a2 + γ a1 a2 + γ a1 a2 ,
dependent and also constituting a model with a rather broad range of practical applicability. The Hamiltonian in Eq. (1) can be written in terms of the Hermitian operators and real parameters as
xi =
1
xi ,
i
H h¯
=
1 2
W jk [xi , x j xk ] =
jk
Ωi j x j ,
(4)
j
where Ω = S W . We can also obtain the operator time evolution from the canonical transformation xi (t ) =
Y i j x j (0).
(5)
j
Using this transformation Eq. (4) can be written in the compact form, dx dt
=
dY dt
x(0) = Ω x(t ) = Ω Y x(0),
(6)
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
leading to, d dt
and
Y = ΩY ,
(7)
with the initial condition Y i j (0) = δi j . Once the canonical transformation Y is obtained the corresponding transformation connecting + the operators a j (0), a+ j (0) to a j (t ), a j (t ), ∗ + a1 (t ) = μ1 a1 (0) + ν1∗ a+ 1 (0) + η1 a2 (0) + ζ1 a2 (0), ∗ +
∗ +
a2 (t ) = μ2 a2 (0) + ν2 a2 (0) + η2 a1 (0) + ζ2 a1 (0),
(8) (9)
solves the system evolution in the Heisenberg picture. The commutation relation [ai , a+ ] = δi j leads to the set of four equations j ∗
∗
∗
∗
1 = μ1 μ1 − ν1 ν1 + η1 η1 − ζ1 ζ1 , ∗
∗
∗
(11)
∗
∗
∗
∗
(12)
0 = μ1 ζ2 − ν1 η2 + η1 ν2 − ζ1 μ2 , 0 = μ1 η2 − ν1 ζ2 + η1 μ2 − ζ1 ν2 , ∗
∗
∗
∗
1 = μ2 μ2 − ν2 ν2 + η2 η2 − ζ2 ζ2 ,
We start from the canonical transformation via the time translation operator Ut connecting a(0) with a(t ), ai (t ) = Ut−1 ai (0)Ut ,
(14)
expressed in the {|n1 , n2 t }-base of the Fock space, constructed with the eigenvectors of the number operator Ni (t ) = a+ i (t )ai (t ), i = 1, 2. We have,
= t n1 , n2 |m1 , m2 0 .
(15)
∗ + t n1 n2 |a1 (t )|m1 m2 0 = t n1 n2 | μ1 a1 (0) + ν1 a1 (0) + + η1 a2 (0) + ζ1∗ a2 (0) |m1m2 0 .
(16)
Since the operator a1 (t ) acts on the base {|n1 , n2 t } whereas a1 (0) acts on the base {|n1 , n2 0 } we get
n1
+
n1 −1,n2 μ1 U m + 1 −1,m2
m1 + 1 ∗ n1 −1,n2 ν1 U m1 +1,m2 n1
m2
n1 −1,n2 η1 U m + 1 ,m2 −1
n1
m2 + 1 ∗ n1 −1,n2 ζ1 U m1 ,m2 +1 . n1
(17)
+ Similarly, repeating the procedure for the operators a+ 1 , a2 and a2 we obtain the three additional relations, n 1, n 2 U m1 ,m2
=
m2 n2
+ n −1,n U m11 ,m2 2
=
n2
n1
+
m2 + 1 ∗ n1 ,n2 −1 ν2 U m1 ,m2 +1 n2
m1
m1
+
μ
n1 ,n2 −1 2 U m1 ,m2 −1
n1 ,n2 −1 2 U m1 −1,m2
η
ν
+
n1 ,n2 1 U m1 −1,m2
+
n1
n1 ,n2 ζ1 U m 1 ,m2 −1
+
+
m2 + 1 ∗ n1 ,n2 μ2 U m1 ,m2 +1 n2
m1 n2
n1 ,n2 ζ2 U m 1 −1,m2
+
m1 + 1 ∗ n1 ,n2 η2 U m1 +1,m2 . n2
(20)
0=
√
0,0
m1 ν1 U m −1,m + 1 2
+ √
√
m1 + 1μ∗1 U m +1,m 1 2
0,0
m2 ζ1 U m ,m −1 + 1 2 0,0
m2 ν2 U m ,m −1 + 1 2
+
√
0,0
m2 + 1η1∗ U m ,m +1 , 1 2
0,0
m2 + 1μ∗2 U m ,m +1 1 2
0,0
m1 ζ2 U m −1,m + 1 2
0,0
m1 + 1η2∗ U m +1,m , 1 2 0,0
(21)
which leads to 0,0
0,0
U 0,1 = U 1,0 = 0,
(22)
for m1 = m2 = 0 and we have taken into account that
μ∗1 μ∗2 − η1∗ η2∗ = 0.
(23)
0=
0=
m1 + 1μ∗1 U m +1,m + 1 2
+
0,0
√
0,0
m1 ν1 U m −1,m + 1 2
+
√
0,0
m1 ζ2 U m −1,m + 1 2
m1 + 1 ∗ n1 ,n2 −1 ζ2 U m1 +1,m2 , n2
m2 + 1 ∗ n1 ,n2 η1 U m1 ,m2 +1 , n1
m2 + 1η1∗ U m ,m +1 1 2
√
0,0
0,0
m2 ζ1 U m ,m −1 , 1 2
m1 + 1η2∗ U m +1,m + 1 2 0,0
√
m2 + 1μ∗2 U m ,m +1 1 2 0,0
0,0
m2 ν2 U m ,m −1 , 1 2
(24)
with the following solutions 0,0 Um 1 ,m2
= =
m2 m1 m1 m2
0,0 δUm 1 −1,m2 −1
+
0,0 δUm 1 −1,m2 −1
+
m1 − 1 m1 m2 − 1 m2
0,0 δ1 U m , 1 −2,m2
(25)
0,0 δ2 U m , 1 ,m2 −2
(26)
in which we have used the constraint relations given by Eq. (12), where
η1∗ ζ2 − ν1 μ∗2 , μ∗1 μ∗2 − η1∗ η2∗ ζ1 η∗ − μ∗1 ν2 δ2 = ∗ 2∗ , μ1 μ2 − η1∗ η2∗ ν1 η∗ − μ∗1 ζ2 δ = ∗ 2∗ . μ1 μ2 − η1∗ η2∗ δ1 =
(27) (28) (29)
Note that Eq. (23) guarantees that these quantities in Eqs. (27)– (29) only assume finite values. By mathematical induction it is 0,0 possible to show that U m1 ,m2 = 0 for odd values of m1 + m2 . 0,0
Eqs. (25)–(36) permit us to determine the elements U m1 m2 when 0,0 U 0,0
(18)
m1 + 1 ∗ n1 ,n2 μ1 U m1 +1,m2 n1
m2
0=
0,0 Um 1 ,m2
ν
Eq. (17) allows us to obtain U n1 ,n2 from U n1 −1,n2 and Eq. (18) determines U n1 ,n2 −1 from U n1 ,n2 . Next taking n1 = 0 and n2 = 0 in Eqs. (19) and (20) we obtain,
The application of Eqs. (8) and (9) furnishes,
m1
+
(Ut )nm11,,nm2 2 = 0 n1 , n2 |U|m1 , m2 0
n ,n
n2
n1 ,n2 2 U m1 ,m2 −1
For m1 = 0 and m2 = 0 Eq. (21) reads,
2. The time translation operator
=
m2
(13)
imposing six constraint relations involving the coefficients μi , νi and ηi , namely, one given by Eq. (10), four given from the real and imaginary parts of Eqs. (11) and (12), with the last one coming from Eq. (13). To obtain the time evolution of the vector state describing our system we must calculate explicitly the time translation operator as detailed in the next section.
n ,n −1 U m11 ,m2 2
(10)
∗
U m11 ,m2 2 =
227
(19)
is known. In resume the recursive relations given by Eqs. (17)–(18) and (25)–(26) permit us to construct the matrix representation of the operator U from the canonical transformation Y obtained from the solution of Eq. (7). Up to here, the method has been developed only for two coupled-subsystem. However it is important to remark that it can be generalized for L interacting subsystems, keeping our basic assumption of a quadratic form involving creation and annihilation
228
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
operators. The starting point is the generalization of Eqs. (8) and (9) given by, ai (t ) =
L
n ,...,n
μik ak (0) + νik∗ ak+ (0); i = 1, . . . , L ,
(30)
k=1
which leads to a set of the L recursive relations, n ,...,n
U m11 ,...,mL L =
L mk k
+
ni
n1 ,...,ni −1,...,n L μik U m 1 ,...,mk −1,...,m L
mk + 1 ∗ n1 ,...,ni −1,...,n L νik U m1 ,...,mk +1,...,m L . ni
(31)
√
0,...,0
mk νik U m ,...,m −1,...,m L 1 k
k
0,...,0 + mk + 1μ∗ik U m , 1 ,...,mk +1,...,m L
(32)
together with the unitariness condition, U † U = U U † = 1. Here we note that Eq. (32) is the analogous of Eq. (21) when dealing with the L-coupled subsystem. From Eqs. (30)–(32) we see that the method can be extended to the case of L-coupled subsystems, via the same procedures for L = 2. For practical use of the method we have to work in a truncated ( N + 1)-dimensional version of a Fock space, where the representation of the U -operator has a total of ( N + 1)2L elements. However, by taking an appropriate time expansion of this operator and considering that n
n
n
2 lim U = δm11 δm2 . . . δmLL ,
(33)
δt →0
it is possible to show (see next section) that in first order approximation in time the only nonnull terms to be determined are the diagonal, the sub- and supra-diagonal elements. To higher order approximation the method offers a systematic procedure to obtain more nonnull off-diagonal elements from the recursive relations (30), which permit us to elaborate a single algorithm for any order of approximation in time. We remark that in the present approximation the number of parameters to be determined in each step of time is similar to that one used in the conventional numerical solutions of the L-coupled differential equations, which is of order ( N + 1) L . As shown in the next sections, the present method is more efficient than the conventional method, requesting a similar machine memory size. For the case of L-coupled subsystems the generalized form of the quadratic Hamiltonian in Eq. (1) is given by, H=
L
ωk ak+ ak +
k=1
+
L
αk∗ 2
(ak+ )2 +
αk 2
ak2
k=1
(36)
in view of the commutation between the total number operator ( N = i a+ a ) and the Hamiltonian. The method makes direct use i i of this condition, allowing us to reach a high computational performance compared with other methods using a similar number of parameters to be determined in each step of time. For simplicity, and to highlight the essential characteristic of the method, hereafter we will mainly deal with the particular case of L = 2 coupled subsystems.
The computational scheme must be established to obtain accurately, step by step, the new values of components of the state vector. As shown in previous section the procedure must hold the 0,0 unitariness of the operator U , determining the values of U 0,0 from the unitary condition U † U = U U † = 1. For a given small value of time step δt one has μ j of order δt 0 and ν j , η j , ζ j just of order δt, as inferred by Eqs. (8), (9). In this case, Eqs. (17), (18), (25) and (26) allow us to write n ,n
U m11 ,m2 2 =
M
δn1 ,m1 +ki δn2 ,m2 +li × O(δt i ),
(37)
i =0
where |k j | + |l j | = 2 j and M is the truncation order for the expansion in δt, with δn1 ,m1 +ki and δn2 ,m2 +li being the Kronecker delta symbol. To estimate the error originated from the numerical method we compare the variance of observable quantities obtained from the xi ), with the correspondnumerical solution of the state vector ( ing analytical solution ( xi ), obtained via the Heisenberg picture. We define the relative error for this observable as
i =
xi | | xi − , xi
(38)
with the xi being the operators appearing in Eq. (2). So, we define the total error by
=
1 4
i .
(39)
i
An illustrative result for this error is calculated here considering the complete form of Hamiltonian in Eq. (1) with a specific choice of parameter values. In Fig. 1 the error is depicted as a function of time, for the time step ω1 δt = 0.5 × 10−2 , the parameter values given in the figure caption and the initial configuration |Ψ (0) = |1 ⊗ |4. The numerical calculation was carried out using a truncated Fock base with a large enough dimension, N = 40. Note that the magnitude of numerical error is affected by the value of M appearing in Eq. (37), as the maximum order in the time expansion of the translation operator. The error is displayed in Fig. 1 for M = 1, 2, 3, 4.
(34)
from which we can construct the generalized recursive relations in Eq. (31). A special subclass of this Hamiltonian models is represented by the following simplified form of Eq. (32), but still significative for physical applications, L
for n1 + · · · + n L = m1 + · · · + m L ,
4. Numerical application
λ i j a+ a j + λ∗i j ai a+ + γi∗j a+ a+ + γi j ai a j , i j i j
i< j
H=
U m11 ,...,mL L = 0,
3. Computational features of the method
Note that Eqs. (17)–(18) are recovered from Eq. (31) for the case L = 2, as it should. These relations allow us to construct the whole time translation operator U t (2L-rank tensor) from the set of el0,...,0 ements U m1 ,...,m L . To determine this set of starting elements we have to solve the following system of L-equations, 0=
In this case the Fock space matrix representation of the operator U reads
ωk ak+ ak +
L i< j
λ i j a+ a j + λ∗i j ai a+ . i j
(35)
It is important to call attention to the convenience of the numerical solution presented here for physical applications. The scheme permits one to discuss the relevant problem of quantum state transfer between coupled HO, dealing with arbitrary initial states. The evolution operator for the vector state associated with the whole system is determined directly in the Fock base and the instant of the state transfer can be also obtained with high precision via numerical solution for general case; an analytical determination of this instant of time is noticeable by using the method for
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
229
Table 1 The values of the quantity (1 − F ), which measures the difference between the numerical and analytical solution of the system state vector calculated at the transfer instant of time, τ = 1/2 (with τ = ωt /π ). M \dτ
10−1
10−2
10−3
1 2 3 4
4.8 × 10−1
6.9 × 10−2
1.7 × 10−5 5.7 × 10−7 4.6 × 10−14
7.4 × 10−3 1.7 × 10−8 5.7 × 10−10 1.0 × 10−14
RK
4.8 × 10−1
1.1 × 10−5
1.1 × 10−10
1.8 × 10−2 6.6 × 10−4 5.2 × 10−8
The results are obtained by using different expansion order in Eq. (37) (M = 1, 2, 3, 4), and for three values of dimensionless time step size. To perform the calculation we have taken ω1 = ω2 = ω = 2λ for the Hamiltonian in Eq. (40). In the last line we display the result obtained with a conventional fourth order Runge– Kutta (RK) method for the same situation.
Fig. 1. Time evolution of the numerical error for the variance of an observable operator calculated by using different expansion order in Eq. (37) ( M = 1, 2, 3, 4). The size of the time step is specified in the figure, and the results are obtained from the general quadratic form of the Hamiltonian in Eq. (1), by taking the arbitrary choice ω1 = 1.835, ω2 = 2.150, α1 = 0.255 + 0.73i, α2 = 0.020 − 1.120i, and γ = 0.290 − 0.210i, λ = 0.200 − 0.550i.
some particular cases. We have discussed previously in Ref. [4] a comparable situation where the statistical distributions characterizing two interacting subsystems are transferred from one to the other. Also, complete transfer of states between two subsystems was studied in a sequential work [3]. Here we are focusing our attention only to the numerical scheme and its implementation to illustrate the accuracy of the numerical method to obtain the time evolution of the state. With this purpose we have chosen a particular quadratic form of the Hamiltonian for which it is possible to get exact solution of the evolving state, in order to compare it with the numerical solution. The following particular Hamiltonian was chosen H h¯
+ + + = ω1 a+ 1 a1 + ω2 a2 a2 + λ a1 a2 + a1 a2 .
(40)
The time evolution of the system is solved analytically and numerically, starting from the initial state,
Ψ (0) = |φ ⊗ |0.
(41)
It is possible to demonstrate that transfer occurs at the instant of time t = π /(2λ), taking ω1 = ω2 = ω = 2λ in the Hamiltonian (40) (see Ref. [3]), when we have
Ψ π /(2λ) = |0 ⊗ |φ,
(42)
where
|φ =
cn |n,
(43)
n
with cn = |cn |e −i αn and αn = ( ω λ + 1)nπ . The efficacy of the method can be tested by calculating the absolute value of the fidelity parameter,
2
F = Ψanal |Ψnum ,
(44)
associated to the analytical and numerical state vector solutions. Table 1 shows the values for the quantity (1 − F ), obtained at the
transfer instant of time, τ = 1/2 (with τ = ωt /π ), using the numerical solution with different order in time expansion Eq. (37) ( M = 1, 2, 3, 4) for three different values of dimensionless time step size, dτ . We remark that one of the significative computational advantages of the method consists in a quite general important characteristic of the method, its algorithm implementation allowing us to increase the precision of the results by only changing the value to the variable-M appearing in Eq. (37) (see also Appendix A). The conventional methods are associated with the solution of the differential equation and, to increase the accuracy, it is necessary to redefine the algorithm, by elaborating a new computational code involving different mathematical expressions and values of parameters. A manifest advantage of using the Fock space, consists in having a computational matrix representation of the evolution operator generated directly by the discrete base of the space. Whereas for the coordinate space or phase space (continuous base) the matrix computational representation of the propagator involves an ad-hoc space discretization procedure with mere numerical meaning. 5. Conclusions and final remarks We have shown the possibility, and comparative advantages, of extending the method for the case involving large number of coupled subsystems. The extended method remains a valuable tool to discuss the quantum state transfer problem between coupled systems in the quite important class of Hamiltonian model described by (34). The state transfer between subsystems depends on the particular adopted Hamiltonian model in this general class, and it constitutes a subject of a future work. The more recent focus of interest concerns the transfer of entangled states in multipartite subsystems. Here we mention a still simple case, with L = 4, which can be used as a groundwork to give some insight to more general situations. Let us consider the system of a linear chain of four coupled harmonic oscillators, with its elements ordered as A 1 , + + + A 2 , A 3 , and A 4 , with (a+ 1 , a1 ), (a2 , a2 ), (a3 , a3 ), and (a4 , a4 ) being their corresponding creation and annihilation operators. Now, assume that one is concerned with the transfer of an initial (entangled) state of the composed subsystem A 1 A 2 , to the composed subsystem A 3 A 4 , by considering the interaction between two adjacent elements, namely: A 1 ⇔ A 2 , A 2 ⇔ A 3 , and A 3 ⇔ A 4 , via the interacting Hamiltonian,
+ + + H int = λ12 a1 a+ 2 + a1 a2 + λ23 a2 a3 + a2 a3
+ + λ34 a3 a+ 4 + a3 a4 .
(45)
Then, for example, if the Hamiltonian leads the initial state
Ψ (0) A
1 A2 A3 A4
= |φ A 1 A 2 ⊗ |0 A 3 A 4 ,
(46)
230
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
of the entire system A 1 A 2 A 3 A 4 , to the final state
Ψ (t )
A1 A2 A3 A4
= |0 A 1 A 2 ⊗ |φ A 3 A 4 ,
(47)
this would represent a (mutual) transfer of states between the composed subsystems A 1 A 2 and A 3 A 4 . The Hamiltonian in Eq. (45) is another particular case of the general form used to formulate the present method (Eq. (34)), for convenient choice of the parameters, e.g., λ12 , λ23 , λ23 = 0 (real numbers); λi j = 0 and γi j = 0 for other i j-pairs; ωk = αk = 0; for all values of k. Consequently, it is possible to discuss, e.g., the important question of the transfer of states using the method in this still simple but illustrative situation, with the additional help of the condition in Eq. (36). In Ref. [3] we have presented a detailed study of state transfer for two coupled subsystems ( L = 2) described by the Hamiltonian in Eq. (40). The detailed treatment for multipartite systems ( L > 2) and its applications to discuss quantum state transfer are in progress, an issue of future work. Acknowledgements
obtain other elements of v, with the constraints m1 + m2 2M. STEP 6: Reset v, through the relation 1 = m1 m2 | v m1 ,m2 |2 . STEP 7: Set n1 = 0, n2 = 1, and u = v. OUT IN STEP 8: Obtain C 00 = m1 m2 u m1 m2 C m . 1 m2 STEP 9: Set u NEW through the relation
ukNEW 1 ,k2
=
+
d Y = ΩY . STEP 2: Solve the differential equation dt STEP 3: Obtain the parameters of the transformation, Eqs. (27)– (29). STEP 4: Set v 0,0 = 1, v m1 ,m2 = 0, for odd values of m1 + m2 . STEP 5: From the relation
v m1 ,m2 =
m1 m2
v m1 −1,m2 −1 δ +
v m1 −1,m2 −1 δ +
n2
η2 uk1 −1,k2 +1 +
k1 + n1 + 1 ∗ ζ2 uk1 +1,k2 +1 . n2
k1 + n1 n1
+
μ1 v k1 ,k2 +
k2 + n2 n1
k1 + n1 + 1 ∗ ν1 v k1 +2,k2 n1
η1 v k1 +1,k2 −1 +
k2 + n2 + 1 ∗ ζ1 v k1 +1,k2 +1 . n1 1
k1 + n1
k2 + n2 + 1 ∗ ν2 uk1 ,k2 +2 n2
STEP 13: Set v = u and n2 = 1 to obtain C nOUT 0 =
INPUT: |Ψ (t = 0) = n1 n2 C nIN1 ,n2 |n1 n2 . STEP 1: Set Y i j (0) = δi j .
m1
μ2 uk1 ,k2 +
uk1 ,k2 =
Sequence of commands:
m2
STEP 10: Set u = u NEW and obtain C nOUT = k1 k2 uk1 k2 C kIN1 +n1 ,k2 +n2 . 1 ,n2 STEP 11: Set n2 = n2 + 1, repeat steps 9 and 10 until n2 = N. STEP 12: Set n1 = n1 + 1 and n2 = 0, to get u through the relation
Appendix A. Computational code scheme
n2
The authors thank the CNPq (S.B.D., B.B.) and FAPERJ (D.P.J.) for the partial supports.
v m1 ,m2 =
k2 + n2
m1 − 1 m1 m2 − 1 m2
v m1 −2,m2 δ1 ,
v m1 ,m2 −2 δ2 ,
k 1 m2
u k 1 m2 ·
C kIN+n ,m . 1 1 2 STEP 14: Repeat steps 9, 10 and 11 until n2 = N. STEP 15: Set n1 = n1 + 1, repeat steps 12, 13 and 14 until n1 = N. OUTPUT: |Ψ (t + δt ) = n1 n2 C nOUT |n1 n2 . 1 ,n2 References [1] A.L. Fetter, J.D. Walecka, Quantum Theory—Many-Body Problems, International Series in Pure and Applied Physics, McGraw-Hill Publishing Co., New York, 1971. [2] L. Mandel, E. Wolf, Quantum Coherence and Quantum Optics, Cambridge Univ. Press, 1995. [3] D. Portes Jr., et al., Eur. Phys. J. D 48 (2008) 145–149. [4] H. Rodrigues, et al., Physica A 311 (2002) 188; D. Portes Jr., et al., Physica A 329 (2003) 391. [5] See, e.g., S.-W. Tsai, A.F.R. de Toledo Piza, Phys. Rev. A 53 (1996) 3683; K.M. Ng, C.F. Lo, Phys. Lett. A 230 (1997).
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Computational procedure to determine quantum state evolution in Fock space D. Portes Jr. a , H. Rodrigues a , B. Baseia b , S.B. Duarte c,∗ a b c
Centro Federal de Educação Tecnológica do Rio de Janeiro, Av. Maracanã 229, CEP 20271-110, Rio de Janeiro, RJ, Brazil Instituto de Física, Universidade Federal de Goiás, PO Box 131, 74.001-970, Goiania, GO, Brazil Centro Brasileiro de Pesquisas Físicas /CNPq, Rua Dr. Xavier Sigaud 150, CEP 22.290-180, Rio de Janeiro, RJ, Brazil
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 13 March 2008 Received in revised form 10 September 2008 Accepted 24 September 2008 Available online 27 September 2008
A new algorithm to solve the quantum state evolution of a system described by a general quadratic Hamiltonian form in creation and the annihilation operators of Fock space is presented. The nonlinear equation for the dynamic operators are obtained in the matrix representation, and by a recursive relation the time evolution operator in the Fock basis is constructed. The method permits to obtain the evolution of entangled quantum states of interacting subsystems when the Hamiltonian of the whole system is in the above mentioned form. Numerical solution with the method is sufficiently accurate to safely analyze the important question of quantum state transfer between the interacting subsystems. A qubits transfer is discussed as an illustrative example when the method is applied to a system described by a particular quadratic Hamiltonian form. © 2008 Elsevier B.V. All rights reserved.
Keywords: Quantum state transfer Fock space Qubits
1. Introduction Many optical and interacting quantum field-matter systems are modeled by a Hamiltonian which can be reduced to a particular quadratic form of the creation and annihilation operators in the Fock space. This happens, for example, when one considers the interaction of a quantized electromagnetic field with a two-level atom inside a cavity. The field states are described by coupled field-modes [2] due to the effect of the atom-field interaction [1]. In the realm of quantum optics and atomic physics this type of model Hamiltonian is widely applied to discuss the transference of quantum state from one subsystem to another along the time evolution of the whole system. The Fock space solution for the time evolution is important to directly recognize when we have a complete‘[3] or only a partial state transfer [4]. This work presents a procedure to solve the time evolution of a state vector of a quantum system described by the most general quadratic form of the creation and annihilation operators, a and a+ [5] for two coupled subsystems, given by H h¯
=
2 j =1
ω j a+j a j +
α ∗j 2
(a+j )2 +
αj 2
(1)
with ω j real and α j , λ, γ being c-numbers. Note that this Hamiltonian involves ten free parameters, which can be even time-
*
Corresponding author. E-mail address: [email protected] (S.B. Duarte).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2008.09.012
©
2008 Elsevier B.V. All rights reserved.
H
=
h¯
1 2
W i j xi x j ,
(2)
ij
where W is a 4 × 4 real matrix and x = {q1 , p1 ; q2 , p2 }, with qi and pi being the canonical coordinate and momentum associated to each oscillator (i = 1, 2). Also, we have [xi , x j ] = i h¯ S i j , with S i j being the elements of matrix,
⎡
0 ⎢ −1 S=⎣ 0 0
1 0 0 0 0 0 0 −1
⎤
0 0⎥ ⎦. 1 0
(3)
The Heisenberg equations of motion are easily obtained and read d dt
a2j
+ ∗ ∗ + + + λa+ 1 a2 + λ a1 a2 + γ a1 a2 + γ a1 a2 ,
dependent and also constituting a model with a rather broad range of practical applicability. The Hamiltonian in Eq. (1) can be written in terms of the Hermitian operators and real parameters as
xi =
1
xi ,
i
H h¯
=
1 2
W jk [xi , x j xk ] =
jk
Ωi j x j ,
(4)
j
where Ω = S W . We can also obtain the operator time evolution from the canonical transformation xi (t ) =
Y i j x j (0).
(5)
j
Using this transformation Eq. (4) can be written in the compact form, dx dt
=
dY dt
x(0) = Ω x(t ) = Ω Y x(0),
(6)
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
leading to, d dt
and
Y = ΩY ,
(7)
with the initial condition Y i j (0) = δi j . Once the canonical transformation Y is obtained the corresponding transformation connecting + the operators a j (0), a+ j (0) to a j (t ), a j (t ), ∗ + a1 (t ) = μ1 a1 (0) + ν1∗ a+ 1 (0) + η1 a2 (0) + ζ1 a2 (0), ∗ +
∗ +
a2 (t ) = μ2 a2 (0) + ν2 a2 (0) + η2 a1 (0) + ζ2 a1 (0),
(8) (9)
solves the system evolution in the Heisenberg picture. The commutation relation [ai , a+ ] = δi j leads to the set of four equations j ∗
∗
∗
∗
1 = μ1 μ1 − ν1 ν1 + η1 η1 − ζ1 ζ1 , ∗
∗
∗
(11)
∗
∗
∗
∗
(12)
0 = μ1 ζ2 − ν1 η2 + η1 ν2 − ζ1 μ2 , 0 = μ1 η2 − ν1 ζ2 + η1 μ2 − ζ1 ν2 , ∗
∗
∗
∗
1 = μ2 μ2 − ν2 ν2 + η2 η2 − ζ2 ζ2 ,
We start from the canonical transformation via the time translation operator Ut connecting a(0) with a(t ), ai (t ) = Ut−1 ai (0)Ut ,
(14)
expressed in the {|n1 , n2 t }-base of the Fock space, constructed with the eigenvectors of the number operator Ni (t ) = a+ i (t )ai (t ), i = 1, 2. We have,
= t n1 , n2 |m1 , m2 0 .
(15)
∗ + t n1 n2 |a1 (t )|m1 m2 0 = t n1 n2 | μ1 a1 (0) + ν1 a1 (0) + + η1 a2 (0) + ζ1∗ a2 (0) |m1m2 0 .
(16)
Since the operator a1 (t ) acts on the base {|n1 , n2 t } whereas a1 (0) acts on the base {|n1 , n2 0 } we get
n1
+
n1 −1,n2 μ1 U m + 1 −1,m2
m1 + 1 ∗ n1 −1,n2 ν1 U m1 +1,m2 n1
m2
n1 −1,n2 η1 U m + 1 ,m2 −1
n1
m2 + 1 ∗ n1 −1,n2 ζ1 U m1 ,m2 +1 . n1
(17)
+ Similarly, repeating the procedure for the operators a+ 1 , a2 and a2 we obtain the three additional relations, n 1, n 2 U m1 ,m2
=
m2 n2
+ n −1,n U m11 ,m2 2
=
n2
n1
+
m2 + 1 ∗ n1 ,n2 −1 ν2 U m1 ,m2 +1 n2
m1
m1
+
μ
n1 ,n2 −1 2 U m1 ,m2 −1
n1 ,n2 −1 2 U m1 −1,m2
η
ν
+
n1 ,n2 1 U m1 −1,m2
+
n1
n1 ,n2 ζ1 U m 1 ,m2 −1
+
+
m2 + 1 ∗ n1 ,n2 μ2 U m1 ,m2 +1 n2
m1 n2
n1 ,n2 ζ2 U m 1 −1,m2
+
m1 + 1 ∗ n1 ,n2 η2 U m1 +1,m2 . n2
(20)
0=
√
0,0
m1 ν1 U m −1,m + 1 2
+ √
√
m1 + 1μ∗1 U m +1,m 1 2
0,0
m2 ζ1 U m ,m −1 + 1 2 0,0
m2 ν2 U m ,m −1 + 1 2
+
√
0,0
m2 + 1η1∗ U m ,m +1 , 1 2
0,0
m2 + 1μ∗2 U m ,m +1 1 2
0,0
m1 ζ2 U m −1,m + 1 2
0,0
m1 + 1η2∗ U m +1,m , 1 2 0,0
(21)
which leads to 0,0
0,0
U 0,1 = U 1,0 = 0,
(22)
for m1 = m2 = 0 and we have taken into account that
μ∗1 μ∗2 − η1∗ η2∗ = 0.
(23)
0=
0=
m1 + 1μ∗1 U m +1,m + 1 2
+
0,0
√
0,0
m1 ν1 U m −1,m + 1 2
+
√
0,0
m1 ζ2 U m −1,m + 1 2
m1 + 1 ∗ n1 ,n2 −1 ζ2 U m1 +1,m2 , n2
m2 + 1 ∗ n1 ,n2 η1 U m1 ,m2 +1 , n1
m2 + 1η1∗ U m ,m +1 1 2
√
0,0
0,0
m2 ζ1 U m ,m −1 , 1 2
m1 + 1η2∗ U m +1,m + 1 2 0,0
√
m2 + 1μ∗2 U m ,m +1 1 2 0,0
0,0
m2 ν2 U m ,m −1 , 1 2
(24)
with the following solutions 0,0 Um 1 ,m2
= =
m2 m1 m1 m2
0,0 δUm 1 −1,m2 −1
+
0,0 δUm 1 −1,m2 −1
+
m1 − 1 m1 m2 − 1 m2
0,0 δ1 U m , 1 −2,m2
(25)
0,0 δ2 U m , 1 ,m2 −2
(26)
in which we have used the constraint relations given by Eq. (12), where
η1∗ ζ2 − ν1 μ∗2 , μ∗1 μ∗2 − η1∗ η2∗ ζ1 η∗ − μ∗1 ν2 δ2 = ∗ 2∗ , μ1 μ2 − η1∗ η2∗ ν1 η∗ − μ∗1 ζ2 δ = ∗ 2∗ . μ1 μ2 − η1∗ η2∗ δ1 =
(27) (28) (29)
Note that Eq. (23) guarantees that these quantities in Eqs. (27)– (29) only assume finite values. By mathematical induction it is 0,0 possible to show that U m1 ,m2 = 0 for odd values of m1 + m2 . 0,0
Eqs. (25)–(36) permit us to determine the elements U m1 m2 when 0,0 U 0,0
(18)
m1 + 1 ∗ n1 ,n2 μ1 U m1 +1,m2 n1
m2
0=
0,0 Um 1 ,m2
ν
Eq. (17) allows us to obtain U n1 ,n2 from U n1 −1,n2 and Eq. (18) determines U n1 ,n2 −1 from U n1 ,n2 . Next taking n1 = 0 and n2 = 0 in Eqs. (19) and (20) we obtain,
The application of Eqs. (8) and (9) furnishes,
m1
+
(Ut )nm11,,nm2 2 = 0 n1 , n2 |U|m1 , m2 0
n ,n
n2
n1 ,n2 2 U m1 ,m2 −1
For m1 = 0 and m2 = 0 Eq. (21) reads,
2. The time translation operator
=
m2
(13)
imposing six constraint relations involving the coefficients μi , νi and ηi , namely, one given by Eq. (10), four given from the real and imaginary parts of Eqs. (11) and (12), with the last one coming from Eq. (13). To obtain the time evolution of the vector state describing our system we must calculate explicitly the time translation operator as detailed in the next section.
n ,n −1 U m11 ,m2 2
(10)
∗
U m11 ,m2 2 =
227
(19)
is known. In resume the recursive relations given by Eqs. (17)–(18) and (25)–(26) permit us to construct the matrix representation of the operator U from the canonical transformation Y obtained from the solution of Eq. (7). Up to here, the method has been developed only for two coupled-subsystem. However it is important to remark that it can be generalized for L interacting subsystems, keeping our basic assumption of a quadratic form involving creation and annihilation
228
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
operators. The starting point is the generalization of Eqs. (8) and (9) given by, ai (t ) =
L
n ,...,n
μik ak (0) + νik∗ ak+ (0); i = 1, . . . , L ,
(30)
k=1
which leads to a set of the L recursive relations, n ,...,n
U m11 ,...,mL L =
L mk k
+
ni
n1 ,...,ni −1,...,n L μik U m 1 ,...,mk −1,...,m L
mk + 1 ∗ n1 ,...,ni −1,...,n L νik U m1 ,...,mk +1,...,m L . ni
(31)
√
0,...,0
mk νik U m ,...,m −1,...,m L 1 k
k
0,...,0 + mk + 1μ∗ik U m , 1 ,...,mk +1,...,m L
(32)
together with the unitariness condition, U † U = U U † = 1. Here we note that Eq. (32) is the analogous of Eq. (21) when dealing with the L-coupled subsystem. From Eqs. (30)–(32) we see that the method can be extended to the case of L-coupled subsystems, via the same procedures for L = 2. For practical use of the method we have to work in a truncated ( N + 1)-dimensional version of a Fock space, where the representation of the U -operator has a total of ( N + 1)2L elements. However, by taking an appropriate time expansion of this operator and considering that n
n
n
2 lim U = δm11 δm2 . . . δmLL ,
(33)
δt →0
it is possible to show (see next section) that in first order approximation in time the only nonnull terms to be determined are the diagonal, the sub- and supra-diagonal elements. To higher order approximation the method offers a systematic procedure to obtain more nonnull off-diagonal elements from the recursive relations (30), which permit us to elaborate a single algorithm for any order of approximation in time. We remark that in the present approximation the number of parameters to be determined in each step of time is similar to that one used in the conventional numerical solutions of the L-coupled differential equations, which is of order ( N + 1) L . As shown in the next sections, the present method is more efficient than the conventional method, requesting a similar machine memory size. For the case of L-coupled subsystems the generalized form of the quadratic Hamiltonian in Eq. (1) is given by, H=
L
ωk ak+ ak +
k=1
+
L
αk∗ 2
(ak+ )2 +
αk 2
ak2
k=1
(36)
in view of the commutation between the total number operator ( N = i a+ a ) and the Hamiltonian. The method makes direct use i i of this condition, allowing us to reach a high computational performance compared with other methods using a similar number of parameters to be determined in each step of time. For simplicity, and to highlight the essential characteristic of the method, hereafter we will mainly deal with the particular case of L = 2 coupled subsystems.
The computational scheme must be established to obtain accurately, step by step, the new values of components of the state vector. As shown in previous section the procedure must hold the 0,0 unitariness of the operator U , determining the values of U 0,0 from the unitary condition U † U = U U † = 1. For a given small value of time step δt one has μ j of order δt 0 and ν j , η j , ζ j just of order δt, as inferred by Eqs. (8), (9). In this case, Eqs. (17), (18), (25) and (26) allow us to write n ,n
U m11 ,m2 2 =
M
δn1 ,m1 +ki δn2 ,m2 +li × O(δt i ),
(37)
i =0
where |k j | + |l j | = 2 j and M is the truncation order for the expansion in δt, with δn1 ,m1 +ki and δn2 ,m2 +li being the Kronecker delta symbol. To estimate the error originated from the numerical method we compare the variance of observable quantities obtained from the xi ), with the correspondnumerical solution of the state vector ( ing analytical solution ( xi ), obtained via the Heisenberg picture. We define the relative error for this observable as
i =
xi | | xi − , xi
(38)
with the xi being the operators appearing in Eq. (2). So, we define the total error by
=
1 4
i .
(39)
i
An illustrative result for this error is calculated here considering the complete form of Hamiltonian in Eq. (1) with a specific choice of parameter values. In Fig. 1 the error is depicted as a function of time, for the time step ω1 δt = 0.5 × 10−2 , the parameter values given in the figure caption and the initial configuration |Ψ (0) = |1 ⊗ |4. The numerical calculation was carried out using a truncated Fock base with a large enough dimension, N = 40. Note that the magnitude of numerical error is affected by the value of M appearing in Eq. (37), as the maximum order in the time expansion of the translation operator. The error is displayed in Fig. 1 for M = 1, 2, 3, 4.
(34)
from which we can construct the generalized recursive relations in Eq. (31). A special subclass of this Hamiltonian models is represented by the following simplified form of Eq. (32), but still significative for physical applications, L
for n1 + · · · + n L = m1 + · · · + m L ,
4. Numerical application
λ i j a+ a j + λ∗i j ai a+ + γi∗j a+ a+ + γi j ai a j , i j i j
i< j
H=
U m11 ,...,mL L = 0,
3. Computational features of the method
Note that Eqs. (17)–(18) are recovered from Eq. (31) for the case L = 2, as it should. These relations allow us to construct the whole time translation operator U t (2L-rank tensor) from the set of el0,...,0 ements U m1 ,...,m L . To determine this set of starting elements we have to solve the following system of L-equations, 0=
In this case the Fock space matrix representation of the operator U reads
ωk ak+ ak +
L i< j
λ i j a+ a j + λ∗i j ai a+ . i j
(35)
It is important to call attention to the convenience of the numerical solution presented here for physical applications. The scheme permits one to discuss the relevant problem of quantum state transfer between coupled HO, dealing with arbitrary initial states. The evolution operator for the vector state associated with the whole system is determined directly in the Fock base and the instant of the state transfer can be also obtained with high precision via numerical solution for general case; an analytical determination of this instant of time is noticeable by using the method for
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
229
Table 1 The values of the quantity (1 − F ), which measures the difference between the numerical and analytical solution of the system state vector calculated at the transfer instant of time, τ = 1/2 (with τ = ωt /π ). M \dτ
10−1
10−2
10−3
1 2 3 4
4.8 × 10−1
6.9 × 10−2
1.7 × 10−5 5.7 × 10−7 4.6 × 10−14
7.4 × 10−3 1.7 × 10−8 5.7 × 10−10 1.0 × 10−14
RK
4.8 × 10−1
1.1 × 10−5
1.1 × 10−10
1.8 × 10−2 6.6 × 10−4 5.2 × 10−8
The results are obtained by using different expansion order in Eq. (37) (M = 1, 2, 3, 4), and for three values of dimensionless time step size. To perform the calculation we have taken ω1 = ω2 = ω = 2λ for the Hamiltonian in Eq. (40). In the last line we display the result obtained with a conventional fourth order Runge– Kutta (RK) method for the same situation.
Fig. 1. Time evolution of the numerical error for the variance of an observable operator calculated by using different expansion order in Eq. (37) ( M = 1, 2, 3, 4). The size of the time step is specified in the figure, and the results are obtained from the general quadratic form of the Hamiltonian in Eq. (1), by taking the arbitrary choice ω1 = 1.835, ω2 = 2.150, α1 = 0.255 + 0.73i, α2 = 0.020 − 1.120i, and γ = 0.290 − 0.210i, λ = 0.200 − 0.550i.
some particular cases. We have discussed previously in Ref. [4] a comparable situation where the statistical distributions characterizing two interacting subsystems are transferred from one to the other. Also, complete transfer of states between two subsystems was studied in a sequential work [3]. Here we are focusing our attention only to the numerical scheme and its implementation to illustrate the accuracy of the numerical method to obtain the time evolution of the state. With this purpose we have chosen a particular quadratic form of the Hamiltonian for which it is possible to get exact solution of the evolving state, in order to compare it with the numerical solution. The following particular Hamiltonian was chosen H h¯
+ + + = ω1 a+ 1 a1 + ω2 a2 a2 + λ a1 a2 + a1 a2 .
(40)
The time evolution of the system is solved analytically and numerically, starting from the initial state,
Ψ (0) = |φ ⊗ |0.
(41)
It is possible to demonstrate that transfer occurs at the instant of time t = π /(2λ), taking ω1 = ω2 = ω = 2λ in the Hamiltonian (40) (see Ref. [3]), when we have
Ψ π /(2λ) = |0 ⊗ |φ,
(42)
where
|φ =
cn |n,
(43)
n
with cn = |cn |e −i αn and αn = ( ω λ + 1)nπ . The efficacy of the method can be tested by calculating the absolute value of the fidelity parameter,
2
F = Ψanal |Ψnum ,
(44)
associated to the analytical and numerical state vector solutions. Table 1 shows the values for the quantity (1 − F ), obtained at the
transfer instant of time, τ = 1/2 (with τ = ωt /π ), using the numerical solution with different order in time expansion Eq. (37) ( M = 1, 2, 3, 4) for three different values of dimensionless time step size, dτ . We remark that one of the significative computational advantages of the method consists in a quite general important characteristic of the method, its algorithm implementation allowing us to increase the precision of the results by only changing the value to the variable-M appearing in Eq. (37) (see also Appendix A). The conventional methods are associated with the solution of the differential equation and, to increase the accuracy, it is necessary to redefine the algorithm, by elaborating a new computational code involving different mathematical expressions and values of parameters. A manifest advantage of using the Fock space, consists in having a computational matrix representation of the evolution operator generated directly by the discrete base of the space. Whereas for the coordinate space or phase space (continuous base) the matrix computational representation of the propagator involves an ad-hoc space discretization procedure with mere numerical meaning. 5. Conclusions and final remarks We have shown the possibility, and comparative advantages, of extending the method for the case involving large number of coupled subsystems. The extended method remains a valuable tool to discuss the quantum state transfer problem between coupled systems in the quite important class of Hamiltonian model described by (34). The state transfer between subsystems depends on the particular adopted Hamiltonian model in this general class, and it constitutes a subject of a future work. The more recent focus of interest concerns the transfer of entangled states in multipartite subsystems. Here we mention a still simple case, with L = 4, which can be used as a groundwork to give some insight to more general situations. Let us consider the system of a linear chain of four coupled harmonic oscillators, with its elements ordered as A 1 , + + + A 2 , A 3 , and A 4 , with (a+ 1 , a1 ), (a2 , a2 ), (a3 , a3 ), and (a4 , a4 ) being their corresponding creation and annihilation operators. Now, assume that one is concerned with the transfer of an initial (entangled) state of the composed subsystem A 1 A 2 , to the composed subsystem A 3 A 4 , by considering the interaction between two adjacent elements, namely: A 1 ⇔ A 2 , A 2 ⇔ A 3 , and A 3 ⇔ A 4 , via the interacting Hamiltonian,
+ + + H int = λ12 a1 a+ 2 + a1 a2 + λ23 a2 a3 + a2 a3
+ + λ34 a3 a+ 4 + a3 a4 .
(45)
Then, for example, if the Hamiltonian leads the initial state
Ψ (0) A
1 A2 A3 A4
= |φ A 1 A 2 ⊗ |0 A 3 A 4 ,
(46)
230
D. Portes Jr. et al. / Computer Physics Communications 180 (2009) 226–230
of the entire system A 1 A 2 A 3 A 4 , to the final state
Ψ (t )
A1 A2 A3 A4
= |0 A 1 A 2 ⊗ |φ A 3 A 4 ,
(47)
this would represent a (mutual) transfer of states between the composed subsystems A 1 A 2 and A 3 A 4 . The Hamiltonian in Eq. (45) is another particular case of the general form used to formulate the present method (Eq. (34)), for convenient choice of the parameters, e.g., λ12 , λ23 , λ23 = 0 (real numbers); λi j = 0 and γi j = 0 for other i j-pairs; ωk = αk = 0; for all values of k. Consequently, it is possible to discuss, e.g., the important question of the transfer of states using the method in this still simple but illustrative situation, with the additional help of the condition in Eq. (36). In Ref. [3] we have presented a detailed study of state transfer for two coupled subsystems ( L = 2) described by the Hamiltonian in Eq. (40). The detailed treatment for multipartite systems ( L > 2) and its applications to discuss quantum state transfer are in progress, an issue of future work. Acknowledgements
obtain other elements of v, with the constraints m1 + m2 2M. STEP 6: Reset v, through the relation 1 = m1 m2 | v m1 ,m2 |2 . STEP 7: Set n1 = 0, n2 = 1, and u = v. OUT IN STEP 8: Obtain C 00 = m1 m2 u m1 m2 C m . 1 m2 STEP 9: Set u NEW through the relation
ukNEW 1 ,k2
=
+
d Y = ΩY . STEP 2: Solve the differential equation dt STEP 3: Obtain the parameters of the transformation, Eqs. (27)– (29). STEP 4: Set v 0,0 = 1, v m1 ,m2 = 0, for odd values of m1 + m2 . STEP 5: From the relation
v m1 ,m2 =
m1 m2
v m1 −1,m2 −1 δ +
v m1 −1,m2 −1 δ +
n2
η2 uk1 −1,k2 +1 +
k1 + n1 + 1 ∗ ζ2 uk1 +1,k2 +1 . n2
k1 + n1 n1
+
μ1 v k1 ,k2 +
k2 + n2 n1
k1 + n1 + 1 ∗ ν1 v k1 +2,k2 n1
η1 v k1 +1,k2 −1 +
k2 + n2 + 1 ∗ ζ1 v k1 +1,k2 +1 . n1 1
k1 + n1
k2 + n2 + 1 ∗ ν2 uk1 ,k2 +2 n2
STEP 13: Set v = u and n2 = 1 to obtain C nOUT 0 =
INPUT: |Ψ (t = 0) = n1 n2 C nIN1 ,n2 |n1 n2 . STEP 1: Set Y i j (0) = δi j .
m1
μ2 uk1 ,k2 +
uk1 ,k2 =
Sequence of commands:
m2
STEP 10: Set u = u NEW and obtain C nOUT = k1 k2 uk1 k2 C kIN1 +n1 ,k2 +n2 . 1 ,n2 STEP 11: Set n2 = n2 + 1, repeat steps 9 and 10 until n2 = N. STEP 12: Set n1 = n1 + 1 and n2 = 0, to get u through the relation
Appendix A. Computational code scheme
n2
The authors thank the CNPq (S.B.D., B.B.) and FAPERJ (D.P.J.) for the partial supports.
v m1 ,m2 =
k2 + n2
m1 − 1 m1 m2 − 1 m2
v m1 −2,m2 δ1 ,
v m1 ,m2 −2 δ2 ,
k 1 m2
u k 1 m2 ·
C kIN+n ,m . 1 1 2 STEP 14: Repeat steps 9, 10 and 11 until n2 = N. STEP 15: Set n1 = n1 + 1, repeat steps 12, 13 and 14 until n1 = N. OUTPUT: |Ψ (t + δt ) = n1 n2 C nOUT |n1 n2 . 1 ,n2 References [1] A.L. Fetter, J.D. Walecka, Quantum Theory—Many-Body Problems, International Series in Pure and Applied Physics, McGraw-Hill Publishing Co., New York, 1971. [2] L. Mandel, E. Wolf, Quantum Coherence and Quantum Optics, Cambridge Univ. Press, 1995. [3] D. Portes Jr., et al., Eur. Phys. J. D 48 (2008) 145–149. [4] H. Rodrigues, et al., Physica A 311 (2002) 188; D. Portes Jr., et al., Physica A 329 (2003) 391. [5] See, e.g., S.-W. Tsai, A.F.R. de Toledo Piza, Phys. Rev. A 53 (1996) 3683; K.M. Ng, C.F. Lo, Phys. Lett. A 230 (1997).