Antoine Petit
email: antoine.a.petit@inria.fr

Stéphane Cotin
email: stephane.cotin@inria.fr

Environment-aware non-rigid registration in surgery using physics-based simulation

This paper presents a system for capturing the deformations of soft objects undergoing elastic deformations and contacts with their environment, using image and point cloud data provided by an RGB-D sensor. We improve upon previous works by integrating environment constraints in the frame-by-frame registration process. The approach combines a physics-based elastic model of the considered objects, computed in real-time using an optimized Finite Element Method (FEM), which is driven by surface constraints on the objects. Additional forces, such as gravity are added. A case study in open surgery on the liver is here described. Yet in this case a major improvement in the accuracy of the registration is provided by the integration of anatomical shape constraints, which are naturally hidden from the RGB-D camera, and that we account for through a registration with the pre-operative CT data. With a comparative study, we demonstrate the relevance of our method in a real world application mimicking an open surgery scenario where the liver has to be tracked to provide an augmented reality view.

Introduction

The problem of capturing the deformations of non-rigid objects is an active research field in computer vision that finds prospective applications in robotic manipulation, motion capture or computer-aided medical interventions. It becomes a challenging issue when considering multiple (soft) objects in interactions in a considered environment. Providing an accurately and reliable mean of registering and tracking these objects is a crucial task in applications such as bin picking or grasping soft objects with a soft manipulator, or providing guidance with augmented reality during surgical interventions.

Although the case of rigid or articulated objects in close interaction has shown an increasing interest in the research community in the recent years, considering deformable objects remains an open problem. The parameter space becoming much larger than 3D rigid poses or articulation joint states, several challenges must be overcome to jointly estimate the deformations of objects interacting with their environment. A main issue concerns inter occlusions between objects and their deformations due to contacts. Tracking independently the different objects and not considering the environment would indeed result, due to inaccuracies and ambiguities in the visual or point cloud information, in inter-penetrations and physically unrealistic behaviours. Instead, a simultaneous registration method relying on modelling the contacts between the objects would inherently handle these issues and reach a reliable behavior. Similarly, the field of non-rigid registration of anatomical structures, typically required to bring preoperative knowledge into the operating room, has been a long time research topic in medical imaging. However, deforming registration but also real-time tracking, while dealing with important occlusions of the organ surface during surgery, make it a real challenge.

In this article we introduce a solution to integrate collision detection and response models into a rigid and non-rigid registration framework combining RGB-D data and physics-based modeling. Collisions are handled by hierarchical bounding volume to contact point detection methods and on penalty forces to compute the contact response. This work extends the approach proposed in [START_REF] Petit | Capturing deformations of interacting nonrigid objects using rgb-d data[END_REF], by handling deformable/rigid contacts, by accounting for gravity as another physical prior, and by integrating patient-specific knowledge about anatomical boundaries.

The study described here shows new experiments in an open surgery context, for which target phantom organs are tracked, considering the surrounding anatomical structures.

The paper is organized as follows: after a review of the related literature, along with the motivations for the proposed system in section 2, descriptions of the elastic and interaction models are provided in section 3 and 4. Section 5 introduces the multiple objects registration framework and some experimental results are reported in section 6 to illustrate and validate the approach.

Related works

The goal is here to continuously estimate, based on successive RGB-D frames, the deformations undergone by one or multiple objects in interaction, given their known initial geometry. When applied to open surgery, the problem can be reformulated as accurately estimating the shape of a target organ, undergoing deformations due to surgical manipulations, using only partial surface information provided by an RGB-D camera. In the literature, various approaches have been proposed to register deformable objects, using vision and/or range data, and they can be classified according to the underlying model of the considered object, the resolution method, and according to the extracted visual and geometric features.

Single object non-rigid registration

Methods in [START_REF] Jordt | Direct model-based tracking of 3d object deformations in depth and color video[END_REF][START_REF] Zollhöfer | Real-time non-rigid reconstruction using an rgb-d camera[END_REF] employ an RGB-D sensor to register the acquired point cloud to a surface mesh by minimizing an error function accounting for geometric or direct depth and color errors, and a stretching penalty function for the mesh. By means of a NURBS parametrization [START_REF] Jordt | Direct model-based tracking of 3d object deformations in depth and color video[END_REF] or an optimized GPU implementation [START_REF] Zollhöfer | Real-time non-rigid reconstruction using an rgb-d camera[END_REF], real-time performance can be achieved. Although these two systems have shown promising and impressive results, they are still limited to isometric or small elastic deformations, by means of regularization functions proportional to squared distances between nodes of the mesh. Here, we wish to model elasticity in a more physically realistic manner, in order to handle volumetric effects and larger strains.

Indeed, another formulation of the problem relies on physics-based deformable models to perform registration. Other sorts (such as non-linear elasticity) and magnitudes of deformations can then be handled, inferring more consistently shape and/or volumetric regularization. The solution can be determined by balancing internal and external forces or, equivalently, minimizing energy functions. Physics-based methods include discrete mass-spring-damper systems [START_REF] Elbrechter | Bi-manual robotic paper manipulation based on realtime marker tracking and physical modelling[END_REF][START_REF] Schulman | Tracking deformable objects with point clouds[END_REF], or more explicit approaches relying on the Finite Element Method (FEM), based on continuum mechanics. [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF] use a linear tetrahedral co-rotational FEM model, coping with larger elastic deformations, external forces being related to correspondences between tracked 3D feature points mapped to the 3D mesh by means of a stereo camera system. Similarly, in [START_REF] Adagolodjo | Markerbased registration for large deformations -application to open liver surgery[END_REF] the authors propose a markers-based method allowing for the tracking and the deformation of a preoperative model in real-time during the surgery. Markers are manually placed on the surface of the organ after opening the abdominal cavity, and tracked in real-time by a set of infrared cameras. Since the use of markers remains an important clinical constraint, alternative methods have proposed to rely on RGB-D data (see [START_REF] Modrzejewski | Soft-body registration of pre-operative 3d models to intra-operative RGBD partial body scans[END_REF] for instance).

Multiple objects and interaction capture

However, previous works do not generally handle interactions between the object of interest and its environment. As we will see later, this can significantly impair the accuracy of the registration, which is essential in medical applications. We can mention the recent system designed in [START_REF] Guo | Real-time geometry, albedo, and motion reconstruction using a single rgb-d camera[END_REF], which introduces a collision detection process in its non-rigid reconstruction pipeline, to deal with self-collisions. However this technique falls in the class of non-rigid surface reconstruction of the whole observed scene. Our method, as a template tracking technique, would instead explicitly consider different known deformable or rigid interacting entities, and this in a volumetric way.

In addition to the generation of occlusions, interactions and contacts also result in deformations on the considered objects. Therefore, we advocate for the use of a physics-based approach to model both the deformations and the interactions, by relying on volumetric FEM-based models, and on collision detection models, as proposed in [START_REF] Petit | Capturing deformations of interacting nonrigid objects using rgb-d data[END_REF]. Based on successive RGB-D data inputs, the main contributions of this work consist in introducing interactions with known fixed rigid structures, while accounting for the impact of gravity. We also propose a way of integrating boundary conditions which are not visible by the RGB-D camera, in the specific context of open surgery. Let us note that the application considered in this work does not involve large rigid motions, so the process to handle these transformations presented in [START_REF] Petit | Capturing deformations of interacting nonrigid objects using rgb-d data[END_REF], through rigid ICPs, is not carried out here.

Overview of the system

As illustrated in Figure 1, our frame-by-frame tracking system can be outlined as follows: Input: the known 3D volumetric meshes of N O objects, a given RGB-D frame, and assuming correct registrations at the previous frame. A fixed rigid object can also be considered in the process.

1. Visually segment each of the structures of interest, using a parallel graph cut-based approach on each tracked object.

2. Derive linear external elastic forces between the current meshes and both the rigidly transformed meshes and the segmented point clouds, using closest point correspondences.

3. Numerically solve the mechanical equations required to compute the deformations on the objects, based on tetrahedral co-rotational FEM models, on collision models, and on gravity.

Elastic deformation model

We assume the structures of interest to have an homogeneous, isotropic elastic behavior. The volumetric displacement field is computed using a Finite Element Method. To handle large displacements under a small strain assumption, we use a co-rotational formulation [10]. It provides physical realism, by relying on continuum mechanics, while being computationally efficient.

Let us note that the modularity of the implemented framework would allow us to use other physical deformation models (hyperelastic, viscous, anisoptropic, ...).

For an exhaustive description of the FEM, the reader can refer to [START_REF] Cook | Finite element modeling for stress analysis[END_REF]. As a concise description, the method consists in tessellating the deformable object into a mesh made of elements connecting a set

X i = {x i,j } n X i
j=1 of 3D vertices. The deformation fields over the elements are approximated as continuous interpolations of the displacements of the vertices. We rely here on a volumetric linear FEM approach with tetrahedral elements. For elasticity modeling, we then resort to the linear elasticity theory, with Hooke's law, and to the infinitesimal strain theory [START_REF] Cook | Finite element modeling for stress analysis[END_REF], modified with the co-rotational approach so as to accommodate to rotational transformations and large displacements.

If the deformations ûe can be written as ûe = x e -x e,0 , we define here ûr e = R -1 e x e -x e,0 , with R -1 e x e the back rotated deformed coordinates of the four vertices of e, stacked into the 12 × 1 vector x e . R e is a 12 × 12 block diagonal matrix containing four copies of the 3 × 3 rotation matrix corresponding the rotational component of the deformations of the element.

The internal elastic forces f e exerted on the vertices of e can then be related to ûr e through:

f e = R e K e ûr e (1)
with K e the stiffness matrix of the element. K e depends on two elastic parameters of the material, the Young modulus E i and the Poisson ratio ν i , which are specific to each involved object.

Interaction model

Modeling the interactions between the considered objects constitutes a central issue to allow a joint registration of the deformations. We propose here, as in [START_REF] Petit | Capturing deformations of interacting nonrigid objects using rgb-d data[END_REF] a state-ofthe-art model of collision detection and collision response strategies. Collisions being performed on triangulated surface, we map a surface mesh of each object on its corresponding volumetric tetrahedral mesh by using the shape functions of the finite element method. These functions are linear and can be pre-computed as long as the object topology is not modified.

Collision detection

Collisions can be handled by employing different approaches to determine the intersections and the contacts between the interacting objects. Reviews of collision detection methods can be found in [START_REF] Jiménez | 3d collision detection: a survey[END_REF][START_REF] Teschner | Collision detection for deformable objects[END_REF]. Intersections are computed by performing successive distance tests, based on a hierarchical set of bounding volumes around the polyhedral triangular collision meshes.

Here, we use hierarchies of axis aligned bounding boxes (AABB) for each of the interacting structures. First, as a broad phase, intersections between the coarsest AABBs are found. Then, as a narrow phase, we descend in the AABB hierarchies of the objects in a coarse-to-fine manner. Within the finest AABBs, a search for intersections between pairs of geometric primitives is performed. The primitives are triangles in our case. The computation of intersections is based on a threshold on the distances between the primitives (triangle-to-triangle), and returns corresponding pairs of collision (or contact) points on the colliding primitives.

Collision response

In order to determine the contact forces between two colliding deformable objects i and j, we use a penalty-based method. For two detected collision points p and q on the collision surfaces of the interacting objects, barycentric mappings on the tetrahedral meshes are carried out to determine two corresponding points x i (p) and x j (q), in both volumetric models. Then a penetration distance δ can be defined as δ = d -D, with d = x i (p) -x j (q) and D a proximity threshold. As a collision normal n i,j we choose n i,j = p-q p-q . A spring-like repulsive penalty force f c j/i (x i ), exerted by x j on x i , can then be computed as:

f c j/i (x i ) = -kδn (2)
with k a stiffness factor. If non-moving (rigid) objects are present in the simulation, the response forces to the rigid/non-rigid contacts are only applied to the deformable body in order to optimize computation time.

Multiple objects registration

The frame-by-frame registration problem that we address consists in fitting the point cloud data, provided by an RGB-D sensor, with the tetrahedral meshes, in terms of nonrigid transformations. The method proposed in [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF] for single objects is here extended to multiple objects. Consequently, point cloud segmentation, matching, non-rigid registration techniques are adapted, as described below.

Let us note that the non-moving rigid object is registered thanks to a rigid ICP algorithm. Since the mounted RGB-D camera is also assumed to be static throughout the manipulation, the rigid transformation is kept fixed from the first frame. This would prevent from the introduction of collision detection issues with the non-rigid entities, due the jittering. Preliminary parallel visual segmentation.

The visual segmentation step presented in [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF], based on a graph cut minimization over foreground/background color models, is first performed. The objective is to restrict the acquired point clouds to the considered objects, so as to avoid ambiguities in the matching process, with the background or with occluding shapes. Dealing with multiple objects, the process is carried out in a parallel manner on each object O i , giving 3D target point clouds {Y i }. We limit the size of each Y i by sampling it on a regular grid in the image plane.

Parallel point cloud matching for non-rigid registration

With the aim of fitting the meshes to the point clouds in a non-rigid manner, similarly to [START_REF] Petit | Capturing deformations of interacting nonrigid objects using rgb-d data[END_REF], we compute external forces that are exerted by the observed segmented point clouds Y i on the different corresponding sets X V,i . We use sets of external forces f d ext,i related to geometrical information, as introduced in [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF]. The method consists first in determining nearest neighbors correspondences, both from the segmented point clouds to the meshes and from the meshes to the segmented point clouds. For each object, based on the two sets of mesh-to-point cloud and point cloud-to-mesh correspondences, we can compute an external elastic force f d ext,i exerted on each x i,j in X V,i as follows:

f d ext,i (x i,j ) = k d ext,i (x i,j -y f i,j ) (3)
where y f i,j is a linear combination of points in the point clouds which are matched to x i,j , either from mesh-to-point cloud and from point cloud-to-mesh correspondence sets. k d ext,i is the stiffness of these external elastic forces, which is here common to each object.

Computation of the solution

To compute the final shape of each object, we need to balance the sets of external forces capturing non-rigid transformations of the objects, the interactions between themselves, and gravity, with the internal forces associated with the different FEM models presented in Section 3. This requires to solve, for each object, a dynamic system of linear ordinary differential equations involving the internal and external forces, based on Lagrangian dynamics:

M i ẍi + C i ẋi + f i = f d ext,i + f c i + g (4)
with

f i = K i x i + f 0,i
where x i is a 3n Xi ×1 vector containing the positions to estimate of the vertices of each object O i . g is gravity, M i and C i are the 3n Xi × 3n Xi mass and damping matrices, K i the 3n Xi × 3n Xi N o global stiffness matrix which sums the 3n Xi × 3n Xi elementwise rotated stiffness matrices K r e = R e K e R -1 e , written with respect to whole set of vertices, and f 0,i the corresponding global offset summing the element-wise ones R e K e x e,0 . f d ext,i and f c i are 3n Xi × 1 vectors containing the external forces and resulting from the interaction model and the collision responses between object O i and the other non-rigid or rigid objects.

An Euler implicit integration scheme is used to obtain an an unconditionally stable time integration, thus allowing for larger time steps than explicit integration methods. To solve the final system of equations with respect to each x i , we use a conjugate gradient method. This iterative method allows for fast (and adjustable) computation times depending on the expected accuracy. We could also consider the static case, with the static equilibrium of the deformations assumed to be reached, so the transient and the dynamic terms of equation ( 4) would be neglected, leading to simply solving the equality between internal elastic forces and external forces with the conjugate gradient method.

Application to open surgery

The results presented here report a preliminary case study in open surgery. Our final objective is to perform augmented reality as a way to overlay internal structures of the organ (tumor, blood vessels) in order to guide the surgeon during the procedure. As a first step, we used two phantoms employed for surgical training and modeling the abdomen of the human body (liver, pancreas, spleen...) with silicon materials, contained in rigid plastic molds in the shape of a human torso (see fig. 2 and5).

Augmented reality for open surgery

For many lesions in the liver, open surgery remains the gold-standard approach, for which augmented reality (AR) appears as an useful tool to assist the surgeon, by improving the perception of the anatomy and avoiding missing lesions (disappeared metastasis after chemotherapy for instance).

Unlike neurosurgery or orthopedic surgery in which organs do not suffer from the deformation, real-time intraoperative AR is very rarely used for routine hepatic surgery.

The previous solution by [START_REF] Adagolodjo | Markerbased registration for large deformations -application to open liver surgery[END_REF] proposed an AR system for liver open surgery, by registering the deformations of the liver using a set of IR reflective markers set on the surface of the liver. The markers are tracked using an OptiTrack system and initially mapped on the surface of the known mesh modelling the liver. An underlying FEM physical model, similar to ours here, is integrated to provide an accurate and physically realistic estimation of the undergone deformations, through a constraints based approach. This technique has been promisingly validated during an in-vivo operation. However, the applicability of this system is limited by the delicate task of installing the IR reflective markers on the patient's liver, and the use of a cumbersome and expensive set-up such as the OptiTrack system around the operating table. This work proposes instead a marker-free approach, with the use of a single consumer RGB-D camera that could be easily mounted above the operating table. An Intel RealSense D435 is here employed, acquiring 640x480 aligned RGB and depth images. The camera is mounted to point vertically.

Patient-specific preoperative data

In the context of open surgery on the liver, the availability of the patient preoperative Computed Tomography (CT) scan can provide segmented 3D segmented of the different organs or tissues of interest. We proceed in the same way with the first phantom (in fig. 2). The 3D surface models of the liver (blue), the plastic mold (red) and the other tissues (greens) have thus been reconstructed from a CT scan. For the second phantom, an off-the-shelf RGB-D based dense 3D reconstruction techniques (InfiniTAM [START_REF] Kähler | Real-time large-scale dense 3d reconstruction with loop closure[END_REF]) was used. The corresponding volumetric tetrahedral meshes are then generated using a Delaunay triangulation, implemented in the CGAL library 2 . As an approximation, we assume the tissues to be isotropic in order to apply the deformation model described in Sect. 3. To implement the physical elastic and collision models, we have employed the Simulation Open Framework Architecture (SOFA) [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF]. Initialization is performed manually, shifted from the actual initial configuration.

A comparative study

Let us discuss the relevance of the approach by comparing the contribution provided by the collision models and gravity, on two scenarios carried out on the two considered phantoms.

In Fig. 2 are reported the results of different tested approaches on the first phantom: the liver only, without collision (third row), the liver with collision with the mold (fourth row), the liver and the registered surrounding tissues in the mold, with collision (fifth row), the sixth row shows the addition of gravity, and finally in the seventh row a lower stiffness for the liver is used. For the second phantom (Fig. 5), we consider the liver only (second row), the liver, the mold and the other tissues without (third row) and with (fourth row) collision.

The advantage of using collision for the mold can be clearly seen in Fig. 6 for the first phamton, where no collision (first row) leads to a large penetration of the liver into the tub. The detected collisions with the mold (second row) or both the mold and the other tissues (third row) are also featured with segments. The percentage of the penetrated volume of the liver phantom into the mold along the sequence is plotted in Fig. 3 (left). The effect of the various compression efforts on the surface can be observed with several peaks. The benefit of adding gravity to the model can be noticed in Fig. 2 (fifth row, second image), where the other modelled tissues tend to be repulsed from the mold when gravity is not active. This is due to the large occlusions undergone by these tissues and the poor visual segmentation of this area (second row), providing few point cloud based forces to fit it and to balance the penalty forces exerted by the mold when collision starts. The point cloud matching errors on the liver visible surface, which are based on nearest neighbor searches and which are used to compute the external forces, are reported in Fig. 3 (right). It shows that the addition of the collision models and gravity, which does not affect the surface error, and even slightly improves it when the other tissues are added. The effect of the lower stiffness for the liver can also be observed, with an expected lower error. It also results in a slightly better behavior around the compressed areas due to the hand pressure (Fig. 2). However, tracking these particular deformed areas is still quite inaccurate and remains a challenge. This is mostly due to the large occlusions by the hands. A more accurate non-linear physical model could also improve the behavior locally around these areas. The benefit of using collision models can also be observed for the second phantom (see Fig. 5), with a better accuracy when considering the other tissues and the mold (especially on the second column), and with penetrations when not considering collisions. More challenging color separation between the organs, occlusions and deformations, also show the robustness of the mehtod. The effect of gravity is here not noticeable. An issue lies in the computation time. For the first phantom for instance, considering the liver, the other tissues and the mold indeed results in a quite slow processing, around 5-6 fps. A more efficient pyramidal segmentation technique, lower mesh resolutions, collision models and point cloud resolutions could increase the framerate.

Conclusion

In this paper we propose a system able to simultaneously register interacting non-rigid and rigid objects, with a focus on an open surgery application. It extends a previous work, by showing its applicability with non-moving interacting rigid objects, and by adding gravity as another physical prior. Co-rotational FEM elastic models are used to accurately and efficiently model elasticity, along with collision detection and response models to cope with contacts between objects. Based on RGB-D input data, point cloud matching and registrations processes are achieved on segmented data. Deformations and contacts of the considered object are then estimated in a common solver, providing a realistic behaviour and coherently constraining registration on the involved entities, given the fixed boundary conditions. Preliminary results an abdominal phantoms show the benefit of introducing collisions and gravity. Future works would aim at improving the capture of local compression deformations, through specific 3D features, and at increasing the framerate. Further experimentations would consist in validation on invivo data.

Fig. 1 :

 1 Fig. 1: Proposed pipeline for multiple deformable objects tracking.

Fig. 2 :

 2 Fig.2: For the first phantom: visual segmentation of both target objects (the liver and the other tissues) on the second row. Registration of the liver (blue visuals) only (third row), colliding with the non-moving mold (red visuals) on the fourth row, colliding with registered other tissues (green visuals) on the fifth row, with the addition of gravity (sixth row), and finally with a lower stiffness (400 Pa against 600 Pa for the Young Modulus) on the last row.

Fig. 3 :

 3 Fig. 3: Penetrated volume of the liver model through the mold (left). Point cloud matching errors on the liver visible surface, for the different tested approaches. E0=600 Pa and E1=400 Pa correspond to two different tested Young Modulus for the liver (right).

Fig. 4 :

 4 Fig. 4: Effect of adding collisions for the first phantom. No collision detection (first row) results in inter-penetrations between the liver and the mold. The detected collisions can be observed with the mold (second row) and with both the mold and the other organs (third row).

Fig. 5 :

 5 Fig. 5: For the second phantom, registration of the liver (green visuals) only (second row), along with the registered other tissues (blue visuals) on the third row, and colliding with both the other tissues (blue visuals) and the mold (red visuals) on the fourth row.

Fig. 6 :

 6 Fig. 6: Effect of adding collisions for the second phantom. No collision detection (first row) results in inter-penetrations between the liver, the mold and the other tissues. The detected collisions can be observed with the mold and the other tissues (second row).

http://www.cgal.org


Anne-Sophie Bonnet-Ben Dhia

Lucas Chesnel
email: lucas.chesnel@inria.fr

Vincent Pagneux
email: vincent.pagneux@univ-lemans.fr

Trapped modes and reflectionless modes as eigenfunctions of the same spectral problem

Keywords: Waveguides, transmission resonances, backscattering, trapped modes, complex scaling, PT symmetry

We consider the reflection-transmission problem in a waveguide with obstacle. At certain frequencies, for some incident waves, intensity is perfectly transmitted and the reflected field decays exponentially at infinity. In this work, we show that such reflectionless modes can be characterized as eigenfunctions of an original non-selfadjoint spectral problem. In order to select ingoing waves on one side of the obstacle and outgoing waves on the other side, we use complex scalings (or Perfectly Matched Layers) with imaginary parts of different signs. We prove that the real eigenvalues of the obtained spectrum correspond either to trapped modes (or bound states in the continuum) or to reflectionless modes. Interestingly, complex eigenvalues also contain useful information on weak reflection cases. When the geometry has certain symmetries, the new spectral problem enters the class of PT-symmetric problems.

Introduction

In bounded domains, eigenvalues and corresponding eigenmodes completely determine the solutions to wave equations. In open systems, it is necessary to take into account also complex resonances, sometimes known as leaky modes or quasi normal modes in literature [START_REF] Aguilar | A class of analytic perturbations for one-body schrödinger hamiltonians[END_REF][START_REF] Balslev | Spectral properties of many-body Schrödinger operators with dilatation-analytic interactions[END_REF][START_REF] Simon | Quadratic form techniques and the Balslev-Combes theorem[END_REF][START_REF] Moiseyev | Quantum theory of resonances: calculating energies, widths and cross-sections by complex scaling[END_REF][START_REF] Aslanyan | Complex resonances in acoustic waveguides[END_REF][START_REF] Sjostrand | Complex scaling and the distribution of scattering poles[END_REF][START_REF] Zworski | Resonances in physics and geometry[END_REF][START_REF] Zworski | Lectures on scattering resonances[END_REF][START_REF] Moiseyev | Non-Hermitian quantum mechanics[END_REF]. Complex resonances provide the backbone of phenomena of wave propagation. For example, they indicate if there is a rapid variation of the scattering coefficients. We emphasize that real eigenvalues and complex resonances are intrinsic and universal objects that can be computed using different approaches. Among them, let us mention methods based on complex scaling or Perfectly Matched Layers (PMLs) [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] or techniques relying on Hardy space infinite elements [START_REF] Hohage | Hardy space infinite elements for scattering and resonance problems[END_REF][START_REF] Hohage | Convergence of infinite element methods for scalar waveguide problems[END_REF]. However complex resonances fail to provide useful information on recurrent questions such as "how large is the backscattered energy?" or "is there good transmission in the system?". Such questions are of tremendous importance in many topics of intense study in wave physics: extraordinary optical transmission [START_REF] Ebbesen | Extraordinary optical transmission through sub-wavelength hole arrays[END_REF], topological states immune to backscattering [START_REF] Wang | Observation of unidirectional backscattering-immune topological electromagnetic states[END_REF][START_REF] Lu | Topological photonics[END_REF], perfect transmission resonances [START_REF] Zhukovsky | Perfect transmission and highly asymmetric light localization in photonic multilayers[END_REF], transmission eigenchannels through disordered media [START_REF] Kim | Maximal energy transport through disordered media with the implementation of transmission eigenchannels[END_REF][START_REF] Sebbah | Scattering media: A channel of perfect transmission[END_REF], reflectionless metamaterials [START_REF] Rahm | Optical design of reflectionless complex media by finite embedded coordinate transformations[END_REF] or metasurfaces [START_REF] Yu | Flat optics with designer metasurfaces[END_REF][START_REF] Asadchy | Broadband reflectionless metasheets: frequency-selective transmission and perfect absorption[END_REF].

In this article, we define a new complex spectrum to quantify these aspects. In particular, our spectrum allows one to identify wavenumbers for which there is an incident field such that the backscattered field is evanescent (Figure 1). For construction of waveguides where such phenomenon occurs at a prescribed frequency, we refer the reader to [START_REF] Bonnet-Ben Dhia | Obstacles in acoustic waveguides becoming "invisible" at given frequencies[END_REF][START_REF] Chesnel | Team organization may help swarms of flies to become invisible in closed waveguides[END_REF][START_REF] Bonnet-Ben Dhia | Perfect transmission invisibility for waveguides with sound hard walls[END_REF][START_REF] Chesnel | Invisibility and perfect reflectivity in waveguides with finite length branches[END_REF]. Our spectrum also contains usual trapped modes [START_REF] Ursell | Trapping modes in the theory of surface waves[END_REF][START_REF] Evans | Trapped acoustic modes[END_REF][START_REF] Evans | Existence theorems for trapped modes[END_REF][START_REF] Davies | Trapped modes in acoustic waveguides[END_REF][START_REF] Linton | Embedded trapped modes in water waves and acoustics[END_REF][START_REF] Nazarov | Sufficient conditions on the existence of trapped modes in problems of the linear theory of surface waves[END_REF] called Bound States in the Continuum (BSCs or BICs) in quantum mechanics [START_REF] Sadreev | Bound states in the continuum in open quantum billiards with a variable shape[END_REF][START_REF] Moiseyev | Suppression of Feshbach resonance widths in two-dimensional waveguides and quantum dots: a lower bound for the number of bound states in the continuum[END_REF][START_REF] González | Bound states in the continuum in graphene quantum dot structures[END_REF][START_REF] Hsu | Bound states in the continuum[END_REF]. To compute this new spectrum, we use complex scaling/PMLs techniques in an original way. Our approach is based on the following basic observation: if for an incident wave, the backscattered field is evanescent, then the total field is ingoing in the input lead and outgoing in the output lead. Therefore, changing the sign of the imaginary part of the usual complex scaling/PML in the input lead should allow one to compute reflectionless modes (see the exact definition after [START_REF] Bonnet-Ben Dhia | Obstacles in acoustic waveguides becoming "invisible" at given frequencies[END_REF]). We show that this simple idea, inspired by the results obtained in [START_REF] Hernandez-Coronado | Perfect transmission scattering as a PT-symmetric spectral problem[END_REF] on a 1D PT-symmetric problem, indeed works. Interestingly, we will observe numerically that complex eigenvalues in this spectrum seem to contain useful information as well, indicating settings where reflection is weak.

u i u s u s + +
Figure 1: Schematic picture of a reflectionless mode. The propagating wave (blue) is not reflected. The backscattered field is purely evanescent (green).

To make the presentation as simple as possible, we stick to a simple 2D wave problem in a waveguide with Neumann boundary conditions. We assume that this waveguide contains a penetrable obstacle (Figure 2). Everything presented here can be generalized to other types of obstacles, to other kinds of boundary conditions (Dirichlet, ...), to higher dimension (3D) and to more complex geometries (one input lead/several output leads).

The paper is organized as follows. In section 2, the governing equations of the scattering problem are presented and the reflectionless case is defined. Section 3 is a reminder of the technique of complex scaling allowing one to compute classical complex resonances. In section 4, we explain how to use conjugated complex scalings to identify reflectionless frequencies. We provide examples of spectra in section 5 dedicated to numerics and we give some proofs of the theoretical results in section 6. Eventually, section 7 contains some concluding remarks. We are interested in the propagation of waves in a waveguide with a bounded penetrable obstacle (Figure 2). We assume that the waveguide coincides with the region Ω := {(x, y) ∈ R 2 | 0 < y < 1} and that the propagation of waves is governed by the Helmholtz equation with Neumann boundary conditions ∆u + k 2 γu = 0 in Ω ∂ y u = 0 on y = 0 and y = 1.

Setting

-L L γ = 5 γ = 1 (a) γ = 5 γ = 1 x y -L L (b)
(

) 1
Here u can represent for instance the acoustic pressure in a compressible fluid, or the transverse component of the displacement field in an elastic solid, or one component of the electromagnetic field in a dielectric medium. ∆ denotes the 2D Laplace operator. Moreover k is the wavenumber such that k = ω/c where c denotes the waves velocity and ω is the angular frequency, corresponding to a time harmonic dependence in e -iωt . We assume that γ is a positive and bounded function such that γ = 1 for |x| ≥ L where L > 0 is given. In other words, the obstacle is located in the region Ω

L := {(x, y) ∈ Ω | |x| < L}. Pick k ∈ (N π; (N + 1)π), with N ∈ N := {0, 1, . . . }. For n ∈ N, set β n = √ k 2 -n 2 π 2 .
In this work, the square root is such that if z = |z| e i arg(z) with arg(z) ∈ [0; 2π), then √ z = |z| e i arg(z)/2 . With this choice, we have m √ z ≥ 0 for all z ∈ C. The function

w ± n (x, y) = (2|β n |) -1/2 e ±iβnx ϕ n (y), (2)
with ϕ 0 (y) = 1, ϕ n (y) = √ 2 cos(nπy) for n ≥ 1, solves Problem (1) without obstacle (γ ≡ 1). For n = 0, . . . , N , the wave w ± n propagates along the (Ox) axis from ∓∞ to ±∞. On the other hand, for n > N , w ± n is exponentially growing at ∓∞ and exponentially decaying at ±∞. For n = 0, . . . , N , we consider the scattering of the wave w + n by the obstacle located in Ω. It is known that Problem (1) admits a solution u n = w + n + u s n with the outgoing scattered field u s n written as

u s n = +∞ p=0 s ± np w ± p for ± x ≥ L (3)
with (s ± np ) ∈ C N . The solution u n is uniquely defined if and only if Trapped Modes (TMs) do not exist at the wavenumber k. We remind the reader that trapped modes are non zero functions u ∈ L 2 (Ω) satisfying [START_REF] Aguilar | A class of analytic perturbations for one-body schrödinger hamiltonians[END_REF], where L 2 (Ω) is the space of square-integrable functions in Ω. We denote by K t the set of k 2 such that TMs exist at the wavenumber k. On the other hand, the scattering coefficients s ± np in (3) are always uniquely defined, including for k 2 ∈ K t . In the following, we will be particularly interested in the features of the reflection matrix (whose size, determined by the number of propagating modes, depends on k)

R(k) := (s - np ) 0≤n,p≤N ∈ C N +1×N +1 . ( 4
)
Definition 2.1. Let k ∈ (0; +∞) \ πN. We say that the wavenumber k is reflectionless if ker R(k) = {0}.

Let us explain this definition. In general, by linearity, for an incident field

u i = N n=0 a n w + n , (a n ) N n=0 ∈ C N +1 , (5)
Problem (1) admits a solution u such that u = u i + u s with

u s = +∞ p=0 b ± p w ± p for ± x ≥ L and b ± p = N n=0 a n s ± np ∈ C. (6)
The above definition says that, if k is reflectionless, then there are (a n ) N n=0 ∈ C N +1 \{0} such that the b - p in (6) satisfy b - p = 0, p = 0, . . . , N . In other words, the scattered field is exponentially decaying for x ≤ -L. Let us notice finally that the corresponding total field u = u i + u s decomposes as

u = N n=0 a n w + n + ũ for x ≤ -L u = N n=0 t n w + n + ũ for x ≥ L ( 7
)
where t n = a n + b + n and where ũ decays exponentially for ±x ≥ L. In other words, the total field is ingoing for x ≤ -L and outgoing for x ≥ L.

In the following, we call Reflectionless Modes (RMs) the functions u admitting expansion [START_REF] Bonnet-Ben Dhia | Non-scattering wavenumbers and far field invisibility for a finite set of incident/scattering directions[END_REF] and we denote by K r the set of k 2 such that the wavenumber k is reflectionless. The main objective of this article is to explain how to determine directly the set K r and the corresponding RMs by solving an eigenvalue problem, instead of computing the reflection matrix for all values of k.

Classical complex scaling

As a first step, in this section we remind briefly how to use a complex scaling to compute trapped modes. Define the unbounded operator A of L 2 (Ω) such that

Au = - 1 γ ∆u with Neumann boundary conditions ∂ y u = 0 on y = 0 and y = 1. It is known that A is a selfadjoint operator (L 2 (Ω) is endowed with the inner product (γ •, •) L 2 (Ω) ) whose spectrum σ(A)
coincides with [0; +∞). More precisely, we have σ ess (A) = [0; +∞) where σ ess (A) denotes the essential spectrum of A. By definition, σ ess (A) corresponds to the set of λ ∈ C for which there exists a so-called singular sequence (u (m) ), that is an orthonormal sequence (u

(m) ) ∈ L 2 (Ω) N such that ((A -λ)u (m)
) converges to 0 in L 2 (Ω). Besides, σ(A) may contain eigenvalues (at most a sequence accumulating at +∞) corresponding to TMs. In order to reveal these eigenvalues which are embedded in σ ess (A), one can use a complex change of variables. For 0 < θ < π/2, set η = e iθ and define the function I θ : R → C such that

I θ (x) = -L + (x + L) η for x ≤ -L x for |x| < L +L + (x -L) η for x ≥ L. (8)
For the sake of simplicity, we will use abusively the same notation I θ for the following map: {Ω → C × (0; 1), (x, y) → (I θ (x), y)}. Note that with this definition, the left inverse I -1 θ of I θ , acting from I θ (Ω) to Ω, is equal to I -θ . One can easily check that for all n ≥ 0, w + n • I θ is exponentially decaying for x ≥ L, while w - n • I θ is exponentially decaying for x ≤ -L. As a consequence, defining from expansion (6) the function v θ = u s • I θ , one has v θ = u s for |x| < L and v θ ∈ L 2 (Ω) (which is in general not true for u s ). Moreover v θ satisfies the following equation in Ω:

α θ ∂ ∂x α θ ∂v θ ∂x + ∂ 2 v θ ∂y 2 + k 2 γv θ = k 2 (1 -γ)u i (9)
with α θ (x) = 1 for |x| < L and α θ (x) = η -1 = η for ±x ≥ L. In particular, for a TM, v θ solves (9) with u i = 0. This leads us to consider the unbounded operator A θ of L 2 (Ω) such that

A θ v θ = - 1 γ α θ ∂ ∂x α θ ∂v θ ∂x + ∂ 2 v θ ∂y 2 (10)
again with homogeneous Neumann boundary conditions. Since α θ is complex valued, the operator A θ is no longer selfadjoint. However, we use the same definition as above for σ ess (A θ ), which is licit for this operator. We recall below the main spectral properties of A θ [START_REF] Simon | Resonances and complex scaling: a rigorous overview[END_REF]:

Theorem 3.1. i) There holds σ ess (A θ ) = n∈N, t≥0 {n 2 π 2 + te -2iθ }. ( 11
)
ii) The spectrum of 1) whose amplitude is exponentially growing at +∞ or at -∞.

A θ satisfies σ(A θ ) ⊂ R - θ with R - θ := {z ∈ C | -2θ ≤ arg(z) ≤ 0}. iii) σ(A θ ) \ σ ess (A θ ) is discrete and contains only eigenvalues of finite multiplicity. iv) Assume that k 2 ∈ σ(A θ ) \ σ ess (A θ ). Then k 2 is real if and only if k 2 ∈ K t . Moreover if v θ is an eigenfunction associated to k 2 such that m k 2 < 0, then v θ • I -θ is a solution of the original problem (
The interesting point is that now TMs correspond to isolated eigenvalues of A θ , and as such, they can be computed numerically as illustrated below. Note that the elements k 2 of σ(A θ ) \ σ ess (A θ ) such that m k 2 < 0, if they exist, correspond to complex resonances (quasi normal modes).

Let us point out that the complex scaling is just a technique to reveal them. Indeed complex resonances are intrinsic objects defined as the poles of the meromorphic extension from {z ∈

C | m z > 0} to {z ∈ C | m z ≤ 0} of the operator valued map z → (∆ + zγ) -1 .
For more details, we refer the reader to [START_REF] Aslanyan | Complex resonances in acoustic waveguides[END_REF].

Conjugated complex scaling

Now, we show that replacing the classical complex scaling by an unusual conjugated complex scaling, and proceeding as in the previous section, we can define a new complex spectrum which contains the reflectionless values k 2 ∈ K r we are interested in. We define the map J θ : Ω → C × (0; 1) using the following complex change of variables

J θ (x) = -L + (x + L) η for x ≤ -L x for |x| < L +L + (x -L) η for x ≥ L, (12)
with again η = e iθ (0 < θ < π/2). Note the important difference in the definitions of J θ and I θ for x ≤ -L: η has been replaced by the conjugated parameter η to select the ingoing modes instead of the outgoing ones in accordance with [START_REF] Bonnet-Ben Dhia | Non-scattering wavenumbers and far field invisibility for a finite set of incident/scattering directions[END_REF]. Now, if u is a RM associated to k 2 ∈ K r , setting w θ = u • J θ , one has w θ = u for |x| < L and w θ ∈ L 2 (Ω) (which is not the case for u). The function w θ satisfies the following equation in Ω:

β θ ∂ ∂x β θ ∂w θ ∂x + ∂ 2 w θ ∂y 2 + k 2 γw θ = 0 ( 13
)
with β θ (x) = 1 for |x| < L, β θ (x) = η for x ≤ -L and β θ (x) = η for x ≥ L. This leads us to define the unbounded operator B θ of L 2 (Ω) such that

B θ w θ = - 1 γ β θ ∂ ∂x β θ ∂w θ ∂x + ∂ 2 w θ ∂y 2 (14)
with homogeneous Neumann boundary conditions. As A θ , the operator B θ is not selfadjoint. Its spectral properties are summarized in the following theorem which is proved in the last section of this article.

Theorem 4.1. i) There holds

σ ess (B θ ) = n∈N, t≥0 {n 2 π 2 + te -2iθ , n 2 π 2 + te +2iθ }. ( 15
)
ii) The spectrum of

B θ satisfies σ(B θ ) ⊂ R θ with R θ := {z ∈ C | -2θ ≤ arg(z) ≤ 2θ}. ( 16
) iii) Assume that k 2 ∈ σ(B θ ) \ σ ess (B θ ). Then k 2 is real if and only if k 2 ∈ K t ∪ K r . Moreover if w θ is an eigenfunction associated to k 2 such that ± m k 2 < 0, then w θ • J -θ is a solution of (1)
whose amplitude is exponentially growing at ±∞ and exponentially decaying at ∓∞.

The important result is that isolated real eigenvalues of B θ correspond precisely to TMs and RMs.

The following proposition provides a criterion to determine whether an eigenfunction associated to a real eigenvalue of B θ is a TM or a RM.

Proposition 4.1. Assume that (k 2 , w θ ) ∈ R×L 2 (Ω) is an eigenpair of B θ such that k ∈ (N π; (N + 1)π), N ∈ N. Set ρ(w θ ) = N n=0 1 0 w θ (-L, y)ϕ n (y) dy 2 ( 17
)
where ϕ n is defined in (2). If ρ(w θ ) = 0 then w θ • J -θ is a TM (k 2 ∈ K t ). If ρ(w θ ) > 0 then w θ • J -θ is a RM (k 2 ∈ K r ).
In this case, the incident field u i defined in [START_REF] Bender | Making sense of non-Hermitian Hamiltonians[END_REF] with

a n = 1 0 w θ (-L, y)ϕ n (y) dy, n = 0, . . . , N,
yields a scattered field which decays exponentially for x ≤ -L.

The next proposition tells that B θ satisfies the celebrated PT symmetry property when the obstacle is symmetric with respect to the (Oy) axis. This ensures in particular the stability of simple real eigenvalues, with respect to perturbations of the obstacle satisfying the same symmetry constraint.

Proposition 4.2. Assume that γ satisfies γ(x, y) = γ(-x, y) for all (x, y) ∈ Ω. Then the operator

B θ is PT-symmetric (PTB θ PT = B θ ) with Pϕ(x, y) = ϕ(-x, y), Tϕ(x, y) = ϕ(x, y) for ϕ ∈ L 2 (Ω). Therefore, we have σ(B θ ) = σ(B θ ).
The proof is straightforward observing that the β θ defined after (13) satisfies β θ (-x, y) = β θ (x, y).

Finally let us mention a specific difficulty which appears in the spectral analysis of B θ . While Theorem 3.1 guarantees that σ(A θ ) \ σ ess (A θ ) is discrete, we do not write such a statement for the operator B θ in Theorem 4.1. A major difference between both operators is that C \ σ ess (A θ ) is connected whereas C \ σ ess (B θ ) has a countably infinite number of connected components. As a consequence, to prove that σ(B θ ) \ σ ess (B θ ) is discrete using the Fredholm analytic theorem, it is necessary to find one λ such that B θ -λ is invertible in each of the components of C \ σ ess (B θ ).

In general, in presence of an obstacle, such a λ probably exists (proofs for certain classes of γ can be obtained working as in [START_REF] Bonnet-Ben Dhia | Non-scattering wavenumbers and far field invisibility for a finite set of incident/scattering directions[END_REF]). But for this problem, we can have surprising perturbation results. Thus, if there is no obstacle (γ ≡ 1 in Ω), then there holds σ(B θ ) = R θ (see ( 16)): all connected components of C \ σ ess (B θ ), except the one containing the complex half-plane e λ < 0, are filled with eigenvalues. To show this result, observe that for k 2 ∈ σ(B θ ) \ σ ess (B θ ), the function u • J θ , with u(x, y) = e ikx , is a non-zero element of ker B θ . Notice that this pathological property is also true when Ω contains a family of sound hard cracks (homogeneous Neumann boundary condition) parallel to the (Ox) axis.

5 Numerical experiments

Classical complex scaling: classical complex resonance modes

We first compute the spectrum of the operator A θ defined in (10) with a classical complex scaling (complex resonance spectrum). For the numerical experiments, we truncate the computational domain at some distance of the obstacle and use finite elements. This corresponds to the socalled Perfectly Matched Layers (PMLs) method. We refer the reader to [START_REF] Kalvin | Analysis of perfectly matched layer operators for acoustic scattering on manifolds with quasicylindrical ends[END_REF] for the numerical analysis of the error due to truncation of the waveguide and discretization. The setting is as follows. We take γ such that γ = 5 in O = (-1; 1) × (0.25; 0.75) and γ = 1 in Ω \ O (see Figure 2 (a)). In the definition of the maps I θ , α θ (see ( 8), ( 9)), we take θ = π/4 (so that η = e iπ/4 ) and L = 1. In practice, we use a P2 finite element method in the bounded domain Ω 12 = {(x, y) ∈ Ω | -12 < x < 12} with Dirichlet boundary condition at x = ±12. Finite element matrices are constructed with FreeFem++1 .

In Figure 3 and in the rest of the paper, we display the square root of the spectrum (k instead of k 2 ). The vertical marks on the real axis correspond to the thresholds (0, π, 2π, ...). In accordance with Theorem 3.1, we observe that σ(A θ ) is located in the region R - θ = {z ∈ C | -θ ≤ arg(z) ≤ 0}. Moreover, the discretisation of the essential spectrum σ ess (A θ ) defined in [START_REF] Chesnel | Invisibility and perfect reflectivity in waveguides with finite length branches[END_REF] and forming branches starting at the threshold points appears clearly. Note that a simple calculation shows that {n 2 π 2 + te -2iθ , t ≥ 0} is a half-line for n = 0 and a piece of hyperbola for n ≥ 1. This is precisely what we get. Eigenvalues located on the real axis correspond to trapped modes (k 2 ∈ K t ). In the chosen setting, which is symmetric with respect to the axis R × {0.5}, one can prove that trapped modes exist [START_REF] Evans | Existence theorems for trapped modes[END_REF]. On the other hand, the eigenvalues in the complex plane which are not the discretisation of the essential spectrum correspond to complex resonances. . The trapped modes are in red, the dashed lines represent the essential spectrum of A θ (see [START_REF] Chesnel | Invisibility and perfect reflectivity in waveguides with finite length branches[END_REF]). The picture on the right is a zoom-in of that on the left.

Conjugated complex scaling: reflectionless modes

Now we compute the spectrum of the operator B θ defined in ( 14) with a conjugated complex scaling. First, we use exactly the same symmetric setting (see Figure 2 (a)) as in the previous paragraph. In Figure 4, we display the square root of the spectrum σ(B θ ). Since γ satisfies γ(x, y) = γ(-x, y), according to Proposition 4.2 we know that B θ is PT-symmetric and that therefore its spectrum is stable by conjugation (σ(B θ ) = σ(B θ )). This is indeed what we obtain. Note that the mesh has been constructed so that PT-symmetry is preserved at the discrete level. PT-symmetry is an interesting property in our case because it guarantees that eigenvalues located close to the real axis which are isolated (no other eigenvalue in a vicinity) are real. Therefore, according to Theorem 4.1, they correspond to trapped modes or to reflectionless modes. Remark that, for the same geometry, the spectrum of B θ (Figure 4) contains more elements on the real axis than the spectrum of A θ (Figure 3): the additional elements (green points in Figure 4) correspond to reflectionless modes. In Figure 5 top, we represent the real part of eigenfunctions associated with seven real eigenvalues of B θ . To obtain these pictures, we take L = 4 in the definition of J θ in . The green and red dots represent respectively the reflectionless modes and the trapped modes computed in Figure 4. We indeed observe that R 00 (k) is null for reflectionless k.

(12) and we display only the restrictions of the eigenfunctions to

Ω L = {(x, y) ∈ Ω | -L < x < L}.
We recognize two trapped modes (images 3 and 5). The other modes are reflectionless modes. In Figure 5 bottom, we provide the value of the indicator function ρ defined in [START_REF] Hernandez-Coronado | Perfect transmission scattering as a PT-symmetric spectral problem[END_REF] for the seven eigenmodes. We have to mention that eigenmodes are normalized so that their L 2 norm is equal to one. The indicator function ρ offers a clear criterion to distinguish between trapped modes and reflectionless modes. Moreover, in order to inspect the scattering coefficient, we remark that for reflectionless modes associated with wavenumbers k smaller than π, the incident field u i in (5) decomposes only on the piston mode w + 0 (x, y) = e ikx / √ 2k (monomode regime). In this case the reflection matrix R(k) in ( 4) is nothing but the usual reflection coefficient. In Figure 6, we thus display the modulus of this coefficient R 00 (k) with respect to k ∈ (0.1; 3.1). As expected, we observe that R 00 vanishes for the values of k obtained in Figure 5 solving the spectral problem for B θ . Of course obtaining the curve k → |R 00 (k)| is relatively costly and it is precisely what we want to avoid by computing the reflectionless k as eigenvalues. Here it is simply a way to check our results. In Figure 7, we represent the modulus of reflectionless mode eigenfunctions of B θ associated with one real eigenvalue and two complex conjugated eigenvalues. We observe, and this is true in general, a symmetry with respect to the (Oy) axis for modes corresponding to real eigenvalues which disappears for complex ones. This is the so-called broken symmetry phenomenon which is well-known for PT-symmetric operators (see e.g. the review [START_REF] Bender | Making sense of non-Hermitian Hamiltonians[END_REF]). 8) for a coefficient γ which is not symmetric in x nor in y. More precisely, we take γ such that γ = 5 in O = (-1; 0] × (0.25; 0.5) ∪ [0; 1) × (0.25; 0.75) and γ = 1 in Ω \ O. We observe that the spectrum is no longer stable by conjugation (σ(B θ ) = σ(B θ )) since the operator B θ is not PT-symmetric, and there is no "help" for the eigenvalues to be real. However, a closer look shows the presence of eigenvalues close to the real axis, in particular for k ≈ 1.0 + 0.13i, k ≈ 1.9 + 0.005i, k ≈ 2.5 + 0.02i, k ≈ 2.8 + 0.08i and k ≈ 3.0 -0.008i. In Figure 9, we represent k → |R 00 (k)| for k ∈ (0.1; 3.1) where there is only one propagating mode in the leads. It is interesting to note that the above computed complex reflectionless modes (located close to the real axis) have an influence on this curve. More precisely, k → |R 00 (k)| attains minima for k ∈ (0.1; 3.1) close to the real part of these complex reflectionless modes. Therefore complex reflectionless modes also have significance for scattering at real frequncies. The dashed lines represent the essential spectrum of B θ (see [START_REF] Evans | Existence theorems for trapped modes[END_REF]). The spectrum is not stable by conjugation. The picture on the right is a zoom-in of that on the left.

Proofs

Finally we give the proofs of Theorem 4.1 and Proposition 4.1.

Proof of Theorem 4.1. i) First, let us explain how to show that for t ≥ 0 and n ∈ N, n 2 π 2 + te -2iθ belongs to σ ess (B θ ). Consider a smooth cut-off function χ defined on R such that χ(x) = 0 for |x| > L and

χ L 2 (R) = 1. For m ≥ 1, set v (m) + (x, y) = m -1/2 χ (x -m 2 L)/m e i √ tx ϕ n (y)
where ϕ n is defined in [START_REF] Asadchy | Broadband reflectionless metasheets: frequency-selective transmission and perfect absorption[END_REF]. One can check that (v

(m)
+ ) is a singular sequence for B θ at n 2 π 2 +te -2iθ (note that the supports of v for m = m do not overlap). On the other hand, we prove that (v

(m) -), with v (m) -(x, y) = v (m)
+ (-x, y), is a singular sequence for B θ at n 2 π 2 + te +2iθ . Summing up, we obtain ∪ n∈N, t≥0 {n 2 π 2 + te -2iθ , n 2 π 2 + te +2iθ } ⊂ σ ess (B θ ). The converse inclusion requires a bit more work. Observe that if λ ∈ σ ess (B θ ), then one can prove by localization the existence of a corresponding singular sequence (v (m) ) supported either in x > L or in x < -L. In the first case, it means that (v (m) ) is also a singular sequence for A θ , so that λ ∈ σ ess (A θ ). In the second case, (v (m) ) is a singular sequence for the complex conjugate of A θ , which implies that λ ∈ σ ess (A θ ). Finally, the result follows from item i) of Theorem 3.1. iii) If k 2 ∈ K t ∪ K r , then k 2 is real and by construction k 2 ∈ σ(B θ ) \ σ ess (B θ ). Conversely, assume that k 2 is real and that k 2 ∈ σ(B θ ) \ σ ess (B θ ). There is one N ∈ N such that k ∈ (N π; (N + 1)π). Consider a non zero w θ ∈ ker (B θ -k 2 ). Setting u = w θ •J -θ , we find that u satisfies ∆u+k 2 γu = 0 in Ω and expands as

ii) Define the form b

θ (•, •) such that b θ (u, v) = Ω β θ ∂ x u ∂ x v + β -1 θ ∂ y u ∂ y v -k 2 β -1 θ u v dxdy. One can check that for k 2 ∈ C \ R θ ,
u = N n=0 a - n w + n + +∞ n=N +1 a - n w - n for x ≤ -L +∞ n=0 a + n w + n for x ≥ L, (18)
with (a ± n ) ∈ C N . If one of the a - n , for n = 0, . . . , N , is non zero, it means that u is a RM. If on the contrary a - n = 0 for all n = 0, . . . , N , then one can prove that a + n = 0 for all n = 0, . . . , N , so that u is a TM. Indeed, integrating by parts, we show that the quantity

F = Σ -L ∪Σ L (∂ n uu -u∂ n u) dy,
with Σ ±L = {±L} × (0; 1) and ∂ n = ±∂ x at x = ±L, satisfies F = 0. On the other hand, a direct calculation using expansion [START_REF] Hohage | Hardy space infinite elements for scattering and resonance problems[END_REF] and the orthonormality of the ϕ n yields F

= N n=0 i(|a + n | 2 - |a - n | 2
) and the result follows. Finally, consider some [START_REF] Hohage | Hardy space infinite elements for scattering and resonance problems[END_REF] for u = w θ • J -θ holds and one of the a - n , n = 0, . . . , N , has to be non zero. Indeed, otherwise u would be exponentially decaying for ±x ≥ L and k 2 would be in σ(A), which is impossible because σ(A) = [0; +∞). Therefore the amplitude of u is exponentially decaying at +∞ and exponentially growing at -∞.

k 2 ∈ σ(B θ ) \ σ ess (B θ ) such that m k 2 > 0 (the case m k 2 < 0 is similar). There is a unique N ∈ N such that 2π-2θ < arg(k 2 -(N π) 2 ) < 2π and π < arg(k 2 -((N +1)π) 2 ) < 2π -2θ. Then if w θ ∈ ker (B θ -k 2 ) \ {0}, expansion
Proof of Proposition 4.1. If (k 2 , w θ ) ∈ R × L 2 (Ω) is an eigenpair of B θ , then u = w θ • J -θ expands as in [START_REF] Hohage | Hardy space infinite elements for scattering and resonance problems[END_REF]. Moreover, we deduce from the orthogonality of the ϕ n that a - n = 1 0 w θ (-L, y)ϕ n (y) dy, n = 0, . . . , N, which gives the result, using the same arguments as in the proof of Theorem 4.1, item iii).

Concluding remarks

It is often desirable to determine frequencies for which a wave can be completely transmitted through a structure, a task usually leading to the tedious work of evaluating the scattering coefficients for each frequency. Here, we have shown that reflectionless frequencies can be directly computed as the eigenvalues of a non-selfadjoint operator B θ (see ( 14)) with conjugated complex scalings enforcing ingoing behaviour in the incident lead and outgoing behaviour in the other lead. The reflectionless spectrum of this operator B θ provides a complementary information to the one contained in the classical complex resonance spectrum associated with leaky modes which decompose only on outgoing waves (see the operator A θ in ( 10)). Note that eigenvalues corresponding to trapped modes belong to both the reflectionless spectrum and to the classical complex resonance spectrum because trapped modes do not excite propagating waves.

Moreover, since ingoing and outgoing complex scalings can be associated, respectively, with gain and loss, one observes that the non-selfadjoint operator B θ leads to consider a natural PT symmetric problem when the structure has mirror symmetry. Interestingly, a direct calculus shows that in the very simple case of a 1D transmission problem through a slab of constant index, reflectionless frequencies are all real. This gives an example of a non-selfadjoint PT symmetric operator with only real eigenvalues.

In this work, we investigated scattering problems in waveguides with N = 2 leads for which two reflectionless spectra exist: one associated with incident waves propagating from the left and another corresponding to incident waves propagating from the right. The more general case with N (N ≥ 2) leads can be considered as well. Among the total of 2 N different spectra with an ingoing or an outgoing complex scaling in each lead, two spectra correspond to eigenmodes which decompose on waves which are all outgoing or all ingoing. As a consequence, there are 2 N -2 reflectionless spectra.

Figure 2 :

 2 Figure 2: Symmetric (a) and non-symmetric (b) obstacles considered in the numerics below.

Figure 3 :

 3 Figure3: Classical complex resonances in the complex k plane corresponding to the spectrum of A θ for a symmetric obstacle (Figure2 (a)). The trapped modes are in red, the dashed lines represent the essential spectrum of A θ (see[START_REF] Chesnel | Invisibility and perfect reflectivity in waveguides with finite length branches[END_REF]). The picture on the right is a zoom-in of that on the left.

Figure 4 : 14 Figure 5 :Figure 6 :

 41456 Figure 4: Reflectionless eigenvalues in the complex k plane corresponding to the spectrum of B θ for a symmetric obstacle (Figure2 (a)). The trapped modes in red are the same as in Figure3. The reflectionless modes are in green and the dashed lines represent the essential spectrum of B θ (see[START_REF] Evans | Existence theorems for trapped modes[END_REF]). The picture on the right is a zoom-in of that on the left.

Figure 7 :

 7 Figure 7: Modulus of eigenfunctions of B θ associated to the eigenvalues k ≈ 5.31 (top), k ≈ 5.29 -0.13i (middle) and k ≈ 5.29 + 0.13i (bottom) obtained in Figure 4. The symmetry x → -x for real k (due to PT-symmetry) disappears for complex k. Now, we use the non-symmetric setting (see Figure 2 (b)), and we display the square root of the spectrum of B θ (in Figure 8) for a coefficient γ which is not symmetric in x nor in y. More

Figure 8 :

 8 Figure 8: Spectrum of B θ in the complex k plane for a non symmetric obstacle (Figure 2 (b)).The dashed lines represent the essential spectrum of B θ (see[START_REF] Evans | Existence theorems for trapped modes[END_REF]). The spectrum is not stable by conjugation. The picture on the right is a zoom-in of that on the left.

Figure 9 :

 9 Figure 9: Curve k → |R 00 (k)| (modulus of the reflection coefficient) for k ∈ (0.1; 3.1) and a non symmetric obstacle. The blue dots and the vertical dashed lines correspond to the real parts of the eigenvalues of B θ located close to the real axis computed in Figure 8. We observe that |R 00 (k)| is minimal for these particular k.

the numbers 1 ,|u| 2 +

 12 e ±iθ , -k 2 , -k 2 e ±iθ are all located in some region {z ∈ C | e (e iκ z) ≥ µ} for some κ ∈ [0; 2π) and µ > 0. We deduce that e (e iκ b θ (u, u)) ≥ µ Ω |∂ x u| 2 + |∂ y u| 2 dxdy which shows that b θ (•, •) is coercive on the Sobolev space H 1 (Ω). As a consequence of the Lax-Milgram theorem, B θ -k 2 is invertible for k 2 ∈ C \ R θ , which guarantees that σ(B θ ) ⊂ R θ .

FreeFem++, http://www.freefem.org/ff++/.