Optimal control of the self-bound dipolar droplet formation process

TL;DR

This paper demonstrates how optimal control algorithms can significantly accelerate the formation of dipolar Bose-Einstein condensate droplets, using advanced numerical methods to improve efficiency and feasibility.

Contribution

It introduces a multilevel B-spline based optimal control approach for dipolar droplet formation, reducing computational effort and incorporating box constraints effectively.

Findings

Optimal control dramatically accelerates droplet formation.

Multilevel B-spline method reduces cost function evaluations.

Efficient evaluation of dipolar interactions maintains spectral convergence.

Abstract

Dipolar Bose-Einstein condensates have recently attracted much attention in the world of quantum many body experiments. While the theoretical principles behind these experiments are typically supported by numerical simulations, the application of optimal control algorithms could potentially open up entirely new possibilities. As a proof of concept, we demonstrate that the formation process of a single dipolar droplet state could be dramatically accelerated using advanced concepts of optimal control. More specifically, our optimization is based on a multilevel B-spline method reducing the number of required cost function evaluations and hence significantly reducing the numerical effort. Moreover, our strategy allows to consider box constraints on the control inputs in a concise and efficient way. To further improve the overall efficiency, we show how to evaluate the dipolar interaction potential in the generalized Gross-Pitaevskii equation without sacrificing the spectral convergence rate of the underlying time-splitting spectral method.

Figures

Equations

$${\hat{\mathcal{J}}}_1(u)={\mathcal{J}}_1({\psi }_u,u),$$

$${\hat{\mathcal{J}}}_2(\mathbf{c})={\mathcal{J}}_2({\psi }_{\mathbf{c}},\mathbf{c})$$

$$u(\mathbf{c},t)=u_0+(u_T-u_0)t/T+\sum _{k=1}^K{c}_k\sin (k\pi t/T),$$

$$u_0=u(t=0),u_T=u(t=T)$$

$$u_{\mathrm{lb}}\le u(\mathbf{c},t)\le u_{\mathrm{ub}},t\in [0,T]$$

$$\displaystyle\imath\hbar\frac{\partial\psi({\bm{r}},t)}{\partial t}=\bigg{[}-\frac{\hbar^{2}\nabla^{2}}{2m}+V_{\mathrm{ext}}({\bm{r}},t)+g(t)|\psi({\bm{r}},t)|^{2}$$

$$\displaystyle+\Phi({\bm{r}},t)+\gamma_{\mathrm{qf}}(t)|\psi({\bm{r}},t)|^{3}-\frac{\imath\hbar L_{3}}{2}|\psi({\bm{r}},t)|^{4}\bigg{]}\psi({\bm{r}},t),$$

$$V_{\mathrm{ext}}(\mathbf{r},t)=\frac{1}{2}m{\omega }_{\rho }^2(t)(x^2+y^2)+\frac{1}{2}m{\omega }_z^2(t)z^2$$

$$\Phi (\mathbf{r},t)=\frac{{\mu }_0{\mu }^2}{4\pi }{\int }_{\text{mathbbm}{R}^3}U(\mathbf{r}-{\mathbf{r}}^{\prime })|\psi ({\mathbf{r}}^{\prime },t){|}^{2}d{\mathbf{r}}^{\prime },$$

$$U(\mathbf{r})=\frac{1-3{\cos }^2\theta }{r^3},$$

$$\gamma_{\mathrm{qf}}=\frac{32}{3}g\sqrt{\frac{a_{s}^{3}}{\pi}}\bigg{(}1+\frac{3}{2}\frac{a_{\mathrm{dd}}^{2}}{a_{s}^{2}}\bigg{)}$$

$$a_s(t=0)$$

$${\omega }_{\rho }(t=0)$$

$${\omega }_z(t=0)$$

$$a_s^{\mathrm{lb}}$$

$${\omega }_{\rho }^{\mathrm{lb}}$$

$${\omega }_z^{\mathrm{lb}}$$

$$a_{s,\mathrm{i}}$$

$$a_s^{\mathrm{lb}}$$

$${\omega }_{\rho ,\mathrm{i}}$$

$${\omega }_{\rho }^{\mathrm{lb}}$$

$${\omega }_{z,\mathrm{i}}$$

$${\omega }_z^{\mathrm{lb}}$$

$$\psi (t_{n+1})\approx e^{-ıB^{+}△t/2}{e}^{-ıA△t}{e}^{-ıB^{-}△t/2}\psi (t_n)$$

$$A=-\hslash {\nabla }^2/(2m),B^{\pm }=V_{\mathrm{eff}}/\hslash$$

$$V_{\mathrm{eff}}=V_{\mathrm{ext}}+g|{\psi }_n^{\pm }{|}^2.$$

$$ı\hslash {\partial }_t\psi =(V_{\mathrm{ext}}+g|\psi (t_n){|}^2)\psi$$

$$ı\hslash {\partial }_t\psi =-{\hslash }^2{\nabla }^2/(2m)\psi ,$$

$$\Omega =[-L_x/2,L_x/2]\times [-L_y/2,L_y/2]\times [-L_z/2,L_z/2].$$

$$V_{\mathrm{eff}}=V_{\mathrm{ext}}+g|{\psi }_n^{\pm }{|}^2+{\Phi }_n^{\pm }+{\gamma }_{\mathrm{qf}}|{\psi }_n^{\pm }{|}^3-\frac{ı\hslash L_3}{2}|{\psi }_n^{\pm }{|}^4$$

$$\Phi=\big{(}\mu_{0}\mu^{2}/(4\pi)\big{)}\mathcal{F}^{-1}\Big{[}\mathcal{F}[U]\,\mathcal{F}[\rho]\Big{]},$$

$$U(\mathbf{r})$$

$$\Phi (\mathbf{r},t)$$

$$\displaystyle\quad-\frac{\mu_{0}\mu^{2}}{4\pi}\int_{\mathbbm{R}^{3}}\partial_{{\bm{n}}{\bm{n}}}\bigg{(}\frac{1}{|{\bm{r}}-{\bm{r}}^{\prime}|}\bigg{)}|\psi({\bm{r}}^{\prime},t)|^{2}\,d{\bm{r}}^{\prime}$$

$$\Phi (\mathbf{r},t)=-g_{\mathrm{dd}}|\psi (\mathbf{r},t){|}^2-3g_{\mathrm{dd}}{\partial }_{\mathbf{n}\mathbf{n}}\phi (\mathbf{r},t),$$

$$\phi (\mathbf{r},t)={\int }_{\text{mathbbm}{R}^3}u(\mathbf{r}-{\mathbf{r}}^{\prime })|\psi ({\mathbf{r}}^{\prime },t){|}^{2}d{\mathbf{r}}^{\prime }$$

$$u(\mathbf{r})=\frac{1}{4\pi |\mathbf{r}|}.$$

$$\Delta \phi (\mathbf{r},t)=-|\psi (\mathbf{r},t){|}^2,\underset{|\mathbf{r}|\to \infty }{\lim }\phi (\mathbf{r},t)=0.$$

Full text

I Introduction

Optimal control has meanwhile become a widely used method for the efficient and robust manipulation of quantum systems in various applications Li_2011 ; glaser_2015 ; goerz_2017 ; machnes_2018 and, as such, is expected to play a key role in many emerging quantum technologies like sensing or computing waldherr_2014 ; dolde_2014 ; noebauer_2015 . There is also particular interest in the optimization of the intricate processes in complex quantum many-body systems hohenester_2007 ; doria_2011 . In this context, the feasibility of optimal control algorithms has been experimentally demonstrated for lattice systems and Bose-Einstein condensates buecker_2011 ; vanfrank_2016 .

The first numerical simulation demonstrating that optimal control provides an efficient way to realize the transfer of a Bose-Einstein condensate (BEC) to a desired target state was presented in a pioneering work by Hohenester et al hohenester_2007 . Mathematical aspects and computational improvements have been published shortly thereafter in winckel_2008 , with numerical implementations freely available hohenester_2014 . While early studies grond_2009 , jaeger_2013 , buecker_2013 , jaeger_2014 focussed on low-dimensional BECs to limit the computational costs, optimal control of the full three-dimensional Gross-Pitaevskii equation (GPE) was presented in mennemann_2015 . Moreover, the original approach was generalized to include multiple control inputs. In the meantime, various computational and technical aspects, as well as different numerical methods have been considered in the literature. These considerations are not limited to numerical simulations but have also been successfully applied in real physical experiments buecker_2011 ; vanfrank_2014 ; vanfrank_2016 .

In the present work we go beyond the isotropic contact interactions typically found in such BECs and study the dynamics in BECs with long-range, anisotropic dipolar interactions lahaye_2008 ; maier_2015 ; baier_2016 ; chomaz_2018 . Moreover, we also go beyond the 3D mean-field treatment and include the Lee-Huang-Yang beyond mean-field correction, which has been shown to lead to the emergence of new and unexpected states of matter, including self-bound quantum droplets baillie_2016 ; schmitt_2016 ; chomaz_2016 and self-organized striped states wenzel_2017 .

As a proof of concept, we consider the formation of a self-bound dipolar droplet state. In particular, we demonstrate that the single droplet formation process studied in baillie_2016 can be significantly accelerated using efficient numerical methods in combination with advanced concepts of optimal control.

Our numerical simulations exploit two novel and effective innovations: First, we show how to solve the generalized GPE describing BECs including dipolar interactions and the Lee-Huang-Yang mean-field correction. More precisely, we show how to apply the popular and highly efficient time-splitting spectral (split-step Fourier) method without sacrificing the spectral convergence rate with respect to the discretization of the spatial operators. Second, we turn the underlying infinite dimensional control problem into a finite-dimensional nonlinear optimization task using a B-spline control vector parameterization. In this context, we also develop a multilevel refinement strategy to enhance the speed of convergence of the underlying optimization algorithm.

Optimal control of systems governed by ordinary or partial differential equations has a long history in mathematics, science and technology lions_1971 ; troeltzsch_2010 ; gerdts_2011 . In order to classify and understand the approach proposed in this article it is useful to give a short review on the two predominant methods frequently found in the literature of quantum optimal control, and in particular, in the field of BEC experiments.

The first approach is based on a reduced cost functional

[TABLE]

where $\mathcal{J}_{1}(\psi,u)$ is a suitably selected cost functional describing the optimal control problem and $u:[0,T]\rightarrow\mathbbm{R}$ denotes a given control input. Furthermore, $\psi_{u}$ denotes the time-evolution of the wave function corresponding to a solution of the state equation using $u$ as the control input. In the context of BEC experiments, the state equation would be given by the GPE. The state equation is considered to be a constraint to the optimization problem. Using the calculus of variations an additional adjoint equation is derived which allows for an efficient computation of the gradient of $\hat{\mathcal{J}}_{1}$ needed to improve the given control input $u$. Until then, no numerical approximations are involved. Finally, the state and adjoint equation are solved using suitable numerical discretizations. From a mathematical point of view, this approach appears to be very elegant, however, the derivation and discretization of an additional adjoint equation is a highly non-trivial task and must be repeated whenever the state equation or the cost functional is modified. Moreover, one should bear in mind that the numerical solution of the adjoint equation is typically more expensive than the solution of the state equation, see, e.g., mennemann_2015 .

In the second approach, the reduced cost functional (1) is replaced by a reduced cost function

[TABLE]

with $\bm{c}=[c_{1},\ldots,c_{K}]$ being a vector collecting the $K$ coefficients of a suitable control vector parameterization (CVP). As an example, we consider a simple sum of sines parameterization

[TABLE]

where

[TABLE]

denote the initial and the final conditions of the control input, respectively. The first argument of $\mathcal{J}_{2}(\psi_{\bm{c}},\bm{c})$ denotes the time evolution of the wave function corresponding to a solution of the GPE using $u(\bm{c},t)$ as the control input. It should be noted that $\hat{\mathcal{J}}_{1}$ in (1) is a functional, whereas $\hat{\mathcal{J}}_{2}$ in (2) is a function. Therefore, the previously infinite-dimensional optimal control problem has been mapped to a finite-dimensional nonlinear Program (NLP) which can be solved using various numerical methods like the popular Nelder-Mead algorithm nelder_mead_1965 .

With these preconditions CVP methods are easy to implement as the only requirement is a solver for the forward problem, i.e., a numerical method to solve the given differential or partial differential equation. Not surprisingly, CVP methods are frequently applied in classical optimal control applications and have been used long before they were applied to quantum optimal control problems, see, e.g., kraft_1985 ; schlegel_2005 ; hadiyanto_2008 .

A CVP method based on a randomized Fourier basis is proposed as a general and versatile quantum optimal control technique in caneva_2011 ; rach_2015 . However, like most randomized algorithms, the proposed method converges extremely slowly as was demonstrated recently in an extensive benchmark problem sorensen_2018 . In fact, the simple ansatz in (3) seems to be more efficient and could be easily applied in a multistart optimization procedure to explore the space of possible solutions in a more rigorous way.

Often, a surprisingly small number of coefficients $K$ is sufficient to obtain a high-quality solution of the underlying control problem. However, in some quantum control applications sorensen_2018 the number of coefficients needed to approximate pulses of very small duration can easily exceed $K\approx 20$. For such large values of $K$, derivative-free methods like the Nelder-Mead algorithm are known to be less efficient.

In principle, the gradient of the reduced cost function (2) with respect to the coefficients $\bm{c}$ can always be approximated using simple finite difference formulas. To this end, the forward problem needs to be solved at least $K+1$ times. More efficient methods to compute the gradient are presented in winckel_2008 and sorensen_2018 for a CVP based on Chebyshev polynomials and a sum of sines ansatz, respectively. Generally speaking, the approach is not limited to a specific CVP, however, it also requires the derivation and discretization of an additional adjoint equation. Once a numerical method to evaluate the gradient is available, the NLP can be solved using most effective implementations of quasi-Newton algorithms (BFGS) or sequential quadratic programming (SQP) methods nocedal_wright_2006 .

While the CVP defined in Eqn. (3) works well in the context of unconstrained optimization problems, it results in suboptimal solutions if the control input is subject to constraints

[TABLE]

naturally present in realistic control applications. This applies in particular if the control horizon $T$ is chosen to be small. The reason for this is that, as the $T$ becomes smaller, the corresponding optimal control tends to switch from one extreme to the other. In minimum-time problems the corresponding control is referred to as a bang-bang solution. In our application, we are not aiming at solving a minimum-time optimal control problem. However, we will choose the time horizon $T$ small enough that the effect of the constraints (5) starts to play an important role. More precisely, we will see that the control inputs resemble bang-bang solutions, i.e., the lower and upper bounds become active for time intervals of non-vanishing lengths. We would like to add that these kind of controls cannot be obtained using the simple ansatz given in Eqn. (3) as this would require an infinite number of sine waves $K$.

As an alternative, we propose to parameterize the control input by a superposition of B-spline basis functions d_boor_1978 . As a result of their compact support, B-spline basis functions represent a much more favorable ansatz with respect to Eqn. (5). Moreover, Dirichlet boundary conditions (4) can be implemented easily and B-spline basis functions are ideally suited for multiresolution function approximations. In fact, we will demonstrate that a multilevel refinement strategy can significantly improve the convergence of the overall optimization procedure.

The paper is organized as follows: Section II covers the generalized GPE, the phase diagram of a trapped dipolar Bose gas and all relevant physical parameters. In the first part of Section III, we review the time-splitting spectral method (TSSM) frequently used to solve the ordinary GPE. Moreover, we introduce the notation concerning the spatial and temporal discretizations and explain why the TSSM can be easily applied to the generalized GPE as well. The second part of Section III is devoted to spectrally accurate numerical methods for the evaluation of the dipolar interaction potential. In Section IV, we introduce the multilevel B-spline method which is finally applied to the dipolar droplet formation process in Section V.

II Generalized Gross-Pitaevskii Equation

Our considerations are based on the description of dipolar Bose-Einstein condensates given by the generalized Gross-Pitaevskii equation baillie_2016

[TABLE]

where $\psi$ is the condensate wave function with initial state $\psi_{0}=\psi(t=0)$. The density of the initial state is assumed to be normalized to the number of atoms $\mathcal{N}_{0}=\|\psi_{0}\|^{2}$. Furthermore, the contact interaction strength is given by $g=4\pi a_{s}\hbar^{2}/m$, where $a_{s}$ denotes the s-wave scattering length and $m$ the mass of the atoms. In the numerical simulations below, the external potential

[TABLE]

is a function of the radial and longitudinal trapping frequencies $\omega_{\rho}$ and $\omega_{z}$, which can vary with time and are considered to be freely adjustable parameters. Dipole-dipole interactions are taken into account via the dipolar interaction potential

[TABLE]

with $r=(x^{2}+y^{2}+z^{2})^{1/2}$. Moreover, $\theta$ denotes the angle between ${\bm{r}}$ and the polarization axis ${\bm{n}}=(n_{1},n_{2},n_{3})^{\top}$ which is assumed to be normalized, i.e., $(n_{1}^{2}+n_{2}^{2}+n_{3}^{3})^{1/2}=1$. Within a local density approximation lima_2011 the correction term $\gamma_{\mathrm{qf}}|\psi|^{3}$ is used to include the effect of quantum fluctuations lima_2011 ; waechtler_2016 ; schmitt_2016 ; baillie_2016 . The prefactor

[TABLE]

depends on the dipolar length $a_{\mathrm{dd}}=m\mu_{0}\mu^{2}/(12\pi\hbar^{2})$ lahaye_2009 . Finally, three-body loss processes are taken into account by the term $-(\imath\hbar L_{3}/2)|\psi|^{4}$ with $L_{3}$ being the loss coefficient.

The validity of this approximation has previously been checked using Monte-Carlo simulations saito_2016 . The effect of quantum fluctuations is usually very small, however it can be dramatically enhanced when the mean-field contributions of contact and dipolar interactions nearly cancel. This is the case, for example, during the collapse following a motional instability of the system waechtler_2016 .

The resulting phase diagram for the system is schematically shown in Fig. 1 baillie_2016 . For high scattering length $a_{s}$ the system is a dipolar BEC, for low $a_{s}$ it forms a single stable quantum droplet. Above a critical aspect ratio $\lambda_{c}\sim 1.87$ the system is bistable ferrier_2018 . Lowering $a_{s}$ into this bistable region triggers a motional instability that leads to the formation of arrays containing several droplets.

Our aim is to optimize the self-bound dipolar droplet formation process presented in baillie_2016 . This work proposed to linearly decrease the scattering length, which - for appropriate trap aspect ratios - connects the BEC and the self-bound droplet by a smooth crossover. This is followed by a linear decrease of the trap frequencies to turn off the trapping potential and reveal the self-bound nature of the droplet. However, adiabatically following such a linear trajectory requires a long timescale and thus leads to non-negligible atom loss. Alternatively, fast linear ramps typically lead to unwanted excitations. As we will show in the following, optimal control can help with both these problems, preparing droplets without observable excitations on very short timescales.

Like in baillie_2016 we consider 164Dy atoms. We therefore choose $m=163.93$ u, $\mu=9.93\,\mu_{\mathrm{B}}$ and $L_{3}=1.2\times 10^{-41}\,\mathrm{m}^{6}s^{-1}$, with $u$ the atomic mass unit. Initially, $\mathcal{N}_{0}=10^{4}$ atoms are trapped in the harmonic potential at $t=0$. Furthermore, $a_{s},\omega_{\rho}$ and $\omega_{z}$ are considered to be control inputs which are subject to the initial and final conditions

[TABLE]

as well as the constraints

[TABLE]

for $t\in[0,T]$ where $T$ denotes the final time of the optimal control problem considered below. All remaining parameters of Eqn. (10) and Eqn. (11) are summarized as follows:

[TABLE]

In experiments $a_{s},\omega_{\rho}$ and $\omega_{z}$ can be changed using a Feshbach resonance and the laser intensities of the trapping lasers, respectively. The initial value of $a_{s}$ in baillie_2016 is given by $a_{s,\mathrm{i}}=130\,a_{0}$, which is linearly decreased to $a_{s,\mathrm{f}}=80\,a_{0}$. In our optimization we restrict ourselves to this range of scattering lengths to take into account the fact that Feshbach resonances in lanthanide atoms are typically very narrow and their distribution highly complex maier_2015 . In contrast to alkali atoms the s-wave scattering length can thus only be tuned over a limited range, which has to be taken into account in the control sequence.

Moreover, we also choose the same initial configuration as in baillie_2016 for the external confinement potential at time $t=0$, i.e., $\omega_{\rho,\mathrm{i}}=2\pi\times 70$ Hz and $\omega_{z,\mathrm{i}}=\lambda\,\omega_{\rho,\mathrm{i}}$ with $\lambda=0.75$ being the trap aspect ratio. As for the scattering length we include limits for these trapping frequencies, which are bounded by the available laser power in experiments. Since we are aiming at creating a self-bound droplet state, the external potential is required to vanish at $t=T$ and thus $\omega_{\rho,\mathrm{f}}=\omega_{z,\mathrm{f}}=0$. In all our calculations, the polarization axis of the dipoles is aligned with the positive $z$-direction, i.e., ${\bm{n}}=(0,0,1)^{\top}$.

III Numerical Discretization

III.1 Time-Splitting Spectral Method

Our strategy to solve the optimal control problem considered in Section V requires hundreds of simulations of the generalized GPE (6). A variety of numerical methods is described in the literature to compute the time-evolution of nonlinear Schrödinger equations. The most straightforward way is to first discretize the spatial operators using finite difference approximations and then to apply explicit Runge-Kutta methods (like the classical Runge-Kutta method of order four) to integrate the resulting system of ordinary differential equations in time. While this approach is simple to implement, very small time step sizes are needed to ensure numerical stability caplan_2013 . Much larger time step sizes can be used if the time-integration is based on implicit methods. A proven method in the context of quantum dynamical simulations is the Crank-Nicolson scheme, see e.g. bao_numerical_solution_2003 . This latter approach represents an excellent choice in one spatial dimension, however, in two or three spatial dimensions significant computational resources are needed to solve the corresponding linear systems of equations. Finally, an efficient and frequently used method to compute the time-evolution of the classical GPE is the time-splitting spectral (TSSM) method bao_numerical_solution_2003 ; bao_numerical_study_2003 ; thalhammer_high-order_2009 . We recall that the most frequently used implementation of the TSSM is based on the Strang time splitting bao_numerical_solution_2003 , i.e.,

[TABLE]

Moreover, $t_{n}=n\triangle t$ with $\triangle t$ being the time step size and $n=0,...,N$ such that $N\triangle t=T$. According to (12a), the $n$th time step consists of three sub-steps. First, we solve

[TABLE]

for a duration of $\triangle t/2$ using $\psi_{n}^{-}=\psi(t_{n})$ as initial value. The result is used as initial value for the free Schrödinger equation

[TABLE]

which is then solved for a duration of $\triangle t$; the result is $\psi_{n}^{+}$. Finally, $\imath\hbar\partial_{t}\psi=(V_{\mathrm{ext}}+g|\psi_{n}^{+}|^{2})\psi$ is solved again with initial value $\psi_{n}^{+}$, again for a duration of $\triangle t/2$. The result of the third sub-step is taken as approximation of $\psi(t_{n+1})$.

The spatial discretization is based on the grid points $x_{j}=-L_{x}/2+j_{x}\triangle x$, $y_{j}=-L_{y}/2+j_{y}\triangle y$ and $z_{j}=-L_{z}/2+j_{z}\triangle z$ with $\triangle x=L_{x}/J_{x}$, $\triangle y=L_{y}/J_{y}$ and $\triangle z=L_{z}/J_{z}$. The parameters $L_{x}$, $L_{y}$ and $L_{z}$ are chosen sufficiently large such that the wavefunction is compactly supported within the computational domain

[TABLE]

Moreover, $j_{x}=0,...,J_{x}-1$, $j_{y}=0,...,J_{y}-1$, $j_{z}=0,...,J_{z}-1$ and $J=J_{x}J_{y}J_{z}$ denotes the total number of spatial grid points.

In a numerical implementation the solution of (13) requires a pointwise evaluation of the exponential function on a real valued array as well as a pointwise multiplication of a real and a complex valued array. The free Schrödinger equation (14) can easily be solved in Fourier space and the numerical implementation requires the application of two discrete Fourier transforms (forward and inverse) and a pointwise multiplication of two complex valued arrays, see e.g. bao_numerical_solution_2003 . For the sake of numerical efficiency, the discrete Fourier transforms are computed using the Fast Fourier Transform (FFT) and hence the numerical costs of a single time step are of order $\mathcal{O}(J\log J)$.

The TSSM outlined above (Strang splitting) is of second order in time. Moreover, the spatial derivatives in the free Schrödinger equation are discretized with spectral accuracy, meaning that the numerical error corresponding to the discretization of the Laplacian operator decreases at an exponential rate with respect to the number of spatial grid points. It should be noted that these convergence rates are based on the assumption that the initial wave function and the external confinement potential are sufficiently smooth. However, in most experimentally relevant simulation scenarios these assumptions are satisfied and thus the number of spatial grid points can be reduced drastically in comparison with finite difference approximations. Moreover, the TSSM allows for the application of large time step sizes $\triangle t$ and as a result the total number of time steps needed to compute the time evolution of the GPE is significantly lower than in the case of the above-mentioned explicit time-integration methods. To summarize, the TSSM represents a very efficient numerical method for the solution of the GPE.

In the following, we show that the TSSM can also easily be applied to the generalized GPE in (6). To this end, we only need to substitute the effective potential energy in (12c) by

[TABLE]

with $\Phi_{n}^{\pm}$ being defined in analogy to (8a) using $\psi_{n}^{\pm}$ instead of $\psi$. The terms $\gamma_{\mathrm{qf}}|\psi^{\pm}|^{3}$ and $-i\hbar L_{3}/2\,|\psi^{\pm}|^{4}$ can be included at almost no additional numerical costs. However, the evaluation of $\Phi$ in (8a) is a non-trivial task if we want to conserve the spectral convergence rate of the TSSM with respect to the discretization of the spatial operators.

III.2 Evaluation of the Dipolar Interaction Potential

A simple way to evaluate the dipolar interaction potential $\Phi$ is based on the convolution theorem which yields

[TABLE]

where $\mathcal{F}$ denotes the Fourier transform and $\rho$ is the density $\rho=|\psi|^{2}$. In a numerical implementation $U$ and $\rho$ are replaced by their discrete representations corresponding to the numerical grid of the computational domain $\Omega$. Moreover, all Fourier transforms are implemented using FFTs in combination with a zero-padding strategy. Due to the singularity of $U({\bm{r}})$ at ${\bm{r}}=(0,0,0)^{\top}$ the accuracy of the approximation is of second order with respect to the number of spatial grid points $J$.

To find a discretization of the dipolar interaction potential that converges at an exponential rate, we first apply the same transformation that has been employed in bao_efficient_2010 . More precisely, we make use of the fact that (8b) can be written (in the sense of distributions) as

[TABLE]

where $\partial_{{\bm{n}}}={\bm{n}}\cdot\nabla=n_{1}\partial_{x}+n_{2}\partial_{y}+n_{3}\partial_{z}$ and $\partial_{\bm{nn}}=\partial_{{\bm{n}}}\partial_{{\bm{n}}}$. Plugging (17) into (8a) yields

[TABLE]

or

[TABLE]

where $g_{\mathrm{dd}}=4\pi\hbar^{2}a_{\mathrm{dd}}/m=\mu_{0}\mu^{2}/3$ and

[TABLE]

The normal derivative operator $\partial_{{\bm{n}}{\bm{n}}}$ can be evaluated using spectral differentiation (at the cost of two additional FFTs of size $J=J_{x}J_{y}J_{z}$). Consequently, the only remaining difficulty is the computation of the convolution in (18), which however is a ubiquitous problem in computational physics and scientific computing in general. In fact, the solution of (18) is equivalent to the solution of the free-space Poisson equation

[TABLE]

A variety of efficient numerical methods has been developed to solve the Poisson equation in unbounded domains ethridge_new_2001 ; aguilar_2005 ; langstron_free-space_2011 ; hejlesen_high_2013 ; malhotra_pvfmm_2015 ; exl_2016 ; vico_2016 ; exl_2017 . While most of them are based on formulation (18) it is also possible to employ (19) directly. This strategy is used by the authors of bao_efficient_2010 in the context of numerical simulations of dipolar Bose-Einstein condensates. By replacing the far-field boundary condition in (19) with a homogeneous Dirichlet boundary condition on the boundary of the computational domain, the Poisson equation can be solved easily using a sine pseudospectral method. The method avoids the application of zero-padding. However, since the potential $\varphi$ in (19) is known to decay slowly, an extremely large computational domain is needed to keep the modeling error small.

A highly efficient numerical method for the evaluation of nonlocal potentials, including the 3D Coulomb potential (18), has been presented recently exl_2016 ; exl_2017 . The method is based on a Gaussian-sum approximation of the singular convolution kernel (combined with a Taylor expansion of the density) and it has been shown that for sufficiently smooth right hand sides (densities) the numerical error decreases at an exponential rate with respect to the number of spatial grid points $J$. Apart from a relatively expensive precomputation, the most time-consuming part of the method is the application of several FFTs in combination with a zero-padding strategy.

Another spectrally accurate algorithm was presented almost simultaneously in vico_2016 , motivated by the truncation technique of fuesti_molnar_2002 ; ronen_2006 . The implementation of the latter method is much simpler but involves a very memory-intensive precomputation. However, in our application these limitations are not an issue which is why we employ the second algorithm. The second algorithm is based on the observation that the solution of (18) is indistinguishable from the solution of

[TABLE]

with the characteristic function

[TABLE]

and $L=\sqrt{L_{x}^{2}+L_{y}^{2}+L_{z}^{2}}$. Since $u_{L}$ is a compactly supported function the Paley-Wiener theorem implies that its Fourier transform

[TABLE]

is an entire function vico_2016 (and thus $C^{\infty}$). Application of the convolution theorem to (20) yields

[TABLE]

which is the starting point for the numerical discretization. However, since $\hat{u}_{L}$ is a highly oscillatory function, the complexity of the algorithm significantly increases. More precisely, $\hat{u}_{L}$ needs to be sampled using a four-fold sampling rate, which finally requires the application of FFTs of size $4J_{x}4J_{y}4J_{z}$. Fortunately, these complications only concern the precomputation phase. During normal operation the algorithm requires the application of two FFTs (forward and inverse) of size $2J_{x}2J_{y}2J_{z}$ and hence the algorithm is of the same complexity as the simple strategy presented at the very beginning of this section.

All considerations in vico_2016 are based on the assumption that $\Omega$ is a unit box with an aspect ratio

[TABLE]

equal to one. In the application considered below the aspect ratio of the computational domain is given by $\zeta=2$. By means of numerical experiments using gaussian density distributions and their known analytical solutions we found that the algorithm still converges exponentially fast if the aspect ratio remains smaller than approximately $\zeta\lessapprox 2.75$. Furthermore, we found that for even larger aspect ratios, the sampling rate in the precomputation phase needs to be increased (six-fold sampling rate for $2.75\lessapprox\zeta\lessapprox 4.5$).

IV Multilevel B-spline method

For the reasons stated in the introduction, we propose to parameterize the control inputs using B-spline functions

[TABLE]

where $B_{k,p}$ denotes the $k$th B-spline basis function of polynomial order $p$. Given a nondecreasing knot vector $\Xi=\{\xi_{1},\xi_{2},\ldots,\xi_{K+p+1}\}$ the $K$ basis functions are defined recursively by means of the Cox-de Boor recursion formula d_boor_1978 ; cottrell_2009

[TABLE]

with

[TABLE]

being piecewise constant functions.

In the simulations below, we employ cubic B-spline functions, i.e., we set $p=3$. Moreover, our optimization algorithm is based on the uniform knot vectors

[TABLE]

corresponding to four different levels of refinement. These knot vectors are said to be open since the first and last knot values appear $p+1$ times. Most applications of B-spline basis functions are based on open knot vectors. This is due to the fact that basis functions formed from open knot vectors are interpolatory at the ends of the interval $[\xi_{1},\xi_{K+p+1}]$ which makes it very simple to implement Dirichlet boundary conditions like the ones given in Eqn. (4).

The B-spline basis functions corresponding to the knot vectors $\Xi_{2}$ and $\Xi_{3}$ are shown in Fig. 2 a) and Fig. 2 b), respectively. Furthermore, Fig. 2 c) shows a given control input using two representations

[TABLE]

corresponding to the basis functions of the second and third B-spline level, respectively. These two representations of a single control input illustrate the basic idea behind the multilevel B-spline method. Instead of solving the optimal control problem using the finest refinement level directly, it is often more promising to solve the optimization problem using a coarse representation of the control input first. Subsequently the solution will be used as an initial guess for the optimal control problem in the next refinement level. In Eqn. (23a), the solution of the optimization problem corresponding to the second B-spline level is represented using $K^{(2)}$ coefficients $c_{1}^{(2)},\ldots,c_{K^{(2)}}^{(2)}$. The same function is represented in Eqn. (23b) using $K^{(3)}>K^{(2)}$ coefficients corresponding to the knot vector $\Xi_{3}$. However, the coefficients $c_{1}^{(3)},\ldots,c_{K^{(3)}}^{(3)}$ do not represent an optimal solution with respect to the third B-spline level. Rather, they serve as an excellent initial guess which will speed up the solution of the next optimization problem significantly.

The solution corresponding to the previous B-spline level can be represented exactly by means of the new basis functions if the knots of the old knot vector are contained in the new knot vector. An efficient numerical method to determine the coefficients in the new basis is called knot insertion prautzsch_2002 . Alternatively, we may determine the coefficients of the next level $c_{1}^{(\ell+1)},\ldots,c_{K^{(\ell+1)}}^{(\ell+1)}$ by solving a linear system corresponding to the equations

[TABLE]

for $m=1,\ldots,K^{(\ell+1)}$ with $\hat{\xi}_{m}$ being the Greville points johnson_2005 ; cottrell_2009 of the knot vector $\Xi_{\ell+1}$.

V Application to the dipolar droplet formation process

The numerical methods introduced in Section II, III and IV will now be applied to optimize the self-bound dipolar droplet formation process.

V.1 Parameterization, Constraints and Cost Function

As mentioned in Section II, the control inputs are given by the s-wave scattering length $a_{s}$ and the trap frequencies $\omega_{\rho}$ and $\omega_{z}$. While the system is initially a Bose-Einstein condensate, we are aiming to find optimal trajectories of $a_{s}$, $\omega_{\rho}$ and $\omega_{z}$ which will smoothly transform this condensate into a single dipolar droplet. The final droplet is stable in the sense that it will exist even after the external potential $V_{\mathrm{ext}}$ has been set to zero at final time $T$.

To this end, the physical control inputs

[TABLE]

are parameterized using three individual B-splines $u_{1}(t)$, $u_{2}(t)$ and $u_{3}(t)$ with $t\in[0,T]$. The initial and final values of $a_{s}$, $\omega_{\rho}$ and $\omega_{z}$ in Eqn. (10) are taken into account using

[TABLE]

As discussed at the end of Section II, we limit the allowed ranges of the control inputs

[TABLE]

which can be easily translated to equivalent constraints for the B-spline functions

[TABLE]

Each of these three B-spline functions are uniquely defined by the basis functions correponding to the knot vectors (22) and the coefficients

[TABLE]

where $\ell$ denotes the level of refinement and $K_{\ell}$ is the number of B-spline basis functions. An important feature of B-spline functions is called the convex hull property d_boor_1978 , which implies that we can implement the constraints (26) by simply replacing $u_{1}$, $u_{2}$ and $u_{3}$ with the corresponding coefficients

[TABLE]

For a given knot vector $\Xi_{\ell}$ of dimension $s_{\ell}$ we have $K_{\ell}=s_{\ell}-(p+1)$ basis functions, see, e.g., Fig. 2. However, as a result of Eqn. (25), the coefficients

[TABLE]

are already fixed. Therefore, the unknown variables of the $\ell$th optimization problem are given by

[TABLE]

It can be easily verified that the total number of unknown coefficients in the $\ell$th B-spline level is given by $6$, $9$, $15$ and $27$, respectively.

For the cost function we chose

[TABLE]

with $\psi(T)$ being implicitly dependent on the coefficients ${\bm{c}}^{(\ell)}$. Furthermore, $\psi_{d}$ denotes the desired single droplet state. In this context, it should be noted that the densities of the initial and desired states $\psi_{0}$ and $\psi_{d}$ are normalized to the same number of particles $\mathcal{N}_{0}$. Furthermore, the norm of $\psi(T)$ is strictly smaller than the norm of $\psi_{d}$, which is a direct consequence of the three-body loss term in the generalized GPE (6).

The effect of the three-body loss term is noticeable during the whole time evolution. Initially, the effect might be small as the density of the initial state is comparatively low. However, as soon as the wavefunction comes close to the target state $\psi_{d}$, a significant number of atoms will inevitably get lost. Regardless of how we select the control inputs, it is therefore impossible to make all atoms contribute to the desired droplet state at final time $T$. Nonetheless, we will see that our optimization algorithm will find solutions minimizing the effect of the three-body loss term such that almost all particles of the initial condensate will be transferred to the desired single droplet state.

The final optimal control problem in the $\ell$th level of refinement will be solved using an advanced implementation of the SQP method, which allows to take into account the box constraints (27) in a highly efficient way. Alternatively, the given constrained optimization problem could be reformulated as an unconstrained optimization problem. To this end, the cost function (29) would have to be extended by additional penalty terms explicitly depending on the coefficients $\bm{c}^{(\ell)}$.

V.2 Implementation Details

The numerical methods presented above have been implemented in Matlab R2018a on a powerful workstation computer. The most time-consuming part in our optimization strategy is the numerical solution of the generalized Gross-Pitaevskii equation (6). Just like for the ordinary GPE mennemann_2015 , the computations can be accelerated significantly by means of a powerful graphics processing unit (GPU). More precisely, we employ a NVIDIA Quadro GV100 (32GB HBM2 GPU Memory) being able to deliver more than $7.4$ TFLOPS of double precision floating-point (FP64) performance.

The numerical solution of the generalized GPE is based on the time-splitting spectral method (TSSM) using a spectrally accurate algorithm to evaluate the dipolar interaction potential (8) as outlined in Section III. In order to take into account the elongated shape of the desired droplet state, we choose the side lengths of the computational domain (15) to be $L_{x}=12\,\mu$m, $L_{y}=12\,\mu$m and $L_{z}=24\,\mu$m. Our highly accurate implementation of the TSSM converges at an exponential rate with respect to the number of spatial grid points $J_{x}$, $J_{y}$ and $J_{z}$. In fact, we employ two sets of discretization parameters. The first set is characterized by the spatial discretization parameters $J_{x}=72$, $J_{y}=72$, $J_{z}=64$ and a time step size of $\triangle t=0.005$ ms. This rather coarse discretization will be used to solve the generalized GPE in the optimization process. As a cross check we will employ a second set of discretization parameters given by $J_{x}=2\times 72$, $J_{y}=2\times 72$, $J_{z}=2\times 64$ and a time step size of $\triangle t=0.0025$ ms. In fact, the graphs of all physical quantities presented in the next section correspond to numerical simulations using the second set of parameters. However, the graphs are indistinguishable from the ones obtained with the set of coarse discretization parameters.

Our aim is to prepare a self-bound dipolar droplet state on a very short time scale of $T=2$ ms corresponding to one-tenth of the duration of the piecewise linear ramps employed in baillie_2016 . Accordingly, only $N=400$ ($N=800$) time steps are needed to compute the wave function $\psi(T)$ at final time $T$. As a result of the small number of required spatial grid points $J=J_{x}J_{y}J_{z}$ and time steps $N$, the computational time to compute $\psi(T)$, and hence to evaluate (29), is less than four seconds in case of the coarse discretization and much shorter than a minute in case of the fine discretization. It seems that the potential of the above mentioned GPU is not fully exploited with the first set of discretization parameters. Therefore, the solution times corresponding to the second set of discretization parameters are significantly shorter than expected from a purely computational complexity based point of view.

Before starting our optimization algorithm we need to compute the initial state $\psi_{0}$ and the desired state $\psi_{d}$. As in the case of the ordinary GPE, these states can be computed using imaginary time propagation chiofalo_2000 , i.e., the time step $\triangle t$ in our implementation of the TSSM is replaced by $-i\triangle t$ and the density of the wave function is normalized to the number of desired atoms after every time step.

V.3 Results

In accordance with Eqn. (28), the initial point of the NLP corresponding to the first B-spline level is given by a $6$-dimensional coefficient vector $\bm{c}_{0}^{(1)}$. The coefficients are chosen such that the corresponding control inputs $a_{s}$, $\omega_{\rho}$ and $\omega_{z}$ are simple linear ramps.

First, we study the convergence history of our optimization algorithm. To this end, we consider the best cost function value

[TABLE]

among the first $k$ cost function evaluations of the underlying SQP method. The normalized convergence history is depicted in Fig. 3 for two different optimization strategies. For normalization we employ the cost function value corresponding to the case where all control inputs are given by simple linear ramps.

In the first strategy, we employ the original multilevel B-spline approach outlined in Section IV. In this context, it should be noted that we employ $15$ iterations of the SQP method before triggering a new B-spline refinement operation. However, the maximum number of iterations in the fourth B-spline level is set to infinity and the optimization is stopped after $2500$ cost function evaluations. The corresponding convergence history is shown in Fig. 3 as a black solid line.

Alternatively, we can start the optimization algorithm using the fourth level of refinement right from the beginning. The initial point of the SQP algorithm is given by a $27$-dimensional coefficient vector $\bm{c}_{0}^{(4)}$. As before, the coefficients are chosen such that the corresponding control inputs are simple linear ramps. The resulting convergence history is depicted as a broken blue line in Fig. 3.

By comparing the results of both strategies it is clear that the multilevel B-Spline approach is far superior to the approach in which the finest B-spline parameterization is used right from the beginning. The successive refinement strategy prevents the optimization algorithm from getting stuck in one of the many very unfavorable local minima existing in a high-dimensional parameter space. In fact, the multilevel B-spline optimization method exhibits some similarities with multigrid preconditioning techniques frequently applied in numerical linear algebra problems such as elliptic partial differential equations hackbusch_2013 .

The time evolution of the control inputs $a_{s}$, $\omega_{\rho}$, $\omega_{z}$ and the corresponding (normalized) cost function as a result of the multilevel B-spline optimization process is shown in Fig. 4 a)–d). It can be seen that $\mathcal{J}(\psi(T))$ decreases with increasing B-spline level and that the optimized control inputs are able to remain on their respective upper or lower bounds for time intervals of non-vanishing lengths. This particular behavior is a result of the compact support of the B-spline basis functions and could not be obtained by more simple control vector parameterizations like the sum of sines ansatz given in Eqn. (3) or Fourier type parameterizations in general.

As an example, we apply the parameterization in Eqn. (3) using $K=9$ sine functions for each of the three controls $a_{s}$, $\omega_{\rho}$ and $\omega_{z}$. In contrast to the multilevel B-spline method, the box constraints on the control inputs cannot be easily transferred to the corresponding coefficients $c_{k}$, $k=1,\ldots,K$. We therefore resort to quadratic penalty functions sanctioning the squares of the constraint violations nocedal_wright_2006 . Since the SQP algorithm does not perform well for the given unconstrained optimization problem, we apply the Nelder-Mead algorithm using $2500$ cost function evaluations as a stopping criterion. Furthermore, the initial guess $\bm{c}_{0}\in\mathbbm{R}^{27}$ corresponds to linear ramps, meaning that all coefficients of the initial control vector are set to zero. The numerical results of this alternative optimization strategy are depicted in Fig. 4 e)–h). As can be seen from Fig. 4 e), the lower bound of the control input $a_{s}$ becomes active only at one or two isolated points in time. Moreover, Fig. 4 h) shows that the corresponding final value of the objective function is almost two orders of magnitude larger than the final value corresponding to the multilevel B-spline method.

The above reported numerical results turned out to be stable with respect to variations of the initial guess used to initialize the sum of sines and multilevel B-spline optimization algorithms. More precisely, having tested a dozen of randomly selected initial conditions we always found that the final cost function value of the multilevel B-spline method was about two orders of magnitude smaller than the corresponding value of the sum of sines approach.

No reduction of the cost function value $\mathcal{J}(\psi(T))$ can be observed in the case where the control inputs are given by simple linear ramps, see Fig. 4 h). In fact, much larger final times $T$ are needed to observe a significant improvement. However, with increasing $T$ more and more particles are lost as a result of the three-body loss term in the generalized Gross-Pitaevskii equation (6). For that reason, and in strong contrast to most quantum control applications, the final desired state cannot be reached by adiabatically changing the control inputs. Finally, we would like to stress that the three-body loss term is of crucial importance for a realistic description of dipolar Bose-Einstein condensates and hence cannot be omitted. This is true in particular if the final state comes close to the desired state which is characterized by a much higher peak density than the initial state. In this context it is important to remember that the three-body loss term in the generalized GPE (6) is proportional to $|\psi|^{4}$.

Several snapshots of the density at $y\equiv 0$ and $z\equiv 0$ are depicted in Fig. 5 and Fig. 6, respectively. They illustrate the time evolution of the density at various times after the final time $T$. In the simulations we set $a_{s}(t)=a_{s,\mathrm{f}}$, $\omega_{\rho}(t)=\omega_{\rho,\mathrm{f}}$ and $\omega_{z}(t)=\omega_{z,\mathrm{f}}$ for $t\geq T$. The density corresponding to the optimized control inputs can be clearly identified as the density of a stable self-bound dipolar droplet state. In contrast, shape and amplitude of the density corresponding to the linearly decreasing control inputs show very strong oscillations.

These oscillations become even more obvious if we consider the time evolution of the peak density as shown in Fig. 7 a). Moreover, it is interesting to note, that the optimization algorithm finds control inputs which keep the peak density low until shortly before the end of the optimization at $T=2$ ms. In this way the effect of the three-body loss term is minimized and less than $100$ atoms are lost during the preparation of the self-bound dipolar droplet state. This corresponds to more than $99\%$ of the initial atoms contributing to the desired final state. However, due to the high density of the final state, the effect of the three-body loss term becomes even more relevant for $t>T$ as can be seen from Fig. 7 b). In case of the linear ramps, the effect of the three-body loss term is rather small, which can be explained by the comparatively small peak (and average) density which deviates strongly from the density of the target state $\psi_{d}$.

To investigate the above mentioned effect in more detail we take a look at Fig. 7 c) showing the time evolution of the number of atoms in the cylindrical domain

[TABLE]

The parameters in (30) are chosen such that the density of the desired droplet state is closely surrounded. As expected, almost all atoms of the initial state are transferred to $\mathcal{Z}$ for the case of the optimized control inputs. However, in case of the linear ramps, strong fluctuations of the density cause the atoms to get distributed over the whole computational domain $\Omega$. While the overall loss of atoms as a result of the three-body loss term is small, a comparatively small fraction of these atoms actually contributes to the desired self-bound droplet state.

With regard to real-world quantum experiments it is important to characterize the stability of the control trajectories found by the optimization algorithm. In this context, we consider the optimized B-spline functions $u_{1}$, $u_{2}$ und $u_{3}$ in combination with small perturbations of the initial and final values

[TABLE]

in Eqn. (24). The result is a systematic perturbation of the control inputs over the whole time interval $[0,T]$. With these modifications we model perturbations typically occurring in experiments, be it uncertainties in the interaction strength or in the atom number. Additionally, we add white Gaussian noise to the control inputs to model fluctuations, e.g. in the magnetic fields to control the scattering length. The variance of the individual noise signals is given by $\sigma^{2}$ with a standard deviation of $\sigma=0.03$. Moreover, the noise signals are scaled with the maximum values of the control inputs, see Fig. 4 e)–g). However, the main instability in BEC experiments is the number of particles $\mathcal{N}_{0}$ trapped in the potential at $t=0$. This effect is taken into account by artificially lowering $\mathcal{N}_{0}$ by $10$ percent, which is a typical value reached in experiments. In other words, the initial state of the perturbed simulation is recalculated using $\tilde{\mathcal{N}}_{0}=9000$ particles, whereas the computations of the multilevel B-spline method where based on the assumption that initially $\mathcal{N}_{0}=10000$ particles are trapped in the quadratic potential at $t=0$. Despite these relatively strong perturbations, we still observe a remarkable reduction of the objective function value at final time $T$, see Fig. 4 h). The time evolution of the corresponding observables is shown in Fig. 7 along with the results of the previous simulations. While the effect of the perturbations is clearly visible, the outcome of the numerical experiment remains largely the same.

VI Conclusion and outlook

We have introduced advanced numerical methods for the optimization of the self-bound dipolar droplet formation process. By means of a multilevel B-spline control vector parameterization we were able to dramatically speed up the transfer of a given Bose-Einstein condensate to this exotic quantum object. Moreover, it has been seen that the optimized control inputs minimize the effect of the three-body loss term in the generalized GPE and hence almost all atoms of the initial state could be transferred to the desired state at final time $T$. In order to reduce the overall numerical costs we have also shown how to evaluate the dipolar interaction potential using recently introduced spectrally accurate algorithms. Our considerations are not limited to the formation of a single self-bound dipolar droplet state. Rather, we are convinced that the presented algorithms can directly be applied to many other interesting parameter regimes, e.g. to efficiently realize the elusive supersolid state of matter wenzel_2017 ; tanzi_2018 .

VII Acknowledgement

We thank Matthias Wenzel for discussions. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 746525. TL acknowledges financial support by the Baden-Württemberg Stiftung through the Eliteprogramme for Postdocs. LE acknowledges support from the Austrian Science Foundation (FWF) under grant No. P31140-N32 (project ’ROAM’). JFM acknowledges support from the Austrian Science Foundation (FWF) under grant No. P32033, (“Numerik für Vielteilchenphysik und Single-Shot Bilder”). NJM acknowledges support from the Austrian Science Foundation (FWF) under grant No. F65 (SFB “Complexity in PDEs”) and grant No. W1245 (DK ”Nonlinear PDEs”). JFM and NJM acknowledge the Wiener Wissenschafts- und TechnologieFonds (WWTF) project No. MA16-066 (“SEQUEX”).