Mahmoud M R Elsawy

Virginie Nazabal

Mathieu Chauvet

Gilles Renversez

Improved nonlinear plasmonic slot waveguide: a full study

Keywords: numbers: 42.65.Wi, 42.65.Tg, 42.65.Hw, 73.20.Mf

We present a full study of an improved nonlinear plasmonic slot waveguides (NPSWs) in which buffer linear dielectric layers are added between the Kerr type nonlinear dielectric core and the two semi-infinite metal regions. Our approach computes the stationary solutions using the fixed power algorithm, in which for a given structure the wave power is an input parameter and the outputs are the propagation constant and the corresponding field components. For TM polarized waves, the inclusion of these supplementary layers have two consequences. First, they reduced the overall losses. Secondly, they modify the types of solutions that propagate in the NPSWs adding new profiles enlarging the possibilities offered by these nonlinear waveguides. In addition to the symmetric linear plasmonic profile obtained in the simple plasmonic structure with linear core such that its effective index is above the linear core refractive index, we obtained a new field profile which is more localized in the core with an effective index below the core linear refractive index. In the nonlinear case, if the effective index of the symmetric linear mode is above the core linear refractive index, the mode field profiles now exhibit a spatial transition from a plasmonic type profile to a solitonic type one. Our structure also provides longer propagation length due to the decrease of the losses compared to the simple nonlinear slot waveguide and exhibits, for well-chosen refractive index or thickness of the buffer layer, a spatial transition of its main modes that can be controlled by the power. We provide a full phase diagram of the TM wave operating regimes of these improved NPSWs. The stability of the main TM modes is then demonstrated numerically using the FDTD. We also demonstrate the existence of TE waves for both linear and nonlinear cases (for some configurations) in which the maximum intensity is located in the middle of the waveguide. We indicate the bifurcation of the nonlinear asymmetric TE mode from the symmetric nonlinear one through the Hopf bifurcation. This kind of bifurcation is similar to the ones already obtained in TM case for our improved structure, and also for the simple NPSWs. At high power, above the bifurcation threshold, the fundamental symmetric nonlinear TE mode moves gradually to new nonlinear mode in which the soliton peak displays two peaks in the core. The losses of the TE modes decrease with the power for all the cases. This kind of structures could be fabricated and characterized experimentally due to the realistic parameters chosen to model them.

I. INTRODUCTION

In recent years, plasmonics has been receiving a great attention from scientists because of its ability to achieve highly confined electromagnetic field and to manipulate light at nano-scale. Nonlinear plasmonics might be of interest because the nonlinearity can be controlled by the power allowing a management of the solution unlike the linear case. Nonlinear plasmonic slot waveguides (NPSWs) have been a flourishing research topic, at least since 2007, [START_REF] Feigenbaum | Plasmon-soliton[END_REF] due to the strong light confinement in the nonlinear dielectric core ensured by the surrounding metal regions and due to its peculiar nonlinear effects [START_REF] Davoyan | Nonlinear plasmonic slot waveguide[END_REF][START_REF] Davoyan | Selffocusing and spatial plasmon-polariton solitons[END_REF][START_REF] Rukhlenko | Exact dispersion relation for nonlinear plasmonic waveguides[END_REF][START_REF] Ferrando | Variational theory of soliplasmon resonances[END_REF][START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF]. One of the key property of symmetric NPSWs is the existence of asymmetric modes that emerge at critical powers through Hopf bifurcations from the first symmetric mode and also from some higher order symmetric * gilles.renversez@fresnel.fr modes [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. Several applications have already been proposed for NPSWs including phase matching in higher harmonic generation processes [START_REF] Davoyan | Quadratic phase matching in nonlinear plasmonic nanoscale waveguides[END_REF], nonlinear plasmonic couplers [START_REF] Salgueiro | Theoretical study of spatial solitons in planar waveguides[END_REF] or switching [START_REF] Nozhat | Switching power reduction in the ultra-compact Kerr nonlinear plasmonic directional coupler[END_REF], and also for biological and chemical sensors [START_REF] Wijesinghe | Low-loss dielectric-loaded graphene surface plasmon polariton waveguide based biochemical sensor[END_REF][START_REF] Charbonneau | Demonstration of surface sensing using long-range surface plasmon waveguides on silica[END_REF]. Nevertheless, the experimental observation of plasmon-soliton waves in these NPSWs is still lacking even if linear slot waveguides have already been fabricated [START_REF] Han | Experimental realization of subwavelength plasmonic slot waveguides on a silicon platform[END_REF]. Similarly to the case of the single nonlinear dielectric/metal interface structures [START_REF] Ariyasu | Nonlinear surface polaritons guided by metal films[END_REF][12][START_REF] Stegeman | Nonlinear surface plasmons guided by thin metal films[END_REF][START_REF] Kushwaha | Nonlinear Surface Polaritons in Thin Films: TM Case[END_REF][START_REF] Walasik | Low-power plasmon-soliton in realistic nonlinear planar structures[END_REF], the modes already studied in the simple NPSWs suffer from high propagation losses that seriously limit the propagation length of the waves. In the present study, we propose and study an improved structure in which buffer linear dielectric layers are added between the nonlinear dielectric core and the two semi-infinite metal regions. For the TM case, these supplementary layers have two main consequences. First, they modify the types of solutions that propagate in the NPSWs adding new profiles, enlarging the possibilities offered by theses nonlinear waveguides. Second, they reduced overall losses and allow the losses to decrease with the power in most of the cases. For the TE case, they reduce the cut-off core thickness for the fundamental linear mode, and increase its propagation length. For the nonlinear case, the losses decrease with the increase of the power for all the TE modes, and the nonlinear symmetric mode moves gradually to new nonlinear mode in which the soliton intensity profile displays two peaks in the middle of the waveguide. We also demonstrate the bifurcation of the asymmetric nonlinear mode from the symmetric one for the TE polarized waves.

The article is organized as follows. In Sec. II, we describe the model and numerical method we use to study the stationary nonlinear waves in the improved NPSW we propose. In Sec. III, we study the properties of the stationary transverse magnetic (TM) waves. In this case, first, we describe the linear solutions of the new structures in order to be able to classify the solutions that will be obtained in the nonlinear case. Second, we give the properties of TM nonlinear stationary solutions and show that, for some linear parameter configurations, new modal spatial transitions as a function of power occur compared to the simple NPSW case. Third, we provide the full nonlinear phase diagram of the main TM modes as a function of the total power and the buffer layer. In Sec. IV, we prove that the added buffer dielectric layers are able to reduce propagation losses and allow them to decrease with power in most of the cases for the TM case. In Sec. V using the FDTD method we study the stability properties of the main nonlinear TM solutions for the most significant cases. Finally, in Sec. VI, we study the existence and the properties of TE modes for both linear and nonlinear core.

II. MODEL AND METHOD

The improved NPSW we study is depicted in Fig. 1, compared to already studied simple NPSW made of a nonlinear dielectric core surrounded by two semi-infinite metal regions [START_REF] Feigenbaum | Plasmon-soliton[END_REF][START_REF] Davoyan | Nonlinear plasmonic slot waveguide[END_REF][START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF], the improved structure contains additional linear dielectric layers between the core and the metal regions. In this study, we will consider only symmetric structures even if asymmetric simple NPSWs have already been considered [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF]. d core is the core thickness, and d buf is the dielectric buffer thickness of each of the supplementary layers with their permittivities denoted by ε buf . We will only consider either monochromatic transverse magnetic (TM) waves or transverse electric (TE) waves propagating along the z direction in a symmetric one-dimensional NPSW. All field components evolve proportionally to exp[i(k 0 n ef f z -ωt)] with k 0 = ω/c, where c denotes the speed of light in vacuum, ω is the angular frequency, and n ef f denotes the effective mode index. The nonlinear Kerr-type dielectric is isotropic with focusing nonlinearity.

For the TM polarized waves, the permittivity in the nonlinear core is defined as being the linear part of the core permittivity, E x being the x component of the electric field), and with α > 0 (focusing nonlinearity). This approximation of the Kerr nonlinearity, taking into account only the transverse component of the electric field has already been used in several models of the NPSWs [START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF]. It provides similar results than the more accurate approach where all the electric field components are considered in the optical Kerr effect [START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Stationary plasmon-soliton waves in metal-dielectric nonlinear planar structures: modeling and properties[END_REF], even if at high power its results differ from the ones obtained with the full model taking into account all the electric field components in the nonlinear term [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. This approximation allows us to use the fixed power algorithm in the finite element method (FEM) to compute the nonlinear stationary solutions and their nonlinear dispersion curves [START_REF] Walasik | Stationary plasmon-soliton waves in metal-dielectric nonlinear planar structures: modeling and properties[END_REF][START_REF] Rahman | Numerical analysis of nonlinear bistable optical waveguides[END_REF][START_REF] Ferrando | Spatial soliton formation in photonic crystal fibers[END_REF][START_REF] Drouart | Spatial Kerr solitons in optical fibres of finite size cross section: beyond the Townes soliton[END_REF]. In the core, the approximated relation between the nonlinear parameter α and the coefficient n 2 , that appears in the definition of an intensity dependent refractive index ℜe(n core ) = ℜe(ε l,core ) 1/2 +n 2 I core , is α = ε 0 cℜe(ε l,core )n 2 , where the intensity inside the core is defined as

ε core = ε l,core +α|E x | 2 (ε l,core
I core = ε 0 c|ε core ||E x | 2 /(2n ef f
), and the vacuum permittivity is denoted by ε 0 . For the TE case, the only nonzero component of the electric field is E y and the permittivity in the core is defined as ε core = ε l,core + α|E y | 2 , with α > 0 (focusing nonlinearity). In general, the TE surface plasmon waves do not exit due to continuity of the component E y . But, the fundamental TE mode exists only if the the core thickness beyond a cut-off [START_REF] Nakano | Insertion loss and extinction ratio of a surface plasmon-polariton polarizer: theoretical analysis[END_REF][START_REF] Sun | Low-loss, compact waveguiding with TE mode in metal/dielectric waveguides for planar lightwave circuit[END_REF].

To compute the linear and nonlinear modes we use our custom-built FEM we developed to study both scalar and vector models for several types of waveguides including nonlinear ones with Kerr regions like microstructured optical fibers and NPSWs [START_REF] Walasik | Stationary plasmon-soliton waves in metal-dielectric nonlinear planar structures: modeling and properties[END_REF][START_REF] Drouart | Spatial Kerr solitons in optical fibres of finite size cross section: beyond the Townes soliton[END_REF][START_REF] Zolla | Foundations of Photonic Crystal Fibres[END_REF][START_REF] Désévédavy | Te-As-Se glass microstructured optical fiber for the middle infrared[END_REF]. In this study, ε l,core = 3.46 2 + i 10 -4 and n 2 = 10 -17 m 2 /W, corresponding to amorphous hydrogenated silicon [START_REF] Lacava | Nonlinear characterization of hydrogenated amorphous silicon waveguides and analysis of carrier dynamics[END_REF][START_REF] Matres | High nonlinear figure-of-merit amorphous silicon waveguides[END_REF], and ε met = -90 + i 10 corresponding to gold at λ = 1.55 µm [START_REF]Handbook of Optical constants of Solids[END_REF][START_REF] Rakić | Optical properties of metallic films for verticalcavity optoelectronic devices[END_REF], wavelength for which all the simulations are done. The propagation losses will be estimated using the method based on the field profiles and imaginary parts of the permittivities described in [START_REF] Ariyasu | Nonlinear surface polaritons guided by metal films[END_REF][START_REF] Snyder | Optical Waveguide Theory[END_REF], that is similar to the ones from [START_REF] Davoyan | Selffocusing and spatial plasmon-polariton solitons[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF]. We stress that the considered parameter values can affect the critical values given in this study but would not affect its main conclusions.

III. RESULTS FOR TM STATIONARY WAVES

III.1. Linear case

For the TM case, we start by the study of the linear case as it is well known that this is a crucial step in the investigation of nonlinear waveguides [START_REF] Snyder | Linear approach for approximating spatial solitons and nonlinear guided modes[END_REF]. Fig. 2(a) provides the linear dispersion curves fo (solid lines) and with dcore = 300 nm (dashed lines) as a function of the thickness of the supplementary buffer dielectric layers d buf for two different layer refractives indices. Typical magnetic field symmetric components Hy(x) for ε buf = 2.5 2 : (b) A cosine-type profile for d buf = 40 nm such that ℜe(n ef f ) < ℜe(ε l,core ) 1/2 . (c) Flat type profile in the core for d buf = d up buf = 30.82 nm, with ℜe(n ef f ) ≈ ℜe(ε l,core ) 1/2 . (d) Plasmonic type profile for d buf = 15 nm such that ℜe(n ef f ) > ℜe(ε l,core )

1/2 . q ǫ buf - r ǫ met q ǫ buf + r ǫ met = q ǫ buf + p ǫ l,core q ǫ buf - p ǫ l,core e 2k0q(d buf ) , (1)
p 2 = (ℜe(n ef f )) 2 -ǫ l,core , (2a)
q 2 = (ℜe(n ef f )) 2 -ǫ buf , (2b)
r 2 = (ℜe(n ef f )) 2 -ǫ met . (2c)
Equations ( 1) and (2) give the linear dispersion relation for TM waves for a linear dielectric buffer layer sandwiched between a semi-infinite metal region and a semi-infinite high-index dielectric region. By letting (ℜe(n ef f )) 2 → ǫ l,core in Eqns. ( 1) and (2), one gets:

d up buf ≈ λ 4π Re(ǫ l,core ) -ǫ buf × ln Re(ǫ met ) √ Re(ǫ l,core )-ǫ buf -ǫ buf √ Re(ǫ l,core )-Re(ǫ met ) Re(ǫ met ) √ Re(ǫ l,core )-ǫ buf +ǫ buf √ Re(ǫ l,core )-Re(ǫ met ) . ( 3
)
One clearly see that, the critical buffer thickness given by Eq. ( 3), does not depend on the core thickness d core which can be inferred from our FEM simulations depicted in Fig. 2(a). Using the parameters in Sec. II, ε buf = 2.5 2 and d buf = 15 nm, the full formula given by Eq. (3) gives d up buf = 31.81 nm, this result agrees quite well with our FEM simulations (30.8 nm) as it is shown for the blue curves in Fig. 2(a). An approximated value for d up buf (28.4 nm) can be estimated for Re(ǫ l,core ) << |Re(ǫ met )| as:

d up buf ≈ λ 2π ǫ buf |Re(ǫ met )| 1 Re(ǫ l,core -ǫ buf ) . (4)

III.2. Nonlinear case

For the nonlinear case, we study the influence of the supplementary buffer dielectric layers on the nonlinear dispersion curves. We consider the case in which ε buf = 2.5 2 and d core = 400 nm, which corresponds to the solid blue curve in Fig. 2(a) with the associated linear field profiles depicted in Fig. 2(b), (c), and(d). We start with the case such that the effective index of the symmetric linear mode is below the core linear refractive index in which the symmetric linear mode is localized in the high-index core with plasmonic tails in the metal regions (see Fig. 2(b)). The real part of effective indices of the first three main modes: symmetric, asymmetric and antisymmetric denoted by; S0-solI, AS1-solI, and AN0-solI, respectively are shown in Fig. 3(b) with the corresponding field profiles depicted in Fig. 4(d), (e), and (f). The results for a simple NPSW made from the improved NPSW by removing the buffer dielectric layers are also shown in Fig. 3(a) for comparison (symmetric S0-plas, asymmetric AS1-plas, antisymmetric AN0-plas) and the associated field profiles are given in Fig. 4(a), (b), and (c). The mode notation used in our study has been chosen to extend coherently the one already introduced for the simple NPSW [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. As it can be seen, in the improved NPSW, the bifurcation of the asymmetric mode from the symmetric mode and the relative positions of the three main modes are preserved. The difference between the nonlinear symmetric mode curve and the asymmetric nonlinear mode curve is smaller in this new structure (see Fig. 3(b)). The reason is that, in this case, the nonlinear symmetric mode profile looks like a symmetric hyperbolic secant spatial soliton and its associated asymmetric mode is simply shifted off the slot center as it is shown in Fig. 4(d) and(e).

We now study the nonlinear dispersion curves of the improved NPSW such that the effective index of the symmetric linear mode is above the core linear refractive index. The symmetric linear mode is of plasmonic type as it is shown in Fig. 2(d). The corresponding nonlinear dispersion curves are given in Fig. 3(c) and (d) in which the values of the buffer layer thicknesses are d buf = 15 nm and d buf = 10 nm, respectively. The global shape of these dispersion curves is similar to the ones already stud- ied for the three main modes of the simple NPSW (see Fig. 3(a)), but the mode field profiles now exhibit a spatial transition from a plasmonic type profile to a solitonic type one. For the asymmetric one: the mode moves gradually from AS1-plas to AS1-solI, and the antisymmetric mode moves from AN0-plas to AN0-solI. The corresponding field profiles for d buf = 15 nm are shown in Fig. 5(b) and (c). For d buf = 10 nm, the modal transitions for the asymmetric and the antisymmetric modes are shown in Fig. 6 (green and red curves, respectively). Nevertheless, depending on the thickness of the buffer layers, the nonlinear symmetric mode S0-plas exhibits spatial modal transition toward two different nonlinear symmetric solitonic modes. In Fig. 3(c): the symmetric S0-plas mode exhibits spatial transition to symmetric solitonic mode S0-solI; one single peak located in the nonlinear core as it is shown in Fig. 5(a). While, in Fig. 3(d) in which d buf = 10 nm, the symmetric mode S0-plas turns gradually to the new nonlinear symmetric mode denoted by S0-solII (dashed blue line to solid dark blue line in Fig. 3(d)); it is composed of two solitons located in the core near each interfaces (see dark blue curves in Fig. 6).

The solid blue curve in Fig. 3(d) above the dark blue curve, represents the symmetric nonlinear mode S0-solI which emerges at specific value of the power as a purely nonlinear mode. The field profile is similar to the solid blue profiles in Figs. 4(d) and 5(a). In general, the symmetric solitonic mode S0-solI looks like the higher order purely nonlinear symmetric mode denoted by SI which has been found in the simple NPSW, but with different boundary conditions at the core interface due to the opposite sign of the real part of the permittivity in the outer regions [START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. Using numerical simulations providing both the dispersion curves and the associated field profiles like the ones shown above, we obtain the phase diagram of the improved NPSW for the main modes. Such diagram is shown in Fig. 7 where the existence and the type of the first TM symmetric and asymmetric modes are given. The buffer thickness d buf and the total power P tot are used as parameters. The core thickness d core = 400 nm and the permittivity of the buffer layers ε buf = 2.5 2 , are fixed which provide the solid blue curve of the linear main symmetric mode shown in Fig. 2, and the nonlinear dispersion curves given in Fig. 3. For these parameters, the critical value of the buffer layers d up buf in which the fundamental linear symmetric mode is flat in the core (see Fig. 2(c)) is ≈ 30.82 nm. The nonlinear phase diagram depicted in Fig. 7, provides us with the type and the existence of the main symmetric and asymmetric modes that can be found in the improved NPSW as a function of the buffer layer thickness d buf . This phase diagram is divided into three regions. The first region is for d buf > d up buf (right most grey region). In this region, the symmetric linear mode is localized in the high-index core with ℜe(n ef f ) < ℜe(n l,core ) as it is given in Fig. 2

(b).

As the power increases, the main symmetric mode tends to be more localized in the core with plasmonic tails in the metal regions, it is of S0-solI type (symmetric solitontype with a single peak as it is shown in Fig. 4(a)). Due to the nonlinearity, the asymmetric nonlinear mode bifurcates from the nonlinear symmetric mode S0-solI at some critical power value (see solid green curve in Fig. 7). It is denoted by AS1-solI (soliton peak shifted-off the slot center as it is depicted in Fig. 4(b)). Above the bifurcation threshold curve, both the nonlinear symmetric mode S0-solI and the nonlinear asymmetric mode AS1-solI exist as it is shown in the nonlinear dispersion diagram obtained for d buf = 40 nm in Fig. 3(b). The second region is for d down buf < d buf < d up buf (yellow region). In this region, the symmetric linear mode is of plasmonic type (see the field profile in Fig. 2(d)). Below the brown curve labelled "S0-plas to S0-solI" in Fig. 7, the nonlinear symmetric mode is of plasmonic type S0plas. With the increase of power, the nonlinear symmetric mode S0-plas exhibits a spatial modal transition from S0-plas to S0-solI (see Fig. 5(a)), in which the mode becomes more localized in the nonlinear core (single spatial soliton peak in the nonlinear core with plasmonic tails in the metal regions as it is given in Fig. 5(a)). Along the brown curve, the nonlinear symmetric mode is flat in the core (the arrow indicates the transition from S0-plas to S0-solI). Above d buf ≃ 24 nm (but below d up buf ), the solid green curve above the brown thin line in Fig. 7, represents the limit curve for the nonlinear asymmetric mode AS1-solI: below the curve, it does not exist. It emerges directly from the symmetric nonlinear mode S0-solI as expected (see Fig. 4(e) for an example of such field profile). For d down buf < d buf < 24 nm, the dashed green curve is the limit curve of the asymmetric plasmonic mode AS0plas: below the curve, it does not exist. It emerges at low powers from the symmetric one, S0-plas (see the dashed green curve in Fig. 3(c)). At higher powers (see the solid green line in Fig. 7), it moves to an asymmetric mode of solitonic type AS1-solI (a peak being located in the core and near one of its interface, see solid line profile in Fig. 5(b)). Fig. 3(c) represent an example of the dispersion curves, and Fig. 5 shows the corresponding field profiles in this second region for d buf = 15 nm. The nonlinear dispersion curves in Fig. 3(c) can be seen as a vertical cut at d buf = 15 nm in the nonlinear phase diagram depicted in Fig. 7.

The third region: for 0 < d buf < d down buf (left most white region). The symmetric linear mode is of plasmonic type as it can be inferred from Fig. 2(a). Below the cyan curve labelled "S0-plas to S0-solII", the nonlinear symmetric mode is of plasmonic type S0-plas. Above the cyan curve, the new symmetric nonlinear solution S0-solII emerges, at high powers, from the S0-plas mode; it is composed of two solitons located in the core near each interfaces (see the dark blue curves in the first row of Fig. 6). The transition curve between the S0-plas region and the S0-solII region is shown in cyan in the phase diagram shown in Fig. 7. It is worth mentioning that, the symmetric nonlinear mode S0-solII (see first row of Fig. 6) is different from the asymmetric mode shown in Fig. 5(b). Furthermore, this mode exists for core thicknesses as small as 100 nm and it can not be considered in this case as equivalent to the mirror-imaged of the mode existing in the single interface configuration [START_REF] Ferrando | Variational theory of soliplasmon resonances[END_REF]. In this third region, the nonlinear symmetric mode S0-solI (single peak in the core) appears as a purely nonlinear mode without any reference in the linear case. The limit curve for the emergence of the S0-solI is shown in thick pink in the phase diagram. It is worth mentioning that, the appearance of the S0-solI is similar to the appearance of the nonlinear symmetric higher order mode denoted by SI for the simple NPSW [START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF], but with different boundary conditions at the core interface due to the opposite sign of the real part of the permittivity in the outer regions. For d buf = 0 nm (simple NPSW), we recover the same type of modes as expected. It is important to point out that in the region above the cyan curve and below the pink one, only the S0-solII mode exists. The asymmetric nonlinear plasmonic mode AS1-plas emerges at low power from the nonlinear symmetric mode S0-plas (dashed green curve), and moves gradually to the nonlinear symmetric mode AS1-solI (solid green curve). For an example of the nonlinear dispersion curves with the corresponding field profiles in this third region, see Fig. 3(d) and Fig. 6, respectively. This example can be seen as a vertical cut in Fig. 7 at d buf = 10 nm.

As we mentioned before, the linear symmetric mode for d buf = d up buf ≈ 30.82 nm is flat in the core (see Fig. 2(c) for an example of such mode profile). For the nonlinear case, the mode is still nearly flat in the core up to some critical value of the power such that above this value a peak emerges in the core, and the symmetric mode becomes of S0-solI type. The critical power is shown as black circle in Fig. 7 which corresponds to P tot ≈ 10 9 [W/m]. The asymmetric mode AS1-solI bifurcates directly from the symmetric mode S0-solI as expected which is indicated by the crossing between the solid green curve and the the vertical thin line at d buf = d up buf in Fig. 7. The results for the simple NPSW correspond to the vertical line at d buf = 0.0 nm. At low power, the symmetric mode is of plasmonic type S0-plas (see Fig. 4(a)). The bifurcation of the asymmetric nonlinear mode AS1-plas from the symmetric mode S0-plas is represented by the beginning of the dashed green curve at d buf = 0.0 nm. Above this point, both the symmetric S0-plas, and the asymmetric mode AS1plas exist as it is shown in Fig. 3(a). The beginning of the pink curve at d buf = 0.0 nm corresponds to the appearance of the symmetric higher order mode denoted by SI [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. This mode looks like the symmetric mode S0-solI (see solid blue profile Figs. 4 and5), but with different boundary conditions at the core interface as we mentioned before. In Fig. 8, we study the influence of the core thickness d core on the lower critical value d down buf . The horizontal line represents the upper critical value d up buf . The permittivities being fixed, we found that, the upper critical value d up buf does not change with the core thickness as expected (see Eq. ( 3)). While, d down buf increases with d core , such that they coincide when the core thickness reaches the critical value d λ core ≃ 1550 nm (one wavelength), above this core thickness, the critical value d down buf will no longer change with d core . We also provide elsewhere the phase diagram of the improved NPSW using d core and P tot as parameters, for a given fixed values of d buf and ε buf [START_REF] Elsawy | Improved nonlinear slot waveguides using dielectric buffer layers: properties of TM waves[END_REF]. We found that the threshold for the Hopf bifurcation, from which the asymmetric mode emerges, strongly decreases with the core thickness as it was already found in the simple NPSW [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. We provided the nonlinear phase diagram describing the existence and the type of the main TM modes in our improved structure. We have seen the importance of the inclusion of the supplementary buffer layers on the types of solutions that propagate in the NPSWs. Nevertheless, these supplementary layers slightly shift the bifurcation threshold toward higher powers as it is depicted in Fig. 7. While, the increase of the core thickness, reduces the bifurcation threshold [START_REF] Elsawy | Improved nonlinear slot waveguides using dielectric buffer layers: properties of TM waves[END_REF] as expected from the simple NPSWs. [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF].

IV. STUDY OF THE LOSSES FOR TM WAVES

In the previous section III.2, we have shown the influence of the supplementary buffer layers on the type of modes (for the TM case) that propagate in the simple NPSW. We demonstrated that, even for a very thin buffer dielectric layers, the type of the nonlinear modes propagating in the simple NPSW is modified (see phase diagram in Fig. 7 and the nonlinear dispersion curves depicted in Fig. 3). In this section, we study the influence of the buffer layers on the overall losses. The imaginary part of the effective indices of the main modes with the corresponding values of the losses are shown in Fig. 9 which corresponds to Fig. 3 for the real part of the effective indices. The losses are estimated using the approach based on the imaginary part of the permittivity and the field profiles. This method is described in the case of linear waveguides [START_REF] Snyder | Optical Waveguide Theory[END_REF] and has also been used for the nonlinear studies [START_REF] Walasik | Low-power plasmon-soliton in realistic nonlinear planar structures[END_REF][START_REF] Walasik | Stationary plasmon-soliton waves in metal-dielectric nonlinear planar structures: modeling and properties[END_REF].

For the simple slot configuration, the losses of the three main modes increase with the power as it is shown in Fig. 9(a). For d buf = 40 nm, in which the symmetric linear mode is core localized and the nonlinear modes are of solitonic type (see Fig. 4), we found that the losses are at least three times smaller than the ones of the simple NPSW, and that they decrease with the power unlike the simple NPSW (see Fig. 9(b)). We now study the losses for the improved NPSW in which the effective index of the symmetric linear mode is above the core linear refractive index. The corresponding curves are given in Fig. 9(c) and(d). As we have seen before in Fig. 3(c (d), the mode field profiles now exhibit a spatial transition from a plasmonic type profile to a solitonic type one. Furthermore, at low powers, the losses are reduced, at least, by a factor two compared to the associated simple NPSW. For d buf = 15 nm, the losses decrease with the power for the main symmetric and antisymmetric modes leading to low loss solutions at high powers (see blue and red curves in Fig. 9(c) and the associated field profiles in Fig. 5(a) and (c)). For the asymmetric mode, loss first increases with power and then decreases. This loss increase is linked to the fact that, in this case, the field first tends to increase near one of the metal interface, increasing the fraction of field in a highly lossy region. Then, the mode moves gradually from a plasmonic type to a solitonic type asymmetric profile, and when the peak of the soliton is far enough from the metal interface the losses start to decrease as it is shown in Fig. 5(b). For d buf = 10 nm, losses of the three main modes first, increase with the power and then decrease. The increase of the losses can be understood by looking at the field profiles in Fig. 6, in which the fields tend to increase near the metal regions. With the increase of the power, the modes move gradually from plasmonic type to solitonic type: S0-plas to S0-solII (for the fundamental symmetric mode), AS1-plas to AS1-solI (for the asymmetric mode), and S0-plas to AN0-solI (for the antisymmetric mode), in which the soliton peaks emerge near the interfaces, and when the soliton peaks are far enough from the metal regions, the losses start to decrease. The symmetric higher order mode S0-solI (see solid blue line in Fig. 9(d)), has lower loss and the loss decrease with the increase of the power as expected due to the localization of the soliton peak in the nonlinear core (see solid blue curve in Fig. 5(a) for an example of such field profile). It is worth mentioning, that the increase of the losses in Fig. 9(d), is similar to the simple NPSW in Fig. 9(a). While, due to the inclusion of very thin buffer layers d buf = 10 nm, the behaviour of the losses change at high power (where the field profiles ex-hibit spatial transitions from plasmonic to solitonic type solutions), which proves the importance of the inclusion of the dielectric buffer layers between the nonlinear core and the metal regions in the simple NPSW.

V. STABILITY STUDY FOR TM WAVES

In this section we study the stability properties for the TM waves in the most significant cases. In a previous work, using two different numerical methods including the FDTD, we have already shown that the results about the stability of the main symmetric and asymmetric modes derived in reference [START_REF] Mitchell | Stability of fundamental nonlinear guided waves[END_REF] from theoretical arguments in the framework of the weak guidance approximation in nonlinear waveguides are still valid, at least numerically, in simple NPSWs except just around the bifurcation point where subtle behaviours occur [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. Here, we have again used the nonlinear capabilities of the FDTD method implemented in the meep software [START_REF] Oskooi | Meep: A flexible freesoftware package for electromagnetic simulations by the FDTD method[END_REF][START_REF] Rodriguez | χ 2 and χ 3 harmonic generation at a critical power in inhomogeneous doubly resonant cavities[END_REF]. The metal permittivity is described by a Drude model to obtain the fixed negative value used at the studied wavelength in the stationnary study. As it can be seen in FIG. 10. Evolution during the propagation of the Hy field profile for dcore = 300nm, d buf = 40 nm. (a) S0-solI symmetric solitonic type mode at a power such that ∆n = 0.043 with ε buf = 2.5 2 . (b) S0-plas symmetric plasmonic type mode at a power such that ∆n = 0.015 with ε buf = 3.2 2 . ∆n is the spatial average along the core of the nonlinear refractive index change [START_REF] Walasik | Symmetric plasmonic slot waveguides with a nonlinear dielectric core: Bifurcations, size effects, and higher order modes[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF], other parameters as in the text.

Fig. 10(a), for the improved NPSW, the symmetric mode S0-solI is stable at low power. The S0-plas mode (plasmonic type) is also stable as it is depicted in Fig. 10(b). The asymmetric mode AS1-solI is also stable for the tested powers above the bifurcation threshold (data not shown). Due to the losses of the metal, and as pointed out in [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF][START_REF] Salgueiro | Complex modes in plasmonic nonlinear slot waveguides[END_REF], the temporal nonlinear solutions are actually not self-coherent [START_REF] Snyder | Linear approach for approximating spatial solitons and nonlinear guided modes[END_REF] as it can be seen in Fig. 10 with the decrease of the field intensity along the propagation.

VI. TE STATIONARY WAVES

In this section we study the transverse electric TE polarized waves in the improved structure depicted in Fig. 1.

In this case the electromagnetic fields have only three nonzero components E y , H x , and H z . The only nonzero component of the electric field is the longitudinal component E y . Consequently, the TE surface plasmon waves do not exist in isotropic planar linear metal/dielectric structures. However, the fundamental transverse electric mode for linear (metal/dielectric/metal) structures exists if the core width is beyond a cut-off thickness [START_REF] Nakano | Insertion loss and extinction ratio of a surface plasmon-polariton polarizer: theoretical analysis[END_REF][START_REF] Sun | Low-loss, compact waveguiding with TE mode in metal/dielectric waveguides for planar lightwave circuit[END_REF]. For this mode, the maximum intensity is located in the core of the waveguide, thus losses in the TE case are much smaller than their counterparts in the TM case in which the field penetrates into the metal.

VI.1. Linear TE case

As we saw in Sec. III.1 for the TM case, the linear case is a fundamental step in understanding the nonlinear behaviour. For this reason, we start by the linear case in our study of the TE modes. In Fig. 11, we compare the fundamental symmetric linear profiles for both TM and TE waves for two different values of d buf . We choose the core thickness d core = 400 nm, which is above the cut-off core thickness of the symmetric linear TE mode (there is no cut-off thickness for the linear TM mode). The linear TE mode being core localized for small values of d buf , even for large d buf it stays core localized (see Fig. 11(b) and (d)). Consequently, unlike the TM mode, it does not exhibit any spatial transition as a function of d buf . While, for the TM case (see Fig. 11(a) and (c)) the buffer layer thicknesses d buf , play an important role in distinguishing between two different types of modes (plasmonic and core localized modes) as we mentioned before in Sec. III.1. The real and the imaginary parts of the effective index for the fundamental symmetric linear TE mode S0-TE (see Fig. 11(b) and (d) for an example of such field profile) are given in Fig. 12 and 13, respectively, as a function of the core thickness d core for five different values of d buf . For each value of d buf , there exist a cut-off thickness for d core such that below this value the fundamental symmetric linear TE mode does not exit. The cut-off core thickness decreases with the increase of the buffer layer thicknesses d buf such that the black curve (which represents the results obtained of the improved structure by removing the supplementary buffer layers i.e the simple NPSW) has the highest cut-off core thickness ≈ 172 nm, while the dark green curve has the lowest cut-off value ≈ 40 nm for d buf = 100 nm. It is worth mentioning that, for d core > 600 nm all the curves merge together regardless the value of the buffer thicknesses, because the core becomes very thick and the S0-TE mode moves away from the metal cladding. In Fig. 13, we show the

(1) S0-TE, d buf =100nm

(2) S0-TE, d buf =40nm

( imaginary parts of the fundamental symmetric linear TE mode S0-TE with the corresponding propagation length L p . The meaning of this figure is similar to Fig. 12. The solid curves for Im(n ef f ) and the dashed ones for the propagation length L p = λ/ (4πIm(n ef f )). For a fixed value of d core , the improved structure supports a fundamental symmetric linear TE mode with longer propagation length than the one obtained from the simple structure. The dark green curve (for d buf = 100 nm), has the longest propagation length and the black curve has the shortest propagation length (for d buf = 0.0 nm).

VI.2. Nonlinear TE case

In this part, we study the nonlinear dispersion curves for the improved NPSW depicted in Fig. 1, for TE polarized waves. In Subsec. VI.1, we have shown that, as soon as the cut-off core thickness is exceeded, the supplementary buffer layers for the TE case have no significant effects on the linear field profiles, unlike the TM case, as it is shown in Fig. 11(b) and(d). Nevertheless, the inclusion of the linear buffer layers between the nonlinear core and the metal regions allows us to reduce the cut-off thickness of d core for the fundamental symmetric linear TE mode (see Fig. 12), and to increase its propagation length (see Fig. 13). For the nonlinear TE case, we consider d core = 400 nm, ε buf = 2.5 2 , and d buf = 15 nm (for the other parameters see Sec. II), these parameters are used for linear field profile shown in Fig. 11(b). We chose the core thickness d core = 400 nm to be above the cut-off thickness which is ≈ 147 nm for these parameters (see green curve in Figs. 12 and13). The only non-zero component of the electric field is the continuous component E y therefore, the nonlinear permittivity of Kerrtype in the core is written as ε core = ε l,core +α|E y | 2 , with α > 0 (focusing nonlinearity). We use the fixed power algorithm in the FEM [START_REF] Walasik | Stationary plasmon-soliton waves in metal-dielectric nonlinear planar structures: modeling and properties[END_REF][START_REF] Drouart | Spatial Kerr solitons in optical fibres of finite size cross section: beyond the Townes soliton[END_REF][START_REF] Zolla | Foundations of Photonic Crystal Fibres[END_REF] to compute the nonlinear dispersion curves and the associated field profiles, in which the input is the total power and the outputs are the effective index and the corresponding field component (E y (x) in this case).

In Fig. 14, we present the nonlinear dispersion diagram for the main nonlinear TE modes in the improved NPSW. The global shape of this diagram is similar to the ones obtained for the TM case such that the relative positions of the three main modes are preserved. However, the field profiles are different. The blue, green, and red curves represent the symmetric, asymmetric, and antisymmetric nonlinear TE modes, respectively. The field profiles for the nonlinear symmetric mode S0-TE are given in Fig. 15(a), (b), and (c) in which the intensity is located in the middle of the waveguide without any plasmonic tails at the interfaces between the linear buffer layers and the metal regions. The antisymmetric field profiles are depicted in Fig. 16(d), (e), and(f).

Due to the nonlinearity, the nonlinear asymmetric mode (denoted by AS1-TE), bifurcates at a critical value of the power through a Hopf bifurcation from the nonlinear symmetric mode S0-TE (see the green and the blue curves in Fig. 14). The nonlinear asymmetric mode AS1-TE profiles are shown in Fig. 16(a), (b), and (c) in which the mode is shifted-off the slot center and moves to one of the interfaces (we remind that the asymmetric mode is doubly degenerated). This kind of Hopf bifurcation is similar to the ones obtained for the TM wave both for the improved NPSW (see Fig. 3) and for the simple NPSW [START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. I. modeling[END_REF][START_REF] Walasik | Plasmon-soliton waves in planar slot waveguides. II. results for stationary waves and stability analysis[END_REF]. If we compare the bifurcation thresholds for the TE case (see green curve in Fig. 14) with its counterpart in the TM case (shown in Fig. 3(c)) using the same parameters, we found that the asymmetric nonlinear mode in the TE case bifurcates from the symmetric nonlinear mode at a higher power than its counterpart in the TM case. Furthermore, in the TM case, the asymmetric mode bifurcates from the symmetric plasmonic mode S0-plas, and turns gradually to the asymmetric solitonic mode AS1-solI as it is shown in Fig. 5(b), while in the TE case, the asymmetric mode AS1-TE, bifurcates from the nonlinear symmetric mode S0-TE and does not exhibit any spatial transition (see Fig. 14). At high power, above the bifurcation threshold in the nonlinear dispersion diagram shown in Fig. 14, the nonlinear symmetric core localized mode S0-TE turns gradually to a nonlinear two-humped mode. We call it SII-TE (see the blue and thick dark blue curves in Fig. 14), in which the soliton intensity profile in the core displays two peaks in the middle of the waveguide. In Fig. 15, we illustrate the modal spatial transition from the nonlinear symmetric mode S0-TE (blue profiles), to the nonlinear symmetric two-humped mode SII-TE (dark blue profiles). It is worth mentioning that such kind of transition between the symmetric modes has already been found in a fully nonlinear dielectric waveguides [START_REF] Trillo | Spatial Solitons[END_REF] at high power.

In Fig. 17, we show the imaginary part of the effective indices of the main TE modes with the corresponding values of the losses. One clearly see that, at low power, the antisymmetric mode has the highest value of the losses (see Fig. 17 to the metal regions as it is shown in the second row of Fig. 16. Whereas, at high power the losses of the three main modes become of the same order (see Fig. In general, the losses decrease with the power increase for all the tested modes. This decrease of the losses with the increase of the power, can be understood from the field profiles in Figs. [START_REF] Walasik | Low-power plasmon-soliton in realistic nonlinear planar structures[END_REF] and 16 in which at high power, all the modes are more localized in the core and move away from the metal regions.

VII. CONCLUSION

We have presented a full study of a nonlinear plasmonic slot waveguide improved by the addition of supplementary linear dielectric buffer layers between the nonlinear core of Kerr-type and the semi-infinite metal regions. We have provided a complete description of its main TM and TE modes for both linear and nonlinear case. For the linear TM case, in addition to the symmetric linear plasmonic profile obtained in the simple structure (without the buffer layers), our improved structure supports two new types of linear modes; the core localized and the flat modes. These new linear modes are not found in the simple symmetric linear structures whatever the opto-geometric parameters used. For the nonlinear TM case, we have provided the nonlinear phase diagram for the existence and the type of the main TM modes. The buffer layer thickness and the total power are used as parameters such that, for a well-chosen opto-geometric parameters, the main TM modes exhibit spatial modal transitions which can be controlled by the power. We have also shown that this structure allows a decrease of the losses compared to the simple nonlinear slot waveguide and that they decrease with the power for most of the cases. The stability properties of the main TM modes has also been studied numerically using the FDTD method. We have also demonstrated the existence (for some configurations), and studied the properties of the main TE modes for both linear and nonlinear cases. This kind of structure could be fabricated and characterized experimentally due to the realistic parameters chosen to model it.

FIG. 1 .

 1 FIG.1. Symmetric improved NPSW geometry with the supplementary linear dielectric buffer layers between the nonlinear dielectric core and the two semi-infinite metal regions.

FIG. 2 .

 2 FIG. 2. (a)Linear dispersion curves of the fundamental symmetric TM mode for the improved NPSW with dcore = 400 (solid lines) and with dcore = 300 nm (dashed lines) as a function of the thickness of the supplementary buffer dielectric layers d buf for two different layer refractives indices. Typical magnetic field symmetric components Hy(x) for ε buf = 2.5 2 : (b) A cosine-type profile for d buf = 40 nm such that ℜe(n ef f ) < ℜe(ε l,core ) 1/2 . (c) Flat type profile in the core for d buf = d up buf = 30.82 nm, with ℜe(n ef f ) ≈ ℜe(ε l,core ) 1/2 . (d) Plasmonic type profile for d buf = 15 nm such that ℜe(n ef f ) > ℜe(ε l,core ) 1/2 .

FIG. 3 .FIG. 4 .

 34 FIG. 3. Nonlinear dispersion curves for the simple NPSW (grey curves) and for the improved one (color curves). The blue, green, and red colors indicate the symmetric, asymmetric and antisymmetric nonlinear TM modes, respectively. dcore = 400 nm, ε buf = 2.5 2 . (b) for d buf = 40 nm, (c) for d buf = 15 nm, and (d) for d buf = 10 nm.

FIG. 5 .

 5 FIG.5. Field profiles Hy(x) for the symmetric (blue curves), asymmetric (green curves), and antisymmetric (red curves) TM modes with d buf = 15 nm, before (dashed line) and after (solid line) the spatial modal transition induced by the power increase. The power is given in W/m.

  FIG. 7. Phase diagram of the TM modes of the improved NPSW as a function of dielectric buffer thickness d buf and of the total power Ptot for dcore = 400 nm and ε buf = 2.5 2 . The green curves refer to the first asymmetric modes. Arrows indicates the corresponding mode transition across the associated curve (thin one).

FIG. 8 .

 8 FIG. 8. Evolution of the critical value d down buf as a function of the core thickness dcore for ε buf = 2.5 2 and d buf = 15 nm. The horizontal line represents the critical value d up buf .

FIG. 9 .

 9 FIG. 9. Losses for the simple NPSW (grey curves) and for the improved one (color curves). The blue, green, and red colors indicate the symmetric, asymmetric and antisymmetric nonlinear TM modes, respectively. dcore = 400 nm and ε buf = 2.5 2 . (b) for d buf = 40 nm, (c) for d buf = 15 nm, and (d) for d buf = 10 nm. For the real part of the effective indices see Fig. 3. For the field profiles see Figs. 4, 5, and 6.

FIG. 11 .

 11 FIG. 11. Distribution of the symmetric linear profiles for fundamental symmetric linear modes with dcore = 400 nm and ε buf = 2.5 2 for two different buffer layers thicknesses d buf = 15 nm and d buf = 40 nm. (a) and (c) for fundamental linear TM modes without cut-off thickness. (b) and (d) for linear fundamental TE mode above the cut-off thickness of dcore.

FIG. 12 .

 12 FIG.12. Linear dispersion relation for the symmetric fundamental transverse electric mode S0-TE as a function of the core thickness dcore for ε buf = 2.5 2 . The color curves represent the linear dispersion relation for different thicknesses of the buffer layer d buf . The cut-off thickness for dcore is given in the top region of the figure for each plotted curve.

FIG. 13 .

 13 FIG. 13. Im(n ef f ) and the propagation distances for the linear symmetric fundamental transverse electric mode S0-TE as a function of the core thickness dcore for ε buf = 2.5 2 . The color curves represent the results for different values of d buf . The solid lines are for Im(n ef f ), and the dashed lines are for the propagation length Lp (see Fig. 12 last sentence).

FIG. 14 .

 14 FIG.[START_REF] Kushwaha | Nonlinear Surface Polaritons in Thin Films: TM Case[END_REF]. Nonlinear dispersion curves of the improved NPSW for the TE polarized waves. The blue, green, and red colors indicate the symmetric, asymmetric and antisymmetric nonlinear TE modes, respectively. dcore = 400 nm, ε buf = 2.5 2 , and d buf = 15 nm.

  (a)) because the peaks are closer

FIG. 15 .

 15 FIG. 15. Spatial transition of the symmetric nonlinear field profiles Ey(x) for the improved NPSW in the TE case. (top) for S0-TE and (bottom) for SII-TE with dcore = 400 nm, ε buf = 2.5 2 , and d buf = 15 nm. The power is given in W/m.

FIG. 16 .

 16 FIG.[START_REF] Lacava | Nonlinear characterization of hydrogenated amorphous silicon waveguides and analysis of carrier dynamics[END_REF]. Nonlinear field profiles Ey(x) for the improved NPSW in the TE case. The top green curves and the bottom red curves for AS1-TE and AN0-TE, respectively, with dcore = 400 nm, ε buf = 2.5 2 , and d buf = 15 nm. The power is given in W/m.

17

 17

  (b)).

FIG. 17 .

 17 FIG.[START_REF] Matres | High nonlinear figure-of-merit amorphous silicon waveguides[END_REF]. Im(n ef f ) and the losses for the TE waves in the improved NPSW. The blue, green, and red colors indicate the symmetric, asymmetric and antisymmetric nonlinear TE modes, respectively. dcore = 400 nm, ε buf = 2.5 2 , and d buf = 15 nm.

ACKNOWLEDGMENTS

This work has been partially funded by the PhD school of Physics and Materials Science, ED352, and the University of Aix-Marseille. It was not funded by the French National Research Agency.


Olivier Aumage

Denis Barthou

Alexandre Honorat

A Stencil DSEL for Single Code Accelerated Computing with SYCL

Keywords: Stencil, DSEL, meta-programming, GPU, parallelization, HPC

Stencil kernels arise in many scientific codes as the result from discretizing natural, continuous phenomenons. Many research works have designed stencil frameworks to help programmer optimize stencil kernels for performance, and to target CPUs or accelerators. However, existing stencil kernels, either library-based or languagebased necessitate to write distinct source codes for accelerated kernels and for the core application, or to resort to specific keywords, pragmas or language extensions.

SYCL is a C++ based approach designed by the Khronos Group to program the core application as well as the application kernels with a single unified, C++ compliant source code. A SYCL application can then be linked with a CPU-only runtime library or processed by a SYCL-enabled compiler to automatically build an OpenCL accelerated application. Our contribution is a stencil domain specific embedded language (DSEL) which leverage SYCL together with expression template techniques to implement statically optimized stencil applications able to run on platforms equipped with OpenCL devices, while preserving the single source benefits from SYCL.

Introduction

Stencil kernels involve concurrently updating every array element with a function of itself and its close neighbourhood cells, using a regular pattern whose width, layout and dimensions characterize a given stencil. Such kernels arise in many scientific codes resulting from the discretization of natural phenomenons with local, short range interactions. While simple in principle, stencil kernels are notoriously hard to code, and even more to optimize, as the arithmetics on array indices and boundary managements are tricky and error prone. Thus many efforts have been put within the high performance computing community on designing stencil programming frameworks and languages, to hide part of this complexity.

Yet, HPC platforms now evolve towards heterogeneity, using nodes equipped with accelerators, resulting in increased complexity and reduced portability. Stencil programming frameworks have This is the author's version of the work, published in the following workshop: SYCL'16, March 13, 2016, Barcelona, Spain. to adapt to heterogeneity but face the dual language issue. A heterogeneous stencil framework has to generate general purpose C/C++/Fortran code for the main application side and CPU kernel instances, and some accelerator code such as OpenCL for the accelerated stencil kernel instances. The SYCL [START_REF] Group | C++ single-source heterogeneous programming for OpenCL[END_REF] language has been designed by the Khronos Group consortiumto address the dual language issue in heterogeneous computing. It leverages regular C++ constructs to enable a SYCL compiler to extract accelerator kernels automatically from a unique valid C++ application source code, without the need for pragmas or keywords.

The contribution of this paper is a C++ Domain Specific Embedded Language (DSEL) based on SYCL to handle stencil computations. It builds on SYCL capabilities and expression template techniques to dramatically reduce the complexity of stencil kernel programming for heterogeneous platforms down to a few lines of code. Our DSEL has been tested with both the triSYCL [START_REF] Keryell | triSYCL[END_REF] CPU-only SYCL library (since SYCL is pure C++, CPU-only versions can indeed be implemented as mere C++ libraries!), and with a preliminary version of Codeplay's ComputeCPP SYCL implementation for OpenCL platforms. A preliminary version of the DSEL is available on the triSYCL repository [START_REF] Keryell | triSYCL[END_REF] inside the tests/jacobi folder, this version only works for triSYCL.

This paper is organized as follows. Section 2 depicts the context and discusses related works on stencil computing frameworks. Section 3 introduces our heterogeneous stencil DSEL on top of SYCL. Section 4 details the implementation and parallelization techniques we used. Section 5 concludes this paper and discusses future work.

Context: Computing Stencils on Heterogeneous Platforms

Stencils are widely used in scientific computing especially for image processing, and for physics simulation. Beyond scientific computing, stencils are also used in other performance demanding domains such as special effects in video games [START_REF] Stam | Real-time fluid dynamics for games[END_REF]. Thus, stencil kernels are key optimization targets for many performance dependent applications. However, the art of stencil writing requires skills in both computer science and in the considered application domain, which may be hard to find in the same programmer, or even in the same scientific team.

To reduce or even remove this impeding dual skill prerequisite, many research efforts have been put into designing stencil computing frameworks that separate the generic "computer science side" stencil execution and optimization work, from the "application specific" mathematical stencil expression, that is to implement the so called separation of concerns. Two fundamental approaches have been followed for doing so: compilers and libraries. A specific compiler can generate highly optimized code, but it forces the user to write the stencil code in a Domain Specific Language (DSL) which may not be practical for large applications. At the opposite, a library allows to write code with the same usual language as the rest of the application, still having several implementations dedicated to the desired targets. We refer to it as a Domain Specific Embedded Language (DSEL) since it provides a set of dedicated methods within and in compliance with the general purpose language used by the application, and as such can be compiled with any compiler available for this general language. C++ is frequently used as a general purpose language for embedding DSELs due to its ability at being specialized through templating.

Dedicated stencil libraries and compilers include the following ones, among many others. Pochoir [START_REF] Tang | Coding stencil computations using the pochoir stencil-specification language[END_REF] uses a stencil declaration language embedded in C++ and provides its own stencil kernel compiler to generate optimized versions, and to provide a level of error checking. Blitz++ [START_REF] Veldhuizen | Arrays in blitz++[END_REF] is a standalone library. By design, its optimizing capabilities are therefore less aggressive than Pochoir's compiler, but it instead offers some predefined stencil patterns (such as the laplacian operator) and a wide panel of functions to process arrays. PATUS [START_REF] Christen | PATUS: A code generation and autotuning framework for parallel iterative stencil computations on modern microarchitectures[END_REF] is a framework for stencil kernel generation and autotuning, which targets both CPUs and accelerators such as nVidia CUDA GPUs, in contrast to Pochoir and Blitz++. It implements a DSL for letting the programmer express the stencil expression and layout. Some general purpose libraries can instead be used to implement stencils such as Kokkos [START_REF] Edwards | Kokkos tutorial, a Trilinos package for manycore performance portability[END_REF]. Kokkos helps with processing arrays in parallel but does not optimize user algorithms. The mapreduce paradigm provides a mean to apply a stencil mapped as possible in SkePU [START_REF] Dastgeer | Skeleton programming for heterogeneous gpu-based systems[END_REF], Bones [START_REF] Nugteren | Bones: An automatic skeleton-based c-to-cuda compiler for gpus[END_REF][START_REF] Nugteren | Algorithmic species: A classification of affine loop nests for parallel programming[END_REF] or Java Lime [START_REF] Dubach | Compiling a high-level language for gpus: (via language support for architectures and compilers)[END_REF], but they cannot be as efficient as dedicated stencil libraries.

Applying a stencil on an array and storing the results on another array is an operation with significant potential for SIMD/vector processing and more generally for accelerator and GPU processing. Exploiting the benefit of such hardware for stencil processing is challenging, especially with respect to portability. Some parameters such as workgroup size and local shared memory size depend on the hardware characteristics. An image filtering stencil C++ embedded language has been developed on the HIPA-CC [START_REF] Membarth | Towards domainspecific computing for stencil codes in HPC[END_REF] (Heterogeneous Image Processing Acceleration) sourceto-source framework to generate optimized OpenCL/CUDA kernels. The Liszt [START_REF] Devito | Liszt: a domain specific language for building portable mesh-based pde solvers[END_REF] stencil compiler and tool-chain also support CPU/GPU code generation as well as distributed MPI-based code generation. The OP2 project [START_REF] Mudalige | OP2: An active library framework for solving unstructured mesh-based applications on multi-core and many-core architectures[END_REF] is dedicated to unstructured mesh computations and can produce OpenMP/AVX vectorized code for CPUs, and CUDA or OpenCL for GPUs. The compiler presented by Holewinski et al. [START_REF] Holewinski | High-performance code generation for stencil computations on gpu architectures[END_REF] generates efficient stencil kernel, but it only targets GPUs. Yet, CPUs may indeed also be used aside or as fall-back for accelerators. Some stencil processing techniques can be leveraged for both CPU and accelerators, such as hierarchical overlapped tiling [START_REF] Zhou | Hierarchical overlapped tiling[END_REF] or 3.5D blocking [START_REF] Nguyen | 3.5dd blocking optimization for stencil computations on modern cpus and gpus[END_REF].

Finally, some frameworks can be specialized to some specific stencil application fields. For instance, quantum chromodynamics (QCD) simulations involve stencils in 4 dimensions, with variable coefficients depending on the direction. Heybrock et al. [START_REF] Heybrock | Lattice QCD with domain decomposition on intel R xeon phi co-processors[END_REF] present a method to efficiently compute QCD kernels on many-cores processors, but the computer science skills required from the programmer to adapt it to other stencil application fields limit its usability by a wider range of non specialists. The Qiral DSL is dedicated to generating an efficient QCD code [START_REF] Barthou | Automated code generation for lattice quantum chromodynamics and beyond[END_REF], unfortunately for CPU only.

In contrast to the related works above, our contribution is instead to build a generic stencil processing framework on top of the recently ratified SYCL standard [START_REF] Group | C++ single-source heterogeneous programming for OpenCL[END_REF] introduced in the next Section, which leverages the unique properties of SYCL to reconcile the library approach and the compiled approach, to target both CPU and accelerator, to require a single source code, written in the C++ general purpose language.

Contribution: a SYCL-Based Stencil DSEL for Heterogeneous Platforms

Our contribution is to develop a stencil DSEL that can produce optimized parallel code for all target types that can be managed by SYCL. SYCL provides a common hardware interface and handles memory management. Our library specializes the code with respect to the stencil parameters (dimensions, order, pattern). This approach enables to achieve stencil kernel parallelization without the need for developing a compiler. More importantly, our approach frees the programmer from writing complicated parallel code. Indeed, he/she only needs to describe the application dependent stencil parameters and the data to be processed.

SYCL in a Nutshell

SYCL [START_REF] Group | C++ single-source heterogeneous programming for OpenCL[END_REF] is a C++ interface standardized by the Khronos Group, currently at revision 1.2. It aims to simplify parallel code writing in modern computers by overlaying the OpenCL parallelism standard. By bridging the gap between C++ and OpenCL, it enables a single application source code to target regular CPU and OpenCL enabled devices such as accelerators and GPUs. Moreover, OpenCL in itself is not easy to use for an end-user programmer, since he/she must handle many initialization and setup steps such as device discovery, kernel loading and compilation, memory management and transfers between the host and the accelerator board, and so on. A typical OpenCL application indeed starts with dozens of such preparatory calls. SYCL hides most of these steps away, and let the user focus on writing application code instead. Yet, the programmer still retains access to most OpenCL methods and can therefore reuse existing OpenCL code. Two parallelization models are supported in SYCL. The global model is a flat "parallel for" loop. The hierarchical model defines workgroups and workitems concepts as in OpenCL, to abstract teams of cooperating hardware threads, and individual threads in such teams respectively. The use of modern C++ functionalities enable to write kernels with lambda functions, and also, for example, to have implicit barriers between two workitem method calls. These capabilities makes SYCL easy to integrate in other libraries. The common OpenCL memory spaces are addressable in SYCL: host, global, local, private and image. Data must be stored in SYCL buffers and can then be referenced in kernels through SYCL accessors. Accessors are used to associate a piece of data with one of the aforementioned memory space.

Currently two SYCL implementations exist. triSYCL [START_REF] Keryell | triSYCL[END_REF] is a standard C++ library designed by R. Keryell which targets CPU exclusively. It is described as a "test-bed to experiment with the provisional specification" of SYCL. ComputeCPP from company Codeplay is a commercial SYCL implementation which targets both CPUs and OpenCL devices. In this work, we experimented with running our own stencil library using both of these SYCL implementations. As a cautionary note, since -at the time of this writing -triSYCL has an experimental status, and since the Com-puteCPP version we used is a preliminary, early access version graciously provided by Codeplay, we will not expose quantitative performance results in this paper.

Stencil DSEL Description

A stencil is an elements' grid pattern which can be applied on every element in order to update its value. The stencil is mapped over the grid and consists of a reduction on the element neighbourhood values. An example of a generic 4-points stencil in two dimensions (involving the two nearest neighbours in each dimension and excluding the central point) is shown in equation 1. A, B and C arrays store the source elements, the destination elements and the coefficients, respectively. Whereas the neighbour elements are by defini-tion close to the target element in the logical grid, they may not be close in memory, thus the DSEL first requires an access function to map logical neighbour indices onto array indices in memory. This is why A and B need the function g to map the element logical grid coordinates i, j onto the actual array indices.

The set of coefficients, provided by the array C, is needed to modulate the respective contribution of the neighbour elements on the reduction result. The function f plays the same index mapping role for C as the function g for A and B, with two extra parameters however, to give the relative position of the coefficient in the stencil pattern.

We then define two operators to link the reduction expression components together. The coefficient modulation operator ⊗ applies a coefficient weight on the neighbour element contributions. The reduction operator ⊕ then accumulates the modulated contributions from each considered neighbour element.

g(i, j, B) =f (i, j, -1, 0, C) ⊗ g(i -1, j + 0, A) ⊕ f (i, j, 1, 0, C) ⊗ g(i + 1, j + 0, A) ⊕ f (i, j, 0, -1, C) ⊗ g(i + 0, j -1, A) ⊕ f (i, j, 0, 1, C) ⊗ g(i + 0, j + 1, A) (1)
Thus in order to describe a stencil, the framework needs the following information:

• Input data, output data, and coefficients arrays;

• Access function for data arrays;

• Coefficient relative coordinates;

• Extended access function for coefficients arrays;

• Operators for linking coefficients and reductions together.

The three input, output and coefficient arrays and the coefficient relative coordinates are directly application dependent, and must therefore always be provided by the programmer. Sensible defaults can be automatically provided for the access functions (C/C++ usual row-major array layout) and for building a default stencil expression together (using operator × for applying coefficients and operator + for computing the reduction by default). Such defaults may be overridden by the programmer if needed, to address more complex situations.

The following example shows how the stencil pattern of the equation 1 is created in the DSEL. It requires the only five following lines declaring the four variable coefficients and then assembling them to form a pattern (not related to any value).

coef_var2D<-1,0> c1; //left coefficient coef_var2D<1,0> c2; //right coefficient coef_var2D<0,-1> c3; //down coefficient coef_var2D<0,1> c4; //up coefficient auto st_pattern = c1+c2+c3+c4;

DSEL Design

The DSEL exposes C++ classes for the various coefficient flavors supported by the DSEL. The other classes used by the DSEL to abstract the stencil pattern (an aggregation of coefficients), the input and output information, and the whole expression resulting from applying a stencil on an input buffer and storing the results on an output buffer do not actually need to be exposed to the programmer. The auto type keyword and the type inference engine of modern C++ compilers are sufficient to properly determine the respective object types automatically.

Coefficients and stencil pattern Supported coefficients currently come in two flavors: variable coefficients such as presented in the example equation 1, and constant coefficients that are independent of the target element coordinates. In the first case, the only template parameter needed is the relative coordinates in the stencil pattern, (0, 0) being the current targeted element. The coefficient value will be defined elsewhere in a specific buffer, thus no type is needed in the constant coefficient template parameter since it only describes a pattern and not values. However internally to the DSEL, the coefficient evaluation will need a type, which is set later, when a stencil is associated to the typed coefficient buffer. In the second case, the constant coefficient value need to be passed as an argument to the constructor, whereas the expression type is needed as the last template parameter (float by default). No coefficient buffer will be needed in this case. coef_var2D<0, 0> c1; //variable coef_fxd2D<-1, 1 /*, float*/ > c2(0.2f); //constant Once several coefficients have been declared, they can be aggregated together to define a stencil pattern using the specially overloaded + operator. It is also possible to add several stencil patterns in a row, provided they are all of the same type, that is with the same number of dimensions and all constant or all variable coefficients.

/* creation of a 4 coefficient stencil*/ auto st1 = c1+c2; auto st2 = c3+c4; auto st = st1+st2; When a stencil is created with a single coefficient, a specific method is called on it to do the transformation. This case arises when building stencils such as copies or scaling transformations. coef_fxd2D<0,0> c_id {1.0f}; //is a coefficient auto st_id = c_id.toStencil(); //is a stencil Access functions Access functions correspond to the functions f and g in the equation 1. f and g abstract the data layout used in the buffers. Function g is for the input and output arrays. Function f is for variable coefficients, and is not required when using constant coefficients. It is declined in two versions, for input and output respectively. The input version may returns the element value directly, while the output version needs to return a writable, L-value compatible reference, to write the computation result into. When using the triSYCL library, these access functions can be passed as function pointers and recorded in the input and output objects. While this approach works for the C++-only triSYCL library, it is not actually compliant with the SYCL 1.2 revision which for-bids function pointers to be used within SYCL kernels -for good reasons such that a host function pointer would not have any meaning for a kernel executed on an accelerator device. Consequently, this illegal use of function pointers in a SYCL kernel is refused by the ComputeCPP compiler. In order to still enable the programmer to specialize the various access functions, we implemented a second, SYCL 1.2 compliant approach (and inspired indeed by the SYCL language itself) using the following template based function specialization scheme. The key aspect of this approach is the empty class passed as the second template parameter. This class does not need to be defined; it only needs to be declared and solely serve the purpose of discriminating multiple function specializations from each others, just like SYCL's own kernel naming classes [START_REF] Group | C++ single-source heterogeneous programming for OpenCL[END_REF]. The DSEL then accepts the naming class as a SYCL compliant alternative everywhere it expects a function pointer. Operation assembly and use To define an operation, input information must first be put together in a specific object containing the type of the element, the input buffer and associated access function.

One may also specify the coefficient buffer and its associated access function when using variable coefficients. The output information is built from the type, the output buffer and its access function.

output_2D<float, &outBuffer, &fdl_out> workload_out; input_var2D<float, &inBuffer, &coefBuffer, &fdl_in, &fac_coef> workload_in;

These pieces of information are assembled with the stencil to form the complete stencil expression. The programmer then only needs to declare a SYCL device queue and pass it to the DSEL for the parallel computation to be instantiated.

/* st is a predeclared stencil pattern */ auto op_work = workload_out << st << workload_in; op_work.doComputation(queue);

The drawback of this convenient syntax is that buffers must be defined statically. However, the DSEL could also support registering buffers through object constructors instead of template parameters, at the cost of disabling some static parameters inference.

DSEL workflow

While the DSEL enables the programmer to describe stencils thanks to the dedicated interface presented in the previous section, parallelization is performed automatically in a transparent manner. Thus, the DSEL workflow summarizes as follows.

Main programmer steps The programmer declares coefficients, aggregates them into a stencil pattern, defines function access to data buffers, assembles all these pieces of information together (including actual data buffers addresses) and eventually calls the computation method. The programmer does not have to specify parallelization-related directives or structures, with the exception of the device queue. All parallelization steps are handled internally.

Automatic parallelization Multiple stencil parallelization techniques exist such as color (e.g. red-black) element partitioning or tiling. Our DSEL is based on the tiling approach. For tiled stencil Figure 1. Stencil tile with a ghost region parallelism to work, a ghost boundary is necessary all around each tile. This ghost region supplies a read-only input replicate of the neighbourhood elements values, for updating the target elements located within the stencil width distance of the tile edges. An example tile and its ghost region are shown in Fig. 1, were the centre blue cells are the tile elements and the red cells are the ghost copies of the neighbour tiles elements. Once the internal computation of each tile is complete, the ghost region of each tile is updated with the newly computed element values of the neighbour tiles.

The DSEL automatically obtains the ghost zone size from the stencil width. After that it computes the parallelization parameters as the size of a workgroup, considering the data type and the number of coefficients. Using these parameters, it builds a specialized parallel SYCL kernel, which can then be further optimized by the SYCL compiler, and then by the back-end OpenCL compiler in accordance with the targeted device.

Implementation Details

The DSEL is written in C++14 and is made of some public and private templated classes. For now only around 800 lines of code are needed to implement a basic DSEL version. The last C++14 revision is used following the triSYCL library approach. It enables useful features such as the auto keyword. Keyword auto leverages the compiler type inference capabilities to relieve the programmer from writing fastidious template types, of which our DSEL makes a large use to automatically compute stencil execution parameters. The present contribution is inspired in this meta-programming approach by Blitz++, developed by T. Veldhuizen. Veldhuizen especially introduced the expression template technique [START_REF] Langer | C++ expression templates[END_REF][START_REF] Veldhuizen | Expression templates[END_REF] used in our DSEL to store and evaluate the stencil pattern. He also made some recommendations in his articles [START_REF] Veldhuizen | Guaranteed optimization for domain-specific programming[END_REF][START_REF] Veldhuizen | Active libraries: Rethinking the roles of compilers and libraries[END_REF], as the active library notion which leads to use the partial specialization technique [START_REF]Partial template specialization[END_REF].

Stencil parallelization through SYCL

The grid where the stencil is applied is partitioned in rectangular tiles to be processed potentially in parallel. Following the usual OpenCL terminology, each block is updated by a workgroup and each workitem in a workgroup updates a single element. Moreover each input element will be read as many times as there are points in the stencil pattern during the stencil computation. Therefore a device local memory area, shared among a workgroup, is used to temporarily cache the input elements to read. First, each workitem loads a subset of the global array of elements, corresponding to the tile elements to be updated and to the tile ghost domain, into the potentially faster local memory (depending on the actual device). Second, this local memory is accessed to compute the stencil kernel on the tile elements. The local memory enables to reuse the tile input values multiple times efficiently during the tile update, while the potentially slower global memory has been read only once for that tile.

To compute an element, the stencil kernel needs the input buffer, as well as the coefficient array when using variable coefficients. Since such coefficients may be different in every direction, no reuse is possible, and the coefficient array is left in global memory. When using constant coefficients, such coefficients are directly inserted in the inlined stencil evaluation function.

Currently the input and output buffers must have the same size and the stencil is applied to the whole data area, writing the results at the same location in the output buffer. Thus, we assume for simplicity, for now, that the output buffer has been allocated with the same ghost domain as the input buffer. Internally the DSEL computes the workgroup sizes independently from the global input/output size. If the global size is not a multiple of the workgroup size, border tiles on the right and at the bottom of the grid are made smaller and a condition disables workitems that would get the out of bound elements to compute.

Static Parameters Computing with Meta-programming

Static parameters being computed at compile time allow highly specialized host binaries and may dramatically improve the run time when used in parallel sections. One can argue that for GPUs the kernel code goes through another device dependent compilation step at run time, beyond the SYCL compilation step, so the effort may not be necessary. However some SYCL implementation such as the CPU-only triSYCL does not involve the use of OpenCL and therefore can benefit from having the native C++ compiler perform full specialization and static parameter computing. Moreover metaprogramming in C++ allows to use these static parameters to raise compilation errors. It helps the developer by detecting structural stencil errors at compile time (if the user attempts to mix fixed and variable coefficients for example), and giving a specific, accurate, and more user-friendly error message.

Static parameters are divided in two categories: those related to the stencil and those related to the parallelization technique. The first one contains for now the stencil width in each dimension direction (and thus the ghost domain width too). The second one is the workgroup size used in SYCL specific code. A distinct metaprogramming algorithm is used for each category: a constructive one based on expression templates for inferring the stencil width, and a recursive one using partial specialization for computing the workgroup size.

Stencil width inference Each 2D coefficient/stencil has four static members: the maximum relative position in the stencil for each direction (two times the two dimensions). Here we will denote them by lef t, right, top and bottom for convenience. For a coefficient with relative position (i, j) these parameters are simply equal to i for the first dimension (so for lef t and right) and j for the second one (top and bottom). When some coefficients are "added" in order to form a stencil pattern, a binary tree is built internally using the expression template technique. This way, each node of the tree is a stencil subpattern and can statically compute the local four parameters depending on their children values. It corresponds to the minimum of the two children' values for the negative directions (i. e. lef t and top), and to the maximum of these values for the positive directions (right and bottom).

The four parameters are needed to compute the ghost size: the needed thickness of the ghost regions in a direction is simply the difference between the dimension maxima ( right -lef t for example). The offset for accessing the updatable (non ghost) values in the input/output buffer is then (max(-lef t, 0), max(-top, 0)). Moreover this technique handles off-center and not symmetric stencils. It also would enable to specialize code with respect to the pattern layout, in a future version of the DSEL. Indeed the pattern building order influences the tree shape. This tree shape controls the input value processing, which are traversed in the same order. For example the two trees and related stencil pattern creation code depicted in Fig. 2 would lead to different performance. Workgroup size computation Under the assumption of some machine dependent environment variables available at compile time, the workgroup size can be statically computed w.r.t. the target machine maximum workgroup size in each dimension, and also w.r.t. the maximum local memory amount shared in a workgroup. This workgroup size corresponds to the tile size excluding the ghost elements. Notice that for now only one target is enabled in the DSEL implementation.

To do so the current naive algorithm increments the workgroup size by one element in each direction until a stop condition is reached. The stop conditions are the local memory overflow (since thanks to the stencil width the amount of accessed elements in a tile is exactly known) and the maximum workgroup size (which stop the workgroup expansion only in the related dimension). This way the workgroup can have a rectangular or square shape in two dimensions. A static_assert checks if the stencil width is larger than the maximum workgroup size at the algorithm beginning. This computations are performed statically using meta-programming with a tail-recursive algorithm.

Discussion

As mentioned above, we implemented two variants of the stencil DSEL API. The first variant enables the programmer to specify access function as function pointers in template parameters, but it only works for triSYCL. The second variant fully complies with the SYCL 1.2 specification (and with both triSYCL and ComputeCPP implementations), and so does not allow these function pointers, since they cannot be used with OpenCL on which SYCL is based. Although the second version is the one that would be chosen if the DSEL implementation continues, we think that the first version may offer some more benefits for the user: it requires less lines of code and is more traditional in terms of programming style.

The current DSEL prototype suffers from other restrictions as the one depicted at the end of section 3.3. Thus the buffers have to be created statically, and the benefit is to enable future versions of the DSEL to detect if the input and output buffers reference the same object. However to impose static buffers may be a too strong restriction, which may be removed in future versions, together with the static workgroup size computation using environment variables available at compile time, which restricts portability.

Finally the implementation details allow us to classify our DSEL among the two following categories: shallow and deep [START_REF] Gibbons | Folding domain-specific languages: Deep and shallow embeddings[END_REF]. Actually the DSEL belongs to both, depending on which aspect we consider in a stencil. The DSEL is shallow regarding the coefficient operators (⊗ and ⊕), which are not stored internally but directly inferred by the compiler (the user can directly override the * and + operators in the stencil elements' type class). On the opposite, the DSEL is deep regarding the stencil pattern, which is stored like a tree and so defines a kind of syntax.

Conclusion and Future Work

The DSEL introduced in this paper allows to express with ease stencil computation, that is automatic parallelized on both CPU and GPU architectures. The stencil concise abstraction has been made possible thanks to modern C++ features and especially to metaprogramming techniques. Moreover, the parallelization is based on the new SYCL interface which offers a unified parallel programming mean for heterogeneous platforms. The benefits of such a DSEL are to ease the developer work, which then neither has to be a parallelization specialist, nor needs to use several different languages/compilers, or be restrained to a single processing unit type.

However, this DSEL is still a proof of concept at this stage, so many improvements and further experiments are needed on production release SYCL compilers/libraries instead of the preliminary version that have been used for this paper. Among our intended improvements, the priority will be put on the DSEL generalization to any n-dimensional stencil (probably using variadic templates [START_REF] Gregor | Variadic templates for c++0x[END_REF]) and on extending the code specialization to both target parameters and modern tiling techniques (e.g. [START_REF] Bandishti | Tiling stencil computations to maximize parallelism[END_REF]).

/

  * output access function, notice the & */ float& fdl_out(int a,int b, accessor<float, 2, access::mode::write, access::target::global_buffer> acc) {return acc[a][b];} /* input access function */ float fdl_in(int a,int b, accessor<float, 2, access::mode::read, access::target::global_buffer> acc) {return acc[a][b];} /* coefficient access function */ float fac_coef(int a,int b, int c, int d, accessor<float, 1, access::mode::read, access::target::global_buffer> acc) {return acc[0];}

  class out_access; //use chosen class name template<> inline float& fdl_out<float, out_access> (...) {...}

Figure 2 .

 2 Figure 2. Two ways to create a stencil pattern from coefficients, and their related internal tree representation.

Acknowledgments

We are especially grateful to the Codeplay company for providing us with a preliminary version of ComputeCPP and with swift, efficient support throughout the study. We are also very much grateful towards Ronan Keryell for implementing the triSYCL library which help kick-start this work, and for contributing valuable comments on our stencil DESL design.