Topological edge states and gap solitons in the nonlinear Dirac model

TL;DR

This paper investigates topological edge states and gap solitons in nonlinear Dirac models, revealing their interactions and tunability in one- and two-dimensional topological photonic systems through analytical and numerical methods.

Contribution

It provides the first detailed analysis of the interaction between topological edge states and gap solitons in nonlinear Dirac systems, including analytical solutions and numerical modeling.

Findings

Strong nonlinear interaction between edge states and gap solitons.

Analytical solutions demonstrate nonlinear tunability of localized states.

Numerical simulations confirm wave transport behaviors in topological systems.

Abstract

Topological photonics has emerged recently as a novel approach for realizing robust optical circuitry, and the study of nonlinear effects in topological photonics is expected to open the door for tunability of photonic structures with topological properties. Here, we study the topological edge states and topological gap solitons which reside in the same bandgaps described by the nonlinear Dirac model, in both one- and two-dimensions. We reveal strong nonlinear interaction between those dissimilar topological modes manifested in the excitation of the topological edge states by scattered gap solitons. Nonlinear tunability of localized states is explicated with exact analytical solutions for the two-component spinor wave function. Our studies are complemented by spatiotemporal numerical modeling of the wave transport in one- and two-dimensional topological systems.

Equations

$$i{\partial }_t\Psi =({\hat{H}}_D(\delta \mathbf{k})+{\hat{H}}_{NL})\Psi ,$$

$${\hat{H}}_D(\delta \mathbf{k})=\delta k_x{\hat{\sigma }}_x+\delta k_y{\hat{\sigma }}_y+M{\hat{\sigma }}_z,$$

$$i\partial_{t}\psi\!=\left(\!{\begin{array}[]{*{22}{c}}{M-g|\psi_{1}|^{2}}&{k-\partial_{y}}\\ {k+\partial_{y}}&{-M-g|\psi_{2}|^{2}}\end{array}}\!\right)\psi.$$

$$\tau=\gamma\left(t-vy\right)\!,\>\xi=\gamma\left(y-vt\right)\!,\>\gamma=1\bigl{/}\sqrt{1-v^{2}}\bigr{.}.$$

$${\psi }_1^s(\xi )=\sqrt{\frac{{\rho }_s}{2}}\left(e^{-i{\alpha }_s}+\gamma (1+v)e^{i{\alpha }_s}\right)e^{i{\beta }_s},$$

$${\psi }_2^s(\xi )=-i\sqrt{\frac{{\rho }_s}{2}}\left(e^{-i{\alpha }_s}-\gamma (1+v)e^{i{\alpha }_s}\right)e^{i{\beta }_s},$$

$${\alpha }_n\left(\xi \right)={\tan }^{-1}\left[\frac{\Omega -\omega }{\sqrt{{\Omega }^2-{\omega }^2}}\tanh \left(\sqrt{{\Omega }^2-{\omega }^2}\xi \right)\right],$$

$${\alpha }_s(\xi )={\alpha }_n-\delta /2,{\varrho }_s(\xi )=\frac{2(\Omega \cos 2{\alpha }_n-\omega )}{\mathcal{B}+\mathcal{A}{\cos }^{2}2{\alpha }_s},$$

$${\beta }_s(\xi )=v\sqrt{\frac{\mathcal{B}}{\mathcal{A}+\mathcal{B}}}{\tan }^{-1}\left(\sqrt{\frac{\mathcal{B}}{\mathcal{A}+\mathcal{B}}}\tan 2{\alpha }_s\right),$$

$$\left(\!{\begin{array}[]{*{20}{c}}{\psi^{s}_{1}(y)}\\ {\psi^{s}_{2}(y)}\end{array}}\!\right)=\sqrt{2\varrho_{s}(y)}e^{-i\omega(\mathcal{P})t}\left(\!{\begin{array}[]{*{20}{c}}{\cos\alpha_{s}(y)}\\ {-\sin\alpha_{s}(y)}\end{array}}\!\right),$$

$${\hat{H}}_D(\delta k)=-\epsilon {\hat{\sigma }}_x+v_F{\hat{\sigma }}_y\delta k,$$

$${\mathcal{P}}_s(\omega ,v)=\frac{4}{g\gamma \sqrt{1+{\gamma }^2}}\left[\frac{\pi }{2}-{\tan }^{-1}\left(\frac{1}{\sqrt{1+{\gamma }^2}}\frac{{\omega }_s}{\sqrt{{\epsilon }^2-{\omega }^2}}\right)\right]$$

$${\mathcal{E}}_s(\omega ,v)={\omega }_s(1-{\gamma }^2){\mathcal{P}}_s+\frac{4{\gamma }^2\epsilon }{g}{\tan }^{-1}\left(\frac{\sqrt{{\epsilon }^2-{\omega }^2}}{\gamma \epsilon }\right),$$

$${\mathcal{P}}^{(1)}(\omega )=\frac{\sqrt{2}}{g}{\tan }^{-1}\frac{|\omega |}{\sqrt{2({\epsilon }^2-{\omega }^2)}},$$

$${\hat{H}}_{K_{\pm }}=\pm \delta k_x{\hat{\sigma }}_x+\delta k_y{\hat{\sigma }}_y+M{\hat{\sigma }}_z,$$

$$\omega =\Delta {\omega }_{NL}\mp k,$$

Full text

The grand vision of robust waveguiding and routing of light motivates studies of topological photonic systems Lu2016 ; Ozawa2019 . Nonlinearities in photonics and related fields such as acoustics and Bose-Einstein condensates are being pursued to combine topological protection with advanced functionalities including active tunability, genuine non-reciprocity, frequency conversion, and entangled particle generation Zhou2017 ; Kartashov2017 ; Segev2018b ; Chen2018 ; Mittal2018 ; Kruk2018 ; Leykam2018 ; Wang2019 . However, rigorously-defined topological invariants are presently restricted to linearized perturbations to nonlinear steady states Bardyn2016 ; Bleu2016 ; Peano2016 and there is a pressing need to better understand nonlinear topological systems in non-perturbative regimes.

Non-perturbative studies of nonlinear edge states and solitons have largely focused on the dynamics along the edge and been limited to numerical simulations, due to the scarceness of exact solutions to nonlinear problems Segev2013 ; Leykam2016 ; Kartashov2017 ; Solnyshkov2017 . A notable exception are models based on nonlinear coupling, which can be understood semi-analytically in terms of nonlinearity-induced domain walls Gerasimenko2016 ; Hadad2016 ; Hadad2017 ; Hadad2018 ; Zhou2017 ; Leykam2018 ; RingSoliton2018 . Although realized in electronic circuits Hadad2018 ; Wang2019 , this class of models does not admit strictly localized nonlinear modes and is challenging to scale to optics, where local on-site Kerr nonlinearities are more feasible. Recent experiments using coupled optical fibre loops have shown that the latter class of nonlinearities can be used to couple between localized topological edge states and bulk modes Bisianov2018 . So far, however, there were no systematic studies of this phenomenon and the nonlinear dynamics transverse to the edge. This paper is aimed at bridging this gap.

Here we study a generic continuum Dirac model with local Kerr nonlinearity describing dimerized chains of optical resonators or photonic graphene with staggered sublattice potentials, illustrated in Fig. 1. We obtain localized nonlinear modes analytically, showing that bulk gap solitons and edge states have a unified origin as heteroclinic orbits of a nonlinear Hamiltonian system, differing only in their boundary conditions. In the nonlinear system, the mid-gap frequency appears a critical point of a phase transition where the topological edge states emerge. Further numerical simulations reveal that scattering of travelling bulk solitons off the edge can efficiently excite topological edge states and there is an optimal soliton velocity maximizing the energy transfer. Our model is highly relevant to recent theoretical and experimental studies Solnyshkov2017 ; Dobrykh2018 , with our exact solutions more informative and accurate than previous approximate variational approaches. This allows us explicate nonlinear tunability and transport properties of spin-polarized localized nonlinear modes residing in the nontrivial band gap.

We begin with the nonlinear Dirac model describing evolution of a spinor wavefunction $\Psi=[\Psi_{1},\Psi_{2}]$:

[TABLE]

where $\delta{\bm{k}}=(\delta k_{x},\delta k_{y})\equiv-i(\partial_{x},\partial_{y})$ is the momentum and {\hat{H}_{N\!L}}=-g\bigl{[}{|\Psi_{1}|^{2}},{|\Psi_{2}|^{2}}\bigr{]} is a local cubic nonlinearity. We reduce Eq. (1) to a quasi-1D form by considering plane wave-like states along the $x$ axis, $\Psi=\psi e^{ikx}$, treating $k$ as a parameter,

[TABLE]

Conserved quantities are the power $\mathcal{P}=\langle{\bm{\psi}}|{\bm{\psi}}\rangle$ and the energy $\mathcal{E}=\langle{\bm{\psi}}|\left({\hat{H}_{D}+\frac{1}{2}\hat{H}_{N\!L}}\right)|{\bm{\psi}}\rangle$, where the inner products denote integration over $y$.

We seek stationary states with time dependence $\sim\exp(-i\omega t)$ and velocity $v$ by applying the Lorentz transformation,

[TABLE]

Recasting the nonlinear eigenvalue problem as a Hamiltonian system governing the transverse mode profiles (Supplemental Material), we find that all localized modes (edge or bulk) must have vanishing Hamiltonian. This leads to the closed form solution,

[TABLE]

where $\rho_{s}$, $\alpha_{s}$, and $\beta_{s}$ are intensity, spin angle, and phase profiles, respectively, which take the form

[TABLE]

with coefficients $\mathcal{A}=\gamma\left(1+v\right)\textsl{g}$, \mathcal{B}=\gamma\textsl{g}\bigl{/}\left(1-v\right), $\Omega=\sqrt{k^{2}+M_{0}^{2}}$, $\delta=\tan^{-1}{\left(\displaystyle\frac{k}{M_{0}}\right)}$, and we have taken $M=M_{0}\geq 0$. The $v=0$ limit returns spinor components of the chiral soliton Solnyshkov2017 ,

[TABLE]

where the frequency $\omega(\mathcal{P})$ implicitly depends on the total power. Remarkably, this exact solution describes both bulk and edge solitons, with the latter interpreted as part of a stationary bulk soliton bound to the edge distinguished by the boundary conditions.

To illustrate the features of these exact solutions, we consider as an example the Su-Schrieffer-Heeger (SSH) model, describing a 1D dimer chain with alternating weak and strong couplings $\kappa_{1,2}=\kappa\pm\varepsilon/2$ between the nearest neighbors. Near the Brillouin zone edge $k_{y}d=\pi$ the bulk spectrum is captured by the continuum Dirac Hamiltonian

[TABLE]

where $v_{F}=\kappa d$ is Fermi velocity and $d$ is the lattice period. Introducing self-focusing nonlinearity, this corresponds to Eq. (3) with $k=-\varepsilon$ and $M=0$.

The phase plane ($\varrho,\bar{\alpha})$ of the Hamiltonian system determining the mode profiles is shown in Fig. 2(b), where $\bar{\alpha}=\alpha_{s}-\pi/4$. There are three equilibrium points: two saddles $S,S^{\prime}$ $\rho=0$, $\omega=\varepsilon\sin(2\alpha)$, and center $C$ $\alpha=\pi/4$, $\rho={(\varepsilon-\omega)}/{g}$. This phase portrait supports a bright soliton as a heteroclinic trajectory at vanishing Hamiltonian, corresponding to a separatrix between two saddles. In contrast to the nonlinear coupling model studied in Refs. Hadad2016 ; Hadad2017 ; Hadad2018 , the bounded trajectory here is unique. The boundary conditions (i) $\alpha_{s}=\pi/2$ and (ii) $\alpha_{s}=0$ describe linear edge states modified by nonlinearity and nonlinearity-induced edge states respectively. We verified these analytical solutions against stationary solutions found numerically using Newton’s method.

Integrating the intensity of the soliton solutions Eq. (5), we obtain their total power and energy

[TABLE]

where $\omega_{s}=\gamma\omega$ is the soliton frequency in the laboratory frame. In the linear limit $\omega\rightarrow\varepsilon$ the energy of nonlinear interaction vanishes and $\mathcal{E}_{s}\rightarrow\omega_{s}{\mathcal{P}}_{s}$. We similarly obtain the type (i) edge state dispersion

[TABLE]

with asymptotic values ${\mathcal{P}}^{(1)}(\omega\rightarrow 0)\rightarrow 0$ and ${\mathcal{P}}^{(1)}(\omega\rightarrow-\varepsilon)=\pi/g\sqrt{2}$. The power of the type (ii) edge state can then be calculated as the difference ${\mathcal{P}}^{(2)}={\mathcal{P}}_{s}(\omega,v=0)-{\mathcal{P}}^{(1)}$ and has a nonzero threshold ${\mathcal{P}}^{(2)}(\omega\rightarrow 0)\rightarrow\sqrt{2}\pi/g$. The modal dispersion plotted in Fig. 2(a) shows that the type (i) edge states bifurcate from the linear midgap topological edge modes at zero intensity, indicating the absence of any threshold necessary for their existence. They exist in the range $-\varepsilon<\omega(\mathcal{P})\leq 0$ that can be controlled by dimerization strength $\varepsilon$. Notably, $\omega=0$ coincides with a bifurcation of the bulk soliton, which loses its stability and produces a type (ii) nonlinear edge state.

To relate these observations to the topological properties of the SSH model it is useful to study the soliton intensity and spin angle profiles ($\rho,\alpha)$, shown in Figs. 2(c-f). The latter determines the $y$-dependent spin polarization density ${\bm{S}}=\frac{1}{2}{\bm{\psi}}^{\dagger}{\hat{{\bm{\sigma}}}}{\bm{\psi}}$. For solitons at rest this yields ${\bm{S}}_{s}(y)=\rho_{s}(-\sin{2\alpha_{s}},0,\cos{2\alpha_{s}})$. The spin angle $\bar{\alpha}$ changes sign at the soliton core and localization requires $\bar{\alpha}$ to asymptotically approach $\cos 2\bar{\alpha}=\omega/\varepsilon$. The bulk solitons hence have total spin $\langle{\bm{S}}\rangle_{s}=\int{\bm{S}}\,\mathrm{d}y=-{2}/{g}\left(\tan^{-1}\left({\sqrt{1-{\omega^{2}}{\varepsilon^{-2}}}}\right),0,0\right)$, which respects the inversion symmetry of the bulk Hamiltonian. On the other hand, the edge solitons are obtained by asymmetrically cutting the bulk solitons into two pieces, such that they individually have nonzero $S_{z}$ and break the sublattice inversion symmetry. This behaviour is strongly reminiscent of the linear limit, in which topological edge states are created when the edge cuts a dimer in half.

Whether a dimer is cut by the edge is determined by the Zak phase $\gamma$, which is quantized by inversion symmetry. When $\gamma=\pi$ (nontrivial phase) the Wannier centres lie at the cell boundary, corresponding to dimers cut by the edge. The Zak phase can be generalized to the nonlinear case by computing the nonlinear Berry phase of the delocalized modes comprising the bulk band structure BerryPhaseNL2010 . For the nonlinear chain, provided the total dimer intensity $I\equiv|\psi_{1}|^{2}+|\psi_{2}|^{2}$ is lower than the critical value $I<2\varepsilon/g$, the bulk band structure features two dispersive pass-bands, $\omega^{N\!L}_{\pm}(k_{y})=\omega_{\pm}^{L}(k_{y})-gI/2$, where linear dispersion $\omega_{\pm}^{L}(k_{y})$ undergoes a uniform shift towards negative frequencies. The corresponding eigenfunctions ${\bf u}_{\pm}=\sqrt{\frac{I}{2}}\left(e^{-i\varphi(k_{y})},\pm 1\right)$ do not change compared to the linear case, so that $\gamma$ retains its same quantized values, i.e. $\gamma_{\pm}=i\int_{-\pi}^{\pi}dk_{y}\langle\mathbf{u}_{\pm}|\frac{\partial\mathbf{u}_{\pm}}{\partial k_{y}}\rangle=\pi$ in the nontrivial case. While a bifurcation occurs above the critical intensity $I>2\varepsilon/g$, creating an additional flat bulk band at frequency $\omega(|k_{y}|<\sqrt{(Ig/2)^{2}-\varepsilon^{2}})=-gI$, this transition does not appear to affect the bulk solitons or the band gap, since it lies outside their frequency range of $[-\varepsilon,\varepsilon]$. Moreover, though localized modes are characterized by their dispersion $\mathcal{P}(\omega)$, which can be related to the local intensity at the edge as $\omega=-gI(y=0)/2$, this cannot be directly compared with the Bloch wave intensity $I$.

These results suggest that the stationary bulk solitons are more relevant than the bulk nonlinear Bloch waves to the properties of edge states. In particular, edge states emerge at a symmetry-breaking bifurcation of the bulk solitons, which lose their stability. This picture is fully consistent with recent studies of topological laser arrays Malzard2018 ; Cancellieri2019 , which proposed nonlinear topological transitions protected by particle-hole symmetry. By contrast, here the relevant symmetry is inversion symmetry Grusdt2013 .

Eq. (5) can additionally describe moving solitons that are capable of exciting edge states by reflecting off the topologically nontrivial edge, as illustrated in Fig. 3. Upon reflection, the conservation of power and energy defines conservation laws for incident and reflected solitons and the excited edge state. The conversion efficiency tends to grow with decreasing frequency of the pump soliton, i.e. as the bifurcation point is approached, and for a given soliton frequency there is an optimal velocity maximizing the energy transfer. This effect holds beyond the continuum limit in finite discrete lattices, which is illustrated further in the Supplementary Material.

The model (3) can also be implemented in nonlinear photonic graphene with staggered sublattice potential. Near the Dirac points of the bandgap structure, the low-energy Hamiltonian assumes the form

[TABLE]

where $M$ is the detuning of the eigenfrequencies (mass) of two resonators in the bipartite unit cell of graphene. Neglecting inter-valley scattering, we restrict our consideration to single valley $K_{+}$. The valley-Hall domain wall is created between two insulators characterized by distinct values of the mass, $M(y>0)=M_{0}$ and $M(y<0)=-M_{0}$. Such PT-symmetric interface with a real-valued mass governs the relation between the components of the edge state wavefunction that can be interpreted as an effective boundary condition Ni2018 . For propagating waves $\propto e^{-i\omega t+ikx}$ bound to the interface $y=0$, the relation $\psi_{1}(0)=\mp\psi_{2}(0)$ holds for $M_{0}>0$ and $M_{0}<0$, respectively. Near the Dirac points, the valley edge states have the linear dispersion $\omega=\mp k$ traversing the gap between the hyperbolic cones of bulk states $\omega^{2}=k^{2}+M_{0}^{2}$. Introducing nonlinearity, with the use of soliton solution (5), we derive the nonlinear dispersion for the edge states

[TABLE]

where the Dirac cross undergoes a shift $\Delta\omega_{N\!L}=-\displaystyle\frac{g}{2}\varrho_{s}\left(\alpha_{s}=\pm\displaystyle\frac{\pi}{4}\right)$ in the band gap $[-M_{0},M_{0}]$. The valley Chern number formally stays the same as its linear counterpart until the upper band of nonlinear Bloch waves forms a self-crossing loop at above-threshold bulk intensities $I\geq 2M_{0}$ Bomantara2017 . Thereby, nonlinearity grants control over the frequency and transverse structure of the edge state as defined by $\varrho_{s}(\alpha_{s})$. We stress that our solution describes the localization transverse to the edge assuming plane wave-like profiles with fixed $k$ parallel to the edge. For finite wavepackets, higher order terms in $k$ will induce diffraction along the edge, with weak nonlinearity inducing edge solitons Ablowitz2014 .

In contrast to the 1D case, 2D bulk solitons at $M\neq 0$ cannot be obtained analytically. Nevertheless, by analogy with the above results we can consider excitation of valley edge sites by scattering of bulk solitons. Stationary bulk semi-vortex solitons with harmonic time dependence and radial symmetry, whose spinor in the polar coordinate system $(r,\varphi)$ is given by $[\Psi_{1},\Psi_{2}]=[u(r),iv(r)e^{i\varphi}]e^{-i\omega t}$ at $M_{0}>0$ and $[\Psi_{1},\Psi_{2}]=[-v(r)e^{-i\varphi},iu(r)]e^{-i\omega t}$ at $M_{0}<0$, are numerically calculated by the shooting technique after Chebyshev discretization on the radial coordinate is performed Maravero2016 ; RingSoliton2018 . At $M_{0}=1$, solitons are found to be stable at $\omega>0.388$ in agreement with Refs. Maraver2017b ; Sakaguchi2018 . Time dynamics is modeled with a custom numerical code employing a split-step scheme and the Fourier spectral method accomplished by means of the fast Fourier transform. Periodic boundary conditions are applied to the rectangular calculation domain with an equispaced grid. In Fig. 4(a), a topological cavity is created by mass inversion in a circular domain. A clock-wise pulse of edge waves is seen to be excited at a closed contour of the cavity by a bulk soliton. The soliton excites persistent topological current and goes away tilted in the opposite direction. In Fig. 4(b), a soliton initially launched in $y$ direction is trapped in a waveguide made of two parallel topological interfaces supporting counter-propagating valley edge states. The soliton moves in a zigzag between the waveguide boundaries. Upon each reflection from the the walls, it emits a pulse of edge waves.

Given the potential to explore the interplay between nonlinear and topological effects in the photonic settings, where the infinite-contrast tight-binding approximation is often less applicable, the perspective of continuum nonlinear edge models and bifurcation analysis offer valuable insights. Our method can be applied to other classes of topological lattices, with the example of the nonlinear distorted Kagome lattice given in Supplemental Material.

In summary, we have studied topological localized states in the nonlinear Dirac model and demonstrated close connections between edge states and self-trapped nonlinear modes in topological band gaps. Both can be inferred from phase portraits of the same nonlinear mapping. The bulk solitons resemble the Wannier functions in the linear limit, with nonlinear edge states emerging precisely when the bulk solitons are cut by the edge. This bifurcation of edge states is accompanied by the bulk soliton destabilizing. We furthermore demonstrated numerically that mutual transformations between edge and bulk states, forbidden in linear limit, can occur in the nonlinear regime in one- and two-dimensional systems. Our findings could have important implications for further developments of nonlinear topological systems not being limited to photonics but spanning through the fields of metamaterials to the topological physics of cold atoms in optical lattices.