Gradient catastrophe of nonlinear photonic valley-Hall edge pulses

We derive nonlinear wave equations describing the propagation of slowly-varying wavepackets formed by topological valley-Hall edge states. We show that edge pulses break up even in the absence of spatial dispersion due to nonlinear self-steepening. Self-steepening leads to the previously-unattended effect of a gradient catastrophe, which develops in a finite time determined by the ratio between the pulse's nonlinear frequency shift and the size of the topological band gap. Taking the weak spatial dispersion into account results then in the formation of stable edge quasi-solitons. Our findings are generic to systems governed by Dirac-like Hamiltonians and validated by numerical modeling of pulse propagation along a valley-Hall domain wall in staggered honeycomb waveguide lattices with Kerr nonlinearity.

More sophisticated effective models such as nonlinear Dirac models explicitly include the nontrivial spinlike degrees of freedom required to create topological band gaps [51][52][53][54][55]. In the bulk, the nonlinear Dirac model supports self-induced domain walls and solitons whose stability and degree of localization are sensitive to the gap size. Infinitely-extended (plane wave-like) nonlinear edge states can be obtained analytically and exhibit similar features. However, the nonlinear dynamics of localized edge pulses within the nonlinear Dirac model framework were not yet considered.
In this paper we study an analytically-solvable nonlinear Dirac model (NDM) describing topological edge pulses, revealing that nonlinear topological edge states' exhibit a self-steepening nonlinearity when the pulse self-frequency modulation becomes comparable to the width of the topological band gap. The selfsteepening nonlinearity leads to the formation of a gradient catastrophe of edge wavepackets within a finite propagation time proportional to the pulse width. Taking the weak spatial dispersion of the topological edge modes into account regularizes the catastrophe and results in the formation of stable edge solitons for sufficiently long pulses. We validate our analysis using numerical simulations of beam propagation in a laser-written valley-Hall waveguide lattice, demonstrating that this effect should be observable even for relatively weak nonlinearities. Our findings suggest valley-Hall photonic crystal waveguides will provide a fertile setting for observing and harnessing higherorder nonlinear optical effects.
As a specific example, the model (1) can be implemented in nonlinear photonic graphene with a staggered sublattice potential as illustrated in Fig. 1(a). A dimerized graphene lattice is composed of singlemode dielectric waveguides with local Kerr nonlinearity. The effective mass M characterises a detuning between propagation constants in the waveguides of two sublattices. The form of Hamiltonian operator (2) assumes normalization of the transverse coordinates arXiv:2102.10569v1 [physics.optics] 21 Feb 2021 x, y and evolution variable in the propagation direction t ∼ z/v D to the Dirac velocity v D = 3κa 0 /2 defined by the lattice parameters, a coupling constant κ and a distance a 0 between two neighboring waveguides [23,55]. This continuum model is valid provided |κ| 2|M | [23].
The valley-Hall domain wall is created between two insulators characterized by different signs of the mass. We take M (y > 0) = M 0 and M (y < 0) = −M 0 and without loss of generality assume M 0 > 0 in the upper half-space. In Ref. [55] we derived exact analytic solutions for the propagating nonlinear valley edge modes bound to the interface y = 0: Based on the close connections between nonlinear edge states and self-trapped Dirac solitons in topological band gaps revealed in Ref. [55], the nonlinear edge mode dispersion can be obtained as with the transverse profile of the edge state determined by the intensity at the interface I 1 = |ψ 1,2 (y = 0)| 2 (see Supplemental Material [56]). In Figs. 1(b,c) we show the plane wave-like profile of the nonlinear edge mode with fixed wavenumber k parallel to the edge.
Using the global parity symmetry with respect to the interface and analytical solutions for the edge states [56], we calculate two characteristics of the edge states via integration in the upper half plane y > 0: the power P and spin S x , and identify a functional relation S x (P) between them: Crucially, this relation is independent of the wavevector k, which allows us to develop a slowly varying envelope approximation to describe the nonlinear dynamics of finite edge wavepackets. Using Eq. (1), it can be shown that the integral characteristics obey the following evolution equation: Next, we assume P(x, t) and S x (x, t) are slowly varying functions of the local frequency and wavenumber, such that Eq. (5) remains valid to a first approximation for smooth x-dependent field envelopes. Plugging Eq. (5) into Eq. (6), and using Eq. (3) assuming weak nonlinearity gI 1 M 0 , we obtain the simple nonlinear wave equation for the longitudinal intensity profile I 1 (x, t): Equation (7) suggests that the evolution of finite wavepackets propagating along the x axis is accompanied by steepening of the trailing wavefront up to the development of a gradient catastrophe. Note, in the linear case g = 0, Eq. (7) shows that the edge wavepacket of any shape does not diffract and propagates along the domain wall with constant group velocity v = −1, being granted with topological robustness. Alternatively, Eq. (7) can be derived using asymptotic methods based on a series expansion of the spinor wavefunction [56] where we have introduced a small parameter µ ∼ gI 1 /M 0 , a hierarchy of time scales: τ n = µ n t, and assumed a smooth dependence of the spinor components on t in the moving coordinate frame (ξ ≡ x + t, y).
To illustrate the key effect of the gradient catastrophe captured by Eq. (7), we model the time dynamics of an edge wavepacket using a custom numerical code, applying a split-step scheme and the fast Fourier transform to solve Eq. (1). Figure 2 shows evolution of the initial distribution set in the form of the edge state across the interface with the Gaussian envelope I 0 Plugging the Gaussian distribution into Eq. (7), we may estimate the pulse breakdown time t * analytically: Thus, pulse breakdown occurs for finite wavepackets when the peak nonlinear frequency shift becomes comparable to the size of the topological band gap. As the pulse propagates its tail becomes increasingly steep, developing a discontinuity (i.e. a shock) in a finite time. The numerical solution of Eq. (1) is fully consistent with our analytical considerations, see Fig. 2.
Weak spatial dispersion effects serve as a possible mechanism regularizing the gradient catastrophe, resulting in the formation of solitons. For honeycomb photonic lattices, dispersion is accounted for by introducing off-diagonal second-order derivatives with the coefficient η = (6κ) −1 into the Dirac model (1): Assuming ηM 0 ∼ µ 2 and developing a perturbation theory with expansion (8), we derive an evolution equation governing the complex-valued amplitude a(ξ, t): which differs from the conventional cubic nonlinear Schrödinger equation by the second higher-order nonlinear term responsible for the phase modulation and self-steepening. This equation enables analysis of both the modulational instability of nonlinear plane-wavelike edge states, and the formation of edge quasisolitons [56].
To verify the validity of the modified NLSE (11) we consider the propagation of a Gaussian pulse in Fig. 3. The conventional NLSE, which lacks the selfsteepening term, only exhibits self-focusing and gradual self-compression of the pulse. On the other hand, the modified Eq. (11) correctly reproduces the growing asymmetry of the edge pulse as it propagates. We show in the Supplemental Material how the selfsteepening leads to the break-up of wide pulses, resulting in the radiation of part of its energy into bulk modes, with the remainder forming a quasi-soliton which continues to propagate along the edge and is capable of traversing sharp bends. We note that even after the initial pulse breakup, self-steepening terms can influence the soliton stability and soliton-soliton interactions [32].
As a possible implementation, we consider a realistic valley-Hall waveguide array of laser-written waveguides with parameters similar to those used in the experimental work Ref. [16]. In this case, the evolution variable t becomes the longitudinal propagation distance z. To simulate the evolution, we solve the paraxial equation numerically in a periodic potential by propagating an optical wavepacket [56]. For realistic laser input powers we observe a notable distortion of the beam and signatures of the catastrophe development at its trailing wavefront [ Fig. 4(a,b)]. agrees with modeling of the corresponding continuum Dirac equations and estimates of the breakdown coordinate z * based on Eq. (9) [Fig. 4(d)].
In conclusion, we have described the gradient catastrophe of the nonlinear edge wavepackets in the spinor-type Dirac equation and the formation of edge solitons at the valley-Hall domain walls. We have derived a higher-order self-steepening nonlinear Schrödinger equation describing these effects. Spatiotemporal numerical modeling confirmed that pulse self-steepening can manifest already in the framework of paraxial optics in weakly nonlinear media, such as topological waveguide lattices, and will likely play a key role in future experiments with topological photonic crystal waveguides. Beyond the specific valley-Hall example we considered, our findings are instructive for other emerging experimental studies of nonlinear dynamic phenomena in topological systems, such as the Chern insulators and their implementations in a variety of physical platforms spanning from metamaterials [4] to optical lattices [43,44] and excitonpolariton condensates [57]. This work was supported by the Australian Research Council (Grant DE190100430), the Russian Foundation for Basic Research (Grant 19-52-12053), the National Research Foundation, Prime Ministers Office, Singapore, the Ministry of Education, Singapore under the Research Centres of Excellence programme, and the Polisimulator project co-financed by Greece and the EU Regional Development Fund. Theoretical analysis of the continuum model was supported by the Russian Science Foundation (Grant No. 20-72-00148). D.A.S. thanks Yuri Kivshar for valuable discussions.