Samuel Nibouche
email: samuel.nibouche@cirad.fr.

Laurent Costet

Raul F Medina

Jocelyn R Holt

Joëlle Sadeyen

Anne-Sophie Zoogones

Paul Brown

Roger L Blackman

Laurent Costet

Morphometric and Molecular Discrimination of the sugarcane aphid, Melanaphis sacchari (Zehntner, 1897) and the sorghum aphid Melanaphis sorghi (Theobald, 1904)

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

Introduction

The species Melanaphis sacchari [START_REF] Zehntner | Overzicht van de Zieken van het Suikerriet op Java[END_REF] and Melanaphis sorghi [START_REF] Theobald | The 'Dura' Aphis or 'Asal Fly[END_REF] were described at the turn of the 20 th century by Zehntner [START_REF] Zehntner | Overzicht van de Zieken van het Suikerriet op Java[END_REF] on sugarcane in Java and by Theobald [START_REF] Theobald | The 'Dura' Aphis or 'Asal Fly[END_REF] on sorghum in Sudan respectively. However, although these two species are commonly treated as synonyms, referred to as the 'sugarcane aphid', no comparative study demonstrating this synonymy has been conducted. When Remaudière & Remaudière [3] considered M. sorghi as a synonym of M. sacchari in their 1997 catalogue, following [START_REF] Eastop | Notes on East African Aphids I-Synonomy[END_REF] [START_REF] Eastop | Notes on East African Aphids I-Synonomy[END_REF], they provided no reference to support this choice. Moreover, Halbert and Remaudière [START_REF] Halbert | A new oriental Melanaphis species recently introduced in North America [Hemiptera, Aphididae][END_REF] referred later to these species as the 'sorghi/sacchari group' and underlined them as 'two very variable species usually regarded as synonyms, but possibly distinct according to [START_REF] Blackman | The biology and taxonomy of the aphids transmitting barley yellow dwarf virus[END_REF]'. The communication by Blackman et al. [START_REF] Blackman | The biology and taxonomy of the aphids transmitting barley yellow dwarf virus[END_REF] indeed supported the separation of M. sorghi and M. sacchari and provided a morphological criterion to separate the species, based on the ratio between the hind tibia length and the antennal processus terminalis length. The same separation was used in Blackman and Eastop's book in 2006 [START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF] and has not been challenged since.

Margaritopoulos et al. [START_REF] Margaritopoulos | Two aphid species newly introduced in Greece[END_REF] stated that 'DNA evidence that might confirm the existence of two species is not yet available, but at this time it seems advisable to recognise that they are probably functioning as distinct taxonomic entities'. Nibouche et al. [START_REF] Nibouche | Low genetic diversity in Melanaphis sacchari aphid populations at the worldwide scale[END_REF] observed genetic structuring in five clonal lineages matching a geographic structure, but they could not separate the two species by using 'universal' COI barcoding.

Regarding host plant association, each species is observed on both sugarcane and sorghum, but M. sorghi is considered preferring sorghum and M. sacchari preferring sugarcane [START_REF] Blackman | The biology and taxonomy of the aphids transmitting barley yellow dwarf virus[END_REF]. However, host plant association is blurred by the existence of biotypes, as shown by Nibouche et al. [START_REF] Nibouche | Host Plant Specialization in the Sugarcane Aphid Melanaphis sacchari[END_REF] who demonstrated the existence of a sorghum and a sugarcane biotype in Reunion populations (within the same multi locus lineage). Interestingly, Melanaphis sorghi is known for a long time to produce very heavy infestations on sorghum, in Africa. Early in the 20 th century, [START_REF] Vuillet | Les pucerons du sorgho au Soudan Français[END_REF] [START_REF] Vuillet | Les pucerons du sorgho au Soudan Français[END_REF] cited this aphid as responsible for famines in West Africa.

The objectives of this study were [START_REF] Zehntner | Overzicht van de Zieken van het Suikerriet op Java[END_REF] to improve the description of the morphometric differences between M. sorghi and M. sacchari, and (2) to delimit the molecular separation of M. sorghi and M. sacchari.

Material and Methods

Material collected

This study is based on 199 samples, collected from 2002 to 2016 in 31 states or countries (S1 Table ). No specific permissions were required when sampling aphids in the locations studied. From these 199 samples, we analyzed 2,409 apterous aphid specimens collected on cultivated sorghum (Sorghum bicolor) (n = 439), Johnson grass (Sorghum halepense) (n = 97), sugarcane (Saccharum officinarum x S. spontaneum) (n = 1,382), Sorghum arundinaceum (≡ Sorghum bicolor verticiliflorum ≡ Sorghum verticiliflorum) (n = 427), maize (Zea mays) (n = 2), perennial sorghum (Sorghum × almum) (n = 8) , Sorghum sp. (n = 16), and pearl millet (Pennisetum glaucum) (n = 38). Most of this material was already analyzed in previous studies [START_REF] Nibouche | Low genetic diversity in Melanaphis sacchari aphid populations at the worldwide scale[END_REF][START_REF] Nibouche | Host Plant Specialization in the Sugarcane Aphid Melanaphis sacchari[END_REF][START_REF] Nibouche | Invasion of sorghum in the Americas by a new sugarcane aphid (Melanaphis sacchari) superclone[END_REF]. Using Blackman & Eastop's ( 2006) key, we identified these specimens as M. sacchari or M. sorghi, but the identification sometimes was ambiguous given the continuous and overlapping nature of the criteria used to separate both species. For convenience, while awaiting this taxonomy study, this material was referred to as M. sacchari in our previous papers.

DNA extraction

DNA of individual aphids was extracted using the 'salting-out' protocol of Sunnucks and Hales [START_REF] Sunnucks | Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae)[END_REF] or using the Qiagen DNeasy Blood & Tissue Kit (Qiagen, Courtaboeuf, France). For slide-fixed specimens, a non-destructive DNA extraction was performed using the Qiagen manufacturer's protocol, but retrieving the insect body from the first elution column [START_REF] Nibouche | Invasion of sorghum in the Americas by a new sugarcane aphid (Melanaphis sacchari) superclone[END_REF].

Microsatellites

Nine microsatellite (Single Sequence Repeat, SSR) markers were selected among 14 previously developed markers [START_REF] Andris | Permanent Genetic Resources added to Molecular Ecology Resources Database 1 June 2010 -31[END_REF]. PCR reactions were performed with labelled primers and multiplexed following previously established protocols [START_REF] Nibouche | Low genetic diversity in Melanaphis sacchari aphid populations at the worldwide scale[END_REF]. Genotyping was carried out using an ABI PRISM 3110 and alleles were identified at each locus by comparison with a size standard using Gene-Mapper version 2.5 software (Applied Biosystems). The total microsatellite genotyping dataset included 2,255 specimens: 2,175 specimens were previously analyzed in [START_REF] Nibouche | Low genetic diversity in Melanaphis sacchari aphid populations at the worldwide scale[END_REF][START_REF] Nibouche | Host Plant Specialization in the Sugarcane Aphid Melanaphis sacchari[END_REF][START_REF] Nibouche | Invasion of sorghum in the Americas by a new sugarcane aphid (Melanaphis sacchari) superclone[END_REF] and 80 additional specimens were genotyped for this study. Single combinations of alleles were characterized and arranged as distinct multilocus genotypes (MLG) and assigned to one of six multilocus lineages (MLL) [START_REF] Nibouche | Low genetic diversity in Melanaphis sacchari aphid populations at the worldwide scale[END_REF][START_REF] Nibouche | Invasion of sorghum in the Americas by a new sugarcane aphid (Melanaphis sacchari) superclone[END_REF]. We carried out a Bayesian clustering analysis with Structure version 2.3.4 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF], results were summarized with Structure Harvester [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF], Clumpp [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF] and Distruct 1.1 [START_REF] Rosenberg | DISTRUCT: a program for the graphical display of population structure[END_REF]. Parameters of the Structure analysis were: admixture, independent allele frequencies, 100,000 iterations after a 25,000 burn-in period, 10 replications for each k value ranging from 1 to 8.

DNA sequencing

Aphids were sequenced for three mitochondrial and three nuclear DNA regions belonging to the cytochrome c oxidase subunit I (COI) [START_REF] Folmer | DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates[END_REF], cytochrome c oxidase II (COII) [START_REF] Kim | A molecular phylogeny of the tribe Aphidini (Insecta: Hemiptera: Aphididae) based on the mitochondrial tRNA/COII, 12S/16S and the nuclear EF1α; genes[END_REF], cytochrome b (CytB) [START_REF] Harry | Molecular evidence for parallel evolution of adaptive syndromes in fig-breeding Lissocephala (Drosophilidae)[END_REF], elongation factor-1α (EF1-α) [START_REF] Von Dohlen | Phylogenetics and evolution of the eastern Asianeastern North American disjunct aphid tribe, Hormaphidini (Hemiptera: Aphididae)[END_REF] and the internal transcribed spacer 1 and 2 (ITS1 and ITS2) [START_REF] Ji | Evolutionary conservation and versatility of a new set of primers for amplifying the ribosomal internal transcribed spacer regions in insects and other invertebrates[END_REF].

Three of the six sequences, COI, COII and EF1-α, produced informative polymorphism and were simultaneously sequenced on a large number of specimens. These sequences were concatenated and the resulting concatenated haplotypes were used to produce a minimum spanning network with PopArt [START_REF] Leigh | POPART: full-feature software for haplotype network construction[END_REF].

Morphometry Material

The morphometry dataset consisted of 89 apterous female specimens that were slidefixed (S1 Fig) after a non-destructive DNA extraction. Among these 89 specimens, 21 were successfully genotyped with SSR and could be assigned to a multilocus lineage (MLL). The remaining 68 were not genotyped, but belonged to a sample whose multilocus lineage (MLL) had been identified from other specimens. Since 186 of the 188 (98.9%) SSR genotyped samples in this study were homogeneous (i.e. contained only a single MLL), we assumed that the non-genotyped slide-fixed samples belonged to the same MLL as the other specimens in the sample that were genotyped. The distribution of the slide-fixed specimens was: 22 MLL-A, 2 MLL-B, 19 MLL-C, 15 MLL-D and 31 MLL-F. Since only two specimens of MLL-B were observed, they were discarded from the statistical analysis. No MLL-E specimen was observed.

We also examined the Theobald type series of Aphis sacchari collected in Sudan in 1902, which is stored in the Natural History Museum of London (NHM). In this series, only one paratype specimen (NHM-1915-81) was an apterous female and could be used for our morphometric analysis. Unfortunately, we were unable to locate the types of Melanaphis sacchari. Hollier & Hollier [START_REF] Hollier | Leo Zehntner, Swiss pioner of tropical applied entomology[END_REF] reported that a fire destroyed the experimental station of Salatiga in Java in 1902, including Zehntner's laboratory and his collections. As Zehntner did not send his types to the Geneva museum, the type of Melanaphis sacchari, collected in 1897, probably was destroyed and should be considered lost.

We also examined some alate M. sorghi (n = 5) and M. sacchari (n = 5) specimens on slides from the MNHN (Muséum National d'Histoire Naturelle, Paris; G. Remaudière collection) and from the GBGP (Centre de Biologie et de Gestion des Populations, Montpellier; F. Leclant collection). The Theobald's type series from NHM contained one alate specimen (NHM-1915-81), but it was of insufficient quality to be included in the morphometric characterization of alatae.

Methods

Twenty-two characters that are used classically in aphid taxonomy [START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF][START_REF] Blackman | Aphids on the world's plants. An online identification and information guide[END_REF] were measured on each slide-fixed specimen, using a binocular lens stereo microscope. Paired appendages (i.e. legs, antenna, siphunculi), were measured on both sides, and the mean value was used for analysis (except for the number of setae on the antenna, which was observed on one side only). Twelve ratios were computed from these characters. The analysis was carried out on a subset of 11 characters (Table 1) that are known to be discriminant between species within the Melanaphis genus [START_REF] Halbert | A new oriental Melanaphis species recently introduced in North America [Hemiptera, Aphididae][END_REF][START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF][START_REF] Blackman | Aphids on the world's plants. An online identification and information guide[END_REF][START_REF] Kim | Taxonomic review and morphometric analysis of the genus Melanaphis van der Goot (Hemiptera: Aphididae) in Korea[END_REF]. Because of missing data, only 50 specimens were used in the multivariate analysis. Because Theobald did not clarify his M. sorghi paratype specimen (NHM-1915-81), we were unable to observe urs and siphBW, causing three missing ratios in the dataset. To include the paratype in the analysis, we replaced these three missing data by zeros after the standardization step (see below). The data were first submitted to a discriminant analysis of principal components (DAPC) with the R package ADEGENET using the find.clusters function [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF]. This method first uses an unsupervised k-means clustering approach to determine the number of clusters without requiring any a priori clustering information. The determination of the number of clusters was based on the Bayesian Information Criteria (BIC). Then in a second step, a principal component analysis (PCA) is carried out, followed by a canonical discriminant analysis (CDA) performed on the coordinates along the principal components. To verify if the clustering was congruent with the separation of M. sorghi and M. sacchari, we compared the HindTibia:pt ratio between the clusters delimited by ADEGENET. Indeed, M. sorghi specimens have a relative length of the processus terminalis shorter than M. sacchari specimens [START_REF] Blackman | The biology and taxonomy of the aphids transmitting barley yellow dwarf virus[END_REF][START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF]. According to [START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF][START_REF] Margaritopoulos | Two aphid species newly introduced in Greece[END_REF][START_REF] Blackman | Aphids on the world's plants. An online identification and information guide[END_REF], the HindTibia:pt ratio range for apterae is (2.0 -3.0) for M. sorghi vs. (1.4 -2.2) for M. sacchari.

Secondly, the complete morphological dataset (i.e. 34 morphological traits on 88 apterous specimens) was submitted to a one-way ANOVA with SAS PROC GLM [29] to detect significant differences between M. sacchari and M. sorghi. Because we carried out multiple analysis, to control the study-wise type-1 error level we used a 5% False Discovery Rate (FDR) approach [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] with SAS PROC MULTTEST [29] to detect significant differences between clusters. To carry out this analysis, slide specimens were assigned to M. sacchari or M. sorghi according to their molecular assignment using SSR and EF1-α data.

Results

Morphometry

According to the unsupervised k-means clustering, the number of inferred morphological clusters was k = 2. The membership probability of each specimen is shown by Fig 1 . The M. sorghi paratype was assigned to the morphological cluster 1, with a 100% membership probability. The loading plot (Fig 3 ) showed that the morphological characters contributing most to the DAPC were pt:cauda, HindTibia:pt and pt:siph. The graphical comparison of these three main contributing traits between both species is presented in Fig 2 . In these plots, the assignment of the specimens to both species was made according to their SSR and EF1-α genotype (see below). For the three traits, both species differ significantly with an uncorrected P-value < 0.0001 (F = 100.3, 78.73 and 62.96 respectively for pt:cauda, HindTibia:pt and pt:siph). However, some overlapping is observed for each trait. The comparison of all morphological traits of M. sorghi and M. sacchari apterous females is presented in S2 Table . There are significant differences between both species in 16 out of 34 traits.

The values for alatae are given in S3 Table. Due to the small number of specimens observed, no statistical comparison was carried out and these values are only suggestive. Because we only observed museum slide-fixed alate specimens, we could not genotype them. In the absence of a key for alatae, species assignment was carried out taking into account geographic origin, according to the worldwide repartition of both species (Fig 5). Specimens from West Africa (Burkina Faso, Ivory Coast and Senegal) were considered as M. sorghi, specimens from Reunion and Brazil were considered as M. sacchari. Although the number of alate specimens is insufficient to form a basis for any statistical comparison, it should be noted that, as in apterae, the pt of alatae assigned to M. sacchari tends to be longer than that of alatae assigned to M. sorghi, and that the ratios involving this parameter, especially pt:siph, may prove to be useful discriminants. 4). The eight concatenated haplotypes (cH1 to cH8) are listed in Table 4. The number of hatch marks represents the number of mutations separating the concatenated haplotypes. Circle sizes are proportional to haplotype frequencies.

Microsatellites

Fifty-nine multilocus genotypes (MLG) have been identified (S6 Table ). Fifty-six were published earlier [START_REF] Nibouche | Low genetic diversity in Melanaphis sacchari aphid populations at the worldwide scale[END_REF][START_REF] Nibouche | Host Plant Specialization in the Sugarcane Aphid Melanaphis sacchari[END_REF][START_REF] Nibouche | Invasion of sorghum in the Americas by a new sugarcane aphid (Melanaphis sacchari) superclone[END_REF] and three new ones (Ms25, Ms26, and Ms58) were observed during this study. The clustering with Structure, followed by the use of the Evanno et al. [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] method, leads to the conclusion that the number of inferred populations is k = 2 (S2 Fig) .    The assignment of each MLG with k = 2 is presented in sacchari and the Structure Bayesian clustering based on SSR genotyping showed unambiguously that the blue cluster is M. sorghi and the orange cluster is M. sacchari (Table 2). Among the 29 specimens belonging to the blue Structure cluster, 27 were morphologically assigned to M. sorghi, one was morphologically assigned to M. sacchari and one was undetermined. Among the 21 specimens belonging to the orange Structure cluster, 19 were morphologically assigned to M. sacchari, one was morphologically assigned to M. sorghi and one undetermined. The resulting accuracy of the congruence of Structure clustering and morphological clustering is 95.8% (46 / 48), when excluding the two undetermined specimens. The specimens with a DAPC membership lower than 0.8 are considered undetermined. b MLL identification was carried out by SSR genotyping on the specimen itself or on other specimens from the same sample. c Specimens whose molecular and morphometric assignments are discordant.

DNA sequencing

A total of 371 aphids were sequenced for at least one of three genes, COI, COII or EF1-α: 340 for COI, 143 for COII, and 163 for EF1-α.

In COI (658 bp), we detected four SNPs defining four haplotypes (S4 Table ). The sequences were deposited in Genbank under accession numbers KJ083108-KJ083215, KX453783-KX453784, MG838208-MG838315, MT813521-MT813656, and MT821458-MT821474.

In COII (763 bp), we detected four SNPs defining five haplotypes (S5 table ). The sequences were deposited in Genbank under accession numbers MT847245-MT847387.

The EF1-α gene portion amplified was 1,014 bp long. The sequences were deposited in Genbank under accession numbers MT847432-MT847594. Prior to analysis, the sequences were trimmed to a 528 bp length (from position 248 bp to 775 bp), to discard low quality 5' and 3' sequence ends in most specimens. In this EF1-α 528 bp sequence portion, we detected 1 unambiguous SNP and 4 ambiguous positions (heterozygous) defining seven haplotypes (S6 Table ). Only two haplotypes were defined when omitting the four ambiguous positions, and we considered only these two haplotypes in further analysis. Haplotype H1 was only observed in M. sacchari and haplotype H2 only in M. sorghi.

In CytB (745 bp), we detected one SNP defining two haplotypes. Both haplotypes were detected in both M. sacchari and M. sorghi. The sequences were deposited in Genbank under accession numbers MT847388-MT847423.

In the ITS1 region, we obtained 445-451 bp length sequences. These sequences included a 3-6 bp indel region and one SNP. The two haplotypes defined by the SNP were present in both M. sacchari and M. sorghi specimens. The sequences were deposited in Genbank under accession numbers MT821305-MT821342.

In the ITS2 region, in a 462 bp-sequence length, we detected one SNP and three indels, but the SNP was located within an indel and therefore could not be used in further analysis. The sequences were deposited in Genbank under accession numbers MT821344-MT821448.

A total of 63 specimens were genotyped at the three genes COI, COII, EF1-α, and were also successfully genotyped with SSR. The relationship between MLL and haplotypes among these specimens is presented in Table 3. EF1-α provided a diagnostic substitution at position 637 separating M. sorghi from M. sacchari. The distinctive base was T for M. sacchari and A for M. sorghi. COI provides an incomplete separation of both species: haplotypes H2 and H3 are diagnostic of M. sacchari, H6 is diagnostic of M. sorghi, but haplotype H1 is present in both species. The number of variable sites was nine on the 1,806 bp of the three concatenated gene sequences (i.e. 0.49%), defining eight haplotypes, with a nucleotide diversity of 0.13%. The minimum spanning network among haplotypes is shown by Fig 5 . M. sorghi specimens defined a star shaped haplogroup centered on haplotype cH1. M. sacchari defined a less homogeneous haplogroup, each haplotype being separated from the other by two or three substitutions. we did not encounter in our study or are H1 haplotype, which is uninformative because it is shared by M. sorghi and M. sacchari. In West and in Southern Africa, M. sorghi is the sole detected species. In East Africa, M. sacchari has been detected in Kenya and Tanzania, and M. sorghi in Uganda and Kenya. In Kenya, both species coexist and were collected in the same sample once. Reunion and Mauritius, in the South West Indian Ocean, are exclusively colonized by M. sacchari. The Neotropical zone is colonized by M. sacchari. The Nearctic zone is colonized by both species, as a result of the recent introduction of M. sorghi in the Americas. In Asia, M. sorghi is present in China and India, and M. sacchari in Cambodia. In Australia and in Hawaii, only M. sacchari was detected.

Discussion

Genetic analyses with SSRs and three gene sequences showed that two genetic clades exist, one grouping MLL-A-E-F and the other grouping MLL-C-D. The multivariate morphometric data analysis separated the specimens in two groups matching the two genetic clades. Comparison with a M. sorghi paratype and comparison of the HindTibia:pt ratio [START_REF] Blackman | The biology and taxonomy of the aphids transmitting barley yellow dwarf virus[END_REF][START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF][START_REF] Blackman | Aphids on the world's plants. An online identification and information guide[END_REF] confirmed that the genetic clade grouping MLL-A-F is M. sorghi while the clade grouping MLL-C-D is M. sacchari. The status of MLL-B (Australian specimens) remains to be confirmed: it is assigned to M. sorghi by SSRs and EF1-α sequence, but we did not confirm this assignment by morphometric means, due to a lack of specimens. Similarly, the lack of specimens prevented us from studying the morphology of MLL-E (all from China) and its taxonomic status remains uncertain, although SSRs and EF1-α sequence data both suggest it belongs to M. sorghi.

Three morphological criteria are useful for species separation: pt:cauda, pt:siph and HindTibia:pt. However, as observed by Blackman and Eastop [START_REF] Blackman | Aphids on the world's herbaceous plants and shrubs[END_REF][START_REF] Blackman | Aphids on the world's plants. An online identification and information guide[END_REF], there are no clear limits between species and values overlap largely (Fig 2 The molecular diagnostic methods for separation of M. sacchari and M. sorghi are summarized in Figs 7 and8. COI can allow for the identification of M. sacchari through two specific SNPs in positions 263 or 294. But, because haplotype H1 (Table 3) is shared by both species, an unambiguous identification of M. sorghi with COI is sometimes not possible. With EF1-α, a specific SNP in position 637 allows the separation of both species. Genotyping with a single SSR locus also allows the separation of both species. For example, the SSR locus CIR-Ms-G01 can be used for this purpose (Fig 8). The two alleles present at this locus in M. sacchari are separated by 25 to 31 bp, whereas M. sorghi genotypes are mostly homozygous or exhibit two alleles separated by 4 to 8 bp only (S7 table). ) and A5 (voucher # SNIB000237_0102). The PCR was carried out according to [START_REF] Andris | Permanent Genetic Resources added to Molecular Ecology Resources Database 1 June 2010 -31[END_REF]. The migration of PCR products was carried out on a Qiagen Qiaxcel electrophoresis analyzer. The image was generated by the Qiaxcel ScreenGel 1.6.0 software. We observed low genetic distances between M. sacchari and M. sorghi when comparing gene sequences that are widely used in aphid taxonomy. A 'borderline' distance between species in the COI barcode region is sometimes used by some authors to confirm species limits [START_REF] Cocuzza | Identification of aphids of Aphis frangulaegroup living on Lamiaceae species through DNA barcode[END_REF][START_REF] Rakauskas | Mitochondrial COI and morphological evidence for host specificity of the black cherry aphids Myzus cerasi (Fabricius, 1775) collected from different cherry tree species in Europe (Hemiptera, Aphididae)[END_REF]. However, it is now well recognized that there is no universal genetic distance separating aphid species and that low genetic difference between species can occur in aphids. Multiple examples of low COI, COII, CytB pairwise distance between species exist [START_REF] Lee | Barcoding aphids (Hemiptera: Aphididae) of the Korean Peninsula: updating the global data set[END_REF][START_REF] Lee | Developing a new molecular marker for aphid species identification: Evaluation of eleven candidate genes with species-level sampling[END_REF][START_REF] Chen | The effectiveness of three regions in mitochondrial genome for aphid DNA barcoding: a case in Lachininae[END_REF]. For example, a situation very similar to ours occurs in the genus Megoura, where Kim and Lee [START_REF] Kim | Molecular systematics of the genus Megoura (Hemiptera: Aphididae) using mitochondrial and nuclear DNA sequences[END_REF] observed an absence of COI, COII, CytB, ITS1, ITS2 divergence between Megoura litoralis Müller and Megoura viciae Buckton, which differ by a 0.2% P-distance on EF1-α.

According to Blackman et al. [START_REF] Blackman | The biology and taxonomy of the aphids transmitting barley yellow dwarf virus[END_REF], M. sorghi is more likely observed on sorghum and M. sacchari more likely on sugarcane, although not absolutely specific to the hosts indicated by their names. In our study, the samples were distributed as follows: 14 from sugarcane vs. 94 from sorghum (all Sorghum species together) for M. sorghi, and 18 samples from sorghum vs. 69 from sugarcane for M. sacchari. Although our sampling plan was not designed to test host plant preference, the difference of preference between the two species appears obvious. The results of Boukari et al. [START_REF] Boukari | Lack of transmission of Sugarcane yellow leaf virus in Florida from Columbus grass and sugarcane to sugarcane with aphids or mites[END_REF] obtained in Florida confirm this preference, showing that sugarcane harbors almost only M. sacchari (COI haplotypes H2 and H3), which is absent from aphids on Sorghum spp. The recent work by Paudyal et al. [START_REF] Paudyal | Evidence of host plant specialization among the US sugarcane aphid (Hemiptera: Aphididae) genotypes[END_REF] in the USA also supports this apparent preference. Indeed, using host transfer experiments, these authors have demonstrated that MLL-F strains collected from Sorghum spp. exhibited a higher fitness on sorghum than on sugarcane, and that an MLL-D strain collected from sugarcane exhibited a higher fitness on sugarcane than on sorghum. According to this host preference difference, we suggest that the common name 'sugarcane aphid' should be used for M. sacchari and 'sorghum aphid' for M. sorghi.

Our study shows that the invasive genotype responsible for outbreaks on sorghum in North and Central America and the Caribbean islands since 2013 is MLL-F [START_REF] Nibouche | Invasion of sorghum in the Americas by a new sugarcane aphid (Melanaphis sacchari) superclone[END_REF], which, belongs to M. sorghi, while the genotype present before 2013 (MLL-D) is M. sacchari. If the hypothesis of a lower fitness of M. sacchari on sorghum compared to M. sorghi is confirmed, this would explain why no damage was observed on sorghum prior to the introduction of M. sorghi to the Americas.

Fig 1 .

 1 Fig 1. DAPC analysis based on 11 morphological traits recorded on 51 slide-fixed

Fig 3 .

 3 Fig 3. Canonical loading plot. The horizontal line is the limit showing the variables

Fig 5 .

 5 Fig 5. Minimum spanning network constructed using the concatenated COI, COII and

Fig 4 .

 4 Fig 4. Assignment of each of the 59 MLG to the two clusters inferred by Structure. MLG

Fig 6 .

 6 Fig 6. Molecular identification with SSRs and sequencing of COI or EF1-α of 2,332

Fig 7 .Fig 8 .

 78 Fig 7. Molecular diagnosis for separation of M. sacchari and M. sorghi using sequencing of COI or EF1-α. Fig 8. Molecular diagnosis for separation of M. sacchari and M. sorghi using the SSR

Table 1 . List of the 11 morphological variables used in the discriminant analysis of principal components (DAPC) analysis.

 1

	Variable name Variable signification
	NsetaeCauda number of setae on the cauda
	pt:cauda	ratio processus terminalis length / cauda length
	HindTibia:pt	ratio hind tibia length / processus terminalis length
	Ant:BL	ratio antenna length / body length
	urs:htII	ratio ultimate rostral segment length / hind tarsa II length
	pt:VIb	ratio processus terminalis length / base length of the 6th antennal segment
	pt:siph	ratio processus terminalis length / siphunculi length
	cauda:urs	ratio cauda length / ultimate rostral segment length
	siph:BL	ratio siphunculi length / body length
	siph:siphBW	ratio siphunculi length / siphunculi basal width
	siph:cauda	ratio siphunculi length / cauda length
	The 11 variables were standardized to obtain a mean of zero and a standard deviation
	of 1 prior to analysis.

Table 2 . Multi Locus Lineage (MLL) identification of the 51 slide fixed specimens morphologically assigned to M. sorghi or M. sacchari by DAPC.

 2

	Structure	MLL	level of MLL	morphological assignment by DAPC
	cluster		identification b
				M. sacchari M. sorghi undetermined a
	orange	MLL-C	specimen	1		1
			sample	9	1 c
		MLL-D	specimen	8
			sample	1
	blue	MLL-A	specimen
			sample		16
		MLL-F	specimen	1 c	7
			sample		4	1
	a

Table 3 . Haplotypes defined by the concatenation of four genes and correspondence with the Multilocus Lineages (MLLs) defined using SSR.

 3

	concatenated	gene / haplotype	species		Structure
	haplotype					MLL	Cluster	n
		CO1 CO2 EF
	cH1	H1	H3	H2	sorghi	A	1	18
	cH1	H1	H3	H2	sorghi	F	1	8
	cH2	H1	H3	H1	sacchari	B	2	2
	cH3	H2	H1	H1	sacchari	C	2	11
	cH4	H2	H5	H1	sacchari	C	2	5
	cH5	H3	H3	H1	sacchari	D	2	10
	cH6	H1	H2	H2	sorghi	E	1	6
	cH7	H6	H3	H2	sorghi	F	1	1
	cH8	H1	H4	H2	sorghi	A	1	2

  ). The ranges for M. sacchari vs. M. sorghi are respectively pt:cauda (1.72-2.35) vs. (2.19-2.8), HindTibia:pt (1.8-2.31) vs. (1.54-1.93) and pt:siph (2.71-4.28) vs. (2.26-3.19). Due to overlap, the use of these morphological criteria can lead to ambiguous results and should be applied quantitatively (i.e. at the population level) rather than qualitatively (i.e. at the individual level).

Acknowledgements

We are grateful to Hughes Telismart, Magali Hoarau and Antoine Franck for technical assistance. We also thank Susan E. Halbert for her helpful comments on this manuscript. We also acknowledge Armelle Coeur d'Acier (CBGP), Laurent Fauvre and Thierry Bourgoin (MNHN) for providing access to the CBGP and MNHN collections.

The within species and between species divergences for the six sequenced gene portions are summarized in Table 4. Table 4. Sequence divergences (pairwise uncorrected P-distances, %) between or within species.

COI

Geographical distribution of M. sacchari and M. sorghi

The geographical distribution of M. sacchari and M. sorghi, is presented in Fig 6 and is based on 2,332 genotyped specimens. Specimens belonging to MLL-A-E-F were assigned to M. sorghi, while specimens belonging to MLL-B-C-D were assigned to M. sacchari. In the absence of SSR genotyping, specimens bearing the EF1-α haplotype H1 or the COI haplotypes H2 or H3 were assigned to M. sacchari while specimens bearing EF1-α haplotype H2 were assigned to M. sorghi (Table 4). The remaining specimens (n = 24) were considered unidentified and therefore not taken into account in Fig 6 . Only one unambiguous (i.e. allowing the species identification) data point from India was available from a public database (see Fig 6 legend). Most COI sequences stored in these public databases are either haplotypes that Supporting information S1 Fig. Melanaphis


Gaia Savino Longo

Micca Longo

Khaled Hassouni

Armelle Michau

Swaminathan Prasanna

Gaia Micca Longo

Stochastic models of systems for Nanotechnology: from micro to macro scale

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

Introduction

Computer models based on the theory of stochastic processes allow to study different systems of interest for Nanotechnology in a quantitative way, and to analyze: the microscopic scale, at the level of single atoms or molecules or ions, taking the quantum effects into account; the mesoscopic scale of objects with dimensions of the order of tens or hundreds of microns, avoiding the detail of the molecular scale; the macroscopic scale, i.e. entire reactors in the treatment of surfaces with special properties and the production of nanostructures on a large scale.

Stochastic systems can be studied by solving one of a limited class of evolution equations for probability.

If w nm are the transitions probabilities between the states n and m, the evolution of most system is governed by a set of ordinary differential equation collectively called Master Equation:

(1)

where Pn(t) is the probability distribution function, such that at time t the system is in the state n. When the states span on a continuum, n is replaced by dynamic variables like r. Such an equation can describe, for example, quantized charge state of a nanoparticle.

As variations of such a continuum formalism, we have the Boltzmann equation, which describe the interaction of classical particles via the related cross sections, and the equations for ray propagations in radiative transfer. When the state variable undergoes small variations in a single microscopic interaction, a partial differential equation, the Fokker-Planck equation, can be applied:

(2)

The diffusion coefficient D accounts for thermal fluctuations, while the vector quantity a represents an advective flow. Such an equation can be applied to the diffusion of particles performing a Brownian motion; with minor variation, it can be applied to calculate the quantum state of a microscopic system. These equations can be solved by techniques based on the so-called realization of the implied stochastic process: this means simulating an ensemble of mathematical dots distributed as P, and simulating events on any single dot in order to reproduce the transition probabilities w in the case of equation (1), and the transport coefficients a and D in equation (2). In the case of equation (1), dots (like particles with a given velocity) are randomly displaced with probability wdt in a short time dt. In the case of equation ( 2), the new random position of any dot is determined based on the short-time propagator (Green function):

(3)

The most important advantage of the stochastic approach is its universality. Either a random jump process or a continuum one, diffusion process can be considered as the basic feature of many problems in molecular structure, gas/solid interaction, transport phenomena. In many cases, the random processes are realized in the computer code by using random numbers, thereby producing one of an endless varieties of Monte Carlo models. These models can describe the quantum wavefunction of a molecule, the interaction of a nanoparticle with its environment, the motion of atoms and ions in the macroscale.

For many years, the scientific activity of our collaboration has involved the development of calculation programs based on stochastic process theory concepts, in some cases including the development of innovative algorithms. Taking advantage of the positive properties of the stochastic simulation, the application of the simulation programs developed to nanotechnology includes both the microscopic and the macroscopic scale. On a microscopic scale, we are able to study the electronic structure of small molecular systems, in particular connected with hydrogen plasma technology. On a macroscopic scale, we are able to perform simulations of plasma reactors and systems in molecular regime for the production and treatment of nanomaterials at low temperature.

On the Italian side, these researches were carried out within the CNR Institute of Nanotechnology in the course of several years of activity, although they are presently continued in association with the newly created Institute of Plasma Science and Technology.

In this short highlight, we will summarize some of these applications by choosing in particular recent case studies, trying to highlight the characteristics of these research themes that have led to their success. Some new calculations are produced to illustrate the kind of results which can be obtained. We will proceed from a decidedly nanoscopic scale of molecules and ions with a small number of electrons confined in cavities of various kinds, then we will analyze the study of nanometer-sized dust grains that exhibit special properties from a thermal and electrical point of view, and we will arrive at the simulation of systems on the laboratory scale and used for the production of nanomaterials in both fundamental and applicative contexts.

Quantum Monte Carlo models of nano-confined hydrogen species

Quantum Monte Carlo (QMC) methods are becoming more and more interesting in the determination of the quantum state of nano-systems, as they allow to minimize the use of approximations and analytical methodologies and do not need heavy adaptations concerning geometry changes. For this reason, they are very useful in the study of systems confined in potential wells of nanometer size and of different shape. Such simple confined physical systems have indeed a very long tradition in physics, not only as basic quantum systems, but also as models for extreme high-pressure states of matter, spectroscopically active defects in solid lattices and chemical species in molecular cages [START_REF] Sabin | The Theory of Confined Quantum Systems, Part I and II[END_REF]Sen, 2014]. In this regard, the confined states of molecular and ionic states of hydrogen, like H 2 + and H 2 have been much studied, since they allow much accurate calculations [Colín-Rodríguez, 2011]. At the same time, these systems have numerous real applications in the fields of nanotechnology, microelectronics, energy management and quantum computing [START_REF] Harrison | Quantum Wells, Wires and Dots: Theoretical and Computational Physics of Semiconductors Nanostructures[END_REF].

Among several Quantum Monte Carlo methods available [Foulkes, 2001], the diffusion Monte Carlo (DMC) method offers the possibility of considering confining walls of any geometry while retaining the use of cartesian coordinates and to include any distribution of positive charge within the cavity. DMC makes use of the similarity between the imaginary time Schrödinger equation and a generalized diffusion equation, which, as shortly discussed in the introduction, can be solved using a stochastic calculus and simulating a random walk. Basically, the diffusion process of numerous points in the h "w k " . Th k gy f h w q with the diffusion process, while the potential energy terms formally produce reactions of destruction or multiplication of the walkers, which are consequently included in the stochastic process. The eigenvalue of the energy is adjusted in the course of the calculation until it reaches its exact value, which is reached when the average number of walkers remains constant during the stochastic process, as demonstrated by the theory. The estimated eigenvalue is then increased if the number of walkers is increasing and decreased if the number of particles is decreasing: this simple technique leads to a statistically exact estimate of the ground state energy.

The method always provides the energy of the ground state of the system, with a very limited number of numerical parameters to be optimized and no need to expand the wave function, even if statistical convergence may require significant calculation times. Excited states can be found by looking for the ground state with restrictions in the form of nodal surfaces, that can be often identified a priori thanks to symmetry considerations.

On this basis, in recent years we have developed a suite of Diffusion Monte Carlo codes designed for the study of hydrogen-based and nano species confined by infinite potential wells of various geometries [Micca Longo, 2015a;Micca Longo, 2015b;Micca Longo, 2015c;Longo, 2018]. In contrast with previous Monte Carlo calculations [Sarsa, 2012], our programs do not employ an analytic guess wavefunctions to reduce variance by the importance sampling technique. Instead, it performs very long imaginary time evolutions, while averages are collected. While this technique leads to somewhat larger oscillations and fluctuations compared to the employ of analytical guesses, it has the advantage to be applicable to different confinement geometries with minimal variations, without having to change the expression of the guess wavefunction, since no such guess is used. This technique allows to study different systems, in the ground and excited states, under different geometries of confinements, using essentially the same program.

We found that the lowest excited states of the H 2 + molecule (the  u ,  g and  u ones) in spherical confinement conditions have potential energy surfaces (PESs) similar to those of the ground state, even if, without confinement, they are repulsive states (Figure 1, a-b) [Micca Longo, 2015a;Micca Longo, 2015b]. The energy reference of any curve in figure 2b is the same and has been chosen in such a way that, asymptotically at , the energy of the unconfined system H(1s) + H + is equal to zero. The same reference is used for the other curves: this is why the dissociation limit energy in figure 2b increases when the confinement become tighter. As concerns the electronic ground state, the compression produces an increase of the curvature of the PES at its minimum due to compression of the electron cloud: accordingly, an increase in the Raman frequencies is predicted, which in principle should be detectable in experiments.

A recent result (Figure 1, c-d) is the first calculation of the PESs of the H 2 + ion confined by an octahedral potential well, which describes the quantum state of this species in the cavities of real crystal lattices much better than a spherical or ellipsoidal potential well [Longo, 2018]. In this case, the PES is quite complex since it depends not only on the internuclear distance but also on the orientation of the molecular axis into the well.

Figure 1. The H 2 + system with different nanoconfinement geometries studied with the DMC code: a) potential energy surfaces of the unconfined molecule, in its ground state and in its three lowest excited electronic states ; b) potential energy surfaces of the excited state confined into a spherical cavity (r 0 is the cavity radius, a.u. = atomic units): the confinement changes the character of the state of symmetry, which develops a pronounced minimum which is absent in the free molecular state.; c) geometry of the tetrahedral confinement of H 2 + ; d) potential energy surfaces of H 2 + for different degrees of octahedral confinement (c is the distance between the center and any vertex of the octahedral hole).

After [Micca Longo, 2015a;Micca Longo, 2015b, Longo 2018].

Another advantage of the DMC method is that, in the case of few-electron systems, quantum correlations are taken into account exactly. Although the results presented in this account refer to a one-electron system, the method can be easily applied to few-electron systems by increasing the dimensionality of the problem: we have recently studied the confinement in different geometries (including nanotubes) of the H 2 molecule and the triatomic ion H 3 + . These calculations are being finalized: preliminary results of the collinear H 3 + confined in a cylindrical well mimicking a nanotube are shown in figure 2 and 3: The universal features of the stochastic description manifest themselves here, not only because of the easy adaptability of the code to different quantum problems, but also in view of the adaptability of a code based on the diffusion equation to completely different problems. For example, DMC was applied to the description of the dissociation of CO 2 molecules in the context of greenhouse effect mitigation [Diomede, 2017]. Despite being developed for a classical problem of chemical kinetics, the last code is essentially the same as a DMC code for solving the quantum problem: same mathematical foundations (but different interpretations of symbols), code structure, management of random events. This universality is a very concrete advantage as to multidisciplinary code development.

Models of the kinetic and thermal behavior of nanoparticles in low temperature plasma

In plasma reactors, nanoparticles of various compositions can be produced as a result of chemical reactions. Arguably, the most important mechanism governing the dynamics of the nanoparticles in dusty plasma is particle charging. The nanoparticles inside plasma acquire a charge due to frequent collisions with the charged species, ions and electrons. In fact, most of the particle dynamics of the nanoparticles such as formation, growth, temperature and transport are strongly influenced by their charging process that is governed by the ion and electron current at the nanoparticle.

Having an analytical expression for the current at a solid nanoparticle is difficult for most cases and molecular dynamics Monte Carlo simulations where the motion of charged particles is integrated using a velocity Verlet algorithm is often used to study the charging of particles. The motion of particles is well described by a Langevin equation with the effects of collisions represented by a statistical Gaussian source term [Gatti, 2008;Gopalakrishnan, 2012;[START_REF] Chahl | [END_REF] with inclusion of various other forces.

The simplest theory for the charging of solid particle immersed in plasma is the orbital motion limited theory (OML) [Allen, 1992] which is valid for most low pressure. The theory highlights that an ion of a given kinetic energy entering a limiting orbit would be eventually captured by the particle. The dimension of the limiting orbit is generally small compared to the mean free path and collisions can be assumed to be non-existent in this small region. In this case, the electron and ion currents can be derived from the straightforward expressions based on Orbit Motion Limited theory (OML) for nocollisional sheath [Allen, 1992]:

(4)

(5)

Where j e /j i , n e /n i and m e /m i are the electron/ion flux densities, and masses, respectively.  is the electrical potential of the particle with respect to the surrounding plasma, and k b is the Boltzmann constant. The heavy species, i.e., neutrals and ions, are assumed to be in thermal equilibrium at a temperature T g , while the electrons are assumed to be at a thermal non-equilibrium with usually a much higher temperature T e . Actually, the RHS of equation ( 4) represents the thermal flux of electrons that can overcome the potential energy barrier imposed by the negative charge of the particle, i. e. only those electrons with a thermal kinetic energy that is above the potential barrier can reach the particles.

The RHS of equation ( 5) represents the population of ions that reach the particles estimated on the basis that an ion dynamic in the vicinity of the particle is governed by the central coulomb potential field (no collision) and that only those ions with an initial velocity that lead to a trajectory that fall on the particle contribute to the ion flux on the particle.

At steady state, the electrical current collected by a particle immerged in a plasma is zero, i.e. the particle is at the floating potential. The particle charge and potential can be estimated by expressing the balance between electron and ion currents as inferred from equations ( 4) and ( 5) and assuming that the particles behave like spherical capacitance, i. e., with q p and r p the charge and radius of the particle.

However, the OML approximation is not valid for various conditions encountered in plasma [Goree, 1994]. For example, at moderate to high pressures, the mean free path decreases and ions and electrons may undergo frequent collisions with other neutral species inside the sheath of the particle. The frequent collisions can reduce the kinetic energy of the charged particles and can be eventually be captured by the particles. In other words, the presence of collisions can enhance the ion and electron currents at the particle surface. Furthermore, there can be other effects such as image potentials which can further alter the current of charged particles. There has been extensive work to improve the charging model incorporating different effects in the literature [Allen, 1992;Khrapak, 2006; D'y hk 2 ; G 2 8; G k h 2 2; h h 2 9].

When the charging process may be assumed continuous, the average charge carried by a solid particle is such that there is a balance between negative and positive currents at the particle. In a typical dusty plasma environment, the particles acquire a net negative charge owing to the higher electron mobility. Further, due to the difference between the ion and electron mobilities, the particle charge is subjected to fluctuations. In fact, description of charge fluctuations has been identified to critical for particle growth and morphology. [Cui, 1994] performed Monte Carlo simulation of particle charging and showed that charge fluctuations can be very significant. Further, the particle charging can be transformed into a Fokker-Planck description to obtain the mean particle charge and distribution [Matsoukas, 1995;Matsoukas, 1997]. It was shown that, for OML theory of particle charging, the particle charge distribution can be represented as a Gaussian [Matsoukas, 1995;Matsoukas1997]. However, it has to be noted that charge is discrete and modeling particle charging using a continuous model may not accurately describe the charging process [Michau, 2019]. This is especially the case when the charge fluctuation is of the same order as the average charge, which usually takes place for small solid particles, i.e., typically less than 10 nm. In this case, we consider the stochastic nature of particle charging as a Poisso ' wh h h g b f g h g states and the charge distribution is obtained by solving the master equations (equation 1) for charge fraction of each charge state q given by: (6)

Where and are the ion and electron currents respectively, and is the surface area of the particle. The solution of the previous equation requires a current law that may be driven from OML or any other theory.

The difference between the continuous model and the discrete model of particle charging is very stark for small particles as seen in the tabulations of average charge of particles using the two models. r p (nm) q ref (number of charges)

q stat (number of charges) 1 -0.25 -0.91 2.5 -0.63 -0.99 5 -1.26 -1.66 50 -12.6 -12.9

Table 1. Average charge of particles calculated using OML equilibrium and statistical approach for n e = n i . q ref is the reference particle charge calculated using the OML theory with the electron and ion current given by equations ( 4) and ( 5).

q stat is the particle charge derived from the statistical model (equation ( 6)). Figure 4 shows the charge distribution for nanoparticles of different sizes obtained using Monte Carlo simulations inside a plasma as described in [Michau, 2019]. While the Gaussian approximation for charge distribution was appropriate for large particles, such an approximation is no more correct for very small particles.

One of the side-effects of charge fluctuation is the stochastic thermal balance of the particle. Due to continuous charging, the average particle temperature can be significantly higher than the surrounding gases [START_REF] Swinkels | [END_REF]Arnas, 2009;Maurer, 2011]. Further it was shown in [Prasanna, 2019] that the average particle temperature calculated using a continuous theory [Maurer, 2011] would severely under predict the temperature averaged from Monte Carlo simulation (Table 2). This is mainly because the stochastic charging process also results in a stochastic heating which is not accounted by the continuous theory. This is even more significant for small particles. Table 2. Average particle temperature using continuous model and Monte-carlo simulations for the plasma conditions described in [Michau, 2019].

In order to accurately describe the stochastic nature of charging and thermal balance, a population balance equation for charge q and temperature T p was developed by [Prasanna, 2019]. The solution of the population balance equations permitted to obtain both the charge and temperature distribution functions for the particles under different plasma conditions. Figure 5 shows the normalized particle cumulative temperature distribution (NCPTD) obtained for a particle size of 1 nm radius for the plasma described in [Michau, 2019] for different electron densities. The value of the NPCTD at a given temperature T corresponds to the fraction of particles with temperature higher than T. It is remarkable that one can find substantial quantities of particles at high temperatures. As a matter of fact, 10% of the particle shows temperature-balues above 1000 K and 1% of the partuicles have a temperatire abve 2000 K for an electron density of 10 11 cm -3 . The fraction of particles having high temperature decreases with the electron densities. The electron densities of the plasma depletes under strong presence of dusts. The results are significant for plasma processing of crystalline nanoparticles. The critical temperatures for obtaining crystalline nanoparticles is a function of particle radius [Hirasawa, 2006] and the possibility of having large fractions of particles at higher temperature permits having higher yield of crystalline particles. In that sense, one needs to design a plasma process such that dilute quantities of particles nucleate in the plasma so that particles can easily reach temperatures greater than crystallization temperature [Mangolini, 2009;Kramer, 2015].

Figure 5. Cumulative particle temperature distribution functions for particle size of 1 nm for the plasma conditions described in [Michau, 2019] for ion densities fixed at 10 11 cm -3 and different electron densities corresponding to different discharge electronegativity.

The other direct consequence of charge fluctuation is growth of particles through coagulation. As most particles are negatively charged, growth of nanoparticles by coagulation is very limited. However, small particles can be neutrally charged due to charge fluctuations (as seen in Figure 5) and there can be coagulation between large negatively charged particles and small neutral or positively charged particles. The importance of charging model discussed earlier is critical in order to obtain a good estimate of the coagulation rates [Le Picard, 2016;Mamunuru, 2017]. Like particle charging, coagulation rates between two particles depends on multitude of parameters such as mean free path (free regime, transition regime or continuous regime). Further, the attractive forces between the charged particles can increase due to polarization [Huang, 1991;Mankelevich, 2002;Mankelevich, 2008]. The velocity Verlet algorithm can also be employed to determine the coagulation rates for different conditions [Gopalakrishnan, 2012].

Finally, the different statistical models described to determine the dynamics of the particles act as inputs to fluid aerosol models which can capture the big picture of the dusty plasma dynamics and determine the space-time distribution of particle density, mean particle diameter and particle size distribution. Several fluid aerosol dynamic studies [Ravi, 2009;Agarwal, 2012;Michau, 2017;Michau, 2019] have stressed the importance of charging models and its other effects discussed in the paper on having a simulation comparable with the experiments.

Models of neutral transport in plasma reactors

Low pressure reactors (0.1-100 millitor range), possibly associated with electrical ionization of the residual gas (plasma reactors), have unique characteristics in the production of nanomaterials and functionalized thin films [START_REF] Chen | Lecture notes on principles of plasma processing[END_REF]Kashtanov, 2007;Meyyappan, 2011;Bianco, 2018], with the possibility to obtain excellent uniformity and reproducibility and much reduced heating of substrate and the film, avoiding thermal reactions that can spoil the final result. Computer modeling allows to obtain useful information on the phenomena occurring in reactors at the macroscopic scale, provides insight for their optimal use avoiding a complex diagnostic, and allows to estimate the effects of geometry changes without actually performing very expensive modification of the real reactor. In this respect, the use of stochastic methodology allows to solve the most complete forms of the transport equation avoiding the use of approximations which are often non tenable in the extreme non equilibrium conditions found in these devices [Longo, 2004;Longo, 2006;Donkó, 2011a].

In the past, our research groups produced several simulation codes for to the prediction of the chemical reactions of transport phenomena and collective phenomena on a macroscopic scale that occur in plasma reactors, with both capacitive and inductive coupling to the electric generator [START_REF] Hassouni | [END_REF][START_REF] Capitelli | Plasma Sources Science and Technology[END_REF]Diomede, 2007;Diomede, 2008].

An example of application where a multiscale approach was adopted is the production of Cs/Mo bimetallic surfaces as catalyst in the Neutral Drive concept for the ITER nuclear fusion experiment [START_REF] Dudnikov | 20 years of cesium catalysis for negative ion production in gas discharges[END_REF]Diomede, 2009, Taccogna 2008]. The macroscopic models for this system, employing Monte Carlo models of the plasma phase and transport processes, were also supported by rigorous ab initio simulations by methods of quantum chemistry [Damone, 2015], in which Cs atoms were launched on a pristine 001 Mo substrate, which allowed to study the partially ordered geometry of a molybdenum surface covered by half a cesium monolayer. This appears to be the most efficient c y y f h h by gg g h h " " f coverage, leaving part of the substrate metal visible to hydrogen molecules, is beneficial to get the highest reaction yield. Our set of nuclear coordinates for Cs-covered Mo (Figure 6), which is available on request to the authors, has been used in several subsequent calculations by other groups, e.g. [Rutigliano, 2017]. At the macroscopic scale, as mentioned, the use of stochastic methodology in the study of the transport of neutral particles in the context of plasma reactors is an extremely important topic, since these chemically activated particles from the plasma (for example, in the form of radicals or excited species) can start high added value chemical reactions that take place on the substrate, while avoiding heating.

Recently, our activity has focused on the development of models for the transport of neutral particles based on the stochastic solution of the transport equation in integral form [Longo, 2006;Longo, 2009;Galtier, 2013]. The method can be defined as a variant of the Ray Tracing technique used in computer graphics: the primary source of neutral particles is assimilated to a source of radiation that produces a large number of rays (from a few thousands to several millions), which simulate both the scattering by inner reactor surfaces and by ambient gas molecules, thanks to a Monte Carlo random selection technique.

A very interesting feature of this approach is the possibility of simply describing, within a handcrafted program, real geometry of even very complex reactors, thus allowing a unique capacity for customization and experimentation, as well as the possibility of adding, with relative ease, additional elements of reactor physics and chemistry, which is often not allowed by commercial programs. In this way the program GAS was recently created in a Fortran90 version, while a Python3 version is being developed. The program uses a space description matrix (Figure 7, left) to reduce strongly the cost of the detection of boundaries. The numerical value associated with each element of this matrix corresponds to a volume occupied by gas or to a solid structure, for which it is also possible to specify the physical properties: for example, the structure can adsorb the traced particles or diffuse them with an isotropic or more complex law. The same matrix is used to sample the particles during the motion and produce a density plot (Figure 7, right). In this way, it is possible to describe reactors with any geometry: an example is shown in figure 8. In order to show the possibilities offered by this type of simulation, figure 9 illustrates a toy system that b " z ": h y h b taken seriously as a geometry in a reactor for the treatment of materials, but it is a demonstration of the achievable geometric complexity with simple instructions of this type, based on the geometric definition of the values of matrix elements.

Figure 7. Uses of the mesh in the GAS model. Left: in this example, the marked cells are used to locate the surface of a cylindrical object as part of a reactor model (the resolution is much better in actual calculations). The figure shows that, while the centers of the mesh cells are scanned, they are compared with the analytical expression of the feature: any cell whose center is included into the feature is marked by a special value of a flag (represented by the yellow color). The simulated particles which enter into such marked cells change their status, in order to describe their interaction with the feature surface. The visible irregularities of the stored shape obviously disappear with a denser mesh. Right: storing particle data on the same matrix to produce a density plot. It is possible to simulate the transport and diffusion of particles introduced from a source into a point of the reactor in a simple way, and obtain any type of potentially measurable quantity in a reactor of this type, such as the collision frequency (hits per unit surface per second) on any of the walls and the average concentration of particles of the species introduced into the volume of the gas in each point of the gas phase.

In the past, we included in models much more simplified geometry, some refined aspects of transport physics: for example, the thermalization of high-energy radicals produced by exothermic reactions [Panarese, 2013] and the transport of atomic particles in transonic and viscous gas flows (including shocks), calculated in turn by non-linear stochastic techniques for the same reactor geometry [Longo, 2009]. These routines are now going rewritten to be included in the GAS code. Thus, we will be able soon to study systems characterized by both chemical-physical and geometric complexity at the macroscopic level.

Even if in the demonstrations here the reactor geometry is two-dimensional, the motion and scattering of the particles is three-dimensional, and there is no difficulty, neither in principle nor in practice, in simulating complex three-dimensional systems using the same method, apart the increased computational cost and some obvious issues of visual rendering.

The setting of the neutrals transport problem as a ray tracing problem also has another advantage: it is possible, and it was demonstrated in a recent publication [Longo, 2019], to take advantage of the formal analogy between the transport problem of neutral particles and Radiative Transfer (RT) theory to reformulate, at least under certain symmetry conditions, the transport problem in a particularly convenient integral form, i.e. a generalization of the Schwarzschild-Milne equation, which include terms describing the chemical reactions in the gas phase and the catalytic properties of the reactor wall. These deterministic methods, which also arise from the theory of stochastic models, allow to obtain considerable accelerations in terms of calculation time with respect to the corresponding Monte Carlo formulation and allow an easier integration of the corresponding module with deterministic multiphysical models of a full reactor.

Conclusions

As can be seen from the study cases presented in this work, the computer simulation based on the ideas of the stochastic process theory concretely allows to study in a quantitative way phenomena and objects in the realms of nanoscience, ranging from the microscopic scale of the smallest confined molecules and nanoparticles to the macroscopic scale of plasma reactors or vacuum systems, used in the production of nanostructured surfaces. It also allows to examine very different aspects of the physics of these systems at their respective scales, ranging from the response to static potentials to the dynamic response to strong non-equilibrium or time-dependent conditions. The universality of the formalism allows to develop custom and native codes, starting from algorithm and sources originally developed for very distant and different applications. Many years of experience in this direction taught us that when, as in our case, there are already considerable computer programming skills in the collaboration, the creation of native codes exploiting these universal features allows to get the highest versatility, prompt answers to current research issues, and also provides a unique first-hand experience in understanding the physical principles.

Figure 2 .

 2 Figure 2. Collinear H 3 + confined in a longitudinally infinite cylindrical potential well, mimicking a nanotube. The two internuclear distances are equal and they are represented by d. c is the confinement dimension.

Figure 3 .

 3 Figure 3. Potential energy curves of a free collinear H 3 + and of a collinear H 3 + confined in a nanotube with confinement dimension c = 2 a. u. .

Figure 4 :

 4 Figure4: Particle charge distributions obtained for different particles sizes in a plasma described in[Michau, 2020] with electron density 10 11 cm -3 and electron temperature of 0.1 eV.

Figure 6 .

 6 Figure 6. The geometry calculated for a half-monolayer of cesium on molybdenum h w g h w g f " f " g f h b ff . g f hy Fusion vol.57, no.3, March 2015 [Damone, 2015].

Figure 8 .

 8 Figure 8. Left: Ray Tracing (RT) in vacuum with scattering surfaces of different shapes. The particles are injected from center-top. The figure is a 2D slice of a 3D geometry infinitely extended vertically: therefore, the two obstacles in the upper region are prisms, the 5 obstacles are cylinders. All surfaces diffuse incident particles isotropically. Right: time-dependent transport of neutral particles, including scattering from the bulk gas, in a model reactor using the GAS code. Two prisms and a cylinder are introduced as isotropic diffusing obstacles.

Figure

  Figure 9. The " z ": g y f h GA by D'A z hy . The traced particles are injected into the bulk gas from the inlet above. Conventions are the same as in Figure 8.