Modal analysis of the ultrahigh finesse Haroche QED cavity

Let us now discuss the geometry considered for the simulations. In order to reduce the number of unknowns, only one quarter of a single mirror is modeled, as shown in figure 3. Furthermore, the air surrounding the quarter mirror is represented by a rectangular parallelepiped, as shown in figure 4(a). As discussed in section 2, a PML is added to model the surrounding space, as depicted in figure 4(b). As a first guess, the PML thickness is taken as 1λ (where $\lambda \simeq 5.9\,\mathrm{mm}$ is the wavelength at $51.1\,\mathrm{GHz}$), and its distance from the mirror is also set to 1λ. More details on the PML can be found in sections 2 and 3.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final geometry is depicted in figure 4(b), and was generated by the Open CASCADE⁸ engine. This library was driven using the python⁹ interface provided by Gmsh¹⁰ [33]. Finally, the computational domain is meshed using the Gmsh mesh engine. Since our objective is to use high-order FE discretizations, curved mesh elements are mandatory [34], not only to achieve a precise representation of the curved mirror, but also to keep the number of mesh elements to an acceptable value.

Our objective is to find the eigenvalues of the photon cavity. In other words, we need to find every possible value of the angular frequency ω, such that

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{\bf{curl}}({{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{\boldsymbol{e}})-{\omega }^{2}\displaystyle \frac{{{\boldsymbol{\varepsilon }}}_{r}}{{c}_{0}^{2}}{\boldsymbol{e}}={\bf{0}}\quad {\rm{on}}\,{\rm{\Omega }},\\ {\rm{boundary}}\,{\rm{conditions}}\,{\rm{of}}\,{\rm{section}}\,{\rm{3.3}},\end{array}\right.\end{eqnarray*}$$

holds, where ${c}_{0}$ is the speed of light in vacuum, ${\boldsymbol{e}}$ the electric field, and where Ω is the computational domain (i.e. the quarter mirror surrounded by the PML). To clarify the presentation, the treatment of the boundary conditions is discussed in the next subsection.

From a FE point of view, this eigenvalue problem writes [35]:

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{\displaystyle \int }_{{\rm{\Omega }}}({{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{\boldsymbol{e}})\cdot ({\bf{curl}}\,{{\boldsymbol{e}}}_{i})\,{\rm{d}}{\rm{\Omega }}-{\omega }^{2}{\displaystyle \int }_{{\rm{\Omega }}}\left(\displaystyle \frac{{{\boldsymbol{\varepsilon }}}_{r}}{{c}_{0}^{2}}\,{\boldsymbol{e}}\right)\cdot {{\boldsymbol{e}}}_{i}\,{\rm{d}}{\rm{\Omega }}=0\quad \forall {{\boldsymbol{e}}}_{i}\in {V}_{0}({\bf{curl}},\,{\rm{\Omega }}),\\ {\rm{boundary}}\,{\rm{conditions}}\,{\rm{of}}\,{\rm{section}}\,{\rm{3.3}},\end{array}\right.\end{eqnarray*}$$

where ${V}_{0}({\bf{curl}},\,{\rm{\Omega }})$ is a finite-dimensional subspace (of size N) of ${H}_{0}({\bf{curl}},\,{\rm{\Omega }})$, the space of square integrable functions with a square integrable curl (these function are also commonly referred to as edge elements). The electrical field ${\boldsymbol{e}}$ is then approximated by:

$$\begin{eqnarray*}&&{\boldsymbol{e}}=\displaystyle \sum _{j=0}^{N-1}{a}_{j}\,{{\boldsymbol{e}}}_{j}\quad {{\boldsymbol{e}}}_{j}\in {V}_{0}({\bf{curl}},\,{\rm{\Omega }}),\end{eqnarray*}$$

where the coefficients a_j are unknown. Thus, the (generalized) eigenvalue problem writes:

$$\begin{eqnarray}&&{\rm{find}}\,{\rm{every}}\,{\omega }_{k}^{2},\,{\rm{such}}\,{\rm{that}}\quad {{\bf{Ba}}}_{k}={\omega }_{k}^{2}{{\bf{Ca}}}_{k}.\end{eqnarray} \tag{ 1 }$$

By identifying the terms in equation (1), we have:

• the vector ${{\bf{a}}}_{k}={[{a}_{k,0},...,{a}_{k,N-1}]}^{T}$ containing the degrees of freedom for ${\boldsymbol{e}};$
• the matrix ${\bf{B}}(i,j)={\int }_{{\rm{\Omega }}}({{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{{\boldsymbol{e}}}_{j})\cdot ({\bf{curl}}\,{{\boldsymbol{e}}}_{i})\,{\rm{d}}{\rm{\Omega }};$
• the matrix ${\bf{C}}(i,j)={\int }_{{\rm{\Omega }}}\left(\tfrac{{{\boldsymbol{\varepsilon }}}_{r}}{{c}_{0}^{2}}\,{{\boldsymbol{e}}}_{j}\right)\cdot {{\boldsymbol{e}}}_{i}\,{\rm{d}}{\rm{\Omega }}$.

In this paper, the eigenvalue problem (1) is solved by the Krylov–Schur algorithm [36, 37] from the SLEPc library [26, 27]. The relative tolerance of the eigenvalue solver is set to 10⁻¹⁵.

As motivated above, only one quarter of a single mirror is modeled. To simulate the actual configuration, appropriate boundary conditions have to be imposed on the geometry of figure 4(b). In what follows, ${\bf{n}}$ is the unit vector outwardly oriented normal to Ω.

• (i)  
  A zero tangential magnetic field on the Oxy-plane:

  $$\begin{eqnarray*}&&{\bf{n}}\times {{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{\boldsymbol{e}}={\bf{0}}\quad {\rm{on}}\,{Oxy}.\end{eqnarray*}$$
• (ii)  
  Depending on whether polarization 1 or 2 is considered.

  • Polarization 1: a zero tangential magnetic field on the Oxz-plane, and a zero tangential electric field on the Oyz-plane:

    $$\begin{eqnarray*}\left\{\begin{array}{ll}{\bf{n}}\times {{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{\boldsymbol{e}}={\bf{0}} & {\rm{on}}\,{Oxz},\\ {\bf{n}}\times {\boldsymbol{e}}\times {\bf{n}}={\bf{0}} & {\rm{on}}\,{Oyz}.\end{array}\right.\end{eqnarray*}$$
  • Polarization 2: a zero tangential magnetic field on the Oyz-plane, and a zero tangential electric field on the Oxz-plane:

    $$\begin{eqnarray*}\left\{\begin{array}{ll}{\bf{n}}\times {\boldsymbol{e}}\times {\bf{n}}={\bf{0}} & {\rm{on}}\,{Oxz},\\ {\bf{n}}\times {{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{\boldsymbol{e}}={\bf{0}} & {\rm{on}}\,{Oyz}.\end{array}\right.\end{eqnarray*}$$
• (iii)  
  A zero tangential electric field on the outer boundary of the PML:

  $$\begin{eqnarray*}&&{\bf{n}}\times {\boldsymbol{e}}\times {\bf{n}}={\bf{0}}\quad {\rm{on}}\,\partial {{\rm{\Omega }}}_{{\rm{PML}}}.\end{eqnarray*}$$

Moreover, since the mirror is made of superconducting material, we assume that it has a zero resistivity:

$$\begin{eqnarray*}&&{\bf{n}}\times {\boldsymbol{e}}\times {\bf{n}}={\bf{0}}\quad {\rm{on}}\,\partial {{\rm{\Omega }}}_{{\rm{Mirror}}}.\end{eqnarray*}$$

For simplicity, we also consider the frame of the mirror to have a zero resistivity:

$$\begin{eqnarray*}&&{\bf{n}}\times {\boldsymbol{e}}\times {\bf{n}}={\bf{0}}\quad {\rm{on}}\,\partial {{\rm{\Omega }}}_{{\rm{Frame}}}.\end{eqnarray*}$$

Let us now specify the PML used for the considered simulations. Since we chose a parallelepiped to surround the mirror, it is a natural choice to use a Cartesian PML [7, 8]. A PML is usually described by its damping functions σ_x(x), σ_y(y) and σ_z(z). In this paper, we chose the following hyperbolic profiles:

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{\sigma }_{x}(x)=\tfrac{{c}_{0}}{{X}_{{\rm{\max }}}+{\rm{\Delta }}X-| x| }-\tfrac{{c}_{0}}{{\rm{\Delta }}X},\\ {\sigma }_{y}(y)=\tfrac{{c}_{0}}{{Y}_{{\rm{\max }}}+{\rm{\Delta }}Y-| y| }-\tfrac{{c}_{0}}{{\rm{\Delta }}Y},\\ {\sigma }_{z}(z)=\tfrac{{c}_{0}}{{Z}_{{\rm{\max }}}+{\rm{\Delta }}Z-| z| }-\tfrac{{c}_{0}}{{\rm{\Delta }}Z},\end{array}\right.\end{eqnarray*}$$

where ΔX (or ΔY, or ΔZ) is the thickness of the PML in the x (or y, or z) direction, and where X_max (or Y_max, or Z_max) is the distance between the center of the mirror and the PML in the x (or y, or z) direction. These profiles are inspired from [38], where a $\tfrac{{c}_{0}}{{\rm{\Delta }}\alpha }$ (with α∈{X, Y, Z}) term is added to remove the jump between the truncated domain and the PML.

From these damping profiles, it is then possible to compute an equivalent relative electrical permittivity and relative magnetic permeability in the PML domain:

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{{\boldsymbol{\varepsilon }}}_{r}^{{\rm{PML}}}=\mathrm{diag}\left[\tfrac{{\gamma }_{y}{\gamma }_{z}}{{\gamma }_{x}},\tfrac{{\gamma }_{x}{\gamma }_{z}}{{\gamma }_{y}},\tfrac{{\gamma }_{x}{\gamma }_{y}}{{\gamma }_{z}}\right],\\ {{\boldsymbol{\mu }}}_{r}^{{\rm{PML}}}=\mathrm{diag}\left[\tfrac{{\gamma }_{y}{\gamma }_{z}}{{\gamma }_{x}},\tfrac{{\gamma }_{x}{\gamma }_{z}}{{\gamma }_{y}},\tfrac{{\gamma }_{x}{\gamma }_{y}}{{\gamma }_{z}}\right],\end{array}\right.\end{eqnarray*}$$

where ${\gamma }_{i}=1-{\rm{i}}\tfrac{{\sigma }_{i}}{\omega }$ for all i ∈ {x, y, z}.

In our cavity problem, it is not required to extract the entire spectrum of (1): indeed, only the TEM₉₀₀ mode is of interest. Therefore, a spectral transform is used. The idea is to modify the eigenvalue problem, so that the eigenvalues around a given target are easier to compute [39]. In practice, the shift-and-invert transform [39, 40] is used. That is, instead of solving (1), the following equivalent (from the eigenvalue point of view) problem is solved:

$$\begin{eqnarray*}&&{({\bf{B}}-\sigma {\bf{C}})}^{-1}{\bf{Ca}}={\theta }_{k}{\bf{a}},\end{eqnarray*}$$

where σ is a well chosen shift and where ${\theta }_{k}={({\omega }_{k}^{2}-\sigma )}^{-1}$. Since the resonance frequency of the cavity is known to be ${f}_{{\rm{Cavity}}}=51.1\,\mathrm{GHz}$ [5], the shift σ is taken as:

$$\begin{eqnarray*}&&\sigma ={(2\pi \times {f}_{{\rm{Cavity}}})}^{2}=1.0309\times {10}^{23}\,{{\rm{r}}{\rm{a}}{\rm{d}}}^{{\rm{2}}}\,{{\rm{s}}}^{-{\rm{2}}}.\end{eqnarray*}$$

It is worth noticing that the spectral transform involves an LU decomposition of $({\bf{B}}-\sigma {\bf{C}})$, which is handled by the MUMPS library.

Based on the computational setup presented previously, a first set of simulations is run in order to assess the numerical convergence of the computed solutions. In practice, five tetrahedral meshes are generated, with 3 to 7 mesh elements per targeted wavelength (i.e. $\lambda \simeq 5.9\,{\rm{mm}}$, as mentioned in section 3). In each case, the mesh exhibits a 3rd-order curvature. From the FE point of view, three discretizations are considered: order 3, order 4 and order 5. Each simulation is carried out with 120 MPI¹¹ processes. Depending on the problem size, these processes are distributed between 5 and 120 computing nodes. On each node, $64\mathrm{GB}$ of memory is available.

For each simulation, three eigenvalues are computed. In each case, only one eigenvalue exhibits a very large damping time and is selected as the resonant one. Figure 5 depicts the numerical convergence of the computed resonance frequencies and damping times for polarization 1.

Zoom In Zoom Out Reset image size

For the resonance frequency, we may conclude that convergence is achieved: the resonance frequency for polarization 1 is found at $51.0847\,\mathrm{GHz}$. However, for the damping time, convergence is more delicate to assess. Clearly, the 3rd-order simulations did not yet converge. The 4th- and 5th-order simulations seem to converge toward a damping time around $5.6\,{\rm{s}}$. More accurate discretizations would of course be welcome, but, due to memory limitations discussed in section 5, these data are not accessible. Nevertheless, a few indicators suggest that the computed damping time will remain in this range. First, the 4th- and 5th-order curves suddenly saturate around $5.6\,{\rm{s}}$. Secondly, at a mesh density of 5 mesh elements per wavelength, the damping time varies by only a factor of 5% between the two highest-order solutions. Finally, a variation of only 1% is recorded between the best 4th-order estimate and the 5th-order one.

Compared to the damping time of $130\,\mathrm{ms}$ determined experimentally, the computed one is 43 times higher. Four major idealizations made in our model could explain this discrepancy. First, our geometry is ideal: the mirrors are exactly toroidal and perfectly aligned. Secondly, our model consists of only the two mirrors floating alone in the void, whereas the mirrors are surrounded by sophisticated equipment in the actual experiment [5], that could introduce additional losses. Thirdly, the residual resistivity of the superconducting material has not been taken into account. Indeed, since we are not operating at absolute zero temperature, some carriers are not in a superconducting state, leading thus to ohmic losses. Furthermore, it can be shown that these losses are increasing with the square of the frequency [41]. And finally, the roughness of the actual mirror is not modeled. As we will see in the following section, this has a dramatic effect on the value of the computed damping time. Let us stress that these issues can be alleviated by technical improvements, while the open nature of the cavity, and the associated computed losses, is an intrinsic property of the device.

To conclude, and for illustration purposes, figure 6 depicts the eigenmode associated to the cavity resonance. For clarity, the vectorial electric field is shown only on the Oxz- and Oxy-planes. It can be directly noticed that this mode exhibits a cylindrical symmetry around the z axis and 4.5 antinodes. Because of the boundary condition on the Oxy-plane, the total number of antinodes is doubled, confirming thus that polarization 1 is a TEM₉₀₀ mode.

Zoom In Zoom Out Reset image size

In the previous convergence test, only polarization 1 was considered. Let us now focus on polarization 2. Based on the previous simulations, only the 4th- and 5th-order FE solutions are considered. Figure 7 depicts the computed frequencies and damping times for polarization 1 and 2. In addition, the frequency difference between the two polarizations is provided in table 2.

Zoom In Zoom Out Reset image size

Table 2.  Frequency difference between polarization 1 and 2: 3rd-order mesh elements, 4th- and 5th-order finite element discretizations.

	Frequency difference (MHz)
Mesh density $({\lambda }^{-1})$	FE order 4	FE order 5
3	1.2631	1.2638
4	1.2659	1.2638
5	1.2636	1.2639
6	1.2639	—
7	1.2639	—

As for polarization 1, we can directly notice that the cavity resonance frequency has numerically converged to $51.0834\,\mathrm{GHz}$. Moreover, we have reached a frequency difference of $1.26\,\mathrm{MHz}$, which matches the difference measured experimentally [5, 32].

For the damping time, we get a result similar to polarization 1: the computed damping seems to converge toward a value around $8.1\,{\rm{s}}$. Again, the computed damping time is significantly larger than the measured one. It is worth noticing that the damping time is larger for polarization 2 than for polarization 1. This phenomenon is quite easy to explain, since polarization 2 is associated to the largest radius.

So far, the geometry has been meshed with 3rd-order mesh elements. Let us now analyze the impact of the mesh curvature on the simulations. To do so, only polarization 1 is considered for simplicity. Moreover, we limit ourselves to FE discretizations of order 4. As done previously, the geometry is meshed with a density of tetrahedra varying between 3 and 7 mesh elements per wavelength. However, this time, the geometrical order of the elements ranges between 1 and 4. Figure 8 shows the computed frequencies and damping times.

Zoom In Zoom Out Reset image size

By analyzing the data of figure 8, we can directly notice that the damping time does not converge with 1st-order mesh elements. On the other hand, there is no significant difference between simulations using high-order meshes. Obviously, straight mesh elements fail to compute the cavity damping time. Since the shape of the mirrors is curved, approximating it with 1st-order elements introduces some kind of numerical roughness on its surface. This roughness does not significantly impact the resonance frequency. However, it can lead to unwanted scattering, destroying the stability of the wave. This problem has been noticed in the previous attempt to simulate the photon cavity [6].

A last question remains: can the lack of curvature of the 1st-order mesh elements be compensated by a better mesh refinement of the mirrors? To answer this question, let us setup the following simulations. A mesh with 5 tetrahedra per wavelength is used everywhere in the computational domain, except on the surface of the mirror, where refinements of 5, 10 and 20 elements per wavelength are used. The geometrical order of the mesh ranges between 1 and 4, and a 4th-order FE discretization is used for the computations. Let us note that only polarization 1 is considered. The computed damping times are reported in figure 9.

Zoom In Zoom Out Reset image size

Based on the results of figure 9, we can directly notice that increasing the mesh density on the surface of the mirror does not help to achieve convergence for the damping time, at least when straight mesh elements are used. Thus, we can conclude that the damping time is highly sensitive to the numerical roughness introduced by the mesh. Therefore, using curved mesh elements is mandatory to simulate the cavity in order to avoid this virtual roughness.

Let us now study the sensitivity of the computed solution with respect to the parameters of the PML: its distance from the mirror and its thickness. For this analysis, the following setup is used:

• a tetrahedral mesh of order 4, with a density of 5 mesh elements per wavelength;
• a FE discretization of order 4;
• a PML thickness ranging from 0.25λ to 2λ;
• two distances between the mirror and the PML, 1λ and 2λ.

The computed damping times are shown in figure 10. The mean resonance frequency and the maximum deviation from this mean are given in table 3.

Zoom In Zoom Out Reset image size

Table 3.  Resonance frequency: polarization 1, 5 mesh elements per wavelength, 4th-order mesh elements and 4th-order finite element discretization.

Distance PML-to-mirror $(\lambda )$	Mean frequency (GHz)	Maximum deviation (GHz)
1	51.08468	2.5 × 10−6
2	51.08468	2.3 × 10−6

Let us start by analyzing the results presented in table 3. Based on the available data, we can conclude that the resonance frequency does not depend on the PML parameters. On the contrary, figure 10 indicates that the damping time is very sensitive to these parameters. Therefore, we cannot conclude that the damping time converged with respect to the PML-to-mirror distance and the PML thickness, since the behavior is too oscillatory. However, the variations of the damping time seem to remain bounded, which indicates that, even though a larger PML (in term of PML-to-mirror distance and thickness) is preferred, the cavity damping lies in the $5.4\,{\rm{s}}$ range.

To conclude this subsection, it is also worth mentioning that two additional numerical experiments were carried out. First, in addition to the boundary condition given in section 3.3 for the PML boundary, we also imposed:

$$\begin{eqnarray*}&&{\bf{n}}\times {{\boldsymbol{\mu }}}_{r}^{-1}{\bf{curl}}\,{\boldsymbol{e}}={\bf{0}}\quad {\rm{on}}\,\partial {{\rm{\Omega }}}_{{\rm{PML}}}.\end{eqnarray*}$$

Secondly, in addition to the PML equivalent materials given in section 3.4, we also tested the following constant PML parameter:

$$\begin{eqnarray*}&&{\gamma }_{x}(x)={\gamma }_{y}(y)={\gamma }_{z}(z)=1-{\rm{i}}.\end{eqnarray*}$$

In both cases, the same conclusions were drawn. Furthermore, our findings are also supported by additional PML convergence results presented in section 6 for a two-dimensional setup. There, the numerical convergence of the PML is clearly observed, and it is shown that a PML-to-mirror distance of 1λ, as well as a PML thickness of 1λ, already give numerically sharp results.

To conclude the analysis of our three-dimensional high-order computations, let us give some order of magnitude about the wall clock time taken by the simulations. The smallest simulation (3 mesh elements per wavelength with a 3rd-order FE discretization) counted 517556 unknowns and took less than 1 minute. It was run using 5 computing nodes. On the other hand, the largest simulation (7 mesh elements per wavelength with a 4th-order FE discretization) counted 11443760 unknowns and took less than 11 hours. It was run using 120 computing nodes.

Let us look at the two simulations hereafter:

• FE discretization of order 5;
• mesh curvature of order 5;
• two mesh densities, 5 and 6 mesh elements per wavelength;
• polarization 1.

These two simulations led to, respectively, 8019588 and 13363722 unknowns. The test cases were launched on 120 computing nodes, with 64 GB of memory available on each node. The smallest simulation ran successfully, whereas the largest simulation failed because of a memory deficiency during the LU factorization stage. This latter case was run again with 240 nodes but failed for the same reason.

Let us now look at the peak virtual memory allocated by each process for the problem with 8019588 unknowns, as shown in figure 11. Statistics are available in table 4. From these data, we can see that the memory usage is not uniformly distributed across the computing processes. On average, we use 17 GB with a quite low standard deviation. However, we have two spikes above 25.5 GB (i.e. 50% above the mean). It is those spikes that limit the memory scaling of our simulations.

Zoom In Zoom Out Reset image size

Let us take another example:

• FE discretization of order 4;
• tetrahedral mesh of order 4, with a density of 7 mesh elements per wavelength;
• polarization 1.

We ran this setup (11443760 unknowns) successfully on 120 and 240 computing nodes. The peak virtual memory distribution is available in figure 12, and statistics are available in table 5.

Zoom In Zoom Out Reset image size

From the above results, we have once again some memory spikes far above the mean value. Let us focus only on the two largest spikes: 46 GB (240 nodes) and 61 GB (120 nodes). We can directly notice that by increasing by $100 \%$ the total available memory, the largest spike is decreased by only 25%, thus limiting the scaling performance.

Since our two-dimensional model exhibits fewer unknowns, let us first use larger mesh densities compared to the three-dimensional case, in order to assess the convergence of the damping time. From our previous analysis, we found that 2nd-order mesh elements were sufficient to represent the geometry precisely. Thus only this kind of element is used. The mesh density varies between 5 and 40 mesh elements per wavelength. For the FE discretization, 3rd-, 4th- and 5th-order FE bases are considered. The PML is located at one targeted wavelength of the structure, and its thickness is set to one wavelength. The computed resonance frequencies and damping times are reported in figure 14.

Zoom In Zoom Out Reset image size

This figure indicates that the computed resonance frequencies are monotonously converging to $49.9804\,\mathrm{MHz}$ for increasing mesh densities and FE discretization orders. Regarding the computed damping times, the calculated results are also monotonously converging to $7.2\,{\rm{s}}$. This value is below the $8.1\,{\rm{s}}$ found for the three-dimensional polarization 2 case. However, let us recall that this two-dimensional case is modeling a similar, but not equivalent, situation. Furthermore, it is worth mentioning that already with 5 mesh elements per wavelength, with a 4th- or a 5th-order FE discretization, numerically valid damping times are computed. This observation agrees with the analysis performed for the three-dimensional case, and is an additional indicator that the results presented in the sections 4.1 and 4.2 are numerically valid. Finally, it is worth noticing that for both two- and three-dimensional situations, the computed damping time is well above the $130\,\mathrm{ms}$ found experimentally. This comforts us in our assumption that the discrepancy between the simulations and the experimental data has to be attributed to our idealizations, rather than in numerical convergence problems.

Let us now analyze the sensitivity of the results with respect to a variation of the PML thickness and position. In this numerical experiment, the following parameters are used: a mesh density of 5 mesh elements per wavelength, 4th-order mesh elements and a 5th-order FE discretization. Based on the previous analysis, we know that these parameters are sufficiently accurate. In this numerical experiment, both the PML thickness and position (with respect to the mirror) are taken equal and in the $\{1\lambda ,2\lambda ,4\lambda ,8\lambda ,16\lambda ,32\lambda \}$ set. This parameter sweep leads to the results presented in figure 15 and table 6.

Zoom In Zoom Out Reset image size

Based on these data, we can conclude that both the computed resonance frequencies and the damping times are not significantly impacted by a change of the PML parameter. It is however interesting to notice that, as expected, the damping time is more sensitive to a PML variation as compared to the resonance frequency. Again, the numerical stability of the results obtained for this two-dimensional case, with a large PML parameter sweep, comforts us that our three-dimensional simulations are correct. In particular, we can directly notice that numerically valid results can be obtained with a PML-to-mirror distance, as well as a PML thickness, of 1λ.

For this numerical experiment, we compared a 1st-order FE method with our high-order simulations. Figure 16 depicts the computed frequencies and damping times for a mesh density ranging from 5 to 320 straight mesh elements per wavelength. As expected (based on our three-dimensional analysis), the virtual roughness introduced by the straight mesh significantly impacts the numerical results. Indeed, we need at least 160 mesh elements per wavelength to reach a value close to the expected value of $7.2\,{\rm{s}}$.

Zoom In Zoom Out Reset image size

It is worth noticing that the largest simulation setup leads to 20363036 unknowns, corresponding to a linear system with twice the size of the largest three-dimensional case. Nevertheless, this simulation does not suffer the same memory limitations. This phenomenon is easily explained by the better non-zero structure of the matrices of the two-dimensional computations, resulting in a lower fill-in.

In this section, let us compare the computation time taken by different mesh densities and FE discretization orders. The results are summarized in table 7. As can be noticed directly, the computation time can be reduced by increasing the FE discretization order and decreasing the mesh density at the same time (while keeping the accuracy constant). This leads to a reduction of the number of unknowns, and thus a decrease of the computation time.

Table 7.  Computation time for different parameters with the two-dimensional model.

Order
FE (–)	Mesh (–)	Mesh density $({\lambda }^{-1})$	Unknowns (–)	Comp. time (s)	Damping (s)
1	1	320	20363036	2682	7.43
2	2	$40$	$799443$	$81$	7.17
3	2	$20$	$371588$	$46$	7.16
4	2	$10$	$148445$	$23$	7.16
5	2	$5$	$56532$	$11$	7.15

For this last part of our two-dimensional simulations, we computed the resonance frequency and the associated damping time for the other 'trapping' modes of the cavity. Based on our previous convergence analysis, we chose to use a 5th-order FE discretization. The mesh exhibits a density of 10 mesh elements per targeted wavelength and a 2nd-order curvature. It is worth stressing that, since the resonance frequencies are different for each mode, the meshes are also different. This is also the case for the actual length of the PML and its distance to the mirror. Let us also note that because of the symmetry conditions, only the odd modes were computed.

By analyzing the results displayed in figure 17, we can directly notice the linear dependence of the resonance frequency on the mode number. This can be also observed for the damping time in a logarithmic scale, with the notable exception of the transition between the 7th and the 9th modes, where the slope suddenly decreases. It is worth noticing that, after a deeper investigation, this phenomenon does not seem to come from a numerical problem. Finally, let us recall that in our simulations, the surface resistivity of the superconducting mirrors was neglected. However, this parameter is expected to limit the damping time increase with the mode number, since it grows with the frequency [41].

Zoom In Zoom Out Reset image size