Efficient optimization of state preparation in quantum networks using quantum trajectories
EEfficient optimization of state preparation inquantum networks using quantum trajectories
Michael H Goerz , , Kurt Jacobs , , E-mail: [email protected] Edward L. Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA. U.S. Army Research Laboratory, Computational and Information SciencesDirectorate, Adelphi, MD 20783, USA. Department of Physics, University of Massachusetts at Boston, Boston, MA02125, USA Hearne Institute for Theoretical Physics, Louisiana State University, BatonRouge, LA 70803, USA
Abstract.
The wave-function Monte-Carlo method, also referred to as the useof “quantum-jump trajectories”, allows efficient simulation of open systems byindependently tracking the evolution of many pure-state “trajectories”. Thismethod is ideally suited to simulation by modern, highly parallel computers. Herewe show that Krotov’s method of numerical optimal control, unlike others, canbe modified in a simple way so that it becomes fully parallel in the pure stateswithout losing its effectiveness. This provides a highly efficient method for findingoptimal control protocols for open quantum systems and networks. We apply thismethod to the problem of generating entangled states in a network consisting ofsystems coupled in a unidirectional chain. We show that due to the existence of adark-state subspace in the network, nearly-optimal control protocols can be foundfor this problem by using only a single pure-state trajectory in the optimization,further increasing the efficiency.
Keywords : quantum optics, quantum networks, optimal control, numericalmethods a r X i v : . [ qu a n t - ph ] A ug fficient optimization using quantum trajectories
1. Introduction
The control of simple quantum systems is central to the future realization of quantumtechnologies [1, 2, 3, 4]. Such systems are usually open, meaning that they experiencea significant source of noise from their environments [5, 6]. Systems that form thenodes in a quantum network, in which the nodes are connected together via traveling-wave fields are also necessarily open due to the interaction with the fields [6, 7, 8, 9].At the simplest level, open quantum systems and networks can be described by afirst-order differential equation for the density matrix called a master equation. Auseful feature of master equations is that the evolving density matrix for the systemor systems can alternatively be written as the average over an ensemble of pure statesin which the evolution of each state is given by an equation with a random (stochastic)component [10, 11, 12, 13, 6]. For a d -dimensional system the representation of thedensity matrix scales as d while that of a pure state scales only as d . For this reasonthe pure-state ensemble approach, or “quantum trajectory” approach, can provide asignificant reduction in numerical overhead. Here we use the trajectory method inwhich the stochastic element is a series of “jumps” at random times, often referred toas the “quantum Monte Carlo” or “quantum-jump trajectories” method [14, 15, 16].The ability to vary one or more parameters of the Hamiltonian of a system overtime is a powerful tool for controlling the system. The problem of determining thetime-dependence, or “pulse shapes” required to realize a given evolution is highlynon-trivial, however, and as a result numerical search methods, often referred to as“optimal control” methods, are very useful for this purpose [17, 18]. To perform anumerical search requires simulating the evolution of a system or network many timesunder differing pulse shapes to sufficiently explore the “control space” of possiblepulses. The use of quantum trajectories for such simulations, especially as the size ofa network increases, is thus desirable.Here we examine the form that gradient-based numerical search methods takewhen written in terms of the ensemble of pure state evolutions (the “quantumtrajectories”) that simulate the master equation. We show that unlike other methods,Krotov’s method can be formulated in such a way that updates to the pulse shapesare obtained either from independent trajectories, or by quadratic cross-referencing oftrajectories. We then show via a numerical example that Krotov’s method remainseffective even for a small number of trajectories, including a single trajectory.In the second part of this work we apply the above “Krotov trajectory method” toa problem of preparing a class of entangled states in a simple network. This networkconsists of a “chain” of systems connected via a single unidirectional field. Eachsystem, or “node” in the network is a cavity containing a system with two stablelevels, and for which the interaction with the cavity mode can be turned on and offvia a local control field. One way in which these nodes can be realized is by having athree-level atom in each cavity, in which the three levels are in a lambda configuration.The two ground states provide the two stable levels, and the interaction between thesetwo levels and the cavity mode is mediated by an off-resonant coupling to the excitedstate. Our problem may be thought of as a generalization of the configuration usedby Cirac et al. in their seminal work in which they showed how the state of a two-levelsystem in a cavity could be transferred to a second two-level system in a second cavityvia a unidirectional field [19]. They found that, so long as the coupling of each cavityto the traveling-wave field was the only source of dissipation, the coupling betweenthe two-level systems and their respective cavities could be controlled in such a way fficient optimization using quantum trajectories et al. is to say that there is a subspace of joint states of the two nodes thatare “dark” in the sense that no photons will be emitted from the systems when inthis subspace [20, 21, 22, 23]. The protocol is able to realize the transfer by evolvingexclusively within this subspace. Protocols that employ the transfer method proposedby Cirac et al. have been realized in a number of experiments [24, 25, 26, 27, 28].Here we consider the preparation of dark states of N nodes in which a singleexcitation, originally prepared in the first node, is shared equally between all thenodes. We find that by controlling only the local coupling between the nodes and thefield link it is possible to prepare these multi-partite entangled states while remainingwithin the dark-state subspace, thus realizing essentially perfect preparation (up tothe infidelity imposed by the finite duration of the protocol and any sources of lossand decoherence additional to the field link). We also find that because the optimalprotocols are able to enter the dark-state subspace quickly on the timescale of thetransfer, the majority of the numerical search can be performed with a single pure-state trajectory, thus greatly increasing the efficiency.Our trajectory-based optimization method is applicable to essentially any controltask in which the effect of the environment can be written as an average over pure-state trajectories. This is true for all Markovian master equations that have theLindblad form, and a few non-Markovian master equations [29, 30, 31, 32, 33, 34, 35].Our method is especially useful for situations in which there are dark-state spacesin which the effect of the bath vanishes. This can be true for quantum networks, asin the example we explore here. Another scenario in which this holds is that of thestorage and retrieval of a photon wave-packet from a collective dark state of an atomicensemble [36, 37, 38, 39].The rest of this paper is structured as follows. In Section 2 we first review theMonte-Carlo wavefunction method and the two preeminent gradient-based methods ofquantum optimal control, Gradient-Ascent-Pulse-Engineering (GRAPE) and Krotov’smethod. Then, in Section 2.3, we develop the central result of this paper, a trajectory-based variant of Krotov’s method. We illustrate the method in a numerical case studyin Section 3. We first introduce the network model for a chain of cavities containingtrapped atoms. Then, in Sections 3.2–3.4, we present the result of applying optimalcontrol to the creation of a dark state in that system, and analyze the convergenceand the noise properties of the optimized pulses. Section 4 concludes.
2. Optimal Control via Quantum Trajectories
We consider the problem of initializing a quantum network from a well-definedseparable state | Ψ( t = 0) (cid:105) to a (generally entangled) state | Ψ tgt ( t = T ) (cid:105) . The networkis characterized by a Hamiltonian ˆH = ˆH + (cid:80) i u i ( t ) ˆH i with one or more control fields u i ( t ), and a set of Lindblad operators { ˆL l } to describe dissipative processes. In the fficient optimization using quantum trajectories ˆ ρ ( t ) = E ( t, ˆ ρ (0) denotes the solution of themaster equation ∂ ˆ ρ∂t = − i ¯ h L ˆ ρ = − i ¯ h (cid:104) ˆH ( { u i ( t ) } ) , ˆ ρ (cid:105) + L D ˆ ρ (1)with the boundary condition ˆ ρ = | Ψ(0) (cid:105)(cid:104)
Ψ(0) | and the dissipator L D .For a network consisting of more than a trivially small number of nodes, in general,the exponential growth of the Hilbert space will quickly render the numerical treatmentof the system dynamics in Liouville space unfeasible. Significantly less resources arerequired to simulate the dynamics through quantum trajectories. In this case, ˆ ρ ( t ) isapproximated as an average of the pure states | Ψ k ( t ) (cid:105)(cid:104) Ψ k ( t ) | resulting from M → ∞ trajectories indexed by k , ˆ ρ ( t ) = lim M →∞ M (cid:88) k | Ψ k ( t ) (cid:105)(cid:104) Ψ k ( t ) | . (2)Each trajectory | Ψ k ( t ) (cid:105) is a possible, statistically independent, pure-state realizationof the system dynamics, taking into account dissipative processes. There are severalways to obtain the trajectories, corresponding to different physical measurements onthe environment (e.g. the field into which the system decays), and governed by differentstochastic equations of motion. Measurements of one or both of the quadratures of thefield result in a quantum diffusion equation (QSDE) [40, 13], while photon countingresults in quantum jumps (MCWF) [15, 16]. For the purposes of evaluating the densitymatrix via (2), the various kinds of trajectories are equivalent. We will focus on thequantum jump method because it is the most straightforward to realize numerically.The MCWF algorithm for propagating the trajectory state | Ψ k (cid:105) is [41] define the non-Hermitian effective Hamiltonian ˆH eff = ˆH + (cid:88) l ˆL † l ˆL l (3) draw a random number r ∈ [0 , propagate until | Ψ k ( t ) | = r apply an instantaneous quantum jump | Ψ k ( t ) (cid:105) → ˆL l | Ψ k ( t ) (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆL l | Ψ k ( t ) (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , choosing ˆL l from the set of all Lindblad operators with probabilityp( ˆL l ) = (cid:68) Ψ k ( t ) (cid:12)(cid:12)(cid:12) ˆL † l ˆL l (cid:12)(cid:12)(cid:12) Ψ k ( t ) (cid:69) (4) draw a new random number r ∈ [0 ,
1) and continue the propagation normalize any resulting | Ψ k ( t ) (cid:105) .In performing the propagation with discrete time-steps, as is essential numerically,it is crucial to determine the jump times in step 3 with high precision (throughinterpolation and bisection), and in this way multiple instantaneous jumps can beincluded within one time step. Using this approach, relatively long time steps arepossible. An alternative “first-order” scheme [16] in which at most one jump mayhappen in each time step is simpler to implement, but requires small time steps to beaccurate. fficient optimization using quantum trajectories We seek optimal solutions u i ( t ) to minimize the error functional J T = 1 − (cid:10)(cid:10) ˆ ρ ( T ) (cid:12)(cid:12) ˆP tgt (cid:11)(cid:11) , (5)for a fixed duration T , where ˆP tgt is the projector onto the target state | Ψ tgt (cid:105) . Weuse the notation (cid:104)(cid:104) a | b (cid:105)(cid:105) ≡ tr[ a † b ] for the Hilbert-Schmidt overlap of two operators.There are well-established gradient-based methods of numerical optimal control tosolve this problem. The two preeminent ones are Gradient-Ascent-Pulse-Engineering(GRAPE) [42], usually as a quasi-Newton method in combination with the L-BFGS-Balgorithm [43, 44], and Krotov’s method [45, 46, 47]. Both of these are iterative : Theystart from “guess” fields u (0) i ( t ) and calculate updates ∆ u i ( t ) = u (1) i ( t ) − u (0) i ( t ) thatdecrease the value of the target functional; for the next iteration, u (1) i ( t ) becomes theguess, and the procedure repeats until convergence is reached.The GRAPE method starts with the assumption that u i ( t ) with t ∈ [0 , T ] isdiscretized into n t time steps, such that the time evolution up to a point t j is ˆ ρ ( t j ) = E ( t j , ˆ ρ (0) = (cid:32) (cid:89) j (cid:48) = j E j (cid:48) (cid:33) ˆ ρ (0) (6)where E j ≡ E ( t j , t j − ) denotes the dynamical map ˆ ρ ( t j − ) → ˆ ρ ( t j ), with t = 0and t n t = T . The indices in the product run backwards to account for timeordering. The control fields are approximated as constant within the time interval, u ij ≡ u i ([ t j − , t j ]). Each time-local value u (0) ij is modified according to the gradient,∆ u ij ∝ ∂J T ∂u ij = − (cid:28)(cid:28) ˆP (0) ( t j ) (cid:12)(cid:12)(cid:12)(cid:12) ∂ E j ∂u ij (cid:12)(cid:12)(cid:12)(cid:12) ˆ ρ (0) ( t j − ) (cid:29)(cid:29) , (7)where ˆP (0) ( t j ) = (cid:32) n t (cid:89) j (cid:48) = j +1 E (0) † j (cid:48) (cid:33) ˆP tgt (8)is the target state backward-propagated with the conjugate Lindbladian. Thesuperscript zero indicates the dynamics under the guess controls.Krotov’s method takes a slightly different approach. Starting from time-continuous guess controls u (0) i ( t ), it analytically constructs pulse updates ∆ u i ( t ) thatminimize an extended functional, commonly J = J T + (cid:88) i λ i S i ( t ) (cid:90) T [ u (1) i ( t ) − u (0) i ( t ) (cid:124) (cid:123)(cid:122) (cid:125) =∆ u i ( t ) ] d t , (9)where λ i is an arbitrary weight determining the overall magnitude of ∆ u i ( t ), and S i ( t ) ∈ [0 ,
1] is a shape function that may be used to enforce boundary conditions. Asthe optimization converges, ∆ u i ( t ) →
0, such that J and J T become equivalent. Thecontrol problem is solved by the update equation∆ u i ( t ) = S i ( t ) λ i Im (cid:28)(cid:28) ˆP (0) ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ∂ L ∂u i ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ˆ ρ (1) ( t ) (cid:29)(cid:29) , (10)where ˆP (0) ( t ) is a state backward-propagated with the conjugate Lindbladian, cf. (8),with the boundary condition ˆP ( T ) = − ∂J T ∂ (cid:104)(cid:104) ρ ( T ) | = ˆP tgt , (11) fficient optimization using quantum trajectories u i ( t ) discretized to a time grid (∆ u i ( t ) → ∆ u ij , t → t j ), resulting in a formula thatsuperficially resembles (7), with at least two crucial differences:First, the forward-propagated state ρ ( t ) uses the updated controls u (1) i ( t ) in (10),and the guess controls u (0) i ( t ) in (7). This makes Krotov’s method sequential (updatesat later times depend on updates at earlier times), whereas GRAPE is concurrent (updates are independent).Second, the local derivative term in each equation differs; The gradient derivative ∂ E j /∂u ij can be evaluated similarly to E j itself, under the assumption that E j is asimple exponential [48]. In comparison, ∂ L /∂u ij is more straightforward to evaluate(and instantaneous). In full generality, L in (10) is the right-hand side of the equationof motion for ˆ ρ , whether or not this is a master equation in Lindblad form [47]. Forthe particular form (1) with a time-independent dissipator, ∂ L /∂u ij = [ ˆH i , ρ ( t j )].The backward-propagated state ˆP (0) ( t ) is the same for Krotov’s method andGRAPE only for the specific case of a direct state-to-state transfer, functional (5).More generally, in Krotov’s method, the boundary condition is determined explicitlyby (11), which for other optimization tasks (e.g., gate optimization) does not yieldstates identical to the target states (but usually, a linear combination of them). Eventhough for GRAPE, the boundary condition also implicitly depends on the choice offunctional, for all optimization tasks in quantum control that we are aware of, thebackwards propagation always starts from the target states. While we have assumedlinear controls and a convex functional, Krotov’s method can be also be extended togo beyond these assumptions [47]. There are two possibilities for rewriting the optimization in terms of quantumtrajectories. First, we can insert the expansion of ˆ ρ (2) into the functional (5) tofind J T = 1 − lim M →∞ M M (cid:88) k =1 (cid:12)(cid:12) (cid:104) Ψ k ( T ) | Ψ tgt (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) ≡ τ k (cid:12)(cid:12) . (12)That is, the fidelity is simply the average of the fidelities from each trajectory. Weseek a single control that simultaneously implements the state-to-state transition forevery | Ψ k (cid:105) .We now wish to ask, given the form of J T , how the update rules for GRAPEand Krotov’s method might be written in terms of the individual trajectories. ForGRAPE this is not so simple; the form of the update rule is based on the fact that theevolution of ˆ ρ in each time step is given by the application of a linear operator E j , andthis is no-longer true for the trajectories. Second, in a given time step the evolutionof a given trajectory may contain one or more jumps at random times, which furthercomplicates the calculation of the derivative with respect to the controls u ij .In contrast, for Krotov’s method the functional (12) leads to the updateequation [46] ∆ u i ( t ) = S i ( t ) M λ i M (cid:88) k =1 Im (cid:68) χ (0) k ( t ) (cid:12)(cid:12)(cid:12) ˆH i (cid:12)(cid:12)(cid:12) Ψ (1) k ( t ) (cid:69)(cid:124) (cid:123)(cid:122) (cid:125) ≡ ∆ u ik ( t ) , (13) fficient optimization using quantum trajectories χ (0) k ( t ) is backward-propagated with the boundary condition χ (0) k ( T ) = − ∂J T ∂ (cid:104) Ψ k | = τ (0) k | Ψ tgt (cid:105) , (14)where τ (0) k is τ k from (12) evaluated for the guess pulse. The crucial difference is thatfor GRAPE, the gradient is evaluated for the time evolution operator, whereas forKrotov’s method, the gradient is with respect to the equation of motion. Only thelatter is immediately compatible with the decomposition into quantum trajectories.From a numerical perspective, the update equation (13) can be evaluated withhigh efficiency: the contributions to the update ∆ u ik ( t ) from every trajectory arecompletely independent. If the numerical method is parallelized, so that differenttrajectories are simulated on different processors, and a Message-Passing-Interface(MPI) is used for interprocessor communication, only these scalar values need to becommunicated. We thus reap the full benefit of the MCWF method: reduction of theoverall dimensionality by a factor of d , the dimension of the total Hilbert space. Therequired memory per compute node is reduced directly by that factor; the worst-caseCPU time reduces by d (as the fundamental operation in the time propagation is amatrix-vector multiplication).We can obtain an alternative formulation of the control problem in terms ofquantum trajectories by leaving the functional (5) in terms of the density matrix, andinstead expanding the update equation in terms of the trajectories. For GRAPE thisapproach has the same issues as before, whereas for Krotov’s method this results inthe update equation∆ u i ( t ) = S i ( t ) M λ i M (cid:88) k,k (cid:48) =1 Im (cid:104) (cid:68) ξ (0) k ( t ) (cid:12)(cid:12)(cid:12) ˆH i (cid:12)(cid:12)(cid:12) Ψ (1) k (cid:48) ( t ) (cid:69) (cid:68) Ψ (1) k (cid:48) ( t ) (cid:12)(cid:12)(cid:12) ξ (0) k ( t ) (cid:69) (cid:105) , (15)using the decomposition ˆP (0) ( t ) = lim M →∞ M M (cid:88) k =1 | ξ k ( t ) (cid:105)(cid:104) ξ k ( t ) | (16)for the backward-propagation with | ξ (0) k ( T ) (cid:105) = | Ψ tgt (cid:105) . The above update nowincorporates cross-trajectory information, including overlaps of the propagated statesfrom each trajectory with those from every other trajectory. Compared to the fulldensity matrix optimization, this alternative method maintains the saving in CPUtime and memory that comes from evolving the trajectories independently, but doesincur significantly more communication overhead.Both trajectory formulations of Krotov’s method, the “independent trajectory”and “cross-trajectory” formulations, guarantee monotonic convergence for M → ∞ .However, the numerical efficiency rests on the assumption that M can be keptreasonably small, ideally smaller than d (and certainly smaller than d ). For asignificantly larger number of trajectories, the total numerical resources will approachor even surpass those required for a full density matrix propagation. Note that evenif there is no benefit in the total numerical resources, the memory per compute node is still reduced by a factor of d .For finite (small) M , there is no guarantee of monotonic convergence. Thisis especially true because the stochastic jumps are different between iterations, aswell as between the forward and the backward-propagation within the same iteration.Moreover, we expect the instantaneous jumps in the trajectory to cause jumps also in fficient optimization using quantum trajectories Ω ( t ) g Node 2 Node N κ . . .. . .. . . | e i | g i| r i Ω ( t ) g ˆ σ rg ˆa ∆ Figure 1.
The structure of our network. A number of cavities, each containinga three-level atom in the Lambda configuration, are connected together by atraveling optical field that enters each cavity in turn via one of its end mirrors.Each cavity is coupled to the field via a circulator (not shown), so that the cavitiesonly couple to field modes that are traveling from left to right. The field that isoutput from the first cavity thus enters the second, but not vice versa, in what isreferred to as a “cascade” connection. the control update. These jumps will average out both with the number of trajectoriesand the number of control iterations, but they may slow down convergence and mayintroduce noise into the optimized control field. To analyze how strongly the resultingoptimized control fields are affected by these two factors, and to what extent thecross-trajectory update equation might be more robust, we now consider a practicalexample.
3. Creation of an entangled dark state in a quantum network
We consider a network in which a number of “nodes” are connected via a traveling-wave field, depicted in Fig. 1. The nodes are connected to the field via circulatorsso that the field propagates from one node to the next in a single direction only, aconfiguration often referred to as a “cascade” network. Each node consists of a cavitycontaining a single atom. The atom encodes a qubit in two degenerate states | g (cid:105) , | e (cid:105) ,driven via a Raman transition through an auxiliary level | r (cid:105) : | g (cid:105) ↔ | r (cid:105) couples througha time-dependent drive Ω( t ) with detuning ∆, while | r (cid:105) ↔ | e (cid:105) is statically coupledwith coupling strength g , see Fig. 1. For large ∆, the level | r (cid:105) can be adiabaticallyeliminated [19], resulting in an effective two-level system whose coupling to the cavityis controlled via the drive Ω( t ).For a single node labeled ( i ), the effective Hamiltonian in the rotating waveapproximation then consists of the drift term ˆH and a driven Jaynes-Cummings term ˆH d . With the cavity detuned from the driving field by g / ∆ to compensate for theStark shift, we find ˆH ( i ) = ˆH ( i )0 + ˆH ( i ) d , (17) ˆH ( i )0 = − g ∆ ˆa † i ˆa i + g ∆ ˆΠ ( i ) g ⊗ ˆa † i ˆa i , (18) ˆH ( i ) d = − iΩ i ( t ) g (cid:16) ˆ σ ( i ) e , g ⊗ ˆa i − c.c. (cid:17) , (19) fficient optimization using quantum trajectories ˆa i is the cavity lowering operator, ˆΠ g is the projector on the qubit ground state,and ˆ σ e , g is the raising operator of the qubit. Leakage of photons out of the cavity isdescribed by the Lindblad operator ˆL ( i ) = √ κ ˆa i . (20)There is also spontaneous decay from level | r (cid:105) of the atom, but the decay is suppressedfor large ∆ (the prerequisite for the adiabatic elimination that yields the effective two-level system), and can be neglected [19].We now couple N cavities to a traveling-wave field via circulators, so that thefield propagates from cavity to cavity in only one direction, as depicted in Fig. 1.The resulting master equation describing the evolution of all the cavities is obtainedusing input-output theory [8, 49, 9]. The explicit application of the theory is throughthe “SLH” formalism of Gough and James [7, 50, 51]. While the calculation can bedone by hand, it is tedious, and we instead use the QNET software package, whichautomates the SLH formalism ‡ .The total Hamiltonian for the entire network is ˆH = N (cid:88) i =1 ˆH ( i ) + N (cid:88) i,j>i ˆH ( i,j ) , (21) ˆH ( i,j ) = i κ ˆa † i ˆa j + c.c. (22)Note that the connection via the field yields the additional static interactionHamiltonian ˆH ( i,j ) between all nodes of the network (not merely between nearestneighbors). The Lindblad operator for the entire network is ˆL = N (cid:88) i =1 ˆL ( i ) . (23)Thus, both the interaction between nodes and the total dissipation of the networkare proportional to the decay rate κ . In order to create an entangled state quickly, κ (the interaction) must be large, which would seem to imply rapid decay of thenetwork. However, the form of the Lindblad operator allows for the individualsummands ˆL i to destructively interfere, canceling overall dissipation. If the dark-statecondition [20, 21, 22, 23] ∀ t : (cid:68) Ψ( t ) (cid:12)(cid:12)(cid:12) ˆL † ˆL (cid:12)(cid:12)(cid:12) Ψ( t ) (cid:69) = 0 (24)is fulfilled, an entangled state may be generated without loss of coherence.Nonetheless, we must consider the dynamics of a quantum network to be inherentlyhighly dissipative, as the dark-state condition is only fulfilled for carefully chosencontrol pulses that evolve the network exclusively within the dark-state manifold.In applying our numerical search method to the problem of generating anentangled dark state in the network we assume that the field link between cavities isthe only significant source of loss. Since the entangled dark state exists only becausethe same field link is coupled to multiple systems, any sources of loss that are notsimilarly collective will necessarily limit the fidelity with which the entangled statecan be prepared. Such losses include spontaneous emission from the atomic excitedstate (which is suppressed by the detuning), internal cavity losses, and loss en route between the cavities (as opposed to loss out of the open end of the field link, which is ‡ The python package QNET allows the user to construct quantum input-output networks, and willcalculate the master equation for any such network [52]. fficient optimization using quantum trajectories
We assume that at time zero, the network is initialized in the state | Ψ( t = 0) (cid:105) = | eg . . . g (cid:105) . That is, the atom in the cavity of the leftmost node is in the excited state,while all other nodes (and all cavities) are in the ground state. The state of thecavities, | . . . (cid:105) is implicit in our notation. We now seek the control fields Ω i ( t ) (onefield for each node) that will bring the system into an entangled dark state at a fixedfinal time T , | Ψ( t = T ) (cid:105) = √ N ( | eg . . . g (cid:105) + | ge . . . g (cid:105) + · · · + | gg . . . e (cid:105) ). The cavity forall nodes again must be in the ground state at time T .The optimization problem is solved by minimizing (5) using Krotov’s method.We first do this optimization using the full density matrix, update (10), and then –for comparison – using a trajectory optimization, both for independent trajectories,functional (12) with update (13), and for the cross-trajectory method, update (15).One reason for choosing this particular physical system as an example to illustrateour proposed optimization method is that the Hamiltonian (21) preserves the totalexcitation in the system. Thus, it is sufficient to model both the cavity and the atom ateach node as a two-level-system. Moreover, when representing the state numerically,we need only include the single-excitation subspace along with the ground state toaccount for the possibility of spontaneous decay (the loss of the excitation out theopen end of the field that links the cavities). The effective dimension of the resultingHilbert space is 2 N + 1, instead of the exponentially scaled 4 N . This makes a fulldensity matrix optimization feasible even for relatively large N , allowing us to compareagainst the trajectory optimizations. We consider first the simplest case of N = 2nodes, with control fields Ω ( t ) and Ω ( t ). In this case, the dark state is the Bell state | Ψ (cid:105) tgt = 1 √ | eg (cid:105) + | ge (cid:105) ) . (25)This allows us to easily illustrate the characteristics of the optimized pulses and nodedynamics. To ensure that our conclusions hold for larger networks, we then consider anetwork of 20 nodes. In principle, even larger networks would be numerically feasible.However (as we will see), the optimization becomes increasingly harder, due to thedrastically larger search space resulting from N independent control pulses. When thenumber of nodes is significantly higher that 20, we do not find good control solutionswithin a reasonable number of iterations (even for the optimization that employs thefull density matrix). We expect that it might be possible to sufficiently reduce thesearch space for larger networks by not treating all pulses as fully independent.For full generality, we use a dimensionless Hamiltonian where energies are writtenin units of g . Consequently, time can be expressed in units of ¯ h/g . For ∆ = 100 g and fficient optimization using quantum trajectories ~ /g )0 . . . . pu l s e a m p li t ud e ˜Ω ind1 ˜Ω cr1 ˜Ω ρ ˜Ω ind2 ˜Ω cr2 ˜Ω ρ guess —0.2— - . - (a) 0 1 2 3 4 5time ( ~ /g )10 − − h ˆL † ˆL i . . . e x c i t a t i o n qubit 1qubit 2cavity 1cavity 2(b) Figure 2.
Optimized pulses and dynamics. (a) Optimized pulses for the creationof the dark state in a two-node network. Ω ρ , are the pulses for node 1, 2 resultingfrom optimizing using the density matrix. Ωind , and Ωcr , are the pulses fromoptimizing using two independent trajectories, respectively two cross-referencedtrajectories. All pulse amplitudes are normalized as ˜Ω ρ, tr , ≡ Ω ρ, tr , / ρi ( t ), as theexpectation values of the qubit excitation, (cid:104) ˆΠ ( i ) e (cid:105) , and the cavity excitation, (cid:104) ˆa † i ˆa i (cid:105) ,with i = 1 ,
2. In the bottom, the value of the dark-state condition (24). κ = g , through systematical variation of T , we find that the two-node dark state cancomfortably be realized for T = 5 ¯ h/g . For 20 nodes, we correspondingly increase theprocess duration to T = 50 ¯ h/g . In general, for a fixed number of nodes, the minimumtime required to entangle the nodes of the network is proportional to the interactionHamiltonian (22), that is, proportional to the value of κ .Fig. 2 shows the result of an optimization for two nodes starting from a simpleBlackman-shaped guess pulse, see the dotted gray line in panel (a). For a directdensity-matrix optimization using Krotov’s method, that is the minimization of (5)with update (10), after 5000 iterations an error of J T = 1 . × − is achieved. Theoptimized pulses Ω ( t ) and Ω ( t ) for the first and second node are shown in Fig. 2 (a).We show the pulses resulting from the full density matrix optimization (superscript“ ρ ”; black, red), from two independent trajectories (superscript “ind”; dark blue, darkorange), and from two cross-referenced trajectories (superscript “cr”; light blue, lightorange). The corresponding dynamics for the pulses resulting from the density matrixoptimization are in panel (b), showing the smooth transition from the initial state | eg (cid:105) to the superposition of | eg (cid:105) and | ge (cid:105) : the excitation of the first and second qubit, theplotted expectation values (cid:104) Π (1) e (cid:105) and (cid:104) Π (2) e (cid:105) , correspond directly to the population of | eg (cid:105) and | ge (cid:105) , as the Hamiltonian preserves excitations ( (cid:104) ee | ee (cid:105) ≡ fficient optimization using quantum trajectories (cid:104) ˆL † ˆL (cid:105) in the bottom of Fig. 2 (b), which is close to zero). The two cavities are phase-shiftedby π , that is (cid:104) ˆa (cid:105) ≈ −(cid:104) ˆa (cid:105) , resulting in destructive interference for the total dissipator.As we have not taken into account any other dissipation channels (such asspontaneous decay of the qubit), the deviation from the dark-state condition (24), isthe only factor that limits the fidelity of the transfer. By continuing the optimizationbeyond 5000 iterations, the error could be made arbitrarily small in principle.Optimizing with independent trajectories, yields very similar results, correspond-ing to the optimized pulses shown in Fig. 2 (a) (dark orange, dark blue) being closebut not identical to the pulses obtained from direct density matrix optimization (red,black). For the cross-trajectory method, the dynamics and optimized pulses (lightcolored) are virtually indistinguishable from the density matrix optimization.In general, running both optimization methods with differing numbers oftrajectories, and the density matrix method, all yield comparable errors. Neverthelessthere are differences in the rate of convergence, as well as in certain details of thepulse features. As shown in the zoomed inset of Fig. 2 (a), there is some noise in theoptimized pulses obtained from trajectory optimization, compared to the perfectlysmooth result of the density matrix optimization. This noise originates from thediscontinuous jumps in the trajectories. The optimized pulses and dynamics for the 20-node network (not shown) have the same qualitative characteristics. In the following,we will analyze how both the convergence rate and the noise in the optimized pulsesdepend on the choice of optimization method and the number of trajectories. Fig. 3 shows the convergence, as measured by the error (5), over the number ofiterations in the control procedure. In each iteration, the pulse is updated accordingto the appropriate update equation. Panel (a) shows the results for the two-nodenetwork, while panel (b) shows the result for the 20-node network.The convergence for the full density matrix optimization, update (10), is shownas the black dashed line. It compares to the convergence for using 1-128 independenttrajectories (“indep. trajs”), update (13), and 2-128 cross-trajectories (“cross-trajs”),update (15).We first observe that the optimization for the 20-node network is much harderthan the optimization for 2 nodes, requiring several thousand iterations to reach errorsbelow 10 − . In particular, the density matrix optimization has a long initial plateaubefore starting to converge. This behavior is exacerbated for more nodes, which iswhy we have limited ourself to considering at most a 20 node network. Remarkably,the use of a trajectory optimization performs significantly better in this initial phase,largely avoiding the plateau. We presume this is because the dissipative jumps affectthe trajectories much more strongly than the smooth decay affects the density matrix,leading to a steeper gradient. Both the two-node network and the 20-node networkhave the same qualitative convergence characteristics; the larger number of nodessimply “stretches” them out.The convergence of the cross-trajectory optimization consistently outperformsthe optimization using independent trajectories. Asymptotically, the independent fficient optimization using quantum trajectories − − − e rr o r ρ indep. trajscross-trajs(a) 2 nodes, T = 5 ~ /g − J T ρ indep. trajscross-trajs (b) 20 nodes, T = 50 ~ /g number of trajectories124, 8, 16, 32, 64, 128 Figure 3.
Comparison of convergence of optimization variants, for the two-nodenetwork, panel (a), and the 20-node network, panel (b). For using between 1and 128 independent trajectories, the error 1 − (cid:104)(cid:104) ˆ ρ ( T ) | ˆP tgt (cid:105)(cid:105) after each iterationis shown. This compares to the optimization using between 2 and 128 cross-trajectories, and the optimization directly in Liouville space (dashed black curve).As the number of trajectories has almost no effect on the convergence and thecurves are not distinguishable, all the curves for 4–128 trajectories are plotted withthe same red dashed line style. All errors are evaluated from the exact densitymatrix evolution under the optimized pulses, which differs from the functional J T actually guiding the optimization. For the example of an optimization using twoindependent trajectories in the two-node network, the value of J T is shown in theinset of panel (a). trajectories eventually lag behind the full density matrix optimization, but onlyslightly so. The cross-trajectory optimization asymptotically reaches the same erroras the full density matrix optimization (clearly for the two-node network, but we mayextrapolate this claim to hold also for 20 nodes).A second remarkable observation is that the number of trajectories has almostno effect on the convergence. Thus, even a single trajectory, or two cross-referencedtrajectories, can yield a satisfactory result, a significant reduction of numerical effortcompared to the use of full density matrices. It is worth pointing out that as longas the number of trajectories does not exceed the number of CPU cores availableon a compute node, it would be easy to parallelize the method not using a message-passing interface, but instead using a shared-memory approach (e.g. OpenMP). In thiscase, the additional numerical overhead of the cross-trajectory optimization relativeto independent trajectories largely disappears.For the trajectory optimizations, the error (10) that is shown in Fig. 3 is notidentical to the functional J T that the optimization directly aims to minimize (12).Because of the finite number of trajectories and the randomness inherent in thequantum jumps, J T itself does not converge monotonically. An example for twoindependent trajectories is shown in the inset of Fig. 3. The loss of a photon fromthe cavity out of the open end of the field link (that is, a quantum jump) in agiven iteration results in a total loss (error 1.0), as visible in the spikes in J T . Inaddition, the lower envelope ( J T in the absence of photon loss) shows a local minimumaround 250 iterations, owing to the two competing requirements of the optimization:achieving the desired target state (in the absence of photon loss), and minimizing decay(the dark-state condition). For cross-trajectory optimization, J T behaves similarly fficient optimization using quantum trajectories number of trajectories M − − − n o i s e m e a s u r e ν (indep.) ν (indep.) ν (cross) ν (cross) Figure 4.
Scaling of noise artifacts in optimized pulses with the numberof trajectories M . For both independent trajectories and cross-trajectoryoptimization, the noise measure according to (26) is shown for the pulse in eachnode. The dashed lines represent a fit to ν i ∝ / √ M (not shown). The non-monotonic convergence of J T should not be held against thetrajectory optimization method: as shown in Fig. 3 the actual physical error doesdecrease monotonically. The use of trajectory optimization results in some noise in the control fields, cf. theinset in Fig. 2 (a). To analyze how the magnitude of this noise scales with the numberof trajectories, we quantify the noise by comparing the optimized control field Ω i ( t )to a smoothed pulse Ω smooth i ( t ), ν i = (cid:90) T (cid:12)(cid:12) Ω i ( t ) − Ω smooth i ( t ) (cid:12)(cid:12) d t . (26)For the purpose of this analysis, a good noise filter should generate an Ω smooth i ( t )that deviates from Ω i ( t ) as little as possible, while achieving a well-defined degreeof smoothness. A standard Savitzky-Golay filter [53] is well-suited to solve thisproblem [54].We find that the noise as defined in (26) asymptotically scales with the number oftrajectories M as ν i ∝ / √ M . This is illustrated in Fig. 4, which compares the noisemeasure for both pulses of the two-node network, and for independent- and cross-trajectory optimization, as a function of the number of trajectories M . The resultsshown here use a five-point cubic convolute in the Savitzky-Golay filter; this choice ofwindow length and polynomial order only affects the proportionality factor, i.e. they-intercept in Fig. 4.For independent trajectories, the 1 / √ M scaling holds for all values of M . Thepulse Ω is noisier than the pulse Ω , which is due to the fact that Ω is closer tothe original guess pulse than Ω , see Fig. 2 (a) – any noise induced by a jump isproportional to the magnitude of the pulse update in a given iteration.The noise for the cross-trajectory optimization is systematically slightly higherthan that for independent trajectories. We may conjecture that this is because theupdate (15) is quadratic in the (discontinuous) states, while the update (13) is only fficient optimization using quantum trajectories < ≈ / √ M scaling, closer to the independent trajectories.In any case, for this specific example, the noise has little effect on the fidelity,as we can see from Fig. 3: the optimization result is robust with respect to the levelof noise introduced by the numerical procedure. If this were not that case, and ifsmooth pulses are required, one may also incorporate the smoothing directly into theoptimization scheme by applying a smoothing filter, or fitting to a parametrized curve,after each iteration (or after some fixed number of iterations). The pulses resultingfrom the optimization in the 20-node network (not shown) have noise characteristicscompatible with the above discussion.
4. Conclusion
We have formulated two variants of Krotov’s method of optimal control based onMonte-Carlo quantum jump trajectories for implementing state-to-state transitions.The fact that it is possible to apply Krotov’s method directly to stochastic trajectoriesis due to the fact that, unlike other methods, Krotov’s method does not require thecalculation of the derivative of the propagator (which for stochastic trajectories doesnot exist). The first variant considers M independent trajectories, and optimizesthe target in each trajectory, averaging pulse updates from the trajectories. It iseasily parallelized to M processes communicating with a message passing interface(MPI). The second “cross-trajectory” variant obtains an update by cross-referencing M trajectories into a coarse approximation of the density matrix evolution. Whenparallelized using MPI, this produces considerable overhead. However, for small M ,a shared-memory parallelization model is efficient.In any case, the use of trajectories provides a considerable numerical advantagefor systems that are strongly dissipative, but for which simulating the full densitymatrix is not feasible due to the large dimension of the Hilbert space. In manycases, the memory to store a density matrix of size d is not available, while a purestate of dimension d can easily be stored and propagated. This is typically the casefor quantum networks, and we have illustrated the method for a simple example ofcreating an entangled dark state for a network of 20 cascaded cavities in which eachcontains a trapped atom.We find that the use of trajectories in the optimization yields results comparableto the full density matrix optimization. In the case of cross-trajectory optimization,the convergence of the optimization even outperforms the density matrix optimization.Remarkably, the number of trajectories has almost no influence in the convergenceor the error that is achieved. The number of trajectories does have an influence in theamount of noise artifacts in the optimized pulses that result from the discontinuous(jump) dynamics. These scale as 1 / √ M in the number of trajectories. At the sametime, they have an only negligible effect on the dynamics and the resulting error, andcould also be removed by applying a smoothing filter to the optimized pulses.Thus, we can conclude that a relatively small number of trajectories (even justa single trajectory) is sufficient to optimize a state-to-state transition in a typicaldissipative system. Further, this trajectory approach is especially efficient whenthe optimal protocol evolves the system within a dark (dissipationless) subspace,a situation which arises, for example, both in quantum networks and atomicensembles [36]. The method naturally extends to optimization problems beyond asimple state preparation, e.g. the realization of quantum gates [55], or any other fficient optimization using quantum trajectories Acknowledgments
The network model used in this paper was derived using the Python QNETpackage [52]. The propagation and optimization were done using the Fortran QDYNlibrary [59]. The pulse smoothing was realized through SciPy’s signal processingtoolbox [60]. This research was supported in part by ASD(R&E) under their QuantumScience and Engineering Program (QSEP), and by the Army High PerformanceComputing Research Center (AHPCRC) (sponsored by the U.S. Army ResearchLaboratory under contract No. W911NF-07-2-0027).
References [1] Blume-Kohout R, Caves C M and Deutsch I H 2002
Found. Phys. Nature Photon. Phys. Rev. A Topological Matter:Topological Insulators, Skyrmions and Majoranas, Lecture notes of the 48th IFF SpringSchool 2017 ed Bluegel S, Mokrusov Y, Schaepers T and Ando Y p 278[5] Breuer H P and Petruccione F 2007
The Theory of Open Quantum Systems (Oxford: OxfordUniversity Press)[6] Jacobs K 2014
Quantum measurement theory and its applications (Cambridge: CambridgeUniversity Press)[7] Combes J, Kerckhoff J and Sarovar M 2017
Adv. Phys. X Phys. Rev. A Quantum Noise (Springer)[10] Carmichael H 1993
An Open Systems Approach to Quantum Optics (Springer-Verlag, Berlin)[11] Diosi L 1986
Phys. Lett. A
Quantum Semiclass. Opt. Phys. Rev. A Classical and Quantum Systems, Proceedings ofthe Second International Wigner Symposium ed Doebner H D, Scherer W and Schroeck F(Singapore: World Scientific) pp 104–115[15] Dum R, Zoller P and Ritsch H 1992
Phys. Rev. A J. Opt. Soc. Am. B Eur. Phys. J. D J. Phys.: Condens. Matter Phys. Rev. Lett. Phys. Rev. A Phys. Rev. Lett. Phys. Rev. Lett.
New J. Phys. Nature fficient optimization using quantum trajectories [25] Yin Y, Chen Y, Sank D, O’Malley P J J, White T C, Barends R, Kelly J, Lucero E, MariantoniM, Megrant A, Neill C, Vainsencher A, Wenner J, Korotkov A N, Cleland A N and MartinisJ M 2013 Phys. Rev. Lett. (10) 107001[26] Wenner J, Yin Y, Chen Y, Barends R, Chiaro B, Jeffrey E, Kelly J, Megrant A, Mutus J, NeillC, OMalley P, Roushan P, Sank D, Vainsencher A, White T, Korotkov A N, Cleland A andMartinis J M 2014
Physical Review Letters
Nature Phys.
URL https://doi.org/10.1038/s41567-018-0115-y [28] Kurpiers P, Magnard P, Walter T, Royer B, Pechal M, Heinsoo J, Salath´e Y, Akin A, Storz S,Besse J C, Gasparinetti S, Blais A and Wallraff A 2018
Nature
Phys. Rev. A Phys. Rev. Lett. Phys. Rev. Lett. Phys. Rev. A Phys. Rev. A Phys. Rev. Lett.
Phys. Rev. Lett.
Phys. Rev. Lett. (22) 5094[37] Liu C, Dutton Z, Behroozi C H and Hau L V 2001 Nature
Phys. Rev. Lett. Phys. Rev. Lett. Quantum State Diffusion (Cambridge University Press)[41] Plenio M B and Knight P L 1998
Rev. Mod. Phys. J. Magnet. Res.
SIAM J. Sci. Comput. ACM Trans. Math. Softw. Global Methods in Optimal Control (New York, NY, USA: Dekker)[46] Palao J P and Kosloff R 2003
Phys. Rev. A J. Chem. Phys.
J. Magnet. Res.
Phys. Rev. Lett. Commun. Math. Phys.
IEEE Transactions on Automatic Control https://github.com/mabuchilab/QNET [53] Savitzky A and Golay M J E 1964 Anal. Chem. Numerical Recipes inFORTRAN; The Art of Scientific Computing (New York, NY: Cambridge University Press)[55] Goerz M H, Reich D M and Koch C P 2014
New J. Phys. J. Chem. Phys.