Finite Difference Time Domain (FDTD) solver introduction
The Finite-Difference Time-Domain (FDTD) method[1,2,3] is a state-of-the-art method for solving Maxwell's equations in complex geometries. Being a direct time and space solution, it offers the user a unique insight into all types of problems in electromagnetics and photonics. In addition, FDTD can also obtain the frequency solution by exploiting Fourier transforms, thus a full range of useful quantities can be calculated, such as the complex Poynting vector and the transmission / reflection of light.
Solver Physics
This section will introduce the basic mathematical and physics formalism behind the FDTD algorithm.
FDTD solves Maxwell's curl equations in non-magnetic materials:
∂→D∂t=∇×→H→D(ω)=ε0εr(ω)→E(ω)∂→H∂t=−1μ0∇×→E
where H, E, and D are the magnetic, electric, and displacement fields, respectively, while ϵr(ω)
is the complex relative dielectric constant (ϵr(ω)=n2
, where n is the refractive index).
In three dimensions, Maxwell equations have six electromagnetic field components: Ex, Ey, Ez and Hx, Hy, and Hz. If we assume that the structure is infinite in the z dimension and that the fields are independent of z, specifically that
εr(ω,x,y,z)=εr(ω,x,y)∂→E∂z=∂→H∂z=0
then Maxwell's equations split into two independent sets of equations composed of three vector quantities each which can be solved in the x-y plane only. These are termed the TE (transverse electric), and TM (transverse magnetic) equations. We can solve both sets of equations with the following components:
TE: Ex, Ey, Hz
TM: Hx, Hy, Ez
For example, in the TM case, Maxwell's equations reduce to:
∂Dz∂t=∂Hy∂x−∂Hx∂yDz(ω)=ε0εr(ω)Ez(ω)∂Hx∂t=−1μ0∂Ez∂y∂Hy∂t=1μ0∂Ez∂x
The FDTD method solves these equations on a discrete spatial and temporal grid. Each field component is solved at a slightly different location within the grid cell (Yee cell), as shown below. By default, data collected from the FDTD solver is automatically interpolated to the origin of each grid point, so the end user does not have to deal with this issue in their analysis.
Dispersive materials with tabulated refractive index (n,k) data as a function of wavelength can be incorporated by using the multi-coefficient material models that automatically generates a material model based on the tabulated data. Alternatively, specific models such as Plasma (Drude), Debye or Lorentz can be used. See the Materials section for details. The FDTD solver supports a range of boundary conditions, such as PML, periodic, and Bloch. See the boundary conditions section here for the complete list. The FDTD solver supports a number of different types of sources such as point dipoles, beams, plane waves, a total-field scattered-field (TFSF) source, a guided-mode source for integrated optical components, and an imported source to interface with external photonic design softwares. Detailed information about each of these sources is contained in the sources section.
You can find more information on Wikipedia.
Meshing
FDTD uses a rectangular, Cartesian style mesh, like the one shown in the following screenshot. It's important to understand that of the fundamental simulation quantities (material properties and geometrical information, electric and magnetic fields) are calculated at each mesh point. Obviously, using a smaller mesh allows for a more accurate representation of the device, but at a substantial cost. As the mesh becomes smaller, the simulation time and memory requirements will increase. The FDTD solver provides a number of features, including the conformal mesh algorithm, that allow you to obtain accurate results, even when using a relatively coarse mesh.
By default, the simulation mesh is automatically generated. To maintain accuracy, the meshing algorithm will create a smaller mesh in high index (to maintain a constant number of mesh points per wavelength) and highly absorbing (resolve penetration depths) materials. In some cases, it is also necessary to manually add additional meshing constraints. Usually, this involves forcing the mesh to be smaller near complex structures (often metal) where the fields are changing very rapidly.
Mesh refinement options, such as conformal mesh, can allow the solver to maintain high accuracy without the need for a very small mesh.
Integrating over lines, surfaces and volumes
The electric and magnetic fields are recorded on the finite-difference mesh, as shown below for a 2D monitor, where the grey dots represent the positions where the fields are recorded. The thick black outline shows the limits of the surface monitor as seen in the Layout Editor. This monitor has an x-span of 12dx and a y-span of 5dy. This monitor records a total of 13x6 data points.
A typical calculation with this monitor might be to integrate the total power flow across the surface of the monitor, or
Power(w)=12∫real{→P(w)}⋅d→S
In order to calculate this quantity, we provide the scripting function integrate. If Pz is a variable of dimension 13x6x1, and x and y are the corresponding position vectors, then the desired quantity is:
power = 0.5*integrate(real(Pz),1:2,x,y);
The integrate script command can be used to integrate over spatial dimensions even when several frequency points have been recorded. For example, if Pz is a variable of dimension 13x6x1x10, representing 10 frequency points, then the following can be used to integrate over y (ie dimension 2) and then x (ie dimension 1)
power = 0.5*integrate(real(Pz),1:2,x,y);
power will be a matrix of dimension 10x1.
Units and normalization
Unless otherwise specified, all quantities are returned in SI units. Please see Units and normalization for more information.
References
1. Dennis M. Sullivan, Electromagnetic simulation using the FDTD method. New York: IEEE Press Series, (2000).
2. Allen Taflove, Computational Electromagnetics: The Finite-Difference Time-Domain Method. Boston: Artech House, (2005).
3. Stephen D. Gedney, Introduction to the Finite-Difference Time-Domain (FDTD) Method for Electromagnetics. Morgan & Claypool publishers, (2011).
Understanding ray vs wave optics with a simple example
This very simple example of a grating demonstrates how the ray approximation breaks down as the structure feature size becomes wavelength scale.
The structure in this example is a simple grating with an index contrast of 1.5:1. If the illumination is at normal incidence, a very simple application of Snell's law suggests that the fields should be completely reflected, since the angle of incidence on the grating surface is 45 degrees which is beyond the critical angle of 41.8 degrees.
The same system setup in FDTD.
Notice that the span is 20um.
Simulation results when the source wavelength is 0.4um. The key point is to recognize that 0.4um is much smaller than the period of the structure (20um). In this situation, the behavior of the system is basically predicted by the simple ray approximation; most of the fields are reflected.
switchtolayout;
save("400nm.fsp");
setnamed("source","center wavelength",400e-9);
run;
Simulation results when the source wavelength is 4um. The key point is to recognize that 4um is much closer in scale to the period of the structure (20um). In this situation, the behavior of the system is quite different. A significant fraction of power is transmitted through the device.
switchtolayout;
save("4000nm.fsp");
setnamed("source","center wavelength",4000e-9);
run;
Units and normalization conventions in Lumerical solvers
FDTD MODE DGTD CHARGE HEAT FEEM
Units for electrical and thermal solvers
Unless otherwise stated, Lumerical uses SI units.
| Quantity | Description | Units | Unit description |
|---|---|---|---|
| L | Length units in semiconductor models This is reflected in the semiconductor device literature, and in the parameter coefficients for the material models. | cm | Centimeter |
| E | Energy The electron energy E is related to the local electrostatic potential (voltage) as E = -qV. All energies (and voltages) are referenced from the (equilibrium) Fermi level of an electrical contact in the system. | eV | Electron volt |
| n | Electron density | 1/cm3 | Per centimeter cube |
| p | Hole density | 1/cm3 | Per centimeter cube |
| Jn | Electron current density | A/cm2 | Ampere per centimeter square |
| Jp | Hole current density | A/cm2 | Ampere per centimeter square |
| N | Net doping density The doping density is positive for p-type (acceptor) dopants and negative for n-type (donor) dopants | 1/cm3 | Per centimeter cube |
| E | Electric field | V/m | Volt per meter |
| V | Electrostatic potential (voltage) | V | Volt |
Units for optical solvers
Unless otherwise stated, Lumerical's optical solvers used SI units at all times.
General
| Quantity | Description | Units | Unit description |
|---|---|---|---|
| f=w/2p | Frequency | Hz | Hertz |
| x,y,z | Position | m | Meter |
| t | Time | s | Seconds |
Time domain electromagnetic fields
| Quantity | Description | Units | Unit description |
|---|---|---|---|
| E(t) | Electric field as function of time | V/m | Volts per meter |
| |E(t)|2 | Electric field intensity as a function of time | (V/m)2 | Volts squared per meter squared |
| H(t) | Magnetic field as a function | A/m | Amperes per meter |
| |H(t)|2 | Magnetic field intensity as a function of time | (A/m)2 | Amperes squared per meter squared |
| P(t) | Poynting vector as a function of time | W/m2 | Watts per meter squared |
| Power(t) | Power as a function of time | W | Watts |
Dipole moments
| Quantity | Description | Units | Unit description |
|---|---|---|---|
| p | Electric dipole in 3D | Cm | Coulomb meters |
| m | Magnetic dipole in 3D | Am2 | Ampere meters squared |
| p | Electric field in 2D | Cm/m | Coulomb meters per meter |
| m | Magnetic dipole in 2D | Am2/m | Ampere meters squared per meter |
Frequency domain electromagnetic fields - Steady state, single frequency, cwnorm data
| Quantity | Description | Units | Unit description |
|---|---|---|---|
| E(w) | Electric field as a function of angular frequency | V/m | Volts per meter |
| |E(w)|2 | Electric field intensity as a function of angular frequency | (V/m)2 | Volts squared per meter squared |
| H(w) | Magnetic field as a function of angular frequency | A/m | Amperes per meter |
| |H(w)|2 | Magnetic field intensity as a function of angular frequency | (A/m)2 | Amperes squared per meter squared |
| P(w) | Poynting vector as a function of angular frequency | W/m2 | Watts per meter squared |
| Power(w) | Power as a function of angular frequency | W | Watts |
| Power(w) | 2D Power as a function of angular frequency | W/m | Watts per meter |
Frequency domain electromagnetic fields - nonorm data
| Quantity | Description | Units | Unit description |
|---|---|---|---|
| E(w) | Electric field as a function of angular frequency | V/m/Hz | Volts per meter per Hertz |
| |E(w)|2 | Electric field intensity as a function of angular frequency | (V/m/Hz)2 | Volts squared per meter squared per Hertz squared |
| H(w) | Magnetic field as a function of angular frequency | A/m/Hz | Amperes per meter per Hertz |
| |H(w)|2 | Magnetic field intensity as a function of angular frequency | (A/m/Hz)2 | Amperes squared per meter squared per Hertz squared |
| P(w) | Poynting vector as a function of angular frequency | W/m2/Hz2 | Watts per meter squared per Hertz squared |
| Power(w) | Power as a function of angular frequency | W/Hz2 | Watts per Hertz squared |
| Power(w) | 2D Power as a function of angular frequency | W/Hz2/m | Watts per Hertz squared per meter |
Source amplitudes
Beam sources
When specifying the amplitude for beam sources, the "amplitude" refers to the peak electric field amplitude in units of V/m. For example, if a Gaussian beam has the following electric field distribution in time and space:
E(x,y,z,t)=E0sin(ω0(t−t0))exp(−(t−t0)22(Δt)2)exp(−(x2+y2)w20)
Then the "amplitude" refers to the value of E0 and has units of V/m. It is worth noting that different beams will inject different amounts of power for a given source amplitude.
Dipole sources
For dipole sources, amplitude refers to the amplitude of the point source whose units are listed below. Base amplitude refers to the amplitude that will generate a radiated CW power of 10 nW/m in 2D simulations and 1 fW in 3D simulations, and total amplitude refers to the amplitude actually used in the simulations which is the product of the amplitude and the base amplitude.
Dipole source amplitude units are
- Cm for 3D electric dipole sources
- Am2 for 3D magnetic dipole sources
- Cm/m for 2D electric dipole sources
- Am2/m for 2D magnetic dipole sources
Understanding frequency domain CW normalization
The frequency domain field monitors in the FDTD and varFDTD solvers record the electric and magnetic fields at a series of user-defined frequencies. These can be returned in either the Continuous Wave Normalization state (cwnorm), or the No Normalization state (nonorm). For most applications, the default cwnorm state is the best choice.
In the nonorm state, the returned fields are simply the Fourier transform of the simulated time domain fields, and we use the subscript sim to refer to these fields in the table below. In the cwnorm state, the fields are normalized by the Fourier transform of the source pulse, thereby yielding the impulse response of the system, and we use a subscript imp to refer to these fields in the table below.
| Quantity | Definition | Normalization state |
|---|---|---|
| Esim(w) | →Esim(ω)=∫exp(iωt)→E(t)dt |
| nonorm | |
| Hsim(w) | →Hsim(ω)=∫exp(iωt)→H(t)dt |
| nonorm | |
| Psim(w) | →Psim(ω)=→Esim(ω)×→H∗sim(ω) |
| nonorm | |
| Eimp(w) | →Eimp(ω)=1s(ω)∫exp(iωt)→E(t)dt=→Esim(ω)s(ω) |
| cwnorm | |
| Himp(w) | →Himp(ω)=1s(ω)∫exp(iωt)→H(t)dt=→Hsim(ω)s(ω) |
| cwnorm | |
| Pimp(w) | →Pimp(ω)=→Eimp(ω)×→H∗imp(ω)=→Esim(ω)×→H∗sim(ω)|s(ω)|2 |
| cwnorm | |
| s(w) | s(ω)=1N∑sources∫exp(iωt)sj(t)dt |
| where sj(t) is the source time signal of the jth source and N is the number of active sources in the simulation volume | N/A |
Understanding Continuous Wave Normalization (cwnorm)
Impulse response
FDTD is a time domain method, ie. the electromagnetic fields are calculated as a function of time. Here, the system being simulated is excited by a dipole, beam, mode or imported source, and the time signal of the source, s(t), is a pulse. For example, this could be
s(t)=sin(ω0(t−t0))exp(−(t−t0)22(Δt)2)
and the Fourier transform of s(t) is s(w)
s(ω)=∫exp(iωt)s(t)dt
Ideally, s(t) would be a dirac delta function (in which case s(w) = 1). This would allow us to obtain the response of the system at all frequencies from a single simulation. For a variety of reasons, it is more efficient and numerically accurate to excite the system with a short pulse such that the spectrum, |s(w)|2, has a reasonably large value over all frequencies of interest.
In the nonorm state, power and profile monitors return the response of the system to the simulated input pulse s(t):
→Esim(ω)=∫exp(iωt)→E(t)dt
The simulated electric field as a function of angular frequency, Esim(w), depends on both the source pulse used, s(t), and the system under study.
In the default cwnorm state, power and profile monitors return the impulse response of the system.
→Eimp(ω)=→Esim(ω)s(ω)
The impulse response of the system is a much more useful quantity because it is completely independent of the source pulse used to excite the system.
Example
Consider a beam source injected into free space at z=z0. The source signal is
s(t)=sin(ω0(t−t0))exp(−(t−t0)22(Δt)2)
The electric field at the source injection plane has the following form:
E(x,y,z0,t)=E0(x,y,z0)s(t)
In the cwnorm state,
E(x,y,z,ω)=E0(x,y,z)
In other words, the field returned in the cwnorm state is the field that would exist if a CW source of amplitude E0 had been used at the angular frequency w. It removes any frequency dependence due to the finite pulse length of the source, and the units of the returned fields are the same as time domain fields.
Spectral averaging introduction and definitions
For some FDTD and varFDTD simulations, it is useful to calculate the incoherent spectral average of the electromagnetic fields or the Poynting vector. For example, we might want to calculate
〈→P〉=∫W(ω)→Pimp(ω)dω
where W(w) is a weighting function and Pimp(w) is the Poynting vector returned in the cwnorm state, in other words, the impulse response of the system. (Please see Frequency domain normalization for an explanation of the impulse response.) We could calculate the Pimp(w) at many different angular frequencies and then perform the spectral average after the simulation. FDTD, however, allows for two methods of more efficient spectral averaging during the simulation which can significantly reduce the memory requirements and increase the speed of the simulation. These methods are called "total spectral averaging" and "partial spectral averaging".
The frequency power and frequency profile monitors can record spectral averages of various quantities during the simulation. Like the other quantities calculated by these monitors, these can be returned in the Continuous Wave Normalization state (cwnorm), or the No Normalization state (nonorm). For most applications, the default cwnorm state is the best choice. In the tables below, we use the subscripts sim and imp to refer to the quantities returned in the nonorm and cwnorm states respectively. Please refer to Frequency domain normalization for more details on the two normalization states, and obtaining the impulse response of the system.
Total spectral averaging
FDTD/Propagator supports a spectral average that uses the source input spectrum as the weighting function. Total spectral averaging can be useful when the source spectrum of the simulation matches the actual illumination conditions. The following table gives the precise definitions of the quantities available using total spectral averaging.
| Quantity | Definition | Normalization state |
|---|---|---|
| 〈∣→Esim∣2〉total |
| 〈∣→Esim∣2〉total=∫+∞0∣→Esim(ω)∣2dω |
| nonorm | |
| 〈∣→Hsim∣2〉total |
| 〈∣→Hsim∣2〉total=∫+∞0∣→Hsim(ω)∣2dω |
| nonorm | |
| 〈→Psim〉total |
| 〈→Psim〉total=∫+∞0→Psim(ω)dω |
| nonorm | |
| 〈∣→Eimp∣2〉total |
| 〈∣→Eimp∣2〉total=∫+∞0∣s(ω)∣2∣→Eimp(ω)∣2dω∫+∞0∣s(ω)∣2dω |
| cwnorm | |
| 〈∣→Himp∣2〉total |
| 〈∣→Himp∣2〉total=∫+∞0∣s(ω)∣2∣→Himp(ω)∣2dω∫+∞0∣s(ω)∣2dω |
| cwnorm | |
| 〈→Pimp〉total |
| 〈→Pimp〉total=∫+∞0∣s(ω)∣2→Pimp(ω)dω∫+∞0∣s(ω)∣2dω |
| cwnorm | |
| s(ω) |
| s(ω)=1N∑sources∫exp(iωt)sj(t)dt |
| , where sj(t) is the source time signal of the jth source and N is the number of active sources in the simulation volume | N/A |
Total spectral averaging can be turned on by selecting total spectral average as shown in the screenshot below.
To see the data available in a monitor, use the script command
?getdata("monitorname");
All data can be obtained in scripting, with commands such as getdata and transmission_avg. The available components for total spectral averaging are shown in the following table.
| Quantity | Data name | Example script command |
|---|---|---|
| 〈∣Ex,sim∣2〉total |
| E2x_avg | nonorm; E2x = getdata("monitor1","E2x_avg"); | |
| 〈∣Ey,sim∣2〉total |
| E2y_avg | nonorm; E2y = getdata("monitor1","E2y_avg"); | |
| 〈∣Ez,sim∣2〉total |
| E2z_avg | nonorm; E2z = getdata("monitor1","E2z_avg"); | |
| 〈∣Hx,sim∣2〉total |
| H2x_avg | nonorm; H2x = getdata("monitor1","H2x_avg"); | |
| 〈∣Hy,sim∣2〉total |
| H2y_avg | nonorm; H2y = getdata("monitor1","H2y_avg"); | |
| 〈∣Hz,sim∣2〉total |
| H2z_avg | nonorm; H2z = getdata("monitor1","H2z_avg"); | |
| 〈Px,sim〉total |
| Px_avg | nonorm; Px = getdata("monitor1","Px_avg"); | |
| 〈Py,sim〉total |
| Py_avg | nonorm; Py = getdata("monitor1","Py_avg"); | |
| 〈Pz,sim〉total |
| Pz_avg | nonorm; Pz = getdata("monitor1","Pz_avg"); | |
| 〈∣Ex,imp∣2〉total |
| E2x_avg | cwnorm; E2x = getdata("monitor1","E2x_avg"); | |
| 〈∣Ey,imp∣2〉total |
| E2y_avg | cwnorm; E2y = getdata("monitor1","E2y_avg"); | |
| 〈∣Ez,imp∣2〉total |
| E2z_avg | cwnorm; E2z = getdata("monitor1","E2z_avg"); | |
| 〈∣Hx,imp∣2〉total |
| H2x_avg | cwnorm; H2x = getdata("monitor1","H2x_avg"); | |
| 〈∣Hy,imp∣2〉total |
| H2y_avg | cwnorm; H2y = getdata("monitor1","H2y_avg"); | |
| 〈∣Hz,imp∣2〉total |
| H2z_avg | cwnorm; H2z = getdata("monitor1","H2z_avg"); | |
| 〈Px,imp〉total |
| Px_avg | cwnorm; Px = getdata("monitor1","Px_avg"); | |
| 〈Py,imp〉total |
| Py_avg | cwnorm; Py = getdata("monitor1","Py_avg"); | |
| 〈Pz,imp〉total |
| Pz_avg | cwnorm; Pz = getdata("monitor1","Pz_avg"); |
Partial spectral averaging
FDTD/Propagator supports a spectral average that uses a Lorentzian weighting function multiplied by the source spectrum. Partial spectral averaging is useful to extract the average response of the system to a variety of different illumination conditions from a single simulation. The following table gives the precise definitions of the quantities available using partial spectral averaging.
| Quantity | Definition | Normalization state |
|---|---|---|
| 〈∣→Esim(ω)∣2〉partial |
| 〈∣→Esim(ω)∣2〉partial=∫+∞−∞∣h(ω,ω′)∣2∣→Esim(ω′)∣2dω′ |
| nonorm | |
| 〈∣→Hsim(ω)∣2〉partial |
| 〈∣→Hsim(ω)∣2〉partial=∫+∞−∞∣h(ω,ω′)∣2∣→Hsim(ω′)∣2dω′ |
| nonorm | |
| 〈→Psim(ω)〉partial |
| 〈→Psim(ω)〉partial=∫+∞−∞∣h(ω,ω′)∣2→Psim(ω′)dω′ |
| nonorm | |
| 〈∣→Eimp(ω)∣2〉partial |
| 〈∣→Eimp(ω)∣2〉partial=∫+∞−∞∣h(ω,ω′)∣2∣s(ω′)∣2∣→Eimp(ω′)∣2dω′∫+∞−∞∣h(ω,ω′)∣2∣s(ω′)∣2dω′ |
| cwnorm | |
| 〈∣→Himp(ω)∣2〉partial |
| 〈∣→Himp(ω)∣2〉partial=∫+∞−∞∣h(ω,ω′)∣2∣s(ω′)∣2∣→Himp(ω′)∣2dω′∫+∞−∞∣h(ω,ω′)∣2∣s(ω′)∣2dω′ |
| cwnorm | |
| 〈→Pimp(ω)〉partial |
| 〈→Pimp(ω)〉partial=∫+∞−∞∣h(ω,ω′)∣2∣s(ω′)∣2→Pimp(ω′)dω′∫+∞−∞∣h(ω,ω′)∣2∣s(ω′)∣2dω′ |
| cwnorm | |
| ∣h(ω,ω′)∣2 |
| ∣h(ω,ω′)∣2=δ(ω−ω′)+(πδ)2 |
∫∣h(ω,ω′)∣2dω′=1
| N/A | |
| s(ω) |
| s(ω)=1N∑sources∫exp(iωt)sj(t)dt |
| , where sj(t) is the source time signal of the jth source and N is the number of active sources in the simulation volume | N/A |
Partial spectral averaging can be turned on by selecting partial spectral average as shown in the screenshot below.
The FWHM of |h(f,f')|2 is called d, and can be modified as shown below.
Please note that when considered as a function of angular frequency, the FWHM of |h(w,w')|2 is 2pd rad/seconds.
To see the data available in a monitor, use the script command
?getdata("monitorname");
All data can be obtained in scripting, with commands such as getdata and transmission_pavg. The available components for total spectral averaging are shown in the following table.
| Quantity | Data name | Example script command |
|---|---|---|
| 〈∣Ex,sim∣2〉partial |
| E2x_pavg | nonorm; E2x = getdata("monitor1","E2x_pavg"); | |
| 〈∣Ey,sim∣2〉partial |
| E2y_pavg | nonorm; E2y = getdata("monitor1","E2y_pavg"); | |
| 〈∣Ez,sim∣2〉partial |
| E2z_pavg | nonorm; E2z = getdata("monitor1","E2z_pavg"); | |
| 〈∣Hx,sim∣2〉partial |
| H2x_pavg | nonorm; H2x = getdata("monitor1","H2x_pavg"); | |
| 〈∣Hy,sim∣2〉partial |
| H2y_pavg | nonorm; H2y = getdata("monitor1","H2y_pavg"); | |
| 〈∣Hz,sim∣2〉partial |
| H2z_pavg | nonorm; H2z = getdata("monitor1","H2z_pavg"); | |
| 〈Px,sim〉partial |
| Px_pavg | nonorm; Px = getdata("monitor1","Px_pavg"); | |
| 〈Py,sim〉partial |
| Py_pavg | nonorm; Py = getdata("monitor1","Py_pavg"); | |
| 〈Pz,sim〉partial |
| Pz_pavg | nonorm; Pz = getdata("monitor1","Pz_pavg"); | |
| 〈∣Ex,imp∣2〉partial |
| E2x_pavg | cwnorm; E2x = getdata("monitor1","E2x_pavg"); | |
| 〈∣Ey,imp∣2〉partial |
| E2y_pavg | cwnorm; E2y = getdata("monitor1","E2y_pavg"); | |
| 〈∣Ez,imp∣2〉partial |
| E2z_pavg | cwnorm; E2z = getdata("monitor1","E2z_pavg"); | |
| 〈∣Hx,imp∣2〉partial |
| H2x_pavg | cwnorm; H2x = getdata("monitor1","H2x_pavg"); | |
| 〈∣Hy,imp∣2〉partial |
| H2y_pavg | cwnorm; H2y = getdata("monitor1","H2y_pavg"); | |
| 〈∣Hz,imp∣2〉partial |
| H2z_pavg | cwnorm; H2z = getdata("monitor1","H2z_pavg"); | |
| 〈Px,imp〉partial |
| Px_pavg | cwnorm; Px = getdata("monitor1","Px_pavg"); | |
| 〈Py,imp〉partial |
| Py_pavg | cwnorm; Py = getdata("monitor1","Py_pavg"); | |
| 〈Pz,imp〉partial |
| Pz_pavg | cwnorm; Pz = getdata("monitor1","Pz_pavg"); |
Using Parseval's theorem to check for energy conservation between the time and frequency domain
Parseval's theorem refers to that information is not lost in Fourier transform. In this example, we verify energy conservation between time and frequency domain results from an FDTD simulation using Parseval's theorem. This is done by evaluating the energy carried by a short pulse both in the time and frequency domain. In this example, we employ the "power" rawdata returned by both planar time and frequency monitor for calculation.
| Note: Energy conservation - nonorm state This example only applies to the nonom state. In nonorm, the source is a pulse by default. The pulse carries a finite amount of energy, in Joule. If there is no loss in Fourier transform, the amount of energy has to be exactly the same in time and frequency domain. Unlike a CW source, the amount of energy accumulated is a function of time. This example does not apply when CWnorm is on. |
Suppose that e(t) and E(f) are the fields in the time and frequency domain respectively, where E(f) is obtained by the fourier transform of e(t). According to Parseval's theorem, the following equation holds,
∫∞−∞|e(t)|2dt=∫∞−∞|E(f)|2df
For a simple plane wave, the Poynting vector(P
) is directly proportional to |E|2
, where E is the electric field. This relation holds in both time and frequency domain. For a short pulse in the nonorm state (without time averaging),
Power –time(t)=∫real(P⊥(t))dS=n√ϵ0μ0∫|e(t)2|dS
Energy –spectrum(f)=∫real(P⊥(t))dS=n√ϵ0μ0∫|E(f)2|dS
In other words, for the same area dS, the theorem can be written as a form of conservation of energy. Note that Energy_spectrum(f) has a unit of Watt/Hz^2, while power_time(t) has a unit of Watt.
Energy=∫∞−∞power(t)dt=∫∞−∞Energy –spectrum(f)df
To verify the above equation, a simple simulation with a plane wave and planar monitors (both time and frequency) is performed in the nonorm state. Download the above associated files and run the script. The script contains two parts: 1) the "power" returned by the monitors and 2) Poynting vector integration.
First, it extracts the "power" recorded by the time and frequency monitor accordingly. The reason for this step is because the returned "power" is not computational expensive since it's just a 1D vector. One can integrate power_time(t) and energy_spectrum(f) to calculate the amount of energy carried by a pulse, in Joule. Note that, the "power" returned by the time monitor is the instantaneous power, while the "power" returned by the frequency monitor is the time averaged power (therefore needs a factor of 2 to compensate the 0.5 set by default, see the sourcepower command). The sourcepower command is also used for comparison, where another factor of 2 is introduced from the integration of the negative frequencies due to Fourier transform.
By integrating the area under the above curves,
Ratio of energy_time to energy_f = 1.00111 Ratio of energy_time to energy_f_sourcepower = 0.99886
Second, an integration of Poynting vector is performed. Note that, this step could be quite computational expensive since the field data is usually a 3D matrix, especially the time domain data. Pz is used to calculate power in the both domains, using the above equation.
Ratio of energy_time_manual to energy_f_manual = 1.00111
This is an example to show that the time and frequency domain monitors record the same data without losing information during Fourier transform.
Far field projections in FDTD overview
The near to far field projection functions allow the EM fields to be calculated anywhere in the far field. The functions require the spatial near field EM data (typically from an FDTD simulation) and the desired far field spatial locations.
A simple way to understand far field projections is to view them as a decomposition of the near field data using a set of plane waves propagating at different angles as the basis for the decomposition. The end result is that the far field projections functions provide a straightforward, accurate, and numerically efficient method for calculating the EM fields anywhere in the intermediate and far field regions.
If your structure is periodic, consider using the related grating projections. If you are using DGTD, see Far field projections in DGTD overview.
Near to far field projection simple example
The example below shows how the far field projection can be used to see the angular distribution of reflected light from a small surface feature. A beam propagating downward at an angle reflects off the surface and bump. The reflected light is recorded in the near field slightly above the surface directly from the FDTD simulation. After the simulation is complete, the far field projection can be used to project the reflected radiation to an arbitrary location. In this case, the fields are calculated on a hemispherical surface located 1m away. The final result is an image of the field intensity on the hemisphere, as viewed from above. The effects of the misalignment of the spot can clearly be seen from the far field projection results.
Near field simulation results
A) Screenshot of simulation, B) Reflected near field |E|^2 data on XY plane above structure
Result from far field projection
A) Hemispherical surface where far field is calculated (1m radius), B) View of hemisphere, looking down
Far field |E|^2 on hemisphere
Understanding when far field projections can be used
Requirements:
- The EM fields must be known everywhere on a plane or on a closed surface.
- The plane or closed surface must be in a single homogeneous material.
- The material must extend out to infinity. There can be no additional structures (or sources) beyond this plane or closed surface.
- The far field projection functions assume that the EM fields are zero beyond the edge of the monitor. This effectively truncates the near fields at the monitor edge. Use spatial filtering to minimize this issue.
- The material can be dispersive as long as the loss is negligible, ie k<<n. However, the loss is not taken into account for the projection to the far field.
There are two types of situations where these conditions are satisfied:
Case 1 - fields known on a plane
When fields are known on a single surface, the projection functions can be used to calculate the fields anywhere beyond that surface. In the above screenshot, a monitor is located above the source. This monitor records all of the reflected fields. There are no additional structures or sources of light above the source. In this situation, the far field projections can be used to calculate the EM fields anywhere above the monitor plane.
Case 2 - fields known on a closed surface
When the fields are known on a closed surface, the projection functions can be used to calculate the fields anywhere beyond that surface. In the above screenshot, a monitor surrounds the source and scattering particle. The monitor records all reflected fields. There are no additional structures or sources of light above the source. In this situation, the far field projections can be used to calculate the EM fields anywhere outside of the monitor object.
Related publications
Allen Taflove, Computational Electromagnetics: The Finite-Difference Time-Domain Method. Boston: Artech House, (2005).
Simple far field projection example
We will calculate the far field distribution of a Gaussian beam propagating at 30 degrees to the z-axis with an azimuthal angle of 15 degrees from a near field FDTD simulation. This example is intended to help new users understand how to calculate and understand results from far field projections.
The simulation file is a simple simulation with a Gaussian source propagating at 30 degrees to the z-axis with an azimuthal angle of 15 degrees.
After running the simulation, right click on the monitor 'z2' and select Visualize -> farfield from the context menu. The following window will appear which allows you to select the frequency to project if more than 1 frequency point was recorded by the monitor.
You can click the "Far field settings" button from the bottom left corner of the frequency selection window to open the following window which allows you to set additional settings.
| NOTE: Changing the frequency and far field settings The window to select frequency and far field settings can always be opened again by right-clicking on the "farfield" result for the monitor under "Result View". This can be used to recalculate the far field projection without rerunning the FDTD simulation. |
Far field settings - details
General
- PROJECTION DIRECTION: This can be set to auto, forward, or backward. Auto will project in the field of maximum power flow, forward will project towards the positive axis direction, and backward will project to the negative axis direction, similarly to the 1 and -1 projection directions respectively when using the farfield3d script command.
- MATERIAL INDEX: This is the refractive index of the medium to use for projection. For more details see Changing the far field refractive index.
- FAR FIELD FILTER: The far field filter alpha parameter. For more details see Far field projections - Spatial filtering.
Resolution
- 2D: Number of points in the far field to project to for 2D simulations.
- 3D: Number of points in the far field to project to in each ux and uy for 3D simulations. For details on the coordinate system see Far field projections - Direction unit vector coordinates.
Assume structure is periodic
The settings in this section are discussed further in Far field projections - Periodic structures.
Near field sampling rate
- OVERRIDE NEAR FIELD MESH: If selected, enables spatial down sampling of the data in the near field that is used in the far field projection. This can be used to speed up far field projection calculations.
- NEAR FIELD SAMPLES PER WAVELENGTH: Target number of samples per wavelength to use for near field spatial down sampling of monitor data.
After clicking OK, the visualizer window will show up. Notice that you can choose to visualize the far field intensity E2 (that is |E|2), or the complex vector field components Es and Ep.
It is also possible to calculate the far field with a few script commands. The script file will perform a projection to the far field using the farfield3d function, which returns |E|2. The plot that is created is shown below. This plot represents the electric field intensity (|E|2) on a hemisphere of radius R = 1m for z > 0.
| Note: To obtain the far field vector components in spherical coordinates, use the farfieldpolar3d script function. Es corresponds to Ephi (Ez in 2D), while Ep corresponds to Etheta. |
| Note: The projection direction is automatically determined by the direction of net power flow through the monitor.
As we expect, the Gaussian beam is propagating at 30 degrees to the z axis with an azimuthal angle of approximately 15 degrees. |
Understanding direction unit vector coordinates in far field projections
This section describes the direction unit vector coordinates used by far field projections and grating projections.
|
| The standard 3D far field and grating projection functions calculate the far field profile on a hemispherical surface 1 meter from the simulation. |
|
| Calculating the field profile on a hemispherical surface creates a minor issue when we try to plot the data on a flat computer screen. The data must be 'flattened' in some way so that we can represent a curved surface on a flat computer screen. By default we plot the far field data as if we look straight down on the hemisphere. |
|
| The position vectors u1, u2 of the data returned by the far field projections use 'direction unit vector' coordinates. Each unit vector goes from -1 to 1 as shown in the figure to the right. The far field is calculated at a linearly spaced set of points as measured in the u1, u2 direction unit vector coordinates. The point u1,u2 = 0,0 corresponds to propagation at normal incidence, in the middle of the hemisphere. |
Coordinate transformations between spherical and direction cosine units are described below.
| Coordinate limits and units radius0<rm |
polar angle0≤θ≤πrad
azimuthal angle0≤ϕ≤2πrad
unit vector−1≤u≤1
|
| |
| Spherical to direction cosine ux=sin(θ)cos(ϕ) |
uy=sin(θ)sin(ϕ)
uz=cos(θ)
uz=√1−u2x−u2y
u2x+u2y+u2z=1
|
| |
| Direction cosine to spherical r=1 |
θ=acos(uz)
ϕ=atan(uyux)
| Note: farfieldspherical The farfieldspherical function can be used to interpolate far field data from ux,uy coordinates to spherical coordinates. |
| Note: Performing integrals We typically want to perform integrals in spherical coordinates such as the following power=∫∫P(θ,ϕ)R2sin(θ)dθdϕ |
where P is the Poynting vector and R is the radius. The far field projections return the electric field as a function of the variables ux and uy which are the x and y components of the unit direction vector. When changing integration variables from (q,j) to (ux,uy), it can be shown that
power=∫∫P(θ,ϕ)R2sin(θ)dθdϕ=∫∫P(ux,uy)R2duxduycos(θ)
Care must be taken to avoid numerical errors due to the cos(q) term. It's best to use the farfield3dintegrate function to evaluate these integrals.
Using spatial filtering to avoid truncating fields in far field projections
This section describes the far field spatial filtering option.
| Note: The descriptions and examples of the far field projection calculation on the following pages are primarily intended for users of FDTD. For users interested in calculating far field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist. |
For far field projections to be accurate, all radiation that will propagate to the far field must pass through the monitor being used for the projection. The far field projection functions assume that the EM fields are zero beyond the edge of the monitor. This effectively truncates the near fields at the monitor edge.
In some cases, the monitor (and simulation region) would have to be impractically large to ensure that all of the radiation passes through the monitor. One such simulation is the angular distribution of a dipole near an interface. To reproduce these results, run the simulation file and then the script.
The blue line in the following figure shows the near field data just above the interface surface. The far field projection functions will use this data to calculate the far field intensity. The important point to notice is that the near field data is non-zero near the edge of the plot (red circles). It is obvious that the these near fields do not immediately go to zero at x=8um, but that is what the far field projections assume.
This truncation will lead to high frequency ripple in the far field projection data. The far field filter option allows us to apply a spatial filter to the near field data (shown in green). When the fields go more smoothly to zero, the high frequency ripple is removed from the projection. Both projections are shown below. Change the far field filter setting in the script file to reproduce each of these results.
|
farfieldfilter(0); |
farfieldfilter(1); |
It is important to realize that the far field filter does not make the projection more accurate. The filter simply removes high frequency ripple caused by the field truncation. To increase the accuracy, the monitor (and simulation region) must be made wider, so the monitor can record more of the fields that will eventually propagate to the far field.
The far field filter has a single input parameter, which is a number between 0 and 1. By default, it is 0, which turns the filter off. This filter is applied to all far field projections. The filter parameter defines the width of the filter by the following formula: α = (a)/(a+b).
| Note: Periodic structures - far field filter The far field filter option should not be used for periodic structures. Set it to zero when using the 'assume periodic' option. |
Understanding field polarization in far field projections
This section discusses how to obtain and understand the vector field components and polarization of far-field projections.
Note: The descriptions and examples of the far-field projection calculation on this page are primarily intended for users of FDTD. For users interested in calculating far-field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist.
For some far-field projections, it is important to calculate the vectorial components of the electric field in order to determine the polarization properties. In this section, we will focus on the script commands farfieldvector3d and farfieldpolar3d. The concepts dealt with here can also be applied to the two-dimensional commands farfieldvector2d and farfieldpolar2d . We will continue to study the far-field example file described in Simple example.
Note: The descriptions and examples of the far-field projection calculation on this page are primarily intended for users of FDTD. For users interested in calculating far-field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist.
These commands return three complex components of the electromagnetic field. All the properties of the polarization of the field can be determined from the amplitudes and phases of these components.
The difference between farfieldvector3d and farfieldpolar3d is the coordinate system for defining the components of the electric field. farfieldvector3d uses a cartesian coordinate system and farfieldpolar3d uses a spherical coordinate system.
|
|
|
The script file performs a projection in both coordinate systems. The results are shown below.
| farfieldvector3d | farfieldpolar3d |
| Ex
| Er
|
| Ey
| Etheta
|
| Ez
| Ephi
|
The results in spherical coordinates are the easiest to interpret. Er is 15 orders of magnitude smaller than the other components. We expect that Er should be zero because in the far field the electric field is perpendicular to the direction of propagation. The small non-zero terms are due to numerical rounding error on double precision numbers. The polarization of the original source is P polarized. Therefore we expect that at the center of the beam, Ephi is zero and this is clearly the case on the image plots above. Due to the diffraction of the gaussian beam Ephi is non zero when the azimuthal angle is different than 15 degrees.
To see the effect of changing to S polarization, modify the property "polarization angle" to 90 degrees from 0, see the polarization angle definition.
After re-running the FDTD Simulation, run the script file again. A comparison of the two source polarizations is shown below. We can see the expected result that at the beam center there is only an Etheta component when the source is P polarized, but only an Ephi component when the source is S polarized.
| P-polarized source (polarization angle = 0 degree) | S polarized source (polarization angle = 90 degree) |
| Etheta
| Etheta
|
| Ephi
| Ephi
|
Calculating magnetic fields in far field projections
This section discusses how to obtain the H field components of far field projections.
Calculating |H|2
in the far field
To calculate ∣H∣2
in the far field, you can simply use the relation ∣H∣2=n2ϵ0μ0∣E∣2
Calculating H in the intermediate and far field
We use the farfieldexact script commands to calculate the components of the electric field at any specific point in the farfield. In this example (based on the nanoslit device shown in Wang et. al.), we can calculate the E field intensity on top of the device to analyze the focusing properties of this nanoslit structure. This is done using the following script commands:
# define desired region of x and y
x = linspace(-5e-6,5e-6,200);
y = linspace(1e-6,50e-6,500);
# do far field projection
E2 = farfieldexact2d('monitor1',x,y);
E2 = sum(abs(E2)^2,3); # calculate E2 from Ex, Ey, Ez
# plot E field intensity profile as a function of x,y
image(x*1e6,y*1e6,E2,'x (microns)','y (microns)','E field intensity');
Next, we can calculate the H fields numerically from the E fields using Maxwell's equation:
∂Ez∂y−∂Ey∂z=iωμ0Hx
∂Ex∂z−∂Ez∂x=iωμ0Hy
∂Ey∂x−∂Ex∂z=iωμ0Hz
For example, in the attached script, Hz is calculated using the following script commands:
delta = 1e-9; # used to calculate the numerical derivative;
f = getdata("T","f");
Ey1 = pinch(farfieldexact2d("T",x-delta,y),3,2);
Ey2 = pinch(farfieldexact2d("T",x+delta,y),3,2);
Ex1 = pinch(farfieldexact2d("T",x,y-delta),3,1);
Ex2 = pinch(farfieldexact2d("T",x,y+delta),3,1);
dEy_dx = (Ey2-Ey1)/(2*delta);
dEx_dy = (Ex2-Ex1)/(2*delta);
Hz = (dEy_dx - dEx_dy) / (1i * 2* pi * f * mu0);
# plot H field intensity profile as a function of x,y
image(x*1e6,y*1e6,abs(Hz)^2,"x (um)","y (um)","Hz");
Related publications
B. Wang and G. P. Wang, "Directional beaming of light from a nanoslit surrounded by metallic heterostructures," Appl. Phys. Lett. 88, 013114 (2006)
Understanding the Assume Periodic option in far field projections
This section describes the "Assume structure is periodic" options of the far field projections tab.
As discussed in previous sections, far field projections are typically used for isolated structures rather than periodic arrays. However, it is possible to use the projection functions to calculate the approximate far field distribution of finite sized periodic arrays. For infinitely periodic structures, the grating functions should be used.
| Note: The descriptions and examples of the far field projection calculation on the following pages are primarily intended for users of FDTD. For users interested in calculating far field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist. |
Non-periodic projection
By default, the far field functions assume that the near fields are zero beyond the monitor boundary (as discussed in the Far field filtering section). This occurs even when using periodic boundary conditions. Physically, this acts like an aperture placed over the structure, only allowing radiation from one period of the structure to propagate to the far field. The aperture causes strong diffraction, which tends to make these projections not very useful.
See the image below. To simulate an infinite plane wave source (red) incident on a periodic array of structures (blue), we model a single period (orange) by using periodic boundary conditions. The default far field projection gives the far field distribution (yellow) from a single period. It is as if an aperture (black) only allows light from one period to propagate to the far field.
The final result is simply the projection from a single period:
→Efar=−→E0far
The following figure shows the reflected far field |E|^2 distribution from the FDTD Blaze grating application example when using the default far field projection settings. To reproduce this figure, run blaze_grating.fsp which you can download from the Blazed Grating Application Gallery example. Using the Analysis window to plot the far field projection of the reflection monitor. Don't select the 'assume periodic' option.
For structures with a finite number of periods, there are two cases:
- The device is uniformly illuminated by the source. In other words, the source beam is much larger than the finite sized device. 'Top hat' illumination is appropriate for this case.
- The device is much larger than the source beam. In this case, the source is only illuminating a portion of the device and 'Gaussian' illumination should be used.
Top hat illumination
You will notice that as the number of periods is increased, the far field distributions will become sharper. Diffraction effects are lessened when more periods are included.
The final result is the phase correct sum of m periods, where kbloch=0 for normal incidence,
→Efar=−→E0far2∑m−2ei(kinplane−kbloch)ma
If we calculate the far field projection of the blaze structure assuming 10 periods with top hat illumination, notice how the features become sharper than the single period illumination. More periods will make the features even sharper.
Gaussian illumination
The difference between Top hat and Gaussian illumination tends to be small. Projections using Gaussian illumination are smoother because the Gaussian weighting function is smoother than the top hat function.
The final result is the phase correct sum of m periods with a Gaussian weighting, where kbloch=0 for normal incidence,
→Efar=−→E0far∞∑m=−∞ei(kinplane−kbloch)mae−(m5/2)2
The Gaussian illumination projection looks very similar to the top hat illumination.
Infinite periodic structures
As the number of periods is increased, the sharp features in the above figures become delta functions. When studying true gratings, you should consider using the Grating projections (GP) functions , rather than the far field projections.
| Note: Periodic structures - far field filter The far field filter option should not be used for periodic structures. Set it to zero when using the 'assume periodic' option. |
Integrating power in far field projections
The section explains how to integrate the far fields on the default hemispherical surface. This is typically done to calculate the amount of far field power within some range of angles.
Note: The descriptions and examples of the far field projection calculation on the following pages are primarily intended for users of FDTD. For users interested in calculating far field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist.
Power Integrals
In general, we want to integrate power over a given solid angle in the far field. There are 2 ways this can be done
- We integrate the fraction of total electric field intensity (|E|2) over the solid angle that we are interested in, and multiply by the normalized power transmission through the monitor in the near field.
far field fraction=∫cone∣E∣2∫hemisphere∣E∣2Power norm=far field fraction∗near field power transmission
- We calculate the Poynting vector in the far field and integrate the power over a given solid angle. We then normalize to the original source power. In the far field, it is easy to calculate the normal component of the Poynting vector directly from the Electric fields by using the plane wave relationships between E and H. We know the tangential components the Poynting vector are zero for a plane wave.
Poynting⊥=n√ϵ0μ0∣E∣2Power=12∫coneRe(Poynting⊥)Power norm=PowerSourcePower
Both methods will give the same result.
The script file will calculate the normalized power in a given cone around the z-axis defined by a half-angle. It will do the calculation by both methods described above. The desired half angle can be set in the first executed line of the script file
# choose the half angle over which we will integrate half_angle = 30; #in degrees
We expect that using a half angle of 30 degrees will account for 50% of the power. When we run with 30 degrees, we find have the following output
> solver_far_field2; The half angle is: 30 degrees at (theta,phi)=(0,0) The normalized transmission by Method 1 is: 45.6613 % The normalized transmission by Method 2 is: 44.5107 %
If we try several angles we find that
| Cone half angle (degrees) | Normalized transmission in cone by method 1 (%) | Normalized transmission in cone by method 2 (%) |
|---|---|---|
| 29 | 38.4055 | 37.4377 |
| 30 | 45.6613 | 44.5107 |
| 31 | 54.7707 | 53.3905 |
As expected, we find that half the power is radiated within a half-angle cone of 30 degrees around the z-axis.
| NOTE: Far field integration The function farfield3dintegrate makes integrating far field data very easy. |
| NOTE: Integration with non-default far field refractive index. Additional normalization is required when using a non-default far field refractive index. See the far field refractive index page or contact Lumerical support for additional information.
|
Creating line plots of far field projections
This section shows how to obtain line plots of 2D far field projection data.
Note: The descriptions and examples of the far field projection calculation on the following pages are primarily intended for users of FDTD. For users interested in calculating far field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist.
The farfield3d functions generally returns the far field projection as a 2D function of ux and uy. This is shown in figure 1 below. Often we want to create a line plot of this data at a specific value of theta or phi, such as the lines marked in red on the following figure. The associated script shows how to use the farfieldspherical script function to interpolate the standard farfield data to an arbitrary location defined by theta/phi.
|
polar image plot of E^2 |
Line plot of E^2 vs theta at phi=0. |
|
Line plot of E^2 vs theta at phi=15. |
Line plot of E^2 vs phi at theta=30. |
Adjusting the projection distance in far field projections
This section describes how to rescale far field projections to distances other than the default of 1m. It also describes how to use the farfieldexact functions to calculate the field distribution at arbitrary positions, including the so called intermediate field (beyond the simulation region boundary, but not yet the far field).
Note: The descriptions and examples of the far field projection calculation on the following pages are primarily intended for users of FDTD. For users interested in calculating far field projections with MODE, these descriptions are basically still correct, although some subtle differences do exist.
The script file first calculates the standard far field distribution. Rather than calculating the distribution on the entire hemisphere, we only get one line at y=0. This data is calculated with both the farfield3d and farfieldexact3d functions. As the following figure shows, both functions return the same result for the field distribution on a hemisphere with a radius of 1m.
In most cases, the standard projection location is sufficient. However, if you wish to know the field amplitude at a different distance (such as a hemisphere with a radius of 1mm), we recognize that the E and |E|^2 scale as shown in the following table. This formula is valid anywhere in the far field.
| Electric field scaling | Electric field intensity scaling | |||
| 2D | →E(R)=−→E0/√R |
| ∣→E(r)∣2=∣→E(1)∣2r |
| 3D | →E(R)=−→E0/R |
| ∣→E(r)∣2=∣→E(1)∣2r2 |
Projections to the intermediate field
If you want to calculate the field distribution quite close to the structures (the intermediate field), then the farfieldexact functions must be used. The intermediate field is the region close to the structure where the fields are not yet propagating like simple plane waves. The farfieldexact functions are able to calculate the field distribution in the far field or intermediate field (to within a few wavelengths of the simulation region).
For example, suppose we wish to know the fields on an a hemisphere with a radius of 10um. One approach (blue line in following figure) might be to simply scale the fields by a factor of (1m/10um)^2. This is NOT correct because 10um is not far enough to be considered as the far field. The EM fields are not yet propagating like simple plane waves. The correct approach is to use the farfieldexact3d function (green line).
As you can see, the two calculations do not give the same results. The blue line is not quite correct because rescaling the far field distribution is not correct in the intermediate field.
| Note: Projections to surfaces other than hemispheres Use the farfieldexact functions when you want the field distribution on a surface other than a hemisphere. The following figure shows the field intensity along a straight line for x=-20:20um at y=0um, z=10um.
|
Grating projections in FDTD overview
The near to far Grating Projections (GP) calculate the far-field profile from a periodic grating structure. The near field data is typically obtained from Lumerical's FDTD. The far-field is then calculated as a post-processing step.
A simple way to understand grating projections is to view them as a decomposition of the near field data using a set of plane waves propagating at different angles as the basis for the decomposition. The end result is that the far-field projection functions provide a straightforward and numerically efficient method for calculating the EM fields in the far-field region.
If your structure is not periodic, see the far-field projections page. If you're using DGTD or FEEM, see grating projections in DGTD.
Grating physics
The grating functions are used to calculate the direction and intensity of light reflected or transmitted through a periodic structure. For example, the grating order directions of a 2D grating can be calculated from the well known grating equation
mλ=ax(sinθm+sinθi)
where m
is an integer, ax is the periodicity of the grating in the x direction, θm is the angle of propagation for the grating order, and θi
is the angle of propagation of the incident light.
For our purposes, it's more convenient to re-write this equation in terms of the wave vector k:
(→km)x=(→kin)x+m2πax
where →km
is the wavevector of the mth grating order and →kin is the wavevector of the incident light, and the x
subscript denotes the component of the wavevector.
In 3D, these equations become:
(→kn,m)x=(→kin)x+n2πax(→kn,m)y=(→kin)y+m2πay
where kn,m
is the wavevector of the (n,m) grating order, ay is the grating period in the y direction and (n,m
) are any integers that satisfy the condition
∣→kn,m∣≤k=2π⋅index/λ0
where index
is the refractive index of the medium in which the grating order is propagating and λ0
is the vacuum wavelength.
It's important to remember that the grating order directions are defined entirely by the device period, the source wavelength and angle of incidence, and the background refractive index. In principle, the grating order directions can be calculated without running a simulation. However, in practice, the simulation must first be meshed in order for these functions to obtain necessary information such as the dimensions and period of the structure. The functions gratingn , gratingm , gratingu1 , gratingu2 , gratingangle can be used to calculate the direction of each order.
After running a simulation, the grating commands can be used to calculate the fraction of power that is scattered in each direction. The grating function uses a technique similar to a far-field projection to calculate what fraction of near field power propagates in each grating order direction. To get polarization and phase information, use gratingpolar and gratingvector.
Related publications
- Allen Taflove, Computational Electromagnetics: The Finite-Difference Time-Domain Method. Boston: Artech House, (2005).
- John B. Schneider, Understanding the Finite-Difference Time-Domain Method, Chapter 14: Near-to-Far-Field Transformation, (2010).
Using grating projections to calculate fields at an arbitrary location
In this example, we show how to decompose the transmitted or reflected fields of a periodic structure into plane waves and propagate them an arbitrary distance through homogeneous material using Grating Projection (GP). We will compare the results of a projection with a direct FDTD simulation of the same near field region, to demonstrate that this method can work well even over small distances. In practice, this allows the user to directly simulate a very small region with FDTD, then use grating projections to obtain the fields in other locations in the near, intermediate or far field.
This technique is useful in some lithography applications.
Setup
The example simulation file contains a thin film of gold on glass with a circular etch hole in it. The simulation volume represents one unit cell of a periodic structure which has a period of 1.5 microns in both the x and y directions. The source covers a wavelength range of 0.4 to 0.6 microns, and is polarized with the E field along the x-axis. We take advantage of the symmetry of the structure to reduce the simulation volume to 1/4 of the unit cell.
Analysis
The script file uses the grating script commands to decompose the field recorded at the monitor called "transmission" onto a basis state of plane waves. The plane waves can then be propagated to any distance and summed up to calculate the E field at any distance. In principle, this approach is straightforward, but in practice, it can be somewhat complicated. It is important to test your analysis against known results to ensure there are no mistakes in your analysis.
This summing is most conveniently done using the chirped z-transform, or czt, script function. There are a number of amplitude and phase correction factors in the script that must be correctly taken into account, however, the idea of the plane wave decomposition and is quite simple and involves 2 steps:
Step 1: The E field at the "transmission" monitor surface is decomposed into a sum of plane waves
→E(x,y,z0)=∑n,m→En,mexp(ikx,nx+iky,my)kx,n=2nπaxky,m=2mπay
where kx,n and ky,m are the in-plane wave vectors associated with the (n,m) diffracted order and ax, ay are the x and y periods respectively. The En,m coefficients are calculated using the gratingvector script command.
Step 2: The E field at any distance z can then be constructed simply by taking the sum
→E(x,y,z)=∑n,m→En,mexp(ikx,nx+iky,my+ikz,n,mz)kz,n,m=√k2−k2x,n−k2y,m
and this final sum is most easily accomplished uzing the czt script function.
Results
The script file propagate_periodic_test.lsf is used to test this projection method by comparing it to simulated FDTD results. After running the file propagate_periodic.fsp, running the script file will generate the following figures.
The near field, where it is recorded by the FDTD monitor
The electric field intensity in the x-z plane for the first 10 microns of propagation, as calculated by FDTD
The electric field intensity in the x-z plane for the first 10 microns of propagation, as calculated by the plane wave decomposition and projection.
A comparison of the electric field intensity vs z at (x,y)=(0,0)
A comparison of the Ex vs z at (x,y)=(0,0)
We see that there is good agreement between the projected fields and the fields as calculated by FDTD. There is some discrepancy which increases as the propagation distance increases. This difference is due to grid dispersion (the speed of light on the FDTD mesh is slightly different than free space as well as being anisotropic) and therefore we can expect that the projected result is in fact more accurate for long propagation lengths. Here, we show a plot based projection and post-processing, for up to 100 um. To obtain this plot below, set test=0 in solvers_propagate_periodic.lsf and run the script.
2703
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
