Aurelie De

Rus Jacquet

Sophie Laye

Frederic Calon

Aurelie De Rus Jacquet

Fre ´de ´ric Calon

How nutrients and natural products act on the brain: Beyond pharmacology

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The tangled tree of phytomolecules and metabolites Knowingly or not, each and every one of us depends on plants for survival, from the food we eat, the air we breathe, to the clothes we wear. Throughout their lives, humans will interact with hundreds of botanical species, and they will rely on the nutritive and medicinal properties of thousands of natural compounds, such as the large family of polyphenolsmore than 8,000 different chemical structures are known to date. To appreciate the medical and scientific importance of phytomolecules, we must put them into context. Health Canada defines natural health products as ''naturally occurring substances that are used to restore or maintain good health.'' However, their potential may be greater than currently appreciated. The plant kingdom, similarly to the animal kingdom, is constantly exposed to environmental stressors. Therefore, plants have evolved to produce thousands of chemically diverse secondary metabolites with potent biological activities, such as antibiotics, antiparasitics, and other biochemicals that repel predators. They also manufacture phytochemical compounds to promote cellular health and survival in situations of environmental stress. In the context of neurodegenerative diseases, it is the latter that we attempt to leverage for neuroprotective purposes. These protective molecules retain their pro-health biological properties when ingested by mammalian organisms, but they undergo addi-tional chemical modifications by the gut microbiota before crossing the epithelial barrier to reach the circulation. These molecules are further modified by firstpass metabolism in the liver and phase II enzymatic conversion (Figure 1). The resulting bioavailable metabolites form a large new pool of bioactive compounds and, together with their parent compounds, further enrich the diversity of plant-derived phytomolecules. These complex metabolic patterns and the resulting molecular diversity support a conceptual framework whereby the individual or synergistic actions of metabolites modulate central and peripheral biological targets to sustain brain health.

Biodistribution of phytomolecules and their metabolites: Where are the cellular targets?

The health benefits of bioactive molecules depend on their ability to reach target cells, tissues, or organs in sufficient concentrations. In the context of neurodegenerative diseases, the primary objective is to increase survival of vulnerable neuronal populations to stop or slow disease progression. Puzzling observations show that dietary interventions result in measurable neuroprotective outcomes in animal and clinical studies despite low apparent brain bioavailability. [START_REF] Scarmeas | Nutrition and prevention of cognitive impairment[END_REF] The concept that bioactive polyphenols and other phytochemicals may not necessarily need to reach the brain tissue to induce neuroprotective benefits remains incredibly relevant and-for the most part-unanswered. Circulating metabolites may affect the central nervous system (CNS) via several mechanisms, from crossing the blood-brain barrier (BBB) and directly acting on the brain parenchyma to modulating the function of peripheral targets that will indirectly benefit brain health 2 (Figure 1). Proposed target 1: Brain microvascular endothelial cells Brain microvascular endothelial cells (BMECs) are unique gatekeepers of the cerebrovascular system that help maintain brain homeostasis and are key components of the neurovascular unit (NVU). Capable of integrating complex environmental cues, they form a dense vascular network that irrigates the brain parenchyma. The biological function of these highly specialized cells is driven by their key positioning at the interface of the blood and the brain, thus acquiring a dual luminal vs. abluminal polarity and accompanying molecular characteristics. [START_REF] Worzfeld | Apicobasal polarity of brain endothelial cells[END_REF] For example, polarized membrane lipid composition, localization of receptors and transporters, release of signaling proteins, and response to extracellular cues differ between the luminal and abluminal sides, [START_REF] Worzfeld | Apicobasal polarity of brain endothelial cells[END_REF] suggesting that blood-brain communication is tightly controlled by the endothelium. Circulating natural compounds and their metabolites naturally come into contact with the vasculature, and one could argue that downstream pro-health effects may result from a BMEC-mediated response triggered by these molecules instead of a direct action on CNS target cells, as recently shown with insulin. [START_REF] Leclerc | Cerebrovascular insulin receptors are defective in Alzheimer's disease[END_REF] Physiological barriers are living cellular units and signaling hubs at the interface of tissues and peripheral circulation. Fascinating new studies revealed the high connectivity existing between cells forming the endothelium, which are organized into synchronized cell clusters that respond to biochemical stimuli and rapidly propagate local signals via network connections. [START_REF] Lee | Small-world connectivity dictates collective endothelial cell signaling[END_REF] signals can be endogenous, like circulating hormones, neurotransmitters, or immunomodulating factors, but we propose that they could also be exogenous molecules such as natural compounds. BMECs are equipped to interact with and process these exogenous factors, as demonstrated by their expression of metabolizing enzymes (e.g., cytochrome P450 family), which suggests an active interplay between circulating factors and the brain vasculature. Overall, much research remains to be accomplished to shed light on how phytomolecules may act on BMECs to improve hallmarks of cerebrovascular health and trigger the release of neuroprotective molecules into the parenchyma. Proposed target 2: Cells of the brain parenchyma A small subset of phytomolecules and their derivatives has the ability to cross the BBB and reach target cells such as neurons and glia. While the specific mechanisms by which these bioactive compounds enter the brain remain to be fully explored, small phenolic metabolites with hydrophobic functional groups and displaying low hydrogen-bonding potency may cross the BBB without specialized transport mechanisms. [START_REF] Vela ´squez-Jime ´nez | Phenolic compounds that cross the blood-brain barrier exert positive health effects as central nervous system antioxidants[END_REF] However, observations regarding the possible stereoselectivity of natural compounds from the polyphenol family (e.g., epicatechin vs. catechin) suggest that some of these molecules may also interact with specialized endothelial transporters. The cellular membranes of BMECs harbor a diversity of uptake and efflux transporters, such as solute carriers (SLCs) and ABC transporters (e.g., p-glycoprotein) that regulate entry of circulating molecules. While the role of the BBB as a functional barrier separating the blood and the brain is conserved among mammals, brain microvessels display species-specific transcriptomic profiles, [START_REF] Song | Transcriptomic comparison of human and mouse brain microvessels[END_REF] and new humanbased in vitro models may better model conditions specific to human diseases. Brain molecular atlases revealed the unique signature of various regions involved in neurodegenerative disorders, and different nutrients may have specific tropisms toward defined brain regions. [START_REF] Lauriola | Anatomical biology guides a search for nutrients for the aging brain[END_REF] While this field of investigation is still at its infancy, better understanding the mechanisms of entry of natural products and metabolites into the brain parenchyma, identifying their tropism toward neuronal and glial targets, and defining the potential impact of regional differences would further support these molecules as sources of neuroprotective agents. Proposed target 3: Peripheral organs Paradigm-shifting discoveries on the dynamics of blood-brain exchanges have recently emerged, and the realization that the peripheral system acts as a significant modulator of brain health is gaining momentum. [START_REF] Pluvinage | Systemic factors as mediators of brain homeostasis, ageing and neurodegeneration[END_REF] It is also increasingly recognized that neurodegenerative, inflammatory, and cardiometabolic disorders share many risk factors and pathophysiological pathways. In this context, it is reasonable to suspect that beneficial phytomolecules and their metabolites induce protective effects on peripheral organs more readily accessible than the brain, and by altering the composition of the systemic circulation, this could further translate into neuroprotective outcomes. The gut-microbiota-brain axis is among the peripheral elements that could play a protective role and possesses intriguing properties. In addition to being the first and only organ to interact with the complete pool of ingested phytomolecules, the gut has a privileged mode of communication with the brain via afferent spinal and efferent vagal nerves (Figure 1). Moreover, the gut microbiome produces neuromodulating molecules, mediates stress re-sponses, releases factors into the systemic circulation, and, depending on the nature of its resident microorganisms, is associated with distinct behavioral outcomes. [START_REF] Liu | Microbiota and the gut-brain-axis: Implications for new therapeutic design in the CNS[END_REF] Since food intake is a significant factor regulating the composition of the gut microbiome, beneficial dietary compounds could induce neuroprotective effects by reshaping the ecosystem of this peripheral organ. [START_REF] Liu | Microbiota and the gut-brain-axis: Implications for new therapeutic design in the CNS[END_REF] Similarly, dietary interventions aiming at improving cardiometabolic health are currently being investigated in various neurodegenerative disorders.

Methodological considerations to study the benefits of natural products in neurodegenerative diseases

Randomized clinical trials Randomized clinical trials (RCTs) are the gold standards to make health recommendations. However, when exploring the benefits of natural health products, we must distinguish between dietary trials involving complex food matrices and purified phytocompounds with defined biochemical activities. Over the years, population-based association and cohort studies have convincingly shown a relationship between brain health and nutrition; however, translation into practice has been hampered by challenges in determining causal relationships between dietary interventions and health outcomes. In contrast to RCTs evaluating the effect of a single molecule, dietary interventions are uniquely affected by confounding factors such as the participants' baseline dietary practices, inherent variations in the nutritional value of food products, individual genetic variations in metabolic responses, and long-term trial durations that reduce compliance to the defined dietary intervention, which are all difficult to control with a placebo group. [START_REF] Mirmiran | Common Limitations and Challenges of Dietary Clinical Trials for Translation into Clinical Practices[END_REF] Nonetheless, efforts aimed at the identification of dietary biomarkers associated with improved brain health could lead the way to rigorous mechanistic experimentation in clinical settings. Once identified, these biomarkers would allow the Ingested compounds are metabolized by the gut microbiome. The parent molecules are retained in the gut, where they can modulate the microbiota-gut-brain axis by communicating through the vagus nerve. The metabolites are transported to the liver via the portal vein, where they undergo first-pass metabolism and enter the systemic circulation to reach the brain and peripheral organs. Metabolites could interact with BMECs through transporters, receptors, and metabolizing enzymes that remain to be confirmed and fully characterized. Abbreviations: BMEC, brain microvascular endothelial cell; CNS, central nervous system.

Cell Reports Medicine 4, 101243, October 17, 2023 3 implementation of depletion-repletion paradigms to validate a direct cause-effect relationship between nutrients and the health outcome under investigation. [START_REF] Brickman | Dietary flavanols restore hippocampal-dependent memory in older adults with lower diet quality and lower habitual flavanol consumption[END_REF] Difficulties in securing intellectual property have also been a strong deterrent to the full development of unpatentable interventions. [START_REF] Calon | Nonpatentable drugs and the cost of our ignorance[END_REF] While these factors may appear as challenges to conduct pricy clinical trials, they actually offer a wealth of opportunities to modulate clinically relevant disease phenotypes, such as the gut-microbiome-brain axis, and pursue nutritional recommendations to promote brain health. Still, many RCTs have been conducted in the last 15 years, suggesting active efforts to improve the well-being of patients through plant-based natural health products. [START_REF] Scarmeas | Nutrition and prevention of cognitive impairment[END_REF][START_REF] Leclerc | Can Natural Products Exert Neuroprotection without Crossing the Blood-Brain Barrier?[END_REF] In vivo models Prior to conducting clinical trials, animal models are key to establishing natural product oral bioavailability, i.e., prediction of absorption, distribution, metabolism, and excretion along with pharmacokinetic (PK) parameters. [START_REF] Leclerc | Can Natural Products Exert Neuroprotection without Crossing the Blood-Brain Barrier?[END_REF] An important aspect to consider for disease-modifying studies involving neurodegenerative disorders is the ability of compounds to cross the BBB to reach neuronal and glial targets. Considering the range of targets and mechanisms of action detailed above, it becomes essential to determine whether a molecule of interest enters the brain parenchyma to reach the minimum effective concentration (MEC) established for this compound. If the brain MEC is not achieved despite clear neuroprotective benefits, biological targets are likely to fall into one of the categories discussed earlier and indirectly improve brain health. To inform on potential neuroprotective mechanisms and CNS penetration, advanced in vivo methodological approaches such as intracarotid in situ cerebral perfusion enable the quantification of kinetic variables including brain uptake and permeability constants. This technique was developed in the 1980s, but it has been underused despite major advantages over other approaches. In this model, the right common carotid artery is catheterized after ligation of the external branch and perfused with the molecule(s) of interest, which directly reach the brain vasculature and bypass the peripheral metabolism to maintain the perfused molecule(s) in their native forms (Figure 2). As a result, transport Advanced experimental approaches used to investigate the BBB in vitro and in vivo will enable to pinpoint the specific mechanisms involved in brain distribution and neuroprotective benefits of natural health products. Abbreviations: 3D, three dimensional; BMEC, brain microvascular endothelial cell; ECM, extracellular matrix.

from the blood to the brain via the BBB can be efficiently quantified while requiring only a small number of animals to complete the study. [START_REF] Leclerc | Cerebrovascular insulin receptors are defective in Alzheimer's disease[END_REF] This direct measure of brain penetration is more informative compared with other approaches that rely on surrogate measures of compound permeability, for example in the cerebrospinal fluid. These BBB-specific techniques are key add-ons to more standard approaches to determine the PK and biodistribution of natural or pharmaceutical compounds. In essence, animal studies of natural product distribution and brain permeability are necessary to inform on the mechanisms by which bioactive molecules improve brain health and function. However, most animals have feeding habits and nutrient needs that can be very different from humans, and interspecies differences in metabolism and endothelial cell molecular composition and function could in vivo findings and limit translation of the results. Therefore, complementary in vitro studies using human-derived cells could be implemented.

In vitro models Advances in microfluidic and organoid technologies have recently allowed the generation of BBB-in-a-dish platforms that recapitulate tissue-level organization and enable the deep investigation of cellular and molecular mechanisms underlying neuroprotection (Figure 2). [START_REF] De Rus | The contribution of inflammatory astrocytes to BBB impairments in a brain-chip model of Parkinson's disease[END_REF] These approaches leverage the differentiation of human iPSCs into BBB-related cell types or utilize human primary cells to reconstruct the three-dimensional organization and blood-brain compartmentalization of this multicellular unit. These models have the distinct advantages of being composed of human cells originating from healthy donors or people with the neurodegenerative disease of interest, being able to include biological sex as a variable, and being suitable to screening experiments. In contrast to two-dimensional models, microfluidic channels mimic capillary shear stress, a mechanical pleiotropic regulator of BMEC function. Together with glial/ neuronal signals originating from the brain compartment, laminar flow promotes physiological BMEC polarization and reorganization into a single-layer three-dimensional vessel with luminal and abluminal surfaces. Therefore, BBB chips are ingenious tools to study signaling at the BBB, in particular how natural compounds and metabolites may cross the endothelial layer or propagate signals from within the vasculature. For instance, studies can estimate the mechanisms of brain penetration by measuring the passage of bioactive molecules across the endothelial barrier via specialized transporters and evaluate how synergistic effects between molecules may promote transport across the BBB. The next key consideration is to define the tropism of natural compounds toward specific cell type(s). While drug discovery for neurodegenerative disorders has mostly focused on neuronal targets, accumulating evidence suggests that glial cells and cellular dysfunctions at the neuro-glia-vascular unit could contribute to neuronal loss. [START_REF] De Rus | The contribution of inflammatory astrocytes to BBB impairments in a brain-chip model of Parkinson's disease[END_REF][START_REF] De Rus | The LRRK2 G2019S mutation alters astrocyte-to-neuron communication via extracellular vesicles and induces neuron atrophy in a human iPSC-derived model of Parkinson's disease[END_REF] In this context, natural products may be at an advantage, as they have been shown to modulate multiple protective pathways in both neuronal and non-neuronal targets. Furthermore, the discovery of new cellular subtypes, particularly at the BBB, is a welcome reminder of the diverse ecosystem in which beneficial compounds may exert their biological properties. However, differences in molecular and functional profiles between in vitro cultures and their in vivo counterparts need to be considered, as changes/adaptations to the in vitro microenvironment could alter native cellular responses.

Ethics, diversity, equity, and inclusion considerations

From an ethical standpoint, natural products and synthetic drugs should be valued equally when searching for neuroprotective molecules. The chemical complexity of natural molecules is often greater than that of synthetic compounds, and they can be manipulated by medicinal chemists to improve desirable biophysical properties. Importantly, botanical sources of phytomolecules and nutrients are found globally from sustainable sources and are not confined within cutting-edge laboratories, thus potentially allowing a greater access of medicinal resources to the global population. Greater collaborations between neuroscientists, pharmacists, and ethnopharmacologists are one way by which we can produce impactful biomedical scientific discoveries. Better understanding the properties of natural products originating from the diet or pharmacopeia of populations generally underrepresented in the sciences could not only unveil new horizons in our search for neuroprotective candidates but also promote ideals of diversity, equity, and inclusion. The North American and European patent systems have greatly favored the discovery of new synthetic scaffolds over botanical sources of natural products, and this legislative ecosystem may explain the relative neglect of phytochemicals in the second half of the 20 th century. In contrast, the blooming multi-billion-dollar market of natural supplements is clear evidence that people long for non-medical approaches to health. We postulate that natural health products are not meant to replace synthetic drugs or become ''alternative medicines.'' Instead, they provide a diverse arsenal of molecules that can work together with synthetic medicines.

Concluding remarks

The complexity of how nutrients and natural products act on the brain calls for the development of innovative conceptual and methodological approaches that consider that several tissues and targets could be affected by the beneficial compound(s). The study of neuroprotective molecules is further complexified by the restrictive BBB, but the vascular network is a complex biological system integrating and disseminating signals from the blood to the brain, and it should be envisioned as a potential key player in mediating CNS benefits. The controversy of where neuroprotective natural compounds act to improve brain health is a reminder of the complexity of body-brain interactions, and it questions conceptual frameworks that tend to uncouple peripheral and central biology.

Cell Reports Medicine 4 ,

 4 101243, October 17, 2023 ª 2023 The Author(s). 1 This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Figure 1 .

 1 Figure 1. Overview of the potential mechanisms involved in the neuroprotective benefits of plant-derived compounds

Figure 2 .

 2 Figure 2. Advanced methodological approaches to study the blood-brain barrier

(legend on next page)2 Cell Reports Medicine 4, 101243, October 17, 2023

Cell Reports Medicine 4, 101243, October 17, 2023

Cell Reports Medicine 4, 101243, October 17, 2023 5

Cell Reports Medicine 4, 101243, October 17, 2023

ACKNOWLEDGMENTS A.d.R.J. is supported by the Parkinson's Foundation (PF-Launch-937655). S.L. is supported by INRAE (AlimH), Universite ´de Bordeaux, Eranet-Neuron Gut2Behave (ANR-19-NEUR-0003-03), JPND SOLID (ANR-21-JPW2-0004-05), and Nouvelle Re ´gion Aquitaine ExoMarquAge (13059720-13062120). F.C. is a Fonds de recherche du Que ´bec-Sante ´research scholar and is supported by the Canadian Institutes of Health Research. This work was also support by internal funds from INAF to A.d.R.J. and F.C.

DECLARATION OF INTERESTS

F.C. is co-inventor on a provisional patent application: CA3002753A (World, PCT: WO2017072219A1).


Wiktor Walasik

Gilles Renversez
email: gilles.renversez@fresnel.fr

Fangwei Ye

Plasmon-soliton waves in planar slot waveguides: II. Results for stationary waves and stability analysis

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

We describe the results of the two methods we developed to calculate the stationary nonlinear solutions in one-dimensional plasmonic slot waveguides made of a finite-thickness nonlinear dielectric core surrounded by metal regions. These two methods are described in detail in the preceding article [Walasik et al., submitted to PRA]. For symmetric waveguides, we provide the nonlinear dispersion curves obtained using the two methods and compare them. We describe the well known low-order modes and the higher-modes that were not described before. All the modes are classified into two families: modes with and without nodes. We also compare nonlinear modes with nodes with the linear modes in similar linear slot waveguides with a homogeneous core. We recover the symmetry breaking Hopf bifurcation of the first symmetric nonlinear mode toward an asymmetric mode and we show that one of the higher modes also exhibits a bifurcation. We study the behavior of the bifurcation of the fundamental mode as a function of the permittivities of a metal and a nonlinear core. We demonstrate that the bifurcation can be obtained at low power levels in structures with optimized parameters. Moreover, we provide the dispersion curves for asymmetric nonlinear slot waveguides. Finally, we give results concerning the stability of the fundamental symmetric mode and the asymmetric mode that bifurcates from it using both theoretical argument and numerical propagation simulations from two different full-vector methods. We investigate also the stability properties of the first antisymmetric mode using our two numerical propagation methods.

I. INTRODUCTION

Nonlinear slot waveguides (NSWs) are structures made of a finite-size nonlinear dielectric layer sandwiched between two semi-infinite metal layers. They have been studied in Refs. [1,2] where it was shown that they allow for sub-wavelength confinement of light and phase matching for the second-harmonic generation. More recently, in the works of Rukhlenko et al. [3,4], analytical formulas for the dispersion relations of these NSWs were presented for symmetric and antisymmetric nonlinear modes only. These dispersion relations were given using integral equations that have to be solved numerically. The study of Davoyan et al. [5] showed, using the numerical shooting method to solve Maxwell's equation in NSWs, that a symmetry breaking bifurcation that generates an asymmetric mode from the fundamental symmetric mode occurs in NSWs. Such bifurcation phenomena are well known in fully dielectric nonlinear structures [6][7][8][9][10][11][12][13]. Recently, higher order modes were also reported in NSWs [14]. Moreover, it was shown that plasmonic coupling and symmetry breaking phenomena can be observed in waveguides built of a linear dielectric core sandwiched by nonlinear metals [15,16]. Nonlinear switching was predicted in NSW-based structures using numerical simulations [17].

In the preceding article [18], we describe two models we developed to study the stationary nonlinear solutions in NSWs where the nonlinear core of the focusing Kerr type was considered. The first model assumes that the nonlinear term depends only on the transverse component of the electric field and that the nonlinear refractive index change is small compared to the linear part of the refractive index. It allows us to describe analytically the field profiles in the whole waveguide. It also provides a closed analytical formula for the nonlinear dispersion relation. This first model is called Jacobi elliptic function based model (JEM). The second model takes into account the full dependency of the Kerr nonlinear term on all electric field components and no assumption is required on the amplitude of the nonlinear term. The disadvantage of this approach is the fact that the field profiles must be computed numerically. This second model is called the interface model (IM). This article is organized in the following way. In Sec. II, we describe the results obtained with our two models for symmetric NSWs. They include a mode classification taking into account the higher order modes we found previously [14] and a detailed study of the field profile transformation as a function of power. We also provide a permittivity contrast study that allows us to decrease by several orders of magnitude the bifurcation threshold at which the first asymmetric mode appears. In Sec. III, we provide the results concerning asymmetric NSWs in which the mode degeneracy is lifted. Finally, in Sec. IV, using both theoretical arguments and numerical propagation simulations from two different full-vector methods we provide results on the stability of the main stationary solutions obtained in symmetric NSWs.

II. RESULTS FOR SYMMETRIC WAVEGUIDES A. Dispersion relations, field profiles and mode classification

In this section, the dispersion relations obtained for the symmetric NSW are presented. The field profiles corresponding to each of the dispersion curves are depicted and allow us to classify the modes according to their symmetry and the number of nodes in the magnetic field profile.

Figure 1 presents dispersion relations for the symmetric NSW obtained using the JEM and the IM. The parameters of the NSW studied are: 1 = 3 = -90 (gold), l,2 = 3.46 2 , α 2 = 6.3 • 10 -19 m 2 /V 2 (hydrogenated amorphous silicon) and d = 400 nm at a wavelength λ = 1.55 µm. The geometry of the structure with its parameters is shown in Fig. 2. The dispersion relations present the dependence of the effective index of the mode β as a function of the power density in the waveguide core P c which is calculated in the following way:

P c = d 0 S z dx, (1)
where S z denotes the z-component of the Poynting vector S = 1/2 e(E × H * ).

We observe a very good qualitative agreement between the dispersion diagrams obtained using our two models. The number and the character of the dispersion curves is very similar in both cases. The qualitative agreement between the results of the two models confirms their validity. Quantitatively speaking, the models agree in the range of low power densities (below 10 9 W/m). Above this value we observe quantitative differences in the results. The origin of the differences is explained by the assumptions made in the JEM (low nonlinearity, only E x component of the electric field contributing to the Kerr nonlinear effect) as described in Ref. [18]. In the following, we will focus on the results obtained using the more accurate IM.

The NSW supports numerous modes with various properties. Firstly, we will discuss the mode classification according to the symmetry of the mode. For the low power region, the NSW studied here supports two modes: a fundamental symmetric mode [blue curve labeled S0 in Fig. 1 and blue field profiles in Fig. 3] and a low-power antisymmetric mode [red curve labeled AN0 in Fig. 1 and red field profiles in Fig. 3]. At P c ≈ 10 9 W/m a symmetry breaking bifurcation occurs that gives birth to an asymmetric mode [5] [green curve labeled AS1 in Fig. 1 and green field profiles in Fig. 3]. This modes and this type of behavior are known in nonlinear waveguides [5-13, 15, 19]. The power density P c of the modes S0, AN0, and AS1 first increases with the increase of the effective index β and decreases for β > ∼ 4.75. Our models allow us to find new, higher order modes in NSWs. The higher order modes can be divided into two families: node-less modes and modes with nodes. Among the node-less modes we find higher order symmetric modes (SI and SII) from which asymmetric modes bifurcate (AS2 and AS3, respectively). Their dispersion curves are labeled with the name of the mode in Fig. 1 and their field profiles are presented in Figs. 4 and 5(e), (f). Higher order node-less modes resemble a single soliton (SI and AS2) or two solitons (SII and AS3) propagating in the NSW core.

All the dispersion curves of the asymmetric modes are FIG. 3. Profiles of magnetic field Hy(x) for the symmetric S0 mode (blue), antisymmetric AN0 mode (red), and the firstorder asymmetric AS1 mode (green). The subplots present the transformation of the field profiles at the points corresponding to the vertical lines labeled a-e indicated in Fig. 6.

doubly degenerate. This means that to one value of the effective index (and the corresponding power density) correspond two solutions localized on one of the two core interfaces [compare green curves in Figs. 3(a) and (e) for the AS1 mode, Figs. 4(a) and (e) for the AS2 mode and the green and the gray profiles in Fig. 5

(f)]

The higher order modes with nodes resemble the modes of a linear slot waveguide with a higher refractive index than the one used here (see Sec. II C). Only symmetric (S1, S2, . . . ) and antisymmetric (AN1, AN2, . . . ) modes with nodes exist. Their dispersion curves are presented in Fig. 1 and their field profiles are shown in Figs. 5(a)-(d). The dispersion curves of the modes with nodes start for β = 1 and their effective index grows with the increase of the power density P c .

In Fig. 6, we present the dispersion relations obtained using the IM in a different coordinate frame. This time we use the total electric field intensity at x = 0 (the interface between the NSW core and the metal cladding) E 0 {see Eq. ([equation] [28][]28) in Ref. [18]}. This quantity is one of the input parameters of the IM. The dispersion diagrams β(E 0 ) have a drastically different character from the β(P c ) diagrams presented in Fig. FIG. 4. Profiles of magnetic field Hy(x) for the symmetric SI mode (blue) and the second-order asymmetric AS2 mode (green). The subplots present the transformation of the field profiles at the points corresponding to the vertical lines labeled v-z indicated in Fig. 6. In each subplot the value of E0 {see Eq. ([equation] [28][]28) in Ref. [18] for its definition} is identical for both modes.

difference is caused by the fact that E 0 is a local quantity, whereas P c is a global quantity, that results from the integration over the core width. In the coordinates of E 0 , the dispersion curves of the asymmetric modes are not degenerate. In Fig. 6, we notice that for asymmetric modes, to a given value of β correspond two values of E 0 , that represent solutions localized on the left and right interface of the waveguide core.

In Fig. 3, the comparison of the field profiles of the three main modes is presented during their transformation along the dispersion curve associated to the increase of E 0 . The field profiles of the S0 and AN0 modes do not change qualitatively. On the contrary, the field profile of the AS1 mode undergoes a qualitative transformation.

For low E 0 values, this mode is highly asymmetric and strongly localized on the right core interface x = d. With the increase of E 0 the asymmetric profile becomes more symmetric, and at the point of bifurcation it perfectly overlaps with the symmetric mode [see Fig. 3(c)]. For E 0 values above the bifurcation point, the mode becomes asymmetric and it tends to localize on the left interface.

In Fig. 4 a similar transformation is shown for the SI and AS2 modes. Here, with the increase of E 0 , the peak amplitude of the soliton H peak = H y (x = d/2) first decreases (it is the lowest at the bifurcation) and then increases, while side lobe peak amplitude H lobe (located at x = 0 and x = d) of the symmetric mode increases monotonously with E 0 . The ratio H peak /H lobe first decreases (up to the bifurcation) and then increases. In case of the asymmetric mode AS2, with the increase of E 0 the soliton peak shifts from left to right. At the same time the amplitude of the left (right) side lobe increases (decreases). More results on the mode transformation, including the results obtained using the JEM are presented in Ref. [START_REF] Walasik | Plasmon-soliton waves in metal-nonlinear dielectric planar structures[END_REF].

- FIG. 5. Typical magnetic field profiles Hy(x), obtained using the IM, corresponding to different dispersion curves indicated in Fig. 6. Abbreviations next to the subfigure labels indicate the dispersion curve to which a given profile corresponds. The color of the profile allows us to distinguish the mode symmetry: symmetric (S -blue), antisymmetric (AN -red), and asymmetric (AS -green). For asymmetric doubly degenerate mode AS3 the second profile is shown in gray.

B. Single interface limit

In Sec.

[section] [4][]IV in Ref. [18], describing the theoretical derivation of the models for NSWs, we mentioned that in the limiting case, where the integration constants c 0 in Eqs. ([equation] [8][]8) and ([equation] [10][]10) or C 0 in Eqs. ([equation] [START_REF] Kivshar | Optical Solitons, From Fibres to Photonic Crystals[END_REF][]26) and ([equation] [33][]33) (all equation numbers correspond to Ref. [18]) are equal to zero, we recover the case of a single interface between a metal and a nonlinear dielectric. Looking at the field profiles of highly asymmetric modes AS1 [see Figs. 3(a), (e)], we see that these modes are mostly localized at one interface only. Therefore, they can be well approximated by a solution of the single-interface problem.

In Fig. 7, we present the dispersion curves for the NSW obtained using the JEM [β(H 0 )] {see Eq. ([equation][1][9]9a) in Ref. [18]} and the IM [β(E 0 )] (compare with Fig. 6). Additionally to the antisymmetric (red), symmetric (blue), and asymmetric (green) dispersion curves, black dispersion curves obtained using single-interface models are presented. In the case of the JEM, the single-interface approximation was obtained using the 'field based model' for configurations with semiinfinite nonlinear medium described in Ref. [START_REF] Walasik | [END_REF]. This model was used for a single interface between a metal and a nonlinear dielectric with the same parameters as our NSW. In case of the IM, the corresponding singleinterface approximation was obtained using the 'exact model' for configurations with a semi-infinite nonlinear medium described in Ref. [START_REF] Walasik | [END_REF].

In Fig. 7, we see that both for the JEM and for the IM, the single-interface dispersion curve always lays between the antisymmetric AN0 curve and the symmetric S0 curve. For high values of E 0 , the asymmetric AS1 curve becomes very close to the black curve, but remains slightly above it. The fact that the black curves overlap with the green AS1 curves confirms that the highly asymmetric AS1 modes (for high effective index β) are well approximated using the single-interface approach. Now, instead of using the corresponding models for configurations with a semi-infinite nonlinear medium, we will use the formulas found in Sec. [section] [4][]IV in Ref. [18] that give us the analytical expressions for the dispersion relations for the single-interface problem. In the case of the JEM, the analytical formula for the dispersion relation of a single metal/nonlinear dielectric interface problem is given by Eq. ([equation] [START_REF] Oskooi | [END_REF][]42) in Ref. [18]. In this equation, as in the entire formulation of the JEM, the primary parameter is the magnetic field amplitude at the interface H 0 . Therefore, we are able to show the dependency described by Eq. ([equation] [START_REF] Oskooi | [END_REF][]42) only in the coordinates where the effective index is presented as a function of the magnetic field amplitude at the interface H 0 [see yellow line in Fig. 7(a)]. We observe that the dispersion relations calculated using the field based model (balck curve) and the yellow curve described by Eq. ([equation] [START_REF] Oskooi | [END_REF][]42) overlap perfectly. The single-interface dispersion curve, which corresponds to the limiting case c 0 = 0 divides the dispersion plot β(H 0 ) into the regions corresponding to the nodeless family and the family with nodes as predicted in Section [subsection] [2][3]III B in Ref. [18] Above the c 0 = 0 curve (for negative values of the integration constant c 0 ), only node-less solutions exist. Below the c 0 = 0 curve (for c 0 > 0), only solutions with nodes exist.

In case of the IM, the analytical formula for the dis-persion relation for the single-interface problem is given by Eq. ([equation] [START_REF] Rodriguez | [END_REF][]46) in Ref. [18]. The effective index of the mode expressed as a function of the material parameters of the structure and the total electric field intensity at the interface E 0 . The curve described by Eq. ([equation] [START_REF] Rodriguez | [END_REF][]46) is plotted in yellow in Fig. 7(b) and it overlaps well with the black curve obtained using the exact model. In case of the IM, the numerical results also show that the dispersion curves are divided in two families: with nodes and node-less. The regions of the dispersion diagram corresponding to these two families are separated by the curve described by the equation C 0 = 0 (black and yellow curves for the single-interface problem). In the frame of the IM, we could not prove this property analytically because the field plots in the IM are calculated numerically.

In Fig. 7, in the region of high effective indices, the dispersion curves of the AS1 mode overlap with the curves obtained using single-interface approximations. This confirms our hypothesis that highly asymmetric modes AS1 can be approximated by solutions obtained using the corresponding single-interface models.

C. Comparison with linear states

In Sec. II A, while discussing the field profiles of the modes belonging to the family with nodes, we noticed that they resemble higher-order modes of linear slot waveguides with parameters similar to the NSW studied here. In this section, we will explain the origin of the similarities between these nonlinear and linear modes.

In Fig. 8, we present the nonlinear dispersion diagram obtained using the IM for our NSW. In this plot, the effective index of the mode β is presented as a function of the averaged nonlinear index modification in the waveguide core ∆n :

∆n = 1 d d 0 ∆ndx = 1 d d 0 n (2) 2 Idx, (2)
where the nonlinear parameter n

(2) 2 = α 2 / 0 c l,2 In addition to this plot we present also a dispersion relation (black curves in Fig. 8) of a linear slot waveguide with a homogeneous and linear core and the following parameters: 1 = 3 = -90, n = n 0 + ∆n lin = 3.46 + ∆n lin and d = 400 nm. The parameters 1 , 3 and n 0 = √ l,2 are identical to these in the case of the nonlinear waveguide studied here.

We notice that this linear dispersion diagram is similar to the dispersion plot of the NSW. For the core with index n = n 0 only two modes are present and they are the linear counterparts of the modes S0 and AN0. With the increase of the core index n, the effective index of these modes increases and they become closer to each other. At ∆n lin ≈ 0.1, a higher-order linear mode appears that is a counterpart of the S1 mode. For ∆n lin ≈ 2 and ∆n lin ≈ 3.5, another two higher-order modes appear. They are the linear counterparts of the AN1 and S2 modes, respectively. The effective index of these modes increases rapidly with the increase of ∆n lin . The only modes not present in the linear dispersion curves are the asymmetric modes AS1, AS2, . . . and the symmetric node-less modes SI, SII, etc. The asymmetric modes can not be observed in the linear case because nothing breaks the symmetry in the symmetric linear slot waveguide. The node-less symmetric modes are not supported by the linear slot waveguide because they have purely nonlinear solitonic character [see Figs. 4 and5(e)].

The dispersion curves of the nonlinear modes AN0 and S0 overlap with the corresponding linear dispersion curves only for small ∆n values. The nonlinear modes increase their effective indices β faster than the linear modes. In case of higher-order modes S1, AN1 and S2, the dispersion curves of the linear modes lay below the corresponding nonlinear modes. There is only one common point per mode for these curves (indicated by an open circle in Fig. 8) and it turns out that at this point, the index distribution induced by the nonlinear mode in the nonlinear core is flat (data not shown).

Figure 9 present the comparison of the field profiles H y (x) and E z (x) for nonlinear S1 and AN1 modes and their linear counterparts, at the points where the index distribution induced by the nonlinear mode in the nonlinear core is flat. We observe that the nonlinear profiles overlap perfectly with the profiles of the linear modes normalized to the same amplitude as the nonlinear modes.

The results presented here prove that the modes with nodes found in the NSW are close to the modes of the linear slot waveguide with similar opto-geometric parameters. We explain the similarities between these nonlinear and linear modes using the self-coherent definition of nonlinear modes. This definition was introduced by Townes and co-workers in Ref. [22] and was used later in other works (e.g., Ref. [23]). It defines a nonlinear mode as a linear mode of a linear (graded refractive index) waveguide that is induced by the light distribution of this mode. According to this definition, there is no difference between the nonlinear modes of the NSW for which the nonlinear index modification has a flat distribution and the linear modes of the waveguide with higher, uniformly distributed refractive index of the linear core.

D. Permittivity contrast

In Ref. [14], we have studied the influence of the width of the NSW core on the nonlinear dispersion for this this structure. Here we will discuss the influence of the permittivity contrast between the dielectric core and the metal cladding on the nonlinear dispersion diagrams of symmetric NSW.

First, we will discuss the influence of the metal cladding permittivity on the nonlinear dispersion diagrams of NSW. We have studied the dispersion plots for the NSWs with identical parameters as these used in Sec. II A but with different values of the metal cladding permittivity. We observe that the cladding with higher permittivity (lower in absolute value) allows us to reduce the ∆n threshold values where the bifurcation of the AS1 mode occurs. For metals with permittivity equal to -40, the bifurcation occurs at ∆n ≈ 0.02, which is 4 times lower than in the case of 1 = 3 = -90. For metals with permittivity equal to -20, the bifurcation occurs at ∆n ≈ 8 • 10 -4 , which corresponds to the reduction of the bifurcation threshold by two orders of magnitude with respect to the configuration with 1 = 3 = -90. The dependency of the AS1 mode bifurcation threshold ∆n th on the metal cladding permittivity is illustrated in Fig. 10(a). Looking at this plot, we conclude that with the increase of the metal cladding permittivity (decrease of its absolute value) the bifurcation threshold of the AS1 mode decreases. This decrease is slow in the range of high index contrast between the metal and the nonlinear dielectric permittivity. Although, for

| 1 | = | 3 | close to l,2
the decrease of the bifurcation threshold is more rapid. Changing the metal permittivity from -20 to -15 allows us to decrease the bifurcation threshold by almost two order of magnitudes. For the metal cladding permittivity 1 = 3 = -15 the bifurcation threshold is at the level of ∆n ≈ 10 -5 . This is four orders of magnitude lower than for the 1 = 3 from the range [-400, -90], for which the bifurcation occurs at ∆n ≈ 0.1.

Next, we will study the influence of the change of the core permittivity on the dispersion diagram of the symmetric NSW. We analyzed the plots of the dispersion curves for the NSWs with different linear parts of the core permittivity l,2 . All the other parameters are identical to these used in Sec. II A. The behavior of the bifurcation threshold expressed as the averaged nonlinear index modification ∆n is presented in Fig. 10(b). The increase of the linear part of the core permittivity l,2 is accompanied by a monotonous decrease of the bifurcation threshold. From Fig. 10(b) we notice that the increase of l,2 from 1 to 25 results in the decrease of the bifurcation threshold by approximately three orders of magnitude.

It is interesting to remind that, in the case of changing the permittivity contrast by varying the metal cladding permittivity [see Fig. 10(a)], we observed a decrease of the bifurcation threshold for the AS1 mode with the decrease of the permittivity contrast between the cladding and the core permittivity. On the contrary, decreasing the permittivity contrast by changing the core permittivty, leads to the increase of the bifurcation threshold [see Fig. 10(b)].

This phenomenon can be explained using the field profiles of the symmetric mode for different values of core and metal permittivities. We observe that increasing the permittivity of the core or increasing the permittivity of the metal (decreasing its absolute value) leads to symmetric modes that are more localized on the waveguide interfaces and look closer to separate plasmon on both metal/dielectric interfaces. Because the overlap and therefore the interaction between the two plasmons is weaker, it is easier to break the symmetry of the mode. This explains the decrease of the bifurcation threshold. We conclude that changing the permittivity contrast by varying the linear part of the nonlinear core permittivity, has opposite effect than changing the permittivity contrast by varying the metal cladding permittivity.

III. RESULTS FOR ASYMMETRIC STRUCTURES

In Sec. II, we have comprehensively discussed dispersion diagrams and mode profiles in symmetric NSW structures. In this section, we will discuss the influence of the NSW asymmetry on the dispersion curves. The asymmetry is introduced by sandwiching the nonlinear core by metals with different values of the permittivity on both sides. Asymmetric NSW structures have not been studied before in literature. Here we present the analysis of these structures for the first time.

A. Dispersion relations

Figure 11 presents the nonlinear dispersion diagram obtained using the IM for the structure with the following parameters: core permittivity l,2 = 3.46 2 , the secondorder nonlinear refractive index n In the asymmetric structure only asymmetric modes are present. However, in the dispersion diagram shown in Fig. 11, we divide the modes in two groups: modes that resemble the antisymmetric modes of the symmetric structure for which sgn[E x,0 ] = -sgn[E x,d ] (red curves labeled AN-like) and modes that resemble the symmetric or asymmetric modes of the symmetric structure for which sgn[E x,0 ] = sgn[E x,d ] (blue curves labeled S-like) {see Eq. ([equation] [28][]28) in Ref. [18] for the notations of the electric field components}.

We compare the nonlinear dispersion curves for the asymmetric structure presented in Fig. 11 with the dispersion curves obtained for the symmetric structure shown in Fig. 1. We notice that the dispersion curves for the symmetric and antisymmetric modes from the family with nodes did not change much. The number of modes and the character of their dispersion curves is conserved. The main difference between the dispersion curves of the asymmetric and symmetric structures can be observed for the symmetric and asymmetric modes of the nodeless family. The asymmetry of the structure lifts the double degeneracy of the asymmetric branch AS1 (see green curve in Fig. 1). This branch splits into two branches (see Fig. 11). One of them (the branch with lower effective indices β) is a continuation of the symmetric-like fundamental mode (blue curve) that starts for small power density P c levels. The second branch lays along the first one but has slightly higher power levels (branch with higher β values). The degeneracy of the higher-order asymmetric modes is also lifted by the asymmetry of the structure. These branches also split into two separate branches, similar to the case of the AS1 mode. It is difficult to observe this effect in Fig. 11, where the power density in the core is used as abscissa (even enlarging the region of interest), because the two dispersion curves into which the dispersion curve of the higher-order asymmetric mode splits lay very close to each other. The degeneracy lift of the AS2 mode can be however observed from the dispersion curve β(E 0 ) presented in Fig. 12, where the effective index is shown as a function of the field intensity at the left core interface. In these coordinates the separation of the SI and AS2 curves reflects the degeneracy lift of the AS2 mode.

B. Permittivity contrast study

To finish our discussion of the asymmetric NSW properties, we directly compare the dispersion diagrams β(P c ) of the symmetric structure with these of the asymmetric structures. In Fig. 13, the dispersion plot of the symmetric structure ( 1 = 3 = -90, see Fig. 1) is compared with the dispersion plot for the asymmetric structure ( 1 = -110, 3 = -90, see Fig. 11). Only a vicinity of the bifurcation point of the AS1 mode is presented. We observe that, for low P c values, the dispersion curves of the two low-power modes are slightly modified due to the waveguide asymmetry. For higher values of P c , the dispersion curve of the fundamental mode (upper blue curve) exactly overlaps with the dispersion curve of the asymmetric mode of the symmetric structure. This is a consequence of the fact that the field profiles corresponding to this upper blue curve are strongly localized on the interface with the metal with higher value of the permittivity. These profiles resemble the profiles of the highly asymmetric modes of the symmetric structure [see Fig. 3(a)]. Therefore, we are not surprised that these two dispersion curves overlap. The second curve that results from the degeneracy lift of the asymmetric mode lays above (in terms of P c ) the dispersion curve of the asymmetric mode AS1 (green curve). In Fig. 14, we present a comparison of the dispersion curves of the symmetric structure ( 1 = 3 = -90, see Fig. 1) and the asymmetric structures, where one of the metal permittivity values is higher than in the case of the symmetric structure. The dispersion curves of the symmetric structure (green curves) are compared with these of the asymmetric structures with 1 = -70, 3 = -90 (blue curves), and 1 = -50, 3 = -90 (red curves).

In the case illustrated in Fig. 14, contrary to the one presented in Fig. 13, it is the higher (in terms of P c ) of the two curves that result from the lift of the degeneracy that overlap with the dispersion curve of the asymmetric modes of the symmetric structure. This higher curve corresponds to the modes that are localized on the interface between the core and the metal with permittivity equal to -90. For the structures studied in Fig. 14,= -90 is the lowest cladding permittivity. For that reason, the dispersion curves corresponding to the mode localized on the interface with metal with lower permittivity, overlap with the dispersion curves of the symmetric structure.

Another effect that can be observed in Fig. 14, is that with the increase of the structure asymmetry | 1 -3 | the separation of the two curves that appear as a result of the degeneracy lift, increases, as expected. In the limiting case 1 → 3 , these two curves merge into one doubly degenerate curve.

IV. STABILITY OF THE MAIN SOLUTIONS FOR SYMMETRIC WAVEGUIDES

In the previous sections, we have studied the stationary properties of plasmon-soliton waves using two different modal approaches. From both theoretical and practical points of view, the issue of the stability of these waves arises. In several works, the general problem of the stability of nonlinear waves was studied [START_REF] Kivshar | Spatial Solitons[END_REF][START_REF] Weilnau | [END_REF][START_REF] Kivshar | Optical Solitons, From Fibres to Photonic Crystals[END_REF]. Despite an enormous interest in the properties of nonlinear waves over the last decades, there in no universal condition on their stability [START_REF] Akhmediev | Solitons, Nonlinear Pulses and Beams[END_REF]23]. In most of the cases, the stability must be studied numerically for each of the cases separately. Stability of nonlinear guided waves in fully dielectric structures was studied numerically in Refs. [8,10,[START_REF] Moloney | [END_REF][28][START_REF] Mihalache | Progress in Optics[END_REF][START_REF] Mihalache | [END_REF][31][32][33].

In structures made of metals and nonlinear dielectrics, due to the presence of media with negative permittivity, the problem of stability of plasmon-solitons is difficult to study even numerically. Only in Refs. [34,35] the stability of plasmon-solitons was analyzed for the single nonlinear dielectric/metal interface case, using numerical algorithms (like finite-difference time-domain -FDTD [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF][START_REF]Advances in FDTD Computational Electrodynamics[END_REF]). The propagation of light in plasmonic couplers was studied using Fourier methods based on mode decomposition in linear [START_REF] Čtyroký | [END_REF] and nonlinear [39] regimes. In this section, we study the stability of the plasmon-soliton waves in symmetric NSWs using two methodologies: (i) the topological criterion for fundamental modes of nonlinear waveguides derived in Ref. [23] and (ii) two numerical full-vector methods (using COM-SOL [START_REF]COMSOL[END_REF] and nonlinear FDTD implemented in meep [START_REF]Meep -Abinitio[END_REF][START_REF] Oskooi | [END_REF]).

A. Theoretical arguments

We use here the topological criterion presented in Ref. [23] that is based on the linear stability analysis [43] and the Vakhitov-Kolokolov criterion [44]. The stability criterion presented in Ref. [23], is based only on the topology of the nonlinear dispersion curves and the stability of the modes can be read by analyzing β(I tot ) diagrams in which I tot ≡ +∞ -∞ I(x)dx where I(x) is the intensity density. The validity of this approach was confirmed in multiple settings dealing with purely dielectric structures [10,[START_REF] Mihalache | [END_REF][31][32][33] First, we will recall the principle used to estimate the stability of nonlinear modes using the criterion from Ref. [23]. Then we will use it to analyze the stability of some of the plasmon-solitons found in NSWs.

The stability criterion derived in Ref. [23] uses several assumptions. It provides stability for the fundamental nonlinear modes in structures composed of arbitrary nonlinear material distributed nonuniformly in the transverse direction. The derivation of the stability criterion from Ref. [23] is obtained in the weak guiding approximation for which the electric field satisfies the scalar wave equation. In our study of the TM polarized waves, we consider the case in which it is the magnetic field component that satisfies the scalar wave equation [Eq. ([equation] [5][]5) in Ref. [18]]. We are fully aware of the fact that the metal/nonlinear dielectric structures studied here, in which plasmon-soliton waves propagate, do not fulfill the weak guiding approximation due to high permittivity contrast between the metal and the nonlinear dielectric. This means that interesting nonlinear effects will occur for quite high nonlinear permittivity modifications. In spite of this fact, we use here the criterion from Ref. [23], because the dispersion diagrams obtained for our structures have similar character to the dispersion plots of the fully dielectric structures where the criterion is applicable and because, as it will be shown below, two different numerical propagation simulations of the full vector nonlinear problem confirm at least partially the theoretical predictions.

In Fig. 15, the rules that will be required here to determine the stability of the modes derived in Ref. [23] are schematically shown. Consider the dispersion relation presented in Fig. 16. It shows a zoom of a dispersion diagram, using I tot as variable, for a region that contains the dispersion curves of the main modes for the same structure as the one presented in Fig. 1. The stability of modes changes only at the bifurcation points [23]. To determine the stability, first we have to identify all the bifurcation points on the dispersion diagram β(I tot ). In Fig. 16, Modes appear from or disappear at the points of bifurcation. The next step is to label the sections between the bifurcation points with numbers. The numbers are assigned in the following way. At first, we arbitrarily choose one section and label it with any number (in Fig. 16 we labeled the low intensity section of the symmetric dispersion curve with a number 0). The numbers of all the other sections of dispersion curves are assigned using the geometric rules given in Fig. 15. Finally, after having numbered all the sections of the dispersion curves, we can read the stability of the modes directly from the β(I tot ) dispersion curves. The topological stability criterion presented in Ref. [23] tells us that only the modes corresponding to the parts of the curves with the largest number are possibly stable. In Fig. 16, only the modes labeled by 0 are possibly stable (ps). All the other modes are unstable (u). The stability of all the possibly stable modes can be specified at once, as soon as the stability of one of them is determined. The stability can be determined either using numerical methods or theoretical arguments. The low-intensity section of the symmetric branch in the linear limit corresponds to a linear plasmon in metal/insulator/metal (MIM) configurations, which is stable. Therefore, the solutions corresponding to this section of the nonlinear dispersion curves should be stable. This hypothesis will be confirmed in Sec. IV B.

The high intensity section of the symmetric branch (above the Hopf bifurcation) corresponds to unstable solutions. On the contrary, the section of the asymmetric branch just above the Hopf bifurcation should correspond to stable solutions, because the stability properties of the sections with the same number are the same [23]. On the asymmetric branch (at β ≈ 10) another bifurcation occurs (fold type bifurcation indicated by an open circle). The high effective index section of the asymmetric branch (above the fold bifurcation point) is unstable.

B. Numerical simulations of nonlinear propagations

In the previous section we provide some results about the stability of the plasmon-solitons of the lowest orders using the topological criterion derived in Ref. [23]. In the NSWs studied here the weak guiding approximation, used in the derivation of this topological criterion, is not fulfilled. This fact makes the conclusions drawn using the criterion not definitive. For this reason we also investigate the stability by full-vector numerical simulations.

First, we have used the capabilities of the FDTD method [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF][START_REF]Advances in FDTD Computational Electrodynamics[END_REF] implemented in the meep software [START_REF]Meep -Abinitio[END_REF][START_REF] Oskooi | [END_REF][START_REF] Rodriguez | [END_REF]. The metal permittivity is described by a Drude model to obtain the fixed negative value used at the studied wavelength in the semi-analytical models we developped. The useful computational domain is surrounded at its four edges by absorber regions that prevent back-reflected fields more efficiently than the perfectly matched layers that have also been tested during our FDTD simulations. An example of the FDTD propagation of the asymmetric plasmon-soliton is presented in Fig. 17 This result provides a confirmation of the stability of the first asymmetric mode in the nonlinear slot waveguides when the intensity is above a critical threshold. It can be noticed in Fig. 17(a) that the plasmon-soliton profile is not fully stationary. This is due to the fact that the used current excitations in our FDTD simulations do not generate perfectly the field profile of the FIG. 17. Evolution of the Hy field profile of an asymmetric plasmon-soliton for the stable case (a) with ∆n = 0.0138 and for the unstable case (b) with ∆n = 0.005. These simulations are realized using the FDTD method implemented in the meep software. The parameters are given in Fig. 19. .

asymmetric plasmon-soliton. The used asymmetric excitations contain components that weakly excite the antisymmetric plasmon-soliton (which is studied later in this section). When the excitations match perfectly with the asymmetric mode profile the observed non-stationary behaviour disappears as shown later when the simulation results from the other used numerical method are described. The main symmetric plasmon-soliton is easier to excite in a simple way due to its symmetry property as it can be seen in Fig. 18 where a stable and stationary propagation is shown.

In order to obtain a more general view on the stability of the main modes of the NSW in the frame of the FDTD method, we systematically studied the propagation properties of the three main modes as a function of the spatially averaged refractive index variation ∆n [see Eq. ( 2)]: the first symmetric, asymmetric, and antisymmetric modes. Typically, three cases occur in the simulation results:

• Case 1: The mode is visible during the entire simulation duration. This case is, for example, the one encountered for the asymmetric mode above a given threshold ∆n as it can be seen for example in Fig. 17(a) or in Fig. 18.

• Case 2: The studied mode is generated at the beginning of the temporal evolution but after some times it does not propagate anymore in a selfsimilar way. This case is the one encountered for the main asymmetric mode below a given threshold ∆n as shown for example in Fig. 17(b) where only the most stable part of the propagation is shown.

• Case 3: The investigated mode is not generated by the chosen current source (symmetric, antisymmetric or asymmetric) used to excite it, even at the beginning of the temporal evolution and in the surrounding of the source. This behaviour is observed for the asymmetric mode below the critical power or critical ∆n associated to the Hopf bifurcation.

It is one of the main advantages of the FDTD method to be able to simulate temporal evolution even in the case of unstable modes unlike the other method used later in this section.

As it is shown in Fig. 19 obtained from the FDTD simulations, we are able to build a dispersion diagram FIG. 18. Evolution of the Hy field profile of a symmetric plasmon-soliton for the stable case with ∆n = 0.0018. These simulations are realized using the FDTD method implemented in the meep software. The parameters are given in Fig. 19.

for the first symmetric and asymmetric modes taking into account their stability properties. The given β values for unstable modes, corresponding to the case 2 in the above paragraph, are the ones extracted from the simulations results in the stable initial part of the evolution. It is evident that for the case 3 above, no dispersion data are obtained.

The stability properties of the asymmetric mode from the FDTD simulations differ from the ones deduced from the topological criterion given in the previous section for the stationary case. The asymmetric mode is not stable just above the bifurcation (see case 2 above) for some range of ∆n (see the thin green curve in Fig. 19), and then it becomes stable when ∆n increases (see the thick green curve in Fig. 19). The instability of the asymmetric mode just after the bifurcation has already been described in the field of the spatial soliton studies [47,48]. In our case, the instability can be observed in a relatively extended range of intensity or equivalently of ∆n . This extension of the instability could be due to the way the asymmetric mode is excited in our FDTD simulations and/or to the fact that the metal permittivity is dispersive due to the used Drude model.

It is worth noting that the FDTD dispersion curve for the asymmetric mode differs at high ∆n from the one computed using the stationary IM: here the β values are smaller and the FDTD curve stays concave while the stationary one is convex. Similar saturation effects in nonlinear full-vector temporal simulations have already been described e.g., in Ref. [49]. From the FDTD implementation we use, we can not conclude about the stability property at higher intensities than the ones shown in Fig. 19 due to the limitations of the nonlinear treatment used (see reference [START_REF] Oskooi | [END_REF]). Consequently, we can not check the stability properties around or above the fold bifurcation point described in Section IV A.

As it was expected from the previous section, the first symmetric mode is stable at low ∆n or equivalently at low intensities (see the thick blue curve in Fig. 19).

Its stability is lost slightly before the birth and the partial propagation of the asymmetric mode (see the thin blue curve in Fig. 19). For all ∆n values tested above this transition region, the first symmetric mode is unstable. It is worth mentionning that the stability of this symmetric mode is recovered numerically as soon as the symmetry is forced in the FDTD simulations prohibiting the appearance of asymmetric behaviour.

The topological criterion given in given in Ref. [23] can not be applied to the first antisymmetric plasmonsoliton mode because it is valid only for fundamental modes. Therefore, the stability of this mode can only be infered from numerical simulations. The first antisymmetric mode starts, in the low-intensity regime, from the stable linear antisymmetric plasmon and there is no bifurcations on its dispersion curve. Therefore, we expect this mode to be stable. An example of this stable propagation is shown in Fig. 20. We observe no change of the field profiles during the propagation (the antisymmetric excitation used does not contain any symmetric component).

In Fig. 21, the dispersion curve of the antisymmetric plasmon-soliton is given together with its stability property. The antisymmetric mode is stable up to the maximum intensity that can be treated with the FDTD implementation we use.

The stability properties of the three main plasmonsoliton modes in nonlinear slot waveguides are also verified using the nonlinear propagation scheme implemented in the lastest version of the RF module of COMSOL Multiphysics [START_REF]COMSOL[END_REF]. This approach was successfully used to study the stability of solitons in lattices built of metals and nonlinear dielectrics [50][51][52]. This method is limited to the cases where the studied mode is stable since the iterative numerical method used to compute the fields do not converge in other cases.

According to the conclusions drawn from Fig. 16 in Sec. IV A and from the FDTD simulations, the lowpower section of the symmetric branch corresponds to FIG. 20. Evolution of the Hy field profile of a stable antisymmetric plasmon-soliton with ∆n = 0.0225. These simulations are realized using the FDTD method implemented in the meep software. The parameters are the same as in Fig. 19 stable solutions. This result is also confirmed by the simulation presented in Fig. 22 obtained for a NSW with d = 400 nm. No stable symmetric solution was found above the ∆n transition region using the method implemented in COMSOL confirming the FDTD results already obtained. The stability of the asymmetric branch above the bifurcation region observed in Fig. 19 is confirmed by these numerical simulations as shown in Fig. 23. Figure 24 presents the evolution of the E x electric field component for the asymmetric solutions in such a case.

Figure 25 shows the transverse profiles of the symmetric and asymmetric plasmons-solitons in the NSW. For each symmetry type, we compare the profiles obtained using the interface model (these profiles are used as input in the COMSOL based propagation simulations) with the cuts of the profiles presented in Figs. 22 and 23(b). These comparisons validate the accuracy of the evolution simu-lations and consequently the results obtained concerning the stability properties of the symmetric and asymmetric modes in the NSW. One can notice that the stationary behaviour is more clearly seen in the COMSOL based simulations than in the FDTD ones. This is due to the fact that in the former case we directly use as input the profiles provided by the interface model while in the latter case we use excitation current sources to generate the fields that mimick the stationary field profiles. Since we are investigating nonlinear phenomena, it is not possible to use in the FDTD simulations a part of a linear waveguide to filter the needed profile as it is usually done in FDTD based linear studies [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF].

V. CONCLUSIONS

We have provided detailed results for the plasmonsoliton waves in planar slot waveguides with a finitethicknees nonlinear dielectric core. In symmetric structures, using the semi-analytical models we developped for stationary states, we haved investigated the properties of the first main modes and reported new higher order modes including asymmetric ones that exist at high intensities only. We have also described complete dispersion diagrams for these different modes as a function of various quantities including the total power, the field value at one interface between the metal and the nonlinear core, and also the spatial average of the nonlinear refractive index change. We have proven that the total intensity or equivalently the spatially averaged nonlinear refractive index change corresponding to the Hopf bifurcation threshold from the first symmetric mode to the first asymmetric mode can be reduced by several orders of magnitude with an increase of the permittivity of the core or of the metal cladding. We have also proven the versatility of our semi-analytical models studying asymmetric structures. For such structures, we have described the impact of the metal permittivity contrast that lifts the degeneracy of the doubly degenerated asymmetric mode providing a more complex dispersion diagram than the one of a symmetric structure.

Concerning, the stability of the main symmetric and asymmetric modes, we have used an already derived topological criterion established only in the weak guidance approximation being fully aware that our structures lay beyond its validity range. Nonetheless, as shown by full-vector simulations, the topological criterion predicts correctly the principal stability properties of the main modes of the studied planar nonlinear slot waveguides. Using this criterion, we have shown that the asymmetric mode emerging through a Hopf bifurcation at a critical intensity is stable between this bifurcation and a fold bifurcation located at higher intensity level. The stability of this asymmetric mode is lost at this fold bifurcation. On the contrary, the symmetric mode is unstable for all intensity levels above the Hopf bifurcation while it is stable below.

Using two different full-vector numerical propagation methods, we have studied the stability of the three main modes: the symmetric, asymmetric, and antisymmetric modes. We have shown that the asymmetric mode is stable above a critical intensity slightly larger than the threshold associated with the Hopf bifurcation computed for the stationary states from our semi-analytical models at least up to the maximum level of tested intensities. The symmetric mode is shown to be unstable slightly below and slightly above the Hopf bifurcation threshold, and to be stable at lower intensities. For all tested intensities, these results confirm qualitatively the results derived from the topological criterion even if quantitative differences exist. Finally, we have also proven numerically that the anti-symmetric mode is stable in the entire range of tested intensities.

These stability results together with those about the decrease of the bifurcation threshold should facilitate the design of specific structures in order to make possible the experimental observation of these plasmon-soliton waves more than thirty years after their theoretical discovery. Future studies should be dedicated to the further reduction of the bifurcation threshold and to the study of more sophisticated configurations.

FIG. 1 .

 1 FIG. 1. Dispersion diagrams β(Pc) for the symmetric NSW obtained using (a) the JEM and (b) the IM.

FIG. 6 .

 6 FIG.6. Dispersion diagram β(E0) for the symmetric NSW obtained using the IM. The vertical black lines indicate to the values of the total electric field amplitude E0 corresponding to the field profiles depicted in Figs.3 and 4. Open circles correspond to the field profiles shown in Fig.5.

FIG. 7 .

 7 FIG. 7. Dispersion diagrams (a) β(H0) obtained using JEM and (b) β(E0) obtained using the IM for the symmetric NSW. Dispersion curves presenting single-interface approximation obtained using models derived in Ref. [21] are shown in black. Additionally, the curves corresponding to the analytical expression for the single-interface dispersion relation {Eqs. ([equation][42][]42) and ([equation][46][]46) in Ref. [18]} are shown in yellow.

  FIG.8. A comparison of the nonlinear (red, blue and green curves) and the linear dispersion plots (black curves) of the symmetric slot waveguides. In case of the linear waveguide ∆n is equivalent to ∆n lin . Points correspond to the modes presented in Fig.9.

  FIG.9. A comparison of (a), (c) Hy(x) and (b), (d) Ez(x) for the nonlinear modes S1 (blue curve) and AN1 (red curve) and the normalized profiles of their linear counterparts (black curves) at the common points of the black (for linear slot waveguide) and blue or red (for NSW) dispersion curves indicated by open circles in Fig.8.

FIG. 10 .

 10 FIG. 10. Average nonlinear index change at the appearance of the asymmetric AS1 modes ∆n th as a function of the absolute value of (a) the metal cladding permittivity of the symmetric waveguide | 1| = | 3| and (b) the linear part of the nonlinear core permittivity l,2 All the other parameters of the NSW are identical to these used in Sec. II A.

( 2 ) 2 = 2 •

 222 FIG. 11. Dispersion diagram obtained using the IM for the asymmetric structure with 1 = -110 and 3 = -90 (for the scheme of the structure see Fig. 2). Blue curves correspond to the modes for which sgn[Ex,0] = sgn[E x,d ] and red curves correspond to the modes for which sgn[Ex,0] = -sgn[E x,d ] {see Eq. ([equation][28][]28) in Ref. [18] for the notations}. Compare this dispersion diagram for the asymmetric structure with the dispersion diagram for the symmetric structure presented in Fig. 1(b).

FIG. 12 .

 12 FIG. 12. Dispersion curves β(E0) for the asymmetric structure with 1 = -110 and 3 = -90. Compare this dispersion diagram for the asymmetric structure with the dispersion diagram for the symmetric structure presented in Fig. 6.

FIG. 13 .

 13 FIG. 13. Dispersion curves of the asymmetric NSW with 1 = -110, 3 = -90 (blue curves) and the symmetric structure 1 = 3 = -90 (green curves).

FIG. 15 .

 15 FIG. 15. Rules of assigning the stability of the modes for two specific cases extracted from Fig. 2 in Ref. [23]: (a) the fold bifurcation (open circle) and (b) the Hopf bifurcation (open square). Thick lines indicate a doubly degenerate branch, whereas thin lines indicate non-degenerate dispersion curves.

10 FIG. 16 .

 1016 FIG. 16. Zoom on the region of the dispersion diagram with the birth of the first order asymmetric mode. Bifurcation points are marked with an open circle for fold bifurcation and an open square for Hopf bifurcation. The numbers facilitating the stability analysis are assigned to the sections of the dispersion curves according to the rules presented in Fig. 15. Labels 'ps' and 'u' denote possibly stable and unstable modes, respectively.

  (a) through the evolution of the electric field component H y , where the sinusoidal phase modifications are visible.

FIG. 19 . 2 = 2 •

 1922 FIG. 19. Dispersion and stability results, obtained from FDTD simulations, for the first symmetric and asymmetric modes of the symmetric NSW. Thick curves denote a stable propagation while thin ones denote an unstable propagation. The parameters are the following: core permittivity l,2 = 3.46 2 , the second-order nonlinear refractive index n (2) 2 = 2 • 10 -17 m 2 /W, core thickness d = 500 nm, metal permittivities 1 = 3 = -6 at a free-space wavelength λ = 1.0 µm.

  FIG.[START_REF] Walasik | Plasmon-soliton waves in metal-nonlinear dielectric planar structures[END_REF]. Evolution of the Hy field profile of a stable antisymmetric plasmon-soliton with ∆n = 0.0225. These simulations are realized using the FDTD method implemented in the meep software. The parameters are the same as in Fig.19.

FIG. 21 . 2 = 2 •

 2122 FIG. 21. Dispersion and stability results for the first antisymmetric modes of the symmetric NSW obtained from FDTD simulations. The thick curve denotes a stable propagation. The parameters are the same as in Fig. 19.

FIG. 23 .FIG. 24 .FIG. 25 .

 232425 FIG. 23. Evolution of the electric field norm during the propagation of asymmetric modes located between the Hopf bifurcation and the fold bifurcation. The average nonlinear index change in the core ∆n induced by these modes is equal to (a) 2 • 10 -3 , (b) 3 • 10 -3 , and (c) 4 • 10 -3 . The shown propagation distance is approximately 12 free-space wavelengths.The parameters are the same as in Fig.22. These simulations are realized using the COMSOL software.

ACKNOWLEDGMENTS

This work was supported by the European Commission through the Erasmus Mundus Joint Doctorate Programme Europhotonics (Grant No. 159224-1-2009-1-FR-ERA MUNDUS-EMJD). This work was not funded by the French Agence Nationale de la Recherche. W. W and G. R. would like to thank Yaroslav Kartashov for helpful comments on this work and Alejandro Rodriguez for useful discussions about the meep software. G. R. would like to thank the Scilab consortium.