BPM is a powerful tool for creating guided-wave optical devices


Rob Scarmozzino

Most commercial optoelectronic CAD software uses the beam propagation method, which provides a practical approach to modeling and simulation of guided-wave devices.

The demands placed on lightwave systems have dramatically increased, spurring development of lightwave components such as fibers, lasers, detectors, modulators, switches, and wavelength multiplexing/demultiplexing devices. All of these components benefit from the development of computer-aided modeling and design software.

Accurate modeling of optoelectronic devices is essential for companies to develop high performance optical components.1 Design tools and device modeling allow engineers to optimize designs, shorten the design cycle for new designs, and evaluate new device concepts.

Among the techniques developed for modeling guided-wave propagation in integrated optics and fiberoptics, the most popular is the beam propagation method (BPM). Unlike the finite-difference time-domain (FDTD) technique, which simulates the device's performance through a rigorous numerical solution of Maxwell's equations, the BPM takes advantage of some simplifications common to many problems (see "FDTD offers high-end solutions," p. 70). If applied appropriately, these simplifications reduce the computational complexity, time to solution, and the hardware requirements with no loss of accuracy.

The BPM is an approach for approximating the wave equation for monochromatic waves and solving the resulting equations numerically. The simplest form of BPM assumes the electric field can be represented as a scalar value rather than vector value (and therefore neglects polarization effects), and assumes that propagation is restricted to a narrow range of angles (that it is paraxial). The scalar field assumption allows the wave equation to be written in the form of the well-known Helmholtz equation for monochromatic waves:

Here the scalar electric field has been written as E(x,y,z,t) = φ(x,y,z)e-iωt, and the notation k(x,y,z) = k0n(x,y,z) has been introduced for the spatially dependent wavenumber, with k0 = 2π/λ being the wavenumber in free space. The geometry of the problem is defined entirely by the refractive index distribution n(x,y,z).

In a typical guided-wave problem the most rapid variation in the field φ is the phase variation due to propagation along the guiding axis, which is predominantly along the z direction. Therefore, this rapid variation can be factored out of the problem by introducing a so-called slowly varying field u, where

Here k is the reference wavenumber, a constant representing the average phase variation of the field φ. When Equation 2 is inserted into Equation 1 the resultant equation contains both first- and second-order partial derivatives of u with respect to z. But since the variation of the field, u, along the propagation axis, z, occurs slowly, the second order derivative can be ignored as negligible, resulting in:

This logic is the heart of BPM.

In practice, this logic can sometimes be simplified to two dimensions by omitting any dependence on y, but for the purposes of this article the 3-D example will be discussed (see Fig. 1).

For a numerical solution one must define:

  • The distribution of the index of refraction throughout the computational domain
  • The size of the computational domain and an appropriate x,y,z computational grid
  • The input field at the z = 0 plane (the amplitude and phase at each of the grid points in the z = 0 plane)
  • A method of handling any energy that impinges upon the outer transverse boundaries of the domain (boundary conditions)

The problem then becomes that of a relatively straightforward numerical integration to determine the field values at the next z plane (see Fig. 2, p. 73). Repeating this, the distribution in each z plane is used as the input field to determine the distribution in the next, until the energy has propagated through the entire structure.

A decision must be made as to the means of handling energy that impinges on the outer transverse boundaries of the computational domain. For example, one can't choose to have the field vanish at the boundary. This decision causes the calculations to assume that the boundary is a perfect reflector. Although some researchers have tried to employ an absorbing material at the edge, setting the correct parameters for such an absorber is complicated, and unintentional reflections may still result.

A commonly used boundary condition is the so-called transparent boundary condition or TBC.2 The basic approach is to assume that near the boundary the field behaves as an outgoing plane wave, with characteristics (amplitude, direction) that are dynamically determined via some heuristic algorithm. In many, but not all cases, the TBC is effective in allowing radiation to freely escape the computational domain.

The simplifications offered by the BPM reduce the computational complexity of most problems, speeding up the computational process for many, and in some cases making otherwise impossible computations realistically solvable. By eliminating the second derivative term in z, the problem is reduced from a second-order boundary value problem requiring iteration or eigenvalue analysis to a first-order initial-value problem that can be solved by simple "integration" along the propagation direction z.

For some problems, adequate results can be obtained even when the distance between calculated z planes (in other words, the longitudinal step size or Dz) is much coarser than the wavelength. However, the choice of step size is of critical importance. On one hand, a smaller step improves accuracy, while a larger step yields a solution more rapidly.

To get this increase in speed, however, the basic method can only model certain types of systems:

  • The fields must be paraxial—in other words, the light must propagate within a limited range of angles about the z axis (approximately ±5°).
  • The rate of change of the index of refraction with respect to z must remain limited if we wish to continue neglecting that second derivative.
  • Because we assumed an average phase variation, devices for which phase variation is critical may not be well served by the basic BPM but can often be handled by a wide-angle BPM approach. Multimode interference devices (MMIs) are one such example.
  • Finally, the method is not suited to calculating backward traveling waves—reflective devices will not be accurately modeled.

The foregoing discussion also assumes isotropic, linear materials, and a scalar field, which eliminates polarization effects. However, none of these are inherent to the BPM and can be included by a rederivation of the basic formula.

The beam propagation method has been extended in a number of ways to overcome specific limitations:

Vectorial BPM. Polarization effects can be included in the BPM by recognizing that the electric field E is a vector, and starting the derivation from the vector wave equation rather than the scalar Helmholtz equation.3

Wide-angle BPM. The paraxiality restriction on the BPM, as well as the related restrictions on index contrast and multimode propagation can be relaxed through the use of extensions that have been referred to as wide-angle BPM.4 The essential idea behind the various approaches is to reduce the paraxial limitations by incorporating, to differing degrees of approximation, the effects of the second order derivative that was neglected in the derivation of the basic BPM (see "BPM simulates performance of AWG," this page).

Bidirectional BPM. While wide-angle BPM allows propagation in a wider cone of angles about the z axis, this cone can only asymptotically approach ±90° from the z axis, and can never be extended to handle simultaneous propagation along the negative z axis (such as 180°). A recent technique that considers multiple interfaces and reflections in a self-consistent and efficient way works quite well (see Fig. 3).5

Additional Extensions. Several additional BPM techniques can handle anisotropic and nonlinear materials. Anisotropic materials can be dealt with by extending the definition of vectorial BPM to include a dielectric tensor.6 Nonlinear materials can be accommodated by allowing the refractive index appearing in the equations to be a function of the optical-field intensity.

The advancement of modeling and simulation software is an integral part of the growth of the optoelectronics industry as it matures toward the size and significance of the current electronics industry. Software based on BPM is a powerful tool for developing guided-wave optical devices. As with most sophisticated design tools, the user benefits from understanding its principles of operation and its limitations.

Rob Scarmozzino is cofounder and chief technical officer at RSoft, 200 Executive Blvd., Ossining, NY 10562. He can be reached at rob@rsoftinc.com.


  1. R. Scarmozzino et. al, IEEE J. Select. Top. Quant. Elect., 6 (1) (January/ February 2000).
  2. G. R. Hadley, Opt. Lett. 16, 624 (1991); G. R. Hadley, J. Quant. Electron. 28, 363 (1992).
  3. W. P. Huang and C. L. Xu, J. Quant. Electron. 29, 2639 (1993).
  4. G. R. Hadley, Opt. Lett.17, 1426 (1992); G. R. Hadley, Opt. Lett. 17, 1743 (1992).
  5. H. Rao, R. Scarmozzino, and R. M. Osgood Jr., Photon. Tech. Lett. 11, 830 (1999).
  6. C. L. Xu, W. P. Huang, J. Chrostowski, and S. K. Chaudhuri, J. Lightwave Tech. 12, 1926 (1994).


FDTD offers high-end solutions
The finite-difference time-domain (FDTD) technique is a widely used and well-established numerical propagation method that has been employed to solve microwave electromagnetic problems. More recently it is being used for optical problems because of, in part, the increased speed of PCs.

On one hand, it is a rigorous numerical solution of Maxwell's equations in the time domain that can, in principle, yield more accurate results than the beam propagation method (BPM). On the other hand, the FDTD technique requires large storage and long computation times, which make it impractical for large-volume 3-D or large-area 2-D simulations. For example, in a 3-D FDTD simulation, a data set of electric and magnetic field vectors must be maintained for every point in a 3-D grid to yield time-domain propagation similar to that along the z direction in BPM. Thus, to propagate a 1 µm wavelength through a 1 mm3 volume with a grid spacing of 0.1 µm, the data set consists of 1012 grid points—a calculation that requires over 50 terabytes of memory.

In practice, for a wide-range of guided-wave design problems, BPM is still the method of choice because it is much more efficient and equally as accurate. At the other end of the spectrum is a class of problems that must use FDTD, such as ring resonators and photonic crystal structures. However, between the two extremes, when pushing the constraints on BPM (for example, a wide-angle problem, but still narrow enough in angle that BPM can handle it), FDTD can be a very effective supplementary tool to BPM; for example, to check the accuracy of the BPM results or to perform hybrid simulations using both BPM and FDTD in problems that have a mix of the two requirements.1


  1. R. Scarmozzino, E. Heller and Z. Huang, SPIE Photonics West 2001, San Jose, CA (January 2001).


BPM simulates performance of AWG

An arrayed waveguide grating (AWG) router1 is an integrated optical wavelength mutiplexer that is an example of a large-scale photonic integrated circuit. An N x N AWG router is composed of N input/output waveguides and two identical star couplers connected by an array of waveguides.

Simulation of such a device in its entirety, while possible, is complicated because of the large domain required by the array and because it entails propagation at large angles from the optical axis. The finite-difference time-domain (FDTD) technique would eliminate the accuracy issues, but this problem is just too large to be practical. A superior method breaks the problem up into sections that can be simulated separately and subsequently combined via a matrix-like approach. Input/output couplers and waveguide sections are simulated via standard beam propagation method (see Fig. A) and the waveguide array is represented simply by a laterally varying phase, derived semianalytically from the arrayed waveguide geometry. Using this approach the entire device is simulated at each wavelength until the full spectrum is produced (see Fig. B), where the individual wavelength responses of each of the eight output ports are overlaid. From this kind of information, the characteristics of this 1 x 8 AWG router can be derived, including channel spacing, passband width and flatness, channel center wavelength, and interchannel crosstalk. This matrix approach is a general method for designing large scale photonic integrated circuits using the beam propagation method.


  1. M. K. Smit and C. van Dam, IEEE J. Select. Top. Quant. Elect. 2(2), 236 (June 1996).
More in Home