## How to Use the Beam Envelopes Method for Wave Optics Simulations

##### Yosuke Mizuyama January 8, 2018

In the wave optics field, it is difficult to simulate large optical systems in a way that rigorously solves Maxwell’s equation. This is because the waves that appear in the system need to be resolved by a sufficiently fine mesh. The beam envelopes method in the COMSOL Multiphysics® software is one option for this purpose. In this blog post, we discuss how to use the *Electromagnetic Waves, Beam Envelopes* interface and handle its restrictions.

### Comparing Methods for Solving Large Wave Optics Models

In electromagnetic simulations, the wavelength always needs be resolved by the mesh in order to find an accurate solution of Maxwell’s equations. This requirement makes it difficult to simulate models that are large compared to the wavelength. There are several methods for stationary wave optics problems that can handle large models. These methods include the so-called diffraction formulas, such as the Fraunhofer, Fresnel-Kirchhoff, and Rayleigh-Sommerfeld diffraction formula and the beam propagation method (BPM), such as paraxial BPM and the angular spectrum method (Ref. 1).

Most of these methods use certain approximations to the Helmholtz equation. These methods can handle large models because they are based on the propagation method that solves for the field in a plane from a known field in another plane. So you don’t have to mesh the entire domain, you just need a 2D mesh for the desired plane.

Compared to these methods, the *Electromagnetic Waves, Beam Envelopes* interface in COMSOL Multiphysics (which we will refer to as the *Beam Envelopes* interface for the rest of the blog post) solves for the exact solution of the Helmholtz equation in a domain. It can handle large models; i.e., the meshing requirement can be significantly relaxed if a certain restriction is satisfied.

*A beam envelopes simulation for a lens with a millimeter-range focal length for a 1-um wavelength beam.*

We discuss the *Beam Envelopes* interface in more detail below.

### Theory Behind the Beam Envelopes Interface

Let’s take a look at the math that the *Beam Envelopes* interface computes “under the hood”. If you add this interface to a model and click the *Physics Interface* node and change *Type of phase specification* to *User defined*, you’ll see the following in the *Equation* section:

Here, \bf E1 is the dependent variable that the interface solves for, called the *envelope function*.

In the phasor representation of a field, \bf E1 corresponds to the amplitude and \phi_1 to the phase, i.e.,

The first equation, the governing equation for the *Beam Envelopes* interface, can be derived by substituting the second definition of the electric field into the Helmholtz equation. If we know \phi_1, the only unknown is \bf E1 and we can solve for it. The phase, \phi_1, needs to be given *a priori* in order to solve the problem.

With the second equation, we assume a form such that the fast oscillation part, the phase, can be factored out from the field. If that’s true, the envelope \bf E1 is “slowly varying”, so we don’t need to resolve the wavelength. Instead, we only need to resolve the slow wave of the envelope. Because of this process, simulating large-scale wave optics problems is possible on personal computers.

A common question is: “When do you want the envelope rather than the field itself?” Lens simulation is one example. Sometimes you may need the intensity rather than the complex electric field. Actually, the square of the norm of the envelope gives the intensity. In such cases, it suffices to get the envelope function.

### What Happens If the Phase Function Is Not Accurately Known?

The math behind the beam envelope method introduces more questions:

- What if the phase is
*not*accurately known? - Can we use the
*Beam Envelopes*interface in such cases? - Are the results correct?

To answer these questions, we need to do a little more math.

#### 1D Example

Let’s take the simplest test case: a plane wave, Ez = \exp(-i k_0 x), where k_0 = 2\pi / \lambda_0 for wavelength \lambda_0 = 1 um, it propagates in a rectangular domain of 20 um length. (We intentionally use a short domain for illustrative purposes.)

The out-of-plane wave enters from the left boundary and transmits the right boundary without reflection. This can be simulated in the *Beam Envelopes* interface by adding a *Matched* boundary condition with excitation on the left and without excitation on the right, while adding a *Perfect Magnetic Conductor* boundary condition on the top and bottom (meaning we don’t care about the *y* direction).

The correct setting for the phase specification is shown in the figure below.

We have the answer Ez = \exp(-i k_0 x), knowing that the correct phase function is k_0 x or the wave vector is (k_0,0) *a priori*. Substituting the phase function in the second equation, we inversely get E1z = 1, the constant function.

How many mesh elements do we need to resolve a constant function? Only one! (See this previous blog post on high-frequency modeling.)

The following results show the envelope function \bf E1 and the norm of \bf E, `ewbe.normE`

, which is equal to |{\bf E1}|. Here, we can see that we get the correct envelope function if we give the exact phase function, constant one, for any number of meshes, as expected. For confirmation purposes, the phase of \bf E1z, ` arg(E1z)`

, is also plotted. It is zero, also as expected.

*The envelope function (red), the electric field norm (blue), and the phase of the envelope function (green) for the correct phase function k _{0}x, computed for different mesh sizes.*

Now, let’s see what happens if our guess for the phase function is a little bit off — say, (0.95k_0,0) instead of the exact (k_0,0). What kind of solutions do we get? Let’s take a look:

*The envelope function (red), the electric field norm (blue), and the phase of the envelope function (green) for the wrong phase function, 0.95 k _{0}x, computed for different mesh sizes.*

What we see here for the envelope function is the so-called *beating*. It’s obvious that everything depends on the mesh size. To understand what’s going on, we need a pencil, paper, and patience.

We knew the answer was Ez = \exp(-i k_0 x), but we had “intentionally” given an incorrect estimate in the COMSOL® software. Substituting the wrong phase function in the second equation, we get \exp(-i k_0 x)={\bf E1z} \exp(-0.95i k_0 x). This results in {\bf E1z} = \exp(-0.05i k_0 x), which is no longer constant one. This is a wave with a wavelength of \lambda_b= 2\pi/0.05k_0 = 20 um, which is called the *beat wavelength*.

Let’s take a look at the plot above for six mesh elements. We get exactly what is expected (red line), i.e., {\bf E1z} = \exp(-0.05i k_0 x). The plot automatically takes the real part, showing {\bf E1z} = \cos(-0.05i k_0 x). The plots for the lower resolutions still show an approximate solution of the envelope function. This is as expected for finite element simulations: coarser mesh gives more approximate results.

This shows that if we make a wrong guess for the phase function, we get a wrong (beat-convoluted) envelope function. Because of the wrong guess, the envelope function is added a phase of the beating (green line), which is -0.05i k_0 x.

What about the norm of \bf E? Look at the blue line in the plots above. It looks like the COMSOL Multiphysics software generated a correct solution for `ewbe.normE`

, which is constant one. Let’s calculate: Substituting both the wrong (analytical) phase function and the wrong (beat-convoluted) envelope function in the second equation, we get {\bf Ez} = \exp(-0.05i k_0 x) \times \exp(-0.95i k_0 x) = \exp(-i k_0 x), which is the correct fast field!

If we take a norm of \bf E, we get a correct solution, constant one. This is what we wanted. Note that we can’t display \bf E itself because the domain can be too large, but we can find \bf E analytically and display the norm of \bf E with a coarse mesh.

This is not a trick. Instead, we see that if the phase function is off, the envelope function will also be off, since it becomes beat-convoluted. However, the norm of the electric field can still be correct. Therefore, it is important that the beat-convoluted envelope function be correctly computed in order to get the correct electric field. The above plots clearly show that. The six-element mesh case gives the completely correct electric field norm because it fully resolves the beat-convoluted envelope function. The other meshes give an approximate solution to the beat-convoluted envelope function depending on the mesh size. They also do so for the field norm. This is a general consequence that holds true for arbitrary cases.

No matter what phase function we use in COMSOL Multiphysics, we are okay as long as we correctly solve the first equation for \bf E1 and as long as the phase function is continuous over the domain. When there are multiple materials in a domain, the continuity of the phase function is also critical to the solution accuracy. We may discuss this in a future blog post, but it is also mentioned in this previous blog post on high-frequency modeling.

#### 2D Example

So far, we have discussed a scalar wave number. More generally, the phase function is specified by the wave vector. When the wave vector is not guessed correctly, it will have vector-valued consequences. Suppose we have the same plane wave from the first example, but we make a wrong guess for the phase, i.e., k_0(x \cos \theta + y \sin \theta) instead of k_0 x . In this case, the wave number is correct but the wave vector is off. This time, the beating takes place in 2D.

Let’s start by performing the same calculations as the 1D example. We have \exp(-i k_0 x)= {\bf E1z}(x,y) \exp(-i k_0 (x \cos \theta+y \sin \theta) ) and the envelope function is now calculated to be {\bf E1z}(x,y) = \exp(-i k_0 (x (1-\cos \theta) -y \sin \theta) ) , which is a tilted wave propagating to direction (1-\cos \theta, -\sin \theta) , with the beat wave number k_b = 2 k_0/\sin (\theta/2) and the beat wavelength \lambda_b=\lambda_0/(2\sin (\theta/2)).

The following plots are the results for *θ* = 15° for a domain of 3.8637 um x 29.348 um for different max mesh sizes. The same boundary conditions are given as the previous 1D example case. The only difference is that the incident wave on the left boundary is {\bf E1z}(0,y) = \exp(i k_0 y \sin \theta) . (Note that we have to give the corresponding wrong boundary condition because our phase guess is wrong.)

In the result for the finest mesh (rightmost), we can confirm that \bf E1z is computed just like we analyzed in the above calculation and the norm of \bf Ez is computed to be constant one. These results are consistent with the 1D example case.

*The electric field norm (top) and the envelope function (bottom) for the wrong phase function k_0(x \cos\theta +y \sin\theta ), computed for different mesh sizes. The color range represents the values from -1 to 1.*

### Simulating a Lens Using the Beam Envelopes Interface

The ultimate goal here is to simulate an electromagnetic beam through optical lenses in a millimeter-scale domain with the *Beam Envelopes* interface. How can we achieve this? We already discussed how to compute the right solution. The following example is a simulation for a hard-apertured flat top incident beam on a plano-convex lens with a radius of curvature of 500 um and a refractive index of 1.5 (approximately 1 mm focal length).

Here, we use \phi_1 = k_0 x, which is not accurate at all. In the region before the lens, there is a reflection, which creates an interference. In the lens, there are multiple reflections. After the lens, the phase is spherical so that the beam focuses into a spot. So this phase function is far different from what is happening around the lens. Still, we have a clue. If we plot \bf E1z, we see the beating.

*Plot of \bf E1z. The inset shows the finest beat wavelength inside the lens.*

As can be seen in the plot, a prominent beating occurs in the lens (see the inset). Actually, the finest beat wavelength is \lambda_0/2 in front of the lens. To prove this, we can perform the same calculations as in the previous examples. The finest beat wavelength is due to the interference between the incident beam and reflected beam, but we can ignore this because it doesn’t contribute to the forward propagation. We can see that the mesh doesn’t resolve the beating before the lens, but let’s ignore this for now.

The beat wavelength in the lens is 3\lambda_0/2 for the backward beam and 2\lambda_0 for the forward beam for *n* = 1.5, which we can also prove in the same way as the previous examples. Again, we ignore the backward beam. In the plot, what’s visible is the 2\lambda_0 beating for the forward beam. The backward beam is only a fraction (approximately 4% for *n* = 1.5 of the incident beam, so it’s not visible). The following figure shows the mesh resolving the beat inside the lens with 10 mesh elements.

*The beat wavelength inside the lens. The mesh resolves the beat with 10 mesh elements.*

Other than the beating for the propagating beam in the lens, the beating in the subsequent air domain is pretty large, so we can use a coarse mesh here. This may not hold for faster lenses, which have a more rapid quadratic phase and can have a very short beat wavelength. In this example, we must use a finer mesh only in the lens domain to resolve the fastest beating.

The computed field norm is shown at the top of this blog post. To verify the result, we can compute the field at the lens exit surface by using the *Frequency Domain* interface, and then using the Fresnel diffraction formula to calculate the field at the focus. The result for the field norm agrees very well.

*Comparison between the* Beam Envelopes *interface and Fresnel diffraction formula. The mesh resolves the beat inside the lens with 10 mesh elements.*

The following comparison shows the mesh size dependence. We get a pretty good result with our standard recommendation, \lambda_b/6, which is equal to \lambda_0/3. This makes it easier to mesh the lens domain.

*Mesh size dependence on the field norm at the focus.*

As of version 5.3a of the COMSOL® software, the Fresnel Lens tutorial model includes a computation with the *Beam Envelopes* interface. Fresnel lenses are typically extremely thin (wavelength order). Even if there is diffraction in and around the lens surface discontinuities, the fine mesh around the lens part does not significantly impact the total number of mesh elements.

### Concluding Remarks

In this blog post, we discuss what the *Beam Envelopes* interface does “under the hood” and how we can get accurate solutions for wave optics problems. Even if we get beating, the beat wavelength can be much longer than the wavelength, which makes it possible to simulate large optical systems.

Although it seems tedious to check the mesh size to resolve beating, this is not extra work that is only required for the *Beam Envelopes* interface. When you use the finite element method, you always need to check the mesh size dependence for accurately computed solutions.

### Next Steps

Try it yourself: Download the file for the millimeter-range focal length lens by clicking the button below.

### References

- J. Goodman,
*Fourier Optics*, Roberts and Company Publishers, 2005.

## Comments

Phillip SpringerJanuary 15, 2018 10:58 amDear Mr. Yosuke Mizuyama,

I am trying to set COMSOL up to solve the Helmholtz equation in 3D. In my attempt, I basically define a cuboid where one plain shall have an initial electric field distribution and I want to know the electric field at the opposing side. My initial electric field shall be arbitrary, but for test purposes I have it as a Gaussian. This works well, as long as there are no focusing terms. When I use a complex-valued Gaussian, e.g. with a focusing term exp(i*k*(x^2+z^2)/2/R_0) [propagation direction is y], then solution does not look ok. I also tried to use a user-defined phase within the Beam Envelopes interface, but no success either. Can you give me some advice on how to propagate an (arbitrary) initial electric field distribution that may contain a focusing phase in 3D?

Yosuke MizuyamaJanuary 16, 2018 10:08 amDear Phillip,

Thank you for reading my blog. It’s a little bit hard to diagnose what went wrong from your problem description but your phase expression seems correct. What I can think of is that you might have used a small domain size compared to the beam width or a small R_0 compared to the wavelength, in which case you need to pay a special attention. I recommend using a large domain size and start from a very large R_0. If you are successful, then you can reduce them.

Best regards,

Yosuke