††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work

# Extreme fragmentation of a Bose gas

Nathan Dupont  Amit Vashisht  CENOLI, Université Libre de Bruxelles, CP 231,
Campus Plaine, B-1050 Brussels, Belgium  Nathan Goldman<EMAIL_ADDRESS>CENOLI, Université Libre de Bruxelles, CP 231, Campus Plaine, B-1050 Brussels,
Belgium  Laboratoire Kastler Brossel, Collège de France, CNRS, ENS-PSL
University, Sorbonne Université, 11 Place Marcelin Berthelot, 75005 Paris,
France

###### Abstract

Fragmentation of an interacting Bose gas refers to the macroscopic occupation
of a finite set of single-particle eigenstates. This phenomenon is related to
the notion of particle-number squeezing in quantum optics, an exquisite
property of quantum states that can offer metrological gain. So far,
fragmentation has only been partially achieved in experiments involving a
large number $N$ of bosons in few modes. Here, we introduce a practical and
efficient scheme to prepare fragmented states in systems realizing the
$L$-mode Bose-Hubbard model. We demonstrate how a large energy detuning
between the modes can be used as a practical control parameter to successfully
fragment a Bose gas over an extremely short preparation time. Applying an
optimal-control approach within realistic experimental constraints, we obtain
total fragmentation at a high filling factor, realizing $\ket{N/L,...,N/L}$
Fock states with hundreds of bosons in very few modes over a few tunneling
times.

Bose-Einstein condensation (BEC) refers to the macroscopic bosonic occupation
of a single-particle ground state [1, 2, 3]. Interestingly, when the single-
particle ground state is degenerate, a rich competition occurs between the
different states in which to condense. In such a scenario, if the order of
degeneracy $D$ is much smaller than the total number of particles $N$, the
system can form a fragmented BEC [4, 5, 6] which is characterized by a
macroscopic occupation $N/D$ of the degenerate states.

Ground state degeneracy, stemming from the single-particle Hamiltonian’s
invariance under a symmetry transformation, is generally lifted, in practice,
by arbitrarily small perturbations (e.g. an external magnetic field) which
break the relevant symmetries [7]. In practice, this strongly complicates the
experimental realization of fragmented BECs with a large $N/D$ ratio.
Nonetheless, strong interactions among the particles, with an energy scale
larger than symmetry-breaking perturbations, can stabilize these elusive
phases; this was recently observed in the three-fold degenerate ground state
of a spin-1 atomic gas [8]. However, the resulting state was far from
exhibiting complete fragmentation, as quantified by the “fragmentation
entropy”, i.e. the von Neumann entanglement entropy computed from the
eigenvalues of the one-body density matrix. In a $L$-mode system, this entropy
is maximized by the $L$-tuple Fock state
$\ket{N/L}\equiv\bigotimes_{i=1}^{L}\ket{N/L}_{i}$, evenly distributing the
bosons in each mode.

A tight-binding lattice model, as realized by cold atoms in optical lattices
[3], naturally offers a set of $L$ orbitals. By activating hopping processes
between the sites, the degeneracy of the single-particle spectrum is lifted
through the formation of a band structure. As a consequence, bosons condense
in the lowest-energy Bloch state, forming a superfluid state over the lattice.
Ramping up strong interactions can then fragment the superfluid state,
reaching the maximally fragmented Fock state $\ket{N/L}$ in the strongly-
correlated limit [9, 3]. However, this transition from the superfluid to the
Mott insulator [10] has only been realized for limited filling factors
$N/L\\!\sim\\!1$ [11, 12, 13, 3, 14].

Considering the case $L\\!=\\!2$, the twin-Fock state $\ket{N/2}$ is the most
resourceful Fock state for quantum sensing and metrology [15, 16, 17, 18, 19,
20, 21]. A multitude of approaches have been proposed and implemented to
perform number squeezing around $\ket{N/2}$ (either in true two-mode systems,
or in a sub-manifold of spinor gases), such as feedback on non-demolition
measurements [22, 23], parametric amplification [24, 18], adiabatic and quasi-
adiabatic transformations [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 8, 35],
shortcuts to adiabaticity [36, 37], optimal control [38, 39], reinforcement
learning [40] and more [41, 42, 43]. Nevertheless, the preparation of a pure
twin Fock state in a two-mode system has not yet been achieved.

In this Letter, we explore and compare schemes for fragmenting an ensemble of
bosons in a $L$-mode Bose-Hubbard system. We first discuss an adiabatic
approach, which consists in slowly ramping up the ratio between Hubbard
interactions and hopping amplitude, either directly [5, 11] or effectively
using a Floquet drive [44, 45]. We then introduce a large energy detuning (or
“tilt”) between the modes, which effectively annihilates the hopping. On the
onset of that regime, we use quantum optimal control [46, 47, 48, 49] to steer
the ground state of the unbiased system towards the $L$-tuple Fock state. Our
results are illustrated for $L\\!=\\!2$ modes; see 111See Supplementary
Material for details and additional results in $L\\!=\\!2$ and $L\\!=\\!3$.
for the $L\\!=\\!3$ case. With this practical approach, we achieve extreme
fragmentation at high filling, two to three orders of magnitude faster than
adiabatic methods, and with significant robustness to parameter fluctuations,
opening up the prospect of experimental realization.

Figure 1:  Fragmentation in a two-mode Bose-Hubbard system. (a) Schematics of
the double-well Hamiltonian (2), with hopping $J$, 2-body contact interactions
$U$ and no tilt $\delta\\!=\\!0$. (b) Many-body spectrum for $N\\!=\\!32$
bosons as a function of $U/J$ (grey lines), with $E_{n}$ the eigenenergy of
eigenstate $\ket{\phi_{n}}$. The projection of $\ket{N/2}$ onto
$\ket{\phi_{n}}$ is encoded in shades of blue on $E_{n}$. (c-e) Husimi
representations [51, 52, 53, 20] of the states identified in (b) (blue areas,
with colorscale ranging from 0 to the maximum value of the distribution) and
corresponding MF trajectories (black lines). (f-j) same as (a-e) as a function
of the tilt $\delta/J$ with $U/J\\!=\\!1$.

#### The system.

We consider a fixed number $N$ of bosons in a two-mode system, such as the
double well illustrated in Fig. 1(a). Under the low-energy single-band
approximation, we describe the system by the two-mode Bose-Hubbard model:

$\displaystyle\hat{H}(t)=$
$\displaystyle-J\left(\hat{a}^{\dagger}_{1}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}_{1}\right)+\dfrac{U}{2}\left(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{1}\hat{a}_{1}\hat{a}_{1}+\hat{a}^{\dagger}_{2}\hat{a}^{\dagger}_{2}\hat{a}_{2}\hat{a}_{2}\right)$
$\displaystyle+\delta\left(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{2}\hat{a}_{2}\right),$
(1)

where $\hat{a}_{i}^{\dagger},\hat{a}_{i}$ are the bosonic creation and
annihilation operators in mode $i$, $J$ is the particle hopping strength, $U$
quantifies the onsite two-body interaction and $\delta$ is an eventual tilt
between the modes. Using the Schwinger (angular-momentum) representation,
$\hat{J}_{x}\\!=\\!(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2$,
$\hat{J}_{y}\\!=\\!(\hat{a}_{1}^{\dagger}\hat{a}_{2}-\hat{a}_{2}^{\dagger}\hat{a}_{1})/2\text{i}$,
$\hat{J}_{z}\\!=\\!(\hat{a}_{1}^{\dagger}\hat{a}_{1}-\hat{a}_{2}^{\dagger}\hat{a}_{2})/2$,
Hamiltonian (The system.) is rewritten as

$\hat{H}(t)=U\hat{J}^{2}_{z}-2J\hat{J}_{x}+2\delta\hat{J}_{z},$ (2)

up to the constant term $U\hat{N}(\hat{N}-2)/4$ (with
$\hat{N}=\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2}$),
describing our system as an effective spin-$1/2$ quantum gas on the Bloch
sphere. For a sufficiently large $N$, intuition can be built on the quantum
dynamics governed by Eq. (2) by taking the semiclassical, mean-field (MF)
limit $\hat{a}_{i}\rightarrow\sqrt{n_{i}}\text{e}^{-\text{i}\theta_{i}}$,
yielding the MF Hamiltonian

$H(\varphi,z)=\frac{\Lambda}{2}z^{2}-\sqrt{1-z^{2}}\cos{(\varphi)}+\frac{\delta}{J}z,$
(3)

where $\Lambda=NU/2J$, and where we have introduced the two canonically
conjugate Bloch-sphere variables: the relative phase
$\varphi=\theta_{1}-\theta_{2}$ and population $z=(n_{1}-n_{2})/N$ [54, 55].

To quantify the fragmentation of an arbitrary state $\ket{\psi}$, we define
the following three metrics: (i) The fragmentation (von Neumann) entropy
$S_{\text{F}}=\sum_{i}-\lambda_{i}\log{\lambda_{i}}$ [56], with $\lambda_{i}$
the normalized eigenvalues of the reduced one-body density matrix with
elements
$\hat{\rho}_{i,j}\\!=\\!\bra{\psi}\hat{a}_{i}^{\dagger}\hat{a}_{j}\ket{\psi}$.
In a $L$-mode system, $S_{\text{F}}$ is maximized to $\ln(L)$ by the $L$-tuple
Fock state $\ket{N/L}$. (ii) The number squeezing parameter
$\xi_{N}^{2}=\Delta\hat{J}^{2}_{z,\psi}/\Delta\hat{J}^{2}_{z,\text{ref}}$,
quantifying the suppression of particle-number fluctuations with respect to
the uncorrelated, binomial distribution of variance
$\Delta\hat{J}^{2}_{z,\text{ref}}=N/4$ for $L\\!=2\\!$ [20, 50]. (iii) The
quantum fidelity to the twin Fock state
$\mathcal{F}=|\innerproduct{N/2}{\psi}|^{2}$.

#### Adiabatic approach.

We first consider the preparation of $\ket{N/2}$ by adiabatically increasing
$U/J$, as traditionally performed to induce a superfluid to Mott-insulator
transition [10, 11]. In Fig. 1(b), we show the many-body spectrum of the
Hamiltonian in Eq. (2) for $N\\!=\\!32$ and $\delta=0$ as a function of $U/J$.
The physics is that of the standard bosonic Josephson junction [57, 58]:
Starting from $U/J\ll 1$, the ground state $\ket{\phi_{0}}$ is a superfluid
coherent spin state [52, 20], pointing along $(\varphi,z)=(0,0)$ on the Bloch
sphere [Fig. 1(c)]; as one increases interactions, this state asymptotically
tends towards $\ket{N/2}$ in the limit $U/J\rightarrow\infty$ [Fig. 1(d,e)].
Crucially, very large values of $U/J$ are needed to achieve high ground-state
fragmentation: considering $N\\!=\\!32$ bosons, we find that reaching a
fidelity $\mathcal{F}>0.99$ requires $U/J\gtrsim 230$, while
$S_{\text{F}}/\ln(2)>0.99$ is reached for $U/J\gtrsim 285$.

When adiabatically following the states (c) to (e) in Fig. 1(b), we find that
the narrowest energy gap that could trigger Landau-Zener transitions is
$\Delta E\\!=\\!4J$ at $U/J\\!=\\!0$, after which $\Delta E$ keeps increasing
monotonically; this observation does not depend on the number of bosons $N$.
An adiabatic preparation of the target state thus requires a preparation time
$t_{\text{f}}\gg\hbar/J$ [50].

As previously noted, total fragmentation requires $U/J\rightarrow\infty$,
which, in practice, would correspond to progressively annihilating the hopping
strength, e.g. by increasing the distance between the sites. In an experiment,
this could be hard to achieve without exciting motional degrees of freedom
[59, 38, 32, 35]. Alternatively, $J$ could be effectively suppressed by means
of Floquet engineering, using a fast modulation of the tilt $\delta$ [44, 45].
We discuss the adiabatic preparation of the $\ket{N/2}$ state using this
Floquet approach in [50]; we demonstrate the applicability of this method deep
in the high-frequency regime, which is found to be only accessible for limited
number of bosons $N$, due to the interacting nature of this periodically-
driven system [50].

#### Adding a tilt.

To overcome the severe limitations of the adiabatic approach, we now add a
strong tilt $\delta$ to the system [Fig. 1(f)]; for the sake of illustration,
we henceforth fix the ratio $U/J$ to a realistic value of 1. In the limit
$\delta\gg J,U$, the tilt effectively annihilates the hopping term, as one can
deduce by moving to a rotating frame generated by the unitary transformation
$\textstyle\hat{R}=\exp\\{\text{i}2\delta\hat{J}_{z}t/\hbar\\}$; see Refs.
[26, 50]. In this frame, the eigenstates of the Hamiltonian (2) become the
Fock states $\ket{n_{1}}\equiv\ket{n_{1},N-n_{1}}$, i.e. the eigenstates of
$\hat{J}_{z}$, with the intuitive lab-frame ground state $\ket{0,N}$ due to
the large tilt. In Fig. 1(g), we show the many-body spectrum associated with
the Hamiltonian (2) as a function of $\delta$, for $N\\!=\\!32$ bosons,
identifying the projection of the target twin Fock state $\ket{N/2=16}$. We
see that for $\delta/J\\!\sim\\!1$, the twin Fock state is projected onto many
eigenstates of the system, while it indeed becomes an eigenstate for
$\delta/J\gg 1$. As $\delta$ increases, the MF fixed-point of lowest energy
migrates towards the South pole $z=-1$ [Fig 1(h-i)]. The remaining
trajectories wraps around the Bloch sphere, as do the eigenstates, turning
into Fock states of definite relative population $\langle\hat{J}_{z}\rangle$.

From this brief study, summarized in Fig. 1, one apprehends how difficult it
is to adiabatically connect the eigenstates of panels (h) to (j), hence
suggesting the development of a more sophisticated approach.

#### Optimal control.

We now consider optimal control as a promising alternative to prepare twin
Fock states $\ket{N/2}$, through an arbitrary variation of the tilt
$\delta(t)$ between the wells. To obtain the optimal $\delta(t)$ that steers
the ground state $\ket{\phi_{0}}$ of (2) (taken at $\delta=0$) towards
$\ket{N/2}$ in a control duration $t_{\text{f}}$ (fixed beforehand, see
below), we employ a first-order gradient-based algorithm [60, 61, 62, 63, 64,
48]. This approach amounts to performing a gradient ascent on a state-
preparation metric, which we choose here to be the fidelity $\mathcal{F}$ to
the target Fock state $\ket{N/2}$ [50]. The quality of the optimal control
fields (in terms of e.g. fidelity maximization, robustness to noise and
experimental feasibility) yielded by such numerical approaches strongly
depends on the initial guess provided as input. For our purposes, we find that
a rather large, constant tilt $\delta_{0}/J\sim\Lambda/2=NU/4J$ is an
efficacious initial guess, as it amounts to a quench into a configuration
where both the initial state and the target Fock state are projected into the
same restricted subspace consisting of only a few quasi-Fock eigenstates [see
the blue states in Fig. 1(g)], effectively reducing the Hilbert space
dimension into which the optimization takes place. In the large $\delta_{0}$
limit, the spacing between two successive Fock eigenstates
$\ket{n_{1}},\ket{n_{1}+1}$ is $\Delta E\sim U(2n_{1}-N+1)+2\delta_{0}\sim
2\delta_{0}$ for $n_{1}=N/2$ [50]. The period $T_{0}$ at which a non-steady
state evolves in the vicinity of $\ket{N/2}$ in Hilbert space is therefore
$T_{0}\sim h/2\delta_{0}$; this establishes a lower bound on the time
associated with state alteration near the equator of the Bloch sphere in the
strongly tilted system. Here, we typically opt for a control duration
$t_{\text{f}}\sim~{}10\,T_{0}\sim 20\,h/NU$ 222Working with finite
$\delta_{0}$, we roughly adjust $t_{\text{f}}$ and $\delta_{0}$ around the
asymptotic expressions given in the text in order to maximize the fidelity
resulting from the initial guess.. This favorable scaling with $N$ stems from
the large tilt $\delta_{0}\sim NU/4$ used in our approach. In a tight-binding
framework, relevant to optical lattices [66], this would ultimately limit the
maximum number of bosons $N$ that one can consider.

Figure 2:  Twin Fock state preparation by optimal control. (a) Optimal control
tilt as a function of time (solid black line) to drive the ground state of Eq.
(2) with $U/J=1$ and $N=32$ (taken at $\delta/J=0$) towards the twin Fock
state $\ket{N/2=16}$, and initial guess $\delta(t)=\delta_{0}$ (dashed grey
line). Inset: norm of the Fourier transform of $\delta(t)$ (solid black line,
$y$-axis going from 0 to $1J$) and $\pm 1/T_{0}$ (dashed red lines, see text).
(b) Fidelity $\mathcal{F}$ to $\ket{N/2=16}$ (blue solid line), fragmentation
entropy $S_{\text{F}}$ (divided by $\ln(2)$, orange dashed) and number
squeezing parameter $\xi^{2}_{N}$ (dB scale, black dash-dotted) as a function
of time. See Table 1 for details. (c) Fock space projection during the
preparation.

Numerical results for the preparation of the twin Fock state, using optimal
control of the tilt, are shown in Fig. 2 for $N=32$ and $U/J=1$. Panel (a)
shows the optimal $\delta(t)$, displaying a fast, non-trivial oscillation
around $\delta_{0}$. This is a general trend observed for different $N$
values, and a Fourier analysis of $\delta(t)$ [inset of Fig. 2(a)] reveals
that its main non-zero frequencies are grouped around the transition frequency
$\Delta E/h=1/T_{0}$ between two adjacent equatorial Fock states in the large
$\delta_{0}$ limit. The algorithm thus naturally converges on a quasi-periodic
modulation that hybridizes neighboring quasi-Fock eigenstates, in order to
then maximize the projection onto the targeted Fock state. The resulting
fidelity $\mathcal{F}$, fragmentation entropy $S_{\text{F}}$ and number
squeezing $\xi_{N}^{2}$ over the preparation are plotted on Fig. 2(b), while
(c) features the Fock-space projection of the prepared state. For those
parameters, our metrics finally reach $\mathcal{F}=0.999$ (breaking condition
of the iterative algorithm), $S_{\text{F}}/\ln(2)=0.9998$ and
$\xi_{N}^{2}=1.33\cdot 10^{-3}=-28.8$ dB, i.e. number-difference fluctuations
suppressed by nearly 3 orders of magnitude compared with coherent spin states
lying on the equator of the Bloch sphere [20]. One observes how, for this
strongly non-adiabatic scheme, the fidelity to the most fragmented state
$\ket{N/2}$ is a much more fragile figure of merit for fragmentation than the
fragmentation entropy. This can be appreciated by considering, for instance,
that the very next Fock states $\ket{N/2\pm 1}$ have high fragmentation
entropy but zero fidelity to $\ket{N/2}$. Results from $N=4$ to 512 are
compiled in Table 1. Apart from the fidelity, which decreases appreciably as
$N$ increases, we see that optimal control allows to reach high fragmentation
entropy and number squeezing up to substantially large $N$. It should be noted
that, as $U/J$ is fixed in this study, the ground state of the unbiased
Hamiltonian (2) gets inherently number squeezed as $N$ increases 333As the
mean-field interaction parameter $\Lambda\propto NU/J$, increasing $N$ at
fixed $U/J$ gets us deeper into the Josephson regime [57, 58, 20] (see Table
1; the number-squeezing gain due to optimal control being the difference
between the last two rows). Overall, we see that optimal control combined with
a strong average tilt $\delta_{0}$ is a promising method to produce the
maximally fragmented Fock states $\ket{N/2}$ in a two-mode Bose-Hubbard
system; see [50] for similar results on the $L\\!=\\!3$ case. Moreover, one
achieves drastic number squeezing, two to three orders of magnitudes faster
than (quasi-)adiabatic methods; see Refs. [25, 26, 27, 28, 31, 33, 34] and
[50].

$N$ |  4 |  8 |  16 |  32 |  64 |  128 |  256 |  512
---|---|---|---|---|---|---|---|---
$\delta_{0}(J)$ |  2.80 |  4.40 |  7.20 |  12.8 |  24.0 |  41.6 |  64.0 |  89.6
$t_{\text{f}}\,(\hbar/J)$ |  5.05 |  3.43 |  2.77 |  2.02 |  1.47 |  0.925 |  0.702 |  0.500
$\mathcal{F}$ |  $0.999$ |  0.987 |  0.951 |  0.923 |  0.888
$S_{\text{F}}/\ln(2)$ |  $>0.999$ |  0.994 |  0.980 |  0.959 |  0.950
$\begin{array}[]{c}\xi^{2}_{N}\\\ \text{(dB)}\end{array}\left(\begin{array}[]{l}\psi_{\text{f}}\\\ \phi_{0}\end{array}\right)$ |  -26.0 |  -22.6 |  -21.6 |  -28.8 |  -19.0 |  -14.7 |  -15.2 |  -19.3
|  -2.30 |  -3.52 |  -4.84 |  -6.23 |  -7.66 |  -9.12 |  -10.6 |  -12.1
Table 1:  Twin Fock state preparation by optimal control ($U/J=1$) for
increasing number of bosons $N$, with guess tilt $\delta_{0}$, control
duration $t_{\text{f}}$, final fidelity $\mathcal{F}$ to $\ket{N/2}$,
fragmentation entropy $S_{\text{F}}$ and number squeezing parameter of
prepared and initial states ($\ket{\psi_{\text{f}}}$ and $\ket{\phi_{0}}$
resp.).

The full quantum dynamics depends on the details of the instantaneous many-
body spectra, over the entire control duration. A priori, a limitation of this
approach, as compared to adiabatic schemes, is that it requires a precise
knowledge of the number of bosons $N$. Quite generically, optimal control
lacks robustness against parameter fluctuations, especially those associated
with the particle number $N$ [49]. However, we find that in this large-
average-tilt configuration (where the eigenenergy spacing between successive
quasi-Fock eigenstates loses its dependence on the average number of particles
in said eigenstates), an optimal $\delta(t)$ computed for a given $N$ remains
efficient to fragment the Bose gas for a close but different $N^{\prime}=N+2k$
(keeping commensurability with $k\in\mathbb{Z}^{*}$), and can be used as a
solid initial guess for another optimization with the correct $N^{\prime}$
[50].

In Fig. 3, we further interpret the algorithm’s “strategy” in performing
number squeezing around $\ket{N/2}$ in this biased two-mode system by drawing
a parallel with spin-squeezing performed by one-axis twisting (OAT) [68, 69,
20]. Through OAT, a coherent spin state initially pointing along $\vec{n}$ on
the Bloch sphere’s equator undergoes spin-squeezing through an evolution
governed by the interaction Hamiltonian $U\hat{J}^{2}_{z}$, due to MF
trajectories flowing in opposite directions on each Bloch hemisphere
($\differential\varphi/\differential t\propto z$). The resulting spin
squeezing, occurring along a tangential direction $\vec{s}$ to the average
spin direction $\vec{n}$, can be quantified by the squeezing parameter
$\xi_{\text{R}}^{2}=N(\Delta\hat{J}_{\vec{s}})^{2}/\langle\hat{J}_{\vec{n}}\rangle^{2}$
[70, 20], assessing the metrological usefulness of the anisotropy of quantum
fluctuations around $\vec{n}$. OAT spin-squeezes an equatorial coherent spin
state for $t\lesssim N^{-1/2}\hbar/U$, while for longer times, it gets over-
squeezed as it wraps around the Bloch sphere, resulting in
$\xi_{\text{R}}^{2}>1$ [68, 20] [Fig. 3(a,c)]. A salient symptom of over-
squeezing is the constellation of zeros riddling the Husimi function [Fig.
3(a)], which are Majorana anti-stars, antipodes of the $N$ individual
$1/2$-spins into which the state can be factored according to the Majorana
stellar representation [71, 72, 73, 74]. However, OAT never number squeezes,
having Fock states as its eigenstates [Fig. 3(c)]. In the case of optimal
control [Fig. 3(b)], the large tilt yields MF trajectories similar to OAT near
the equator, resulting in a slight short-term spin squeezing (panel (c)). The
control on $\delta(t)$ then constrains the state towards the equator while it
wraps around the Bloch sphere, gradually sending the Majorana anti-stars to
its poles. Our optimal control approach, operating around a large tilt value,
can thus be viewed as number squeezing through optimized axis twisting.

Figure 3:  Optimized axis twisting ($N=32$). (a,b) Husimi state
representations (colored area), their zeros (Majorana anti-stars, red crosses)
and instantaneous MF trajectories (solid black lines) on the Bloch sphere. (a)
correspond to standard one-axis twisting of a coherent spin state, with
initial condition (left), maximum spin-squeezing time (middle) and
oversqueezing at a longer time (right). (b) is the optimal preparation of
$\ket{N/2}$ featured in Fig. 2 for three different increasing times (left to
right). (c) Number squeezing $\xi^{2}_{N}$ (solid lines) and spin squeezing
$\xi^{2}_{\text{R}}$ (dashed lines) of (a,b) (resp. orange, blue) over the
protocols (of independent duration $t_{\text{f}}=2N^{-1/2}\hbar/U$ (a) and
$t_{\text{f}}=2.02\,\hbar/J$ (b)).

#### Conclusion.

In this work, we introduced a practical and efficient approach to achieve
total fragmentation of a Bose gas in the $L$-mode Bose-Hubbard model, i.e. to
prepare the evenly distributed $L$-tuple Fock state $\ket{N/L,...,N/L}$ with
$N$ bosons. We illustrated our results for systems with $L\\!=\\!2$ and 3
modes. We first considered the conventional adiabatic ramping of the ratio
between interaction $U$ and hopping $J$ (either directly varied, or
effectively via Floquet engineering), finding strong limitations and
constraints: long preparation time, low particle number and high Floquet
frequency. We then applied quantum optimal control with the energy bias
$\delta$ between the modes as our control parameter. We found that by
operating around a large tilt $\delta\approx NU/4$, one is able to
diabatically prepare the twin Fock state in a duration $t_{\text{f}}\approx
20\,\hbar/NU$, with appreciable robustness against atom number fluctuations
[50]. This work offers realistic approaches to experimentally achieve
$L$-tuple Fock state preparation, the most resourceful number state for
metrological purposes [20]. We emphasize that our scheme relies on modulating
the tilt (or detuning between the modes) around a large mean value. This can
be advantageous for certain physical settings that cannot operate close to
resonance (e.g. schemes using two coupled internal states operating in the
large-detuning regime).

#### Acknowledgements.

The authors thank Jean Dalibard for insightful discussions, and for his
careful reading of the manuscript. We also acknowledge discussions with Markus
Oberthaler, Helmut Strobel, Maximilian Prüfer, Leon Zaporski and Shai Tsesses,
as well as support from the ERC (LATIS project) and the EOS (CHEQS project).

## References

  * Cornell and Wieman [2002] E. A. Cornell and C. E. Wieman, Nobel Lecture: Bose-Einstein condensation in a dilute gas, the first 70 years and some recent experiments, Rev. Mod. Phys. 74, 875 (2002).
  * Ketterle [2002] W. Ketterle, Nobel lecture: When atoms behave as waves: Bose-Einstein condensation and the atom laser, Rev. Mod. Phys. 74, 1131 (2002).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Nozières and Saint James [1982] P. Nozières and D. Saint James, Particle vs. Pair Condensation in Attractive Bose Liquids, J. Phys. France 43, 1133 (1982).
  * Spekkens and Sipe [1999] R. W. Spekkens and J. E. Sipe, Spatial fragmentation of a Bose-Einstein condensate in a double-well potential, Phys. Rev. A 59, 3868 (1999).
  * Mueller  _et al._ [2006] E. J. Mueller, T.-L. Ho, M. Ueda, and G. Baym, Fragmentation of Bose-Einstein condensates, Phys. Rev. A 74, 033612 (2006).
  * Beekman  _et al._ [2019] A. Beekman, L. Rademaker, and J. van Wezel, An Introduction to Spontaneous Symmetry Breaking, SciPost Physics Lecture Notes , 011 (2019).
  * Evrard  _et al._ [2021] B. Evrard, A. Qu, J. Dalibard, and F. Gerbier, Observation of fragmentation of a spinor Bose-Einstein condensate, Science 373, 1340 (2021).
  * Gerbier  _et al._ [2005] F. Gerbier, A. Widera, S. Fölling, O. Mandel, T. Gericke, and I. Bloch, Phase Coherence of an Atomic Mott Insulator, Phys. Rev. Lett. 95, 050404 (2005).
  * Jaksch  _et al._ [1998] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Cold Bosonic Atoms in Optical Lattices, Phys. Rev. Lett. 81, 3108 (1998).
  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature 415, 39 (2002).
  * Widera  _et al._ [2005] A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, Coherent Collisional Spin Dynamics in Optical Lattices, Phys. Rev. Lett. 95, 190405 (2005).
  * Gerbier  _et al._ [2006] F. Gerbier, S. Fölling, A. Widera, O. Mandel, and I. Bloch, Probing Number Squeezing of Ultracold Atoms across the Superfluid-Mott Insulator Transition, Phys. Rev. Lett. 96, 090401 (2006).
  * Lepoutre  _et al._ [2019] S. Lepoutre, J. Schachenmayer, L. Gabardos, B. Zhu, B. Naylor, E. Maréchal, O. Gorceix, A. M. Rey, L. Vernac, and B. Laburthe-Tolra, Out-of-equilibrium quantum magnetism and thermalization in a spin-3 many-body dipolar lattice system, Nature Communications 10, 1714 (2019).
  * Holland and Burnett [1993] M. J. Holland and K. Burnett, Interferometric detection of optical phase shifts at the Heisenberg limit, Phys. Rev. Lett. 71, 1355 (1993).
  * Kim  _et al._ [1998] T. Kim, O. Pfister, M. J. Holland, J. Noh, and J. L. Hall, Influence of decorrelation on Heisenberg-limited interferometry with quantum correlated photons, Phys. Rev. A 57, 4004 (1998).
  * Huang and Moore [2008] Y. P. Huang and M. G. Moore, Optimized Double-Well Quantum Interferometry with Gaussian Squeezed States, Phys. Rev. Lett. 100, 250406 (2008).
  * Lücke  _et al._ [2011] B. Lücke, M. Scherer, J. Kruse, L. Pezzé, F. Deuretzbacher, P. Hyllus, O. Topic, J. Peise, W. Ertmer, J. Arlt, L. Santos, A. Smerzi, and C. Klempt, Twin Matter Waves for Interferometry Beyond the Classical Limit, Science 334, 773 (2011).
  * Zou  _et al._ [2018] Y.-Q. Zou, L.-N. Wu, Q. Liu, X.-Y. Luo, S.-F. Guo, J.-H. Cao, M. K. Tey, and L. You, Beating the classical precision limit with spin-1 Dicke states of more than 10,000 atoms, Proceedings of the National Academy of Sciences 115, 6381 (2018).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).
  * Barbieri [2022] M. Barbieri, Optical Quantum Metrology, PRX Quantum 3, 010202 (2022).
  * Kuzmich  _et al._ [1998] A. Kuzmich, N. P. Bigelow, and L. Mandel, Atomic quantum non-demolition measurements and squeezing, Europhysics Letters 42, 481 (1998).
  * Mølmer [1999] K. Mølmer, Twin-correlations in atoms, The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics 5, 301 (1999).
  * Bücker  _et al._ [2011] R. Bücker, J. Grond, S. Manz, T. Berrada, T. Betz, C. Koller, U. Hohenester, T. Schumm, A. Perrin, and J. Schmiedmayer, Twin-atom beams, Nature Physics 7, 608 (2011).
  * Orzel  _et al._ [2001] C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, and M. A. Kasevich, Squeezed States in a Bose-Einstein Condensate, Science 291, 2386 (2001).
  * Li  _et al._ [2007] W. Li, A. K. Tuchman, H.-C. Chien, and M. A. Kasevich, Extended Coherence Time with Atom-Number Squeezed States, Phys. Rev. Lett. 98, 040402 (2007).
  * Jo  _et al._ [2007] G.-B. Jo, Y. Shin, S. Will, T. A. Pasquini, M. Saba, W. Ketterle, D. E. Pritchard, M. Vengalattore, and M. Prentiss, Long Phase Coherence Time and Number Squeezing of Two Bose-Einstein Condensates on an Atom Chip, Phys. Rev. Lett. 98, 030407 (2007).
  * Sebby-Strabley  _et al._ [2007] J. Sebby-Strabley, B. L. Brown, M. Anderlini, P. J. Lee, W. D. Phillips, J. V. Porto, and P. R. Johnson, Preparing and Probing Atomic Number States with an Atom Interferometer, Phys. Rev. Lett. 98, 200405 (2007).
  * Rodríguez  _et al._ [2007] M. Rodríguez, S. R. Clark, and D. Jaksch, Generation of twin Fock states via transition from a two-component Mott insulator to a superfluid, Phys. Rev. A 75, 011601 (2007).
  * Rodríguez  _et al._ [2008] M. Rodríguez, S. R. Clark, and D. Jaksch, Adiabatic melting of two-component Mott-insulator states, Phys. Rev. A 77, 043613 (2008).
  * Estève  _et al._ [2008] J. Estève, C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler, Squeezing and entanglement in a Bose-Einstein condensate, Nature 455, 1216 (2008).
  * Maussang  _et al._ [2010] K. Maussang, G. E. Marti, T. Schneider, P. Treutlein, Y. Li, A. Sinatra, R. Long, J. Estève, and J. Reichel, Enhanced and Reduced Atom Number Fluctuations in a BEC Splitter, Phys. Rev. Lett. 105, 080403 (2010).
  * Zhang and Duan [2013] Z. Zhang and L.-M. Duan, Generation of Massive Entanglement through an Adiabatic Quantum Phase Transition in a Spinor Condensate, Phys. Rev. Lett. 111, 180401 (2013).
  * Luo  _et al._ [2017] X.-Y. Luo, Y.-Q. Zou, L.-N. Wu, Q. Liu, M.-F. Han, M. K. Tey, and L. You, Deterministic entanglement generation from driving through quantum phase transitions, Science 355, 620 (2017).
  * Zhang  _et al._ [2024] T. Zhang, M. Maiwöger, F. Borselli, Y. Kuriatnikov, J. Schmiedmayer, and M. Prüfer, Squeezing Oscillations in a Multimode Bosonic Josephson Junction, Phys. Rev. X 14, 011049 (2024).
  * Juliá-Díaz  _et al._ [2012] B. Juliá-Díaz, E. Torrontegui, J. Martorell, J. G. Muga, and A. Polls, Fast generation of spin-squeezed states in bosonic Josephson junctions, Phys. Rev. A 86, 063623 (2012).
  * Yuste  _et al._ [2013] A. Yuste, B. Juliá-Díaz, E. Torrontegui, J. Martorell, J. G. Muga, and A. Polls, Shortcut to adiabaticity in internal bosonic Josephson junctions, Phys. Rev. A 88, 043647 (2013).
  * Grond  _et al._ [2009] J. Grond, J. Schmiedmayer, and U. Hohenester, Optimizing number squeezing when splitting a mesoscopic condensate, Phys. Rev. A 79, 021603 (2009).
  * Pichler  _et al._ [2016] T. Pichler, T. Caneva, S. Montangero, M. D. Lukin, and T. Calarco, Noise-resistant optimal spin squeezing via quantum control, Phys. Rev. A 93, 013851 (2016).
  * Guo  _et al._ [2021] S.-F. Guo, F. Chen, Q. Liu, M. Xue, J.-J. Chen, J.-H. Cao, T.-W. Mao, M. K. Tey, and L. You, Faster State Preparation across Quantum Phase Transition Assisted by Reinforcement Learning, Phys. Rev. Lett. 126, 060401 (2021).
  * Dunningham  _et al._ [2001] J. A. Dunningham, K. Burnett, and M. Edwards, Relative number squeezing in Bose-Einstein condensates, Phys. Rev. A 64, 015601 (2001).
  * Ebert  _et al._ [2014] M. Ebert, A. Gill, M. Gibbons, X. Zhang, M. Saffman, and T. G. Walker, Atomic fock state preparation using rydberg blockade, Phys. Rev. Lett. 112, 043602 (2014).
  * Xin  _et al._ [2023] L. Xin, M. Barrios, J. T. Cohen, and M. S. Chapman, Long-Lived Squeezed Ground States in a Quantum Spin Ensemble, Phys. Rev. Lett. 131, 133402 (2023).
  * Eckardt  _et al._ [2005a] A. Eckardt, C. Weiss, and M. Holthaus, Superfluid-Insulator Transition in a Periodically Driven Optical Lattice, Phys. Rev. Lett. 95, 260404 (2005a).
  * Lignier  _et al._ [2007] H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zenesini, O. Morsch, and E. Arimondo, Dynamical Control of Matter-Wave Tunneling in Periodic Potentials, Phys. Rev. Lett. 99, 220403 (2007).
  * Glaser  _et al._ [2015] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Training Schrödinger’s cat: quantum optimal control, The European Physical Journal D 69, 279 (2015).
  * Koch  _et al._ [2022] C. P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S. J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Quantum optimal control in quantum technologies. Strategic report on current status, visions and goals for research in Europe, EPJ Quantum Technology 9, 19 (2022).
  * Ansel  _et al._ [2024] Q. Ansel, E. Dionis, F. Arrouas, B. Peaudecerf, S. Guérin, D. Guéry-Odelin, and D. Sugny, Introduction to Theoretical and Experimental aspects of Quantum Optimal Control (2024), arXiv:2403.00532 [quant-ph] .
  * Lapert  _et al._ [2012] M. Lapert, G. Ferrini, and D. Sugny, Optimal control of quantum superpositions in a bosonic Josephson junction, Phys. Rev. A 85, 023611 (2012).
  * Note [1] See Supplementary Material for details and additional results in $L\\!=\\!2$ and $L\\!=\\!3$.
  * Husimi [1940] K. Husimi, Some Formal Properties of the Density Matrix, Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 22, 264 (1940).
  * Agarwal [1981] G. S. Agarwal, Relation between atomic coherent-state representation, state multipoles, and generalized phase-space distributions, Phys. Rev. A 24, 2889 (1981).
  * Haas  _et al._ [2014] F. Haas, J. Volz, R. Gehr, J. Reichel, and J. Estève, Entangled States of More Than 40 Atoms in an Optical Fiber Cavity, Science 344, 180 (2014).
  * Smerzi  _et al._ [1997] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Quantum Coherent Atomic Tunneling between Two Trapped Bose-Einstein Condensates, Phys. Rev. Lett. 79, 4950 (1997).
  * Goldman  _et al._ [2023] N. Goldman, O. Diessel, L. Barbiero, M. Prüfer, M. Di Liberto, and L. Peralta Gavensky, Floquet-Engineered Nonlinearities and Controllable Pair-Hopping Processes: From Optical Kerr Cavities to Correlated Quantum Matter, PRX Quantum 4, 040327 (2023).
  * Raventós  _et al._ [2017] D. Raventós, T. Graß, M. Lewenstein, and B. Juliá-Díaz, Cold Bosons in Optical Lattices: A Tutorial for Exact Diagonalization, J. Phys. B: At. Mol. Opt. Phys. 50, 113001 (2017).
  * Leggett [2001] A. J. Leggett, Bose-Einstein condensation in the alkali gases: Some fundamental concepts, Rev. Mod. Phys. 73, 307 (2001).
  * Gati and Oberthaler [2007] R. Gati and M. K. Oberthaler, A bosonic Josephson junction, Journal of Physics B: Atomic, Molecular and Optical Physics 40, R61 (2007).
  * Schumm  _et al._ [2005] T. Schumm, S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer, and P. Krüger, Matter-wave interferometry in a double well on an atom chip, Nature Physics 1, 57 (2005).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172, 296 (2005).
  * Werschnik and Gross [2007] J. Werschnik and E. K. U. Gross, Quantum optimal control theory, Journal of Physics B: Atomic, Molecular and Optical Physics 40, R175 (2007).
  * Boscain  _et al._ [2021] U. Boscain, M. Sigalotti, and D. Sugny, Introduction to the Pontryagin Maximum Principle for Quantum Optimal Control, PRX Quantum 2, 030203 (2021).
  * Dupont  _et al._ [2021] N. Dupont, G. Chatelain, L. Gabardos, M. Arnal, J. Billy, B. Peaudecerf, D. Sugny, and D. Guéry-Odelin, Quantum State Control of a Bose-Einstein Condensate in an Optical Lattice, PRX Quantum 2, 040303 (2021).
  * Dupont  _et al._ [2023] N. Dupont, F. Arrouas, L. Gabardos, N. Ombredane, J. Billy, B. Peaudecerf, D. Sugny, and D. Guéry-Odelin, Phase-space distributions of Bose–Einstein condensates in an optical lattice: optimal shaping and reconstruction, New Journal of Physics 25, 013012 (2023).
  * Note [2] Working with finite $\delta_{0}$, we roughly adjust $t_{\text{f}}$ and $\delta_{0}$ around the asymptotic expressions given in the text in order to maximize the fidelity resulting from the initial guess.
  * Eckardt [2017] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017).
  * Note [3] As the mean-field interaction parameter $\Lambda\propto NU/J$, increasing $N$ at fixed $U/J$ gets us deeper into the Josephson regime [57, 58, 20].
  * Kitagawa and Ueda [1993] M. Kitagawa and M. Ueda, Squeezed spin states, Phys. Rev. A 47, 5138 (1993).
  * Micheli  _et al._ [2003] A. Micheli, D. Jaksch, J. I. Cirac, and P. Zoller, Many-particle entanglement in two-component Bose-Einstein condensates, Phys. Rev. A 67, 013607 (2003).
  * Wineland  _et al._ [1994] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. Heinzen, Squeezed atomic states and projection noise in spectroscopy, Phys. Rev. A 50, 67 (1994).
  * Majorana [1932] E. Majorana, Atomi orientati in campo magnetico variabile, Il Nuovo Cimento (1924-1942) 9, 43 (1932).
  * Bloch and Rabi [1945] F. Bloch and I. I. Rabi, Atoms in Variable Magnetic Fields, Rev. Mod. Phys. 17, 237 (1945).
  * Bruno [2012] P. Bruno, Quantum Geometric Phase in Majorana’s Stellar Representation: Mapping onto a Many-Body Aharonov-Bohm Phase, Phys. Rev. Lett. 108, 240402 (2012).
  * Evrard  _et al._ [2019] A. Evrard, V. Makhalov, T. Chalopin, L. A. Sidorenkov, J. Dalibard, R. Lopes, and S. Nascimbene, Enhanced Magnetic Sensitivity with Non-Gaussian Quantum Fluctuations, Phys. Rev. Lett. 122, 173601 (2019).
  * Rahav  _et al._ [2003] S. Rahav, I. Gilary, and S. Fishman, Effective Hamiltonians for periodically driven systems, Phys. Rev. A 68, 013820 (2003).
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Periodically Driven Quantum Systems: Effective Hamiltonians and Engineered Gauge Fields, Phys. Rev. X 4, 031027 (2014).
  * Eckardt and Anisimovas [2015] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a Floquet-space perspective, New Journal of Physics 17, 093039 (2015).
  * Novičenko  _et al._ [2017] V. Novičenko, E. Anisimovas, and G. Juzeliūnas, Floquet analysis of a quantum system with modulated periodic driving, Phys. Rev. A 95, 023615 (2017).
  * Eckardt  _et al._ [2005b] A. Eckardt, T. Jinasundera, C. Weiss, and M. Holthaus, Analog of Photon-Assisted Tunneling in a Bose-Einstein Condensate, Phys. Rev. Lett. 95, 200401 (2005b).
  * Pu  _et al._ [2000] H. Pu, S. Raghavan, and N. Bigelow, Manipulating spinor condensates with magnetic fields: Stochastization, metastability, and dynamical spin localization, Phys. Rev. A 61, 023602 (2000).
  * Pontryagin [1987] L. S. Pontryagin,  _Mathematical Theory of Optimal Processes_ (CRC Press, 1987).
  * Bellman [1966] R. Bellman, Dynamic Programming, Science 153, 34 (1966).
  * Sussmann and Willems [1997] H. J. Sussmann and J. C. Willems, 300 Years of Optimal Control: From The Brachystochrone to the Maximum Principle, IEEE Control Systems Magazie , 32 (1997).
  * Liberzon [2012] D. Liberzon,  _Calculus of Variations and Optimal Control Theory_ (Princeton University Press, 2012).
  * Dupont [2022] N. Dupont,  _Control and transport of matter waves in an optical lattice_ , Ph.D. thesis, Université Toulouse 3 - Paul Sabatier (2022).
  * Note [4] For instance with $L=2$ and $N=128$, the first optimization was performed in a Hilbert space of dimension 19, only considering the basis Fock states $\ket{n_{1}}$ with $55\leq n_{1}\leq 73$.

Supplementary Material

##  S1 A. Fragmentation from adiabatic ground-state following

We consider the fragmentation of the many-body ground state in the $L$-mode
Bose-Hubbard Hamiltonian

$\hat{H}=-J\sum_{\ell=1}^{L-1}\left(\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell+1}+\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}\right)+\dfrac{U}{2}\sum_{\ell=1}^{L}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell},$
(S1)

as one tunes the ratio $U/J$; see Fig. 1 of the main text. In this work, we
focus on the two-mode ($L=2$) and three-mode ($L=3$) cases.

Figure S1:  Adiabatic ground state following as $U/J$ increases. (a) Energy
gap between the ground state and the first accessible excited state as a
function of $U/J\in[0,(U/J)_{\text{f}}]$ (see text) and $N\\!\in\\![2,32]$
(only considering even integers for commensurability). (b) Final fidelity
$\mathcal{F}_{0}=|\langle\phi^{(U/J)_{\text{f}}}_{0}|\psi(t_{\text{f}})\rangle|^{2}$
(see text) after a linear ramping of $U/J$ from 0 to $(U/J)_{\text{f}}$ and as
a function of the ramping rate $1/t_{\text{f}}$ and $N$. The initial state is
$|\phi^{(0)}_{0}\rangle$.

First for $L\\!=\\!2$, one can rewrite the Hamiltonian
$\hat{H}=U\hat{J}^{2}_{z}-2J\hat{J}_{x}$ with
$\hat{J}_{x}=(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2$
and
$\hat{J}_{z}=(\hat{a}_{1}^{\dagger}\hat{a}_{1}-\hat{a}_{2}^{\dagger}\hat{a}_{2})/2$.
Considering $N$ bosons, we denote $E^{(U/J)}_{n}$ and
$|\phi^{(U/J)}_{n}\rangle$ the $N+1$ eigenenergies and eigenstates of (S1) as
a function of $U/J$. Starting from $U/J\ll 1$, the hopping term dominates, and
we have

$\displaystyle-2J\hat{J}_{x}$
$\displaystyle=-J(\hat{a}^{\dagger}_{1}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}_{1})$
$\displaystyle=-J(\hat{a}^{\dagger}_{+}\hat{a}_{+}-\hat{a}^{\dagger}_{-}\hat{a}_{-}),$
(S2)

with $\hat{a}_{\pm}=(\hat{a}_{1}\pm\hat{a}_{2})/\sqrt{2}$. In the limit $U=0$,
the eigenstates are thus of the form

$\ket{\phi^{(0)}_{n}}=\dfrac{\left(\hat{a}_{+}^{\dagger}\right)^{N-n}}{\sqrt{(N-n)!}}\dfrac{\left(\hat{a}_{-}^{\dagger}\right)^{n}}{\sqrt{(n)!}}\ket{\varnothing},$
(S3)

with $E_{n}^{(0)}=J(2n-N)$ and where $\ket{\varnothing}$ is the vacuum state.
The total Hamiltonian (hopping terms and interactions) is invariant under the
swapping of modes $1\leftrightarrow 2$, which is associated with the parity
operator $\hat{P}$ acting as $\hat{P}\hat{a}_{\pm}=\pm\hat{a}_{\pm}$. From Eq.
(S3), we see that the eigenstates of $\hat{J}_{x}$ are directly sorted in
subspaces of opposite parity, through the parity of $n$:
$\hat{P}|\phi^{(U/J)}_{n}\rangle=(-1)^{n}|\phi^{(U/J)}_{n}\rangle$. Starting
from $|\phi^{(0)}_{0}\rangle$ and given that $\hat{H}$ does not couple the two
subspaces of opposite parity $\forall\,U/J$, the first possible excitation at
$U/J=0$ is towards $|\phi^{(0)}_{2}\rangle$, with energy $\Delta
E^{(0)}=E^{(0)}_{2}-E^{(0)}_{0}=4J$ (Fig. S1(a)). This gap enlarges as $U/J$
increases, up to $\Delta E\sim U$ between $\ket{\phi_{0}}=\ket{N/2}$ and
$\ket{\phi_{2}}=(\ket{N/2+1}+\ket{N/2-1})/\sqrt{2}$ in the strongly
interacting regime $U/J\gg 1$ (with the notation
$\ket{n_{1}}\equiv\ket{n_{1}}_{1}\otimes\ket{N-n_{1}}_{2}$).

Figure S2:  Fragmentation in a three-mode Bose-Hubbard system. (a) Many-body
spectrum for $N\\!=\\!9$ bosons as a function of $U/J$ (grey lines), with
$E_{n}$ the eigenenergy of eigenstate $\ket{\phi_{n}}$. The projection of
$\ket{N/3}$ onto $\ket{\phi_{n}}$ is encoded in shades of blue on $E_{n}$.
(b-d) Fock state projection of the states identified in (a) (blue areas, with
independent maximum values $v_{\text{max}}=0.0514$ (b), 0.177 (c) and 0.999
(d)). (e-h) same as (a-d) as a function of the tilt $\delta/J$ with
$U/J\\!=\\!1$, and $v_{\text{max}}=0.177$ (b), 0.592 (c) and 0.884 (d).

The maximally fragmented state $\ket{N/2}$ can be reached from
$|\phi^{(0)}_{0}\rangle$ through adiabatic ground-state following by slowly
ramping-up $U/J$. The adiabaticity condition on the ramping time can be
estimated in first approximation by the inverse of the frequency associated
with the narrowest gap in the spectrum, i.e. $t_{\text{f}}\gg h/4J$. To go
beyond this estimate, we perform simulation of adiabatic ground state
following as a function of the ramping time $t_{\text{f}}$ and $N$. Starting
from $|\phi^{(0)}_{0}\rangle$, we linearly ramp up $U/J$ until a value
$(U/J)_{\text{f}}$ (increasing with $N$) such that $|\langle
N/2|\phi_{0}^{(U/J)_{\text{f}}}\rangle|^{2}\approx 0.99$. Figure S1(b) shows
the ground state fidelity
$\mathcal{F}_{0}=|\langle\phi^{(U/J)_{\text{f}}}_{0}|\psi(t_{\text{f}})\rangle|^{2}$
of the final state as a function of $N$ and $t_{\text{f}}$. While the
narrowest relevant gap in the spectrum is $4\,J$ $\forall\,N$, we see that the
minimum $t_{\text{f}}$ required to reach a satisfactory fidelity
$\mathcal{F}_{0}$ increases with $N$. This is due to the increasing number of
levels in the spectrum that provide additional channels towards which the
evolved state can propagate. We see that the lowest ramping time to get
$\mathcal{F}_{0}\approx 1$ is already of the order or $10^{3}\hbar/J$ for
$N\\!=32$.

Regarding the $L\\!=\\!3$ case: we display in Fig. S2 a plot analogous to Fig.
1 (main text for $L=2$). Similarly to the case $L\\!=\\!2$, adiabatically
increasing $U/J\rightarrow\infty$ allows one to reach the maximally fragmented
state $\ket{N/3}\equiv\ket{N/3,N/3,N/3}$ from the ground state of the system
with finite interactions (panels (a-d). There, the scaling of the Hilbert
space dimension $d_{\text{H}}$ with $N$ is quadratic:
$d_{\text{H}}=(N+1)(N+2)/2$, which, as discussed previously, is expected to
hinder adiabaticity, even more than for $L\\!=\\!2$. In panels (e-h), we show
that, as for $L\\!=\\!2$, an additional constant tilt reveals the target state
in the spectrum. In Sec. D, we utilize this tilt as a parameter to prepare
$\ket{N/3}$ by optimal control.

##  S2 B. Floquet control of $U/J$

One can introduce a fast periodic modulation of the tilt around $\delta=0$ to
dynamically control the effective tunneling between the two modes [44]. In
this way, one can drive a superfluid to Mott insulator transition [45] by
slowly ramping down the effective tunneling of the particles.

More precisely, we consider a strong sinusoidal modulation of the tilt
$\delta(t)=A\cos(\omega t)$ and the Hamiltonian of the periodically driven
system is written as

$\displaystyle\hat{H}(t)=$
$\displaystyle-J\sum_{\ell=1}^{L-1}\left(\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell+1}+\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}\right)+\frac{U}{2}\sum_{\ell=1}^{L}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell}+A\cos{(\omega
t)}\sum_{\ell=1}^{L}\ell\,\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell},$ (S4)

Changing to a frame of reference rotating with the same frequency $\omega$,
generated by the unitary transformation

$\displaystyle\hat{R}$
$\displaystyle=\exp\left\\{\text{i}\frac{A}{\hbar}\int_{0}^{t}\mathrm{d}t^{\prime}\,\cos{(\omega
t^{\prime})}\sum_{\ell=1}^{L}\ell\,\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\right\\},$
(S5)

we obtain a transformed Hamiltonian

$\displaystyle\tilde{H}(t)$
$\displaystyle=-J\left(\exp\left\\{\text{i}\frac{A}{\hbar\omega}\sin{(\omega
t)}\right\\}\sum_{\ell=1}^{L-1}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell+1}+\mathrm{h.c.}\right)+\frac{U}{2}\sum_{\ell=1}^{L}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell}.$
(S6)

In the high-frequency approximation, the dynamics of the system can be
captured at stroboscopic times — integer multiples of the modulation period
$T$ — by a time-independent effective Hamiltonian $\hat{H}_{\mathrm{eff}}$,
which can be perturbatively approximated using an infinite series in powers of
$\omega^{-1}$ based on the Magnus expansion [75, 76, 66]. At zeroth order (in
the energy frame moving with the periodic tilt [50]), the Hamiltonian sees its
hopping strength renormalized as [44, 45]:

$\displaystyle\hat{H}_{\text{eff}}^{(0)}=\frac{1}{T}\int_{0}^{T}\differential
t\,\tilde{H}(t)=\frac{U}{2}\sum_{\ell=1}^{L}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell}-J\mathcal{J}_{0}\left(\frac{A}{\hbar\omega}\right)\sum_{\ell=1}^{L-1}\left(\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell+1}+\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}\right).$
(S7)

The first order correction is given by

$\displaystyle\hat{H}_{\text{eff}}^{(1)}$
$\displaystyle=\frac{1}{2T\text{i}}\int_{0}^{T}\differential
t_{1}\int_{0}^{t_{1}}\differential t_{2}\,[\tilde{H}(t_{1}),\tilde{H}(t_{2})]$
$\displaystyle=\frac{UJ\pi}{2\hbar\omega}\mathcal{H}_{0}\left(\frac{A}{\hbar\omega}\right)\sum_{\ell=1}^{L-1}\left(\hat{a}^{\dagger}_{\ell}\hat{n}_{\ell}\hat{a}_{\ell+1}-\hat{a}^{\dagger}_{\ell}\hat{n}_{\ell+1}\hat{a}_{\ell+1}+\hat{a}^{\dagger}_{\ell+1}\hat{n}_{\ell}\hat{a}_{\ell}-\hat{a}^{\dagger}_{\ell+1}\hat{n}_{\ell+1}\hat{a}_{\ell}\right)$
$\displaystyle\quad+\frac{J^{2}\pi}{\hbar\omega}\mathcal{H}_{0}\left(\frac{A}{\hbar\omega}\right)\mathcal{J}_{0}\left(\frac{A}{\hbar\omega}\right)\left(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{L}\hat{a}_{L}\right),$
(S9)

where $\mathcal{J}_{0}(x)$ is the zeroth order Bessel function and
$\mathcal{H}_{0}(x)$ is the zeroth order Struve function. In the limit
$\hbar\omega\gg J,U$, one can negelct the effect of the first (and higher)
order term and approximate
$\hat{H}_{\text{eff}}\approx\hat{H}_{\text{eff}}^{(0)}$. Under such a
scenario, one can engineer a Floquet quasi-adiabatic [77, 78] state
preparation scheme [44]: starting from the ground state of
$\hat{H}(t)=\hat{H}^{(0)}_{\text{eff}}$ at $A=0$, one slowly ramps up
$A/\hbar\omega$ (at fixed $\omega$) up to the first zero $x_{0}\approx 2.4$ of
$\mathcal{J}_{0}(x)$. Following the ground state of
$\hat{H}^{(0)}_{\text{eff}}$, the resulting state is the maximally fragmented
state $\ket{N/L}$, ground state of the purely interacting Hamiltonian.

Figure S3:  Fragmented state preparation by fast modulation of the tilt,
starting from the ground state of (S4) with $A=0$ and $U/J=1$. (a) Tilt
$\delta(t)=A(t)\cos{(\omega t)}$ (grey line, $\omega=30J/\hbar$) with
$A(t)=mt$ (blue envelope, $m=0.2J^{2}/\hbar$) until $t_{\text{f}}\approx
2.4\,\hbar\omega/m$ and constant thereafter. $\mathcal{J}_{0}(x)$ (green line)
with $x=A/\hbar\omega$. (b) Fidelities to the twin Fock state of
$\ket{\psi(t)}$ ($\mathcal{F}^{\text{eff}}$, solid line) and of the ground
state of (S7) (dotted line), and fragmentation entropy $S_{\mathrm{F}}$ of
$\ket{\psi(t)}$ for $N=16$ (blue) and 32 (red).

Numerical results for $U/J=1$ and $L=2$ and $3$ are shown in Fig. S3. The
protocol is drawn in panel (a), showing $\delta(t)$ with frequency
$\omega=30\,J/\hbar$, and envelope $A(t)=mt$ with slope $m=0.2J^{2}/\hbar$,
reaching $A(t)=x_{0}$ at time $t_{\text{f}}\approx
2.4\,\hbar\omega/m=360\,\hbar/J$. Upon reaching its target value, the
amplitude is maintained to observe the realization. The second row features
results of numerical simulations for $N=18$ and $30$ at $L=2$ (Panel (b)) and
$L=3$ (Panel (c)). We show, as a function of time, the fragmentation entropy
$S_{\mathrm{F}}$ of the state propagated with the exact Hamiltonian (S4), as
well as two fidelities to $\ket{N/L}$: the fidelity $\mathcal{F}^{\text{ex}}$
of the state $\ket{\psi(t)}$ propagated according to the exact Hamiltonian,
and the fidelity $\mathcal{F}^{\text{eff}}$ of the ground state of
$\hat{H}_{\text{eff}}^{(0)}$ (S7) (expected to reach 1 $\forall N$). For
$N=18$, we achieve exact preparation of $\ket{N/2=9}$, maximizing the final
values of $\mathcal{F}^{\text{ex}}$ and $S_{\text{F}}$. For $N=30$ however,
while the fragmentation entropy reaches almost $\ln(2)$, the final exact
fidelity is not maximized which indicates a breakdown of the effective
Hamiltonian approximation.

Figure S4: (a) Instantaneous Floquet quasi-energy spectrum $\varepsilon_{n}$
as a function of the ramping time $t$ of the modulation amplitude (black dots)
and quasi-energy of the Floquet state maximizing the overlap with $\ket{N/2}$
(red dashed line) for $N=32$. The projection of $\ket{\psi(t)}$ on the
instantaneous Floquet states is encoded in the radius and colorbar shade. (b)
Final exact fidelity $\mathcal{F}^{\text{ex}}$ as a function of the ramp
duration $t_{\text{f}}$ for different $N$ at $\hbar\omega=30J$. (c)
$\mathcal{F}^{\text{ex}}$ as a function of $\omega$ for a fixed ramp slope
$m=0.2J^{2}/\hbar$

To understand the reason behind the failure of the effective model, we plot
the instantaneous Floquet quasi-energy spectra [66] of (S4) as a function of
$t$ (and thus $A$) for $L=2$ and $N=30$ in panel (a) of Fig. S4 . Following
how the prepared state projects in the spectrum, we see that the failure of
the effective Hamiltonian to model the exact dynamics corresponds to a leakage
towards other, untargeted Floquet states [see purple eigenenergies in Panel
(a)]. Indeed, many narrow avoided crossings are to be seen, resulting from the
Floquet folding of the static spectrum (a realistic necessity considering the
width of the many body spectrum [66, 78]). Hence, as we consider more bosons,
the Hilbert space dimension increases (linearly for $L=2$), and with it the
number of fortuitous avoided crossings that the prepared state has to cross
diabatically [66]. To achieve better fidelity in many-body “adiabatic” Floquet
engineering [77, 78], one must embrace diabaticity and moderately increase the
slope of parameter variation to favor Landau-Zener transitions through the
inevitable narrow crossings. One can see in Fig. S3(b) that there exists an
optimum fidelity achievable with this quasi-adiabatic Floquet driving, however
quite limited in $N$ at fixed $\omega$.

We finally study the requirement on $\omega$ as a function of $N$ (at fixed
slope $m$ and $L=2$). In Fig. S3(c), one can see that the driving frequency
$\omega$ required to maintain a high fidelity (i.e. proficiency of the high
frequency approximation) scales almost linearly with $N$. This scaling puts a
limit on the maximum number of atoms that one can use to create a fragmented
state in a realistic cold-atom setup using this method. Indeed, for
experimentally realistic parameters, the single-band Bose-Hubbard model is
valid for typical energy scales $E\ll 100J$ [79, 66]; since the maximum
$\omega$ is limited to around $30J$ before the driving strength
$A(t_{\text{f}})=2.4\hbar\omega$ (zero of the Bessel function) reaches $100J$,
limiting $N$ to approximately 32 atoms.

##  S3 C. Negating the hopping by a large constant tilt

We show how a strong tilt between the modes of the Bose-Hubbard model
effectively annihilates the hopping, detailing the derivation for $L=2$ and
$L=3$ modes.

###  S3.1 a. $L=2$

We consider the Hamiltonian

$\hat{H}=U\hat{J}^{2}_{z}-2J\hat{J}_{x}+2\delta\hat{J}_{z},$ (S10)

in the rotating frame accessed by the unitary transformation
$|\widetilde{\psi}(t)\rangle=\hat{R}|\psi(t)\rangle$ with
$\hat{R}=\text{exp}\\{\text{i}2\delta\hat{J}_{z}t/\hbar\\}$. In this frame,
the Hamiltonian reads

$\displaystyle\hat{\widetilde{H}}(t)$
$\displaystyle=\hat{R}\hat{H}\hat{R}^{\dagger}-\text{i}\hbar\hat{R}\dfrac{\differential\hat{R}^{\dagger}}{\differential
t}$
$\displaystyle=U\hat{J}^{2}_{z}-2J\,\text{e}^{\text{i}2\delta\hat{J}_{z}t/\hbar}\hat{J}_{x}\,\text{e}^{-\text{i}2\delta\hat{J}_{z}t/\hbar}$
$\displaystyle=U\hat{J}^{2}_{z}-2J\cos\left(\dfrac{2\delta
t}{\hbar}\right)\hat{J}_{x}+2J\sin\left(\dfrac{2\delta
t}{\hbar}\right)\hat{J}_{y},$ (S11)

using the algebra
$[\hat{J}_{i},\hat{J}_{j}]=\text{i}\epsilon_{ijk}\hat{J}_{k}$. If $\delta\gg
U,J$, the Hamiltonian $\hat{\widetilde{H}}$ can be approximated by its Magnus
expansion at zeroth order, namely its average over one period
$\textstyle\hat{\widetilde{H}}\approx\hat{\widetilde{H}}_{\text{eff}}^{(0)}=U\hat{J}^{2}_{z}$.
Finally, going back to the initial frame gives

$\displaystyle\hat{H}$
$\displaystyle=\hat{R}^{\dagger}\hat{\widetilde{H}}\hat{R}-\text{i}\hbar\hat{R}^{\dagger}\dfrac{\differential\hat{R}}{\differential
t}$ $\displaystyle\approx U\hat{J}^{2}_{z}+2\delta\hat{J}_{z},$ (S12)

whose eigenstates are the Fock states, among which the maximally fragmented
state $\ket{N/2}$; see Fig. 1 of the main text.

###  S3.2 b. $L=3$

We consider the titled Bose-Hubbard Hamiltonian,

$\displaystyle\hat{H}(t)$
$\displaystyle=-J\sum_{\ell=1}^{2}\left(\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell+1}+\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}\right)+\frac{U}{2}\sum_{\ell=1}^{3}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell}+\delta\left(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{3}\hat{a}_{3}\right).$
(S13)

Again, moving to a frame rotating with angular frequency $\delta/\hbar$ using
the unitary operator $\hat{R}=\mathrm{exp}\\{\text{i}\delta
t(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{3}\hat{a}_{3})/\hbar\\}$,
we get the transformed Hamiltonian

$\displaystyle\hat{\widetilde{H}}(t)$
$\displaystyle=J\sum_{\ell=1}^{2}\left(\text{e}^{\text{i}\delta
t/\hbar}\,\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell+1}-\text{e}^{-\text{i}\delta
t/\hbar}\,\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}\right)+\frac{U}{2}\sum_{\ell=1}^{3}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell}.$
(S14)

Using a similar procedure as for $L\\!=2$ for $\delta\gg U,J$, one can write
an effective Hamiltonian which till the second order of Magnus expansion is
given as

$\displaystyle\hat{\widetilde{H}}_{\text{eff}}$
$\displaystyle\approx\left(\frac{U}{2}-\frac{J^{2}U}{\delta^{2}}\right)\sum_{\ell=1}^{3}\hat{a}^{\dagger}_{\ell}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}_{\ell}-\frac{J^{2}U}{\delta^{2}}\hat{a}^{\dagger}_{2}\hat{a}^{\dagger}_{2}\hat{a}_{2}\hat{a}_{2}+\frac{4J^{2}U}{\delta^{2}}\sum_{\ell=1}^{2}\hat{a}^{\dagger}_{\ell}\hat{a}_{\ell}\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell+1}-\frac{J^{2}U}{\delta^{2}}\left(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{3}\hat{a}_{2}\hat{a}_{2}+\mathrm{h.c.}\right)$
$\displaystyle\quad+\frac{J^{2}}{\delta}\left(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{3}\hat{a}_{3}\right).$
(S15)

As a matter of fact, the terms obtained here resemble the Hamiltonian
describing rotationally invariant spin-1 gases [80, 8],

$\displaystyle\frac{\hat{H}_{\mathrm{S}=1}}{U_{s}}$
$\displaystyle=2\hat{N}+\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{1}\hat{a}_{1}\hat{a}_{1}+\hat{a}^{\dagger}_{3}\hat{a}^{\dagger}_{3}\hat{a}_{3}\hat{a}_{3}+2\left(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{2}\hat{a}_{1}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}^{\dagger}_{3}\hat{a}_{2}\hat{a}_{3}-\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{3}\hat{a}_{1}\hat{a}_{3}\right)+2\left(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{3}\hat{a}_{2}\hat{a}_{2}+\mathrm{h.c.}\right),$
(S16)

with $U_{s}$ the spin interaction energy. However, a one-to-one mapping
between the two Hamiltonians does not exist, as $U,J$ and $\delta$ in Eq.
(S3.2) cannot be tuned to match the corresponding terms in Eq. (S16); in this
sense, the 3-mode tilted Bose-Hubbard model is not a direct analogue of spin-1
Bose gases.

##  S4 D. Quantum optimal control of $L$-tuple Fock states

Given a dynamical system described by the state vector $\vec{x}(t)$ ($0\leq
t\leq t_{\text{f}}$), and whose evolution is parameterized by control
parameters $\vec{u}(t)$, optimal control [81, 82, 83, 84] is the mathematical
formalism that allows to determine $\vec{u}(t)$ in order to optimally control
$\vec{x}(t)$. In the framework of optimal control, optimality is captured by a
cost function $\mathcal{C}$ to be minimized. By virtue of the Pontryagin
Maximum Principle, the optimal control $\vec{u}(t)$ is the one that maximizes
the Pontryagin Hamiltonian $H_{\text{P}}(t),\;\forall t\in[0,t_{\text{f}}]$
(see below, and Refs. [81, 82, 84, 48]). Quantum optimal control [46, 62, 47,
48] is the application of optimal control theory to quantum systems, with
typically $\vec{x}(t)\rightarrow\ket{\psi(t)}$.

###  S4.1 a. Implementation

In this work, we use the energy bias $\delta(t)$ between the $L$-modes of a
Bose-Hubbard system as our control parameter to optimally steer the ground
state of the system (taken at $\delta=0$ with $U/J$ fixed) towards the
$L$-tuple Fock state $\ket{N/L}$ in a fixed time $t_{\text{f}}$. We do not
place any constraint on $\delta(t)$. The cost function that we minimize is
$\mathcal{C}=1-\mathcal{F}$, with
$\mathcal{F}=|\innerproduct{N/L}{\psi(t_{\text{f}})}|^{2}$, leading to the
following definition of the Pontryagin Hamiltonian [84, 63, 85]:

$H_{\text{P}}(t)=\Re{\bra{\chi(t)}\dfrac{\partial}{\partial
t}\ket{\psi(t)}}=\dfrac{1}{\hbar}\Im{\bra{\chi(t)}\hat{H}(\delta(t))\ket{\psi(t)}},$
(S17)

using Schrödinger equation with the Hamiltonian $\hat{H}$ given in Eqs. (S10)
and (S13) for $L\\!=2$ and 3 respectively. In Eq. (S17), $\ket{\psi(t)}$ is
state during its preparation and $\ket{\chi(t)}$ is a non-normalized adjoint
state, also solution of Schrödinger equation and defined by the final
condition:

$\innerproduct{u_{n}}{\chi(t_{\text{f}})}=\dfrac{\partial\mathcal{F}}{\partial\left(\innerproduct{\psi(t_{\text{f}})}{u_{n}}\right)},$
(S18)

in any basis $\left\\{\ket{u_{n}}\right\\}$, giving
$\ket{\chi(t_{\text{f}})}=\innerproduct{N/L}{\psi(t_{\text{f}})}\ket{N/L}$.
Given an non-optimal control $\delta(t)$, the first-order correction
$\Delta\delta(t)$ required to increase $H_{\text{P}}$ (and thus $\mathcal{F}$)
is proportional to the derivative of $H_{\text{P}}$ with respect to the
control parameter [62, 63, 48], that is

$\Delta\delta(t)=\dfrac{\epsilon}{\hbar}\Im{\bra{\chi(t)}\dfrac{\partial\hat{H}(\delta(t))}{\partial(\delta(t))}\ket{\psi(t)}},$
(S19)

with $\epsilon$ a small adjustable line-search parameter [62, 48]. We
iteratively correct $\delta(t)\rightarrow\delta(t)+\Delta\delta(t)$ until we
reach the desired fidelity, following the algorithm given in Refs. [62, 63,
48].

As discussed in the main text, working around a large average tilt allows to
restrict the quantum dynamics of the evolved state to a narrow subset of
Hilbert space. We use it to our advantage by first computing an approach
control field in a Hilbert space strongly truncated “around” $\ket{N/L}$
444For instance with $L=2$ and $N=128$, the first optimization was performed
in a Hilbert space of dimension 19, only considering the basis Fock states
$\ket{n_{1}}$ with $55\leq n_{1}\leq 73$.. We then use this first control
field as a strong guess for a second optimization in the complete Hilbert
space.

###  S4.2 b. Details of the controls featured in the main text

We detail in Fig. S5 the twin Fock state preparations presented in Table 1 of
the main text, i.e. as a function of $N$. For each number of bosons
considered, we show the optimal-control tilt $\delta$, the three metrics of
state preparation considered in the main text (namely the fidelity
$\mathcal{F}$ to the twin Fock state, the fragmentation entropy $S_{\text{F}}$
and the number squeezing parameter $\xi^{2}_{N}$), and the projection over the
Fock basis as a function of the preparation time. In particular, panel (d),
with $N=32$, corresponds to Fig. 2 of the main text.

Figure S5:  Details of twin Fock state preparation by optimal control. (a) Top
row: Optimal control tilt as a function of time (solid black line) to drive
the ground state of the unbiased system with $U/J=1$ and $N=4$ towards the
twin Fock state $\ket{N/2}$, and initial guess $\delta(t)=\delta_{0}$ (dashed
grey line). Middle row: Fidelity $\mathcal{F}$ to $\ket{N/2}$ (blue solid
line), fragmentation entropy $S_{\text{F}}$ (divided by $\ln(2)$, orange
dashed) and number squeezing parameter $\xi^{2}_{N}$ (dB scale, black dash-
dotted) as a function of time. Bottom row: Fock space projection during the
preparation. (b-h) Same as (a) for $N=8$, 16, 32, 64, 128, 256, 512
respectively. For (e-h), the displayed evolution in Fock space is cropped over
$n_{1}=N/2\pm 10$ (bottom rows).

###  S4.3 c. Robustness to number fluctuations

To demonstrate the significant robustness of our optimal-control scheme to
number fluctuations, which results from inducing fragmentation around a large
tilt (see main text), we now apply the control fields that were optimized for
a given $N_{0}$ to other settings with $N\neq N_{0}$. Figure S6 shows those
results for $N_{0}=16$, 32, 64 and 128, corresponding to the controls of Fig.
S5(c-f). The fidelity remains above 90% of its maximum over the large interval
$N\in[10,18]$ for $N_{0}=16$, $N\in[20,46]$ for $N_{0}=32$, $N\in[44,80]$ for
$N_{0}=64$ and $N\in[96,160]$ for $N_{0}=128$. Figure S6 further illustrates
that the fidelity is a fragile metric compared to the fragmentation entropy,
which remains very high around $N_{0}$, furthermore over a larger interval as
$N$ increases.

We mention that dedicated approaches exist to enhance robustness in optimal
control theory (see [48] and references therein), however these go beyond the
scope of this study.

Figure S6:  Robustness of the optimal control preparations to number
fluctuations. Fidelity (top row), fragmentation entropy (divided by $\ln(2)$,
middle row) and number squeezing (bottom row) as a function of $N$. Columns
(a-d) correspond to control fields optimized for $N_{0}=16$ , 32, 64, 128
respectively (marked by the vertical black lines). Solid lines shows the three
metrics for the optimized states and dashed lines for the initial state
(ground state at $U/J=1$ and $\delta=0$). The blue areas mark the intervals
around $N_{0}$ over which the fidelity stays above 90% of its maximum. The
corresponding control fields are plotted in Fig. S5(c-f).

###  S4.4 d. Quantum optimal preparations of triple Fock states

Figure S7 shows results of triple Fock state preparation by optimal control
for $N=3$, 6, 9, 15, 21 and 27 bosons. Our approach is the same as for
$L\\!=\\!2$ (see main text and previous sections). Details of the preparations
are given in Table S1.

Figure S7:  Triple Fock state preparations by optimal control. (a) Top row: Optimal control tilt as a function of time (solid black line) to drive the ground state of the unbiased system with $U/J=1$ and $N=3$ towards the triple Fock state $\ket{N/3}$, and initial guess $\delta(t)=\delta_{0}$ (dashed grey line). Bottom row: Fidelity $\mathcal{F}$ to $\ket{N/3}$ (blue solid line), fragmentation entropy $S_{\text{F}}$ (divided by $\ln(3)$, orange dashed) and number squeezing parameter $\xi^{2}_{N}$ (dB scale, black dash-dotted) as a function of time. (b-f) Same as (a) for $N=6$, 9, 15, 21, 27 respectively. See Table S1 for details.  $N$ |  3 |  6 |  9 |  15 |  21 |  27
---|---|---|---|---|---|---
$\delta_{0}(J)$ |  2.63 |  7.50 |  13.5 |  22.5 |  26.3 |  40.5
$t_{\text{f}}\,(\hbar/J)$ |  15.3 |  11.8 |  5.93 |  4.12 |  4.25 |  2.99
$\mathcal{F}$ |  0.999 |  0.990 |  0.980 |  0.971 |  0.962 |  0.937
$S_{\text{F}}/\ln(3)$ | $>0.999$ |  0.996 |  0.997 |  0.999 |  0.987
$\begin{array}[]{c}\xi^{2}_{N}\\\ \text{(dB)}\end{array}\left(\begin{array}[]{l}\psi_{\text{f}}\\\ \phi_{0}\end{array}\right)$ |  -26.9 |  14.4 |  -14.8 |  -15.7 |  -14.9 |  -14.9
|  -0.994 |  -2.21 |  -2.97 |  -3.99 |  -4.68 |  -5.20
Table S1:  Triple Fock state preparations by optimal control ($U/J=1$) for
increasing number of bosons $N$, with guess tilt $\delta_{0}$, control
duration $t_{\text{f}}$, final fidelity $\mathcal{F}$ to $\ket{N/3}$,
fragmentation entropy $S_{\text{F}}$ (divided by $\ln(3)$) and number
squeezing parameter of prepared and initial states ($\ket{\psi_{\text{f}}}$
and $\ket{\phi_{0}}$ resp.).

##  S5 E. Number fluctuations of reference

We derive, for $L\\!=\\!2$ and 3, the uncorrelated number fluctuations
$\Delta\hat{J}^{2}_{z,\text{ref}}$ used as the reference value when computing
the number squeezing parameter
$\xi^{2}_{N}(\psi)=\Delta\hat{J}^{2}_{z,\psi}/\Delta\hat{J}^{2}_{z,\text{ref}}$.
In each case considered, the reference state that is taken corresponds to the
ground state of the hopping Hamiltonian with $N$ bosons.

###  S5.1 a. Uncorrelated number fluctuations for $L\\!=\\!2$

The reference state is the coherent spin state given in Eq. (S3) for $n=0$,
explicitly

$\ket{\phi^{(0)}_{0}}=\dfrac{\left(\hat{a}^{\dagger}_{1}+\hat{a}^{\dagger}_{2}\right)^{N}}{2^{N/2}\sqrt{N!}}\ket{0,0}.$
(S20)

Number fluctuations for that state are given by the variance on $\hat{J}_{z}$:
$\Delta\hat{J}^{2}_{z,\text{ref}}=\langle\hat{J}_{z}^{2}\rangle_{\text{ref}}-\langle\hat{J}_{z}\rangle^{2}_{\text{ref}}$,
with
$\hat{J}_{z}=(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{2}\hat{a}_{2})/2$.
As the reference state lies symmetrically on the equator of the Bloch sphere,
the first moment of $\hat{J}_{z}$ is zero:

$\displaystyle\langle\hat{J}_{z}\rangle_{\text{ref}}$
$\displaystyle=\dfrac{1}{2^{N+1}N!}\bra{0,0}\left(\hat{a}_{1}+\hat{a}_{2}\right)^{N}\left(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{2}\hat{a}_{2}\right)\left(\hat{a}^{\dagger}_{1}+\hat{a}^{\dagger}_{2}\right)^{N}\ket{0,0}$
$\displaystyle=\dfrac{1}{2^{N+1}N!}\left(\sum_{n=0}^{N}\binom{N}{n}\sqrt{(N-n)!n!}\bra{N-n,n}\right)\left(\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{2}\hat{a}_{2}\right)\left(\sum_{m=0}^{N}\binom{N}{m}\sqrt{(N-m)!m!}\ket{N-m,m}\right)$
$\displaystyle=\dfrac{1}{2^{N+1}}\sum_{n=0}^{N}\binom{N}{n}(N-2n)=0,$ (S21)

where the last equality is evaluated using the following identities on the sum
of the binomial coefficients:

$\sum_{n=0}^{N}\binom{N}{n}=2^{N}\qquad\text{and}\qquad\sum_{n=0}^{N}n\binom{N}{n}=2^{N-1}N.$
(S22)

Number fluctuations are thus given by the second moment of $\hat{J}_{z}$, that
is

$\Delta\hat{J}^{2}_{z,\text{ref}}=\langle\hat{J}_{z}^{2}\rangle_{\text{ref}}=\dfrac{1}{2^{N+2}}\sum_{n=0}^{N}\binom{N}{n}(N-2n)^{2}=\dfrac{N}{4},$
(S23)

using Eqs. S22 and

$\sum_{n=0}^{N}n^{2}\binom{N}{n}=2^{N-2}N(N+1).$ (S24)

The result in Eq. (S23) is well known, see e.g. Refs. [68, 5]

###  S5.2 b. Uncorrelated number fluctuations for $L\\!=\\!3$

For $L=3$, we evaluate the number fluctuations partially, via the variance of
the population imbalance between modes 1 and 3:
$\hat{J}_{z}=\hat{a}^{\dagger}_{1}\hat{a}_{1}-\hat{a}^{\dagger}_{3}\hat{a}_{3}$
[6, 20]. The hopping Hamiltonian (without coupling between modes 1 and 3) is

$\displaystyle\hat{H}$
$\displaystyle=-J\left(\hat{a}^{\dagger}_{1}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}_{3}+\hat{a}^{\dagger}_{2}\hat{a}_{1}+\hat{a}^{\dagger}_{3}\hat{a}_{2}\right)$
$\displaystyle=-\sqrt{2}J\left(\hat{a}^{\dagger}_{+}\hat{a}_{+}-\hat{a}^{\dagger}_{-}\hat{a}_{-}\right),$
(S25)

with $\hat{a}_{\pm}=(\hat{a}_{1}\pm\sqrt{2}\hat{a}_{2}+\hat{a}_{3})/2$. The
reference state for number fluctuations is the ground state of Hamiltonian
(S5.2), namely

$\displaystyle\ket{\phi^{(0)}}$
$\displaystyle=\dfrac{\left(\hat{a}^{\dagger}_{1}\pm\sqrt{2}\hat{a}^{\dagger}_{2}+\hat{a}^{\dagger}_{3}\right)^{N}}{2^{N}\sqrt{N!}}\ket{0,0,0}$
$\displaystyle=\dfrac{1}{2^{N}\sqrt{N!}}\sum_{n=0}^{N}\binom{N}{n}2^{(N-n)/2}\sqrt{(N-n)!}\sum_{r=0}^{n}\binom{n}{r}\sqrt{(n-r)!r!}\ket{n-r,N-n,r},$
(S26)

which is not a coherent spin state in the absence of coupling between modes 1
and 3. For that state (having equal populations in modes 1 and 3), the first
moment of $\hat{J}_{z}$ is zero:

$\langle\hat{J}_{z}\rangle_{\text{ref}}=\dfrac{1}{2^{N}}\sum_{n=0}^{N}\binom{N}{n}\dfrac{1}{2^{n}}\sum_{r=0}^{n}\binom{n}{r}(n-2r)=0,$
(S27)

as the second sum is null from Eqs. (S22). The variance of $\hat{J}_{z}$ is
thus

$\Delta\hat{J}^{2}_{z,\text{ref}}=\langle\hat{J}_{z}^{2}\rangle_{\text{ref}}=\dfrac{1}{2^{N}}\sum_{n=0}^{N}\binom{N}{n}\dfrac{1}{2^{n}}\sum_{r=0}^{n}\binom{n}{r}(n-2r)^{2}=\dfrac{N}{2},$
(S28)

using Eqs. (S22) and (S24).

For completeness, we give that the uncorrelated number fluctuations for
$L\\!=\\!3$ with coupling between modes 1 and 3 (i.e. considering a coherent
spin state, e.g. the ground state of a ring of three spatial sites) are given
by $\Delta\hat{J}^{2}_{z,\text{ref}}=2N/3$.