Bloch oscillations in a Bose-Hubbard chain with single-particle losses

We theoretically investigate Bloch oscillations in a one-dimensional Bose-Hubbard chain, with single-particle losses from the odd lattice sites, described by the Lindblad equation. For a single particle the time evolution of the state is completely determined by a non-Hermitian effective Hamiltonian. We analyse the spectral properties of this Hamiltonian for an infinite lattice and link features of the spectrum to observable dynamical effects, such as frequency doubling in breathing modes. We further consider the case of many particles in the mean-field limit leading to complex nonlinear Schr\"odinger dynamics. Analytic expressions are derived for the generalised nonlinear stationary states and the nonlinear Bloch bands. The interplay of nonlinearity and particle losses leads to peculiar features in the nonlinear Bloch bands, such as the vanishing of solutions and the formation of additional exceptional points. The stability of the stationary states is determined via the Bogoliubov-de Gennes equation and is shown to strongly influence the mean-field dynamics. Remarkably, even far from the mean-field limit, the stability of the nonlinear Bloch bands appears to effect the quantum dynamics. This is demonstrated numerically for a two-particle system.


Introduction
In recent years there has been great interest in the dynamics of ultracold atoms and Bose-Einstein condensates (BECs) trapped in optical lattices with a static tilt [1][2][3][4]. The application of a tilt can cause the atoms to undergo Bloch oscillations, as has been demonstrated, for example, in remarkable experiments with one and two bosonic atoms [5]. A BEC of weakly-interacting particles also exhibits Bloch oscillations, although this behaviour can dramatically change when the particle interaction strength increases [6]. In a full many-particle treatment, a system of ultracold atoms in a tilted one-dimensional optical lattice with M = 2L + 1 lattice sites can be described by the Bose-Hubbard Hamiltonian Theâ † j (â j ) create (annihilate) particles in the ground state of the j th site and satisfy the canonical commutation relations [â j ,â † k ] = δ jk , [â j ,â k ] = 0; (2) arXiv:2002.12727v1 [quant-ph] 28 Feb 2020 J is a positive coupling constant; U describes pairwise on-site interactions; and F describes the lattice tilt. Throughout this paper we use appropriate dimensionless units with J = 1 and = 1. Open quantum systems have become a focus of intense research, motivated in part by the drive towards engineering quantum effects into new technologies. While dissipation and decoherence were once viewed solely as unwanted features, engineered dissipative systems have opened up new opportunities for the control of quantum dynamics. For example, the ability to remove single atoms from specific sites of an optical lattice [7,8] may be used to control the dynamics of many-body quantum systems. Theoretically, the single-particle losses in such a setup can be described by the Lindblad equation [9][10][11][12][13], an approximate master equation that is often valid in atomic, molecular and optical systems.
The time evolution given by the Lindblad equation can be related to non-Hermitian Schrödinger dynamics in various ways. For example, the quantum jump approach consists of an effective non-Hermitian evolution interrupted by irreversible transformations [14,15]. Effective non-Hermitian Hamiltonians have long been used to model open quantum systems, and constitute an active research area in their own right. The class of non-Hermitian PT-symmetric Hamiltonians, the eigenvalues of which are real or else come in complex conjugate pairs, exhibits interesting properties and has attracted much attention over the last few decades [16]. There has been a growing interest in non-Hermitian PT-symmetric quantum dynamics, with many experiments performed in analogous wave systems, such as microwave resonators and optical waveguides [17]. The current consensus is that it is not possible to implement PT-symmetry in a bona fide open quantum system. This would require a perfect balance between gain and loss, which cannot be achieved due to the inevitable noise that accompanies gain. However, passively PT-symmetric systems can be implemented, which are equivalent to PT-symmetric ones up to an overall exponential decay. The first experimental realisations of passively PT-symmetric quantum systems have recently been acheived [18][19][20].
Cold atoms in optical lattices provide a natural platform for investigating PTsymmetry in many-body systems. Here we study the interplay of Bloch oscillations with PT-symmetric particle losses in a Bose-Hubbard chain. Specifically, we consider the effects of localised single-particle losses from the odd sites of the lattice, which can be described by the Lindblad equation with Lindblad operators of the form L j ∝â j [11,13]. We assume that the decay rate γ is the same on each dissipative lattice site and takeL j = √ 4γâ j , where the factor of four is included to simplify subsequent expressions. The Lindblad equation then has the form where the effective Hamiltonian is defined aŝ withĤ being the Bose-Hubbard Hamiltonian (1). The effective Hamiltonian has passive PT-symmetry, where the parity operator interchanges even and odd sites and the time-reversal operator is simply the complex conjugation. Some aspects of a single-particle version of this model have been investigated in [21][22][23][24]. In the mean-field limit of many particles this model leads to a complex discrete nonlinear Schrödinger equation. Nonlinear PT-symmetric systems can have properties quite unlike their linear counterparts. They appear naturally in models of BECs in PTsymmetric potentials [25][26][27], as well as in the context of classical wave systems [28] and nonlinear optics [29,30]. The model arising here could also be experimentally realised in coupled wave guides with absorption and Kerr nonlinearity.
In what follows we shall analyse the combined effect of particle losses and interaction/nonlinearity on Bloch oscillations. We start by considering the single particle case in Section 2. In this case the time evolution of the stateρ is completely determined byĤ eff . The spectral properties ofĤ eff are analysed for an infinite lattice and features in the spectrum are linked to observable dynamical effects. In Section 3 we consider the mean-field limit of a many-particle system, and derive analytic expressions for the generalised nonlinear stationary states and the nonlinear Bloch bands. The interplay of nonlinearity and particle losses is found to lead to peculiar features in the nonlinear Bloch bands, such as the vanishing of solutions and the formation of additional exceptional points. The stability of the stationary states is determined via the Bogoliubov-de Gennes equation and connected to features in the mean-field dynamics. Surprisingly, the stability of the nonlinear Bloch bands appears to influence the quantum dynamics even far from the mean-field limit. This is demonstrated numerically for a two-particle system. We end with a brief summary in Section 4.

State evolution
Let us first consider the time evolution of the quantum stateρ for a single particle. In this case the time evolution ofρ is completely determined by the effective Hamiltonian (4). To see this, consider the expansion of the density operator in the basis {|j , |v }, where |j is a Fock state corresponding to a single particle on the j th lattice site and |v is the vacuum state. Inserting this expansion into the Lindblad equation (3) yields the equations of motion for the matrix elements By defining the components of the vector ψ via ρ jk = ψ jψk it follows from the first equation above that the evolution of ψ is governed by the non-Hermitian Schrödinger equation where H eff is the effective Hamiltonian (4) represented in the basis {|j }. In what follows it is convenient to express the single-particle effective Hamiltonian in terms of the basis states {|j } aŝ The non-Hermitian Schrödinger equation (7) yields the equations of motion for the components For a single particle state we have ρ jv (0) = 0 and thus ρ jv (t) = 0 for all time. Therefore, the time evolution ofρ is completely determined by ψ, where the expression for ρ vv (t) follows from the trace conservation of the Lindblad equation.

Spectral features
Let us now examine the spectrum of the effective Hamiltonian (8) with an infinite lattice. We will use algebraic methods to show that the spectrum generally consists of two ladders of eigenvalues. The analysis in this section closely follows [31], in which a Hermitian one-dimensional Bose-Hubbard Hamiltonian with a variation of the on-site energy of every other site is analysed. Consider the complex conjugation operatorT and the translation operatorŜ m = j |j − m j|, which shifts the lattice by m sites to the left. A short calculation shows that this pair of operators and the Hamiltonian satisfy the commutation relations Ŝ m ,T = 0.
Let |λ be an eigenvector ofĤ eff with the eigenvalue λ (F, γ), that is, To simplify the notation we temporarily drop the functional dependence of λ(F, γ) on the system parameters. Using the commutation relations above it can be shown that where l ∈ Z. An application of the operatorsŜ 2l andŜ 2l+1T to the eigenvector |λ then yields the ladders of eigenvectorŝ where α = 0, 1 labels the ladder, and the states are associated with the energies E 0,n = Re λ + 2nF + i Im λ, The expressions for these energies can be simplified by noting that the effective Hamiltonian possesses a chiral symmetrŷ whereX is a unitary and Hermitian operator defined aŝ The chiral symmetry ofĤ eff , together with the (passive) PT-symmetry, implies that the eigenvalues λ have real parts symmetric around Re λ = 0, and complex conjugate imaginary parts (up to an overall imaginary energy shift). Note that chiral symmetries in non-Hermitian systems have attracted considerable attention recently, see, e.g., [32][33][34][35]. By acting on the eigenvalue equation (16) withX one can show thatX|E α,n is an eigenvector ofĤ † eff with eigenvalue −E α,n and the spectrum ofĤ † eff may be written as modulo 2F . However, the eigenvalues ofĤ † eff can also be put into correspondence with the complex conjugate eigenvalues ofĤ eff , that is, modulo 2F . Equality of the sets (22) and (23) implies that Re λ = 0, and the energies in (18) and (19) simplify to Thus, the spectrum ofĤ eff generally consists of two energy ladders. The real part has values centred around Re E = 0, with equal spacings of size F , while the imaginary part of each ladder is fixed and independent of n. The real and imaginary parts of the spectrum are shown in Figure 1 in dependence on γ for two different values of F . The eigenvalues were obtained from numerical simulations in a lattice of size M = 201 and eigenvalues due to the finite lattice size were discarded.
In the unitary case (γ = 0) it is well known that the dynamics are periodic with Bloch period T = 2π/F (see, e.g., [36]). From the above spectral analysis we see that when γ = 0 the differences in the real parts of the energies are now given by 2F within each ladder. But the ladders are stacked such that the resulting energy difference of the ordered real parts is F , just as in the unitary case. This gives rise to periodic dynamics with the same Bloch period T = 2π/F . However, due to the difference in the imaginary parts of the two ladders, the components from one ladder may decay faster than those from the other. This can lead to the emergence of a motion dominated by only the stable ladder, which has half the period. For example, when γ = 0 a single particle initialised on a single lattice site |Ψ(0) = |j performs left-right symmetric Bloch oscillations with a time period T [36]. When γ = 0 the decay of one of the ladders eventually leads to an effective frequency doubling of the system, resulting in an additional structure with half the period. This effect is illustrated in Figure 2, where the renormalised density |ψ j | 2 /|ψ| 2 on each lattice site j is shown for both the unitary and dissipative time evolution.

Dispersion relation
In contrast to the breathing motion of a particle initialised on a single lattice site, a broad initial beam performs spatial Bloch oscillations, an example of which is depicted in the left panel in Figure 4. In order to understand this behaviour we first derive the dispersion relation for the field-free (F = 0) effective Hamiltonian (8) with an infinite lattice. To this end, we consider the effective Hamiltonian (8) in the basis of Bloch states where the quasimomentum k is confined to the region −π ≤ k ≤ π. The Bloch states are orthogonal and normalised to the 2π-periodic delta comb k|k = δ 2π (k − k). The effective Hamiltonian is not diagonal in this basis and for F = 0 one finds the equations of motion Following [37] we introduce the two-component function The time evolution can then be written as the twocomponent Schrödinger equation where σ i are Pauli matrices and the Bloch Hamiltonian is defined as The k-dependent eigenvalues of h(k) define the dispersion relation of the two-band system [21,22,24] E ± (k) = −iγ ± 4 cos 2 k − γ 2 .
When γ < 2 there are exceptional points at 2| cos k| = γ, and for values of the quasimomentum 2| cos k| > γ the energy has a real part. At the exceptional points the passive PT-symmetry breaks and for 2| cos k| < γ the energy becomes purely imaginary. The band structure is illustrated for an example of a small γ = 0.05 in Figure 3. For γ ≥ 2 the bands are purely imaginary for all values of the quasimomentum.

Beam dynamics
Consider an initial Gaussian wave packet where x 0 is the centre of the Gaussian; k 0 is the initial momentum; σ is the width parameter; and N is a normalisation constant. In [24] it was shown that a broad Gaussian beam in position space is an approximate eigenstate of the Bloch Hamiltonian (29), provided the Bloch states are close to the standard basis vectors. Furthermore, a static tilt can result in Bloch oscillations, together with a splitting of the beam in real space due to a transfer of population between the Bloch bands. This is illustrated in the left frame of Figure 4 for γ = 0.05 and F = 0.2, where the Gaussian wave packet effectively maps out the real part of the band structure depicted in Figure  3. Numerical simulations show that this behaviour does not require the decay rate to be exactly the same on each dissipative site. This can be seen in the middle panel of Figure 4, where the decay rates on the dissipative sites γ j are randomly selected from the interval (0.025, 0.125). Thus, a robust experimental realisation should be feasible.
Let us now attempt to gain some insight into the oscillatory behaviour of the beam in position space, using a quasiclassical argument similar to [31,36,38]. In the present case there are two energy bands that degenerate at the exceptional points. Therefore, in a semiclassical treatment we can consider the dynamics on each of the branches of the band structure separately, and account for transitions between the bands using a Landau-Zener type formula. The right frame shows the total particle number dynamics for γ = 0.05 (black solid) and the randomised decay rates (red dashed).

The application of a static tilt F introduces a term F x into the Bloch Hamiltonian
where where (q, p) ∈ R × R are a pair of canonical phase-space coordinates. The eigenvalues of (33) are readily found to be where E ± (p) is the dispersion relation (30), i.e., E ± (p) are the eigenvalues of the fieldfree Bloch Hamiltonian (29). In the single-band approximation the population in each band ± is assumed to evolve independently according to the semiclassical dynamics generated by the non-Hermitian Hamiltonian ξ ± (q, p). Thus, as the coupling between the bands is ignored, this approximation cannot account for any population transfer between the bands. Note that the q and p are different in the ± bands, so that technically we should write ξ ± (q ± , p ± ). However, for notational simplicity we do not explicitly indicate this, as it should be clear from the context which band is being referred to. The semiclassical dynamics of a Gaussian wave packet evolving under a non-Hermitian phase-space Hamiltonian H(q, p) = Re H + i Im H are given by [39] where Σ encodes the (co)variances of the canonical coordinates according to and the determinant of Σ is one. The covariances are also time dependent, following the dynamical equationṡ where Ω is the symplectic matrix and H is the Hessian of H. In our approximation we have assumed that the semiclassical non-Hermitian dynamics above can be applied to each band individually. A more careful treatment would derive the semiclassical limit of non-Hermitian Diractype equations, thus extending the semiclassical dynamics identified in [40] for the unitary case. This is an interesting problem for future investigations. We must distinguish between the regions where the square root 4 cos 2 p − γ 2 appearing in E ± (p) is real or imaginary. In the region where it is real we find and an application of the semiclassical equations (35) and (36) yieldṡ Therefore, the acceleration theorem p(t) = −F t + p(0) is recovered and Similarly, in the region where the square root is imaginary we have The general semiclassical equations of motion then yielḋ andΣ qq = Σ 2 qp − 1 ∂ pp Im ξ ± ,Σ qp = Σ qp Σ pp ∂ pp Im ξ ± ,Σ pp = Σ 2 pp ∂ pp Im ξ ± , (45) where Im ξ ± is defined in (43). Now observe that in the approximation of an initially infinitely narrow momentum wave packet with Σ pp (0) = 0, the (co)variances remain constant for all times, i.e., Σ pp (t) = 0 and Σ qp (t) = Σ qp (0). In this case the equations of motion simplify toq = Σ qk (0)∂ k Im ξ ± ,k = −F, and the acceleration theorem p(t) = −F t + p(0) is again recovered. In the special case where Σ qp (0) = 0 we find that q(t) = q(0). The position of each wave packet follows the real dispersion relation in the unbroken passive PT-symmetry phase (42) and remains stationary in the broken phase. Therefore, as observed in Figure 4, the position maps out the real part of the band structure of the Bloch Hamiltonian. Close to the exceptional point, at which the two branches of the band structure intersect, a partial transition to the other branch of the band structure takes place. This can be described by a Landau-Zener-type scenario. Details of the calculation can be found in [24]. Here we simply quote the estimate of the transition probability between the bands, which fits the numerical results with surprising accuracy for a large range of parameters. This phenomenon has already been observed experimentally, for example in PT-symmetric photonic lattices [41].

Mean-field limit of many particles
Let us now consider the time evolution of a pure BEC, defined as, Here |v is the vacuum state, N 0 is the number of particles and the coefficients ψ j are not necessarily normalised to unity, The many-body state (48) is then normalised to Ψ(0)|Ψ(0) = N N0 . A far richer set of dynamics are expected to be possible for many-particle systems, owing to the interplay between particle interactions and dissipation. However, for many particles on a large lattice the size of the Hilbert space prohibits numerical investigations. Thus, to obtain insights into the case of large particle numbers we consider the meanfield limit N 0 → ∞. For unitary time evolution it is well known that the mean-field dynamics provide a good approximation, provided that the state remains close to a pure BEC for all time [1]. For a large number of weakly-interacting particles, the meanfield approximation is also valid for dissipative processes described by the Lindblad equation (see, e.g., [42]). One way to construct the mean-field approximation is by working with the elements of the single-particle density matrix (SPDM) σ ij = â † iâ j /N 0 . Due to the interaction term, the equations of motion for the SPDM elements do not form a closed system of equations. Rather, they form the so-called Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy. In the mean-field approximation the hierarchy is truncated by approximating the four-point correlation functions as The mean-field evolution of the SPDM elements is then found to be where g = U N 0 is fixed as N 0 → ∞. By identifying σ jk =ψ j ψ k the mean-field equations may equivalently be described by the complex discrete nonlinear Schrödinger equation The initial conditions ψ j (0) are the parameters ψ j appearing in the pure BEC state (48). The same equations of motion can also be obtained directly from a general dissipative phase-space dynamics derived from Gaussian states in [43]. We note that (52) could also be implemented in classical optics using an array of waveguides with absorption and Kerr nonlinearity (see, e.g. [29,30]).

Nonlinear Bloch bands
For a single particle we found that a knowledge of the band structure was crucial for understanding the Bloch oscillations of a broad beam. In the mean-field limit it is possible to define nonlinear Bloch states and their corresponding nonlinear generalised eigenvalues. In closed systems it has been shown that the presence of particle interactions can lead to novel features in the nonlinear Bloch bands. For instance, the formation of loops at the band edges [44] for a Bose-Einstein condensate in a bichromatic lattice [45]. Here we shall examine the nonlinear band structure when the odd lattice sites are dissipative, by deriving a nonlinear complex Bloch Hamiltonian, and following the procedure in [25] to deduce its nonlinear eigenvalues. The nonlinear Bloch states are defined as the stationary states φ j of the field-free version of equation (52) where µ is a generalised nonlinear complex eigenvalue. The translational symmetry motivates the ansatz for the stationary states, with the quasimomentum −π ≤ k ≤ π and the normalisation convention |A| 2 + |B| 2 = 1. The coefficients A and B are then determined by the nonlinear two-level system To simplify the calculations we subtract an overall complex energy g 2 − iγ from the diagonal, leading to the PT-symmetric system In [25] it was shown that the eigenvalues of a nonlinear complex two-level system of the form (56) can be obtained analytically. Here we closely follow the procedure in [25] to obtain the nonlinear Bloch bands. Multiplying the first equation byĀ, the second byB, and making use of the normalisation condition |A| 2 + |B| 2 = 1, yields the generalised eigenvalues For k = ± π 2 equation (58) reduces tõ and (56) provides the two conditions In this case the only solutions are z = ±1, corresponding to (A = 1, B = 0) and (A = 0, B = 1) with generalised eigenvalues µ = g and µ = g − 2iγ respectively. When k = ± π 2 it can be shown that the possible values of the population imbalance z are given by the real roots of the polynomial which has the four solutions Let us now parameterise the nonlinear Bloch states as Following the analysis in [25] it can be shown that when k = ± π 2 there are up to four possible solutions with where the phase coordinate for the second solution in (64) is given by π/2 − v. The nonlinear eigenvalues corresponding to the two solutions in (64) are given by respectively. Those corresponding to the two solutions (65) may be written in the form µ = 2c 1 − 2 cos 2 k c 2 + γ 2 (1 + sgn(cos k)) + iγ(z − 1).
In the latter expression the real part of µ is constant when cos k is negative, resulting in nonlinear Bloch bands that are 2π-periodic. This should be contrasted with the non-interacting (g = 0) Bloch bands (Figure 3), which were found to be π-periodic.
Only pairs (z, v) that are real valued correspond to solutions. Thus, there are only two solutions for |2 cos k| < γ and two for 4 cos 2 k > c 2 + γ 2 . There are three solutions at the exceptional points 2| cos k| = γ, where the two solutions (64) coalesce in an exceptional point (EP2). The nonlinearity introduces additional exceptional points at the quasimomentum values 2| cos k| = c 2 + γ 2 . At these points there can be three or two solutions. There are three solutions for sgn (cos k) > 0, as the two solutions in (65) coalesce (EP2). Surprisingly, there are only two solutions for sgn (cos k) < 0, as both the solutions in (65) coalesce with one solution in (64) (EP3). Otherwise there are four solutions. The existence of an EP2, EP3 and even an EP4 has previously been demonstrated for a BEC in a PT-symmetric double well potential [26,27,46].
The real and imaginary parts of µ are plotted in Figure 5 for weak (g = 2) and strong (g = 7) interaction strengths, when the decay rate is γ = 0.1. The Bloch bands for the non-interacting (g = 0) case are also shown in the top row for comparison. The combination of interaction and dissipation leads to significant changes in the band structure compared to the single particle (non-interacting) case. When g is switched on a pair of solutions vanish just after the exceptional points at 2| cos k| = γ, visible around k = ±π/2, and a new band emerges (coloured blue). The interaction introduces additional exceptional points at 2| cos k| = c 2 + γ 2 that vanish when c 2 + γ 2 > 2. Away from exceptional points the solutions in the blue band have the same real part but different imaginary parts. The solutions in the black band have different real parts but degenerate imaginary parts.

Stability analysis
In contrast to the Bloch states of the single-particle system, the nonlinear Bloch states can become dynamically unstable. The stability of the stationary states can be determined by examining the effect of a small time-dependent perturbation of the stationary solution of the form [45] ψ j (t) = e −iµt φ j + e ikj−iµt e −iωt δχ A − + e iωt δχ A + , j even, e ikj−iµt e −iωt δχ B − + e iωt δχ B + , j odd.
By inserting the ansatz (68) into the field-free version of (52) the problem reduces to the two-mode system analysed in [25], and the stability analysis obtained therein can immediately be applied. In particular, it is sufficient to consider the nonlinear Schrödinger equation where H nl (|Ψ 1 | 2 , |Ψ 2 | 2 ) is the two-level nonlinear Hamiltonian In the previous section we found that the stationary states of this system are given by χ = (A, B), with corresponding nonlinear eigenvalues µ. Therefore, the stability analysis (68) reduces to analysing the linearised effect of a small time-dependent perturbation to the stationary solution χ of the form where the components of δχ ± = (δχ A ± , δχ B ± ) are defined in (68). In this case it can be shown [25] that the Bogoliubov-de Gennes equations yield the eigenvalue equation where M is the 4 × 4 matrix Here H 0 nl is the nonlinear Hamiltonian (70) evaluated at the stationary solution χ = (A, B) and P is the projection matrix orthogonal to χ, The imaginary parts of the eigenvalues ω determine the stability of the stationary solution χ. When Im ω ≤ 0 the perturbations are decaying or oscillating, and the solution is stable. On the other hand, if Im ω > 0 the perturbations grow exponentially and the solution is unstable. A short calculation reveals that at k = ±π/2 the solution χ = (1, 0) with µ = g is stable, while the second solution χ = (0, 1) with µ = g − 2iγ is unstable. In general the stability of the stationary solutions changes as the system parameters are varied. This is illustrated in Figure 6 for two different values of g and varying quasimomentum values, with γ = 0.1.
It follows that for sgn(cos k) > 0 the solution corresponding to µ + becomes unstable for quasimomentum values satisfying 4 cos 2 k < c 2 + γ 2 , while the solution corresponding to µ − remains stable. The behaviour of these solutions is reversed when sgn(cos k) < 0, as illustrated in the top row of Figure 6. The stability of these solutions is completely determined by the interaction strength g. On the other hand, changing g and/or γ can change the stability of the solutions (65). The change in stability as g is increased is shown in the bottom row of Figure 6. It is important to note that the stability obtained from the Bogoliubov-de Gennes equation can change during the time evolution of the system. Due to particle losses the norm (|A| 2 + |B| 2 ) is no longer conserved. This effectively leads to a timedependent interaction strength and, therefore, different stability behaviour during the time evolution. We shall now show that the stability of the nonlinear Bloch bands can have a strong influence on the mean-field dynamics.

Mean-field and two-particle dynamics
We now examine the dynamics of an initial pure BEC state (48). The parameters ψ j of the BEC state are taken to be the components of a nonlinear stationary state φ j (54), weighted by a Gaussian envelope ψ j ∝ φ j e −(j−x0) 2 /2σ 2 . In particular, we assume that k 0 = 0 and consider the two stationary solutions in (64). The renormalised density dynamics are shown in Figure 7 for both initial states, with γ = 0.1 and F = 0.2. For a weak interaction strength (g = 2) and a state initialised in the lower band, the particles perform Bloch oscillations accompanied by a splitting of the beam in position space. As expected, the dynamics are similar to the single-particle system. However, if the state is initialised in the upper band then the dynamics deviate from the single-particle case. This can be attributed to the change in the stability of the upper band, which is depicted in the top left panel of Figure 6. For a strong interaction strength (g = 7) the lower band is initially stable, while the upper band is initially unstable (Figure 6). This has a dramatic effect on the meanfield dynamics, as can clearly be seen in Figure 7. In particular, the strong instability of the upper band results in a focussing of the density towards a single lattice site. The subsequent dynamics are not dissimilar to the frequency doubled single-particle Bloch oscillations depicted in Figure 2.
Surprisingly, even far from the mean-field limit, the stability of the density dynamics has a dependence on the initial state. This is illustrated in Figure 8 for a twoparticle system with γ = 0.1 and F = 0.2. For a weak interaction strength (U = 1) and a state initialised in the lower band, the density behaves like a single particle. However, if the state is initialised in the upper band, then a squeezing and broadening of the density beam is observed, as well as additional structure from t ≈ 0.5T onwards. Increasing the interaction strength (U = 7) leads to additional many-body features, yet there is still a clear difference between the density dynamics for different initial states. A thorough investigation of the two-particle quantum dynamics is left for future work.

Summary
We have investigated Bloch oscillations in a Bose-Hubbard chain with single-particle losses from the odd lattice sites. We demonstrated that the single-particle system may be viewed as a quantum implementation of a PT-symmetric tight-binding Hamiltonian. The spectrum of the single-particle effective Hamiltonian was analysed and found to consist of two ladders of eigenvalues, which we linked to dynamical effects. In the mean-field limit of a many-particle system we derived analytic expressions for the generalised nonlinear stationary states and the nonlinear Bloch bands. The combination of nonlinearity and particle losses led to unusual features in the nonlinear Bloch bands, such as the vanishing of solutions and the formation of additional exceptional points. The stability of the stationary states was determined via the Bogoliubov-de Gennes equation, and was shown to strongly influence the mean-field dynamics. Finally, for a two-particle system far from the mean-field limit, the stability of the density dynamics was shown to have a dependence on the initial state.