Nonlinear one-way edge-mode interactions for frequency mixing in topological photonic crystals

TL;DR

This paper explores nonlinear frequency mixing in topological photonic crystals, demonstrating robust edge-mode interactions for harmonic generation and slow-light enhancement, opening new avenues for topologically protected nonlinear photonics.

Contribution

It introduces the study of nonlinear wave mixing in topological photonic edge modes, a largely unexplored area, and demonstrates practical frequency conversion processes in such systems.

Findings

Second- and third-harmonic generation via edge modes

Slow-light enhancement of nonlinear interactions

Harmonic generation with backward-propagating edge modes

Abstract

Topological photonics aims to utilize topological photonic bands and corresponding edge modes to implement robust light manipulation, which can be readily achieved in the linear regime of light-matter interaction. Importantly, unlike solid state physics, the common test bed for new ideas in topological physics, topological photonics provide an ideal platform to study wave mixing and other nonlinear interactions. These are well-known topics in classical nonlinear optics but largely unexplored in the context of topological photonics. Here, we investigate nonlinear interactions of one-way edge-modes in frequency mixing processes in topological photonic crystals. We present a detailed analysis of the band topology of two-dimensional photonic crystals with hexagonal symmetry and demonstrate that nonlinear optical processes, such as second- and third-harmonic generation can be conveniently implemented via one-way edge modes of this setup. Moreover, we demonstrate that more exotic phenomena, such as slow-light enhancement of nonlinear interactions and harmonic generation upon interaction of backward-propagating (left-handed) edge modes can also be realized. Our work opens up new avenues towards topology-protected frequency mixing processes in photonics.

Figures

c therein) it has been demonstrated the existence of topological edge modes within different band gaps of a graphene plasmonic crystal. Moreover, in another study [Yeung_NanoL14 ](#bib.bib51) (see Fig. 2a therein) the existence of two Dirac points at the frequencies of $3.4\text{\,}\mathrm{THz}$ and $6.8\text{\,}\mathrm{THz}$ has been demonstrated. These Dirac points can be gapped out to form two topological band gaps (ideal for SHG) under an external magnetic field as shown in [Jin_PRL17 ](#bib.bib52); [Pan_NC17 ](#bib.bib53). In fact, recently, we have studied the topologically protected fou

c therein) it has been demonstrated the existence of topological edge modes within different band gaps of a graphene plasmonic crystal. Moreover, in another study [Yeung_NanoL14 ](#bib.bib51) (see Fig. 2a therein) the existence of two Dirac points at the frequencies of $3.4\text{\,}\mathrm{THz}$ and $6.8\text{\,}\mathrm{THz}$ has been demonstrated. These Dirac points can be gapped out to form two topological band gaps (ideal for SHG) under an external magnetic field as shown in [Jin_PRL17 ](#bib.bib52); [Pan_NC17 ](#bib.bib53). In fact, recently, we have studied the topologically protected fou

Equations

$$\displaystyle\mu=\left(\begin{array}[]{ccc}\mu_{0}&i\mu_{i}&0\\ -i\mu_{i}&\mu_{0}&0\\ 0&0&\mu_{0}\\ \end{array}\right).$$

$$C_n=\frac{1}{2\pi }\iint {\mathcal{F}}_n(\mathbf{k})\cdot d\mathbf{k},$$

$${\mathbf{E}}_1^f(\mathbf{r},{\omega }_0)$$

$${\mathbf{H}}_1^f(\mathbf{r},{\omega }_0)$$

$${\mathbf{E}}_1^s(\mathbf{r},{\Omega }_2)$$

$${\mathbf{H}}_1^s(\mathbf{r},{\Omega }_2)$$

$$\frac{1}{4}{\int }_{-\infty }^{\infty }\left({\mathbf{e}}_{f,s}^{\ast }\times {\mathbf{h}}_{f,s}+{\mathbf{e}}_{f,s}\times {\mathbf{h}}_{f,s}^{\ast }\right)\cdot \hat{\mathbf{x}}dy=P_{f,s},$$

$$P_{f,s}=\frac{W_{f,s}}{a}{v}_g^{f,s}=\frac{W_{f,s}^{\text{el}}+W_{f,s}^{\text{mag}}}{a}{v}_g^{f,s},$$

$$W_{f,s}^{\text{el}}$$

$$W_{f,s}^{\text{mag}}$$

$$w_{f,s}^{\text{eff}}=\frac{{\left({\int }_{-\infty }^{\infty }|{\mathbf{e}}_{f,s}\times {\mathbf{h}}_{f,s}|dy\right)}^2}{{\int }_{-\infty }^{\infty }|{\mathbf{e}}_{f,s}\times {\mathbf{h}}_{f,s}{|}^{2}dy}.$$

$$\nabla \times {\mathbf{E}}_1=i\omega \hat{\mu }{\mathbf{H}}_1,$$

$$\nabla \times {\mathbf{H}}_1=-i\omega \widehat{ϵ}{\mathbf{E}}_1,$$

$${\mathbf{E}}_2^f(\mathbf{r},{\omega }_0)$$

$${\mathbf{H}}_2^f(\mathbf{r},{\omega }_0)$$

$${\mathbf{E}}_2^s(\mathbf{r},{\Omega }_2)$$

$${\mathbf{H}}_2^s(\mathbf{r},{\Omega }_2)$$

$$\nabla \times {\mathbf{E}}_2=i\omega \hat{\mu }{\mathbf{H}}_2,$$

$$\nabla \times {\mathbf{H}}_2=-i\omega \widehat{ϵ}{\mathbf{E}}_2-i\omega {\mathbf{P}}_{NL}^{\omega },$$

$${\mathbf{P}}_{NL}^f(\mathbf{r},{\omega }_0)$$

$${\mathbf{P}}_{NL}^s(\mathbf{r},{\Omega }_2)$$

$${\int }_A\nabla \cdot \mathbf{F}dA=\frac{\partial }{\partial x}{\int }_A\mathbf{F}\cdot \hat{\mathbf{x}}dA+{\oint }_{\partial A}\mathbf{F}\cdot \hat{\mathbf{n}}dl,$$

$${\int }_{-\infty }^{\infty }\nabla \cdot \mathbf{F}dy=i\omega {\int }_{-\infty }^{\infty }{\mathbf{E}}_1^{\ast }\cdot {\mathbf{P}}_{NL}^{\omega }dy.$$

$${\int }_{-\infty }^{\infty }\nabla \cdot \mathbf{F}dy=\frac{2i{\omega }_0{A}_f^{\ast }{A}_s{e}^{i\Delta kx}}{P_f\sqrt{P_s}}{\int }_{-\infty }^{\infty }{\mathbf{e}}_f^{\ast }\cdot {\hat{\chi }}^{(2)}({\omega }_0;-{\omega }_0,{\Omega }_2):{\mathbf{e}}_f^{\ast }{\mathbf{e}}_{s}dy,$$

$$\frac{\partial }{\partial x}$$

$$=$$

$$\frac{d{A}_f(x)}{dx}=i{\gamma }_f^{(2)}(x)A_f^{\ast }(x)A_s(x)e^{i\Delta kx},$$

$${\gamma }_f^{(2)}(x)=\frac{{\omega }_0}{2P_f\sqrt{P_s}}{\int }_{-\infty }^{\infty }{\mathbf{e}}_f^{\ast }\cdot {\hat{\chi }}^{(2)}({\omega }_0;-{\omega }_0,{\Omega }_2):{\mathbf{e}}_f^{\ast }{\mathbf{e}}_{s}dy.$$

$${\int }_{-\infty }^{\infty }\nabla \cdot \mathbf{F}dy=\frac{i{\Omega }_2{A}_f^2{e}^{-i\Delta kx}}{P_f\sqrt{P_s}}{\int }_{-\infty }^{\infty }{\mathbf{e}}_s^{\ast }\cdot {\hat{\chi }}^{(2)}({\Omega }_2;{\omega }_0,{\omega }_0):{\mathbf{e}}_f{\mathbf{e}}_{f}dy,$$

$$\frac{\partial }{\partial x}{\int }_{-\infty }^{\infty }$$

$$=$$

$$\frac{d{A}_s(x)}{dx}=i{\gamma }_s^{(2)}(x)A_f^2(x)e^{-i\Delta kx},$$

$${\gamma }_s^{(2)}(x)=\frac{{\Omega }_2}{4P_f\sqrt{P_s}}{\int }_{-\infty }^{\infty }{\mathbf{e}}_s^{\ast }\cdot {\hat{\chi }}^{(2)}({\Omega }_2;{\omega }_0,{\omega }_0):{\mathbf{e}}_f{\mathbf{e}}_{f}dy.$$

$${\gamma }_f^{(2)}(x)$$

$${\gamma }_s^{(2)}(x)$$

$${\chi }_{\mathrm{eff},f}^{(2)}(x)$$

Full text

I Introduction

One of the most important developments in condensed matter physics in the past decades is the discovery of topological insulating materials Kane_RMP10 ; Qi_RMP11 . These materials feature gapped bulk but gapless edge modes, which propagate unidirectionally along the system edge and are immune to local disorder, thus opening a promising avenue towards robust wave manipulation protected by topology. Inspired by this development, the emerging field of topological photonics aims to extend these topology related ideas to the realm of photonics Lu_NP14 ; Shvets_NP17 ; Chen_PQE17 ; Gong_AOM17 ; Chen_OE18 ; Rider_JAP19 ; Ozawa_RMP19 , which holds great promises for innovative optical devices by exploiting robust, scattering-free light propagation and manipulation. As the concept of energy band exists at the single particle level both in condensed matter physics and photonics, the goal of realizing photonic topological insulators can be readily achieved in the linear regime of light matter interaction. Indeed, topological phenomena of electromagnetic waves in a linear medium can be understood by mapping Maxwell equations to the Schrödinger equation Haldane_PRL08 ; Nittis_AP14 .

Photonics, however, has several features not present in solid-state physics. For example, optical gain and loss can be utilized to implement non-Hermitian photonics based on parity-time symmetry Ge_NP17 . The recently realized topological insulator laser demonstrates the power of this new ingredient and could deepen our understanding of the interplay between non-hermiticity and topology in active optical systems TIL_Science1 ; TIL_Science2 . Another well-known feature is the existence of nonlinearity in many optical materials. In fact, optical nonlinear effects play a key rôle in modern photonic applications, giving rise to a variety of important phenomena, including the formation of solitons, modulation and all-optical switching of optical signals, and frequency conversion for the generation of ultrashort pulses Boyd_Book08 . Thus one expects new physics to emerge when adding nonlinearity to photonic systems with nontrivial topological properties. Indeed, it has been shown that when a photonic topological insulator is embedded in an optical medium with Kerr nonlinearity, lattice edge solitons could arise Lumer_PRL13 ; Chong_PRL16 . The possibility to enhance the conversion efficiency of harmonic generation in the presence of topological edge states has also been studied Qian_OE18 ; Chong_NC19 ; yuri_natnano_19 . Moreover, traveling-wave amplifiers Peano_prx16 , topological sources of quantum light Hafezi_nature18 , nonlinear control Yuri_PRL2018 and mapping Yuri_PRL2019 of photonic topological edge states have also been achieved. Despite these important advances, the feasibility of achieving nonlinear optical mixing of edge states of topological photonic crystals via phase matching, which is one of the most fundamental nonlinear optical processes, has not been explored yet. We would also like to highlight the key differences between our work and previous works Qian_OE18 ; Chong_NC19 ; yuri_natnano_19 ; Yuri_PRL2019 on nonlinear optics pertaining to topological edge states: the works Qian_OE18 ; Chong_NC19 ; yuri_natnano_19 are based on one-dimensional (1D) systems, so that the topological edge states are non-propagating optical modes localized at the edges of the 1D system, whereas the system investigated in Yuri_PRL2019 is 2D, i.e., the same as ours, but there is only one topological band gap at the fundamental frequency, and thus there are no nonlinear optical interactions between topological modes.

In this work, we study nonlinear optical interactions of edge modes in topological photonic crystals (PhCs), as per Fig. 1. In particular, we present a detailed study of the band topology of 2D photonic crystals with hexagonal symmetry by mapping out the Chern-number-graded gap phase diagrams. Interestingly, we find that most gaps of the phase diagrams have exactly one edge state in each gap, thus providing a convenient configuration to study the nonlinear interaction of these modes. To this end, by properly tailoring the edge configuration to achieve phase matching, we show that key nonlinear optical processes, such as second- and third-harmonic generation (SHG, THG) can be readily realized in this setup. Beyond this proof-of-principle demonstration of these nonlinear optical processes, we further show that some more exotic nonlinear optical phenomena can also be observed in these topological PhCs, including slow-light enhanced frequency conversion efficiency and higher-harmonic generation upon interaction of so-called backward-propagating (left-handed) modes. All these novel ideas open up new avenues towards active photonic devices with novel functionalities for photonic applications.

The article is organized as follows. In the next section we present and discuss the linear optical properties of the topological photonic crystal, whereas in Sec. III we describe the nonlinear optical interaction between one-way topological modes and the coupled-mode theory that governs the nonlinear propagation of interacting topological modes. Moreover, in Sec. IV we briefly discuss possible experimental implementations of the ideas presented in this study, whereas in the last section we summarize the main conclusions of this work.

II Linear optical properties of the investigated topological photonic crystal

In this section we describe the geometry and material parameters of the topological PhC investigated in this work, as well as the topological properties of the bulk frequency bands and edge topological modes.

II.1 The system

We begin by describing the system setup. First, 2D PhCs possessing topological frequency gaps around frequencies of $\omega_{0}$, $\Omega_{2}=2\omega_{0}$, and $\Omega_{3}=3\omega_{0}$ are designed in order to study SHG and THG via the corresponding edge modes located inside these gaps. As such, in principle any PhC satisfying this condition could be employed. Nonetheless, it would be beneficial if the first gap is topological since typically, the spectral separation among frequency bands and the gap widths decrease as the frequency increases. In view of this, employing the transverse magnetic (TM) modes of a PhC with hexagonal symmetry lattice, which features Dirac cones at $K$ and $K^{\prime}$ points of the first Brillouin zone (FBZ) between the first and second bands pc_book08 , is a natural choice. More specifically, one expects that for this configuration the first gap becomes topological when gapping the Dirac cones by breaking the time-reversal symmetry. Consequently, we consider triangular PhCs whose unit cell contains only one cylinder with radius, $r$, as depicted in Fig. 2a. Lattice structures with hexagonal symmetry but having more cylinders in each unit cell, like honeycomb and Kagome lattices with two and three cylinders, respectively, could potentially be employed, too. The second step of our design procedure is to include magnetic and nonlinear materials. To guide potential experimental implementations and for the sake of specificity, we consider cylinders with low-permittivity ($\epsilon_{1}$), non-magnetic nonlinear material immersed in a magnetic background material with high-permittivity ($\epsilon_{2}$). Note that the permittivity of the cylinders has to be lower than that of the background to ensure that Dirac cones exist.

II.2 Topological properties of the bulk frequency bands

We now move on to the topological properties of the bulk frequency bands of the proposed non-magnetic ($\mu=\mu_{0}$) PhCs whose unit cell and FBZ are shown in Figs. 2a and 2b. In the following, we use normalized frequency and momentum, $\overline{\omega}=\omega a/2\pi c$ and $\overline{k}=ka/\pi$, respectively, where $c$ is the speed of light and $a$ is the lattice constant. Figure 2c shows the photonic band structure of the PhC with $r=0.4a$, $\epsilon_{1}=3$, and $\epsilon_{2}=18$ from which one can see the Dirac cone between the first and second bands at $K$ and $K^{\prime}$ points. All band structures presented in this work were calculated using COMSOL Multiphysics 5.3 COMSOL , a commercial software package based on the finite-element method, and validated using Synopsis’s BandSOLVE software BandSOLVE .

As known, the Chern number of each band is zero in systems with time-reversal symmetry Wang_PRL08 . A common way to break time-reversal symmetry and generate bands with nonzero Chern number is to use magnetic materials wang_Nature09 ; Chan_PRL11 ; Skirlo_PRL15 ; Bahari_Science17 ; feifei_NC18 , where the permeability tensor of the material under an external magnetic field along the $z$-axis possesses off-diagonal components in the $x-y$ plane, i.e.,

[TABLE]

Here, we set $\mu_{0}=1$ and take $\mu_{i}$ as a parameter to quantify the effect of time-reversal symmetry breaking. Figure 2d shows the photonic band structure for $\mu_{i}=0.8$, where one can see that the Dirac cone is now gapped.

To characterize the topology of the frequency bands, we calculate the Chern number of the $n$th band, defined as Haldane_PRL08 ; Lu_NP14 :

[TABLE]

where $\mathcal{A}_{n}(\mathbf{k})=\expectationvalue{i\nabla_{\mathbf{k}}}{\mathbf{E}_{n}({\mathbf{k}})}$ and $\mathcal{F}_{n}(\mathbf{k})=\nabla_{\mathbf{k}}\times\mathcal{A}_{n}(\mathbf{k})$ are the Berry connection and Berry curvature, respectively, with $\mathbf{E}_{n}({\mathbf{k}})$ being the electric field of the $n$th band mode with momentum, $\mathbf{k}$. The momentum-space integral is performed over the FBZ, whereas the inner product of the Berry connection is defined as $\bra{\mathbf{E}_{\alpha}}\ket{\mathbf{E}_{\beta}}=\iint\epsilon(\mathbf{r})\mathbf{E}_{\alpha}(\mathbf{r})\cdot\mathbf{E}_{\beta}(\mathbf{r})d\mathbf{r}$, with the real-space integral performed over the unit cell. The Chern number is calculated using the algorithm described in Fukui_JPSJ05 (for more details about this algorithm, see Appendix A).

The calculated Chern numbers of the photonic bands are indicated in Fig. 2d on top of each band, and the gap Chern number, defined as the sum of the Chern numbers of the bands below the gap, is also given for each gap. The gap Chern number characterizes the topology of the gap in the sense that its sign determines the propagation direction of the edge modes and its value indicates the number of edge states located inside the gap Lu_NP14 . An interesting feature revealed by Fig. 2d is that the first few gaps have Chern number $C=1$, which means that there is one edge mode in each gap and all propagate in the same direction. Therefore, this configuration provides a convenient platform to study nonlinear optical processes, e.g., SHG and THG.

To understand intuitively what regimes can be achieved with this setup, it is instructive to map out the Chern-number-graded gap phase diagrams, defined as the variation of the gap Chern numbers with the system parameters $r$, $\epsilon_{1}$, $\epsilon_{2}$, and $\mu_{i}$. We show these gap phase diagrams in Fig. 2e for $r=0.4a$, $\epsilon_{1}=3$, $\epsilon_{2}=18$, and $\mu_{i}=0.8$, when one parameter is varied while keeping fixed the others. The results show that most domains of the phase diagrams have $C=1$. Moreover, one can also see gaps with $C=2$, $C=3$ and, importantly, even gaps with negative Chern numbers, $C=-1$ and $C=-2$. As we will demonstrate, this variety of values of the gap Chern numbers leads to particularly rich physics when nonlinear interactions of topological modes are considered.

II.3 Topological properties of the edge modes

Guided by the phase diagrams in Fig. 2e, we choose the suitable parameters to create photonic gaps suitable to study SHG and THG. According to the principle of bulk-edge correspondence in systems with finite size, when the gap has nonzero Chern number, one-way edge modes will emerge in the gap. We present in Fig. 3a the photonic band structure of a PhC strip with 30 unit cells along the $y$-axis and periodic along the $x$-axis. This figure illustrates the emergence of various edge states across bulk photonic gaps in a range of frequencies. For the sake of clarity, we mark in red and blue the edge states formed on the top and bottom edges of the PhC strip, respectively. The field profiles of the edge states at frequencies $\overline{\omega}=0.2$, $\overline{\Omega}_{2}=2\overline{\omega}=0.4$, and $\overline{\Omega}_{3}=3\overline{\omega}=0.6$, presented in Fig. 3b, highlight the key feature of the edge state – exponential decay of the field away from the edge.

A well-known prerequisite for achieving efficient frequency conversion processes is to phase match the interacting waves Boyd_Book08 . In the current context, this requires a method to tune the wave vectors of the edge modes. As far as we know this issue has not been previously discussed, perhaps due to the irrelevance of phase matching in other (linear) physics involving edge modes. We find that the wave vector of edge modes can be readily tuned by simply changing the configuration of the edge termination. Figure 3c shows how the edge-mode band of the top edge changes when varying the location of the edge termination (see the sketch at the top of Fig. 3c). In general, we find that the edge-mode band shifts by about one reciprocal lattice vector, $G=2\pi/a$, as one increases the width of the PhC strip by one unit cell. Note that due to the periodicity of the system along the $x$-axis, one can always shift the wave vector of the edge mode to the region of $[-\pi/a,~{}\pi/a]$ by adding a momentum of $nG$, with $n$ a suitable integer.

III Nonlinear optical interaction between one-way topological modes

In this section we first introduce the coupled-mode theory (CMT) for SHG mediated by the topological edge modes (a similar derivation for THG is given in the Appendix B). We then present proof-of-concept results, both numerical and analytical, on topics, such as SHG and THG upon edge-mode interaction, SHG in the slow-light regime, and SHG via interaction between forward- and backward-propagating edge modes.

III.1 Coupled-mode theory describing second-harmonic generation

The derivation of the CMT for SHG upon nonlinear interaction of topological edge-modes follows the formulation of the CMT governing nonlinear pulse interactions in nonlinear PhC slab wavegudes pmw10jstqe ; lp16prb . Thus, to begin with, we note that our system contains magneto-optic materials, i.e., $\hat{\mu}\neq\hat{\mu}^{T}$, but rather $\hat{\mu}=\hat{\mu}^{\dagger}$, so that we will use the conjugated form of the Lorentz reciprocity theorem LPT_book for the vector field $\mathbf{F}=\mathbf{E}_{1}^{*}\times\mathbf{H}_{2}+\mathbf{E}_{2}\times\mathbf{H}_{1}^{*}$, where $\{\mathbf{E}_{1},\mathbf{H}_{1}\}$ are the electric and magnetic fields of an eigenmode of the topological PhC, whereas $\{\mathbf{E}_{2},\mathbf{H}_{2}\}$ are the electric and magnetic fields under the effect of nonlinear polarization $\mathbf{P}_{NL}$. Here, $\hat{\mu}$ is the magnetic permeability tensor and the symbols “$T$” and “$\dagger$” indicate the transpose and Hermitian transpose operations, respectively. Note also that as we consider 2D PhCs, we denote by $\hat{\mathbf{z}}$ the out-of-plane direction and all the physical quantities only depend on the $x$ and $y$ coordinates.

To start with, we write the eigenmodes of the topological PhC, both at the fundamental ($f$) and second-harmonic ($s$) frequencies, denoted by $\omega_{0}$ and $\Omega_{2}=2\omega_{0}$, respectively, in their Bloch forms, i.e.,

[TABLE]

and

[TABLE]

In these equations, $\{\mathbf{e}_{f},\mathbf{h}_{f}\}$ and $\{\mathbf{e}_{s},\mathbf{h}_{s}\}$ are the lattice-periodic electric and magnetic Bloch fields of the fundamental and second-harmonic waves, respectively, $k_{f,s}$ are the corresponding Bloch wavevectors, assuming that the waves propagate along the $x$-axis, and the vector $\mathbf{r}$ lies in the transverse $(x,y)$-plane. If we choose the normalization constants $P_{f,s}$ such that

[TABLE]

the modes $\{\mathbf{E}_{1}^{f},\mathbf{H}_{1}^{f}\}$ and $\{\mathbf{E}_{1}^{s},\mathbf{H}_{1}^{s}\}$ in Eqs. (6) and (7), respectively, carry $1\text{\,}\mathrm{W}$ per unit length along the longitudinal $z$-axis.

The mode power is related to the mode energy per unit length contained in one unit cell, $W$, and the group velocity, $v_{g}$, via the relation:

[TABLE]

where, for non-dispersive media, the electric and magnetic energies are given by the following formulae:

[TABLE]

One can also define an effective width of the mode, $w^{\textrm{eff}}$, in terms of the Poynting vector of the field lp16prb ,

[TABLE]

From the Maxwell equations, one can readily infer that, in the frequency domain, the optical modes satisfy the following equations:

[TABLE]

where $\omega=\omega_{0}$ and $\omega=\Omega_{2}$ for the fundamental and second-harmonic waves, respectively.

Following the CMT, for the perturbed problem, we write the electric and magnetic fields of the fundamental and second-harmonic waves as:

[TABLE]

and

[TABLE]

where $A_{f}(x)$ and $A_{s}(x)$ are the slowly-varying envelopes of the fundamental and second-harmonic waves, respectively, under the effect of the nonlinear polarization $\mathbf{P}_{NL}$. The power per unit length carried by the fundamental and second-harmonic fields are $|A_{f}(x)|^{2}$ and $|A_{s}(x)|^{2}$, respectively, meaning that these field amplitudes are measured in $\sqrt{\mathrm{W/m}}$.

The Maxwell equations for the perturbed fields are:

[TABLE]

where $\omega=\omega_{0}$ and $\omega=\Omega_{2}$ for the fundamental and second-harmonic waves, respectively. Moreover, starting from the equation for the second-order nonlinear polarization in the time domain, $\mathbf{P}_{NL}(\mathbf{r},t)=\hat{\chi}^{(2)}(\mathbf{r},t):\mathbf{E}(\mathbf{r},t)\mathbf{E}(\mathbf{r},t)$, where $\hat{\chi}^{(2)}$ is the second-order susceptibility tensor, one can easily show that the nonlinear polarization at the fundamental and second-harmonic frequencies can be written as:

[TABLE]

We now use the 2D form of the divergence theorem, which states that for any general function $\bf{F}$,

[TABLE]

where $A$ is an arbitrary cross-section perpendicular to the direction of wave propagation, $\hat{\mathbf{x}}$, and $\hat{\mathbf{n}}$ is the unit vector outwardly normal onto $\partial A$ in the plane of $A$. If we take $A$ to extend to infinity along the $y$-axis, the line integral vanishes for fields that decay exponentially to infinity. Moreover, the l.h.s. of Eq. (19) can be written as:

[TABLE]

In deriving this equation, we used the vector identity $\nabla\cdot(\mathbf{A}\times\mathbf{B})=\mathbf{B}\cdot(\nabla\times\mathbf{A})-\mathbf{A}\cdot(\nabla\times\mathbf{B})$ and the fact that, since $\hat{\epsilon}$ and $\hat{\mu}$ are Hermitian, the identities $\mathbf{E}_{1}^{*}\cdot\hat{\epsilon}\mathbf{E}_{2}=(\hat{\epsilon}\mathbf{E}_{1})^{*}\cdot\mathbf{E}_{2}$ and $\mathbf{H}_{1}^{*}\cdot\hat{\mu}\mathbf{H}_{2}=(\hat{\mu}\mathbf{H}_{1})^{*}\cdot\mathbf{H}_{2}$ hold.

Let us now consider Eq. (20), written for the fundamental frequency:

[TABLE]

where $\Delta k=k_{s}(\Omega_{2})-2k_{f}(\omega_{0})$ is the wavevector mismatch. Moreover, the r.h.s. of Eq. (19) can be cast as:

[TABLE]

Comparing Eqs. (21) and (III.1), we arrive at the equation describing the slowly-varying mode amplitude $A_{f}(x)$:

[TABLE]

where the nonlinear coefficient at the fundamental frequency is

[TABLE]

The equation governing the evolution of the slowly-varying mode amplitude of the second-harmonic, $A_{s}(x)$, is derived in a similar way. Thus, when $\omega=\Omega_{2}$, the l.h.s. of Eq. (19) becomes:

[TABLE]

and the r.h.s. of Eq. (19) can be expressed as:

[TABLE]

Finally, from Eqs. (25) and (III.1), we obtain the governing equation for the slowly-varying amplitude $A_{s}(x)$:

[TABLE]

where the nonlinear coefficient at the second-harmonic is

[TABLE]

In order to better understand the effects of the slow light on the strength of the nonlinear interaction, we use Eq. (9) in conjunction with Eq. (10) to express the nonlinear coefficients as follows:

[TABLE]

where $Z_{0}$ is the vacuum impedance, $n_{g,f/s}=c/v_{g,f/s}$ are the group indices of the two interacting modes, and the effective second-order susceptibilities $\chi_{\mathrm{eff},f/s}^{(2)}(x)$ are defined by the following relations:

[TABLE]

It should be noted that the nonlinear coefficients and nonlinear effective susceptibilities vary with the $x$-coordinate over a characteristic length equal to the lattice constant, $a$, whereas the characteristic length over which the field amplitudes vary is equal to $1/\Delta k$. When the two interacting waves are nearly phase-matched, $a\ll 1/\Delta k$. Therefore, in order to describe the nonlinear mode interactions, it is convenient to introduce the averaged physical quantities

[TABLE]

where $x_{0}$ is arbitrary. Then, by averaging Eqs. (23) and (27), we arrive at the system of CMT governing the nonlinear dynamics of the interacting modes:

[TABLE]

where $\overline{A}_{f}(x)$ and $\overline{A}_{s}(x)$ are the averaged mode amplitudes at the fundamental and second-harmonic frequencies, respectively.

III.2 Description of the full-wave numerical simulations approach

The full-wave dynamics of the nonlinear interaction of topological edge-modes (SHG and THG) were determined numerically using the module “Electromagnetic Waves, Frequency Domain” of COMSOL Multiphysics. Thus, to simulate the nonlinear frequency mixing processes in COMSOL, we defined two “Electromagnetic Waves, Frequency Domain” models, one for the fundamental frequency $\omega_{0}$ and one for the second (third) harmonic frequency $\Omega_{2}$ ($\Omega_{3}$). The two models are coupled using a “Polarization” feature added to each of the models. We assumed that for both the SHG and THG cases the nonlinear susceptibilities are diagonal tensors, the diagonal elements being $\chi_{2}$ and $\chi_{3}$, respectively.

For the study of the SHG (where we consider TM-polarized modes), the nonlinear polarizations at the FF and SH are:

[TABLE]

For the study of THG, the corresponding nonlinear polarizations are:

[TABLE]

III.3 Second-harmonic generation and third-harmonic generation upon edge-mode interaction

We now investigate nonlinear frequency conversion processes via the edge modes indicated in Fig. 3a using full-wave numerical simulations of Maxwell equations (for simulation details, see Appendix A), with the results being summarized in Fig. 4. In the following, we mainly focus on the discussion of SHG, as the results of THG can be understood similarly.

We consider cylinders made of homogeneous and isotropic nonlinear material with a scalar nonlinear second-order susceptibility of $\chi^{(2)}=${10}^{-21}\text{,}\mathrm{C}\text{,}{\mathrm{V}}^{-2}$$ (typical value of $\chi^{(2)}$ varies from ${10}^{-24}\text{\,}\mathrm{C}\text{\,}{\mathrm{V}}^{-2}{10}^{-21}\text{\,}\mathrm{C}\text{\,}{\mathrm{V}}^{-2}$ Saleh_Book07 ). The pump electric field $\mathbf{E}_{1}$ is induced by an external source whereas $\mathbf{E}_{2}$ by the nonlinear polarization at the SH, generated by $\mathbf{E}_{1}$. The amplitude of $\mathbf{E}_{1}$ is chosen such that the undepleted pump approximation holds, i.e. the amplitude of $\mathbf{E}_{1}$ is much larger than that of $\mathbf{E}_{2}$ and thus $\mathbf{E}_{1}$ is roughly constant during the frequency conversion process. Note, however, that our analysis remains valid when this condition is not fulfilled, too, our choice being chiefly guided by a more facile comparison between numerical and theoretical results, which is possible in this propagation regime.

From Figs. 4a and 4b, one can observe that the field profiles of $\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ are indeed the same as the profiles of the edge modes shown in Fig. 3b, indicating that the two edge modes are indeed nonlinearly interacting via the SHG – a key result of our work. The physics of this nonlinear process can be accurately captured by the CMT (see Appendix C). In particular, the period of spatial oscillations of the SH field $\mathbf{E}_{2}$ is determined by the wave-vector mismatch $\Delta k=k_{SH}-2k_{FF}$ ($\Delta k=k_{TH}-3k_{FF}$ for THG). As a result, we can straightforwardly compare the numerically extracted oscillation period of $\mathbf{E}_{2}$ with the theoretical prediction of $\Lambda_{2}=2\pi/\Delta\overline{k}$, thus confirming that the key physics of nonlinear frequency conversion processes is validated by our simulations.

We further validate these conclusions using a much larger simulation domain, with the corresponding results being presented in Figs. 4d–4f. In Fig. 4d, where $\overline{\omega}_{0}=0.2$, $\overline{\Omega}_{2}=0.4$, and $\Delta\overline{k}=0.054$, we have $\Lambda_{2}=37$. The agreement between the predictions the CMT and direct numerical simulations, of both the period and amplitude of SH power oscillation along the propagation distance, is excellent. We calculate $\Delta\overline{k}$ for all the frequencies of the interacting edge modes from Fig. 4e and present the theoretically calculated and numerically extracted oscillation period $\Lambda_{2}$ in Fig. 4f. An excellent agreement between the two results can be observed, which confirms the key physics of SHG, namely phase matching is indeed at work in our photonic system and SHG purely via nonlinear interaction of edge modes occurs in our setup. We also confirm the edge mode mediated THG as shown in Figs. 4g–4i, where the discrepancy in Fig. 4i between the numerical and theoretical results is due to inherent limitations of numerical simulations at very small $\Delta k$.

While the SHG and THG of one-way edge modes are governed by the usual mechanism of topological protection from the chiral nature of the edge modes, where the modes can bypass structural defects, such as the sharp bends shown in Figs. 4a-c, without undergoing backscattering, the effect of structural defects on the coherence properties of nonlinear processes is an important but less understood problem. In order to answer this key question, and taking SHG as an example, we introduce a structural defect at a certain location along the edge where the fundamental and second-harmonic modes co-propagate and present the corresponding simulation results in Fig. 4j. Thus, it can be seen that, as expected, both the fundamental and second-harmonic modes bypass the defect without experiencing backscattering. More importantly, by comparing the results in Figs. 4j and 4d, one can infer that the coherence length of the SHG, namely the oscillation period of the amplitude of the interacting modes, is not altered by the interaction with the structural defect. In other words, the coherence of the nonlinear interaction is preserved in the presence of defects, meaning that the phase-matched nature of the nonlinear mode interaction process is immediately regenerated after the mode interaction with each defect. As the coherence length crucially determines the conversion efficiency of the nonlinear processes, topological protection of the coherence length in nonlinear frequency mixing processes demonstrates a new area where topology can boost the performance of photonic devices based on nonlinear optical processes.

Another relevant question, which we do not intend to fully answer here, is how other types of perturbations affect the phase-matching condition of the nonlinear optical interaction of topological modes. From a practical point of view, the most relevant such perturbation is structural disorder, and in this context we consider two cases, namely weak and strong disorder. The case of weak disorder can be analyzed using perturbative methods. Thus, let us assume that in the absence of disorder the frequency-dependent wave-vector mismatch, $\Delta k_{i}(\omega)$, vanishes at a certain frequency $\omega_{0}$, that is $\Delta k_{i}(\omega_{0})=0$. Then, adding weak disorder to the unperturbed system will change the modal dispersion of the interacting modes, such that the wave-vector mismatch varies, say, by a small quantity $\delta k(\omega)$. Furthermore, since one expects that $\delta k(\omega)$ has constant sign around the frequency $\omega_{0}$, it can be seen that there is a frequency $\omega_{0}+\delta\omega_{0}$ at which the wave-vector mismatch in the presence of disorder, $\Delta k_{f}(\omega)=\Delta k_{i}(\omega)+\delta k(\omega)$, vanishes, that is $\Delta k_{f}(\omega_{0}+\delta\omega_{0})=0$. To validate this argument, one would have to employ full-wave simulations of the nonlinear optical interaction of topological modes and average the results over an ensemble of disorder realizations large enough to achieve convergence of the ensemble average. This computational analysis would be rather unfeasible as just a single full-wave simulation requires several days to complete. On the other hand, for large values of disorder strength, concepts such as optical modes, topological properties, and photonic band gaps cease to exist, and therefore we do not investigate this case any further.

III.4 Second-harmonic generation in the slow-light regime

The slow-light regime, characterized by a significantly reduced group velocity, $v_{g}=d\omega/dk$, can be particularly effective in enhancing the efficiency of nonlinear wave interactions. In the context of SHG, this can be achieved when $v_{g}$ is reduced at one or both interacting waves. As Fig. 3c suggests, when one varies the location of the edge termination, the shift of the dispersion curve of the edge modes is accompanied by a change of its shape. For example, we find that for $t=0.24a$, the dispersion curve of the edge mode at the FF has a plateau (see Fig. 5a) leading to a peak of the group index, $n_{g}=c/v_{g}$, at $\overline{\omega}_{g}$ (see Fig. 5b, where we define the slow-light regime by the condition $n_{g}>20$ lp16prb ). On the other hand, the fact that the phase-matching spectral point $\overline{\omega}_{k}\equiv\overline{\omega}(\Delta\overline{k}=0)\neq\overline{\omega}_{g}$ provides a valuable approach to study the interplay between slow-light and phase-matching effects in the enhancement of SH conversion efficiency.

To this end, we fix the power of FF at $P_{1}=$1\text{,}\mathrm{W}$$ and compare the conversion efficiencies at two frequencies $\overline{\omega}_{l}$ and $\overline{\omega}_{h}$, lower and higher as compared to $\overline{\omega}_{k}$, respectively, chosen such that $\Delta\overline{k}(\overline{\omega}_{l})=-\Delta\overline{k}(\overline{\omega}_{h})$ (see the sketch in Fig. 5c). As such, phase-matching has the same influence on the wave interaction in the two cases and the conversion efficiency enhancement, $\eta=P_{2}(\overline{\omega}_{l})/P_{2}(\overline{\omega}_{h})$, is purely due to slow-light effects. The results in Fig. 5c show that slow-light contribution to $\eta$ can be larger than $20\times$. Alternative scenarios to enhance the SH conversion efficiency can also be devised, e.g., by using PhCs for which $\overline{\omega}_{k}=\overline{\omega}_{g}$ or PhCs for which the edge modes at both the FF and SH are slow-light edge modes.

III.5 Second-harmonic generation via interaction between forward- and backward-propagating edge modes

We now move on to an important class of nonlinear processes, which are challenging to achieve in regular optical media, and demonstrate SHG via interaction of backward-propagating edge modes. To this end, we exploit the existence of photonic gaps with negative Chern number in our system, as per Fig. 2e. In particular, we explore a case where the gap at the FF has $C=1$, while the gap at the SH has $C=-1$. The corresponding edge modes are shown in Fig. 6c, which illustrates the signature of modes with negative Chern number, i.e., the slope of the mode dispersion curve (at $\Omega_{2}$) is negative. The simulation results of the field profiles are presented in Figs. 6a and 6b. It clearly shows that whereas the mode at the FF propagates clockwise, as it is a forward-propagating mode, the SH wave propagates counterclockwise because in the left-half region of the simulation domain there is no field at the FF and consequently no nonlinear polarization at the SH (note that we placed the source of the FF wave at the middle of the top edge and used absorbing boundary condition for the bottom edge). We also confirm that the phase matching mechanism is involved in this unusual mode interaction regime, as per the results in Fig. 6d. This backward-propagating mode regime is promising for practical applications where one requires to separate the FF mode from the mode generated at the SH.

IV Experimental considerations

The key idea of this work involves coupling the edge states within different topological frequency band gaps via optical nonlinearity. This central idea of our study is rather general and can potentially be implemented in different experimental platforms available for topological photonics. In this section, for the sake of completeness, we first present a set of materials and parameters that can be used to experimentally implement our model system. Then we discuss further possible experimental platforms that can be used to implement the theoretic ideas and results presented in our work.

Since we separate the magnetic and nonlinear material components as the background and cylindrical regions of the PhC, respectively, possible experimental implementations could be readily conceived considering the fact that at microwave frequencies, magnetic materials to realize topological band gaps wang_Nature09 ; Chan_PRL11 ; Skirlo_PRL15 ; Bahari_Science17 ; feifei_NC18 and nonlinear materials to realize frequency mixing Boyd_PRL71 ; Bridges_APL72 ; Boyd_APL72 ; Ahn_JAP83 are routinely used. For example, to demonstrate the topology-protected SHG, one can use as the background medium yttrium iron garnet (YIG) microwave_book , a magnetic material widely used in recent experiments in topological photonics wang_Nature09 ; Chan_PRL11 ; Skirlo_PRL15 ; Bahari_Science17 ; feifei_NC18 . This material has relative permittivity $\epsilon=15$ and saturation magnetization of $4\pi M_{s}=$1780\text{,}\mathrm{G}\mathrm{a}\mathrm{u}\mathrm{s}\mathrm{s}. At frequency of $10\text{\,}\mathrm{GHz}$ and an external magnetic field of $H_{0}=$1000\text{\,}\mathrm{O}\mathrm{e} along the $z$-axis, the components of the permeability of YIG in the $x-y$ plane are $\mu_{\text{diag}}=0.85$ and $\mu_{\text{off-diag}}=0.54$. For the nonlinear material, one could use NaNO2, which has relative permittivity of $\epsilon=4.18$ and $\chi^{(2)}=$3.2\text{\times}{10}^{-22}\text{,}\mathrm{C}\text{,}{\mathrm{V}}^{-2} at frequency of tens of $\mathrm{GHz}$ [Boyd_PRL71 ](#bib.bib39). Setting the radius of the cylinders to $r=0.35a$, with $a\approx$3\text{\,}\mathrm{mm}, the system has topological band gaps at $\overline{\omega}_{0}\in[0.195,~{}0.23]$ and $\overline{\Omega}_{2}\in[0.385,~{}0.41]$, so that SHG is achieved in the frequency interval of $\overline{\omega}_{0}\in[0.195,~{}0.205]$. The effects of material losses and frequency dispersion in a ferrite at microwave frequencies have been discussed in Ref. Wang_PRL08 , showing that for YIG, the decay length is around $1300a$, thus far exceeding practical structural dimensions, and the band gap width slightly decreases by about $6\text{\,}\mathrm{\char 37\relax}$.

Furthermore, recent advances in different experimental platforms for topological photonics Lu_NP14 ; Shvets_NP17 ; Chen_PQE17 ; Gong_AOM17 ; Chen_OE18 ; Rider_JAP19 ; Ozawa_RMP19 provide a variety of choices available to implement the idea of topologically protected nonlinear frequency mixing processes in diverse photonic systems. More specifically, as nonlinearity exists in many optical materials and is easy to incorporate in an experimental platform, the key task reduces to designing photonic systems with several topological frequency band gaps that can satisfy the frequency and phase matching conditions. Regarding this requirement, we stress that several recent experiments have demonstrated the existence of topological band gaps using coupled waveguides or resonators Rechtsman_Nature13 ; Hafezi_NP13 ; Noh_PRL18 . In these tight-binding systems, further topological band gaps at higher frequencies can readily be created by considering waveguide or resonance modes at higher frequencies.

In what follows, we outline and briefly discuss several such experimental configurations in which this can potentially be achieved. i) Nonlinear and magnetic metamaterials: magnetism and nonlinearity can be easily implemented using metamaterials Smith_PRL2000 ; Linden_Science06 ; Lapine_RMP13 ; Gao_PRL15 , which provides an alternative to using regular materials, such as YIG and NaNO2 suggested above. For example, to create nonlinear metamaterials, approaches, such as insertion of nonlinear elements, nonlinear host medium, and nonlinear transmission lines are commonly used Lapine_RMP13 . Artificial magnetism at optical frequencies can be created using metamaterials based on split-ring resonators Smith_PRL2000 . Harmonic generation and topological edge states have also been studied in such types of metamaterials Linden_Science06 ; Gao_PRL15 . ii) Graphene plasmonic crystals Yeung_NanoL14 ; Jin_PRL17 ; Pan_NC17 ; You_arxiv19 : in Ref. Jin_PRL17 (see Fig. 3c therein) it has been demonstrated the existence of topological edge modes within different band gaps of a graphene plasmonic crystal. Moreover, in another study Yeung_NanoL14 (see Fig. 2a therein) the existence of two Dirac points at the frequencies of $3.4\text{\,}\mathrm{THz}$ and $6.8\text{\,}\mathrm{THz}$ has been demonstrated. These Dirac points can be gapped out to form two topological band gaps (ideal for SHG) under an external magnetic field as shown in Jin_PRL17 ; Pan_NC17 . In fact, recently, we have studied the topologically protected four-wave mixing process in a graphene metasurface You_arxiv19 . iii) Mimicking time-reversal breaking and synthesizing magnetic fields for photons Rechtsman_Nature13 ; Hafezi_NP11 ; Hafezi_NP13 : using an array of evanescently coupled helical waveguides, topological one-way edge states similar to that presented in Fig. 3 can be created without the need of an external magnetic field (see, e.g., Fig. 2b in Ref. Rechtsman_Nature13 ). Alternatively, one can also create one-way edge states using synthetic magnetic fields rather than a real magnetic field in an array of coupled optical-ring resonators (e.g., see Fig. 4 in Ref. Hafezi_NP13 ). Furthermore, additional topological band gaps at higher frequencies can be readily created by considering waveguide modes at corresponding frequencies. iv) Photonic quantum valley Noh_PRL18 ; Shalaev_NN19 ; He_NC19 ; Chen_AOM19 or spin Yuri_PRL2019 ; Hu_PRL15 ; Barik_Science18 ; Peng_PRL19 ; Noh_NP18 ; Mittal_PRL19 ; Shi_PNAS17 Hall (QVH or QSH) crystals with time-reversal symmetry: recent experiments on valley-Hall-like photonic systems have demonstrated the existence of one Noh_PRL18 ; Shalaev_NN19 ; He_NC19 or two Chen_AOM19 (see Fig. 2b therein) topological band gaps and corresponding edge modes within the gaps. As for the photonic QSH systems, the experiments in Yuri_PRL2019 ; Barik_Science18 ; Peng_PRL19 have demonstrated the existence of one topological band gap. Moreover, the proposal in Shi_PNAS17 (see Figs. 3 and 4 therein) using optically-passive elements to realize the optical version of the QSH insulator has shown the existence of multiple topological band gaps.

V Conclusions

In this work, we have demonstrated topology-protected nonlinear frequency conversion processes via one-way edge modes of topological photonic crystals. Apart from the proof-of-concept implementations, such as SHG and THG, we also showed that more complex behaviors, such as slow-light effects and counter-propagating mode interaction, can also be realized within the setup. A special aspect of nonlinear processes, i.e., the phase-matching condition, requires a new level of control of the edge modes, which has not been discussed previously. This condition requires a method to tune the dispersion of the edge modes, which we found can be conveniently achieved by tailoring the geometry of the edge termination. Our work reveals that the coherence length characterizing the nonlinear optical interactions considered in this study, which crucially determines the efficiency of the topologically protected nonlinear frequency mixing processes, is robust against structural defects. Our proposed setup provides a platform for studying additional phenomena, e.g., when the frequency gap has large Chern number ($|C|\neq 0,1$ Skirlo_PRL04 ) one can explore how to excite one of the several edge modes in the gap or how to couple edge modes belonging to different gaps via the nonlinearity of the medium.

Importantly, nonlinear interactions of topological modes, such as sum- and difference-frequency generation, high-harmonic generation, and four-wave-mixing, can be readily implemented within our setup. Our work may also stimulate the search for other lattice geometries or setups where one can optimize the gap properties for specific applications. For example, in the Chern number graded gap phase diagrams of Fig. 2e, apart from the gap of $C=1$, other gaps with $C=-1,-2,2,3$ are typically narrow and appear at high frequencies, so designing setups where these gaps are wide and are formed at low frequencies is particularly relevant from experimental point of view. Beyond the experimental implementation of our model system at microwave frequencies using magnetic and nonlinear materials, we also discussed several different possible implementations using diverse experimental platforms available for topologic photonics. Last but not least, the concept of topology-protected nonlinear frequency mixing is very general in that it applies not only to photonics, but also to plasmonics Jin_PRL17 ; Pan_NC17 ; Kauranen_NP12 ; Panoiu_JO18 , phononics Zhang_PRL15 ; Bojahr_PRL15 ; Ge_NSR18 , magnonics Chumak_NP15 ; Sebastian_PRL13 ; Wang_PRAPP18 and exciton-polariton systems Karzig_PRX15 ; Klembt_Nature18 ; Warkentin_PRB18 ; Knuppel_Nature19 , thus we expect that our work will have a broad impact.

Acknowledgements

This work was supported by the European Research Council (ERC) (Grant no. ERC-2014-CoG- 648328). We acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL) and associated support services in the completion of this work.

Appendix A Computation of Chern numbers

The Chern number of the $n$th-band of the photonic crystal is defined by Eq. (5) based on Refs. Haldane_PRL08 ; Lu_NP14 . An efficient numerical algorithm to compute it was introduced in Ref. Fukui_JPSJ05 . In this method, from the eigenmode $|{\mathbf{E}_{n}(\mathbf{k})}\rangle$ of the $n$th band, which can be accessed directly from the eigenfrequency solver of COMSOL, one can define a $U(1)$ link variable,

[TABLE]

where $\mathbf{k}_{l}$ is a lattice point in the discretized Brillouin zone and $\mathbf{e}_{\alpha}$ is the lattice displacement in the direction $\alpha$ $(\alpha=1,2)$. Furthermore, a lattice field strength can be defined by the link variable:

[TABLE]

where the lattice field strength is defined as the principal branch of the logarithm $-\pi<F_{12}((\mathbf{k}_{l})/i\leq\pi$. Finally, the Chern number of the $n$th band can be calculated from the lattice field strength according to:

[TABLE]

where the sum is taken over all the lattice points in the discretized Brillouin zone. The $C_{n}$ defined above is manifestly gauge-invariant and strictly an integer for arbitrary lattice parameters. This is because if we introduce a gauge potential

[TABLE]

one can get

[TABLE]

where $\Delta_{\alpha}f(\mathbf{k}_{l})=f(\mathbf{k}_{l}+\mathbf{e}_{\alpha})-f(\mathbf{k}_{l})$ and $n_{12}(\mathbf{k}_{l})$ is an integer valued field, chosen in such a way that $F_{12}(\mathbf{k}_{l})/i$ takes a value within the principal branch. Consequently, we can get

[TABLE]

which shows that $C_{n}$ is an integer. Certainly, this does not mean that any coarse discretization of the first Brillouin zone will ensure a converged Chern number; nevertheless, asymptotic convergence of the Chern number requires only a moderately dense discretization. In our calculations of the data presented in Fig. 2e, we find that a $20\times 20$ discretization of the first Brillouin zone suffices.

Appendix B Coupled-mode theory for the third-harmonic generation

The coupled-mode equations for the third-harmonic generation can be derived in a way similar to that for the second-harmonic generation presented in Section III.1. Here we outline the main steps pmw10jstqe ; lp16prb . First, one can write down the expressions for the fundamental ($f$) and third ($t$) harmonic waves, $(\mathbf{E}_{1}^{f},\mathbf{H}_{1}^{f},\mathbf{E}_{2}^{f},\mathbf{H}_{2}^{f},\mathbf{P}_{NL}^{f})$ and $(\mathbf{E}_{1}^{t},\mathbf{H}_{1}^{t},\mathbf{E}_{2}^{t},\mathbf{H}_{2}^{t},\mathbf{P}_{NL}^{t})$, in the following Bloch forms,

[TABLE]

and

[TABLE]

where $\Omega_{3}=3\omega_{0}$ and $\mathbf{P}_{NL}^{f},\mathbf{P}_{NL}^{t}$ can be derived from $\mathbf{P}_{NL}(\mathbf{r},t)=\hat{\chi}^{(3)}(\mathbf{r},t)\vdots\mathbf{E}(\mathbf{r},t)\mathbf{E}(\mathbf{r},t)\mathbf{E}(\mathbf{r},t)$ with $\mathbf{E}(\mathbf{r},t)=\mathbf{E}_{2}^{f}(\mathbf{r},t)+\mathbf{E}_{2}^{t}(\mathbf{r},t)$ and $\hat{\chi}^{(3)}$ the third-order susceptibility tensor. After some simple algebra, the nonlinear polarizations at the fundamental and third-harmonic frequencies can be written as:

[TABLE]

If we now consider the fields $(\mathbf{E}_{1}^{f},\mathbf{H}_{1}^{f},\mathbf{E}_{2}^{f},\mathbf{H}_{2}^{f},\mathbf{P}_{NL}^{f})$, the l.h.s. of Eq. (19) gives:

[TABLE]

where $\Delta k=k_{t}(\Omega_{3})-3k_{f}(\omega_{0})$, and the r.h.s. of Eq. (19) leads to:

[TABLE]

Comparing Eqs. (44) and (45), we obtain the coupled-mode equation for the slowly-varying mode amplitude $A_{f}(x)$:

[TABLE]

where the nonlinear coefficient at the fundamental frequency is

[TABLE]

Similarly, if we consider the fields $(\mathbf{E}_{1}^{t},\mathbf{H}_{1}^{t},\mathbf{E}_{2}^{t},\mathbf{H}_{2}^{t},\mathbf{P}_{NL}^{t})$, the l.h.s. of Eq. (19) gives:

[TABLE]

and the r.h.s. of Eq. (19) can be written as:

[TABLE]

Finally, from Eqs. (48) and (49), we get the coupled-mode equation for the slowly-varying amplitude $A_{t}(x)$ as:

[TABLE]

where the nonlinear coefficient at the third-harmonic frequency is:

[TABLE]

We can also express the nonlinear coefficients in terms of the group indices and the effective third-order susceptibilities, $\chi_{\mathrm{eff},f/t}^{(3)}(x)$, of the two interacting modes:

[TABLE]

where

[TABLE]

Introducing the averaged physical quantities,

[TABLE]

and then averaging Eqs. (46) and (51) over one lattice constant, we obtain the coupled-mode equations governing the nonlinear dynamics of the interacting modes:

[TABLE]

where $\overline{A}_{f}(x)$ and $\overline{A}_{t}(x)$ are the averaged mode amplitudes at the fundamental and third-harmonic frequencies, respectively.

Appendix C Comparison between rigorous numerical simulations and the coupled-mode theory

In this Appendix, we present a computational analysis that illustrates how the coupled-mode theory derived above can be used to explain the full-wave dynamics obtained using COMSOL.

C.1 Second-harmonic generation described by the coupled-mode theory

The key quantities that characterize the coupled-mode theory describing the second-harmonic generation are $\chi^{(2)}_{\textrm{eff},f}(x)$ and $\chi^{(2)}_{\textrm{eff},s}(x)$, as defined by Eqs. (30a) and (30b), respectively. As we consider photonic crystals that are periodic in space, $\chi^{(2)}_{\textrm{eff},f}(x)$ and $\chi^{(2)}_{\textrm{eff},s}(x)$ are periodic functions of $x$, and we only need to show their $x$-dependence in one unit cell. Thus, in Fig. 7 we depict the $x$-dependence of $\chi^{(2)}_{\textrm{eff},f}(x)$ and $\chi^{(2)}_{\textrm{eff},s}(x)$ in one unit cell. As one can see, while the real parts of $\chi^{(2)}_{\textrm{eff},f}(x)$ and $\chi^{(2)}_{\textrm{eff},s}(x)$ are even functions, with respect to the center of the unit cell, the imaginary parts of $\chi^{(2)}_{\textrm{eff},f}(x)$ and $\chi^{(2)}_{\textrm{eff},s}(x)$ are odd functions. In fact, it can be easily demonstrated that $\chi^{(2)}_{\textrm{eff},f}(x)={\chi^{(2)}_{\textrm{eff},s}}^{*}(x)$. Therefore, the averages of $\chi^{(2)}_{\textrm{eff},f}(x)$ and $\chi^{(2)}_{\textrm{eff},s}(x)$ in one unit cell are real numbers and are equal to each other.

After one computes $\overline{\gamma}^{(2)}_{f,s}$ and $\overline{\chi}^{(2)}_{\mathrm{eff},f/s}$ by averaging $\gamma^{(2)}_{f,s}(x)$ and $\chi^{(2)}_{\mathrm{eff},f/s}(x)$, respectively, in one unit cell, one can straightforwardly solve the coupled-mode equations expressed as Eqs. (32). Thus, we present in Fig. 8 the evolution of the generated second-harmonic wave as a function of propagation distance. It can be seen in this figure that there is a good agreement between the coupled-mode theory and rigorous numerical simulations both in regard of the oscillation period and the amplitude of the power of the second-harmonic wave.

C.2 Third-harmonic generation described by the coupled-mode theory

Similar to the case of SHG, the key quantities that characterize the coupled-mode theory of third-harmonic generation are $\chi^{(3)}_{\textrm{eff},f}(x)$ and $\chi^{(3)}_{\textrm{eff},t}(x)$, as defined by Eqs. (53). Both these physical quantities are periodic functions of the $x$-coordinate, their $x$-dependence being presented in Fig. 9. Also similar to the case of SHG, while the real parts of $\chi^{(3)}_{\textrm{eff},f}(x)$ and $\chi^{(3)}_{\textrm{eff},t}(x)$ are even functions, the imaginary parts of $\chi^{(3)}_{\textrm{eff},f}(x)$ and $\chi^{(3)}_{\textrm{eff},t}(x)$ are odd functions. Therefore, the averages of $\chi^{(3)}_{\textrm{eff},f}(x)$ and $\chi^{(3)}_{\textrm{eff},t}(x)$ in one unit cell are real numbers and are equal to each other.

We have solved Eqs. (55) using $\overline{\gamma}^{(3)}_{f,t}$ obtained by averaging $\gamma^{(3)}_{f,t}(x)$ in one unit cell and present the results in Fig. 10. We again find a good agreement between the coupled-mode theory and rigorous numerical simulations regarding both the oscillation period and amplitude of the power of the third-harmonic wave.

C.3 Coupled-mode theory of slow-light nonlinearity enhancement

The group index $n_{g}$ of the edge modes in a generic case, e.g., Fig. 4, is limited to 5$\sim$7. However, we have shown in the main text that the shape of the dispersion curve of the one-way edge modes can be tailored, leading to much larger $n_{g}$, as per Fig. 5. As it is well known, this will enhance the efficiency of the nonlinear process. We plot in Fig. 11a the function $\chi^{(2)}_{\textrm{eff},f/s}(x)$ at two different frequencies, $\bar{\omega}_{\textrm{slow}}=0.197$ and $\bar{\omega}_{\textrm{fast}}=0.21$ with $n_{g}^{f}(\bar{\omega}_{\textrm{slow}})\simeq 174$ and $n_{g}^{f}(\bar{\omega}_{\textrm{fast}})\simeq 12$ [$n_{g}^{s}(\bar{\omega}_{\textrm{slow}})\simeq 5$ and $n_{g}^{s}(\bar{\omega}_{\textrm{fast}})\simeq 4$ have similar values at these two frequencies]. While the amplitudes of $\chi^{(2)}_{\textrm{eff},f/s}(x)$ at these two frequencies are comparable, from Fig. 11b, which shows the nonlinear coefficient $\gamma_{f/s}^{(2)}(x)$ of Eqs. (24)/(28), one can see that $\gamma_{f/s}^{(2)}(x)$ at $\bar{\omega}_{\textrm{slow}}$ is significantly enhanced as compared to $\gamma_{f/s}^{(2)}(x)$ at $\bar{\omega}_{\textrm{fast}}$. The ratio of the averaged nonlinear coefficient in one unit cell according to Eq. (31) and Fig. 5b is $\bar{\gamma}^{(2)}_{f,s}(\bar{\omega}_{\textrm{slow}})/\bar{\gamma}^{(2)}_{f,s}(\bar{\omega}_{\textrm{fast}})\simeq 18$, which roughly agrees with $n_{g}^{f}(\bar{\omega}_{\textrm{slow}})/n_{g}^{f}(\bar{\omega}_{\textrm{fast}})\simeq 15$. This enhancement of the nonlinear coefficient leads to the enhancement of the efficiency of the SHG, see, e.g., Figs. 5c and 5d.