Study of plasmonic slot waveguides with a nonlinear metamaterial core: semi-analytical and numerical methods

# Study of plasmonic slot waveguides with a nonlinear metamaterial core: semi-analytical and numerical methods

Mahmoud M. R. Elsawy    Gilles Renversez Aix–Marseille Univ, CNRS, Ecole Centrale Marseille, Institut Fresnel, 13013 Marseille, France
July 5, 2019
###### Abstract

Two distinct models are developed to investigate the transverse magnetic stationary solutions propagating in one-dimensional anisotropic nonlinear plasmonic structures made from a nonlinear metamaterial core of Kerr-type embedded between two semi-infinite metal claddings. The first model is semi-analytical, in which we assumed that the anisotropic nonlinearity depends only on the transverse component of the electric field and that the nonlinear refractive index modification is small compared to the linear one. This method allows us to derive analytically the field profiles and the nonlinear dispersion relations in terms of the Jacobi elliptical functions. The second model is fully numerical, it is based on the finite-element method in which all the components of the electric field are considered in the Kerr-type nonlinearity with no presumptions on the nonlinear refractive index change. Our finite-element based model is valid beyond weak nonlinearity regime and generalize the well-known single-component fixed power algorithm that is usually used. Examples of the main cases are investigated including ones with strong spatial nonlinear effects at low powers. Loss issues are reduced through the use of gain medium in the nonlinear metamaterial core.

Nonlinear waveguides, optical, Optical solitons, Kerr effect: nonlinear optics, Plasmons on surfaces and interfaces / surface plasmons
###### pacs:
42.65.Wi, 42.65.Tg, 42.65.Hw, 73.20.Mf

## I Introduction

Nonlinear optical properties play a crucial role in all-optical integrated circuits due to the different control functionalities they offer Tien:71-integrated-optics (); stegeman1985nonlinear-integrated-optics (); stegeman1988third-nl-integrated (). Utilizing plasmonics as a part of nonlinear structures may be a promising choice because of the reduced footprint achievable compared with all-dielectric structures, and because of the enhancement of the field intensities, which can be used to boost the nonlinearity kauranen_nonlinear_2012 (); rukhlenko_nonlinear_2011 (). Several nonlinear plasmonic waveguides have already been studied Ariyasu85 (); Agranovich80 (); Walasik14 (); ajith2015linear-Metal-NL (). The structures composed of a nonlinear isotropic core of Kerr-type sandwiched between two semi-infinite metal claddings, have received a great attention since their study in 2007 Feigenbaum07 (), due to the strong light confinement obtained and due to its peculiar nonlinear effects Davoyan08 (); rukhlenko2011exact-dispersion (); Ferrando13 (); salguiero14complex-modes-plasmonic-nonlinear-slot-waveguides (); Walasik15b (). These nonlinear plasmonic slot waveguides (NPSWs) promise a family of exciting applications such as phase matching in higher harmonic generation processes Davoyan09a (), nonlinear plasmonic couplers sammut_theoretical_1993 () or switching Nozhat12 (). In order to study light propagating in such structures, different methods have been developed Davoyan08 (); rukhlenko2011dispersion (); rukhlenko2011exact-dispersion () to describe its main modes. Recently, a full description of the solutions including the higher-order ones was introduced Walasik15a (); Walasik15b (). In addition, an improvement of such structures by the inclusion of two linear dielectric buffer layers between the nonlinear core and the two metal claddings has also been proposed Elsawy_OL_16 (); Elsawyspie15 (). Nevertheless, all these studies deal only with the standard nonlinear isotropic core with focusing Kerr-type. It appears that the power needed to observe the symmetry breaking of the fundamental symmetric mode is in the range of GW/m which is still too high. Otherwise, during the last few years, it was demonstrated that the nonlinear effects can be extremely enhanced using epsilon-near-zero (ENZ) materials ciattoni2010extreme (); campione2013electric-field-enhancement-ENZ (); ciattoni2016enhanced (); neira2015eliminating (); alam2016large-indium-tin-oxide-ENZ (). However, even if recent experimental results have demonstrated a large enhancement of third harmonic generation in the epsilon-near-zero (ENZ) regime Capretti15OL-enhanced-third-harmonic-generation (), the usual modelling approaches are no longer applicable in this high nonlinear regime since the spatial nonlinearity cannot be treated as a perturbation. We have shown in a recent work that in order to fully take advantage of the ENZ nonlinearity enhancement in nonlinear waveguides, it is important to introduce anisotropy Elsawy16-OL-anisotropic-metamaterial-waveguides (). Even if, plasmonic waveguides have already been studied with metamaterial layers either in the cladding ishii2014plasmonic-HMM-cladding () or in the core avrutsky2007-HMM-core-linear (); rukhlenko2012guided-plasmonic-anisotropic (), these last studies focused only on linear guided waves. In ciattoni2010extreme (), nonlinear guided waves were studied in anisotropic media, but with an effective nonlinear response identical to the standard isotropic one for transverse magnetic (TM) polarized waves. In our study, we consider a nonlinear metamaterial with an anisotropic effective response for the TM polarized waves. To allow this approach, we generalize two of our previous methods to tackle anisotropic nonlinear cores and not only isotropic ones. For the first method, this implies both a new classification of the possible cases for the effective nonlinearity and more complicated nonlinear dispersion relations compared with the isotropic case. The second method is based on the finite element method and the fixed power algorithm where the nonlinear stationary solutions are computed numerically as a function of the fixed input power. Usually, it is assumed that only the transverse component of the electric field is taken into account in the nonlinear form. Here, this is not the case, since we take into account all the electric field components in the Kerr-type nonlinearity. In order to achieve this result, we introduce and solve two coupled scalar nonlinear equations, on per the continuous tangential component of the electromagnetic field.

In this study, we present the full derivation of two different methods to study nonlinear TM stationary solutions propagating in one-dimensional plasmonic slot waveguides with an anisotropic metamaterial nonlinear core. The first one is semi-analytical while the second model is fully numerical and is based on the finite-element method (FEM).

These models have been briefly introduced in Elsawy16-OL-anisotropic-metamaterial-waveguides (), while here we present their detailed derivations. The article is organized as follows, in section II the statement of the problem is presented. The semi-analytical method and the numerical FEM are described in section III. In section IV the validation of our models is given by adapting them to study the isotropic case and the comparison with the results from previously published work is described. In addition, we present several examples to show the influence of the anisotropic nonlinearity on the nonlinear dispersion curves and on the field profiles.

## Ii Problem formalism

We consider a structured metamaterial nonlinear core, formed by bulk layers, embedded between two semi-infinite isotropic metal claddings (see figure 1(a)). The core layers consist of two different isotropic media with two different permittivities , and two different thicknesses , , respectively as shown in figure 1(b). We propose two different models to study light propagation in such anisotropic nonlinear structures. The first model is based on the approach presented in Walasik15a (); Walasik15b () for isotropic nonlinear core, extending it to structures containing an anisotropic nonlinear core. In this method, the nonlinearity is treated in an approximated way such that all the components of the permittivity tensor depend only on the transverse component of the electric field. Another limitation is that the nonlinear refractive index change is small compared to the linear refractive index. Based on the above assumptions, analytical formulas for the field profiles in terms of the Jacobi elliptical functions Abramowitz-Handbook-Mathematical-functions (); Byrd54-Handbook-elliptic-integrals () are obtained; with the continuity of the tangential electromagnetic field components, they allow us to derive analytical formulas for the nonlinear dispersion relations. This method will be called the extended Jacobi elliptical model (EJEM). Our EJEM, is different from the simple model used in references Ariyasu85 (); Stegeman85 (), in which only two components of the permittivity tensor depend on the transverse component of the electric field, while in our case we consider the dependence of the transverse component of the electric field for all the components of the permittivity tensor. In addition, the methods used in Ariyasu85 (); Stegeman85 (), developed to study stationary solutions propagating in structures containing a semi-infinite nonlinear region, whereas in our case we consider a finite-size nonlinear anisotropic core. The second model described in the current study is fully numerical and based on the FEM to solve the stationary TM problem in nonlinear layered structures. This numerical FEM does not require any of the assumptions used in the semi-analytical EJEM, and all the components of the effective nonlinear permittivity tensor depend on both the transverse and the longitudinal components of the electric field, moreover the nonlinear term does not need to be small, which means that it is valid beyond the weak nonlinearity regime. In order to treat the nonlinearity in the FEM, we generalize the fixed power algorithm presented in Drouart08 (); Walasik14 () and we consider a coupled nonlinear eigenvalue problem to take into account all the electric field components in the Kerr-type nonlinearity instead of the single scalar eigenvalue problem solved previously. Furthermore, in our new model, the nonlinearity is treated without any assumptions relative to its amplitude. It is worth mentioning that the first semi-analytical model (EJEM), provides more insight and understanding into the nature of the of finding stationary solutions in the structure than the second, more numerical model. Nevertheless, the second model treats the nonlinearity in a proper way without any assumptions on the Kerr-nonlinearity amplitude, but the field profiles are computed numerically.

Our models are written for TM light polarization in which the magnetic field has only one component such that and the electric field has two components . We consider only monochromatic TM waves propagating along the -direction in one-dimensional symmetric structures which are invariant along and directions (see figure 1(a)). In these structures, all the field components evolve proportionally to : ,

where represents the wave number in vacuum, denotes the angular frequency of the wave, and c is the speed of light in vacuum. Here and denote the effective index and the propagation constant of the wave, respectively. The magnetic field induction vector is defined as , where the permeability is constant and is denoted by , that of vacuum. In our models, we assume that the relative permittivity tensor is complex and diagonal: . We consider only the real part of the permittivity in the derivation of the nonlinear dispersion relations, such that the displacement vector is defined as , with being the vacuum permittivity, the imaginary parts will be used to estimate the losses, in which we use the same method as in Davoyan09 (); Walasik14 (); Stegeman85 (), but extended to the anisotropic case (see section IV). Depending on the chosen orientation of the compound layers relative to the Cartesian coordinate axes, different anisotropic permittivity tensors can be build for the core. Due to the -invariance hypothesis required by our modelling, only two types where the -axis belongs to the layers have to be considered. For the first one where the layers are parallel to the -axis (figure 1(b)), one has for the diagonal terms of : . For the second case where the layers are parallel to the -axis, one gets: (see figure 1(c)). The first one allows us to obtain an anisotropic effective nonlinear response for TM polarized waves, unlike the second case in which the effective nonlinear response coincides with the standard isotropic one. We will focus on the first orientation of the layers as depicted in figure 1(b). The nonlinearity considered in this study is the usual Kerr-type, where all the components of the relative permittivity tensor depend on the electric field intensity:

 ϵj=ϵjj+αjj|E(x)|2\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ∀j∈{x,y,z}, (1)

in which is the j-th component of the real linear permittivity tensor and is the corresponding nonlinear parameter, . This kind of nonlinearity has already been used in many nonlinear waveguide studies Feigenbaum07 (); Agranovich80 (); Maradudin88 (); Walasik15a (), as a first step we do not need to consider the complex Kerr-type nonlinearity used in ciattoni2010extreme (); rizza2011two-peaked-flat-top ().

Finally, using the definition of the magnetic field induction , the displacement vector we can express Maxwell’s equations for the TM polarized waves as:

 Ex(x)=neffHy(x)ϵ0ϵx(x)c, (2a) Ez(x)=1ϵ0ϵz(x)ωdHy(x)dx, (2b) k0neffEx(x)−dEz(x)dx=ωμ0Hy(x). (2c)

## Iii Derivation of the Methods

### iii.1 Extended Jacobi elliptical model (EJEM)

We begin with the derivation of the field profiles and the nonlinear dispersion relation in the frame of the EJEM, in which strong assumptions on the form of the Kerr-type nonlinearity are required to establish it. This method is a generalization and an extension to the Jacobi elliptical model (JEM) developed to study the stationary nonlinear solutions in isotropic plasmonic slot waveguides Walasik15a (); Walasik15b (), and based on the same assumptions: (i) the nonlinearity depends only on the transverse component of the electric field and (ii) the nonlinear permittivity modifications are small compared to the linear refractive index such that :

 ϵj=ϵjj+αjj|Ex|2, (3) αjj|Ex|2<<ϵjj.

These assumptions are valid only at low power as it was shown in Walasik15a (); Walasik14 (); Maradudin88 (); Stegeman84 (); Agranovich80 (). However, they allow us to derive analytical formulas for the field profiles inside the anisotropic nonlinear core in terms of the Jacobi elliptical functions Abramowitz-Handbook-Mathematical-functions (); Byrd54-Handbook-elliptic-integrals () which will be used together with the field in the linear metal claddings and the continuity of the tangential components at the core interface to acquire analytical expressions for the nonlinear dispersion relations.

In order to derive the nonlinear wave equation in terms of the magnetic field component in the frame of the EJEM assumptions, we use equations 2 and proceeding as in the isotropic case Walasik14 (), we find

 d2Hydx2−k20q2(x)Hy(x)+k20a\tiny{% EJEM}nl(x)H3y(x)=0, (4)

where

 q2(x)=⎧⎪⎨⎪⎩q2core=([Re(neff)]2ϵzzϵxx−ϵzz)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak $inthecore$,q2m=([Re(neff)]2−ϵm)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak $inthecladdings$, (5)

in which being the real part of the permittivity in the metal claddings. In equation equation (5), we consider only the real part of the effective index as used in Walasik15a (); Ariyasu85 (); Stegeman85 (). The nonlinear coefficient is null in the linear metal claddings while in the nonlinear anisotropic core it is given by:

 a\tiny{EJEM}nl=−[Re(neff)]2ϵ4xxc2ϵ20([Re(neff)]2(αzzϵxx−αxxϵzz)−αzzϵ2xx). (6)

Our model can be reduced to the isotropic state Walasik14 (); Walasik15a (); Walasik15b () by setting and in equations (5) and (6) from which we recover the same expressions for and as in the JEM (the same isotropic response can be obtained with the orientation shown in figure 1(c)). It is worth mentioning that the nonlinear coefficient given by equation (6) is different from what has been developed for the anisotropic case presented in Ariyasu85 (); Stegeman85 (). In these two references (using to our notations to facilitate the comparison), it was assumed that and depend on the transverse component of the electric field, while is assumed not to depend on this component (it is assumed to be constant). In the present study, we consider the dependence of the transverse component of the electric field for all the components of the permittivity tensor. As a special case, in our recent work Elsawy16-OL-anisotropic-metamaterial-waveguides (), we presented the results considering isotropic nonlinear term () with an anisotropic linear one () in order to show the influence of the anisotropic linear part of the permittivity on the nonlinear dispersion diagrams and on the field profiles. In Elsawy16-OL-anisotropic-metamaterial-waveguides (), we distinguished between two different cases; the elliptical case, in which and , and the hyperbolic case, where and . While, in the present study, we present a general and full treatment for anisotropic linear and nonlinear permittivity terms with more general classifications based on the sign of , which depends on both the effective linear and nonlinear permittivity terms (see equation (6)). In order to solve equation (4), we use the first integral approach Maradudin88 (); Mihalache87 (); Walasik14 () and integrate in each of the structure layers separately. Consequently, we can write the nonlinear wave equation as:

 (dHydx)2−k20q2(x)H2y(x)+k20a\tiny{EJEM}nl(x)2H4y(x)=C0, (7)

where is a constant of integration. In the semi-infinite metal cladding, is null since the magnetic field and its derivative tends to zero as . Moreover, in the linear cladding, the nonlinear parameter and equation (4) reduces to the standard linear wave equation and the magnetic field matches the usual decaying exponential in the metal regions whose solutions are given by:

 \leavevmode\nobreak \leavevmode\nobreak Hy=H[−dcore/2]ek0qm(x+dcore2)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak −∞≤x<−dcore2, (8) Hy=H[dcore/2]e−k0qm(x−dcore2)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak dcore2≤x<+∞.

Here, and are the values of the magnetic field amplitudes at the left and the right interfaces, respectively. We are searching for guided wave solutions in the anisotropic waveguides depicted in figure 1, consequently, we will look only for the solutions with positive attenuation coefficient in the metal claddings (see equation (8)). While, the quantity in the core can be either real or imaginary leading to positive or negative values of (see equation (5)).

In the nonlinear core, and this constant can be obtained in terms of the magnetic field amplitudes at the core interfaces using the continuity of the longitudinal components at the core interface such that

 C0=k20H[−dcore/2]⎡⎣(ϵzzϵm)2q2m−q2core+a\tiny{% EJEM}nl2H2[−dcore/2]⎤⎦. (9)

A similar expression can be obtained for at the right interface. It is important to mention that the value of at the left interface is replaced by in equation (9) since we assumed, in the EJEM, that the nonlinear refractive index change is small compared to the linear one. The expression of is a generalization to what have been obtained in the isotropic case Walasik15a () with different values of and , such that we recover the same expression in the isotropic case by setting and . Due to the anisotropy, the sign of the nonlinear parameter could be positive or negative depending on the values of , , , and (see equation (6)). Therefore, we will consider a general classification of the nonlinear solutions in the core according to the sign of the nonlinear parameter .

#### iii.1.1 Case a\tiny{EJEM}nl<0

In this case, we set in equation (9) to write the integration constant in the nonlinear core as

 C0=k20H[−dcore/2]⎡⎣(ϵzzϵm)2q2m−q2core−|a\tiny{% EJEM}nl|2H2[−dcore/2]⎤⎦. (10)

Here, we cannot determine directly the sign of according to the sign of , hence the magnetic field will be classified according to the signs of and . Nevertheless, we have found that the only bounded solutions in the nonlinear core in case of a negative value of can be obtained only when and with the following criterion

 q4core≥2C0|a\tiny{EJEM}nl|k20. (11)

In order to clarify this point, we write equation (7) in the reduced from

where the reduced parameters are given by

 A=2/(k20|a\tiny{EJEM}nl|),\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak Q=k20|q2core|. (13)

Now, equation (12) can be written as

with

 γ2=(AQ+√(AQ)2−4AC0)/2, (15) δ2=(AQ−√(AQ)2−4AC0)/2,

where . We are looking for real values of the magnetic field , thus and must be real, which means that we need to consider the case , this gives the condition shown in equation (11). In order to ensure that the quantity under the square root in the denominator of the left-hand side of equation (14) is positive, we must have . It is worth mentioning that we have found an unbounded solution in the core when the condition in equation (11) is not satisfied (see appendix A for more details). Integrating  equation (14) in the nonlinear core and using formula 17.4.45 from Abramowitz-Handbook-Mathematical-functions (), we can express the magnetic field profile in the nonlinear core (for , , and ) in terms of the bounded Jacobi elliptical function sn with argument and parameter  Abramowitz-Handbook-Mathematical-functions (); Byrd54-Handbook-elliptic-integrals (); Salas14-duff-eq (); David-mathematica07 () as:

 Hy(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩δsn[−√γ2A(x+dcore2)+X0∣∣∣m]\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ϵzz>0,δsn[+√γ2A(x+dcore2)+X0∣∣∣m]\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ϵzz<0. (16)

Where,

 m=δ2γ2\leavevmode\nobreak \leavevmode\nobreak andX0=sn−1[H[−dcore/2]δ∣∣∣m]. (17)

Here, sn is the inverse Jacobi elliptical function sn Abramowitz-Handbook-Mathematical-functions (). In order to demonstrate the choice of the sign in  equation (16), we need to look at equation (12). If the condition given by equation (11) is satisfied, the denominator of equation (12) is positive and the sign in the right-hand side depends only on the infinitesimal change of the magnetic field . In order to be able to study the sign of the magnetic field derivative, we consider the continuity condition of the longitudinal component at the left core interface

 −k0qmϵmHy∣∣∣x=(−dcore/2)−=1ϵzz[dHydx]x=(−dcore/2)+. (18)

The real part of the permittivity in the metal cladding is negative , and we are looking for solutions with positive attenuation coefficient in the cladding ; additionally, we assume that . Now, since the magnetic field changes its sign at the interface ( being continuous) the derivative of the magnetic field inside the nonlinear core and near the interface depends on the sign of such that: for the sign of the magnetic field derivative is positive, while if the correct choice of the magnetic field derivative at the interface is the negative sign. Consequently, according to the sign of , we choose the sign in the right-hand side of equation (12). It is important to notice that, the solution shown in equation (16) for , cannot be obtained in the isotropic case with positive Kerr-nonlinearity. Nevertheless, due to the full anisotropic treatment of the effective nonlinearity in the present study, we can obtain a nonlinear solution with , starting from a positive Kerr-nonlinearity (see section IV for an example of this case), which can also be seen by looking at equation (6). This is one of the consequences of the anisotropic nonlinear core considered in the current study.

The procedure of the derivation of the nonlinear dispersion relation is similar to what have already been used for the isotropic case Walasik15a (); Walasik15b (): we use the magnetic field profile in the nonlinear core equation (16), the magnetic field in the linear metal claddings (equation (8)), and Maxwell’s equations (equations 2). Using the continuity conditions for the tangential electromagnetic field components and at the right core interface (), we can write the nonlinear dispersion relation in terms of the bounded Jacobi elliptical functions sn, cn, and dn with argument and parameter as

In equation (19), the sign in front of the square roots is related to the sign of as we discussed before (see equation (16) and the text after). We remind the reader that appears in equation (19) through that is used to define and consequently and . Since we are searching for nonlinear bounded solutions with finite energy in the core, we will consider only the subcase with and for negative (which satisfies the condition shown in equation (11)) in the derivation of the nonlinear dispersion relation and the other subcases which provide unbounded solutions in the nonlinear core will be summarized in appendix A and table 1.

#### iii.1.2 Case a\tiny{EJEM}nl>0

Now we consider the case in which the effective parameters , , , and provide positive values for the effective nonlinearity (see equation (6) ). The situation is quite similar to the isotropic case with focusing Kerr-type nonlinearity Walasik15a (); Walasik15b (), while the present anisotropic case is more general since it is not necessary to use focusing Kerr-type nonlinearity to get positive values of , as it can be inferred from equation (6). The isotropic case studied previously Walasik15a (), can be seen as a special case of the current study. For , the nonlinear wave equation given by equation (7) can be written as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩dHy√AC0+AQH2y−H4y(x)=±√1Adx\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak q2core>0,dHy√AC0−AQH2y−H4y(x)=±√1Adx\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak q2core<0, (20)

with the parameters and such that

 A=2(k20a\tiny{EJEM}nl),\leavevmode\nobreak and\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak Q={k20q2core\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak %for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak q2core>0,k20|q2core|\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak q2core<0. (21)

This means that equation (20) takes different forms according to the signs of the integration constant and the quantity .

For : in this subcase, the integration constant can only take positive value as it can be inferred from equation (9) for . In this case, the magnetic field profile will be written in terms of the bounded Jacobi elliptical function cn as

 Hy(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩δcn[+√γ2+δ2A(x+dcore2)+X0∣∣∣% m]\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak ϵzz>0,δcn[−√γ2+δ2A(x+dcore2)+X0∣∣∣m]\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak ϵzz<0, (22)

where,

 γ2=+AQ+√(AQ)2+4AC02, (23) δ2=−AQ+√(AQ)2+4AC02, m=δ2γ2+δ2,and\leavevmode\nobreak \leavevmode\nobreak X0=cn−1[H[−dcore/2]δ∣∣∣m].

The choice of the sign in front of the square roots in equation (22) is related to the sign of and it is treated as in the previous case (see equation (18) and the text after together with equation (20)). In order to derive the nonlinear dispersion relation, in this case, we use the magnetic field profile in the nonlinear core (equation (22)) and the field profile in the linear metal cladding (equation (8)) together with the continuity condition of the tangential components of the electromagnetic field at the right core interface. The final expression of the nonlinear dispersion relation reads

We choose the lower positive sign for and the upper negative sign for .

For : unlike the former subcase, the sign of the integration constant cannot be determined directly.
First, we consider . Again, the magnetic field profile in the nonlinear core can be obtained by integrating equation (20) using formula 17.4.43 from Abramowitz-Handbook-Mathematical-functions () with the reduced parameters and defined in equation (21) such that

 Hy(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩δnd[−√γ2A(x+dcore2)+X0∣∣∣m],\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak ϵzz>0,δnd[+√γ2A(x+dcore2)+X0∣∣∣m],\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak \leavevmode\nobreak ϵzz<0, (25)

where

 γ2=AQ+√(AQ)2−4A|C0|2, (26) δ2=AQ−√(AQ)2−4A|C0|2, m=γ2−δ2γ2,and\leavevmode\nobreak \leavevmode\nobreak X0=nd−1[H[−dcore/2]δ∣∣∣m].

It is worth mentioning that according to the parameters used, the quantities under the square roots in the definitions of and shown in equation (26) are positive which ensures that both and are real quantities and the magnetic field profile provided by equation (25) is real. Using the magnetic field profile shown in equation (25) and proceeding as for the previous case , we can write the nonlinear dispersion relation as

In equation (27) we choose the upper negative sign in front of the square roots for , and we choose the bottom positive sign for  (see equation (18), and the derivative of the Jacobi elliptical functions nd in Abramowitz-Handbook-Mathematical-functions ()).

Second, for , we can write the magnetic field profile in the nonlinear core as:

 Hy(x)= (28) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩γδ√γ2+δ2sd[−√γ2+δ2A(x+dcore2)+X0∣∣∣m]\leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak ϵzz>0,γδ√γ2+δ2sd[+√γ2+δ2A(x+dcore2)+X0∣∣∣m]\leavevmode\nobreak \leavevmode\nobreak for\leavevmode\nobreak ϵzz<0,

where

 γ2=−AQ+√(AQ)2+4AC02, (29) δ2=+AQ+√(AQ)2+4AC02,

The associated nonlinear dispersion relation gives

It is important noticing that an alternative formula to equation (28) can also be obtained in terms of the Jacobi elliptical function cn as it was already used in the isotropic case Walasik15a (). The choice of the sign of the nonlinear dispersion relation in equation (30) is linked to the sign of as it is shown in equation (28) in which we choose the upper negative sign for and the lower positive one for .

The nonlinear dispersion relations depicted in equations (19), (24), (27), and (30) represent the full nonlinear dispersion relations for bounded solutions propagating in the anisotropic nonlinear plasmonic slot waveguide depicted in figure 1 under the EJEM assumptions described in subsection III.1. For a fixed opt-geometrical parameters, the above nonlinear dispersion relations are solved for the real part of the effective index to obtain the nonlinear dispersion diagrams.

### iii.2 Finite-element method

In this part, we use the finite-element based model (FEM) to solve the nonlinear TM eigenvalue problem in the one-dimensional anisotropic layered structure depicted in figure 1. It is well-known that the FEM is generally versatile and can be applied to complex nonlinear waveguide problems, including the two-dimensional ones with arbitrary shape and field profiles Rahman90 (); li1992variational (); polstyanko1996nonlinear (); desevedavy2009Te-As-Se-optical-fiber (). Generally speaking, in the frame of the FEM, the initial physical problem is transformed into a variational form (weak formulation) by multiplying the initial partial differential equation by chosen test functions that belong to a particular functional space. The next step is the discretization of the problem in which the waveguide cross section is first divided into a patchwork of elements. The unknown fields are expanded in terms of interpolation polynomials over each element. The expansion coefficients that define the values of the fields at the nodal elements can then be obtained by solving a standard matrix eigenvalue problem. For a general and recent review of the finite-element method in the frame of optical waveguides, the reader can refer to chapter in reference livre12-FPCF ().

In this article, in order to treat the anisotropic nonlinearity in the frame of our FEM, the fixed power algorithm Rahman90 (); Ferrando03_fixed_power (); Drouart08 () will be used, in which the input is the total power and the outputs are the field profile (the eigenfunction) and the corresponding effective index (the eigenvalue). This algorithm uses a simple iterative scheme based on a sequence of linear modal solutions which converge to the nonlinear solution after few steps. The fundamental issue is that the amplitudes of the eigenmode are irrelevant and the numerical solution of the intermediary eigenvalue problem used at each iteration has an uncontrolled amplitude. But, since the nonlinear problem depends on the amplitude of the field, the numerical eigenvector which is computed at each step as a solution of the eigenvalue problem has to be scaled by a scalar factor. This process allows the eigenvector to be used as input for the next iteration. The scaling factor is computed at each iteration for a fixed value of the input power. For the single-component eigenvalue problem in terms of the magnetic field component , an approximated formula has already been used Drouart08 (); Walasik14 (); Elsawy_OL_16 (); Elsawyspie15 () in order to compute the scaling factor using only the linear part of the permittivity and neglecting the nonlinear term (see equation (38) in appendix B). Moreover, only the transverse component of the electric field was used in the isotropic Kerr-type nonlinearity. It is worth mentioning, that it is not possible to take into account all the electric field components in the Kerr-type nonlinearity using the single-component eigenvalue problem Walasik14 () within the FEM implementation of fixed power algorithm, as it is demonstrated in equation (39) and the text after it.

In the current study, we present a new and generalized formalism for the fixed power algorithm in arbitrary nonlinear layered structures where the nonlinearity will be treated in a more rigorous way such that all the components of the electric field are taken into account in the Kerr-type nonlinearity and no assumptions are needed to compute the scaling factor. Additionally, unlike the previous isotropic cases Rahman90 (); Drouart08 (); Walasik14 (); Elsawy_OL_16 (), we consider a fully anisotropic nonlinear treatment for the permittivity tensor as it is shown in equation (1).

Our approach is based on the solution of a coupled nonlinear eigenvalue problem in terms of the continuous tangential components of the electromagnetic field and . In order to obtain the correct weak formulation, we must consider the full TM wave equations with both the inhomogeneous permittivity term induced by the nonlinearity and the structure interfaces. The coupled weak formulation reads:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−1k20∫Γ1ϵz(x)dhydxdh′ydxdx+∫Γhy(x)h′y(x)dx−n2eff∫Γ1ϵx(x)hy(x)h′y(x)dx=0\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ,∫Γez(x)e′z(x)dx−1ϵ0k0c∫Γ1ϵz(x)dhydxe′zdx=0. (31)

Here, and stand for the tangential component of the magnetic and electric field, respectively, and is the domain of integration (in the present case the full cross section of the waveguide). In equation (31), we look for , where is Sobolev space of order with null Dirichlet boundary conditions on the domain of integration . It is important to point out that for the magnetic field , we need to pick a test function from a functional space (Sobolev space) at least of order , while for the tangential component of the electric field we can choose a test function from a functional space of order or as it can be understood from the relation between the tangential components shown in equation (2b). We use the fixed power algorithm described in algorithm 1 to solve the coupled eigenvalue problem shown in equation (31). Our FEM is implemented using the free and open-source software GMSH as a mesh generator and GETDP as a solver dular-getdp-1998general (); geuzaine2007getdp (); geuzaine2009gmsh (). These software programs have already been used to solve both one-dimensional and two-dimensional isotropic nonlinear electromagnetic waveguide problems Drouart08 (); Elsawy_OL_16 (); Walasik14 (). It is worth mentioning that our FEM with its fixed power algorithm shown in the algorithm 1 can be reduced to take into account all EJEM assumptions as described in appendix B.

## Iv Numerical results

In this section, we present several numerical results to validate our models and to demonstrate the influence of the anisotropy on the nonlinear dispersion diagrams and on the effective nonlinearity. We build nonlinear cores from stacks of realistic bulk materials, and we use the effective medium theory to retrieve their linear and the nonlinear effective parameters ciattoni2010extreme (); ciattoni2011all (). Moreover, we assume that both and (see figure 1(b)) are much smaller than the operating wavelength, so we can derive the effective linear permittivities and the effective nonlinear susceptibilities in the frame of the simple effective medium theory (EMT) ciattoni2010extreme (); ciattoni2011all ()

 ~ϵzz=rϵ2+(1−r)ϵ1,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ~ϵxx=ϵ1ϵ2rϵ1+(1−r)ϵ2, (32) ~αzz=3(rχ(3)2+(1−r)χ(3)1),\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak ~αxx=3rχ(3)2ϵ21+(1−r)χ(3)1ϵ22((1−r)ϵ2+rϵ1)2,

where, is the ratio of the second material in the core, and , , are the linear permittivities and the third order nonlinear susceptibilities of the constituent materials in the core, respectively. Generally, and are complex, the real parts are used in the derivation of the models whereas the imaginary parts will be used to estimate the losses. In this study we will consider only the linear losses, nonlinear processes like two-photon absorption are neglected