Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Theory of soliton propagation in nonlinear photonic crystal waveguides

Vakhtang Jandieri, Ramaz Khomeriki, Jamal Berakdar, and Daniel Erni

Open Access Open Access

Abstract

Propagation of the temporal soliton in Kerr-type photonic crystal waveguide is investigated theoretically and numerically. An expression describing the evolution of the envelope of the soliton based on the full-wave modal analysis, taking into account all space-harmonics, is rigorously obtained. The nonlinear coefficient is derived, for the first time, based on a modification of the refractive indices for each space-harmonic due to the Kerr-type nonlinearity. For illustrating the general formulation and results, we performed extensive computational electromagnetics simulations for the propagation of gap solitons in an experimentally feasible photonic crystal waveguides, endorsing the correctness and usefulness of the proposed formalism.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

After John Scott Russell’s experimental discovery and the subsequent mathematical formulations by Rayleigh and Boussinesq, the generality of solitonic behaviour has been recognized in many branches of physics, electronics and biology [1–3]. Several theoretical and numerical techniques have been developed so far to analyze the formation and the propagation of solitons (kink-solitons, spatial and spatio-temporal solitons, gap solitons, etc.) in different kind of media [4]. In general, the main requirement for the propagation of solitons is the one- or quasi-one dimensionality of the medium except for the spatial soliton [5] for which a beam propagation direction plays the role of time. In photonics a quasi-one dimensionality can be achieved either by sharp refractive index change between the guiding layer and the surrounding medium [6], or by the photonic crystal (PhC) along (PhC fiber [7]) or perpendicular (PhC waveguide [8]) to the direction of the soliton propagation.

The transmission of localized energy in a structure with a periodically varying refractive index, such as 1D Bragg grating [9] or 2D PhC [10], occurs due to the modulation of corresponding Bloch modes. In weakly nonlinear systems, the problem could be solved by multiple-scale perturbation method and the soliton propagation is described by the nonlinear Schrödinger equation for a slowly-varying envelope of the soliton [11,12]. Planar PhC waveguides, which are designed by removing one or a few rows from the original PhCs, are considered as very promising candidates for various realizations of optical devices on a chip. Moreover, PhC waveguides guarantee a very good confinement of the modes in the slow light regime [13].

The main tool to deal with the nonlinearities is variational Ansatz, from which a soliton shape is obtained [14,15]. Here, for the first time, a theoretical framework for the propagation of temporal solitons in PhC waveguides based on a multi-harmonic treatment of the nonlinear setup is presented. The method is based on a full-wave modal analysis using the Floquet-Bloch theorem. A rigorous expression for the nonlinear coefficient (due to Kerr-type nonlinearity) based on a modification of the refractive indices of each space-harmonic is derived. For the sake of confirmation and illustration we performed numerical simulations for the propagation of gap-solitons. The gap-solitons are formed if the power of the launched signal exceeds a certain threshold value [16–18]. Their formation time is relatively short and the registered gap soliton inside the guiding structure can be used for comparison with the analytical solution of the nonlinear Schrödinger equation. A very good agreement between the numerical and theoretical results is demonstrated.

2. Formulation of the problem

Without loss of generality we analyze three coupled symmetric PhC waveguides composed of a hexagonal lattice made of circular air-holes that are periodically distributed along the ${x}$-axis with a common period ${h}$, as shown in Fig. 1. The structure as such is already well studied [19]. The guiding regions ${(a)}$, ${(b)}$ and ${(c)}$ having the same width ${w}$ are separated by the barrier layers of two PhCs with a barrier thickness according to the number of layers $N_B$. The radius of the air-holes is ${r}$ and $\epsilon _s=\ n_s^{2}$ is the relative dielectric permittivity of the background material. The number of layers ${N}$ of the upper and lower PhCs is taken large enough in order to minimize the leakage of the power along the transverse ${y}$ - axis. The modulated ${x}$ component of the electric field $E_x(x,y,t)$ propagating in the nonlinear coupled PhC waveguides can be written in the following form:

$$E_x= \frac{\Psi}{2} \sum_{m} (u^{+}_m e^{ik_{ym}y} + u^{-}_m e^{{-}ik_{ym}y} ) e^{i\beta_m x -i\omega t} + c.c.$$
where $\Psi \equiv \Psi (x,t)$ is a slowly varying amplitude of its arguments, $k_{ym} = \sqrt {k_{s}^{2} - \beta _{m}^{2}}$, $\beta _m = k_{x0} + \frac {2m\pi }{h}$, $k_s = \omega n_s\sqrt {\epsilon _0 \mu _0}$, $\omega$ is the angular frequency, $k_{x0}$ is the mode propagation constant along the ${x}$-axis, $k_{ym}$ is the transverse wavenumber of the ${m}$ - th space-harmonic and ’$c.c.$’ stands for complex conjugate. We symbolize by $\mathbf {u^{+}}$ and $\mathbf {u^{-}}$ the vectors with components $\left \{u^{+}_m\right \}$ and $\left \{u^{-}_m\right \}$, respectively. Thus, $\mathbf {u^{+}}$ and $\mathbf {u^{-}}$ are the amplitude vectors of the up-going and down-going space-harmonics along the ${y}$-axis. Let us denote by $\mathbf {a^{+}}$ and $\mathbf {a^{-}}$ the amplitude vectors of the Floquet modes defined at the upper and lower interfaces of the guiding region ${(a)}$, by $\mathbf {b^{+}}$ and $\mathbf {b^{-}}$ those of the guiding layer ${(b)}$ and by $\mathbf {c^{+}}$ and $\mathbf {c^{-}}$ those of the guiding layer ${(c)}$, as shown in Fig. 1. From Maxwell’s equations one infers the expression for the ${y}$-component of the electric field and finds the corresponding harmonics as:
$$E^{{\pm}}_{xm}=\Psi u_m^{{\pm}}; \qquad E_{ym}^{{\pm}} ={\mp} \Psi \frac{\beta_m}{k_{ym}} u^{{\pm}}_m$$
The ${x}$-component of the electric displacement field ${D_x}$ reads:
$$D_x = E_x \left[ \chi^{(1)} + \chi^{(3)} ( E_x^{2} + E_y^{2}) \right]$$
where $\chi ^{(i)}$ is the i-th order optical susceptibility (particularly, $\chi ^{(1)}\equiv \epsilon _s=n_s^{2}$). In our analysis we assume that the second-order susceptibility term $\chi ^{(2)} = 0$, which is true in centrosymmetric crystals such as silicon [20]. Substituting the electric field components defined from (1) and (2) into (3) and collecting terms for the $x$-component of the electric displacement field $D_{xm}^{\pm }$ for the up-going and down-going m-th space harmonic containing the same $e^{i(\beta _m x-\omega t)} e^{\pm i k_{ym}y}$, we obtain:
$$\begin{aligned} D_{xm}^{{\pm}} = \Psi u_m^{{\pm}} \left[ n_s^{2} + \chi^{(3)}|\Psi|^{2} \left( \frac{3}{4} \left[ 1 + \frac{\beta_m^{2}}{k_{ym}^{2}} \right] |u_m^{{\pm}}|^{2} + \frac{1}{2}\left[ 3 - \frac{\beta_m^{2}}{k_{ym}^{2}} \right] |u_m^{{\mp}}|^{2} \right. \right.\\ + \sum_{\mu \neq m} \frac{1}{2} \left[ 3 + \frac{\beta_{\mu}^{2}}{k_{y\mu}^{2}} \right] \left(|u_{\mu}^{{\pm}}|^{2} + |u_{\mu}^{{\mp}}|^{2}\right) + \sum_{\mu \neq m} \frac{\beta_m \beta_{\mu}}{k_{ym}k_{y\mu}} \left(|u_{\mu}^{{\pm}}|^{2} - |u_{\mu}^{{\mp}}|^{2}\right) \Biggr) \Biggr] \end{aligned}$$
Note that in (4) all the higher harmonic terms are neglected with respect to the leading one. If the guided modes are well confined by the upper and lower PhCs with a number of layers $N$ (in our studies $N$ = 5), $u^{-}_m\approx u^{+}_m$ and the last term in the square brackets of (4) can be considered as negligibly small. Thus, a refractive index for the up-going and down-going m-th space harmonic in each guiding region can be written in the following form:
$$\begin{aligned} n_{m,\nu}^{{\pm}} = n_s + \frac{\chi^{(3)} |\Psi|^{2}}{4n_s} \Biggl[ \frac{3}{2} \Biggl[ 1 + \frac{\beta_m^{2}}{k_{ym}^{2}} \Biggr] |u_{m,\nu}^{{\pm}}|^{2} + \Biggl[ 3 - \frac{\beta_m^{2}}{k_{ym}^{2}} \Biggr] |u_{m,\nu}^{{\mp}}|^{2} \Biggr. \Biggl.\\ + \sum_{\mu \neq m} \Biggl[ 3 + \frac{\beta_{\mu}^{2}}{k_{y\mu}^{2}} \Biggr] \left( |u_{\mu,\nu}^{{\pm}}|^{2} + |u_{\mu, \nu}^{{\mp}}|^{2} \right) \Biggr] \end{aligned}$$
where $\nu = a,b,c$. Here, the scattering amplitude for the guiding region “$a$” is $u_{m,a}^{\pm } = a_{m}^{\pm }$, for the guiding region “$b$” is $u_{m,b}^{\pm } = b_{m}^{\pm }$ and for the guiding region “$c$” is $u_{m,c}^{\pm } = c_{m}^{\pm }$ (Fig. 1). Note that $D_{ym}^{\pm }$ is expressible through $D_{xm}^{\pm }$ exploiting the Maxwell’s equations $D_{ym}^{\pm } = \mp \frac {\beta _m}{k_{ym}}D_{xm}^{\pm }$. Besides, from (5) it follows that the refractive indices for the up-going and down-going m-th space harmonics are very close to each other $n_{m,\nu }^{+}\approx n_{m,\nu }^{-}$ and are proportional to the wave intensity $|\Psi |^{2}$, which leads to a corresponding shift in the angular frequency $\omega '= \omega + \delta \omega$ at a given propagation constant $k_{x0}$ due to the involved Kerr-type nonlinearity. Applying a ray tracing of the space-harmonics along the transverse ${y}$-axis, the relations between the scattering amplitudes in three guiding layers are given in terms of the reflection and transmission matrices of the PhCs as:
$${\mathbf{a^{-}}} = \mathbf{W}(\omega', k_{x0},n_{m,a}^{{\pm}} ) \mathbf{R}(\omega', k_{x0}, n_s)\cdot {\mathbf{a^{+}}}$$
$${\mathbf{a^{+}}} = {\mathbf{W}}(\omega',k_{x0},n_{m,a}^{{\pm}}) \left[ \mathbf{ R_B } (\omega', k_{x0}, n_s) \cdot{ \mathbf{a^{-}}} + {\mathbf{ T_B }} (\omega', k_{x0}, n_s) \cdot{ \mathbf{b^{+}}} \right]$$
$${\mathbf{b^{+}}} = {\mathbf{W}}(\omega',k_{x0},n_{m,b}^{{\pm}}) \left[ \mathbf{ T_B } (\omega', k_{x0}, n_s) \cdot{ \mathbf{c^{+}}} + {\mathbf{ R_B }} (\omega', k_{x0}, n_s) \cdot{ \mathbf{b^{-}}} \right]$$
$${\mathbf{b^{-}}} = \mathbf{W}(\omega',k_{x0},n_{m,b}^{{\pm}}) \left[ \mathbf{ T_B } (\omega', k_{x0}, n_s) \cdot{ \mathbf{a^{-}}} + {\mathbf{ R_B }} (\omega', k_{x0}, n_s) \cdot { \mathbf{b^{+}}} \right]$$
$${ \mathbf{c^{+}}} = \mathbf{W}(\omega', k_{x0},n_{m,c}^{{\pm}}) \mathbf{R}(\omega', k_{x0}, n_s) \cdot {\mathbf{ \mathbf{c^{-}}}}$$
$${\mathbf{c^{-}}} = \mathbf{W}(\omega',k_{x0},n_{m,c}^{{\pm}}) \left[ \mathbf{ T_B } (\omega', k_{x0}, n_s) \cdot{ \mathbf{b^{-}}} + \mathbf{ R_B } (\omega', k_{x0}, n_s) \cdot {\mathbf{c^{+}}} \right]$$
with
$$\mathbf{W}(\omega', k_{x0},n_{m,\nu}^{{\pm}}) = [ e^{ik_{ym} w} ], \hspace{25pt}\nu=a,b,c$$
where ${\mathbf {W}}(\omega ', k_{x0},n_{m,\nu }^{\pm })$ is a diagonal matrix, which defines the phase shift of each m-th up-going and down-going space harmonics within the guiding region(s), ${\mathbf {R}}(\omega ', k_{x0}, n_s)$ is the generalized reflection matrix for the upper and lower N-th layered PhCs, whereas the generalized reflection and transmission matrices of the PhCs barriers are ${\mathbf { R_B}}(\omega ', k_{x0}, n_s)$ and ${\mathbf { T_B }} (\omega ', k_{x0}, n_s)$, respectively. Note that the reflection and transmission matrices are functions of $\ n_s$, assuming that the shift of the refractive index due to the Kerr effect (according to the injected intensity $|\Psi |^{2}$) occurs only in the guiding region having the thickness $\ w-2r$ (Fig. 2). Then, we use our original formulation to accurately and efficiently calculate the generalized reflection and transmission matrices of the PhCs [21]. The formalism is based on the lattice sums (LSs) technique [22] combined with the transition matrix (T-matrix) and the generalized reflection matrix approach [23]. It is numerically slim and applicable to a wide class of periodic and bandgap structures from microwave to optical regions. The formulation is briefly discussed in Appendix A.

 figure: Fig. 1.

Fig. 1. (Left) Schematic view of the three symmetric PhC waveguides with guiding regions ${(a)}$, ${(b)}$ and ${(c)}$ having the same width ${w}$ and separated by the barrier layers of two PhCs with number of layers $N_B$. Planar PhCs are composed of hexagonal lattice of the circular air holes periodically distributed along $x$-axis with a period $h$. The air-holes are infinitely long along the ${z}$-axis. The radius of the air holes is $r$, $n_s$ is the refractive index of the background medium. Here ${{a^{+}}}$ and ${{a^{-}}}$, ${{b^{+}}}$ and ${{b^{-}}}$, ${{c^{+}}}$ and ${{c^{-}}}$ define the up-going and down-going space-harmonics in the guiding regions ${(a)}$, ${(b)}$ and ${(c)}$, respectively. (Right) Dispersion curves of the symmetric (blue curve and green curve) and the antisymmetric (red curve) super-modes for the $H$-polarized field. The operating frequency for the realization of functional all-optical logic gates is $\frac {h\omega }{2\pi c}=0.232$, and it is marked by a red dot. The distributions of the magnetic field $H_z$ for the super-modes, as well as the geometry of the setup are shown in the corresponding insets.

Download Full Size | PDF

 figure: Fig. 2.

Fig. 2. Schematic view of PhC waveguide with a guiding region ${(a)}$ having a refractive index $\tilde {n}_s=n_s + \delta n$, whereas the refractive index of the other regions is the same as in Fig. 1 and it is equal to $n_s$.

Download Full Size | PDF

The linear system of Eqs. (6)–(11) can be re-written in the following compact form:

$$\mathbf{\Omega}(\omega', k_{x0}, n_{m,\nu}^{{\pm}})\cdot \mathbf{A}^{T}(\omega', k_{x0}, n_{m,\nu}^{{\pm}}) = 0$$
where ${\mathbf {A}} = \left \{ \left [ \cdots a^{+}_{-m},\ldots , a^{+}_m\cdots \right ],\ldots , \left [\cdots b^{-}_{-m},\ldots , b^{-}_m\cdots \right ], ..\right \}$, “$T$” stands for the transpose of the vector, ${\mathbf {\Omega }}$ is a block matrix and the size of each sub-matrix is $(2M + 1) \times (2M + 1)$. Expanding (13) in terms of the perturbation $\delta \omega$ and $|\Psi |^{2}$, and vector multiplying with ${\mathbf {A'}}$ from the left while considering that ${\mathbf {A'}}\cdot \left [ {\mathbf {\Omega }}(\omega , k_{x0}, n_s) \right ]=0$, the following expression can be obtained:
$$\delta\omega\mathbf{A'}\frac{\partial\mathbf{\Omega}(\omega,k_{x0},n_s)}{\partial\omega}\mathbf{A}^{T} =\left.{-}|\Psi|^{2}\mathbf{A'}\frac{\partial\mathbf{\Omega}(\omega,k_{x0},n_{m,\nu}^{{\pm}})}{\partial\left(|\Psi|^{2}\right)}\mathbf{A}^{T}\right|_{\substack{n_{m,\nu}^{{\pm}}=n_s}}$$
It should be emphasized that the eigenvectors ${\mathbf {A}}'$ and ${\mathbf {A}}$ are numerically calculated for the truncated block matrix ${\mathbf {\Omega }}(\omega , k_{x0},n_s)$ after confirming the convergence of the solutions that have to satisfy $\det \left[ {\mathbf {\Omega }}(\omega , k_{x0}, n_s) \right] = 0$. The second term in (14) is rewritten as:
$$\frac{{\partial}\mathbf{\Omega}(\omega, k_{x0},n_{m,\nu}^{{\pm}})} {\partial\left(|\Psi|^{2}\right)}=\left. \displaystyle\sum_{m,\nu}\frac{\partial\left(n_{m,\nu}^{{\pm}}\right)}{\partial\left(|\Psi|^{2}\right)}\frac{\partial\mathbf{\Omega}(\omega, k_{x0}, n_{m,\nu}^{{\pm}})}{\partial\left(n_{m,\nu}^{{\pm}}\right)}\right|_{\substack{n_{m,\nu}^{{\pm}}=n_s}}$$
Note that $\frac {\partial \left (n_{m,\nu }^{\pm }\right )}{\partial \left (|\Psi |^{2}\right )}$ can be inferred directly from (5), whereas $\left.\frac {\partial {\mathbf {\Omega }}(\omega , k_{x0},n_{m,\nu }^{\pm })}{\partial \left (n_{m,\nu }^{\pm }\right )}\right |_{\substack {n_{m,\nu }^{\pm }=n_s}}$ should be numerically calculated for each ${m}$-th space-harmonic in each ${\nu }$-th guiding region. The calculation procedure follows the formulation [21,22], albeit a slight modification is needed of (6)–(11) by implementing the Fresnel matrices. Note, Fresnel matrices are diagonal matrices. Figure 2 displays a schematic view of only the guiding region ${(a)}$ having a refractive index $\tilde {n}_s=n_s + \delta n$, whereas the refractive index of the other regions is $n_s$. The calculation technique is briefly described in Appendix B and is applicable to the other guiding regions having different refractive indices.

Once the non-linear dispersion relation $\omega (k_{x0},|\Psi |^{2})$ for weak nonlinearities is defined from (13) and (14), there exists a well-established procedure to reduce the problem to a nonlinear Schrödinger equation using the multiple-scale treatment of the wave equations [11]. Following the analysis described in [24], we find the following expression:

$$i\left( \frac{\partial\Psi}{\partial t} + v_g\frac{\partial\Psi}{\partial x} \right) + \frac{\omega^{''}}{2}\frac{\partial^{2}\Psi}{\partial x^{2}} + \gamma |\Psi|^{2}\Psi = 0$$
where $\gamma = -\frac {\partial (\delta \omega )}{\partial (|\Psi |^{2})}$ is a nonlinear coefficient and can be directly deduced from (14), the group velocity $v_g = \frac {\partial \omega }{\partial k_{x0}}$ and the group velocity dispersion $\omega ^{''} = \frac {\partial ^{2}\omega }{\partial k_{x0}^{2}}$ can be retrieved from the dispersion diagram. The solution for the slowly varying amplitude $\Psi (x,t)$ is then written in the following form:
$$\Psi(x,t) = \frac{Fe^{{-}i(\Delta\omega)t}}{\cosh \left[(x-v_g t)/\Lambda \right]}$$
where
$$\Lambda = \sqrt{\frac{\omega^{''}}{\gamma F^{2}}}, \qquad \Delta\omega ={-}\frac{\gamma F^{2}}{2} = \frac{F^{2}}{2}\frac{\partial(\delta\omega)}{\partial(|\Psi|^{2})}$$
Note that $\Lambda$ stands for the width of the temporal soliton, $F$ denotes the soliton’s amplitude, and $\Delta \omega$ is the shift of the angular frequency due to Kerr nonlinearities. Finally, substituting (16) into (1) and (2), we arrive at the expressions for the electric fields from which the magnetic field is retrieved.

3. Numerical results and discussions

We analyze a 2D model system of three coupled PhC waveguides composed of experimentally feasible planar hexagonal lattice of air holes formed in a dielectric nonlinear background medium with a linear refractive index $n_s = 2.95$ (crystalline silicon) in conjunction with a Kerr-type nonlinearity [19]. This 2D model has proven to be a very good approximation of the original 3D structure. The thickness of the upper and lower PhCs is taken $N=5$ and the radius of the air-holes is $r = 0.32h$, where $h$ is a period of the PhCs. The barrier layers are composed of 1-layered structures, $N_B=1$, and the length of the PhC is $30h$. The dispersion diagram of the structure is shown in Fig. 1. Here, we are interested only in the symmetric mode described by the blue line, since this mode is responsible for the formation of the gap soliton. The dispersion diagram of this mode can be well approximated by the parabola taking into account the terms up to the square of the angular frequency and the group velocity dispersion amounts to $\frac {\partial ^{2}\omega }{\partial k_{x0}^{2}} \approx 1.9 \frac {hc}{2\pi }$. A full-wave computational electromagnetics analysis is conducted based on the finite-difference time-domain (FDTD) method [25] together with uniaxial (perfectly matched layer) PML at the operating normalized frequency $\frac {h\omega }{2\pi c} = 0.232$ (Fig. 1).

A continuous wave (CW) signal with $(H_z,E_x,E_y)$ is launched through the middle waveguide of the coupled PhC waveguides [guiding region ${(b)}$]. The injected peak power of the CW signal is chosen as $\chi ^{(3)}E_0^{2} = 0.1410$, which means that for silicon with the nonlinear refractive index $n_2={3}\cdot {10^{-18}}[m^{2}\cdot W^{-1}], E_0^{2} = {2.7}\cdot {10^{18}} [V^{2}\cdot m^{-2}]$ (we excite the waves very close to the left edge inside the PhCs). A nonlinear medium due to the Kerr effect results in a dispersion shift of the symmetric mode (blue line in Fig. 1) to the lower frequencies yielding the formation of the gap soliton. The shift of the angular frequency due to the Kerr-type nonlinearity can be calculated theoretically based on (17) and it amounts to $\frac {h\Delta \omega }{2\pi c}=0.0016$. Figure 3(a) depicts the magnetic field distribution of the gap soliton propagation. From the numerical simulations it follows that the maximum $F$ and the width $\Lambda$ of the gap soliton exhibiting a stable profile propagating inside the waveguide, amounts to $F=300 A/m$ and $\Lambda =4.10 h$, respectively.

 figure: Fig. 3.

Fig. 3. (a) Schematic distribution of the magnetic field of the gap soliton inside the coupled PhC waveguides, when a continuous signal with the peak power $\chi ^{(3)}E_0^{2} = 0.1410$ is injected into the guiding region ${(b)}$. The operating frequency is chosen $\frac {h\omega }{2\pi c}=0.232$ (Fig. 1). (b) Magnetic field $H_z$ versus the non-dimensional parameter $x/h$ at $y=0$ [Fig. 3(a)] where the blue line shows the numerical results based on FDTD analysis while the red dashed line stands for the theoretical result using (1) together with (6)–(11) and (16), (17). (c) Magnetic field $H_z$ versus $y/h$ at $x=0$ [Fig. 3(a)]. The blue line represents the numerical results based on FDTD analysis and the red dashed line shows the theoretical analysis based on (1) together with (6)–(11). Here [a.u] stands for the arbitrary units.

Download Full Size | PDF

Figure 3(b) illustrates the longitudinal dependence of the magnetic field $H_z$ versus the dimensionless parameter $x/h$ at $y=0$ [cf. Figure 3(a)], where [a.u] stands for arbitrary units. Numerical results from the FDTD analysis are shown by the blue line, whereas the theoretical result is indicated by a dashed red line. The latter we find by substituting (16) and (17) together with (6)–(11) into (1) and (2) from which the magnetic field $H_z$ is easily retrieved. Calculating the nonlinear coefficient $\gamma$ from (14), the width of the bandgap soliton $\Lambda$ is directly obtained from (17). Based on our theoretical analysis, the width of the bandgap soliton amounts to $\Lambda =3.95 h$, which is in a very good agreement with the numerical result ($\Lambda =4.10 h$). It is worth emphasizing that only the scattering amplitudes of the space harmonics $m = - 1$ and $m = 0$ gives rise to the formation of the electromagnetic field (1), while the scattering amplitudes of other space-harmonics are very small and can be neglected. All scattering amplitudes are calculated by solving for the eigenvalue problem in the linear regime defined as ${\mathbf {\Omega }} (\omega , k_{x0}, n_s) \cdot {\mathbf {A}}^{T} = 0$.

Figure 3(c) depicts the the transversal dependence of the magnetic field $H_z$ versus $y/h$ at $x=0$ [cf. Figure 3(a)]. The numerical results of the FDTD calculations are depicted as blue line, whereas the theoretical results are indicated by the dashed red line. The latter is obtained by substituting (6)–(11) into (1) from which the magnetic field $H_z$ is inferred. A very good agreement is observed between the results obtained from our theory and those that follow from the corresponding full-wave computational electromagnetics simulations.

4. Conclusion

A rigorous theoretical formulation describing the evolution of the envelope of temporal solitons propagating in Kerr-type nonlinear PhC waveguides has been proposed. The formalism is based on the full-wave modal analysis and is very general in nature and can be applied to various configurations of planar PhCs including the most challenging ones, such as plasmonic crystals with intrinsic losses. We investigated an experimentally feasible hexagonal lattice formed by air-holes in crystalline silicon. The expression for the nonlinear coefficient has been rigorously derived while accounting for all the space-harmonics and the interactions between them. The width of the soliton and the shift of the angular frequency due to Kerr nonlinearities have been analytically calculated. Extensive computational electromagnetics simulations based on the FDTD were performed demonstrating the correctness of the proposed formalism.

Appendix A

In case of a multilayered periodic structure such as a PhC (Fig. 1), the scattered space harmonics are conceptualized as new incident waves on the neighbouring PhC arrays and thus scattered into another set of space harmonics, which then impinges back on the original array. The scattering process from each layer is described by reflection and transmission matrices, which relate a set of the incident space harmonics to a set of reflected and transmitted ones. Reflection ${\mathbf {r}_{j}}$ and transmission ${\mathbf{f}_{j}}$ matrices for ${j}$-th layer are derived as follows [22]:

$${\mathbf{r}_{j}}=\mathbf{U}^{+}(k_{x0})[\mathbf{I} - \mathbf{T}(k_s)\mathbf{L}(k_{x0}h,k_sh)]^{{-}1} \mathbf{T}(k_s)\mathbf{P}(k_{x0})$$
$$\mathbf{f}_{j}={\mathbf{I}} + \mathbf{U}^{-}(k_{x0})[\mathbf{I} - \mathbf{T}(k_s)\mathbf{L}(k_{x0}h,k_sh)]^{{-}1} \mathbf{T}(k_s)\mathbf{P}(k_{x0})$$
with
$$\mathbf{P}(k_{x0})=\Biggl[i^{s} e^{i s \cos^{{-}1}(k_{xq}/k_s)}\Biggr]$$
$$\mathbf{U}^{{\pm}}(k_{x0})=\Biggl[\frac{2({-}i)^{s}}{k_{yn}h}e^{{\pm} i s \cos^{{-}1}(k_{xn}/k_s)}\Biggr]$$
where $\mathbf {P}(k_{x0})$ is a matrix that transforms the down-going ${q}$-th incident space harmonic wave to the ${s}$-th scattered cylindrical harmonic wave, $\mathbf {U}^{+}$ and ${\mathbf {U}}^{-}$ are the matrices that transform the ${s}$-th scattered cylindrical harmonic wave back to the up-going and down-going ${n}$-th space-harmonic waves (with $n, s, q = -M, -M+1, \ldots , 0, \ldots , M-1, M$), $\mathbf {T}(k_s)$ is the T-matrix describing the nature of the circular scatterer per unit cell, ${\mathbf {I}}$ is unit matrix and $\mathbf {L}(k_{x0}h, k_sh)$ is the matrix called the lattice sums (LSs). All the matrices are given as functions of the complex wavenumber $k_{x0}$, whereas the T-matrix does not depend on $k_{x0}$, but depends only on the geometrical and material parameters of the circular scatterer per unit cell and the wavenumber $k_s$ characterizing the background medium. The LSs whose elements are $L_m$, with $m = q - s$, are defined as [21,22]:
$$L_m(k_{x0}h,k_sh) = \sum_{n=1}^{\infty} H^{(1)}_m (k_snh) \Biggl[ e^{ink_{x0}} +e^{i\pi m-ink_{x0}h}\Biggr]$$
where $m= 0,1,2,\ldots$ and $H_m^{(1)}$ is the ${m}$-th order Hankel function of the first kind. The LSs only characterize the periodic arrangements of the scatterers. Recently we have used the formulation of the Ewald method for a fast and accurate calculation of the LSs applicable for real and complex-valued phase shifts $k_{x0}$ [21]. The latter is needed when analyzing complex guided and leaky waves. Since the reflection and transmission matrices for each array have been rigorously defined based on (18) and (19) together with (20) and (21), the generalized reflection matrix ${\mathbf {R}} (\omega , k_{x0}, n_s)$ for the $N$-layered periodic structure array can be easily obtained using a recursive algorithm [22,23].

Appendix B

A schematic view of a guiding region ${(a)}$ having a refractive index $\tilde {n}_s=n_s + \delta n$, whereas the refractive index of other regions is equal to $\ n_s$ is depicted in Fig. 2. In this case we need a slight modification only of (6) and (7) by implementing the diagonal Fresnel matrices. Other expressions are the same as in (8)–(11). Expressions in (6) and (7) are re-written in the following form:

$$\mathbf{\mathbf{a^{-}}} = \mathbf{W}_1 \mathbf{F}_{\tilde{n}_s,n_s} \mathbf{W}_2 \mathbf{F}_{n_s,\tilde{n}_s} \mathbf{W}_1 \mathbf{R}\cdot \mathbf{\mathbf{a^{+}}}$$
$$\mathbf{\mathbf{a^{+}}} = \mathbf{W}_1 \mathbf{F}_{\tilde{n}_s,n_s}\mathbf{W}_2 \mathbf{F}_{n_s,\tilde{n}_s}\mathbf{W}_1 \mathbf{R}_B \cdot \mathbf{\mathbf{a^{-}}}+\mathbf{W}_1 \mathbf{F}_{\tilde{n}_s,n_s}\mathbf{W}_2 \mathbf{F}_{n_s,\tilde{n}_s} \mathbf{W}_1 \mathbf{T}_B \cdot \mathbf{\mathbf{b^{+}}}$$
where ${\mathbf {W}}_1=[\exp (ik_{yn}r)]$ and ${\mathbf {W}}_2=\{\exp (i\tilde {k}_{yn}[w-2r])\}$ are the diagonal matrices which represent the phase shift of the up-going and down-going space-harmonics, $k_{yn}$ and $\tilde {k}_{yn}$ are the wavenumbers along the $y$-axis of the $n$-th space harmonic in media with the refractive indices $\ n_s$ and $\tilde {n}_s$, respectively; ${\mathbf {F}}_{\tilde {n}_s, n_s}$ denotes the Fresnel matrix, which connects the transmitted space harmonics to the medium with a refractive index $n_s$ and the space-harmonics impinging from the medium with a refractive index $\tilde {n}_s$, ${\mathbf {F}}_{n_s,\tilde {n}_s}$ stands for the Fresnel matrix, which connects the transmitted space harmonics to the medium with a refractive index $\tilde {n}_s$ and the space-harmonics impinging from the medium having a refractive index $n_s$. It can be checked that when the refractive indices of both media are the same, the Fresnel diagonal matrices become unit matrices and (22) and (23) are reduced to (6) and (7). The calculation procedure for the other guiding regions with different refractive indices is the same.

Acknowledgments

V. Jandieri would like to express his sincere gratitude to Professor Kiyotoshi Yasumoto from Kyushu University, Japan for fruitful discussions about the theory of photonic crystals.

References

1. M. Remoissenet, Waves Called Solitons: Concepts and Experiments, (Springer, 1996).

2. T. Dauxois and M. Peyrard, Physics of Solitons, (Cambridge University, 2010).

3. Y. Kivshar and G. Agrawal, Optical Solitons: From Fibers to Photonic Crystals, (Academic Press, 2003).

4. A. Scott, Encyclopedia of Nonlinear Science, (Taylor & Francis Group, 2005).

5. G. I. Stegeman and M. Segev, “Optical spatial solitons and their interactions: universality and diversity,” Science 286(5444), 1518–1523 (1999). [CrossRef]  

6. A. Hasegawa and Y. Kodama, “Guiding-center soliton in optical fibers,” Opt. Lett. 15(24), 1443–1445 (1990). [CrossRef]  

7. P. St. J. Russell, “Photonic crystal fibers,” Science 299(5605), 358–362 (2003). [CrossRef]  

8. H. Gersen, T. J. Karle, R. J. P. Engelen, W. Bogaerts, J. P. Korterik, N. F. van Hulst, T. F. Krauss, and L. Kuiper, “Real-Space observation of ultraslow light in photonic crystal waveguides,” Phys. Rev. Lett. 94(7), 073903 (2005). [CrossRef]  

9. B. J. Eggleton, R. E. Slusher, C. M. de Sterke, P. A. Krug, and J. E. Sipe, “Bragg grating solitons,” Phys. Rev. Lett. 76(10), 1627–1630 (1996). [CrossRef]  

10. N. A. R. Bhat and J. E. Sipe, “Optical pulse propagation in nonlinear photonic crystals,” Phys. Rev. E 64(5), 056604 (2001). [CrossRef]  

11. T. Taniuti and N. Yajima, “Perturbation method for a nonlinear wave modulation,” J. Math. Phys. 10(8), 1369–1372 (1969). [CrossRef]  

12. M. Malishava and R. Khomeriki, “All-Phononic digital transistor on the basis of gap-soliton dynamics in an anharmonic oscillator ladder,” Phys. Rev. Lett. 115(10), 104301 (2015). [CrossRef]  

13. A. Blanco-Redondo, C. Husko, D. Eades, Y. Zhang, J. Li, T. F. Krauss, and B. J. Eggleton, “Observation of soliton compression in silicon photonic crystals,” Nat. Commun. 5(1), 3160 (2014). [CrossRef]  

14. B. B. Baizakov, B. A. Malomed, and M. Salerno, “Multidimensional solitons in periodic potentials,” EPL 63(5), 642–648 (2003). [CrossRef]  

15. A. Trombettoni and A. Smerzi, “Discrete solitons and breathers with dilute Bose-Einstein condensates,” Phys. Rev. Lett. 86(11), 2353–2356 (2001). [CrossRef]  

16. F. Geniet and J. Leon, “Energy transmission in the forbidden band gap of a nonlinear chain,” Phys. Rev. Lett. 89(13), 134102 (2002). [CrossRef]  

17. R. Khomeriki, “Nonlinear bandgap transmission in optical waveguide arrays,” Phys. Rev. Lett. 92(6), 063905 (2004). [CrossRef]  

18. V. Jandieri, R. Khomeriki, D. Erni, and W. C. Chew, “Realization of all-optical digital amplification in coupled nonlinear photonic crystal waveguides,” Prog. Electromagn. Res. 158, 63–72 (2017). [CrossRef]  

19. V. Jandieri, R. Khomeriki, and D. Erni, “Realization of true all-optical AND logic gate based on the nonlinear coupled air-hole type photonic crystal waveguide,” Opt. Express 26(16), 19845–19853 (2018). [CrossRef]  

20. J. Leuthold, C. Koos, and W. Freude, “Nonlinear silicon photonics,” Nat. Photonics 4(8), 535–544 (2010). [CrossRef]  

21. V. Jandieri, P. Baccarelli, G. Valerio, and G. Schettini, “1-D periodic lattice sums for complex and leaky waves in 2-D structures using higher-order Ewald formulation,” IEEE Trans. Antennas Propag. 67(4), 2364–2378 (2019). [CrossRef]  

22. K. Yasumoto, H. Toyama, and R. Kushta, “Accurate analysis of two-dimensional electromagnetic scattering from multilayered periodic arrays of circular cylinders using lattice sums technique,” IEEE Trans. Antennas Propag. 52(10), 2603–2611 (2004). [CrossRef]  

23. V. Jandieri and K. Yasumoto, “Electromagnetic scattering by layered cylindrical arrays of circular rods,” IEEE Trans. Antennas Propag. 59(6), 2437–2441 (2011). [CrossRef]  

24. J. W. Boyle, S. A. Nikitov, A. D. Boardman, J. G. Booth, and K. Booth, “Nonlinear self-channeling and beam shaping of magnetostatic waves in ferromagnetic films,” Phys. Rev. B 53(18), 12173–12181 (1996). [CrossRef]  

25. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, (Artech House, 1995).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (3)
Fig. 1.
Fig. 1. (Left) Schematic view of the three symmetric PhC waveguides with guiding regions ${(a)}$, ${(b)}$ and ${(c)}$ having the same width ${w}$ and separated by the barrier layers of two PhCs with number of layers $N_B$. Planar PhCs are composed of hexagonal lattice of the circular air holes periodically distributed along $x$-axis with a period $h$. The air-holes are infinitely long along the ${z}$-axis. The radius of the air holes is $r$, $n_s$ is the refractive index of the background medium. Here ${{a^{+}}}$ and ${{a^{-}}}$, ${{b^{+}}}$ and ${{b^{-}}}$, ${{c^{+}}}$ and ${{c^{-}}}$ define the up-going and down-going space-harmonics in the guiding regions ${(a)}$, ${(b)}$ and ${(c)}$, respectively. (Right) Dispersion curves of the symmetric (blue curve and green curve) and the antisymmetric (red curve) super-modes for the $H$-polarized field. The operating frequency for the realization of functional all-optical logic gates is $\frac {h\omega }{2\pi c}=0.232$, and it is marked by a red dot. The distributions of the magnetic field $H_z$ for the super-modes, as well as the geometry of the setup are shown in the corresponding insets.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. Schematic view of PhC waveguide with a guiding region ${(a)}$ having a refractive index $\tilde {n}_s=n_s + \delta n$, whereas the refractive index of the other regions is the same as in Fig. 1 and it is equal to $n_s$.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. (a) Schematic distribution of the magnetic field of the gap soliton inside the coupled PhC waveguides, when a continuous signal with the peak power $\chi ^{(3)}E_0^{2} = 0.1410$ is injected into the guiding region ${(b)}$. The operating frequency is chosen $\frac {h\omega }{2\pi c}=0.232$ (Fig. 1). (b) Magnetic field $H_z$ versus the non-dimensional parameter $x/h$ at $y=0$ [Fig. 3(a)] where the blue line shows the numerical results based on FDTD analysis while the red dashed line stands for the theoretical result using (1) together with (6)–(11) and (16), (17). (c) Magnetic field $H_z$ versus $y/h$ at $x=0$ [Fig. 3(a)]. The blue line represents the numerical results based on FDTD analysis and the red dashed line shows the theoretical analysis based on (1) together with (6)–(11). Here [a.u] stands for the arbitrary units.

View in Article | Download Full Size | PDF

Equations (25)

Equations on this page are rendered with MathJax. Learn more.

E x = Ψ 2 ∑ m ( u m + e i k y m y + u m − e − i k y m y ) e i β m x − i ω t + c . c .
E x m ± = Ψ u m ± ; E y m ± = ∓ Ψ β m k y m u m ±
D x = E x [ χ ( 1 ) + χ ( 3 ) ( E x 2 + E y 2 ) ]
D x m ± = Ψ u m ± [ n s 2 + χ ( 3 ) | Ψ | 2 ( 3 4 [ 1 + β m 2 k y m 2 ] | u m ± | 2 + 1 2 [ 3 − β m 2 k y m 2 ] | u m ∓ | 2 + ∑ μ ≠ m 1 2 [ 3 + β μ 2 k y μ 2 ] ( | u μ ± | 2 + | u μ ∓ | 2 ) + ∑ μ ≠ m β m β μ k y m k y μ ( | u μ ± | 2 − | u μ ∓ | 2 ) ) ]
n m , ν ± = n s + χ ( 3 ) | Ψ | 2 4 n s [ 3 2 [ 1 + β m 2 k y m 2 ] | u m , ν ± | 2 + [ 3 − β m 2 k y m 2 ] | u m , ν ∓ | 2 + ∑ μ ≠ m [ 3 + β μ 2 k y μ 2 ] ( | u μ , ν ± | 2 + | u μ , ν ∓ | 2 ) ]
a − = W ( ω ′ , k x 0 , n m , a ± ) R ( ω ′ , k x 0 , n s ) ⋅ a +
a + = W ( ω ′ , k x 0 , n m , a ± ) [ R B ( ω ′ , k x 0 , n s ) ⋅ a − + T B ( ω ′ , k x 0 , n s ) ⋅ b + ]
b + = W ( ω ′ , k x 0 , n m , b ± ) [ T B ( ω ′ , k x 0 , n s ) ⋅ c + + R B ( ω ′ , k x 0 , n s ) ⋅ b − ]
b − = W ( ω ′ , k x 0 , n m , b ± ) [ T B ( ω ′ , k x 0 , n s ) ⋅ a − + R B ( ω ′ , k x 0 , n s ) ⋅ b + ]
c + = W ( ω ′ , k x 0 , n m , c ± ) R ( ω ′ , k x 0 , n s ) ⋅ c −
c − = W ( ω ′ , k x 0 , n m , c ± ) [ T B ( ω ′ , k x 0 , n s ) ⋅ b − + R B ( ω ′ , k x 0 , n s ) ⋅ c + ]
W ( ω ′ , k x 0 , n m , ν ± ) = [ e i k y m w ] , ν = a , b , c
Ω ( ω ′ , k x 0 , n m , ν ± ) ⋅ A T ( ω ′ , k x 0 , n m , ν ± ) = 0
δ ω A ′ ∂ Ω ( ω , k x 0 , n s ) ∂ ω A T = − | Ψ | 2 A ′ ∂ Ω ( ω , k x 0 , n m , ν ± ) ∂ ( | Ψ | 2 ) A T | n m , ν ± = n s
∂ Ω ( ω , k x 0 , n m , ν ± ) ∂ ( | Ψ | 2 ) = ∑ m , ν ∂ ( n m , ν ± ) ∂ ( | Ψ | 2 ) ∂ Ω ( ω , k x 0 , n m , ν ± ) ∂ ( n m , ν ± ) | n m , ν ± = n s
i ( ∂ Ψ ∂ t + v g ∂ Ψ ∂ x ) + ω ″ 2 ∂ 2 Ψ ∂ x 2 + γ | Ψ | 2 Ψ = 0
Ψ ( x , t ) = F e − i ( Δ ω ) t cosh ⁡ [ ( x − v g t ) / Λ ]
Λ = ω ″ γ F 2 , Δ ω = − γ F 2 2 = F 2 2 ∂ ( δ ω ) ∂ ( | Ψ | 2 )
r j = U + ( k x 0 ) [ I − T ( k s ) L ( k x 0 h , k s h ) ] − 1 T ( k s ) P ( k x 0 )
f j = I + U − ( k x 0 ) [ I − T ( k s ) L ( k x 0 h , k s h ) ] − 1 T ( k s ) P ( k x 0 )
P ( k x 0 ) = [ i s e i s cos − 1 ⁡ ( k x q / k s ) ]
U ± ( k x 0 ) = [ 2 ( − i ) s k y n h e ± i s cos − 1 ⁡ ( k x n / k s ) ]
L m ( k x 0 h , k s h ) = ∑ n = 1 ∞ H m ( 1 ) ( k s n h ) [ e i n k x 0 + e i π m − i n k x 0 h ]
a − = W 1 F n ~ s , n s W 2 F n s , n ~ s W 1 R ⋅ a +
a + = W 1 F n ~ s , n s W 2 F n s , n ~ s W 1 R B ⋅ a − + W 1 F n ~ s , n s W 2 F n s , n ~ s W 1 T B ⋅ b +
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use