# Study of a degenerate non-elliptic equation to model plasma heating

Patrick Ciarlet, Maryna Kachanovska, Etienne Peillon

► To cite this version:

Patrick Ciarlet, Maryna Kachanovska, Etienne Peillon. Study of a degenerate non-elliptic equation to model plasma heating. 2023. hal-04152267

HAL Id: hal-04152267

https://inria.hal.science/hal-04152267

Preprint submitted on 5 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Study of a degenerate non-elliptic equation to model plasma heating

Patrick CIARLET*

Maryna KACHANOVSKA*

Étienne PEILLON*

July 5, 2023

Abstract

In this manuscript, we study solutions to resonant Maxwell's equations in heterogeneous plasmas. We concentrate on the phenomenon of upper-hybrid heating, which occurs in a localized region where electromagnetic waves transfer energy to the particles. In the 2D case, it can be modelled mathematically by the partial differential equation $-\text{div}(\alpha\nabla u) - \omega^2 u = 0$, where the coefficient $\alpha$ is a smooth, sign-changing, real-valued function. Since the locus of the sign change is located within the plasma, the equation is non-elliptic, and degenerate. On the other hand, using the limiting absorption principle, one can build a family of elliptic equations that approximate the degenerate equation. Then, a natural question is to relate the solution of the degenerate equation, if it exists, to the family of solutions of the elliptic equations. For that, we assume that the family of solutions converges to a limit, which can be split into a regular part and a singular part, and that this limiting absorption solution is governed by the non-elliptic equation introduced above. One of the difficulties lies in the definition of appropriate norms and function spaces in order to be able to study the non-elliptic equation and its solutions. As a starting point, we revisit a prior work [13] on this topic by A. Nicolopoulos, M. Campos Pinto, B. Després and P. Ciarlet Jr., who proposed a variational formulation for the plasma heating problem. We improve the results they obtained, in particular by establishing existence and uniqueness of the solution, by making a different choice of function spaces. Also, we propose a series of numerical tests, comparing the numerical results of Nicolopoulos et al to those obtained with our numerical method, for which we observe better convergence.

# 1 Introduction

When an electromagnetic wave is sent inside a plasma, it can transfer energy to the particles to produce plasma heating in a localized region. This phenomenon is related to the so-called resonant waves. In a magnetized plasma set in a region of $\mathbb{R}^3$ under an exterior constant magnetic field $B_0 = (0,0,B_0)$, the electric field $E$ and the magnetic induction $B$ are governed by the following Maxwell system in the time-harmonic regime ($\omega > 0$):

$$ \begin{cases} \text{curl } E = i\omega B, \\ \text{curl } B = -\frac{i\omega}{c^2} \varepsilon E, \end{cases} \quad \text{with} \quad \varepsilon = \begin{pmatrix} \alpha & i\delta & 0 \\ -i\delta & \alpha & 0 \\ 0 & 0 & \beta \end{pmatrix} \qquad (1.1) $$

<footnote>*POEMS, CNRS, INRIA, ENSTA Paris, Institut Polytechnique de Paris, 828 Bd des Maréchaux, 91762 Palaiseau Cedex, France. E-mail addresses: patrick.ciarlet@ensta-paris.fr, maryna.kachanovska@inria.fr, etienne.peillon@ensta-paris.fr</footnote>
where $\epsilon$ is the cold plasma dielectric tensor [15, 11]. For a single species plasma, the coefficients $\alpha$, $\delta$ and $\beta$ have the form:

$$ \alpha(\mathbf{x}) = 1 - C_{\alpha,\omega}\mathcal{N}_e(\mathbf{x}), \qquad \delta(\mathbf{x}) = C_{\delta,\omega}\mathcal{N}_e(\mathbf{x}), \qquad \beta(\mathbf{x}) = 1 - C_{\beta,\omega}\mathcal{N}_e(\mathbf{x}), \qquad (1.2) $$

where $\mathcal{N}_e(\mathbf{x})$ is the particle density (assumed to be a continuous function) and $C_{\alpha,\omega}$, $C_{\delta,\omega}$ and $C_{\beta,\omega}$ are three real-valued frequency-dependent constants which do not vary in space. We consider an upper hybrid resonance in the plasma, see [15, Chapter 2-6] and recent works [11, 7, 4, 5, 8, 6, 12, 14, 13], which is characterized by the fact that $\alpha = 0$ on some curve inside the region. Like in the above cited works, we will be particularly interested in the cases when the density $\mathcal{N}_e(\mathbf{x})$ is s.t. the sign of $\alpha$ changes continuously between subregions separated by an interface.

We will consider a simplified model, which we derive in detail below. Let $(e_1, e_2, e_3)$ be an orthonormal basis of $\mathbb{R}^3$, with $(x_1, x_2, x_3)$ the normalized orthogonal coordinates. With obvious notations, $E = E_1e_1 + E_2e_2 + E_3e_3$, etc. We will assume in this manuscript that all quantities are independent of $x_3$ (a variable corresponding to the direction of the exterior magnetic field), i.e., $E_1 = E_1(x_1, x_2)$, etc. Then, because of the block diagonal structure of $\epsilon$, we observe a decorrelation between $E_3$, $B_\perp = B_1e_1 + B_2e_2$ on one hand, and $E_\perp = E_1e_1 + E_2e_2$, $B_3$ on the other. As a matter of fact, the Maxwell system (1.1) can be split into two independent systems, referred to as the system for the Ordinary mode:

$$ \left\{ \begin{array}{l} \mathrm{curl}_{\perp} E_3 = i\omega \mathbf{B}_{\perp}, \\ \mathrm{curl}_{\perp} \mathbf{B}_{\perp} = -\frac{i\omega\beta}{c^2} E_3, \end{array} \right. \qquad (\mathrm{O-mode}) $$

and the eXtraodinary mode:

$$ \left\{ \begin{array}{l} \mathrm{curl}_{\perp} \mathbf{E}_{\perp} = i\omega B_3, \\ \mathrm{curl}_{\perp} B_3 = -\frac{i\omega}{c^2} \epsilon_{\perp} \mathbf{E}_{\perp}, \end{array} \right. \qquad (\mathrm{X-mode}) $$

where $\epsilon_\perp = (\frac{\alpha}{-i\delta}, \frac{i\delta}{\alpha})$. The differential operators are defined by $\mathrm{curl}_\perp \mathbf{v} = \partial_1 v_2 - \partial_2 v_1$ and $\mathrm{curl}_\perp \mathbf{v} = \partial_2 v e_1 - \partial_1 v e_2$.

For the discussion that follows, we will need to introduce auxiliary notation. Let $\mathbf{x}_\perp = x_1e_1 + x_2e_2$, $\Delta_\perp v = \partial_{11}v + \partial_{22}v$, $\mathrm{div}_\perp \mathbf{v} = \partial_1 v_1 + \partial_2 v_2$ and $\nabla_\perp v = \partial_1 v e_1 + \partial_2 v e_2$. In this case we evidently have that $\mathrm{curl}_\perp = -R_{\pi/2} \nabla_\perp$ and $\mathrm{curl}_\perp = -\mathrm{div}_\perp R_{\pi/2}$, where $R_{\pi/2} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$ is the $\pi/2$ rotation matrix.

Let us now focus on the equations governing the scalar unknowns $E_3$ and $B_3$. The second-order PDE derived from the system for the Ordinary mode is $-\Delta_\perp E_3 = \frac{\omega^2\beta}{c^2} E_3$. In the case when the sign of $\beta$ changes continuously, this equation is reminiscent of an Airy equation, cf. [11].

On the other hand, the second-order PDE derived from the (X-mode) is

$$ \mathrm{div}_{\perp} (R_{\pi/2} \epsilon_{\perp}^{-1} R_{\pi/2} \nabla_{\perp} B_3) = \frac{\omega^2}{c^2} B_3. $$

We will assume that the tensor $\epsilon_\perp$ is invertible everywhere in the region, more precisely, that $\alpha^2(\mathbf{x}_\perp) - \delta^2(\mathbf{x}_\perp) \neq 0$, and thus the above expression is well-defined. Let us define the two-by-two tensor

$$ \omega(\mathbf{x}_\perp) := c^2 R_{\pi/2} \epsilon_{\perp}^{-1} R_{\pi/2} = \frac{c^2}{\delta^2 - \alpha^2} \begin{pmatrix} \alpha & -i\delta \\ i\delta & \alpha \end{pmatrix}. $$
Then, since it holds that $\alpha$ and $\delta$ depend on the space variable $\mathbf{x}_{\perp}$ only via the density of the plasma $N_e(\mathbf{x}_{\perp})$, $\alpha$ and $\delta$ have the same level curves.

As discussed before, we assume that the coefficient $\alpha(\mathbf{x}_{\perp})$ vanishes on some interface $I$. In view of the last remark, the tensor $\omega(\mathbf{x}_{\perp})$ is constant on $I$, and it holds that $\omega = iA$ there, with $A$ being a real-valued skew-symmetric matrix. Denoting by $\delta^{+}$ the value of $\delta$ on the interface $I$, we remark that $\delta(\mathbf{x}_{\perp}) - \delta^{+} = -\delta^{+}\alpha(\mathbf{x}_{\perp})$. From this it follows that we can decompose $\omega(\mathbf{x}_{\perp})$ as

$$ \omega(\mathbf{x}_{\perp}) = -\alpha_0(\mathbf{x}_{\perp})\mathbb{H}(\mathbf{x}_{\perp}) + iA, $$

where $\alpha_0 = \frac{c^2\alpha}{\alpha^2 - \delta^2}$, and $\mathbb{H}(\mathbf{x}_{\perp})$ is a Hermitian matrix given by

$$ \mathbb{H} = \begin{pmatrix} 1 & -i(\delta + \alpha/\delta^{+}) \\ i(\delta + \alpha/\delta^{+}) & 1 \end{pmatrix}. $$

In what follows, we will assume that $\mathbb{H}(\mathbf{x}_{\perp})$ is positive definite in the whole computational region, which in particular requires that

$$ |\delta(\mathbf{x}_{\perp}) + \alpha(\mathbf{x}_{\perp})/\delta^{+}| < 1. $$

Since $A$ is skew-symmetric, $\mathrm{div}_{\perp}(A\nabla_{\perp}B_3) = 0$, so the second-order PDE governing $B_3$ becomes

$$ -\mathrm{div}_{\perp}(\omega_0 \nabla_{\perp} B_3) - \omega^2 B_3 = 0. \qquad (1.3) $$

We suppose that the electron density $N_e$ is $C^2$-regular, so that $\alpha_0$ and $\mathbb{H}$ are also $C^2$-regular. We assume that the interface $I = \{\alpha_0(\mathbf{x}_{\perp}) = 0\}$ is a $C^1$-loop (without self-intersections), and that $|\alpha_0(\mathbf{x}_{\perp})|$ behaves like $\mathrm{dist}(\mathbf{x}_{\perp}, I)$ in its neighborhood.

Hence, considering the model derived from the eXtraodinary mode with unknown $B_3$ in the neighborhood of the interface leads to a degenerate elliptic PDE. How to solve this equation will be the manuscript's goal. This problem was originally investigated by A. Nicolopoulos, M. Campos Pinto, B. Després and P. Ciarlet Jr. in [13], where a numerical method based on a mixed variational formulation was proposed. In this manuscript we revisit this method and propose some improvements. The outline is as follows. We introduce in section 2 the geometrical and functional settings used in this paper. We also discuss the a priori regularity assumptions on the solution using the limiting absorption principle. Then in section 3, we investigate the construction and the discretization of the method introduced in [13]. Numerical experiments indicate that the discrete solution does not converge to the exact solution. Then, in section 4, we propose a new method, designed to overcome the limitations of the original one. Its properties are studied in section 5. Finally, we discretize the new method in section 6: numerical experiments show drastic improvement over the original method.

## 2 Mathematical setting

We consider a bounded Lipschitz domain $D$ in $\mathbb{R}^2$. Let $\lambda > 0$ and $f \in L^2(\partial D)$. We study the following boundary value problem:

$$ \left\{ \begin{array}{ll} \mathrm{div}_{\perp}(\omega_0 \nabla_{\perp} B_3) - \omega^2 B_3 = 0 & \text{in } D, \\ (\omega_0 \nabla_{\perp} B_3) \cdot \mathbf{n}_{\perp} + i\lambda B_3 = f & \text{on } \partial D, \end{array} \right. \qquad (2.1) $$
where $\mathbf{n}_{\perp}$ denotes the outward unit vector field to $\partial D$. According to our model, it holds that $\alpha_0(\mathbf{x}_{\perp}) = \alpha_0(\mathbf{x}_{\perp})\mathbb{H}(\mathbf{x}_{\perp})$, where the scalar field $\alpha_0$ and the hermitian matrix field $\mathbb{H}$ are $\mathcal{C}^2(\bar{D})$-regular. We set $D_p = \{\mathbf{x}_{\perp} \in D : \alpha_0(\mathbf{x}_{\perp}) > 0\}$, $D_n = \{\mathbf{x}_{\perp} \in D : \alpha_0(\mathbf{x}_{\perp}) < 0\}$, and recall that the interface $I = \{\mathbf{x}_{\perp} \in D : \alpha_0(\mathbf{x}_{\perp}) = 0\}$ is a $\mathcal{C}^1$-loop (without self-intersections). We assume here that meas$(D_{p,n}) > 0$, and that $I$ does not intersect $\partial D$. Observe that outside every neighborhood of $I$ we are solving a classical second-order elliptic PDE with smooth coefficients. Hence, following the classical theory, we shall look for a solution that belongs to $H^1$ outside this neighborhood. To fix ideas, we consider the case where that $D$ is a tubular neighborhood of $I$. Finally, we recall that $|\alpha|$ behaves like $\text{dist}(\cdot, I)$ in a neighborhood of the interface.

Like in [13], we focus on the problem (2.1) posed in the neighborhood of the interface. Let $\Omega = (-a, a) \times (0, L)$ be a subset of $\mathbb{R}^2$, with the normalized orthogonal coordinates $(x, y)$. Introduce the volume preserving bijective transform $\psi : (x, y) \to \mathbf{x}_{\perp}$ (more precisely, $|\det(D\psi)| = 1$) which maps $\bar{\Omega}$ to $\bar{D}$ with the following properties, see Figure 1:

* the preimage of the interface $I$ is the straight line $\Sigma = \{0\} \times [0, L]$;

* the preimage of the subregion $D_n$ is the rectangle $\Omega_n = (-a, 0) \times (0, L)$;

* the preimage of the subregion $D_p$ is the rectangle $\Omega_p = (0, a) \times (0, L)$;

* the preimage of $\partial D_n \setminus I$ is the straight line $\{-a\} \times [0, L)$;

* the preimage of $\partial D_p \setminus I$ is the straight line $\{a\} \times [0, L)$;

* the image of $(-a, a) \times \{0\}$ is equal to the image of $(-a, a) \times \{L\}$.

<figure><img src="image_2.png" /><figcaption>Figure 1: [Left] The tubular neighborhood $D$ of $I$. [Right] The domain $\Omega = (-a, a) \times (0, L)$. [Center] The transform $\psi : \bar{\Omega} \to \bar{D}$ with $\psi(\Sigma) = I$, $\psi(\Omega_{p,n}) = D_{p,n}$ and $\psi(\Gamma_{p,n}) = \partial D_{p,n} \setminus I$.</figcaption></figure>

We split the boundary of $\Omega$ into 4 components:

$$ \Gamma_p = \{a\} \times [0, L), \quad \Gamma_n = \{-a\} \times [0, L), \quad \Gamma_1 = (-a, a) \times \{0\}, \quad \Gamma_2 = (-a, a) \times \{L\}. $$

Via the mapping $\psi^{-1}$, the problem is recast onto the rectangle $\Omega = (-a, a) \times (0, L)$. Let us introduce

$$ \underline{\alpha}(x, y) = [D\psi(x, y)]^{-1}\alpha_0(\mathbf{x}_{\perp})[D\psi(x, y)]^{-t} \text{ with the correspondence } \mathbf{x}_{\perp} = \psi(x, y). $$
Because the transform is volume preserving (see e.g., [2, §2.1.3]), we find that the BVP (2.1) is equivalently reformulated with the unknown $u = B_3 \circ \psi$ (and data still denoted by $f$) as

$$ egin{cases} -\operatorname{div}(\underline{\omega}\nabla u) - \omega^2 u = 0 & \text{in } \Omega, \\ \underline{\omega}\nabla u \cdot \mathbf{n} + i\lambda u = f & \text{on } \Gamma_n \cup \Gamma_p, \\ u(x, 0) = u(x, L), \quad (\underline{\omega}\nabla u \cdot \mathbf{e}_y)(x, 0) = (\underline{\omega}\nabla u \cdot \mathbf{e}_y)(x, L), \quad x \in (-a, a). \end{cases} \quad (2.2) $$

Above, the divergence and gradient operators are the classical 2D operators, while $\mathbf{n}$ denotes the outward unit vector field to $\partial\Omega$. The last conditions account for periodicity.

The requirement that the transform is volume preserving does not reduce the scope of the study. Indeed, if the transform is not volume preserving, then the definition of $\underline{\omega}(x, y)$ in (2.2) changes to $\underline{\omega}(x, y) = [D\psi(x, y)]^{-1}\omega_0(\mathbf{x}_{\perp})[D\psi(x, y)]^{-t}\det(D\psi(x, y))$. This applies for both the second-order term in $\Omega$, and for the flux-like term on $\partial\Omega$. Regarding the zero-order terms, we note first that, in $\Omega$, the zero-order term becomes $-\omega^2 u(x, y) \det(D\psi(x, y))$. But, since there exist two constants $C_{\min}, C_{\max} > 0$ such that $C_{\min} \le |\det(D\psi(x, y))| \le C_{\max}$ for all $(x, y) \in \bar{\Omega}$, one simply needs to add a smooth, bounded above and away from below, weight to the zero-order term in the domain $\Omega$. For the zero-order term on the boundary $\partial\Omega$, one has to make a similar modification, ending up with a second weight with the same properties.

In what follows, we make two simplifying technical assumptions. First, that $\underline{\omega}$ is pointwise proportional to the identity matrix, that is $\underline{\omega}(x, y) = \alpha(x, y) (\begin{smallmatrix} 1 & 0 \\ 0 & 1 \end{smallmatrix})$ everywhere in $\bar{\Omega}$. Here we use, with an abuse of notation, the same letter for the new coefficient $\alpha(x, y)$ and the coefficient $\alpha(\mathbf{x}_{\perp})$ in (1.1); those however are not to be confused. In this situation, the model can be recast as

$$ egin{cases} -\operatorname{div}(\alpha\nabla u) - \omega^2 u = 0 & \text{in } \Omega, \\ \alpha\partial_n u + i\lambda u = f & \text{on } \Gamma_n \cup \Gamma_p, \\ u(x, 0) = u(x, L), \quad (\alpha\partial_y)u(x, 0) = (\alpha\partial_y)u(x, L), \quad x \in (-a, a). \end{cases} \quad (2.3) $$

Second, we assume that $\partial_x\alpha(0, y) > 0$ for every $y \in (0, L)$. Because the coefficient $\alpha$ is a scalar, we observe that the interface $\Sigma$ is now described by $\{(x, y) : \alpha(x, y) = 0\}$, while the two subdomains are respectively described by $\Omega_p = \{(x, y) : \alpha(x, y) > 0\}$ and $\Omega_n = \{(x, y) : \alpha(x, y) < 0\}$.

There remains to specify the requested regularity of $u$, so as to allow for the modelling of plasma heating. In this manuscript, we look for limiting absorption solutions of the above problem, namely, we look for $u$ being an $L^2$-weak limit of $u^{\nu}$, as $\nu \to 0^{+}$, where $u^{\nu}$ is the unique solution of the (coercive) limiting absorption problem:

$$ egin{cases} \text{find } u^{\nu} \in H^1(\Omega) \text{ s.t.} \\ -\operatorname{div}((\alpha + i\nu)\nabla u^{\nu}) - \omega^2 u^{\nu} = 0 & \text{in } \Omega, \\ \alpha\partial_n u^{\nu} + i\lambda u^{\nu} = f & \text{on } \Gamma_n \cup \Gamma_p, \\ u^{\nu}(x, 0) = u^{\nu}(x, L), \quad (\alpha\partial_y)u^{\nu}(x, 0) = (\alpha\partial_y)u^{\nu}(x, L), \quad x \in (-a, a). \end{cases} \quad (2.4) $$

**Remark 2.1.** *Remark that we study the family of problems (2.4), where the absorption appears in the principal part of the operator, rather than in zero-order terms, as in more classical cases. Indeed, in real plasmas, due to the presence of dissipation, the frequency-dependent cold plasma dielectric tensor $\epsilon$ has a small imaginary part, cf. [9, Chapters 15.4.4, 15.5] (which in the model (2.4) we study is modelled by $\nu > 0$).*
Up to our knowledge, the limiting absorption principle can be justified in 1D, as well as for particular values of $\alpha(x, y)$ in slab geometries, cf. [7, 6]. Let us provide an illuminating example whose goal is two-fold. On one hand, we will show how the limiting absorption principle leads to the occurrence of a logarithmic singularity in the solution. On the other hand, we will highlight the difficulty in the choice of the functional framework that would accommodate such singular solutions. Consider the 1D boundary-value problem: given $c_1, c_2 \in \mathbb{R}$, find $u$ solving

$$-(xu')' = 0 \text{ on } \mathcal{I} := (-a, a), \quad u(-a) = c_1, u(a) = c_2.$$

We could have looked for the solution to the above problem in the space $H_{1/2}^1(\mathcal{I}) = \overline{C^\infty(\bar{\mathcal{I}})}^{\|\cdot\|_{|\alpha|^{1/2}}}$, where

$$\|v\|_{|\alpha|^{1/2}}^2 = \int_{\mathcal{I}} |v|^2 + \int_{\mathcal{I}} |x||v'|^2.$$

In the definition of the norm $\|\cdot\|_{|\alpha|^{1/2}}$, the choice of the weight in front of $|v'|^2$ is motivated by the ODE itself, since after multiplying it by any admissible function $v$ supported away from 0 in either $\mathcal{I}_+ = (0, a)$ or $\mathcal{I}_- = (-a, 0)$ and integrating by parts, one gets a volume term like $\pm \int_{\mathcal{I}_\pm} |x|u'v'$. In this space, the bilinear form associated to the above equation is obviously continuous. A straightforward computation shows that in this case $u$ is a piecewise constant function

$$u = c_1 \text{ on } (-a, 0), \quad u = c_2 \text{ on } (0, a).$$

We see that the above solution does not contain any singularity except for the jump at the origin.

On the other hand, we can have a look at the limiting absorption solution to the above equation, where the absorption solution solves

$$-( (x + i\nu)(u^{\nu})' )' = 0 \text{ on } \mathcal{I}, \quad u^{\nu}(-a) = c_1, u^{\nu}(a) = c_2.$$

In particular, for each $\nu > 0$ the $H^1(\mathcal{I})$-solution to this problem is unique. With $z \mapsto \log z$ defined by its principal value (i.e., $\log(z) = \log|z| + i \operatorname{Arg}(z)$, $\operatorname{Arg} z \in (-\pi, \pi])$, we then compute the solution to the above equation

$$u^{\nu} = a_{\nu} \log(x + i\nu) + b_{\nu}, \quad \text{with } a_{\nu} = \frac{c_2 - c_1}{\log(a + i\nu) - \log(-a + i\nu)} \quad \text{and} \quad b_{\nu} = c_2 - a_{\nu} \log(a + i\nu).$$

Since $\log(x + i\nu) \xrightarrow{\nu \to 0^+} \log|x| + i\pi 1_{x<0}$ for $x \in \mathbb{R}^*$, the limiting absorption solution $u^+(x) = \lim_{\nu \to 0^+} u^\nu$ is given by the pointwise limit

$$u^+(x) = a_+(\log|x| + i\pi 1_{x<0}) + b_+, \quad a_+ = \frac{c_1 - c_2}{i\pi}, \quad b_+ = c_2 - a_+ \log a.$$

Note that $\int_{\mathcal{I}} |x| |(u^+)'|^2 = +\infty$ as soon as $a_+ \neq 0$. We thus see the difference between the two solutions $u \in H_{1/2}^1(\mathcal{I})$ and $u^+$: the first one has a jump singularity only, while the second one has both a logarithmic and a jump singularities. Interestingly, from the physics viewpoint, we know that it is the logarithmic singularity that is responsible for plasma heating phenomenon [7, 6]. Therefore, for the 2D model (2.3), we focus on using the latter one, that is solutions that include the jump and the logarithmic singularities ($\log|x| + i\pi 1_{x<0}$). From now on, we use the notation

$$S(x) := \log|x| + i\pi 1_{x<0}$$
to describe those singularities.

Let $C_{per,y}^{\infty}(\bar{\Omega}_j) = \{v \in C^{\infty}(\bar{\Omega}_j) : \partial_y^m v(x, 0) = \partial_y^m v(x, L), \forall m\}, j \in \{p, n\}$. Introduce the two spaces $H_{1/2}^1(\Omega_j) = \overline{C_{per,y}^{\infty}(\bar{\Omega}_j)}^{\|\cdot\|_{|\alpha|^{1/2}}}$, $j \in \{p, n\}$, with associated norm

$$ \|v\|_{|\alpha|^{1/2}}^2 = \int_{\Omega_j} |v|^2 + \int_{\Omega_j} |\alpha||\nabla v|^2, \quad j \in \{p, n\}. $$

Defining the above spaces is motivated by the same observation as above: multiplying the second-order PDE by any admissible function $v$ supported either in $\Omega_p$ or $\Omega_n$ and integrating by parts, one gets a volume term like $\pm \int_{\Omega_{p,n}} |\alpha| \nabla u \cdot \nabla v$. So, in light of the 1D example, we now consider that one can recast the 2D model with solution $u$ as follows.

**Assumption 2.2.** The family of solutions $(u^\nu)_\nu > 0$ of (2.4) converges in $L^2(\Omega)$ to the limiting absorption solution $u^+ \in L^2(\Omega)$

$$ u^\nu \xrightarrow{\nu \to 0+} u^+. \qquad (2.5) $$

Moreover, $u^+$ can be represented as

$$ u^+ = u_{reg}^+ + u_{sing}^+ $$

where the pair $(u_{reg}^+, u_{sing}^+)$ is such that $u_{reg|\Omega_{p,n}}^+ \in H_{1/2}^1(\Omega_{p,n})$ and $u_{sing}^+(x, y) = g^+(y) S(x)$ with $g^+ \in H_{per}^1(\Sigma)$.

Starting from (2.5), it is easy to verify that the limiting absorption solution $u^+$ is governed by the 2D model

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\alpha \nabla u^+) - \omega^2 u^+ = 0 & \text{in } \Omega_{p,n}, \\ \alpha \partial_n u^+ + i\lambda u^+ = f & \text{on } \Gamma_{p,n}, \\ u^+(x, 0) = u^+(x, L), \quad (\alpha \partial_y) u^+(x, 0) = (\alpha \partial_y) u^+(x, L), & x \in (-a, a). \end{array} \right. \qquad (2.6) $$

In Appendix A, we prove further results for $u^+$.

We identify the function $u_{reg}^+$ with a pair

$$ \mathbf{u}^+ = (u_{reg|\Omega_p}^+, u_{reg|\Omega_n}^+) \in H_{1/2}^1(\Omega_p) \times H_{1/2}^1(\Omega_n). $$

For generic $g(y)$, we use the notation

$$ s_g(x, y) := g(y) S(x). $$

As noticed in the 1D example, when the singular coefficient $g^+$ does not vanish, $s_{g^+}$ does not belong to the space $H_{1/2}^1(\Omega_p) \times H_{1/2}^1(\Omega_n)$, hence the notation $s_{g^+}$, with $s$ for “singular”.

Given $\mathbf{u}^+ = (u_p^+, u_n^+) \in H_{1/2}^1(\Omega_p) \times H_{1/2}^1(\Omega_n)$ and $g^+ \in H_{per}^1(\Sigma)$ a solution of the system (2.6), no transmission condition through $\Sigma$ is imposed a priori between $u_p^+$ and $u_n^+$. On the other hand, the convergence assumption (2.5) contains a hidden transmission condition through $\Sigma$, as we will see in section 5.

Let us conclude by the simple result below.
**Proposition 2.3.** Let $u$ be governed by (2.3).

In the framework of Assumption 2.2, there holds that $u \in H^1(\Omega \setminus V_\Sigma)$, for every neighborhood $V_\Sigma$ of the interface.

Conversely, if $u$ can be decomposed as $u = u_{reg} + u_{sing}$ where $u_{reg}|_{\Omega_{p,n}} \in H_{1/2}^1(\Omega_{p,n})$ is periodic in $y$-direction and $u_{sing}(x, y) = g(y)S(x)$, and if there holds that $u \in H^1(\Omega \setminus V_\Sigma)$ for every neighborhood $V_\Sigma$ of the interface, then $g \in H_{per}^1(\Sigma)$.

# 3 The numerical method from Nicolopoulos et al [13]

We now recall the main ingredients that were used by Nicolopoulos et al in [13, 14] to build a numerical approximation to problem (2.3), with its solution split into a regular and a singular part. A stronger assumption on the singular part is used in [13], namely that $g \in H_{per}^2(\Sigma)$. This choice has strong consequences on the variational reformulation of the model, see §3.1, and also on the numerical approximation of the model, see §3.2.

## 3.1 Main ideas and formulation of the method

In order to explain the method of [13], let us introduce the following functions, which we describe as "singularities with absorption"

$$s_g^{\nu}(x, y) := g(y) \log \left( x + \frac{i\nu}{\mathbf{r}(y)} \right) \quad \text{with} \quad \nu > 0, g \in H_{per}^1(\Sigma), \qquad (3.1)$$

and where $\mathbf{r}(y) = \partial_x \alpha(0, y)$. The absorption parameter scaled by $1/\mathbf{r}(y)$ will ensure some nice convergence properties on second order derivatives, see lemma 5.11. We also introduce the weighted $L^2(\Sigma)$-norm $\|g\|_{\mathbf{r}} := (\int_{\Sigma} |g(y)|^2 \mathbf{r}(y) dy)^{1/2}$ and its associated inner product is denoted $(\cdot, \cdot)_{\mathbf{r}}$.

Note that for any $g \in L^2(\Sigma)$, it is easily checked that $s_g^{\nu} \to s_g$ in $L^2(\Omega)$ as $\nu \to 0^+$. We then have the following lemma, whose proof is left to the reader.

**Lemma 3.1.** Given $g \in H_{per}^1(\Sigma)$, the following limits hold in $L^2(\Omega)$ as $\nu \to 0+$:

$$\begin{array}{lll} s_g^{\nu} \to s_g, & \partial_y s_g^{\nu} \to \partial_y s_g, & \\ (\alpha + i\nu) \partial_x s_g^{\nu} \to \alpha \partial_x s_g, & \partial_x ((\alpha + i\nu) \partial_x s_g^{\nu}) \to \partial_x (\alpha \partial_x s_g). & \end{array}$$

Let $\varphi$ be a truncation function satisfying

**Definition 3.2.** Given $\varphi_1 \in C_0^1((-a, a); \mathbb{R})$ and $\varphi_1 = 1$ in the vicinity of $x = 0$, let $\varphi(x, y) = \varphi_1(x)$.

The method of [13] relies on the observation that, for $g \neq 0$, the singular ansatz $s_g$ does not belong to $H_{1/2}^1(\Omega_n) \times H_{1/2}^1(\Omega_p)$. One can check through simple computations that, as soon as $g \in H_{per}^1(\Sigma)$, and in particular for $g \in H_{per}^2(\Sigma)$, it holds that

$$\lim_{\nu \to 0+} \int_{\Omega} \nu |\nabla s_g^{\nu}|^2 \varphi \, d\mathbf{x} = \pi \|g\|_{\mathbf{r}}^2 > 0. \qquad (3.2)$$

Physically, the above identity is related to the plasma heating [14, 7]. The above observation serves as a basis to construct a functional to minimize; this minimization procedure will yield a
variational formulation. Let $u^{\nu} \in H^1(\Omega)$ be the unique solution of (2.4) (see [13], Proposition 6). By assumption 2.2, $u^{\nu} \to u_{reg}^{+} + s_{g^{+}}$ in $L^2(\Omega)$, as $\nu \to 0^+$. We then split $u^{\nu} = u_{reg}^{\nu} + s_{g^{+}}^{\nu}$; evidently, $s_{g^{+}}^{\nu} \to s_{g^{+}}$ in $L^2(\Omega)$, and $u_{reg}^{\nu} \to u_{reg}^{+}$ in $L^2(\Omega)$. Recall that $u_{reg}^{+}$ is identified with a pair of functions $\mathbf{u}^{+} = (u_{p}^{+}, u_{n}^{+})$.

Next we introduce the function $s_h^{\nu}$, with $h \in H_{per}^2(\Sigma)$ being an artificial variable, and define the following energy functional:

$$ \mathcal{J}^{\nu}(u^{\nu}, h) := \operatorname{Im} \left( \int_{\Omega} (\alpha(x, y) + i\nu) |\nabla (u^{\nu} - s_h^{\nu})|^2 \varphi \, d\mathbf{x} \right) = \int_{\Omega} \nu |\nabla (u_{reg}^{\nu} + s_{g^{+}-h}^{\nu})|^2 \varphi \, d\mathbf{x}. \quad (3.3) $$

Clearly, $\mathcal{J}^{\nu}(u^{\nu}, h)$ is non-negative.

To be able to carry out the derivation of the model, one needs some a priori convergence results. One can prove the following result (see lemma A.4, appendix A).

**Proposition 3.3.** Let $(u^{\nu})_{\nu}$ be a family governed by (2.4) fulfilling assumption 2.2. Then,

$$ \lim_{\nu \to 0^+} \int_{\Omega} \nu |\nabla u_{reg}^{\nu}|^2 \varphi \, d\mathbf{x} = 0. \quad (3.4) $$

As a consequence of the above, $\mathcal{J}^{\nu}(u^{\nu}, h)$ converges to the limit $\pi \|g^{+} - h\|_{\mathbf{r}}^2$ when $\nu \to 0^+$. Obviously, $\lim_{\nu \to 0^+} \mathcal{J}^{\nu}(u^{\nu}, g^{+}) = 0$. Hence, the limit of $(u^{\nu}, g^{+}) = (u_{reg}^{\nu}, s_{g^{+}}^{\nu}, g^{+})$ as $\nu \to 0^+$ should solve the limiting minimization problem $\min_{(\mathbf{u}, g^{+}, h)} \mathcal{J}^{+}(\mathbf{u}, g^{+}, h)$, where $\mathcal{J}^{+} := \lim_{\nu \to 0^+} \mathcal{J}^{\nu}(u^{\nu}, h)$, subject to the constraint of $\mathbf{u}, g^{+}$ satisfying assumption 2.2. It remains to write the respective saddle-point problem.

To do so, we follow [13], where it is suggested to use a stronger assumption, namely that $g^{+}$ (and, consequently, $h$) belongs to the space $H_{per}^2(\Sigma)$. Introducing $Q := H_{1/2}^1(\Omega_p) \times H_{1/2}^1(\Omega_n)$ and $V^{(2)} := Q \times H_{per}^2(\Sigma) \times H_{per}^2(\Sigma)$, we find that the limit is governed by the following mixed variational formulation (cf. appendix B for a detailed derivation):

$$ \begin{cases} \text{Find } (\mathbf{u}, g, h) \in V^{(2)} \text{ and } \boldsymbol{\lambda} \in Q \text{ such that} \\ \left\{ \begin{array}{l} a^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) - \overline{b^{(2)}((\mathbf{v}, k, l), \boldsymbol{\lambda})} = 0, \quad \forall (\mathbf{v}, k, l) \in V^{(2)}, \\ b^{(2)}((\mathbf{u}, g, h), \boldsymbol{\mu}) = \ell^{(2)}(\boldsymbol{\mu}), \quad \forall \boldsymbol{\mu} \in Q. \end{array} \right. \end{cases} \quad (3.5) $$

First, the form $a^{(2)}: V^{(2)} \times V^{(2)} \to \mathbb{C}$ is s.t. Im $a^{(2)} = d\mathcal{J}^{+}$, and can be recast as

$$ \begin{aligned} a^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) := & \sum_{j \in \{p, n\}} \int_{\Omega_j} \alpha(u_j + s_{g-h}) \overline{\partial_x (v_j + s_{k-l})} \partial_x \varphi \, d\mathbf{x} \\ & - \int_{\Omega_j} \alpha \overline{(v_j + s_{k-l})} \partial_x (u_j + s_{g-h}) \partial_x \varphi \, d\mathbf{x} \\ & - \int_{\Omega_j} (-\operatorname{div} \alpha \nabla s_h - \omega^2 s_h) \overline{(v_j + s_{k-l})} \varphi \, d\mathbf{x} \\ & + \int_{\Omega_j} (-\operatorname{div} \alpha \nabla s_l - \omega^2 s_l) (u_j + s_{g-h}) \varphi \, d\mathbf{x}. \end{aligned} \quad (3.6) $$

The sesquilinear form $b^{(2)}: V^{(2)} \times Q \to \mathbb{C}$ is, in its turn, given by

$$ b^{(2)}((\mathbf{u}, g, h), \mathbf{v}) = b_{reg}^{(2)}(\mathbf{u}, \mathbf{v}) + b_{sing}^{(2)}(g, \mathbf{v}), \quad (3.7) $$
where, for all $\mathbf{u}, \mathbf{v} \in Q, g \in H_{per}^2(\Sigma)$

$$b_{reg}^{(2)}(\mathbf{u}, \mathbf{v}) := \sum_{j \in \{p,n\}} \int_{\Omega_j} (\alpha \nabla u_j \cdot \overline{\nabla v_j} - \omega^2 u_j \overline{v_j}) d\mathbf{x} + \int_{\Gamma_j} i\lambda u_j \overline{v_j} ds, \quad (3.8)$$

$$b_{sing}^{(2)}(g, \mathbf{v}) := \sum_{j \in \{p,n\}} \int_{\Omega_j} (-\operatorname{div}(\alpha \nabla s_g) - \omega^2 s_g) \overline{v_j} d\mathbf{x} + \int_{\Gamma_j} (\alpha \partial_n s_g + i\lambda s_g) \overline{v_j} ds. \quad (3.9)$$

Remark that neither $b_{reg}^{(2)}$ nor $b_{sing}^{(2)}$ do not depend on $h$. Finally, the antilinear form $\ell^{(2)}: Q \to \mathbb{C}$ is defined as

$$\ell^{(2)}(\boldsymbol{\mu}) = \sum_{j \in \{p,n\}} \int_{\Gamma_j} f \overline{\boldsymbol{\mu}} \, ds.$$

To stabilize the method in order to guarantee uniqueness, two stabilization terms are added to (3.5). Namely, one considers

$$\text{Find } (\mathbf{u}, g, h) \in V^{(2)} \text{ and } \boldsymbol{\lambda} \in Q \text{ such that} \\ \left\{ \begin{array}{l} a_{\rho, \mu}^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) - \overline{b^{(2)}((\mathbf{v}, k, l), \boldsymbol{\lambda})} = 0, \quad \forall (\mathbf{v}, k, l) \in V^{(2)}, \\ b^{(2)}((\mathbf{u}, g, h), \boldsymbol{\mu}) = \ell^{(2)}(\boldsymbol{\mu}), \quad \forall \boldsymbol{\mu} \in Q, \end{array} \right. \quad (3.10)$$

where

$$a_{\rho, \mu}^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = a^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) + i(-\rho(g, k)H^2(\Sigma) + \mu(\partial_{yy}h, \partial_{yy}l)_{L^2(\Sigma)}),$$

with $\rho, \mu > 0$. It is shown in [13, Theorem 16] that for $\rho, \mu > 0$, and $f \in L^2(\Gamma_p \cup \Gamma_n)$, the problem (3.10) is well-posed. Up to our knowledge, there exists no proof that the solution to (3.10) is a limiting absorption solution of the original problem.

## 3.2 Numerical experiments and comments

In [13], a conforming discretization of (3.10) was proposed, with $V_{h_1, h_2}^{(2)} = Q_{h_1} \times H_{h_2}^2 \times H_{h_2}^2$,

$$Q_{h_1} = \{v_{h_1} \in Q : v_{h_1}|_K \in P_1(K), \text{ for all } K \in \mathcal{T}_{h_1}^{\Omega}\},$$

$$H_{h_2}^2 = \{p_{h_2} \in H_{per}^2(\Sigma) : p_{h_2}|_K \in H_m(K), \text{ for all } K \in \mathcal{T}_{h_2}^{\Sigma}\},$$

where $H_m(K)$ is Hermite finite element of order $m$, $\mathcal{T}_{h_1}^{\Omega}$ is a triangulation of $\Omega$ with meshsize $h_1$ that is conforming with respect to the interface $\Sigma$ (for all $K \in \mathcal{T}_{h_1}$, $\mathrm{int}(K) \cap \Sigma = \emptyset$), and $\mathcal{T}_{h_2}^{\Sigma}$ is a triangulation of $\Sigma$ with meshsize $h_2$. Notice that the restriction to $\Omega_p$ (respectively $\Omega_n$) of elements of $Q_{h_1}$ belongs to $H^1(\Omega_p)$ (resp. $H^1(\Omega_n)$). On the other hand, there is no matching condition at the interface for elements of $Q_{h_1}$.

In the original paper [13], the numerical experiments were done for a single discretization. Structured meshes were used for both regular and singular parts, with $h_2 = 4h_1$. In particular, the question of the convergence of the discrete solution to the continuous one was not addressed. The goal of this section is to provide insight into this question, by letting $h_1$ vary and keeping $h_2 = 4h_1$.

We use the code provided by A. Nicolopoulos written in FreeFem++ [10]. The singular coefficient $g^+$ is discretized with the 2D HCT finite elements penalized along the normal direction.
We consider the case $\alpha(x, y) = x, \omega = 0$ and perform two experiments with $L = 2$ on the domain $\Omega = (-1, 1) \times (-1, 1)$ with known exact solutions. We choose the boundary data so that, in the first case, the exact solution is purely regular and equal to $u^+(x, y) = 1$, and in the second case it is given by $u^+(x, y) = -K_0(\pi x)e^{i\pi y}$, where $K_0$ is a modified Bessel function of the second kind $^1$. In the first case $g^+(y) = 0$, while in the second case, there is a non-zero singular part $s_{g^+}(x, y)$ with $g^+(y) = e^{i\pi y}$.

The stabilization parameters $\rho_2, \mu_2$ are taken equal.

We denote by $e_{L^2}(\mathbf{u}^+)$, resp. $e_Q(\mathbf{u}^+)$, the relative error of the regular part in $L^2(\Omega)$-norm, resp. in $\|\cdot\|_Q$ norm. And we denote by $e_{L^2}(g^+)$ the relative error of the singular coefficient in $L^2(\Sigma)$-norm. Note that, when measuring volume errors, we do not take into account the cells that touch the interface.

Although $u^+(x, y) = 1$ belongs to the discrete space, the computed solution does not seem to converge, see figure 2. This phenomenon happens regardless of the value of the stabilization parameters $\rho_2, \mu_2$. The situation improves somewhat for the singular solution, where one observes

<figure><img src="image_3.png" /><figcaption>Figure 2: Relative errors for $u^+(x, y) = 1$, with $\rho_2 = \mu_2 = 10^{-5}$.</figcaption></figure>

a monotonic decrease of the relative error in $L^2(\Sigma)$-norm for the singular coefficient $g^+$, see figure 3a. However, convergence for the regular part is not obvious in $\|\cdot\|_Q$ norm nor in $L^2(\Omega)$-norm, see again figure 3a. In figure 3b, we provide convergence curves depending on the choice of the parameters $\rho_2 = \mu_2$. As expected, the convergence seems to stagnate for larger values of $\rho_2$, and finally we see that decreasing $\rho_2$ from $10^{-5}$ to $10^{-6}$ has no visible effect on the convergence curves.

These experiments seem to indicate that the numerical method of [13] does not converge numerically. We do not know whether the source of the instability is intrinsic to the numerical variational formulation itself, or is due to the penalization of the HCT elements in the normal direction, used in the implementation. In this manuscript, we will not dwell on the precise reason for this instability. Instead, we propose in section 4 an alternative method, which we study mathematically in section 5. Interestingly, this method can be discretized in a conforming manner without using Hermite elements on the interface: it relies on the more common $P_1$ finite elements.

<footnote n="1">¹For $x < 0$, we use the convention $K_0(x) = \lim_{\nu \to 0+} K_0(x + i\nu)$.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Relative errors for $u^+(x, y) = -K_0(\pi x)e^{i\pi y}$.</figcaption></figure>

## 4 An alternative method

One of the drawbacks of the formulation (3.10) is that it requires the extra regularity on $g$ ($g \in H_{per}^2(\Sigma)$), which renders the underlying numerical method quite cumbersome to implement. Also, a priori, there is no reason to assume such a regularity, especially from the variational viewpoint. As a matter of fact, this regularity assumption can be weakened. To see that, we observe that the form $a^{(2)}$ can be rewritten by integration by parts of the last two integrals along the $y$-direction (recall that $\varphi$ is independent of $y$):

$$ 
\begin{aligned}
a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) & \\
&:= \sum_{j \in \{p, n\}} \int_{\Omega_j} \alpha \left[ (u_j + s_{g-h}) \overline{\partial_x (v_j + s_{k-l})} - \partial_x (u_j + s_{g-h}) \overline{(v_j + s_{k-l})} \right] \partial_x \varphi \, d\mathbf{x} \\
&\quad - \int_{\Omega_j} \left( \alpha \partial_y s_h \overline{\partial_y (v_j + s_{k-l})} + (-\partial_x (\alpha \partial_x s_h) - \omega^2 s_h) \overline{(v_j + s_{k-l})} \right) \varphi \, d\mathbf{x} \\
&\quad + \int_{\Omega_j} \left( \overline{\alpha \partial_y s_l} \partial_y (u_j + s_{g-h}) + (-\partial_x (\alpha \partial_x s_l) - \omega^2 s_l) \overline{(u_j + s_{g-h})} \right) \varphi \, d\mathbf{x}.
\end{aligned} 
\quad (4.1) $$

The periodic boundary condition on $\partial_y g$ is now a natural boundary condition.

Remark that the expression $\int_{\Omega_j} \alpha \partial_y s_h \overline{\partial_y (v_j + s_{k-l})} d\mathbf{x}$ is well-defined for all $h, k, l \in H_{per}^1(\Sigma)$ and $v_j \in H_{1/2}^1(\Omega_j)$. Indeed,

$$ 
\begin{aligned}
\left| \int_{\Omega_j} \alpha (\log |x| + i\pi \mathbb{1}_{x<0}) \partial_y h \overline{\partial_y v_j} d\mathbf{x} \right| &< +\infty, \\
\left| \int_{\Omega_j} \alpha |\log |x| + i\pi \mathbb{1}_{x<0}|^2 \partial_y h \overline{\partial_y (k-l)} d\mathbf{x} \right| &< +\infty.
\end{aligned} 
$$

Thus, we consider the larger function space $V^{(1)} \times Q$ with $V^{(1)} := Q \times H_{per}^1(\Sigma) \times H_{per}^1(\Sigma)$. The form $a^{(1)}$ is naturally continuous and sesquilinear on $V^{(1)} \times V^{(1)}$.
**Remark 4.1.** The form $a^{(1)}$ can also be obtained as the differential of a limit functional, like $a^{(2)}$, by repeating almost verbatim the derivation in appendix B, and the limit functional writes:

$$ \mathcal{J}^{+}(\mathbf{u}, g, h) := \frac{1}{2i} a^{(1)}((\mathbf{u}, g, h), (\mathbf{u}, g, h)). \qquad (4.2) $$

While for $g, h \in H_{per}^2(\Sigma)$ one has also:

$$ \mathcal{J}^{+}(\mathbf{u}, g, h) = \frac{1}{2i} a^{(2)}((\mathbf{u}, g, h), (\mathbf{u}, g, h)). $$

**Remark 4.2.** It does not seem possible to further reduce the regularity of $g, h, k, l$ because their first derivative appears in (4.1).

The same integration by parts along the $y$-direction can be applied to the form $b_{sing}^{(2)}$: we define, for all $\mathbf{v} \in Q, g \in H_{per}^1(\Sigma)$,

$$ b_{sing}^{(1)}(g, \mathbf{v}) := \sum_{j \in \{p, n\}} \int_{\Omega_j} (\alpha \partial_y s_g \overline{\partial_y v_j} + (-\partial_x(\alpha \partial_x s_g) - \omega^2 s_g) \overline{v_j}) d\mathbf{x} + \int_{\Gamma_j} (\alpha \partial_n s_g + i\lambda s_g) \overline{v_j} ds. \quad (4.3) $$

We now define $b^{(1)}: V^{(1)} \times Q \to \mathbb{C}$, which is, in its turn, given by

$$ b^{(1)}((\mathbf{u}, g, h), \mathbf{v}) := b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{sing}^{(1)}(g, \mathbf{v}), \qquad (4.4) $$

with

$$ b_{reg}^{(1)} = b_{reg}^{(2)}: Q \times Q \to \mathbb{C}. \qquad (4.5) $$

Finally, we introduce $\ell^{(1)} = \ell^{(2)}: Q \to \mathbb{C}$.

**Remark 4.3.** Compared to the definition of $b_{sing}^{(2)}$ (3.9), only the first order derivatives of $g$ appear. Hence, if $g \in H_{per}^1(\Sigma)$, then $b_{sing}^{(1)}(g, \cdot)$ defines a continuous antilinear form on $Q$, and, moreover, $\|b_{sing}^{(1)}(g, \cdot)\|_{Q'} \lesssim \|g\|_{H^1(\Sigma)}$ (to be compared with $\|b_{sing}^{(2)}(g, \cdot)\|_{Q'} \lesssim \|g\|_{H^2(\Sigma)}$). While for $g \notin H^1(\Sigma)$, $b_{sing}^{(1)}(g, \cdot)$ cannot be defined. On the other hand, the regular parts of $b^{(2)}$ and $b^{(1)}$ are identical, and so are the right-hand sides.

Hence, we end up with the following variational formulation:

$$ \begin{array}{ll} \text{Find } ((\mathbf{u}, g, h), \boldsymbol{\lambda}) \in V^{(1)} \times Q \text{ such that} & \\ \left\{ \begin{array}{ll} a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) - \overline{b^{(1)}((\mathbf{v}, k, l), \boldsymbol{\lambda})} = 0, & \forall (\mathbf{v}, k, l) \in V^{(1)}, \\ b^{(1)}((\mathbf{u}, g, h), \boldsymbol{\mu}) = \ell^{(1)}(\boldsymbol{\mu}), & \forall \boldsymbol{\mu} \in Q. \end{array} \right. & (4.6) \end{array} $$

With this method, a stabilized counterpart of (4.6) reads:

$$ \begin{array}{ll} \text{Find } ((\mathbf{u}, g, h), \boldsymbol{\lambda}) \in V^{(1)} \times Q \text{ such that} & \\ \left\{ \begin{array}{ll} a_{\rho}^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) - \overline{b^{(1)}((\mathbf{v}, k, l), \boldsymbol{\lambda})} = 0, & \forall (\mathbf{v}, k, l) \in V^{(1)}, \\ b^{(1)}((\mathbf{u}, g, h), \boldsymbol{\mu}) = \ell^{(1)}(\boldsymbol{\mu}), & \forall \boldsymbol{\mu} \in Q, \end{array} \right. & (4.7) \end{array} $$
where

$$a_{\rho}^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) - i\rho((g, k)_{H^1(\Sigma)} - (h, l)_{H^1(\Sigma)}), \quad \text{with } \rho > 0.$$

The stabilization terms involve the $H_{per}^1(\Sigma)$ scalar product, which corresponds to the fact that $g$ does belong to $H_{per}^1(\Sigma)$.

Retracing the steps of [13], we can prove the well-posedness result below regarding the stabilized variational formulation. We use a classical approach to the well-posedness of the mixed formulations. According to the Babuška-Brezzi theory, it is sufficient to prove a surjectivity property of the operator $\mathbf{B}^{(1)}: V^{(1)} \to Q'$ associated to the form $b^{(1)}$ and an inf-sup condition for the form $a_{\rho}^{(1)}$ on the kernel of $\mathbf{B}^{(1)}$. The kernel of the operator $\mathbf{B}^{(1)}$ can be characterized as

$$\mathrm{Ker} \mathbf{B}^{(1)} = \{(\mathbf{u}, g, h) \in V^{(1)} : b_{\mathrm{reg}}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{\mathrm{sing}}^{(1)}(g, \mathbf{v}) = 0, \text{ for all } \mathbf{v} \in Q\}. \quad (4.8)$$

We observe that, since the third variable $h$ does not appear in the characterization of the kernel, it can take any value. We use this property on several occasions throughout the manuscript. To study the surjectivity of the operator $\mathbf{B}^{(1)}$, we introduce the following problem: given $\ell \in Q'$,

$$\begin{gathered} \text{Find } \mathbf{v} \in Q \text{ s. t.} \\ b_{\mathrm{reg}}^{(1)}(\mathbf{v}, \boldsymbol{\mu}) = \ell(\boldsymbol{\mu}), \quad \forall \boldsymbol{\mu} \in Q. \end{gathered} \quad (4.9)$$

Problem 4.9) is well-posed. The result can be proven like in proposition 13 of [13]. As a straightforward consequence, one finds that

**Proposition 4.4.** The operator $\mathbf{B}^{(1)}$ is surjective.

We can now state the well-posedness result of the stabilized variational formulation.

**Theorem 4.5.** Let $\rho > 0$. For all $f \in L^2(\Gamma_p \cup \Gamma_n)$, the stabilized mixed formulation (4.7) admits a unique solution.

*Proof.* One needs to verify an inf-sup condition for the form $a_{\rho}^{(1)}$ on the kernel of $\mathbf{B}^{(1)}$. Again, the proof mimics the one in [13], and is based on the following two equalities for $(\mathbf{u}, g, h) \in \mathrm{Ker} \mathbf{B}^{(1)}$:

$$\begin{align} a^{(1)}((\mathbf{u}, g, 0), (\mathbf{u}, g, 0)) &= -2i\lambda \sum_{j \in \{p, n\}} \|u_j + s_g\|_{L^2(\Gamma_j)}^2, \tag{4.10} \\ a^{(1)}((\mathbf{0}, 0, h), (\mathbf{0}, 0, h)) &= 2i\pi \|h\|_{\mathbf{r}}^2. \notag \end{align}$$

For completeness, we give the proof of these equalities in corollary D.5.

Regarding the inf-sup condition, we note that for any $(\mathbf{u}, g, h) \in \mathrm{Ker} \mathbf{B}^{(1)}$,

$$\begin{align} \mathrm{Im} a_{\rho}^{(1)}((\mathbf{u}, g, h), (-\mathbf{u}, -g, h)) &= 2\pi \|h\|_{\mathbf{r}}^2 + 2\lambda \sum_{j \in \{p, n\}} \|u_j + s_g\|_{L^2(\Gamma_j)}^2 \tag{4.11} \\ &\quad + \rho \|g\|_{H^1(\Sigma)}^2 + \rho \|h\|_{H^1(\Sigma)}^2 \notag \\ &\ge \tilde{C} \left( \|\mathbf{u}\|_Q^2 + \|g\|_{H^1(\Sigma)}^2 + \|h\|_{H^1(\Sigma)}^2 \right), \tag{4.12} \end{align}$$

since the norm of $\mathbf{u}$ is controlled by the norm of $g$. In the statement of proposition 13 (with $f = 0$) in [13], one shows that $\|\mathbf{u}\|_Q \le \|b_{\mathrm{sing}}^{(2)}(g, \cdot)\|_{Q'} \lesssim \|g\|_{H^2(\Sigma)}$, while in our case the same argument yields a stability bound $\|\mathbf{u}\|_Q \lesssim \|b_{\mathrm{sing}}^{(1)}(g, \cdot)\|_{Q'} \lesssim \|g\|_{H^1(\Sigma)}$, according to the discussion after (4.3). $\square$
# 5 Existence and uniqueness of the solution of the mixed variational formulation (4.6)

We now study the mixed problem (4.6) (without stabilization) in more details. Namely, we address the uniqueness and existence of its solution. First, one can prove that the solution to the mixed problem (4.6) is unique. Moreover, its existence can be ensured through a direct construction: indeed, under assumption 2.2, the limiting absorption solution satisfies (4.6) with an appropriate choice of the Lagrange multipliers $h, \lambda$.

In order to prove these results, we introduce a notion of weak jump of the regular part of the solution (section 5.1). This notion will serve to prove, on one hand, the uniqueness of the solution to the mixed formulation (4.6) (section 5.2) and, on the other hand, its consistency with the original limiting absorption problem (2.4) (section 5.3).

## 5.1 Weak jump of a regular part

It is possible to define a notion of a jump in a weak sense for functions $\mathbf{u} \in Q$ which satisfy the constraint in the mixed formulation (4.6). This definition is formalized in the statement of the following lemma.

For the sake of conciseness, $\langle \cdot, \cdot \rangle_{\Sigma}$ denotes $\langle \cdot, \cdot \rangle_{(H_{per}^{1}(\Sigma,r))'}$, $H_{per}^{1}(\Sigma,r)$ until the end of the manuscript.

**Lemma 5.1.** Let $\mathbf{u} \in Q$, $g \in H_{per}^{1}(\Sigma)$ and $\ell \in Q'$, $\text{supp} \ell \cap \Sigma = \emptyset^{2}$ related by

$$b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{sing}^{(1)}(g, \mathbf{v}) = \ell(\mathbf{v}), \quad \forall \mathbf{v} \in Q. \qquad (5.1)$$

Let $\varphi$ be a cutoff function as in definition 3.2, and $\varphi_{\varepsilon}(x) = \varphi(\frac{x}{\varepsilon})$ with $\varepsilon > 0$. The jump $[\mathbf{u}]_{\Sigma}$ of the regular part is defined as

$$\forall k \in H_{per}^{1}(\Sigma), \quad \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma} := \lim_{\varepsilon \to 0} \sum_{j \in \{p,n\}} \int_{\Omega_j} u_j \overline{(-\partial_x (\alpha \partial_x (s_k \varphi_{\varepsilon})))} d\mathbf{x}, \qquad (5.2)$$

and it is finite for all $k$. The limit is independent of the choice of $\varphi$.

Before proving the lemma, let us make a few comments. For piecewise regular $\mathbf{u} \in H^{1}(\Omega_{p}) \times H^{1}(\Omega_{n})$, the above definition of the jump coincides with the classical definition $[\mathbf{u}]_{\Sigma} = \gamma_{0}(u_{p}) - \gamma_{0}(u_{n})$, $\gamma_{0}$ denotes the trace on $\Sigma$, seen either as a part of $\Omega_{p}$, or of $\Omega_{n}$. Indeed, in this case (5.2) yields

$$\sum_{j \in \{p,n\}} \int_{\Omega_j} u_j \overline{(-\partial_x (\alpha \partial_x (s_k \varphi_{\varepsilon})))} d\mathbf{x} = \int_{\Sigma} [\gamma_0(u_p) - \gamma_0(u_n)] \overline{(\alpha \partial_x (s_k \varphi_{\varepsilon}))} |_{\Sigma} dy + \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha \partial_x u_j \overline{\partial_x (s_k \varphi_{\varepsilon})} d\mathbf{x}.$$

The first term in the right-hand side of the above can be made more explicit. Indeed, by the regularity assumption on $\alpha$ and using an explicit form of $s_k$, we have $\alpha \partial_x (s_k \varphi_{\varepsilon}) \in H^1(\Omega)$ and

<footnote n="2">²By writing $\text{supp} \ell \cap \Sigma = \emptyset$, we mean that there exists $\delta > 0$ and a $\delta$-vicinity $\Omega^{\delta}$ of $\Sigma$, $\Omega^{\delta} := \{(x, y) \in \Omega : |x| < \delta\}$, s.t. $\ell(v) = 0$ for all $v \in H_{1/2}^{1}(\Omega) : \text{supp } v \subset \Omega^{\delta}$.</footnote>
$\alpha \partial_x(s_k \varphi_\varepsilon)|_{\Sigma} = k(y)\mathbf{r}(y)$. Observing that $\|\alpha \partial_x(s_k \varphi_\varepsilon)\|_{L^2(\Omega_j)} \lesssim \|k\|_{L^2(\Sigma)}$, the second term is bounded with the help of Cauchy-Schwarz inequality:

$$ \left| \int_{\Omega_j} \alpha \partial_x u_j \overline{\partial_x(s_k \varphi_\varepsilon)} d\mathbf{x} \right| \le \|\mathbf{u}\|_{H^1(\Omega_j \cap \{|x|<a\varepsilon\})} \|\alpha \partial_x(s_k \varphi_\varepsilon)\|_{L^2(\Omega_j)} \\ \lesssim \|\mathbf{u}\|_{H^1(\Omega_j \cap \{|x|<a\varepsilon\})} \|k\|_{L^2(\Sigma)} \xrightarrow{\varepsilon \to 0} 0, $$

so that

$$ \sum_{j \in \{p,n\}} \int_{\Omega_j} u_j \overline{(-\partial_x (\alpha \partial_x (s_k \varphi_\varepsilon)))} d\mathbf{x} = \int_{\Sigma} [\mathbf{u}]_{\Sigma} \bar{k} r dy + o_{\varepsilon \to 0}(1). $$

On the other hand, the jump defined by (5.2) is not finite for all $\mathbf{u} \in Q$; take e.g., $\mathbf{u} = (\log |\log |x||, 0)$. Finally, we note that, even for $\mathbf{u}$ satisfying (5.1), the jump $[\mathbf{u}]_\Sigma$ is defined in a very weak sense, since it is taken in the dual space of $H_{per}^1(\Sigma, r)$.

Before proceeding to prove lemma 5.1, let us define, for $\ell \in Q'$, s.t. $\text{supp } \ell \cap \Sigma = \emptyset$, any $\psi \in C^1(\bar{\Omega})$ and $k \in H_{per}^1(\Sigma)$, the quantity

$$ \ell_{\infty}(s_k \psi) := \lim_{\varepsilon \to 0} \ell(s_k \psi(1 - \varphi_{\varepsilon})), \qquad (5.3) $$

where the sequence $(\varphi_\varepsilon)_{\varepsilon>0}$ is defined in the statement of lemma 5.1. The above quantity is well-defined and independent of the choice of $\varphi$, since $\text{supp } \ell \cap \Sigma = \emptyset$.

*Proof of lemma 5.1.* We test (5.1) with $\mathbf{v} = s_k(1 - \varphi_\varepsilon) \in Q$. On one hand, we have

$$ b_{reg}^{(1)}(\mathbf{u}, s_k(1-\varphi_\varepsilon)) \xrightarrow{\varepsilon \to 0} \ell_\infty(s_k) - b_{sing}^{(1)}(g, s_k), \qquad (5.4) $$

where $\ell_\infty(s_k)$ is defined like in (5.3), with $\psi \equiv 1$. On the other hand, integrating by parts in the $x$-direction $b_{reg}^{(1)}(\mathbf{u}, s_k(1-\varphi_\varepsilon))$ gives:

$$ \begin{aligned} b_{reg}^{(1)}(\mathbf{u}, s_k(1-\varphi_\varepsilon)) &= \sum_{j \in \{p,n\}} \int_{\Omega_j} \left\{ u_j \overline{(-\partial_x (\alpha \partial_x (s_k(1-\varphi_\varepsilon))))} + \alpha \partial_y u_j \overline{\partial_y (s_k(1-\varphi_\varepsilon))} - \omega^2 u_j \overline{s_k(1-\varphi_\varepsilon)} \right\} d\mathbf{x} \\ &\quad + \sum_{j \in \{p,n\}} \int_{\Gamma_j} \left( u_j \overline{(\alpha \partial_n s_k)} + i\lambda u_j \overline{s_k} \right) ds \\ &= \overline{b_{sing}^{(1)}(k, \mathbf{u})} + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} ds - \sum_{j \in \{p,n\}} \int_{\Omega_j} \left\{ \alpha \partial_y u_j \overline{\partial_y s_k} \varphi_\varepsilon - \omega^2 u_j \overline{s_k} \varphi_\varepsilon \right\} d\mathbf{x} - J_\varepsilon, \end{aligned} $$

with

$$ J_\varepsilon = \sum_{j \in \{p,n\}} \int_{\Omega_j} u_j \overline{(-\partial_x \alpha \partial_x (s_k \varphi_\varepsilon))} d\mathbf{x}, $$

cf. the definition of the jump (5.2). Next, using Lebesgue's dominated convergence theorem yields

$$ \lim_{\varepsilon \to 0} b_{reg}^{(1)}(\mathbf{u}, s_k(1-\varphi_\varepsilon)) = \overline{b_{sing}^{(1)}(k, \mathbf{u})} + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} ds - \lim_{\varepsilon \to 0} J_\varepsilon. $$
Replacing the left-hand side of the above by (5.4) shows that $\lim_{\epsilon \to 0} J_\epsilon$ is finite and, with the definition (5.2), the jump of $\mathbf{u}$ is expressed as

$$ \langle [\mathbf{u}]_{\Sigma}, k \rangle_{(H_{\text{per}}^{1}(\Sigma, r))', H_{\text{per}}^{1}(\Sigma, r)} = \overline{b_{\text{sing}}^{(1)}(k, \mathbf{u})} + b_{\text{sing}}^{(1)}(g, s_k) + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} ds - \ell_{\infty}(s_k). \quad (5.5) $$

As claimed, the last expression does not depend on the chosen cutoff function $\varphi$. $\square$

**Remark 5.2.** We can relax the condition on the support of $\ell$ by requiring only that $\lim_{\epsilon \to 0} \ell(s_k(1 - \varphi_\epsilon))$ exists and is independent of $(\varphi_\epsilon)_\epsilon$. In particular, if it holds that $b_{\text{reg}}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{\text{sing}}^{(1)}(g, \mathbf{v}) = \ell(\mathbf{v})$ for any $\mathbf{v} \in Q$ with

$$ \ell(\mathbf{v}) = \sum_{j \in \{p,n\}} \int_{\Omega_j} f_j \overline{v_j} d\mathbf{x} \quad \text{where} \quad f_j \in L^2(\Omega_j), \quad j \in \{p,n\}, $$

then the trace can be defined as in (5.5).

**Remark 5.3.** The definition of the jump (5.2) does not depend on the jump part of the singularity $i\pi k(y)\mathbb{1}_{x<0}$. Indeed, in (5.2), $s_k$ can be replaced by $k(y) \log|x|$ or $k(y)(\log|x| - i\pi\mathbb{1}_{x<0})$. This holds because, given $\varphi$ as in lemma 5.1, we have, by the Cauchy-Schwarz inequality and after integration by parts in the $x$-direction,

$$ \begin{aligned} \left| \int_{\Omega} u \partial_x (\alpha \partial_x (\mathbb{1}_{x<0} \varphi_{\varepsilon} k(y))) d\mathbf{x} \right| &= \left| \int_{\Omega_n} u_n \partial_x (\alpha \partial_x (\varphi_{\varepsilon} k(y))) d\mathbf{x} \right| \\ &= \left| \int_{-a\varepsilon < x < 0} \alpha \partial_x u_n \partial_x \varphi_{\varepsilon} k(y) d\mathbf{x} \right| \lesssim \|u_n\|_{H_{1/2}^{1}(\{|x|<a\varepsilon\})} \|k\|_{L^2(\Sigma)} \xrightarrow{\varepsilon \to 0} 0. \end{aligned} $$

Therefore, the jump can also be computed as

$$ \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma} = \lim_{\varepsilon \to 0} \sum_{j \in \{p,n\}} \int_{\Omega_j} u_j \overline{(-\partial_x (\alpha \partial_x (k(y) (\log|x| - i\pi \mathbb{1}_{x<0}) \varphi_\varepsilon)))} d\mathbf{x}. \quad (5.6) $$

This identity will be useful later, see section 5.3.2.

## 5.2 Uniqueness of the solution

The goal of this section is to prove the following result.

**Theorem 5.4.** The solution to (4.6), if it exists, is unique.

**Remark 5.5.** Let us remark that the uniqueness result does not follow directly from the identities (4.10). Indeed, applying these identities allows to conclude that $h = 0$, and $(u_j + s_g)|_{\Gamma_j} = 0$; the latter, however, does not imply that $\mathbf{u} = 0$ and $g = 0$.

The proof of this theorem hinges on some technical results, which will additionally allow us to understand some properties of solutions to (4.6). First, given a solution $((\mathbf{u}, g, h), \mathbf{\lambda}) \in V^{(1)} \times Q$ of (4.6), the condition (5.1) is obviously satisfied for $\mathbf{u}, g$ with $\ell = \ell^{(1)}$, and consequently the jump $[\mathbf{u}]_{\Sigma}$ is well-defined. Therefore, we can reexpress $a^{(1)}$ using the jump $[\mathbf{u}]_{\Sigma}$.
**Lemma 5.6.** Let $(\mathbf{u}, g, h) \in V^{(1)}$ be such that $\ell: \mathbf{v} \mapsto b^{(1)}((\mathbf{u}, g, h), \mathbf{v}) \in Q'$ with $\text{supp} \ell \cap \Sigma = \emptyset$. Let $\varphi$ satisfy definition 3.2 and $\text{supp} \ell \cap \text{supp} \varphi = \emptyset$. For any $(\mathbf{v}, k, l) \in \text{Ker } \mathbf{B}^{(1)}$, we have

$$a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = 2i\pi (g - h, k - l)_{\mathbf{r}} - \langle [\mathbf{u}]_{\Sigma}, k - l \rangle_{\Sigma} + \overline{\langle [\mathbf{v}]_{\Sigma}, g - h \rangle_{\Sigma}}. \quad (5.7)$$

A proof of lemma 5.6 can be found in appendix D, see Lemma D.1.

**Remark 5.7.** In the context of lemma 5.6 and in the light of remark 4.1, one can rewrite the minimization functional $\mathcal{J}^{+}(\mathbf{u}, g, h) = \frac{1}{2i}a^{(1)}((\mathbf{u}, g, h), (\mathbf{u}, g, h))$ on the kernel of $\mathbf{B}^{(1)}$:

$$\mathcal{J}^{+}(\mathbf{u}, g, h) = \pi \|g - h\|_{\mathbf{r}}^{2} - \text{Im} \langle [\mathbf{u}]_{\Sigma}, g - h \rangle_{\Sigma}, \quad \forall (\mathbf{u}, g, h) \in \text{Ker } \mathbf{B}^{(1)}.$$

We now have the necessary ingredients to prove theorem 5.4.

**Proof of theorem 5.4.** Let $((\mathbf{u}, g, h), \boldsymbol{\lambda}) \in V^{(1)} \times Q$ be such that

$$\left\{ \begin{array}{l} a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) - \overline{b^{(1)}((\mathbf{v}, k, l), \boldsymbol{\lambda})} = 0, \quad \forall (\mathbf{v}, k, l) \in V^{(1)}, \\ b^{(1)}((\mathbf{u}, g, h), \boldsymbol{\mu}) = 0. \quad \forall \boldsymbol{\mu} \in Q. \end{array} \right.$$

First, we recall that $(\mathbf{0}, 0, h) \in \text{Ker } \mathbf{B}^{(1)}$, cf. (4.8). Using the identities (4.10) yields $h = 0$. So, one has in particular

$$a^{(1)}((\mathbf{u}, g, 0), (\mathbf{v}, k, l)) = 0, \quad \forall (\mathbf{v}, k, l) \in \text{Ker } \mathbf{B}^{(1)}. \quad (5.8)$$

We now prove that $(\mathbf{u}, g) = (\mathbf{0}, 0)$ by choosing ad hoc test functions in $\text{Ker } \mathbf{B}^{(1)}$. As a matter of fact, since $(\mathbf{0}, 0, g), (\mathbf{u}, g, 0) \in \text{Ker } \mathbf{B}^{(1)}$, from (5.7) it follows that

$$\begin{aligned} a^{(1)}((\mathbf{u}, g, 0), (0, 0, g)) &= -2i\pi \|g\|_{\mathbf{r}}^{2} + \langle [\mathbf{u}]_{\Sigma}, g \rangle_{\Sigma}, \\ a^{(1)}((\mathbf{u}, g, 0), (\mathbf{u}, g, 0)) &= 2i\pi \|g\|_{\mathbf{r}}^{2} - 2i \text{Im} \langle [\mathbf{u}]_{\Sigma}, g \rangle_{\Sigma}, \end{aligned}$$

so that

$$\text{Im } a^{(1)}((\mathbf{u}, g, 0), (\mathbf{u}, g, 2g)) = -2\pi \|g\|_{\mathbf{r}}^{2}.$$

Recalling (5.8), the above implies that $g = 0$.

Finally, $(\mathbf{u}, 0, 0) \in \text{Ker } \mathbf{B}^{(1)}$ implies that $b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) = 0$ for any $\mathbf{v} \in Q$. According to the well-posedness of the problem (4.9), this implies that $\mathbf{u} = 0$.

The above allows to handle the part of the solution that belongs to $V^{(1)}$. To conclude, one has to handle the Lagrange multiplier. According to the above, it is governed by

$$b^{(1)}((\mathbf{v}, k, l), \boldsymbol{\lambda}) = 0, \quad \forall (\mathbf{v}, k, l) \in V^{(1)}.$$

But the operator $\mathbf{B}^{(1)}: V^{(1)} \to Q'$ is surjective thanks to proposition 4.4: one has indeed $\boldsymbol{\lambda} = \mathbf{0}$. $\square$

## 5.3 Further properties and existence of the solution

Below, we denote by $\mathbf{u}$ the pair $(u_p, u_n) \in Q$.
### 5.3.1 Consistency and jump of the regular part

Theorem 5.4 shows that the mixed variational formulation (4.6) has at most one solution. Therefore, if we can construct a solution to this formulation, it will be unique. It is thus reasonable to look for $(\mathbf{u}, g)$ as the limiting absorption solution of the original formulation (2.4), as $\nu \to 0^+$. From the content of section 3.1, we should expect that $h = g$. Moreover, again, by explicit computations, in this case one can show that the Lagrange multiplier $\lambda$ is equal to $\mathbf{u}\varphi$. However, this is not straightforward when comparing the section 3.1 and the remark 5.7, due to the presence of the extra term involving the jump $[\mathbf{u}]_{\Sigma}$ in the functional $\mathcal{J}^+$, which did not seem to occur in the original functional $\mathcal{J}^\nu$. This term becomes more apparent in the following proposition.

**Proposition 5.8.** Let $(\mathbf{u}, g) \in Q \times H_{per}^1(\Sigma)$ be such that $b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{sing}^{(1)}(g, \mathbf{v}) = \ell^{(1)}(\mathbf{v})$ for any $\mathbf{v} \in Q$. Then, for any $(\mathbf{v}, k, l) \in V^{(1)}$,

$$a^{(1)}((\mathbf{u}, g, g), (\mathbf{v}, k, l)) - \overline{b^{(1)}((\mathbf{v}, k, l), \mathbf{u}\varphi)} = -\langle [\mathbf{u}]_{\Sigma}, k-l \rangle_{\Sigma}. \quad (5.9)$$

As a consequence, $((\mathbf{u}, g, g), \mathbf{u}\varphi)$ is the solution of (4.6) if and only if $[\mathbf{u}]_{\Sigma} = 0$.

The above proposition shows that $h = g$ and $\lambda = \mathbf{u}\varphi$ if and only if the jump of the regular part of the solution $\mathbf{u}$ vanishes. Notice that the above proposition does not ensure the existence of a solution to the mixed variational formulation (4.6), nor that $(\mathbf{u}, g, g)$ is the solution, since it may happen that the solution of (4.6) satisfies $[\mathbf{u}]_{\Sigma} \neq 0$.

Nonetheless, the above shows that the question of the consistency of the mixed variational formulation with the original limiting absorption problem reduces to the question of vanishing of the jump of the regular part $[\mathbf{u}]_{\Sigma}$, where we seek $(\mathbf{u}, g)$ to be the limiting absorption solution. Precisely, it will be proved that its jump vanishes in proposition 5.13.

The proof of proposition 5.8 again relies on two auxiliary lemmas, which we present below and which are proven in appendix D. First, let us introduce the form

$$C_{\varphi}(U, V) := \int_{\Omega} \alpha [U \overline{\partial_x V} - \partial_x U \overline{V}] \partial_x \varphi d\mathbf{x}, \quad (5.10)$$

with $\varphi$ as in definition 3.2, so that $\partial_x\varphi$ vanishes in the vicinity of $\Sigma$. Using the definitions of $C_\varphi$ and $b_{sing}^{(1)}$ (cf. (4.3)), we rewrite the sesquilinear form $a^{(1)}$ originally defined in (4.1). Below there is no contribution on $\Gamma_n$ or $\Gamma_p$ because $\varphi$ also vanishes in their vicinity:

$$\begin{aligned} a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = & \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha(u_j + s_{g-h}) \overline{\partial_x (v_j + s_{k-l})} \partial_x \varphi d\mathbf{x} \\ & - \int_{\Omega_j} \alpha(\overline{v_j + s_{k-l}}) \partial_x (u_j + s_{g-h}) \partial_x \varphi d\mathbf{x} \\ & - \int_{\Omega_j} \left( \alpha \partial_y s_h \overline{\partial_y (v_j + s_{k-l})} + (-\partial_x (\alpha \partial_x s_h) - \omega^2 s_h) \overline{(v_j + s_{k-l})} \right) \varphi d\mathbf{x} \\ & + \int_{\Omega_j} \left( \alpha \overline{\partial_y s_l} \partial_y (u_j + s_{g-h}) + (-\partial_x (\alpha \partial_x s_l) - \omega^2 s_l) (u_j + s_{g-h}) \right) \varphi d\mathbf{x} \\ = & C_{\varphi}(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l}) - b_{sing}^{(1)}(h, (\mathbf{v} + s_{k-l})\varphi) \\ & + b_{sing}^{(1)}(l, (\mathbf{u} + s_{g-h})\varphi). \end{aligned} \quad (5.11)$$

We will need the following two identities, whose proofs can be found in appendix C.
**Lemma 5.9.** For $\mathbf{u}, \mathbf{v} \in Q$, $C_\varphi(\mathbf{u}, \mathbf{v}) = \overline{b_{reg}^{(1)}(\mathbf{v}, \mathbf{u}\varphi)} - b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}\varphi)$.

*Proof.* This is a direct application of proposition C.2, where we used $\mathrm{supp} \, \ell^{(1)} \cap \mathrm{supp} \, \varphi = \emptyset$. $\square$

**Lemma 5.10.** Let $(\mathbf{u}, g) \in Q \times H_{per}^1(\Sigma)$ be like in proposition 5.8. Then, for any $k \in H_{per}^1(\Sigma)$,

$$C_{\varphi}(\mathbf{u}, s_k) = \overline{b_{sing}^{(1)}(k, \mathbf{u}\varphi)} + b_{sing}^{(1)}(g, s_k\varphi) - \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma}.$$

*Proof.* This is a direct application of proposition C.6, where we used $\mathrm{supp} \, \ell^{(1)} \cap \mathrm{supp} \, \varphi = \emptyset$. $\square$

With these identities, we can prove proposition 5.8.

*Proof of proposition 5.8.* By (5.11),

$$a^{(1)}((\mathbf{u}, g, g), (\mathbf{v}, k, l)) = C_{\varphi}(\mathbf{u}, \mathbf{v} + s_{k-l}) - b_{sing}^{(1)}(g, (\mathbf{v} + s_{k-l})\varphi) + \overline{b_{sing}^{(1)}(l, \mathbf{u}\varphi)}. \quad (5.12)$$

Replacing in the above $C_\varphi$ from lemmas 5.9 and 5.10 yields

$$\begin{align*}
a^{(1)}((\mathbf{u}, g, g), (\mathbf{v}, k, l)) &= \overline{b_{reg}^{(1)}(\mathbf{v}, \mathbf{u}\varphi)} - b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}\varphi) + \overline{b_{sing}^{(1)}(k-l, \mathbf{u}\varphi)} + b_{sing}^{(1)}(g, s_{k-l}\varphi) \\
&\quad - \langle [\mathbf{u}]_{\Sigma}, k-l \rangle_{\Sigma} - b_{sing}^{(1)}(g, (\mathbf{v} + s_{k-l})\varphi) + \overline{b_{sing}^{(1)}(l, \mathbf{u}\varphi)} \\
&= \overline{b_{reg}^{(1)}(\mathbf{v}, \mathbf{u}\varphi)} + \overline{b_{sing}^{(1)}(k, \mathbf{u}\varphi)} - b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}\varphi) - b_{sing}^{(1)}(g, \mathbf{v}\varphi) \\
&\quad - \langle [\mathbf{u}]_{\Sigma}, k-l \rangle_{\Sigma}.
\end{align*}$$

The sum of the first two terms above is by definition $\overline{b^{(1)}((\mathbf{v}, k, l), \mathbf{u}\varphi)}$. The sum of the next two terms gives $\ell^{(1)}(\mathbf{v}\varphi)$ by the assumption of the proposition, and vanishes because $\mathrm{supp} \, \ell^{(1)} \cap \mathrm{supp} \, \varphi = \emptyset$. This allows us to conclude. $\square$

### 5.3.2 Properties of the limiting absorption solution

Given $g \in H_{per}^1(\Sigma)$, let us introduce some “artificial singularities” with non-zero absorption $\lambda$

$$s_g^{\lambda}(x, y) = g(y) \log \left( x + \frac{i\lambda}{\mathbf{r}(y)} \right).$$

For positive $\lambda$, one recovers the "singularities with absorption" of (3.1). We remark that one has convergence almost everywhere as $\nu \to 0+$:

$$s_g^{\nu} \xrightarrow{a.e.} s_g^{+} := g(\log|x| + i\pi\mathbb{1}_{x<0}), \quad \text{and} \quad s_g^{-} \xrightarrow{a.e.} s_g^{-} := g(\log|x| - i\pi\mathbb{1}_{x<0}).$$

We then have the following lemma (which generalizes lemma 3.1 to the case of artificial singularities), whose proof is left to the reader.

**Lemma 5.11.** Given $g \in H_{per}^1(\Sigma)$, the following limits hold in $L^2(\Omega)$ as $\nu \to 0+$:

$$\begin{align*}
s_g^{\pm\nu} &\to s_g^{\pm}, & \partial_y s_g^{\pm\nu} &\to \partial_y s_g^{\pm}, \\
(\alpha \pm i\nu)\partial_x s_g^{\pm\nu} &\to \alpha\partial_x s_g^{\pm}, & \partial_x((\alpha \pm i\nu)\partial_x s_g^{\pm\nu}) &\to \partial_x(\alpha\partial_x s_g^{\pm}).
\end{align*}$$
Note that $s_g = s_g^+$. We adopt this convention from now on. From the above lemma, it follows in particular that for any $\psi \in C^\infty(\bar{\Omega})$,

$$ \partial_x((\alpha \pm i\nu)\partial_x(s_g^{\pm\nu}\psi)) \xrightarrow[\nu \to 0+]{L^2(\Omega)} \partial_x(\alpha\partial_x(s_g^{\pm}\psi)) \quad \text{in } L^2(\Omega). \qquad (5.13) $$

Going back to our problem, one can prove the result (see lemma A.1)

$$ (\alpha + i\nu) \nabla u^{\nu} \xrightarrow[\nu \to 0+]{L^2(\Omega)} \alpha \nabla (u_{reg}^{+} + u_{sing}^{+}). \qquad (5.14) $$

Then, one has the result below, which follows directly from assumption 2.2 on the limiting absorption solution and lemma 5.11. As argued in proposition 5.8, in order to construct a solution to (4.6), it is sufficient to

find $\mathbf{u} \in Q, g \in H_{per}^1(\Sigma)$ such that

$$ b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{sing}^{(1)}(g, \mathbf{v}) = \ell^{(1)}(\mathbf{v}), \quad \forall \mathbf{v} \in Q, $$

$$ [\mathbf{u}]_{\Sigma} = 0. $$

Let us show that the limiting absorption solution verifies the sufficient conditions above. The problem with absorption (2.4) can be written in the variational form:

$$ \text{find } u^{\nu} \in H^1(\Omega) \text{ such that} \qquad (5.15) \\ b^{\nu}(u^{\nu}, v) = \ell^{(1)}(v), \quad \forall v \in H^1(\Omega), $$

$$ \text{with } b^{\nu}(u, v) = \int_{\Omega} [(\alpha + i\nu) \nabla u \cdot \overline{\nabla v} - \omega^2 u \overline{v}] \, dx + i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u \overline{v} ds. $$

First, we remark that the limiting absorption solution $(\mathbf{u}^+, g^+)$ verifies the desired variational formulation.

**Lemma 5.12.** Let $(\mathbf{u}^+, g^+) \in Q \times H_{per}^1(\Sigma)$ be as in assumption 2.2. Then, for any $\mathbf{v} \in Q$,

$$ b_{reg}^{(1)}(\mathbf{u}^+, \mathbf{v}) + b_{sing}^{(1)}(g^+, \mathbf{v}) = \ell^{(1)}(\mathbf{v}). $$

*Proof.* See lemma A.3 (Appendix A). $\square$

Next, we show that the limiting absorption solution $(\mathbf{u}^+, g^+)$ has a vanishing jump.

**Proposition 5.13.** Let $(\mathbf{u}^+, g^+) \in Q \times H_{per}^1(\Sigma)$ be as in assumption 2.2. Then $[\mathbf{u}^+]_{\Sigma} = 0$.

*Proof.* Let $\mathbf{u}^+, g^+$ be like in assumption 2.2. To prove that $[\mathbf{u}^+]_{\Sigma} = 0$, we will use the identity (5.6) defining the jump with $s_k^-$, for a given $k \in H_{per}^1(\Sigma)$. More precisely, let $\varphi$ be a truncation function as in the definition 3.2 and, for $\varepsilon > 0$, $\varphi_{\varepsilon}(x, y) = \varphi(\frac{x}{\varepsilon}, y)$. We will show that the quantity below is well-defined and converges to 0 as $\varepsilon \to 0+$:

$$ J_{\varepsilon}(k) = \sum_{j \in \{p,n\}} \int_{\Omega_j} u_j^{+} \overline{(-\partial_x(\alpha \partial_x (s_k^{-} \varphi_{\varepsilon})))} d\mathbf{x}. $$
We reexpress $J_\varepsilon(k)$ with the help of (5.13) and the convergence $(u_{reg}^\nu)_\nu$ to $u^+$ of Lemma A.2:

$$J_\varepsilon(k) = \lim_{\nu \to 0^+} J_\varepsilon^\nu(k), \quad \text{with} \quad J_\varepsilon^\nu(k) = \sum_{j \in \{p,n\}} \int_{\Omega_j} u_{reg}^\nu \overline{(-\partial_x ((\alpha - i\nu) \partial_x (s_k^{-\nu} \varphi_\varepsilon)))} d\mathbf{x} \quad \text{for } \varepsilon > 0.$$

The main idea of the proof consists in reexpressing $J_\varepsilon^\nu$ via $b^\nu(u^\nu, s_k^{-\nu} \varphi_\varepsilon)$. By (5.15) and the decomposition $u^\nu = u_{reg}^\nu + s_{g^+}^\nu$, defined in Lemma A.2,

$$b^\nu(u_{reg}^\nu, s_k^{-\nu} \varphi_\varepsilon) + b^\nu(s_{g^+}^\nu, s_k^{-\nu} \varphi_\varepsilon) = \ell^{(1)}(s_k^{-\nu} \varphi_\varepsilon) = 0, \quad (5.16)$$

since $\supp \ell^{(1)} \cap \supp \varphi_\varepsilon = \emptyset$. One has, by integrating by parts in the $x$-direction,

$$\begin{aligned} b^\nu(u_{reg}^\nu, s_k^{-\nu} \varphi_\varepsilon) &= \int_{\Omega} \partial_x u_{reg}^\nu \overline{(\alpha - i\nu) \partial_x (s_k^{-\nu} \varphi_\varepsilon)} d\mathbf{x} + \sum_{j \in \{p,n\}} \int_{\Omega_j} \left[ (\alpha + i\nu) \partial_y u_{reg}^\nu \overline{\partial_y (s_k^{-\nu} \varphi_\varepsilon)} - \omega^2 u_{reg}^\nu \overline{s_k^{-\nu} \varphi_\varepsilon} \right] d\mathbf{x} \\ &= J_\varepsilon^\nu(k) + \sum_{j \in \{p,n\}} \int_{\Omega_j} \left[ (\alpha + i\nu) \partial_y u_{reg}^\nu \overline{\partial_y s_k^{-\nu}} - \omega^2 u_{reg}^\nu \overline{s_k^{-\nu}} \right] \varphi_\varepsilon d\mathbf{x}. \end{aligned}$$

Indeed, the boundary terms vanish due to the choice of $\varphi$, which is supported away from $\Gamma_n$ and $\Gamma_p$. As $\nu \to 0^+$, by lemmas A.2, 5.11 and the limit (5.13), it holds that

$$b^\nu(u_{reg}^\nu, s_k^{-\nu} \varphi_\varepsilon) \to J_\varepsilon(k) + I_\varepsilon, \quad I_\varepsilon = \sum_{j \in \{p,n\}} \int_{\Omega_j} \left[ \alpha \partial_y u_j^+ \overline{\partial_y s_k^-} - \omega^2 u_j^+ \overline{s_k^-} \right] \varphi_\varepsilon d\mathbf{x}. \quad (5.17)$$

Next let us consider the second term in (5.16). Performing once again integration by parts in the $x$-direction, one finds

$$b^\nu(s_{g^+}^\nu, s_k^{-\nu} \varphi_\varepsilon) = \int_{\Omega} \left[ (\alpha + i\nu) \partial_y s_{g^+}^\nu \overline{\partial_y (s_k^{-\nu} \varphi_\varepsilon)} + \left( \left( -\partial_x \left( (\alpha + i\nu) \partial_x s_{g^+}^\nu \right) \right) - \omega^2 s_{g^+}^\nu \right) \overline{(s_k^{-\nu} \varphi_\varepsilon)} \right] d\mathbf{x},$$

and by lemma 5.11, as $\nu \to 0^+$,

$$b^\nu(s_{g^+}^\nu, s_k^{-\nu} \varphi_\varepsilon) \to b_{sing}^{(1)}(g^+, s_k^- \varphi_\varepsilon) = \int_{\Omega} \left[ \alpha \partial_y s_{g^+}^+ \overline{\partial_y (s_k^- \varphi_\varepsilon)} + \left( \left( -\partial_x (\alpha \partial_x s_{g^+}^+) \right) - \omega^2 s_{g^+}^+ \right) \overline{(s_k^- \varphi_\varepsilon)} \right] d\mathbf{x}. \quad (5.18)$$

Finally, by Lebesgue's dominated convergence theorem, as $\varepsilon \to 0^+$, $I_\varepsilon$, $b_{sing}^{(1)}(g^+, s_k^- \varphi_\varepsilon)$ both go to 0. Therefore, combining (5.17) and (5.18) in (5.16), and taking $\varepsilon \to 0^+$, we obtain that

$$\lim_{\varepsilon \to 0^+} J_\varepsilon(k) = 0,$$

which leads to the conclusion thanks to the alternate definition of the jump (5.6). $\square$

### 5.3.3 Existence and uniqueness of the solution

We start with the following proposition about the existence of the solution to (4.6).

**Proposition 5.14.** Let $(u^+, g^+)$ be like in assumption 2.2. Then $(u^+, g^+, g^+, u^+)$ satisfies (4.6).
Proof. It suffices to verify that the limiting absorption solution $(\mathbf{u}^+, g^+)$ as defined in assumption 2.2 satisfies the assumptions of proposition 5.8, with $[\mathbf{u}^+]_{\Sigma} = 0$. This follows from lemma 5.12 and proposition 5.13. □

The principal result of this section is then summarized below.

**Theorem 5.15.** Given $(\mathbf{u}^+, g^+)$ satisfying assumption 2.2, $(\mathbf{u}^+, g^+, g^+, \mathbf{u}^+\varphi)$ is the unique solution of (4.6).

*Proof.* The solution to the mixed formulation (4.6) is unique by theorem 5.4 and exists by proposition 5.14. □

**Remark 5.16.** The mixed formulation takes its origin in the minimization of a functional $\lim_{\nu \to 0+} \mathcal{J}^{\nu}$, see (3.3). The minimum of this functional is achieved in particular when $h = g^+$, and thus it is unsurprising that the Lagrange multiplier $h$ is chosen as $g^+$ in the above. As for the explicit form $\lambda = \mathbf{u}^+\varphi$, it follows from the computations.

# 6 Numerical experiments

Below, we study the numerical approximation of (4.6), or of its stabilized version (4.7). In order to test the accuracy of the method described in section 4, we reproduce the experiments conducted in section 3.2, replacing the discrete space $V_{h,h}^{(2)}$ with the space $V_h^{(1)} = Q_h \times H_h^1 \times H_h^1$,

$$ 
\begin{aligned}
Q_h &= \{v_h \in Q : v_h|_K \in P_1(K), \text{ for all } K \in \mathcal{T}_h^{\Omega}\}, \\
H_h^1 &= \{g_h \in H_{per}^1(\Sigma) : g_h|_K \in P_1(K), \text{ for all } K \in \mathcal{T}_h^{\Sigma}\}.
\end{aligned}
$$

Above, $\mathcal{T}_h^{\Omega}$ is a triangulation of $\Omega$ that is conforming with respect to the interface $\Sigma$, and $\mathcal{T}_h^{\Sigma}$ is a triangulation of $\Sigma$, both with meshsize $h$, however we do not impose that $\mathcal{T}_h^{\Sigma}$ is the trace of $\mathcal{T}_h^{\Omega}$ on $\Sigma$. Different triangulations $\mathcal{T}_h^{\Omega}$ are used, which are all symmetric with respect to the interface $\Sigma$, and we choose equidistant triangulations $\mathcal{T}_h^{\Sigma}$. Like in section 3.2, elements of $Q_h$ have no matching condition at the interface. The relative errors $e_{L^2}$ and $e_Q$ are the same as those of section 3.2. Furthermore, we denote by $e_{H^s}$ is the relative error in a discrete $H^s(\Sigma)$-norm, with $s \in \{-1, -1/2, 1\}$.

The code is written in **FreeFem++** [10]. Whereas 2D HCT finite elements were used to discretize the singular part $g_h^+$ in [13], we now use $P_1$ Lagrange finite elements on the interface $\Sigma$.<fnref n="3" />

We consider the same setting as in section 3.2, where $\alpha(x, y) = x$ and $\omega = 0$, with a purely regular solution $u^+(x, y) = 1$, and with a singular solution $u^+(x, y) = -K_0(\pi x)e^{i\pi y}$ where the singular coefficient is equal to $g^+(y) = e^{i\pi y}$.

The cutoff function we use is the $C^1(\Omega)$ function $\varphi(x, y) = \frac{1}{2}(1 + \cos(2\pi x)) \mathbb{1}_{|x|<0.5}$. Notice that $\varphi$ is prescribed equal to 1 only on the interface in the experiments (compare with definition 3.2). The approach in section 4 and its theoretical justification in section 5 remain valid also for this choice of $\varphi$. It has been also checked numerically that the results presented below do not depend on the choice of $\varphi$ provided that $\varphi \in C^1(\Omega)$, $\partial_y\varphi = 0$, $\varphi|_{\Sigma} = 1$ and is compactly supported in $x \in (-a, a)$.

<footnote n="3">³To our knowledge, in **FreeFem++**, it is not possible to combine 1D and 2D discretizations. So, in practice, to represent elements of $H_h^1$, we use $P_1$ Lagrange finite elements on a single elongated cell in the $x$-direction, and as many cells in the $y$-direction as there are in $\mathcal{T}_h^{\Sigma}$, with periodic conditions in $x$.</footnote>
**Influence of the triangulation.** The design of $\mathcal{T}_h^\Omega$ has a noticeable influence on the numerical stability of the method. In particular, we observe that the method is unstable with an unstructured triangulation $\mathcal{T}_h^{\Omega,unstr}$, see figure 4a, even though it is symmetric. On the same figure, we see that one can stabilize the method by using a structured triangulation $\mathcal{T}_h^{\Omega, str}$, as long as structuring occurs on the geometrical support of $\varphi$ (see figure 4b).

<figure><img src="image_5.png" /><figcaption>Figure 4: Influence of structuring $\mathcal{T}_h^\Omega$ on the stability of the method.</figcaption></figure>

**Numerical convergence and stabilization parameter.** We observe on figure 5 that the method using $H_1^1$ performs significantly better than the one using $H_{h2}^2$, compare with figures 2 and 3. Let us remark that, as before, when computing $e_{L^2}(\mathbf{u}^+)$ and $e_Q(\mathbf{u}^+)$, we exclude cells that are adjacent to the interface.

In figure 5a, we notice that the errors increase when decreasing $h$: this is likely due to the fact that already at the most coarse discretization the machine precision had been reached, and for finer discretizations we can observe the effects of the round-off errors in cells close to the interface (note that the total area of the excluded cells diminishes linearly with the meshsize).

We observe on figure 6 that the relative error on $\mathbf{u}_{\rho_1,h}$ decreases proportionally to the stabilization parameter $\rho_1$. Moreover, one can still compute the discrete solution for $\rho_1 = 0$ and, for the range of mesh sizes that we use, it gives the same results as those obtained for $\rho_1 = 10^{-6}$. The latter is due to the fact that, for the chosen mesh sizes, the error due to stabilization is negligible, even for "small" values of $\rho_1$. On the other hand, we have proved that the solution to the (non-stabilized) variational formulation (4.6) exists and is unique, see theorem 5.15. So uniqueness also holds for its conforming discretization: as a result, the discrete solution exists.

**Experiments.** Now, we take the same geometry, with $\alpha(x, y) = x(1 + \frac{1}{2}\cos(\pi y)) + \frac{x^2}{2}\cos(\pi y)$, $\omega = 0$ and data $f(x, y) = i(\mathbb{1}_{\{1\}}\times(0,1) - \mathbb{1}_{\{-1\}}\times(-1,0))$. Remark that $\alpha$ depends on $y$ non-trivially, and thus the exact solution is not known. According to section 5.3.2, given the limiting absorption solution $(\mathbf{u}^+, g^+)$, which has a vanishing jump according to proposition 5.13, $(\mathbf{u}^+, g^+, g^+, \mathbf{u}^+\varphi)$ is equal to the solution $(\mathbf{u}, g, h, \boldsymbol{\lambda})$ of (4.6). Therefore, we expect that $\mathbf{u}_h\varphi - \boldsymbol{\lambda}_h$, $g_h - h_h$ and finally $[\mathbf{u}_h]_\Sigma$ go to zero in the appropriate norm $\|\cdot\|_{\bullet}$ when the triangulations $\mathcal{T}_h^\Sigma$ and $\mathcal{T}_h^{\Omega, str}$ are refined. We will refer to the norms of these quantities as indicators. First, we observe that the value of each
<figure><img src="image_6.png" /><figcaption>Figure 5: Relative errors with $T_h^{\Omega, str}$. Notice the difference in the scale of the two figures.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: $e_{L^2}(\mathbf{u}^+)$ errors with $T_h^{\Omega, str}$ for $u^+(x, y) = -K_0(\pi x)e^{i\pi y}$.</figcaption></figure>

norm $||\lambda_h||$, $||g_h||$ and $||h_h||$ stabilizes quickly with respect to the meshsize $h$. Hence, in figure 7, we can report the relative errors defined by

$$d_{\bullet}(\mathbf{u}\varphi, \lambda) = \frac{\|\mathbf{u}_h\varphi - \lambda_h\|_{\bullet}}{\|\lambda_h\|_{\bullet}}, \quad \text{or} \quad d_{\bullet}(g, h) = \frac{\|g_h - h_h\|_{\bullet}}{\|g_h\|_{\bullet}}.$$

In figures 7a and 7b, we see that the first two indicators converge nicely to 0. Regarding the last indicator (the norm of the jump $[\mathbf{u}_h]_{\Sigma}$), we observe in figure 7c that it does not converge in $L^2(\Sigma)$-norm (or, if it does, the convergence is extremely slow), but it clearly does so in discrete $H^{-1}(\Sigma)$- and $H^{-1/2}(\Sigma)$-norms<fnref n="4" />. This indicates that the jump must indeed be handled carefully.

<footnote n="4">⁴The discrete $H^{-s}$-norm of $[\mathbf{u}_h]_{\Sigma}$ is evaluated as $h^s \|[\mathbf{u}_h]_{\Sigma}\|_{L^2(\Sigma)}$, see [1].</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 7: Experiment with $\alpha(x, y)$ that depends non-trivially on $y$.</figcaption></figure>

# 7 Conclusion and perspectives

In this paper, we have discussed, analyzed and improved the method introduced by Nicolopoulos et al in [13] for solving a degenerate PDE with a continuously sign-changing coefficient $\alpha$. This has led to the design of an improved mixed variational formulation. Analysis has shown that this new mixed variational formulation is coherent with the limiting absorption principle under suitable assumption on the splitting of its solution. On the other hand, when it is stabilized, this formulation can be fully analyzed via classical techniques of functional analysis. Numerical examples support the analysis.

Nevertheless, several issues need to be addressed. The first one is the limiting absorption principle assumption 2.2. Of special interest is the case where the support of the data intersects with the interface. Another one consists in simplifying the mixed formulation (4.6), the analysis of which is somewhat cumbersome due to its mixed nature.

Regarding the model, the degeneracy is caused by the continuous sign-change of the coefficient $\alpha$ through an interface $\Sigma$. One could also consider another family of degeneracy, namely the sign of $\alpha$ changes but $\alpha$ also blows up through an interface. This case has been investigated in [11] for instance.
# Acknowledgements

We thank Bruno Després and Martin Campos Pinto for their support, and we would also like to show our gratitude to Anouk Nicolopoulos for providing the code she initially developed.

# References

[1] N. ALSHEHRI, D. BOFFI, AND L. GASTALDI, *Unfitted mixed finite element methods for elliptic interface problems*, (2022).

[2] D. BOFFI, F. BREZZI, AND M. FORTIN, *Mixed finite element methods and applications*, vol. 44, Springer, 2013.

[3] H. BREZIS, *Functional Analysis, Sobolev Spaces and Partial Differential Equations*, Springer, 2011.

[4] C. CALDINI-QUEIROS, B. DESPRÉS, L.-M. IMBERT-GÉRARD, AND M. KACHANOVSKA, *A numerical study of the solution of $X$-mode equations around the hybrid resonance*, in CEMRACS 2014—numerical modeling of plasmas, vol. 53 of ESAIM Proc. Surveys, EDP Sci., Les Ulis, 2016, pp. 1–21.

[5] M. CAMPOS PINTO AND B. DESPRÉS, *Constructive formulations of resonant Maxwell's equations*, SIAM J. Math. Anal., 49 (2017), pp. 3637–3670.

[6] B. DESPRÉS, L.-M. IMBERT-GÉRARD, AND O. LAFITTE, *Solutions to the cold plasma model at resonances*, J. Éc. polytech. Math., 4 (2017), pp. 177–222.

[7] B. DESPRÉS, L.-M. IMBERT-GÉRARD, AND R. WEDER, *Hybrid resonance of Maxwell's equations in slab geometry*, J. Math. Pures Appl. (9), 101 (2014), pp. 623–659.

[8] B. DESPRÉS AND R. WEDER, *Hybrid resonance and long-time asymptotic of the solution to Maxwell's equations*, Physics Letters A, 380 (2016), pp. 1284–1289.

[9] J. P. FREIDBERG, *Plasma physics and fusion energy*, Cambridge university press, 2008.

[10] F. HECHT, *New development in FreeFem++*, J. Numer. Math., 20 (2012), pp. 251–265.

[11] L.-M. IMBERT-GÉRARD, *Mathematical and numerical problems of some wave phenomena appearing in magnetic plasmas*, phd thesis, Université Pierre et Marie Curie - Paris VI, 2013.

[12] A. NICOLOPOULOS, M. CAMPOS-PINTO, AND B. DESPRÉS, *A stable formulation of resonant Maxwell's equations in cold plasma*, J. Comput. Appl. Math., 362 (2019), pp. 185–204.

[13] A. NICOLOPOULOS, M. CAMPOS PINTO, B. DESPRÉS, AND P. CIARLET, JR., *Degenerate elliptic equations for resonant wave problems*, IMA J. Appl. Math., 85 (2020), pp. 132–159.

[14] A. NICOLOPOULOS-SALLE, *Variational formulations of resonant Maxwell equations and corner problems in wave propagation*, PhD thesis, Sorbonne university, 2019.

[15] T. H. STIX, *Waves in plasmas*, Springer Science & Business Media, 1992.
In the appendices, we write $|| \cdot ||$ for the $L^2(\Omega)$-norm, and $|| \cdot ||_\infty$ for the $L^\infty(\Omega)$-norm. For the convenience of the reader, we also make an abuse of notations by writing $\mathbf{u}^+$ for $(u_{reg|\Omega_p}^+, u_{reg|\Omega_n}^+) \in Q$ and vice versa when referring to the regular part of the limiting absorption solution.

Moreover, where applicable, we will denote by $\mathbf{n} = (n_x, n_y)$ the unit normal to $\Sigma$ directed from $\Omega_n$ towards $\Omega_p$, and by $\partial_n$ the respective normal derivative.

## A Auxiliary results for the limiting absorption solution

The goal of this appendix is to prove several auxiliary results on the limiting absorption solution under assumption 2.2.

**Lemma A.1.** Let $(u^\nu)_{\nu>0}$ be a family of the solutions of (2.4). Assuming (2.5), there holds

$$ (\alpha + i\nu)\nabla u^\nu \frac{L^2(\Omega)}{\nu \to 0+} \to \alpha \nabla u^+. $$

*Proof.* The variational formulation for (2.4) reads: find $u^\nu \in H^1(\Omega)$ s.t. for all $v \in H^1(\Omega)$,

$$ \int_{\Omega} [(\alpha + i\nu) \nabla u^\nu \cdot \overline{\nabla v} - \omega^2 u^\nu \overline{v}] dx + i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u^\nu \overline{v} ds = \sum_{j \in \{p,n\}} \int_{\Gamma_j} f \overline{v} ds. \quad (A.1) $$

We split the proof into two steps.

**Step 1.** Proof that $\nu||\nabla u^\nu|| \to 0$ as $\nu \to 0+$. We test the equation (A.1) with $v^\nu = u^\nu$ and take the imaginary part of the resulting expression. This yields

$$ \nu ||\nabla u^\nu||^2 + \lambda ||u^\nu||^2_{L^2(\Gamma_p \cup \Gamma_\nu)} \le ||f||_{L^2(\Gamma_p \cup \Gamma_n)} ||u^\nu||_{L^2(\Gamma_p \cup \Gamma_\nu)}. $$

One obtains that $\lambda ||u^\nu||_{L^2(\Gamma_p \cup \Gamma_\nu)} \le ||f||_{L^2(\Gamma_p \cup \Gamma_n)}$, so $\nu ||\nabla u^\nu||^2 \le \lambda^{-1} ||f||^2_{L^2(\Gamma_p \cup \Gamma_n)}$, hence the claim.

**Step 2.** Proof that $\alpha \nabla u^\nu \to \alpha \nabla u^+$ in $L^2(\Omega)$. We will show that $(\alpha \nabla u^{\nu_n})_{n \in \mathbb{N}}$ is a Cauchy sequence for any $(\nu_n)_{n \in \mathbb{N}} \subset \mathbb{R}^+$, s.t. $\lim_{n \to +\infty} \nu_n = 0$.

For this we consider the difference of (A.1) written for $\nu = \nu_n$ and $\nu = \nu_m$, namely

$$ \int_{\Omega} [\alpha (\nabla u^{\nu_n} - \nabla u^{\nu_m}) \cdot \overline{\nabla v} + i (\nu_n \nabla u^{\nu_n} - \nu_m \nabla u^{\nu_m}) \cdot \overline{\nabla v} - \omega^2 (u^{\nu_n} - u^{\nu_m}) \overline{v}] dx \\ + i\lambda \int_{\Gamma_p \cup \Gamma_n} (u^{\nu_n} - u^{\nu_m}) \overline{v} ds = 0. \quad (A.2) $$

Let us denote $e_{nm} = (u^{\nu_n} - u^{\nu_m})$. We test the equation (A.2) with $v = \alpha e_{nm}$, which yields

$$ \int_{\Omega} [ |\alpha \nabla e_{nm}|^2 + \alpha \nabla e_{nm} \cdot \nabla \alpha \overline{e_{nm}} ] dx + \int_{\Omega} i (\nu_n \nabla u^{\nu_n} - \nu_m \nabla u^{\nu_m}) \cdot (\nabla \alpha \overline{e_{nm}} + \alpha \overline{\nabla e_{nm}}) dx \\ - \int_{\Omega} \omega^2 \alpha |e_{nm}|^2 dx + i\lambda \int_{\Gamma_p \cup \Gamma_n} \alpha |e_{nm}|^2 ds = 0. $$
Taking the real part of the above, and using the Cauchy-Schwarz inequality to bound sign-indefinite terms yields:

$$ \int_{\Omega} |\alpha \nabla e_{nm}|^2 dx \leq \|\nabla \alpha\|_{\infty} \|\alpha \nabla e_{nm}\| \|e_{nm}\| + \|\nu_n \nabla u^{\nu_n} - \nu_m \nabla u^{\nu_m}\| (\|\nabla \alpha\|_{\infty} \|e_{nm}\| + \|\alpha \nabla e_{nm}\|) \\ + \omega^2 \|\alpha\|_{\infty} \|e_{nm}\|^2. $$

With the help of the Young inequality, we obtain the following bound:

$$ \|\alpha \nabla e_{nm}\| \leq C (\|e_{nm}\| + \|\nu_n \nabla u^{\nu_n} - \nu_m \nabla u^{\nu_m}\|), $$

where the constant C depends on $\|\alpha\|_{W^{1,\infty}}$ and $\omega^2$ only. Because $\nu\|\nabla u^{\nu}\| \to 0$ as $\nu \to 0$ and $u^{\nu}$ converges in $L^2(\Omega)$ as $\nu \to 0$, we conclude that $(\alpha \nabla u^{\nu_n})_{n \in \mathbb{N}}$ is an $L^2(\Omega)$-Cauchy sequence, and thus converges. Evidently, its limit is $\alpha \nabla u^+$; this follows from the following expression (which allows to define the distribution $\alpha \nabla v$ for $v \in L^2(\Omega)$ and $\alpha \in C^1(\bar{\Omega})$):

$$ \alpha \nabla u^{\nu} = \nabla(\alpha u^{\nu}) - u^{\nu} \nabla \alpha. $$

We have that $u^{\nu} \nabla \alpha \to u^+ \nabla \alpha$ in $L^2(\Omega)$; similarly, $\alpha u^{\nu} \to \alpha u^+$, thus, in the sense of distributions, $\nabla(\alpha u^{\nu}) \to \nabla(\alpha u^+)$, and by the uniqueness of the distributional limit we conclude with the desired result. □

Let $\nu > 0$, and $u^{\nu}$ satisfy (2.4). Introducing $u_{reg}^{\nu} := u^{\nu} - s_{g^{+}}^{\nu}$ and using the previous lemma together with lemma 3.1, one has the improved result below.

**Lemma A.2.** Let $\nu > 0$ and $u^{\nu}$ satisfy (2.4). Under assumption 2.2, the unique solution $u^{\nu} \in H^1(\Omega)$ of (2.4) can be decomposed as

$$ u^{\nu} = u_{reg}^{\nu} + s_{g^{+}}^{\nu}, $$

where $s_{g^{+}}^{\nu}$ is defined as in (3.1), and it holds that

$$ u_{reg}^{\nu} \xrightarrow{\frac{L^2(\Omega)}{\nu \to 0+}} u_{reg}^{+}, \quad \text{and} \quad (\alpha + i\nu) \nabla u_{reg}^{\nu} \xrightarrow{\frac{L^2(\Omega)}{\nu \to 0+}} \alpha \nabla u_{reg}^{+}. $$

The result that follows relies on the definitions (4.5), (4.3) of the forms $b_{reg}^{(1)}$ and $b_{sing}^{(1)}$.

**Lemma A.3.** Let $(u^+, g^+) \in Q \times H_{per}^1(\Sigma)$ be governed by (2.6). Then, for any $v \in Q$,

$$ b_{reg}^{(1)}(u^+, v) + b_{sing}^{(1)}(g^+, v) = \sum_{j \in \{p,n\}} \int_{\Gamma_j} f \bar{v} \, ds. \quad (A.3) $$

*Proof.* Let us proceed as follows.<fnref n="5" /> First of all, we know that

$$ -\operatorname{div}(\alpha \nabla u^+) - \omega^2 u^+ = 0 \text{ in } \mathcal{D}'(\Omega). \quad (A.4) $$

<footnote n="5">⁵The main difficulty in the proof lies in the fact that a priori it is not clear whether $H^1(\Omega)$ functions are dense in the space $Q$, and thus passing from the variational formulation (5.15) with absorption to the variational formulation in the statement of the lemma is not completely straightforward.</footnote>
Moreover, using the decomposition of assumption 2.2 and the fact that $\partial_x(\alpha\partial_x s_{g^+}) \in L^2(\Omega)$, cf. lemma 5.11, we conclude that

$$ \mathbf{d}^+ := \operatorname{div}(\alpha \nabla u_{reg}^+) + \partial_y (\alpha \partial_y s_{g^+}) \in L^2(\Omega). \quad (A.5) $$

Testing the equation (A.4) with $v_p \in C_{per,y}^{\infty}(\bar{\Omega}_p)$ on $\Omega_p$, using the boundary conditions of (2.6) and integrating by parts yields (cf. (A.5), lemma A.2 and the periodicity of $g^+$):

$$ \int_{\Omega_p} \left\{ \alpha \partial_x u_{reg}^+ \overline{\partial_x v_p} - \partial_x (\alpha \partial_x s_{g^+}) \overline{v_p} + \alpha \partial_y (u_{reg}^+ + s_{g^+}) \overline{\partial_y v_p} - \omega^2 (u_{reg}^+ + s_{g^+}) \overline{v_p} \right\} d\mathbf{x} \\ + i\lambda \int_{\Gamma_p} (u_{reg}^+ + s_{g^+}) \overline{v_p} ds + \int_{\Gamma_p} \alpha \partial_n s_{g^+} \overline{v_p} ds - \underbrace{\langle \alpha \partial_x u_{reg}^+, \gamma_0 v_p \rangle_{(H_{per}^{1/2}(\Sigma))', H_{per}^{1/2}(\Sigma)}}_{I_\Sigma} = \int_{\Gamma_p} f \overline{v_p} ds. $$

Finally, it remains to show that $I_\Sigma$ vanishes. Let $\varphi$ be a truncation function as in the definition 3.2, and, for given $\varepsilon > 0$, $\varphi_\varepsilon(x, y) = \varphi(\frac{x}{\varepsilon}, y)$. By integration by parts, we have for any $\varepsilon > 0$,

$$ I_\Sigma = \int_{\Omega_p} \left\{ \mathbf{d}^+ \overline{v_p \varphi_\varepsilon} + \alpha \partial_y (u_{reg}^+ + s_{g^+}) \overline{\partial_y (\varphi_\varepsilon v_p)} \right\} d\mathbf{x} + \int_{\Omega_p} \alpha \partial_x u_{reg}^+ \overline{\partial_x (v_p \varphi_\varepsilon)} d\mathbf{x}. $$

It suffices to consider the case $\varepsilon \to 0$. The convergence of the first integral follows from Lebesgue's dominated convergence theorem and the fact that $\partial_y \varphi_\varepsilon = 0$. As for the second integral, we can estimate it as follows (where $\Omega_p^\varepsilon = \Omega_p \cap \operatorname{supp} \varphi_\varepsilon$):

$$ \left| \int_{\Omega_p} \alpha \partial_x u_{reg}^+ \overline{\partial_x (v_p \varphi_\varepsilon)} d\mathbf{x} \right| \le \left| \int_{\Omega_p^\varepsilon} \alpha \partial_x u_{reg}^+ \overline{\partial_x v_p} \varphi_\varepsilon d\mathbf{x} \right| + \left| \int_{\Omega_p^\varepsilon} \alpha \partial_x u_{reg}^+ \overline{v_p} \partial_x \varphi_\varepsilon d\mathbf{x} \right| \\ \le C \left( \|u_{reg}^+\|_{H_{1/2}^1(\Omega_p^\varepsilon)} \|v_p\|_{H_{1/2}^1(\Omega_p^\varepsilon)} + \|u_{reg}^+\|_{H_{1/2}^1(\Omega_p^\varepsilon)} \|v_p\|_{L^\infty(\Omega_p^\varepsilon)} \|\varphi_\varepsilon\|_{H_{1/2}^1(\Omega_p^\varepsilon)} \right). $$

Remark that there exists $C_1 > 0$, s.t. for all $\varepsilon > 0$, $\|\varphi_\varepsilon\|_{H_{1/2}^1(\Omega_p^\varepsilon)} \le C_1$. Since, additionally, $u_{reg}^+ \in H_{1/2}^1(\Omega_p)$, we conclude that, as $\varepsilon \to 0+$,

$$ \left| \int_{\Omega_p} \alpha \partial_x u_{reg}^+ \overline{\partial_x (v_p \varphi_\varepsilon)} d\mathbf{x} \right| \to 0. $$

Therefore, $I_\Sigma = 0$, and hence for all $v_p \in C_{per,y}^{\infty}(\bar{\Omega}_p)$, it holds that

$$ \int_{\Omega_p} \left\{ \alpha \partial_x u_{reg}^+ \overline{\partial_x v_p} + \alpha \partial_y (u_{reg}^+ + s_{g^+}) \overline{\partial_y v_p} - \partial_x (\alpha \partial_x s_{g^+}) \overline{v_p} - \omega^2 (u_{reg}^+ + s_{g^+}) \overline{v_p} \right\} d\mathbf{x} \\ + i\lambda \int_{\Gamma_p} (u_{reg}^+ + s_{g^+}) \overline{v_p} ds + \int_{\Gamma_p} \alpha \partial_n s_{g^+} \overline{v_p} ds = \int_{\Gamma_p} f \overline{v_p} ds. $$

Repeating the argument for $v_n \in C_{per,y}^{\infty}(\bar{\Omega}_n)$, we conclude that a similar identity holds true in $\Omega_n$. By density of the functions $C_{per,y}^{\infty}(\bar{\Omega}_n) \times C_{per,y}^{\infty}(\bar{\Omega}_p)$ in $Q$, we arrive at the formulation (A.3). $\square$
**Lemma A.4.** Let $(u^\nu)_\nu$ be a family governed by (2.4) fulfilling assumption 2.2, and $\varphi \in C^1(\Omega)$ be such that $\text{supp} \varphi \cap \overline{\Gamma_p} \cup \overline{\Gamma_n} = \emptyset$. Then

$$ \lim_{\nu \to 0+} \int_{\Omega} \nu |\nabla u_{reg}^\nu|^2 \varphi \, dx = 0. $$

*Proof.* Firstly, remark that $\int_{\Omega} \nu |\nabla u_{reg}^\nu|^2 \varphi \, dx = \text{Im} \mathcal{E}_{reg}^\nu$ with

$$ \mathcal{E}_{reg}^{\nu} = \int_{\Omega} \{ (\alpha + i\nu) |\nabla u_{reg}^{\nu}|^2 \varphi - \omega^2 |u_{reg}^{\nu}|^2 \varphi \} d\mathbf{x}. $$

Therefore, using that $\nabla u_{reg}^\nu \varphi = \nabla(u_{reg}^\nu \varphi) - u_{reg}^\nu \nabla \varphi$ and that $u^\nu = u_{reg}^\nu + s_g^\nu$ is a weak solution of the problem with absorption (2.4), one has

$$ \begin{aligned} \mathcal{E}_{reg}^{\nu} &= \int_{\Omega} \{ (\alpha + i\nu) \nabla u_{reg}^{\nu} \cdot \overline{\nabla(u_{reg}^{\nu} \varphi)} - \omega^2 u_{reg}^{\nu} \overline{u_{reg}^{\nu} \varphi} \} d\mathbf{x} - \int_{\Omega} (\alpha + i\nu) \nabla u_{reg}^{\nu} \cdot \overline{u_{reg}^{\nu} \nabla \varphi} d\mathbf{x} \\ &= - \int_{\Omega} \{ (\alpha + i\nu) \nabla s_g^{\nu} \cdot \overline{\nabla(u_{reg}^{\nu} \varphi)} - \omega^2 s_g^{\nu} \overline{u_{reg}^{\nu} \varphi} \} d\mathbf{x} - \int_{\Omega} (\alpha + i\nu) \nabla u_{reg}^{\nu} \cdot \overline{u_{reg}^{\nu} \nabla \varphi} d\mathbf{x} \\ &= - \int_{\Omega} \{ (\alpha + i\nu) \partial_y s_g^{\nu} \overline{\partial_y(u_{reg}^{\nu} \varphi)} + [-\partial_x ((\alpha + i\nu)\partial_x s_g^{\nu})] \overline{u_{reg}^{\nu} \varphi} - \omega^2 s_g^{\nu} \overline{u_{reg}^{\nu} \varphi} \} d\mathbf{x} \\ &\quad - \int_{\Omega} (\alpha + i\nu) \nabla u_{reg}^{\nu} \cdot \overline{u_{reg}^{\nu} \nabla \varphi} d\mathbf{x}, \end{aligned} $$

where an integration by parts in the x-direction is made in the last equality. According to lemma 3.1 and lemma A.2, the following convergences hold in $L^2(\Omega)$ as $\nu \to 0+$:

$s_g^\nu \to s_g$, $\partial_y s_g^\nu \to \partial_y s_g$, $-\partial_x ((\alpha + i\nu)\partial_x s_g^\nu) \to -\partial_x (\alpha \partial_x s_g)$, $u_{reg}^\nu \to u_{reg}^+(\alpha + i\nu)\nabla u_{reg}^\nu \to \alpha \nabla u_{reg}^+$, so that $\mathcal{E}_{reg}^\nu$ converges to (see (4.5), (4.3) for the definitions of the forms $b_{reg}^{(1)}$ and $b_{sing}^{(1)}$):

$$ \begin{aligned} \mathcal{E}_{reg}^{+} &= -b_{sing}^{(1)}(g^{+}, \mathbf{u}^{+}\varphi) - \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha \nabla u_j^{+} \cdot \overline{u_j^{+} \nabla \varphi} d\mathbf{x} \\ &= b_{reg}^{(1)}(\mathbf{u}^{+}, \mathbf{u}^{+}\varphi) - \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha \nabla u_j^{+} \cdot \overline{u_j^{+} \nabla \varphi} d\mathbf{x} = \int_{\Omega} \{ \alpha |\nabla u_j^{+}|^2 \varphi - \omega^2 |u_j^{+}|^2 \varphi \} d\mathbf{x}, \end{aligned} $$

where are used the facts that $(\mathbf{u}^+, g^+)$ is a weak solution of (2.6) by lemma A.3, $\nabla(u_j^+ \varphi) - u_j^+ \nabla \varphi = \nabla u_j^+ \varphi$ and the condition on the disjoint supports. Finally, considering $\text{Im} \mathcal{E}_{reg}^+$ gives the desired result. $\square$

## B From the energy functional to the mixed formulation

Let $\varphi$ be as in definition 3.2, $(u^\nu)_\nu$ such that $u^\nu \in H^1(\Omega)$ and $-\text{div}((\alpha + i\nu)\nabla u^\nu) - \omega^2 u^\nu = 0$ in $\text{supp} \varphi$, and $h \in H_{per}^2(\Sigma)$. Recall the energy functional defined in (3.3):

$$ \mathcal{J}^{\nu}(u^{\nu}, h) = \int_{\Omega} \nu |\nabla (u^{\nu} - s_{h}^{\nu})|^2 \varphi \, d\mathbf{x} = \text{Im} \mathcal{E}^{\nu}, \text{ where} $$

$$ \mathcal{E}^{\nu} = \int_{\Omega} \left[ (\alpha(x, y) + i\nu) |\nabla (u^{\nu} - s_h^{\nu})|^2 - \omega^2 |u^{\nu} - s_h^{\nu}|^2 \right] \varphi \, d\mathbf{x}. $$
Our goal is to minimize the limit $\mathcal{J}^+$ of $\mathcal{J}^\nu(u^\nu, h)$ when $\nu$ goes to $0^+$.

In order to compute the limit of $\mathcal{J}^\nu(u^\nu, h)$, we will integrate by parts the expression for $\mathcal{E}^\nu$. First, using definition 3.2 of $\varphi$, i.e., $\partial_y\varphi = 0$, and the identity $\nabla U\varphi = \nabla(U\varphi) - U\nabla\varphi$, one can rewrite $\mathcal{E}^\nu$ as

$$ \mathcal{E}^\nu = \int_{\Omega} \left[ (\alpha(x, y) + i\nu)\nabla (u^\nu - s_h^\nu) \cdot \overline{\nabla ((u^\nu - s_h^\nu)\varphi)} - \omega^2 (u^\nu - s_h^\nu) \overline{((u^\nu - s_h^\nu)\varphi)} \right] d\mathbf{x} \\ - \int_{\Omega} \left[ (\alpha(x, y) + i\nu)\partial_x (u^\nu - s_h^\nu) \overline{\overline{(u^\nu - s_h^\nu)}\partial_x\varphi} \right] d\mathbf{x}. $$

Then we separate

$$ \mathcal{E}^\nu = e^\nu(u^\nu, (u^\nu - s_h^\nu)\varphi) - e^\nu(s_h^\nu, (u^\nu - s_h^\nu)\varphi) - c^\nu(u^\nu - s_h^\nu, u^\nu - s_h^\nu), $$

where

$$ e^\nu(u, v) = \int_{\Omega} \left[ (\alpha(x, y) + i\nu)\nabla u \cdot \overline{\nabla v} - \omega^2 u\bar{v} \right] d\mathbf{x}, $$

$$ c^\nu(u, v) = \int_{\Omega} (\alpha(x, y) + i\nu) \partial_x u \bar{v} \partial_x \varphi d\mathbf{x}. $$

Since $u^\nu$ is a weak solution of the problem with absorption (2.4), $(u^\nu - s_h^\nu)\varphi \in H^1(\Omega)$ and $\text{supp}\,\varphi \subset (-a, a)$ (cf. definition 3.2 for $\varphi$), we have that for all $\nu > 0$, $e^\nu(u^\nu, (u^\nu - s_h^\nu)\varphi) = 0$. It remains to integrate by parts the term $e^\nu(s_h^\nu, (u^\nu - s_h^\nu)\varphi)$, which is allowed since $h \in H_{per}^2(\Sigma)$. Using the boundary conditions, we find that:

$$ e^\nu(s_h^\nu, \varphi(u^\nu - s_h^\nu)) = \int_{\Omega} \left[ -\text{div}((\alpha(x, y) + i\nu)\nabla s_h^\nu) - \omega^2 s_h^\nu \right] \overline{(u^\nu - s_h^\nu)} \varphi d\mathbf{x}. $$

Given that $s_h^\nu = h \log\left(x + \frac{i\nu}{\bar{r}(y)}\right)$, according to lemma 3.1 and noting that $\partial_y((\alpha + i\nu)\partial_y s_h^\nu) \to \partial_y(\alpha\partial_y s_h)$ for $h \in H_{per}^2(\Sigma)$, the following limits hold in $L^2(\Omega)$:

$$ s_h^\nu \to s_h, \quad (\alpha(x, y) + i\nu) \partial_x s_h^\nu \to \alpha \partial_x s_h, \quad \text{div}((\alpha(x, y) + i\nu)\nabla s_h^\nu) \to \text{div}(\alpha \nabla s_h). $$

Therefore, using an assumption 2.2 and lemma A.1, according to the limits above and [3, Prop.3.5], we have that

$$ e^\nu(s_h^\nu, \varphi(u^\nu - s_h^\nu)) \to e(s_h, \mathbf{u}^+ + s_{g^+-h}) \quad \text{and} \quad c^\nu(u^\nu - s_h^\nu, u^\nu - s_h^\nu) \to c(\mathbf{u}^+ + s_{g^+-h}, \mathbf{u}^+ + s_{g^+-h}), $$

where the sesquilinear forms are respectively given by

$$ e(s_h, v) = \sum_{j \in \{p, n\}} \int_{\Omega_j} \left[ -\text{div}(\alpha \nabla s_h) \bar{v} - \omega^2 s_h \bar{v} \right] \varphi d\mathbf{x}, $$

$$ c(u, v) = \sum_{j \in \{p, n\}} \int_{\Omega_j} \alpha(x, y) \partial_x u \bar{v} \partial_x \varphi d\mathbf{x}. $$

We observe that the limit depends on the triple $(\mathbf{u}^+, g^+, h)$. Then we define

$$ \mathcal{J}^+(\mathbf{u}^+, g^+, h) := -\operatorname{Im} e(s_h, \mathbf{u}^+ + s_{g^+-h}) - \operatorname{Im} c(\mathbf{u}^+ + s_{g^+-h}, \mathbf{u}^+ + s_{g^+-h}). $$
From this point on, one finds by integrating by parts that $\mathcal{J}^+(\mathbf{u}^+, g^+, g^+) = 0$ (this is reminiscent of the proof of lemma A.4). Since we know that the limit of $\mathcal{J}^\nu(u^\nu, h)$ is non-negative, we conclude that $(\mathbf{u}^+, g^+, g^+)$ is a minimizer of $\mathcal{J}^+$. Finally, it is straightforward to compute the differential of $\mathcal{J}^+$:

$$d\mathcal{J}^+((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = \operatorname{Im} a^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l))$$

where $a^{(2)}$ is a sesquilinear form defined on $V^{(2)} \times V^{(2)}$ by

$$a^{(2)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = \overline{e(s_l, \mathbf{u} + s_{g-h})} - e(s_h, \mathbf{v} + s_{k-l}) + \overline{c(\mathbf{v} + s_{k-l}, \mathbf{u} + s_{g-h})} - c(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l}),$$

cf. the definition (3.6).

## C Auxiliary Green's identities

Below, we denote by $\mathbf{u}$ the pair $(u_p, u_n) \in Q$. Recall that the term $C_\varphi(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l})$ appears in the expression (5.11) of the form $a^{(1)}$, with

$$C_{\psi}(U, V) := \int_{\Omega} \alpha [U \overline{\partial_x V} - \partial_x U \overline{V}] \partial_x \psi d\mathbf{x}. \qquad (C.1)$$

The goal of this appendix is to express $C_\psi(U, V)$ using the sesquilinear forms $b_{reg}^{(1)}$ and $b_{sing}^{(1)}$.

The first step is the following manipulation, which will be used elsewhere:

$$\begin{aligned} [u \overline{\partial_x v} - \partial_x u \overline{v}] \partial_x \psi &= \partial_x(u\psi)\overline{\partial_x v} - \partial_x u \overline{\partial_x v} \psi - \partial_x u \overline{\partial_x (v\psi)} + \partial_x u \overline{\partial_x v} \psi \quad (C.2) \\ &= \partial_x(u\psi)\overline{\partial_x v} - \partial_x u \overline{\partial_x (v\psi)}. \end{aligned}$$

Therefore, given a Lipschitz domain $\mathcal{O}$ and $u, v$ smooth in $\overline{\mathcal{O}}$, we have

$$\int_{\mathcal{O}} \alpha [u \overline{\partial_x v} - \partial_x u \overline{v}] \partial_x \psi d\mathbf{x} = \int_{\mathcal{O}} \alpha [\partial_x(u\psi)\overline{\partial_x v} - \partial_x u \overline{\partial_x (v\psi)}] d\mathbf{x}. \quad (C.3)$$

From now on, unless stated otherwise, we will assume that $\psi$ satisfies the following assumption.

**Assumption C.1.** Let $\psi \in C^1(\overline{\Omega}; \mathbb{R})$ be such that $\partial_y \psi = 0$.

With this assumption, we obviously have that

$$\partial_y(u\psi) \partial_y v = \partial_y u \partial_y (v\psi) = \partial_y u \partial_y v \psi. \qquad (C.4)$$

We observe that, depending on whether $U, V$ are regular, i.e., belonging to $H_{1/2}^1(\Omega_p) \times H_{1/2}^1(\Omega_n)$, or singular, i.e., of the form $s_g$, with $g \in H_{per}^1(\Sigma)$, the expression (C.1) of $C_\psi(U, V)$ will obviously change. There are three different cases:

* $U, V$ are both regular, in $Q$, see proposition C.2,

* $U, V$ are both singular, i.e., $U = s_g$ and $V = s_k$, see proposition C.4,

* $U$ is regular and $V$ is singular, see proposition C.6.
The simplest case is when $U, V \in Q$. According to the above, using (C.4) with $u = U$ and $v = V$ allows us reexpress the right-hand side of the identity (C.3) with $\mathcal{O} = \text{int}(\bar{\Omega}_p \cup \bar{\Omega}_n)$. Namely,

$$ \int_{\mathcal{O}} \alpha [U \overline{\partial_x V} - \partial_x U \overline{V}] \partial_x \psi dx = \int_{\mathcal{O}} [\alpha \nabla(U\psi) \cdot \overline{\nabla V} - \omega^2(U\psi)\overline{V}] dx \\ - \int_{\mathcal{O}} [\alpha \nabla U \cdot \overline{\nabla(V\psi)} - \omega^2 U \overline{(V\psi)}] dx. \qquad (C.5) $$

Recalling the definition (4.5) of $b_{reg}^{(1)}$:

$$ b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) = \sum_{j \in \{p,n\}} \int_{\Omega_j} (\alpha \nabla u_j \cdot \overline{\nabla v_j} - \omega^2 u_j \overline{v_j}) dx + i\lambda \int_{\Gamma_j} u_j \overline{v_j} ds, $$

one has the following proposition.

**Proposition C.2.** Let $\mathbf{u}, \mathbf{v} \in Q$ and $\psi$ satisfy assumption C.1. Then

$$ \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha (u_j \overline{\partial_x v_j} - \partial_x u_j \overline{v_j}) \partial_x \psi dx = \overline{b_{reg}^{(1)}(\mathbf{v}, \mathbf{u}\psi)} - b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}\psi) + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{v_j} \psi ds. $$

Let us now consider the second case, when $U$ and $V$ are both singular, that is $U = s_g$ and $V = s_k$. Evidently, we cannot apply (C.5) for $\mathcal{O} = \Omega_{p,n}$ and $\psi$ non-vanishing on the interface, since the terms $\int_{\Omega_{p,n}} \alpha \nabla(s_g\psi) \cdot \overline{\nabla s_k} dx$ and $\int_{\Omega_{p,n}} \alpha \nabla s_g \cdot \overline{\nabla(s_k\psi)} dx$ are not defined. This difficulty can be overcome by integrating by parts in the $x$-direction. Let $u, v$ be sufficiently smooth in $\overline{\mathcal{O}}$, then

$$ \int_{\mathcal{O}} \alpha [u \overline{\partial_x v} - \partial_x u \overline{v}] \partial_x \psi dx = \int_{\mathcal{O}} [(u\psi)(\overline{-\partial_x(\alpha\partial_x v)}) - (-\partial_x(\alpha\partial_x u))\overline{v\psi}] dx \quad (C.6) $$

$$ + \int_{\partial \mathcal{O}} [u(\overline{\alpha \partial_n v}) - (\alpha \partial_n u) \overline{v}] \psi ds $$

$$ = \int_{\mathcal{O}} [\alpha \partial_y (u\psi) \overline{\partial_y v} + (u\psi)(\overline{-\partial_x(\alpha\partial_x v) - \omega^2 v})] dx \quad (C.7) $$

$$ - \int_{\mathcal{O}} [\alpha \partial_y u \overline{\partial_y (v\psi)} + (-\partial_x(\alpha\partial_x u) - \omega^2 u) \overline{v\psi}] dx $$

$$ + \int_{\partial \mathcal{O}} [u(\overline{\alpha \partial_n v}) - (\alpha \partial_n u) \overline{v}] \psi ds. $$

**Lemma C.3.** For $g, k \in H_{per}^1(\Sigma)$ and $\psi \in C^1(\bar{\Omega}; \mathbb{R})$ (not necessarily satisfying assumption C.1), for $j \in \{p, n\}$, it holds that

$$ \int_{\Omega_j} \alpha [s_g \overline{\partial_x s_k} - \partial_x s_g \overline{s_k}] \partial_x \psi dx = \int_{\Omega_j} [s_g \overline{-\partial_x(\alpha \partial_x s_k)} - (-\partial_x(\alpha \partial_x s_g)) \overline{s_k}] \psi dx \\ + \int_{\Gamma_j} [s_g \overline{(\alpha \partial_n s_k)} - (\alpha \partial_n s_g) \overline{s_k}] \psi ds + \sigma_j \int_{\Sigma} g(y) \overline{k(y)} \mathbf{r}(y) \psi(0, y) dy, $$

where $\sigma_p = 0$ and $\sigma_n = 2i\pi$.
Proof. Applying (C.6) in $O = \Omega_j^\varepsilon = \{\mathbf{x} \in \Omega_j : \mathrm{dist}(\mathbf{x}, \Sigma) > \varepsilon\}$, $\varepsilon > 0$, $j \in \{p, n\}$, with $u = s_g$ and $v = s_k$ yields

$$ \int_{\Omega_j^\varepsilon} \alpha \left[ s_g \overline{\partial_x s_k} - \partial_x s_g \overline{s_k} \right] \partial_x \psi d\mathbf{x} = \int_{\Omega_j^\varepsilon} \left[ s_g \overline{-\partial_x(\alpha \partial_x s_k)} - \overline{s_k} \overline{-\partial_x(\alpha \partial_x s_g)} \right] \psi d\mathbf{x} \\ + \int_{\Gamma_j} \left[ s_g \overline{(\alpha \partial_n s_k)} - (\alpha \partial_n s_g) \overline{s_k} \right] \psi ds - a_j I_j^\varepsilon, $$

with $a_p = 1$ and $a_n = -1$ and

$$ I_j^\varepsilon = \int_{\{x=a_j\}} \left[ s_g \overline{(\alpha \partial_x s_k)} - (\alpha \partial_x s_g) \overline{s_k} \right] \psi dy. $$

As $\varepsilon \to 0+$, the volume integrals over $\Omega_j^\varepsilon$ converge to the volume integrals over $\Omega_j$, since the integrands are obviously in $L^1(\Omega_j)$. Let us compute the remaining limit $\lim_{\varepsilon \to 0+} I_j^\varepsilon$. Recall that $s_g(x, y) = g(y) S(x)$ with $S(x) = \log|x| + i\pi 1_{x<0}$. As $\varepsilon \to 0+$,

$$ I_j^\varepsilon = \int_{\Sigma} g(y) \overline{k(y)} \frac{\alpha(a_j\varepsilon, y)}{a_j\varepsilon} \left[ S(a_j\varepsilon) - \overline{S(a_j\varepsilon)} \right] \psi(a_j\varepsilon, y) dy \to \sigma_j \int_{\Sigma} g(y) \overline{k(y)} \mathbf{r}(y) \psi(0, y) dy, $$

where $\sigma_p = 0$, $\sigma_n = 2i\pi$.

<figure><img src="image_9.png" /></figure>

The proposition below is a rewriting of formula (C.7), using the above lemma and the definition (4.3) of $b_{sing}^{(1)}$, recalled here for convenience:

$$ b_{sing}^{(1)}(g, \mathbf{v}) = \sum_{j \in \{p, n\}} \int_{\Omega_j} \left( \alpha \partial_y s_g \overline{\partial_y v_j} + (-\partial_x(\alpha \partial_x s_g) - \omega^2 s_g) \overline{v_j} \right) d\mathbf{x} + \int_{\Gamma_j} (\alpha \partial_n s_g + i\lambda s_g) \overline{v_j} ds. $$

**Proposition C.4.** Let $g, k \in H_{per}^1(\Sigma)$ and $\psi$ satisfy assumption C.1. It holds that

$$ \int_{\Omega} \alpha \left[ s_g \overline{\partial_x s_k} - \partial_x s_g \overline{s_k} \right] \partial_x \psi d\mathbf{x} = \overline{b_{sing}^{(1)}(k, s_g \psi)} - b_{sing}^{(1)}(g, s_k \psi) \\ + 2i\lambda \sum_{j \in \{p, n\}} \int_{\Gamma_j} s_g \overline{s_k} \psi ds + 2i\pi \int_{\Sigma} g(y) \overline{k(y)} \mathbf{r}(y) \psi(0, y) dy. $$

Applying the last proposition with $\psi = 1$ yields immediately the following counterpart of Green's third formula.

**Corollary C.5.** For each $g, k \in H_{per}^1(\Sigma)$,

$$ b_{sing}^{(1)}(g, s_k) - \overline{b_{sing}^{(1)}(k, s_g)} = 2i\pi(g, k)\mathbf{r} + 2i\lambda \sum_{j \in \{p, n\}} \int_{\Gamma_j} s_g \overline{s_k} ds. \quad (C.8) $$

The third and last case consists in taking $U$ regular and $V$ singular, namely $U$ belonging to a certain subspace of $H_{1/2}^1(\Omega_{p,n})$ and $V = s_k$. Let us introduce $\psi_\Sigma(x, y) = \psi(0, y)$ in $\bar{\Omega}$. According to assumption C.1, $\psi_\Sigma \in C^1(\bar{\Omega}; \mathbb{R})$ with $\partial_x \psi_\Sigma = \partial_y \psi_\Sigma = 0$, so $\psi_\Sigma$ is actually a constant.
**Proposition C.6.** Let $\mathbf{u} \in Q$, $g \in H_{per}^1(\Sigma)$ and $\ell \in Q'$ with $\text{supp} \ell \cap \Sigma = \emptyset$, be such that

$$b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{sing}^{(1)}(g, \mathbf{v}) = \ell(\mathbf{v}), \quad \forall \mathbf{v} \in Q. \qquad (C.9)$$

Let $k \in H_{per}^1(\Sigma)$ and $\psi$ satisfy assumption C.1. Then

$$\sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha [u_j \overline{\partial_x s_k} - \partial_x u_j \overline{s_k}] \partial_x \psi d\mathbf{x} = \overline{b_{sing}^{(1)}(k, \mathbf{u}\psi)} - \ell_{\infty}(s_k\psi) + b_{sing}^{(1)}(g, s_k\psi) \qquad (C.10) \\ + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} \psi ds - \psi_{\Sigma} \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma},$$

where the jump $[\mathbf{u}]_{\Sigma}$ is defined in the statement of lemma 5.1 and $\ell_{\infty}(s_k\psi)$ as in (5.3).

The proof of the previous proposition relies on the following technical lemma, most of the proof of which is left to the reader.

**Lemma C.7.** Let $\mathbf{u} \in Q$, $g \in H_{per}^1(\Sigma)$, and $\psi$ satisfy assumption C.1. Then, $(\psi - \psi_{\Sigma}) \partial_x u_j \in L^2(\Omega_j)$, and $(\psi - \psi_{\Sigma}) \partial_x s_g \in L^2(\Omega)$. As a consequence, $u \in L^2(\Omega)$, s.t. $u|_{\Omega_{p,n}} = u_{p,n}$, satisfies $(\psi - \psi_{\Sigma}) u \in H^1(\Omega)$. Moreover, the trace $\gamma_0 [(\psi - \psi_{\Sigma})u] = (\psi - \psi_{\Sigma})u|_{\Sigma}$ vanishes.

*Proof.* We will prove the statement about the trace of $h_p := (\psi - \psi_{\Sigma})u_p \in H^1(\Omega_p)$ only. With the standard density argument, it suffices to prove the result for $u_p \in C^{\infty}(\Omega_p)$. We start with the expression

$$h_p(0, y) = h_p(x, y) - \int_0^x \partial_s h_p(s, y) ds, \quad (x, y) \in \Omega_p.$$

Applying the Cauchy-Schwarz inequality in $\mathbb{R}^2$ and $L^2(\Omega_p)$ yields

$$|h_p(0, y)|^2 \leq 2|h_p(x, y)|^2 + 2 \left| \int_0^x \partial_s h_p(s, y) ds \right|^2 \\ \leq 2|h_p(x, y)|^2 + 2x \int_0^x |\partial_s h_p(s, y)|^2 ds.$$

Remark that a priori $(x, y) \mapsto h_p(x, y)/x \in L^2(\Omega_p)$. Integrating both sides of the above inequality in the strip $\Omega_p^\varepsilon = \{(x, y) \in \Omega_p : |x| < \varepsilon\}$, $\varepsilon \in (0, 1)$, allows to obtain the following inequality, where all terms in the right-hand side are finite:

$$\int_{\Sigma} |h_p(0, y)|^2 dy \leq \varepsilon^{-1} \left( 2 \int_{\Omega_p^{\varepsilon}} |h_p(x, y)|^2 d\mathbf{x} + \varepsilon^2 \|\nabla h_p\|_{L^2(\Omega_p^{\varepsilon})}^2 \right) \\ \leq \varepsilon \left( 2 \left\| \frac{h_p(x, y)}{x} \right\|_{L^2(\Omega_p^{\varepsilon})}^2 + \|\nabla h_p\|_{L^2(\Omega_p^{\varepsilon})}^2 \right).$$

The above is valid for all $\varepsilon > 0$, hence taking $\varepsilon \to 0$ in the above shows that $\|\gamma_0 h_p\|_{L^2(\Sigma)} = 0$. $\square$
Proof of proposition C.6. We start by using (C.3), with $\mathcal{O} = \Omega_{p,n}$:

$$ \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha [u_j \overline{\partial_x s_k} - \partial_x u_j \overline{s_k}] \partial_x \psi d\mathbf{x} = I_1 - I_2, \quad \text{with} $$

$$ I_1 = \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha \partial_x (u_j (\psi - \psi_{\Sigma})) \overline{\partial_x s_k} d\mathbf{x}, \quad \text{and} \quad I_2 = \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha \partial_x u_j \overline{\partial_x (s_k (\psi - \psi_{\Sigma}))} d\mathbf{x}. $$

Remark that the above two integrals are well-defined by lemma C.7. On one hand, integrating by parts in $\Omega_j$, $j = p, n$, and noting that, according to proposition C.6, $(\psi - \psi_{\Sigma}) u \in H^1(\Omega)$ with vanishing trace on $\Sigma$, and $\alpha \partial_n s_k|_{\Sigma} = k(y) \mathbf{r}(y)$ yields

$$ \begin{aligned} I_1 &= \sum_{j \in \{p,n\}} \int_{\Omega_j} u_j (\psi - \psi_{\Sigma}) \overline{(-\partial_x (\alpha \partial_x s_k))} d\mathbf{x} + \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j (\psi - \psi_{\Sigma}) \overline{\alpha \partial_n s_k} ds \\ &= \overline{b_{\text{sing}}^{(1)}(k, \mathbf{u}(\psi - \psi_{\Sigma}))} - \sum_{j \in \{p,n\}} \int_{\Omega_j} (\alpha \partial_y (u_j (\psi - \psi_{\Sigma})) \overline{\partial_y s_k} - \omega^2 (u_j (\psi - \psi_{\Sigma})) \overline{s_k}) d\mathbf{x} \\ &\quad + i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j (\psi - \psi_{\Sigma}) \overline{s_k} ds. \end{aligned} $$

On the other hand,

$$ I_2 = b_{\text{reg}}^{(1)}(\mathbf{u}, s_k(\psi - \psi_{\Sigma})) - \sum_{j \in \{p,n\}} \int_{\Omega_j} \left( \alpha \partial_y u_j \overline{\partial_y (s_k(\psi - \psi_{\Sigma}))} - \omega^2 u_j \overline{(s_k(\psi - \psi_{\Sigma}))} \right) d\mathbf{x} \\ - i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j (\psi - \psi_{\Sigma}) \overline{s_k} ds. $$

Hence, using $\partial_y \psi = 0$,

$$ \begin{aligned} \sum_{j \in \{p,n\}} \int_{\Omega_j} \alpha [u_j \overline{\partial_x s_k} - \partial_x u_j \overline{s_k}] \partial_x \psi d\mathbf{x} = & \overline{b_{\text{sing}}^{(1)}(k, \mathbf{u}(\psi - \psi_{\Sigma}))} - b_{\text{reg}}^{(1)}(\mathbf{u}, s_k(\psi - \psi_{\Sigma})) \\ & + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} (\psi - \psi_{\Sigma}) ds. \end{aligned} \quad (C.11) $$

Comparing the above with the statement of the proposition, it remains to rewrite $b_{\text{reg}}^{(1)}(\mathbf{u}, s_k(\psi - \psi_{\Sigma}))$, using the identity (C.9) and the fact that $\psi_{\Sigma}$ is constant:

$$ \begin{aligned} b_{\text{reg}}^{(1)}(\mathbf{u}, s_k(\psi - \psi_{\Sigma})) &= \ell(s_k(\psi - \psi_{\Sigma})) - b_{\text{sing}}^{(1)}(g, s_k(\psi - \psi_{\Sigma})) \\ &= \ell_{\infty}(s_k\psi) - b_{\text{sing}}^{(1)}(g, s_k\psi) - \psi_{\Sigma} \left(\ell_{\infty}(s_k) - b_{\text{sing}}^{(1)}(g, s_k)\right), \end{aligned} \quad (C.12) $$

where $\ell_{\infty}(s_k\psi)$ and $\ell_{\infty}(s_k)$ are well-defined because $\text{supp } \ell \cap \Sigma = \emptyset$. Notice that $b_{\text{sing}}^{(1)}(g, s_k\psi)$ and $b_{\text{sing}}^{(1)}(g, \mathbf{v})$ are also well-defined since $b_{\text{sing}}^{(1)}(g, \mathbf{v})$ is well-defined as soon as $\mathbf{v} \in L^2(\Omega)$, $\partial_y \mathbf{v} \in L^2(\Omega)$ and the trace of $\mathbf{u}$ on $\Gamma_j$ belongs to $L^2(\Gamma_j)$ for $j = p, n$. Recall that the jump $[\mathbf{u}]_{\Sigma}$ satisfies (5.5), namely

$$ \ell_{\infty}(s_k) - b_{\text{sing}}^{(1)}(g, s_k) - \overline{b_{\text{sing}}^{(1)}(k, \mathbf{u})} = 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} ds - \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma}. $$

Combining (C.11), (C.12) and the above identity results in the desired expression. $\square$

<figure><img src="image_10.png" /></figure>
**Remark C.8.** Let $\varphi$ be as in definition 3.2. With this particular regular function, the previous propositions are respectively summarized as, with $\mathbf{u}$, $\mathbf{v}$, $g$ and $k$ satisfying the assumptions of the corresponding propositions,

$$ (\text{prop. C.2}) \quad C_{\varphi}(\mathbf{u}, \mathbf{v}) = \overline{b_{\text{reg}}^{(1)}(\mathbf{v}, \mathbf{u}\varphi)} - b_{\text{reg}}^{(1)}(\mathbf{u}, \mathbf{v}\varphi), \qquad (\text{C.13}) $$

$$ (\text{prop. C.4}) \quad C_{\varphi}(s_g, s_k) = \overline{b_{\text{sing}}^{(1)}(k, s_g\varphi)} - b_{\text{sing}}^{(1)}(g, s_k\varphi) + 2i\pi(g, k)\mathbf{r}, \qquad (\text{C.14}) $$

$$ (\text{prop. C.6}) \quad C_{\varphi}(\mathbf{u}, s_k) = \overline{b_{\text{sing}}^{(1)}(k, \mathbf{u}\varphi)} - \ell_{\infty}(s_k\varphi) + b_{\text{sing}}^{(1)}(g, s_k\varphi) - \langle[\mathbf{u}]_{\Sigma}, k\rangle_{\Sigma}. \qquad (\text{C.15}) $$

## D Expressions of $a^{(1)}$

This section is dedicated to the study of $a^{(1)}$ on $\text{Ker } \mathbf{B}^{(1)}$. For the convenience of the reader, we recall the expression (5.11) of the form $a^{(1)}$, written using the respective definitions (5.10) and (4.3) of the forms $C_{\varphi}$ and $b_{\text{sing}}^{(1)}$:

$$ a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = C_{\varphi}(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l}) - b_{\text{sing}}^{(1)}(h, (\mathbf{v} + s_{k-l})\varphi) + \overline{b_{\text{sing}}^{(1)}(l, (\mathbf{u} + s_{g-h})\varphi)}. \qquad (\text{D.1}) $$

The following technical lemma allows to reexpress the form $a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l))$ for $(\mathbf{u}, g)$ satisfying the assumption of proposition C.6 and $(\mathbf{v}, k, l) \in \text{Ker } \mathbf{B}^{(1)}$.

**Lemma D.1.** Let $\mathbf{u} \in Q$, $g \in H_{\text{per}}^{1}(\Sigma)$ and $\ell \in Q'$ with $\text{supp } \ell \cap \Sigma = \emptyset$, be such that

$$ b_{\text{reg}}^{(1)}(\mathbf{u}, \mathbf{v}) + b_{\text{sing}}^{(1)}(g, \mathbf{v}) = \ell(\mathbf{v}), \quad \forall \mathbf{v} \in Q. \qquad (\text{D.2}) $$

Let $\varphi$ in the definition of $a^{(1)}$ satisfy additionally to definition 3.2 $\text{supp } \ell \cap \text{supp } \varphi = \emptyset$. Then, for any $(\mathbf{v}, k, l) \in \text{Ker } \mathbf{B}^{(1)}$, it holds that

$$ a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = 2i\pi(g-h, k-l)\mathbf{r} - \langle[\mathbf{u}]_{\Sigma}, k-l\rangle_{\Sigma} + \overline{\langle[\mathbf{v}]_{\Sigma}, g-h\rangle_{\Sigma}}. \qquad (\text{D.3}) $$

*Proof.* We start by developing the first term in the definition of $a^{(1)}$ given by (D.1). Our goal is to rewrite it in terms of the forms $b^{(1)}$ and $b_{\text{sing}}^{(1)}$. Since $(\mathbf{v}, k, l) \in \text{Ker } \mathbf{B}^{(1)}$, $\mathbf{u} = \mathbf{v}$, $g = k$ satisfy the assumption of proposition C.6 with $\ell = 0$. Then, the form $C_{\varphi}$ is skew-hermitian, which yields, together with (C.13), (C.14), (C.15) the following expression:

$$ \begin{aligned} C_{\varphi}(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l}) &= C_{\varphi}(\mathbf{u}, \mathbf{v}) + C_{\varphi}(\mathbf{u}, s_{k-l}) + C_{\varphi}(s_{g-h}, \mathbf{v}) + C_{\varphi}(s_{g-h}, s_{k-l}) \\ &= \overline{b_{\text{reg}}^{(1)}(\mathbf{v}, \mathbf{u}\varphi)} - b_{\text{reg}}^{(1)}(\mathbf{u}, \mathbf{v}\varphi) \\ &\quad + \overline{b_{\text{sing}}^{(1)}(k-l, \mathbf{u}\varphi)} - \ell(s_{k-l}\varphi) + b_{\text{sing}}^{(1)}(g, s_{k-l}\varphi) - \langle[\mathbf{u}]_{\Sigma}, k-l\rangle_{\Sigma} \\ &\quad - b_{\text{sing}}^{(1)}(g-h, \mathbf{v}\varphi) - \overline{b_{\text{sing}}^{(1)}(k, s_{g-h}\varphi)} + \overline{\langle[\mathbf{v}]_{\Sigma}, g-h\rangle_{\Sigma}} \\ &\quad + \overline{b_{\text{sing}}^{(1)}(k-l, s_{g-h}\varphi)} - b_{\text{sing}}^{(1)}(g-h, s_{k-l}\varphi) + 2i\pi(g-h, k-l)\mathbf{r}. \end{aligned} $$
Remark that in the above the term $\ell(s_{k-l}\varphi)$ is well-defined and vanishes, since $\text{supp} \ell \cap \text{supp} \varphi = 0$. Rearranging the terms in the above yields

$$ C_{\varphi}(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l}) = \overline{b_{reg}^{(1)}(\mathbf{v}, \mathbf{u}\varphi) + b_{sing}^{(1)}(k, \mathbf{u}\varphi)} - \left[ \begin{aligned} & b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}\varphi) + b_{sing}^{(1)}(g, \mathbf{v}\varphi) \\ & - \overline{b_{sing}^{(1)}(l, (\mathbf{u} + s_{g-h})\varphi)} + b_{sing}^{(1)}(h, (\mathbf{v} + s_{k-l})\varphi) \\ & - \overline{\langle[\mathbf{u}]_{\Sigma}, k-l\rangle_{\Sigma} + \langle[\mathbf{v}]_{\Sigma}, g-h\rangle_{\Sigma}} + 2i\pi(g-h, k-l)_{\mathbf{r}} \end{aligned} \right]. $$

Using the definition (4.4) of the form $b^{(1)}$, namely $b^{(1)} = b_{reg}^{(1)} + b_{sing}^{(1)}$, and the assumptions of the lemma on $(\mathbf{v}, k, l)$ and $(\mathbf{u}, g, h)$ we rewrite the above as follows:

$$ C_{\varphi}(\mathbf{u} + s_{g-h}, \mathbf{v} + s_{k-l}) = -\ell(\mathbf{v}\varphi) - \overline{b_{sing}^{(1)}(l, (\mathbf{u} + s_{g-h})\varphi)} + b_{sing}^{(1)}(h, (\mathbf{v} + s_{k-l})\varphi) \\ - \overline{\langle[\mathbf{u}]_{\Sigma}, k-l\rangle_{\Sigma} + \langle[\mathbf{v}]_{\Sigma}, g-h\rangle_{\Sigma}} + 2i\pi(g-h, k-l)_{\mathbf{r}}. $$

Again, since $\text{supp} \ell \cap \text{supp} \varphi = \emptyset$, the first term in the above vanishes. Plugging in the resulting expression into the definition (D.1) of $a^{(1)}$ yields the desired expression (D.3). $\square$

The results that follow lead to an alternative expression to $a^{(1)}$ on $\text{Ker } \mathbf{B}^{(1)} \times \text{Ker } \mathbf{B}^{(1)}$.

**Proposition D.2.** Let $(\mathbf{u}, g, h)$, $(\mathbf{v}, k, l) \in \text{Ker } \mathbf{B}^{(1)}$. Then we have the following identity:

$$ \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma} - \overline{\langle [\mathbf{v}]_{\Sigma}, g \rangle_{\Sigma}} = 2i\pi(g, k)_{\mathbf{r}} + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} (\mathbf{u} + s_g) \overline{(\mathbf{v} + s_k)} ds. \quad (D.4) $$

*Proof.* Let $J := \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma} - \overline{\langle [\mathbf{v}]_{\Sigma}, g \rangle_{\Sigma}}$. According to the jump formula (5.5), we have

$$ J = \overline{b_{sing}^{(1)}(k, \mathbf{u})} + b_{sing}^{(1)}(g, s_k) + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{s_k} ds \\ - b_{sing}^{(1)}(g, \mathbf{v}) - \overline{b_{sing}^{(1)}(k, s_g)} + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} s_g \overline{v_j} ds. $$

Making use of the fact that $b^{(1)}((\mathbf{u}, g, h), \mathbf{v}) = 0$, $b^{(1)}((\mathbf{v}, k, l), \mathbf{u}) = 0$ and using the definition (4.4) of $b^{(1)}$ yields

$$ J = \overline{b_{reg}^{(1)}(\mathbf{u}, \mathbf{v}) - b_{reg}^{(1)}(\mathbf{v}, \mathbf{u})} + \overline{b_{sing}^{(1)}(g, s_k) - b_{sing}^{(1)}(k, s_g)} + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} (u_j \overline{s_k} + s_g \overline{v_j}) ds. $$

From the definition (4.5) of $b_{reg}^{(1)}$, it follows that $J_1 = 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} u_j \overline{v_j} ds$. Applying (C.8) to reformulate $J_2$, we readily arrive at (D.4). $\square$

**Remark D.3.** The above proposition can easily be extended to any $(\mathbf{u}, g)$, $(\mathbf{v}, k) \in Q \times H_{per}^{1}(\Sigma)$ as in lemma 5.1. Indeed, reproducing the computations in the proof leads to:

$$ \langle [\mathbf{u}]_{\Sigma}, k \rangle_{\Sigma} - \overline{\langle [\mathbf{v}]_{\Sigma}, g \rangle_{\Sigma}} = 2i\pi(g, k)_{\mathbf{r}} + 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} (\mathbf{u} + s_g) \overline{(\mathbf{v} + s_k)} ds \\ + \overline{b_{reg}^{(1)}(\mathbf{v}, \mathbf{u} + s_g)} + \overline{b_{sing}^{(1)}(k, \mathbf{u} + s_g)} - b_{reg}^{(1)}(\mathbf{u}, \mathbf{v} + s_k) - b_{sing}^{(1)}(g, \mathbf{v} + s_k). \quad (D.5) $$
The above proposition yields immediately the following property.

**Corollary D.4.** Let $(\mathbf{u}, g, h) \in \text{Ker } \mathbb{B}^{(1)}$ with $g \neq 0$. Then $\text{Im} \langle [\mathbf{u}]_{\Sigma}, g \rangle_{\Sigma} > 0$.

*Proof.* It is a direct application of previous proposition with $\mathbf{v} = \mathbf{u}$ and $k = g$, so that

$$ \text{Im} \langle [\mathbf{u}]_{\Sigma}, g \rangle_{\Sigma} = \pi \|g\|_{\mathbf{r}}^{2} + \lambda \sum_{j \in \{p,n\}} \|\mathbf{u} + s_{g}\|_{L^{2}(\Gamma_{j})}^{2}. $$

<figure><img src="image_11.png" /></figure>

Finally, proposition D.2 and lemma D.1 allow us to prove the following result, the second part of which is a reformulation of proposition 23 from [13], now in $V^{(1)}$.

**Corollary D.5.** Let $(\mathbf{u}, g, h), (\mathbf{v}, k, l) \in \text{Ker } \mathbb{B}^{(1)}$. Then

$$ a^{(1)}((\mathbf{u}, g, h), (\mathbf{v}, k, l)) = 2i\pi(g-h,k-l)_{\mathbf{r}} - 2i\pi(g,k)_{\mathbf{r}} - 2i\lambda \sum_{j \in \{p,n\}} \int_{\Gamma_j} (\mathbf{u} + s_g)(\overline{\mathbf{v} + s_k})ds. $$

In particular, it holds that

$$ a^{(1)}((\mathbf{u}, g, 0), (\mathbf{u}, g, 0)) = -2i\lambda \sum_{j \in \{p,n\}} \|u_j + s_g\|_{L^2(\Gamma_j)}^2, \quad \text{and} \quad a^{(1)}((\mathbf{0}, 0, h), (\mathbf{0}, 0, h)) = 2i\pi \|h\|_{\mathbf{r}}^2. \qquad (D.6) $$