## Abstract

Multiple-scale analysis is employed for the analysis of plane-wave refraction at a nonlinear slab. It will be demonstrated that the perturbation method will lead to a nonuniformly valid approximation to the solution of the nonlinear wave equation. To construct a uniformly valid approximation, we will exploit multiple-scale analysis. Using this method, we will derive the zeroth-order approximation to the solution of the nonlinear wave equation analytically. This approximate solution clearly shows the effects of self-phase modulation (SPM) and cross-phase modulation (XPM) on plane-wave refraction at the nonlinear slab. As will be shown, the proposed method can be generalized to the rigorous study of nonlinear wave propagation in one-dimensional photonic band-gap structures.

©2005 Optical Society of America

## 1. Introduction

In recent years, there have been many studies on nonlinear wave propagation in one-dimensional photonic band-gap structures as a result of which a number of phenomena such as bistability, gap and Bragg solitons has been discovered in these structures [1–4]. There is an increasing interest among researchers for theoretical study and experimental observation of these phenomena as well as their applications to optical signal processing [5]. The coupled mode theory or the nonlinear Shrödinger equation (NLSE) have been conventionally used for investigation of nonlinear wave propagation in these structures near the Bragg wavelength and where the grating is very shallow [5]. However, rigorous study of nonlinear wave propagation in these structures opens up a way to interpret the phenomena which cannot be accurately explained using the conventional methods. A rigorous study of wave refraction in a nonlinear slab along with a method applicable to nonlinear photonic band-gap structures will provide a reliable approach for analyzing nonlinear wave propagation in these structures.

In a pioneering paper [6], Armstrong *et al*. investigated the interaction of light waves at different frequencies with a nonlinear dielectric having second-order nonlinearity. They proposed a Quasi-Phase Matching (QPM) technique for efficient second-harmonic generation. Later, Bloembergen and Pershan [7] studied the propagation of harmonic waves emanating from the interface of a linear medium with a nonlinear one, and derived the general laws of reflection and refraction of harmonic waves at the boundary. In addition, Carroll [8] studied the propagation of a circularly polarized *monochromatic* plane electromagnetic wave at the interface of a nonlinear half-space. However, to our best knowledge, there exist only limited works on the problem of plane-wave refraction at dielectric slabs with Kerr-type nonlinearity. For example, in [9] by means of Jacobian elliptic functions, optical tunneling through a nonlinear film has been studied in the special case where the background refractive index of the film is less than that of the adjacent media for angles near total internal reflection. The majority of authors have concentrated on the analysis of nonlinear slab-guided waves where the nonlinear propagation modes of a nonlinear slab sandwiched between two or more layers of homogeneous linear dielectrics have been studied [10–12].

It is the aim of this work to study the refraction of normally incident plane wave in a dielectric slab with Kerr-type nonlinearity. As a first step, plane-wave refraction in the slab is modeled by wave propagation in a nonlinear transmission line. A multiple-scale analysis will then be exploited for derivation of the solution to the resulting nonlinear wave equation.

## 2. Formulation

Fig. 1 depicts a plane wave incident on a nonlinear slab. It is assumed that the slab is made of Kerr-type nonlinear material for which the refractive index is given by *n*= *n*_{o}
+*n*
_{2}|*E*_{x}
|^{2}. Therefore, the relative permittivity of the slab is [13]

where *n*_{o}
is the linear refractive index of medium which is dimensionless, and *n*
_{2} is the nonlinear refractive index with a unit of *m*
^{2}/*V*
^{2}. Since the problem is one-dimensional and the incident plane wave is normal to the slab, only *E*_{x}
and *H*_{y}
field components exist. Maxwell’s equations in the slab are given as follows

Hence, the nonlinear wave equation in the slab is

in which ${k}_{o}=\omega \sqrt{{\mu}_{o}{\epsilon}_{o}}$ and *α*=2*n*_{o}*n*
_{2}.

According to Eq. (2) and Eq. (3), we can model wave propagation in the slab with wave propagation in a nonlinear transmission line with an inductance per unit length of *L* = *μ*_{o}
and a capacitance per unit length of *C* = *ε*_{o}
(${n}_{o}^{2}$ + *α*|*E*_{x}
|^{2}) where the electric field *E*_{x}
and magnetic field *H*_{y}
play the role of voltage and current on this line, respectively. Since the region *z* < 0 is free space and a plane wave is incident on the slab, the equivalent nonlinear transmission line is connected to a linear transmission line with inductance per unit length of *L*_{o}
= *μ*_{o}
and capacitance per unit length of *C*_{o}
= *ε*_{o}
with a known forward wave which resembles the incident wave on the slab. Since the region *z* > *L* is free space, we terminate the equivalent nonlinear transmission line to free-space characteristic impedance *Z*_{o}
= *μ*_{o}
/*ε*_{o}
. Hence, we concentrate on the derivation of the voltage and current on the equivalent nonlinear transmission line with the mentioned boundary conditions.

It should be noted that *α* is typically very small compared with the linear relative permittivity of nonlinear materials. This allows us to use perturbation method for solving the nonlinear wave equation. In this method, the electric field *E*_{x}
is expanded to power series of *α* [14], i.e.

If this expansion is substituted in the nonlinear wave equation and the terms with equal powers of *α* are balanced, the differential equations for determining the coefficients of this expansion are derived. For example, the differential equations for the zeroth- and first-order perturbation approximation are given by

Note that the boundary conditions for Eq. (6) are those of the original problem whereas the boundary conditions for Eq. (7) and other higher-order approximations are zero. In other words, we assume *E*_{m}
= 0 and *dE*_{m}
/*dz* = 0 at *z* = 0 and *z* = *L* for every *m* ≥ 1. From Eq. (6), it is obvious that the zeroth-order approximation *E*_{o}
is given by

in which *a* and *b* are complex numbers representing the amplitude and phase of the forward and backward waves and *β* = *k*_{o}*n*_{o}
. It can be easily shown that

Now, from Eq. (9), it is obvious that if |*E*_{o}
|^{2} is multiplied by *E*_{o}
, some terms proportional to exp (-*jβz*) and exp (*jβz*) will appear in the right-hand side of Eq. (7). These terms are in fact the solutions of the homogeneous differential equation. This means that some nonphysical terms (secular terms [14]) of the form *z* exp (*jβz*) and *z* exp (-*jβz*) will appear in the solution of *E*
_{1}. This implies that if the slab thickness approaches infinity, the solution will be unbounded. In fact, the solution given by perturbation series Eq. (5) does not converge uniformly to the solution of the original nonlinear wave equation, i.e. with increasing *z* the error of the approximate solution increases rapidly. Note that although the individual terms of the perturbation series are secular, the secularity disappears when the series is summed up [14]. For removing this secularity, use will be made of multiple-scale analysis to be explained in the next section.

#### 2.1 Multiple-scale analysis for removing secularity

For eliminating the most secular terms to all orders, a new variable ζ = *αz* is introduced. Even though the exact solution *E*_{x}
(*z*) is a function of *z* alone, multiple-scale analysis seeks solutions which are functions of both ζ and *z* treated as *independent* variables. We wish to emphasize that expressing *E*_{x}
as a function of two variables is merely a mathematical technique to remove secularity; the actual solution has *z* and ζ related by ζ = *αz* so that *z* and ζ are ultimately not independent.

The formal procedure consists of assuming a perturbation expansion of the form

Chain rule for partial differentiation is to be used for computing derivatives of *E*_{x}
(*z*). Hence, for the first-order derivative with respect to *z*, we have

However, since ζ = *αz*,

Similarly, it can be shown that for the second-order derivative with respect to *z*, we have

If this relation is substituted in the nonlinear wave equation and the terms with equal powers of *α* are balanced, we will obtain the following differential equations for *E*_{o}
(*z*,ζ) and *E*
_{1}(*z*,ζ)

It is obvious that *E*_{o}
is given by

*a*(ζ) and *b*(ζ) will be determined under the condition that secular terms do not appear in the solution to Eq. (15). From Eq. (16), the right-hand side of Eq. (15) is

If the coefficients of exp (-*jβz*) and exp (*jβz*) are nonzero, then the solution to *E*
_{1} would be secular. To preclude the appearance of secularity, we require that *a*(ζ) and *b*(ζ) satisfy

For solving these equations, we write *a*(ζ) and *b*(ζ) in their polar representation, i.e.

It can be easily shown that *dR*
_{1}/*d*ζ = 0 and *dR*
_{2}/*d*ζ = 0 and the differential equations for determining the angles are

Hence, *a*(ζ) and *b*(ζ) are given by

Using complex representation of *a*(ζ) and *b*(ζ), we will arrive at the following relations

in which *a*
_{1} = *R*
_{1}(0)exp (*jθ*
_{1}(0)) and *b*
_{1} = *R*
_{2}(0)exp (*jθ*
_{2}(0)) are complex constants which are found by satisfying the boundary conditions. This solution clearly shows the simultaneous effects of SPM and XPM on wave refraction. In another word, the nonlinearity of the medium has imposed a nonlinear phase shift on the forward and backward waves. As can be seen, the forward and backward waves impose nonlinear phase shift on themselves (SPM) and on each other (XPM).

For the sake of simplicity, the higher-order expansion coefficients in Eq. (10) are neglected. From Eq. (10), it can be easily seen that the approximate solution of the nonlinear wave equation given by Eqs. (24) and (25) is accurate to *O*(*α*). This approximation is completely consistent with the standard approximations (see e.g. [1,5,6]) for the derivation of the solution of the nonlinear wave equation. In the majority of the conventional methods, the solution of the nonlinear wave equation is assumed as *E*_{x}
(*z*) = *A*(*z*) exp (-*jβz*)+*B*(*z*) exp (*jβz*) This approximate solution is then substituted in Maxwell’s equations and the differential equations for the determination of *z*-dependent amplitudes of forward and backward waves are obtained in the framework of slowly varying amplitude functions, i.e., under the condition of |{*βd*
^{2}
*A*/*dz*
^{2}|<<|*dA*/*dz*| and |*βd*
^{2}
*B*/*dz*
^{2}|<<|*dB*/*dz*|. Obviously, the substitution of the mentioned approximate solution in Maxwell’s equations produces a number of terms proportional to exp (-3*jβe*) and exp (3*jβz*), but these terms are usually neglected. In these methods, the aforementioned differential equations are usually solved using numerical methods such as shooting method. In contrast to the above mentioned methods, utilizing the multiple-scale analysis, we have derived analytical expressions for these *z*-dependent amplitudes; only the amplitudes *a*
_{1} and *b*
_{1} are to be determined numerically.

Now, we concentrate our attention on the calculation of *a*
_{1} and *b*
_{1} from the boundary conditions. The required boundary condition at *z* = 0 is the continuity of tangential electric and magnetic fields, i.e., the continuity of the equivalent voltage and current at this interface. Therefore, we need to find an expression for the magnetic field *H*_{o}
(*z*,ζ). It can be simply shown that

in which *Z*_{c}
= *Z*_{o}
/*n*_{o}
, *k*
_{1} = *k*_{o}
(|*a*
_{1}|^{2} + 2|*b*
_{1}|^{2})/(2*n*_{o}
), and *k*
_{2} = *k*_{o}
(2|*a*
_{1}|^{2} + |*b*
_{1}|^{2})/(2*n*_{o}
). From Eq. (16) and Eq. (26), the boundary condition at *z* = 0 or the condition of continuity of voltage and current at this interface is as

where *a*_{o}
and *b*_{o}
, as can be seen in Fig. 1, are the complex amplitudes of the forward and backward waves of the equivalent linear transmission line connected to the nonlinear transmission line at *z* = 0. It is assumed that the nonlinear transmission line is terminated to impedance *Z*_{L}
at *z* = *L*. Therefore, the boundary condition at this interface will be

Using Eqs. (16) and (26), it can be shown that the boundary condition at *z*= *L* is given by

in which *Z̄* = *Z*_{L}
/*Z*_{c}
is the load impedance normalized to the characteristic impedance of the
line in the linear regime. Note that *a*_{o}
being the complex amplitude of the incident wave is assumed known. Hence, Eqs. (27), (28), and (30) form a system of nonlinear equations for determination of *b*_{o}
, *a*
_{1}, and *b*
_{1}. By eliminating *b*_{o}
from these equations and using the normalized complex amplitudes *ā*
_{1} = *a*
_{1}/*a*_{o}
and *b̄*
_{1} = *b*
_{1}/*a*_{o}
, one obtains

This system of nonlinear equations is to be numerically solved for determination of *ā*
_{1} and *b̄*
_{1} for a known amplitude of the incident wave *a*_{o}
. These quantities are the only unknowns that characterize the solution. After their evaluation, we would be able to derive the electric and magnetic fields. These equations can be easily solved using the conventional methods for finding the roots of a nonlinear system of equations. The fact that facilitates the derivation of the roots of this system is that there exists a good initial guess for the solution when *a*_{o}
is relatively weak. In this case, *ā*
_{1} and *b̄*
_{1} are nearly equal to their counterparts in the linear regime. Hence, for the derivation of the roots of the system for a specific value of the incident-wave amplitude, *a*_{o}
can initially be assumed to be weak, so that the solution in the linear regime can be regarded as the initial guess for the derivation of *ā*
_{1} and *b̄*
_{1} . Later, this parameter is gradually increased to the desired value so that the solution in the last step can be regarded as the initial guess for the derivation of the solution in the present step until the solution for desired value of *a*_{o}
is found.

## 3. Numerical results

Using the method outlined in the last section, we have studied the refraction of an incident plane wave on a nonlinear slab made of Type-RN Corning glass (*n*_{o}
= 2.46 and *n*
_{2}= 1.25×10^{-18}
*m*
^{2}/*V*
^{2}) at *λ*_{o}
= 1.53μ*m*. As a specific example, we have assumed *L* = *λ*_{o}
and *Z*_{L}
= *Z*_{o}
. The variation of the magnitude of reflection coefficient, i.e. |*b̄*_{o}
|, at *z* = 0 with respect to the magnitude of the incident wave *a*_{o}
has been computed and depicted in Fig. 2 (a). This Fig. shows that as *a*_{o}
is increased, the reflection coefficient decreases and approaches zero for *a*_{o}
= 1.927×10^{8}
*V*/*m* and increases for larger values of the amplitude of the incident wave. The variation of the magnitude of transmission coefficient, i.e. |*ā*
_{2}|, seen at *z* = *L* has been depicted in Fig. 2 (b). It also shows that as the intensity of the incident electric field is increased, the magnitude of transmission coefficient is increased and approaches unity for *a*_{o}
= 1.927×10^{8}
*V*/*m* and decreases for larger amplitudes of the incident wave. In addition, the reflected and transmitted intensities have been depicted in Fig. 2 (c) with blue and red lines, respectively. As can be seen, these quantities sum up to unity which ensures the conservation of power in our proposed solution.

As shown in Fig. 2 (a) and (b), there is a resonance in the magnitude of the reflection and transmission coefficient at *a*_{o}
=1.927×10^{8}
*V*/*m*. The nonlinear phase shift imposed by SPM and XPM is responsible for this behavior of reflection and transmission coefficients. The variation of the magnitude of the reflection coefficient at *z* = 0 for wavelengths in the neighborhood of the operating wavelength *λ*_{o}
= 1.53μ*m* has been shown in Fig. 3 for the linear regime. As can be seen, there is a resonance in the reflection coefficient of the structure for wavelengths shorter than the operating wavelength. Since Type-RN Corning glass is a self-focusing medium (*n*
_{2} > 0), its refractive index increases with increasing intensity. Consequently, as the incident-wave magnitude is increased, the resonance shifts toward the operating wavelength *λ*_{o}
= 1.53μ*m*. For *a*_{o}
=1.927×10^{8}
*V*/*m*, the resonance is shifted to *λ*_{o}
= 1.53μ*m*. If the magnitude of the incident wave is increased further, the resonance shifts away the operating wavelength and transmission coefficient decreases.

The normalized magnitude and real part of the electric field phasor inside the slab for three different values of *a*_{o}
have been shown in Fig. 4 (a) and Fig. 4 (b), respectively. These Figs. clearly illustrate the simultaneous effects of SPM and XPM on the wave refraction in the nonlinear slab. As *a*_{o}
is increased, the distribution of the field in the slab deviates rapidly from the field distribution in the linear regime.

## 4. Conclusions

A perturbation approach based on a multiple-scale analysis has been proposed for the study of plane-wave refraction at a slab with Kerr-type nonlinearity. A number of alternative analytical and numerical techniques could have been used for this analysis. For example, a Jacobian elliptic function can be used as the analytical solution of the nonlinear wave equation [9,15,16]. The latter suffers from a number of disadvantages. Firstly, since Jacobian elliptic functions does not have explicit representations, satisfaction of complicated boundary conditions such as impedance boundary condition, which was considered in the present article, is very difficult [9,16]. Secondly, in this approach, the solution cannot be easily subdivided into forward and backward waves and the reflection and transmission coefficients cannot be easily defined. Thirdly, the generalization of the Jacobian method to the study of nonlinear wave propagation in one-dimensional photonic band-gap structures is very involved. A number of numerical methods such as FDTD or BPM can also be regarded as candidates for the study of the present problem, but these methods are computationally inefficient and suffer from high memory requirements. In contrast to the aforementioned methods, our proposed method clearly shows the simultaneous effects of SPM and XPM on the wave refraction and can be readily generalized to the rigorous analysis of nonlinear wave propagation in one-dimensional photonic band-gap structures. As was demonstrated, the proposed approximate solution is completely consistent with the standard approximations for derivation of solution of nonlinear wave equation. In addition, it can be simply shown that there is a perfect agreement between our results for *L*→ ∞ and the ones reported in [8] for the case of plane wave refraction at the interface of a half-space of nonlinear medium with a linear medium. As shown, the only numerical step in our method is the solution of a system of nonlinear equations which can be simply accomplished using existing efficient methods such as, Newton-Raphson method.

## Acknowledgments

The authors would like to thank both the Research Council and the Center of Excellence for Applied Electromagnetics Systems at the University of Tehran.

## References and links

**1. **H. G. Winful, J. H. Marburger, and E. Garmire, “Theory of bistability in nonlinear distributed feedback structures,” Appl. Phys. Lett. **35**, 379–381 (1979). [CrossRef]

**2. **D. N. Christodoulides and R. I. Joseph, “Slow Bragg solitons in nonlinear periodic structures,” Phys. Rev. Lett. **62**, 1746–1749 (1989). [CrossRef] [PubMed]

**3. **N. G. R. Broderick, D. J. Richardson, and M. Ibsen, “Nonlinear switching in a 20-cm long fiber Bragg grating,” Opt. Lett. **25**, 536–538 (2000). [CrossRef]

**4. **K. Senthilnathan, P. Malathi, and K. Porsezian, “Dynamics of nonlinear pulse propagation through a fiber Bragg grating with linear coupling,” J. Opt. Soc. Am. B **20**, 366–372 (2003). [CrossRef]

**5. **R. E. Slusher and B. J. Eggleton, eds., *Nonlinear Photonic Crystals* (Springer-Verlag Berlin Heidelberg, Berlin, 2003).

**6. **J. A. Armstrong, N. Bloembergen, J. Ducuing, and P. S. Pershan, “Interaction between light waves in a nonlinear dielectric,” Phys. Rev. **127**, 1918–1939 (1962). [CrossRef]

**7. **N. Bloembergen and P. S. Pershan, “Light waves at the boundary of nonlinear media,” Phys. Rev. **128**, 606–622 (1962). [CrossRef]

**8. **M. M. Carroll, “Plane waves of constant amplitude in nonlinear dielectrics,” Phys. Rev. A **6**, 1977–1980 (1972). [CrossRef]

**9. **Th. Peschel, P. Dannberg, U. Langbein, and F. Lederer, “Investigation of optical tunneling through nonlinear films,” J. Opt. Soc. Am. B **5**, 29–36 (1988). [CrossRef]

**10. **K. Hayata, M. Nagai, and M. Koshiba, “Finite-element formalism for noninear slab-guided waves,” IEEE Trans. Microwave Theory Tech. **36**, 1207–1215 (1988). [CrossRef]

**11. **S. V. Polstyanko, R. Dyczij-Edlinger, and J. F. Lee, “A full vectorial analysis of a nonlinear slab waveguide based on the nonlinear hybrid vector finite-element method,” Opt. Lett. **21**, 98–100 (1996). [CrossRef] [PubMed]

**12. **V. Van and S. K. Chaudhuri, “A hybrid implicit-explicit FDTD scheme for nonlinear optical waveguide modeling,” IEEE Trans. Microwave Theory Tech. **47**, 540–545 (1999). [CrossRef]

**13. **R. M. Joseph and A. Taflove, “FDTD Maxwell’s equations models for nonlinear electrodynamics and optics,” IEEE Trans. Microwave Theory Tech. **45**, 364–374 (1997).

**14. **C. M. Bender and S. A. Orszag, *Advanced Mathematical Methods for Scientists and Engineers* (McGraw-Hill, Singapore, 1978).

**15. **K. Ogusu, “Self-switching in hollow waveguides with a Kerrlike nonlinear permittivity,” IEEE J. Lightwave Technol. **8**, 1541–1547 (1990). [CrossRef]

**16. **U. Trutschel, F. Lederer, and M. Golz, “Nonlinear guided waves in multilayer systems,” IEEE J. Quantum Electron. **25**, 194–200 (1989). [CrossRef]