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 ₂|E_x |². Therefore, the relative permittivity of the slab is [13]

where n_o is the linear refractive index of medium which is dimensionless, and n ₂ is the nonlinear refractive index with a unit of m ²/V ². 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

 

Fig. 1. A plane wave incident on a slab with Kerr-type nonlinearity.

Download Full Size | PDF

Hence, the nonlinear wave equation in the slab is

in which $k_o=\omega \sqrt{{\mu }_o{\epsilon }_o}$ and α=2n_on ₂.

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 |²) 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_on_o . It can be easily shown that

Now, from Eq. (9), it is obvious that if |E_o |² 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 ₁. 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 ₁(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 ₁ 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 ₁/dζ = 0 and dR ₂/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 ₁ = R ₁(0)exp (jθ ₁(0)) and b ₁ = R ₂(0)exp (jθ ₂(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 ² A/dz ²|<<|dA/dz| and |βd ² B/dz ²|<<|dB/dz|. Obviously, the substitution of the mentioned approximate solution in Maxwell’s equations produces a number of terms proportional to exp (-3jβe) and exp (3jβ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 ₁ and b ₁ are to be determined numerically.

Now, we concentrate our attention on the calculation of a ₁ and b ₁ 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 ₁ = k_o (|a ₁|² + 2|b ₁|²)/(2n_o ), and k ₂ = k_o (2|a ₁|² + |b ₁|²)/(2n_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 ₁, and b ₁. By eliminating b_o from these equations and using the normalized complex amplitudes ā ₁ = a ₁/a_o and b̄ ₁ = b ₁/a_o , one obtains

This system of nonlinear equations is to be numerically solved for determination of ā ₁ and b̄ ₁ 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, ā ₁ and b̄ ₁ 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 ā ₁ and b̄ ₁ . 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.

 

Fig. 2. (a) Variation of magnitude of reflection coefficient at z=0 as a function of a_o . (b) Variation of magnitude of transmission coefficient at z=L as a function of a_o . (c) Normalized reflected (blue line) and transmitted intensities (red line) as a function of a_o . The green line shows the summation of normalized reflected and transmitted intensities.

Download Full Size | PDF

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 ₂= 1.25×10^-18 m ²/V ²) 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⁸ V/m and increases for larger values of the amplitude of the incident wave. The variation of the magnitude of transmission coefficient, i.e. |ā ₂|, 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⁸ 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.

 

Fig. 3. Magnitude of reflection coefficient at z=0 versus wavelength in linear regime.

Download Full Size | PDF

 

Fig. 4. (a) Normalized magnitude of electric field phasor, |E_x /a_o |, for three different values of a_o . (b) Normalized electric field at t=0 (real part of E_x /a_o ) for three different values of a_o . Blue line: a_o =7×10⁶ V/m. Green line: a_o =12×10⁷ V/m . Red line: ao=22×10⁷ V/m.

Download Full Size | PDF

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⁸ 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 ₂ > 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⁸ 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.