{'fraction_non_alphanumeric': 0.04081128638903467, 'fraction_numerical': 0.016268946994721083, 'mean_word_length': 4.256206392990427, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Elastic interior transmission eigenvalues and their computation via the method of fundamental solutions

A stabilized version of the fundamental solution method to catch ill-conditioning effects is investigated with focus on the computation of complex-valued elastic interior transmission eigenvalues in two dimensions for homogeneous and isotropic media without voids. Its algorithm can be implemented very shortly and adopts to many similar PDE-based eigenproblems as long as the underlying fundamental solution function can be easily generated. We develop a corroborative approximation analysis which also implicates new basic results for transmission eigenfunctions and present some numerical examples which together prove successful feasibility of our eigenvalue recovery approach.


Introduction
Interior transmission eigenvalues (ITEs) arise primarily in the study of inverse scattering theory, cf. [2,3,8]. Their major role comes along with conventional reconstruction algorithms for which incident wave frequencies of that specific order need to be excluded to fully justify the feasibility of qualitative methods for the target recovery process. While the acoustic and electromagnetic cases have been thoroughly investigated in the last years, see [9,19,12,7,11,14,13,16], there are comparably only a few papers covering the numerical study of elastic ITEs so far, e.g. [10,18,20,21], especially when focusing on the overall complex-valued spectrum. It is the purpose of this paper to introduce the approved algorithm from [13] and [14] also for their computation. As for the latter, we confine to ideally isotropic, homogeneous and planar-symmetric scatterers to resemble the 2D case, but allow for less regular shapes within our theoretical foundation now which are assumed to be known in advance. ‡ Author to whom any correspondence should be addressed.
Numerical challenges in computing ITEs origin from the fact that the underlying interior transmission eigenvalue problem (ITP) is non-self-adjoint, non-elliptic and nonlinear in the eigenvalue parameter. Its straightforward discretization would therefore result in non-Hermitian matrices whose pseudospectra are generally harder to capture, especially when the matrix size becomes large. In this context, the advantage of Ritztype methods such as the method of fundamental solution (MFS) whose prototype goes back to Kupradze in the 1960th for approximating solutions of well-posed boundary value problems, see [15], is that surprisingly accurate results can already be obtained for many regular domains while only a relatively small number of trial functions is used. Due to this observation our MFS-based eigenvalue algorithm similar to [4] will be very efficient at remarkably low numerical costs and is still easy-to-implement as being meshand integration-free unlike the usual competitive methods, cf. [10].
In the course of substantiating our approach mathematically, we prove some novel findings concerning norm relations of ITP eigenfunction pairs which surprisingly differ for purely real and complex-valued ITEs with non-vanishing imaginary part. Since the field of eigenfunction properties is quite undiscovered unlike their eigenvalues themselves in the ITP context, see [5], we want to point out that most of our results presented here should be adaptable for the acoustic and electromagnetic case, too.
The remainder of this paper is structured as follows: In section 2 we will setup the modeling framework of the ITP. Section 3 will then consider the eigenvalue approximation problem from an abstract and thus more general perspective that also fits other boundary control techniques such as the boundary element method. Definitions and prerequisite results that are specific for our MFS ansatz will be given in section 4 and are followed by numerical examples in section 5. A final conclusion will summarize the merits of our proposed method at the end.

Setup of the elastic interior transmission eigenvalue problem
We model time harmonic vibrations of a bounded and elastically penetrable solid D with homogeneous mass density ρ = const > 0 which is assumed to differ from its normalized background material through ρ = 1. Mathematically, we think of D as a domain D ⊂ R 2 of class C 1,1 whose regional displacement, when absent or present in comparison with its surrounding, will be modeled by vector fields u, v : D → C 2 , respectively. For both the generalized but isotropic version of Hooke's law shall relate their gradients in a linear way to the internal stress σ : C 2×2 → C 2×2 via σ(z) = 2µǫ(z) + λtr(ǫ(z))I .
Here, I ∈ R 2×2 is the identity matrix and ǫ(z) = (∇z + (∇z) T )/2 is the symmetric part of the gradient for any displacement field z measuring its corresponding strain. Further, µ, λ are Lamé parameters which coincide for simplicity for both the scatterer and the background and which are constrained to µ > 0 as well as 2µ + λ > 0 to guarantee strong ellipticity of the governing Navier system, see [17] div(σ(z)) + ̺ω 2 z = µ∆z + (λ + µ)∇(div z) + ̺ω 2 z =: ∆ * z + ̺ω 2 z = 0 . (1) These equations describe from a physical perspective the spatial part of the elastic wave propagation, where ̺ represents any mass density under consideration. With this notation at hand we may now formulate the elastic interior transmission problem (ITP): where u, v ∈ L 2 (D, C 2 ) and the scattered part (u − v) is supposed to be an element of the more regular Sobolev space : ϕ = 0 and σ(ϕ)ν = 0 on ∂D} with outer normal ν along the boundary ∂D. Frequencies ω ∈ C\{0} which admit nontrivial solutions (u, v) to the above ITP will be called interior transmission eigenvalues (ITEs). Apparently, they refer to those pairs of harmonic waves in D whose behavior along the scattering boundary coincides and which thus lock the possibility of detecting the scatterer on the basis of close-by data.
Note that the co-normal derivative appearing in the Neumann boundary condition of (2) is connected naturally to the highest order term of their PDEs via integration by parts, also known as Betti's first formula which holds by duality for any ψ ∈ L 2 (D, C 2 ) and ϕ ∈ H 2 (D, C 2 ). Here, the colon symbol denotes the Frobenius inner product given by A : B = tr(AB H ), whereas the single dot refers to the scalar-product-like operation a·b = a H b for a, b ∈ C 2 . Having thus set the mathematical framework of our eigenproblem, in the sequel we try to recover ITEs ω by solving (2) approximately in the sense that we allow for small deviations within the boundary data that are assumed to vanish in some limiting procedure.

Approximation analysis of elastic interior transmission eigenvalues via boundary control
We define our relaxed space of trial functions for approximating solutions of (2) subject to boundary optimization by where for any ω ∈ C\{0} By definition, any pair in H(ω) fulfills the required PDE conditions from (2) automatically and it is the choice of ω that determines in how far their boundary data are compatible in the sense of the ITP. Note that the above restrictions on ω in H are due to the following theorem on the overall locations of complex ITEs and the fact that both ω, −ω and −ω are each ITEs if and only if some ω from the first quadrant in the complex plane is. We drop the proof of the former since it would exactly follow the lines from its acoustic analogon in [6] with the obvious operator adaptions.
Then it holds that where λ(D) is the smallest interior Dirichlet eigenvalue for the Navier problem (1) with ̺ = ρ. In particular, if ω lies in the first quadrant of the complex plane, then 0 ≤ arg(ω) < π/4.
We aim to extract those ω which allow for approximate ITP eigenfunctions in H with a relatively small ratio of boundary misfit to interior norm. The next theorem states that if these residual quotients can be made arbitrarily small while the corresponding ω accumulate, their limit is indeed an ITE. Its proof slightly refines the technique from [14,Theorem 2] to encompass also complex-valued eigenvalues now.
Proof. By rescaling and redefining (u m , v m ), if necessary, we can assume without loss of generality that C = 1 and aim to apply weak compactness in order to construct a solution candidate which will indeed meet all the required ITE criteria. By assumption (ii) we know (modulo the extraction of subsequences which we will not relabel in m) that u m ⇀ u and v m ⇀ v in L 2 (D, C 2 ) which implies that (u, v) ∈ L 2 (D, C 2 )×L 2 (D, C 2 ) fulfills the interior conditions of the ITP, or equivalently the PDE system (1), according to for any bump function ϕ ∈ C ∞ c (D, C 2 ) and with a similar calculus for v. In order to prove that (u − v) also has the correct ITP boundary data, it suffices to prove that the differences (u m − v m ) are bounded in H 2 (D, C 2 ) in combination with assumption (iii) because w → w |∂D and w → σ(w) |∂D ν are continuous as operators H 2 (D, C 2 ) → H 3 2 (∂D, C 2 ) and H 2 (D, C 2 ) → H 1 2 (∂D, C 2 ), respectively, due to our regularity assumptions on ∂D. Noting that (u m −v m ) are solutions to the inhomogeneous Navier system elliptic estimates like in [17] tell us that which therefore gives the desired uniform bound with respect to m. It remains to show that (u, v) is non-trivial. For this we recall that the embedding > 0 so we will assume contrarily that (u m − v m ) → 0 in L 2 (D, C 2 ). Then, on the one hand, the bounded sequence a m := D u m · v m dx may be singled out to converge to some a ∈ C from which we know by (ii) with C = 1 that i.e. a = 1/2. Since arg(ω) < 1/4 and ρ = 1, we can even conclude |Re a(ω 2 − ρω 2 )| > 0. On the other hand, (3) with exchanged roles of its test functions ϕ, ψ ∈ C ∞ (D) yields the analogon of Green's second identity Duality shows that the separated boundary contributions of the approximate eigen- are even uniformly bounded in m thanks to (i)-(iii). Therefore we may further compute which is a contradiction.
In order to establish for fixed m some more qualitative relation between the boundary misfit of some sufficiently approximate eigenfunction pair and its absolute eigenvalue deviation from the actual one, we derive an interconnecting estimate in the following which can also be seen as practical guiding principle for numerical calculations. Similar to the acoustic case, see [14,Lemma 5], but now improved for even less regular shapes, its rigorous validity relies on a critical integral expression which must not vanish.
Lemma 3. Let (u, v) be an ITP eigenfunction pair with ITE ω and assume that ( u, v) ∈ H( ω) for some aritrary frequency w. If the integral constraint is fulfilled, then there exists a constant C > 0 determined only by the boundary data of u (or equivalently v) such that Proof. Applying Betti's formula twice, using the identical ITP boundary data for u and v, yields and after rearranging we obtain Taking absolute values gives .
Duality shows that C < ∞ since u solves the Navier equation.
We state a direct consequence with respect to the L 2 (D, C 2 )-norm of eigenfunctions for frequencies ω ∈ C\R.
→ 0 and the left hand side of (7) converges because when setting u := u m and v := v m . Since |ω 2 − ω 2 | > 0 by assumption, the right hand side of (8) forces any uniform bound on ε to vanish in the limit m → ∞.
The situation is different for ω ∈ R and especially when it is the ITE of smallest magnitude. In this case, with the additional assumption that ρ > 1 or 0 < ρ < 1 is sufficiently large or small to be made precise next, respectively, we can show that = 0 which also guarantees a uniform bound ε > 0 for approximate eigenfunction pairs ( u, v) in the vicinity of (u, v) (cf. Corollary 6 later). To avoid misleading confusion at that point, we want to emphasize that Corollary 4 does not imply the non-existence of some positive constant ε when dealing with complex-valued ITEs ω admitting a non-vanishing imaginary part.
Theorem 5. Let ω be the smallest real-valued ITE for the scatterer D with constant density ρ = 1. If ρ > 1 is large or 0 < ρ < 1 is small enough, where the corresponding thresholds depend only on D and the Lamé parameters µ, λ, we have the relations Proof. For the sake of presentation we will assume that ρ > 1 since the case 0 < ρ < 1 works structurally similar. Because u and v can be expressed each in terms of their difference w := u−v by u = (∆ * w+ρω 2 w)/(ω 2 (ρ−1)) and v = (∆ * w+ω 2 w)/(ω 2 (ρ−1)), the basic idea of our proof will be to exploit isometry of the Fourier transform with respect to the single field w to obtain a well-behaved algebraic integrand in terms of ρ. The fact w ∈ H 2 0 (D) then shows that u, v, w extend naturally by zero outside of D so Plancherel's identity gives In the last step we employed the pointwise estimate with c := min{µ, λ + 2µ} > 0 inherited from strong ellipticity. As will be shown later, we have for the smallest ITE ω = ω(ρ) with density ρ that ρω 4 → 0 and being the disc centered at the origin of radius R(ρ), we split the integral above and exploit monotonicity of p to The first summand can be made arbitrarily small in terms of ρ because of where L 2 (D) denotes the two-dimensional Lebesgue measure of D. Putting everything together, we finally obtain for ρ large enough due to the decay of R(ρ).
As announced for the latter, it remains to show that ρω 4 → 0 as ρ → ∞, where ω is the smallest real-valued ITE of D with density ρ. According to Corollary 1 from [2], the magnitude of ω can be bounded from above by the smallest ITE from any included disc B r ⊂ D which thus amounts to show the asymptotics assertion just for the unit disc as scatterer. In this case, as was derived in the appendix of [10] with an ansatz for a purely compressional wave instead of its orthogonal shear part that we will consider now for completion, ω is some ITE corresponding to a radial symmetric eigenfunction if and only if where J 1 is the first Bessel function of order one. This condition can be restated as finding roots ω of the piecewise continuous function which can be recursively decomposed into .
Let j 1 < j 2 be the two smallest positive roots of J 1 and choose ρ > 1 large enough to have j 2 < j 1 √ ρ. Then set ω 1 := j 1 µ/ρ as well as ω 2 := j 2 µ/ρ and observe that g is singular at those points, but continuous in between. Also, those poles have different signs according to which follows from the basic facts that J 1 < 0 in (j 1 , j 2 ), J ′ 1 (j 1 ) < 0 but J ′ 1 (j 2 ) > 0 and that both f (ω 1 ), f (ω 2 ) are finite. Therefore we can make use of the intermediate value theorem which guarantees for any large ρ a root ω of g fulfilling ω 1 ≤ ω ≤ ω 2 or equivalently the uniform bound j 2 1 µ ≤ ρω 2 ≤ j 2 1 µ as ρ → ∞. In particular, ρω 4 → 0 for the same limit procedure which finally proves our lemma.
Finally, we cite a result from [14,Corollary 6] that summarizes the previous findings for our conceptual ITE recovery approach.

Corollary 6. Let the conditions of Theorem 2 hold for
Then, for sufficiently large m ∈ N, we have (modulo the relabeling of the weakly convergent subsequence) where C > 0 depends on the boundary data of u (or equivalently v) and the magnitude of (10). (10) that ε > 0 in Lemma 3 uniformly for m large enough and thus the existence of C > 0.
Altogether, we have proved that the conditioned process from Theorem 2 for detecting ITEs is spurious free as limiting procedure and provided a constrained a posteriori estimate for the eigenvalue approximation accuracy at each step m. To meet the conditions therein, we will construct "simple" trial functions in the next section that are unspecified in (4) so far.

The abstract setting
In this section we want to focus on how to generate interior solutions to a given PDE on the basis of its fundamental solution. A fundamental solution for the free space Navier system (1) with arguments x = y ∈ R 2 , constant coefficients µ, λ, ̺ and unknown eigenfrequency parameter ω is given by with wave numbers uniformly in all directions of the radial distance r from the origin. In particular, freezing one of the arguments in (11) such as y without loss of generality, these properties apply column-wise (and by symmetry of the fundamental matrix also row-wise) to x → Φ ̺ω 2 (x, y) in R 2 \{y} for any source point y ∈ R 2 and density ̺ = const. Choosing some simply closed and sufficiently smooth contour Γ ∈ R 2 \D, called the artificial or source boundary, we easily see that any coefficient function c : Γ → C 2 would generate a smooth solution of (1) in D by the continuous superposition which is the starting point for the MFS. We now refine from (4) and prove in the following that the resulting approximation space is still sufficiently dense to finally recover exact eigenfunctions via boundary control. For this, we need some prerequisites concerning the construction of fundamental solutions for higher order PDEs. It also shows that the MFS ansatz in (12) for approximating u and v separately is equivalent to recovering their difference u − v via its synthesized fundamental solution.
Proof. Let ϕ ∈ C ∞ c (D, C 2 ) be an arbitrary bump function and set Φ y ρω 2 ,ω 2 (x) := Φ ρω 2 ,ω 2 (y − x) for some fixed y ∈ R 2 (and likewise for its generating fundamental solutions). Then we can easily check which proves the lemma.
We are now ready to prove that for any eigenfunction pair (u, v) of (2) with eigenfrequency ω lying in the first complex quadrant we can find approximations in H(ω) with arbitrarily small boundary misfits. In particular, all the recovery criteria from Theorem 2 can be satisfied for corresponding ITEs.
Theorem 8. Let w ∈ H 2 (D, C 2 ) be any distributional solution to the fourth order equation (∆ * + ω 2 )(∆ * + ρω 2 )w = 0 for some frequency ω such that 0 ≤ arg(ω) < π/4 (holds even more generally for Im ω ≥ 0). Then there exists a sequence of elements Proof. Assume w ∈ H −2 (D, C 2 ), where the latter denotes the negative Sobolev space with compact support in D identifying the dual space of H 2 (D, C 2 ), is chosen such that for all (c 1 , c 2 ) ∈ L 2 (Γ, C 2 ) × L 2 (Γ, C 2 ). The integral expression is hereby overloaded in notation with the corresponding duality pairing and by definition of Φ ρω 2 ,ω 2 from the previous lemma we see that the kernel of w contains functions of the form (u − v) with (u, v) ∈ H(ω). Therefore, if we can show that (14) implies for every distributional solution w * ∈ H 2 (D, C 2 ) of (∆ * + ω 2 )(∆ * + ρω 2 )w * = 0, the Hahn-Banach theorem would yield the desired density claim since no other extension is possible. For this, we define the auxiliary functions w : , where the gain of regularity in D results from the fact that the convolution with Φ ρω 2 ,ω 2 or Φ ρω 2 are pseudodifferential operators of order −2 due to their logarithmic singularity type. Rewriting (14), we obtain Γ c 1 (y) · w(y) + c 2 (y) · v(y) ds(y) = 0 for all (c 1 , c 2 ) ∈ L 2 (Γ, C 2 ) × L 2 (Γ, C 2 ) and setting one of the coefficient functions to zero, respectively, we may conclude that w |Γ = v |Γ = 0. Using the pointwise estimate , and similarly for v s , while employing standard differentiation properties and decay estimates for the resulting Hankel functions expansion within the first norm on the right, we deduce that Kupradze's radiation conditions are completely inherited by v as convolution of the correspondingly radiating fundamental solution and some compactly supported distribution. By uniqueness of the exterior Navier problem for Im ω ≥ 0 and (∆ * + ρω 2 )v = w in the sense of distributions, we may then conclude that v = 0 outside of Γ. Due to analyticity, v even needs to vanish completely in R 2 \D because the right hand side of the Navier equation is zero here by assumption. Similarly, we want to prove in the following that w ∈ H 2 0 (D) for justifying its role as a valid test function later: Using the definition of Φ ρω 2 ,ω 2 , direct calculations show the distributional relations The same uniqueness and analyticity reasoning as for v above but with frequency ω 2 now yields w = 0 in R 2 \D again and by a bootstrap argument due to v ∈ L 2 (D, C 2 ) we then conclude w ∈ H 2 loc (R 2 , C 2 ) or equivalently w ∈ H 2 0 (D). Therefore we can find a sequence of bump functions Taking then any distributional solution w * ∈ H 2 (D, C 2 ) of (∆ * + ω 2 )(∆ * + ρω 2 )w * = 0 as mentioned in (15), we may finally compute where we especially incorporated in the first step that (∆ * + ρω 2 )(∆ * + ω 2 )w = w in the sense of distributions according to our original definition of w. Since w * was an arbitrary homogeneous solution, the desired density result for the interior domain is thereby proven. An application of standard trace theorems eventually takes the approximation result over to the boundary of D in corresponding norms.
In the next section we will focus on the numerical implications of these abstract findings and especially present its associated recipe of how to approximate ITEs in practice.

The numerical setting
The MFS aims to discretize (12) in terms of a Riemann sum with unknown integration weights determined in some optimization process which we derive in the following. For arbitrary frequencies ω and some approximation order m ∈ N we focus on trial functions defined on D of the form where {y 1 , . . . , y m } ⊂ Γ are preselected source points and c u , c v ∈ C 2×m are up to now unspecified block coefficient vectors. Note (u m , v m ) are actually pseudo-elements of H(ω) by using delta-distribution-type coefficient functions c u and c v along Γ, respectively. However, it can be readily seen that they are dense in H(ω) for m → ∞ in any interior Sobolev norm, cf. [14] where this was shown for the acoustic framework. In order to characterize the optimal pair ( c u , c v ) ∈ C 2×m × C 2×m , we will now develop an MFS-based scheme that is supposed to meet the criteria listed in Theorem 2 from a numerical perspective. It uses the collocation concept as a discretized representative for the Sobolev-norm-conditioned quantities involved. First, we treat the vanishing boundary misfit condition (iii). Therefore we pick m collocation points {x 1 , . . . , x m } ⊂ ∂D and consider for fixed ω the constrained minimization of This optimization must be subject to criterion (ii) in a coherently discretized way for which we will take quadrature-like sample points { x 1 , . . . , x m I } ⊂ D with m I being large but fixed and demand That is, we scale all eigenfunction candidate pairs to some almost constant but nonvanishing discretized interior norm to guarantee (u m , v m ) = 0 for all ω. Both sums above can be restated in compact matrix form, for which we define where the upper dense rows correspond to the boundary control matrices with 1 ≤ i, j ≤ m and co-normal derivatives affecting only the first arguments on the right hand side of involved definitions. A similar pattern applies for σ( Φ ̺ω 2 )ν. The diagonal lower part of (16) then similarly embodies the interior samples where now 1 ≤ i ≤ m I and again 1 ≤ j ≤ m. The benefit of this matrix reformulation is that we can now perform a QR-factorization with Q ∈ C 4m×4m and Q I ∈ C 4m I ×4m , which conveniently reflects the above minimalboundary-to-maximal-interior coupling through the unitary property of Q M . This observation goes back to Betcke and Trefethen, see [4], who analyzed Dirichlet eigenvalues with a slightly different ansatz but still via boundary control. Our common ingredient left is then to find local minima of Those solutions that are (almost) roots for the smallest singular value σ 1 (ω), see Figure  1 (left), will be denoted by ω m to trace back to the underlying trial space dimension of (u m , v m ). Our described solution algorithm will be called modified MFS and its output, ω m , approximate ITEs. In the spirit of (i) from Theorem 2 we hope them to converge for m → ∞ also in practice and Lemma 3 would then provide a measure for their convergence speed in terms of the continuous analogon of σ 1 (ω m ). The next section will demonstrate this for some exemplary scatterers.

Numerical results
We use the modified MFS from the previous section to compute some ITEs for a disc-, ellipse-, kite-and square-shaped scatterer D whose (collocation) boundaries are given by respectively. Note that D s is actually not covered by our ITE approximation analysis due to its corners, but still a very typical scatterer feasible for our algorithm since collocation is invisible with respect to the boundary regularity except for singular points that need to be excluded. Throughout we fix the constitutive parameters which were also used in the context of [10] and therefore serve as independent reference values for our exemplary MFS-based findings listed in Figure 2. Those approximate the first four real-valued ITEs from each scatterer and were individually obtained via the  2. Approximations of the first four real-valued ITEs (counting without multiplicity) for some exemplary scatterers obtained via the modified MFS with material parameter µ = 1/16, λ = 1/4, ρ = 4: One clearly sees that the more advanced the shape of the boundary becomes, the less ITE digits can be effectively recovered.
source point boundaries Γ e := 0.95 cos(t) 1.9 sin(t) = 1.9 · ∂D e , t ∈ [0, 2π) , The extracted scaling factors {2, 1.9, 1.6, 1.3} were preselected to follow the conclusions from [14] and approach unity the more scattering shapes seem to deviate from the disc that is considered as the most promising state. Having thus set all the necessary computational contours for the modified MFS, both the m-dependent collocation points {x 1 , . . . , x m } as well as the auxiliary sources {y 1 , . . . , y 2m } were distributed equiangular on their corresponding boundaries, e.g. equidistant with respect to t if our representation above allows. However, it is well known that the individual optimization of Γ itself and its source point distribution do have a noticeable impact on the MFS approximation quality as was shown in [1], but resulting in a more advanced non-linear problem for each scatterer in total. In contrast, concerning the interior points, m I = 10 of them were fixed randomly in a centered disc with radius 0.5 as their overall location and number turns out not to affect the ITE output significantly.
With these input arrangements and the focus on real-valued ITEs first, the modified MFS can be successfully exploited by incrementing m within 40 ≤ m ≤ 80, where the lower bound just provided an averaged value for when to expect good results. However, exceeding the given D-specific threshold for our setup, the graph of (17) always began to suffer drastically from the intrinsic ill-conditioning effects of the discretization matrix M(ω) in (16) via impeding oscillations and finally lead to unreliable approximations apart, see Figure 1 (right). Obtained via the regime in between, we believe our given cut-off results from Figure 2 to be correct up to that point modulo round-off errors compared to the exact ITE mantissa since they correspond to the nonfluctuating digits within the modified MFS output ω m when m ր 80. Some exact reference values for D d corresponding to rotational-symmetric eigenfunctions can be obtained by computing roots of (9) which even confirms our approximation for its smallest ITE to be exact up to machine precision, see Figure 3 for the convergence history. In general, the recoverable accuracy strongly correlates with the scattering shape and especially with the existence of corners. These technical observations mostly agree with the ones from the acoustic case analyzed in [14,13] and only differ in the retarded yet D-specific convergence regime in m for starting the actual ITE approximations which was about m ≥ 20 before. It manifests the fact that the eigenfunctions from elasticity are vector-valued and thus numerically slightly more expensive. As in the particular scope of this paper, we also want to present some approximations from the complex-valued eigenvalue spectrum. Therefore Figure 4 plots (17) in the complex range 0 ≤ Re ω ≤ 2.5, −2 ≤ Im ω ≤ 2 for the disc D d and the square D s since only their distribution of ITEs seems not that dense among our analyzed scatterers according to Figure 2 and is thus better suited for broad contour visualizations. We can immediately see that ITEs are correctly displayed in conjugated pairs and as an example the closest of them with respect to the positive imaginary axis are computed for the two scatterers to be 1.8624 + 0.3104i and 1.987178187576699 + 0.283125784408650i, respectively. The correctness of all given ITE digits for D d can even be confirmed again by comparing with the correspondingly coinciding root of (9). Regarding accuracy, the modified MFS does not show any remarkable difference for its extension to the upper half space of the complex plane. However, although the plot indicates a certain axial symmetry with below, one clearly sees that inaccuracies due to large condition numbers of the underlying matrix (16) dominantly propagate from the lower half space upwards for increasing m, giving a first hint where the oscillations in Figure 1 (right) originate from. Fortunately, ITEs always arise in conjugated pairs and because of Theorem 8 we even know that restricting to complex numbers with non-negative imaginary part within our investigations, as already implemented in our initial ansatz (4) and (13), is not limiting at all.  To get a rough impression about how competitive our method is compared to others, we use the approximation results for real-valued k 2 from the paper [10] in which the authors also analyzed the scatterers D d and D s but via special finite element methods. With a mesh size of h ≤ 0.0125, they obtained the approximations k 2 d ≈ 2.110723 and k 2 s ≈ 1.9428781 for the smallest ITE squared of the disc and the square, respectively. Allowing for possible round-off deviations in the last digit within our corresponding cutoff values, we obtain the tolerance ranges k 2 d ∈ [2.106283380546506, 2.106283380546512] and k 2 s ∈ [1.9423, 1.9429] by setting at most m = 80 which may be converted via a boundary partitioning to a collocation point distance of the same order as the alternative mesh size. While the results regarding D s are of the same order for both solution approaches, the modified MFS clearly dominates in accuracy for D d and is therefore believed to do so for slightly perturbed scatterers D, too.

Conclusion
We analyzed the method of fundamental solution in a stabilized version for the computation of complex-valued elastic transmission eigenvalues in two dimensions. Our theoretical studies show that the short algorithm shall produce spurious-free results in the limit whose approximation error per step we could quantify in terms of some discretized residual output value. Our numerical experiments confirm these expectations in practice within some scatterer-specific collocation point regime for which the method's feasibility is unaffected from ill-conditioning pollutions. This was tested for a collection of scatterers whose boundaries are easily-parametrizable. In accordance with the conclusions from previous works that analyzed our algorithm in the context of related transmission problems but restricted to real-valued eigenvalues so far, the best results, including also the complex spectrum from now, can still be obtained for the unit disc. The more the actual scattering shape then deviates from the disc, the less accurate approximations are finally achievable while even more collocation points are needed in total. However, depending on the choice of fundamental solution for generating the trial functions, different areas of the complex spectrum (in our case the lower half space) should be avoided as initial guess input for the algorithm due to ill-conditioning effects which are, however, not that restrictive because all eigenvalues arise in conjugated pairs. Generally, since the source boundary necessary for the MFS setup was individually preselected by intuition for simplicity, even more promising approximations can be achieved by its optimization which was, however, not covered in the scope of this paper. As a conclusion, the modified MFS is especially effective for quite regular domains and dominates here over many competitive methods in the general context of transmission eigenvalue problems.