## Abstract

Nonlinear evolution of coupled forward and backward fields in a multi-layered film is numerically investigated. We examine the role of longitudinal and transverse modulation instabilities in media of finite length with a homogeneous nonlinear susceptibility *χ ^{(3)}*. The numerical solution of the nonlinear equations by a beam-propagation method that handles backward waves is described.

© 2002 Optical Society of America

## 1. Introduction

Photonic band gap (PBG) structures have opened new pathways for harnessing light. The PBG optoelectronic devices are designed to mold the properties of light by manipulating the dispersion, symmetry and amplitude of the field. Nonlinear optics in PBG’s has emerged as a way to make compact and efficient photonic devices. Many applications have been envisioned in one-dimensional systems consisting of multi-layer, dielectric stacks or patterned waveguides. Gap solitons was a term coined by Chen and Mills [1] to describe the hyperbolic-secant shaped stationary structures they observed in numerical simulations of periodic, nonlinear dielectrics. Their results were thoroughly analyzed by a number of researchers using envelope function formalisms; a summary can be found in the review of deSterke and Sipe [2]. Potential device applications include nonlinear optical limiters [3]; optical diodes [4]; a photonic band edge laser [5]; and a high-gain optical parametric amplifier for nonlinear frequency conversion [6].

In this paper we numerically solve the nonlinear, coupled equations for counter-propagating waves in a finite, periodic medium. An initial pulse is launched outside the medium and the ensuing interaction leads to reflections and pulse break up. The transverse and longitudinal modulation instability occurs at a threshold value of the field amplitude. Modulation instability phenomena occur in a variety of nonlinear optical systems. For instance, they are identified in a homogeneous optical media [7–9] and optical fibers [10–12]. Recently Lichinitzer et al. [12] examined the regions of spatio-temporal instabilities in a Kerr medium with a Bragg grating. Their linear stability analysis of an infinite medium predicts the appearance of transverse and longitudinal modulations of the field amplitude. We present a computational method, which is similar in spirit to the method introduced by Scalora and Crenshaw [13] to describe beam propagation with simultaneous forward- and backward-propagating waves.

## 2. Mode-coupled equations and computational approach

#### 2.1 Mode-coupled Equations

A weakly periodic index profile is assumed in order to simplify the calculations for the electromagnetic modes of the structure. In this paper the nonlinear coefficient is assumed to be homogeneous throughout the medium. The field evolution is described by a set of coupled equations for the forward- and backward propagating waves [10]:

where *δ*=(*k*- *k _{0}*) is the detuning of the laser wavenumber in the medium,

*k*, from half the wavenumber corresponding to the lattice periodicity, i.e.

*k*=π/ʌ, where ʌ is the lattice constant ;

_{0}*κ*=

*ω*/

^{2}Δε*2 k c*, is the grating coupling coefficient;

^{2}*η*=

*12πω*)/

^{2}χ(3)*2kc*, is the nonlinear coefficient (related to the average Kerr coefficient of the medium),

^{2}*F*is a diffraction parameter, and

*v*is the group velocity . Our boundary conditions are applied at opposite ends of the sample:

The function *S*(*x*,*y*,*t*) is the pulse input to the sample. The linear equations are analytically solved by taking Fourier transforms of the time and transverse coordinates, using *exp(-iωt)* for the time-harmonic component. The solutions are

where Ω = *δ* + *ω*/*v*-*q*
^{2}/*F*, and Δ^{2} =Ω^{2} -*κ*
^{2}. The solutions have the transversal and temporal coordinate dependence, albeit in Fourier transformed variables.

#### 2.2 Numerical Computations for Coupled Mode Equations

Equations (1) with boundary conditions (2) are solved by a variation of the Slowly Varying Envelope in Time (SVEAT) method used for deep gratings [13]. Eqs. (1a) and (1b) are integrated by stepping forward in time; a split step operator form is used. The propagation is separated into a 2×2 diagonal matrix of a linear operator, whose components are

The second operator, *V*, is a matrix with diagonal components we denote by *N*

and off-diagonal components, which are denoted by *K*

Eqs. (1) are rewritten in the following vector form (defining *U*= (*E _{f}*,

*E*)

_{b}^{T}, where the superscript

*T*denotes the transpose of the vector) is

The solution for the time increment t+Δt is expressed in the split-step operator form as:

The operator *L* is diagonalized by Fourier transforming the (*x*,*z*) variables. The split-step method therefore uses fast Fourier transforms to simplify the linear operator matrices. The central operator in Eq. (7) can be subsequently decomposed into a linear portion and a nonlinear contribution, namely

The decomposition of the equations may be done in several ways. Here Eqs. (7–8) provide a simply implemented algorithm. The operator L is solved by spectral methods; fast Fourier transforms (FFT) diagonalize the spatial derivatives. In the longitudinal direction we use 2048 points. We use only one transverse direction; in this paper the number of FFT points is 128.

Since this method steps the solution forward in time, the first boundary condition in Eq. (2) is modified to an initial condition by displacing the function from the boundary at *z*=*0*. The relation between the two conditions is *F*
*(x,z,0)*=*S(x,0,*
*(z*-*z _{0}F)v)*, because the pulse propagates without dispersion in the medium outside the nonlinear medium. We also assume that the boundary is impedance matched, so that in the absence of the coupling

*κ*no reflected wave is generated. The initial pulse has been displaced by

*z*

_{0}to a position outside the medium. The accuracy of the code was verified by comparing results with the analytic solutions of the linear equations and by verifying the conservation of energy at each time step.

## 3. Results

Neglecting transverse coupling, the linear dispersion relation is a solution of Eqs. (3) for the form *E _{x}* =

*A*exp(

_{x}*iKz*), where the subscript is

*x*=

*f*,

*b*. The result is a quadratic equation, which can be expressed it in the form

The detuning parameter takes the role of the fields’ angular frequency. The wavenumber is real as long as the magnitude of *δ* is larger than *κ*. A band gap separates the two branches. Tuning the laser close to one band or another determines the characteristics of the wave’s propagation. The transmission curve is determined from the coefficient of F(*q*,*ω*) on the right hand side of Eq.(3a) when *z*=*L*. In our numerical equations the space, time and field amplitude are scaled so that the following parameter values are used: *κ*= *1*,*v* = *1* and *η* = *1*. Only real, positive values of the nonlinear coefficient are discussed here and for simplicity we chose *F*=*100*. For a medium length *L*=6.*1*, the transmission curve is shown in Figure 1.

The equations of motion are solved by stepping forward in time. The Gaussian-shaped pulse is initially displaced to a position outside the medium:

The initial displacement, *z _{0}*, of the pulse is chosen to be several times longer than its pulse width

*σ*= 1/8

_{z}*π*. This prevents any significant overlap between the medium and the leading edge of the wave. The transverse width was chosen as

*σ*=

_{x}*π*/

*64*. The z-axis spans 20π and the x-axis span is 8π. In our computations

*F*=

*100*so that the transverse coupling is weak and does not noticeably affect the pulse tuning at a selected frequency.

Figure 2 is an animation of a low-power pulse amplitude propagating through the sample (*A*= *0.2*). The absolute amplitude of the field at each pixel is represented by a color. The vertical axis displays the transverse coordinate and the horizontal axis is the propagation direction. The evolution is stepped through a time sequence beginning with the pulse to the left of the nonlinear PBG medium. The time increment between frames is *Δt*=*1*. Three subfigures displayed are displayed in the animation. The top figure is the input Gaussian field. The initial pulse is displaced in the window and propagates to the right. The nonlinear medium begins at the far left of the top panel and is not shown. The middle figure is the continuation of the forward-propagating pulse as it passes into and through the medium. The boundaries of the medium are denoted by the green boxes positioned to the left side on the top and bottom of the figure. The bottom panel shows the propagation of the backward field in the medium. The green boxes on the right side of the figure denote the boundaries of the medium. The pulse width is shorter than the medium; it is reflected inside the medium as contributes to pulse break up along the longitudinal axis. This behavior can be expected by comparing the pulse spectrum with the transmission curve in Figure 1. The pulse is short enough to sample a large portion of the dispersive spectrum in the medium.

A similar animation is shown in Figure 3; however, now the input field amplitude has been increased to *A* = *0.6*. In this case we do not plot the input field displayed in the top panel of Figure 2. The laser is still detuned from the center of the gap to the first transmission maximum on the high frequency side, *δ* = *1.12*. As the field interacts with the medium the pulse begins to break up in both the transverse and longitudinal directions. There is also strong focussing of the pulse as it enters the medium.

The front end of the forward-propagating pulse is strongly altered by the passage through the medium, while the backward pulse has a relatively smooth leading edge. The tailing portion of both pulses has strong oscillations in both longitudinal and transverse portions. The pulse spread is increased over the small amplitude case. For higher intensities the pulse break up becomes stronger; the spatial frequencies observed in the transverse direction become higher and the pulse spread is stronger. Increasing diffraction effects by lowering *F* also leads to a larger spread of the pulse in passage through the medium and greater divergence of the nonlinearly modulated beam in free space.

The break up of the pulses is illustrated in the three-dimensional plots in Fig. 4. The top frame is the forward-propagating pulse amplitude is shown for the last frame in Fig. 3. The extent of the pulse’s transverse spreading is more obvious in this figure as the smooth Gaussian envelope break up into multiple peaks. The longitudinal instability is obvious in the lower frame showing the backward-propagating field amplitude. The backward-propagating pulse in Fig. 3 displays a rib structure that is a characteristic of the longitudinal instability.

The modulation instability is readily apparent in our simulations for *A*=*0.4*. But a closer look at the transverse spatial spectra of the forward propagating field, defined by

shows deviations from the initial transverse Gaussian. The spectrum *H*(*q*) is plotted for three amplitudes in Figure 5. The pulses have propagated through the medium. The spectrum for *A*=*0.2* in Figure 5 is Gaussian. However, the spectrum for *A*=*0.3* is already deformed from a Gaussian shape; the wings of the initial Gaussian profile acquire a broad shoulder. At threshold the transverse coordinate has relatively low spatial frequencies. Above threshold the spatial frequencies become higher; we find the maximum develops in the nonlinear regime around *q*=*2.2* and moves towards larger values.

The pulse break up depends on the laser tuning. The normal dispersion regime occurs for tuning frequencies below the band gap, so its ability to shape pulses is changed [12]. The animation in Figure 6 shows pulse evolution for the case where the laser detuning from the gap center is *δ* = -*1.12*. A large field amplitude is chosen for this example *A*= *0.7*. The pulse does undergo a pulse focussing effect in the medium due to the dominance of the nonlinearity. However, the pulse break up does not occur. For our chosen parameters the absence of modulations instability is expected based on the linear stability analysis of Lichtinitzer et al. [12]. The gain curves are peaked at zero wavenumbers in both longitudinal and transverse directions. The one distinctive difference between the upper and lower branches is the opposite signs of their dispersion. The lower branch dispersion is normal, while the upper branch is anomalous. In the normal dispersion regime simulations run to amplitudes of unity have not resulted in pulse break up.

## 4. Conclusions

We have presented a computational algorithm to numerically solve simultaneous forward- and backward-propagating, time-dependent fields in a nonlinear media. In a periodic medium the modulation instability is manifested in the transverse and longitudinal directions. The modulation instabilities lead to a rapid spread of the pulse in the transverse direction and temporal break up of the pulse at threshold intensity. Our threshold modulation instability amplitude is around *A*=*0.3*. This compares favorably with the value *A*=*0.14* we deduce from the paper of Lichinitzer et al. [12] by using the 1/e widths of the initial Gaussian profiles in Fourier transformed variables, i.e. transverse wavenumber, *k _{t}*=

*16*/

*π*

^{2}, and longitudinal wavenumber,

*k*=

_{p}*0.4*. This comparison is only qualitative though, since the stability analysis is made for an infinite medium. In a finite medium we chose the central frequency at the first transmission maximum, because the field is locally enhanced. Our pulses are short and the spectral width of the input pulse samples a spread of dispersion values. Also the modulation instability gain needs a sufficient path length to grow a local perturbation and our sample was kept short. A stability analysis for the finite system is required for a quantitative comparison with our simulations, but the problem may not be manageable by analytic techniques.

The instability in the transverse direction can be applied to designing optical limiters. On the high frequency branch the low intensity transmission is about 0.4 when tuned to the first transmission maximum; it decreases to 0.19 for *A*=*0.7*. On the high frequency branch and tuned to the first transmission maximum the transmission coefficient increases with the input intensity from about 0.4 at low intensity to about 0.6 for *A*=*0.6*; however, the diffraction reduces the on-axis energy, so that by using an apertured detector the off-axis portion is eliminated from the system. To put values on the threshold intensity for the modulation instability consider a sample made from two dielectrics with dielectric contrast, 0.01 and a nonlinear optical coefficient n_{2}= 10^{-10} cm^{2}/W, the corresponding threshold intensity in finite systems is about 3 MW/cm^{2}. At threshold the nonlinear change of the index is about 0.0003.

## Acknowledgements

The work of JWH was supported by NSF grant ECS-9630068.

## References and links

**1. **W. Chen and D. L. Mills, “Gap solitons and the nonlinear optical response of superlattices,” Phys. Rev. Lett. **58**, 160–163 (1987). [CrossRef] [PubMed]

**2. **C. M. de Sterke and J. E. Sipe, “Gap Solitons,” Progress in Optics **33**, 203–260 (1994).

**3. **M. Scalora, J.P. Dowling, C.M. Bowden, and M. J. Bloemer,“ Optical limiting and switching of ultrashort pulses in nonlinear photonic band gap materials,” Phys. Rev. Lett. **73**, 1368–1371 (1994). [CrossRef] [PubMed]

**4. **M. Scalora, J. P. Dowling, M. J. Bloemer, and C. M. Bowden, “The photonic band edge optical diode,” J. Appl. Phys. **76**, 2023–2026 (1994). [CrossRef]

**5. **J. P. Dowling, M. Scalora, M. J. Bloemer, and C. M. Bowden, “The photonic band edge laser: A new approach to gain enhancement,” J. Appl. Phys. **75**, 1896–1899 (1994). [CrossRef]

**6. **M. Scalora, R. J. Flynn, S. B. Reinhardt, R. L. Fork, M. D. Tocci, M. J. Bloemer, C. M. Bowden, H. S. Ledbetter, J. M. Bendickson, J. P. Dowling, and R. P. Leavitt, “Ultrashort pulse propagation at the photonic band edge: Large tunable group delay with minimal distortion and loss,” Phys. Rev. E **54**, R1078–1081 (1996). [CrossRef]

**7. **G. G. Luther and C. J. McKinstrie, “ Transverse modulational instability of collinear waves,” J. Opt. Soc. Am. B **7**, 1125–1141 (1990). [CrossRef]

**8. **L.W. Liou, X.D. Cao, C.J. McKinstrie, and G.P. Agrawal, “Spatiotemporal instabilities in dispersive nonlinear media,” Phys. Rev. A **46**, 4202–4208 (1992). [CrossRef] [PubMed]

**9. **Y. Silberberg, “Collapse of optical pulses,” Opt. Lett. **15**, 1282–1284 (1990). [CrossRef] [PubMed]

**10. **B. J. Eggleton, C. M. deSterke, R. E. Slusher, and J. E. Sipe, “Distributed feedback pulse generator based on nonlinear fiber grating,” Electron. Lett. **32**, 2341–2342 (1996). [CrossRef]

**11. **B. J. Eggleton, C. M. deSterke, and R. E. Slusher, “Nonlinear pulse propagation in Bragg gratings,” J. Opt. Soc. Am. B **14**, 2980–2993 (1997). [CrossRef]

**12. **N. M. Lichinitzer, C. J. McKinstrie, C. M. de Sterke, and G. P. Agrawal, “Spatiotemporal instabilities in nonlinear bulk media with Bragg gratings,” J. Opt. Soc. Am. B **18**, 45–54 (2001). [CrossRef]

**13. **M. Scalora and M. Crenshaw, “Beam propagation method that handles reflections,” Optics Commun. **108**, 191–196 (1994). [CrossRef]