Preparation of ordered states in ultra-cold gases using Bayesian optimization
Rick Mukherjee, Frederic Sauvage, Harry Xie, Robert Löw, Florian Mintert
PPreparation of ordered states in ultra-cold gasesusing Bayesian optimization
Rick Mukherjee , Fr´ed´eric Sauvage , Harry Xie , Robert L¨ow and Florian Mintert Department of Physics, Imperial College London, SW7 2AZ, London, UK
5. Physikalisches Institut and Center for Integrated Quantum Science andTechnology, Universit¨at Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, GermanyE-mail: [email protected]
September 2019
Abstract.
Ultra-cold atomic gases are unique in terms of the degree of controllability,both for internal and external degrees of freedom. This makes it possible to usethem for the study of complex quantum many-body phenomena. However in manyscenarios, the prerequisite condition of faithfully preparing a desired quantum statedespite decoherence and system imperfections is not always adequately met. To pavethe way to a specific target state, we implement quantum optimal control based onBayesian optimization. The probabilistic modeling and broad exploration aspects ofBayesian optimization are particularly suitable for quantum experiments where dataacquisition can be expensive. Using numerical simulations for the superfluid to Mott-insulator transition for bosons in a lattice as well as for the formation of Rydbergcrystals as explicit examples, we demonstrate that Bayesian optimization is capableof finding better control solutions with regards to finite and noisy data compared toexisting methods of optimal control.
1. Introduction
In this paper, the focus is on creating spatially ordered states in two different ultra-coldsystems, atoms in optical lattice and highly excited Rydberg atoms, as testbeds forBayesian optimization. Optical lattices [1, 2] provide a natural and versatile platformto study strongly correlated condensed matter systems [3, 4], with implementationsin quantum information processing [5] and quantum metrology [6, 7, 8]. Likewise theextraordinary properties of Rydberg atoms [9] in combination with ultra-cold gasesestablish a facile way to realize strongly interacting many-body systems [10, 11, 12].Recently Rydberg atoms trapped in tweezers have been used to produce some ofthe largest non-classical states [13] and have emerged as a serious competitor for therealization of quantum computer [14]. While adiabatic preparation still remains anintuitive and straightforward way to create these interesting states, it requires longtimescales especially for large systems with diminishing energy gaps. This makes it a r X i v : . [ qu a n t - ph ] M a y reparation of ordered states in ultra-cold gases using Bayesian optimization reparation of ordered states in ultra-cold gases using Bayesian optimization M iterations Q(cid:88)an(cid:87)(cid:88)me(cid:91)pe(cid:85)imen(cid:87)Op(cid:87)imi(cid:93)e(cid:85)
ExperimentNumerical
Simulation
Q(cid:88)an(cid:87)(cid:88)me(cid:91)pe(cid:85)imen(cid:87)Op(cid:87)imi(cid:93)e(cid:85)
Quantum system
Initialization
Next
Building the modelUpdating the modelDecision rule
Bayesian Optimisation F ( ✓ m ) ✓ ✓ m Figure 1: Quantum optimal control scheme: F ( θ m ) is the FoM to maximize whichdepends on input parameters θ m that describes the control fields of the quantumsystem. F ( θ m ) is evaluated either numerically or experimentally. The optimizationtask is initiated by evaluating this FoM for a set of random control parameters θ . Atevery stage, the Bayesian optimizer suggests the next set of parameters for thequantum system based on previous evaluations of the FoM. This process is repeatedover M iterations.
2. Quantum Optimal Control formalism
Historically, quantum optimal control was first applied to manipulate chemical reactionsby selectively breaking bonds in molecules using lasers [28, 29, 30, 31, 32]; more recentexamples involve quantum metrology [33, 34], quantum computing [35, 36, 37, 38],control over few qubits [39, 40, 41, 42] as well as many-body systems [43, 44, 45].Within the framework of optimal control, the entire optimization task involvesmaximizing the FoM, F ( θ ) with respect to the set of parameters θ until an optimumsolution, θ opt is reached, as expressed below, θ opt = arg max θ F ( θ ) . (1) θ = [ θ , θ , . . . , θ P ] is a vector of size P that represents the parametrization of theinput control fields. For time-dependent control, this parametrization can be done usingspecial functions [46], Fourier series [45] as well as piece-wise constant functions [47].As illustrated in Fig. 1, the FoM is evaluated for some random set of initial parameters θ . This is passed onto the optimizer which in turn suggests the next set of parameters,for which the FoM is again evaluated thus repeating the optimization task iteratively.At iteration m , F ( θ m ) is obtained either from a real experiment or through numerical reparation of ordered states in ultra-cold gases using Bayesian optimization θ opt well before exhausting theexperimental/numerical resources. Although experimental constraints are not explicitlyincluded in Eq. (1), they can be incorporated either in the FoM directly or in the choiceof the parametrization.A variety of optimization algorithms are available and are broadly classified basedon the ability to evaluate the gradient of F ( θ ) with respect to θ . Gradient basedalgorithms such as the Krotov [48, 49] and the gradient ascent pulse engineering (GRAPE) [47] methods have been successfully applied to numerical simulations.Other gradient (or approximate gradient) based methods such as finite-differences [50], simultaneous perturbation stochastic approximation (SPSA) [38, 51] and an hybridquantum-classical approach [37, 38] can be applied to experiments. These gradientbased methods are considered as local optimization methods whose convergence is fastas long as the optimization landscape is well-behaved. However, these methods getcompromised by the presence of any local minima and plateaus in the optimizationlandscape [52, 53, 54]. Alternatively, gradient-free algorithms are able to explore theoptimization landscape more globally than gradient based methods making them lessvulnerable to local minima. Among them, Nelder-Mead has been extensively usedin the context of the chopped random basis (CRAB) method [45, 55, 56] as well asindependently [34, 35, 57]. It has the advantage of simplicity but can still get trappedin local minima and its convergence is limited by the presence of noise in the observations[51]. Other popular non-gradient methods include evolutionary algorithms [30, 52, 58]and the recently introduced reinforcement learning techniques [53, 59, 60]. Thesemethods usually require large number of iterations to find the optimal solution.Bayesian Optimization (BO) is a non-gradient based method offering an appealingalternative as it cleverly selects the next set of parameters to evaluate at each step of theoptimization. This can lead to significant reduction in the number of iterations neededfor convergence while performing global optimization. BO has the additional advantageof integrating probabilistic elements of data acquisition which can be exploited to furtherincrease its efficiency [61]. These characteristics make it especially well suited to performoptimal control on quantum experiments but also when numerical simulations of thesystem are time-consuming.
3. Bayesian Optimization
The essential steps of BO used for quantum optimal control are shown in Fig. 1 andare discussed in this section: BO relies on an approximate model of the optimizationlandscape, which is updated at each iteration, and is leveraged to choose the next setof parameters for which the FoM is evaluated.As highlighted in the previous section, the FoM can be obtained eitherexperimentally or numerically and is passed to the optimization routine. Each of theseevaluations can be time consuming, and in experimental setups, subject to noise. Forthese reasons, one cannot expect to resolve perfectly the true optimization landscape and reparation of ordered states in ultra-cold gases using Bayesian optimization f is specifiedin terms of a distribution, p ( f ) taken to be a Gaussian process [62] which is detailedin Appendix A.1. In essence, p ( f ) allows to favor smooth and regular functions todescribe the unknown control landscape F . As this distribution does not yet incorporateevaluations of the FoM it is therefore referred as the prior distribution .The next step in BO is to update the probability distribution p ( f ) based on thevalues of the FoM already collected. At an arbitrary step D of the optimization, D of such evaluation have been obtained and these set is denoted by a vector, y = [ y ( θ ) , . . . , y D ( θ D )]. The posterior distribution is the probability distribution p ( f | y ) defined as the distribution for f conditioned on all the evaluations y obtained sofar and is specified using Bayes rule, p ( f | y ) = p ( f ) p ( y | f ) p ( y ) , (2)where p ( y ) is the probability distribution for the set of observations y and p ( y | f ) isthe likelihood for the set of observations y to occur based on a given model f . Specificdetails regarding the evaluation of the predictive distribution p ( f ( θ ) | y ) for any controlparameters θ are described in Appendix A.2.Finally it remains to decide which set of control parameters θ D +1 to use in thenext step of the iterative optimization. One could naively choose it where the model f takes its maximal mean value. However, this model is only approximative and thuslikely to miss some interesting features of the optimization landscape, especially whenthe number of observations is low. Therefore it is also of interest to evaluate the FoMwhere uncertainty in the model is high. In BO, these two conflicting aspects, sometimesreferred respectively as exploitation and exploration, are captured by an acquisitionfunction α ( θ ) and the next set of parameters is chosen where this function reaches itsmaximum. For example it could be taken as α ( θ ) = µ f ( θ ) + kσ f ( θ ) , (3)with a positive scalar k and where both the mean value µ f ( θ ) and the standard deviation σ f ( θ ) are obtained according to the model predictive distribution given in Eq. (2). Themaximum of this function is reached when both its mean value and standard deviation,which quantifies the uncertainty in the modeling approach, are high (more or lessemphasis on one or the other can be modulated with k ). More details are given inAppendix A.3.To illustrate the points discussed above, a simple example of one dimensionaloptimization is considered in Fig. 2. The maximum value of an arbitrary function of asingle parameter θ ∈ [0 ,
4] is searched for using BO, implemented using the numericalpackage [63]. This function represents the true FoM which is shown in dashed red andthe noisy evaluations of the function, which serve as elements of y , are shown as reddots. Ten of such evaluations have been obtained and the posterior distribution p ( f | y )is plotted in the top panel of Fig. 2(a). The posterior distribution is defined by its mean reparation of ordered states in ultra-cold gases using Bayesian optimization θ after (a) 10, (b) 11 and (c) 20 iterations. The top panel depicts the underlyingfunction F to maximize (dashed red) along with its actual data points (red circles).The mean of the probabilistic model that approximates F is represented by a solidblue line and a 95% confidence interval (shaded region delimited by the solid graylines). The lower panel depicts the acquisition function, α ( θ ) (red solid line) given byEq. (3) with k = 4 (a-b) and 0 (c), whose maxima locations correspond to theselection of the next set of parameters in the optimization loop. For example, in (b)the new parameter was chosen at θ = 2 . α ( θ ) in (a).(shown in solid blue line) and its variance which are used to compute a 95% confidenceinterval (depicted as shaded blue across the mean). The width of this interval resultsfrom the finite number of evaluations and the noisy evaluations.As seen from Fig. 2(a), on comparing the model with the true FoM, one finds that itdoes not adequately reproduce the underlying control landscape due to the small numberof evaluations. This lack of knowledge of the true underlying landscape is captured bythe variance of the distribution resulting in a large confidence interval. The acquisitionfunction given in Eq. (A.7) is plotted in the lower panel of Fig. 2(a). In this case,the maximum of the acquisition function is at θ ≈ . θ . Fig. 2(b) incorporatesthis new evaluation (green square) and exhibits the updated model based on a total of11 iterations. This cycle of updating the model and suggesting the next parameter toevaluate is repeated another nine times and resulting in Fig. 2(c). In this figure, themodel is in close agreement with the true landscape and has successfully identified aparameter θ close to the true global maximum.Before moving to the next section showcasing the results of BO applied to specificproblems, we comment on the limitations and possible variations of this approach.Gaussian processes are a flexible tool to model well-behaved functions but are inadequateto model functions exhibiting discontinuities or varying degree of smoothness. However,this may be dealt with more complicated models as discussed in [64]. Similarly, reparation of ordered states in ultra-cold gases using Bayesian optimization p ( f | y ) in Eq. (2) isassumed to be Gaussian in a standard BO implementation. There may be scenarioswhere this is not the case, nonetheless, it is possible to incorporate non-Gaussian noiseinto the framework [61]. Another possible limitation is the computational complexityassociated to evaluating the posterior distribution given in Eq. (2) which grows as D with D as the number of accumulated evaluations. Alternatives to an exact treatmentof Gaussian processes are discussed in the outlook which may help to alleviate thiscomputational burden.
4. Preparation of ordered states in ultra-cold systems beyond adiabaticmethods
Controlled preparation of many-body states with spatial correlations resembling solid-state matter is of great interest to the condensed matter community. For example, theSF-MI (Superfluid-Mott insulator) transition realized with bosonic atoms trapped inoptical lattices [65, 66] simulates the transition of a conducting state to an insulatorstate in solid state physics. Similarly creating highly correlated ordered states withlong range interactions in a cold gas of Rydberg atoms [67, 68, 69] is akin to crystals insolids with long range Coulomb interactions between the electrons. In the Bose-Hubbardmodel, the only form of interaction is short-range while in the Rydberg system, it canbe relevant over a range longer than a typical lattice spacing [70, 71]. Both of theseexamples are well studied in the context of ultra-cold physics and will serve as perfecttestbeds for establishing the merits of BO techniques.In Sec. 4.1, using the Bose-Hubbard model as a toy example, BO’s ability to obtainan optimal protocol that drives the SF-MI transition is tested within the context ofquantum speed limit [49]. Furthermore, BO’s performance is benchmarked with otheroptimization methods in terms of its efficiency towards convergence. In Section 4.2, BOis applied to the optimization of laser pulse dynamics to create Rydberg crystals withlarge fraction of Rydberg excitations in different lattice geometries; one-dimensional(1D), two-dimensional (2D) as well as three-dimensional (3D).
The Bose-Hubbard model is widely studied in solid state physics and is simulated withbosonic atoms in an optical lattice. Although conceptually simple, this many-bodysystem cannot be mapped onto a single particle problem and contains the interestingphenomenon of transitioning from a SF state to a MI state [65]. Experimental realizationwas first achieved [66] by slowly varying the depth of the optical lattice potential. reparation of ordered states in ultra-cold gases using Bayesian optimization V ( t = ) V ( t = T ) V (t = T ) V (t = ) (a)(b) (c) Figure 3: Setup of Bose-Hubbard model: The goal of the optimization is to drive thesystem from an initial superfluid state illustrated in (a) to a Mott insulating phase asshown in (b) by dynamically changing the depth of the optical lattice V ( t ). (c) Energyspectrum for the Hamiltonian in Eq. (5) for different values of the control Γ. For the Bose-Hubbard model, a homogeneous system of N bosonic atomswith repulsive interactions in a lattice with L sites is considered. The Hamiltonian forthis model is given asˆ H ( t ) = − J ( t ) L (cid:88) i =1 (ˆ b i ˆ b † i +1 + ˆ b i +1 ˆ b † i ) + U ( t )2 L (cid:88) j =1 ˆ n j (ˆ n j − . (4)The first term in the Hamiltonian describes the tunnelling of bosons between neighboringpotential sites whose strength is given by J ( t ). The annihilation (creation) of a bosonat site i is defined by operators ˆ b i (ˆ b † i ) which follow the usual commutation relation,[ˆ b i , ˆ b † i ] = δ ij . The tunneling term tends to delocalize each atom over the lattice andthe repulsive interaction between two bosonic atoms at each site is quantified by U ( t ).Owing to the short range nature of the interactions compared to the lattice spacing,they are referred to as on-site interactions where ˆ n i = ˆ b † i ˆ b i is the number operator for agiven site i .As shown schematically in Fig. 3(a-b), by varying the potential depth of the opticallattice in time, both terms J ( t ) and U ( t ) are affected. J ( t ) changes depending on thetunnelling barrier between neighboring lattice sites, while the change in the on-siteinteraction is due to the variation in atomic wave function confinement. Introducinga single dimensionless quantity, Γ( t ) = U ( t ) / ( U ( t ) + J ( t )) such that Γ( t ) ∈ [0 , H ( t ) = − (1 − Γ( t )) L (cid:88) i =1 (ˆ b i ˆ b † i +1 + ˆ b i +1 ˆ b † i ) + Γ( t )2 L (cid:88) j =1 ˆ n j (ˆ n j − . (5)The Bose-Hubbard Hamiltonian has two distinct ground states depending on thestrength of the interactions U relative to the tunneling J . In the limit where Γ( t ) = 0, reparation of ordered states in ultra-cold gases using Bayesian optimization | ψ SF (cid:105) ∝ (cid:16) L (cid:88) i b † i (cid:17) N | (cid:105) . (6)Here | (cid:105) is the many-body vacuum state in the Fock basis representation. This SF stateis characterized by a non-zero variance of the number operator and a constant globalphase across the lattice.In the opposite limit with Γ( t ) = 1, the hopping between adjacent sites issuppressed and the ground state of the system consists of localized atomic wavefunctions that minimize the interaction energy. For a homogeneous system, one canhave commensurate filling of n atoms per lattice site and the many-body state is thena product of local Fock states in the atom number for each lattice site. For the groundstate where n = N/L , the Mott insulator is given as | ψ MI (cid:105) ∝ L (cid:89) i ( b † i ) N/L | (cid:105) . (7)Global phase coherence of the matter wave field is lost at the expense of perfect atomnumber correlations between lattice sites. In the next section we proceed to analyze theenergy spectrum of the Bose-Hubbard model. Interestingly, in the limits of Γ( t ) = 0or Γ( t ) = 1, the Hamiltonian in Eq. (5) is completely integrable. However, for theintermediate values of Γ, the system is non-integrable [72]. For a chain of L = 5 sites with periodic boundaryconditions and unit filling ( L = N = 5) of bosons, the many-body energy spectrum asa function of Γ is provided in Fig. 3(c). This is solved using exact diagonalization inthe Fock basis representation and is implemented with QuSpin [73]. The energies E n ofthe many-body eigenstates | n (cid:105) are plotted where n = 0 , , . . . enumerates each statein ascending order of energy at every value of Γ. The energies at Γ = 0 and Γ = 1can be expressed analytically [72]. The ground state is highlighted in red ranging fromsuperfluid ground state (Γ = 0) to Mott insulating ground state (Γ = 1) and takesspecific values, with energies E (Γ = 0) = − N and E (Γ = 1) = 0 respectively. Thereare multiple avoided crossings for all the other values of Γ. The aim is to find an optimal control function Γ opt ( t ) thatwill efficiently drive an initial superfluid state as defined in Eq. (6)) to a Mott-insulatingstate (see Eq. (7)) in time T . The general set of parameters θ described in Sec. 2 takespecific meaning in this context, as the time-dependent control Γ( t ) is parametrized with10 parameters representing the control values taken at equidistant time, Γ( t = i × T / i ∈ { , . . . , } . With the boundary conditions, Γ(0) = 0 and Γ( T ) = 1, the reparation of ordered states in ultra-cold gases using Bayesian optimization P n ( t ) (defined in the main text)of the instantaneous eigenstates for the linear (b) and optimized protocol (d) are alsoreported. For clarity they are grouped by: ground state only (blue), first to fifth(green) and higher (turquoise) eigenstates.values of the control for intermediate times are obtained by fitting a cubic spline.For the FoM, we will first consider the fidelity to the target state. It is defined as F MI = |(cid:104) ψ MI | ψ ( T ) (cid:105) | , where | ψ ( T ) (cid:105) is the realized state at time t = T . While thisfidelity is a natural choice for state-preparation task in numerical simulations, it is notreadily accessible in an experiment. A second choice of FoM, F exp , takes into accountexperimental constraints such as the effect of finite sampling. It is defined as F exp = (cid:104) V i (cid:105) i ,where V i is the variance in the occupation number estimated by averaging over for agiven lattice site i .In all cases the dynamics is solved numerically by exact diagonalization, and wedenote the eigenstates of H ( t ) from Eq. (5) as | n ( t ) (cid:105) . Thus, | n (0) = 0 (cid:105) and | n ( T ) = 0 (cid:105) are the ground state SF and MI states respectively. Using fidelity as the FoM, the duration of the protocolis fixed to the quantum speed limit as defined in [16, 49], T QSL = π/ ∆, where ∆ is theminimum energy gap between the ground and first excited state within the bounds ofΓ (see Fig. 3(c)). This limit theoretically bounds the minimum time to reach perfectfidelity by the system. Fig.4(a) shows the drive with a simple linear increase of Γ( t )and is compared to the optimal driving found with BO which is plotted in Fig.4(c).Fig.4 (b) and (d) exhibit the squared overlap of the state | ψ ( t ) (cid:105) with | n ( t ) (cid:105) , defined as P n ( t ) = |(cid:104) ψ ( t ) | n ( t ) (cid:105)| for the linear and optimized drives. As Γ( t ) is varied from 0 to 1,we retrieve the ground state fidelity at the end of the protocol with almost 90% in thecase of the optimized protocol while the linear ramp has a fidelity of less than 30%. For reparation of ordered states in ultra-cold gases using Bayesian optimization (a) (e)(b) (c)(f)(d) T = 0.5 T
QSL
T = T
QSL
T = 1.2 T
QSL
IterationsIterationsIterations I n fi d e li t y I n fi d e li t y Figure 5: The infidelity is reported as a function of the number of iterations fordifferent optimization routines (colors in legend). All the results obtained are based on30 runs of each setup, for which the median (dots) and the first to third quartileinterval (shaded region) are reported. Optimizations for different protocol times areshown: T = 0 . × T QSL (first column), T QSL (second) and 1 . × T QSL (last), where T QSL is defined in the text. Two figures of merit are used: the optimizer has access to eitherexact evaluations of the fidelity (first row) or to a noisy (hence experimentallyrealistic) FoM (second row) as described in the main text.the linear protocol to reach the same fidelity as the optimized one, it would take aboutan order of magnitude longer in time thus establishing the non-adiabatic character ofthe optimum protocol.To benchmark BO’s overall performance for the SF-MI task, results obtained usingBO are compared with other optimization routines (some of which were mentionedin Sec. 2) in Fig.5. These include both gradient as well as non-gradient methods.For comparison with gradient based method, the spontaneous perturbation stochasticapproximation (SPSA) algorithm was chosen [74] as it was shown to be competitivein the number of iterations for convergence and robust to noise [51]. For comparisonwith non-gradient based methods, the differential evolution(DE) and Nelder-Mead(NM)optimizers were implemented using Scipy [75]. Finally random search (RANDOM) wasalso included in the list of optimizers. It relies on randomly sampling a new set ofparameters at each iteration.In Fig.5, the three columns correspond to optimization results for different protocoldurations, T = 0 . T QSL , T
QSL , . T QSL , as these regimes can have different influenceson the complexity of the optimization task [52, 53]. Fig.5(a-c) show the infidelity,defined as (1 − F MI ), as a function of the number of iterations. However for Fig.5(d-f), reparation of ordered states in ultra-cold gases using Bayesian optimization
12a more experimentally realistic FoM is used, denoted by F exp . As mentioned before,this particular choice of FoM is motivated by experimental constraints such that theimpossibility to project onto the target state and finite size sampling effects. Theinfidelity in this case is directly given by F exp which is related to the variance in theoccupation number at each lattice site. F exp is minimized (vanishes) when a MI phasehas been reached, thus providing an appropriate guidance for the optimization. Ratherthan the true expected value of the variance, an estimation of its value based on a finitenumber of repetitions (here 1000) is used.As shown in Fig.5, over all these different configurations, BO (orange curve) exhibitsfaster convergence towards low infidelity. Both SPSA and DE also converge to goodparameters but at the expense of more iterations compared to BO. NM in most of thecases get quickly stuck in local minima with poor fidelities. When the time allowed forevolution is too short (see first column in Fig. 5), the best final infidelity reached isof the order of only 50%. All optimizers are quite competitive in this case. However,when the time allowed for the transition to occur is more or equal to T QSL (second andthird columns in Fig.5), an order of magnitude less iterations are required by BO toconverge to high fidelity control. Using the experimentally motivated FoM F exp , thefinal fidelities dropped by 5% for all the optimizers indicating additional difficulties toperform optimization in the noisy scenario. Interestingly SPSA performs better withthe noisy FoM Fig.5(d) while getting stuck when it has access to the true fidelity (c).This highlights the ability of SPSA to leverage noise to escape local minima.Based on these results, BO seems to provide a decisive advantage when it comes tonumber of iterations. After the first observations of interacting Rydberg gases [76, 77, 78, 79, 80, 81], morecontrol over the nature, form and structure of the quantum many-body state has beenachieved by shaping the excitation sequence [82] and the arrangement of the atoms inspecific patterns [83]. Early proposals [67, 68] for creating ordered phases in Rydberggases were adiabatic in nature, in the sense that the overall dynamics was slower than theminimum energy gap of the many-body system. Unfortunately, this energy gap decreasesexponentially with the number of Rydberg excitations preventing the preparation ofcrystals with large number of Rydberg excitations.Other proposals rely on admixing the ground state atoms with small amount ofRydberg state to create super-solid droplet crystals [84, 85]. This was experimentallyimplemented in [69], but only for a very small system where the energy gap was ratherlarge compared to the coherence time. Moreover Rydberg admixing of ground stateatoms requires considerable fine tuning [86] especially to avoid spontaneous scatteringand molecular resonances. More recently optimal control has been applied to optimizethe preparation of Rydberg crystals [87] for specific geometries (quasi one-dimensional)and very low Rydberg fraction. Having proved BO’s efficiency in the last example, there reparation of ordered states in ultra-cold gases using Bayesian optimization | g i i | g j i r ij | e i i | e j i C r ij ⌦ | e i| g i Figure 6: Schematic diagram of the setup: Ground state atoms denoted by | g (cid:105) aretrapped in a deep optical lattice. The ground state and Rydberg state ( | e (cid:105) ) of eachatom is optically coupled with an effective Rabi frequency Ω and detuning ∆. TwoRydberg atoms with a relative distance r ij interact via the van der Waals interaction.is a strong motivation to apply BO for preparing Rydberg ordered states by findingbetter protocols that give enhanced number of Rydberg excitations even in the presenceof experimentally realistic parameter noise and possible lattice filling imperfections. Fig. 6 illustrates the setup under consideration with Rb atoms in theirground state | g (cid:105) = 5 S coupled to a particular Rydberg state | e (cid:105) = 50 S with the helpof a laser, thus treating each atom effectively as a two-level system. The atoms can beeither trapped in a deep optical lattice [82, 88] or in an array of optical dipole traps[89] with uniform unit filling and lattice spacing l = 1 . µ m. The advantage of trappingatoms at suitable spacings is that unwanted molecular resonances and scattering can beavoided. Moreover, if needed, simultaneous trapping of ground and Rydberg state atomsis possible in magic wavelength lattices [90, 91, 92]. Two Rydberg atoms with positions r i and r j interact via the van der Waals interaction given by V ij = C ( n ) / | r i − r j | ,where C > n = 50 S takes the specific value, C = 1 . × − Hz m [93]. The Hamiltonian describing the full setup in the rotating wave approximation isgiven asˆ H = − (cid:126) ∆ (cid:88) i | e i (cid:105) (cid:104) e i | + (cid:126) Ω2 (cid:88) i ( | g i (cid:105) (cid:104) e i | + h . c . ) + (cid:88) i 10 GHz5 GHz E n E n (a)0 2 0 02 2[GHz] [GHz][GHz] [ G H z ][ G H z ] Figure 7: Panels (a-c) correspond to zero laser intensity, many-body energy spectrum(given in Eq. (9)) with Rydberg state, n = 50 (for which C = 1 . × − Hz m )shown for different lattice geometries (with lattice spacing l = 1 . µ m). Themany-body ground and excited states are shown explicitly as green and red linesrespectively. Panels (d-f) depict the many-body energy spectrum for the sameparameters as (a-c) but now with a coupling of Ω = 100 MHz.the system. Any motional dynamics is neglected as their typical timescales are muchlonger than the fast excitation dynamics that is considered here. Spontaneous decay ofthe Rydberg state including black-body radiation for the chosen state, 50S is about 65 µ s [14]. A gas of N atoms with n e Rydberg excitationscan be represented by the many-body eigenstate | n e , k (cid:105) where k labels the differentcombinatorial possibilities to distribute n e excitations over N atoms. The many-bodyenergy spectrum with Ω = 0 for different geometries is shown in Fig. 7(a-c). In the zerointensity field scenario, the energy of a many-body state is given as E n e ,k = − n e ∆ + (cid:15) n e ,k (9)where (cid:15) n e ,k are the interaction energies. The many-body ground state ( n e = 0) haszero energy and is highlighted in green in Fig. 7(a-c). States with the same numberof Rydberg excitations form a manifold with the same slope with respect to detuning.Thus the state that has the highest slope (shown in red) corresponds to all atoms in thelattice excited, n e = N . For non-zero laser intensity, anti-crossings occur between themany-body states as shown in Fig. 7(d-f). reparation of ordered states in ultra-cold gases using Bayesian optimization ν = n e /N (cid:28) , ∆) control landscape. The dynamics for Eq. (8) is numerically solved using alinear multi-step integration method implemented using the QuTiP package [94]. Forthis system, the relevant figures of merit are target state fidelity, spatial correlationfunction and number of Rydberg excitations, all of which can be potentially measuredin a typical Rydberg experiment with spatially resolved detection methods [95]. Thecontrol protocol is optimized for a specific number of Rydberg excitations n e whichis represented by the many-body eigenstates | n e , k (cid:105) . By denoting the fidelity as F n e ( t ) = (cid:80) k |(cid:104) Ψ( t ) | n e , k (cid:105)| , the FoM provided at each iteration to BO is given by F n e ( T ). The sum is taken over all the configurations containing n e Rydberg excitations.The total duration of the dynamics is fixed to T = 1 µ s which is well below theoverall lifetime of the many-body Rydberg system. The control field is parametrizedusing six parameters, [Ω( t ) , Ω( t ) , Ω( t ) , ∆( t ) , ∆( t ) , ∆( t )] which correspond to Rabifrequencies and detunings at three different times. For intermediate times, the valuesof (Ω( t ) , ∆( t )) are interpolated using quadratic polynomial. Additionally, boundaryconditions are imposed such that the intensity of the laser vanishes at t = 0 and t = T .This is done by multiplying Ω( t ) with a Tukey window function w ( t ) [96]. The boundsfor Ω ∈ [0 , . ∈ [ − . , The total number of atoms for 1D and 2D latticesis nine and the protocol is optimized for five Rydberg excitations. The 3D latticeconsists of eight atoms with the protocol optimized for four Rydberg excitations. Thefidelities F n e ( t ) over time for the optimal control protocols found by BO are depicted inFig. 8 (a-c) and Fig. 8 (d-f) show the optimized control parameters for different latticegeometries. In all cases, BO achieves high fidelity for the selected number of excitationsreaching convergence in the optimization within 10 iterations. In order to reproduce realexperiments and test BO’s resilience, 5% relative noise is introduced in Ω( t ) and ∆( t )respectively. The standard deviation for Ω( t ) and ∆( t ) is shown in shaded blue and red reparation of ordered states in ultra-cold gases using Bayesian optimization (a) (c) [ G H z ] 1D 2D 3D [ µ s] t [ µ s] t [ µ s] t site i P n ( T )site j s i t e k (c) (f ) s i t e i s i t e j s i t e k . . . . . . . . . . . . . . . . . . . . . P n ( T ) (a) site i . . . . . . site i P n ( T ) s i t e j (b) (e) site i s i t e j 11 22 33 0 . . . . . . P ( T ) P ( T ) F n e (b) (g) (h) (i) Figure 8: Results for optimized protocol dynamics using BO targeting a specificRydberg excitations with different lattice geometries: (a-c) show the plots of fidelities F n e ( t ) (defined in the text) for a given number of Rydberg excitations n e . It wasoptimized in (a) to maximize for F ( t ) in a 1D lattice, in (b) to maximize for F ( t ) ina 2D lattice and in (c) to maximize F ( t ) in a 3D lattice. Panels (d-f) are thecorresponding optimized protocols given in terms of Ω( t ) and ∆( t )). The three panels(g-i) show the final Rydberg excitation probability at any given lattice sitecorresponding to (a), (b) and (c) respectively.colors respectively across their mean values (which are depicted with solid lines). It ismost perceptible for Ω in Fig. 8 (e). These noises are usually negligible but could arisedue to experimental drifts or Doppler broadening.Inspecting the protocols optimized in Fig. 8, it seems like BO naturally selectsprotocols with the general trend where the detuning starts from large negative valuesand then transitions to positive values. However, the temporal profiles of the optimizedRabi frequencies are non-trivial especially for the 1D and 2D lattice models. Fig. 8(g-i) show the probability of Rydberg excitation P ( T ) for every site at the end of theoptimum protocol. It verifies the formation of Rydberg crystalline states as the atomicexcitations are maximally separated for any given lattice model. For 1D and 2D lattices,there is only one unique spatial configuration that corresponds to ordered state withfive Rydberg excitations. However in 3D lattice, there are two configurations that reparation of ordered states in ultra-cold gases using Bayesian optimization (a) 1D 3D2D [ G H z ] [ µ s] t [ µ s] t [ µ s] t (a) 1D 2D [ G H z ] (d) (b) (c) (e) (f ) F F F F n e F n e Figure 9: Results of optimizations in the presence of imperfect lattice fillings fordifferent lattice geometries. The solid (dashed) line refers to the protocol optimized fora perfect (imperfect) lattice filling which is denoted by OP1 (OP2). In (a-c), solid anddashed orange curves represent the fidelities for these specific protocols OP1 and OP2both applied with imperfect lattice filling. (d-f) show the details of these protocols interms of the Rabi frequency and detuning over time.equally contribute to crystal with four Rydberg excitations. In Fig. 8 (i), one of the twoconfigurations is rotated onto the other.Errors originating from inaccurate readouts and imperfect lattice fillings areprevalent in lot of the current cold atom experiments [13, 88]. Both these imperfectionsare included in the simulations by associating a probability for Rydberg detection foreach atom to be 0.9 when in Rydberg state. In one case, the optimum protocol isobtained for a perfect lattice and is labeled as OP1 while the other protocol is obtainedfor imperfect lattice filling which is labeled as OP2. The chosen FoM is the number ofRydberg excitations which was optimized for in Fig. 8.Fig. 9 (a-c) show the FoM calculated by applying both control dynamics, OP1(shown in solid orange) and OP2 (shown in dashed orange) on an imperfect lattice fill-ing averaged over 50 realizations. The profiles of the drivings that account for latticeimperfections are vastly different from the ones obtained in the idealized case. Takinginto account imperfections during the optimization results in better performances as canbe seen in Fig. 9 (a-c) where OP2 exhibits higher final fidelities compared to OP1 acrossall lattice geometries. This shows the importance of performing optimizations directlyonto the experiments rather than relying on idealized models or theoretical protocols forcharacterizing the noise sources and incorporating them in numerical simulations. Alsothe number of iterations required when noise is incorporated is five times higher whichreiterates the need for optimization routines to converge fast in order to be of practicaluse. reparation of ordered states in ultra-cold gases using Bayesian optimization ν = n e /N = 0 . 55 in T = 1 µ s. To create Rydbergcrystals with such large Rydberg fraction would take much longer in an adiabatic scheme.For example, the adiabatic preparation of Rydberg crystals in one dimension wererealized for ν = 0 . , T = 1 µ s in [67] and ν = 0 . , T = 4 µ s in [68]. The experiment[82] was done for ν = 0 . , T = 4 µ s while the optimal control found using Nelder-Mead[87] was obtained for a quasi-1D system with ν = 0 . , T = 4 µ s. In summary, BO isable to efficiently obtain high fidelity protocols despite the parameter noise and latticefilling imperfections that were included in the simulations. 5. Discussion and Outlook Quantum optimal control aims at designing fast high fidelity schemes to preparedesirable target states in a quantum system. An optimization routine that is ableto converge even when there is limited availability of data is advantageous for both,numerical simulations as well as for experiments. Along with this, the ability of theoptimization routine to handle noisy data is highly beneficial for experimental setups.In this context, BO seems to be a promising optimization tool and to make this case, itwas successfully applied to the creation of spatially ordered phases of ultra-cold gasestrapped in lattices.In the SF-MI transition task, the performance of BO in comparison to otheroptimization routines was found to be the most efficient with regards to the numberof iterations needed for convergence. This is primarily due to its ability to cleverlyselect the next set of parameters to probe. In the second example, notwithstanding thecomplexity in the many-body energy spectrum and the vanishing energy gap for largeRydberg systems, BO obtained control protocols to create Rydberg crystalline states inlattice models of arbitrary dimension with much larger Rydberg fraction than previousmethods. More interestingly, BO successfully navigated the optimization landscape evenwith the inclusion of noisy parameters and imperfect lattice fillings which is evidence ofthe advantage in accounting for noisy data in the modeling approach.In the examples studied here, the parametrization of the control field was takento be either spline or quadratic functions. Other types of smooth parametrizations, forexample with Fourier components or even with Gaussian processes, could also be tested.Furthermore in all the optimizations performed in this work no initial guess functionwas assumed for the control. Such a guess could be directly taken into account in theparametrization of the control function [55] and may result in optimized controls withhigher fidelity. Similarly the number of parameters used for the time-varying control fieldcould be increased to allow for finer control. However, in practice optimizing over higherdimensional parameter spaces may be challenging. For a given number of iterations, atrade-off between flexibility in the control field and realistic number of iterations has tobe found.Alternatively it would be interesting to study iterative refinements of the control reparation of ordered states in ultra-cold gases using Bayesian optimization 6. Acknowledgments The authors acknowledge H. Weimer for insightful and productive discussions. BothImperial and Stuttgart have received funding from the QuantERA ERANET Cofundin Quantum Technologies implemented within the European Unions Horizon 2020 Pro-gramme under the project Theory-Blind Quantum Control TheBlinQC and from EP-SRC under the grant EP/R044082/1. This work was supported through a studentshipin the Quantum Systems Engineering Skills and Training Hub at Imperial College Lon-don funded by the EPSRC(EP/P510257/1). [1] Bloch I, Dalibard J and Nascimb`ene S 2012 Nat. Phys. https://doi.org/10.1038/nphys2259 [2] Gross C and Bloch I 2017 Science https://science.sciencemag.org/content/357/6355/995 [3] Moses S, Covey J, Miecnikowski M, Jin D and Ye J 2017 Nat. Phys. https://doi.org/10.1038/nphys3985 [4] Bernien H, Schwartz S, Keesling A, Levine H, Omran A, Pichler H, Choi S, Zibrov A S, EndresM, Greiner M, Vuletic V and Lukin M D 2017 Nature https://doi.org/10.1038/nature24622 [5] Deutsch I H, Brennen G K and Jessen P S 2000 Fortschritte der Physik https://doi.org/10.1002/1521-3978(200009)48:9/11<925::AID-PROP925>3.0.CO;2-A [6] Pezz`e L, Smerzi A, Oberthaler M K, Schmied R and Treutlein P 2018 Rev. Mod. Phys. (3)035005 URL https://link.aps.org/doi/10.1103/RevModPhys.90.035005 [7] Hosten O, Engelsen N J, Krishnakumar R and Kasevich M A 2016 Nature https://doi.org/10.1038/nature16176 [8] Katori H 2011 Nat. Photonics https://doi.org/10.1038/nphoton.2011.45 reparation of ordered states in ultra-cold gases using Bayesian optimization [9] Gallagher T F 1988 Rep. Prog. Phys. https://doi.org/10.1088%2F0034-4885%2F51%2F2%2F001 [10] Robicheaux F and Hern´andez J V 2005 Phys. Rev. A (6) 063403 URL https://journals.aps.org/pra/abstract/10.1103/PhysRevA.72.063403 [11] Honer J, Weimer H, Pfau T and B¨uchler H P 2010 Phys. Rev. Lett. (16) 160404 URL https://link.aps.org/doi/10.1103/PhysRevLett.105.160404 [12] Zeiher J, Schauß P, Hild S, Macr`ı T, Bloch I and Gross C 2015 Phys. Rev. X (3) 031015 URL https://link.aps.org/doi/10.1103/PhysRevX.5.031015 [13] Omran A, Levine H, Keesling A, Semeghini G, Wang T T, Ebadi S, Bernien H, Zibov A S, PichlerH, Choi S, Cui J, Rossignolo M, Rembold P, Montangero S, Calarco T, Endres M, Greiner M,Vuleti´c V and Lukin M D 2019 Science https://science.sciencemag.org/content/365/6453/570 [14] Saffman M, Walker T G and Mølmer K 2010 Rev. Mod. Phys. (3) 2313–2363 URL https://link.aps.org/doi/10.1103/RevModPhys.82.2313 [15] Rosi S, Bernard A, Fabbri N, Fallani L, Fort C, Inguscio M, Calarco T and Montangero S 2013 Phys. Rev. A (2) 021601 URL https://link.aps.org/doi/10.1103/PhysRevA.88.021601 [16] van Frank S, Bonneau M, Schmiedmayer J, Hild S, Gross C, Cheneau M, Bloch I, Pichler T,Negretti A, Calarco T and Montangero S 2016 Sci. Rep. https://doi.org/10.1038/srep34187 [17] Wigley P B, Everitt P J, van den Hengel A, Bastian J W, Sooriyabandara M A, McDonald G D,Hardman K S, Quinlivan C D, Manju P, Kuhn C C N, Petersen I R, Luiten A N, Hope J J,Robins N P and Hush M R 2016 Sci. Rep. https://doi.org/10.1038/srep25890 [18] Tranter A D, Slatyer H J, Hush M R, Leung A C, Everett J L, Paul K V, Vernaz-Gris P,Lam P K, Buchler B C and Campbell G T 2018 Nat. Commun. https://doi.org/10.1038/s41467-018-06847-1 [19] Brochu E, Cora V M and De Freitas N 2010 arXiv:1012.2599 URL https://arxiv.org/abs/1012.2599 [20] Snoek J, Larochelle H and Adams R P 2012 In Proc. Advances in neural in-formation processing systems https://papers.nips.cc/paper/4522-practical-bayesian-optimization-of-machine-learning-algorithms.pdf [21] Shahriari B, Swersky K, Wang Z, Adams R P and De Freitas N 2015 Proc. IEEE https://doi.org/10.1109/JPROC.2015.2494218 [22] Cully A, Clune J, Tarapore D and Mouret J B 2015 Nature https://doi.org/10.1038/nature14422 [23] Calandra R, Seyfarth A, Peters J and Deisenroth M P 2016 Annals of Mathematics and ArtificialIntelligence https://doi.org/10.1007/s10472-015-9463-9 [24] Zhu D, Linke N M, Benedetti M, Landsman K A, Nguyen N H, Alderete C H, Perdomo-Ortiz A,Korda N, Garfoot A, Brecque C, Egan L, Perdomo O and Monroe C 2019 Sci. Adv. URL https://advances.sciencemag.org/content/5/10/eaaw9918 [25] Henson B M, Shin D K, Thomas K F, Ross J A, Hush M R, Hodgman S S and Truscott A G2018 Proc. Natl. Acad. Sci [26] Nakamura I, Kanemura A, Nakaso T, Yamamoto R and Fukuhara T 2019 Opt. Express [27] Weimer H, L¨ow R, Pfau T and B¨uchler H P 2008 Phys. Rev. Lett. (25) 250601 URL https://link.aps.org/doi/10.1103/PhysRevLett.101.250601 [28] Rabitz H, de Vivie-Riedle R, Motzkus M and Kompa K 2000 Science https://science.sciencemag.org/content/288/5467/824 [29] Judson R S and Rabitz H 1992 Phys. Rev. Lett. (10) 1500–1503 URL https://link.aps.org/doi/10.1103/PhysRevLett.68.1500 [30] Assion A, Baumert T, Bergt M, Brixner T, Kiefer B, Seyfried V, Strehle M and Gerber G 1998 reparation of ordered states in ultra-cold gases using Bayesian optimization Science https://doi.org/10.1126/science.282.5390.919 [31] Bartels R, Backus S, Zeek E, Misoguti L, Vdovin G, Christov I P, Murnane M M and KapteynH C 2000 Nature https://doi.org/10.1038/35018029 [32] Herek J L, Wohlleben W, Cogdell R J, Zeidler D and Motzkus M 2002 Nature https://doi.org/10.1038/417533a [33] N¨obauer T, Angerer A, Bartels B, Trupke M, Rotter S, Schmiedmayer J, Mintert F andMajer J 2015 Phys. Rev. Lett. (19) 190801 URL http://link.aps.org/doi/10.1103/PhysRevLett.115.190801 [34] Poggiali F, Cappellaro P and Fabbri N 2018 Phys. Rev. X (2) 021059 URL https://link.aps.org/doi/10.1103/PhysRevX.8.021059 [35] Kelly J, Barends R, Campbell B, Chen Y, Chen Z, Chiaro B, Dunsworth A, Fowler A G, HoiI C, Jeffrey E, Megrant A, Mutus J, Neill C, O’Malley P J J, Quintana C, Roushan P, SankD, Vainsencher A, Wenner J, White T C, Cleland A N and Martinis J M 2014 Phys. Rev. Lett. (24) 240504 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.240504 [36] Lu D, Li K, Li J, Katiyar H, Park A J, Feng G, Xin T, Li H, Long G, Brodutch A, BaughJ, Zeng B and Laflamme R 2017 npj Quantum Inf. 45 ISSN 2056-6387 URL https://doi.org/10.1038/s41534-017-0045-z [37] Li J, Yang X, Peng X and Sun C P 2017 Phys. Rev. Lett. (15) 150503 URL https://link.aps.org/doi/10.1103/PhysRevLett.118.150503 [38] Feng G, Cho F H, Katiyar H, Li J, Lu D, Baugh J and Laflamme R 2018 Phys. Rev. A (5)052341 URL https://link.aps.org/doi/10.1103/PhysRevA.98.052341 [39] Sp¨orl A, Schulte-Herbr¨uggen T, Glaser S J, Bergholm V, Storcz M J, Ferber J and Wilhelm F K2007 Phys. Rev. A (1) 012302 URL https://link.aps.org/doi/10.1103/PhysRevA.75.012302 [40] Bartels B and Mintert F 2013 Phys. Rev. A (5) 052315 URL https://journals.aps.org/pra/abstract/10.1103/PhysRevA.88.052315 [41] Montangero S, Calarco T and Fazio R 2007 Phys. Rev. Lett. (17) 170501 URL https://link.aps.org/doi/10.1103/PhysRevLett.99.170501 [42] Banchi L, Pancotti N and Bose S 2016 npj Quantum Inf. https://doi.org/10.1038/npjqi.2016.19 [43] Sklarz S E and Tannor D J 2002 Phys. Rev. A (5) 053619 URL https://link.aps.org/doi/10.1103/PhysRevA.66.053619 [44] Hohenester U, Rekdal P K, Borz`ı A and Schmiedmayer J 2007 Phys. Rev. A (2) 023602 URL https://link.aps.org/doi/10.1103/PhysRevA.75.023602 [45] Doria P, Calarco T and Montangero S 2011 Phys. Rev. Lett. (19) 190501 URL https://link.aps.org/doi/10.1103/PhysRevLett.106.190501 [46] Machnes S, Ass´emat E, Tannor D and Wilhelm F K 2018 Phys. Rev. Lett. (15) 150401 URL https://link.aps.org/doi/10.1103/PhysRevLett.120.150401 [47] Khaneja N, Reiss T, Kehlet C, Schulte-Herbr¨uggen T and Glaser S J 2005 J. Mag. Res. http://dx.doi.org/10.1016/j.jmr.2004.11.004 [48] Krotov V 1995 Global methods in optimal control theory vol 195 (CRC Press) URL https://doi.org/https://doi.org/10.1007/978-1-4612-0349-0_3 [49] Caneva T, Murphy M, Calarco T, Fazio R, Montangero S, Giovannetti V and Santoro G E 2009 Phys. Rev. Lett. (24) 240501 URL https://link.aps.org/doi/10.1103/PhysRevLett.103.240501 [50] Dive B, Pitchford A, Mintert F and Burgarth D 2018 Quantum 80 ISSN 2521-327X URL https://doi.org/10.22331/q-2018-08-08-80 [51] Ferrie C and Moussa O 2015 Phys. Rev. A (5) 052306 URL https://link.aps.org/doi/10.1103/PhysRevA.91.052306 [52] Zahedinejad E, Schirmer S and Sanders B C 2014 Phys. Rev. A (3) 032310 URL https://link.aps.org/doi/10.1103/PhysRevA.90.032310 reparation of ordered states in ultra-cold gases using Bayesian optimization [53] Bukov M, Day A G R, Sels D, Weinberg P, Polkovnikov A and Mehta P 2018 Phys. Rev. X (3)031086 URL https://link.aps.org/doi/10.1103/PhysRevX.8.031086 [54] McClean J R, Boixo S, Smelyanskiy V N, Babbush R and Neven H 2018 Nat. Commun. https://doi.org/10.1038/s41467-018-07090-4 [55] Caneva T, Calarco T and Montangero S 2011 Phys. Rev. A (2) 022326 URL https://link.aps.org/doi/10.1103/PhysRevA.84.022326 [56] Rach N, M¨uller M M, Calarco T and Montangero S 2015 Phys. Rev. A (6) 062343 URL https://link.aps.org/doi/10.1103/PhysRevA.92.062343 [57] Egger D J and Wilhelm F K 2014 Phys. Rev. Lett. (24) 240503 URL https://link.aps.org/doi/10.1103/PhysRevLett.112.240503 [58] Palittapongarnpim P, Wittek P, Zahedinejad E, Vedaie S and Sanders B C 2017 Neurocomputing [59] F¨osel T, Tighineanu P, Weiss T and Marquardt F 2018 Phys. Rev. X (3) 031084 URL https://link.aps.org/doi/10.1103/PhysRevX.8.031084 [60] Niu M Y, Boixo S, Smelyanskiy V N and Neven H 2019 npj Quantum Information 33 ISSN2056-6387 URL https://doi.org/10.1038/s41534-019-0141-3 [61] Sauvage F and Mintert F 2019 arXiv:1909.01229 URL https://arxiv.org/abs/1909.01229 [62] Williams C K and Rasmussen C E 2006 Gaussian processes for machine learning vol 2 (MITpress Cambridge, MA) URL [63] 2016 Gpyopt: A bayesian optimization framework in python URL http://github.com/SheffieldML/GPyOpt [64] Neal R M 1998 Bayesian statistics 475 URL [65] Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 Phys. Rev. B (1) 546–570URL https://link.aps.org/doi/10.1103/PhysRevB.40.546 [66] Greiner M, Mandel O, Esslinger T, H¨ansch T W and Bloch I 2002 Nature https://doi.org/10.1038/415039a [67] Pohl T, Demler E and Lukin M D 2010 Phys. Rev. Lett. (4) 043002 URL https://link.aps.org/doi/10.1103/PhysRevLett.104.043002 [68] van Bijnen R M W, Smit S, van Leeuwen K A H, Vredenbregt E J D and KokkelmansS J J M F 2011 J. Phys. B:At. Mol. Opt. Phys. https://doi.org/10.1088%2F0953-4075%2F44%2F18%2F184008 [69] Zeiher J, van Bijnen R, Schauß P, Hild S, Choi J y, Pohl T, Bloch I and Gross C 2016 Nat. Phys. https://doi.org/10.1038/nphys3835 [70] Bohlouli-Zanjani P, Petrus J A and Martin J D D 2007 Phys. Rev. Lett. (20) 203005 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.203005 [71] B´eguin L, Vernier A, Chicireanu R, Lahaye T and Browaeys A 2013 Phys. Rev. Lett. (26)263201 URL https://link.aps.org/doi/10.1103/PhysRevLett.110.263201 [72] Kolovsky A R and Buchleitner A 2004 Europhys. Lett. https://doi.org/10.1209%2Fepl%2Fi2004-10265-7 [73] Weinberg P and Bukov M 2017 SciPost Phys. (1) 003 URL https://scipost.org/10.21468/SciPostPhys.2.1.003 [74] Spall J C 1998 IEEE Transactions on Aerospace and Electronic Systems https://ieeexplore.ieee.org/document/705889 [75] Virtanen P and others SciPy 1 0 2019 arXiv:1907.10121 URL https://arxiv.org/abs/1907.10121 [76] Tong D, Farooqi S M, Stanojevic J, Krishnan S, Zhang Y P, Cˆot´e R, Eyler E E and Gould P L 2004 Phys. Rev. Lett. (6) 063001 URL https://link.aps.org/doi/10.1103/PhysRevLett.93.063001 [77] Vogt T, Viteau M, Zhao J, Chotia A, Comparat D and Pillet P 2006 Phys. Rev. Lett. (8) reparation of ordered states in ultra-cold gases using Bayesian optimization https://link.aps.org/doi/10.1103/PhysRevLett.97.083003 [78] Reetz-Lamour M, Deiglmayr J, Amthor T and Weidemller M 2008 New J. Phys. https://doi.org/10.1088%2F1367-2630%2F10%2F4%2F045026 [79] Mohapatra A K, Jackson T R and Adams C S 2007 Phys. Rev. Lett. (11) 113003 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.113003 [80] Raitzsch U, Bendkowsky V, Heidemann R, Butscher B, L¨ow R and Pfau T 2008 Phys. Rev. Lett. (1) 013002 URL https://link.aps.org/doi/10.1103/PhysRevLett.100.013002 [81] L¨ow R, Weimer H, Krohn U, Heidemann R, Bendkowsky V, Butscher B, B¨uchler H P and Pfau T2009 Phys. Rev. A (3) 033422 URL https://link.aps.org/doi/10.1103/PhysRevA.80.033422 [82] Schauß P, Zeiher J, Fukuhara T, Hild S, Cheneau M, Macr`ı T, Pohl T, Bloch I and Gross C 2015 Science https://science.sciencemag.org/content/347/6229/1455 [83] Labuhn H, Barredo D, Ravets S, de L´es´eleuc S, Macr`ı T, Lahaye T and Browaeys A 2016 Nature https://doi.org/10.1038/nature18274 [84] Henkel N, Nath R and Pohl T 2010 Phys. Rev. Lett. (19) 195302 URL https://link.aps.org/doi/10.1103/PhysRevLett.104.195302 [85] Cinti F, Jain P, Boninsegni M, Micheli A, Zoller P and Pupillo G 2010 Phys. Rev. Lett. (13)135301 URL https://link.aps.org/doi/10.1103/PhysRevLett.105.135301 [86] Balewski J B, Krupp A T, Gaj A, Hofferberth S, Lw R and Pfau T 2014 New J. Phys. https://doi.org/10.1088%2F1367-2630%2F16%2F6%2F063012 [87] Cui J, van Bijnen R, Pohl T, Montangero S and Calarco T 2017 Quantum Science Technol. https://doi.org/10.1088%2F2058-9565%2Faa7daf [88] Schauß P, Cheneau M, Endres M, Fukuhara T, Hild S, Omran A, Pohl T, Gross C, KuhrS and Bloch I 2012 Nature https://doi.org/10.1038/nature11596 [89] Nogrette F, Labuhn H, Ravets S, Barredo D, B´eguin L, Vernier A, Lahaye T and Browaeys A 2014 Phys. Rev. X (2) 021034 URL https://link.aps.org/doi/10.1103/PhysRevX.4.021034 [90] Zhang S, Robicheaux F and Saffman M 2011 Phys. Rev. A (4) 043408 URL https://link.aps.org/doi/10.1103/PhysRevA.84.043408 [91] Topcu T and Derevianko A 2013 Phys. Rev. A (5) 053406 URL https://link.aps.org/doi/10.1103/PhysRevA.88.053406 [92] Mukherjee R, Millen J, Nath R, Jones M P A and Pohl T 2011 J. Phys. B:At. Mol. Opt. Phys. https://doi.org/10.1088%2F0953-4075%2F44%2F18%2F184010 [93] Reinhard A, Liebisch T C, Knuffman B and Raithel G 2007 Phys. Rev. A (3) 032712 URL https://link.aps.org/doi/10.1103/PhysRevA.75.032712 [94] Johansson J, Nation P and Nori F 2013 Comput. Phys. Commun. https://doi.org/10.1016/j.cpc.2012.11.019 [95] Ott H 2016 Rep. Prog. Phys. https://doi.org/10.1088%2F0034-4885%2F79%2F5%2F054401 [96] Harris F J 1978 Proceedings of the IEEE https://ieeexplore.ieee.org/document/14551063 [97] Ghahramani Z 2015 Nature https://doi.org/10.1038/nature14541 [98] Hensman J, Fusi N and Lawrence N D 2013 In Proc. Twenty-Ninth Conference on Uncertaintyin Artificial Intelligence [99] Snoek J, Rippel O, Swersky K, Kiros R, Satish N, Sundaram N, Patwary M, Prabhat M andAdams R 2015 In Proc. 32nd International Conference on Machine Learning http://proceedings.mlr.press/v37/snoek15.html [100] Swersky K, Snoek J and Adams R P 2013 In Proc. Advances in neural in- reparation of ordered states in ultra-cold gases using Bayesian optimization formation processing systems https://papers.nips.cc/paper/5086-multi-task-bayesian-optimization.pdf [101] 2012 Gpy: A gaussian process framework in python URL http://github.com/SheffieldML/GPy Appendix A. Bayesian optimization: techniques The problem of quantum optimal control defined in Section.2 consists on finding optimal(or good enough) control parameters θ such that the figure of merit F ( θ ), is maximized.BO relies on an approximative model f , called the surrogate model , of the figure of meritwhich is used to guide the optimization. Appendix A.1. Building the surrogate model: Prior distribution In the context of BO the surrogate model is taken to be a random function , alsoknown as stochastic process . Random functions extend the notion of a finite set ofrandom variables to an infinite one. Thus, considering f to be a random functionmeans that for any of the infinitely many input parameters θ m , the value f ( θ m )taken by the function is itself a random variable. More precisely, f is taken to bea Gaussian Process (GP) which is a specific type of random functions [62]. A GPis the natural extension of a Multi-Variate (MV) Gaussian distribution: while a MVGaussian distribution is entirely specified by a mean vector and a covariance matrix ,a GP is given in terms of a mean function m and a covariance function k . Thesetwo functions define respectively the mean (cid:104) f ( θ m ) (cid:105) = m ( θ m ), and the covariances (cid:104) f ( θ m ) f ( θ n ) (cid:105) − (cid:104) f ( θ m ) (cid:105)(cid:104) f ( θ n ) (cid:105) = k ( θ m , θ n ) of the values taken by the model at anyinput parameters θ m and θ n . Considering the model f to be a GP with mean function m and covariance function k is denoted as p ( f ) = GP ( m, k ) and is defined such that any finite collection of random variables f ( θ i ), of arbitrary length N , follows a multivariateGaussian distribution: p (cid:16) f ( θ )... f ( θ N ) (cid:17) = N (cid:16) m = m ( θ )... m ( θ N ) , K = k ( θ , θ ) .. ( θ , θ N )... ... k ( θ N , θ ) .. ( θ N , θ N ) (cid:17) (A.1)where N denotes a Gaussian distribution here with a mean vector m of length N anda covariance matrix K of dimension N × N . The distribution over functions p ( f ) orequivalently the joint distributions over finite sets of function values as given in Eq. A.1,are referred as the prior distribution as they do not incorporate any observations yet.A specific choice of the mean and covariance functions defines the global propertiesof the model. The mean is often taken to be vanishing, m ( θ ) ≡ f are entirely delegated to the choice of the covariancefunction k . A common choice, followed here, for the covariance function is the Mat´ern5/2 function imposing the model to be twice differentiable [62]. This Mat´ern 5/2reparation of ordered states in ultra-cold gases using Bayesian optimization (a) (b) Figure A1: (a) The figure shows three functions sampled from a GP with a Mat´erncovariance (Eq.A.2) for different values of the length-scale l and the varianceparameter σ . (b) Predictive distribution obtained according to Eq.A.5 for a 1-dparameter θ ∈ [0 , 1] with 10 noisy observations (red dots). The distribution isdisplayed in terms of its mean function µ f ( θ ) (blue thick line) and a 95% confidenceinterval (grey light envelope) with boundaries µ f ( θ ) ± . × σ f ( θ ). Additionallyexamples of the full density for two specific inputs θ = 0 . 25 and θ = 2 . k / ( θ m , θ n ) = k / ( x = | θ m − θ n | ) = σ (1 + xl + x l ) e − x/l , (A.2)where the variance σ and the length-scale l , are free hyper-parameters. Essentially thevariance σ = k / ( x = 0) = (cid:104) f ( θ m ) (cid:105) specifies to which extend any variable f ( θ m )is expected to deviate from its mean value (here 0 with our choice of mean function m ), while the length-scale l , appearing in the exponentially decaying term, scales thedistance x between the parameters θ m and θ n which is indicative of the extent ofthe correlations between different values of the function. To illustrate the impact ofthese hyper-parameters, Fig.A1 (a) shows 3 functions with have been sampled withdifferent hyper-parameters values (given in legend). Each sample is obtained accordingto Eq. (A.1) for a vector of 1-d input parameters [ θ , . . . , θ N ] with N large enough suchthat each sample with finite size N effectively looks like a function.In general, it is not possible to have a precise idea of the values of these hyper-parameters beforehand but they can be fitted to the observations at any stage of theoptimization. This fitting is often done by minimizing the log marginal-likelihood [62]and is implemented in any GP library such as [101]. In summary the prior distributionis entirely defined by the choice of a GP, a mean function m , a covariance function k and the hyper-parameters σ and l . In particular, under the choice of a mean anda kernel function presented here, a priori each value f ( θ m ) taken by the model has adistribution: : p ( f ( θ m )) = N (0 , σ ) . (A.3) reparation of ordered states in ultra-cold gases using Bayesian optimization Appendix A.2. Updating the surrogate model: Predictive distribution Starting with the prior distribution described in the previous section, it remains toincorporate observations to the model to obtain the posterior distribution. In thepresent context this posterior is also called the predictive distribution as it allowsto make prediction for unseen parameters. Observations are denoted by a vector y = [ y ( θ ) , . . . , y ( θ M )] where each element y ( θ m ) corresponds to an evaluation of thefigure of merit at iteration m for parameters θ m , over a total of M iterations.In the case of noiseless evaluations of the figure of merit one directly records thevalues taken by the model for the different parameters tried. Writing this vector ofvalues f = [ f ( θ ) , . . . , f ( θ M )], it follows the equality f = y . Applying the definition ofa GP in Eq. (A.1) to the vector of random variables [ f ( θ ) , ..., f ( θ M ) , f ( θ ∗ )] and usingthe conditioning properties of Gaussian distributions [62], the predictive distribution isgiven by p ( f ( θ ∗ ) | f ) = N ( µ f ( θ ∗ ) = k T K − f , σ f ( θ ∗ ) = σ − k T K − k ) , (A.4)where the column vector k has entries k m = k ( θ m , θ ∗ ), and elements of the covariancematrix K are given by K m,n = k ( θ m , θ n ). The symbols µ f ( θ ∗ ) and σ f ( θ ∗ ) denote themean and standard deviation of this predictive distribution. They both are function ofthe parameter θ ∗ as each element of k depends on it. Compared to the prior distribution p ( f ( θ ∗ )) in Eq. (A.3), the mean of the predictive distribution has been shifted from 0 to k T K − f and its variance has decreased from σ by a positive quantity k T K − k , resultingfrom the incorporation of the observations. The most computational demanding partof evaluating this mean and variance comes from the inversion of the matrix K whichhas complexity O ( M ).In an experimental scenario, the set of measurements y does not directly revealthe values taken by the model, still model and observations can be related by positinga noise model. This noise, originating from both experimental imperfections and alsoquantum fluctuations, can be approximated as an additive constant Gaussian noise suchthat y ( θ j ) = f ( θ j )+ ε j , where the noise terms ε j are assumed to be independently drawnfrom the same distribution p ( ε j ) = N (0 , σ N ). The extra parameter σ N , capturing theamount of noise, can be fixed or could also be fitted to the data in the same way as thehyper-parameters σ and l . Under the independence assumption for the noise terms, thelikelihood of recording the full data set of observations y for a given set of values f takenby the model can be written as p ( y | f ) = N ( f , σ N I ). With this simple phenomenologicalmodel of the noise the sought-after predictive distribution can be derived: p ( f ( θ ∗ ) | y ) = (cid:90) df p ( f ( θ ∗ ) | f ) p ( f | y ) = (cid:90) df p ( f ( θ ∗ ) | f ) p ( y | f ) p ( f ) /p ( y )= N (cid:16) µ f ( θ ∗ ) = k T ( K + σ N I ) − f , σ f ( θ ∗ ) = σ − k T ( K + σ N I ) − k (cid:17) , (A.5)where the first line shows how it can be decomposed, using Bayes rule andmarginalization over the vector of values f , as a product of the prior distribution p ( f ) reparation of ordered states in ultra-cold gases using Bayesian optimization p ( y | f )previously defined and the probability of obtaining a given set of observations p ( y ). Thenext line provides the analytical result of this integral [62]. This final result is quitesimilar to Eq. (A.4) except that the inverse of the covariance matrix K − has now beenreplaced by ( K + σ N I ) − thus incorporating noise in the observations.A practical example of evaluating this predictive distribution for a one-dimensionalparameter θ ∈ [0 , 4] based on a set of M = 10 observations is provided in Fig. A1(b).The mean function µ f and a confidence interval with bounds µ f ± . σ f have beenrepresented with, in addition, the explicit densities p ( f ( θ ) | y ) for two specific values ofthe parameter. Far away from the observations, for example for θ ≈ . 75, the confidenceinterval is wider indicating uncertainty in the predictive distribution, while closer toobservations, for example for θ ≈ . 5, it shrinks but does not vanish as the modelidentified noise in the observations (more precisely the value of the hyper-parameter σ N representing the strength of noise after fitting is non 0). Appendix A.3. Decision rule: Acquisition function The last step in the BO framework is the choice of the next parameter θ M +1 toobserve based on the surrogate model. This new parameter could be picked suchthat it maximizes the mean of the predictive distribution given in Eq. (A.5): θ M +1 =arg max θ µ f ( θ ). However, with only a limited amount of observations, it is likely thatthe model does not capture the full range of variations in the figure of merit and thisapproach is likely to end up in a local minima. This can be seen for example in Fig.2(a)where the model misses one of the peak. Alternatively an explorative strategy wouldconsist on taking measurements where the uncertainty in the model is maximal. Thiscan be done by choosing θ M +1 = arg max θ σ f ( θ ), where the standard deviation σ f ( θ )quantifies the uncertainty in the model. Balancing these two objectives, which arefinding the location of the maximum and also exploring region of the parameter spacewhere not enough observations have been made, is often referred as the exploration-exploitation trade-off. In BO this trade-off is dealt with by introducing an acquisitionfunction α ( θ ) aiming at capturing these two aspects; and the choice of the nextparameter is taken such that: θ M +1 = arg max θ α ( θ ) . (A.6)Among the most popular acquisition functions we recall the definition of the UpperConfidence Bound (UCB) function already provided in the main text: α UCB ( θ ) = µ f ( θ ) + kσ f ( θ ) , (A.7)where the positive scalar k balances the bias toward exploration (for high value of k)or exploitation (small values of k). This acquisition function was used to produce thegraphs in the bottom panels of Fig.2. In Fig.2(a-b) a value of k = 4 was used, while in reparation of ordered states in ultra-cold gases using Bayesian optimization k = 0 in order to showthe optimal parameter found according to the model.Another acquisition function often used is the Expected Improvement (EI) function,where the improvement, defined with regards to the best observation of the figure ofmerit obtained so far y max , is averaged under the predictive distribution: α EI ( θ ) = (cid:90) max (0 , f ( θ ) − y max ) p ( f ( θ ) | y ) df ( θ ) , (A.8)which has analytic solution [20]. The choice of the acquisition function can have insome cases significant impact on the final results and both the UCB and EI functionswere considered in this study. In practice finding the maximum of these acquisitionfunctions is often performed by gradient ascent (possibly repeated over several randomstarting parameters). For both the UCB and EI, the gradients of the acquisitionfunction with regards to the parameters have analytical expressions and can be efficientlyevaluated. These two acquisition functions and the routines for finding their maximaare implemented in most of the BO libraries such as [63]. Appendix B. Bayesian optimization: implementation details The specific details of the BO runs used to obtain the results presented in Sec.4 aregiven here. Appendix B.1. SF - MI transition For the SF-MI transition example, all the optimizations shown in Fig. 5 were performedwith the same configuration. In each case, N init = 100 observations were initially takenfor randomly chosen parameters followed by M = 1750 iterations of BO. To decrease thecomputational burden associated to this large number of iterations we found it efficientnot to update the hyper-parameters (noise and parameters of the kernel functions) ateach iteration but rather every 10 iterations. This resulted in a saving of almost 40%in computational time without affecting the quality of the convergence results. For thechoice of the next parameter to probe the UCB acquisition function in Eq.A.7 was usedwith a value of k linearly decreasing from k (1) = 5 at the first iteration to k (1750) = 0 atthe end of the optimization. For this problem this acquisition function performed betterthan the EI one. Finally other types of kernels [62] were compared and we found thatthe Mat´ern 5/2 kernel performed almost always better than the Radial Basis Function(RBF) kernel which is also a popular choice [17]. Appendix B.2. Rydberg crystalline states