Optimal control of interacting particles: a multi-configuration time-dependent Hartree–Fock approach

We combine optimal control theory with the multi-configuration time-dependent Hartree–Fock method to control the dynamics of interacting particles. We use the resulting scheme to optimize state-to-state transitions in a one-dimensional (1D) model of helium and to entangle the external degrees-of-freedom of two rubidium atoms in a 1D optical lattice. Comparisons with optimization results based on the exact solution of the Schrödinger equation show that the scheme can be used to optimize even involved processes in systems consisting of interacting particles in a reliable and efficient way.


Introduction
Controlling the dynamics of quantum systems lies at the heart of many of the most exciting scientific developments today. Quantum computation and information processing [1], the design of molecular switches for electronic circuits [2] and the use of attosecond pulses to study the dynamics of electrons on their natural timescale [3] are just a few examples where quantum control is a prerequisite for success. Fundamental tests of quantum mechanics that require precise control of state preparation, e.g. [4], and the control of chemical reactions [5]- [7] are further examples of the importance of quantum control. All of these examples have one central question in common: how to achieve a certain objective given some external parameters that can be controlled, e.g. a laser's frequency, amplitude and polarization?
In order to answer this question, several techniques have been developed over the years. Examples are stimulated Raman adiabatic passage schemes [8], learning algorithms [9], coherent control via interfering pathways [6] and schemes based on optimal control theory (OCT) [10,11]. Because of the generality of the formulation of OCT, it is arguably the most universal approach and, neglecting technical considerations, the resulting control is by construction the optimal one. In addition, the huge success of OCT in 'classical' engineering [12] makes OCT a promising tool for advanced quantum engineering.
Despite this promise, applying OCT to quantum systems entails the practical issue that it requires solving the time-dependent Schrödinger equation. Although this can be done in a straightforward way for one particle, solving the time-dependent Schrödinger equation for two interacting particles in three dimensions (3D) is already a severe challenge; especially in situations involving ionization processes as, e.g. in strong laser fields. For more particles, the exact solution is in general unfeasible because of the large size of the underlying Hilbert space. Thus, one has to resort to approximations. In the context of OCT, these approximations should be as unbiased as possible because it is not known a priori what the optimal control and what the optimal dynamics of the system look like. For fermions, an approximation that satisfies this desired feature perfectly is the multi-configuration time-dependent Hartree-Fock (MCTDHF) method [13]- [15]. This ab initio approach has been used successfully to describe 3 the correlated dynamics of fermions [16]- [19] and thus, the approach is a natural candidate for the optimization of processes involving several interacting fermions. If the particles are distinguishable one can use the multi-configuration time-dependent Hartree (MCTDH) method [20,21] and for bosons the multi-configuration time-dependent Hartree for bosons (MCTDHB) is a natural candidate [22].
The purpose of the present paper is to investigate the combination of the MCTDHF method with OCT. In section 2, we present the MCTDHF method and discuss two different approaches how it can be combined with OCT. In section 3, we investigate one of the two approaches in more detail by comparing the optimization results to the results obtained from the exact solution of the Schrödinger equation. The processes that we optimize in this section are a state-to-state transition in a model of helium and the entanglement of two rubidium atoms in a 1D optical lattice. Finally, section 4 gives a summary and conclusion.

The MCTDHF method
In general, the dynamics of any N -particle system is governed by the time-dependent Schrödinger equation where x k labels both the spin and position of the kth particle (we are using atomic units if not stated otherwise). The Hamiltonian H is given in the following by the usual form In this equation V pp is the particle-particle interaction and the one-particle Hamiltonian h( p k , x k , t) is given by Here, α(t) = α 1 (t), . . . , α M (t) are the time-dependent parameters that can be used to control the dynamics. The basic idea of the MCTDHF method is to approximate the state ψ(x 1 , . . . , x N , t) by the ansatz where both the coefficients c k 1 ... k N (t) and the single-particle spin orbitals ϕ k (t) are allowed to be time dependent. Due to the Pauli principle, the coefficients c k 1 ... k N (t) are antisymmetric under particle exchange. As explained in more detail below, the number of spin orbitals N S controls the accuracy. The crucial idea now is that the time dependence of the coefficients and the spin orbitals is not determined by the Schrödinger equation, but by the Dirac-Frenkel variational principle which allows the spin orbitals to adapt to the solution [21]. As explained in detail 4 in [21], the resulting equations-of-motion for the coefficients and spin orbitals can be cast in several different, but essentially equivalent, forms. The set we are using is where K = {k 1 , . . . , k N }, | K = N l=1 |ϕ k l (t) , and P = N S k=1 |ϕ k ϕ k | have been used. ρ −1 is the inverse density matrix and V pp is the 'mean-field' matrix with elements depending on x and t. The detailed definition of ρ −1 and V pp can be found in, e.g. [15]. Important for the following is that ρ −1 depends on the coefficients and V pp depends on the coefficients and spin orbitals, i.e.
. Several properties of (5) and (6) are worth noting. First, the projector (1 − P) guarantees that the spin orbitals, which are initially orthogonal, remain orthogonal. Second, (6) is, in contrast to the Schrödinger equation, a nonlinear equation because of the dependence of V pp on the spin orbitals. This is the price one has to pay for the self-adaption of the spin orbitals. Thirdly, for N S → ∞ the spin orbitals become a complete basis and P becomes the identity. Thus, the nonlinearity vanishes and the exact time-dependent Schrödinger equation is recovered. Since in practice a finite, time-independent basis of size N g is required to represent the spatial part of the orbitals ϕ k (x, t), the 'exact' Schrödinger equation (in the finite, time-independent basis) is recovered in the case N o = N g , where N o is the number of spatial orbitals (N o = N S /2). Finally, for N S = N , one obtains the time-dependent Hartree-Fock approximation.
In our implementation the spatial orbitals ϕ k are represented on a 1D real-space grid with N g grid points and equidistant spacing dx. We evaluate the derivatives via Fast-Fourier transformation and use a fifth-order Runge-Kutta scheme [23] to integrate the coupled equations (5) and (6). The 'exact' Schrödinger equation is solved on a 2D real-space grid with N g × N g grid points using the split-operator technique [24]. The grid spacing dx is the same as in the MCT-DHF case.

MCTDHF and OCT for state-to-state transitions
Given a set of control parameters α(t), our aim is to steer a quantum system from a fixed initial state |ψ I to a target state |ψ T at fixed final time T . Mathematically, this can be achieved by maximizing the functional where the second term can be used to implement constraints on the values of the controls at different times via the penalties λ 1 (t), . . . , λ M (t) [25]. For instance, a possible constraint can be a smooth switch-on and switch-off of the controls. In the case that the state |ψ(t) obeys the Schrödinger equation, the control equations derived from OCT are [10,11] id t |ψ(t) = H |ψ(t) , |ψ(0) = |ψ I , where |χ(t) is the 'conjugate' state to |ψ(t) . In order to maximize J , these equations must be solved self-consistently. For this purpose, several monotonically convergent algorithms have been developed over the years [26]- [28]. The one we are using for solving (8)-(10) is the Krotov iteration method [26,29].
Since the MCTDHF state |ψ MF (t) satisfies the nonlinear equations (5) and (6), the control equations (8)- (10) in general cannot be used in the context of the MCTDHF approach. Instead it is necessary to derive control equations based on (4)- (6). The resulting equations, however, are very much involved [30]. This can be best understood from calculating the time-derivative of (4). One obtains where the time derivatives are given by (5) and (6). If the spin orbitals do not constitute a complete basis, the right-hand side of (11) cannot be expressed explicitly and solely in terms of the state |ψ MF (t) . Thus, the functional derivative of the rhs of (11) with respect to |ψ MF (t) , which is required for the equation-of-motion of the 'conjugate' state |χ(t) [31], cannot be evaluated in a straightforward way. This problem can be avoided by introducing 'conjugate' variables for the coefficients c K and the spin orbitals ϕ k . For these 'conjugate' variables, the required derivatives can be evaluated directly from (5) and (6). However, as mentioned above, the resulting control equations are extremely involved because of the strong nonlinear dependence of the second term in (6) on the coefficients and spin orbitals. The above problems can be avoided by an alternative approach that has been used earlier by Wang et al in the context of the MCTDH method for nuclei [32]. This approach is based on the observation that if one has a method to accurately solve the Schrödinger equation one can revert to the linear optimization equations (8)-(10) because the conjugate equation (9) has the same structure as the original Schrödinger equation and therefore can itself be accurately solved. In other words, in this approach the MCTDH/HF method is regarded as an efficient tool for solving the Schrödinger equation and the equation-of-motion for |χ(t) . The conceptual difference between this approach and the one discussed in the previous paragraph is schematically shown in figure 1. It is clear that the second approach is justified rigorously only in the case that the MCTDHF solution is close enough to the exact solution, i.e. N S must be chosen large enough so that the MCTDHF state is close to the exact state. In the limit N S = N , i.e. the time-dependent Hartree-Fock limit, the second approach will almost certainly fail because the nonlinearities in the Hartree-Fock propagation are not properly taken into account.
From a practical point of view the second approach is very attractive because one can avoid solving the highly involved control equations of the first approach arising from the 'conjugate' variables for c K and ϕ k . In addition, the second approach allows one to take advantage of the vast knowledge that has been gathered over the years about solving (8)-(10). However, a disadvantage of this approach is that the MCTDHF state must be a sufficiently accurate approximation to the exact state, i.e. the approach may require a large number of Figure 1. Schematic of the two different approaches for combining OCT and the MCTDHF method. The approach on the left-hand side is the one described first in the text, whereas the other approach is described second. In the present paper, we are using the approach on the right-hand side.
spin orbitals N S . One can argue that in this case an optimization using fewer spin orbitals in combination with nonlinear OCT would be unphysical anyway; however, this argument does not take into account situations in which the objective can be well described by the MCTDHF method even if the MCTDHF wavefunction is not highly accurate, e.g. if one is interested in creating a specific one-particle density distribution instead of a fully correlated N -particle wavefunction.
In the following, we are interested in optimizing state-to-state transitions using the second approach. In particular, we are interested in the performance of the approach, i.e. how many spatial orbitals N o are required for a successful optimization (monotonic increasing objective) and how the resulting controls compare to the controls obtained from an optimization based on the exact Schrödinger equation.

Results and discussion
We have chosen two different systems for testing the optimization approach presented above. In both cases, the system consists of two interacting particles in one spatial dimension. Since for such systems the interacting Schrödinger equation can still be solved exactly, these systems allow us to assess the quality and performance of the optimization based upon the MCTDHF method.

Results for a model of helium
We first consider a 1D model of helium described by the single-particle Hamiltonian and the soft-Coulomb interaction [33] V pp (x 1 ,  The second and the third term in (12) are, respectively, the soft-Coulomb interaction with the nucleus and the laser field in the dipole approximation. The electric field (t) is our control, and the objective is, similar to [34], to steer the system from the ground state to the first excited singlet state at T = 300 a.u. Since the Hamiltonian is spin-independent and the ground state is a singlet, the spatial wavefunction is always symmetric during the process.
As mentioned earlier, we use a real-space grid as a primitive basis to represent the orbitals. The grid spacing is d x = 0.2 a.u. and the grid size is ±200 a.u. The large box size guarantees that ionization channels are properly taken into account. For the penalty, we use with λ 0 = 1 and 2σ 2 = (0.232T ) 2 to ensure a smooth switch-on and switch-off of the pulse. To guarantee that the target state can be completely represented by the MCTDHF state with N o spatial orbitals, we use the ground state and excited state from an imaginary time propagation with N o spatial orbitals for the optimization at level N o . Figure 2 shows the value of J 1 = ψ(T )|ψ T ψ T |ψ(T ) during the iteration process for different numbers N o of spatial orbitals in the MCTDHF ansatz.    In order to clarify the mechanism of the state-to-state transition, we plot the power spectrum of the N o = 4 pulse in figure 4. The arrows in the figure indicate the transition 9 frequency and its odd multiples. The large peak at the transition frequency indicates that the main mechanism of the transition is a resonant process from the ground state directly to the first excited state. Other processes are also present, as one can see, e.g. from the energetically higherlying peaks in the spectrum, but their contribution is small. This observation is also supported by the fact that a resonant π-pulse with a sin 2 -envelope and the same maximal amplitude as the optimal pulse already leads to an overlap of J 1 = 0.918.
The previous findings suggest that the optimization in the present case corresponds to a large extent to an optimization in a two-level system. Consequently, the MCTDHF optimization procedure works well in this context if two conditions are satisfied: first, the two involved states must be represented well by the MCTDHF method and, second, their dynamics, as well as the dynamics of their superpositions, must be described accurately. Since for N o = 2 the states are very similar to the states for N o = 4 (their energy differs by <0.1 eV and the overlap is >0.997), we conclude that the failure of the optimization procedure for N o = 2 is caused by the wrong dynamics of the intermediate superposition states. This conclusion is also backed by the observation that a resonant π -pulse does not induce any significant transition from the ground state to the first excited state for N o = 2. One obtains a value of only J 1 = 0.064. Thus, in the present situation, the failure of the MCTDHF method for N o = 2 ultimately originates from the violation of the superposition principle which is caused by the nonlinearity of the MCTDHF method. A related observation was reported in [35] and has been called 'time-dependent state averaging' by the authors. For larger N o , the influence of the nonlinearity is reduced, and the superposition principle is, at least in the relevant subspace, restored.

Results for rubidium atoms in an optical lattice
The optimization and the dynamics in the previous subsection can be well understood in terms of a two-level system coupled to a laser field. To investigate the performance of the optimization method in a more challenging setting, we investigate the control of the external degrees-offreedom of two 87 Rb atoms in an optical lattice 2 . The external potential is given by [36,37] V ext (x k , V 0 (t), β(t), θ (t)) = −V 0 (t) cos 2 β(t) 2 (1 + cos 2 (k L x k − π/2)) and the atom-atom interaction is a 1D contact interaction, i.e.
In (15) k L = 2π/λ L is the wavevector of the laser with wavelength λ L = 810 nm. α(t) = {V 0 (t), β(t), θ (t)} are the controls and g/ h ≈ 400 kHz is the effective interaction strength. By changing the control parameters, one can switch the optical lattice from a double-well configuration with periodicity λ L /2 to a single-well configuration with periodicity λ L (see figure 5 and [37] for a more detailed description of the influence of V 0 (t), β(t) and θ(t)). Initially, the lattice is in the double-well configuration and the system is in the state where φ L/R 0 are states localized in the left/right well (see figure 5). To update the control, we use in combination with (14) for λ(t). Forα k (t), we choose the value of the control in the previous iteration. In contrast to (10), this update procedure does not penalize the magnitude of the controls, but only the change from one iteration to the next [38].

Process 1: merging two wells.
The objective is to merge the two wells in such a way that the atoms end up in two orthogonal spatial single-particle orbitals in the resulting well at T = 0.15 ms. In more detail, we choose the target state where φ M 0/1 are spatial single-particle orbitals very similar to the vibrational ground state and the first excited state in the well defined by V 0 / h = 87.5 kHz, β/π = 0.5 and θ/π = −0.474 (see figure 5). For the initial guess, we use the control sequence which is very close to the one used in [37], where practically the same process has been optimized based on the Schrödinger equation. This optimization provides a good performance test for the MCTDHF approach since, due to the short time of 0.15 ms, non-adiabatic processes are important for achieving a large value of J 1 and, consequently, several eigenstates of the instantaneous Hamiltonian are occupied during the process [37]. Figure 6 shows the value of the functional J 1 during the optimization. The target state is given by (19). The inset shows the complete optimization, whereas the main graph shows the last four iteration steps in more detail. As in the previous case, the optimization works very well for N o = 4. The maximal difference from the exact result is 0.01 in the first iteration step.
quality of the control sequence is practically the same for N o = 4 and N o = ∞. As shown in figure 7, the optimized controls in the two cases are also very similar, as in the previous subsection.

Process 2: creating entanglement.
The objective in the final optimization is to entangle the motion of the two atoms in the double-well configuration. For this purpose, we choose the entangled state as our target state at final time T = 0.30 ms. The states φ L/R 1 are given by where x L/R 0 are the centres of the left/right well and C is a normalization factor. The initial state is the same as in the previous optimization, i.e. equation (17). For the initial control sequence, we use (20) to switch from the double-well configuration to the single-well configuration and back. Although creating a vibrationally excited entangled state as in (21) is usually not desired in quantum information processing [39], we have chosen the state (21) since optimizing the resulting transition is a challenge for the MCTDHF approach for at least two reasons. First, as in process 1, the Hamiltonian changes considerably during the process and, second, the particle-particle interaction V pp plays a crucial role in the optimization because it is required for entanglement production. The second feature especially provides a stringent test for the performance of the MCTDHF optimization since it is the contribution from V pp to the time-derivative of |ψ MF (t) that is approximated in the MCTDHF method, and not the contribution from the single-particle Hamiltonian. As one can see from (5) and (6), the MCTDHF equations are exact for systems without a particle-particle interaction. The value of the functional J 1 during the optimization is plotted in figure 8. In order to avoid any differences in the results due to the extremely small initial value of J 1 (∼10 −6 ), we use the control sequence from the third iteration of the N o = 4 optimization (J 1 ∼10 −2 ) for the initial control sequence in the N o = ∞ optimization. For comparison, we also plot the results of an optimization with vanishing particle-particle interaction. In this case, the controls can change only the spatial shape of the orbitals in (17), but the number of configurations cannot increase during the process. As a consequence, the largest possible value of J 1 is 0.5 (see the appendix), which is reproduced well by the optimization. For the cases with a particle-particle interaction, the value of J 1 = 0.75 after 55 iterations is still small in comparison with the values of J 1 = 0.94 for helium and J 1 = 0.98 for process 1. Intuitively, it is not surprising that the optimization in the present case is more involved since the particle-particle interaction is essential for the transition, but can only be controlled in an indirect way. A direct control over the two-particle interaction, e.g. via Feshbach resonances, would certainly help to improve the optimization by allowing more flexibility. By how much it would speed up the optimization process is difficult to judge because little is known in general about the influence of different controls on the iterative optimization algorithm; a detailed study is beyond the scope of the present paper. From a control-theory perspective a control-dependent interaction would be desirable because it would remove the 'drift Hamiltonian' and thus, in principle, allow one to perform the entanglement operation arbitrarily fast (assuming that the system is controllable) [40].
In addition to the slower convergence of the iteration process, one observes a violation of the monotonic increase in some iteration steps for N o = 4. However, the total difference from the optimization with N o = ∞ never exceeds 3% and the maximal decrease from the previous iteration is 1%. These small violations are within the expected accuracy of the MCTDHF solution for N o = 4; moreover, after the short periods of decrease the N o = 4 results return to track the convergence curve of the exact solution. The N o = 4 results continue to rise smoothly and monotonically for at least 25 iterations beyond those shown in the figure, reaching J 1 = 0.80. Thus, we can conclude that for all processes investigated in this paper, the optimization based on the MCTDHF method produces useful results already for N o = 4.

Summary and conclusion
We have combined the MCTDHF method with OCT to optimize state-to-state transitions in quantum systems containing interacting particles. As we have explained in section 2, there are two different approaches to how the combination can be achieved. One approach is based directly on the MCTDHF equations-of-motion and uses nonlinear OCT. The other approach is based on the fact that any method that allows one to solve the Schrödinger equation sufficiently accurately can be used also to solve the corresponding linear OCT problem. In order to investigate the performance of the second approach in different settings, we have used the approach in section 3 to optimize state-to-state transitions in a 1D model of helium and in a system containing two rubidium atoms in a 1D optical lattice. For comparison, we have also optimized the same processes by solving the control equations exactly, which is still feasible for the systems discussed in this paper.
In all our test cases, the optimization based on the MCTDHF method works very well already for N o = 4 spatial orbitals. For instance, during the complete optimization the values of the objective never differ by more than J 1 = 0.03 from the exact values. In addition, the resulting controls for N o = 4 are also very similar to the ones obtained from the exact Schrödinger equation. For the optimization with N o = 6, we find perfect agreement between the results, while for N o = 2, the optimization procedure does not work. The case N o = 2 shows that in the present case more spatial orbitals are required than would be expected from the 'rule of thumb', which states that in general as many spatial orbitals as electrons are required for an accurate solution.
Since the optimization approach always works for a sufficiently large number of spin orbitals N S , the crucial question is how large N S must be in general. Clearly, for large values of N S the approach cannot be used in practice because the MCTDHF method scales as N S N with the number of particles and spin orbitals used. Judging in advance the required number of spin orbitals is difficult because it depends on the system and the process that one wants to optimize. However, it is possible to state at least two necessary conditions that must be satisfied for a successful optimization based on the MCTDHF method. First, the initial state and the target state must both be represented well by the MCTDHF method. Based on the numerical evidence in [15], we thus expect that if the states are eigenstates of an N -particle system at least N spatial orbitals are required. However, as shown by the failure of the optimization for N o = 2, this condition is not necessarily enough for a successful optimization. The second necessary condition is that it must be possible to steer the system with the given controls from the initial state to the target state in such a way that the exact time-derivative of the state |ψ(t) can always be approximated well by the time-derivative of the corresponding MCTDHF state. In other words, the initial state and the final state must be 'connected' by a low-dimensional subspace, which is accessible with the given controls and in which the MCTDHF dynamics approximates the exact dynamics well. Since the iterative optimization can take the state out of this subspace (or it already starts outside), the second condition again is not a sufficient condition for a successful optimization. In the final analysis, for each new application one must determine anew how large N S must be for a successful optimization. Depending on the optimization objective, the dimensionality of the problem and the system's Hamiltonian, we expect the approach to be useful for systems with N = 5 to N = 10 fermions. Despite the mentioned drawback that N S must be determined anew for different applications, the present applications clearly show how promising the scheme is for controlling interacting few-particle systems in a non-perturbative way.
Appendix. Upper bound for J 1 in process 2 with vanishing particle-particle interaction In this appendix, we show that the maximal value of J 1 in process 2 without a particle-particle interaction is 0.5. We start by noting that without a particle-particle interaction the number of configurations cannot change during the process as can be seen from (5). Thus, the state at final time T has the general form where ϕ a and ϕ b are some spatial orbitals resulting from the time evolution. We write ϕ a in the form with a L 0 = φ L 0 |ϕ a , . . . , and ϕ orth a labelling the part of ϕ a that is orthogonal to the four states φ L/R 0/1 . The orbital ϕ b is written in the same way with coefficients b L/R 0/1 . The overlap of the state at time T with the target state is given by Since the orbitals ϕ a/b are normalized, the absolute square value of the bracket on the rhs is 1 (Cauchy-Schwartz inequality). Thus, J 1 = | ψ T |ψ(T ) | 2 can never exceed a value of 0.5.