Gaëtan Lévêque
email: gaetan.leveque@univ-lille.fr

Yan Pennec

Pascal Szriftgiser

Alberto Amo

Alejandro Martínez

Scattering-matrix

come

protection takes place, in particular at sharp corners, which results in transmission minima and split-resonances, otherwise non-existent. In order to evaluate the small coupling coefficients associated to this breakdown, a phenomenological model based on an exact parameterization of scattering matrices at splitters and corners of the resonators is then introduced. By comparison with the numerical simulations, we are able to quantify the loss of topological protection at sharp bends and splitters. Finally, we use the obtained set of phenomenological parameters to compare the predictions of the phenomenological model with full numerical simulations for fractal-inspired cavities based on the Sierpiński triangle construction. We show that the agreement is overall good, but shows more differences for the cavity composed of the smallest triangles.

Our results suggest that even in a system exempt of geometrical and structural defects, topological protection is not complete at corners, sharp bends and splitters. However, simpler but predictive calculations can be realized with a phenomenological approach, allowing simulations of very large devices beyond the reach of standard simulation methods, which is crucial to design photonic devices which gather compactness and low losses through topological conduction of electromagnetic waves.

Topological photonics has recently become a disruptive paradigm enabling exotic ways to manipulate light propagation. [1][2][START_REF] Segev | Topological photonics: Where do we go from here?[END_REF][START_REF] Iwamoto | Recent progress in topological waveguides and nanocavities in a semiconductor photonic crystal platform[END_REF][START_REF] Price | Roadmap on topological photonics[END_REF] Amongst the different platforms to implement photonic structures relying on topological effects, two-dimensional (2D) high-index photonic crystal slabs display interesting features such as lossless propagation and large bandwidth, compatible with standard microfabrication tools. [START_REF] Barik | A topological quantum optics interface[END_REF][START_REF] Shalaev | Robust topologically protected transport in photonic crystals at telecommunication wavelengths[END_REF][START_REF] He | A silicon-on-insulator slab for topological valley transport[END_REF][START_REF] Parappurath | Direct observation of topological edge states in silicon photonic crystals: Spin, dispersion, and chiral routing[END_REF][START_REF] Arora | Breakdown of Spin-to-Helicity Locking at the Nanoscale in Topological Photonic Crystal Edge States[END_REF][START_REF] Barczyk | Interplay of Leakage Radiation and Protection in Topological Photonic Crystal Cavities[END_REF][START_REF] Kumar | Phototunable chip-scale topological photonics: 160 Gbps waveguide and demultiplexer for THz 6G communication[END_REF] An interesting proposal to build a topologicallyprotected waveguide in a 2D photonic crystal was presented in Wu et al. [START_REF] Wu | Scheme for Achieving a Topological Photonic Crystal by Using Dielectric Material[END_REF] Essentially, the idea is to design the unit cell of a honeycomb lattice so that it shows a Dirac cone at the Γ point at a given frequency. Then either by shrinking or expanding the motif inside unit cells, a topological band-gap arises. The interface between two semi-infinite shrunken and expanded lattices supports topologically protected modes showing a certain pseudo-spins for a given propagation directions. [START_REF] Wu | Scheme for Achieving a Topological Photonic Crystal by Using Dielectric Material[END_REF] Remarkably, when this approach is applied to 2D photonic crystal slabs, the guided modes are always over the light line, meaning that they are always radiative, a property that has been used to identify the pseudo-spin of the guided modes via far-field measurements. [START_REF] Parappurath | Direct observation of topological edge states in silicon photonic crystals: Spin, dispersion, and chiral routing[END_REF][START_REF] Arora | Breakdown of Spin-to-Helicity Locking at the Nanoscale in Topological Photonic Crystal Edge States[END_REF] The realization of large-scale photonic integrated circuits requires, however, waveguides that do not radiate. In photonic crystal slabs, this means that the guided modes should be below the light line to ensure perfect confinement by total internal reflection. The realization of topological waveguides supporting fully guided modes would require thus a honeycomb lattice showing -when undeformed -a Dirac point at symmetry points different to Γ in the first Brillouin zone. In contrast to the shrunken-expanded configuration, which mimics the spin Hall effect for photons, it has been proposed [START_REF] Ma | All-Si valley-Hall photonic topological insulator[END_REF][START_REF] Noh | Observation of Photonic Topological Valley Hall Edge States[END_REF] and experimentally realized [START_REF] Shalaev | Robust topologically protected transport in photonic crystals at telecommunication wavelengths[END_REF][START_REF] He | A silicon-on-insulator slab for topological valley transport[END_REF][START_REF] Wu | Direct observation of valley-polarized topological edge states in designer surface plasmon crystals[END_REF] photonic analogs of the valley Hall effect. Since experimental works use standard silicon technology, valley Hall photonic waveguides show a huge potential to become key elements in silicon photonics. One of the great advantage over light waveguiding along line defects in trivial photonic crystals is the ability, thanks to the topological protection, to conduct light even along sharp corners, with angles as small as 60 o . [START_REF] Shalaev | Robust topologically protected transport in photonic crystals at telecommunication wavelengths[END_REF][START_REF] He | A silicon-on-insulator slab for topological valley transport[END_REF][START_REF] Ma | Topological Photonic Integrated Circuits Based on Valley Kink States[END_REF][START_REF] Arora | Direct quantification of topological protection in symmetry-protected photonic edge states at telecom wavelengths[END_REF] This would allow for increasing the compacity and decreasing the footprint of future devices for information and communication technologies. However, and unlike in other topological photonic systems in which the time-reversal symmetry is broken (for instance by applying an external magnetic field [START_REF] Wang | Observation of unidirectional backscattering-immune topological electromagnetic states[END_REF][START_REF] Bahari | Nonreciprocal lasing in topological cavities of arbitrary geometries[END_REF] ), in photonic crystals the band structure is symmetric with respect to the wave vector along the propagation direction. This means that for any topologically-protected guided mode having a certain helicity relying its pseudo-spin and wave vector, there will be an identical state with opposite helicity. Albeit the topology of the system provides certain robustness to the propagation, 7,8 backward scattering is not prohibited by nature. Few works concern the assessment of the robustness of the topological protection, [START_REF] Arregui | Quantifying the Robustness of Topological Slow Light[END_REF][START_REF] Rosiek | Observation of strong backscattering in valley-Hall photonic topological interface modes[END_REF] and they essentially concern the back-reflection induced by inhomogeneities inside the photonic crystal itself, which as a result introduces propagation losses distributed over the whole length of the topological edge. To our knowledge, a quantitative evaluation of the pseudo-spin conversion in perfect topological crystals in presence of corners or other guiding elements like splitters is still lacking, and is a mandatory step to assess the applicability of such waveguides to more complex photonic systems, where those processes can be highly detrimental. This is particularly the case of ring-like cavities in which the finesse is highly sensitive to small losses due to the continuous recirculation of light.

In this work we analyze via numerical simulations the properties of valley topological edge-modes built in 2D photonic crystals. In particular, we focus on the loss of topological protection at sharp corners and splitters of triangular resonators coupled to a linear waveguide, [START_REF] Siroki | Topological photonics: From crystals to particles[END_REF][START_REF] Barik | Chiral quantum optics using a topological resonator[END_REF][START_REF] Jalali Mehrabad | Chiral topological photonics with an embedded quantum emitter[END_REF][START_REF] Devi | Topological edge states in an alldielectric terahertz photonic crystal[END_REF] either from a corner or the middle of an edge. In such systems, in absence of absorption, no resonant features are expected in the transmission spectra if perfect topological protection is realized: any deviation from a flat transmission band can in principle be traced back to a breakdown of topological protection somewhere along the path of light. We show that, in triangular resonators with different coupling conditions, transmission spectra present minima and split-resonances due to the coupling between counter-propagating waves, as in the case of other non-topological traveling-wave resonators. [START_REF] Zhang | Resonance-splitting and enhanced notch depth in SOI ring resonators with mutual mode coupling[END_REF] In order to elucidate and quantify the origin of the phenomenon, we then introduce a phenomenological model relying on the description of corners and splitters by scattering matrices whose exact expressions are derived. The numerical evaluation of the eight real parameters describing those matrices allow reproducing very precisely the simulations and assessing the small, frequency-dependent, coupling coefficients corresponding to pseudo-spin conversions at corners and splitters. Our results evidence that, even in the case of propagation of light in topological circuits free from geometrical and structural imperfections, reflection at sharp corners dominates the overall shape of transmission spectra, whose finer details are attributed to weaker ruptures of topological protection at the splitter. We conclude our study by demonstrating that our phenomenological model allows predictive and faster numerical simulations of complex circuit, taking the example of fractal-inspired resonators based on a Sierpiński triangle construction.

Numerical approach

Our topological photonic crystal (TPC), see Fig. 1(a), is based on the well-known hexagonallattice geometry made of circular holes with radii r 1 and r 2 , investigated for example in He et al 8 and other works. [START_REF] Barik | Chiral quantum optics using a topological resonator[END_REF][START_REF] Jalali Mehrabad | Chiral topological photonics with an embedded quantum emitter[END_REF] The study is restricted to TE (in-plane) polarization. For our bi-dimensional system, we chose a lattice constant b 0 = 385 nm, an average hole radius r 0 = 130 nm, and a refractive index n = 2.7. Those values allow matching the band-gap of the silicon membrane described in He et al. [START_REF] He | A silicon-on-insulator slab for topological valley transport[END_REF] Unless specified, all numerical simulations have been performed using the finite-elements-method (FEM) software Comsol Multiphysics.

When r 1 = r 2 = r 0 , the band diagram of the resulting honeycomb lattice presents Dirac cones at the six K points at the edge of the first Brillouin zone (FBZ), close to f 0 = 216 THz (Fig. 1(b), black line). A band-gap is then opened around f 0 for r 1 ̸ = r 2 (Fig. 1(b), red line), where 6-fold rotational symmetry of the lattice point group is lowered to 3-fold due to the breaking of inversion symmetry. For r 1 = 180 nm and r 2 = 80 nm, the band-gap corresponds to the frequency window [205.7 THz, 220.7 THz], indicated by the gray area.

The bands surrounding this bandgap present a non trivial local topology as expected from the valley Hall effect. This can be readily seen by computing the Berry curvature of those bands using a planewave expansion method [START_REF] Blanco De Paz | Tutorial: Computing Topological Invariants in 2D Photonic Crystals[END_REF] (see details in section 1 of the supplementary informations file). Figure 1(c) shows the calculated Berry curvature for the lowest band (lower right corner of the figure) and for the ensemble of bands 2 and 3 (upper right corner), which are touching. In both cases, the Berry curvature is concentrated at the K and K ′ points. For a given band, it has opposite sign at the K and K ′ , as expected for a time reversal symmetric system, and at a given K/K ′ point, each set of bands presents Berry curvatures of opposite sign. This configuration of opposite signs at opposite K/K ′ points and different bands is at the origin of the interface topological modes when two mirror symmetric photonic crystals are pasted together. This is confirmed by the non-zero values of valley Chern numbers, calculated by integrating the Berry curvature around the K and K ′ points: we obtain ∓0.090 for the first band and ±0.076 for the second and third bands.

The computed Chern numbers are low compared to the usually expected values of ±1/2.

Actually, as mentioned by several authors, [START_REF] He | A silicon-on-insulator slab for topological valley transport[END_REF][START_REF] Wong | Gapless unidirectional photonic transport using all-dielectric kagome lattices[END_REF] the Chern number is a half integer in the limit of weak perturbations, which corresponds here to small dissymmetries of the holes radii. In order to open an appreciable band-gap, r 1 must be significantly different of r 2 , which leads to an overlap of the Berry curvatures with opposite signs in each half of the 1BZ, finally resulting in lower Chern numbers. [START_REF] He | A silicon-on-insulator slab for topological valley transport[END_REF] We show on Fig. 1(d) the dispersion relation of the topological edge mode propagating along a bearded ΓK edge, in between two semi-infinite TPCs with glide mirror symmetry (in orange and blue on Fig. 1(e)). The system being now uni-dimensional with period b 0 , its FBZ is a segment with length 2π/b 0 along the ΓK direction (see inset), fully included inside the FBZ of the TPC. The breakdown of periodicity along the direction perpendicular to the edge direction induces a projection of the bulk modes of the infinite TPC onto the linear FBZ, which correspond to the light blue background, and of K and K ′ points onto the points indicated by the blue dotted lines. The shape of a typical unit cell used in the numerical simulation is plotted on Fig. 1(e), left panel, together with the distribution of the magnetic field amplitude for a frequency of 215.6 THz : the field is concentrated at the interface between the two crystals with a penetration length of about 1.5 unit cells into the bulk. The frequency of the topological mode has a local minimum at the Γ point, with a frequency of about 217.1 THz, lower than the bottom of the bulk band-gap. For that reason, the effective band-gap for the topological edge mode is [205.7 THz, 217.1 THz], underlined in dark gray on Fig. 1(d). From this simulation, we can extract the evolution of the wavevector or, equivalently, the effective index of the topological edge mode with frequency.

Triangular resonators

To investigate the robustness of the topological protection, we now characterize the properties of triangular edge mode resonators coupled to waveguides built on ΓK bearded edges. Such cavities can be coupled to a waveguide in two different ways, either from the edge, Fig. 2(a), or from the corner, Fig. 2(b). In both situations, the injection of the topological mode inside the resonator is realized through a splitter with four branches and labeled by S, while the corners of a triangular cavity will be noted C. As shown on Fig. 2(a) and (b), if topological protection is perfect in the system (which means at S and C points), all topological modes propagate along each edge with the same helicity, and then along the same direction. As a consequence, the system cannot reflect waves in the excitation guide, which means that the reflection coefficient R in power is zero, and through energy conservation the transmission The numerical simulations (see section 2 in the supplementary informations file) of tri-angular cavities with edge length L ≈ 28b 0 are presented on Fig. 2(c) to (f). Notice that the TPC containing the resonators and the coupling guide has been surrounded by Perfectly Matched Layers (PML's), whose role is to absorb the field along the outer edge of the simulation domain to simulate an infinite system. The spectra show, in contrast to the previous analysis, narrow transmission dips regularly separated in frequency, both for edge-and corner-addressed resonators, however with different profiles and frequencies. The simulation domain has been taken large enough to minimize the coupling of evanescent fields emanating from the structures (for example corners) with PMLs: we can then verify numerically that T ≈ 1 -R, as black (T ) and dashed-red (1 -R) curves overlap in the band-gap. In order to fully explain the transmission and reflection spectra, we need to suppose the existence of losses of topological protection in the system. For that reason, the topological edge mode can travel along the directions corresponding to the same helicity as for the excitation mode (blue arrows on Fig. 2(d) and (f)) or the opposite helicity indicated by the red arrows. This breakdown of topological protection can originate from the splitter or the triangle corners.

Both spectra present split resonances with low transmission (T ≈ 0, two and four respectively for the edge-and corner-addressed resonators inside the effective gap) and profiles characterized by a quasi-unity transmission (T ≈ 1, three and one respectively) in between two transmission minima above 0.5 in power, called anti-resonances below. Split resonances can be related to mode splitting arising from back-scattering in standard ring resonators. [START_REF] Zhang | Resonance-splitting and enhanced notch depth in SOI ring resonators with mutual mode coupling[END_REF][START_REF] Zhu | On-chip single nanoparticle detection and sizing by mode splitting in an ultrahigh-Q microresonator[END_REF] Typical distributions of the magnetic field amplitude are presented on Fig. 2(d) and(f) for split resonances and anti-resonances. For both transmission minima of the split resonance at F ≈ 209 THz (edge-addressed cavity), a clear interference pattern is obtained inside the cavity and in the excitation guide. The difference between both distributions is visible along the bisector plane crossing each corner: along those planes, the magnetic field is maximum or minimum respectively for the lowest and highest frequencies. The frequency difference is about 160 GHz, and the full width at half maximum (FWHM) is 100 GHz and 50 GHz respectively, corresponding to quality factors of Q = 2090 and 4180. At the anti-resonance (F ≈ 211 THz, blue arrow), the transmission is close to unity and no interference pattern is observed accordingly in the coupling waveguide, but a small intensity modulation is visible along the triangular cavity. Notice that the colorscale is the same for all distributions, except for the anti-resonances where the maximum value of the magnetic field amplitude is three times higher. Similar observations can be made on the field distributions of the corneraddressed cavity, see Fig. 2(f): the shape of the input waveguide does not modify the field distribution at the split-or anti-resonances, despite of their frequency shift as compared to the corner-addressed cavity.

If the occurrence of transmission split resonances and anti-resonances is a signature of a loss in topological protection, it is difficult from the numerical simulation to quantify this breakdown and find its origin. To this purpose, we propose in the next section a phenomenological approach, where both the splitter and triangle corners are described by a scattering matrix, allowing a simplified description of the systems.

Phenomenological model

Our phenomenological model of the triangular resonators relies on the description of the splitter S and corners C by a scattering matrix, which expresses a linear relation between outcoming waves with amplitudes A -, B -, C -... and incoming waves with amplitudes A + , B + , C + ... (see Fig. 3). The matrices elements can be partly extracted from numerical simulations. The topological nature of those modes implies that several of those elements are expected to be zero or much smaller than unity. The details of the calculations are given in the section 3 of the supplementary informations file, but we outline the main results below.

A simple corner, in the most general situation, behaves like a coupler between the two incident edge modes described by the complex vector

X + = [A + , B + ]
T with the corresponding transmitted and back-scattered modes, X -= [A -, B -] T , see Fig. 3(a). The shape of the scattering matrix, defined by the relation X -= M C X + , is constrained by energy conserva- tion, which implies that M C is unitary, and time-reversal symmetry, which, combined with unitary, implies that t 1 = t 2 . Finally, M C have the following general form:

M C =    t 1 r 1 r 2 t 2    = e iτ    cos σ i sin σ i sin σ cos σ   
The phase τ and angle σ can take arbitrary values. For convenience, we note in the following

a = t 1 = t 2 and b = r 1 = r 2 .
Concerning the splitter, four outputs,

X + = [A + , B + , C + , D + ] T , are now connected via the scattering matrix M S to four inputs, X -= [A -, B -, C -, D -]
T , see Fig. 3(b), with a priori 16 complex coefficients defined on the figures. However, the symmetry of the system implies that α 1 = α 2 and α 3 = α 4 , where α = r, t, ϵ, ϵ ′ , and we will note

t 1 = t 2 = t, t 3 = t 4 = t ′ , ϵ ′ 1 = ϵ ′ 2 = ϵ ′ , ϵ ′ 3 = ϵ ′ 4 = ϵ ′′ .
Energy conservation implies the unitarity of M S , and time reversal symmetry allows showing that all the r i and ϵ i coefficients are equal. Finally, additional symmetry considerations on the geometry of the splitter lead to the following parameterization of M S :

M S =          r ϵ ϵ ′′ t ′ ϵ r t ′ ϵ ′′ ϵ ′ t r ϵ t ϵ ′ ϵ r          = e iα 2 ×          (c ϕ e iρ + c ϕ ′ e -iρ ) (c ϕ e iρ -c ϕ ′ e -iρ ) i s ϕ e -iδ -s ϕ ′ e iδ e -iβ i s ϕ e -iδ + s ϕ ′ e iδ e -iβ (c ϕ e iρ -c ϕ ′ e -iρ ) (c ϕ e iρ + c ϕ ′ e -iρ ) i s ϕ e -iδ + s ϕ ′ e iδ e -iβ i s ϕ e -iδ -s ϕ ′ e iδ e -iβ i(s ϕ e -iγ -s ϕ ′ e iγ )e iβ i(s ϕ e -iγ + s ϕ ′ e iγ )e iβ (c ϕ e iρ + c ϕ ′ e -iρ ) (c ϕ e iρ -c ϕ ′ e -iρ ) i(s ϕ e -iγ + s ϕ ′ e iγ )e iβ i(s ϕ e -iγ -s ϕ ′ e iγ )e iβ (c ϕ e iρ -c ϕ ′ e -iρ ) (c ϕ e iρ + c ϕ ′ e -iρ )          (1)
with the additional constraint γ + δ + 2ρ = 0, and the definitions c x = cos x and s x = sin x.

The M S matrix is then parameterized by six free parameters : the two angles ϕ and ϕ ′ , and the four phases α, β, γ and δ.

We now need to numerically evaluate eight real parameters: two for M C and six for M S . To that purpose, reference points have to be defined in order to evaluate the phases of the matrices coefficients. Concerning the corner, the reference is taken as the intersection of the average lines of the bearded edges, as shown on Fig. 4 For the splitter, we have evaluated both the amplitudes and phases of the reflection (r) and transmission (t and t ′ ) coefficients defined in the semi-analytical model, see Eq. 1.

As indicated on Fig. 5(a), the four-branches splitter can be addressed either from one of the oblique (left panel) or the vertical (right panel) branches. Following the definition of the coupling coefficients and taking the example of the oblique excitation, the wave can be transmitted with a coefficient t ′ to the second oblique branch (called "transmission" branch), reflected with a coefficient r into a topological mode with same helicity propagating along the downward vertical branch ("reflection" branch), coupled to a mode with opposite helicity with a coefficient ϵ propagating along the upward vertical branch ("forbidden" branch), or back-scattered into the mode with opposite helicity into the excitation branch with coefficient ϵ ′′ . Similarly in the second case, and with equivalent denomination, the wave coming from a vertical branch can be coupled to the transmission branch (t), the reflection branch (downward oblique branch, r), the forbidden branch (upward branch, ϵ), or back-scattered (ϵ ′ ). As shown above, the coefficients r and ϵ must be the same in both cases (excitation from a vertical or an oblique branch), but the transmission and back-scattering coefficients can be different. The evolutions with frequency of the squared modulus of the transmission, reflection and forbidden transmission coefficients are plotted on Fig. 5(b), in solid lines for the vertical excitation, and with colored crosses for the oblique excitation. First, we can verify that the transmission and the reflection coefficients have nearly equal values in both configurations. Second, it appears clearly that the coefficient ϵ is weak (see green solid line).

For this reason, it cannot be evaluated by this method because the average value, below 1% in intensity, could be related to the reflection of the transmitted and reflected waves on the PMLs surrounding the simulation domain. As a consequence, |ϵ| 2 can have any value between 0 and about 1%. As numerically |t| 2 + |r| 2 ≈ 1, ϵ ′ and ϵ ′′ are confirmed to be, as ϵ, much smaller than 1. In that case, as the coefficient r is supposed to be the same in oblique or vertical excitation, energy conservation implies that |t| ≈ |t ′ |, which is correctly reproduced by the numerical simulation. This tends to show that topological protection is mostly conserved at the splitter. The dashed line is a second order fit of |r| 2 , which varies between 84 and 76%. The distribution of the magnetic field amplitude is plotted on Fig. 5(c) and (d) respectively for oblique and vertical excitation, at a frequency of 212 THz. We can visually verify the equality of the transmission and reflection coefficients. A very faint field can be distinguished along the forbidden channel, which may be again attributed to weak reflection of the transmitted and reflected fields on PMLs. As a conclusion, it is reasonable to consider as a first approximation that the splitter conserves topological protection: an incident topological mode can be either coupled to the transmission channel (characterized by coefficients t of t ′ ), the reflection channel (coefficient r), but not to the forbidden channel, neither being back-reflected. Those results agree well with those obtained by Ma et al [START_REF] Ma | Topological Photonic Integrated Circuits Based on Valley Kink States[END_REF] for a similar splitter but with triangular holes. In equation 1, it appears that perfect topological protection corresponds to ϕ = ϕ ′ and γ = δ = 0, which we suppose in the next paragraph.

In order to evaluate the phases of t, t ′ and r, reference points have been defined as indicated on Fig. 5(e): the phase will be for each coefficient the phase difference of the topological mode between two points linked by the corresponding arrow. Figure 5(f) shows the frequency evolution of the four phases ϕ t , ϕ t ′ , ϕ r and ϕ r ′ , normalized to π, corresponding to coefficients t, t ′ , r and r ′ . We can first verify numerically that, as predicted by energy conservation and time-inversion symmetry, ϕ r = ϕ r ′ , and finally r = r ′ . The second point is that the phase ϕ t for the transmission coefficient along the straight edge of the connection is equal, within the numerical uncertainties, to the propagation phase of the wave along the distance 3b 0 /2 between the two reference points of t. No additional phase is introduced by the presence of the nearby oblique edges. Next, the three phases ϕ t , ϕ t ′ and ϕ r are related to α and β (γ = ρ = 0) through:

α + β + π/2 = ϕ t α -β + π/2 = ϕ t ′ α = ϕ r which leads to π[2π] = 2ϕ r -ϕ t -ϕ t ′ .
We can see on Fig. 4(f) that this relation is very well verified numerically. As we suppose for now that the splitter preserves topological protection, the three remaining coefficients ϵ, ϵ ′ and ϵ ′′ are taken as 0.

Comparison with FEM simulations

The semi-analytical modelization of the triangular resonators is realized using a coupled wave approach, whose details are given in section 5 of the supplementary informations file.

The principle is summarized on Fig. S6 The transmission spectra can be computed for triangular resonators either addressed from the corner or the edge. For the resonators of Fig. 2, the spectra obtained by FEM simulations (Fig. 6(a,b)) compare very correctly with the semi-analytical model (bottom red lines) even if the splitter is supposed to perfectly preserve topological protection, Fig. 6(c,d). Note that Figure S7 shows the estimated frequency-evolution of the squared amplitude of the eight coefficients corresponding to the corner and splitter scattering matrices. As supposed, the largest coefficient at the origin of the pseudo-spin conversion is b, the reflection on a corner, but the coefficients ϵ, ϵ ′ and ϵ ′′ are finally comparable. For example, at F = 212 GHz, we have The semi-analytical model starts to differ more significantly from the numerical simulations Two main explanations can be proposed. First, as the interference paths are more complex, small errors in the estimations of the matrices coefficients have more significant impact on the transmission. Second, as the edges are shorter, we can expect that the model becomes less valid, as the edge-mode, which has a lateral extend, can directly tunnel laterally to adjacent edges across the lattice. An example of field distribution is shown on Fig. S8, for the resonance at F=214.61 GHz. Despite those limits, our results are encouraging and show that such a semi-analytical approach can be employed to predict, with low computer resources, the propagation of the electromagnetic signal along a complex and extended topological circuit.

Conclusion

In this work, we have investigated, using a full numerical approach based on the finite elements method and a semi-analytical model relying on an exact parameterization of scattering matrices at splitters and corners, the breakdown of topological protection evidenced by the resonance properties of valley topological triangular resonators. In resonators, perfect topological protection implies a flat transmission band, then any resonance feature must result from a back-scattering or forbidden transmission between waves with opposite helicity occurring at particular points of the system (corners, splitters...). In our system, we have demonstrated that the split-resonances, together with anti-resonances, must be mainly attributed to back-scattering at corners of the triangular cavity with a lower but comparable contribution of the splitter. Quantitatively, and for the considered valley topological crystal, the amount of power back-scattered at corners is about 1% of the incident edge-mode, while the backscattering and forbidden transmission at the splitter is lower than 0.5%. We have then demonstrated, by simulations of fractal-inspired larger resonators, that our semianalytical approach can be employed for fast and reliable predictive simulations of larger and more complex topological systems, if however the length of the edges composing the circuit are not too short (larger than 13 periods in our study) in order to avoid unwanted tunneling between close edges through the photonic lattice. We believe that the proposed methodology can be applied to different geometry of topological photonic devices in order to evaluate the quality of the topological protection depending on the shape (triangular vs circular) of air holes, edges (bearded or zigzag holes), or configuration of splitters (four or six branches) and corners, which is a crucial point in order to design photonic devices which gather compactness and low losses made possible through topological conduction of 20 electromagnetic waves.

Comsol simulations

Numerical simulations were realized using the finite elements software Comsol Multiphysics. 3 Semi-analytical model: coupling matrices

We derive in this section the analytical expressions of the scattering matrices between in and out-coming topological edge modes at either the corner of the linear waveguide (two inputs and two outputs) or a splitter (four inputs and four outputs). The calculation is based on energy conservation, time-reversal invariance and mirror symmetries of the configuration.

Scattering matrix of a corner

A simple corner, in the most general situation, behaves like a coupler between one of the two incident edge modes A + and B + with the corresponding transmitted and reflected modes,

A -and B -, see Fig. S3. The process is linear and obeys the following relation:

   A - B -    = M C    A + B +    =    t 1 r 2 r 1 t 2       A + B +   
Figure S3: Definition of the wave amplitudes and coupling coefficients for a corner.

Energy conservation implies that:

A -2 + B -2 = A + 2 + B + 2 , ∀A + , B +
The matrix M C is then unitary:

M † C M C = 1.
Additionally, time-reversal symmetry implies that :

   B + A +    * = M C    B - A -    * ⇔ S C    A + B +    * = M C S C    A - B -    * with S C =    0 1 1 0   
which, combined with unitarity, implies M C = S C M T C S C , and t 1 = t 2 :

M C =    t r 1 r 2 t   
The unitarity condition imposes:

|r 1 | 2 + |t| 2 = 1 |r 2 | 2 + |t| 2 = 1 r 1 t * + tr * 2 = 0 (1)
The two first equations impose

|r 1 | = |r 2 | = R. Let us write t = T exp (iτ ), r 1 = R exp (iρ 1 )
and r 2 = R exp (iρ 2 ). We have R 2 + T 2 = 1, and:

RT e i(ρ 1 -τ ) + e i(τ -ρ 2 ) = 0 =⇒ ρ 1 -τ = τ -ρ 2 + π + m2π, m ∈ Z or: ρ 1 -τ -π/2 + ρ 2 -τ -π/2 + m2π = 0
Writing:

r 1 = R exp (i(ρ 1 -τ -π/2 + τ + π/2)) r 2 = R exp (i(ρ 2 -τ -π/2 + τ + π/2)) = R exp (i(-(ρ 1 -τ -π/2) + τ + π/2)),
we pose ρ = ρ 1τπ/2 and we obtain r 1 = iR exp (i(ρ + τ )) and r 2 = iR exp (i(-ρ + τ )).

Finally, the M C matrix reads:

M C = e iτ    cos σ i sin σe iρ i sin σe -iρ cos σ    (2)
where we have introduced the angle σ such that T = cos σ and R = sin σ. We end up with three degrees of freedom. Note that the eigenvalues of M C are λ ± = exp{i(τ ± σ)}.

Assuming symmetry of the junction, we can as well impose r 1 = r 2 ⇒ ρ = 0, and, with σ ∈ [0, π] such as R = sin σ and T = cos σ:

M C = e iτ    cos σ i sin σ i sin σ cos σ   
For perfect topological protection, σ = 0 and M C is diagonal.

Scattering matrix of a splitter

Let us consider now the splitter. Following Fig. S4, four outputs, X + = [A + , B + , C + , D + ] T , are now connected to four inputs, X -= [A -, B -, C -, D -] T . We must then extend the 2 × 2 to a 4 × 4 matrix with the following shape:

X -= M S X + ⇔          A - B - C - D -          =          r 1 ϵ 2 ϵ ′ 3 t 4 ϵ 1 r 2 t 3 ϵ ′ 4 ϵ ′ 1 t 2 r 3 ϵ 4 t 1 ϵ ′ 2 ϵ 3 r 4                   A + B + C + D +         
The mirror-symmetry of the corner, with respect to the dashed line on Fig. 4, leads to

         C + D + A + B +          * = M S          C - D - A - B -          * ⇔ S S X + * = M S S S X - * , S S =          0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0         
Combined with unitarity of M S , that last relation leads to:

M S = S S M T S S S .
By comparison of both members of that equation, we finally see that r 1 = r 2 = r 3 = r 4 and

ϵ 1 = ϵ 2 = ϵ 3 = ϵ 4 .
The matrix M S can then be expressed under the following form:

M S =          r ϵ ϵ ′′ t ′ ϵ r t ′ ϵ ′′ ϵ ′ t r ϵ t ϵ ′ ϵ r         
As M S is unitary, it has four eigenvalues with unit modulus. We can exploit this property in order to find the full parameterization of the matrix. For that, we notice that the eigenvectors must be either symmetric or anti-symmetric with respect to the line of symmetry of the splitter (dotted line on Fig. 4).

A set of symmetric vectors is A 1 = [1, 1, 0, 0] T and A 2 = [0, 0, 1, 1] T , in which basis the matrix M S reduces to:

M A S =    r + ϵ t ′ + ϵ ′′ t + ϵ ′ r + ϵ   
As this matrix is unitary, it has the form given by equation 2, and :

r + ϵ = cos ϕe iα , t ′ + ϵ ′′ = i sin ϕe i(α+β) , t + ϵ ′ = i sin ϕe i(α-β)
Similarly, a set of anti-symmetric vectors is B 1 = [1, -1, 0, 0] T and B 2 = [0, 0, -1, 1] T , in which basis the matrix M S reduces to:

M B S =    r -ϵ t ′ -ϵ ′′ t -ϵ ′ r -ϵ   
which leads to:

r -ϵ = cos ϕ ′ e iα ′ , t ′ -ϵ ′′ = i sin ϕ ′ e i(α ′ +β ′ ) , t -ϵ ′ = i sin ϕ ′ e i(α ′ -β ′ )
Finally, we obtain:

r = 1 2 c ϕ e iα + c ϕ ′ e iα ′ ϵ = 1 2 c ϕ e iα -c ϕ ′ e iα ′ t ′ = i 2 s ϕ e i(α+β) + s ϕ ′ e i(α ′ +β ′ ) ϵ ′′ = i 2 s ϕ e i(α+β) -s ϕ ′ e i(α ′ +β ′ ) t = i 2 s ϕ e i(α-β) + s ϕ ′ e i(α ′ -β ′ ) ϵ ′ = i 2 s ϕ e i(α-β) -s ϕ ′ e i(α ′ -β ′ )
with c ϕ = cos ϕ, c ϕ ′ = cos ϕ ′ , s ϕ = sin ϕ and s ϕ ′ = sin ϕ ′ .

After some change of notation, we obtain:

M S = e iα 2 ×          (c ϕ e iρ + c ϕ ′ e -iρ ) (c ϕ e iρ -c ϕ ′ e -iρ ) i s ϕ e -iδ -s ϕ ′ e iδ e -iβ i s ϕ e -iδ + s ϕ ′ e iδ e -iβ (c ϕ e iρ -c ϕ ′ e -iρ ) (c ϕ e iρ + c ϕ ′ e -iρ ) i s ϕ e -iδ + s ϕ ′ e iδ e -iβ i s ϕ e -iδ -s ϕ ′ e iδ e -iβ i(s ϕ e -iγ -s ϕ ′ e iγ )e iβ i(s ϕ e -iγ + s ϕ ′ e iγ )e iβ (c ϕ e iρ + c ϕ ′ e -iρ ) (c ϕ e iρ -c ϕ ′ e -iρ ) i(s ϕ e -iγ + s ϕ ′ e iγ )e iβ i(s ϕ e -iγ -s ϕ ′ e iγ )e iβ (c ϕ e iρ -c ϕ ′ e -iρ ) (c ϕ e iρ + c ϕ ′ e -iρ )          (3)
with ρ = -(γ + δ)/2. There are then 6 degrees of freedom.

If the topological protection is perfect (ϵ = ϵ ′ = ϵ ′′ = 0), we have ϕ ′ = ϕ and γ = δ = ρ = 0, which gives:

M 0 S = e iα          cos ϕ 0 0 i sin ϕe -iβ 0 cos ϕ i sin ϕe -iβ 0 0 i sin ϕe iβ cos ϕ 0 i sin ϕe iβ 0 0 cos ϕ         
Finally, starting from the expressions of M A S and M B S , we show easily that the four eigenvalues of M S are:

λ ± A = e i(α+ρ±ϕ) λ ± B = e i(α+ρ±ϕ ′ )
superposition H c of the incident and reflected fields with amplitude H 0 reads:

H c (y) = H 0 [exp{ik(y -y 0 )} + b exp{-ik(y -y 0 )}] ⇒ |H c |(y) = |H 0 | 1 + |b| 2 + 2|b| cos {2k(y -y 0 ) -ϕ b }
where y 0 is the position of the corner and k the wavevector of the edge mode. We can verify that the contrast in amplitude is |b|. Numerically, the simulation data corresponding to In both cases, the incoupling field have the form:

X + =       A + B + Y +       A + = 1, B + = 0
where Y + contains the components of the input field inside the triangular cavity.

We introduce as well a matrix S which selects the cavity components Y + :

      0 0 Y +       = SX + ⇒ S ij = δ ij (1 -δ 11 )(1 -δ 22 )
where δ ij is the Kronecker symbol. We can then express:

X + = SX + + X 0 , X 0 = [1, 0, • • • , 0] T
The two relations between X + and X -are:

X -= CX + X + = PX -
The coupling matrix C is:

C =          M S M C 0 0 . . . M C         
where there are as many M C in the matrix as there are corners in the triangular cavity, not counting the splitter. Note that in the edge-connected cavity, the (t, ϵ ′ ) and (t ′ , ϵ ′′ ) coefficients must be exchanged as compared to the tip-connected cavity. The matrix P describes the phase relations between the field components inside he cavity.

Finally, we end up solving the following linear system of equations: For the tip-coupled cavity, we have the following definitions (see Fig. 6(a)):

X -= [1 -CSP]
X ± = A ± B ± C ± D ± E ± F ± G ± H ± T C =       M S 0 M C 0 M C       P =                      
0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 e iϕ 1 0 0 0 0 0 0 e iϕ 1 0 0 0 e iϕ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 e iϕ 0 0 0 0 0 e iϕ 0 0 0 0 e iϕ 1 0 0 0 0 0 0 0

                     
with ϕ i = nk 0 L i , where n is the effective index of the edge mode and k 0 = 2πF/c with F the frequency of the wave and c the speed of light in vacuum.

For the edge-coupled cavity (see Fig. S6(b)):

X ± = A ± B ± C ± D ± E ± F ± G ± H ± I ± J ± T C =          M ′ S 0 M C 0 0 M C 0 M C          P =                              0
                            
.

Figure 1 :

 1 Figure 1: (a) Geometry of the TPC and representation of the reciprocal space and associated high-symmetry points. (b) Dispersion diagram of the TPC for equal (black) or different (red) radii r 1 and r 2 of the two holes inside each primitive cell. (c) Berry curvature and valley Chern numbers simulated for the disymmetric TPC (r 1 = 180 nm and r 2 = 80 nm). (d) Dispersion curves (solid blue lines) of the edge modes propagating along the bearded interface in between two semi-infinite mirror-symmetric TPCs, parallel to the ΓK direction (the light blue background indicates the projected bulk modes). The solid red lines show the dispersion curves of the infinite TPC. The inset compares the FBZ of the interface (thick blue line with length 2π/b 0 ) and the FBZ of the infinite TPC. (e) Typical unit cell used in the simulation (left panel) and distribution of magnetic field amplitude of the edge mode (right panel).

Figure 2 :

 2 Figure 2: (a) and (b), respective representations of an edge-or corner-addressed topological triangular cavity. (c) Comparison between the transmission (solid black line) and reflection (dashed red line) spectra obtained from numerical simulations for the edge-addressed cavity. (d) Distribution of the magnetic field amplitude for frequencies indicated by a colored arrow on (c). For F = 211.42 THz, the maximum value is three times larger as for the two first frequencies. (e) and (f): same as (c) and (d) but for the corner-addressed cavity.

Figure 3 :

 3 Figure 3: Definition of the coupling coefficients for the corner, (a), and the splitter, (b), as employed in the scattering matrix semi-analytical model.

  (a). The transmission, a = exp{iτ } cos σ, and reflection, b = i exp{iτ } sin σ, coefficients have been numerically evaluated simulating the propagation of the edge mode along a single corner, whose magnetic field amplitude distribution is shown on Fig.4(b). The phase ϕ a = τ of the transmission coefficient a is obtained by comparing the phase of the transmitted mode to the phase of a mode propagating along a straight waveguide (see section 4 of the supplementary informations file).The numerical phase and a second order polynomial fit have been plotted on Fig.4(c). The amplitude |b| of the reflection coefficient is directly given by the contrast of the interference pattern in the injection guide, which is observable, despite being weak, on the magnetic field distribution of Fig.4(b). Figure4(d)shows the contrast as a function of the frequency, which presents a noticeable oscillation due to the interference with the wave which is slightly backscattered on the PML. However, we can obtain a correct approximation of |b| by evaluating the average contrast, on the order of 8-9% in amplitude. As ϕ b = ϕ a + π/2 and |a| = 1 -|b| 2 , the matrix M C is fully determined, within the simulation uncertainties.

Figure 4 :

 4 Figure 4: (a) Transmission, a, and reflection, b, coefficients on a corner. (b) Distribution of the magnetic field amplitude at F = 212 THz. (c) Transmission phase ϕ a normalized to π computed from numerical simulations. (d) Contrast of the interference pattern along the vertical input edge of figure (a), black solid line, and average contrast, solid red line.

Figure 5 :

 5 Figure 5: (a) Geometry of the splitter and coupling coefficients when addressed from an oblique (left) or a vertical (right) branch. (b) Evolution of the intensity of the waves propagating along each branch of the splitter. Solid lines (resp. crosses) correspond to the illumination along the vertical (resp. oblique) branch. (c) and (d) Distribution of the magnetic field amplitude respectively for oblique and vertical incidence. (e) Transmission and reflection coefficients across the splitter for vertical (t, r) and oblique (t ′ , r ′ ) incidence. (f) Phases of the different coefficients across the splitter.

  , which shows how C and S points are connected through segments of different lengths. Due to the choice of reference points for the splitter, the lengths of the different triangle edges are slightly different (L 0 = 28.5 b 0 , L 1 = 27 b 0 and L 2 = 13.5 b 0 ).

aFigure 6 :

 6 Figure 6: Comparison between the transmission spectra obtained from finite elements simulations and the semi-analytical model for the, (a,c,e), corner-addressed and, (b,d,f), edgeaddressed cavity. The splitter is either supposed to be perfectly topology-protected, (c,d), or to present small breakdown of topological protection, (e,f).

|a| 2 =Fig. 7 (

 27 Fig. 7(h): single and split-resonances are recovered at the same frequencies, but small differences in amplitude are observed for shallower resonances, close to 209, 210, 213 and 217 THz.

Figure 7 :

 7 Figure 7: Comparison between the full numerical simulation and the semi-analytical model for the three first iteration of the Sierpiński triangle construction, (a,b,c). (d,e,f) Corresponding transmission spectra computed with the finite-element-method. (g,h,i) Transmission spectra calculated with the semi-analytical model.

A

  Figure S2: Comsol configuration.

Figure S4 :

 S4 Figure S4: Definition of the wave amplitudes and coupling coefficients for a splitter.

Fig. S5(

  Fig.S5(c) have been fitted using a function f (y, r, s, t, u) = (r + s sin {ty + u}) 1/2 , leading to:

Figure S6 :

 S6 Figure S6: Simplified description of tip-and edge-coupled cavities in the semi-analytical model. The waves crossing the splitters and corners are connected by scattering matrices at nodes and propagation phases along the segments.

- 1

 1 CX 0 from which we extract the transmission and reflection coefficients T = |D -| 2 and R = |C -| 2 .

Acknowledgements This work was supported by the Horizon-RIA action project "Magnific" (101091968).

The phase ϕ a can then simply be extracted from the complex value of H z along the curved edge by comparison with the straight edge. 7 Field distribution in a fractal-type resonator


J D Béguin

D Samélor

C Vahlas

A Gleizes

J A Petit

B Sheldon

Evaluation of Al 2 O 3 MOCVD Coating for Titanium Alloys Protection Under Severe Conditions At High Temperature

Keywords: MOCVD, Al 2 O 3 coating, NaCl-induced corrosion, corrosion protection

The corrosion resistance of MOCVD Al 2 O 3 coating system was investigated to protect a TA6V Alloy under hot salt corrosion conditions: This coating was corroded with a salt deposit without mechanical loading at 480°C during 100 h. Corrosion products formed in salted areas were studied by Energy Dispersive Spectroscopy (EDS). Although all coated specimens were damaged with corrosion products presence in salted area, Al 2 O 3 coatings showed the lowest salt damage on titanium substrate after a metallographic cross section observation compared to uncoated ones. As well as these interesting experimental results, coated specimens exhibit a good adherence on titanium substrate

Introduction

Titanium alloys have been largely developed for aeronautic components like turboshaft engines, because of their good specific mechanical resistance, high tenacity and structural stability. However, under severe conditions in presence of aerated salted atmosphere and high temperature, corrosion damage can reduce the mechanical properties and lead to a premature failure. Over recent years, a variety of protective coatings has been developed for hot sections engines. It is well-known that the best salt corrosion barrier is an oxide scale such as alumina, silica or nickel oxide. There is 2 ways to obtain an oxide barrier: a deposition of an oxide forming coating or direct deposition of an oxide layer. These coatings must be thermomechanically and chemically compatible with the substrate, with respect to its microstructural, hence mechanical integrity. Previous studies showed the beneficial effect of Al 2 O 3 forming alloys for the hot corrosion resistance in chlorine-containing environements [1]. In this work, Al 2 O 3 coating was investigated to evaluate his protection performance for titanium alloys under hot salt corrosion. This alumina coating was applied by metal organic chemical vapor deposition (MOCVD) at 420°C. Such processing temperature allows maintaining the microstructural integrity of the Ti alloys. Deposited films are 180 nm thick and are transparent and amorphous. Their moderate thickness allows sustaining stresses originated from CTE mismatch between the coating and the substrate.

Experimental

MOCVD processing of the alumina coating

Deposition of the coating was performed from aluminium tri-isopropoxide under nitrogen flow in a horizontal, hot-wall reactor. Processing temperature and pressure were 420°C and 5 torr, respectively. Such processing temperature deposition was fixed not only in order to obtain the amorphous alumina phase but also to maintain the microstructural integrity of the substrate. Flat TA6V substrate samples were ground using a 600 grit abrasive SiC paper before surface treatment. This roughness value is similar to the roughness of aeronautical components in gas turbine engines. Deposition time of 2.5 min provided 180 nm thick films as was measured by reflectometry. Films are transparent, present a smooth appearance and do not contain any heteroatom, namely carbon.

coating properties

The surface morphology of the coating was examined with an AFM (Veeco Nanoscope IIIa ) in contact mode with a nitride silicon microlever under air atmosphere, and its structural homogeneity by lateral force in parallel (LFM). Because LFM is very sensitive to differences in the chemistry of the surface [3], it would indicate any differences in the chemistry between the coating and the substrate (Fig. 1). Although the initial surface roughness of the substrate, the alumina film seems to be present all over the substrate surface, on peaks and valleys of surface topography (Fig. 1A). A low friction variation on LFM image seems to confirm a chemical homogeneity of the surface (Fig. 1B).

Hot salt corrosion tests

Hot salt corrosion tests were performed on coated samples after NaCl deposit. The salt deposit was obtained by applying a calibrated droplet of saline solution on sample surface. NaCl solution was prepared with natural sea salt and deionised water. The salt crystallisation was achieved in an oven for 5 minutes at 115° C. This procedure allows obtaining a reproducible deposit. The deposited mass of NaCl was 120µg with a local salt concentration of 3 mg/cm 2 (Fig. 2). This local salt concentration is 10 times greater than that present in engines exposed to severe service conditions A B

Peaks Valley [4]. The salted sample was then placed in a furnace to 480° C during 100h. After the experiment, the degradation mechanism was observed by optical microscopy after immersion in deionised water to eliminate residual salt. In order to compare the effect of NaCl on coated sample, uncoated TA6V specimen was evaluated in the same conditions for reference. Fig. 2. Calibrated salt deposit after crystallisation on alumina coating.

Results and discussion

Hot salt corrosion behaviour of coated samples. The observation of corroded area for the latter sample by Scanning Electron Microscopy (SEM) indicates that pitting corrosion is the mode of degradation under the present experimental conditions (Fig. 4A). However it appears less severe compare to the uncoated sample (Fig. 5). At higher magnification, this attack seems to be due to local imperfections in the oxide film because of substrate Roughness (Ra=47 nm, Rt=300nm), and/or low thickness of oxide scale (Fig. 4A). In the of corroded area, the salt attack is not significant. These results are in agreement with previous studies [5]. It is well known that the presence of chloride ions leads to pitting attack, which generally initiates at imperfections in the oxide film. These chlorides penetration can facilitate the reaction with the substrate.

A B Corroded area

Residual corrosion products

In order to study the extend of salt attack, the corrosion depth was measured after a metallographic section (Fig. 5). The results clearly show that the depth of salt attack is higher on the uncoated sample than alumina coated sample (26 µm and 3 µm respectively).

Residual corrosion products Salt degradation

Degradation mechanism of alumina coating. The nature of corrosion products was revealed by Energy Dispersive X-ray spectroscopy (EDS) analysis and morphological examination (Fig. 6). The major chemical elements detected near the location pits were Ti (52 wt%), O (31 wt %), Al (7 wt%), Na and Cl less than 2 wt%. It is concluded that corrosion products are composed mainly of phases with the presence of remanent sea salt.

6. EDS analysis on corrosion products near a pit in salted area.

surface morphology of localised salt attack suggests that the salt deposit reacts with the titanium substrate to form corrosion products under the alumina coating. This volume change leads to its breakdown (Fig. 7). It's generally believed that the NaCl-induced hot corrosion is dominated by the chlorination/oxydation reactions [1,6]. Exposure of TA6V alloys to sea salt air atmosphere results in heavily corroded materials, with compromised mechanical performance. Application of protective coatings is a potential solution to meet this problem.

Fig. 1 .

 1 Fig. 1. Alumina film on TA6V substrate. (A) is a contact mode topology image, while (B) is the LFM image of the same area ("chemical image"). In LFM image, lower friction is darker shade.

Fig. 3

 3 presents optical surface micrographs of uncoated TA6V sample (3A) and (3B) Al 2 O 3 coated sample. A significant amount of small corrosion products (revealed in white) can be observed on the reference (uncoated) sample. In contrast, only a very small amount of such corrosion products are present on the corroded area of the Al 2 O 3 coated TA6V alloy.

Fig. 3 .

 3 Fig. 3. Optical microscopy of water washed samples after salt attack at 480°C/100h (A) uncoated TA6V sample (B) Al 2 O 3 Coated sample.

Fig. 4 .

 4 Fig. 4. SEM observations of corroded damage on alumina coating (A) low salt attack (B) salt attack on superficial defect (C) No visible salt attack at the center.

Fig. 5 .

 5 Fig. 5. Quantitative evaluation of corrosion damage after salt attack at 480°C/100h. TA6V SEM observation and (B) Al 2 O 3 SEM observation.

Fig. 7 .

 7 Fig. 7. Surface morphology on localised salt attack.

MOCVD alumina coatings on TA6V strongly improve the resistance to salt attack. The corrosion mechanism changes from uniform attack to pitting attack. Initiation of this localised corrosion is attributed to the imperfections of the coating. It is concluded that this is a promising coating for hot corrosion protection. Future work will aim the fine tuning of the MOCVD process in order to obtain adherent coatings with thickness equivalent to or higher than the surface roughness of the samples; i.e. 300 nm.


Martial Phélippé-Guinvarc'

Jean Cordier

Economic Viability of companion plants to lower the spread of Virus Yellows in Sugar Beets

Following EU's decision to ban neonics in 2018, several alternative solutions have emerged from the dedicated French national research program initiated in 2020. Companion plants, such as barley or oats, sown at the same time as beets, proved to have a repellent effect and reduce aphids' presence in sugar beet fields. This paper evaluates the economic viability of this particular solution and assesses how it reduces the risk of sugar beet infections by virus yellows. We show that although companion plants may be effective in repelling aphids, the strict economic impact associated with the reduction of aphids' population is mainly negative under cost-benefit assumptions currently adopted. Our analysis is based on the incidence model developed by Qi et al. It factors in costs induced by the planting and destruction of companion plants at a [4-6] leaf stage and factors in a 'competition' effect between companion plants and beets. As a consequence of our analysis, we highlight the need to run further research work in two directions, (1) analyse the potential of additional benefits offered by the solution to complement its strict impact on the number of aphids in sugar beet fields and (2) analyse potential benefits offered by synergies and combinations of companion plants with other viable solutions used to mitigate the spread of virus yellows in sugar beets.

Introduction

Following European Union's (EU) 1 decision to ban neonics in 2018, and given the catastrophic sugar beet campaign witnessed in 2020, the French government has implemented -since January 2021 -a Program for Research and Innovation (PNRI). This program aims at exploring alternative solutions to neonics for the sugarbeet sector. One notable approach involves sowing (concurrently with sugarbeets) companion plants like barley and oats in order to repel aphids and therefore mitigate their impacts on sugar beet production (aphids being the main vector that spreads yellow viruses in sugar beet fields). [START_REF] Jaworski | Crop diversification to promote arthropod pest management: A review[END_REF] explain that intercropping refers to the simultaneous or relay cultivation of multiple crops during a growing season within the same field. This includes companion plants, strip cropping, and relay intercropping. As for crop diversification strategies, it may include trap cropping, agroforestry, and landscape diversification. They show that recent literature has addressed the benefits of crop diversification for arthropod pest control, examining essential ecological mechanisms, tools, and scales for implementing diverse cropping systems.

As illustrated in Figure 1, [START_REF] Sarkar | Application of trap cropping as companion plants for the management of agricultural pests: A review[END_REF] elucidate the functioning of companion plants. In their comprehensive literature review, authors recognize that the development of an effective trap crop system is a costly process which often entails a secondary product but often with minimal market value. While biocontrol and natural repelling solutions may offer favourable outcomes, reaching, the desired results is typically a gradual and long-term process requiring much more time than the straightforward application of pesticides on trap crops (once pests have been trapped).

Our focus on companion plants is prompted by recent findings from 84 trials run by the French Technical Institute for Sugarbeet in collaboration with Agronomic Sugar Services from sugar companies as well as the participation of several agricultural colleges. This work has therefore paved the way for additional information regarding the feasibility and sustainability of implementing the 'companion plants' solution should it be used to mitigate the economic risk associated with virus yellows infections in sugar beet fields. Using trial results coming from the project mentioned above as well as economic results originating from the DEFT project run by the Association de Recherche Technique Betteravière (ARTB), as well as the model described in previous papers from [START_REF] Phélippé-Guinvarc'h | Actuarial Implications and Modeling of Yellow Virus on Sugar Beet After the EU's Ban on Neonicotinoids and Climate Change[END_REF], itself based on papers by Dewar and Qi (2021), by [START_REF] Qi | Decision making in controlling virus yellows of sugar beet in the UK[END_REF], by [START_REF] Werker | Modelling the incidence of virus yellows in sugar beet in the UK in relation to numbers of migrating myzus persicae[END_REF] and by [START_REF] Luquet | Predicting the seasonal flight activity of myzus persicae, the main aphid vector of virus yellows in sugar beet[END_REF] amongst others, a specific work has been performed on this solution.

While the efficiency of companion plants to repel aphids is acknowledged, our modelling work suggests that, based on current cost-benefit assumptions, there is hardly any economic incentive to use such solution to reduce aphids' presence in the fields. After defining the break-even loss above which beet losses are sufficiently big to make the use of companion plants economically viable, we show that the impact of a reduction in aphids' number is not linear. Consequently, a 23% decrease in aphids' number, as triggered by the companion plants, does not result in a 23% reduction of beet losses.

1 Two reasons for economic non-viability

Efficiency and costs of the solution

The primary effect to be factored in our study is the reduction in aphids' number triggered by companion plants which is expected to lower sugarbeet infections by virus yellows. The second aspect refers to the costs of such solution which includes planting and destruction of companion plants as well as the inclusion of an estimated beet yield loss resulting from competition between companion plants and sugar beets in the fields. The net cost-benefit value BC P C could then be written as follows :

BC P C = (ℓ -ℓ CP ) -(C ICP + C CCP )
where ℓ is the virus yellows beet loss without companion plants, ℓ CP is the virus yellows beet loss with companion plants, C ICP is the combined planting and destruction cost of companion plants, and C CCP is the reduction in beet yield loss due to competition between companion plants and sugar beets. Obviously, if there is no virus yellows beet loss, there is no companion plants' benefit. BC P C can be equivalently expressed in euros (e ) per hectare (ha) and as a percentage of the total expected beet revenue. Assuming a beet yield of 80 tons (at 16% sugar content) per hectare together with a price of e 45/ton of beets and considering an estimated average cost for the planting and destruction of companion plants amounting to e 97/ha (results from the project DEFT run by ARTB), the total expected beet revenue is e 3600/ha while the companion cost for planting and destruction (CICP) represents 2.7% of this expected beet revenue per ha).

As reported by [START_REF] Fabarez | Projet « fermes pilotes d'expérimentation » -rapport pour le comité de suivi scientifique[END_REF], pilot trials run at the farm allowed to test various solutions. Three crops: oats, spring barley and faba beans were tested as companion plants. Companion plants were tested with and without pesticides, paraffin oil, and Lecanicillium muscarium. Amongst those crops, oats showed the most promising results, especially when destroyed once sugar beets at a 6-leaf stage, with an average 23% reduction in aphids' number (noted ξ). In the meantime, beet yield losses due to competition with companion plants represented about 3% or e 108/ha (Figures 12 and9 of Fabarez ( 2024)). Using a straight assumption (further discussed and challenged in the paper) that 23% fewer aphids means a 23% reduction in beet yield losses,it is then possible to calculate the minimum beet loss level that equals the cost of using the companion plants (later refereed to as the break-even loss). For BC P C = 0 or ℓ -ℓ CP = C ICP + C CCP and total costs (aggregating competition effect costs and planting/destruction costs) of 108 + 97 = 205 e /ha, the break-even loss is e 891/ha or 24.8% of the expected beet revenue per hectare (e 205/ha being equivalent to 5.7% of the expected beet revenue per ha divided by 23%).

Run with MPI climatic data -RCP 8.5 scenario, our model results show that estimated beet losses never exceed the 24.8% break-even loss for the Polerovirus (BMYV). In other words and under our cost-benefit assumptions, the companion plants have no positive economic impact when sugar beets are infected by Polerovirus.

The model also determined that the 24.8% break-even loss threshold can be reached for BYV with a 1.3% probability of occurrence over the period (weighted by expected crop value).2 Although assuming a 23% reduction in winged aphids directly translates into a 23% reduction in beet losses is questionable, the rationale holds if we can establish that the rate of beet loss reduction rate, based on the aphids' reduction rate, can be encapsulated by a linear function. This conclusion becoming even more valid if the loss reduction rate is equal to or less than the rate of abundance reduction, as explained in the following sub-section.

The impact of a reduction in aphids' abundance 1.2.1 Main assumptions

This sub-section only recalls the elements of the article of [START_REF] Phélippé-Guinvarc'h | Actuarial Implications and Modeling of Yellow Virus on Sugar Beet After the EU's Ban on Neonicotinoids and Climate Change[END_REF] where the model of [START_REF] Werker | Modelling the incidence of virus yellows in sugar beet in the UK in relation to numbers of migrating myzus persicae[END_REF] is the reference. The observed proportion of plants exhibiting symptoms of virus yellows infections (denoted as Y ) increases over time and can be attributed to two sources of inoculum, namely primary infection (P ) and secondary infection (S). This increase is a result of the interplay between the proportion of disease-free plants (1 -Y ) and these two sources of infection.

Considering that the total proportion of diseased plants is given by Y = Y p + Y s , the combined rate of change of the total proportion of diseased plants with respect to time is given by

dY /dt = (r p P + r s Y )(1 -Y )
The primary inoculum, represented by the variable P , is defined as the number of migrant viruliferous aphids that lands on beets and feeds on it. The determination of P is based on the total number of aphids migrating in spring, which is assessed using suction traps from the Rothamsted Insect Survey.

The relationship between P and N , the total number of migrating aphids in spring, is set by the equation:

P = 1 -e -pN (1)
This equation corresponds to the multiple infection transformation and involves the constant of proportionality denoted as p. The use of this transformation is necessary to stabilize the values of P since the number of migrating aphids (N ) exhibits significant variations from one year to another.

In this investigation the infection cycle is initiated at time t 0 , when Y 0 = 0 and the solution is:

Y = 1 -e -(rpP +rs)(t-t 0 ) 1 + rs rpP e -(rpP +rs)(t-t 0 )
(2) Time t 0 coincides with the start of the migration of the aphids' (Myzus persicae) migration in spring, and at this point, the proportion of infected plants, denoted as Y 0 , is set to be zero. Last but not least, losses depend on the beet phenological stage when the virus is inoculated by aphids. The more mature the plant is, the more resistant it becomes and the lower beet losses will be. Furthermore and as highlighted in Figure 3, in the case of early infections, a plant infected by BMYV will suffer a loss of about 30%, whereas in the case of BYV the loss will be close to 50%. As a consequence, the whole beet field is not contaminated the same way and losses may differ from one place to another within the same field even if it is eventually infected and beets present similar yellow-leaf symptoms later in the campaign.

Abundance and first aphids' flights impact on incidence and losses

The equation 2 shows that incidence is a non-linear function of abundance and first aphids' flights. In addition, while the number of new infections increases, the rate of beet losses decreases as the plant matures and becomes more resistant. It is therefore hard to determine how a 23% reduction in aphids' number will exactly result when expressed in beet yield losses at the end. The aim of this section is to clarify the relationship between the reduction in aphids' numbers and the reduction in beet yield losses. From a methodological point of view and as stated on Figure 5, the number of days required by beets to reach a certain phenological stage are based on the number of average days calculated over the future 2025-2100 period. As recommended by ARTB (Phélippé-Guinvarc'h and Cordier, 2023), the assumption of Qi for the northern region is, to date, the most appropriate in the French context and was then used for our modeling work.

Figure 5 shows the relationship between beet losses and aphids' abundance with and without companion plants. The top-left graph shows that the beet losses increase with the aphids' number, and that the reduction in aphids' number due to the use of companion plants (dotted line) also reduces losses. The top-right graph shows the efficiency of the companion plants. What's most remarkable is that the reduction in incidence at the 6-leaf stage is stable regardless of abundance, while abundance increases as the crop progresses. Thus, although a decrease of 23% in aphids' abundance has an effect on beet yield losses close to 23% when abundance is close to 0, this effect decreases asymptotically to 0 at high abundances. In other words, the greater the loss due to aphids' abundance is, the less effective the companion plants are.

If we carry out the same exercise on the first aphids' flight in Figure 5 only day values between 109 and 128 are plotted. This has to do with the fact that if the first aphids' flight occurs earlier than that they cannot land on beets and if it occurs after that, companion plants have already been destroyed (the [6]leaf stage having been reached already). The bottom-left graph shows that the earlier first aphids' flight happens, the higher beet losses are. The bottom-right graph on the other hand, shows that companion plants are less effective when first aphids' flight happens early. The conclusion is therefore the same, the greater beet losses are due to early aphids' flights, the less effective companion plants are.

For the entire MPI scenario (IPCC 8.5), we calculated the beet losses with companion plants ℓ CP for BYV virus using abundance and stage days provided by the model, as detailed in Phélippé-Guinvarc'h and Cordier (2023) (450450 scenarios = 429 Agricultural Regions × 7 planting dates × 150 years). The results showed that the probability of ℓ CP < ℓ is equal to 0.

Discussion

Our cost-benefit analysis of companion plants as a way to contain the economic risk of virus yellows infections on sugar beets is largely negative under assumptions currently adopted and generally accepted. No circumstances of infections may justify the use of companion plants and result in a positive economic result when compared to do nothing. We found that the use of companion plants is too costly (5.7% of income per hectare) when compared to expected virus yellow beet losses (2.621% for BMYV and 6.711% for BYV) and to efficiency (23%). In addition, we have shown that significant losses occur when aphids are abundant and land early on the beets: a time at which the estimated efficiency of companion plants is the lowest.

Of course, these results are conditional to the quality of the Qi model and to the quality of the cost-benefit assumptions that have been used. As for the Qi model, its use as a reference for many projects of the French national research program surely calls for its use by the GREcoS project project Figure 5: At the top : the effect of abundance on the loss ℓ and the incidence Y . On the left, the figure shows that the beet losses that depends on abundance without (lined) and with (dashed) companion plants. On the right, the efficiency of reducing incidence (green) and losses (oranges) is shown. At the bottom: the effect of first aphid's flight day on the loss ℓ and the incidence Y . On the left, the figure shows the loss that depends on first aphids' flight day without (lined) and with (dashed) companion plants. On the right, the efficiency of reducing incidence (green) and losses (oranges) is shown.

(this publication arising from the work performed in this project) although some theoretical aspects of it may be questionable. Regarding the estimation of the full cost of companion plants, it has been set to an optimal level as we only considered a 3% beet yield losses due to competition between plants, the reality in the fields being slightly higher than that.

It is also worth noting that our results are in line with previous studies. In particular, [START_REF] Chauvin | The virus yellows epidemic in sugar beet in the UK in 2020 and the adverse effect of the EU ban on neonicotinoids on sugar beet production[END_REF] who used machine learning to identify risk factors for virus yellows spreading in sugar beet fields. It was surprising to find that the presence of aphids was only the 11th most significant factor (out of 154). The authors argue that since winter severity variables are excellent predictors of aphids' abundance, the role of the abundance variable is overshadowed by other factors in the statistical analysis. Although this argument is relevant, our analysis also highlights that the non-linear relationship between abundance, first aphid flight time and losses may also explain Chauvin and Soubeyrand's findings.

Meanwhile, one notable aspect of the work done by Fabarez is: the variability in companion plants' performance or efficiency. This is also the case for other alternatives to neonics proposed by Laurent et al.. Although introducing such variability could have allowed us to identify certain cases where the reduction in beet yield losses due to the companion plants was better than average, this would also have included 22% cases where companion plants would have been far less efficient than average. Integrating this variability would therefore have meant adding a risk premium to the full cost of the measure, thus increasing the break-even loss level of the companion plants.

Although the results may seem rather disappointing, companion plants can be considered in synergy with other solutions, for example by using odor repellent measures. Synergy between two solutions could indeed improve efficiency or generate significant positive externalities, such as soil enrichment or improved plant resistance to other risk factors. But while all the literature consulted unanimously stresses the need to combine technical solutions to combat the development of beet yellows, they should be subject to an economic and operational cost-benefit analysis. As Fabarez explained: 'It has been difficult to select and massively deploy trials that combine multiple levers (or solutions) because of the lack of efficiency from the levers alone.' At this point, these solutions have not shown tangible economic benefits (see for example [START_REF] Francis | Complementary strategies for biological control of aphids and related virus transmission in sugar beet to replace neonicotinoids[END_REF], [START_REF] Sarkar | Application of trap cropping as companion plants for the management of agricultural pests: A review[END_REF], [START_REF] Jaworski | Crop diversification to promote arthropod pest management: A review[END_REF]).

Conclusions

The originality of the article is to provide a cost-benefit analysis based on modelling work. Although the use of companion plants is almost never a viable economic solution within the framework of our cost-benefit assumptions and given its efficiency in reducing aphid numbers, this study work calls for further simulations to be done using other technical or combined solutions. These new and/or combined solutions may significantly reduce the number of aphids landing on sugar beets and transmitting virus yellows, or have other effects such as reducing implementation costs or providing other valuable positive externalities which could be valued economically and therefore change the economic equation of the solution.

Acknowledgment

Thanks to Luc Boucher, director of Risques Agricoles (DiagoRisk), to Azilis Lesteven and Alexis Patry, Association de Recherche Technique Betteravière (ARTB), for their valuable comments and suggestions.

Funding

The National Research and Innovation Program (PNRI) finances 23 projects and mobilizes many different players and institutions. This work was performed within the fourth axis of PNRI called 'Transition towards a sustainable economic model' and as a part of the Grecos project led by the Association de Recherche Technique Betteravière (ARTB). The GREcoS project aims at prefiguring what risk management system would be the most appropriate to cover beet yellows risk using alternative technical solutions to the neonics which are developed by other projects of the PNRI.

Figure 1 :

 1 Figure 1: From Sarkar et al. (2018). 'Role of trap crop to attract insect pests and natural enemies in a farming system.'

Figure 2 :

 2 Figure 2: BMYV (top) and BYV (bottom) beet loss density over 2025-2100 (weighted by expected crop value).

Figure 3 :

 3 Figure 3: Expected beet yield losses ℓ due to Polerovirus (dark orange) or BYV (orange) infections relative to this beet phenological stages. On this figure, the number of days required to reach the different phenological stages -emergency (blue), [6]-leaf, [12]-leaf and [18]-leaf (green) -have been determined in accordance with results obtained from simulations run over the (2025-2100) period.

Figure 4 :

 4 Figure 4: Abundance density 2025-2100 (weighted by expected crop value and where Abundance < 800)

Considering all oat trials, the efficiency of companion plants allowed a

36% reduction in aphids' number, but with an average beet yield loss of over 9%(Figures 2 and 9 of Fabarez (2024)). The costs of companion plants would therefore be e 421/ha (e 324+97/ha or 11.7% of the expected beet revenue per ha and the break-even loss would represent 1169/ha (or 32.48% of the expected beet revenue per ha). Assuming that in the future we can optimize the use of the companion plants to get a maximum efficiency of 36% and a minimal loss due to competition of 3%, the break-even of loss to reach may eventually be 15.82% and could be reached with a 10.8% probability of occurrence over the period.


Huaibin Tang
email: hutang@irisa.fr

Qinghua Zhang
email: zhang@irisa.fr

An Efficient Numerical Inverse Scattering Algorithm for Generalized Zakharov-Shabat Equations with Two Potential Functions

The inverse scattering problem for coupled wave equations has various applications such as waveguide filter design and electric transmission line fault diagnosis. In this paper, an efficient numerical algorithm is presented for solving the inverse scattering problem related to generalized Zakharov-Shabat equations with two potential functions. Inspired by the work of Xiao and Yashiro on Zakharov-Shabat equations with a single potential function, this new algorithm considerably improves the numerical efficiency of an algorithm proposed by Frangos and Jaggard, by transforming the original iterative algorithm to a one-shot algorithm.

Introduction

The inverse scattering problem for coupled wave equations [START_REF] Song | Inverse scattering problem for the coupled-wave equations when the reflection coefficient is a rational function[END_REF] has various applications, such as waveguide filter design [START_REF] Song | Design of corrugated waveguide filters by the gelfandlevitanmarchenko inverse-scattering method[END_REF][START_REF] Papachristos | Design of corrugated optical waveguide filters through a direct numerical solution of the coupled gelfandlevitanmarchenko integral equations[END_REF] and electric transmission line fault diagnosis [START_REF] Zhang | Inverse scattering for soft fault diagnosis in electric transmission lines[END_REF]. The solution to such problems can be computed, in principle, by solving the Gel'fand-Levitan-Marchenko (GLM) linear integral equations. In practice, numerical algorithms are developed to solve the GLM equations. After discretization, the GLM equations lead to a set of linear algebraic equations, typically with a large number of equations and unknowns. Though in principle such a large set of linear algebraic equations can be solved by standard algorithms with digital computers, in most practical situations, due to the large size of the equations, it is important to develop efficient numerical algorithms by taking into account the particularities of the inverse scattering problem. Such a numerical algorithm was proposed in [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF] for the inverse scattering problem related to the Zakharov-Shabat equations with a single potential function, then extended in [START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF] to the generalized Zakharov-Shabat equations with two potential functions. In this approach, by combining the GLM equations with some second order partial differential equations derived from Zakharov-Shabat equations, iterative algorithms are designed for solving discretized GLM equations. These iterative algorithms can be initialized by approximative solutions of discretized GLM equations. The result of [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF] is later improved in [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] where the second order partial differential equations combined to the GLM equations are replaced by first order partial differential equations. This modification has the advantage of triangularizing the set of linear algebraic equations resulting from the discretized differential and integral equations, hence the solution to the large set of linear algebraic equations becomes much more efficient.

The algorithm after the improvement of [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] is much faster than that of [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF]. However, in [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] this improvement has only been made to the the case of single-potential Zakharov-Shabat equations. In the present paper, the method of [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] is extended to the case of generalized Zakharov-Shabat equations with two-potential functions, resulting in a similar fast algorithm. The two-potential generalized Zakharov-Shabat equations have a larger scope of applications than the single-potential Zakharov-Shabat equations. For instance, when applied to electric transmission lines, the single-potential equations correspond to lossless lines, whereas the two-potential equations cover the case of lossy lines [START_REF] Jaulent | The inverse scattering problem for LCRG transmission lines[END_REF].

After a brief formulation of the inverse scattering problem for the generalized Zakharov-Shabat equations with two potential functions and a short presentation of their theoretic solution through GLM equations, the new numerical algorithm, extending the results of [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF][START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF][START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF], will be presented.

The inverse scattering problem

Consider the following generalized Zakharov-Shabat equations with two potential functions q + (x) and q -(x) which will be referred to as ZS + :

dν 1 (k, x) dx + ikν 1 (k, x) = q + (x)ν 2 (k, x) (1a) dν 2 (k, x) dx -ikν 2 (k, x) = q -(x)ν 1 (k, x) (1b)
ZS + :

    
where q ± (x) are the unknown coupling potential functions and k is the wave number. Denote with r l (k), r r (k) and t(k) respectively the left reflection coefficient, the right reflection coefficient and the transmission coefficient of ZS + , or of (1), the inverse scattering problem considered in this paper is to determine the coupling potential functions q ± (x) from the scattering data. Assume that ZS + has no bound state (square integrable solution for x ∈ R), then the scattering data are composed of r l (k), r r (k) and t(k) .

For a given value of k, the left and right Jost solutions ν = f l (k, x) and ν = f r (k, x) of ZS + are particular solutions satisfying the limiting conditions

lim x→+∞ f l1 (k, x) e -ikx f l2 (k, x) = 0 1 , and lim x→-∞ e ikx f r1 (k, x) f r2 (k, x) = 1 0 .
The three coefficients constituting the scattering data of ZS + can be expressed as

r l (k) = lim x→-∞ f l1 (k, x) f l2 (k, x) e 2ikx r r (k) = lim x→+∞ f r2 (k, x) f r1 (k, x) e -2ikx t(k) = lim x→+∞ f l2 (k, x)e -ikx lim x→-∞ f l2 (k, x)e -ikx = lim x→-∞ f r1 (k, x)e ikx lim x→+∞ f r1 (k, x)e ikx
To solve the inverse scattering problem as formulated above, it is useful to introduce the following auxiliary system, referred to as ZS -:

dν - 1 (k, x) dx + ikν - 1 (k, x) = q -(x)ν - 2 (k, x) (3a) dν - 2 (k, x) dx -ikν - 2 (k, x) = q + (x)ν - 1 (k, x). (3b) ZS -:     
which is obtained by interchanging the two potential functions q ± (x) of ZS + . The left reflection coefficient, the right reflection coefficient and the transmission coefficient of ZS -will be respectively denoted by r - l (k), r - r (k) and t -(k) in the following.

The solution to the inverse scattering problem of ZS + summarized below will retrieve the two potential functions q ± (x) from the two left reflection coefficients r l (k) and r - l (k) which are respectively related to ZS + and ZS -. In most applications, only the reflection and transmission coefficients of ZS + are available, not those of the auxiliary system ZS -. Fortunately, in the case of real potential functions q ± (x), r - l (k) can be computed from the reflection and transmission coefficients of ZS + based on the following lemmas.

Lemma 1. If the potential functions q ± (x) are real, the following equalities hold, where " * " represents the complex conjugate:

r r (-k) = [r r (k)] * , r l (-k) = [r l (k)] * , t(-k) = [t(k)] * .
This result can be proved by simply taking the complex conjugate at both sides of ZS + . Lemma 2. The reflection and transmission coefficients of ZS + and ZS -are related by

t(k) r r (k) r l (k) t(k) -1 = t -(-k) r - l (-k) r - r (-k) t -(-k) (4)
See [START_REF] Eckhaus | The inverse scattering transformation and the theory of solitons: an introduction[END_REF] for a proof of this result.

By combining the above two lemmas, we have

r - l (k) = r r (-k) r r (-k)r l (-k) -[t(-k)] 2 = r r (k) r r (k)r l (k) -[t(k)] 2 * (5)
Hence r - l (k) can be computed from r l (k), r r (k) and t(k), the scattering data of ZS + , which are thus sufficient for retrieving the potential functions q ± (x) through the inverse scattering transform.

To avoid infinite integral intervals in the GLM equations, we assume that q ± (x) = 0 for all x < 0. This is the case if r - l (k), and r l (k) can be extended analytically to the upper half of the complex k-plane without any pole. Then, as derived in [START_REF] Eckhaus | The inverse scattering transformation and the theory of solitons: an introduction[END_REF], the inverse scattering transform consists of the following steps for computing the potential functions q ± (x) from r l (k), r r (k) and t(k).

1. Compute r - l (k) from r l (k), r r (k) and t(k) through (5). 2. Compute the Fourier transforms,

R l (y) 1 2π +∞ -∞ r l (k)e -iky dk (6a) R - l (y) 1 2π +∞ -∞ r - l (k)e -iky dk (6b)
3. Solve the following GLM integral equations for the kernel functions A r1 (x, y), A r2 (x, y), A - r1 (x, y), and A - r2 (x, y),

A - r1 (x, y) + x -y dsR l (y + s)A r2 (x, s) = 0 (7a) A - r2 (x, y) + R l (y + x) + x -y dsR l (y + s)A r1 (x, s) = 0 (7b) A r1 (x, y) + x -y dsR - l (y + s)A - r2 (x, s) = 0 (8a) A r2 (x, y) + R - l (x + y) + x -y dsR - l (y + s)A - r1 (x, s) = 0 (8b)
4. Compute the potential functions q ± (x) from

q + (x) = 2A - r2 (x, x) (9a) q -(x) = 2A r2 (x, x) (9b)
For more details about the inverse scattering theory, we refer the readers to [START_REF] Eckhaus | The inverse scattering transformation and the theory of solitons: an introduction[END_REF].

Numerical algorithm

Inspired by the results of [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF][START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF][START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF], here we present an efficient numerical algorithm to solve the inverse scattering problem formulated in the previous section.

As the computation of r - l (k) from r l (k), r r (k) and t(k) is trivial, and the Fourier transforms R l (y) and R - l (y) can be computed with the aid of FFT from discretized values of r l (k) and r - l (k), the main numerical computations concern the solution of the GLM equations ( 7) and ( 8) for the kernel functions A r1 (x, y), A r2 (x, y), A - r1 (x, y) and A - r2 (x, y).

In addition to the GLM equations ( 7) and ( 8), the following two pairs of equations respectively derived from ZS + and ZS -will also be used:

∂A r1 (x, y) ∂x + ∂A r1 (x, y) ∂y = q + (x)A r2 (x, y) (10a) ∂A r2 (x, y) ∂x - ∂A r2 (x, y) ∂y = q -(x)A r1 (x, y) (10b) ∂A - r1 (x, y) ∂x + ∂A - r1 (x, y) ∂y = q -(x)A - r2 (x, y) (11a) ∂A - r2 (x, y) ∂x - ∂A - r2 (x, y) ∂y = q + (x)A - r1 (x, y) (11b)
These equations are redundant to the GLM equations ( 7) and [START_REF] Jaulent | The inverse scattering problem for LCRG transmission lines[END_REF]. Their combination with the GLM equations will make easier the solution of the large set of linear algebraic equations resulting from the discretization of the GLM equations.

For the case of single-potential Zakharov-Shabat equations (q + (x) = ±q -(x)), only the pair (10) is used in [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF], where it was proposed as an alternative to the similar second order equations used in [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF]. This improvement made in [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] leads to an algorithm much faster than the original algorithm of [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF].

Similar to the improvement made in [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] to the algorithm of [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF], the use of the two pairs of first order equations ( 10) and (11) in this paper is in replacement of the second order equations used in the algorithm of [START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF].

As in [START_REF] Frangos | A numerical solution to the Zakharov-Shabat inverse scattering problem[END_REF][START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF][START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF], we introduce the following coordinate change which will transform the region of x ≥ |y| in the x-y plane to the first quadrant in the ξ-η plane:

then the potential functions can be expressed in terms of ξ and η as q + (x) = q + (ξ + η), q -(x) = q -(ξ + η)

and the A-kernels will be rewritten as B-kernels after the coordinate change:

Then the GLM equations ( 7), [START_REF] Jaulent | The inverse scattering problem for LCRG transmission lines[END_REF] and equations (9) become

By setting ξ = 0 in the first four equations of (13), the following boundary conditions are derived:

After the coordinate change (12), equations (10) and (11) become

Now discretize all these equations over the grid on the ξ-η plane defined by ξ = id and η = jd, with i, j = 0, 1, 2, . . . , and d being the discretization step size 1 . The discrete counterparts of equations ( 13) and ( 15) then write

The numerical inverse scattering algorithm should solve for B 1 (i, j), B 2 (i, j), B - 1 (i, j), B - 2 (i, j), and hence for q ± (i+j), from the above discretized equations. Usually the discretization grid is defined for i and j ranging from 0 to a large integer value (typically hundreds or thousands). Hence these equations constitute a large system 1 The index notations i, j used here should not be confused with the imaginary unit number.

of linear algebraic equations in the unknowns B 1 (i, j), B 2 (i, j), B - 1 (i, j), B - 2 (i, j). Fortunately, these equations can be efficiently solved by observing the following particular structure of the discretized equations.

In what follows, the L-th diagonal of the discretization grid will refer to the set of all the grid nodes (i, j) such that i + j = L and i, j ≥ 0. For any node (i, j) on diagonal L = i + j of the grid, we said the kernel value B 1 (i, j) is on diagonal L, so are the other similar kernel values. Moreover, q + (i + j) and q -(i + j) are also said to be on diagonal L = i + j.

The kernel value B 1 (i, j) on diagonal L = i + j can be computed from kernel and potential functions values on diagonal L -1, following (16g), since the indexes of each quantity appearing at the right hand side of this equation sum to i + j -1. If the kernel values on diagonal L-1 are already computed, and so is q + (i+j -1), then the kernel value B 1 (i, j) can be easily obtained through (16g). The same observation can be made for the other kernel values and q -(i + j -1) by examining (16h), (16i) and (16j). Therefore, the kernel and potential function values can be computed diagonal by diagonal, by repeatedly increasing the value of L. This reasoning does not apply to nodes at the border of the first quadrant where i = 0 or j = 0, for which the boundary conditions (14) or equations (16b) or (16d) provide the kernel values. Following this algorithm, at each step a single unknown kernel value is computed with a new equation where the other unknowns are already computed in previous steps, without requiring any algebraic operation for eliminating unknonws. This simplicity indicates that the set of linear algebraic equations to be solved is in a triangular form.

The complete algorithm is summarized as follows.

• At diagonal L = 0, apply the boundary conditions B - 1 (0, 0) = 0, B - 2 (0, 0) = -R l (0), B 1 (0, 0) = 0, B 2 (0, 0) = -R - l (0), q + (0) = 2B - 2 (0, 0) = -2R l (0), q -(0) = 2B 2 (0, 0) = -2R - l (0). • For L = 1, 2, 3, . . . , the nodes (i, j) on diagonal L are such that i + j = L.

0). 2. Use (16g) and (16i) to find all the other values of B 1 (i, j) and B - 1 (i, j), and use (16h) and (16j) to find all the other values of B 2 (i, j) and B - 2 (i, j), except B 2 (L, 0) and B - 2 (L, 0). 3. Compute B 2 (L, 0) and B - 2 (L, 0) through (16b) and (16d). 4. Compute q + (L) and q -(L) through (16e) and (16f):

Conclusion

In this paper an efficient numerical algorithm is presented for solving the inverse scattering problem related to the generalized Zakharov-Shabat equations with two potential functions. Inspired by the work of [START_REF] Xiao | An efficient algorithm for solving Zakharov-Shabat inverse scattering problem[END_REF] on Zakharov-Shabat equations with a single potential function, the second order differential equations derived from Zakharov-Shabat equations used in the algorithm of [START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF] are replaced by first order differential equations. This modification has the advantage of triangularizing the large set of linear algebraic equations resulting from the discretization of the GLM equations, hence these equations can be efficiently solved by a one-shot algorithm, instead of the iterative algorithm used in [START_REF] Frangos | Analytical and numerical solution to the two-potential Zakharov-Shabat inverse scattering problem[END_REF].