Nonlinearity-induced topological phase transition characterized by the nonlinear Chern number

As first demonstrated by the characterization of the quantum Hall effect by the Chern number, topology provides a guiding principle to realize robust properties of condensed matter systems immune to the existence of disorder. The bulk-boundary correspondence guarantees the emergence of gapless boundary modes in a topological system whose bulk exhibits nonzero topological invariants. Although some recent studies have suggested a possible extension of the notion of topology to nonlinear systems such as photonics and electrical circuits, the nonlinear counterpart of topological invariant has not yet been understood. Here, we propose the nonlinear extension of the Chern number based on the nonlinear eigenvalue problems in two-dimensional systems and reveal the bulk-boundary correspondence beyond the weakly nonlinear regime. Specifically, we find the nonlinearity-induced topological phase transitions, where the existence of topological edge modes depends on the amplitude of oscillatory modes. We propose and analyze a minimal model of a nonlinear Chern insulator whose exact bulk solutions are analytically obtained and indicate the amplitude dependence of the nonlinear Chern number, for which we confirm the nonlinear counterpart of the bulk-boundary correspondence in the continuum limit. Thus, our result reveals the existence of genuinely nonlinear topological phases that are adiabatically disconnected from the linear regime, showing the promise for expanding the scope of topological classification of matter towards the nonlinear regime.

Topology is utilized to realize robust properties of materials that are immune to disorders [1,2].A prototypical example of topological materials is the quantum Hall effect [3,4], which was discovered in a two-dimensional semiconductor under a magnetic field.In such a two-dimensional system, the Chern number characterizes the topology of the band structure and the corresponding gapless boundary modes.This bulk-boundary correspondence lies at the heart of the robustness of topological devices utilizing boundary modes.Recent studies have also explored topological phenomena in a variety of platforms, such as photonics [5], electrical circuits [6], ultracold atoms [7], fluids [8], and mechanical lattices [9].
While band topology has been well-explored in linear systems, nonlinear dynamics is ubiquitous in classical [10][11][12][13][14][15] and interacting bosonic systems [16,17].For example, nonlinear interactions can naturally emerge in the mean-field analysis of bosonic many-body systems, as is seen in the Gross-Pitaevskii equations of ultracold atoms.Recent studies have also analyzed the nonlinear effects on topological edge modes [11,12,[18][19][20][21][22][23][24][25][26][27] and revealed unique topological phenomena intertwined with solitons [28][29][30][31][32][33][34] and synchronization [35][36][37].Nonlinearity can further modify the conventional notion of topological phase; recent studies have revealed that onedimensional systems can exhibit nonlinearity-induced topological phase transitions, where the existence of topological edge modes depends on the amplitude of the oscillatory modes [38][39][40][41][42].While these previous studies have indicated the existence of topological edge modes in nonlinear systems, one cannot straightforwardly extend the topological invariants to nonlinear systems because they have no band structures.In addition, despite the advantage that nontrivial topology in twodimensional systems requires no additional symmetries other than the U (1) and translation symmetries as is the case for the quantum Hall effect, nonlinear topology in two-dimensional systems [18,29,30,34] is much less understood than that in one-dimensional systems.
In this paper, we introduce the notion of the nonlinear Chern number and reveal its relation to the bulk-boundary correspondence.To define the nonlinear Chern number of two-dimensional systems, we consider the nonlinear extension of the eigenvalue problem [40,42,43] and make an analogy to band structures.While it is not obvious that the nonlinear eigenvalue problem elucidates the bulk-boundary correspondence in nonlinear topological insulators, we theoretically prove the bulk-boundary correspondence of the nonlinear Chern number in weakly nonlinear systems.Furthermore, in stronger nonlinear regimes where nonlinear terms are larger compared to the linear band gap, we find that the nonlinearity-induced topological phase transition can occur in two-dimensional systems as depicted in Fig. 1.We construct a minimal model of a nonlinear Chern insulator, for which we can obtain exact bulk solutions and thus can analytically show the amplitude dependence of the associated nonlinear Chern number without any approximations.In the continuum limit of the model, we analytically show the nonlinear bulk-boundary correspondence, which states that the nonlinear Chern number predicts the existence of edge modes even under strong nonlinearity.In addition, we numerically confirm the bulk-boundary correspondence in the lattice model at the parameter where the lattice model well-approximates the behavior of the continuum model.We also discuss that the nonlinearity-induced edge modes can be observed in the quench dynamics in experimental setups such as nonlinear topological photonics.Since the nonlinearity-induced topological phases are disconnected amplitude

𝐻 𝜓
FIG. 1. Schematic of the nonlinear Chern insulators and the nonlinearity-induced topological phase transition.While topology of a noninteracting linear system can be characterized by the Chern number that is computed from its eigenvectors, topology of a nonlinear system is classified by the nonlinear Chern number, which utilizes the nonlinear extension of the eigenequation.In weakly nonlinear regions (i.e., small amplitude), the nonlinear Chern number predicts the existence of edge modes corresponding to those in linear systems.Specifically, when nonlinear systems exhibit edge-localized steady states, both the nonlinear and linear Chern numbers are nonzero as shown in the upper part of the figure.If we inject higher energy into the system and consider the eigenmodes with large amplitudes, the nonlinear band structure can become gapless.At such a gapless point, nonlinearity-induced topological phase transition can occur, where topological boundary modes appear with the nonzero nonlinear Chern number.The nonlinearity-induced topological phases exhibit boundary modes that cannot be predicted from the linear Chern number.Therefore, such topological phases are genuinely unique to nonlinear systems.
from the linear limits under adiabatic deformations, our results show the existence of genuinely nonlinear topological phases.
The scope of this paper applies to a broad class of twodimensional systems with U (1)-gauge and spatial translation symmetries, which can be realized in a variety of experimental setups.As is the case that the Thouless-Kohmoto-Nightingaleden Nijs formula [4] has triggered the research of a variety of topological materials, the nonlinear Chern number is expected to open up the research stream of nonlinear topological materials including their systematic classification.From the experimental point of view, one can realize nonlinear Chern insulators with the U (1)-gauge symmetry in, e.g., photonics [11,12,18,20,22,25,27,[29][30][31][32][33], ultracold atoms [16,17,27], and electrical circuits [6,35,36], where both linear band topology and nonlinear effects have been investigated.In particular, since Kerr nonlinearity [10] is fairly common in photonic systems, it should be possible to extend the current topological photonic devices to nonlinear ones.
Nonlinear eigenvalue problem and nonlinear Chern number.-Theconventional band topology is based on eigenvalue problems in linear systems.To define a topological invariant, one should consider eigenvalues E j (k) and eigenvectors |ψ j (k)⟩ of the Bloch Hamiltonian H(k) corresponding to the Fourier component of the real-space Hamiltonian of a periodic system.Then, one can define the Chern number as C j = (1/2πi) ∇ k × ⟨ψ j (k)|∇ k |ψ j (k)⟩d 2 k.The Chern number is unchanged under the continuous deformation of the Hamiltonian unless the Bloch Hamiltonian is degenerate.In this sense, the band topology is robust against perturbations, which is of practical use to realize topological edge modes in experimental setups.Thus, to define the topological invariant in nonlinear systems, it is essential to extend the notion of band structures.
We here consider the nonlinear extension of the eigenvalue equations [40,42,43] and define the nonlinear Chern number by using the nonlinear eigenvectors.We start from the general nonlinear dynamics, where Ψ j (r) is the state variable and f j (•; r) is a nonlinear function of the state vector Ψ.In lattice systems, we use the notation that r is a representative point in each unit cell of the lattice as shown in Fig. 2a.Then, j represents the internal degrees of freedom that include, e.g., sublattices and effective spin degrees of freedom.When we consider continuum systems, r should simply represent the location and j corresponds to the internal degrees of freedom such as spins.For example, the Gross-Pitaevskii equation in the continuous space is given by f (Ψ; r) = −∇ 2 Ψ(r)/(2m)+V Ψ(r)+(4πa/m)|Ψ(r)| 2 Ψ(r) with V being a potential and m and a being constants, and the nonlinear function f (•; r) depends on the Ψ(r) and its derivative and has no internal degrees of freedom.Since the quantum Hall system has the U (1) and translation symmetries, we impose them on the nonlinear equation to study the analogy of such a prototypical topological insulator.Concretely, the U (1) symmetry is represented as f j (e iθ Ψ; r) = e iθ f j (Ψ; r), which is satisfied in, e.g., the Kerr-like nonlinearity κ|Ψ j (r)| 2 Ψ j (r) [10].The translational symmetry in lattice systems is defined as f j (Ψ; r + a) = f j (Φ; r) with a being a lattice vector and Φ being translated state variables Φ j (r) = Ψ j (r + a).The translational symmetry in continuum systems is also defined in the same equations, while f j (Φ; r) still remains r dependence due to, e.g., the periodic potential.We also focus on conservative dynamics analogous to Hermitian Hamiltonians where the sum of squared amplitudes j,r |Ψ j (r)| 2 is preserved.
Corresponding to the nonlinear dynamical system in Eq. ( 1), the nonlinear eigenvector and eigenvalue are defined as the state vector Ψ with components Ψ j (r) and the constant E that satisfy We term the equation as a nonlinear eigenequation and analyze its bulk-boundary correspondence below.We note that we can regard the nonlinear eigenvector as a periodically oscillating steady state Ψ j (r; t) = e −iEt Ψ j (r) of the nonlinear system when the eigenvalue is real.Since the sum of the squared amplitude is conserved under the U (1) symmetry, we here focus on special solutions with the fixed sums of squared amplitudes j,r |Ψ j (r)| 2 (resp.j d n r|Ψ j (r)| 2 ) in lattice (resp.continuum) systems, where we take the sum both on the location and the internal degrees of freedom.
To extend the Chern number to nonlinear systems, we introduce the eigenvalue problems in the wavenumber space, which is analogous to the linear eigenequation of the Bloch Hamiltonian.In a lattice system with the translation symmetries, we assume an ansatz state [40] that we name the Bloch ansatz: Ψ j (r) = e ik•r ψ j (k).We note that in linear systems the Bloch theorem guarantees that every eigenvector is given by the form of the Bloch ansatz.On the other hand, in nonlinear systems, there can be nonlinear eigenvectors out of the description of the Bloch ansatz, including bulk-localized ones.Despite the existence of such localized modes, we here only focus on nonlinear bulk eigenvectors described by the Bloch ansatz and show that even such periodic bulk solutions can exhibit unique topological phenomena to nonlinear systems, i.e., the nonlinearity-induced topological phase transition.Under this ansatz and the U (1) symmetry, one can rewrite the nonlinear eigenequation as f j (k, ψ(k)) = E(k)ψ j (k) parametrized by k (see Supplementary Information for the detailed derivation).
To capture the nonlinearity-induced topological phase transition depending on amplitudes, we focus on a special so-lution of the nonlinear eigenvector at each k whose sum of the squared amplitudes j |ψ j (k)| 2 = w is fixed independently of the wavenumber k.We note that the assumption of such fixed-amplitude Bloch-ansatz solutions is consistent with the perturbation calculation of the nonlinear eigenvectors (see Supplementary Information for the detail).By using fixed-amplitude nonlinear eigenvectors, we define the nonlinear Chern number as (3) We note that this definition reduces to the conventional linear Chern number if f defined in Eq. ( 2) is a linear function.It is also noteworthy that since special solutions of nonlinear eigenvectors should exist at arbitrary w in ordinary nonlinear systems, we can define the nonlinear Chern number at any positive w except for gap-closing points.One can prove that the nonlinear Chern number is an integer by embedding the nonlinear eigenvectors into an eigenspace of a linear Bloch Hamiltonian (see also Supplementary Information).Since the eigenvector can be changed by the amplitude w, the nonlinear Chern number also depends on w as shown in Fig. 1.Therefore, the nonlinear Chern number can predict the nonlinearity-induced topological phase transition by the change of the amplitude in nonlinear systems, which is a qualitatively novel phenomenon absent in linear systems [22,24,[38][39][40][41][42].
In continuum systems with a periodic potential, the Bloch ansatz should read Ψ j (r) = e ik•r ψ j (k, r) with ψ j (k, r) being a periodic function of r whose period is equal to that of the periodic potential.Then, the wavenumberspace representation of the nonlinear eigenequation becomes f j (k, ψ(k); r) = E(k)ψ j (k, r), and the nonlinear Chern number can be written as , where S represents the unit cell of the periodic system.The squared amplitude is defined as w = j S d 2 r|ψ j (k, r)| 2 in this continuum case.
Compared to a previous study of the nonlinear topological invariant in one-dimensional systems [42], the nonlinear Chern number has no higher-order correction terms.This is because we assume the U (1) symmetry which guarantees that one can describe periodic steady states by using a single frequency mode.Therefore, there is no need to consider the multi-mode expansion of nonlinear eigenvectors that leads to the nonlinear collection terms.In fact, the one-dimensional topological invariant defined in the previous study is also reduced to the conventional form under the U (1) symmetry, which is consistent with our results.
We note that the Bloch ansatz does not describe bulklocalized solutions that can be obtained in strongly nonlinear systems.Since such strong nonlinearity can also generate edge-localized modes [34], it is not straightforward to identify whether the edge modes originate from the bulk topology or the nontopological nonlinear effects, which makes the bulk-boundary correspondence unclear in strongly nonlinear regimes.Therefore, in the following sections, we mainly fo- cus on weakly and more strongly nonlinear systems where the nonlinear terms are smaller than the linear terms (see also Supplementary Information for the Bloch-wave-like solutions of the bulk modes in this parameter region).
Nonlinear Chern number calculated from exact solutions of a lattice model.-Toinvestigate the bulk-boundary correspondence, i.e., the correspondence between the nonzero nonlinear Chern number and the existence of the edge-localized steady state, we propose and analyze the nonlinear extension of the Qi-Wu-Zhang (QWZ) model [44] defined in the real space (see Supplementary Information for the detail of the real-space description of the model).By using the Bloch ansatz, we rigorously obtain its wavenumber-space description, where w is the squared amplitude u and κ are dimensionless parameters of linear and nonlin- ear mass terms.We here introduce the staggered Kerr-like nonlinearity ±κw to the linear Chern-insulator model [44].
To calculate the nonlinear Chern number, we focus on special solutions where the squared amplitude w is fixed independently of the wavenumber k.Then, we can regard Eq. ( 4) as a linear equation and thus can diagonalize it as the linear QWZ model.Analytically diagonalizing the right-hand side of Eq. ( 4), we obtain the exact bulk solutions of the nonlinear eigenvalues and eigenvectors.Then, using the exactly obtained nonlinear eigenvectors, we calculate the nonlinear Chern number and obtain the phase diagram in Fig. 2b (see also Supplementary Information for the detail of the calculation and Supplementary Information for the numerical confirmation).The amplitude dependence of the nonlinear Chern number indicates the existence of the nonlinearity-induced topological phase transition in the nonlinear QWZ model.We note that the above nonlinear eigenvectors are exact solutions of the nonlinear QWZ model under the periodic boundary conditions, and thus our result reveals the existence of the nonlinearity-induced topological phase transition without any approximations.Such an analytical demonstration of the nonlinearity-induced topological phase transition is achieved by considering nonlinear equations with the form, (5) with ψ j and k being the state variables and the wavenumber (see also Supplementary Information), while the addition of off-diagonal or non-uniform diagonal nonlinear terms prevents us to obtain the exact solutions.It is also noteworthy that the obtained exact solutions are consistent with the results of the perturbation or self-consistent calculations of the nonlinear eigenvalue problem (see also Supplementary Information).
Bulk-boundary correspondence in weakly nonlinear systems.-Wefirst numerically confirm the bulk-boundary correspondence in weakly nonlinear systems.We simulate the dynamics of the nonlinear QWZ model (Eq.( 4)) with weak nonlinearity where the nonlinear Chern number is the same as that in the linear limit κw → 0. In the topological phase C NL = ±1 (Fig. 2d), we find a long-lived localized state that corresponds to a topological edge mode in the QWZ model.Meanwhile, in the case of C NL = 0 (Fig. 2c), the edgelocalized initial state is spread to the bulk, which indicates the absence of edge modes.We also confirm the bulk-boundary correspondence from the perspective of the nonlinear band structure as shown in Fig. 3 (see Supplementary Information for the details of the numerical method).
In fact, the bulk-boundary correspondence between the nonlinear Chern number and the gapless edge modes can be established in general weakly nonlinear systems.We mathematically show the bulk-boundary correspondence under weak nonlinearity compared to the linear band gap.We describe the detail of the theorem and its proof in Supplementary Information.
Nonlinearity-induced topological phase transition by stronger nonlinearity.-Wenext show that nonlinearityinduced topological phase transitions occur in the stronger nonlinear regime, where the nonlinear Chern number becomes nonzero and topological edge modes appear at a critical amplitude.To analyze the behavior of topological edge modes, we derive the effective theory of the low-energy dispersion of the nonlinear QWZ model around the critical amplitude.For example, if we focus on the critical amplitude w c = (−2 − u)/κ, the nonlinear band structure of the model closes the gap at (k x , k y ) = (0, 0).Then, around the critical amplitude and the gap closing point, w ∼ w c and (k x , k y ) ∼ (0, 0), we remain the leading order term of the wavenumber-space description of the nonlinear QWZ model.Finally, by substituting the wavenumbers with the derivative, we derive the real-space description of the continuum model m is positive and κ is negative, the nonlinear Dirac Hamiltonian exhibits the nonlinearity-induced topological phase transition.We obtain an unconventional localized mode if the amplitude is larger than a critical value.In this localized mode, there exist nonvanishing amplitudes even in the limit of x → ∞.Therefore, the localized mode can be unphysical in a truly semi-infinite system because we cannot normalize such a nonvanishing mode, while it still exists and is physically relevant in experimentally realizable finite systems.We set m, κ, and D as m = 0. i∂ t Ψ(r) = Ĥ(Ψ(r))Ψ(r), with m being m = u + 2 and ψ(x, y) = (ψ 1 (x, y), ψ 2 (x, y)) T being the state-vector function at the location (x, y).This statedependent Hamiltonian has a similar structure to the Dirac Hamiltonian except for the nonlinear mass term κ(|Ψ , and thus we term it the nonlinear Dirac Hamiltonian.Starting from other critical amplitudes w c = −u/κ and w c = (2 − u)/κ, one can derive similar state-dependent Hamiltonians (see also Supplementary Information).In general, the nonlinear Dirac Hamiltonian should describe the lowenergy dispersion of nonlinear topological insulators, and thus its localized modes unveil the existence of topological edge modes in various continuum systems.By using the Bloch ansatz Ψ i (x, y) = ψ i (k) exp(i(k x x + k y y)) (without explicit (x, y) dependence because of the continuous translational symmetry), one can determine the nonlinear Chern number of the nonlinear Dirac Hamiltonian as where w = S dS[|Ψ 1 (x, y)| 2 +|Ψ 2 (x, y)| 2 ]/|S| is the average squared amplitude of plain waves in this nonlinear system.We note that the nonlinear Chern number of the nonlinear QWZ model corresponds to the sum of those of the nonlinear Dirac Hamiltonians obtained from the expansion around the gapclosing points (k x , k y ) = (0, 0), (0, π), (π, 0), (π, π).
We analytically show that the Chern number predicts the existence of localized modes in continuum systems by calculating the localized modes of the nonlinear Dirac Hamiltonian.We assume the ansatz Ψ i (x, y) = e ikyy Ψ ′ i (x) that is periodic in the y direction and consider the x > 0 region with the open boundary at x = 0. We calculate the localized mode with the wavenumber k y and the eigenvector E = k y .Constructing an analogy to the linear case, one can use an ansatz We can analytically ob- The legend shows the lattice constants h.We can confirm that the size of the gap decreases as the system size becomes larger.However, if the lattice constant is so large that the discretized system cannot imitate the behavior of the continuum nonlinear Dirac Hamiltonian, we find a sudden decrease in the size of the band gap.This sudden decrease contradicts the estimation that the band gap is proportional to 1/(L + 1/2) and indicates the existence of irregular gapless modes.Since such gapless modes appear at smaller amplitudes than the critical point of the nonlinear Chern number, the numerical result indicates that the bulkboundary correspondence can be ensured only after taking both the continuum and thermodynamic limits.
tain the solution of this equation, where D is the integral constant and −(κ/m) + De −2mx must be positive.Combining w and D by the averaged squared amplitude w ), we find the bulkboundary correspondence, i.e., correspondence between the positive nonlinear Chern number in Eq. ( 8) and the existence of a left-localized mode in Eq. ( 9).We note that L can take an arbitrary value because the local amplitude w(x) = 1/ −(κ/m) + De −2mx of the left-localized mode also exists in the topological parameter region, m + κw(x) < 0.
Figure 4 summarizes the behaviors of gapless modes in the nonlinear Dirac Hamiltonian at different parameters.In the case of m < 0, where the Chern number is C = 1/2 in the linear limit, we obtain the localized states at the left side as in the linear case.These localized states are consistent with the bulk-boundary correspondence in weakly nonlinear systems shown in the previous section.If we consider m > 0, κ < 0, we still obtain the localized state, while the amplitude remains to be |m/κ| at the limit of x → ∞.The residual amplitude |m/κ| corresponds to the transition point of the nonlinear Chern number that satisfies m + κw = 0. Therefore, the localized state at m > 0, κ < 0 indicates the nonlinearity-induced topological phase transition associated with the amplitude-dependent Chern number.We note that the nonvanishing amplitude at the limit of x → ∞ indicates that it is impossible to normalize the edge mode.Meanwhile, in finite systems, such a nonvanishing localized mode can be normalized and thus can robustly emerge.We can also check that no localized modes appear in the case of positive m and κ, where the nonlinear Chern number is C NL = −1/2 at any amplitudes.
In lattice systems, we numerically validate the existence of the amplitude-dependent gapless modes and their correspondence to the nonlinear Chern number.We reveal that the finite-size effect leads to nonzero gaps of the localized modes, while such band gaps converge to zero in the thermodynamic limit L → ∞.We also find that discontinuity of edge modes can alter the phase boundary, and thus lattice systems exhibit the bulk-boundary correspondence at the parameters where they well-approximate the continuum models.To show that, we consider the rediscretization of the continuum nonlinear Dirac Hamiltonian in Eq. ( 7) and reveal the correspondence between the parameters in the nonlinear QWZ model (Eq.( 4)) and the lattice constant (see Supplementary Information for the detail).Using the small lattice constant h, one can assume that the discretized Hamiltonian reproduces the behavior of the continuum nonlinear system.
To numerically confirm the bulk-boundary correspondence in the lattice model, we calculate the minimum of the absolute values of the nonlinear eigenvalues at different sizes.We fix the parameters m = −1 and κw/L = 1, where the nonlinear band structure is gapless in the corresponding continuum model.Figure 5 shows the size and lattice-constant dependencies of the energy gaps.We obtain nonzero gaps even at the topological parameters C NL ̸ = 0 due to the finite-size effect.However, the gap becomes smaller as the system size becomes larger.Specifically, if the lattice constant h is small enough to reproduce the behavior of the continuum Hamiltonian, we confirm that the gap is proportional to 1/(L + 1/2) with L being the system size, which corresponds to the analytical solution of the eigenvalue equation of the linear QWZ model (see Supplementary Information).Thus, the gap will be closed in the continuum and thermodynamic limit.Meanwhile, we find sudden decreases in the sizes of the band gaps at h > 0.14.The inconsistency between the numerical results at large lattice constants and the analytical estimation in the continuum limit indicates that the bulk-boundary correspondence can be modified by strong nonlinearity (see Supplementary Information for additional numerical calculations).
Observation protocol of nonlinear edge modes via the quench dynamics.-Inrealistic setups, it is difficult to directly prepare a single nonlinear edge mode because it has a complicated amplitude distribution and thus one needs finetuning of the strength of the excitation of each site.Instead, one can observe the topological properties via quench dynamics.In quench dynamics, one only has to excite the edge sites at homogeneous amplitudes and observe the dynamics without any other external interactions as depicted in Fig. 6a.
To confirm the correspondence between the existence of the nonlinear edge modes and the localized states in quench dynamics, we numerically simulate quench dynamics of the nonlinear QWZ model (Eq.( 4)) at various parameters.Figures 6b and c show the time evolution of the quench dynamics with and without nonlinear edge modes.We also obtain the phase diagrams in Figs.6d and e, which are classified by the amplitude at the edge of the sample in the long-time limit.Here, we consider two initial states equivalent to the nonlinear edge modes at u+κw = ±1 (see also Supplementary Information).In weakly nonlinear regimes, the localized states remain in the topological cases and vanishes in the trivial cases.We can also confirm the shift of the phase boundary by the change of the amplitude, which indicates the existence of the nonlinearityinduced topological phase transitions.We note that strong nonlinearity induces trap phases [34,41,45] where the localized initial states remain without topological origins (see also Supplementary Information).It is also noteworthy that the phase boundary obtained from quench dynamics does not exactly match that of the nonlinear Chern number because the perfectly-localized initial condition is quite different from the Bloch wave and thus quench dynamics cannot fully reproduce the topological properties obtained from the Bloch ansatz.Specifically, previous studies have argued topological photonics using ring-resonator arrays [46,47].By focusing on a pseudo-spin sector in light propagating through ring-resonator arrays, one can realize the photonic counterpart of a Chern insulator.If we utilize different materials in the ring resonators at different sublattices, one can possibly tune the sign of the Kerr nonlinear effect [48] at each site.Thus, topological photonic metamaterials using different resonators can be a candidate system to realize the nonlinear topological insulators analyzed in the present paper.
Discussion.-We introduced the nonlinear Chern number in two-dimensional systems by considering the nonlinear eigenvalue problem (Eq.( 2)).We theoretically proved the bulk-boundary correspondence of the nonlinear Chern number in weakly nonlinear regimes, which guarantees that if the nonlinearity is small enough compared to the bulk band gap, the bulk topology can predict the existence of the edge modes.More importantly, we investigated the bulk-boundary correspondence in the stronger nonlinear regime where the nonlinearity is larger than the linear band gap while it is smaller than the linear couplings.We showed the existence of the nonlinearity-induced topological phase transition that depends on the amplitude and thus has no counterparts in linear systems.We proposed a minimal model of a nonlinear Chern insulator, whose exact bulk solutions indicate the existence of the nonlinearity-induced topological phase transition detected by the nonlinear Chern number.We analytically show that the nonlinear Chern number exactly predicts the nonlinearityinduced topological phase transition in the nonlinear Dirac Hamiltonian, implying the nonlinear bulk-boundary correspondence..In a nonlinear lattice system, we numerically checked that the bulk-boundary correspondence is recovered in the continuum and thermodynamic limit.
Our results indicate the existence of unique topological phenomena beyond the weakly nonlinear regime.There remain intriguing issues to fully establish the topological classification of nonlinear systems including further strongly nonlinear cases where nonlinearity is even stronger than the linear couplings.Since such strong nonlinearity can induce bulklocalized modes [18,34,41,45] that are out of the description of the Bloch ansatz, it is unclear whether or not the nonlinear Chern number still fully works.
From the perspective of many-body quantum physics, nonlinear systems are regarded as the result of the mean-field approximation of interacting many-body systems [16,17].The nonlinear Chern number should correspond to the topology of the excited states of interacting systems, which we can observe by applying external fields oscillating at the corresponding frequency.In fact, the nonlinear eigenvalues correspond to the energies per particle of the many-body wavefunctions in the mean-field analysis (see Supplementary Information), which are not the ground state except for that of the minimum eigenvalue.Therefore, the nonlinear topology would reveal the unexplored topological phases in quantum many-body systems at the level of mean-field approximation.

Supplementary Materials Justification of the Bloch ansatz via the perturbation analysis.
If the nonlinear term is small compared to the linear band gap, one can regard the nonlinear effect as a perturbation to the linear band structure.Under such an assumption, one can perturbatively calculate the nonlinear eigenvectors.To show the perturbation-calculation protocol of the nonlinear eigenvalue problem we rewrite the nonlinear eigenequation as We consider the perturbation expansion by κ, One can determine Ψ (0) and E (0) from the eigenvalue and eigenvector of the linear Hamiltonian H 0 Then, the first-order perturbation is calculated from the eigenequation of (H 0 + κH NL (ψ (0) )) as ))(Ψ (0) + κΨ (1) ) = (E (0) + κE (1) )(Ψ (0) + κΨ (1) ).
In translation-invariant systems, the nonperturbed eigenvector Ψ (0) is described by a Bloch wave Ψ (0) = e ikx ψ (0) , due to the Bloch theorem of linear systems.Since the Bloch wave exhibits no site-dependence of the amplitude, one can assume that the nonlinear term κH NL (Ψ) is also uniform, and thus the whole effective Hamiltonian H 0 + κH NL (Ψ) still has the translational symmetry.Therefore, in weakly nonlinear systems, one can believe that all of the nonlinear eigenvectors are described by the Bloch ansatz.We note that in strongly nonlinear regimes, there can be localized modes that cannot be described by the Bloch ansatz [34,41,45].However, the periodic solutions obtained from the Bloch ansatz are still exact nonlinear eigenvectors under the periodic boundary conditions.

Real-space description of the nonlinear QWZ model.
To investigate the existence or absence of topological edge modes in lattice systems, we construct a minimal lattice model of a nonlinear Chern insulator, which we term the nonlinear QWZ model (see Eq. ( 4) for the wavenumber-space description).Its real-space dynamics is described as where j, l and x, y represent the internal degree of freedom and the location, respectively, and Ψ j (x, y) is the j-th component of the state vector at the location (x, y).(σ i ) jl 's are the (j, l)-component of the ith Pauli matrices.This lattice model introduces the staggered Kerr-like nonlinearity −(−1) to the linear QWZ model [44].
Exact bulk solutions of the nonlinear QWZ model.
To obtain the phase diagram in Fig. 2b, we analytically solve the nonlinear eigenequation in Eq. ( 4).If we focus on special solutions where the squared amplitude |Ψ 1 (k)| 2 + |Ψ 2 (k)| 2 = w has no k-dependence, the nonlinear eigenequation (Eq.( 4)) exactly corresponds to a linear one.Therefore, by solving the corresponding linear eigenequation, we obtain the following exact bulk eigenvalues and eigenvectors, where c ± (k) = (u + κw + cos k x + cos k y + E ± (k x , k y )) 2 + sin 2 k x + sin 2 k y is a normalization constant.By using these nonlinear eigenvectors, we analytically obtain the nonlinear Chern number, as summarized in Fig. 2b.We note that one can generally obtain exact bulk solutions if the nonlinear equation has the form in Eq. ( 5) (see also the following section).

Derivation of exact bulk solutions.
If the nonlinear dynamics is described as one can obtain its exact bulk nonlinear eigenvectors.Here, ψ n (k) is the state variable in the momentum space k, and n, m and j label the unit cell or the internal degree of freedom such as the spin satisfying The basic idea is to construct a set of special exact solutions self-consistently by requiring that the quantity w defined by has no k dependence.Then, the nonlinear eigenequation becomes a linear equation which is diagonalizable by using eigenfunctions ψ together with the eigenenergy E (α) .We note that is also an eigenfunction of the same eigenequation in the linear case, where c (k) is an arbitrary function.When we set c (k) independently of k as we obtain It follows from Supplementary Eqs.(S12) and (S16) that which indicates that Supplementary Eq. (S14) is the nonlinear eigenvector of Supplementary Eq. (S12).Hence, we have solved Supplementary Eq. (S12) self-consistently by requiring Supplementary Eq. (S11).Supplementary Equation (S14) presents a set of exact solutions parametrized by the real number w in terms of solutions of the linear equation (S12).
Perturbation analysis of the bulk modes of the nonlinear QWZ model.
While we calculate the exact solutions of the bulk modes of the nonlinear QWZ model in the main text, we can also obtain the same bulk modes from the perturbation analysis or the self-consistent calculations.Specifically, if we conduct the perturbation analysis described in the previous section, the calculation stops at the first-order perturbation and derives the same bulk modes as the exact solutions.
The zeroth-order solutions, i.e., the linear solutions of the bulk modes of the QWZ model are where we fix the norm as Then, substituting these solutions into the state-dependent Hamiltonian of the nonlinear QWZ model, one can obtain the first-order-perturbation solutions.Due to the nonlinear terms only depending on the norm of the nonlinear eigenvector, the substituted effective Hamiltonian is described as independently of the wavenumber.The eigenvalues and eigenvectors of this Hamiltonian are the same as those of the nonlinear QWZ model (Eqs.(S7) and (S8)), and thus the first-order perturbation calculation is consistent with the exact solutions.The self-consistent calculation is equivalent to the higher-order perturbation calculation (see also the following section).However, if one substitutes the first-order solutions into the state-dependent Hamiltonian of the nonlinear QWZ model, one obtains the same effective Hamiltonian as Eq.(S21).Therefore, the nonlinear eigenvectors and eigenvalues obtained from the self-consistent calculation are the same as those obtained from the first-order perturbation and the exact solutions.
Numerical simulations of the real-space dynamics of the nonlinear QWZ model.
In Fig. 2, we numerically calculate the dynamics of the nonlinear QWZ model by using the fourth-order Runge-Kutta method.We consider the 20 × 20 square lattice, where each lattice point has two internal degrees of freedom.We impose the open boundary condition in the x direction and the periodic boundary condition in the y direction.We set the time step dt = 0.005.The simulation starts from the localized initial states where We use the parameters u = 3, κ = 0.1, w = 1 in Fig. 2c and u = −1, κ = 0.1, w = 1 in Fig. 2d.In these figures, we plot the square root of the sum of the square of absolute values of the first and second components.

Self-consistent calculation of nonlinear band structures.
To obtain the nonlinear band structures in Fig. 3, we numerically calculate the nonlinear eigenvalue problem by using the self-consistent method.We first rewrite the nonlinear dynamics as i∂ t Ψ = f (Ψ) = H(Ψ)Ψ, where H(Ψ) is a state-dependent effective Hamiltonian.Then, we conduct the self-consistent calculation in the following procedure: (i) We numerically diagonalize H( ⃗ 0) and set the initial guess of the eigenvalue and eigenvector Ψ 0 , E 0 by adopting a pair of the obtained eigenvalue and eigenvector of H( ⃗ 0).We fix the norm of Ψ 0 to be ||Ψ 0 || 2 = w.(ii) We substitute the guessed eigenvector Ψ i after i iterations into H(Ψ) and diagonalize H(Ψ i ).(iii) We choose the obtained eigenvalue that is the closest to the previous guess E i , and the corresponding eigenvector as the next guess E i+1 , Ψ i+1 .(iv) We iterate the steps (ii) and (iii) until the distance between Ψ i and Ψ i+1 becomes smaller than the threshold, ||Ψ i+1 − Ψ i || < ϵ, or the iteration reaches a fixed number.We also perform these calculations starting from all the eigenvectors of H( ⃗ 0) and obtain a set of nonlinear eigenvectors and eigenvalues of f (Ψ).
In Fig. 3, we consider the parameterized state-dependent Hamiltonian that corresponds to the nonlinear Chern insulator under the assumptions of the y-periodic Bloch ansatz and the open boundary condition in the x direction.The state-dependent Hamiltonian is described as where ∆ x and ∆ 2 x are the difference operators defined as ∆ x Ψ(x) = [Ψ(x + 1) − Ψ(x)]/2 and ∆ 2 x Ψ(x) = Ψ(x + 1) + Ψ(x − 1) − 2Ψ(x), respectively.Then, we calculate the nonlinear eigenvalues of this state-dependent Hamiltonian at k y = n∆k, n = −N, −N + 1, • • • , N − 1, N , and ∆k = π/N with N being N = 50.
To calculate the nonlinear band gap in Fig. 5, we also use the self-consistent method.However, to stably obtain the band structure of topological gapless modes, we start from the eigenvectors of H(( w/2L, w/2L, • • • ) T ) instead of those of H( ⃗ 0), where L is the system size.This is because H( ⃗ 0) has no topological gapless modes, and thus the initial guess corresponding to the gapless modes in H(Ψ) cannot be obtained from H( ⃗ 0).We also note that the self-consistent calculation is only stable at weak nonlinearity.To calculate the band structures in strongly nonlinear systems than those analyzed in the main text, one should instead use the Newton method, whose details are described in the following section.
Quasi-Newton method to solve strongly nonlinear eigenvalue problems.
Here, we explain the (quasi-)Newton method to solve nonlinear eigenvalue problems in strongly nonlinear regimes.While we use the self-consistent calculation in weakly nonlinear systems, it often fails to converge in strongly nonlinear regimes.Instead, we should use the (quasi-)Newton method [51], which approximately calculates the roots of algebraic equations.To apply the Newton method, we reconsider the set of the nonlinear eigenvalue problem and the normalization condition, If we consider the n-component vector Ψ, there are n + 1 unknown variables consisting of the set of the components of the eigenvector Ψ i and the eigenvalue E.Then, the nonlinear eigenequation provides n algebraic equations and the normalization condition does one algebraic equation.Therefore, the total number of algebraic equations is equal to the number of unknown variables, and thus we can solve the eigenvalue problem by calculating the roots of the algebraic equations.The Newton method solves the algebraic equations by using the information of the residual and the local gradient.We rewrite the equations as F (Ψ, E) = 0, where 2 −w.We first set an initial guess of the eigenvalue E 0 and eigenvector Ψ 0 .In our calculations, we randomly determine the initial guess of the eigenvalue from the uniform distribution [−2, 2] and the real part and the imaginary part of each component of the initial guess of the eigenvector from the uniform distribution [−1, 1].Then, we renormalize the initial guess of the eigenvector.If the guess is close to the genuine solution Ψ g , E g , we can approximate the residual as where J F is the Jacobian matrix of F .Thus, (Ψ 1 , E 1 ) = (Ψ 0 , E 0 )+αJ −1 F F (Ψ 0 , E 0 ) with α being the update step becomes better approximations of the nonlinear eigenvector and eigenvalue.Iterating this update procedure of the approximated eigenvector and eigenvalue, we can numerically calculate the nonlinear eigenvector and eigenvalue.In our calculations, we fixed the update step as α = 0.01.If we start from different initial guesses, we can obtain other sets of eigenvectors and eigenvalues, which can cover all the eigenvectors of the nonlinear equation.
Since the Jacobian matrix of F is tough to calculate in our model, we also numerically obtain the approximation of the Jacobi matrix, which is known as the quasi-Newton method.Specifically, we use an update algorithm called the Broyden method [52].We denote the approximated Jacobian matrix at the ith step of the quasi-Newton method by B i .We set B 0 as the identity matrix, B 0 = I.Then, the approximated Jacobian matrix is updated as where y i and ∆Ψ ′ i are Finally, we substitute B i into J F in Supplementary Eq. (S26) and update the guess of eigenvectors and eigenvalues.

Theorem of the bulk-boundary correspondence in weakly nonlinear systems.
We mathematically show the bulk-boundary correspondence in weakly nonlinear systems.We here use a simple notation , where ⃗ Ψ is the nonlinear eigenvector whose components correspond to the state variables at the locations r and the internal degrees of freedom j.The claim of the theorem is as follows: Ψ is a nonlinear eigenvalue problem on a two-dimensional lattice system that satisfies the following assumptions: (1) When we rewrite the nonlinear function f as f ( ⃗ Ψ) = H( ⃗ Ψ) ⃗ Ψ, there exists a positive real number c < 1 that satisfies ||H( ⃗ Ψ) − H( ⃗ 0)|| < gc/2 for any complex vector ⃗ Ψ with g being the bandgap of H( ⃗ 0) and || • || being the operator norm.(2) There exists a positive real number c ′ < 1 such that for any pairs of complex vectors Ψ and Ψ + ∆Ψ with the norm w, they satisfy (3) For any complex vector ⃗ Ψ, one can rewrite the nonlinear function f where H is a Hermitian matrix.(4) The nonlinear equation satisfies the U (1) symmetry, f (e iθ ⃗ Ψ) = e iθ f ( ⃗ Ψ).We also assume that the number of nonlinear eigenvectors is equal to that of the linear eigenvectors of H( ⃗ 0).Then, the nonlinear eigenequation f ( ⃗ Ψ) = E ⃗ Ψ exhibits robust gapless boundary modes if and only if its nonlinear Chern number is nonzero.
To prove Theorem, we first show the following proposition: Proposition 1: In the nonlinear eigenvalue problem that satisfies the assumptions in Theorem, the self-consistent calculation converges, Furthermore, there exists an eigenvector ⃗ Ψ 0 and eigenvalue E 0 of H( ⃗ 0) that satisfy To show this proposition, we utilize the perturbation theorem of the eigenvectors in linear systems: Suppose H as a nondegenerate Hermitian matrix and g as the minimum difference between its two eigenvalues.If ||A|| < g/2 in terms of the operator norm, for an arbitrary eigenvector ⃗ Ψ of H + A, there exists an eigenvector ⃗ . We iteratively use this theorem and evaluate the distance between the guess of eigenvectors ⃗ Ψ i and ⃗ Ψ i+1 at each step.
We also show that the resulting eigenvalue and eigenvector, E ∞ and ⃗ Ψ ∞ are indeed a pair of nonlinear eigenvalue and eigenvector of f ( ⃗ Ψ), which is summarized in Proposition 2: The converged solutions of the self-consistent calculation, E ∞ and ⃗ We show this proposition by using simple inequalities and limit evaluations.By using these propositions, we finally show the following lemma to prove Theorem.
Lemma: For an arbitrary eigenvector of a nonlinear eigenvalue problem , where w is the norm of ⃗ Ψ and ⃗ Ψ ϵ , and C is the constant independent of the eigenvector ⃗ Ψ and the constant ϵ.The lemma indicates that in weakly nonlinear systems, one can map the nonlinear eigenvalues and eigenvectors onto those of a linear eigenvalue problem.Thus, we can show the bulk-boundary correspondence in weakly nonlinear systems.

Proof of the main theorem.
In this section, we prove Theorem in the main text.As discussed in the previous section, we use two propositions and a lemma to prove the theorem.In Propositions 1 and 2, we show that the nonlinear eigenvectors are obtained as the convergent values of the self-consistent calculation.We also show Lemma by using a similar technique in Proposition 1. Lemma indicates the existence of a continuous path that connects nonlinear and linear eigenvectors.Therefore, one can conclude the bulk-boundary correspondence in weakly nonlinear systems via that in linear systems.In this section, we use Ψ to denote the nonlinear eigenequation, where E and ⃗ Ψ are the nonlinear eigenvalue and eigenvector, and H( ⃗ Ψ) represents the effective Hamiltonian that depends on the nonlinear eigenvector.

Proof of the perturbation theorem of linear eigenvectors.
We utilize the perturbation theorem of linear eigenvectors to prove Proposition 1.The perturbation theorem is as follows: Suppose H as a nondegenerate Hermitian matrix and g as the minimum difference between its two eigenvalues.If ||A|| ∞ < g/2 in terms of the operator norm, for an arbitrary eigenvector ⃗ Ψ of H + A, there exists an eigenvector ⃗ . The perturbation theorem is proved as follows.Let Ψ and E be the eigenvector and eigenvalue of H + A. From Weyl's perturbation theorem [53], there exists only one pair of an eigenvector and an eigenvalue, ⃗ where we rewrite the difference of eigenvectors and eigenvalues as ∆ ′ ⃗ Ψ = ⃗ Ψ − ⃗ Ψ ′ 0 and ∆E = E − E 0 , respectively.Then, we separate ∆ ′ ⃗ Ψ into the components parallel and perpendicular to ⃗ with (H − E 0 ) − being the generalized inverse matrix of H − E 0 .The operator norms of ∆E − A and Therefore, we obtain the bound of the norm of ∆ ⃗ Ψ as Then, we determine and its norm is bounded by which indicates the existence of the eigenvector ⃗ Proof of Proposition 1.
Next, we prove Proposition 1.We assume that the nonlinear eigenvalue problem f ( ⃗ Ψ) = E ⃗ Ψ satisfies the following conditions in Theorem: (1) When we rewrite the nonlinear function Then, Proposition 1 reads as follows: In the nonlinear eigenvalue problem that satisfies the assumptions in Theorem, the self-consistent calculation converges, Furthermore, there exists an eigenvector ⃗ Ψ 0 and eigenvalue We denote the eigenvector obtained in the ith step of the self-consistent calculation by ⃗ Ψ i .Mathematically, ⃗ Ψ i is defined as follows: whose eigenvalue is the closest to the eigenvalue of ⃗ Ψ i−1 and that minimizes the distance between ⃗ Ψ i−1 and ⃗ Ψ i .We also denote the band gap of H( ⃗ Ψ i ) as g i+1 .
To prove Proposition 1, we iteratively use the perturbation theorem of linear eigenvectors and evaluate the distance between ⃗ Ψ i−1 and ⃗ Ψ i .In the following, we explain how to evaluate the distance at each step.From the first condition, we obtain an inequality ||H( ⃗ Ψ 0 ) − H( ⃗ 0))|| < gc/2.Then combining the inequality and the perturbation theorem, we show the upper bound of the distance between ⃗ Ψ 0 and ⃗ The linear perturbation theorem [53] also indicates g 1 > g(1 − c).We next use the second condition, which indicates Combining these inequalities and the perturbation theorem of linear eigenvectors, we obtain where we use 0 < c < 1 in the last inequality.From the conditions c < 1 and c ′ < (1 − c)/2 < 1 and the linear perturbation theorem, we also obtain We again use the second condition and derive From the perturbation theorem, we obtain Then, we obtain the bounds To continue the calculations, we utilize the inequality g By iterating the evaluation of the upper bound of the band gap g i and the distance between ⃗ Ψ i−1 and ⃗ Ψ i , we obtain The above inequality indicates that the distance between ⃗ Ψ i and ⃗ Ψ 0 is also bounded by a constant Therefore, this infinite series is a Cauchy sequence, and thus ⃗ Ψ i converges to lim i→∞ ⃗ Ψ i .Supplementary Equation (S38) also presents the upper bound of the distance between lim i→∞ ⃗ Ψ i and ⃗ Ψ 0 .We can also obtain the bound of the difference between the nonlinear eigenvalue and the corresponding linear eigenvalue from half of the right-hand side of Supplementary Eq. (S36).
We note that we have loosely (but rigorously) evaluated the bound in Supplementary Eqs.(S34) and (S35).To obtain the tight inequality, we need to solve a complicated recurrence relation and thus we leave it to future works.While the obtained bound is loose, Proposition 1 still guarantees that the algorithm of the self-consistent calculation can solve the nonlinear eigenvalue problem in weakly nonlinear systems.It is also noteworthy that one can also prove the convergence of the self-consistent calculation in the case (1 − c)/2 ≤ c ′ < 1 in the completely same way.The condition of c ′ < (1 − c)/2 is needed only for proving Lemma.

Proof of Proposition 2.
While Proposition 1 guarantees the convergence of the self-consistent calculation, there still remains a possibility that the convergent solution is not a pair of genuine nonlinear eigenvalue and eigenvector.Thus, we need to prove Proposition 2, which guarantees that the obtained solution satisfies the nonlinear eigenequation.
Proposition 2: The converged solutions of the self-consistent calculation, Ψ i and E i be the obtained solution of an eigenvalue and an eigenvector at the ith step of the self-consistent calculation, respectively.We also consider their limits, ⃗ Ψ = lim i→∞ ⃗ Ψ i and E = lim i→∞ E i .Then, Proposition 2 is equivalent to The norm of the left-hand side is upper-bounded by Since we assume that H( ⃗ Ψ) is always a bounded operator, the first term converges to zero.By using the second condition in Theorem, we obtain the upper bound of the second part as 2) with g n being the minimum difference of the nonlinear eigenvalues of H( ⃗ second term also converges to zero.Lastly, one can also check the convergence of the third term to zero from the definitions of E and ⃗ Ψ.Therefore, the right-hand side of Supplementary Eq. (S39) converges to zero, which deduces that the left-hand side also converges to zero in the limit of n → ∞.

Proof of Lemma.
The proof of Lemma is similar to that of Proposition 1.The major difference is that we use the nonlinear eigenvector ⃗ Ψ of H + f ( ⃗ Ψ) as the initial guess of the nonlinear eigenvector of H + (1 − ϵ)f ( ⃗ Ψ).In this subsection, we define ⃗ Ψ ′ i as follows: The initial guess ⃗ The obtained solution at the ith step ⃗ We also describe the minimum difference of the eigenvalues of H + f ( ⃗ Ψ) by g ′ , which is lower-bounded by Supplementary Eq. (S36).The minimum difference of the eigenvalues of , we obtain the following inequality, Then, we obtain the bound of As in the proof of Proposition 1, we iteratively evaluate the bound of ||H( ⃗ for i ≥ 3. From the condition of c ′ < (1 − c)/2, we obtain gc ′ /2g ′ < 1.Therefore, Supplementary Equation (S45) indicates the convergence of the self-consistent calculation.We can also derive the upper bound of the distance between ⃗ Ψ and ⃗ Ψ ϵ , whose right-hand side is the product of ϵ and a constant that is independent of the eigenvector ⃗ Ψ.
Proof of Theorem.
To prove the bulk-boundary correspondence in weakly nonlinear systems that satisfies the conditions in Theorem, we show that the nonlinear eigenvectors can be continuously connected to linear eigenvectors.Here, we consider the nonlinear eigenvectors whose nonlinear eigenvalues are the closest to a certain eigenvalue E of H of all the nonlinear eigenvalues.Since we assume that the number of nonlinear bands is equal to that of linear bands of H, the nonlinear eigenvector of H + ϵf ( ⃗ Ψ) should match that obtained from the self-consistent calculation starting from the nonlinear eigenvector of H + ϵf ( ⃗ Ψ)/(1 − δ), which is considered in Lemma.Therefore, Lemma indicates that ⃗ Ψ ϵ is a continuous vector function of ϵ and thus provides a continuous path connecting the nonlinear eigenvector of H + f ( ⃗ Ψ) and the linear eigenvector of H.
Since the nonlinear Chern number is a topological invariant that cannot change its value under the continuous deformation of the nonlinear eigenvectors, the nonlinear Chern number coincides with the linear Chern number of H.One can also construct adiabatic connections between the edge modes of linear systems and nonlinear eigenvectors in nonlinear systems with open boundaries and nonzero nonlinear Chern numbers.Therefore, the bulk-boundary correspondence in weakly nonlinear systems is proved from that in linear systems.

Discretization of the nonlinear Dirac Hamiltonian.
To confirm the bulk-boundary correspondence in the nonlinear QWZ model with stronger nonlinearity, we discretize the nonlinear Dirac Hamiltonian in Eq. (7).Naively taking the lattice constant h, one seems to achieve the discretization of the Dirac Hamiltonian, where ∆ i are the difference operators defined as . However, at the gapless parameter 2 ) = 0, the obtained lattice Hamiltonian closes gaps at four points in the wavenumber space (k x , k y ) = (0, 0), (0, π), (π, 0), (π, π), which have different signs of topological charges around them.Thus, we cannot obtain the expected topological features such as the nonzero Chern number and the localized modes in this naive discretization.To avoid this inconsistency, we introduce the wavenumber-dependent mass terms ∆ 2 /h, where ∆ 2 is the discretized Laplacian, The added mass terms correspond to the Wilson fermion action [49,50] in the lattice field theory.We finally obtain the discretized Hamiltonian, which is denoted as in the wavenumber space, where w is the squared amplitude w As expected from the derivation of the nonlinear Dirac Hamiltonian in the main text, this discretized Hamiltonian is equivalent to that of the nonlinear QWZ model in Eq. ( 4) except for the existence of the lattice constant h.In the h → 0 limit, the obtained lattice Hamiltonian reproduces the behavior of the nonlinear Dirac Hamiltonian.

Numerical calculations of the quench dynamics.
We have numerically solved the nonlinear Schrödinger equation (Eq.(S6)) with the initial condition localized at the left edge as shown in Fig. 6a.As the initial conditions, we study the two cases.One is which corresponds to the solution with u = −1 of the linear model (κ = 0).The other is which corresponds to the solution with u = 1 of the linear model (κ = 0).The derivation is shown in Supplementary Note 10.
For the numerical calculation, we used the "NDSolve" function in Mathematica.We consider a 10 × 10 lattice under the open boundary condition in the x direction and the periodic boundary condition in the y direction.We solve the nonlinear Schrödinger equation (Eq.(S6)) under the two initial conditions in Eqs.(S49) and (S50).The time evolution of 6d and e.We define the phase indicator P by with T = 10.We have P ∼ 0 in the trivial phase as shown in Fig. 6b, which implies the absence of localized edge states.On the other hand, we have P ∼ 1 in a topological phase as shown in Fig. 6c, which implies the presence of a localized edge mode.In order to elucidate the topological phase diagram, we plot P in the u-κw plane in Fig. 6d and e.We find that the phase indicator is 1 along the blue line, where the exact solution is valid.It shows that the exact solution is realized by the quench dynamics.
Derivation of a nonlinear eigenvalue problem in the wavenumber space.
We use the wavenumber-space description of a nonlinear eigenequation to introduce the nonlinear Chern number.Here, we clarify the relationship between the real-space and wavenumber-space description of the nonlinear eigenequation.For clarity, we consider a one-dimensional lattice system with translation symmetry and a periodic boundary.In general, one can describe the eigenequation of such a one-dimensional nonlinear system as where ) represent the nonlinear eigenvalue and eigenvector, respectively.Due to the periodicity, ⃗ F does not explicitly depend on the location x.
To derive the wavenumber-space description of the nonlinear eigenequation, we have assumed the Bloch ansatz, with We have also assumed the U (1) symmetry denoted by where θ is a real phase factor.Substituting the Bloch wave ansatz in Supplementary Eq. (S52), we obtain Then, we finally derive the wavenumber-space description of the nonlinear eigenequation where both sides are independent of the location x.
One can also derive the wavenumber-space description of the nonlinear eigenequation corresponding to a continuum nonlinear system (cf.the nonlinear Dirac Hamiltonian in Eq. ( 7) in the main text).When a nonlinear system has continuous translation symmetry, the real-space description of an eigenequation is in a general continuum nonlinear system.To derive its wavenumber-space description, one should assume the Bloch-wave ansatz, We also impose the U (1) symmetry described as Then, the eigenequation (Supplementary Eq. (S57)) reads Therefore, one can obtain the wavenumber-space description of the eigenequation In continuum systems with discrete translation symmetry, the real-space description of an eigenequation is We assume the Bloch ansatz where ⃗ ψ ′ (x) is a periodic function with the same period as that of ⃗ F (x; •).The U (1) symmetry is described as Then, we rewrite the eigenequation (Supplementary Eq. (S62)) as Therefore, we obtain the wavenumber-space description of the eigenequation In both lattice and continuum cases, the size of the eigenvector in the wavenumber space is smaller than that in the real space, which is of practical advantage in the calculation of the nonlinear Chern number.We also note that one can derive the wavenumber-space description of nonlinear eigenequations in higher dimensions via the same calculations as one-dimensional cases.
In linear systems, one can derive the wavenumber-space descriptions of Hamiltonians from the Fourier transformation of the linear equations.In contrast, the Fourier transformation of a nonlinear equation reveals the interaction between eigenmodes with different wavenumbers.For example, in the Fourier transformation of a Kerr-type nonlinearity, |ψ(x)| 2 ψ(x) becomes However, since we assume the Bloch-wave ansatz ⃗ ψ(x) = e ikx ⃗ ψ ′ or equivalently the twisted boundary condition in a unit cell in the derivation of the wavenumber-space description of the nonlinear eigenequation, the integrant in Supplementary Eq. (S69) contributes only when the wavenumbers are

Quantization of the nonlinear Chern number.
The nonlinear Chern number defined in Eq. ( 3) in the main text are quantized.We can show the quantization by embedding the nonlinear eigenvectors into eigenspaces of a linear Bloch Hamiltonian.For clarity, we consider two-band cases in the following.
Let ψ(k) be the nonlinear eigenvectors of the nonlinear eigenvalue equation F (ψ) = Eψ derived from the Bloch ansatz of the wavenumber k.Then, one can uniquely determine ϕ(k) so that ϕ(k) is perpendicular to ψ(k) (we here ignore the phase ambiguity).By using ϕ(k) and ϕ(k), one can construct the following linear Bloch Hamiltonian This Bloch Hamiltonian has the eigenvector ψ(k) and the corresponding eigenvalue E(k) = −1.Therefore, the Chern number of the lower band of this linear Hamiltonian is which is equivalent to the nonlinear Chern number of F (ψ) = Eψ.Therefore, the nonlinear Chern numbers must be integers as linear ones.We note that in many-band cases, one can also construct a linear Bloch Hamiltonian that has the same eigenvectors, and thus the nonlinear Chern number must be quantized.

Amplitude distributions of nonlinear bulk modes.
Here, we numerically calculate the nonlinear eigenstates of the nonlinear Qi-Wu-Zhang (QWZ) model (Eq.( 4) in the main text) of nonlinear Chern insulators as in Fig. 3 in the main text.Supplementary Figure S1 shows an amplitude distribution of a nonlinear bulk eigenstate.The bulk eigenstate at k y = π/5 resembles sine curve, which can justify the use of the Bloch-wave ansatz to approximate the bulk modes and derive the nonlinear Chern number in this parameter region.In contrast, the bulk modes around k y = 0 are localized as shown in Supplementary Figs.S1c,f because of the nonlinear interactions between linear bulk modes with different wavenumbers.These results are consistent with Proposition 1, which indicates that the nonlinear eigenstates are almost unchanged from the corresponding linear modes in weakly nonlinear regimes, while stronger nonlinearity than the linear band gap can lead to drastic changes in nonlinear eigenvectors.
While numerical techniques cannot stably achieve periodic solutions around k y = 0, periodic solutions still exist in such a stronger nonlinearity region.To capture the topological properties of nonlinear systems, we focus only on such periodic solutions.We also note that weak nonlinearity has a negligible impact on the existence of gapless edge modes because their eigenvalues are isolated from those of bulk modes.
Numerical calculation of the Chern number in the nonlinear QWZ model.
In this section, we numerically demonstrate the quantization and the amplitude dependence of the nonlinear Chern number.We calculate the nonlinear eigenvectors of the wavenumber-space description of the nonlinear QWZ model (Eq.( 4) in the main text) as in Fig. 3 in the main text.Then, applying the Fukui-Hatsugai-Suzuki method [54] to the obtained nonlinear eigenvectors, we calculate the nonlinear Chern number.
Supplementary Figure S2a shows the result of the numerical calculation of the nonlinear Chern number.We confirm that the nonlinear Chern numbers always take integer values as in linear cases.The figure also indicates the existence of the nonlinearityinduced topological phase transitions where the nonlinear Chern numbers are altered at the critical amplitudes w = 0.5, 2.5, 4.5.While these nonlinear Dirac Hamiltonians have different signs in their wavenumber-dependent terms, the nonlinear Chern numbers are determined by the signs of mass terms similar to Eq. ( 8) in the main text.The sum of the nonlinear Chern numbers of four nonlinear Dirac Hamiltonians is equal to that of the nonlinear QWZ model.
System-size dependence of nonlinear band gaps in the nonlinear QWZ Hamiltonian in the continuum limit.
We analytically estimate the system-size dependence of nonlinear band gaps in the nonlinear QWZ model (Eq.( 4) in the main text) in the continuum limit at the critical amplitude where the nonlinearity-induced topological phase transition occurs.We impose the open boundary condition in the x direction and assume the uniform eigenvector in the y direction.Then, the eigenvalue equation of the discretized system is equivalent to that of the nonlinear Su-Schriefer-Heeger (SSH) model [23,28,38,[40][41][42]55], under a properly chosen unitary transformation of the effective Hamiltonian.We use the same parameters m, κ, and h as those in Eq. ( 7) in the main text.At the critical amplitude | ⃗ ψ(2x − 1)| 2 + | ⃗ ψ(2x)| 2 = |m|/κ, the strengths of intercell and intracell hoppings of the nonlinear SSH model become the same on average, we can estimate the band gap of the discretized nonlinear Dirac Hamiltonian with L unit cells from that of a simple one-dimensional chain with 2L lattice points.The eigenvector of the one-dimensional chain is described by sine curves ψ(x) = sin(kx).To satisfy the open boundary conditions ψ(0) = 0 and ψ(2L + 1) = 0, k must be a multiple of 2π/(2L + 1), i.e., k = 2nπ/(2L + 1) with n being an integer n = 1, • • • , 2L.The corresponding eigenvalues are E = a(2 + cos(2nπ/(2L + 1))) with a being a real constant, and thus we estimate the band gap around E = 0 as 2πa/(2L + 1) + O((2L + 1) −2 ).This estimation indicates that the nonlinear Dirac Hamiltonian is gapless at the critical amplitude in the thermodynamic limit, which is consistent with the bulk dispersion of the nonlinear Dirac Hamiltonian.
We also numerically estimate the nonlinear band gaps in the thermodynamic limit at different amplitudes.We calculate the minimums of the absolute values of eigenvalues at different amplitudes and system sizes.Then, we fit the squares of the minimums by a/(2L + 1) 2 + b, where a and b are the fitting parameters.As we discussed in the previous paragraph, b should become zero at the critical amplitude, while b should be positive below the critical amplitude, which indicates the existence of the band gap in the thermodynamic limit.Supplementary Figure S3 shows the obtained fitting parameters b at different amplitudes w.As we expect, we obtain positive b when the amplitudes are smaller than the critical one, and b is zero at the critical amplitude.Since the fitting function does not completely capture the finite-size scaling of spectrum gaps at noncritical amplitudes, b becomes negative at larger amplitudes than the critical one.However, this numerical result still indicates that the nonlinear QWZ model is gapped (gapless) when the nonlinear Chern number is zero (nonzero).
We note that the nonlinear Dirac Hamiltonian obtained from the low-energy effective theory of the nonlinear QWZ model is different from that in previous studies [19,56] because they have considered the nonlinear terms with the same sign at the first  4) in the main text) in a 10 × 1 supercell structure at u = −2.5 and κ = 10.The transition point is w = 0.5 with w being the squared amplitude, while gapless edge modes appear around w ≃ 0.15 as is pointed out by the red arrow.We note that some eigenvalues disappear at large amplitudes due to the limitation of the numerical technique.b, Amplitude distribution of an anomalous gapless mode.We confirm the localization of the gapless mode at the edge of the sample.The eigenvalue corresponds to that pointed by the red arrow on the left side of panel a. c, Amplitude distribution of a topological gapless mode.We confirm the localization of the topological mode at both edges of the system.One can understand this localized edge mode as the superposition of two edge modes localized at the left and right edges, respectively.The eigenvalue corresponds to that pointed by the green arrow on the right side of panel a.
main text, we can expect the emergence of the nonlinearity-induced topological phase transition, where topological edge modes appear at the critical amplitude.In fact, numerically calculating the dynamics of the nonlinear QWZ model above the critical amplitude, we obtain a long-lived localized state shown in Supplementary Fig. S4.Here, we conduct the Runge-Kutta simulation of the lattice model as in Fig. 2 in the main text, setting the time step as dt = 0.005 and parameters as u = −2.1 and κ = 1.
We also numerically calculate the nonlinear band structure corresponding to the lattice model.Since the self-consistent calculation cannot stably obtain the eigenvalues in the strongly nonlinear regimes, here we instead use the quasi-Newton method to solve the nonlinear eigenvalue problems.We impose the open boundary condition in the x direction and the periodic boundary condition in the y direction and calculate the nonlinear eigenvalues at each phase parameter k y to obtain the nonlinear band structure.Supplementary Figure S5 shows the nonlinear band structure at the same parameters in Supplementary Fig. S4.One can confirm the existence of the gapless modes.We note that if the amplitude is smaller than that considered in Supplementary Figs.S4, S5, there are no gapless edge modes, which indicates that the strongly nonlinear lattice system exhibits the nonlinearityinduced topological phase transition.
The numerical results in Supplementary Figs.S4, S5 indicate the existence of topological edge modes in stronger nonlinear lattice systems, while one can see the inconsistency between the parameters where edge modes appear and the nonlinear Chern number is changed.We confirm such inconsistency in the numerical calculations of the eigenvalues at k y = 0 and different amplitudes.Supplementary Figure S6a shows the numerical results at the parameters u = −2.5 and κ = 10 and various amplitudes.In this case, the nonlinear Chern number is changed at the amplitude w = 0.5, while the nonlinear band becomes gapless around w = 0.15, which implies the breakdown of the bulk-boundary correspondence.We check that the anomalous gapless modes in Supplementary Fig. S6b are also localized at the edge of the sample as in topological gapless modes in Supplementary Fig. S6c.
The emergence of the anomalous gapless modes is induced by the discontinuous change of the amplitude.At the amplitude w = 0.15, the local on-site term at x = 0 can become u + 1 + κ|Ψ(x = 1)| 2 = 0, which leads to the existence of perfectly localized modes exhibiting a sudden decrease in their amplitude.Comparing the nonlinear QWZ model and the nonlinear Dirac Hamiltonian (Eq.( 7) in the main text), a sudden change of the amplitude is unphysical in the nonlinear Dirac Hamiltonian.However, the effect of the momentum term becomes much stronger than that of the mass term at the corresponding parameter region, we need to take a small lattice constant to approximate the nonlinear Dirac Hamiltonian by the nonlinear QWZ model.Therefore, using the same lattice constant as in the weakly nonlinear regime, one cannot accurately calculate the high-frequency components of the wavefunction in this parameter region.
Exact solutions of the perfectly localized edge modes in the nonlinear QWZ model.
As discussed in the previous section, there can be perfectly localized modes unique to lattice systems.We here derive the exact solutions of such perfectly localized modes in the nonlinear QWZ model (Eq.( 4) in the main text).We consider the localized modes described in the following ansatz, We impose the open boundary condition in the x direction and the twisted boundary condition in the y direction.We obtain the gapless edge modes corresponding to the nonzero nonlinear Chern number as in the case of the fixed norm.We use the parameters u = −1 and κw = 0.1.We note that some bulk modes are not obtained due to the limitation of the numerical method.b, We show an example of an eigenvector of an edge mode, which corresponds to the red circle in panel a.We can check its localization at the edge.previous papers [34,41] that have revealed strong nonlinearity induces the bulk-localization and thus eliminates topological edge modes.

Effects of definitions of normalization conditions.
While we fix the L 2 norm of the nonlinear eigenvector as the normalization condition, one can consider various types of normalization conditions to restrict the solutions of a nonlinear eigenequation.Specifically, some previous studies [40,42] have fixed the maximum absolute value of the components of the nonlinear eigenvector.In this section, we compare the nonlinear band structures under such a normalization condition to those obtained when we fix the L 2 norm.Under both normalization conditions, one can confirm the bulk-boundary correspondence in weakly nonlinear systems and in strongly nonlinear continuum systems.
We numerically calculate the eigenvalues of the nonlinear QWZ model (Eq.( 4) in the main text) under the normalization condition of max x (|ψ 1 (x)| 2 +|ψ 2 (x)| 2 ) = w.Supplementary Figure S8 shows the obtained nonlinear band structure under weak nonlinearity.We confirm that the existence and absence of gapless edge modes correspond to the nonzero and zero nonlinear Chern number for each.Under the weak nonlinearity compared to the band gap in the linear limit (κ = 0 in Eq. ( 4) in the main text), the nonlinear terms can be considered as the perturbation and thus do not alter the topology of the system independently of the normalization conditions.Therefore, weakly nonlinear systems also exhibit the bulk-boundary correspondence under the normalization condition fixing the maximum of the amplitudes.
We also numerically calculate the nonlinear eigenvalues at k y = 0 under the stronger nonlinearity.Supplementary Figure S9a shows the obtained spectra in a 10 × 1 supercell structure under the open boundary condition in the x direction and the twisted boundary condition in the y direction.We use the parameters u = −2.5 and κ = 1.Since the mean of the amplitude is less than the maximum (see Supplementary Fig. S9b), the nonlinear band is still gapped at w = 0.5, where the bulk bands close the gap and the nonlinear Chern number becomes nonzero.Instead, we obtain gapless modes around w = 1.5 and confirm its localization at the edge of the system as shown in Supplementary Fig. S9c.The emergence of the gapless edge modes indicates the nonlinearity-induced topological phase transition, while the bulk-boundary correspondence is broken under the maximum normalization condition due to similar mechanisms to that under the norm-fixed normalization condition.
Exact bulk solutions of the modified nonlinear QWZ model.
While we focus on the nonlinear QWZ model with nonlinear terms −(−1) j κ(|Ψ 1 (x, y)| 2 + |Ψ 2 (x, y)| 2 )Ψ i (x, y), replacing the nonlinear terms with on-site ones −(−1) j |Ψ j | 2 ψ j has no significant effects on the emergence of topological edge modes and the nonlinearity-induced topological phase transitions.We here analytically calculate the bulk solutions of the modified nonlinear QWZ model,

FIG. 2 .
FIG. 2. Phase diagram and dynamics of the nonlinear Qi-Wu-Zhang (QWZ) model.a, The figure shows a schematic of the nonlinear QWZ model.The model has two sublattices (the black circles) at each lattice point encircled by the blue ellipse.The green lines represent the linear couplings.We use the notation Ψi(x, y) to represent the state variable at each sublattice, where (x, y) is the location of the representative point of each lattice point denoted by the red cross.b, We analytically obtain the phase diagram of the prototypical model of nonlinear topological insulators.The horizontal axis represents the parameter of the mass term and the vertical axis corresponds to the strength of the nonlinearity.The color of each separated region represents the difference in the nonlinear Chern number.c, We numerically show the absence of the edge modes in the topologically-trivial parameter region.We simulate the dynamics of the prototypical model of a nonlinear Chern insulator starting from an initial state localized at the left edge.We impose the open boundary condition in the x direction and the periodic boundary condition in the y direction.The figure shows the snapshot at t = 1.The color shows the absolute value of the components of the state vector at each site.The parameters used are u = 3, κ = 0.1, and w = 1, which corresponds to the red square in panel b.d, We numerically show the existence of the long-lived localized state in the weakly nonlinear topological insulator.The figure shows the snapshot of the simulation at t = 1.The sites at the left edge show large amplitudes, which indicates the existence of the edge-localized state.The parameters used are u = −1, κ = 0.1, and w = 1, which corresponds to the blue circle in panel b.

FIG. 3 .
FIG. 3. Nonlinear band structure and edge modes of the weakly nonlinear QWZ model.a, We numerically calculate the nonlinear band structure of the topologically trivial system under the open boundary condition in the x direction and the periodic boundary condition in the y direction.One can confirm the absence of gapless modes.The parameters used are u = 3, κ = 0.1, and w = 1.b, Nonlinear band structure of the topologically nontrivial system is numerically calculated.There are gapless modes that connect the upper and lower bulk bands.The parameters used are u = −1, κ = 0.1, and w = 1.c, The spatial distribution of the gapless mode is presented.The eigenvalue of the localized mode corresponds to the red circle in the band structure in panel b.

FIG. 4 .
FIG. 4. Phase diagram and localized modes of the nonlinear Dirac Hamiltonians.a, We illustrate the phase diagram of the nonlinear Dirac Hamiltonian, which demonstrates the nonlinear bulk-boundary correspondence in the continuum model.The vertical axis represents the amplitude, and the horizontal axes correspond to the parameters of the nonlinear Dirac Hamiltonian.The blue curved surface shows the phase boundary that separates a trivial phase without boundary modes and a topological phase exhibiting localized modes at the left boundary.The red surfaces present the boundaries where the sign of the parameters of the Dirac Hamiltonian changes.The red lines show the phase boundaries in the surfaces of w = 1 and w = 2.In the linear limit (w = 0), the topological phases are separated by the m = 0 axis, while the nonlinearity modifies the boundary of the topological phases.b-d, Each panel shows the representative shape of the localized mode in each of topologically nontrivial parameter regions.b, When m is negative and κ is positive, we obtain a localized mode in the small-amplitude region, which is regarded as a counterpart of a conventional topological edge mode.We set m, κ, and D as m = −0.5, κ = 1, and D = 3. c, When m is positive and κ is negative, the nonlinear Dirac Hamiltonian exhibits the nonlinearity-induced topological phase transition.We obtain an unconventional localized mode if the amplitude is larger than a critical value.In this localized mode, there exist nonvanishing amplitudes even in the limit of x → ∞.Therefore, the localized mode can be unphysical in a truly semi-infinite system because we cannot normalize such a nonvanishing mode, while it still exists and is physically relevant in experimentally realizable finite systems.We set m, κ, and D as m = 0.5, κ = −2, and D = −3.d, When both m and κ are negative, a localized mode appears independently of the amplitude as in linear topological insulators.We set m, κ, and D as m = −0.5, κ = −1, and D = 3.
FIG. 4. Phase diagram and localized modes of the nonlinear Dirac Hamiltonians.a, We illustrate the phase diagram of the nonlinear Dirac Hamiltonian, which demonstrates the nonlinear bulk-boundary correspondence in the continuum model.The vertical axis represents the amplitude, and the horizontal axes correspond to the parameters of the nonlinear Dirac Hamiltonian.The blue curved surface shows the phase boundary that separates a trivial phase without boundary modes and a topological phase exhibiting localized modes at the left boundary.The red surfaces present the boundaries where the sign of the parameters of the Dirac Hamiltonian changes.The red lines show the phase boundaries in the surfaces of w = 1 and w = 2.In the linear limit (w = 0), the topological phases are separated by the m = 0 axis, while the nonlinearity modifies the boundary of the topological phases.b-d, Each panel shows the representative shape of the localized mode in each of topologically nontrivial parameter regions.b, When m is negative and κ is positive, we obtain a localized mode in the small-amplitude region, which is regarded as a counterpart of a conventional topological edge mode.We set m, κ, and D as m = −0.5, κ = 1, and D = 3. c, When m is positive and κ is negative, the nonlinear Dirac Hamiltonian exhibits the nonlinearity-induced topological phase transition.We obtain an unconventional localized mode if the amplitude is larger than a critical value.In this localized mode, there exist nonvanishing amplitudes even in the limit of x → ∞.Therefore, the localized mode can be unphysical in a truly semi-infinite system because we cannot normalize such a nonvanishing mode, while it still exists and is physically relevant in experimentally realizable finite systems.We set m, κ, and D as m = 0.5, κ = −2, and D = −3.d, When both m and κ are negative, a localized mode appears independently of the amplitude as in linear topological insulators.We set m, κ, and D as m = −0.5, κ = −1, and D = 3.

FIG. 5 .
FIG.5.The bulk-boundary correspondence for the lattice model in the continuum limit.We calculate the band gaps of the nonlinear QWZ model at different sizes and lattice constants.We impose the open boundary condition in the x direction and the periodic boundary condition in the y direction.We fix the parameters m = −1, and κw = 1.0,where the nonlinear band gap is closed in the corresponding continuum model.The band gaps versus 10 2 /(L + 1/2) (L is the size of the system) are plotted.The legend shows the lattice constants h.We can confirm that the size of the gap decreases as the system size becomes larger.However, if the lattice constant is so large that the discretized system cannot imitate the behavior of the continuum nonlinear Dirac Hamiltonian, we find a sudden decrease in the size of the band gap.This sudden decrease contradicts the estimation that the band gap is proportional to 1/(L + 1/2) and indicates the existence of irregular gapless modes.Since such gapless modes appear at smaller amplitudes than the critical point of the nonlinear Chern number, the numerical result indicates that the bulkboundary correspondence can be ensured only after taking both the continuum and thermodynamic limits.

FIG. 6 .
FIG. 6. Phase diagram of the quench dynamics of the nonlinear QWZ model.a, The schematics of the experimental protocol of the quench dynamics is shown.First, one excites the edge sites by, e.g., applying lasers to the edge resonators in nonlinear topological photonic insulators.Then, one observes the nonlinear dynamics without external fields and confirms the existence or absence of a long-lived localized state.b, We plot the time evolution of the quench dynamics in a trivial phase.We use the parameters u = −2.5, κ = 1.5, and w = 1, which correspond to the white square in panel d.We confirm the absence of localized edge modes.c, We plot the time evolution of the quench dynamics of nonlinear edge modes.We use the parameters u = −2.5, κ = 0.25, and w = 1, which correspond to the white circle in panel d.We confirm that a localized state remains for a long time, which indicates the existence of edge modes.d and e, We simulate the quench dynamics and plot the amplitude remaining at the edge sites in the long-term limit after a quench in the u-κw plane.The initial configuration is taken as the edge modes in the linear limit at u = −1 (u = 1) in panel d (e).The light-blue lines indicate the parameters where we can obtain exact edge-localized solutions.The gray lines represent the phase boundary derived from the nonlinear Chern number.These lines agree with the right (left) boundary of the topological phase in panel d (e), and thus the quench dynamics shows the shift of the phase boundary by the nonlinearity.

6 FIG
FIG. S1.Bulk nonlinear eigenvectors of the nonlinear Qi-Wu-Zhang (QWZ) model.a, Nonlinear band structures under the open boundary condition in the x direction and the twisted boundary condition in the y direction.The red circle and the green square correspond to the nonlinear eigenvectors in panels b and c, respectively.The parameters used are u = −0.1 and κw = 0.1.b, Wave function of a sine-wave bulk mode under the open boundary condition in the x direction.The red (blue) lines show the real part of the first (second) component of the wave function, and the red (blue) dashed lines show the imaginary part of the first (second) component (these correspondences are the same in panels c, e, and f).c, Wave function of a localized bulk mode under the open boundary condition in the x direction.The localization indicates that there can be eigenstates out of the Bloch-wave ansatz.d, Nonlinear band structures under the periodic boundary condition in the x direction.We consider a 10 × 1 supercell structure and impose the twisted boundary condition in the y direction.The red circle and the green square correspond to the nonlinear eigenvectors in panels e and f, respectively.e, Wave function of a sine-wave bulk mode under the periodic boundary condition in the x direction.f, Wave function of a localized bulk mode under the open boundary condition in the x direction.The localization and periodicity are not altered by the boundary condition.
FIG.S3.Recovery of the bulk-boundary correspondence at the continuum and thermodynamic limit.a, Band gaps at each size and strength of nonlinearity.The squared gaps versus 10 4 /(L+1/2) 2 (L is the size of the system) are plotted.The legend shows the correspondence between the color and the strength of nonlinearity κw.The other parameters used are h = 0.01 and m = −1.We can confirm that the gaps decrease as the size become larger.b, Squared gaps at the thermodynamic limit estimated from the least-square fitting.We fit the squared gaps with the function a/(L + 1/2) 2 + b and plot the obtained b's at different amplitudes.W/L = 1.0 is the transition point, where the band gap is closed and the gapless edge modes appear.
FIG.S6.Nonlinear eigenvalues and the existence of anomalous gapless modes in a stronger nonlinear regime.a, Amplitude dependence of nonlinear eigenvalues at ky = 0. We calculate the nonlinear eigenvalues of the nonlinear QWZ model (Eq.(4) in the main text) in a 10 × 1 supercell structure at u = −2.5 and κ = 10.The transition point is w = 0.5 with w being the squared amplitude, while gapless edge modes appear around w ≃ 0.15 as is pointed out by the red arrow.We note that some eigenvalues disappear at large amplitudes due to the limitation of the numerical technique.b, Amplitude distribution of an anomalous gapless mode.We confirm the localization of the gapless mode at the edge of the sample.The eigenvalue corresponds to that pointed by the red arrow on the left side of panel a. c, Amplitude distribution of a topological gapless mode.We confirm the localization of the topological mode at both edges of the system.One can understand this localized edge mode as the superposition of two edge modes localized at the left and right edges, respectively.The eigenvalue corresponds to that pointed by the green arrow on the right side of panel a.

Ψ 1
FIG.S7.Recovery of the bulk-boundary correspondence in the infinite-amplitude limit.We calculate the nonlinear eigenvalues of the modified model of a nonlinear Chern insulator in Supplementary Eq. (S98).We impose the open boundary condition in the x direction and the twisted boundary condition in the y direction and plot the nonlinear eigenvalues at ky = 0. We use the parameters u = −3, κ = 2, and c = 10.The nonlinear Chern number becomes CNL = −1 around the squared amplitude w ≃ 1 while the band gap is not completely closed at that point.However, one can confirm that the band gap becomes smaller as the amplitude becomes larger, which indicates the existence of genuinely gapless modes and the recovery of the bulk-boundary correspondence in the infinite-amplitude limit.
FIG. S9.Nonlinear eigenvalues and eigenvectors under the maximum-based normalization.a, We calculate the nonlinear eigenvalues at various maximum amplitudes.We impose the open boundary condition in the x direction and the twisted boundary condition in the y direction.The transition point of the nonlinear Chern number should correspond to w = 0.5 with w being the maximum of the squared amplitude, while the nonlinear band has a gap at that parameter.We instead obtain gapless modes at w = 0.5, which indicates the breakdown of the bulk-boundary correspondence under the maximum-based normalization.b, Eigenvector of a bulk mode at w = 0.5, which is pointed by the red arrow on the left side of panel a.The red (blue) curves show the real part of the first (second) component of the wave function, and the red (blue) dashed curves show the imaginary part of the first (second) component.We can confirm that the eigenvector exhibits sine curves, which indicates the absence of edge modes at this amplitude.c, Eigenvector of an edge mode at w = 1.5, which is pointed by the red arrow on the left side of panel a.As in panel b, the red (blue) curves show the real part of the first (second) component of the wave function, and the red (blue) dashed curves show the imaginary part of the first (second) component.We confirm the localization at the edge of the sample.