Laurent Guillier
email: laurent.guillier@anses.fr

Philippe Fravalo
email: philippe.fravalo@lecnam.net

Alexandre Leclercq

Anne Th Ébaut

Pauline Kooh

Vasco Cadavez
email: vcadavez@ipb.pt

Ursula Gonzales-Barron

Anne Thébaut

Risk factors for sporadic Yersinia enterocolitica infections: a systematic review and meta-analysis

Keywords: Research synthesis, Case-control studies, Meta-regression, Yersinia enterocolitica

Yersinia enterocolotica is an important causative agent of diarrheal illness. A systematic review and meta-analysis of case-control studies were performed to determine the main risk factors associated with sporadic Y. enterocolitica infections. Suitable scientific articles published between 1987 and 2017 were identified through a systematic literature search and subject to methodological quality assessment. From each study, odds ratios (OR) were extracted or calculated, as well as study characteristics such as population type, design, type of model used and risk factor categorization. Mixed-effects meta-analytical models were adjusted by population type to appropriate data partitions. From 807 identified references, the quality assessment stage was passed by 14 case-control studies focusing on sporadic Y. enterocolitica infections which provided 165 ORs for meta-analysis. All studies considered Y. enterocolitica as the cause of sporadic infections and are mainly located in Europe.

The meta-analysis identified host-specifics factors, animal and food exposures as significant risk factors. The meta-analysis confirms the predominant role of the pig reservoir. The occupational contact with pigs and the consumption of pork meat are significantly associated with sporadic Y. enterocolitica infections occurrence. The consumption of raw or undercooked pork meat is also a very important risk factor. Untreated drinking water was also identified as risk factor. Further studies with other enteropathogenic Yersinia species, especially Y. pseudotuberculosis, and/or from other continents would help to refine conclusions of the meta-analysis of the risk factors of yersiniosis.

Highlights

-Meta-analysis of 14 case-control studies about sporadic human yersiniosis -Available studies focus on Y. enterocolitica and mainly from Europe -Consumption of pork meat and untreated drinking water were the main determinant of sporadic Y. enterocolitica infections -The consumption of undercooked or raw pork meat considerably increases risk levels

Introduction

The genus Yersinia is composed of 19 species, among which two are enteropathogenic to humans (Y. enterocolitica, and Y. pseudotuberculosis) [START_REF] Savin | Genus-wide Yersinia core-genome multilocus sequence typing for species identification and strain characterization[END_REF]. Enteropathogenic Yersinia are mainly found in temperate or cold regions, such as Central and Northern Europe, New-Zealand and North America. Strains of Yersinia are ubiquitous and occur in soil, surface water, food and in the digestive tract of various animal species [START_REF] Le Guern | Yersiniosis in France: overview and potential sources of infection[END_REF].

Y. enterocolitica is the main species in the genus associated with yersiniosis, which can be defined as a mild-moderate self-limiting gastroenteritis [START_REF] Galindo | Pathogenesis of Y. enterocolitica and Y. pseudotuberculosis in human yersiniosis[END_REF]. Y. enterocolitica is subdivided into 6 biotypes (1A/1B, 2, 3, 4, 5) based on biochemical tests and more than 70 serotypes [START_REF] Le Guern | Yersiniosis in France: overview and potential sources of infection[END_REF][START_REF] Wauters | Revised biogrouping scheme of Yersinia enterocolitica[END_REF]. Y. enterocolitica infection is usually characterized by diarrhea, fever and abdominal pain (Savin et al., 2008). Patients do not always present all three symptoms. Complications such as reactive arthritis or, more rarely, sepsis can also be observed in individuals [START_REF] Hoffmann | Acute Complications and Sequelae from Foodborne Infections: Informing Priorities for Cost of Foodborne Illness Estimates[END_REF][START_REF] Rosner | Clinical aspects and self-reported symptoms of sequelae of Yersinia enterocolitica infections in a population-based study, Germany 2009-2010[END_REF]. Symptoms develop 4 to 7 days after exposure and persist for 5 to 14 days, or several weeks in the case of diarrhea [START_REF] Laukkanen-Ninios | Epidemiology, virulence genes, and reservoirs of enteropathogenic Yersinia species, Foodborne and waterborne bacterial pathogens: epidemiology, evolution and molecular biology[END_REF].

Y. enterocolitica infections are the fourth most frequently reported bacteriologically related foodborne zoonosis in Europe [START_REF] Efsa | The European Union One Health 2018 Zoonoses Report[END_REF]. Reported incidence varies according to country, the method used for diagnostic and the performance of the reporting system. It ranges from less than 1 case per 100,000 in the US [START_REF] Tack | Preliminary incidence and trends of infections with pathogens transmitted commonly through food -Foodborne diseases active surveillance network, 10 U.S. Sites, 2015-2018[END_REF] up to 19.2 in New Zealand [START_REF] Strydom | Evaluating sub-typing methods for pathogenic Yersinia enterocolitica to support outbreak investigations in New Zealand[END_REF]. Yet the number of illness could be higher when considering underreporting and underdiagnosis [START_REF] Van Cauteren | Estimated Annual Numbers of Foodborne Pathogen-Associated Illnesses, Hospitalizations, and Deaths, France, 2008-2013[END_REF]. As an example, in France, the corrected annual community incidence rate was estimated at 36 cases per 100,000 for yersiniosis, about 30 times higher than the reported incidence rate [START_REF] Van Cauteren | Estimated Annual Numbers of Foodborne Pathogen-Associated Illnesses, Hospitalizations, and Deaths, France, 2008-2013[END_REF]. Microbiological and epidemiological investigations make it possible to identify outbreaks and to trace back the food at the origin of the human cases [START_REF] Espenhain | Cross-border outbreak of Yersinia enterocolitica O3 associated with imported fresh spinach, Sweden and Denmark[END_REF]. Although outbreaks are identified [START_REF] Konishi | Two outbreaks of Yersinia enterocolitica O: 8 infections in Tokyo and the characterization of isolates. Kansenshogaku zasshi[END_REF][START_REF] Macdonald | Yersinia enterocolitica outbreak associated with ready-to-eat salad mix, Norway[END_REF], a majority of yersiniosis cases caused by Y. enterocolitica are not associated with any known outbreak [START_REF] Marimon | Thirty years of human infections caused by Yersinia enterocolitica in northern Spain: 1985-2014[END_REF] and are classified as sporadic cases. The question of specific environmental or food exposures and their respective weight for these sporadic cases arise.

A large variety of methodological approaches for source attribution of sporadic cases is available [START_REF] Mughini-Gras | Critical orientation in the jungle of currently available methods and types of data for source attribution of foodborne diseases[END_REF]). Yet, very few source attribution studies have been carried out for yersinosis. Expert elicitation was the main method applied so far [START_REF] Batz | Ranking the disease burden of 14 pathogens in food sources in the united states using attribution data from outbreak investigations and expert elicitation[END_REF][START_REF] Zanabria | Source attribution at the food sub-product level for the development of the Canadian Food Inspection Agency risk assessment model[END_REF].

Another way of identifying the sources of sporadic cases is case-control studies. In casecontrol studies, the association of cases with various food exposures is usually measured through odds ratios (ORs) (Pires et al., 2009). Meta-analyses of these studies can provide information on exposure pathways of interest [START_REF] Devleesschauwer | Associating sporadic, foodborne illness caused by Shiga toxin-producing Escherichia coli with specific foods: a systematic review and meta-analysis of case-control studies[END_REF].

The aim of this study is to perform a systematic review of case-control studies for human sporadic Y. enterocolitica infections, and, subsequently, to perform a meta-analysis to synthesize data on factors associated with sporadic infection, combining the odds ratio from a selection of relevant studies [START_REF] Gonzales-Barron | Strategy for systematic review of observational studies and meta-analysis modelling of risk factors for sporadic foodborne diseases[END_REF].

Methods

The systematic review process as well as the meta-analysis model are described in depth in a methodological paper [START_REF] Gonzales-Barron | Strategy for systematic review of observational studies and meta-analysis modelling of risk factors for sporadic foodborne diseases[END_REF].

Systematic review

The keywords were defined attending the review question which is to evaluate the association between a (risk) factor and sporadic Y. enterocolitica infections risk in a population exposed to it. This review question was identified to have a typical PECO structure (Population, Exposure, Comparator and Outcome as key elements) [START_REF] Efsa | Application of systematic review methodology to food and feed safety assessments to support decision making[END_REF]..The literature search was conducted in March 2017 using a combination of keywords related to (1) "Yersinia enterocolitica OR yersiniosis", (2) "case-control OR risk factor OR cohort" (3) "infection OR disease", joined by the logical connector AND. Systematic searches using a combination of suitable keywords were conducted using five bibliographic search engines (ISI Web of Science, PubMed, Scielo, Science Direct and Scopus). The literature search was conducted for English, French, Portuguese and Spanish languages.

The screening criteria followed the PRISMA (Preferred Reporting Items for Systematic Reviews and Meta-analyses) method [START_REF] Moher | Preferred reporting items for systematic reviews and meta-analyses: the PRISMA statement[END_REF]. Each reference record was screened by at least two persons for relevance for inclusion in the meta-analysis study, and subsequently, the methodological quality of the "candidate" studies was assessed using preset quality criteria. The first criterion for inclusion was related to the definition of the disease. The included studies should have used a laboratory-based definition of a case, based on a detection or isolation from cultures of the cases. The second criterion refers to the study design. Only the case-control study design was considered.

Moreover, the methodological quality of each "relevant" primary case-control study was assessed using a checklist comprised of six areas of concern. It concerns the appropriate selection of the controls in order to avoid selection bias; the adjustment to correct for confounders, the criteria for the comparability between cases and controls; the data analysis appropriate to the type of design, matched or unmatched; the responses rates for the exposed and control groups; and the provision of crude data and/or adjusted ORs and either confidence interval or p-value. All these points permitted to assess the overall quality of underlying methods, appropriate statistics, sensible data and quality of reporting/interpretation.

Finally, the quality and completeness of statistical analysis were assessed. Primary studies that passed the screening for relevance were marked as having the potential for bias if they failed to meet at least one of the methodological quality assessment criteria. After careful data extraction, meta-analysis models were adjusted within appropriate risk factor data partitions in order to estimate overall ORs, extracting variability due to primary studies and type of statistical analysis. Diagnostics based on Cook's distance was assessed for every meta-analysis model in order to remove any influential OR originated from studies deemed as having some potential for bias.

Data synthesis

The joint meta-analytical data was first described using basic descriptive statistics. The Meta-analysis models were then fitted to each of the data partitions or subsets in order to estimate the overall OR due not only for food vehicles but also to travel, host-specific factors and transmission pathways related to person-to-person transmission, animal contact and environmental exposures. In accordance with Gonzales-Barron et al. ( 2019), the meta-analytical models were fitted separately by population type, which are children and mixed population. For some food classes, the effects of handling (i.e., eating raw, undercooked)

and setting (i.e., eating out) were assessed by dividing the mean ORs when food was mishandled (or, alternatively, when food was prepared outside the home) by the base ORs.

The statistical analysis was designed to assess the effect of study period and analysis type

Results

Descriptive statistics of the case-control studies

In the systematic review of risk factors pertaining to sporadic Y. enterocolitica infections, a total of 807 bibliographic sources were identified using the defined keywords in the five search engines, from which only 33 passed the full assessment for eligibility comprising case-control/case and cohort studies from both sporadic illnesses and outbreaks (Figure 1).

A total of 19 fully-documented case-control studies investigated the source(s) of outbreaks and were thus not included in the study. The meta-analysis was undertaken using 14 primary studies focused on sporadic disease (Figure 1). These published studies were conducted between 1987 and 2010. Table 1 compiles a list of the case-control studies along with their main features.

The eligible studies jointly provided 165 categorized odds-ratios for meta-analysis. A total of 66 ORs were retrieved from 8 case-control studies performed before the year 2000, while 99

ORs were excerpted from 6 case-control studies undertaken after 2000. The majority of primary studies investigated sporadic Y. enterocolitica infections caused by undetermined serogroups (9 case-control studies) representing 72% of the ORs. The countries whose case-control studies contributed the largest body of results for sporadic Y. enterocolitica infections were Finland (2 studies, 55 ORs), Norway (4 studies, 30 ORs), New Zealand (2 studies, 25 ORs) and Germany (1 study, 23 ORs).

Twelve case-control studies investigated pathways of exposure in adult or mixed population, while only three case-control studies had children as the target population. The only primary study investigating determinants of disease in both children and mixed population was that of [START_REF] Rosner | Risk factors for sporadic Yersinia enterocolitica infections, Germany 2009-2010[END_REF] (Table 1). Whereas 80% of the ORs originated from exposures evaluated in the mixed/adult population, 20% of the ORs were quantified from ill children cases. As a rule, because of their distinct routes of exposure, the ORs for children and mixed populations were not joined in a single meta-analysis model, but in separate meta-analyses.

Data from both populations were merged only when the ORs belonging to the children population were too few to run a separate meta-analysis model.

Meta-analysis

For every data partition, the meta-analyzed risk factors are presented in summary tables only when significant (Tables 2, 3 and 4) and when more than one study inform the risk factor.

Pooled ORs were considered significant when the lower bound of the 95% CI was equal or greater than 1. Appendix-2 provides values for non-significant factors as well as significant risk factors that arise from a single study.

Risk factors associated with non food-related transmission pathways

Among travel, host-related factors, contact with the environment and animals, only a few categories were found as significantly associated with sporadic Y. enterocolitica infection cases (Table 2).

Among the main risk categories that were amenable to be meta-analysed, the host-specific factors, associated with chronic diseases in the mixed population, represented as a whole the most important risk factor for acquiring yersiniosis (overall OR=3.132; 95% CI: 1.191 -8.235). Malnutrition was also associated to yersiniosis in the children population but it arises from a single study (Appendix 2).

Within "contact with animals" categories, the occupational animal contact (overall OR=1.434; 95% CI: 1.164 -1.766) was the only significant risk factor of sporadic Y. enterocolitica infections. It is worthy to mention that, within the mixed population, the significance is mainly explained for occupational exposure data by exposure to pigs (Figure 2). Contacts with pets or occasional contact with farm animals were not associated with sporadic Y. enterocolitica infections. A single study explore the contact with wild animals [START_REF] Rosner | Risk factors for sporadic Yersinia enterocolitica infections, Germany 2009-2010[END_REF] and it does not allow us to draw any conclusions for this risk factor.

For the mixed population and children, sporadic Y. enterocolitica infections occurred more frequently among people exposed to untreated drinking water (overall OR=1.804; 95% CI:

1.230 -2.645) (Figure 3). For children, activities in playground were associated with sporadic Y. enterocolitica infections (OR=1.580; 95% CI: 1.343 -1.857).

In all the meta-analytical data partitions, there was no effect of the study period (before and after 2000) on the measured ORs. The meta-analysis revealed that who travelled recently had no significantly higher (OR=2.947; 95% CI: 0.920 -9.438) probability of getting infected with sporadic Y. enterocolitica infections that people who did not travel.

Risk factors associated with food-related transmission pathways

The studies explore a wide range of foods and food preparation/consumption practices (e.g. hygiene, cooking, and place of consumption). For the mixed population and children, within main large food categories, meat was the only significant risk factor (Table 2). Among meats, pork presents the highest risk factors (Table 3). The risk appears to be higher for children 4).

(
For most of the meta-analytical models reported in Table 2, 3 and 4, the statistical tests indicated the absence of potential significant publication bias is above 5%. Exceptions were observed in partitions related to meat and pork in both mixed and children population (Tables 2 and3), fast food in children (Table 3), and the effect of meat handling (Table 4). For better assessing the publication bias (above 5%), the funnel plot for those models is given in Figure 4. For all of them, they were an asymmetry, and a lack of non-significant studies with smaller ORs, that could lead to OR overestimation. Moreover, the intra class correlation I 2 indicates, in all Tables, a high heterogeneity (<75%). However, residual between-study heterogeneity (p-value often below 0.05 for Q or QE) was observed for the data partitions.

Discussion

In this study, the aim was to synthesize data produced by published case-control studies.

The weight on common factors associated with sporadic Y. enterocolitica infections, among them the relative importance of host-related factors, contact with animals and the environment, as well as food-related factors have been characterized. The results of the case-control surveys confirm the importance of the pig reservoir and of the practices of eating pork meat (raw or undercooked). Indeed pigs are the main carriers of Y. enterocolitica [START_REF] Drummond | Yersinia enterocolitica: a brief review of the issues relating to the zoonotic pathogen, public health challenges, and the pork production chain[END_REF]. Moreover, the main biotype isolated in pork, that is biotype 4 [START_REF] Raymond | Diversity of Yersinia enterocolitica isolated from pigs in a French slaughterhouse over 2 years[END_REF] matches with the most common biotype isolated in sporadic cases.

Despite this association, it should be noticed that consumption of raw, undercooked pork, is rarely associated with outbreaks [START_REF] Grahek-Ogden | Outbreak of Yersinia enterocolitica serogroup O:9 infection and processed pork, Norway[END_REF][START_REF] Self | Outbreaks attributed to pork in the United States, 1998-2015[END_REF]. The situation is the opposite for vegetables and milk products. Both food categories were not identified as risk factors in the meta-analysis (see Appendix 2), although several outbreaks have been associated to these food products, especially for vegetables, in recent years [START_REF] Espenhain | Cross-border outbreak of Yersinia enterocolitica O3 associated with imported fresh spinach, Sweden and Denmark[END_REF][START_REF] Konishi | Two outbreaks of Yersinia enterocolitica O: 8 infections in Tokyo and the characterization of isolates. Kansenshogaku zasshi[END_REF][START_REF] Macdonald | Yersinia enterocolitica outbreak associated with ready-to-eat salad mix, Norway[END_REF].

Untreated drinking water was also identified by the meta-analysis as a risk factor for sporadic Y. enterocolitica infections. This finding is also confirmed by epidemiological analysis of some outbreaks that were related to the contamination of food products by untreated drinking water contaminated by Y. enterocolitica [START_REF] Ackers | An outbreak of Yersinia enterocolitica O:8 infections associated with pasteurized milk[END_REF][START_REF] Tacket | An outbreak of Yersinia enterocolitica infections caused by contaminated tofu (soybean curd)[END_REF].

Playground attendance was found to be a risk factor for children by the meta-analysis. This finding is somewhat in contradiction with the results from literature where Y. enterocolitica strains isolated in the environment are mainly non-pathogenic strains [START_REF] Le Guern | Yersiniosis in France: overview and potential sources of infection[END_REF]. enterocolitica and is regularly implicated in foodborne outbreaks [START_REF] Pärn | Outbreak of Yersinia pseudotuberculosis O:1 infection associated with raw milk consumption, Finland, spring 2014[END_REF][START_REF] Rimhanen-Finne | Yersinia pseudotuberculosis causing a large outbreak associated with carrots in Finland, 2006[END_REF]. Specific case-control studies would help to confirm that risk factors also hold for other enteropathogenic Yersinia species.

A potential difficulty to interpret the results of the meta-analysis is related to the relative importance of biotypes/serotypes in sporadic cases according to the different continents or the different period. Pathogenic Y. enterocolitica strains belong to biotypes 1B, 2, 3, 4 and 5, while biotype 1A is currently considered non-pathogenic, although its association with diarrheal forms is still debated [START_REF] Huovinen | Symptoms and sources of Yersinia enterocolitica-infection: a case-control study[END_REF]. Although bioserovar 4/O:3 is the most commonly identified type world-wide in human cases [START_REF] Hunter | Identification and typing of Yersinia enterocolitica and Yersinia pseudotuberculosis isolated from human clinical specimens in England between 2004 and 2018[END_REF][START_REF] Le Guern | Yersiniosis in France: overview and potential sources of infection[END_REF][START_REF] Strydom | Evaluating sub-typing methods for pathogenic Yersinia enterocolitica to support outbreak investigations in New Zealand[END_REF], some particularities can be identified. For example, biotype 2 has recently emerged as the most common biotype causing yersiniosis in New Zealand [START_REF] Strydom | Evaluating sub-typing methods for pathogenic Yersinia enterocolitica to support outbreak investigations in New Zealand[END_REF]. In the same way, biotype 1B was the most common biotype in the 1980s for US while now biotype 4 is the most commonly isolated biotype [START_REF] Savin | Isolation of a Yersinia enterocolitica biotype 1B strain in France, and evaluation of its genetic relatedness to other European and North American biotype 1B strains[END_REF].

Another limitation of the meta-analysis is related to the origin of the studies. Most of the case-control studies included in the meta-analysis came from Europe, whereas this pathogen is also present in other continents [START_REF] Duan | Prevalence of Yersinia enterocolitica bioserotype 3/O: 3 among children with diarrhea, China, 2010-2015[END_REF][START_REF] Lucero-Estrada | An overview of Yersinia enterocolitica and related species in samples of different origin from San Luis, Argentina[END_REF]. Casecontrol studies from other continents would help to refine our conclusions.

Several sub-typing methods have been developed for outbreak investigation [START_REF] Strydom | Evaluating sub-typing methods for pathogenic Yersinia enterocolitica to support outbreak investigations in New Zealand[END_REF], such as biotyping, serotyping, PFGE, MLST and MLVA (Mughini-Gras et al., 2019).

These methods have been used in the context of outbreak investigations, but they have never been used as input of source attribution models for sporadic cases of yersiniosis.

Frequency-matching models for source attribution could infer the most likely sources of human sporadic cases by comparing their subtype frequencies, weighted by factors like prevalence in these sources and the human exposure to them.

In the same way, the application of WGS methods was until recently limited to outbreak investigations [START_REF] Espenhain | Cross-border outbreak of Yersinia enterocolitica O3 associated with imported fresh spinach, Sweden and Denmark[END_REF][START_REF] Inns | First use of whole-genome sequencing to investigate a cluster of Yersinia enterocolitica, Liverpool, United Kingdom[END_REF]. The recent publication of the Yersinia cgMLST schema [START_REF] Savin | Genus-wide Yersinia core-genome multilocus sequence typing for species identification and strain characterization[END_REF] will help to apply population structure models or machine learning approaches for estimating the attribution of sporadic yersiniosis cases to the potential sources.

Conclusion

The results of this meta-analysis confirm the importance of pork and untreated drinking water for sporadic Y. enterocolitica infections. The study also stresses the role of other routes than food, such as environmental and animal contact. The published case-control studies included this meta-analysis are stored in a database, which will be updated with other relevant studies published in the future.

To refine these results, risk factors should be investigated by Y. enterocolitica biotypes. Such investigations, together with specific studies dedicated to Y. pseudotuberculosis, could allow to identify species or sub-type source-specific risk factors and to infer the underlying transmission pathways.

  number of ORs per country, per year, and per type were calculated. Next, data was partitioned into subsets of categories of risk factors. The used source categorization scheme included travel, host-specific factors and hierarchical pathways of exposurecomprising person-to-person, animal, environment and food routes of transmission. Gonzales-Barron et al. (2019) provides the full list of risk factors, categories and subcategories. Food class and handling are also used for exploring in depth food related factors. The full hierarchy of risk factors is provided in Appendix 1.

(

  univariate/multivariate) on the final result. All meta-analytical models were essentially weighted random-effects linear regression models. A random effect term allow taking into account the study effect ("dependence between studies") nested into the risk factors categories[START_REF] Gonzales-Barron | Strategy for systematic review of observational studies and meta-analysis modelling of risk factors for sporadic foodborne diseases[END_REF]. Once a meta-analytical model was fitted, influential diagnostics statistics were assessed in order to remove any influential observation originating from studies marked as having potential-for-bias. Publication bias was assessed by funnel plots and statistical tests. Heterogeneity between studies was assessed by three indicators, the between-study variability (τ 2 ), the QE test investigating residual heterogeneity, the variance of residuals and the intra-class correlation I 2 (Gonzales-[START_REF] Gonzales-Barron | Strategy for systematic review of observational studies and meta-analysis modelling of risk factors for sporadic foodborne diseases[END_REF]. All analyses were conducted in the R software implemented with the metafor package (Viechtbauer, 2010).

  It's important to note that all the 14 studies focused on Y. enterocolitica. None of them integrated cases of sporadic Y. enterocolitica infection associated to other species, especially Y. pseudotuberculosis. This latter is thought to have close ecology with Y.

Figure 1 .

 1 Figure 1. Flow chart of the literature search for case-control studies of sporadic Y.

Figure 2 .

 2 Figure 2. Forest plot of the associations of yersiniosis with the category of

Figure 3 .

 3 Figure 3. Forest plot of the associations of sporadic Y. enterocolitica infections with

Figure 4 :

 4 Figure 4: Funnel plots of studies investigating categorized risk factors. A) Food in

Table 2 . Results of the meta-analysis on the main risk factors of sporadic Y. enterocolitica infections

 2

	Population	Risk factor	Pooled OR [IC95%]	N/n* p-value of risk factor	Publication bias p-value	Points removed **	Heterogeneity analysis***
				Host specific
								Tau2=0.777
			3.132					Q(df = 4) = 35.551,
	Mixed	Host specific	[1.191 -	4/5	0.021	0.956	0	p-val < .0001
			8.235]					S2=1.241
								I2=38,491
				Animal

Acknowledgments

The authors would like to thank Anses staff and the members of the Anses Working Group on Source Attribution of Foodborne Diseases: Moez Sanaa, Jean Christophe Augustin, Frédéric Carlin, Julie David, Nathalie Jourdan da-Silva, Lapo Mughini-Gras, Nicole Pavio, Isabelle Villena, and Laurence Watier. U. Gonzales-Barron and V. Cadavez are grateful to the Foundation for Science and Technology (FCT, Portugal) for financial support through national funds FCT/MCTES to CIMO (UIDB/00690/2020). U. Gonzales-Barron acknowledges FCT, P.I., for the institutional scientific employment programme. Author Statement Laurent Guillier: Writing -original draft. Philippe Fravalo: Writing -Review & Editing Alexandre Leclercq: Writing -Review & Editing Anne Thébault: Methodology, Formal analysis, Visualization Pauline Kooh: Methodology, Project administration Vasco Cadavez: Methodology, Investigation, Formal analysis. Ursula Gonzales-Barron: Methodology, Investigation, Formal analysis, Writing-reviewing & editing.

Declarations of interest:

The authors declare no conflict of interest.

Appendix 1: List of risk factors used for meta-analysis of food-borne pathogen ).


Gérard Granet
email: gerard.granet@uca.fr

A new spectral element method with modified Legendre polynomials for modal analysis of lamellar gratings

We report on the derivation of a spectral element method the originality of which comes from the use of a hierarchical basis built with modified Legendre polynomials. We restrict our work to the TM polarization which is the most challenging. The validation and the convergence is carefully checked for metallic dielectric gratings. The method is shown to be highly efficient and remains stable for huge truncation numbers. All the necessary information is given so that non-specialists can implement the method.

I. Introduction

During the last decade, plasmonics has taken an increasingly important part in optics. In this context, dielectric metallic structures, whether periodic or not, play a fundamental role because they support plasmons which correspond to electromagnetic fields with a very fast spatial variation and local huge enhancement. These latter basic properties make the design of plasmonic structures very challenging for optical engineers regardless of the numerical method used. For translation invariant structures, modal methods and mode matching techniques are very interesting in that they avoid discretization along the axis of invariance. Numerical modal methods differ in the choice of expansion and test functions. One can distinguish full domain functions or sub-domain functions according to whether they are defined on the whole domain or only on some part of it. For periodic problems, the most popular method is certainly the well-known Fourier Modal Method (FMM) [START_REF] Moharam | Diffraction analysis of dielectric surface-relief gratings[END_REF] in which expansion and test functions are pseudo periodic functions. On the contrary, using sub-domain expansions within each domain that constitute the structure allows to express rigorously the different continuity relations that determine the eigenvalue problem and thus leads to exponential convergence for the eigenvalues and eigenvectors. The polynomial modal method [START_REF] Morf | Exponentially convergent and numerically efficient solution of Maxwell's equations for lamellar gratings[END_REF] is one such method. It was revisited in [START_REF] Edee | Modal method based on sub-sectional Gegenbauer polynomial expansion for lamellar grating[END_REF], [START_REF] Edee | Modal method based on subsectional Gegenbaueur polynomial expansion for lamellar gratings: weighting function, convegence and stability[END_REF] and [START_REF] Randriamihaja | Polynomial modal analysis of lamellar diffraction gratings in conical mounting[END_REF]. It was shown to outperform all other numerical modal methods for dielectric metallic structures. However, its implementation may scare those who are not familiar * Corresponding author with numerical methods. Due to locality of basis functions, the numerical polynomial modal methods look quite similar to finite-element based methods. The aim of the present paper is precisely to reformulate it in the frame work of the Spectral Element Method (SEM) and to make it easy for every one. The detailed numerical theory and implementation are given in section II. Numerical results are shown in section III. We restrict our analysis to TM polarized waves and apply it to metallic gratings in the optical domain. This limited area of investigation is actually not a limitation. Indeed all the experience accumulated on the grating theory shows that if a method works well in this particular case, it also works in all other cases [START_REF] Popov | Gratings:Theory and Numeric Applications[END_REF]. The converse is not true.

II. Presentation of the problem

The lamellar diffraction problem is depicted in Fig. 1. The grating is illuminated from a homogeneous medium with optical index n 1 by a monochromatic plane wave with a wavelength λ and an angular frequency ω. The wave vector k forms an angle θ with respect to the z-axis. Hence its Cartesian coordinates are k x = kα 0 = kn 1 sin(θ) and k z = -kn 1 cos(θ) with k = 2π/λ the wavenumber. The exp(iωt) time dependence will be omitted throughout this paper. The grating layer with thickness h is characterized by a piecewise homogeneous permittivity function ε(x) that is periodic in the x direction with period d. The general solution is a linear combination of two fundamental polarizations: the TE polarization and the TM polarization whose non-null components of the electromagnetic field are H z ,H x ,E y and E z ,E x ,H y respectively. We restrict our analysis to the TM polarization which is known to be the most challenging. The difficulty comes from the fact that for some combinations of wavelengths and material, some field components happen to be very confined to the interfaces and others undergo huge discontinuities. In each region, which we refer to by the subscript l, l = 1, 2, 3, H y and E x can be represented by the superposition of forward and backward eigenmodes:

H yl (x, z) = N ∑ n=1 A ± ln exp(∓iγ ln z)ψ ln (x) (1a)
E xl (x, z) = 1 ωε l (x) N ∑ n=1 ±γ l A ± ln exp(∓iγ ln z)ψ ln (x) (1b)
where A + ln and A - ln are constant amplitudes and γ ln and ψ ln are modal eigenvalues and modal eigenfuntions that are determined by the boundary-value problem:

d dx 1 ε(x) d dx + k 2 ψ ln (x) = γ 2 ln 1 ε(x) ψ ln (x) (2)
ψ ln (c -) = ψ ln (c + )

ψ ln (x + d) = τψ ln (x) p - l dψ ln dx (c -) = p + l dψ ln dx (c + ) p + l dψ ln dx (d) = τ p - l dψ ln dx (0) (3)
with τ = exp(-ikα 0 d) and p - l = p + l = 1 except in the grating region where we have p - 2 = 1/n 2 21 and p + 2 = 1/n 2 22 . The eigenvalues are deduced from their square number by:

γ ln =          γ 2 lq if γ 2 lq ∈ R + -i -γ 2 lq if γ 2 lq ∈ R - γ 2 lq with negative imaginary part ifγ 2 lq ∈ C (4)
We have to compute the reflection and transmission amplitude coefficients A + 1n and A - 3n which in turn determine the diffraction efficiencies. Assuming that the unit amplitude incident plane wave corresponds to eigenvector ψ 10 , the reflection and transmission efficiencies R q , T q are given by:

R q = γ 1q γ 10 |A + 1q | 2 T q = n 2 1 n 2 3 γ 3q γ 10 |A - 3q | 2 (5)
The determination of the constant coefficients requires the knowledge of H y and E x in every region. H yl is obtained as a linear combination of the eigenvectors of ( 2) whereas E xl is deduced from the H yl thanks to (1b). The next section is devoted to the numerical derivation of H y and E x .

III. Numerical solution i. Eigenvalue eigenvector system

i.1 Weak form of the Helmholtz equation

The numerical solution follows the well-known steps of the finite elements method [8] [9]. Firstly, the differential equation (2) is turned into an integral equation by multiplying it by some continuous function φ(x) and integrating it on computational domain Ω:

Ω φ(x) d dx 1 ε(x) d dx ψ(x)dx + k 2 Ω φ(x)ψ(x)dx = γ 2 Ω φ(x) 1 ε(x) ψ(x)dx (6)
Secondly, in order to deal with the discontinuity of ε(x), the whole domain is partitioned into a finite number Q of subdomains Ω j . Let the position of the vertical interfaces be denoted by

x j , j = 0, • • • , Q, that is Ω j = x | x ∈ [x j-1 , x j ] .
The integral over domain Ω is obtained by summation of integrals over elements Ω j . Let us consider the first term of the left hand-side of ( 6) calculated on element Ω j and let us integrate it by parts. We get:

Ω j φ(x) d dx 1 ε(x) d dx ψ(x)dx = I j+1 -I j - Ω j d dx φ(x) 1 ε(x) d dx ψ(x)dx (7)
with

I j+1 = φ(x j+1 ) 1 ε(x j+1 ) d dx ψ(x) x=x j+1 I j = φ(x j ) 1 ε(x j ) d dx ψ(x) x=x j (8)
We recognize E z in

1 ε(x) d dx ψ(x) which is continuous at x = x j , j = 1, • • • , Q -1. Hence (6) writes: Q ∑ q=1 U q + V q + I N -I 0 = Q ∑ q=1 W q (9)
with

U q = x q x q-1 d dx φ(x) 1 ε(x) d dx ψ(x)dx (10) V q = - x q x q-1 φ(x)ψ(x)dx (11) W q = x q x q-1 φ(x) 1 ε(x) ψ(x)dx (12)

i.2 Algebraic eigensystem

In finite elements methods, the algebraic eigen-system from which eigenvalues and eigenvectors are sought is obtained by assembling element specific matrices. The present subsection is devoted to the building of such matrices. In element Ω j , H j y and E j x are represented as

H j y (x) = ∑ N j +1 n=1 H j y,n B n (u) E j x (x) = ∑ N j n=1 E j x,n B n (u) (13)
Generally, spectral finite elements are based on a set of Lagrange polynomials utilizing a specific nodal distribution on the interval [-1, 1] [START_REF] Pozirikidis | Introduction to Finite and Spectral Element Methods Using Matlab[END_REF]. The originality of our approach lies in the chosen set of basis functions B m (u), which are build from Legendre polynomials in the following way:

B m (u) =    0.5(L 0 -L 1 ), m = 1 0.5(L 0 + L 1 ), m = 2 L m (u) -L m-2 (u), m >= 2 (14)
where L m is the Legendre polynomial of degree m. The reduced variable

u = 2 ∆ l (x -∆ l ) , ∆ l = x l -x l+1 (15)
results from the mapping of the element Ω j to the unit interval [-1, 1] where the Legendre polynomials are defined. The above basis sketched in Fig. 2 was already used in [START_REF] Jorgensen | Higher Order Hierarchical Legendre Basis Functions for Electromagnetic Modeling[END_REF] . It allows a straightforward assembly of elementary matrices in view of the building of the global matrix. Please, note that H y expansion contains one more basis function than that of E x . This comes from the need to enforce the continuity of H y at the border of elements during the assembly process. Let us denote by P q the transformation matrix that gives the new defined basis in terms of the Legendre polynomials. In order to make the further assembly process more readable, we swap the second and the last element in the set of B m . For a basis with N q + 1 elements, P q is a square matrix of size N q + 1 which writes:

P q =                .5 -1 0 • • • 0 .5 -.5 0 -1 . . . 0 .5 0 1 0 . . . 0 0 0 0 1 . . . -1 0 0 0 0 . . . 0 0 0 0 0 0 1 0                (16)
Now, we have to choose the testing functions φ m (x). A common choice is to take also the B m (u). Thus, we are now able to express the elements of matrices U q , V q and W q :

U q mn = 2 ∆ j 1 -1 dB m du 1 ε(u) dB n du du V q mn = ∆ j 2 1 -1 B m B n du W q mn = ∆ j 2 1 -1 B m 1 ε(u) B n du (17)
The last task is to compute the above elements. We restrict our analysis to the case where ε(u) is constant and we denote by ε j the value of ε(u) in element Ω j . Recall that the Legendre polynomials are orthogonal and that

1 -1 L m L n du = 2 2m + 1 δ mn ( 18
)
with δ mn the Kronecker symbol. Let us denote by D q the diagonal matrix whose elements are the 2/(2m + 1), m = 0, • • • Q. The calculation of elementary matrices V q and W q is straightforward and gives.

V q = ∆ q 2 (P q ) T D q P q , W q = 1 ε q V q (19)
where the letter T designates the transposition. The calculation of U q is slightly more complicated. It makes use of the matrix Dr formed by the expansion coefficients of the derivative of the Legendre polynomials onto the Legendre polynomials themselves, that is:

dL m du = N ∑ n=0 D rmn L n (u) (20)
The derivation of matrix Dr is carried over to appendix A.

Let us introduce the matrix Ṗq corresponding to the expansion of the derivative of the B m (u) onto the Legendre polynomials. We have: Ṗq = Dr q P q (21)

The above described matrices allow the calculation of elementary matrix U q in the form:

U q = ∆ q ε q ( Ṗq ) T D q Ṗq (22) i.

Assembly

The next step in the finite element procedure is the assembly of elemental matrices in order to satisfy [START_REF] Volakis | Finite Element Method for Electromagnetics,Antennas, Microwave circuits, and scattering applications[END_REF] with the matrices given by ( 19) and ( 22). Let us assume that the computational domain is divided into Q elements with nodes labeled from 0 to Q giving rise to elemental matrices with dimension N q + 1, q = 1, 2, • • • Q. For periodic structures, the dimension of the global matrix is N = ∑ Q q=1 N q . Let us consider a node q, q = 1, • • • , Q -1 which is shared by elements Ω q and Ω q+1 . The value of H y at node q is simultaneously the (N q + 1)th expansion coefficients in element Ω q and the first one in element Ω q+1 . We go from local numbering p in the P th element Ω p to global numbering n in the whole computational domain by:

n = p + Q ∑ q=1 (N q + 1) -(Q -1) (23)
Hence, in the global matrix, the N p + N p+1 non-null elements on nth line corresponding to node p writes:

X p N p +1,1 , X p N p +1,2 , • • • , X p N p +1,N p +1 + X p 11 , • • • , X p+1 1,2 , X p+1 1,3 , X p+1 1,N p+1 +1 (24)
where the X q , q ∈ (p, p + 1) refers to any of elemental matrices U q ,V q or W q .

i.4 Bloch periodic boundary conditions

For the grating under consideration, the Bloch periodic boundary conditions stipulate that H y satisfies

H y (x + d) = τH y (x), τ = exp(-iα 0 d) (25)
The above particular boundary condition can be easily enforced. Indeed, it is enough to replace in element Ω Q the last expansion basis function

B Q N Q +1 by τB 1 1 .
Here, the exponents 1 and Q refer to the element to which the basis function belongs. The set of test functions has to be modified accordingly: the test function B Q N Q +1 is replaced by τ * B 1 1 where * designates the complex conjugation. One can see that I N and I 0 in (9) are equal and thus disappear in the summation. Moreover, only the last line and the last column of X Q are modified. Let us denote by F the global matrix built from the matrices X q . Enforcing the Bloch periodic boundary conditions consists in doing the three following operations

F 11 = X 1 11 + X Q N Q +1 N Q +1 F 1 N-N Q , F 1,N-N Q +1 , • • • , F 1,N = τ * X Q N Q 1 , X Q N Q 2 , • • • , X Q N Q N Q F N-N Q 1 , F N-N Q +1 1 , • • • , F N 1 = τ X Q 1 N Q , X Q 2 N Q , • • • , X Q N Q N Q
(26) Fig. 3 illustrates the filling of the global matrices U,V and W for three elements and N q = 5, 6, 4 respectively. Finally, we have transformed (2) into a generalized algebraic eigenvector equation; ii

(U + V) H y = γ 2 WH y (27)

. E x Component

In a simple spectral problem, we would not necessarily compute the E x component. In the present problem, we have to determine the constant coefficients A ± ln , which requires to enforce Dirichlet boundary conditions for H y and E x at z = 0 and z = h. Let us consider an eigen vector in element Ω q . For the sake of conciseness, in what follows, we do not make any reference neither to the region nor to the number of the eigenvector. From (1b) and ( 13) we have:

N q ∑ n=1 E q xn B n (u) = γ ωε q N q +1 ∑ n=1 H q yn B n (u) (28)
In a second step, (28) is projected onto the set of basis function B m , m = • • • , N q which gives:

N q ∑ n=1 E q xn 1 -1 B m (u)B n (u)du = γ ωε q N q ∑ n=0 H j yn 1 -1 B m (u)B n (u)du (29)
Let us introduce the matrix Pq obtained by removing the N q th column from Pq . Hence Pq has N q + 1 rows and N q columns. Similarly let us define matrix D by removing the (N q + 1)th column and the (N q + 1)th line from matrix D. In matrix form,(29) writes:

E x q = Z q H y q (30) with E x q = [E N q x1 , E N q x2 , • • • , E N q xN q ] T (31) Z q = 1 ωε q ( Pq ) T Dq Pq -1 ( Pq ) T D q D q (32)
Matrix Z q has N q rows and N q + 1 columns. Let us denote by Z the global matrix that results from the assembly of elemental matrices Z q . It takes into account the continuity of H y at the borders of elements in the following way: the nth column in matrix Z corresponding to the connection between two consecutive elements Ω p and Ω p+1 is the concatenation of the last column of Ω p and the first column of Ω p+1 . Lastly, the quasi-periodicity of H y amounts to put the N Q entries of the last column of matrix Z Q multiplied by the quasi periodic coefficient τ at the end of the first column of matrix Z. This process is illustarted in Fig. 4.

iii. Special case of homogeneous media

As discussed in [START_REF] Granet | Fourier-matching pseudospectral modal method for diffraction gratings:comment[END_REF], we use computed eigenvectors in every region, which allows matching the different fields at the boundary of the grating without additional projection. In homogeneous media 1 and 3, exact eigen-solutions are the wellknown Rayleigh expansions. In those media, the number of basis functions should be large enough so that any computed real eigenvalues, when they exist, can be associated to a diffracted order. In that case, provided that the eigenvalue is not degenerate, up to a multiplicative constant coefficient, the associated eigenvector tends to the vector formed by the expansion coefficient of the corresponding plane wave. Indeed, in a similar way to what was done in [START_REF] Granet | Fourier-matching pseudospectral modal method for diffraction gratings:comment[END_REF], we replace the computed eigenvectors associated to real eigenvalues, and only these ones, by their Rayleigh expansion counterpart. In other words, for any α q = α 0 + q2π/d such that -1 ≤ α q ≤ 1, we expand exp(-iα q x) onto B m basis. Since the P m are themselves linear combination of Legendre polynomials, we first expand the exponential functions onto the latter:

exp(-iα q x) = N p +1 ∑ m=1 ψp mq L m (u(x)) x p-1 ≤ x ≤ x p (33)
Using the following relation

1 -1 L n (x) exp(iαx)dx = i n 2π α J n+1/2 (α) (34)
and the orthogonality of Legendre polynomials we may derive:

ψp mq = (-i) m 2m + 1 2 2π ζ q J m+1/2 ζ p q exp(-iζ p q ) (35)
with The ψ p mq such that:

ζ p q = α q ∆ p 2 (36)
exp(-iα q x) = N p +1 ∑ m=1 ψp mq B m (u(x)) x p-1 ≤ x ≤ x p (37)
are linked to the ψp mq by:

ψ p q = (P p ) -1 ψp q ψ p q =    ψ p 1q
. . .

ψ p N q    ψp q =    ψp 1q . . . ψp N q    (38)
IV. Results

In this section, we validate our code by comparing our results with already published data. The parameters of the considered gratings are collected in table 1. In the past those cases were investigated because convergence problems or instabilities occurred when the Fourier Modal Method was used to analyse them. Let us first consider Grating 1 which was studied in [START_REF] Peng | Efficient implementation of rigorous coupled-wave analysis for surface relief gratings[END_REF], [START_REF] Lalanne | Highly improved convergence of the coupled-wave method for TM polarization[END_REF] and [START_REF] Chiou | Analysis of Hightly Conducting Lamellar Gratings With Multidomain Pseudospectral Method[END_REF]. We have first to verify that our method converges to the same value as others methods.

For that purpose we have plotted computed efficiencies as a function of the inverse of the truncation number which is defined as the dimension of the matrix from which eigenvalues are sought. Fig. 5 illustrate the convergence of the zeroth transmitted order computed with the SEM and with the Fourier Matching Pseudo Spectral Modal Method (PSMM) [START_REF] Song | Fourier-matching pseudospectral modal method for diffraction gratings[END_REF], [START_REF] Granet | Fourier-matching pseudospectral modal method for diffraction gratings:comment[END_REF]. We found that both converged to T 0 = .698305 which becomes the assumed "exact" value. The same figure also shows the error in the computation of T 0 defined as:

Err = log(|T 0 -T 0exact |) (39)
We made the same kind of simulations for a lossless metallic grating for which numerical data was available [START_REF] Chiou | Analysis of Hightly Conducting Lamellar Gratings With Multidomain Pseudospectral Method[END_REF], [START_REF] Foresti | Modal method in deep metal-dielectric gratings:the decisive role of hidden modes[END_REF]. By running our codes, including Fourier Modal Method with Adaptive Spatial Resolution, (ASR-FMM), [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF], with huge truncation orders, we found that the "exact" zeroth reflected order was R 0 = .127698. The computed values are summarized in table 2 and table 3. From a practical view point, it is of course unnecessary to use such large truncation numbers. On the one hand, our aim is to verify the stability and the robustness of the method and on the other hand to determine a numerical "exact" value that allows to calculate an error.

The convergence of any modal methods applied to gratings depends on the accuracy with which eigenvalues and eigen modes are computed on the one hand and on the matching of plane waves at the border of the grating region on the other hand. Or, for lamellar grating, the junction between the grating region itself and the homogeneous regions in which it is embedded creates a wedge with three different materials and a right angle. It is well-known that in Tm polarization the transverse components of the electric field and magnetic field may be singular at the sharp edge of the wedge. When the singularity is not incorporated into the algorithms as is the case in this paper, the accuracy of the numerical method can be tested by verifying whether computed data behaves as predicted in the vicinity of singular points. With reference to Fig. 1, it can be shown that at z = h, x > c, E x has the following asymptotic form:

E x = A × (x -c) τ-1 (40)
where xc is the distance from the sharp edge. The singularity exponent τ has the following expression in the special case where n 1 = n 21 [START_REF] Li | Field singularities at lossless metaldielectric right-angle edges and their ramifications to the numerical modeling of gratings[END_REF]:

τ = 2 π arcsin 1 - 1 4 n 2 22 -n 2 1 n 2 22 + n 2 1 1/2 (41)
For grating 1, we have τ = .666 + .0335i and for grating 2 τ = .635. The smaller the τ value, the larger the truncation number. Indeed the comparison of Fig. 5 and Fig. 6 shows that the numerical convergence of the efficiencies is faster for grating 1 than for grating 2. However, the difference in the value of τ is probably not the main reason that explains the speed of convergence between the two gratings configurations. In Fig. 7, the magnitude of E x is plotted against the distance from the wedge in logarithmic scale, for a truncation order of 300. The distance is expressed in multiple of the wavelength. Fitting the plot with a straight line of the form y = a × x + b between -4.5 and -3 gives a = .33 in very good accordance with the theoretical value .3340. Other simulations, not shown here, with lower truncation orders give the same behavior but the parameter a is smaller and varies around 0.29. One can conclude that the method reproduces with a satisfactory accuracy the near field.

V. Conclusion

We have developed a code based on the spectral element method in which hierarchical basis functions derived from Legendre polynomials were introduced. The code has been validated by comparison with other well established methods on cases for which obtaining numerical convergence was considered challenging. On these examples, the SEM demonstrated accuracy and stability. The formulation which we have just presented enriches the broad range of the numerical modal methods. However compared to the pioneering work of Morf [START_REF] Morf | Exponentially convergent and numerically efficient solution of Maxwell's equations for lamellar gratings[END_REF], who was the first to use Legendre polynomials expansions in numerical modal methods, it is our opinion that the SEM is simpler to implement. In the future, it would be interesting to benchmark honestly the many numerical methods with appropriate metrics in order to better understand their optimal field of use.

Figure 1 :

 1 Figure 1: Typical lamellar grating.

Figure 2 :

 2 Figure 2: Illustration of the first four modified Legendre polynomials

Figure 3 :

 3 Figure 3: Example of assembly process of three elementary matrices for the solution of eigen-vector equation 6 with pseudo-periodic boundary conditions. The cardinality of elements is 6,7 and 5 respectively.

Figure 4 :

 4 Figure 4: Example of assembly process of three elementary matrices for the derivationn of the E x electric field component given by Eq32. H y has 6, 7 and 5 spectral coefficient in elements Ω 1 , Ω 2 and Ω 3 , respectively whereas E x has one less.

Figure 5 :

 5 Figure 5: Zeroth transmitted order of grating 1 against the inverse of the truncation order computed with SEM and with PSFM. It is observed that both methods converge to the same value. The inset gives the error of T 0 .

Figure 7 :

 7 Figure 7: Magnitude of E x at the upper surface of grating 1. The straight line fits the data between -4.5 and -3, that is close to the wedge. Its equation is y = (.3291 × x -.48).

Table 1 :

 1 Parameters of investigated dielectric metallic gratings

	Parameter Grating 1 Grating 2
	θ	0		35
	λ	.55		.45
	c d	.3 × .25 .25		.5 .2
	h	.2		.521
	n 3	1.5		1.5
	n 21 n 22	1 3.18-4.41j	-	1 √ 25j

Table 2 :

 2 Specular Reflected Efficiency for the lossless grating 2 computed with different modal methods.TN is the dimension of the matrix from which eigenvalues are sought.

	TN	SEM	PSMM[14]	AMM[15]
	33	1.2955593e-01 1.3059319(-1) 1.2544877(-1)
	65	1.2812373e-01 1.2835611(-1) 1.2694599(-1)
	129	1.2779410e-01 1.2787523(-1) 1.2745205(-1)
	257	1.2771968e-01 1.2771968(-1) 1.2761832(-1)
	513	1.2770294e-01 1.2769345(-1) 1.2767238(-1)
	1025 1.2769926e-01 1.2768770(-1) 1.2768988(-1)

Table 3 :

 3 Specular Reflected Efficiency for the lossless grating 2. Benchmark between three in-house codes with huge truncation order. TN is the dimemsion of the matrix from which eigenvalues are sought.

		0.138					PSMM
	Reflected Efficiencie R 0	0.13 0.132 0.134 0.136	Error of R 0	-6 -5 -4 -3 -2	0 400 800 1200 1600 Truncation Number SEM PSMM	SEM
		0.128
		0.126
		0			0.01	0.02	0.03	0.04	0.05
					Inverse of Truncation Number

Figure 6: Specular reflected order of grating 2 against the inverse of the truncation order computed with SEM and with PSFM. It is observed that both methods converge to the same value. The inset gives the error of R 0 .

Appendix A: derivation of the derivative operator associated with Legendre polynomial expansions

Consider the expansion of functions F(u) and its derivative G(u) onto the set of N + 1 Legendre polynomials :

We seek the elements of the D matrix which relate the expansion coefficients of the derivative of the F(u) function to the coefficients of the function itself that is:

by deriving F(u), we are led to:

The projection of the above equation onto Legendre polynomial L m (u) gives:

Taking into account the two following properties of Legendre polynomials:

1

where δ mn is the Kronecker symbol we get:

with

VI. Acknowledgments

This work was sponsored by a public grant overseen by the French National Agency as part of the "Investissements d'Avenir" through the IDEX-ISITE initiative CAP 20-25 (ANR-16-IDEX-0001).

VII. Disclosures

The authors declare no conflicts of interest.