1. Introduction

A number of interesting nonlinear optical phenomena occur due to the interaction between matter and electromagnetic radiation [1]. Because of this interaction, the optical properties of matter are noticeably changed in the presence of strong external radiation and these changes, in turn, affect the propagation of radiation inside the matter. Charged particles in matter are influenced by the electric and magnetic fields associated with electromagnetic waves. In nonlinear optics, only dielectric responses to the electric field are usually considered, because magnetic responses are often negligible in many problems interested in that field. In some magnetic materials, however, nonlinear magnetic responses are important in certain frequency ranges and their nonlinear optical properties have been studied [2,3].

In recent years, there has been a great interest in the so-called left-handed media, or equivalently negative refractive index media [4–7]. In these media, a prototypical example of which is an array of wires and split-ring resonators, both the dielectric permittivity and magnetic permeability are negative in a certain range of frequencies and the effective refractive index is also considered to be negative. Zharov, Shadrivov and Kivshar have suggested that in a split-ring type structure made of a nonlinear dielectric, the nonlinearity in the dielectric response can cause strong nonlinearity in the magnetic response, too [8]. This suggestion gives a motivation for studying the propagation of waves in generalized nonlinear media, where the dielectric and magnetic responses are both nonlinear [9–12].

We are interested in the application of nonlinear media in nano- or micro-structured photonic devices, where the length scale of inhomogeneity is comparable to or smaller than the wavelength. In these devices, the interference effects between the forward and backward propagating waves are important and the usual approximate theoretical methods such as the slowly varying envelope approximation, which ignores them, are often inadequate [13,14]. Therefore it is highly desirable to have a theoretical method that solves the nonlinear wave equation accurately without using any approximation for the analysis of photonic devices.

In this paper, we apply a theoretical method called the invariant imbedding method [15–18] to the study of the propagation of electromagnetic waves in nonlinear inhomogeneous media. This method has been used extensively in the study of wave propagation in various kinds of linear stratified media. The main idea of this method is to transform the boundary value problem of original wave equations to an initial value problem of coupled first-order ordinary differential equations in an exact manner. Since the numerical integration of an initial value problem of coupled first-order ordinary differential equations can be done with extremely high accuracy, this transformation allows very efficient numerical calculations of various wave propagation characteristics in arbitrarily inhomogeneous stratified media.

Although most applications of the invariant imbedding method have been limited to the study of wave propagation phenomena in stratified media, where the wave propagation is essentially one-dimensional, this method still has great advantages over other methods and is capable of solving otherwise difficult problems in a numerically precise manner. For instance, using the invariant imbedding method, one can solve a large number of problems where several coupled waves propagate in arbitrarily inhomogeneous media [19]. Furthermore, it is possible to obtain disorder-averaged wave propagation characteristics for any kind of waves propagating in stratified random media using the method [20]. When resonances occur at discrete points in inhomogeneous media, the wave function shows a singularity at resonance points. This singular problem can also be handled using the invariant imbedding method without much difficulty [21,22].

Many years ago, Babkin and Klyatskin developed an invariant imbedding method for solving the wave equation for electromagneticwaves in stratified nonlinear dielectric media, where only the dielectric response was nonlinear [13]. In the present paper, we generalize their method to the case where both the dielectric permittivity and the magnetic permeability are simultaneously nonlinear and derive the corresponding invariant imbedding equations. We demonstrate the validity and utility of our equations by applying them to several uniform and nonuniform cases, including the propagation of waves in a uniform nonlinear slab. We find that in the presence of external radiation, the uniform slab can change into a complex inhomogeneous medium where regions of positive or negative refractive index and metallic regions appear. We also study the wave transmission properties of periodic nonlinear media and the influence of nonlinearity on the mode conversion phenomena in inhomogeneous plasmas. We argue that our theory is very useful in the study of the optical properties of a variety of nonlinear media including nonlinear negative index media fabricated using wires and split-ring resonators.

2. Wave equation

We consider the propagation of a plane electromagnetic wave of vacuum wave number k ₀=ω/c in inhomogeneous nonlinear media. We focus on the stationary self-action of waves and ignore the generation of harmonics. The wave is assumed to be incident from a uniform region onto a stratified medium, where the dielectric permittivity ε and the magnetic permeability µ vary only along the z axis. The inhomogeneous medium lies in 0 ≤ z ≤ L and the wave propagates in the xz plane. The x component of the wave vector, q, is then a constant.

For obliquely incident waves, we need to distinguish two cases of linear polarization. In the s wave case, the electric field vector is perpendicular to the xz plane and the complex amplitude of the electric field, E=E(z), satisfies

In the p wave case, the magnetic field vector is perpendicular to the xz plane. There exists a simple symmetry between the s and p waves. All results for the p wave case can be obtained by replacing ε and E with μ and -H and vice versa in the equations for the s wave case, where H is the complex amplitude of the magnetic field. In the present paper, we will mainly focus on the s wave. We will also restrict our attention to the cases where both ε and μ depend only on the intensity of the wave field, not on the propagation direction.

We assume that the wave is incident from the region where z>L and is transmitted to the region where z<0. In the inhomogeneous region 0≤z≤L, ε and μ are given by

where ε_L, α, μ_L and β are arbitrary complex functions of z. In the incident (transmitted) region, their values are ε ₁ (ε ₂) and μ ₁ (μ ₂) respectively. The functions f and g are also arbitrary functions of the wave intensity. When θ is defined as the angle of incidence, q is equal to $\sqrt{{\epsilon }_1{\mu }_1}{k}_0\sin \theta$ .

3. Invariant imbedding equations

We consider an s-polarized plane wave Ẽ(x, z)=νu(z)exp(iqx)=νexp[ip(L-z)+iqx], where $p=\sqrt{{\epsilon }_1{\mu }_1}{k}_0\cos \theta$ , incident on the medium from the right. |ν|²(≡w) is the electric field intensity of the incident wave. The quantities of main interest are the reflection and transmission coefficients, r=r(L) and t=t(L), defined by the wave functions outside the medium:

where $p\prime =\sqrt{{\epsilon }_2{\mu }_2}{k}_0\cos \theta \prime$ and θ’ is the angle that outgoing waves make with the negative z axis.

Inside the inhomogeneous medium, the normalized electric field is represented by u(z)=E(z)/ν. Let us consider the u field as a function of z, L and the electric field intensity of the incident wave, w. We start from the crucial observation that the boundary value problem of the wave equation, Eq. (1), can be transformed into an integral equation

where I(z,L,w)≡ |u(z,L,w)|² and J(z,L,w) (≡ |H(z,L,w)|²/w) is given by

In the next step, we take partial derivatives of Eq. (4) with respect to L and w respectively and obtain the quasi-linear partial differential equation [13]

where

Using the relationships

we derive the partial differential equations satisfied by r and t:

where we have changed the variable L to l. One can obtain r(L,w) and t(L,w) by integrating these equations from l=0 to l=L.

Finally, we use the method of characteristics and transform the quasi-linear partial differential equations, Eq. (9), to a set of equivalent ordinary differential equations [18]:

where

and the variables r, t and w are considered to be functions of l. The functions ε(l) and μ(l) are expressed in terms of r(l) and w(l) using

where

In the simplest case where f=|E|² and g=|H|², Eq. (12) takes the form

We note that in a nonlinear system, the variables r and t are coupled to w. The invariant imbedding equations, Eq. (10), are supplemented with the initial conditions for r and t, which are obtained using Fresnel’s formulas:

The initial condition for w is w(0)=w ₀, where w ₀ is chosen such that the final solution for w(L) is the same as the physical input intensity. In general, there are several w ₀ values corresponding to a given w(L) value. For given values of ε ₁, ε ₂, μ ₁, μ ₂, k ₀ and θ and for arbitrary functions ε_L(l), α(l), μ_L(l), β(l), f(|E|²) and g(|H|²), we integrate Eq. (10) using Eqs. (11), (12) and (13) and the initial conditions, Eq. (15), from l=0 to l=L and obtain the reflection and transmission coefficients r(L) and t(L).

We can also derive the invariant imbedding equation for the normalized electric field amplitude u(z)=E(z)/ν inside the inhomogeneous medium. Using the method of characteristics, we transform Eq. (6) to

where u is considered to be a function of z and l. For a given z (0<z<L), u(z,L) is obtained by integrating this equation, together with the equations for r(l) and w(l) in Eq. (10), from l=z to l=L using the initial condition u(z,z)=1+r(z).

4. Application

The invariant imbedding equations presented above, Eqs. (10) and (16), can be used in studying wave propagation phenomena in a wide variety of inhomogeneous nonlinear media. As already stated, the numerical integration of an initial value problem of coupled first-order ordinary differential equations can be done very efficiently with extremely high accuracy. In fact, the numerical error in all of the figures presented in this paper is unnoticeably small.

Before applying our method to generalized media with simultaneously nonlinear dielectric permittivity and magnetic permeability, we test its validity by applying it to ordinary nonlinear dielectric media with nonlinearity only in the dielectric response and comparing the results with previous works. We consider the simplest case where s- and p-polarized plane waves are incident obliquely at an incident angle θ on a uniform dielectric slab with the Kerr-type nonlinearity. The magnetic response is linear and the magnetic permeability is always equal to 1 in this case. In Fig. 1, we plot the transmittance T (=|t|²) obtained by solving Eq. (10) with the initial conditions r(0)=0, t(0)=1 and w(0)=w ₀ as a function of the nonlinearity parameter αw/ε_L, where w (=w(L)) is the input intensity of incident light and ε_L is the linear part of the dielectric permittivity of the slab. The parameters used are k ₀ L=2.2π, ε_L=16, θ=45°, μ_L=ε ₁=ε ₂=μ ₁=μ ₂=1 and β=0. We find that the curves shown here agree perfectly with Fig. 3 and Fig. 7 presented in [23]. We have also reproduced the curves precisely the same as those shown in [24]. In addition, we have applied our method to nonlinear multilayer structures considered in [25] and obtained results which agree very well with those shown there.

 

Fig. 1. (a) Transmittance of a uniform slab with the Kerr-type nonlinearity only in ε plotted versus the nonlinearity parameter αw/ε_L, in the (a) s and (b) p wave cases. The parameters used are k ₀ L=2.2π, ε_L=16, μ_L=ε ₁=ε₂=μ₁=μ₂=1, θ=45° and β=0.

Download Full Size | PDF

The main purpose of this work is to study generalized nonlinear media with simultaneously nonlinear dielectric and magnetic responses. We first consider the simplest case where a plane wave is incident perpendicularly on a uniform slab of thickness L. We consider only the Kerrtype nonlinearity, though our theory can be applied equally easily to more general cases. In order to focus on the effect of nonlinearity, we choose very simple parameter values ε_L=μ_L=ε ₁=ε ₂=μ ₁=μ ₂=1 and assume that α and β are real constants. Then ε and μ inside the slab are given by ε=1+α|E|² and μ=1+β|H|².

It turns out that the most interesting effects are observed when the signs of α and β are different from those of ε_L and μ_L respectively. Therefore, we limit our interest to the cases where α<0 and β<0 here. In the case of perpendicular incidence, the transmittance of our nonlinear slab is always equal to 1 if α=β. We make α and β slightly different by choosing β=0.8α.

In Fig. 2(a), we plot the transmittance T of a slab of thickness k ₀ L=20 versus the nonlinearity parameter |α|w. We note that the transmittance is almost equal to 1 except in the parameter region 0.42<|α|w<1.13, where strong multistability is observed. In Fig. 2(b), we show the dependence of w on the initial value w ₀. Multistability occurs when there are several w ₀ values for a given value of w. For instance, the points A, B, C and D correspond to the same value of αw (≈-1.12579). The distance between A and B is so small that they cannot be distinguished in (a) and (b) (see the inset of (b)). There is no multistability when |α|w<0.42.

 

Fig. 2. (a) Transmittance of a nonlinear slab plotted versus the nonlinearity parameter |α|w. The parameters used are k ₀ L=20, ε_L=μ_L=ε ₁=ε ₂=μ ₁=μ ₂=1, θ=0 and β=0.8α (< 0). (b) Dependence of w on the initial value w ₀. The points A, B, C and D correspond to the same value of |α|w (≈ 1.12579). The distance between B and C is too small for them to be distinguished in (a) and (b) (see the inset of (b)).

Download Full Size | PDF

 

Fig. 3. Spatial distributions of (a) the normalized electric field intensity, (b) the effective dielectric permittivity and magnetic permeability, and (c) the effective refractive index, corresponding to the point A in Fig. 2 (|α|w ₀=0.42191892). The dashed line in (c) is a plot of $\sqrt{-\epsilon \mu }$ when εμ<0.

Download Full Size | PDF

Even in uniform media, the electric and magnetic fields will be nonuniform in the presence of external radiation. Therefore, according to the formulas ε=1+α|E|² and μ=1+β|H|², both ε and μ become nonuniform inside strongly nonlinear media. We have calculated the spatial distributions of the electric field intensity inside the nonlinear slab normalized to the input intensity, |E|²/w (=|u|²), for the parameter values corresponding to the points A, B, C and D in Fig. 2 using Eq. (16). The effective dielectric permittivity is obtained easily from ε=1+αw|u|². Equivalently, we can calculate the spatial distribution of ε using Eq. (14), which becomes ε(l)=1+αw(l)|1+r(l)|² in the present case. r(l) and w(l) are obtained by solving Eq. (10). The effective magnetic permeability is calculated most easily using Eq. (14), which is equivalent to μ(l)=1+βw(l)|1-r(l)|² in the present situation.

Based on the signs of ε and μ, we can distinguish four different cases. If both ε and μ are positive, the effective refractive index n defined as $\sqrt{\epsilon \mu }$ is positive. If both ε and μ are negative, n is defined as $-\sqrt{\epsilon \mu }$ and takes a negative value. In the cases where only ε or μ is negative, the wave becomes evanescent and the medium behaves like a metal. The refractive index becomes a pure imaginary number, the size of which is equal to $\sqrt{-\epsilon \mu }$ . One very interesting consequence of our calculation is that an initially uniform dielectric medium can spontaneously change into a highly inhomogeneous medium where regions of positive or negative refractive index as well as metallic regions appear in the presence of external radiation.

In Fig. 3, we plot the spatial distributions of the normalized electric field intensity, the effective dielectric permittivity and magnetic permeability, and the effective refractive index, corresponding to the point A in Fig. 2 where |α|w ₀=0.42191892. Near this point, w and T change extremely rapidly as w ₀ varies. We find that the electric field distribution is highly nonuniform and the field is strongest near the boundary on the incoming side. This causes both ε and μ to be negative near that boundary and there appears a thin region of negative refractive index in 19.722<k ₀ z<20. In the region 18.142<k ₀ z<19.722, only ε is negative and the medium becomes metallic. The dashed line in (c) is a plot of $\sqrt{-\epsilon \mu }$ in this case. Finally, in the region 0<k ₀ z<18.142, both ε and μ are positive and the medium behaves as a normal dielectric. The complicated spatial dependence of ε and μ is the main reason why the transmittance at the point A (≈ 0.375) is very small.

We have found that the spatial distributions corresponding to the point B (|α|w ₀=0.42195175), which we do not show here, are substantially different from those in Fig. 3, though w ₀ is changed only very slightly. The slab is divided into five layers and there appear a couple of negative index layers near the boundary on the incoming side.

In Fig. 4, we show the result for the point C where |α|w ₀=0.500849. We notice that the nonuniformity is quite strong and the medium is separated into 9 layers. There appear three layers of negative index, two layers of positive index, three layers where only ε is negative and one layer where only μ is negative. The transmittance is a little larger (≈ 0.445) than the values for the points A and B.

In Fig. 5, we show the result for the point D where |α|w ₀=1.0559815. Though the input intensity is the same as the A, B and C cases, the transmittance corresponding to the point D (≈0.938) is quite different. We find that the nonuniformity is weak and the absolute values of ε and μ are very small in this case. We observe that there is no positive index region and the medium is divided into three layers of negative index, two layers where only ε is negative and two layers where only μ is negative.

In addition, we have found that in the weakly nonlinear region where |α|w<0.42, the nonuniformity is rather weak and the refractive index is always positive and smaller than 1. We expect the slab to behave as a normal self-defocusing medium. In the strongly nonlinear region where |α|w>1.13, the nonuniformity is also weak and the refractive index is negative in the entire slab.

 

Fig. 4. Spatial distributions of (a) the normalized electric field intensity, (b) the effective dielectric permittivity and magnetic permeability, and (c) the effective refractive index, corresponding to the point C in Fig. 2 (|α|w ₀=0.500849).

Download Full Size | PDF

 

Fig. 5. Spatial distributions of (a) the normalized electric field intensity, (b) the effective dielectric permittivity and magnetic permeability, and (c) the effective refractive index, corresponding to the point D in Fig. 2 (|α|w ₀=1.0559815).

Download Full Size | PDF

 

Fig. 6. Transmittance of a nonlinear slab plotted versus the nonlinearity parameter |α|w in the s wave case. The parameters used are k ₀ L=20, ε_L=μ_L=ε ₁=ε ₂=μ ₁=μ ₂=1, θ=1° and β=0.8α(<0). The point P corresponds to αw ₀=-0.5982.

Download Full Size | PDF

Next we study the propagation of plane waves incident obliquely on the nonlinear slab considered in Figs. 2, 3, 4 and 5. It turns out that when the slab is sufficiently thick, wave propagation phenomena depend very strongly on the incident angle and the polarization. In Fig. 6, we plot the transmittance of an s wave incident at a very small incident angle (θ=1°) versus the nonlinearity parameter |α|w. We find that this figure is substantially different from Fig. 2(a) in the region 0.42<|α|w<1.13.

 

Fig. 7. Spatial distributions of (a) the normalized electric field intensity, (b) the effective dielectric permittivity, and (c) the effective magnetic permeability, corresponding to the point P in Fig. 6. Note that μ changes discontinuously wherever it passes through zero.

Download Full Size | PDF

In Fig. 7, we show the spatial distributions of the normalized electric field intensity and the effective dielectric permittivity and magnetic permeability, corresponding to the point P in Fig. 6, where αw ₀=-0.5982. We observe that μ changes discontinuously wherever it passes through zero. Similar discontinuities are observed in the spatial dependence of ε in the p wave case. This phenomenon, occurring due to resonant enhancement of the z component of the magnetic (electric) field near the μ=0 (ε=0) points in the s (p) wave case, has been studied originally in the context of nonlinear plasmas in [26] and [27].

 

Fig. 8. (a) Scatter plot of the transmittance of a nonlinear slab with sinusoidal spatial variations of ε_L (=1+0.3sin(2πz/Λ)) and μ_L (=1-0.2sin(2πz/Λ)) versus the nonlinearity parameter |α|w. The parameters used are k ₀Λ=2.85, L/Λ=20, ε ₁=ε ₂=μ ₁=μ ₂=1, θ=0 and β=0.8α (< 0). (b) Scatter plot of the dependence of w on the initial value w ₀. The point Q corresponds to the third resonant transmission peak with αw ₀=-0.03547.

Download Full Size | PDF

 

Fig. 9. Spatial distribution of the normalized electric field intensity corresponding to the point Q in Fig. 8.

Download Full Size | PDF

So far we have considered only the case of uniform slabs. The main advantage of the invariant imbedding method, however, is in its efficiency in solving strongly inhomogeneous cases. We illustrate this by applying our method to two kinds of nonuniform media. The first case we study is a nonlinear slab with sinusoidal spatial variations of the linear parts of ε and μ, which are given by ε_L=1+0.3sin(2πz/Λ) and μ_L=1-0.2sin(2πz/Λ). We note that ε_L and μ_L vary periodically with opposite phases and the same period Λ. The slab thickness is L=20Λ and we consider only the case of perpendicular incidence. In the absence of nonlinearity, there exists a photonic band structure in the frequency spectrum and the lowest forbidden band gap appears in the frequency region 2.8<k ₀Λ< 3.6. The nonlinear coefficients are assumed to satisfy β=0.8α<0.

In Fig. 8(a), we plot the transmittance versus the nonlinearity parameter |α|w. We choose the wave frequency to be slightly above the lower edge frequency of the first forbidden gap, k ₀Λ=2.85. In the linear case with αw=0, there is no transmission. We observe that as the strength of nonlinearity increases, the influence of periodicity becomes weaker and the transmittance increases. We find weak multistability and a series of resonant transmission peaks where T=1. In the present case with α<0 and β<0, as the strength of nonlinearity increases above |α|w=0.37, an intriguing phenomenon occurs. In the parameter region where 0.37<|α|w<2, extremely strong multistability with a huge number of multiple solutions sets in. Since it is difficult to plot curves in this extreme case, we calculate T for α ₀ w=0.00005n (n=1,2,…) and plot the result in a scatter plot. In Fig. 8(b), we plot |α|w versus |α|w ₀ in a similar scatter plot. Extremely strong multistability is clearly seen in the region 0.37<|α|w ₀<1.34.

In the presence of weak nonlinearity, transmission of electromagnetic waves through a photonic crystal, when the frequency is within the forbidden gap, is mediated by the formation of gap solitons [25]. It is well-known that the n-th resonant transmission peak corresponds to the formation of an n-soliton state. In Fig. 9, we show the spatial distribution of the normalized electric field intensity corresponding to the point Q in Fig. 8, where αw ₀=-0.03547. This point corresponds to the third resonant transmission peak and the electric field distribution shows the penetration of three gap solitons as expected. peak and the electric field distribution shows the penetration of three gap solitons as expected.

 

Fig. 10. (a) Absorptance of a nonlinear slab with the Kerr-type nonlinearity only in ε and with a parabolic spatial variation of ε_L (=1+8z(z-L)/L ²+0.00001i), plotted versus the nonlinearity parameter αw (> 0) in the p wave case. The parameters used are k ₀ L=20, μ_L=ε ₁=ε ₂=μ ₁=μ ₂=1, θ=1° and β=0. The point R corresponds to αw ₀=0.0001156.

Download Full Size | PDF

As the second example of strongly inhomogeneous cases, we consider a nonuniform plasma with nonlinearity only in ε and study the influence of nonlinearity on the phenomenon called mode conversion or resonant absorption, which is the resonant conversion of transverse electromagnetic waves into longitudinal plasma oscillations at the points where ε=0 in the case of unmagnetized plasmas [21,22,28]. We assume that the linear part of ε is given by ε_L=1+8z(z-L)/L ²+0.00001i, the real part of which becomes zero at z/L=0.146 and 0.854. We also assume a self-focusing nonlinearity with α>0. Other parameters used are k ₀ L=20, μ_L=ε ₁=ε ₂=μ ₁=μ ₂=1, θ=1° and β=0. We consider only the p wave case, since mode conversion cannot occur in the s wave case in the present situation.

 

Fig. 11. Spatial distributions of (a) the normalized magnetic field intensity and (b) the effective dielectric permittivity, corresponding to the point R in Fig. 10. ε changes discontinuously wherever it passes through zero. The discontinuities at k ₀ z=2.95 and 17.05 are too small to be noticed.

Download Full Size | PDF

In Fig. 10, we plot the absorptance, A=1-|r|²-|t|², versus the nonlinearity parameter αw when θ (=1°) is very small. In the linear case, the absorptance is small (A=0.0017). As the nonlinearity parameter increases above αw=0.00016, optical bistability sets in and the absorptance can be either small or large. The maximum absorption of A=0.255 occurs when αw=0.00016 and αw ₀=0.0001156 and is designated by the point R.

In Fig. 11, we show the spatial distributions of the normalized magnetic field intensity and the effective dielectric permittivity, corresponding to the point R in Fig. 10. We find that the magnetic field becomes extremely strong near k ₀ z=9.02 and 10.98, where ε passes through zero. This phenomenon can be understood from the behavior of ε shown in Fig. 11(b), where a sharp peak is formed at the center of the slab. This plays a role of a narrow potential well for electromagnetic waves and a standing wave is formed there. Due to the interplay between the standing wave formation and the mode conversion at the points where ε=0, the field gets resonantly enhanced and the absorptance becomes very large. Similarly to Fig. 7(c), we note that ε changes discontinuously wherever it passes through zero. The discontinuities at k ₀ z=2.95 and 17.05 are too small to be noticed.

5. Conclusion

In summary, we have developed a new version of the invariant imbedding method, which allows us to solve the electromagnetic wave equations in stratified media with nonlinearity in both dielectric and magnetic responses in a numerically precise manner. We have applied our theory to a uniform nonlinear slab and found that in the presence of external radiation, an initially uniform medium of positive refractive index can change into a highly inhomogeneous medium where regions of positive or negative index as well as metallic regions can appear. We have also studied the wave transmission properties of periodic nonlinear media and the influence of nonlinearity on the mode conversion phenomena in inhomogeneous plasmas.

In the present paper, we have concentrated on the simplest kind of nonlinearity, namely the Kerr-type nonlinearity. In [8], the authors derived a somewhat complicated functional dependence of μ on the magnetic field intensity in nonlinear negative index media fabricated using wires and split-ring resonators. This suggests that one has to be careful in applying the Kerrtype model to nonlinear negative index media. Since our method can handle equally easily the cases where f and g in Eq. (2) are complicated functions of the electric and magnetic field intensities respectively, it can be directly applied to the model proposed in [8]. Research in this direction will be pursued in the future.