Global optimization for quantum dynamics of few-fermion systems
Xikun Li, Daniel Pęcak, Tomasz Sowiński, Jacob Sherson, Anne E. B. Nielsen
GGlobal optimization for quantum dynamics of few-fermion systems
Xikun Li, Daniel Pęcak, Tomasz Sowiński, Jacob Sherson, and Anne E. B. Nielsen ∗ Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark Institute of Physics, Polish Academy of Sciences, Aleja Lotnikow 32/46, PL-02668 Warsaw, Poland Max-Planck-Institut für Physik komplexer Systeme, D-01187 Dresden, Germany
Quantum state preparation is vital to quantum computation and quantum information processingtasks. In adiabatic state preparation, the target state is theoretically obtained with nearly perfectfidelity if the control parameter is tuned slowly enough. As this, however, leads to slow dynamics, itis often desirable to be able to do processes faster. In this work, we employ two global optimizationmethods to estimate the quantum speed limit for few-fermion systems confined in a one-dimensionalharmonic trap. Such systems can be produced experimentally in a well controlled manner. Wedetermine the optimized control fields and achieve a reduction in the ramping time of more than afactor of four compared to linear ramping. We also investigate how robust the fidelity is to smallvariations of the control fields away from the optimized shapes.
I. INTRODUCTION
Quantum optimal control is essential to manipulateand engineer complex quantum systems in quantum in-formation processing and quantum computation [1–3].The operations in experiments are often executed adi-abatically to guarantee the transition to the target statewith almost perfect fidelity [4]. The adiabatic process,however, needs to be done slowly, and it is therefore in-teresting to look for ways to achieve a speed up which isthe topic of the field of quantum optimal control [5, 6].The minimal allowed time for driving such transitionswith perfect fidelity is known as the quantum speed limit(QSL) [7, 8]. The quantum speed limit is a lower boundfor the duration in which the quantum system can becompletely steered to the target state [8–13]. For du-rations shorter than the quantum speed limit, defectsemerge that lead to a drop in fidelity between the tar-get state and the obtained state. The quantum optimalcontrol theory is important to obtain the quantum speedlimit [14] and has been applied using certain numericalmethods in many quantum systems like the NMR [15],Bose-Einstein Condensates [16–18] and spin chain mod-els [19, 20].Except for a few special cases in which analytical re-sults are available, one has to perform numerical calcula-tions, which are highly non-trivial, due to the high dimen-sionality of the Hilbert space. Generally, quantum con-trol theory relies on numerical techniques including thelocal optimization algorithms, such as Krotov, GRAPEand CRAB [21–24], as well as the global optimizationmethods like Differential Evolution (DE) [25–27] and co-variance matrix adaptation evolution strategy (CMA-ES) [28, 29]. In Ref. [30], it is proposed that numericaloptimization relies on an appropriate balance betweenlocal and global optimization approaches and problemrepresentation. When the quantum system is fully con- ∗ On leave from: Department of Physics and Astronomy, AarhusUniversity, DK-8000 Aarhus C, Denmark trollable and free of constraints, there are no traps in theform of suboptimal local extrema [31]. In such cases, thelocal algorithms are preferred as the computational costof local optimization methods is lower than that of theglobal ones. When the duration of the process is short orif there are constraints on the control field, the local al-gorithms are often stuck in the local suboptimal traps inthe quantum control landscape. For the low-dimensionalquantum system, the computational cost of multistart-ing the local optimization algorithms, which is able togive sufficiently good results, is comparable with that ofglobal ones. In Ref. [25], the local optimization meth-ods fail to obtain a satisfactory result determined by cer-tain threshold infidelity for quantum gates, though theglobal optimization methods succeed. The superiorityof global optimization methods are also highlighted forhigh-dimensional Hamiltonians studied in the Ref. [26].The ultimate goal is to fully control any many-bodyquantum system. In cold atom experiments, one caninfluence the interparticle interaction with an externalmagnetic field thanks to Feshbach resonances. Due tothe adiabatic change of the interaction it is possible toobtain for example a highly correlated state known as theTonks-Girardeau gas starting from the non-interactingstate [32–34]. Unfortunately, the full control of systemswith a large number of particles is very challenging. Apossible way to overcome the difficulty of such complexsystems is to fully control smaller physical systems anduse them to build the real many-body systems. A possi-ble candidate to serve this purpose are quantum systemsof a few ultra-cold atoms [35–41]. In the two-componentmixtures of fermions one can deterministically prepare asystem confining a well-established number of atoms withastonishing precision. The properties of few-body ultra-cold systems were also studied recently theoretically, in-cluding energy spectra and density profiles [42–50]. Thetwo different flavors in the mixture of same-mass fermionsare realized experimentally by using two different hyper-fine states of ultra-cold lithium Li. A natural way togeneralize this idea is to change the hyperfine states totwo completely different species, for example lithium andpotassium [51, 52]. Such an experiment on a two-flavor a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a r mixture of lithium and potassium on a many-body scalehas already been performed [53]. Recently, few-bodymass-imbalanced systems were broadly explored theoret-ically [54–60].In this paper, we employ two global optimizationalgorithms, CMA-ES [28, 29] and self-adaptive DE(SaDE) [26, 27], to numerically estimate the quantumspeed limit for few-fermion mass-imbalanced systems,and show the optimized control field for various dura-tions. We consider the fidelity of the final state withrespect to the target state as the fitness function to beoptimized. As a proof of concept we show how to fullycontrol a system of a few fermions and drive it from thenon-interacting state to the strongly correlated one. II. THE MODEL
We consider a two-flavor system of a few ultracoldfermions confined in a one-dimensional harmonic trap.Here we assume that the frequencies ω of the harmonictrap are the same for both flavors. Fermions of oppositeflavors interact in the ultra-cold regime via short-rangeforces modeled by the delta-like potential U ( x − x (cid:48) ) = gδ ( x − x (cid:48) ) , where g is the interaction strength [61]. Inthis approximation fermions of the same type do not in-teract as a consequence of the Pauli exclusion principle.Fermions of opposite flavors are fundamentally distin-guishable and may have the same or different masses (inthe following we denote the mass ratio µ = m ↑ /m ↓ ).The Hamiltonian of the mass-imbalanced system reads(see [54–56]): ˆ H = N ↓ (cid:88) i =1 (cid:20) − ∂ ∂x i + 12 x i (cid:21) + N ↑ (cid:88) j =1 (cid:34) − µ ∂ ∂y j + µ y j (cid:35) + g N ↓ ,N ↑ (cid:88) i,j =1 δ ( x i − y j ) . (1)All quantities are measured in appropriate harmonicoscillator units, i.e., positions are measured in unitsof (cid:112) (cid:126) / ( m ↓ ω ) , time in units of /ω , energies in (cid:126) ω ,and the interaction strength g is measured in units of ( (cid:126) ω/m ↓ ) / .We employ the exact diagonalization approach tostudy the dynamics of the few-fermion system (seeRef. [54] for details of the numerical method). The inter-action g is controlled experimentally with the help of themagnetic field B , so by changing the magnetic field intime, one can also modify the interaction g ( B ( t )) . Thus,it is convenient to treat the interaction strength g ( t ) asthe control field. The many-body spectrum for the sys-tem of three fermions is shown in Fig. 1. In this paper,we focus on the transformation from the ground stateof the non-interacting Hamiltonian ˆ H ( g = 0) to that of ˆ H ( g = 10) , where strong correlations are present. Byincreasing the interaction adiabatically, one can trans-fer the non-interacting state Ψ g =0 to the interacting E ne r g y Interaction g Ψ g=0 Ψ g=10 E ne r g y Interaction g Ψ g=0 Ψ g=10 E ne r g y Interaction g Ψ g=0 Ψ g=10 Figure 1. The energy spectrum of the system of three fermionsfor three different scenarios, namely 1Li-2Li, 2K-1Li, and 1K-2Li. The thick orange line gives the ground-state energy. Thenon-interacting ground state g = 0 and the strongly corre-lated target state g = 10 are marked with the red dots. Forthe 1Li-2Li and the 2K-1Li systems, there are, respectively, 2and 1 other states that have energies close to the ground stateenergy in the strong coupling limit. These states do, however,have different symmetries than the ground state, and the pop-ulation of these states remain zero throughout the dynamics.The energies and the interaction strengths are measured inunits of (cid:126) ω and ( (cid:126) ω/m ↓ ) / , respectively. one Ψ g =10 . Note that depending on whether the mass-imbalance is present in the system, and on the specificconfiguration, the ground state of the system might bequasi-degenerated in the strong interaction limit. For theequal-mass system there is a three-fold degeneracy, whilefor the system 1K-2Li the ground-state is not degener-ated. For the dynamics considered below, there is nocoupling to these additional low-energy states, becausethey have a different symmetry. The time scale for adi-abatic ramping is hence not determined by the energyof these states relative to the ground state, but by theenergy of the higher lying states relative to the groundstate.For a typical quantum control problem, the Hamil-tonian depends on a time-dependent control field H = H ( g ( t )) . We wish to optimize the fidelity as a fitnessfunction F ( g ( t ) , T ) = |(cid:104) Ψ g =10 |T exp( − i (cid:90) T H ( g ( t ))d t ) | Ψ g =0 (cid:105)| (2)with the time evolution driven by the control field g ( t ) ,where T is the time-ordering operator and the duration T is discretized with time-step size δt = 10 − for numer-ical evaluation of the time evolution. It is worth notingthat the most time-consuming part during the numer-ical calculation is the time evolution which consists ofa long sequence of evaluations of exponentials of largeHamiltonian matrices (see Sec. V). Therefore the max-imal number of iterations is set to be 500 for CMA-ESand SaDE.To perform this optimization systematically, we choose g ( t ) to be decomposed into a truncated Fourier basis(CRAB method, see [16, 21]) g ( t ) = g ( t ) (cid:34) N ( t ) N c (cid:88) n =1 ( A n cos( ω n t ) + B n cos( ω n t )) (cid:35) , (3)where g ( t ) is the initial guess of the control field and N ( t ) = T / t ( T − t ) is a time-dependent function to fixthe initial and final control field value to be g ( t = 0) = 0 and g ( t = T ) = 10 . { A n , B n } are Fourier coefficients and ω n = 2 πn (1 + r n ) /T are “randomized” Fourier harmon-ics, and r n ∈ [0 , . The choice of the cut-off number N c of the Fourier basis may vary from one case to another:it may depend on the Hamiltonian, the fitness functionand the optimization algorithm. The parameter space(search space) consists of the Fourier coefficients and har-monics { A n , B n , ω n } , which can be numerically obtainedusing an optimization method, e.g., the simplex method,gradient-based strategies and global optimization algo-rithms. In this paper, we restrict to non-negative inter-actions g ( t ) ≥ by simply putting the negative valuesof g ( t ) to be zero, in which case the local optima in thequantum control landscape are usually not global optima. III. GLOBAL OPTIMIZATION
We employ two evolutionary computation techniques,CMA-ES and SaDE, as global optimization methods. Wecompare CMA-ES with SaDE for three different systems:( i ) the mixture of three Li atoms with two different hy-perfine states ( N ↑ = 1 , N ↓ = 2 , µ = 1 ); ( ii ) the mixtureof one K atom and two Li atoms ( N ↑ = 1 , N ↓ = 2 , µ = 40 / ); ( iii ) the mixture of two K atoms and one Li atom ( N ↑ = 2 , N ↓ = 1 , µ = 40 / ). We numeri-cally estimate the duration T QSL for which the fidelity F ( T QSL ) = 0 . . We then compare the control field g ( t ) for various durations T and depict the deviations betweenthe optimized control field and non-optimal ones. Forsimplicity, we will present results on the 1K-2Li systemunless stated otherwise.CMA-ES and SaDE are variants of evolution strategy(ES) and DE, respectively. Both ES and DE belong tothe class of evolutionary algorithms and are stochastic,derivative-free algorithms for global optimization of fit-ness functions. An evolutionary algorithm works througha loop of variations (including recombination and muta-tion) and selection in each iteration (also called gener-ation). New candidates are generated by variation ofcurrent parent individuals in each iteration. Then somecandidates are selected, based on their fitness, to be par-ents for the next generation. In this way, search pointswith better and better values of the fitness function aregenerated over the sequence of iterations. Table I. Fidelity for various combinations of cut-off number N c and population size N p with duration T = 0 . using theCMA-ES method for the 1K-2Li system. The maximal valueof fidelity is indicated by the bold font, and the correspondingpair of values is N c = 15 and N p = 60 , respectively. N c N p
20 0.5995 0.5869 0.5933 0.585740 0.6021 0.5984 0.5924 0.584860 0.5968 0.6015
In CMA-ES, new search points (parameter vectors) aresampled according to a multivariate normal distributionin the parameter space. CMA-ES begins with a randomlyinitiated population of search points in the parameterspace with initial mean and covariant matrix. The popu-lation size N p is the number of search points in each itera-tion. In the selection and recombination step, the searchpoints with the best m fitness, where m is the parent sizeand not larger than the population size, are chosen as theparents to update the new mean, step-size and covariantmatrix. Recombination amounts to selecting a new meanvalue for the multivariate normal distribution. In the mu-tation step, the parameter vectors are further added byrandom vectors with zero mean and updated covariancematrix. The fitness function evolves iteratively towardsits optimal state. In contrast to most other evolutionaryalgorithms, CMA-ES is quasi parameter-free: one needsonly to randomly choose an initial value of the step-size.In addition, the population size N p does not depend onthe dimension of the parameter space and can hence bechosen freely (which is in contrast to DE). In general,large population sizes help to circumvent local optima,while small population sizes usually lead to faster conver-gence. Therefore, a trade-off between the computationalcost and performance needs to be carefully determined ifthe computational time for each iteration is considerablylong. See [29] for a review of the CMA-ES method.In Table I we show the values of fidelity for variouscombinations of N c and N p with duration T = 0 . ob-tained using the CMA-ES method for the 1K-2Li sys-tem. The maximal fidelity in Table I is F = 0 . with( N c = 15 , N p = 60 ). It is hardly possible to infer the opti-mal combination ( N c , N p ) to obtain the maximal fidelityfor arbitrary durations, nor reasonable to try all possi-ble combinations of ( N c , N p ) as the computation cost ishuge. Therefore we fix the population size to N ES p = 60 and the cut-off number of the Fourier basis N c = 15 forall durations in the three different few-fermion systems.In SaDE, an initial population of parameter vectors(called genome or chromosome) is randomly sampled.Then the mutant chromosomes are obtained from thedifferential mutation operation (origin of the term “DE”).In the mutation step, three mutually exclusive parametervectors are generated randomly. The new set of param-eter vectors are generated by adding one of those threevectors to the difference between the other two vectorswith the mutation scale factor S which controls the differ-ential variation. In the recombination step, an offspringis formed by recombining the original and those mutantchromosome in a stochastic way, where the crossover rate Cr controls the probability of recombination. In the se-lection step, comparison of values of the fitness functiondetermines whether the offspring or the original chromo-some survives to the next generation. See Ref. [63] for areview of DE.In the conventional DE algorithm, there are two freeparameters: S , Cr which are fixed through the itera-tions. In SaDE, however, S and Cr are adapted in eachiteration to enhance the convergence rate for the high-dimensional optimization problem, and to obtain betterquality solutions more efficiently, compared with the con-ventional DE. A reasonable value of N p for DE and itsvariants is usually chosen between D and D ( D = 3 N c is the dimension of the parameter space), as suggested inthe field of evolutionary computation science [62]. Note,however, that this is not tested in great detail for phys-ically motivated quantum systems. In this work, we fixthe population size N p = 5 D and set the cut-off numberto N c = 5 , thus N DE p = 75 , for SaDE to reduce the com-putational cost. Details of SaDE can be found in [26, 27]. IV. RESULTS
First, we numerically compute the fidelity obtained foroptimized and non-optimal control fields, when the dura-tion T of the time evolution is fixed to a certain value. Weuse the constraint that the interactions must be nonneg-ative at all times. We consider two typical non-optimalrampings: linear ramping and exponential ramping. The Figure 2. Fidelity versus duration for the 1K-2Li system usingdifferent control fields: optimized (blue, solid line), exponen-tial (red, dashed) and linear (black, dash-dot). The curve forthe optimized control field is obtained using the CMA-ES al-gorithm. The stars mark the shortest durations for which afidelity of F = 0 . is reached: T Opt = 2 . , T Exp = 6 . and T Lin = 11 (the time unit is ω − ). latter is described by g ( t ) = g max − e t/τ − e T/τ , (4)where τ = T / (In [17], the exponential ramping withparticular values of T and τ is representative of quasi-adiabatic ramping of the lattice depth used in experi-ments of optical lattices). In Fig. 2, we take the 1K-2Lisystem and compare the fidelity (as a function of dura-tion) for the optimized ramping obtained using the CMA-ES method with that for exponential ramping and linearramping. The shortest duration with fidelity F = 0 . is T Exp = 6 . for the exponential ramping and T Lin = 11 forthe linear ramping, whereas T Opt = 2 . for the optimizedramping (note the time unit is /ω ). Thus the shortestduration obtained by the optimized control is approxi-mately one fourth of that using linear ramping, and twofifth of that using the exponential ramping, in the 1K-2Lisystem. It means that by controlling the interaction inthe optimized way, one can significantly reduce the timeof preparing the system in the strongly correlated state.In Fig. 3, we show the optimized control field, the con-trol field of linear ramping and exponential ramping fordifferent durations. For very short duration ( T = 0 . ),the optimized control field shows a few large oscillations,while g ( t ) remains to be zero for most of the time (Notethe requirement that g ( t ) ≥ ). As the duration ap-proaches T Opt = 2 . , which is an estimate of the quantumspeed limit T ≈ T QSL (Fig. 3b), more oscillations emergeto make the transformation as fast as possible and thefidelity as large as possible. The deviation between theoptimized control field and the linear and exponentialones is thus large for such durations. When the durationis much larger than T Opt , e.g., T = 5 (Fig. 3c) and T = 15 (Fig. 3d), the highly oscillating components are no longer Figure 3. Control field for different durations in the 1K-2Lisystem. (a) T = 0 . , (b) T = 2 . , (c) T = 5 and (d) T = 15 (the time unit is ω − ). The blue solid lines are the optimizedcontrol fields obtained using the CMA-ES method, the reddashed lines are exponential rampings, and the black dash-dot lines are linear rampings. necessary and the deviations between different rampingsare much smaller than the cases in which T < T
Opt . Alsothe fidelities do not depend for the longer durations thatmuch on the way we approach the strong interaction.This observation means that for long enough times, theexact shape of the control field does not matter. This iswhy slow enough processes are quasi-adiabatic and thequantum state can be transferred when carefully man-aged. Since we require the interaction strength to benon-negative g ≥ , the lower bound of the control field iszero, but there is no upper bound. As mentioned above,if g ( t ) (as given in Eq. (3)) is negative in some time inter-val, we put g ( t ) equal to zero in that interval, as shown inFig. 3. Such processes introduce the possibility of sharppeaks in the control fields and thereby high Fourier com-ponents. Given the particular experimental constraintsof the system, one may perform an analogue optimizationwith appropriate constraints included.The comparisons between CMA-ES and SaDE forthree different few-fermion systems are shown in Fig. 4.To reduce the computational cost, we set the maximalnumber of iterations to 500. In addition, we set up haltcriteria for the CMA-ES and SaDE methods. For both ofthem, the calculations stop if the distance between theminimal and the maximal fidelity in the population issmaller than a threshold value ( Error = 10 − ). The topleft panel depicts the infidelity IF = 1 − F as a functionof duration T for three different systems using the CMA-ES and SaDE methods. For short durations T < . , theinfidelities obtained using CMA-ES and SaDE are veryclose. For long durations T > . , however, the infideli-ties obtained by CMA-ES are smaller than SaDE (apartfrom T = 4 for the 1Li-2Li system), which means the -4 -2 Figure 4. Comparing CMA-ES and SaDE. The top left panelshows the infidelity ( IF = 1 − F ) as a function the duration T using CMA-ES (square) and SaDE (asterisk) for 1Li-2Li(blue, solid), 1K-2Li (red, dashed) and 2K-1Li (black, dash-dot), where the scale of vertical axis is logarithm. In the topright panel, the red squares are numerical data, which areobtained using CMA-ES method, of infidelity IF as a functionof duration T for the 1K-2Li system, while the blue solid linesare the fitting obtained using the function k cos ( Tπ β ) for thenumerical results, where k and β are free fitting parameters.In the bottom panel, the comparisons of log -infidelity versusiterations in the 1K-2Li system are shown between CMA-ES (cyan, solid) and SaDE (magenta, dashed) for durations T = 2 . (bottom left) and T = 5 (bottom right) with the timeunit ω − . performance of CMA-ES is better than SaDE in terms ofthe best infidelity with the specific parameters ( N p , N c )used for both global optimization methods in this work.It is worth noting that the estimate of QSL is approx-imately proportional to the inverse of energy gap, i.e., T QSL ∼ π/ ∆ , where ∆ ≈ is the energy gap to the near-est coupled excited state for all three few-fermion sys-tems (See Fig. 1). This fact agrees with the conclusionobtained in Ref. [21]. In the top right panel, we show thenumerical results of infidelity IF as a function of dura-tion T for the 1K-2Li system using the CMA-ES methodand the curve fitting using the cosine square function k cos ( T π β ) with two free fitting parameters ( k, β ) . Thesingle cosine square function does not fit well for the nu-merical results. Such deviations or discrepancies are alsofound in different quantum systems [13, 18, 64].In the bottom panel of Fig. 4, we demonstrate the log -infidelity versus iterations (also called generation in evo-lutionary computation) for duration T = 2 . (bottomleft) and T = 5 (bottom right) in the 1K-2Li system.We observe that the log -infidelity starts converging af-ter a certain number of iterations (several tens to severalhundreds) for the CMA-ES method, while the “staircase”pattern emerges in most cases of the SaDE method. Fromthe lower panel in Fig. 4, the CMA-ES method performs Figure 5. Effects of deviations from the optimized controlfields, which is obtained mixing the optimized control fieldand linear ramping with certain weight. In the top panel,the deviation between the fidelity of the modified optimizedcontrol field F w and that of the optimized control field F o areshown as a function of weight for T = 2 . (left) and T = 5 (right) with the time unit ω − . In the bottom panel, both theoptimized control field (blue, solid line) and its modificationwhose weight is (red, dash-dot) are shown for T = 2 . (left) and T = 5 (right). The fidelity of the optimized controlfield F o and that of the modified one F w are also given, and thedifference between them is of order − for both durations. better than the SaDE method in terms of the best in-fidelity and convergence rate. The reason why CMA-ES performs better than SaDE might be that the cut-offnumber of the Fourier basis of CMA-ES N ES c = 15 islarger than that of SaDE N DE c = 5 , such that the searchspace of CMA-ES is larger than that of SaDE. Note,however, the population size of CMA-ES ( N ES p = 60 )is smaller than that of SaDE ( N DE p = 75 ).Since the numerical calculations studied in this workare considerably time-consuming, the convergence rateof the optimization method is one of the most impor-tant considerations. Therefore, CMA-ES is more prefer-able than SaDE because CMA-ES does not need to use alarge population size [29]. When the computational costis small, the primary consideration is whether a satisfac-tory result, e.g., a certain threshold value of the fitnessfunction, is achieved by the optimization method [25].From the experimental point of view a very importantquestion arises: how sensitive is the final fidelity withrespect to small changes in the optimal control field? Todemonstrate the robustness of the optimized control fieldobtained using global optimization methods, we depict inFig. 5 the comparisons between optimized control fieldsand the corresponding modification for the durations T =2 . and T = 5 . The modified optimized control fields areobtained by mixing the optimized control field and linearramping with different values of weight w ∈ [0 . , . , i.e., g ( t ) = wg Lin ( t ) + (1 − w ) g Opt ( t ) . As shown in the toppanel in Fig. 5, the differences between the fidelity of theoptimized control field F o and that of the modification F w increases as the weight of linear ramping grows (Notethat F w < F o ). For a special case where the weight is , as shown in bottom panel in Fig. 5, the differencesbetween the fidelities are of order − for T = 2 . and T = 5 . It means that the optimized control fields obtainedare robust to the imperfection or external noise which isnaturally present in experiments. Thus, the scheme maybe especially relevant in future experiments. V. DISCUSSION
As two of the most promising evolutionary algorithms,the CMA-ES method and the SaDE method outperformother evolutionary computational algorithms and localoptimization algorithms for high-dimensional optimiza-tion problems in certain quantum systems [28, 29, 63].Both algorithms, however, have its own advantagesand disadvantages, and the preferences may vary fromone case to another. The CMA-ES method is quasiparameter-free, while the SaDE method requires moreinitial parameters to be determined by the user. Forthe SaDE method, as mentioned in Sec. III, the popu-lation size N p is fixed to be D which is suggested as alower bound and tested in great detail in the field of com-putational science (but not as extensively for physicallymotivated quantum systems). For the CMA-ES method,there is no guide for choosing the value of N p , thus wechoose N p = 60 which is large enough to guarantee thatthe fidelity is saturated. Therefore, in general, the pop-ulation size of SaDE is larger than that of CMA-ES,especially in high-dimension parameter space, thus thecomputational time of SaDE is longer than CMA-ES. Asfor the convergence rate, in general, CMA-ES convergesfaster than SaDE. The slow convergence of SaDE is de-picted in Fig. 4 (bottom right) where the width of thestaircase indicates the stagnation of the SaDE method.Note, however, if the N c of SaDE is the same as thatof CMA-ES, the final fidelity obtained using the SaDEmethod is generally larger than that using the CMA-ES method, though the computational time of SaDE ismuch longer than CMA-ES. For instance, suppose that( N c = 15 , N p = 5 D ) is taken for the SaDE method, whichis the same as CMA-ES, then the computational time ofSaDE is approximately three times larger than that ofCMA-ES.The calculations are performed in parallel using Mat-lab
R2017a on cluster (Intel Xeon E5-2680 CPU with28 cores and 251 GB RAM). Take CMA-ES for instance,the computational time of the CMA-ES method over 500iterations is about 31h for T = 1 . and 71h for T = 2 . .For the same process duration T and number of itera-tions, the computational time ratio of CMA-ES to SaDEis roughly
60 : 75 , which is the ratio of population size.This is because the maximal number of cores in the clus-ter is 28. If the number of cores is larger than the popula-tion size, then the computational time ratio of CMA-ESto SaDE is approximately 1:1.
VI. CONCLUSION
We have given first numerical estimates of the quan-tum speed limit for three different few-fermion systemsconfined in a one-dimensional harmonic trap using theCMA-ES and the SaDE methods, and shown that theshortest duration obtained employing optimized, nonadi-abatic processes is much faster than in the case of linearramping and exponential ramping. One can achieve atleast double speed-up in obtaining the target three-bodyground state by using our optimized approach comparedto the exponential ramping (see Fig. 2). Since the Hilbertspace increases greatly with the number of particles, thespeed-up might increase even further for the systems withmore than three particles. We observed that for dura-tions shorter than the estimate of the quantum speedlimit the optimized fields are of oscillation type, whilefor longer times, the optimized fields do not change dras-tically in time, which is analogous to the linear rampingand the exponential ramping. We have compared theperformance of the CMA-ES and the SaDE methods,and found that the performance of CMA-ES is betterthan SaDE in terms of the best infidelity and conver-gence rate for the parameters considered in this paper. In addition, we have explained the advantages and dis-advantages of the CMA-ES method, as well as the SaDEmethod, and the preference varies from one case to an-other. Moreover, we also demonstrated the robustnessof the optimized control field to minor variation. Thisstability of the above scheme on small variations let usbelieve that the obtained optimized fields are not justa purely numerical prediction, but can be useful in thenoisy laboratory environment.The present work shows the encouraging result thatcontrol theory can be used to obtain a significant speedup in producing a target state with the same symmetryas the initial state. As a next step, it would be very inter-esting to investigate how one can design control protocolsto produce any of the low energy states in the spectrumwith high fidelity starting from the ground state of thenoninteracting system.
ACKNOWLEDGMENTS
We would like to thank Nikolaj T. Zinner for discus-sions. This work has in part been supported by theVillum Foundation. DP and TS acknowledge supportfrom the (Polish) National Science Center Grants No.2016/21/N/ST2/03315 (DP) and 2016/22/E/ST2/00555(TS). JS acknowledges support from ERC. XL thanks theMax Planck Institute for the Physics of Complex Systemsfor hospitality during visits to the institute. [1] D. D’Alessandro,
Introduction to Quantum Control andDynamics , (Chapman & Hall/CRC 2007).[2] C. Brif, R. Chakrabarti, and H. Rabitz, New. J. Phys. , 075008 (2010).[3] V. F. Krotov, Global Methods in Optimal Control Theory (Marcel Dekker, New York, 1996).[4] T. Gericke, F. Gerbier, A. Widera , S. Fölling, O. Mandeland I. Bloch, J. Mod. Opt. , 735 (2007).[5] R. Modak, L. Vidmar, and M. Rigol, Phys. Rev. E ,042155 (2017).[6] S. Deng, P. Diao, Q. Yu, A. del Campo, and H. Wu,Phys. Rev. A , 013628 (2018).[7] L. I. Mandelshtam and I. E. Tamm, J. Phys. (USSR) ,249 (1945).[8] M. Gajdacz, K. K. Das, J. Arlt, J. F. Sherson, andT. Opatrný. Phys. Rev. A , 062106 (2015).[9] L. B. Levitin and T. Toffoli, Phys. Rev. Lett. , 160502(2009).[10] S. Deffner and E. Lutz, Phys. Rev. Lett. , 010402(2013).[11] A. del Campo, I. L. Egusquiza, M. B. Plenio andS. F. Huelga, Phys. Rev. Lett. , 050403 (2013).[12] M. G. Bason, et al. , Nat. Phys. , 147 (2012).[13] J. J. W. H. Sørensen, et al. , Nature , 210 (2016).[14] S. Lloyd and S. Montangero, Phys. Rev. Lett. ,010502 (2014).[15] S. J. Glaser, et al. , Eur. Phys. J. D , 279 (2015).[16] P. Doria, T. Calarco, and S. Montangero, Phys. Rev. Lett. , 190501 (2011).[17] S. Rosi, et al ., Phys. Rev. A , 021601 (2013).[18] S. van Frank, et al ., Scientific Reports , 34187 (2016).[19] T. Caneva, T. Calarco, and S. Montangero, New J. Phys. , 093041 (2012).[20] D. Burgarth, et al ., Phys. Rev. A , 040303 (2010).[21] T. Caneva, T. Calarco, and S. Montangero, Phys. Rev.A , 022326 (2011).[22] S. E. Sklarz and D. J. Tannor, Phys. Rev. A , 053619(2002).[23] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen,and S. J. Glaser, J. Magn. Reson. , 296 (2005).[24] S. Machnes, D. J. Tannor, F. K. Wilhelm, and E. Assé-mat, arXiv:1507.04261 (2015).[25] E. Zahedinejad, S. Schirmer, and B. C. Sanders, Phys.Rev. A , 032310 (2014).[26] P. Palittapongarnpim, P. Wittek, E. Zahedinejad,S. Vedaie, B. C. Sanders, Neurocomputing , 116(2017).[27] J. Brest, S. Greiner, B. Boškovićc, M. Mernik, and V.Žumer, IEEE Trans. Evol. Comput., , 646 (2006).[28] O. M. Shir, J. Roslund, D. Whitley, and H. Rabitz,arXiv:1112.4454 (2011).[29] N. Hansen, "The CMA evolution strategy: a comparingreview", Towards a new evolutionary computation. Ad-vances on estimation of distribution algorithms , pp 75-102. (Springer 2006).[30] J. J. W. H. Sørensen et al ., arXiv:1802.07521 (2018). [31] H. A. Rabitz, M. Hsieh, and C. M. Rosenthal, Science , 1998 (2004).[32] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling,I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch,Nature , 277 (2004).[33] T. Kinoshita, T. Wenger, and D. S. Weiss, Science ,1125 (2004).[34] E. Haller, M. Gustavsson, M. J. Mark, J. G. Danzl,R. Hart, G. Pupillo, and H. C. Nägerl, Science , 1224(2009).[35] F. Serwane, G. Zürn, T. Lompe, T. B. Ottenstein,A. N. Wenz and S. Jochim, Science , 6027 (2011).[36] G. Zürn, F. Serwane, T. Lompe, A. N. Wenz, M. G. Ries,J. E. Bohn, and S. Jochim, Phys. Rev. Lett. , 075303(2012).[37] G. Zürn, A. N. Wenz, S. Murmann, A. Bergschneider,T. Lompe, and S. Jochim, Phys. Rev. Lett. , 175302(2013).[38] A. N. Wenz, G. Zürn, S. Murmann, I. Brouzos, T. Lompe,and S. Jochim, Science , 457 (2013).[39] S. Murmann, A. Bergschneider, V. M. Klinkhamer,G. Zürn, T. Lompe, and S. Jochim, Phys. Rev. Lett. , 080402 (2015).[40] S. Murmann, F. Deuretzbacher, G. Zürn, J. Bjerlin, S. M.Reimann, L. Santos, T. Lompe, and S. Jochim, Phys.Rev. Lett. , 215301 (2015).[41] A. M. Kaufman, B. J. Lester, M. Foss-Feig, M. L. Wall,A. M. Rey, and C. A. Regal, Nature , 208 (2015).[42] T. Sowiński, M. Gajda, and K. Rz¸ażewski, Europhys.Lett. , 26005 (2015).[43] E. J. Lindgren, J. Rotureau, C. Forssén, A. G. Volosniev,and N. T. Zinner, New J. Phys. , 063003 (2014).[44] P. D’Amico, and M. Rontani, Phys. Rev. A , 043610(2015).[45] S. E. Gharashi, and D. Blume, Phys. Rev. Lett. ,045302 (2013).[46] T. Grining, M. Tomza, M. Lesiuk, M. Przybytek, M. Mu-sial, P. Massignan, M. Lewenstein, and R. Moszynski,New J. Phys. , 115001 (2015). [47] J. Decamp, P. Armagnat, B. Fang, M. Albert, A. Min-guzzi, and P. Vignolo, New J. Phys. , 055011 (2016).[48] F. Deuretzbacher, D. Becker, J. Bjerlin, S. M. Reimann,and L. Santos Phys. Rev. A , 013611 (2014).[49] Li Yang, Liming Guan, and Han Pu, Phys. Rev. A ,043634 (2015).[50] L. Yang, and X. Cui, Phys. Rev. A , 013617 (2016).[51] E. Wille, F. M. Spiegelhalder, G. Kerner, D. Naik,A. Trenkwalder, G. Hendl, F. Schreck, R. Grimm,T. G. Tiecke, J. T. M. Walraven, S. J. J. M. F. Kokkel-mans, E. Tiesinga, and P. S. Julienne, Phys. Rev. Lett. , 053201 (2008).[52] T. G. Tiecke, M. R. Goosen, A. Ludewig, S. D. Gensemer,S. Kraft, S. J. J. M. F. Kokkelmans, and J. T. M. Wal-raven, Phys. Rev. Lett. , 053202 (2010).[53] M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. M. Wal-raven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt,M. Knap, E. Demler, Science , 96 (2016).[54] D. Pęcak, M. Gajda, and T. Sowiński, New J. Phys. ,013030 (2016).[55] D. Pęcak, and T. Sowiński, Phys. Rev. A , 042118(2016).[56] D. Pęcak, M. Gajda, and T. Sowiński, Few-Body Syst. , 159 (2017).[57] M. A. García-March, A. S. Dehkharghani, N. T. Zinner,J. Phys. B: At. Mol. Opt. Phys. , 075303 (2016).[58] N. J. S. Loft, A. S. Dehkharghani, N. P. Mehta,A. G. Volosniev, and N. T. Zinner, Eur. Phys. J. D ,65 (2015).[59] N. L. Harshman, M. Olshanii, A. S. Dehkharghani, A. G.Volosniev, S. G. Jackson, and N. T. Zinner, Phys. Rev.X , 041001 (2017).[60] A. G. Volosniev, Few-Body Systems . 54 (2017).[61] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[62] R. Storn, and K. V. Price, J. Global Optimiz., , 341(1997).[63] S. Das, and P.N. Suganthan, IEEE. Trans. Evol. Comput. :4-31 (2011).[64] M. Bukov, et al.et al.