Introduction to subpixel averaging

We'll continue our tutorial on the FDTD method. In particular, we would be talking about subpixel averaging.

How to capture geometric features below grid size in FDTD?

In the previous lecture, we highlighted the importance of understanding how the dielectric function is distributed over the grids. As a simple illustration, we simulate the transmission of light through a dielectric slab. Naturally, the electric constant is assigned to reflect the thickness of the slab. Now suppose that the discretization level is 25 nm. In the simplest algorithm, you examine every single grid point (electric grid point here), and if the electric field grid point is inside the dielectric slab, the dielectric constant of the grid is assigned to be the dielectric value in the slab. For example, in the case of silicon, it will be 12.25. If the grid point is outside the slab, as indicated as blue dots here, the dielectric constant is assigned to be that of the air, which is 1.

However, because the grid size is always larger than 0, the variation of the slab thickness below the discretization level will not be captured. Here is an example: suppose you intend to simulate a slab of thickness of 150 nm, in the upper figure, the red dots are assigned to be silicon; however, when you increase the thickness to 175 nm where the change is right at the discretization level, the number of grid points assigned to silicon remains the same. So even though you intend to simulate two slabs of different thickness, the results are actually identical. Thus you don’t have access to the information about the thickness beyond the spatial discretization level. It will be nicer to develop a scheme that takes into account the dielectric function distribution beyond integer multiples of the grid size.

Basic idea of subpixel averaging

That's where the idea of subpixel averaging comes in. In subpixel averaging, for the grid point in the vicinity of the material boundary, one assigns the dielectric constant in a way that takes into account where the boundary lies within the unit cell.

Simplest scheme of subpixel averaging

Here we illustrate a very simple idea of what one can do: as shown in the left figure, the open region is air and the solid region is the dielectric, and there is a grid point in the vicinity of this boundary. Now let’s focus on the computational unit cell around the grid point. In the top figure where the slab thickness is 150 nm, the grid point lies right at the interface between silicon and air. Therefore, in the computational unit cell around the grid point, half is air and half is silicon. As a simplest way to average it, you can assign the dielectric constant at this grid point to be the average of silicon and air.

On the other hand, in the bottom figure where the thickness is 175 nm, by repeating the same process, you can see that in the computational unit cell around the grid point in the vicinity of the boundary, it is entirely silicon. Thus, the dielectric can be assigned to be silicon.

So with a simple formula, the dielectric constant assigned to the grid point has information about where the interface between the two materials really is. In both these cases, the permittivity assignment depends on the thickness of the layer.

Numerical demonstration of the simplest scheme of subpixel averaging

This is a very simple idea, but it actually works quite well. Here we illustrate a numerical demonstration. We consider the light propagating along the z-axis, which is perpendicular to the slab. The electric field is assumed to polarize along the x-direction which is parallel to the slab.

The simulation results are shown in the figure on the right. The blue and the green curves correspond to the two cases where we change the thickness from 175 nm to be 150 nm. As you can see, the spectrum from the FDTD simulation agrees actually quite well with the analytical results. Also the spectrum from the FDTD simulation for the two thicknesses are different even though the change in the thickness is far below the length scale of the grid size. Thus this simple formula actually works quite well.

Now perhaps fortunately or unfortunately, this simple idea in fact only works in this particular setup where the electric field is parallel to the interface. To illustrate that, it will be helpful to go through a bit of a rigorous derivation of where subpixel averaging formulas come from.

Field update procedure during time stepping

In FDTD, a typical field update procedure is shown in the figure. One updates both the D-field and the B-field using the Maxwell equation that contains the curl operations; and then converts the D-field to the E-field through dielectric permittivity, as well as B-field to H-field through magnetic permeability.

The role of subpixel averaging in time stepping

To do subpixel averaging, one would need to look carefully at both of these curl equations in order to derive the effective dielectric function near the interface. This would allow a conversion from the D-field to the E-field that takes into account the actual geometry of the structure.

Derivation of subpixel averaging for tangential components with Ampere’s law

In our simple case with one-dimensional configuration, the electric field is polarized along the x-direction, and propagating along the z-direction. In this case, the only relevant E-grid points are those for the $$E_x$$-fields as indicated by the red dots. Likewise, the grids for the $$H_y$$-fields are indicated by the blue dots. The curl equation that relates the update of D-field to the curl of H-field can be converted to an integral form. In this form, the D-field is now spatially integrated along the line that connects the two H-field components.

Up to this point at the FDTD level, this is an exact representation of the curl equation. The next step is to convert the D-field to the E-field in order to update the B-field in the next step.

Interface condition in deriving subpixel averaging formula

An important observation is that the $$E_x$$-field is continuous because it's parallel to the interface. Consequently, to carry out the spatial integration of the $$D_x$$-field in a small enough unit cell, one can assume $$E_x$$ to be uniform, and take it out of the integration. This leads to the effective dielectric formula: it is the spatial average of the dielectric function within the unit cell. That is exactly the formula that we have naively guessed and tested.

Limitation of the simple subpixel averaging formula

To repeat, the formula here is a simple spatial average of the dielectric function within the unit cell surrounding the grid point. Note that the derivation here strongly relies upon the fact that the $$E_x$$-field is continuous parallel to the interface. In other words, it works well only in the case where the dominant electric field is parallel to the interface. As you can see in the figure, it does work well.

Subpixel averaging formula is more complicated in 2D and 3D

The derivation here points to the fact that in general, subpixel averaging in fact is quite a bit more complicated. In the two dimensional case, the subpixel averaging would be different depending on the polarization: the s-polarization where the electric field is parallel to the interface, and the p-polarization where the electric field has both the parallel and perpendicular field components to the interface. The subpixel formulas are different for the two polarizations. In the three dimensional case, the interface and the field in general can have arbitrary orientation with respect to each other, so a tensorial average is needed in order to do correct subpixel averaging. This is the algorithm that we have implemented in Tidy3D.

To summarize the discussion, I hope to give you an illustration of the basic idea behind subpixel averaging. I also show perhaps the simplest example where you can take a simple spatial average of the dielectric function. But I also want to point to some of the subtleties associated with subpixel averaging that are related to the field component and the surface orientation. When implemented correctly, it is a very powerful algorithm in the FDTD code because it allows you to simulate certain geometric features at the length scale smaller than the discretization level.