This page describes how to calculate the quality factor (Q) of resonance peaks in a resonant cavity.

There are two classes of cavities for Q factor calculations, low Q cavities and high Q cavities. The 2D example file includes both the standard (high) Q analysis object, and the low Q analysis object. The 3D example only includes the standard analysis object because the quality factor of this cavity is too high for the low Q object. Both simulation files contain a cavity that supports two modes, and a Q analysis group to compute the Q factor. The high and low Q analysis groups can also be inserted into your simulation from the object library in the resonators section.

## Low Q cavities

A cavity is called a low Q cavity when the electromagnetic fields decay completely from the simulation in a time that can be simulated reasonably by FDTD. In this case, the quality factor can be determined from the Fourier transform of the field by finding the resonance frequencies of the signal and measuring the full width half maximum (FWHM) of the resonant peaks. We can then use Q = f_{R}/f where f_{R} is the resonant frequency and f is the FWHM.

### Example:

The low Q analysis group contained in the quality_factor_2D.fsp simulation uses this method to find the Q factor of resonance frequencies. Run the simulation file. Then choose to edit the low Q analysis object and press the RUN ANALYSIS button. The analysis script output will contain the location of the resonance frequencies and their corresponding Q factors.

Resonance 1: frequency = 230.727THz, or 1299.34 nm Q = 156.883 +/- 0.82616 Resonance 2: frequency = 205.814THz, or 1456.62 nm Q = 77.498 +/- 0.226738

The analysis script also creates two plots. The plot shown below to the left contains one of the field components (Hz). You can see that the fields have decayed by the end of the simulation time. The second plot shows the location and relative amplitude of the resonance peaks.

Note that the initial transients of the source are neglected by setting the "start time" for the time monitors to 200fs. The "start time" for the time monitors is the time at which the monitors begin recording data. This setting can be changed in the user properties for the analysis group. Also, note that in the analysis group, it is possible to use one time monitor or an array of time monitors for the Q factor calculation. The problem with using one time monitor is that if the one monitor is placed at or near a null of the cavity mode, then due to the fact that the field intensity is very low, the Q factor can have a large uncertainty (if it is even possible to obtain a meaningful result).

The quality_factor_3D.fsp simulation file contains a 3D version of the low Q analysis object. This object can also be added from the object library.

## High Q cavities

A cavity is considered to be a high Q cavity when the electromagnetic fields cannot completely decay from the simulation in a time that can be simulated reasonably by FDTD. In this case, we cannot determine Q from the frequency spectrum because the FWHM of each resonance in the spectrum is limited by the time of simulation, T_{sim}, by FWHM ~ 1/T_{sim}. Instead, the quality factor should be determined by the slope of the envelope of the decaying signal using the formula

$$Q=\frac{-2 \pi_f{R} \log _{10}(e)}{2 m}$$

where f_{R} is the resonant frequency of the mode, and m is the slope of the decay in SI units.

### Derivation of Q factor formula:

The quality factor (Q) is defined as

$$Q=\frac{\omega_{r}}{F W H M}$$

where w_{r} is the resonant frequency ( ω_{r}=2π f_{R}) and FWHM is the full width half max of the resonance intensity spectrum. The time domain signal of the resonance is described by

$$E(t)=e^{-r\left(\alpha-i \omega_{r}\right)} u(t)$$

where α is the decay constant. The Fourier transform of E(t) is easy to calculate

$$|E(\omega)|^{2}=\frac{1}{\alpha^{2}+\left(\omega-\omega_{r}\right)^{2}}$$

The maximum value of |E( ω)|^2 is clearly 1/α^2, at ω= ω_{r}. With a little more work, we can determine that the half max frequencies occurs at ω= ω_{r} + α and ω= ω_{r} - α. Therefore, FWHM = 2α. Substituting this value into the original Q formula and solving for α gives

$$\alpha=\frac{\omega_{r}}{2 Q}$$

Now that we know how to relate α to Q, we must determine how the slope of the time signal decay is related to Q. We must take the log of the time signal to make the envelope a linear function.

$$\log _{10}(|E(t)|)=\frac{-\omega_{r} t}{2 Q} \log _{10}(e)=m t$$

where m is the slope of the log of the time signal envelope. Solving for Q, we get

$$Q=\frac{-\omega_{r} \log _{10}(e)}{2 m}$$

### Example:

Calculation of the Q factor for high Q cavities is complicated because

- separating the decay of the envelope from the underlying sinusoidal signal is difficult since the fields are typically real-valued
- if there are multiple resonant modes, they will interfere with each other in the time domain, making it hard to estimate the decay rate.

By opening the edit dialog box for the Q factor analysis object located in quality_factor_3D.fsp, you can see that the analysis object solves these problems by

- accurately calculating the envelope of the time-domain field signal
- isolating each resonance peak in the frequency domain using a Gaussian filter, and then taking the inverse Fourier transform to calculate the time decay separately for each peak. The slope of the time decay is then used to calculate the Q factor and obtain an error estimate.

In addition, note that:

- the Q analysis object has setup variables that allow you to choose how many time monitors to use to calculate the Q factor. It is often a good idea to add a few point monitors at different locations to reduce the chances that a monitor is placed at a node in the mode profile of a cavity mode yielding a weak signal.
- in the analysis tab, there is a parameter that can be set to choose how many resonant peaks to look for
- all the field components that are available are used to calculate the Q factor
- it is possible to change other parameters, such as the Gaussian filter width and resolution in the frequency domain. These parameters are set in the analysis script.
- in the script, only the part of the time signal lying in 40-60% of the time signal collected is used for the slope calculation. These percentages can easily be changed. However, setting the upper limit to anything greater than 90% can lead to errors due to the fact that Fourier transforms, and inverse transforms were used when the Gaussian filter was used to isolate the peak. The Fourier transforms introduce errors to the end of the time signal due to the fact that discrete Fourier transforms assume periodicity of the signal.

Next, run the simulation. When the simulation is complete, choose to edit the analysis object and press RUN ANALYSIS button. The analysis script output will contain the location of the resonance frequencies and their corresponding Q factors.

Resonance 1: frequency = 178.786THz, or 1676.82 nm Q = 306.279 +/- 1.41318 Resonance 2: frequency = 227.307THz, or 1318.89 nm Q = 274.874 +/- 4.50921

The analysis object also produces the following plots (by modifying make_plots=0 to make_plots=1 before running the analysis).

The time decay of the field components and their envelopes. Note that not all the field components have decayed by the end of the simulation time. The image shown at the top of this topic shows a close-up of the data for the field drawn in red (Ey).

The spectrum and the Gaussian filters used to isolate the resonant peaks.

The spectrum of resonances. Each resonant peak appears in a different color.

The time decay of the sum of squared field components on a logarithmic scale. The slope of these lines is used in the Q factor calculation.

The quality_factor_2D.fsp simulation file contains a 2D version of the Q analysis object.

# Advanced notes

### Resonator response to a given source excitation

In a linear system, the response of a resonator is an intrinsic property of the device, independent of the type the source excitation used. In the frequency domain, the response of the system can be related to the source by

$$t(\omega)=h(\omega) \cdot i(\omega)$$

where t is the response (for example, the electric field at a point in space), i is the source signal (for example, the dipole moment) and h is the transfer function that relates both quantities for each angular frequency w. The quantity h is an intrinsic property of the resonator and is the impulse response of the system, ie the response when i(w) = 1, or i(t) is a delta function.

When calculating quality factors, we are interested in the case where h contains a single pole response

$$h(\omega)=\frac{r}{i\left(\omega-\omega_{r}\right)-\alpha}+b(\omega)$$

where b(w) are additional contributions to h, typically representing much broader features in the frequency domain that decay very rapidly in the time domain. We know α is related to Q by

$$\alpha=\frac{\omega_{r}}{2 Q}$$

**Low Q cavities**

The ideal method for calculating Q is to calculate h(w) and extract Q. This is what we do when Q is small enough that we can run the simulation until the fields have decayed. We then calculate h(w) by

$$h(\omega)=t(\omega) / i(\omega)$$

In FDTD, with the cw normalization turned on, we are calculating h(w) and therefore obtain the intrinsic quality factor, independent of the type of the source excitation pulse or spectrum. Repeating Q measurements with different lengths of source pulses, or even with extremely complicated source pulses will not change the function h(w) measured except for some numerical error which can become large if i(w) is very small near the resonant frequency.

**High Q cavities**

When working with high Q cavities, we want to find a method that is more efficient to calculate Q because it is impractical to run FDTD simulations long enough to calculate h(w) correctly. As a result, we recognize that **if and only if** the source pulse is a delta function, we have

$$t(\omega)=h(\omega)$$

$$t(t)=r \exp \left(-\alpha t+i \omega_{r} t\right) u(t)$$

We can therefore extract the value of Q by measuring α from the exponential decay in the time domain. It is critical to note that we only have exponential decay of the time signal if the source pulse is a delta function. If we use any other function to excite the the resonator, we cannot use this method to extract the quality factor as we cannot disentangle the contribution of the source pulse and the resonator itself without running the simulation until the fields have completely decayed. In practice, we cannot use a delta function for the source, but we can make the approximation that the source is a delta function as long as the function i(w) is nearly constant over the bandwidth of the resonator. As a result, it is not relevant to try and normalize away the source spectrum as we do in the low Q case, because we must assume that the source is approximately a delta function to perform this calculation. Using other source spectra will result in time domain fields that do not necessarily decay exponentially.

By using a near delta function as the source input, we can assume that the time domain fields will decay exponentially, and we can try to extract the quality factor from that decay. In reality, h(w) is more complex than a single pole response, due to the contributions of b(w), and we must wait for the other contributions of h to disappear in the time before we can measure the decay rate. Fortunately, the contributions of b(t) generally decay much faster. For this reason, we must remove some of the early time signal to measure the exponential decay accurately. The amount of time that must be removed is somewhat problem dependent.

The majority of the calculation in the high Q cavity calculation is associated with extracting the exponential decay from a time signal of the form

$$\exp (-\alpha t) \cos \left(\omega_{r} t\right)$$

and eliminating the early part of the time signal where h(t) has other terms that make it difficult to extract the decay. The more advanced parts of the calculation involve extracting several decay rates from several single pole resonances at different resonant frequencies by doing some filtering in the frequency domain before returning to the time domain to calculate the decay.