On the Quantization of AB Phase in Nonlinear Systems

Self-intersecting energy band structures in momentum space can be induced by nonlinearity at the mean-field level, with the so-called nonlinear Dirac cones as one intriguing consequence. Using the Qi-Wu-Zhang model plus power law nonlinearity, we systematically study in this paper the Aharonov–Bohm (AB) phase associated with an adiabatic process in the momentum space, with two adiabatic paths circling around one nonlinear Dirac cone. Interestingly, for and only for Kerr nonlinearity, the AB phase experiences a jump of π at the critical nonlinearity at which the Dirac cone appears and disappears (thus yielding π-quantization of the AB phase so long as the nonlinear Dirac cone exists), whereas for all other powers of nonlinearity, the AB phase always changes continuously with the nonlinear strength. Our results may be useful for experimental measurement of power-law nonlinearity and shall motivate further fundamental interest in aspects of geometric phase and adiabatic following in nonlinear systems.


Introduction
The dynamics depicted by a nonlinear discretized Schödinger equation (NDSE) can be extremely rich, including the emergence of many-dimensional chaos, solitons, and breathers, etc. The problem can be much reduced by assuming the translational invariance of a wave under consideration. With this assumption, the main physics is about the features of Bloch waves, the associated energy bands, and how they respond to changes in the parameters of a nonlinear system. Interestingly, the nonlinear Bloch bands of NDSE can induce gapless band structures absent in linear systems, such as 2-dimensional (2D) nonlinear Dirac cones [1] induced by Kerr nonlinearity [2]. Even more peculiarly, such nonlinear Dirac cones are formed by exotic nonlinear energy bands in a subregime of the Brillouin zone [1,[3][4][5][6][7][8].
As a close analog to a setting in real space to measure the Aharonov-Bohm (AB) phase around a singularity point with magnetic flux, let us now imagine two adiabatic paths, in the momentum space, circling around a band-crossing point. If we adiabatically change the Bloch momentum, so as to guide the Bloch wave to evolve along the two adiabatic paths, the final phase difference thus generated between the two adiabatic paths is termed the nonlinear AB phase [1]. One may naïvely think of the following: provided that the dynamical phases between the two adiabatic paths are identical and hence have zero contribution to the phase difference of interest, the obtained AB phase would be just the Berry phase associated with the band degeneracy point. The actual physics turns out to be more interesting than just a Berry phase. As a result of nonlinearity, any small deviation of the adiabatically following state from the instantaneous Bloch wave causes a tiny correction to the dynamical phase, and accumulation of such tiny corrections over the entire adiabatic protocol yields an unfamiliar geometrical phase on top of the expected Berry phase. Remarkably, as a possible means of topological characterization of nonlinear

Hamiltonian and Energy Spectrum
The momentum-space Hamiltonian is composed of a QWZ model with power law nonlinearity characterized by a parameter p: where σ i are Pauli matrices and ψ a are two components of the wavefunction, ψ = ψ 1 ψ 2 .
The normalization of the wavefunction means that |ψ 1 | 2 + |ψ 2 | 2 = 1. The nonlinearity parameter p is a non-negative real number. The Kerr nonlinearity corresponds to p = 1. The parameters k 1 and k 2 are two quasimomenta, whose values will be adiabatically tuned in order to implement an actual adiabatic protocol to generate the nonlinear AB phase.
To solve the nonlinear eigenvalue problem, we introduce a real parameter x as We will see later that the angular variable ϕ is the same as in the Figure 1. It turns out that x is the central quantity for expressing energy, dynamical phase, Berry phase, and nonlinear AB phase. It can be shown that the instantaneous eigenenergy is where x satisfies the following algebraic equation, with γ := J 1 sin k 1 − iJ 2 sin k 2 . In order to have a Dirac point in the energy spectrum, the energy must be doubly degenerate at k 1 = k 2 = 0. Since γ = 0 at this point, x must satisfy For simplicity, we choose β(k 1 , k 2 ) = B(−1 + cos k 1 + cos k 2 ).
Hence, β(0, 0) = B. It is clear that the nonlinearity strength g and energy E can be scaled in terms of B. Energy spectra with p = 1, 1.5, 2 and g = 2.5B are shown in Figure 2, where the Dirac cone is clearly visible around the origin. A perturbative analysis of energy spectrum near the Dirac cone can be found in Appendix A.1.

Dynamics of Adiabatic Following
To obtain the nonlinear AB phase, let us consider two adiabatic paths along a small circle around the origin k 1 = k 2 = 0. As shown in Figure 1, starting at the same point S, along each path, the system is guided to move along one half of the perimeter of the circle using the same amount of time. The two adiabatic paths are "recombined" at the end of the evolution at point N. As introduced in Section 1, the phase difference acquired by the system between two adiabatic paths is called the nonlinear AB phase. Clearly, the nonlinear AB phase here is the sum of the dynamical phase difference and the Berry phase associated with the closed loop around the band-degeneracy point. We shall study below the possible AB phase quantization for a varying nonlinearity strength g and for different nonlinear parameters p. The quasimomenta k 1 and k 2 associated with two spatial dimensions are parameterized by ϕ and will be made to adiabatically change.
At the starting point S, the system is assumed to be prepared in the Bloch eigenstate at momentum space location S. As the system adiabatically evolves along the path SEN or SWN, the time-evolving state deviates from the instantaneous eigenstate along the path, with the tiny deviation at the order of the adiabatic parameter ε. The slower the rate of adiabatic change is, the smaller ε is, and the smaller the deviation. Here, nonlinearity plays a key role. That is, the dynamical phase also obtains a correction at the order of ε. Since the total evolution time is of order O(ε −1 ), the O(ε) term in this phase correction will contribute an ε-independent term through accumulation, yielding a geometric phase term out of the dynamical phase. This will not occur in linear terms because such correction accumulated over the entire adiabatic process is at most of the order of ε, which vanishes for sufficiently slow adiabatic protocols.
The dynamics of the states is governed by the time-dependent Schrödinger equation, where the Hamiltonian is given by Equation (1) with ψ being replaced by Ψ. Here, the overhead dot denotes the time derivative. We will solve this equation up to the order of ε as described above. Through the lengthy computation, as illustrated in Appendix A.2, we obtain the instantaneous change rate of the overall phase of a time-evolving state aṡ with We recognize that the circular integration of the second term in Equation (10) is nothing but the Berry phase θ B , because it assumes the same form as in the linear limit. The rest of the phase is from the dynamical phase θ D , which contains two parts: the first part comes from the instantaneous eigenenergy E and the second part from the third term in Equation (10) as a new contribution from the nonlinearity. Specifically, In the event that the Dirac cone does exist at the point k 1 = k 2 = 0, the obtained phase difference between the two adiabatic paths described in Figure 1 then becomes the nonlinear AB phase θ AB . Since the two adiabatic paths are symmetric by construction and they take the same amount of time, the leading term in Equation (13) contributes the same in each of the two paths. Thus, the difference of the dynamical phases between two paths comes from the second term of Equation (13) only. Thus, the total nonlinear AB phase is Note that we take into account that the paths are chosen to be close to the Dirac cone (so that the cones indeed have linear dispersion relations), namely, |k 1 | and |k 2 | are small at all times. The leading behavior of the dynamical phase difference term is then found to be where ∆ 0 is ∆ evaluated at x = x 0 and k 1 = k 2 = 0, x 0 is the solution of Equation (6), and For the Berry phase, the leading behavior is As detailed in Appendix A.1, For |g| > 2B, a nonlinear Dirac cone is located at the origin. For |g| < 2B, the only possible solutions to Equation (5) are x = ±1 and there is no Dirac cone. For g ∈ (0, 2B), we can hence assign x 0 = −1, and for g ∈ (−2B, 0), we may assign x 0 = 1. With this convention, it is clear to see that θ B is constantly 0 (mod 2π) for g ∈ (−2B, 2B). The Berry phase θ B becomes nonzero and changes continuously for |g| > 2B. For each p, as we continuously tune g, x 0 can be easily solved numerically using Equation (6), thus obtaining the theoretical values of the leading terms of the dynamical phase, Berry phase and AB phase around the origin. We also numerically solve the evolution using the Schrödinger equation Equation (9) along the two paths, and compute the dynamical phase, AB phase and Berry phase using numerical solutions of the evolution. The evolution is computed using an operator-splitting algorithm. The results are presented in Figure 3.  > 1, the jump of dynamical phase at g = ±2B is ±2π; for p = 1, the jump of dynamical phase at g = ±2B is ±π; and for p < 1, the dynamical phase changes continuously. (b) The Berry phases change continuously for any p. (c) Only for Kerr nonlinearity p = 1, the AB phase has a quantized jump of π at the critical value g = ±2B and stays at π for |g| > 2B.
In each plot, solid lines are theoretical values, while dots on the solid lines are computed from numerical evolutions. In Figure 3a, for any p = 0.5, 1, 1.5, 2, 2.5, 3, the dynamical phase around the origin is 0 for g ∈ (−2B, 2B). At the critical value g = ±2B where the Dirac cone appears, for p = 0.5, the Dirac cone changes continuously with respect to g. For p = 1, there is a quantized jump of ±π at g = ±2B. For p = 1.5, 2, 2.5, 3, there is a quantized jump of ±2π at the critical value g = ±2B (so this is equivalent to no change). In Figure 3b, the Berry phase (modulo 2π) is identically 0 for g ∈ (−2B, 2B), and changes continuously with respect to g. In Figure 3c, the AB phase (modulo 2π) is the sum of the dynamical phase in Figure 3a and the Berry phase in Figure 3b. Only for p = 1, the AB phase has a quantized jump of π at the critical value g = ±2B and stays at π for |g| > 2B, as discovered by Ref. [1]. For all other values of p, the AB phase changes continuously with respect to g. The special behavior of p = 1 is a result of the fact that p = 1 is a critical value for the limit lim g→±2B ± δθ D , as will be explained in the next section.

Mechanism of the Jump of AB Phase at g = ±2B for Kerr Nonlinearity
For p > 1, we can factor out a factor (1 − x 2 0 ) from ∆ 0 which cancels the same factor in the numerator of δθ D , which equals ∓2π or equivalently zero since x 0 = ±1, for |g| = 2B or when the Dirac cone starts to appear. Likewise, for p = 1, we have which equals ∓π since x 0 = ±1, for |g| = 2B. Finally, for 0 < p < 1, which vanishes for x 0 = ±1, for |g| = 2B. The calculations above make it clear that the nonlinear AB phase associated with Kerr nonlinearity (p = 1) is most special as the extra nonlinearity-induced correction to dynamical phase experiences a π jump when the Dirac cone appears. What is intriguing for Kerr nonlinearity is that the nonlinear AB phase stays quantized at π for |g| > 2B, as θ B and δθ D happen to be complementary to each other, as shown in Equations (17) and (19). For all other forms of power-law nonlinearity, there is no such jump, π-quantization is thus absent, and consequently, the nonlinear AB phase only changes continuously with respect to g. This finally explains why in Figure 3 only the nonlinear AB phase for Kerr nonlinearity (p = 1) displays a quantization plateau for |g| > 2B.

Conclusions
In this paper, we analytically and computationally examined the so-called nonlinear AB phase around Dirac cones induced by power-law nonlinearity added to the QWZ model often used for studies of topological band structures. With our analytical results, we are able to explain why the nonlinear AB phase has a quantized jump of π when Dirac cone starts to appear or disappear, for and only for Kerr nonlinearity. In the context of nonlinear AB phase that can be in principle measured in experiments, Kerr nonlinearity is thus identified as a critical form of nonlinearity. As seen from our theoretical considerations above, our result will not be restricted to the QWZ model alone since it is based only on the asymptotic dispersion relation in the vicinity of the nonlinear Dirac cone. It is thus of considerable interest to investigate the generality of our results in other models with nonlinearity.
where χ is at least in the first order in k 1 and k 2 , and ρ and η are at least in the second order in k 1 and k 2 . Plugging Equations (A12)-(A14) into Equation (5), we have To this order, we obtain the correction to the parameter x, Plugging this into Equation (4), we find the expression for the eigenenergy, where the nonperturbed eigenenergy is We can see clearly from the expansion of E that there is a Dirac cone structure at the origin, provided |x 0 | < 1, which corresponds to |g| > 2B. For |g| < 2B, the system contains two smooth energy bands. At the critical value g = 2B (g = −2B), a kink will develop on the lower (upper) band at k 1 = k 2 = 0. Once g > 2B (g < −2B), a 2D self-intersection structure, i.e., a nonlinear Dirac cone, will appear from the lower (upper) band, whose vertex is at k 1 = k 2 = 0. This is true for any p > 0. We show five plots with different values of nonlinearity in Figure A1, along section k 1 = 0 and with −0.1π ≤ k 2 ≤ 0.1π. In each plot, the red dots are perturbative eigenenergies around the Dirac point (or at the origin for |g| ≤ 2B), while the blue lines are numerical solutions by solving Equation (5) exactly. We can see that the perturbative solutions perfectly match the numerical solutions for sufficiently small |k 2 |.