Transmission operators for the non-overlapping Schwarz method for solving Helmholtz problems in rectangular cavities

In this paper we discuss different transmission operators for the non-overlapping Schwarz method which are suited for solving the time-harmonic Helmholtz equation in cavities (i.e. closed domains which do not feature an outgoing wave condition). Such problems are heavily impacted by back-propagating waves which are often neglected when devising optimized transmission operators for the Schwarz method. This work explores new operators taking into account those back-propagating waves and compares them with well-established operators neglecting these contributions. Notably, this paper focuses on the case of rectangular cavities, as the optimal (non-local) transmission operator can be easily determined. Nonetheless, deviations from this ideal geometry are considered as well. In particular, computations of the acoustic noise in a three-dimensional model of the helium vessel of a beamline cryostat with optimized Schwarz schemes are discussed. Those computations show a reduction of 46% in the iteration count, when comparing an operator optimized for cavities with those optimized for unbounded problems.


Introduction
It is well known that large-scale time-harmonic Helmholtz problems are hard to solve because of i) the pollution effect [1] and ii) the indefiniteness of the discretized operator [2]. While the pollution effect can be alleviated by using higher order discretization schemes [3], the indefiniteness is an intrinsic property of time-harmonic wave problems, at least with standard variational formulations [4,5], and significantly limits the performance of classical iterative solvers, such as the generalized minimal residual method (GMRES) for instance. Of course, as an alternative to iterative algorithms, direct solvers can be used. However, because of the fill-in effect, whose minimization is known to be a NP-complete problem [6], the amount of memory needed to treat large-scale systems can become prohibitively high (see for instance [7]).
As an alternative to direct and (unpreconditioned) iterative methods for solving large-scale, highfrequency time-harmonic Helmholtz problems, domain decomposition (DD) algorithms, and optimized Schwarz (OS) techniques [8][9][10][11] in particular, have attracted a lot of attention during the last decades.
The key idea thereof is: i) to decompose the computational domain into (possibly overlapping) subdomains, creating thus new subproblems, ii) to solve each subproblem independently, iii) to exchange data at the interfaces between the subdomains via an appropriate transmission operator and iv) to repeat this "solve and exchange" procedure until convergence of the solution. Since all subproblems are solved independently, domain decomposition methods are parallel by nature 1 and are thus very well suited for the treatment of large-scale problems. Furthermore, as the subproblems are of reduced size, direct solvers can be used. Let also note that DD methods are rarely used as stand-alone solvers, but most of the time as a preconditioner for a Krylov subspace method such as GMRES. The design of such preconditioners for time-harmonic Helmholtz problems remains an active and challenging topic [13].
The convergence rate of an OS scheme strongly depends on its transmission operator. It is well known that the optimal operator is the Dirichlet-to-Neumann (DtN) map at the interface between two subdomains [14] (i.e. the operator relating the trace of the unknown field to its normal derivative at a given interface). However, the DtN map is rarely employed, as it is a non-local operator which leads to a numerically expensive scheme. Instead, in practice, local approximations of the DtN map are used, which lead to many different computational schemes [8][9][10][11]. To the best of our knowledge, those OS techniques share a common drawback: they ignore the impact of back-propagating waves. While this assumption is legitimate in many cases (antenna arrays [15], medical imaging reconstruction [16] or photonic waveguides [17] just to cite a few), it becomes questionable when the geometry allows resonances (even if the source does not oscillate exactly at a resonance frequency), as found for instance in lasers [18], accelerator cavities [19] or quantum electrodynamic devices [7].
The objective of this work is to develop new transmission conditions taking into account the effect of back-propagating waves, and to compare them with well-established operators neglecting these contributions. To this end, we will study a rectangular cavity, determine the DtN map and localize it by following different strategies. We will then apply the resulting new transmission operators to more general geometries. This paper is organized as follows. In sections 2 and 3 the model problem with Dirichlet boundary conditions and the associated DtN map are presented for both overlapping and non-overlapping decompositions. New transmission operators are afterwards presented in section 4 and generalized (i.e. multiple subdomains and Neumann boundary conditions) in section 5. This is followed by section 6 showing a comparison with the classical DtN map related to unbounded problems and the use of transmission operators optimized for unbounded problems as an approximation of the cavity DtN map. The new transmission operators are then validated and compared with numerical experiments involving the reference rectangular cavity in section 7. The case of geometries deviating from this last model problem is further discussed in section 8 and an engineering problem involving a model of the helium vessel of a beamline cryostat is analyzed in section 9. Finally, conclusions and final remarks are drawn in section 10.
Let us solve the following Helmholtz problem on Ω: div grad p + k 2 p = g on Ω, where p(x, y) is the unknown function, g(x, y) is a known source term and k ∈ R is the fixed wavenumber of the Helmholtz problem. Because of its boundary condition, it is obvious that (1) models a cavity problem  exhibiting both forwardand back-propagating waves. It is important to stress that for this problem to be well defined, we must assume that k 2 is not an eigenvalue of (1).
Let us now set up the following optimized non-overlapping Schwarz iterative scheme, indexed by n, to solve the cavity Helmholtz problem (1): where S ij is the transmission operator of the optimized Schwarz algorithm at Σ ij and p n i (x, y) is the solution at iteration n on domain Ω i . Let us stress that, since the subdomains do not overlap, the following holds true: n 01 = −n 10 . Once the Schwarz algorithm has converged, the solution p(x, y) of the original problem (1) is recovered by concatenating the solutions p 0 (x, y) and p 1 (x, y).
In practice, let us note that the above fixed-point scheme is usually recast into the linear system [14]: where one application of the operator A amounts to one iteration of the fixed-point method with homogeneous Dirichlet boundary conditions, where I is the identity operator, where d concatenates all n · grad p + S(p) at the interface between the subdomains and where the right hand side vector b results from the nonhomogeneous Dirichlet boundary conditions. This linear system can then be solved with a matrix-free Krylov subspace method such as GMRES.

Optimal transmission operator for the rectangular cavity problem with homogeneous Dirichlet boundary conditions
In this section, we will first determine the optimal transmission operators at Σ 01 and Σ 10 of the Schwarz scheme (2) involving the non-overlapping subdomains in Figure 1b. While this work focuses on nonoverlapping decompositions, the impact of an overlap on the optimal transmission operator is also discussed in this section.

Non-overlapping case
In order to further simplify the problem at hand, let us now assume that the source term g is zero. Obviously, by not imposing a source in our problem the solution p(x, y) is trivially p = 0 since k is not an eigenvalue. This however does not jeopardize the generality of the results derived in this section.
Let us start by taking the sine Fourier series of p n i (x, y) along the y-axis: where the functions p n i (x, s) are the Fourier coefficients and where s is the Fourier variable, whose values are restricted to the set Indeed, by restricting s to the set S, the boundary conditions are automatically satisfied. Then, by exploiting decomposition (4), the partial differential equation (2) becomes the following ordinary differential equation: , γ and ∀s ∈ S, where λ ij is the Fourier symbol of S ij . Furthermore, and for simplicity, let us define P n i (s) as: In order to find the best symbol λ ij , we need to determine the convergence radius of the iterative scheme (6). This objective can be achieved by following the strategy discussed in [10], which can be summarized as follows.

Compute
∂ p n i (x, s) ∂x at x = γ from the solutions p n i (x, s) found in the previous step.
3. The convergence radius is obtained by simplifying the transmission conditions (6b) and (6e) with the expressions found above.
By following this approach, it can be shown (see Appendix A) that the convergence radius ρ of (6) satisfies 2 : where The best transmission operator S c, dtn ij , that is the Dirichlet-to-Neumann map, is thus

Overlapping case
Let us now assume a partitioning of the domain in Figure 1a into two overlapping rectangles, as shown in Figure 2. As suggested by this figure, we define 01 (resp. 10 ) as the length of Ω 0 (resp. Ω 1 ) including the overlap and 01 (resp. 10 ) as the length of Ω 0 (resp. Ω 1 ) without the overlap. By following the same strategy as in the non-overlapping case, but taking into account that Σ 01 and Σ 10 have now different locations, the convergence radius ρ overlap of the overlapping variant of (6) reads (see Appendix B): A few conclusions can be drawn from the above expressions. Firstly, it is clear from (11) that the optimal transmission operator writes that is λ c, dtn, overlap ij (s) is equal to λ c, dtn ij (s), up to the substitution ji → ji . Secondly, it is easy to notice that when there is no overlap, it holds that ji = ji and the non-overlapping case is recovered. Thirdly and most importantly, it is obvious that an overlap does not necessarily improve the convergence radius ρ overlap unlike for the unbounded case [14], since when s 2 < k 2 (i.e. propagating modes) the term involving the sin functions in (11a) exhibits poles and oscillates with respect to the size of the overlap. Nonetheless, the term with the sinh functions in (11c) introduces a damping proportional to the overlap, as in the unbounded case, but only when s 2 > k 2 (i.e. evanescent modes).

Some local transmission operators for the cavity problem
In this section, we discuss some local transmission operators based on the Fourier symbol (9).

Zeroth-order transmission condition
A zeroth-order transmission condition (OO0 c ) can easily be constructed by approximating the symbol of the DtN map with the constant term of its Taylor series expansion around s = 0. For the considered cavity setting, we obtain: and the OO0 c transmission condition reads As such, this operator exhibits a rather poor behavior. Indeed, the denominator of the convergence radius (8) involves terms of the form which can change their sign multiple times when s 2 < k 2 , since λ c, dtn ji is nothing but a cotangent and λ c, oo0 ij is constant. Consequently, it is possible that d ij ≈ 0 for some s ∈ S, leading to a very large convergence radius. In the worst case scenario, one can also have d ij = 0 for some s ∈ S and the problem becomes illposed. Regularization procedures for preventing this behavior are further discussed in sections 4.5 and 4.6.

Truncated Mittag-Leffler expansion based transmission condition
In order to improve the performance of the above OO0 c transmission condition for cavity problem, we need a condition whose symbol is a better approximation of λ c, dtn ij . To this end, an option is to exploit the Mittag-Leffler [20] expansion of cot(z) according to its poles, leading to the following partial fraction decomposition [20]: that can be exploited to expand the symbol of the DtN map as This symbol can hence be localized by truncating the series up to the N th term, enabling us to form a N -term truncated Mittag-Leffler expansion based (ML c ) transmission condition: As increasing the number of terms in this expansion makes the transmission condition arbitrarily precise, the poor convergence rate of the previous OO0 c operator can be alleviated. This comes however at the cost of N auxiliary computations to account for the N inverse operations appearing in (19) [11]. For this reason, it is desirable to devise an approximation of λ c, dtn ij (s) with a limited number of auxiliary terms. Let us also mention that while the above operator has been developed for a one-dimensional interface, it can be used as well with a two-dimensional one, as shown in section 9 for instance.

Padé approximant based transmission condition
A Padé rational approximation exhibits usually a good convergence rate with respect to its order [M/N ], where M (resp. N ) denotes the order of the numerator (resp. denominator). One can construct it by exploiting the continued fraction expansion of the function to approach [21]. By taking the reciprocal of the continued fraction expansion of the tangent function [22], we have: where b n = 2n + 1 ∀n ≥ 0, a n = −z 2 ∀n ≥ 1.
The Padé approximant can then be determined from the following recurrence formula [21]: with That is, we have for the z cot(z) function: Starting from this recurrence formula and choosing l = m, we can devise a N -term decomposition of the form However, compared with the unbounded case where the coefficients of the Padé approximant of the DtN map are known analytically and exploited to construct the PADE u operator [11], no closed form formulae were found for the coefficientsC 0 ,Ã i andB i of (24). Nevertheless, those can be computed numerically by The numerically demanding part of this approach is the calculation of the poles of R/B N , i.e. the zeros of B N , which requires arbitrary precision arithmetic as the coefficients of the monomials appearing in B N can be vary large.
In this work, this is achieved with the MPSolve library 3 [23]. Within that framework, it takes less than 5 minutes 4 to compute theC 0 ,Ã i andB i coefficients in the very large case of N = 1024 for instance. Of course, these coefficients can be pre-computed and tabulated for various values of N and the actual transmission condition recovered with the change of variable z = ji √ k 2 − s 2 (see paragraph below). For illustration purposes, the Padé coefficients are presented in Table 1 for N ∈ [1,4].
by exploiting the change of variable z = ji √ k 2 − s 2 . The operator associated with this symbol then reads: As in the truncated Mittag-Leffler expansion case, the above operator has been developed for a onedimensional interface, but can be used as well with a two-dimensional one, as shown in section 9 for instance. Before concluding this subsection, let us mention that the truncated Mittag-Leffler expansion of z cot(z) is related to its Padé approximant in the following way: the square-root of the i th pole appearing in (24) is converging towards (i + 1)π, i.e. the (i + 1) th pole of z cot(z) appearing in (17), for a fixed value of i and as N increases. Formally, we have that lim N →∞ as show in Figure 3 for the first six poles.

Quality of the transmission operators in the one-dimensional case
As the DtN and the three transmission operators discussed above are purely real, they reduce to a simple real function in the one-dimensional case, allowing a simple graphical comparison, as suggested in Figure 4. For clarity reasons, this figure is restricted to a relatively low frequency problem with a ratio /λ w = 4.3, where λ w = 2π/k is the wavelength, and two subdomains of equal size. Nonetheless, the discussion below remains general.   To start with, it is clear that OO0 c is a quite poor approximation of the DtN, as it obviously cannot capture its oscillations and its poles. It is also apparent that changing the expansion point of the Taylor series (here s = 0, as a recall) will not improve the situation.
On the other hand, both ML c and PADE c are able to capture the oscillation and poles of the DtN, at least for a sufficiently high value of N . Those approximations exhibit however major differences, which are summarized in what follows.
• For the evanescent modes (i.e. when s/k > 1), the convergence of λ c, ml(N ) ij towards λ c, dtn ij as N increases is quite slow, at least when compared with λ c, pade(N ) ij . Let us mention that in the unbounded case, λ u, pade(N ) ij is known to perform well when s/k > 1. This aspect will be further discussed in section 6.
• For the non-evanescent modes (i.e. when s/k < 1), each pole introduced by λ To conclude this subsection, let us also note that the approximation is better as s/k tends to 1 for both λ

Regularization with a constant imaginary part
As already mentioned earlier, the zeroth-order transmission condition optimized for the cavity problem can lead to an ill-posed problem when one of the d ij (s) terms in the denominator of the convergence radius (8) equals (or is sufficiently close to) zero. This problem is however not peculiar to the OO0 c condition and affects the ML c and PADE c ones as well, requiring therefore a regularization procedure.
The most straightforward and simple strategy to regularize the OO0 c , ML c and PADE c conditions is to exploit the fact that those operators are purely real-valued. Consequently, by adding a purely imaginary part, the d ij (s) terms can be pushed away from zero and the convergence radius ρ(s) can be guaranteed to be well defined. Formally, the regularized OO0 c , ML c and PADE c operators read 5 where χ ∈ R is the regularization parameter. Numerical experiments showing the impact of this regularization approach are further discussed in sections 7 and 8.

Regularization by mixing operators optimized for cavity and unbounded problems
The above regularization is in some sense suboptimal as it acts on all s ∈ S, while regularization is required only in the s 2 < k 2 case. A more selective approach can be achieved by exploiting the PADE u operator. Indeed, assuming a sufficiently high value of N and an appropriate rotation of the branch cut [24] of PADE u , the latter exhibits the following properties [24]: it is approximately real when s 2 > k 2 and 3. it is a good approximation of the DtN map when s 2 > k 2 (see discussion in section 6).
Therefore, a regularized operator can be constructed by combining either the OO0 c , the ML c or the PADE c operator with the PADE u one in a convex way. These new operators will be further referred to as mixed operator and write 6 where ε ∈]0, 1[ denotes the regularization parameter of the mixed formulation.

Estimate for the minimum number of auxiliary unknowns
Given the ML c and PADE c transmission conditions, one question naturally arises: how many auxiliary terms should be selected? In order to answer this question, one could opt for the following criterion: the number of auxiliary terms N should at least be equal to the number of poles of λ c, dtn ij , as it seems legitimate to assume that the behavior of λ c, dtn ij is driven by its poles in the range s 2 < k 2 . Given the Mittag-Leffler expansion of λ c, dtn ij , it is clear that the poles must satisfy which implies that and thus since k, n and ij are positive by construction. Therefore, according to the pole criterion stated above, the minimum number of terms N poles min for localizing λ c, dtn and depends on the size of the subdomains and on the wavelength. As discussed further in section 7.4, this criterion seems however pessimistic, as lower values of N pole min provide already acceptable results.

Generalizations
The transmission operators detailed in the pervious section are restricted to the very particular case of a cavity with Dirichlet boundary conditions divided into two subdomains. In order to generalize this setting, we discuss in section the cases of many subdomains in a one-dimensional partitioning and of Neumann boundary conditions.

One-dimensional partitioning with more than two subdomains
Let us start with the one-dimensional partitioning of the computational domain into D subdomains, as shown in Figure 5a. In such a one-dimensional domain decomposition, the physical meaning of the ji coefficient appearing in λ c, dtn ij (s) must be clarified. In the two subdomains case, the ji coefficients represents the distance between Σ ij and the reflecting wall located in the n ij direction. This interpretation can be directly applied to the D subdomains case to define the ji coefficients, as shown in Figure 5b for the D = 3 case.
x (a) Partitioning into D subdomains.

Neumann boundary conditions
So far, we considered only the situation where the reflecting wall associated with S ij is implemented with an homogeneous Dirichlet boundary condition, i.e. a soft-wall condition. In the case of an homogeneous Neumann boundary condition, i.e. a hard-wall condition, the following DtN map is obtained: by following the same strategy as in section 3. It is worth stressing the similarities between (9) and (34). This DtN map can then be localized using the previously presented approaches and a OO0 c , ML c or PADE c transmission condition can be devised. In this regard, let us note that the Mittag-Leffler expansion of tan(z) reads [25] z tan(az and its continued fraction expansion has the following form [22]: These results are given here for the sake of completeness and will not be further discussed in this work.

Operators optimized for unbounded problems without obstacles in a cavity context
In this section we derive some estimates on the performance of operators optimized for unbounded problems without obstacles when used in a cavity problem. In particular, i) we first compare the DtN operator related to an unbounded problem without obstacles S u, dtn ij with its cavity counter part S c, dtn ij , then ii) we discuss the convergence radius of the OS scheme when using S u, dtn ij as a transmission operator for the rectangular cavity problem (1), as well as iii) the particular case of the optimized order 0 operator (OO0 u ) [8].

Dirichlet-to-Neumann operators
Let us consider the following Helmholtz problem without obstacles: where r 2 = x 2 + y 2 . In this case, it can be shown that the optimal transmission operator S u, dtn ij for solving this problem with an OS scheme is [11]: where By comparing (39) and (9), it is easy to realize that Interestingly, by exploiting the definition of the hyperbolic cotangent [22], the case s 2 > k 2 can be further simplified into which yields: lim In other words, for the case s 2 > k 2 , the symbol λ u, dtn ij (s) is converging towards λ c, dtn ij (s) as s grows. Furthermore, as the difference between both symbols is decreasing exponentially, λ u, dtn ij (s) is an excellent Regarding the case s 2 < k 2 , as the codomains of λ u, dtn ij (s) (which is purely imaginary) and λ c, dtn ij (s) (which is purely real) do not match, the expression in (40a) cannot be further simplified. For illustration purposes, the graphs of λ u, dtn ij and λ c, dtn ij are depicted in Figure 6 for different values of k ( and are respectively denoting the real and imaginary part functions).

Best convergence radius
We already know from the previous section that λ u, dtn ij (s) is a good approximation of λ c, dtn ij (s) when s 2 > k 2 , that is for evanescent modes. For this reason, local approximations of λ u, dtn ij are legitimate, yet suboptimal, candidates for approximating λ c, dtn ij .
In terms of convergence radius, as defined in (8), it is easy to show that and . For this reason, the transmission operators that are good approximations of λ u, dtn ij , such that the optimized order 2 (OO2 u ) [10] or the N -term Padé-localized (PADE u ) [11] operators, should exhibit a convergence radius close to (44). In other words, those local operators should exhibit a slow convergence for the propagating modes (s 2 < k 2 ) and a fast convergence for the evanescent ones (s 2 > k 2 ).

Particular case of the optimized order 0 operator
Before concluding this subsection, it is worth mentioning that in the case of the OO0 u operator, we have that λ u, oo0 ij = −k and therefore ρ(s) = 1. (45)

Numerical validation and comparison between the different operators
In this section we analyze the performance of the different transmission conditions developed in section 4 and compare them with the operators of section 6. To this end, we consider the rectangular cavity shown in Figure 1a and solve the time-harmonic Helmholtz equation over Ω with the following boundary conditions imposed on Γ: where K designates the number of modes used to excite the cavity. This Helmholtz problem is then decomposed into D subdomains of equal size and solved with an optimized Schwarz scheme combined with a GMRES algorithm without restart. Let us mention as well that the software implementation relies on the GmshDDM and GmshFEM [26] frameworks 7 and exploits a finite element (FE) discretization of the subproblems.
In the case of the S c, pade(N ) ij and S c, ml(N ) ij operators, the FE variational formulations involve auxiliary unknowns for the treatment of the inverse operation, as proposed in [11]. However, let us note that since the new transmission operators are not symmetric, i.e. S c, oo0/ml(N )/pade(N ) ij = S c, oo0/ml(N )/pade(N ) ji , the amount of auxiliary fields must be doubled.
As the novel transmission conditions involve operators oscillating rapidly in a wide range, the finite precision arithmetic aspects of the linear solver must be treated with care. In this regard, the orthogonalization step of the GMRES solver is critical and the modified version of the Gram-Schmidt algorithm [27,28] is required for convergence. Concerning the software implementation, the linear solvers of the PETSc [29] (GMRES) and MUMPS [30] (LU factorization) libraries are used.
Unless stated otherwise, the cavity has an aspect ratio of /h = 2 and a length-to-wavelength ratio of /λ w = 157.085 2π ≈ 25.001. This configuration allows 25 non-evanescent modes and is excited with the K = 50 first modes. The geometry is discretized with a mesh consisting of 8 triangular elements per wavelength and the subproblems are discretized with an FE method of order 4. The stopping criterion of the GMRES solver is set to a relative tolerance decrease of r i / r 0 = 10 −6 , where r i is the residual vector at iteration i and r 0 is the residual vector of the first guess, which was chosen equal to zero.

Two subdomains case
Let us start by studying the performance of the transmission conditions developed in section 4 when applied to a rectangular cavity divided into two subdomains. The convergence history of the GMRES solver is shown in Figure 7. From these data, it is clear that all transmission conditions converge. In addition, it appears clearly that S  As the closed form solution p CF of this canonical problem is known, i.e.
with k y (m) = m π h and k x (m) = k 2 − k 2 y (m), let us determine the accuracy of the above simulations by computing the relative L 2 error E between the FE solution p FE and p CF : where the above integrals are evaluated using a quadrature rule with twice the amount of integration points than the number used for the FE computations. The L 2 errors associated with the different OS schemes of Figure 7 are gathered in Table 2 together with the error associated with the MUMPS direct solver. These data show clearly that all transmission conditions and the direct solver lead to errors of the same order of magnitude.  Let us now focus on the wall-clock time required to complete the above computations 8 , as reported in Table 3. In order to analyze these values, let us stress that they heavily depend on the actual software implementation of the FE and OS tools. Nonetheless, some general remarks can be drawn.
1. The OO0 u , OO2 u and OO0 c operators are the computationally cheapest to apply, as they do not involve auxiliary unknowns.
2. The ML c and PADE c operators are computationally more expensive than PADE u , since they require two sets of auxiliary unknowns as they are not symmetric, i.e.  3. The regularization procedure involving PADE u further increases the computational cost, as additional auxiliary unknowns are introduced.
In the current software implementation, the subproblems are solved with the direct solver MUMPS and the resulting LU factorization is reused in the subsequent OS iterations. Therefore, the first iteration is more time consuming than the other ones and the data in Table 3 are thus split in different subquantities. By defining T i as the share of the total wall clock time T tot dedicated to the i th iteration, i.e. T tot = I i=1 T i with I the total number of GMRES iterations required for convergence, we report in Table 3 the following values: T tot , T 1 , T 2,I , where T i,j is the mean value of the sequence [T i , . . . , T I ], and I.
From the data gathered in Table 3, it is clear that S c, pade(32) ij leads to both the minimal amount of iterations and the fastest computation with respect to the wall clock time. Nonetheless, it is evident that the cost of T 1 and T 2,I is higher for S c, pade(32) ij and S c, ml(32) ij than for the other transmission operators, in accordance with the remarks drawn above. The novel operators will therefore lead to the fastest computations only when the reduction of the iteration count is sufficiently high, when compared with the transmission conditions optimized for unbounded problems (see related discussion in sections 9).  Table 3: Wall-clock times -rectangular cavity with D = 2 subdomains (T i is the share of the total wall-clock time, in seconds, taken by the i th iteration, I is the total number of GMRES iterations and T i,j is the mean value of the sequence [T i , . . . , T j ]).

Spectrum of the iteration operators
In this section, let us briefly discuss the spectra of the discretized iteration operator F = I − A, which we will refer to as F from now on 9 . As the explicit construction of F is computationally heavy, only the case of D = 2 subdomains is considered here. It is also important to stress that F is not a normal matrix 10 , as it can be directly seen from the formal expression of F in the case of two subdomains (see [14, section 2.3.1] for instance). As a consequence, the behavior of the GMRES cannot be predicted from the spectrum Eig(F ) of the system matrix F [31]. Nonetheless, eigenvalues that are clustered near 1 are a good indicator that those modes are well captured by a given transmission operator.
The spectra in Figure 8 are associated with the novel transmission operators optimized for the rectangular cavity. Let us stress that the eigenvalues are all real in the present example. However, the spectrum of F can be complex, even when its coefficients are all real, as it is non-normal. With the S  When using a transmission operator optimized for unbounded problem in a cavity setting, the spectra in Figure 9 are obtained. As expected from the analysis of section 6, the S u, oo0 ij operator leads to eigenvalues 9 We use a calligraphic (resp. non-calligraphic) typeface for continuous (resp. discrete) operators. 10 That is F F H = F H F , where F H is the conjugate transpose of F . lying on the unit circle C centered around 1. In addition, the behavior of the PADE u operators has been well anticipated by our previous discussions as well. It is indeed clear that the eigenvalues of F fall into two categories: 1. those associated with evanescent modes, which form a cluster similar to the PADE c one and 2. those associated with non-evanescent modes, that lie on C.
To conclude this subsection, let us also stress that despite a clustering significantly better for the PADE u operators than for OO2 u , both operators exhibit rather similar GMRES convergence curves, as shown in Figure 7. This shows the difficulty in predicting the convergence of GMRES from Eig(F ), since F is nonnormal.

Regularized and mixed transmission conditions
As discussed in sections 4.5 and 4.6, a regularization term or a mixed transmission condition can be used to prevent the convergence radius from becoming very large (or, in the worst case, ill-defined). The performance of those transmission conditions is shown in Figures 10a and 10b for S c, oo0, r(χ) ij and S m(ε, M ), oo0 ij respectively with D = 8 subdomains. It is clear from the figures that regularizing the OO0 c operator by adding an imaginary part comes at the cost of an increased number of iterations. On the other hand, mixing OO0 c with PADE u leads to an improvement in terms of iteration count. Nonetheless, in all cases, the regularized operators do not lead to any improvement with respect to the PADE u conditions. Let us also note that the regularized operators tend to the original OO0 c as χ becomes small (resp. ε becomes large).
A similar numerical experiment can also be carried out for the regularization of the ML c condition. From the GMRES convergence histories depicted in Figure 11, it can be noticed that while the performance of ML c remains better than PADE u , the regularization increases the number of iterations required to reach convergence. Nonetheless, this increase declines as the regularization becomes lighter, i.e. when χ (resp. ε) becomes small (resp. large).
An analogous numerical experiment is performed once more for the regularization of the PADE c condition. This last scenario shows a behavior similar to ML c , as it can be directly observed in Figure 12.      Figure 13 for a numerical experiment involving D = 8 subdomains of equal size. It is clear from the data that every operator converges with values of N as low as N = 1. Nonetheless, for the novel operators to outperform the convergence of the PADE u operator, a minimum value of N data min = 8 is required in this example. It is also evident in this example that S c, pade(N ) ij always performs better than S c, ml(N ) ij for a fixed value of N . Last but not least, it is evident from Figure 13 that N has only a mild effect on the performance of  It is also interesting to compare the value of N data min found experimentally and the value predicted with the pole criterion discussed in section 4.7. In this numerical experiment, the largest ij is max ij = 7 8 and thus N pole min = 44 according to the pole criterion (33). Consequently, N pole min is a pessimistic estimation in this case, since N data min < N pole min .

Increase in the number of rectangular subdomains
In this section we focus on the impact of the number of subdomains D onto the GMRES iteration count, as shown in Figure 14. From this plot, is clear that S    of GMRES iterations with the optimal slope of 2D, at least for the considered range of D. Let us mention also that the unbounded transmission operators exhibit a significantly larger slope, motivating thus the use of transmission conditions specifically devised for cavity problems. Furthermore let us note that S u, pade(64) ij and S u, oo2 ij present the same slope, as they are both excellent localization of the unbounded DtN map.
In addition, while the S c, ml(64) ij operator shows a suboptimal scaling with respect to D, it outperforms S u, pade(64) ij and S u, oo2 ij . Last but not least, it is evident from Figure 14 that the scaling behavior of S c, oo0 ij is the worst of all.

Impact of the length-to-wavelength ratio
The dependence of the GMRES iteration count on the wavenumber is a major performance indicator of an OS scheme, and a numerical experiment analyzing this is therefore carried out. In this investigation the values of /λ w are chosen close to an integer value, which corresponds to a cavity driven at a frequency close to one of its resonance frequencies. The computational domain is partitioned into D = 8 subdomains and the cavity is excited with double the amount of non-evanescent modes (this number depending on /λ w ). It is clear from the data shown in Figure 15 that the iteration count increases rapidly with /λ w with the unbounded transmission operators. This growth in the iteration count is, nonetheless, not as fast as for S c, oo0 ij . On the other hand, the situation is significantly improved with S c, ml(N ) ij and the iteration count becomes almost independent from /λ w with S c, pade(N ) ij , at least in the considered /λ w range. Let us also note that similar results are obtained when /λ w is selected away from a resonance.

Discussion
From the data gathered in the previous numerical examples, it is obvious that the S c, oo0 ij operator leads to a poor scaling when D > 2. For this reason, this operator and its regularized variants are not further considered in this work.
It is also clear that for rectangular cavities the It is also worth noticing that the unbounded operator S u, pade(N ) ij and S u, oo2 ij lead to very similar convergence profiles, which is not surprising as they both are excellent approximations of the unbounded DtN map. For this reason, only the S u, pade(N ) ij operator will be further considered. Another interesting result concerns the impact of a change of the wavenumber on the convergence profile of the above transmission operators. For the considered range of length-to-wavelength ratios /λ w , we observed that the S c, pade(N ) ij behaves quasi-independently from /λ w , while the other operators exhibit a (quasi) linear increase of the iteration count with /λ w . Nonetheless, the slope of this increase remains small for S c, ml(N ) ij when compared with S u, oo0 ij , S u, oo2 ij and S u, pade(N ) ij . Let us finally note that similar experiments were carried out on three-dimensional rectangular parallelepipedic cavities with a square cross-section of h × h and a length of . Results analogous to the above two-dimensional test cases were obtained.

Sensitivity to geometrical parameters
The previous section considered only one geometry (a rectangular cavity and its three-dimensional equivalent), which refers to the specific configuration used to devise the transmission operators discusses in section 4. Nonetheless, relevant simulations involve computational domains that differ from this original setting. We therefore discuss in this section two geometries deviating from the canonical one.

Trapezoidal geometry -modified Gram-Schmidt variant
For this first numerical experiment, let us consider an isosceles trapezoid whose base is characterized with a parameter δ, as the one depicted in Figure 16. This computational domain is partitioned into D = 16 subdomains and, as with the previous rectangular cavity, the aspect ratio /h equals 2 and the length-to-wavelength ratio /λ w is approximately 25.001. The number of GMRES iterations associated with various operators are shown in Figure 17 for different values of δ. It is clear from the depicted results that the novel rational operators outperform, in all considered cases, the S u, pade(64) ij one in terms of iteration count, even when no regularization is applied. Additionally, the S c, pade(64) ij is systematically the best choice and its regularization always leads to a slight increase in the iteration count. In the case of the ML c variant, regularization slightly improves the convergence for the highest value of δ. Furthermore, there is an evident trend in the displayed data showing an increase in the iteration count as δ increases. However, as this increase also affects S u, pade(64) ij , it is hard to predict if and when this operator will overtake the novel ones, at least in terms of iteration count.

Trapezoidal geometry -classical Gram-Schmidt variant
In the introduction of section 7, we mentioned that using the modified version of the Gram-Schmidt (GS) algorithm is critical for the convergence, as the novel transmission conditions involve operators oscillating rapidly in a wide range. In order to assess the importance of this choice, we carry out once more the previous numerical experiments using the classical version of the GS orthogonalization procedure [27].
As it would be impractical to show the data for all possible combinations, only the case with PADE u and PADE c (with and without regularization) at δ = 32 is shown in Figure 18, as it contains the different behavior we want to highlight. First of all, it is clear from the displayed data that the modified GS does not impact the behavior of the PADE u operator. In addition, it is also evident that it prevents stagnation and improves the overall convergence with PADE c . Regarding the other combinations that are not shown in Figure 18, we observed that: i) a relative residual r i / r 0 smaller than 10 −6 is not always reached when using the classical Gram-Schmidt procedure, ii) this relative residual is always reached with the modified GS, iii) the convergence history of the PADE u operator is identical for both GS variants and iv) the modified GS is never worse than the classical GS.

Rectangular cavity involving obstacles -impact of the number of obstacles
Let us now study again a rectangular cavity, but let us introduce O circular obstacles in the domain. As in the previous case, we have /h = 2 and /λ w ≈ 25.001. Furthermore, we assume that the obstacles exhibit a hard-wall behavior and that each subdomain includes at most one obstacle located in its center. A sketch of a possible configuration is shown in Figure 19, with O = 3 obstacles and D = 3 subdomains.

Rectangular cavity involving obstacles -impact of the size of the obstacles
In the previous subsection we assumed a rather large radius for the obstacles compared with the wavelength and focused only on the number of obstacles. Let us now reverse the study and let us determine the impact of the obstacle size on the performance of the transmission operators. To this end, let us again consider a rectangular cavity with D = 17 subdomains and O = 7 obstacles of radius R and let us gather in Figure  in all cases and ii) that regularization leads to an increased iteration count.

Discussion
The above numerical experiments clearly show that, while the performance of the novel S c, ml(N ) ij and S c, pade(N ) ij operators deteriorates as the geometry of the problem deviates from the reference rectangular case, this deterioration does not exhibit sudden jumps as the geometry is modified. In addition, the sensitivity of the novel operators with respect to numerical errors stemming from finite precision arithmetic is also clearly shown by our data, motivating the use of the modified Gram-Schmidt algorithm in the orthogonalization step of GMRES. When compared to the PADE u operator, PADE c and ML c clearly outperform the former for small deviations. On the other hand, in all considered cases, the reduction of the iteration count brought by PADE c and ML c can become modest (with respect to PADE u ) for large deviations. Nonetheless, we did not encounter situations where PADE u leads to clear gain, at least in terms of iteration count and within the scope of the numerical experiment considered in this work. This suggests a rather robust behavior of PADE c and ML c .

Engineering test case: acoustic noise in a three-dimensional model of the helium vessel of a beamline cryostat
Going back to the experiment with a rectangular cavity filled with obstacles, it has been shown that the PADE c operator shows a very fast convergence when the cavity exhibits only one obstacle, as shown in Figure 20. Additionally, the same behavior is observed for other configurations exhibiting a unique circular obstacle. A such characteristic can be very well suited when, for instance, simulating the acoustic noise in the helium vessel of a beamline cryostat, which basically consists of a rectangular cavity interrupted by a circular obstacle (e.g. see the cryostat discussed in [32]). Such acoustic noise analyses can become critical, for instance, when designing a cryogenic current comparator (CCC) with a large bore [33]. A CCC is one of the most sensitive instrument for measuring very low electric currents with high accuracy and can be used e.g. in particle accelerators for the non-destructive monitoring of slowly extracted charged particle beams (current intensities below 1µA) [33].
Within this context, let us compare the behavior of the different transmission operators on a helium vessel model consisting in a rectangular parallelepipedic cavity interrupted with a cylindrical obstacle and partitioned into D = 8 subdomains. The model is excited with its first spatial mode, the length-to-wavelength ratio is chosen as /λ w 12.5004 and the cross-section is h × h with h = /2. This numerical experiment leads to the convergence history displayed in Figure 22 for each transmission operator. It is clear from these data that the S c, pade(24) ij leads to the quickest convergence and requires only 50 iterations to converge. Compared with the 93 iterations required by S u, pade(24) ij , a reduction of the iteration count of approximately 46% is achieved. Figure 22 shows as well the impact of a light regularization of the PADE c and ML c operators. As expected, no performance gain is obtained, as the unregularized counterparts already converge well. For illustration purposes, the computed field map is available in Figure 23.
The above analysis would not be complete without considering the wall clock time 11 . As already dis- 11 The simulations of section 9 were carried out on the NIC5 cluster hosted at the University of Liège, Belgium. cussed in section 7, the novel transmission operators optimized for cavity problems are associated with an increased computational cost, when compared with their unbounded counterparts. As a result, in the context of the considered acoustic study, the OO0 u and OO2 u operators lead to the fastest computations (i.e. approximately 3.5 and 3 hours respectively), despite their rather slow convergence. In comparison with OO2 u , the rational operators see their wall clock time increased by a factor of i) 2 for PADE u , ii) 2.7 for ML c and iii) 2.5 for PADE c . Before concluding this section, it is important to stress that the current implementation does not treat the auxiliary unknowns in the most efficient way. For instance, it assembles all auxiliary unknowns into a single large system. However, the extra unknowns stemming from the asymmetric nature of S c, ml(N )/pade(N ) ij could be pulled out into smaller auxiliary systems, which could improve the overall computational time. Let us also mention that the OS iterative scheme can also be altered, such that all auxiliary unknowns are decoupled from the subproblems, as proposed in [34]. Such improvements are left for a future work.

Conclusion and final remarks
In this work we presented new transmission operators for the Schwarz method optimized for timeharmonic Helmholtz problems in a rectangular cavity. Those operators rely on different localizations of the Dirichlet-to-Neumann map of this reference geometry. Three different strategies were considered for devising local approximations of the Fourier symbol of the DtN, namely: i) a zeroth-order Taylor approximation, ii) a truncated Mittag-Leffler partial fraction expansion and iii) a Padé approximant. As these operators do not necessarily lead to a well defined scheme, two regularization procedures where proposed to correct this problem. Let us however mention that, when combined with the modified variant of the Gram-Schmidt procedure in the orthogonalization step of GMRES (without restart), the proposed operators led to convergent iterative schemes in all the numerical experiments considered in this work, even in the unregularized case.
The new operators optimized for cavity problems were also compared with operators optimized for unbounded geometries, yet applied to cavities. In the case of the reference rectangular geometry, the gain in the iteration count brought by the novel rational operators is clear when compared with the unbounded ones. On the other hand, when considering geometries deviating from the canonical one, this iteration count gain decreases. Nonetheless, even if the gain can be modest for large deviations, the new rational operators always led to a lower iteration count (than their unbounded counterparts) in the considered numerical experiments. Of course, this does not guarantee that the novel rational operators will always perform better, but this suggests that they are sufficiently robust with respect to geometrical changes. The computational cost of the different transmission operators was also briefly discussed, as well as some possible improvements for reducing it.
To conclude this paper, let us draw some final remarks regarding the selection of an operator. While a detailed analysis is out of the scope of this work, some general trends can be extracted from this study. First of all, when considering the wall-clock time, as the computational cost of the novel rational operators is high, a simple OO0 u operator might be a good default choice. Nonetheless, as already mentioned, alterations of the Schwarz scheme could improve this aspect. On the other hand, when focusing only on the iteration count, the PADE c operator (with a small regularization) can be recommended as a first guess. Indeed, since evanescent modes in cavities and unbounded problems converge quickly to each other, both Padé based approximations should behave similarly for these modes if N is sufficiently large. Let us note that this is coherent with the spectra discussed in section 7.2. Thus, the difference between PADE c and PADE u lies mainly in the non-evanescent part, for which the PADE u is clearly a poor approximation. Therefore, in the worst case scenario, it is legitimate to expect that both PADE c and PADE u will behave similarly, at least with a GMRES without restart, as shown in our numerical experiments. Let us note that this expectation still needs to be investigated more formally, which is left for a future work.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Koch, Mr. Achim Wagner, Mr. Dragos Munteanu, Mr. Christian Schmitt, Dr. Wolfgang F.O. Müller and Dr. David Colignon for the administrative and technical support. Finally, the authors would like to thank the anonymous Reviewers, whose comments improved significantly the quality of this work.
Appendix A. Calculation of the convergence radius -non-overlapping case In this appendix, let us briefly discuss the main steps of the calculation leading to the Fourier symbol of the Dirichlet-to-Neumann map shown in equation (9). To this end, we start with case s 2 = k 2 , and write the solutions of (6a) and (6d), together with the boundary conditions (6c) and (6f) and the definition (7). Formally, we obtain: This solution can then be derived with respect to x and the result evaluated at the interface between the subdomains, i.e. at x = γ(t) = t − /2 with t ∈ [0, 1] as discussed in section 2. We thus have that The convergence radius ρ(s) is then obtained by simplifying the transmission conditions (6b) and (6e) with the above expressions. In particular, we can write revealing thus the convergence radius ρ(s = k). Again, equation (9) is recovered, for the case s 2 = k 2 , by identifying the above terms and by exploiting the definition of 01 and 10 .