Fast initialization of a high-fidelity quantum register using optical superlattices

We propose a method for the fast generation of a quantum register of addressable qubits consisting of ultracold atoms stored in an optical lattice. Starting with a half filled lattice we remove every second lattice barrier by adiabatically switching on a superlattice potential which leads to a long wavelength lattice in the Mott insulator state with unit filling. The larger periodicity of the resulting lattice could make individual addressing of the atoms via an external laser feasible. We develop a Bose-Hubbard-like model for describing the dynamics of cold atoms in a lattice when doubling the lattice periodicity via the addition of a superlattice potential. The dynamics of the transition from a half filled to a commensurately filled lattice is analysed numerically with the help of the time evolving block decimation algorithm and analytically using the Kibble–Zurek theory. We show that the timescale for the whole process, i.e. creating the half filled lattice and subsequent doubling of the lattice periodicity, is significantly faster than adiabatic direct quantum-freezing of a superfluid into a Mott insulator for large lattice periods. Our method therefore provides a high-fidelity quantum register of addressable qubits on a fast timescale.


Introduction
Systems of cold atoms trapped in optical lattices provide the unique opportunity to coherently manipulate a large number of atoms [1]- [3]. The remarkable degree of experimental control offered by these systems, as well as the possibility to use the internal hyperfine states of the atoms to encode qubits, make them particularly suited for quantum information processing (QIP). In this context, optical lattices in a Mott-insulating (MI) state with unit filling can be viewed as the realization of a quantum register, and it is possible to collectively manipulate the qubits stored in such a register experimentally [4,5]. However, in many quantum computing schemes based on neutral atoms stored in optical lattices the application of single qubit gates [3] or single qubit measurements [6] requires the ability to address single atoms with a focused laser beam. This remains experimentally challenging since these operations have to be performed without perturbing the state of other atoms in their vicinity.
A number of strategies have been proposed to circumvent this problem by using global operations [7]- [9], e.g. via 'marker atoms'which are moved to a particular lattice site and interact with the corresponding register qubit such that an external laser affects only that qubit [10]. Another way is simply to use a quantum register in which the atoms are distant enough such that they can be addressed individually by a laser. This method requires an optical lattice with filling factor n = 1 in the MI state with a sufficiently large distance between the atoms 2 . The initialization time of a MI state is proportional to the tunnelling time of the atoms between neighbouring sites. Therefore, by using the conventional method of quantum-freezing a superfluid (SF) state [4], [12]- [14], it scales exponentially with the lattice spacing [15]- [17].
In this paper, we propose an alternative method to generate a long wavelength lattice with one atom per site in the MI state. Starting with a one-dimensional (1D) lattice with a short period-and The periodicity of the lattice is then progressively doubled (b) until it reaches its final profile (c) associated with the parameter s = 0. The number of sites M corresponds to the number of unit cells in the large lattice limit. By starting with a filling factor of n = 1/2, this procedure leads to a lattice with filling factor of n = 1.
hence a short initialization time-and filling factor n = 1/2 we remove every second potential barrier by adiabatically turning on a superlattice. This superlattice potential has already been experimentally realized [18,19]. This procedure eventually leads to a long wavelength lattice in which the periodicity has been doubled and where the atoms are in a MI state with n = 1 (see figure 1). This scheme does not require changing the angles of the intersecting laser beams. Furthermore, we show that our method allows for the initialization of a MI state with unit filling factor on a timescale which, although scaling exponentially with the final lattice spacing, is approximately one order of magnitude smaller than the direct quantum-freezing method. Although we only consider the case of a homogeneous lattice, the results presented in this paper extend to the case of weak harmonic confinements, quartic [20] and box traps [21]. This paper is structured as follows. In section 2, we introduce the model used to describe the system dynamics. In section 3.1, we discuss ground state properties, particularly two-site correlation functions and quasi-momentum distributions. Also, we present and discuss numerical results for the probability of staying in the ground state during the transition depending on the speed of the ramping. In section 3.2, we apply the analytical Kibble-Zurek (KZ) theory and compare it with our numerical results. Finally, we summarize and conclude in section 4.

Model
We consider a gas of interacting ultracold bosonic atoms loaded into a 3D optical lattice. The lattice is formed by pairwise orthogonal standing wave laser fields and its optical potential is given by [22] Here V T is the depth of the potential in the y-and z-directions created by pairs of lasers with wavenumber k = 2π/λ, wave length λ and period a T = λ/2. The extension to the case of different optical potentials in the y-and z-directions is straightforward. In the x-direction two pairs of laser beams with a long wavelength λ L and short wavelength λ S = λ L /2 are applied. The potential in the x-direction is thus given by

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
with V 0 the depth of the lattice, k L = 2π/λ L and k S = 2π/λ S . The depths of the potentials will be expressed in units of the recoil energy E R = k 2 L /2m with m the mass of the atoms (takinḡ h = 1 throughout). The parameter s ∈ [0, 1] is determined by the relative intensities of the two pairs of lasers. By changing s from 1 to 0 the lattice profile is continuously transformed from a sinusoidal potential with a small period a S = λ S /2 to one with a long period a = λ L /2, thus halving the number of lattice sites per unit length along the x-direction (see figure 1). The lattice constant a corresponds to the size of a unit cell for s < 1. 3 We refer to the lattice profile with parameters s = 1 and s = 0 as to the small lattice limit and the large lattice limit, respectively.
The Hamiltonian of the system in second quantization readŝ whereĥ 0 (r) = −(1/2m)∇ 2 + V OL (r) is the one-particle Hamiltonian. The symbolr = (r, s) represents the position variable r and lattice parameter s. The interaction between the atoms is modelled by s-wave scattering with g = 4πa s /m where a s is the s-wave scattering length. The bosonic field operators obey the usual commutation relations [ψ(r), ψ † (r )] = δ(r − r ) with δ denoting the Dirac delta function. We restrict our considerations to the case where V T is sufficiently large so that motion along the y-and z-directions is frozen. The dynamics of the system is then effectively 1D along the x-direction. As shown in figure 2(a) the two lowest bands of the Hamiltonianĥ 0,x = −(1/2m) (d/dx) 2 + V S (x, s) are separated in energy by much less than the typical motional excitation energy E ex = √ 4V 0 E R [22] for values s ≈ 1. 4 Therefore, despite assuming that atoms loaded into the lattice are ultracold, we have to consider the two lowest Bloch bands in x-direction to obtain an accurate description of the atomic dynamics. However, excitations to higher bands in the y-and z-directions are neglected in our investigations since we assume that the temperature of the atomic cloud is k B T E ex . We then expand the bosonic field operator aŝ where M is the number of lattice sites andα † i (α = a, b) creates a particle in the mode associated with the localized function φ α,i (r) centred at site i. The mode functions φ α,i (r) = w α,i (x, s)W i,0 (y)W i,0 (z) are factorized into a product of well-localized Wannier functions (WFs) W i,0 of the lowest Bloch band in the y-and z-directions [23] and mode functions w α,i (x, s) in x-direction.
The aim of the next section is to describe the single particle dynamics in the tight binding (TB) approximation. If we were to use WFs for w α,i (x, s) this approximation would restrict our model to sinusoidal Bloch bands [24]. Because of the deviation of the lowest two bands from a sinusoidal dispersion relation (see figure 2(a)) when s ≈ 1 we instead use generalized WF (GWFs) for w α,i (x, s) (see appendix A for a detailed definition and a description of their properties) [25]. By exploiting the optimization procedure described in appendix A, we calculate GWFs w α,i (x, s) which are well localized at lattice sites i. Typical shapes of GWFs and the effects of optimizing their localization are shown in figure 3. We note that these GWFs are in

Single-particle Hamiltonian
Inserting the approximate field operator equation (4) into the first term of equation (3) yields the single-particle part of the Hamiltonian in terms ofâ † i andb † i . Applying the TB approximation, which amounts to keeping only nearest-neighbour hopping terms, the single-particle Hamiltonian can be approximated bŷ where is the hopping matrix element between neighbouring sites along the x-axis and is the local on-site energy of a particle in mode α. Note that hopping between modes a and b within one site is not allowed by the symmetry properties of the GWFs. However, the inclusion of nonzero hopping matrix elements J ab and J ba is essential to accurately reproduce the single particle behaviour of the full Hamiltonian (3). The symmetry properties of WFs would not allow the inclusion of these terms [24].  For periodic boundary conditions, the parameters V α and J α,β can be found from the band structure without explicit calculation of the mode functions (see appendix B). The numerical values of V α and J α,β for V 0 = 30E R are shown in figure 4(a). Using these parameters, we find thatĤ 0 (s) very accurately reproduces the band structure of the exact Hamiltonian for all values of s, thus justifying the utilization of the TB approximation and the corresponding GWFs.

Interaction Hamiltonian
To calculate the interaction matrix elements ofĤ the explicit form of the localized GWFs is needed. We find that for V 0 > 10E R , off-site interaction terms are at least two orders of magnitude smaller than on-site interactions. We therefore only keep the dominant on-site terms and find the interaction HamiltonianĤ Î wheren a i =â † iâ i andn b i =b † ib i are the site-occupation number operators. The on-site interaction matrix elements are given by Parameters of the effective HamiltonianĤ eff as functions of s for a lattice depth of V 0 = 30E R and V T = 60E R withã = a a 2 T . The hopping matrix elements (a) have been calculated using the method described in appendix B, while the on-site interaction energies (b) are calculated using optimized GWFs.
The numerical values of U α,β as a function of the lattice profile s are shown in figure 4(b). Figure 4(b) shows that their values become equal as s → 1 for sufficiently large V 0 .
Combining the single-and two-particle contributions the effective Hamiltonian describing the system dynamics is given bŷ By usingĤ eff (s) for s varying in time, we implicitly assume that the system adiabatically follows changes in the mode functions w α,i . For all dynamical calculations carried out in this work we have carefully chosen the time dependence of s so that such non-adiabatic contributions can safely be neglected.

Limiting cases
In the small lattice limit (s = 1), the superpositions ϕ L,i = (w a,i − w b,i )/ √ 2 and ϕ R,i = (w a,i + w b,i )/ √ 2 correspond to mode functions localized in the left and in the right well of site i respectively (see figures 3(c) and (d)). The associated bosonic operators are defined bŷ Given that the new mode functions are sufficiently localized within each well, the parameters of H eff for s = 1 can be written as Notice that the parameters shown in figure 4 are consistent with these equations. Expressing the Hamiltonian (10) in terms of the operatorsĉ † α,i (α = L, R) and using the parameters (12) we find that Therefore, as expected,Ĥ eff (s = 1) is equivalent to the one-band Bose-Hubbard model (BHM) [22] with 2M sites.
The mode functions keep their symmetry with the two peaks moving towards the centre of the cell when s is decreased. When s ≈ 0 is reached w a,i (x, s) and w b,i (x, s) are equivalent to the WFs of the first and the second Bloch band, respectively. In this limit we obtain a standard two band BHM for M sites.

Timescale for the preparation of the quantum register
In this section, we present numerical as well as analytical results characterizing the ground state properties of the system and the timescale necessary to initialize the quantum register.

Numerical results
The numerical calculations have been carried out using both the exact matrix representation of the HamiltonianĤ eff and the time-evolving Block decimation (TEBD) algorithm (see appendix C for details).
We have evaluated the ground state and dynamical properties of our system for two different values of g corresponding to different interaction regimes. These values have been chosen such that in the large lattice limit (with filling factor n = 1) the ground state is a Mott-insulator state for g 1 and g 2 with J aa /U aa < 0.3 [26]. In the small lattice limit, g 1 and g 2 produce ground states corresponding to a strongly interacting Tonks-Girardeau (TG) gas (U/J = 215) and a SF (U/J = 7), respectively.

Ground states properties.
In the large lattice limit, the ground state |ψ 0,L of the system is populated exclusively by particles in the lowest Bloch band, i.e. a-mode particles. Therefore, in this limit, we only consider the one-particle density matrix given by where |ψ is the state of the system. The quasi-momentum distribution of particles in the a-mode is given by [27] n As expected, in this limit and for commensurate filling n = 1 the ground state is a Mott-insulator for both values of g (see figures 5(a)-(d)). The quasi-momentum distributions show that in this limit particles are uniformly distributed in the first band (see figures 5(a) and (c)).  Figure 5. Ground state one-particle density matrix and quasi-momentum distribution for M = 12 lattice sites and the parameters of panels parameters. The panels (a)-(d) correspond to the large lattice limit and the panels (e)-(h) to the small lattice limit. The quasi-momentum distribution (a), (e) and one-particle density matrix (b), (f) are for g = g 1 . The quasi-momentum distribution (c), (g) and one-particle density matrix (d), (h) are for g = g 2 .

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
In the small lattice limit, the one-particle density matrix is given by where the operatorsĉ † i are constructed from theâ † i andb † i operators using the transformations (11). The quasi-momentum distribution is given by In this limit (with filling factor n = 1/2), the characteristics of the ground state |ψ 0,S depend on the value of g. For g = g 1 , the ratio U/J = 215 and the ground state's correlations as well as the quasi-momentum distributions (see figures 5(e) and (f)) are characteristic of a TG gas (see e.g. [28] and references therein). The value g = g 2 yields the ratio U/J = 7 and the system exhibits the behaviour of a SF (see figures 5(g) and (h)), with particles occupying mainly the q = 0 momentum state. Notice that 1D systems described by the BHM with a fixed filling factor n = 1 cross the MI-SF phase boundary at the critical point (J/U) c ≈ 0.3 [26]. Thus, in our case the SF behaviour of the system in the small lattice limit is due to the fractional filling of the lattice.

3.1.2.
Simulation of the dynamics. Starting from a system with half-filling and a lattice profile s = 1, we investigate the timescale required to obtain a nearly perfect MI state-or quantum register-with filling factor n = 1 by ramping the lattice profile down to s = 0. The quality of the register is determined by the fidelity defined as the overlap between the state of the system |ψ at the end of the ramping process and the ground state in the large lattice limit |ψ 0,L . Furthermore, we calculate the fluctuations of the number of particles in the a-mode at site i where • = ψ| • |ψ . Since particle-number fluctuations are suppressed in a MI state, nonzero fluctuations indicate the presence of excitations, such as double occupancies or particles in the second band, in the final state. We test different ramps by simulating the system dynamics between t i = 0 to t f = τ Q where τ Q is the ramping time (the time required to complete the ramping process from s = 1 to s = 0). Here, each ramp σ corresponds to a function s = s σ (t) varying from s σ (0) = 1 to s σ (τ Q ) = 0.
For linear ramping we use s lin (t) = (τ Q − t/E R )/τ Q . The fidelity and particle-number fluctuations obtained using this strategy for different quench times and g = g 1 and g = g 2 are shown in figures 6(a)-(d). The linear ramp is shown in figure 7(a).
Another ramping strategy we use consists of adapting the velocity of the ramp proportionally to the energy gap between the ground and the first excited state. This ramp is denoted by s gap (t).  (e) and (f) the ramp s gap and g = g 1 ; (g) and (h) the ramp s man and g = g 2 . The vertical lines indicate the value of the particle tunneling time (t tun ∼ 1/J aa ) in the large lattice limit.
sites and V 0 = 30E R . The ramp s gap (t) for g = g 1 is shown in figure 7(a). We expect this ramping strategy to be more efficient than the linear one, since accelerating the ramp when the gap is large while slowing it down when the gap is small should suppress transitions of particles to excited levels. Our numerical calculations have shown that when g = g 1 , the utilization of this ramp does indeed reduce the quench time needed to obtain a nearly perfect fidelity to a half of the tunnelling time in the large lattice limit (see figures 6(e) and (f)). For g = g 2 , we find that compared to s lin (t), this strategy only marginally improves the fidelity. In order to reduce the time required to obtain a given fidelity for g = g 2 , a more sophisticated ramping strategy is needed. In the following, we provide a simple method to estimate the efficiency of different ramps without running a complete dynamical simulation of the system. The transition probability between the ground and some excited state |k at time t for a ramp σ is approximately given by [29] P σ where |k = |k(σ, t) is the kth instantaneous eigenstate 5 ofĤ eff (s σ (t)) and ω 0k is the transition frequency between the ground state and |k . Therefore, the assessment of the efficiency of a ramp can be done by evaluating the functional where the index k runs over all the values associated with levels connected to the ground state. The functional A(σ, τ Q ) corresponds to the average transition probability per unit time for a given strategy σ and quench time τ Q . We calculate the value of the functional A(σ, τ Q ) numerically via exact diagonalization ofĤ eff for a small system. This method allows to optimize ramps by minimizing the value of A(σ, τ Q ). The optimized ramp for a small system is then used in the simulation of larger systems. For instance, the strategy s man (t) shown in figure 7(a) was designed and optimized manually using this method. For g = g 2 , we find that A(lin, τ Q )/A(man, τ Q ) ≈ 2.3 for τ Q = 20/E R (see figure 7(b)), thus showing the better efficiency of the strategy s man (t) compared to s lin (t). As shown in figure 6, dynamical simulations of the system with g = g 2 confirm that this strategy reduces the time required to obtain a given fidelity. For systems initially in the SF regime (g = g 2 ), the fidelity curves exhibit small oscillations (see figure 6(e)). These can be understood from time-dependent perturbation theory as oscillations occurring when some of the frequencies involved in the Fourier decomposition of the perturbation Hamiltonian enter into resonance with system frequencies. This is expected since SF have a dense spectrum at low energies and are therefore likely to enter into resonance with one of the frequencies of the perturbation Hamiltonian [30]. Hence, the amplitude of the oscillations in the fidelity curve associated with a ramping strategy s 1 (t) should be smaller than those associated with a ramping strategy s 2 (t) if A(1, τ Q ) < A(2, τ Q ) for all τ Q . This is what is observed from our numerical simulations (see figures 6(f)-(h)).
The quasi-momentum distribution of the particles in the a-mode at the end of the different ramping processes for a system with g = g 1 are shown in figure 8. For the linear ramp, the quasimomentum distribution of particles shown in figure 8(a) does not correspond to that of a MI state for the quench times considered. In contrast, for the ramp s gap the quasi-momentum distribution becomes approximately flat for quench times of τ Q > 200/E R . Even for the fastest ramps, we find that the occupation of the b-mode is less than 2% of the total number of particles. Thus, the experimental measure of the register fidelity can be made by comparing the quasi-momentum distribution of particles in the final state with, e.g. that shown in figure 5(c).

Discussion of the numerical results.
In the BHM with U/J 1, the tunnelling time 1/J determines the adiabatic timescale of the system. However, as soon as many-body interactions are sufficiently large, this timescale often becomes very non-adiabatic [31]. The main observation that can be drawn from the numerical calculations presented in the last section is that by preparing the system in a TG state (g = g 1 ), and using an efficient ramping strategy, it is possible to initialize a very deep MI state on a timescale equivalent to half the tunnelling time in the large lattice limit (see figure 6(b)). The time required to initialize a MI state with unit filling as well as a TG state with half filling is approximately ten times the tunnelling time of the final system [15,16,31,32]. Since the tunnelling time in the large lattice limit is two orders of magnitude larger than in the small lattice limit, the total time required to initialize a MI state using our procedure is an order of magnitude faster than the direct quantum-freezing method. In this estimation we assume that the initial BEC has zero temperature, i.e. we do not take the effect of defects present in the initial state into account.
The experimental realization of initial states with g = g 1 and g = g 2 can be achieved using Feshbach resonances. For a magnetic Feshbach resonance fluctuations in the magnetic field result in fluctuations of the size of the gap between the ground and the first excited state which will affect the performance of our scheme. For instance, in the case g = g 1 magnetic field fluctuations of 10 mG will change the gap by approximately 0.5% for 85 Rb or 133 Cs atoms [33,34]. We assume adiabatic evolution of the system and thus these fluctuations will have negligible repercussions on the fidelity of the final state. To realize the SF regime with g = g 2 very stable magnetic fields are required. Magnetic field fluctuations of e.g. 1 mG will lead to gap fluctuations of approximately 1% in 23 Na and 85 Rb. We finally remark that our scheme could also be used without employing a Feshbach resonance. This case corresponds to an intermediate value of g between g 1 and g 2 . While a detailed analysis of the intermediate regime is beyond the scope of the present work, we do not expect qualitative differences compared to the interaction strengths considered here.

Analytical results
In this section we derive an approximate expression for the quench time required to obtain a given fidelity in the case of a linear ramp.
The energy spectrum of systems in the TG and SF regime is gapless 6 (see e.g. [35]), while in the MI regime, the gap between the ground and first excited state is proportional to the on-site interaction energy [30]. Since the relaxation time τ(t) of the system-the time required by the system to adjust to a change of parameters at time t-is inversely proportional to the gap between the ground and the first excited state, the relaxation time in the small lattice limit is large, while it is small and finite in the large lattice limit. This observation suggests that the adiabatic-impulse (AI) assumption from KZ theory can be used to evaluate the adiabaticity of a ramp with respect to the quench time [36]- [38].
The AI approximation is based on the following considerations: (i) when the gap between the ground and the first excited state is large, the relaxation time of the system is short and thus a system starting its evolution in the ground state remains in the ground state, i.e. its evolution is adiabatic. (ii) when the gap between the ground and the first excited state is small, the system's relaxation time is large and the system no longer adapts to changes of the Hamiltonian's parameters and its state becomes effectively frozen. The system is then in the impulse regime. The instantt at which the system passes from the impulse to the adiabatic regime, and inversely, is defined by the equation [36,37,39] where α = O(1) is a constant. Note thatt is a time and not an operator. In the AI approximation, the time-evolution of the system is either adiabatic or impulse. Thus, the density of defects D, which corresponds to the density of excitations caused by a change of parameters which drive the system from the impulse to the adiabatic regime can be approximated by [37] where | g (0) and | e (t) are the ground and first excited states at the initial time t = 0 and at time t =t, respectively. Hence, without solving the time-dependent Schrödinger equation it is possible to make predictions for the density of defects (24) resulting from a given dynamical process.
In order to apply the KZ theory to our problem, we develop an effective model describing the system dynamics. A similar model was recently examined by Cucchietti et al [38]. We find from numerical calculations (see figure 7(b)) that most of the excitations created in the system are caused by transitions from the ground to the first excited state. Furthermore, we examined the form of the eigenvectors of the Hamiltonian (10) in both limits for different system sizes via exact diagonalization. This revealed that both the ground state and the first accessible excited state can be approximated by an expansion in only two basis states. The elements of this reduced basis set are given by |2 = (|2; 0; 1; · · · ; 1 + |2; 1; 0; 1; · · · ; 1 + |2; 1; 1; 0; 1; · · · ; 1 + |1; 2; 1; 0; 1; · · · ; 1 + · · ·)/ M(M − 1), where, e.g. |2; 0; 1; · · · ; 1 = |2 The basis state |2 corresponds to a superposition of all possible states of a system of M particles in the a-modes with (M − 2) singly occupied sites, one doubly occupied site, and one empty site. In the limit M → ∞, the matrix representationĤ R of Hamiltonian (10) in the reduced basis {|1 , |2 } reads, up to a constant energy V â The instantaneous eigenstates of equation (27) associated with the energies of the ground and first excited levels are given by |g(t) = − sin(θ(t)/2)|1 + cos(θ(t)/2)|2 and |e(t) = cos(θ(t)/2)|1 + sin(θ(t)/2)|2 , respectively, with cos(θ(t)) = −U aa (t)/[U aa (t) 2 + 8J aa (t) 2 ] 1/2 , θ ∈ [0, π]. Furthermore, we approximate the parameters J aa and U aa as linear functions of time where J aa = |J aa (s = 1) − J aa (s = 0)|, U aa = |U init − U aa (s = 0)| and U init = min[U aa ]. We further simplify the HamiltonianĤ R by replacing the hopping term J aa (t) by its time average. Setting J aa (t) =J aa withJ aa = (1/τ Q ) τ Q 0 dt J aa (t) and rescaling the time as t → t − (U init τ Q / U aa ) yields the transformation which turnsĤ R into the Landau-Zener formĤ R =Ât +B, whereÂ andB are Hermitian matrices andÂ is diagonal [40]- [42]. The energy spectrum of the Hamiltonian (29) reproduces approximately the features of the spectrum ofĤ eff except that the gap is overestimated in the small lattice limit. Starting at t = 0 in the small lattice limit where the system is impulse, we use the AI approximation to derive the density of defects at the end of a linear ramp which drives the system to the large lattice limit, where the system is adiabatic. Defining the relaxation time as the inverse of the energy gap 1/ E between the levels ofĤ R , with E = [( U aa t/τ Q ) 2 + (2 √ 2J aa ) 2 ] 1/2 , equation (23) can be solved analytically and the instantt at which our system exits the impulse Analytical formula for Analytical formula for regime is given bŷ where η = 4J 2 aa / U aa . We evaluate the density of defects D using equation (24) with | g (t i ) = |g(0) and | e (t) = |e(t) . In order to simplify the expression of D, we have set θ(0) = π/2, which is a greater value than the one we would obtain using the real parameters of the system. This has no significant physical consequences in our case as it is actually equivalent to considering a smaller energy gap in the small lattice limit. The density of defects is then given by whereθ = θ(t). Inverting equation (31), we obtain which gives an approximate analytical expression of the quench time τ Q required to obtain a density of defects D.
In our system, defects correspond mainly to doubly occupied sites. Thus, the number of defects is approximately measured by the operator Hence, the density of defects is given by D = K /M = (n a i ) 2 − n a i . Numerical calculations show that the density of defects D has the same scaling behaviour as n 2 a,i . A comparison between equation (31) and the numerical values of the density of defects in the large lattice limit for different values of the ramping time τ Q is shown in figure 9. For M = 12 particles and g = g 2 , we find that the analytical formula for D fits the numerical data well for α = 1.41. For g = g 1 the fit is less accurate. This is expected since for this value of g, the system has a less distinct separation between the adiabatic and impulse regime than for g = g 2 . For the number of particles we have been able to simulate, the fit improves as we increase the number of particles for both values of g. In addition to this, we see from equation (32) that the time required to initialize a register with a given fidelity-and thus the adiabatic time for small D-scales with the ratio U aa /J 2 aa .

Conclusion
We have shown that the dynamics of an optical lattice whose periodicity is doubled via superlattice potentials is very well described by a two-mode Hubbard-like Hamiltonian. The parameters of this Hamiltonian have been evaluated in the TB approximation using optimally localized GWFs. The doubling of the period removes half of the lattice sites and doubles the filling factor. We have shown that this doubling can be used for the fast initialization of a quantum register. By starting from a half filled lattice in the small lattice limit filled by either a TG (g = g 1 ) gas or a SF (g = g 2 ), a commensurate MI state corresponding to an atomic quantum register can be obtained on timescales shorter than those achieved by direct quantum freezing of a SF with same lattice spacing. Furthermore, we derived an analytical expression for the density of defects as a function of the quench time for linear ramping of the superlattice. We found that the time required to achieve a given density of defects is proportional to the ratio U aa /J 2 aa . Our numerical calculations of ground state properties suggest that doubling the lattice period drives the system through a quantum phase transition for large lattices M → ∞. The eventual abrupt change in the ground state properties might be observable by time-of-flight measurements. An investigation of whether such a quantum phase transition indeed exists is beyond the scope of the current work but will be investigated in future work.
In this work we concentrated on the transition from filling factor n = 1/2 to n = 1. We finally note that the idea developed in this paper may be extended by considering lattices with an initial filling factor of n = 1/2 (where is an integer). Subsequently removing every second barrier will create a lattice with period 2a and filling factor 1/2 −1 . This procedure could be repeated times providing a lattice with filling factor n = 1 and large lattice spacing 2 a.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The elements of the matricesε(µ) are given by [43] ε αβ (µ) = w α,i |ĥ 0,x |w β,i+µ . (A.2) For µ = 0 and 1, the elements of the matricesε(µ) correspond to the local site energy and hopping matrix elements between neighbouring sites, respectively. For the TB approximation to be accurate, we need the eigenvalues E α,q ofĤ 0,q to reproduce very closely the band structure of the exact single particle Hamiltonianĥ 0,x for all values of the lattice profile s. If we were to use WFs as mode functions w α,i , the matricesε would be diagonal [24] and the dispersion relations of the two Bloch bands sinusoidal. In order to obtain a more accurate description, we use GWFs as mode functions. The definition of GWFs is given by β,q is called the generalized-Bloch orbital with α,q the Bloch function associated with the band α and R i is the centre of site i [25,44]. The rows of the 2 × 2 matrix U (q) contain the (real) normalized eigenvectors ofĤ 0,q associated with the eigenvalues E α,q , that is µ=a,b (Ĥ 0,q ) n,µ U (q) α,µ = E α,q U (k) α,n [43]. Inserting GWFs in equation (A.2), we recover the elements ε αβ (µ) and, hence, GWFs correspond to the mode functions associated with the effective HamiltonianĤ 0,q [43]. Notice that the definition of GWFs reduces to that of WFs for U (q) = 1 1.

A.1. Localization properties of GWFs
Given a valid set of GWFs, another equally valid set of GWFs can be obtained by applying the following transformation on the U (q) matrices where φ α (q) are (real) functions of q which can be chosen freely as long as they do not introduce discontinuities in the generalized Bloch function [45]. The gauge transformation (A.4) is equivalent to re-phasing each Bloch function as α,q → e iφ α (q) α,q . Notice that gauge transformations do not affect the value of the parameters V α and J α,β calculated using the relations (7) and (6), respectively, but they alter the localization properties-the spread-of the GWFs. Following the convention suggested by Blount [45], we set the phase functions φ α (q) in a manner that leads to maximally localized GWFs. That is, we choose the phase functions such that the resulting GWFs minimize the spread functional where in our case α = a, b while x 2 α = w α,i |x 2 |w α,i and x α = w α,i |x|w α,i correspond to the centre of a GWF and its second moment, respectively. where K j = 2πj/L with L = Ma. Invoking the translational symmetry of the lattice and the convolution theorem [45], the value of the functional is minimized when the expansion coefficients G α,j (q, s) are chosen real, which is always possible when the lattice possesses mirror symmetry [43]. This is equivalent to the choice of purely real GWFs for even generalized-Bloch functions and purely imaginary ones for odd generalized-Bloch functions. Notice that this conclusion is in agreement with a conjecture of Marzari et al [25] on the real nature (up to a global phase) of maximally localized WFs.
We have numerically evaluated the phases φ α (q) using the algorithm described in [25] for the special case of 1D WFs. This procedure minimized the functional in the limit of very fine sampling of the q-space. The effect of this localization procedure is illustrated in figure 3(a).
zone (see figure 2). The numerical values of the parameters obtained via this procedure are shown in figure 4(a).
The accuracy of the TB approximation can be tested by evaluating the standard deviation between the exact and approximated band structure. That is, taking N q different points q i on each band, we define the standard deviation between the exact and approximate band structure for a given lattice profile by σ 2 s = (1/2N q ) N q i=1 ( ε 2 a,q i + ε 2 b,q i ), with ε α,q i = ε α,q i − E α,q i . Averaging over N s different lattice profiles s i , we obtain σ s = [(1/N s ) N s i=1 σ 2 s i ] 1/2 = 3.4 × 10 −2 E R for a lattice depth of V 0 = 10 E R . This excellent agreement improves further as we increase the value of V 0 , and hence fully justifies the TB approximation.
operations which scales as O(d 4 χ 3 ). The resulting tensors are then systematically truncated back to a maximum rank of χ.
Dynamical simulations can be performed by decomposing the time evolution operator exp(−iHδt), for small time step δt, into a sequence of pairwise unitaries via a Suzuki-Trotter expansion. Given the properties outlined such a calculation is likely to be accurate for a practical value of χ if both the initial state and the states generated by the dynamics remain in the lowenergy manifold of the system. To determine the appropriate χ calculations are repeated with increasing values of χ until the final result converges and are unaffected by further increases. For practical purposes the convergence is usually quantified by the robustness of the expectation values calculated. The accuracy of a calculation is also gauged by the sum of the discarded Schmidt coefficients at each time step-a quantity which should necessarily be small-and the deviation of normalization of the final state from unity which indicates the accumulated effect of truncation.
Finally, initial states are typically taken to be the ground state of the system which are found either by applying the DMRG procedure or, as in this work, by simulating imaginary time evolution through the repeated application of exp(−Hδt) and subsequent renormalization of the state. In our simulations, we have used χ = 40.