Breakdown of quantization in nonlinear Thouless pumping

The dynamics of solitons driven in a nonlinear Thouless pump and its connection with the system's topology were recently explored for both weak and strong nonlinear strength. This work uncovers the fate of nonlinear Thouless pumping in the regime of intermediate nonlinearity, thus establishing a fascinating crossover from the observation of nonzero and quantized pumping at weak nonlinearity to zero pumping at strong nonlinearity. We identify the presence of critical nonlinearity strength at which quantized pumping of solitons breaks down regardless of the protocol time scale. Such an obstruction to pumping quantization is attributed to the presence of loop structures of nonlinear topological bands. Our results not only unveil a missing piece of physics in nonlinear Thouless pumping, but also provide a means to detect loop structures of nonlinear systems investigated in real space.


Introduction
Thouless pumping is one of the most prominent examples of quantized topological transport [1,2,3].
The behaviour of the above-mentioned nonlinear systems is best described by the nonlinear Schrödinger equation known as the Gross-Pitaevskii (GP) equation, where the nonlinear terms come from the mean-field treatment of the interacting bosons [24,25]. This type of nonlinear problem has been extensively studied in recent years, and has been shown to exhibit several exotic features, two of which are especially relevant for the present work. First, it is known that the energy spectrum of some nonlinear systems (to be referred to as nonlinear energy spectrum hereafter) with respect to some parameter may develop self-crossing bands such as looped-bands or swallowtail structures [26,27,28,21,29,23,30,31,32]. Second, certain nonlinear systems are known to support solitons [33,34,35,36,37,38,39,40,41,42]. In this manuscript, we refer to solitons as localized wavepackets that arise at finite nonlinearity.
In exploring physics arising from the interplay of nonlinearity and topology, the Thouless pumping of such solitons has been experimentally observed [43], and theoretically studied in both the weak nonlinear regime [44,45,46] and the strong nonlinear regime [46,47]. In the fermionic case, the breakdown of quantized Thouless pumping in the presence of Hubbard interactions has been experimentally demonstrated [48]. Among more recent developments, nonlinearity is known to fundamentally modify the system's adiabatic following dynamics, leading to pathdependent dynamical phases that modify the Berry connection and, consequently, the Aharanov-Bohm phase [49] and the Zak phase [50]. This motivated us to reexamine the nonlinear Thouless pumping especially if its quantization sustains due to nonlinearity. Using a typical model system plus nonlinearity, here we demonstrate that a dynamically evolved state during a pumping protocol may differ from the instantaneous nonlinear stationary state no matter how slow the protocol is, even if the corresponding linear model remains gapped throughout the whole process. This puzzle is explained by finding the presence of self-crossing bands (to be further elaborated below) in the energy spectrum as a function of some adiabatic parameter. The emergence of selfcrossing bands in the energy spectrum at moderate nonlinearity marks the breakdown of adiabaticity, which then leads to non-quantized pumping. We hence discover a mechanism accounting for a definite breakdown of quantized soliton pumping, serving as a crossover from quantized nonzero pumping at weak nonlinearity [44,45,46] to zero pumping at strong nonlinearity [46].
To gain another perspective on the role of self-crossing bands in altering the pumping dynamics, we consider a different but analytically simpler type of pumping involving Bloch states, then explicitly derive the nonlinear correction to the displacement that in turn also leads to the breakdown of quantization. In this case, we find that the critical nonlinearity strength at which loop structures start to appear corresponds to an extremum for the nonlinear drift and a singularity in both the Berry curvature and its nonlinear correction.
To gain analytical insight into the role of loop structures in altering the pumping dynamics, we explicitly derive the nonlinear correction in the particle displacement for a Bloch state. Specifically, we find that the critical nonlinearity strength at which loop structures start to appear corresponds to an extremum for the nonlinear drift and a singularity in both the Berry curvature and its nonlinear correction. This paper is organized as follows. In section 2, we introduce our model of a nonlinear Thouless pump, and present the dynamical evolution of solitons for different regimes of nonlinearity strength. Major results include the failure of adiabatic following no matter how slow the protocol is, with the main peak of the soliton still being displaced in spite of some spreading all over the lattice. In section 3, we study the nonlinear energy spectrum of our model. We highlight the presence of self-crossing bands for a certain range of nonlinearity, and show how it correlates with three very distinguishable behaviors for the pumped soliton. In section 4 we study the pumping of a Bloch state in our model and derive an additional contribution to the average displacement of a particle due to nonlinear dynamics. In section 5 we highlight the crucial differences between the two types of pumping considered, and discuss how self-crossing bands developed in nonlinear energy bands break the quantization of both pumpings in different ways. In section 6 we summarize the main findings of this work and suggest possibilities for future studies.

Nonlinear Thouless pump
In this work, we consider a nonlinear Thouless pump, consisting of a chain of N dimers. It is made up of a time dependent Rice-Mele model to which we add focusing nonlinearity [51]. The resulting system is described by the following set of nonlinear Schrödinger equations for j = 0 . . . N − 1, where v(t) = −[J +δ sin ωt] is the intracell coupling, w(t) = −[J −δ sin ωt] is the intercell coupling, u(t) = −∆ cos ωt is the staggered on-site potential, g > 0 is the strength of the focusing Kerr-like nonlinearity, and ω is the modulation frequency, taken small enough to allow adiabatic evolution. Here, adiabaticity is defined in a similar manner as its linear counterpart. Namely, it refers to the condition under which an initially prepared stationary state evolves into another stationary state of a set of slowly varying nonlinear Schrödinger equations. We further refer to the latter as adiabatic following. The parameters of the system are taken to be periodically driven sufficiently slowly, so that in the linear limit g = 0 the model realizes adiabatic charge pumping. In the linear case, the lower band has a Chern number of 2, meaning that a wavepacket uniformly occupying the lower band will be displaced by 2 sites over one adiabatic pumping cycle. The corresponding nonlinear model is a type of a time-dependent GP equation, and although it was historically used to describe the dynamics of Bose-Einstein condensates of ultracold bosonic matter [24,25], it is equally useful to describe the propagation of pulsed light through an array of waveguides, the most popular experimental platform to realize nonlinear pumping of solitons [43]. In this case, the nonlinear term is caused by the optical Kerr effect. For all calculations below, we use J = 1, δ = 0.5 and ∆ = 1, with N = 100 unit cells (2N = 200 sites), and we impose open boundary conditions (OBC) with Ψ −1 = Ψ 2N = 0. All physical variables presented in this work are assumed to be scaled, and therefore are in dimensionless units.
In the following, we will compare the dynamics of a soliton throughout the pumping cycle obtained from two different methods. In the first method, the soliton at time t is obtained by finding the stationary state of equation (1) via an iterative procedure detailed in Appendix A, using the instantaneous soliton at the previous time step as the trial state. In the second method, the soliton at time t is obtained by directly solving equation (1) using 4th order Runge Kutta method under a given initial soliton. Note that particular care must be taken in the implementation of the latter as the nonlinear term that can strongly vary even for small time steps. To this end, we updated the rate of change in the method at each step based of the current state of the system and additionally checked the convergence of the method by comparing solutions at different time steps. The first method corresponds to a perfect adiabatic following, meaning that at any given time the soliton obtained is an instantaneous stationary state of the nonlinear Hamiltonian. On the other hand, the second method corresponds to the actual dynamical evolution of the soliton and depends on the modulation frequency ω. In experiments, the soliton dynamics is then expected to follow the description of the second method. Under the same initial soliton and assumption of adiabaticity, it is naïvely expected that the soliton dynamics obtained from both methods coincide at all time. We will however see that this is not the case for some range of nonlinearity, hinting for a breakdown of the adiabatic process. In our numerical studies, we prepare our initial soliton to be localized around site n = 100, which is found using an iterative, self-consistent algorithm presented in Appendix A under a trial hyperbolic secant profile wavepacket.
In figure 1, we show the position expectation value ⟨X⟩ = j j |Ψ j | 2 of the soliton, hereafter to be equivalently referred to as the position of the center of mass (COM) of the soliton, over one adiabatic cycle following the two methods above. Figure 1a) shows that for weak nonlinear strength, both methods yield identical results at all time; the evolution of the soliton remains adiabatic, and the total displacement over one cycle is quantized to the Chern number of a band in the linear limit. This quantized displacement was explained in [44,45,46] in terms of the Wannier states expansion of the soliton. On the other hand, in the strongly nonlinear regime, shown in figure 1-c), the soliton merely oscillates around the same site and yields zero displacement over one cycle. This absence of pumping can be explained by the nonlinearity induced Rabi oscillations between two of the lowest bands, whose dynamical Chern numbers sum up to zero [46]. Another qualitative explanation for this zero displacement at high nonlinearity strength is that as when the nonlinearity strength g increases to the point that all other terms become negligible with respect to the constant onsit potential, the pumping path becomes effectively flat and the dynamically evolved solitons appear almost constant. Note that both methods once again yield identical results at all time, showing that the evolution of the soliton is also adiabatic in the strong nonlinearity regime.
It is important to note that even if the displacements of the soliton in the weak and strong nonlinearity regimes are quantized, the nonlinear pumping of solitons is a completely different phenomenon than linear Thouless pumping. While linear Thouless pumping assumes uniform occupation of the pumped band, this is definitely not the case when considering the pumping of a soliton, for which the localization of the wavepacket can oscillate between different sites [46].
As a main observation of this work, there exists an intermediate regime where pumping is still occurring, but no longer quantized, as shown in figure 2-b). Over one period, the trajectory of the COM of instantaneous solitons itself is discontinuous at two points. It is precisely at these points that the dynamically evolved soliton starts to deviate from the instantaneous soliton, and its position expectation value begins to oscillate irregularly. These discontinuities in the trajectory of the COM of instantaneous solitons are tied to the observation that, as nonlinearity increases, these trajectories tend to remain close to integer positions ⟨X⟩ = n for longer durations, swiftly moving from one near-integer position to the next one over a short period of time, as shown later in figure 6-a). This tendency culminates in the opening of gaps in the trajectory around the half-integer position ⟨X⟩ = n+ 1 2 as instantaneous stable solitons with such positions cease to exist.
In figures 2-c,d,e,f), we inspect the soliton's profile and its Inverse Participation Ratio (IPR) at different time steps in the intermediate nonlinearity regime. The IPR of a state |Ψ⟩ is defined by and is small (large) for localized (delocalized) states. Despite the irregular variations in the position expectation value, the dynamically evolved state actually remains strongly localized throughout the pumping cycle and ends with a strong peak at site n = 102 after one period. Figure 2-c) indeed shows that the IPR remains close to the minimum value of 1 during the entire cycle. The seemingly random variations in the position of the COM are in fact due to the apparition of noisy perturbation as the wavepacket evolves from a soliton, instantaneous stationary state of the nonlinear Hamiltonian, to a noisier profile (imperfect soliton) as it slightly disperses to the left side of the lattice. This spreading first occurs at the time of the first jump as shown in figure 2-e). This perturbation ends up spreading all over the lattice as shown in figure 2-f) and leads to the nonquantized pumping despite the wavepacket still exhibiting a strong peak at site n = 102 as in the weak nonlinearity regime. Therefore, despite the nonquantized value in the position of the COM, a quantized result is still obtained when solely viewing the dynamics of the wavepacket's peak. However, given that the position of the COM is directly linked to the Chern number in the linear limit, it remains a more important quantity that provides insights into the topology of the system. That this quantity is no longer quantized thus suggests a hidden mechanism that leads to a nonlinearity induced topological phase transition in the intermediate nonlinearity regime. We will uncover such a mechanism in the next section.
Before ending this section, it is worth investigating the fate of the soliton over more than one pumping cycle. For many pumping cycles, the position of the COM of the dynamically evolved state follows a rightward trend, although it fails to follow the instantaneous solitons due to the perturbations generated to the left at each discontinuous jump, as shown in figure 3-a). Figure 3-b) shows that the IPR remains low with regular peaks at each discontinuity where the spreading occurs. One may notice that in between jumps the IPR decreases, meaning the wavepacket tends to relocalize around its peak. This is due to the instantaneous solitons being dynamically stable. Hence, when moved slightly away from an instantaneous soliton, the wavepacket will tend to go back to a dynamically stable instantaneous soliton. Such a small perturbation away from a stable instantaneous stationary states happens at each jump, and tends to correct itself due to the dynamical stability, until another perturbation occurs at the next jump. This is why we can finally see that even after several pumping cycles, the final state shown in figure 3-c) still has a strong peak on the expected site, i.e., N = 120 when corresponding to 10 pumping cycles.

Breakdown of the adiabatic path due to the emergence of self-crossing bands
The breakdown of adiabaticity in the intermediate nonlinearity regime must be investigated. To that end we examine the nonlinear energy spectrum of the system, using an iterative approach described in Appendix A. In doing so we take as trial states at any given time t the eigenstates of the corresponding linear model, along with the nonlinear eigenstates found at the previous time step, in order to obtain a large Note that the energy bands corresponding to soliton states are highly degenerate, since solitons can live at any unit cell in the bulk.
In the weakly nonlinear regime of figure 4-a), the lowest band corresponds to the pumped solitons, while the third band is another soliton band that is not involved in the pumping process. All the other bands are delocalized bulk states (with IPR(|Ψ⟩) > 100). The quantized pumping by one unit cell is clear from the color-coding of the band and is made possible by the continuity in the energy band that can be followed adiabatically.
However, as nonlinearity strength increases, the energy spectrum of the system will undergo radical transformations, as shown in figure 5. The first well documented effect of the focusing nonlinear term on the system is the apparition of soliton states bifurcating from the linear bulk bands [50,44,47]. In the system under study, one soliton band bifurcates from each bulk band, corresponding to the two soliton bands shown in figure 4-a), the lowest one being the pumped band. As nonlinearity strength increases, the delocalized bulk bands then start to disappear due to the action of the focusing nonlinear term, that progressively favors strongly localized soliton states over delocalized states [50]. The next interesting observation that can be made is that as their energy reached the level of the original linear lower bulk band, the solitons bifurcating from the higher band disappear, only leaving the pumped solitons bifurcating from the lower band. However, for even higher nonlinear strength, another soliton band appears below the energy level of the lower band. This soliton band is different in nature from the one that disappeared earlier, as it emerges from the self-crossing band, whose apparition is explained later in figure 6.
Hence, the energy spectrum has already drastically changed in the intermediate nonlinearity regime of figure 2-a). Delocalized bulk states have already completely disappeared due to the focusing nonlinearity and the entire energy spectrum is constituted of soliton states. More strikingly, the soliton band now develops selfcrossings typically observed in the energy bands of nonlinear Bloch systems [26]. In this system, the self-crossing bands result in two incomplete V-shaped bands in the vicinity of t = 0.25T and t = 0.75T in figure 2-a). These self-crossing structures are attached to the original soliton band that spans the entire time domain and exists for all values of the nonlinearity strength, termed the main band. These incomplete Vshaped bands, along with looped-bands and swallowtail structures described in previous literature [26,27,28,29,31] collectively form a unique feature of nonlinear systems with no linear counterpart.
To our knowledge, the existence of such self-crossing bands depicting the energetics of real-space localized states, let alone their effect on pumping, has not been reported before. In the strongly nonlinear regime of figure 4-b), the two self-crossing structures merge, which causes the energy spectrum to consist of two gapless bands, each corresponding to solitons at even-and odd-numbered sites respectively. In particular, the crossing points between the two bands is responsible for the observation of Rabilike slight oscillations in [46] and figure 1-b). Note that although the existence of band crossing is known to threaten adiabaticity, the solitons belonging to the red and blue bands are strongly localized at different sublattice and hence cannot hybridize, preserving the adiabatic following throughout the cycle in figure 1-c).
We will now highlight the correlation between the emergence of self-crossing bands and the breakdown of the adiabatic following. For this purpose, we study in figure 6 the transition between weak and intermediate nonlinearity regime, by considering the energy spectra for values of g just lower and higher than the critical value g crit ≈ 3.2 at which self-crossing first appears, along with the corresponding dynamical evolution of the pumped solitons. As the system is approaching the critical value g crit in figure 6a) and figure 6-b), the lowest energy band starts to develop cusps around which the solitons evolve at a faster rate while remaining perfectly adiabatic. However, as soon as the nonlinearity strength exceeds the critical value as in figure 6-c) and figure 6d), the cusps evolve into self-crossing structures, around which the trajectory of the instantaneous solitons becomes discontinuous and adiabaticity is broken.
The breakdown of adiabaticity in the presence of self-crossing bands can also be intuitively understood from figure 2. Figure 2-b) shows that a soliton initially (at t = 0) peaked on site n = 100 evolves by adiabatically following the path formed by the blue part of the energy band in figure 2-a) up to t ⪆ 0.25T , slightly moving towards site n = 101. Ultimately, the blue band reaches a "dead-end", beyond which there are no dynamically stable soliton solutions to adiabatically follow into. The soliton is then forced to perform a non-adiabatic jump to the green part of the energy band, which corresponds to the dynamically stable soliton with the most overlap. This in turn leads to the jump in the trajectory of instantaneous solitons, which cannot be followed by the dynamically evolved soliton. In the latter case, most of intensity of the dynamically evolved soliton will make the jump to the next suitable stable nonlinear eigenstate, while the remaining intensity will spread linearly through the lattice as shown in figure 3. A similar dead-end is observed again at t ⪆ 0.75T , which is also accompanied by another discontinuity in the trajectory of instantaneous solitons.
On the other hand, the transition from intermediate to strong nonlinear regime can also be explained by the closing of the self-crossing bands. As nonlinear strength increases, the self-crossing structures will expand in time domain. As long as there remain incomplete bands, the "dead-end" phenomenon leading to non-adiabatic jumps will take place leading to the non-quantized, perturbed pumping shown in figures 2-a,b). However, beyond a second critical value for g, the incomplete bands will finally spread over the entire time domain, forming the second complete band seen in figure 4-b), marking the beginning of the strong nonlinear regime characterized by zero pumping and a recovered adiabatic following throughout the entire pumping cycle, as shown in figure 1-b).
The effect of self-crossing bands on soliton pumping can be clearly understood by considering all possible trajectories of instantaneous solitons over one pumping cycle as shown in figure 7. Figures 7-a,c) respectively correspond to the weak and strong nonlinear regime, where quantized pumping is possible and was experimentally observed [43] thanks to the contiguous path that can be adiabatically followed, leading respectively two sites away or to the same site after one pumping cycle. In the intermediate nonlinear regime shown in figure 7-b), the apparition of self-crossing bands coincide with the disappearance on any contiguous path making adiabatic following impossible. Although this intermediate regime only occurs over a very small window for the system considered in reference [43], it is an important transition regime showing the failure of adiabaticity due to nonlinear effects.
The analysis of the intermediate nonlinearity regime, characterized by the presence of self-crossing bands in the real space energy spectrum and the breakdown of the adiabatic following due to the dead-ends in the adiabatic path, plays the role of missing link between the two other, more extensively studied, regimes. It allows us to draw the bridge between the weakly nonlinear regime where adiabatic following is possible and pumping of the soliton remains quantized [44,45,46], and the strongly nonlinear regime, where solitons oscillate around the same state, amounting to a zero total displacement [46]. Moreover, the correspondence between the presence of self-crossing bands and the breakdown of quantized pumping established here provides a means of detecting the former via an adiabatic pumping experiment.

Nonlinear adiabatic pumping of Bloch states
To gain a different perspective into the breakdown of adiabaticity and nonquantized pumping due to self-crossing bands, we will now consider a different type of pumping involving Bloch states rather than a localized soliton. To this end, we consider the model under Periodic Boundary Conditions (PBC), and assume Bloch state solutions At any given time, we can solve the system for instantaneous stationary states by solving the nonlinear eigenvalue problem H(t, k, |Φ(t, k)⟩) |Φ(t, k)⟩ = E |Φ(t, k)⟩, where |Φ⟩ = [Φ e , Φ o ] T is a pseudo-spinor and H(t, k, |Φ(t, k)⟩) is the instantaneous two-band Gross-Pitaevskii (GP) Hamiltonian By defining Σ = |Φ o | 2 − |Φ e | 2 as the population difference between the two pseudospinor components, we can rewrite the Hamiltonian in a more compact form where h(t, Σ) = h z (t) + g 2 Σ and I 2 is a 2 × 2 identity matrix that will be ignored in later calculations as it merely contributes a shift of the energy.
In the corresponding linear model, the system is gapped at all times for all quasimomenta, allowing for adiabatic pumping of particles. As shown in figure 8a), a similar gapped band structure is observed at weak nonlinearity. However, at intermediate nonlinearity, i.e., g > g c = 2, the nonlinear energy spectrum starts to sport self-crossing bands, which this time form fully looped structures made of a pair of incomplete energy bands that do not span the whole two-dimensional Brillouin Zone, see e.g. the orange-and green-colored bands in figure 8-b). We will now analytically show that the development of loop structures plays a significant role in modifying the pumping of such a Bloch state, this time without breaking the adiabaticity of the process.
Recall that the displacement after one pumping cycle can be written as where ⟨...⟩ represents the average over the state at a given (k, t). During adiabatic evolution, this state is known to only slightly deviate from the instantaneous stationary state by a quantity whose average yields the Berry curvature. However, we show in Appendix B that the presence of nonlinear terms causes the state to deviate further from the stationary state, resulting in an additional contribution which averages in a nonlinear correction to the Berry curvature. For the system given in equation (7) and a state Φ = (cos [θ(k, t)/2], sin [θ(k, t)/2] exp[iϕ(k, t)]) T , the average displacement over one adiabatic cycle is given by where B(k, t) is the Berry curvature and is the nonlinear correction to the Berry curvature causing the average displacement to drift away from the integral of the conventional Berry curvature. This shows that unlike linear Thouless pumping, the displacement of the nonlinear adiabatic pumping is not simply proportional to the Chern number due to this extra second term. It is however much closer conceptually to the linear Thouless pumping than the soliton pumping, as we also assume an averaged displacement over all possible Bloch waves, so that it reduces to the standard linear pumping in the linear limit g = 0.
We then compute the average displacement over one adiabatic cycle for the lowest band of our system, for different values of the nonlinear strength. Results are shown on figure 9. Note that although the result presented above is derived for the highest band of the two band system by choosing |+⟩ = (cos (θ/2), sin (θ/2)e iϕ ) T as our initial state, it can easily be adapted to the lowest band with initial state |−⟩ = (sin (θ/2), − cos (θ/2)e iϕ ) T by simply replacing θ by θ ′ = θ − π in all expressions.
The variations of the displacement with the nonlinearity strength shown in figure 9a) present many interesting features. First for weak nonlinearity strength, the pumping occurs but is no longer quantized to the Chern number due to the additional nonlinear drift. Moreover, in the strongly nonlinear limit, the average displacement goes to zero, with or without the additional nonlinear drift. Interestingly, the transition between a quantized displacement and a zero displacement starts exactly at g = g c = 2, when the loop structures first appear. Beyond this critical value, the displacement is no longer quantized even without taking the nonlinear drift into account, and starts continuously decreasing to zero. Intuitively, this is attributed to the fact that in the presence of loop structures, the nonlinear energy bands are no longer gapped and their Chern number is consequently no longer quantized.
Considering the additional nonlinear drift alone also offers some interesting observations as shown in figure 9-b). Starting small at weak nonlinearity, it reaches a peak value, then becomes small again at strong nonlinearity. Specifically, the peak value occurs exactly at the critical value g c where the loop structures start to appear. Coincidentally, the critical value g c corresponds to the first apparition of singular points in the expression of both the Berry curvature, and its nonlinear correction. These results are shown in detail in Appendix C.

Discussion
While the analysis presented in section 4 highlights the direct correlation between selfcrossing bands and nonquantized pumping, it is to be stressed that the results therein may not have a direct correlation with our numerical results presented in section 2 and 3. Indeed, as elucidated before, the former represents pumping of Bloch states, whereas the latter describes pumping of localized solitons.
The two pumping scenarios above only reconcile in the absence of nonlinearity due to the momentum integral in equation (8), which effectively turns the various Bloch states into a localized Wannier state. Since superposition principle no longer holds in the nonlinear setting, equation (8) is in general inequivalent to the displacement of the particle over one pumping cycle at finite nonlinearity. This in turn explains the order of magnitude difference between the nonlinear correction analytically derived in section 4 and that numerically observed in section 2 and 3. Nevertheless, it is remarkable to highlight that the pattern of increasing discrepancy between dynamically evolved solitons and instantaneous solitons as the self-crossign structures grow holds in both types of pumping. This demonstrates that our finding that self-crossing bands, which are a purely nonlinear phenomenon, can lead to the breakdown of quantized pumping in spite of the topological stability of the system holds for a general type of pumping.
It is also worth noting the difference in normalization prescription for the two pumping treatments. That is, in the soliton pumping of sections 2 and 3, the wave function is normalized over the entire lattice, i.e., 2N n=1 |Ψ n | 2 = 1. By contrast, in the Bloch states pumping of section 4, the wave function is only normalized over a unit cell, i.e., |Φ e | 2 + |Φ o | 2 = 1. This leads to different scaling of the nonlinearity strength g which, among others, contributes to the discrepancy between the critical nonlinearity values in the two pumping treatments. Establishing a quantitative comparison between the two pumping treatments thus remains a challenging question and will be left for future work.

Concluding remarks
In this work, we have studied the correlation between the pumping of solitons and the nonlinear energy bands of a system for a wide range of nonlinearity strength. We established the missing link between earlier studies of [44,45,46,47] by identifying an intermediate nonlinearity regime in which the pumping is nonzero but no longer quantized. This phenomenon is attributed to the presence of self-crossing bands in the nonlinear energy band, disrupting the trajectory of the evolved solitons. Finally, we analytically quantify the nonlinear correction to displacement in the case of Bloch state pumping, which further quantifies the effect of such self-crossing bands. In another recent work, we explicitly utilize this nonlinear correction to probe nonlinear effects on Weyl semimetals [52].
In a future study, it would be interesting to modify the adiabatic driving, e.g., in the spirit of the driven nonlinear Landau-Zener process presented in [53,30], to bypass the self-crossings. This would enable a successful (quantized) adiabatic pumping for a larger range of nonlinearity, extending it to the intermediate nonlinearity regime, and possibly even in the strongly nonlinear regime. Appendix A. Self-consistent iterative method for instantaneous stationary states In order to find instantaneous stationary states of the nonlinear lattice at any given time, we use an iterative method. For a given nonlinear, state dependent Hamiltonian H, at a given time t, the iterative process from a state |Ψ n (t)⟩ to the state |Ψ n+1 (t)⟩ is as follows: • We first compute H n = H(|Ψ n ⟩ , t), the nonlinear state-dependent Hamiltonian evaluated at the state |Ψ n (t)⟩.
we repeat this iterative process until the overlap between old and new state is close enough to 1 by an arbitrary ϵ, i.e. |⟨Ψ n+1 (t)|Ψ n (t)⟩| > 1 − ϵ. Throughout this work, we take ϵ = 10 −10 . In order to execute this iteration method, one also needs to choose the initial state used as a starting point of the iteration. Note that the stationary state obtained is quite dependent on the chosen starting point (trial state). In order to obtain instantaneous solitons, we found that using sech states as trial states gives the best results [54,55,56,57]. To obtain the nonlinear energy spectra, we also used eigenstates of the corresponding linear model as trial states, allowing us to obtain delocalized bulk states in the weakly nonlinear regime.

Appendix B. Displacement correction in Bloch state pumping
We consider a general two level GP Hamiltonian We assume that H, h x , h y , h are all functions of a quasimomentum k and a periodic time modulation ωt where ω is the frequency of said modulation, always appearing with t. We define a state Ψ a = e −if Φ a with a = 1, 2, which corresponds to an element of a projective Hilbert space [58]. The total phase f is taken to capture both dynamical and geometric phases of the state |Φ⟩. With this notation, the nonlinear Schrödinger equation now reads (summation of repeated indices being implied) Applying a Ψ * a ... to both sides, we obtain Doing a perturbative expansion of f and Ψ a under an adiabatic parameter ϵ gives df dt = α 0 + α 1 ϵ + ...
where we notice that the first term in the right hand side of the bottom line corresponds to the conventional Berry connection, and the second term is the geometric contribution coming from the dynamical phase, due to nonlinear dynamics. In our case, we have The general formula for α 0 and α 1 given in equation(B.6) becomes then and for a = 2 For a two-level system, the stationary state |+⟩ can be written without loss of generality in the form so that we can simplify equation (B.9) by taking its real part, and making use again of the normalization condition cos θ 2 Re(Ψ 2 ) = 0 to get (B.14) Using the normalization condition, we also get Re(e −iϕ Ψ 2 ) as In order to determine their imaginary part, we now take the imaginary part of equation (B.9) and equation (B.10) to get, respectively, which can be more compactly written as Now using the eigenstate Ψ (0) = |+⟩ with energy E and the hidden eigenstate Ψ (0)⊥ = + ⊥ with energy E ⊥ = −E, we can rewrite H (0) as and using the fact that Ψ Since Ψ (1) is defined as a small correction to Ψ (0) we have a freedom to make it orthogonal to Ψ (0) (components of Ψ (1) parallel to Ψ (0) can be separated and ‡ + ⊥ is however not a stationary state of the system, as H (0) is state dependent, and + ⊥ is an eigenstate of H (0) (|+⟩), but not necessarily of H (0) + ⊥ combined with Ψ (0) , which only affects the global phase of Ψ (0) ). In such a case, the projection of Ψ (1) onto Ψ (0)⊥ will not change anything, and we finally find On the other hand, equation (B.14) and equation (B.15) give which allows us to write In the case of the Hamiltonian presented in equation (7) in the main text, where h(Σ) = h z + g 2 Σ, we have . (B.27) The first term on the right hand side is responsible for the Berry curvature appearing in the expression of the average displacement, while the second term on the right hand side is due to nonlinear effects, and will lead to an additional displacement. We then move on to the evaluation of the displacement of a particle, whose expression is given by equation (8) in the main text, while taking the average over the state |Ψ⟩ = Ψ (0) + Ψ (1) with Ψ (0) = |+⟩ with eigenenergy E. In equation (B.27), we rewrite the second term on the right hand side as The first two terms are also present in the linear case. In particular, we can identify the second term as the Berry curvature B(k, t), the integral of which yields a quantized Chern number. We can further simplify the third to sixth lines by noticing that to finally obtain the expression in equation (9) of the main text.
Using Σ = − cos θ and ∂Σ = sin θ ∂θ (where ∂ can be any partial derivative, ∂ t or ∂ k ), we can rewrite D(k, t) in equation (10) of the main text as The Berry curvature can also be simplified in a similar fashion. First, we express the Berry connection as where α can be t or k. Then, using B(k, t) = ∂ t A k − ∂ k A t we get At points (k, t) = (π, π 2 ) and (k, t) = (π, 3π 2 ) using the self-consistency equation given in equation (C.4), we also get . (C.8) By considering the solution Σ = g 2 −4(h 2 x +h 2 y ) g 2 in equations (C.5), (C.7) and (C.8), we get where ∂ t ϕ and ∂ k ϕ do not depend on g. We can then see that both quantities exhibit a singularity point at the critical value g c = 4δ.