Stability of topologically protected edge states in nonlinear quantum walks: Additional bifurcations unique to Floquet systems

TL;DR

This paper investigates how increasing nonlinearity causes new bifurcations in topologically protected edge states of Floquet quantum walks, revealing unique dynamical behaviors not seen in static systems.

Contribution

It uncovers additional bifurcations in nonlinear Floquet quantum walks where edge states transition from stable to unstable, a phenomenon absent in non-Floquet systems.

Findings

Identification of bifurcations caused by nonlinearity in Floquet systems

Edge states can switch stability with increasing nonlinearity

Unique bifurcation phenomena specific to Floquet dynamics

Abstract

Recently, effects of nonlinearity on topologically nontrivial systems have attracted attention and the stability of topologically protected edge states has been studied for a quantum walk with nonlinear effects, which is akin to time-periodically driven systems (Floquet systems). In the previous work, it has been found that the edge states can be stable attractors or unstable repellers depending on their intrinsic topological property, while the stability is not affected by the strength of nonlinearity. In the present work, we find additional bifurcations at which edge states change from stable attractors to unstable repellers with increasing the strength of nonlinearity in nonlinear quantum walks, for the first time. The new bifurcations are unique to Floquet systems, since we take dynamical properties of Floquet systems into consideration by directly applying the time-evolution operator of the quantum walks to the linear stability analysis. Our results shed new light on nonlinear effects on topological edge states in Floquet systems.

Equations

$$|\psi (t)⟩=\sum _{x,s=L,R}{\psi }_{x,s}(t)|x⟩\otimes |s⟩,$$

$$|\psi (t+1)⟩=U|\psi (t)⟩.$$

$$C(\theta )=\sum _x|x⟩⟨x|\otimes \mathcal{C}[\theta (x)],\mathcal{C}[\theta (x)]=e^{-i\theta (x){\sigma }_2},$$

$$S=\sum_{x}\left(\begin{array}[]{cc}\ket{x-1}\bra{x}&0\\ 0&\ket{x+1}\bra{x}\end{array}\right),$$

$$U_1=C({\theta }_0/2)SC({\theta }_0/2).$$

$$U_2=U_{2b}{U}_{2a},$$

$$U_{2a}=C({\theta }_2/2)SC({\theta }_1/2),U_{2b}=C({\theta }_1/2)SC({\theta }_2/2).$$

$$|{\Phi }_{\epsilon ,\eta }⟩=\sum _{x,s}{\Phi }_{x,s}^{(\epsilon ,\eta )}|x⟩\otimes |s⟩,$$

$$U_{1/2}|{\Phi }_{\epsilon ,\eta }⟩=e^{-i\epsilon }|{\Phi }_{\epsilon ,\eta }⟩,\Gamma |{\Phi }_{\epsilon ,\eta }⟩=\eta |{\Phi }_{\epsilon ,\eta }⟩.$$

$${\Phi }_{x,L}^{(\epsilon ,\eta )}=\eta {\Phi }_{x,R}^{(\epsilon ,\eta )}.$$

$$|{\Phi }_x^{(\epsilon )}|\equiv |{\Phi }_{x,L}^{(\epsilon ,\eta )}|=|{\Phi }_{x,R}^{(\epsilon ,\eta )}|,$$

$$C(\kappa \Theta )=\sum _x|x⟩⟨x|\otimes \mathcal{C}[\kappa \Theta (x,t)]$$

$$\Theta (x,t)=|{\psi }_{x,L}(t){|}^2-|{\psi }_{x,R}(t){|}^2.$$

$$|\psi (t+1)⟩=U_{1}C(\kappa \Theta )|\psi (t)⟩.$$

$$i\frac{\partial }{\partial t}|{\psi }_x(t)⟩=-i{\sigma }_3\frac{\partial }{\partial x}|{\psi }_x(t)⟩+[{\theta }_0(x)+\kappa \Theta (x,t)]{\sigma }_2|{\psi }_x(t)⟩,$$

$$|\psi (t)⟩=|{\Phi }_{\epsilon ,\eta }⟩+|\delta \psi (t)⟩,|\delta \psi (t)⟩=\sum _{x,s}\delta {\psi }_{x,s}(t)|x⟩\otimes |s⟩,$$

$$i\frac{\partial }{\partial t}|\delta {\psi }_x(t)⟩=\Omega |\delta {\psi }_x(t)⟩,\Omega =-i{\sigma }_3\frac{\partial }{\partial x}+{\theta }_0{\sigma }_2-2i\kappa |\Phi {|}^2[\eta \text{mbox}1\text{mbox}{l}_s-{\sigma }_1]$$

$$\omega =-2i\kappa \eta |\Phi {|}^2\pm \sqrt{q^2+{\theta }_0^2-4{\kappa }^2|\Phi {|}^4},$$

$${\psi }_{x,L}(0)={\psi }_{x,R}(0)=N\exp (-x^2/2{\Delta }^2),$$

$$|\delta \psi (t+1)⟩=V_1^{(\epsilon ,\eta )}|\delta \psi (t)⟩,V_1^{(\epsilon ,\eta )}=U_1{D}_{\kappa ,\eta }(|{\Phi }_x^{(\epsilon )}|),$$

$$D_{\kappa ,\eta }(|{\Phi }_x^{(\epsilon )}|)=\sum _x|x⟩⟨x|\otimes {\mathcal{D}}_{\kappa ,\eta }(|{\Phi }_x^{(\epsilon )}|),$$

$$\mathcal{D}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)=\mbox{1}\mbox{l}_{s}-2\kappa|\Phi_{x}^{(\varepsilon)}|^{2}\left(\begin{array}[]{cc}\eta&-1\\ -1&\eta\end{array}\right)=(1-2\kappa\eta|\Phi_{x}^{(\varepsilon)}|^{2})\mbox{1}\mbox{l}_{s}+2\kappa|\Phi_{x}^{(\varepsilon)}|^{2}\sigma_{1}.$$

$$\max (|{\lambda }_1^{(\epsilon ,\eta )}|)\le \underset{x}{\max }(1,|1-4\kappa \eta |{\Phi }_x^{(\epsilon )}{|}^2|),$$

$${\kappa }_c=\frac{1}{2{\max }_x(|{\Phi }_x^{(\epsilon )}{|}^2)},$$

$${\Phi }_{x,L}^{(\epsilon ,\eta )}=\eta {\Phi }_{x,R}^{(\epsilon ,\eta )}=N_1(-1{)}^{\frac{\epsilon }{\pi }x}{e}^{-\gamma |x-\eta m^{\prime }|},$$

$$\gamma =\log (\frac{1+\sin {\theta }_0}{\cos {\theta }_0}),N_1=\sqrt{\frac{1-e^{-2\gamma }}{4}},$$

$$F_{\epsilon ,\eta }(t)=|⟨{\Phi }_{\epsilon ,\eta }|\psi (t)⟩|.$$

$$|\psi (t+1)⟩=U_{2b}C[-\kappa \tilde{\Theta }]|\tilde{\psi }(t)⟩,|\tilde{\psi }(t)⟩=U_{2a}C[\kappa \Theta ]|\psi (t)⟩.$$

$$\tilde{\Theta }(x,t)=|{\tilde{\psi }}_{x,L}(t){|}^2-|{\tilde{\psi }}_{x,R}(t){|}^2.$$

$$\Phi_{x,L}^{(0,+)}=\left\{\begin{array}[]{ll}N_{2}e^{-\gamma_{0}|x-m|}&(x\leq m)\\ N_{2}pe^{-\gamma_{0}|x-(m+2)|}&(x>m)\end{array}\right.$$

Full text

1 Introduction

The study of topological phases of matter has been very active in broad research fields, from original condensed matter physics to photonic, acoustic, and exitonic systems, and so on. In the latter new systems, the research of nonlinear effects has attracted theoretical and experimental interests [1, 2, 3, 4, 5, 6, 7] since nonlinear effects can cause interesting phenomena related to topologically protected edge states, such as emergence of solitons [2, 5], frequency shift of edge modes [4], and non-reciprocity [7]. While the nonlinear effects have been actively studied for topological phases in static systems, those for topological phases in time-dependent systems have been unclear, because of difficulty for treating nonlinear effects in time-dependent systems, even in periodically driven systems.

Floquet systems, time-periodically driven systems, have attracted great deal of attention [8, 9] to explore non-trivial topological phases induced by periodically driving external fields [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], analyze the stability of systems itself or limit cycles [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31], to name a few. Floquet systems with high tunability have been realized in various experimental setups. One example is a discrete time quantum walk (hereafter, quantum walk) [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. Quantum walks are described by time-evolution operators and have been exploited to explore Floquet topological phases [38, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. There are several advantages in quantum walks in comparison to condensed matter systems, such as real space observations of edge states [38] and existence of edge states with long lifetime in the presence of dissipation [46]. In addition to these unique features, nonlinear effects in quantum walks using a feed-forward scheme have been studied [58, 59].

In 2016, it was shown that, in a nonlinear quamtum walk, topologically protected edge states become a stable attractor or an unstable repeller after long time evolution, depending only on their topological properties [60], but not on the strength of nonlinearity. To this end, the linear stability analysis was applied using an effective Hamiltonian in a continuum limit. However, taking the continuum limit results in losing an exact description of dynamical properties of the system, which plays a key role in Floquet systems.

In this work, we study the stability of edge states in time-periodically driven systems with nonlinear effects, by focusing on nonlinear quantum walks. To this end, we apply the linear stability analysis of edge states in nonlinear quantum walks, by directly treating time-evolution operators to take the dynamical properties into account. As a consequence, in addition to obtaining consistent results in Ref. [60], we find additional bifurcations where edge states change from stable attractors to unstable repellers with increasing the strength of nonlinearity. The bifurcations are typical examples of phenomena unique to Floquet systems, and cannot be predicted by the approach of the effective Hamiltonians in the continuum limit. We discuss the origin of this bifurcation at the end of this paper.

This paper is organized as follows. In Sec. 2, we explain linear quantum walks and edge states. The stability analysis done in Ref. [60] is reviewed in Sec. 3. Section 4 is devoted to a linear stability analysis for edge states of single (4.1) and two (4.2) step nonlinear quantum walks, by directly treating time-evolution operators. Summary and discussion are given in Sec. 5.

2 Quantum walks without nonlinear effects

Before dealing with nonlinear quantum walks, we explain standard quantum walks without nonlinearity. We consider quantum walks in which walkers move in one dimensional position space $\ket{x}$ and have two internal states $\ket{L}=(1,0)^{\rm T}$ and $\ket{R}=(0,1)^{\rm T}$, where the superscript T denotes the transpose. Using these bases, a state in a time step $t$ is written as

[TABLE]

where $\psi_{x,s}(t)$ denotes the wave function amplitude. By using a time-evolution operator $U$, the state at time step $t+1$ is described as

[TABLE]

In the standard quantum walks, $U$ is composed of a coin operator $C(\theta)$ and a shift operator $S$, which are defined as

[TABLE]

[TABLE]

respectively, where $\sigma_{2}$ is one of Pauli matrices $\sigma_{i}\ (i=1,2,3)$. The position-dependent rotation angle $\theta(x)$ is a parameter which determines how internal states are mixed, and the shift operator changes the position of wave function amplitudes depending on internal states. In quantum walks, quasienergy $\varepsilon$ is defined as $\mu=e^{-i\varepsilon}$, where $\mu$ is the eigenvalue of $U$. The quasienergy $\varepsilon$ has $2\pi$ periodicity.

We consider two types of quantum walks. The first one is a single-step quantum walk, whose time-evolution operator is given by

[TABLE]

The second one is a two-step quantum walk defined by

[TABLE]

where $U_{2a}$ and $U_{2b}$ are

[TABLE]

Note that, since $U_{2}$ shifts wave function amplitudes on even (odd) sites to even (odd) sites as it includes two shift operators, $U_{2}$ is decomposed to a block matrix structure for even sites and odd sites. As defined in a symmetry time frame [50], time-evolution operators of single-step and two-step quantum walks, $U_{1}$ and $U_{2}$, have time-reversal symmetry $\mathcal{T}U\mathcal{T}^{-1}=U^{-1}$, particle-hole symmetry $\Xi U\Xi^{-1}=U$, and chiral symmetry $\Gamma U\Gamma^{-1}=U^{-1}$. Each symmetry operator is defined as $\mathcal{T}=\mbox{1}\mbox{l}_{x}\otimes\sigma_{1}\mathcal{K},\ \Xi=\mbox{1}\mbox{l}_{x}\otimes\mbox{1}\mbox{l}_{s}\mathcal{K},$ and $\Gamma=\mbox{1}\mbox{l}_{x}\otimes\sigma_{1}$, where $\mbox{1}\mbox{l}_{x}=\sum_{x}\ket{x}\bra{x},\ \mbox{1}\mbox{l}_{s}=\sum_{s=L,R}\ket{s}\bra{s}$, and $\mathcal{K}$ is the complex conjugation operator. Thus, $U_{1}$ and $U_{2}$ are classified into BDI class [50, 51, 53, 60, 61]. It is known that, if we consider the system depicted in Fig. 1, there exist topologically protected edge states with $\varepsilon=0$ and $\varepsilon=\pi$ near $x=\pm m$ where $\theta_{i}(x)$ changes its sign. We describe edge states in single and two-step quantum walks as

[TABLE]

where the labels $\varepsilon=0,\pi$ and $\eta=\pm$ represent the quasienergy and chirality of edge states, respectively. Thereby, the edge states $\ket{\Phi_{\varepsilon,\eta}}$ satisfy

[TABLE]

Since the chiral symmetry operator is $\Gamma=\mbox{1}\mbox{l}_{x}\otimes\sigma_{1}$, wave function amplitudes of edge states satisfy the relation

[TABLE]

By taking Eq. (10) into account, we use a simplified notation

[TABLE]

when we focus on absolute values of wavefunction amplitudes and $\varepsilon$.

3 Nonlinear quantum walks and stability of edge states: previous work

It has been proposed to introduce nonlinear effects in quantum walks, using feed-forward control in optical systems. In the proposed feed-forward scheme [58, 59], the intensity of light is measured during the propagation and the conditions of optical elements corresponding to coin operators through which light passes after the measurement are changed depending on results of the measurement. In Ref. [60], a nonlinear coin operator

[TABLE]

was introduced, where $\Theta(x,t)$ corresponding to the rotation angle of the coin operator is defined as

[TABLE]

Since $\Theta(x,t)$ depends on wave function amplitudes at every time step $t$, the coin operator $C(\kappa\Theta)$ introduces nonlinear effects. The parameter $\kappa$ denotes the strength of nonlinearity. For convenience, we fix the sign of $\kappa$ as $\kappa\geq 0$ in the following. From Eqs. (10), (12), and (13), we can understand that the nonlinear coin operator has no effect on edge states, i.e. $C(\kappa\Theta)\ket{\Phi_{\varepsilon,\eta}}=C(0)\ket{\Phi_{\varepsilon,\eta}}=\ket{\Phi_{\varepsilon,\eta}}$, because of $|\Phi_{x,L}^{(\varepsilon,\eta)}|^{2}-|\Phi_{x,R}^{(\varepsilon,\eta)}|^{2}=0$. Therefore, edge states are stationary states in the nonlinear quantum walks we consider below. The time evolution of a single-step nonlinear quantum walk is described by $U_{1}$ and $C(\kappa\Theta)$,

[TABLE]

Since $\Theta$ introduces nonlinear effects, Eq. (14) is a nonlinear equation. As mentioned above, edge states $\ket{\Phi_{\varepsilon,\eta}}$ are stationary states in this dynamics. Taking the continuum limit in $x$ and $t$, Eq. (14) is transformed into a nonlinear Dirac equation

[TABLE]

where $\ket{\psi_{x}(t)}=[\psi_{x,L}(t),\psi_{x,R}(t)]^{\rm T}$ and only the first order in $\theta_{0}(x)$ and $\kappa\Theta(x,t)$ is taken into account. Using Eq. (15), the stability of edge states is studied in Ref. [60].

Suppose that $\ket{\psi(t)}$ consists of the stationary edge state with quasienergy $\varepsilon$ and chirality $\eta$, $\ket{\Phi_{\varepsilon,\eta}}$, and an infinitesimally weak fluctuating state around the edge state $\ket{\delta\psi(t)}$,

[TABLE]

where $|\delta\psi_{x,s}(t)/\Phi_{x,s}^{(\varepsilon,\eta)}|\ll 1$ is assumed. Substituting Eq. (16) into Eq. (15) and expanding the equation up to the first order in $\delta\psi_{x,s}(t)$, the time-evolution equation for $\ket{\delta\psi(t)}$

[TABLE]

is obtained, where $\ket{\delta\psi_{x}(t)}=[\delta\psi_{x,L}(t),\delta\psi_{x,R}(t)]^{\rm T}$. In Eq. (17), for simplicity, $x$ dependence of $\theta_{0}(x)$ is ignored. In addition, $x$ and $\varepsilon$ dependences of $|\Phi_{x}^{(\varepsilon)}|$ are also ignored, and we write it as $|\Phi|$. From Eq. (17), we can understand that the time evolution of a plane wave state with wave number $q$, $\ket{\delta\psi_{x}(t)}=e^{i(qx-\omega t)}\ket{\delta\psi_{0}(0)}$, is determined by the complex frequency

[TABLE]

which is the eigenvalue of $\Omega$ (see A for details). On one hand, if all of Im($\omega$) are negative, $\ket{\delta\psi_{x}(t)}$ decays with time steps. Then, $\ket{\Phi_{\varepsilon,\eta}}$ remains as a stable state. Therefore, Eq. (18) means that an edge state with chirality $\eta=+$ is always stable when $\kappa>0$. On the other hand, if there is $\omega$ whose imaginary part is positive, then, $\ket{\delta\psi_{x}(t)}$ grows with time steps and overhelms the edge state $\ket{\Phi_{\varepsilon,\eta}}$. In this case, $\ket{\Phi_{\varepsilon,\eta}}$ is unstable. So, when $\kappa>0$, $\ket{\Phi_{\varepsilon,-}}$ are always repellers. From the analysis done in Ref. [60], it is concluded that the strength of nonlinearity $\kappa\ (>0)$ is irrelevant to determining whether an edge state becomes an attractor or a repeller.

It is worth noting that, $\Omega$ in Eq. (17) can be seen as a non-Hermitian Hamiltonian and possesses so-called $\mathcal{PT}$ symmetry. However, $\mathcal{PT}$ symmetry does not have any effect on the stability of edge states, while $\mathcal{PT}$ symmetry breaking occurs as the value of $\kappa$ is varied. See A for details.

4 Nonlinear quantum walks and stability of edge states: present work

In this section, we explore the stability of edge states in single and two-step nonlinear quantum walks, without taking the continuum limit, which differs from the analysis in Ref. [60] explained in Sec. 3. We will demonstrate below that this scheme is essential to find additional bifurcations unique to Floquet systems.

In both single and two-step nonlinear quantum walks, there are cases where soliton-like states appear when an initial state $\ket{\psi(0)}$ is localized at a single site, which we do not focus on in the present work. Since we empirically know that soliton-like states do not appear if $\ket{\psi(0)}$ is a Gaussian wave packet, we employ the following initial state

[TABLE]

where the standard deviation $\Delta$ is a parameter and $N$ is a real normalization constant. As wave function amplitudes of the initial state are all real, $\psi_{x,s}(t)$ are always real during the nonlinear time evolution.

4.1 single-step nonlinear quantum walks

First, we consider the single-step nonlinear quantum walk explained in Sec. 3, the same model in Ref. [60]. Substituting Eq. (16) into Eq. (14) we obtain the time-evolution equation for $\ket{\delta\psi(t)}$,

[TABLE]

up to the first order in $\delta\psi_{x,s}(t)$, where $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ is

[TABLE]

[TABLE]

We remark that, if $\kappa=0$, $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ is an identity operator and $V_{1}^{(\varepsilon,\eta)}=U_{1}$ is a unitary operator. Then, the eigenvalues of $V_{1}^{(\varepsilon,\eta)},\ \lambda_{1}^{(\varepsilon,\eta)}$, are complex numbers on a unit circle in a complex plane. For $\kappa>0$, since $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ and then $V_{1}^{(\varepsilon,\eta)}$ are nonunitary operators, $\lambda_{1}^{(\varepsilon,\eta)}$ are not on the unit circle in general. We note that $\lambda_{1}^{(\varepsilon,\eta)}$ appears as quartets because of symmetries of the non-unitary time-evolution operator $V_{1}^{(\varepsilon,\eta)}$. Due to the bipartite structure of $V_{1}^{(\varepsilon,\eta)}$, $V_{1}^{(\varepsilon,\eta)}$ has sublattice symmetry ${\mathcal{S}}V_{1}^{(\varepsilon,\eta)}{\mathcal{S}}^{-1}=-V_{1}^{(\varepsilon,\eta)}$, where ${\mathcal{S}}$ is a unitary operator. In addition, by following the symmetry classification of non-Hermitian systems in Ref. [62], the non-unitary operator $V_{1}^{(\varepsilon,\eta)}$ possesses AZ† particle-hole symmetry $[\tilde{\Xi}(V_{1}^{(\varepsilon,\eta)})^{*}\tilde{\Xi}^{-1}]=V_{1}^{(\varepsilon,\eta)}$, where $\tilde{\Xi}$ is a unitary operator. These two symmetries guarantee that eigenvalues appear as quartets, $\pm\lambda_{1}^{(\varepsilon,\eta)}$ and $\pm(\lambda_{1}^{(\varepsilon,\eta)})^{*}$.

In the linear stability analysis, eigenstates which have the largest value of $|\lambda_{1}^{(\varepsilon,\eta)}|$ dominate the dynamics of the nonunitary time evolution in Eq. (20). On one hand, if max$(|\lambda_{1}^{(\varepsilon,\eta)}|)\leq 1$, $|\delta\psi_{x,s}(t)|$ decays with time steps or remains small. Then, the edge state $\Ket{\Phi_{\varepsilon,\eta}}$ becomes a stable attractor. On the other hand, if max$(|\lambda_{1}^{(\varepsilon,\eta)}|)>1$, $|\delta\psi_{x,s}(t)|$ grows with time steps and $\ket{\Phi_{\varepsilon,\eta}}$ becomes an unstable repeller. Taking it into account that $U_{1}$ is a unitary operator and $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ is Hermitian, the upper bound of $\max(|\lambda_{1}^{(\varepsilon,\eta)}|)$ is determined by

[TABLE]

where $1$ and $1-4\kappa\eta|\Phi_{x}^{(\varepsilon)}|^{2}$ are eigenvalues of $\mathcal{D}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ (see B for derivations). While Eq. (23) gives only the upper bound of $\max(|\lambda_{1}^{(\varepsilon,\eta)}|)$, as we numerically demonstrate later, this inequality can correctly estimate the criteria of the linear stability analysis, namely, $\max(|\lambda_{1}^{(\varepsilon,\eta)}|)\leq 1$ or $\max(|\lambda_{1}^{(\varepsilon,\eta)}|)>1$, at least for the present model. Thereby, here we summarize the results derived from Eq. (23) under the facts that $\kappa>0$ and $|\Phi_{x}^{(\varepsilon)}|>0$:

In the case of $\eta=-$, since max$(|\lambda_{1}^{(\varepsilon,-)}|)\leq 1+4\kappa\max_{x}(|\Phi_{x}^{(\varepsilon)}|^{2})>1$ is satisfied for any $\kappa$, $\ket{\Phi_{\varepsilon,-}}$ can always be unstable repellers. 2. 2.

In the case of $\eta=+$, $\ket{\Phi_{\varepsilon,+}}$ are inevitably stable attractors when $\kappa$ is smaller than a threshold value $\kappa_{c}$, which guarantees max$(|\lambda_{1}^{(\varepsilon,+)}|)\leq 1$. For $\kappa>\kappa_{c}$, however, an additional bifurcation can occur since max$(|\lambda_{1}^{(\varepsilon,+)}|)$ can be larger than one. 3. 3.

The threshold value $\kappa_{c}$ is given by

[TABLE]

which is obtained from the condition $1-4\kappa\eta\max_{x}(|\Phi_{x}^{(\varepsilon)}|^{2})=-1$.

The first result is consistent with the conclusion in Ref. [60]. However, the second and third results claim that $\ket{\Phi_{\varepsilon,+}}$ can be unstable for $\kappa>\kappa_{c}$ and an additional bifurcation occurs with increasing the strength of nonlinearity $\kappa$.

In order to confirm the validity of Eq. (24) by numerically calculating eigenvalues of $V_{1}^{(\varepsilon,\eta)}$, we derive the analytical solution of $\ket{\Phi_{\varepsilon,\eta}}$. For simplicity, we consider the case in which $\theta_{0}>0$ is satisfied. Under an assumption that left and right boundaries in Fig. 1 locate far away and the overlap of edge states on each boundary is negligible, $\ket{\Phi_{\varepsilon,\eta}}$ becomes

[TABLE]

where an inverse of a localization length $\gamma$ and a normalization constant $N_{1}$ are

[TABLE]

and $m^{\prime}=m$ for $|x|\leq m$ and $m^{\prime}=m+1$ for $|x|\geq m+1$. The derivation of Eqs. (25) and (26) is given in C. Edge states with $\eta=+,\ \ket{\Phi_{\varepsilon,+}}$, are localized at the right boundary, and $\ket{\Phi_{\varepsilon,-}}$ are localized at the left boundary. From Eqs. (24) and (25), we obtain the expected threshold value $\kappa_{c}$. Alternating sign changes on the position space, $(-1)^{\frac{\varepsilon}{\pi}x}$, do not influence $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ and $V_{1}^{(\varepsilon,\eta)}$, which can be understood from Eq. (22).

Thus, in the single-step nonlinear quantum walk, $V_{1}^{(0,\eta)}=V_{1}^{(\pi,\eta)}$ and $\lambda_{1}^{(0,\eta)}=\lambda_{1}^{(\pi,\eta)}$ hold, since $\ket{\Phi_{0,\eta}}$ and $\ket{\Phi_{\pi,\eta}}$ have the same localization length $\gamma^{-1}$. Therefore, we focus only on $V_{1}^{(0,\eta)}$ and $\lambda_{1}^{(0,\eta)}$ hereafter. Because of this equivalence, however, there are always two stable states $\ket{\Phi_{0,\eta}}$ and $\ket{\Phi_{\pi,\eta}}$ when $\max(|\lambda_{1}^{(0,\eta)}|)=\max(|\lambda_{1}^{(\pi,\eta)}|)\leq 1$, and we cannot predict which state is realized after long time evolution from the linear stability analysis we shall explain. Substituting Eq. (25) into Eqs. (20)-(22), we calculate $\lambda_{1}^{(\varepsilon,\eta)}$ by numerical diagonalizations. Figure 2 shows $\kappa$ dependence of max$(|\lambda_{1}^{(0,\eta)}|)$. In the case of $\eta=-$, max$(|\lambda_{1}^{(0,-)}|)>1$ is satisfied for all $\kappa>0$ as shown in Fig. 2 (a). Therefore, $\ket{\Phi_{\varepsilon,-}}$ are always unstable repellers. In the case of $\eta=+$, max$(|\lambda_{1}^{(0,+)}|)$ remains to be $1$ when $\kappa$ is smaller than a threshold value $\kappa_{c}\approx 3.61$ (when $\theta_{0}=\pi/8$), as shown in Fig. 2 (b). In this case, $\ket{\Phi_{\varepsilon,+}}$ are stable attractors. However, for $\kappa>\kappa_{c}$, max$(|\lambda_{1}^{(0,+)}|)$ becomes larger than one, which means that $\ket{\Phi_{\varepsilon,+}}$ also become unstable. These numerical results completely agree with the theoretical predictions from Eq. (23). In Fig. 2 (b), the absolute value of $1-4\kappa\max_{x}(|\Phi_{x}^{(\varepsilon)}|^{2})$, one of the eigenvalues of $D_{\kappa,+}(|\Phi_{x}^{(\varepsilon)}|)$, is also plotted to support this conclusion.

Here, we verify the additional bifurcation predicted from the stability analysis by numerically simulating time evolution of the nonlinear quantum walk in Eq. (14) with the initial state in Eq. (19). To begin with, we employ the same parameters as shown in Fig. 2 where the stability analysis predicts the additional bifurcation at a threshold value $\kappa_{c}\approx 3.61$, from Eq. (24). Below the threshold, the probability distribution (normalized intensity) is accumulated around a boundary near $x=150$, forming a stationary state as shown in Fig. 3 (a-1) and (a-2). Figure 3 (a-3) compares the probability distribution of the stationary state at $T=6\times 10^{5}$ and that of an edge state $|\Phi_{0,+}\rangle$ of the quantum walk without nonlinearity. Since both states are almost identical, we can see that the edge state $|\Phi_{0,+}\rangle$ is a stable attractor for $\kappa<\kappa_{c}$. Above the threshold, the dynamics is similar to that with $\kappa<\kappa_{c}$ up to a certain time step; i.e. the probability distribution is accumulated around the boundary as shown in Fig. 3 (b-1). Figure 3 (b-2) clarifies that the shape of the probability distribution near the boundary abruptly shrinks around $3.8\times 10^{5}$ time steps. Remarkably, the state after the abrupt shrink is not stationary as the state slightly fluctuates. From Fig. 3 (b-3), we can see that the edge state $|\Phi_{0,+}\rangle$ is unstable for $\kappa>\kappa_{c}$, since the final state localizes stronger than the edge state $|\Phi_{0,+}\rangle$. As the probability distribution of the other edge state $|\Phi_{\pi,+}\rangle$ localized near $x=150$ is the same as that of $|\Phi_{0,+}\rangle$, we can also see that $|\Phi_{\pi,+}\rangle$ is unstable form Fig. 3 (b-3).

In order to quantitatively study the stability of edge states, we calculate a fidelity which is defined as

[TABLE]

If the edge state $|\Phi_{\varepsilon,\eta}\rangle$ is stable, the fidelity is close to one after many time steps.

In the case of $\kappa=3.6<\kappa_{c}$, as shown in Fig. 4 (a), the fidelity for the edge state $|\Phi_{0,+}\rangle$, $F_{0,+}(t)$, remains to be almost one after long time steps. In the case of $\kappa=3.9>\kappa_{c}$ as shown in Fig. 4 (b), $F_{0,+}(t)$ and $F_{\pi,+}(t)$ cannot reach one and abruptly decrease around $3.8\times 10^{5}$ time steps, which is consistent with the observation in Fig. 3 (b-2). Fidelities for other edge states with $\eta=-$ are almost zero during the time evolution. These observations clearly validate the prediction from the linear stability analysis. We note that the abrupt shrink of the probability distribution in Fig. 3 (b-2) or the sharp drop of $F_{0,+}(t)$ in Fig. 4 (b) cannot be predicted from the linear stability analysis. This is because the linear stability analysis cannot predict behaviours in the unstable region where fluctuations become dominant, since the analysis assumes weak fluctuations around the edge states, while it can predict stable-unstable transitions.

For completeness, we explain how the additional bifurcation occurs, by showing all eigenvalues of $V_{1}^{(0,+)},\ \lambda_{1}^{(0,+)}$, for various values of $\kappa$ in the complex plane as shown in Fig. 5. Here, the inside (outside) of the unit circle in the complex plane corresponds to a stable (unstable) region of the edge states. First of all, we remark that almost all of $\lambda_{1}^{(0,+)}$ are identical or very close to eigenvalues of $U_{1}$ which are on the unit circle due to unitarity of $U_{1}$, except a quartet of eigenvalues by following reasons. Firstly, the edge states of $U_{1},\ \ket{\Phi_{\varepsilon,\eta}}$ are also eigenstates of $V_{1}^{(\varepsilon,\eta)}$ with the same eigenvalue $e^{-i\varepsilon}$, because $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)\ket{\Phi_{\varepsilon,\eta}}=\ket{\Phi_{\varepsilon,\eta}}$ is satisfied, which can be understood from Eqs. (21) and (22). Secondly, ${\mathcal{D}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)}\simeq\mbox{1}\mbox{l}_{s}$ unless $x$ is near the boundaries at $x=\pm m$, since $|\Phi_{x}^{(\varepsilon)}|^{2}$ takes exponentially small values far away from the boundaries. Then, extended bulk states of $U_{1},\ \ket{\Psi_{b}}$, whose eigenvalues are $\mu_{b}$ satisfy $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)\ket{\Psi_{b}}\simeq\ket{\Psi_{b}}$ and $V_{1}\ket{\Psi_{b}}\simeq\mu_{b}\ket{\Psi_{b}}$, as the bulk states have small wave function amplitudes at $x=\pm m$. The only exception is topologically trivial localized states, i.e. impurity states, which are localized near boundaries but whose quasienergy is neither [math] nor $\pi$. These states take important roles for us to understand the additional bifurcation as we explain in the next paragraph.

We find that the absolute values of these eigenvalues are not equal to one and strongly depend on $\kappa$, the strength of nonlinearity. Increasing $\kappa$, the four eigenvalues, plotted as light blue squares, flow toward inside of the unit circle, and reach on the real axis [Fig. 5 (a)]. Further increasing $\kappa$, two eigenvalues (blue squares) flow toward the origin, while other two eigenvalues (purple triangles) move towards the opposite direction [Fig. 5 (b)]. After two eigenvalues (blue squares) collide and pass through each other at the origin [Fig. 5 (c)], four eigenvalues (purple triangles and blue squares) flow toward outside of the unit circle, being on the real axis [Fig. 5 (c-2),(d)]. Again increasing $\kappa$, the eigenvalues (plotted as red triangles and red squares) go out from the unit circle [Fig. 5 (e)], which makes $\ket{\Phi_{\varepsilon,+}}$ unstable. As we have explained, the motion of $\lambda_{1}^{(0,+)}$ in the complex plane changes $\ket{\Phi_{\varepsilon,+}}$ from stable to unstable as $\kappa$ is increased. This additional bifurcation is related to that the stable region of $\lambda_{1}^{(\varepsilon,\eta)}$ is bounded, which is one of peculiar features of Floquet systems (see also Sec. 5). Furthermore, as mentioned in Sec. 3, the linear stability analysis in the continuum limit does not predict the additional bifurcation. These facts suggest that the bifurcation we have shown is unique to Floquet systems.

Before closing this subsection, we extensively check the validity of the additional bifurcation. Figures 6 (a) and (b) show $\max(|\lambda_{1}^{(0,+)}|)$ and $\max[F_{0,+}(T),F_{\pi,+}(T)]$ at $T=3\times 10^{6}$, respectively, for various values of $\theta_{0}$ and $\kappa$. In Fig. 6 (b), as the linear stability analysis cannot predict which state is realized in the nonlinear dynamics, $\ket{\Phi_{0,+}}$ or $\ket{\Phi_{\pi,+}}$, we show $\max[F_{0,+}(T),F_{\pi,+}(T)]$. We also plot analytically predicted values of $\kappa_{c}$ in Eq. (24). We remark that the predicted values agree well with those obtained by numerical results. We also note that, as shown in Fig. 6, $\kappa_{c}$ is a decreasing function of $\theta_{0}$ in $0<\theta_{0}<\pi/2$. This is because, when $\theta_{0}$ and the gap size for $\varepsilon=0,\pi$ are large, $|\Phi_{x}^{(\varepsilon)}|^{2}$ have a large value at localization positions $\pm m$ and $\pm(m+1)$, since the localization length $\gamma^{-1}$ is small (see B for details).

4.2 two-step nonlinear quantum walks

Next, we consider a two-step nonlinear quantum walk whose time evolution is described as

[TABLE]

Here $\Theta(x,t)$ is defined from wave function amplitudes of $\ket{\psi(t)}$ as shown in Eq. (13). $\tilde{\Theta}(x,t)$ is also defined from $\ket{\tilde{\psi}(t)}$ in a simillar way:

[TABLE]

For the same reason with the single-step nonlinear quantum walk, edge states $\ket{\Phi_{\varepsilon,\eta}}$ are stationary states in this dynamics because of $\Theta(x,t)=\tilde{\Theta}(x,t)=0$ for $\ket{\psi(t)}=\ket{\Phi_{\varepsilon,\eta}}$. We emphasize that we cannot apply the analysis in terms of the effective Hamiltonian in Ref.  [60] to the two-step nonlinear quantum walk. This is because two nonlinear coin operators $C(\kappa\Theta)$ and $C(-\kappa\tilde{\Theta})$ in Eq. (28) cancel out each other up to the first order in $\kappa$ by taking the continuum limit in $x$ and $t$. Therefore, we have to directly analyze the time-evolution operator in order to check the stability of edge states in the two-step nonlinear quantum walk.

Before analyzing the stability of edge states $\ket{\Phi_{\varepsilon,\eta}}$, we explain properties of $\ket{\Phi_{\varepsilon,\eta}}$ for the two-step quantum walk without nonlinearlity in Eq. (6). We assume that two boundaries in Fig.1 are separated enough and edge states localized at different boundaries have no overlap each other. See D for details. Note that, in the two-step quantum walk, edge states have finite wave function amplitudes only at even sites or odd sites because of the decoupling between even and odd sites. While we use $\ket{\Phi_{\varepsilon,\eta}}$ at even sites for the linear stability analysis, our conclusion does not depend on the even-odd parity. Figure 7 shows how the values of $(\nu_{0},\nu_{\pi})$ depend on the coin parameters $\theta_{1}$ and $\theta_{2}$, where $\nu_{0}$ and $\nu_{\pi}$ are topological numbers for zero-energy gap and $\pi$-energy gap, respectively. In the following, we consider that the coin parameters in the inner region shown in Fig. 1 sweep the thick line (green) in regions A and B as shown in Fig. 7 for simplicity.

In both regions A and B, $\nu_{0}$ has the same value and the analytical form of $\ket{\Phi_{0,\eta}}$ is

[TABLE]

[TABLE]

where $\gamma_{0},\ p$, and $N_{2}$ are

[TABLE]

See D for derivations. Irrespective of the parameter regions, $\ket{\Phi_{0,+}}$ and $\ket{\Phi_{0,-}}$ are localized at the right and left boundaries, respectively. On the other hand, for $\varepsilon=\pi$, the value of $\nu_{\pi}$ and localization centers of $\ket{\Phi_{\pi,\eta}}$ depend on which parameter region $\theta_{1}$ and $\theta_{2}$ exist in. In the region A, $\ket{\Phi_{\pi,\eta}}$ is

[TABLE]

[TABLE]

where

[TABLE]

Equations (33) and (34) mean that $\ket{\Phi_{\pi,+}}$ is localized at the left boundary and $\ket{\Phi_{\pi,-}}$ is vice versa, in the region A. In the region B,

[TABLE]

[TABLE]

where

[TABLE]

Localization centers of $\ket{\Phi_{\pi,+}}$ and $\ket{\Phi_{\pi,-}}$ in the region B are the opposite of those in the region A. Note that we can obtain $\ket{\Phi_{\varepsilon,\eta}}$ in odd sites by changing $m$, $m+2$, and $p$ to $m-1$, $m+1$, and $p^{-1}$, respectively, in Eqs. (30)-(37). In regions A and B, the normalization constant $N_{2}$ for $\ket{\Phi_{\pi,\eta}}$ is obtained by substituting $\gamma_{\pi}^{A}$ and $\gamma_{\pi}^{B}$ into $\gamma_{0}$ in Eq. (32), respectively.

Now, we consider the linear stability analysis for edge states of the two-step nonlinear quantum walk. Using the analytical form of $\ket{\Phi_{\varepsilon,\eta}}$, Eqs. (30)-(37), we linearize the time-evolution equation (28). In the same way with the single-step nonlinear quantum walk, we assume the infinitesimally weak fluctuating state in Eq. (16) and substitute it into Eqs. (28). Ignoring higher order terms in $\delta\psi_{x,s}(t)$, we can obtain the time-evolution equation of $\ket{\delta\psi(t)}$,

[TABLE]

where $|\tilde{\Phi}_{x}^{(\varepsilon)}|$ denotes absolute values of wave function amplitudes, in the same way with Eq. (11). Edge states $\ket{\tilde{\Phi}_{\varepsilon,\tilde{\eta}}}$ are defined by

[TABLE]

As explained in E, $\ket{\tilde{\Phi}_{\varepsilon,\tilde{\eta}}}$ is also the eigenstate of the same chiral symmetry operator $\Gamma$ and the following relation is satisfied:

[TABLE]

Here, we emphasize that chiralities $\tilde{\eta}$ depend on $\varepsilon$. In the same way as we explained in Sec. 4.1, the largest value of $|\lambda_{2}^{(\varepsilon,\eta)}|$, where $\lambda_{2}^{(\varepsilon,\eta)}$ are eigenvalues of $V_{2}^{(\varepsilon,\eta)}$, dominates the time evolution in Eq. (39). If $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)>1\ (\leq 1)$ is satisfied, where $\ket{\delta\psi(t)}$ grows (does not grow) with time steps and $\ket{\Phi_{\varepsilon,\eta}}$ is unstable (stable). As derived in B, the upper bound of $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$ is determined by a product of the maximum value of eigenvalues of $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ and that of $D_{-\kappa,\tilde{\eta}}(|\tilde{\Phi}_{x}^{(\varepsilon)}|)$. The eigenvalues of $\mathcal{D}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ and $\mathcal{D}_{-\kappa,\tilde{\eta}}(|\tilde{\Phi}_{x}^{(\varepsilon)}|)$ are

[TABLE]

respectively. Taking these eigenvalues into account, the upper bound of $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$ is expressed as

[TABLE]

Contrary to the single-step nonlinear quantum walk in Sec. 4.1, $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$ depends on the quasienergy $\varepsilon$, since the localization length of $\ket{\Phi_{0,\eta}}$ and $\ket{\Phi_{\pi,\eta}}$ are different.

Although Eq. (43) gives only the upper bound of $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$, numerical calculations that we will show later guarantee that the inequality properly predicts the bifurcation points for the present model. Thereby, we again summarize the results obtained from Eq. (43) here.

The edge states $\ket{\Phi_{0,\eta=\pm}}$ and $\ket{\Phi_{\pi,-}}$ can be unstable for any $\kappa>0$, since the right hand side of Eq. (43) is always larger than one and $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)>1$ can be satisfied. 2. 2.

In the case of $\ket{\Phi_{\pi,\eta=+}}$, the edge state is inevitably stable as long as both $\max_{x}(|\delta_{2}(x)|)$ and $\max_{x}(|\tilde{\delta}_{2}(x)|)$ are smaller than or equal to one, corresponding to $\max(|\lambda_{2}^{(\pi,+)}|)\leq 1$. This is satisfied when $\kappa$ is smaller than a threshold value $\kappa_{c}$. For $\kappa>\kappa_{c}$, $\max(|\lambda_{2}^{(\pi,+)}|)>1$ can be satisfied since $\max_{x}(|\delta_{2}(x)|)$ and/or $\max_{x}(|\tilde{\delta}_{2}(x)|)$ is larger than one. Therefore, a transition from a stable attractor to an unstable repeller occurs. 3. 3.

The threshold is given by

[TABLE]

which comes from the condition that one of $\max_{x}(|\delta_{2}(x)|)$ and $\max_{x}(|\tilde{\delta}_{2}(x)|)$ becomes one.

Now, we numerically confirm the validity of the above results. To this end, we calculate eigenvalues of $V_{2}^{(\varepsilon,\eta)}$ by numerical diagonalizations after substituting the analytical solutions of edge states in Eqs. (30)-(38) into Eq. (39). Figure 8 (a)-(d) show $\kappa$ dependence of the maximum value of $|\lambda_{2}^{(\varepsilon,\eta)}|$ for four kinds of edge states under specific conditions. We clearly observe that $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$ for edge states $\ket{\Phi_{x}^{0,+}}$, $\ket{\Phi_{x}^{0,-}}$, and $\ket{\Phi_{x}^{\pi,-}}$ become larger than one for $\kappa>0$. However, $\max(|\lambda_{2}^{(\pi,+)}|)$ for the edge state $\ket{\Phi_{x}^{\pi,+}}$ remains to be one as long as $\kappa<\kappa_{c}$, and starts to increase with increasing $\kappa$ further. When $\theta_{1}=-\pi/4$ and $\theta_{2}=2\pi/5$, the value of $\kappa_{c}$ is estimated as $\kappa_{c}\approx 1.03$ from Eq. (44). These observations are consistent with the analytically derived predictions from Eq. (43).

We verify the above result of the stability analysis by numerically calculating the time evolution of the nonlinear quantum walk in Eq. (28). The fidelity $F_{\varepsilon,\eta}(T)$ at a time step $T=10^{7}$ is shown in Fig. 8 (e). As expected from $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$ in Fig. 8 (a)-(d), $F_{\pi,+}(T)$ is almost one (much smaller than one) for $\kappa\leq\kappa_{c}\ (\kappa>\kappa_{c})$, while $F_{0,+}(T)$, $F_{0,-}(T)$, and $F_{\pi,-}(T)$ are always much smaller than one.

Here, we explain details for $\kappa$ dependences of $\lambda_{2}^{(\pi,+)}$, eigenvalues of $V_{2}^{(\pi,+)}$, in the complex plane shown in Fig. 9. Eigenvalues near the unit circle (gray circles) are inside or on the unit circle. Increasing $\kappa$, one real eigenvalue inside the unit circle (green triangle) flows toward the origin, while two complex eigenvalues (light blue squares) approach to the real axis [Fig. 9 (a) and (b)]. Further increasing $\kappa$, two eigenvalues (green triangle and red empty triangle) flow toward the negative side and go beyond the unit circle [Fig. 9 (c) and (d)]. This establishes $\kappa$ dependences of $\max(|\lambda_{2}^{(\pi,+)}|)$ in Fig. 8.

As the behavior of $\lambda_{2}^{(\pi,+)}$ in the complex plane is crucial, the bifurcation of $\ket{\Phi_{\pi,+}}$ in the nonlinear two-step quantum walk is one of the typical phenomena in Floquet systems.

Finally, we present a comprehensive result to verify the additional bifurcation in the two-step nonlinear quantum walk. Figure 10 (a) and (b) show $\max(|\lambda_{2}^{(\pi,+)}|)$ and $F_{\pi,+}(T)$ at $T=3\times 10^{6}$, respectively, for various $\theta_{1}$, accordingly $\theta_{2}=-0.6\theta_{1}+\pi/5$, and $\kappa$. Comparing Fig. 10 (a) with (b), it is obvious that the stability analysis based on $\max(|\lambda_{2}^{(\pi,+)}|)$ gives the correct prediction, excepting near $\theta_{1}=\theta_{2}=\pi/8$ where the topological number $\nu_{\pi}$ is not defined as shown in Fig. 7 and the edge states $\ket{\Phi_{\pi,\eta}}$ do not exist. The reason is as follows. By putting $\theta_{1}=\theta_{2}$ into Eqs. (35) and (38), we obtain $\gamma_{\pi}^{A/B}=0$. Therefore, the localization length of the edge state diverges. Even near this point, $\ket{\Phi_{\pi,\eta}}$ have a huge localization length. Then, edge states localizing at left and right boundaries could largely overlap each other, and the assumption that $\ket{\Phi_{\pi,+}}$ and $\ket{\Phi_{\pi,-}}$ are independent is not suitable. Therefore, we can shrink this exceptional region by making the system size larger. Note that, when $\ket{\Phi_{\pi,+}}$ is stable and the initial state $\ket{\psi(0)}$ in Eq. (19) has amplitudes at both even and odd sites, $F_{\pi,+}(T)$ for large $T$ is almost $0.5$, since $\ket{\Phi_{\pi,+}}$ have amplitudes only at even or odd sites. Then, we need to replace $\ket{\Phi_{\pi,+}}$ with $\ket{\Phi_{\pi,+}}/\sqrt{2}$ in order to calculate $\lambda_{2}^{(\pi,+)}$, and this makes $\kappa_{c}$ become $2\kappa_{c}$.

5 Summary and discussion

In this work, we have shown that the stability of topologically protected edge states of the single and two-step nonlinear quantum walks can be predicted by the linear stability analysis, precisely taking the dynamical properties of Floquet systems into consideration. Contrary to the previous work [60], we have analytically found the additional bifurcations depending on the strength of nonlinearity for both nonlinear quantum walks by directly applying time-evolution operators of the quantum walks to the linear stability analysis. Then, we have confirmed the validity of the analysis by numerical calculations.

Finally, we discuss the origin of the additional bifurcation from the viewpoint of eigenvalues. In single and two-step nonlinear quantum walks, as shown in Figs. 5 and 9 respectively, the flow of eigenvalues in the complex plane as a function of the strength of nonlinearity is of importance for the additional bifurcations of edge states. Here, we discuss further on these bifurcations, by pointing out that there is a fundamental difference for the stability analysis between time-independent static systems and time-periodically driven systems, in which linearized equations are described by effective static Hamiltonians and time-evolution operators, respectively. On one hand, in the static systems, eigenvalues of non-Hermitian Hamiltonians such as $\omega$ in Eq. (17) are crucially important for the stability of stationary states. On the other hand, in the Floquet systems, eigenvalues of nonunitary time-evolution operators, such as $\lambda_{1}^{(\varepsilon,\eta)}$ and $\lambda_{2}^{(\varepsilon,\eta)}$ in Eqs. (20) and (39) respectively, determine the stability of stationary states. In the former case, a spectral region which makes a stationary state stable is unbounded as shown in Fig. 11 (a), while the stable region is bounded in the latter case, which is surrounded by the unstable region as shown in Fig. 11 (b). This difference makes stationary states of Floquet systems more fragile. To clearly see this, we consider a simplified situation of the two-step nonlinear quantum walk in Sec. 4.2 by focusing on one eigenvalue, say $\lambda$, as shown by the green triangle in Fig. 9. Since $\lambda$ is real and monotonically decreases with increasing the strength of nonlinearity as shown in Fig. 11 (b), We can easily understand that $\lambda$ can flow into the unstable region in this case. This peculiarity is more highlighted by introducing quasienergy defined from $\lambda$ as $\xi=i\log(\lambda)$. The quasienergy $\xi$ plays a similar role of $\omega$, as the stability of stationary states depends on the sign of Im($\omega$) and Im($\xi$) as shown in Figs. 11 (a) and (c), respectively. Figure 11 (c) shows the trajectory of $\xi$ corresponding to $\lambda$ in Fig. 11 (b). When $\lambda$ passes through the origin, the value of $\xi$ jumps from $(0,-\infty)$ to $(\pi,-\infty)$ and the direction of motion is reversed, and then $\xi$ enters in the unstable region. The discontinuity of $\xi$ originates from the phase singularity of $\lambda$ at the origin, and such a discontinuity never occurs in time-independent systems or systems described by effective static Hamiltonians which are approximately derived. While we need to consider more complicated situations for the single-step nonlinear quantum walk, the boundedness of the stable region is important in both nonlinear quantum walks. Since the boundedness is a consequence of treating time-evolution operators for linear stability analysis, the additional bifurcations of edge states that we have shown in the present work are phenomena unique to Floquet systems. It should be interesting to discuss the same type of bifurcations in which stationary states become unstable due to the boundedness of the stable region in other Floquet systems. To this end, we note that quantum walks are defined by simple time-evolution operators, and this advantage makes it possible to derive the bifurcation points analytically. Therefore, it may be reasonable to employ Floquet systems in which time-evolution operators can be derived without approximation, such as models treated in Refs. [19, 63], while the derivation is difficult in Floquet systems described by Hamiltonians which are smooth functions of time, e.g. electronic systems under the application of light [9, 10, 12].

In the single-step nonlinear quantum walk, the effective Hamiltonian in the continuum limit is non-Hermitian (A), although our analysis is based on time-evolution operators. Recently, systems described by non-Hermitian Hamiltonians with $\mathcal{PT}$ symmetry have been explored enthusiastically [64, 65], where $\mathcal{P}$ and $\mathcal{T}$ represent parity and time-reversal, respectively. As shown in A, the effective Hamiltonian has $\mathcal{PT}$ symmetry. Although $\mathcal{PT}$ symmetry does not influence the stability of edge states in the single-step quantum walk (A), generally, it can have the large effect on the stability of stationary states in nonlinear systems. While large number of $\mathcal{PT}$ symmetric non Hermitian systems are experimentally realized in linear optical systems, it may be also interesting to study the relation between stability of stationary states and $\mathcal{PT}$ symmetry in other nonlinear systems.

An open problem is the stronger localization of probability distributions in the unstable region as shown in Fig.3 in Sec. 4. Since the linear stability analysis is not useful in this region, other methods would lead to better understandings of stronger nonlinear effects in quantum walks.

Acknowledgement

We thank Y. Asano, H. Hirori, R. Okamoto, and K. Yakubo for helpful discussions. This work was supported by KAKENHI (Grants No. JP18J20727, No. JP19H01838, No. JP18H01140, No. JP18K18733, and No. JP19K03646) and a Grant-in-Aid for Scientific Research on Innovative Areas (KAKENHI Grant No. JP15H05855 and No. JP18H04210) from the Japan Society for the Promotion of Science.

Appendix A $\mathcal{PT}$ symmetry of the effective Hamiltonian in the continuum limit of the single-step nonlinear quantum walk

We argue the non-Hermitian Hamiltonian in Eq. (17) from the viewpoint of $\mathcal{PT}$ symmetry. $\mathcal{PT}$ symmetry is a combined symmetry of parity and time-reversal symmetries. When a non-Hermitian operator has $\mathcal{PT}$ symmetry, there are two phases, that is, a $\mathcal{PT}$-symmetry unbroken phase and $\mathcal{PT}$-symmetry broken phase [64, 65]. While all eigenvalues are real in $\mathcal{PT}$-symmetry unbroken phase, in the $\mathcal{PT}$-symmetry broken phase, eigenvalues are partially or fully complex. The transition between the two phases is called $\mathcal{PT}$-symmetry breaking.

Assuming a homogeneous system and applying Fourier transformation to the infinitesimally weak fluctuating state $\ket{\delta\psi_{x}(t)}$, Eq. (17) is expressed by a wave number $q$ as,

[TABLE]

where $\Omega(q)$ corresponds to a non-Hermitian effective Hamiltonian

[TABLE]

Here, we focus only on the second term in Eq. (46), and identify that $\tilde{\Omega}(q)$ has $\mathcal{PT}$ symmetry

[TABLE]

where the $\mathcal{PT}$ symmetry operator is

[TABLE]

Since the eigenvalue of $\tilde{\Omega}(q)$ is

[TABLE]

the eigenvalue of $\Omega(q)$ is given by

[TABLE]

If the imaginary part of the eigenvalue of $\Omega(q)$ is negative (positive), the edge state $\ket{\Phi_{\varepsilon,\eta}}$ is stable (unstable). The first term in Eq. (46), $-2i\kappa\eta|\Phi|^{2}$, makes $\ket{\Phi_{\varepsilon,+}}$ ($\ket{\Phi_{\varepsilon,-}}$) stable (unstable). When $\kappa<|\theta_{0}|/4|\Phi|^{2}$, $\tilde{\omega}$ is entirely real because the $\mathcal{PT}$ symmetry is unbroken. Therefore, the stability of the edge states is determined by the first pure imaginary term in Eq. (50). Increasing $\kappa$, $\tilde{\omega}$ can become pure imaginary with positive and negative signs since the $\mathcal{PT}$ symmetry is broken when $\theta_{0}<2\kappa|\Phi|^{2}$. However, the $\mathcal{PT}$ symmetry breaking does not affect the stability of edge states since the first imaginary term in Eq. (50) always dominates the sign of ${\rm Im}(\omega)$ due to the following relation

[TABLE]

Therefore, the $\mathcal{PT}$ symmetry breaking is irrelevant to the stability of $\ket{\Phi_{\varepsilon,\eta}}$ for the analysis in terms of the effective Hamiltonian.

Appendix B Derivation of the bifurcation points : upper bounds of the maximum

for absolute values of eigenvalues for nonunitary operators

Given a square matrix $A$ which is diagonalizable, its spectral radius $\max(|\lambda_{A}|)$ and spectral norm $\sigma(A)$ always satisfy

[TABLE]

where $\lambda_{A}$ ($\lambda_{A^{\dagger}A}$) is the eigenvalue of $A$ ($A^{\dagger}A$). Here, $A^{\dagger}$ denotes Hermitian conjugation of $A$. We can derive the bifurcation points in single and two-step nonlinear quantum walks, where edge states become unstable, using Eq. (52).

First, we derive the bifurcation points in the single-step nonlinear quantum walk. As mentioned in the main text, the stability of edge states $\ket{\Phi_{\varepsilon,\eta}}$ is determined by the maximum value of $|\lambda_{1}^{(\varepsilon,\eta)}|$, which is equivalent to the spectral radius of $V_{1}^{(\varepsilon,\eta)}=U_{1}D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ in Eq. (52). Since $U_{1}$ is unitary, $V_{1}^{(\varepsilon,\eta)}$ satisfies

[TABLE]

Therefore, from Eq. (52), the upper limit of $\max(|\lambda_{1}^{(\varepsilon,\eta)}|)$ is obtained by diagonalizing $D^{\dagger}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$. Since $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ defined in Eq. (21) is Hermitian, the eigenvalues of $D^{\dagger}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ correspond to $\delta_{1}^{2}$ and $\delta_{2}^{2}(x)$, where $\delta_{1}$ and $\delta_{2}(x)$ are the eigenvalues of $\mathcal{D}_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$, given by

[TABLE]

From Eqs. (52)-(54), the spectral norm $\sigma[V_{1}^{(\varepsilon,\eta)}]$ is determined by

[TABLE]

On one hand, in the case of $\eta=-$, $1+4\kappa|\Phi_{x}^{(\varepsilon)}|^{2}$ becomes larger than $1$ and

[TABLE]

is satisfied for nonzero $\kappa\ (>0)$. Thereby, $\ket{\Phi_{\varepsilon,-}}$ can be always unstable. On the other hand, in the case of $\eta=+$, the stability of $\ket{\Phi_{\varepsilon,+}}$ depends on $\kappa$. When $\max_{x}(|1-4\kappa|\Phi_{x}^{(\varepsilon)}|^{2}|)\leq 1$ is satisfied for small $\kappa$, $\ket{\Phi_{\varepsilon,+}}$ is stable because of

[TABLE]

When $1-4\kappa\max_{x}(|\Phi_{x}^{(\varepsilon)}|^{2})<-1$ is satisfied, $\ket{\Phi_{\varepsilon,+}}$ can be unstable due to

[TABLE]

The threshold for the stable to unstable transition is derived as

[TABLE]

from the condition that $1-4\kappa\max_{x}(|\Phi_{x}^{(\varepsilon)}|^{2})=-1$. When the value of $\max_{x}(|\Phi_{x}^{(\varepsilon)}|^{2})$ is small, large $\kappa$ is needed to make $\ket{\Phi_{\varepsilon,+}}$ unstable. Therefore, $\ket{\Phi_{\varepsilon,+}}$ tends to be stable (unstable) if the localization length $[\log(1/\cos\theta_{0}+\tan\theta_{0})]^{-1}$ is large (small).

Second, we derive the bifurcation points in the two-step nonlinear quantum walk, where the edge states $\ket{\Phi_{\varepsilon,\eta}}$ become unstable. In addition to Eq. (52), we also use an inequality for square matrices $B$ and $C$,

[TABLE]

Substituting $V_{2}^{(\varepsilon,\eta)}=U_{2b}D_{-\kappa,\tilde{\eta}}(|\tilde{\Phi}_{x}^{(\varepsilon)}|)U_{2a}D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ into $A$ in Eq. (52) and using Eq. (60), we obtain

[TABLE]

since the spectral norms of the unitary operators, $U_{2a}$ and $U_{2b}$, are always equal to $1$. The eigenvalues of $D_{\kappa,\eta}(|\Phi_{x}^{(\varepsilon)}|)$ are given in Eq. (54), and those of $D_{-\kappa,\tilde{\eta}}(|\tilde{\Phi}_{x}^{(\varepsilon)}|)$ are given by

[TABLE]

From Eqs. (54), (61), and (62), the upper bound of $\max(|\lambda_{2}^{(\varepsilon,\eta)}|)$ is given by

[TABLE]

From Eq. (63), we can understand that, for $(\varepsilon,\eta)=(0,+),\ (0,-),$ and $(\pi,-)$, $\max_{x}(|\delta_{2}(x)|)$ and/or $\max_{x}(|\tilde{\delta}_{2}(x)|)$ are larger than one for nonzero $\kappa\ (>0)$. Thereby, the edge states $\ket{\Phi_{x}^{(0,\pm)}}$ and $\ket{\Phi_{x}^{(\pi,-)}}$ can be unstable. Regarding the edge state with $(\varepsilon,\eta)=(\pi,+)$, it is inevitably stable as long as $\max_{x}(|\delta_{2}(x)|)\leq 1$ and $\max_{x}(|\tilde{\delta}_{2}(x)|)\leq 1$ are satisfied, which means

[TABLE]

When one of $\max_{x}(|\delta_{2}(x)|)$ and $\max_{x}(|\tilde{\delta}_{2}(x)|)$ is larger than one, the upper bound of $\max(|\lambda_{2}^{(\pi,+)}|)$ becomes

[TABLE]

and $\ket{\Phi_{\pi,+}}$ can be unstable, since the right hand side of Eq. (65) is larger than one. Therefore, the threshold of the stable to unstable transition is given by

[TABLE]

Appendix C Edge states in the single-step quantum walk

In C, we derive Eqs. (25) and (26). Since edge states satisfy Eq. (10), in the inner region where $-m\leq x\leq m$ and $\theta_{0}(x)=\theta_{0}$, the time-evolution equation of $\ket{\Phi_{\varepsilon,\eta}}$ becomes

[TABLE]

[TABLE]

where $\alpha_{\eta}$ and $\beta_{\eta}$ are

[TABLE]

Considering edge states localized near the right boundary, we assume

[TABLE]

where $N_{1}$ is the normalization constant. While $\Phi_{x,s}^{(0,\eta)}$ have the same sign in any position, the sign of $\Phi_{x,s}^{(\pi,\eta)}$ is opposite to that of $\Phi_{x\pm 1,s}^{(\pi,\eta)}$. Substituting Eq. (70) into Eqs. (67) and (68), we obtain

[TABLE]

[TABLE]

for both $\varepsilon=0$ and $\varepsilon=\pi$. Since Eqs. (71) and (72) are quadratic equations of $e^{\gamma}$, we obtain two solutions of $e^{\gamma}$ from each equation. We employ one solution which is obtained from both Eqs. (71) and (72). Then, $e^{\gamma}$ becomes

[TABLE]

From Eqs. (69) and (73), $\eta$ must be $+$ for the solution of the localized state ($\gamma>0$). In the outer region, where $x\geq m+1$ and $\theta_{0}(x)=-\theta_{0}$, we obtain the same solution in the same way, while Eq. (70) is changed to $\Phi_{x,s}^{(\varepsilon,\eta)}=N_{1}(-1)^{\frac{\varepsilon}{\pi}x}e^{-\gamma[x-(m+1)]}$. In the case $\ket{\Phi_{\varepsilon,\eta}}$ is localized near the left boundary, $\eta$ must be $-$, while the localization length $\gamma^{-1}$ is the same. The boundary condition of $\ket{\Phi_{\varepsilon,+}}$ is

[TABLE]

[TABLE]

near the right boundary. Assuming

[TABLE]

eqs. (74) and (75) are satisfied. For edge states with $\eta=-$ and localized near the left boundary, the boundary condition is the same with Eq. (76), changing $m$, $m+1$, and $+$ into $-m$, $-m-1$, and $-$, respectively.

Appendix D Edge states in the two-step quantum walk

In D, we derive Eqs. (30)-(38). The edge states have amplitudes only in even sites or odd sites, since wave funtions at even (odd) sites always shift to even (odd) sites after one time step, in the two-step quantum walk described by $U_{2}$. Although we derive edge states at even sites, edge states in odd sites are obtained in the same way. In the inner region where $-m\leq x\leq m$ and $\theta_{i}(x)=\theta_{i}$ ($i=1,2$), Eq. (9) can be wriiten as

[TABLE]

[TABLE]

where $\alpha_{i,\eta}$ and $\beta_{i,\eta}$ are defined as

[TABLE]

Considering the edge states localized near the left boundary, we assume

[TABLE]

where $N_{2}$ is the normalization constant. Substituting Eq. (80) into Eqs. (77) and (78), we obtain

[TABLE]

[TABLE]

Solving Eqs. (81) and (82) for $e^{2\gamma_{\varepsilon}}$, we obtain four solutions. Among them, we choose one which satisfies both Eqs. (81) and (82). Then, $\gamma_{0}$ and $\gamma_{\pi}$ are

[TABLE]

respectively. In the outer region where $x\leq-m-1$ and $\theta_{i}(x)=-\theta_{i}$, $\gamma_{0}$ and $\gamma_{\pi}$ become the same with Eq. (83), while Eq. (80) is changed to $\Phi_{x,s}^{(\varepsilon,\eta)}=N_{2}(-1)^{\frac{\varepsilon x}{2\pi}}e^{-\gamma_{\varepsilon}(-m-2-x)}$. As mentioned in the main text, we focus only on parameter regions A and B in Fig. 7, for simplicity. From Eqs. (79) and (83), it is known that chirality must be minus for $\varepsilon=0$, as $\gamma_{0}>0$. For $\varepsilon=\pi$, $\gamma_{\pi}>0$ is satisfied when $\eta=+$, in the parameter region A where $0<\theta_{2}<\pi/2$ and $-\theta_{2}<\theta_{1}<\theta_{2}$. In the parameter region B where $0<\theta_{1}<\pi/2$ and $-\theta_{1}<\theta_{2}<\theta_{1}$, $\eta=-$ for $\gamma_{\pi}$ to be positive. Edge states localized near the right boundary have the opposite chirality, while $\gamma_{0}$ and $\gamma_{\pi}$ have the same values with Eq. 83). Near the left boundary, $\Phi_{-m,s}^{(\varepsilon,\eta)}$ satisfies

[TABLE]

Substituting Eqs. (80) and (83) into Eq. (84), the boundary condition of $\Phi_{x,s}^{(0,-)}$ is

[TABLE]

in both parameter regions. For $\ket{\Phi_{\pi,\eta}}$, taking it into account the dependence of $\eta$ and $\gamma_{\pi}$ on parameter regions, boundary conditions become

[TABLE]

in the parameter region A, and

[TABLE]

in the parameter region B, from Eqs. (80), (83), and (84). The boundary conditions near the right boundary are obtained in the same way.

Appendix E Derivation of Eq. (41)

In order to derive Eq. (41), we define two time evolution operators with chiral symmetry $\Gamma U\Gamma^{-1}=U^{-1}$, $U=U^{\prime}$ and $U=U^{\prime\prime}$. In the symmtery time frame [50], they are defined as

[TABLE]

where $U_{a}$ and $U_{b}$ satisfy

[TABLE]

Note that, in the main text, $U_{a}$, $U_{b}$, and $U^{\prime}$ correspond to $U_{2a}$, $U_{2b}$, and $U_{2}$, respectively. The edge states in the system described by $U^{\prime}$ are distinguished by the quasienergy $\varepsilon=0,\pi$

[TABLE]

and chirality $\eta=\pm$,

[TABLE]

Acting $U_{a}$ from left on both sides of Eq. (90) and taking Eq. (88) into account, we obtain

[TABLE]

From Eq. (92), it is seen that edge states with quasienergy $\varepsilon$ in the system described by $U^{\prime\prime}$, $\ket{\Phi_{\varepsilon,\tilde{\eta}}^{\prime\prime}}$, are obtained from $\ket{\Phi^{\prime}_{\varepsilon,\eta}}$,

[TABLE]

ignoring a global phase. Chirality of $\ket{\Phi_{\varepsilon,\tilde{\eta}}^{\prime\prime}},\ \tilde{\eta},$ is determined by $\varepsilon$ and chirality of $\ket{\Phi^{\prime}_{\varepsilon,\eta}}$, $\eta$. Acting the chiral symmetry operator $\Gamma$ from left on both sides of Eq. (93), $\Gamma\ket{\Phi_{\varepsilon,\tilde{\eta}}^{\prime\prime}}$ becomes

[TABLE]

since Eqs. (89) and (91) hold. Using Eqs. (88) and (90), the right hand side of Eq. (94) can be written as

[TABLE]

From Eqs. (93), (94), and (95), $\ket{\Phi_{\varepsilon,\tilde{\eta}}^{\prime\prime}}$ satisfies

[TABLE]

which means that Eq. (41) is satisfied, where $\ket{\Phi^{\prime}_{\varepsilon,\eta}}$ and $\ket{\Phi_{\varepsilon,\tilde{\eta}}^{\prime\prime}}=U_{a}\ket{\Phi^{\prime}_{\varepsilon,\eta}}$ correspond to $\ket{\Phi_{\varepsilon,\eta}}$ and $\ket{\tilde{\Phi}_{\varepsilon,\eta}}$, respectively, in the main text. Therefore, chirality of $\ket{\Phi_{0,\tilde{\eta}}^{\prime\prime}}$ is the same as that of $\ket{\Phi^{\prime}_{0,\eta}}$, $\eta$. On the other hand, $\ket{\Phi_{\pi,\tilde{\eta}}^{\prime\prime}}$ has chirality opposite to that of $\ket{\Phi^{\prime}_{\pi,\eta}}$, $-\eta$.

Reference