Preparation of ordered states in ultra-cold gases using Bayesian optimization

2020·Arxiv

Abstract

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 use them for the study of complex quantum many-body phenomena. However in many scenarios, the prerequisite condition of faithfully preparing a desired quantum state despite decoherence and system imperfections is not always adequately met. To pave the way to a specific target state, we implement quantum optimal control based on Bayesian optimization. The probabilistic modeling and broad exploration aspects of Bayesian optimization are particularly suitable for quantum experiments where data acquisition can be expensive. Using numerical simulations for the superfluid to Mottinsulator transition for bosons in a lattice as well as for the formation of Rydberg crystals as explicit examples, we demonstrate that Bayesian optimization is capable of finding better control solutions with regards to finite and noisy data compared to existing methods of optimal control.

1. Introduction

In this paper, the focus is on creating spatially ordered states in two different ultra-cold systems, atoms in optical lattice and highly excited Rydberg atoms, as testbeds for Bayesian optimization. Optical lattices [1, 2] provide a natural and versatile platform to study strongly correlated condensed matter systems [3, 4], with implementations in quantum information processing [5] and quantum metrology [6, 7, 8]. Likewise the extraordinary properties of Rydberg atoms [9] in combination with ultra-cold gases establish 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 of the largest non-classical states [13] and have emerged as a serious competitor for the realization of quantum computer [14]. While adiabatic preparation still remains an intuitive and straightforward way to create these interesting states, it requires long timescales especially for large systems with diminishing energy gaps. This makes it prone to accumulation of errors due to prolonged interaction with the environment. Furthermore, there is no guarantee that a particular target state of interest is adiabatically connected to the initial state. Thus in recent years, there has been a growing interest in applying optimal control to ultra-cold systems [13, 15, 16, 17, 18].

At the core of any quantum optimal control framework, lie optimization algorithms that aim at maximizing the figure of merit of the experiment with respect to input control fields. A figure of merit (FoM) measures the success of the experiment in achieving a specific task. This could be the preparation of a highly desirable quantum state, implementation of a gate or even achieving a certain phase transition. The input control fields of a typical quantum experiment can be tunable variables such as external electric or magnetic fields, amplitude and phase of lasers or microwave pulses. The explicit dependence of the FoM on these input control fields defines the control (optimization) landscape. To obtain an optimal solution, an optimization scheme generally requires to probe several points in the control landscape. Undoubtedly the more points the optimization algorithm can probe the better is the optimal solution. However, if the FoM is obtained from a quantum experiment, probing several points in the control landscape would imply repeating the experiment several times. Depending on the complexity of the quantum task at hand and the optimization routine, the experiment may have to be performed for unusually high number of times, which is highly inefficient. This constraint applies equally well to numerical calculations where the Hilbert space grows exponentially with system size. Most conventional methods of optimization find fast protocols but not necessarily efficiently with regards to the number of data points needed from the control landscape. Furthermore, there is always some uncertainty in evaluating the FoM due to experimental noise and imperfections which makes the optimization process even more challenging.

In this regard, Bayesian Optimization (BO) [19, 20, 21] is an attractive option and has been successfully used in robotics [22, 23], machine-learning [20] and has recently made its way to quantum optimal control problems [17, 24, 25, 26]. Section 2 discusses how BO fits into the general framework of quantum optimal control followed by a discussion outlining the essential steps of BO in Section 3. In Section 4, BO is implemented for two different examples. In the first example, BO’s efficiency with regards to the minimum number of evaluations required to drive a superfluid (SF) to Mott insulator (MI) state transition is established when compared to several other optimization routines. In the second example, BO is applied to a strongly interacting gas of trapped Rydberg atoms. Although initial studies of Rydberg gases revealed the possibility of creating ordered structures via a phase transition [27], the presence of a quantum critical point makes the realization of crystalline ordering with large number of Rydberg excitations still challenging. Section 5 summarizes these results and provides an outlook for future work.

Figure 1: Quantum optimal control scheme: ) is the FoM to maximize which depends on input parameters that describes the control fields of the quantum system. ) is evaluated either numerically or experimentally. The optimization task is initiated by evaluating this FoM for a set of random control parameters every stage, the Bayesian optimizer suggests the next set of parameters for the quantum system based on previous evaluations of the FoM. This process is repeated over M iterations.

2. Quantum Optimal Control formalism

Historically, quantum optimal control was first applied to manipulate chemical reactions by selectively breaking bonds in molecules using lasers [28, 29, 30, 31, 32]; more recent examples 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 involves maximizing the FoM, ) with respect to the set of parameters until an optimum solution, is reached, as expressed below,

] is a vector of size P that represents the parametrization of the input control fields. For time-dependent control, this parametrization can be done using special 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 ) is obtained either from a real experiment or through numerical calculations. The hope is to converge close enough to well before exhausting the experimental/numerical resources. Although experimental constraints are not explicitly included in Eq. (1), they can be incorporated either in the FoM directly or in the choice of the parametrization.

A variety of optimization algorithms are available and are broadly classified based on the ability to evaluate the gradient of ) with respect to Gradient based algorithms 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 hybrid quantum-classical approach [37, 38] can be applied to experiments. These gradient based methods are considered as local optimization methods whose convergence is fast as long as the optimization landscape is well-behaved. However, these methods get compromised by the presence of any local minima and plateaus in the optimization landscape [52, 53, 54]. Alternatively, gradient-free algorithms are able to explore the optimization landscape more globally than gradient based methods making them less vulnerable to local minima. Among them, Nelder-Mead has been extensively used in the context of the chopped random basis (CRAB) method [45, 55, 56] as well as independently [34, 35, 57]. It has the advantage of simplicity but can still get trapped in 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]. These methods usually require large number of iterations to find the optimal solution.

Bayesian Optimization (BO) is a non-gradient based method offering an appealing alternative as it cleverly selects the next set of parameters to evaluate at each step of the optimization. This can lead to significant reduction in the number of iterations needed for convergence while performing global optimization. BO has the additional advantage of integrating probabilistic elements of data acquisition which can be exploited to further increase its efficiency [61]. These characteristics make it especially well suited to perform optimal control on quantum experiments but also when numerical simulations of the system are time-consuming.

3. Bayesian Optimization

The essential steps of BO used for quantum optimal control are shown in Fig. 1 and are discussed in this section: BO relies on an approximate model of the optimization landscape, which is updated at each iteration, and is leveraged to choose the next set of parameters for which the FoM is evaluated.

As highlighted in the previous section, the FoM can be obtained either experimentally or numerically and is passed to the optimization routine. Each of these evaluations can be time consuming, and in experimental setups, subject to noise. For these reasons, one cannot expect to resolve perfectly the true optimization landscape and probabilistic modeling becomes convenient. A suitable probabilistic model f is specified in terms of a distribution, p(f) taken to be a Gaussian process [62] which is detailed in Appendix A.1. In essence, p(f) allows to favor smooth and regular functions to describe the unknown control landscape F. As this distribution does not yet incorporate evaluations 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 the values 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, is the probability distribution p(f|y) defined as the distribution for f conditioned on all the evaluations y obtained so far and is specified using Bayes rule,

where p(y) is the probability distribution for the set of observations y and p(y|f) is the likelihood for the set of observations y to occur based on a given model f. Specific details regarding the evaluation of the predictive distribution ) for any control parameters are described in Appendix A.2.

Finally it remains to decide which set of control parameters to use in the next 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 thus likely to miss some interesting features of the optimization landscape, especially when the number of observations is low. Therefore it is also of interest to evaluate the FoM where uncertainty in the model is high. In BO, these two conflicting aspects, sometimes referred respectively as exploitation and exploration, are captured by an acquisition ) and the next set of parameters is chosen where this function reaches its maximum. For example it could be taken as

with a positive scalar k and where both the mean value ) and the standard deviation ) are obtained according to the model predictive distribution given in Eq. (2). The maximum 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 less emphasis on one or the other can be modulated with k). More details are given in Appendix A.3.

To illustrate the points discussed above, a simple example of one dimensional optimization is considered in Fig. 2. The maximum value of an arbitrary function of a single parameter 4] is searched for using BO, implemented using the numerical package [63]. This function represents the true FoM which is shown in dashed red and the noisy evaluations of the function, which serve as elements of y, are shown as red dots. 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

Figure 2: Example of Bayesian optimization for a function of a single parameter after (a) 10, (b) 11 and (c) 20 iterations. The top panel depicts the underlying function 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 solid blue line and a 95% confidence interval (shaded region delimited by the solid gray lines). The lower panel depicts the acquisition function, ) (red solid line) given by Eq. (3) with k = 4 (a-b) and 0 (c), whose maxima locations correspond to the selection of the next set of parameters in the optimization loop. For example, in (b) the new parameter was chosen at 7 (shown in green box) which corresponds to the maximum of

(shown in solid blue line) and its variance which are used to compute a 95% confidence interval (depicted as shaded blue across the mean). The width of this interval results from 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 it does not adequately reproduce the underlying control landscape due to the small number of evaluations. This lack of knowledge of the true underlying landscape is captured by the variance of the distribution resulting in a large confidence interval. The acquisition function 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 7 (highlighted as dashed vertical line), where the model has both, a high mean value (but not necessarily maximal) and high uncertainty. The FoM is then evaluated for this value of . Fig. 2(b) incorporates this new evaluation (green square) and exhibits the updated model based on a total of 11 iterations. This cycle of updating the model and suggesting the next parameter to evaluate is repeated another nine times and resulting in Fig. 2(c). In this figure, the model is in close agreement with the true landscape and has successfully identified a parameter close to the true global maximum.

Before moving to the next section showcasing the results of BO applied to specific problems, we comment on the limitations and possible variations of this approach. Gaussian processes are a flexible tool to model well-behaved functions but are inadequate to model functions exhibiting discontinuities or varying degree of smoothness. However, this may be dealt with more complicated models as discussed in [64]. Similarly, approximating the true landscape well-enough may fail in this approach as could be the case with any model-based approach. The underlying noise model associated with the evaluation of the FoM, which enters in the likelihood term p(f|y) in Eq. (2) is assumed to be Gaussian in a standard BO implementation. There may be scenarios where this is not the case, nonetheless, it is possible to incorporate non-Gaussian noise into the framework [61]. Another possible limitation is the computational complexity associated to evaluating the posterior distribution given in Eq. (2) which grows as

with D as the number of accumulated evaluations. Alternatives to an exact treatment of Gaussian processes are discussed in the outlook which may help to alleviate this computational burden.

4. Preparation of ordered states in ultra-cold systems beyond adiabatic methods

Controlled preparation of many-body states with spatial correlations resembling solidstate matter is of great interest to the condensed matter community. For example, the SF-MI (Superfluid-Mott insulator) transition realized with bosonic atoms trapped in optical lattices [65, 66] simulates the transition of a conducting state to an insulator state in solid state physics. Similarly creating highly correlated ordered states with long range interactions in a cold gas of Rydberg atoms [67, 68, 69] is akin to crystals in solids with long range Coulomb interactions between the electrons. In the Bose-Hubbard model, the only form of interaction is short-range while in the Rydberg system, it can be relevant over a range longer than a typical lattice spacing [70, 71]. Both of these examples are well studied in the context of ultra-cold physics and will serve as perfect testbeds for establishing the merits of BO techniques.

In Sec. 4.1, using the Bose-Hubbard model as a toy example, BO’s ability to obtain an optimal protocol that drives the SF-MI transition is tested within the context of quantum speed limit [49]. Furthermore, BO’s performance is benchmarked with other optimization methods in terms of its efficiency towards convergence. In Section 4.2, BO is applied to the optimization of laser pulse dynamics to create Rydberg crystals with large fraction of Rydberg excitations in different lattice geometries; one-dimensional (1D), two-dimensional (2D) as well as three-dimensional (3D).

4.1. Superfluid to Mott insulator transition in Bose-Hubbard model

The Bose-Hubbard model is widely studied in solid state physics and is simulated with bosonic atoms in an optical lattice. Although conceptually simple, this many-body system cannot be mapped onto a single particle problem and contains the interesting phenomenon of transitioning from a SF state to a MI state [65]. Experimental realization was first achieved [66] by slowly varying the depth of the optical lattice potential.

Figure 3: Setup of Bose-Hubbard model: The goal of the optimization is to drive the system from an initial superfluid state illustrated in (a) to a Mott insulating phase as shown in (b) by dynamically changing the depth of the optical lattice V (t). (c) Energy spectrum for the Hamiltonian in Eq. (5) for different values of the control Γ.

4.1.1. Setup: For the Bose-Hubbard model, a homogeneous system of N bosonic atoms with repulsive interactions in a lattice with L sites is considered. The Hamiltonian for this model is given as

The first term in the Hamiltonian describes the tunnelling of bosons between neighboring potential sites whose strength is given by J(t). The annihilation (creation) of a boson at site i is defined by operators ˆ) which follow the usual commutation relation, [ˆ. The tunneling term tends to delocalize each atom over the lattice and the 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 ˆis the number operator for a given site i.

As shown schematically in Fig. 3(a-b), by varying the potential depth of the optical lattice in time, both terms J(t) and U(t) are affected. J(t) changes depending on the tunnelling barrier between neighboring lattice sites, while the change in the on-site interaction is due to the variation in atomic wave function confinement. Introducing a single dimensionless quantity, Γ(t) = U(t)/(U(t) + J(t)) such that Γ(Hamiltonian is re-scaled as:

The Bose-Hubbard Hamiltonian has two distinct ground states depending on the strength of the interactions U relative to the tunneling J. In the limit where Γ(t) = 0, the interaction term vanishes and the ground state of each atom is delocalized over the entire lattice. There is uncertainty in the number of atoms per site and the many-body state can be described as a superposition of different atom number states,

Here is the many-body vacuum state in the Fock basis representation. This SF state is characterized by a non-zero variance of the number operator and a constant global phase across the lattice.

In the opposite limit with Γ(t) = 1, the hopping between adjacent sites is suppressed and the ground state of the system consists of localized atomic wave functions that minimize the interaction energy. For a homogeneous system, one can have commensurate filling of n atoms per lattice site and the many-body state is then a product of local Fock states in the atom number for each lattice site. For the ground state where n = N/L, the Mott insulator is given as

Global phase coherence of the matter wave field is lost at the expense of perfect atom number correlations between lattice sites. In the next section we proceed to analyze the energy spectrum of the Bose-Hubbard model. Interestingly, in the limits of Γ(t) = 0 or Γ(t) = 1, the Hamiltonian in Eq. (5) is completely integrable. However, for the intermediate values of Γ, the system is non-integrable [72].

4.1.2. Many-body energy spectrum: For a chain of L = 5 sites with periodic boundary conditions and unit filling (L = N = 5) of bosons, the many-body energy spectrum as a function of Γ is provided in Fig. 3(c). This is solved using exact diagonalization in the Fock basis representation and is implemented with QuSpin [73]. The energies the many-body eigenstates are plotted where n = 0, 1, 2 . . . enumerates each state in ascending order of energy at every value of Γ. The energies at Γ = 0 and Γ = 1 can be expressed analytically [72]. The ground state is highlighted in red ranging from superfluid ground state (Γ = 0) to Mott insulating ground state (Γ = 1) and takes specific values, with energies (Γ = 1) = 0 respectively. There are multiple avoided crossings for all the other values of Γ.

4.1.3. Optimization task: The aim is to find an optimal control function Γwill efficiently drive an initial superfluid state as defined in Eq. (6)) to a Mott-insulating state (see Eq. (7)) in time T. The general set of parameters described in Sec. 2 take specific meaning in this context, as the time-dependent control Γ(t) is parametrized with 10 parameters representing the control values taken at equidistant time, Γ(with . With the boundary conditions, Γ(0) = 0 and Γ(T) = 1, the

Figure 4: Two controls are presented: a linear driving in (a) and an optimized one found by BO in (c). The evolution of the population ) (defined in the main text) of the instantaneous eigenstates for the linear (b) and optimized protocol (d) are also reported. 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 is the realized state at time t = T. While this fidelity is a natural choice for state-preparation task in numerical simulations, it is not readily accessible in an experiment. A second choice of FoM, , takes into account experimental constraints such as the effect of finite sampling. It is defined as where is the variance in the occupation number estimated by averaging over for a given lattice site i.

In all cases the dynamics is solved numerically by exact diagonalization, and we denote the eigenstates of H(t) from Eq. (5) as are the ground state SF and MI states respectively.

4.1.4. Results and discussions: Using fidelity as the FoM, the duration of the protocol is fixed to the quantum speed limit as defined in [16, 49], ∆, where ∆ is the minimum 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 perfect fidelity 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 , defined as 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 the case of the optimized protocol while the linear ramp has a fidelity of less than 30%. For

Figure 5: The infidelity is reported as a function of the number of iterations for different optimization routines (colors in legend). All the results obtained are based on 30 runs of each setup, for which the median (dots) and the first to third quartile interval (shaded region) are reported. Optimizations for different protocol times are shown: (first column), (second) and 1(last), where is defined in the text. Two figures of merit are used: the optimizer has access to either exact evaluations of the fidelity (first row) or to a noisy (hence experimentally realistic) FoM (second row) as described in the main text.

the linear protocol to reach the same fidelity as the optimized one, it would take about an order of magnitude longer in time thus establishing the non-adiabatic character of the optimum protocol.

To benchmark BO’s overall performance for the SF-MI task, results obtained using BO are compared with other optimization routines (some of which were mentioned in Sec. 2) in Fig.5. These include both gradient as well as non-gradient methods. For comparison with gradient based method, the spontaneous perturbation stochastic approximation (SPSA) algorithm was chosen [74] as it was shown to be competitive in the number of iterations for convergence and robust to noise [51]. For comparison with non-gradient based methods, the differential evolution(DE) and Nelder-Mead(NM) optimizers were implemented using Scipy [75]. Finally random search (RANDOM) was also included in the list of optimizers. It relies on randomly sampling a new set of parameters at each iteration.

In Fig.5, the three columns correspond to optimization results for different protocol durations, , as these regimes can have different influences on the complexity of the optimization task [52, 53]. Fig.5(a-c) show the infidelity, defined as (1 ), as a function of the number of iterations. However for Fig.5(d-f), a more experimentally realistic FoM is used, denoted by . As mentioned before, this particular choice of FoM is motivated by experimental constraints such that the impossibility to project onto the target state and finite size sampling effects. The infidelity in this case is directly given by which is related to the variance in the occupation number at each lattice site. is minimized (vanishes) when a MI phase has been reached, thus providing an appropriate guidance for the optimization. Rather than the true expected value of the variance, an estimation of its value based on a finite number of repetitions (here 1000) is used.

As shown in Fig.5, over all these different configurations, BO (orange curve) exhibits faster convergence towards low infidelity. Both SPSA and DE also converge to good parameters but at the expense of more iterations compared to BO. NM in most of the cases get quickly stuck in local minima with poor fidelities. When the time allowed for evolution is too short (see first column in Fig. 5), the best final infidelity reached is of 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 (second and third columns in Fig.5), an order of magnitude less iterations are required by BO to converge to high fidelity control. Using the experimentally motivated FoM final fidelities dropped by 5% for all the optimizers indicating additional difficulties to perform optimization in the noisy scenario. Interestingly SPSA performs better with the 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 to number of iterations.

4.2. Creation of Rydberg crystalline states

After the first observations of interacting Rydberg gases [76, 77, 78, 79, 80, 81], more control over the nature, form and structure of the quantum many-body state has been achieved by shaping the excitation sequence [82] and the arrangement of the atoms in specific patterns [83]. Early proposals [67, 68] for creating ordered phases in Rydberg gases were adiabatic in nature, in the sense that the overall dynamics was slower than the minimum energy gap of the many-body system. Unfortunately, this energy gap decreases exponentially with the number of Rydberg excitations preventing the preparation of crystals with large number of Rydberg excitations.

Other proposals rely on admixing the ground state atoms with small amount of Rydberg state to create super-solid droplet crystals [84, 85]. This was experimentally implemented in [69], but only for a very small system where the energy gap was rather large compared to the coherence time. Moreover Rydberg admixing of ground state atoms requires considerable fine tuning [86] especially to avoid spontaneous scattering and molecular resonances. More recently optimal control has been applied to optimize the 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

Figure 6: Schematic diagram of the setup: Ground state atoms denoted by trapped in a deep optical lattice. The ground state and Rydberg state (atom is optically coupled with an effective Rabi frequency Ω and detuning ∆. Two Rydberg atoms with a relative distance interact via the van der Waals interaction.

is a strong motivation to apply BO for preparing Rydberg ordered states by finding better protocols that give enhanced number of Rydberg excitations even in the presence of experimentally realistic parameter noise and possible lattice filling imperfections.

4.2.1. Setup: Fig. 6 illustrates the setup under consideration with Rb atoms in their ground state coupled to a particular Rydberg state with the help of a laser, thus treating each atom effectively as a two-level system. The atoms can be either trapped in a deep optical lattice [82, 88] or in an array of optical dipole traps [89] with uniform unit filling and lattice spacing m. The advantage of trapping atoms at suitable spacings is that unwanted molecular resonances and scattering can be avoided. Moreover, if needed, simultaneous trapping of ground and Rydberg state atoms is possible in magic wavelength lattices [90, 91, 92]. Two Rydberg atoms with positions interact via the van der Waals interaction given by where 0 is the van der Waals coefficient which for our chosen Rydberg state with principal quantum number n = 50S takes the specific value,

[93]. The Hamiltonian describing the full setup in the rotating wave approximation is given as

where is the projection operator onto the excited state for the iis the transition operator. The first two terms in the Hamiltonian represent the atom-laser interaction with the laser detuning ∆ and the two-photon effective Rabi frequency Ω of the system. For small system sizes of less than ten atoms, a laser profile with uniform excitation is assumed such that all atoms experience the same Rabi coupling. There are different energy scales at play here where one can increase interaction strength by selecting higher principal quantum number n, which in turn provides more bandwidth for Ω and can be tuned with respect to the overall lifetime of

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 shown for different lattice geometries (with lattice spacing many-body ground and excited states are shown explicitly as green and red lines respectively. Panels (d-f) depict the many-body energy spectrum for the same parameters as (a-c) but now with a coupling of Ω = 100 MHz.

the system. Any motional dynamics is neglected as their typical timescales are much longer than the fast excitation dynamics that is considered here. Spontaneous decay of the Rydberg state including black-body radiation for the chosen state, 50S is about 65

4.2.2. Many-body energy spectrum: A gas of N atoms with Rydberg excitations can be represented by the many-body eigenstate labels the different combinatorial possibilities to distribute excitations over N atoms. The many-body energy spectrum with Ω = 0 for different geometries is shown in Fig. 7(a-c). In the zero intensity field scenario, the energy of a many-body state is given as

where are the interaction energies. The many-body ground state (zero energy and is highlighted in green in Fig. 7(a-c). States with the same number of 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 the lattice excited, . For non-zero laser intensity, anti-crossings occur between the many-body states as shown in Fig. 7(d-f).

Since repulsive interactions are assumed for Rydberg atoms, the overall energy of the system is minimized by maximizing the spacing between Rydberg atoms and hence forming an ordered phase of Rydberg excitations. For a 1D system, sweeping the detuning from large negative values to specific positive values, the state with all atoms in the ground state gets connected to states with one, two or more Rydberg excitations depending on the final value of the detuning. One way to ensure that the final state is a Rydberg crystalline state, is to have the entire pulse dynamics performed slower than the inverse of the relevant energy gap. It is non-trivial to calculate the energy gap as it depends not only on the Rabi frequency but also on the interactions [67]. As the number of excitations increases, the energy gap decreases exponentially preventing the formation of Rydberg crystals with large Rydberg excitation fraction, Moreover for 2D and 3D lattice geometries, the many-body energy spectrum is far more complicated as there are more degeneracies and requires further careful optimization over the (Ω, ∆) control landscape.

4.2.3. Optimization task: The dynamics for Eq. (8) is numerically solved using a linear multi-step integration method implemented using the QuTiP package [94]. For this system, the relevant figures of merit are target state fidelity, spatial correlation function and number of Rydberg excitations, all of which can be potentially measured in a typical Rydberg experiment with spatially resolved detection methods [95]. The control protocol is optimized for a specific number of Rydberg excitations is represented by the many-body eigenstates By denoting the fidelity as , the FoM provided at each iteration to BO is given by ). The sum is taken over all the configurations containing Rydberg excitations.

The total duration of the dynamics is fixed to s which is well below the overall lifetime of the many-body Rydberg system. The control field is parametrized using six parameters, [Ω()] which correspond to Rabi frequencies and detunings at three different times. For intermediate times, the values of (Ω()) are interpolated using quadratic polynomial. Additionally, boundary conditions 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 bounds for Ω 5 GHz] and ∆

4.2.4. Results and discussion: The total number of atoms for 1D and 2D lattices is nine and the protocol is optimized for five Rydberg excitations. The 3D lattice consists of eight atoms with the protocol optimized for four Rydberg excitations. The fidelities ) over time for the optimal control protocols found by BO are depicted in Fig. 8 (a-c) and Fig. 8 (d-f) show the optimized control parameters for different lattice geometries. In all cases, BO achieves high fidelity for the selected number of excitations reaching convergence in the optimization within 10 iterations. In order to reproduce real experiments and test BO’s resilience, 5% relative noise is introduced in Ω(respectively. The standard deviation for Ω() is shown in shaded blue and red

Figure 8: Results for optimized protocol dynamics using BO targeting a specific Rydberg excitations with different lattice geometries: (a-c) show the plots of fidelities ) (defined in the text) for a given number of Rydberg excitations optimized in (a) to maximize for ) in a 1D lattice, in (b) to maximize for a 2D lattice and in (c) to maximize ) in a 3D lattice. Panels (d-f) are the corresponding optimized protocols given in terms of Ω()). The three panels (g-i) show the final Rydberg excitation probability at any given lattice site corresponding to (a), (b) and (c) respectively.

colors respectively across their mean values (which are depicted with solid lines). It is most perceptible for Ω in Fig. 8 (e). These noises are usually negligible but could arise due to experimental drifts or Doppler broadening.

Inspecting the protocols optimized in Fig. 8, it seems like BO naturally selects protocols with the general trend where the detuning starts from large negative values and then transitions to positive values. However, the temporal profiles of the optimized Rabi 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 the optimum protocol. It verifies the formation of Rydberg crystalline states as the atomic excitations 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 with five Rydberg excitations. However in 3D lattice, there are two configurations that

Figure 9: Results of optimizations in the presence of imperfect lattice fillings for

different lattice geometries. The solid (dashed) line refers to the protocol optimized for a perfect (imperfect) lattice filling which is denoted by OP1 (OP2). In (a-c), solid and dashed orange curves represent the fidelities for these specific protocols OP1 and OP2 both applied with imperfect lattice filling. (d-f) show the details of these protocols in

terms of the Rabi frequency and detuning over time.

equally contribute to crystal with four Rydberg excitations. In Fig. 8 (i), one of the two configurations is rotated onto the other.

Errors originating from inaccurate readouts and imperfect lattice fillings are prevalent in lot of the current cold atom experiments [13, 88]. Both these imperfections are included in the simulations by associating a probability for Rydberg detection for each atom to be 0.9 when in Rydberg state. In one case, the optimum protocol is obtained for a perfect lattice and is labeled as OP1 while the other protocol is obtained for imperfect lattice filling which is labeled as OP2. The chosen FoM is the number of Rydberg 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 lattice imperfections are vastly different from the ones obtained in the idealized case. Taking into account imperfections during the optimization results in better performances as can be seen in Fig. 9 (a-c) where OP2 exhibits higher final fidelities compared to OP1 across all lattice geometries. This shows the importance of performing optimizations directly onto the experiments rather than relying on idealized models or theoretical protocols for characterizing the noise sources and incorporating them in numerical simulations. Also the number of iterations required when noise is incorporated is five times higher which reiterates the need for optimization routines to converge fast in order to be of practical use.

BO found protocols to create Rydberg crystals in all three different lattice geometries with a Rydberg fraction of s. To create Rydberg crystals with such large Rydberg fraction would take much longer in an adiabatic scheme. For example, the adiabatic preparation of Rydberg crystals in one dimension were realized for s in [67] and s in [68]. The experiment [82] was done for s while the optimal control found using Nelder-Mead [87] was obtained for a quasi-1D system with s. In summary, BO is able to efficiently obtain high fidelity protocols despite the parameter noise and lattice filling imperfections that were included in the simulations.

5. Discussion and Outlook

Quantum optimal control aims at designing fast high fidelity schemes to prepare desirable target states in a quantum system. An optimization routine that is able to 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 the optimization 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, it was successfully applied to the creation of spatially ordered phases of ultra-cold gases trapped in lattices.

In the SF-MI transition task, the performance of BO in comparison to other optimization routines was found to be the most efficient with regards to the number of iterations needed for convergence. This is primarily due to its ability to cleverly select the next set of parameters to probe. In the second example, notwithstanding the complexity in the many-body energy spectrum and the vanishing energy gap for large Rydberg systems, BO obtained control protocols to create Rydberg crystalline states in lattice models of arbitrary dimension with much larger Rydberg fraction than previous methods. More interestingly, BO successfully navigated the optimization landscape even with the inclusion of noisy parameters and imperfect lattice fillings which is evidence of the advantage in accounting for noisy data in the modeling approach.

In the examples studied here, the parametrization of the control field was taken to be either spline or quadratic functions. Other types of smooth parametrizations, for example with Fourier components or even with Gaussian processes, could also be tested. Furthermore in all the optimizations performed in this work no initial guess function was assumed for the control. Such a guess could be directly taken into account in the parametrization of the control function [55] and may result in optimized controls with higher fidelity. Similarly the number of parameters used for the time-varying control field could be increased to allow for finer control. However, in practice optimizing over higher dimensional parameter spaces may be challenging. For a given number of iterations, a trade-off between flexibility in the control field and realistic number of iterations has to be found.

Alternatively it would be interesting to study iterative refinements of the control

field keeping the number of parameters fixed and small at each stage. For example, during the optimization, it may be realized that some temporal parts of the control field are more important than others and a new parametrization could be tested accordingly. Another option would be to repeatedly optimize over additive perturbations of the control function, each with a randomized parametrization, as is done in [56]. Such iterative procedures could be made experimentally feasible if used in conjunction with the low data requirement provided by BO.

Finally, we would like to remark that the field of probabilistic machine learning [97] is advancing at a fast pace and has a lot to offer for the characterization and optimization of quantum systems. In particular several alternatives for the model used here for BO have been suggested [98, 99]. They could be easily integrated into the framework and decrease the computational burden associated with Gaussian processes, especially when the number of iterations becomes really large. Another interesting direction of research for even more efficient optimizations could originate from the field of transfer-learning [100]. For example, it could allow to leverage results of optimizations performed on an (approximative) numerical simulator to speed up optimizations performed on the real experiments.

6. Acknowledgments

The authors acknowledge H. Weimer for insightful and productive discussions. Both Imperial and Stuttgart have received funding from the QuantERA ERANET Cofund in Quantum Technologies implemented within the European Unions Horizon 2020 Programme under the project Theory-Blind Quantum Control TheBlinQC and from EPSRC under the grant EP/R044082/1. This work was supported through a studentship in the Quantum Systems Engineering Skills and Training Hub at Imperial College London funded by the EPSRC(EP/P510257/1).

[1] Bloch I, Dalibard J and Nascimb`ene S 2012 Nat. Phys. 8 267–276 ISSN 1745-2481 URL https://doi.org/10.1038/nphys2259

[2] Gross C and Bloch I 2017 Science 357 995–1001 URL https://science.sciencemag.org/ content/357/6355/995

[3] Moses S, Covey J, Miecnikowski M, Jin D and Ye J 2017 Nat. Phys. 13 13–20 ISSN 1745-2481 URL https://doi.org/10.1038/nphys3985

[4] Bernien H, Schwartz S, Keesling A, Levine H, Omran A, Pichler H, Choi S, Zibrov A S, Endres M, Greiner M, Vuletic V and Lukin M D 2017 Nature 551 579–584 ISSN 1476-4687 URL https://doi.org/10.1038/nature24622

[5] Deutsch I H, Brennen G K and Jessen P S 2000 Fortschritte der Physik 48 925–943 URL 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. 90(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 529 505–508 ISSN 1476-4687 URL https://doi.org/10.1038/nature16176

[8] Katori H 2011 Nat. Photonics 5 203–210 ISSN 1749-4893 URL https://doi.org/10.1038/ nphoton.2011.45

[9] Gallagher T F 1988 Rep. Prog. Phys. 51 143–188 URL https://doi.org/10.1088% 2F0034-4885%2F51%2F2%2F001

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. 105(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 5(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, Pichler H, 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 365 570–574 URL https://science.sciencemag. org/content/365/6453/570

[14] Saffman M, Walker T G and Mølmer K 2010 Rev. Mod. Phys. 82(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 88(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. 6 34187 URL 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. 6 25890 URL 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. 9 4360 ISSN 2041-1723 URL 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 information processing systems 2951–2959 URL 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 104 148–175 URL https://doi.org/10.1109/JPROC.2015.2494218

[22] Cully A, Clune J, Tarapore D and Mouret J B 2015 Nature 521 503–507 ISSN 1476-4687 URL https://doi.org/10.1038/nature14422

[23] Calandra R, Seyfarth A, Peters J and Deisenroth M P 2016 Annals of Mathematics and Artificial Intelligence 76 5–23 URL 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. 5 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 G 2018 Proc. Natl. Acad. Sci 115 13216–13221 ISSN 10916490 URL https://www.pnas.org/ content/115/52/13216

[26] Nakamura I, Kanemura A, Nakaso T, Yamamoto R and Fukuhara T 2019 Opt. Express 27 20435– 20443 URL http://www.opticsexpress.org/abstract.cfm?URI=oe-27-15-20435

[27] Weimer H, L¨ow R, Pfau T and B¨uchler H P 2008 Phys. Rev. Lett. 101(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 288 824–828 ISSN 0036-8075 URL https://science.sciencemag.org/content/288/5467/824

[29] Judson R S and Rabitz H 1992 Phys. Rev. Lett. 68(10) 1500–1503 URL https://link.aps.org/ doi/10.1103/PhysRevLett.68.1500

[31] Bartels R, Backus S, Zeek E, Misoguti L, Vdovin G, Christov I P, Murnane M M and Kapteyn H C 2000 Nature 406 164–166 ISSN 1476-4687 URL https://doi.org/10.1038/35018029

[32] Herek J L, Wohlleben W, Cogdell R J, Zeidler D and Motzkus M 2002 Nature 417 533–535 ISSN 1476-4687 URL https://doi.org/10.1038/417533a

[33] N¨obauer T, Angerer A, Bartels B, Trupke M, Rotter S, Schmiedmayer J, Mintert F and Majer J 2015 Phys. Rev. Lett. 115(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 8(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, Hoi I C, Jeffrey E, Megrant A, Mutus J, Neill C, O’Malley P J J, Quintana C, Roushan P, Sank D, Vainsencher A, Wenner J, White T C, Cleland A N and Martinis J M 2014 Phys. Rev. Lett. 112(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, Baugh J, Zeng B and Laflamme R 2017 npj Quantum Inf. 3 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. 118(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 98(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 K 2007 Phys. Rev. A 75(1) 012302 URL https://link.aps.org/doi/10.1103/PhysRevA.75. 012302

[40] Bartels B and Mintert F 2013 Phys. Rev. A 88(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. 99(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. 2 16019 ISSN 2056-6387 URL https://doi.org/10.1038/npjqi.2016.19

[43] Sklarz S E and Tannor D J 2002 Phys. Rev. A 66(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 75(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. 106(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. 120(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. 172 296 URL 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. 103(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 2 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 91(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 90(3) 032310 URL https: //link.aps.org/doi/10.1103/PhysRevA.90.032310

[53] Bukov M, Day A G R, Sels D, Weinberg P, Polkovnikov A and Mehta P 2018 Phys. Rev. X 8(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. 9 4812 ISSN 2041-1723 URL https://doi.org/10.1038/s41467-018-07090-4

[55] Caneva T, Calarco T and Montangero S 2011 Phys. Rev. A 84(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 92(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. 112(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 268 116–126 ISSN 0925-2312 URL http://www.sciencedirect.com/science/article/pii/ S0925231217307531

[59] F¨osel T, Tighineanu P, Weiss T and Marquardt F 2018 Phys. Rev. X 8(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 5 33 ISSN 2056-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 (MIT press Cambridge, MA) URL http://www.gaussianprocess.org/gpml/

[63] 2016 Gpyopt: A bayesian optimization framework in python URL http://github.com/ SheffieldML/GPyOpt

[64] Neal R M 1998 Bayesian statistics 6 475 URL https://www.cs.toronto.edu/~radford/ftp/ val6gp.pdf

[65] Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 Phys. Rev. B 40(1) 546–570 URL 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 415 39–44 ISSN 1476-4687 URL https://doi.org/10.1038/415039a

[67] Pohl T, Demler E and Lukin M D 2010 Phys. Rev. Lett. 104(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 Kokkelmans S J J M F 2011 J. Phys. B:At. Mol. Opt. Phys. 44 184008 URL 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. 12 1095–1099 ISSN 1745-2481 URL https://doi.org/10.1038/nphys3835

[70] Bohlouli-Zanjani P, Petrus J A and Martin J D D 2007 Phys. Rev. Lett. 98(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. 110(26) 263201 URL https://link.aps.org/doi/10.1103/PhysRevLett.110.263201

[72] Kolovsky A R and Buchleitner A 2004 Europhys. Lett. 68 632–638 URL https://doi.org/10. 1209%2Fepl%2Fi2004-10265-7

[73] Weinberg P and Bukov M 2017 SciPost Phys. 2(1) 003 URL https://scipost.org/10.21468/ SciPostPhys.2.1.003

[74] Spall J C 1998 IEEE Transactions on Aerospace and Electronic Systems 34 817–823 ISSN 2371-9877 URL 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. 93(6) 063001 URL https://link.aps.org/doi/10.1103/PhysRevLett.93. 063001

[78] Reetz-Lamour M, Deiglmayr J, Amthor T and Weidemller M 2008 New J. Phys. 10 045026 URL 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. 98(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. 100(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 T 2009 Phys. Rev. A 80(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 347 1455–1458 ISSN 0036-8075 URL 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 534 667–670 ISSN 1476-4687 URL https://doi.org/10.1038/nature18274

[84] Henkel N, Nath R and Pohl T 2010 Phys. Rev. Lett. 104(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. 105(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. 16 063012 URL 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. 2 035006 URL 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, Kuhr S and Bloch I 2012 Nature 491 87–91 ISSN 1476-4687 URL 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 4(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 84(4) 043408 URL https://link. aps.org/doi/10.1103/PhysRevA.84.043408

[91] Topcu T and Derevianko A 2013 Phys. Rev. A 88(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. 44 184010 URL 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 75(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. 184 1234 – 1240 ISSN 0010-4655 URL https://doi.org/10.1016/j.cpc.2012.11.019

[95] Ott H 2016 Rep. Prog. Phys. 79 054401 URL https://doi.org/10.1088%2F0034-4885%2F79% 2F5%2F054401

[96] Harris F J 1978 Proceedings of the IEEE 66 51–83 ISSN 1558-2256 URL https://ieeexplore. ieee.org/document/14551063

[97] Ghahramani Z 2015 Nature 521 452–459 ISSN 1476-4687 URL https://doi.org/10.1038/ nature14541

[98] Hensman J, Fusi N and Lawrence N D 2013 In Proc. Twenty-Ninth Conference on Uncertainty in Artificial Intelligence 282–290 URL http://www.auai.org/uai2013/prints/papers/244. pdf

[99] Snoek J, Rippel O, Swersky K, Kiros R, Satish N, Sundaram N, Patwary M, Prabhat M and Adams R 2015 In Proc. 32nd International Conference on Machine Learning 37 2171–2180 URL http://proceedings.mlr.press/v37/snoek15.html

formation processing systems 2004–2012 URL 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 ), is maximized. BO relies on an approximative model f, called the surrogate model, of the figure of merit which 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, also known as stochastic process. Random functions extend the notion of a finite set of random variables to an infinite one. Thus, considering f to be a random function means that for any of the infinitely many input parameters , the value taken by the function is itself a random variable. More precisely, f is taken to be a Gaussian Process (GP) which is a specific type of random functions [62]. A GP is the natural extension of a Multi-Variate (MV) Gaussian distribution: while a MV Gaussian 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. These two functions define respectively the mean ), and the covariances ) of the values taken by the model at any input parameters . 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 ), of arbitrary length N, follows a multivariate Gaussian distribution:

where N denotes a Gaussian distribution here with a mean vector m of length N and a covariance matrix K of dimension . The distribution over functions p(f) or equivalently 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 properties of the model. The mean is often taken to be vanishing, 0 and the prior assumptions on the function f are entirely delegated to the choice of the covariance function k. A common choice, followed here, for the covariance function is the 5/2 function imposing the model to be twice differentiable [62]. This

Figure A1: (a) The figure shows three functions sampled from a GP with a Mat´ern covariance (Eq.A.2) for different values of the length-scale l and the variance parameter . (b) Predictive distribution obtained according to Eq.A.5 for a 1-d parameter 1] with 10 noisy observations (red dots). The distribution is displayed in terms of its mean function ) (blue thick line) and a 95% confidence interval (grey light envelope) with boundaries ). Additionally examples of the full density for two specific inputs 6 are explicitly shown (vertical grey shaded slices).

function is defined as

where the variance and the length-scale l, are free hyper-parameters. Essentially the variance specifies to which extend any variable 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 the distance x between the parameters which is indicative of the extent of the correlations between different values of the function. To illustrate the impact of these hyper-parameters, Fig.A1 (a) shows 3 functions with have been sampled with different hyper-parameters values (given in legend). Each sample is obtained according to Eq. (A.1) for a vector of 1-d input parameters [large enough such that 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 the optimization. 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 distribution is entirely defined by the choice of a GP, a mean function m, a covariance function k and the hyper-parameters . In particular, under the choice of a mean and a kernel function presented here, a priori each value ) taken by the model has a distribution: :

Appendix A.2. Updating the surrogate model: Predictive distribution

Starting with the prior distribution described in the previous section, it remains to incorporate observations to the model to obtain the posterior distribution. In the present context this posterior is also called the predictive distribution as it allows to make prediction for unseen parameters. Observations are denoted by a vector )] where each element ) corresponds to an evaluation of the figure of merit at iteration m for parameters , over a total of M iterations.

In the case of noiseless evaluations of the figure of merit one directly records the values taken by the model for the different parameters tried. Writing this vector of values )], it follows the equality f = y. Applying the definition of a GP in Eq. (A.1) to the vector of random variables [)] and using the conditioning properties of Gaussian distributions [62], the predictive distribution is given by

where the column vector k has entries ), and elements of the covariance matrix K are given by ). The symbols ) denote the mean and standard deviation of this predictive distribution. They both are function of the parameter as each element of k depends on it. Compared to the prior distribution )) in Eq. (A.3), the mean of the predictive distribution has been shifted from 0 to and its variance has decreased from by a positive quantity , resulting from the incorporation of the observations. The most computational demanding part of evaluating this mean and variance comes from the inversion of the matrix K which has complexity

In an experimental scenario, the set of measurements y does not directly reveal the values taken by the model, still model and observations can be related by positing a noise model. This noise, originating from both experimental imperfections and also quantum fluctuations, can be approximated as an additive constant Gaussian noise such that , where the noise terms are assumed to be independently drawn from the same distribution ). The extra parameter , capturing the amount of noise, can be fixed or could also be fitted to the data in the same way as the hyper-parameters . Under the independence assumption for the noise terms, the likelihood of recording the full data set of observations y for a given set of values f taken by the model can be written as ). With this simple phenomenological model of the noise the sought-after predictive distribution can be derived:

where the first line shows how it can be decomposed, using Bayes rule and marginalization over the vector of values f, as a product of the prior distribution p(f) appearing in Eq. A.1, the conditional distribution in Eq. (A.4), the likelihood p(y|f) previously defined and the probability of obtaining a given set of observations p(y). The next line provides the analytical result of this integral [62]. This final result is quite similar to Eq. (A.4) except that the inverse of the covariance matrix has now been replaced by (thus incorporating noise in the observations.

A practical example of evaluating this predictive distribution for a one-dimensional parameter 4] based on a set of M = 10 observations is provided in Fig. A1(b). The mean function and a confidence interval with bounds represented with, in addition, the explicit densities ) for two specific values of the parameter. Far away from the observations, for example for 75, the confidence interval is wider indicating uncertainty in the predictive distribution, while closer to observations, for example for 5, it shrinks but does not vanish as the model identified noise in the observations (more precisely the value of the hyper-parameter 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 observe based on the surrogate model. This new parameter could be picked such that it maximizes the mean of the predictive distribution given in Eq. (A.5): arg max). However, with only a limited amount of observations, it is likely that the model does not capture the full range of variations in the figure of merit and this approach 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 would consist on taking measurements where the uncertainty in the model is maximal. This can be done by choosing ), where the standard deviation quantifies the uncertainty in the model. Balancing these two objectives, which are finding the location of the maximum and also exploring region of the parameter space where not enough observations have been made, is often referred as the explorationexploitation trade-off. In BO this trade-off is dealt with by introducing an acquisition ) aiming at capturing these two aspects; and the choice of the next parameter is taken such that:

Among the most popular acquisition functions we recall the definition of the Upper Confidence Bound (UCB) function already provided in the main text:

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 the graphs in the bottom panels of Fig.2. In Fig.2(a-b) a value of k = 4 was used, while in 2(c) in the last iteration of the optimization it was taken to be k = 0 in order to show the 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 of merit obtained so far , is averaged under the predictive distribution:

which has analytic solution [20]. The choice of the acquisition function can have in some cases significant impact on the final results and both the UCB and EI functions were considered in this study. In practice finding the maximum of these acquisition functions is often performed by gradient ascent (possibly repeated over several random starting parameters). For both the UCB and EI, the gradients of the acquisition function with regards to the parameters have analytical expressions and can be efficiently evaluated. These two acquisition functions and the routines for finding their maxima are 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 are given here.

Appendix B.1. SF - MI transition

For the SF-MI transition example, all the optimizations shown in Fig. 5 were performed with the same configuration. In each case, = 100 observations were initially taken for randomly chosen parameters followed by M = 1750 iterations of BO. To decrease the computational burden associated to this large number of iterations we found it efficient not to update the hyper-parameters (noise and parameters of the kernel functions) at each 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 the choice of the next parameter to probe the UCB acquisition function in Eq.A.7 was used with a value of k linearly decreasing from k(1) = 5 at the first iteration to k(1750) = 0 at the end of the optimization. For this problem this acquisition function performed better than the EI one. Finally other types of kernels [62] were compared and we found that the 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

For the creation of Rydberg crystalline states, the optimization problem was found to be simpler and the number of initial observations and iterations was significantly smaller than in the previous case. This can be attributed, in part, to the smaller dimensionality of the parameters vector. For Fig. 8 BO routines were initialized with observations followed by M = 10 iterations. In the noisy case presented in Fig. 9 a larger number of iterations was needed: all the optimizations were run starting with = 6 initial observations and M = 50 iterations. The EI acquisition function was used for all the optimizations included in this part.