Baptiste Suchéras-Marx

Gilles Escarguel

Jorge Ferreira

Øyvind Hammer

Statistical confidence intervals for relative abundances and abundance-based ratios: simple practical solutions for an old overlooked question

Micropaleontologists often consider relative abundances of taxa to infer past ecological, environmental and climate conditions and dynamics. However, most published micropaleontological studies involving relative abundance data still do not routinely consider the counting uncertainty inherent to any sample, and thus simply ignore the statistical confidence interval (CI) related to a relative abundance or abundance-based ratio value. In an attempt to make this rather classic computation freely and easily available to the scientific community, we highlight here the calculation of binomial proportion CIs based on the 'exact' Clopper-Pearson method as implemented in the user-friendly PAST freeware. We also introduce a general solution for the computation of the CI related to any abundance-based ratio. In all cases, we strongly recommend that future studies involving taxonomic abundance-based data should systematically display the CIs associated to sample estimates, the only way to integrate sampling uncertainties into result interpretation.

Introduction

Micropaleontology is a diversified field within paleontology, focusing on different fossil groups from centimeter to micrometer-size and addressing a large array of scientific questions. It is commonly used in biostratigraphy, geochemistry, paleoecology, paleoclimatology and paleoenvironmental reconstructions, as well as evolutionary studies (e.g., [START_REF] Armstrong | Microfossils, Second Edition[END_REF][START_REF] Saraswati | Micropaleontology[END_REF]. Many works consider the relative abundances of taxa (expressed as proportions or percentages, based on counted fossil specimens per sample) as the starting point from which all subsequent analyses derive. Hence, many conclusive interpretations in micropaleontological studies are dependent upon the reliability of such relative abundance results. Nannofossil micropaleontology (which focuses on the study of micrometer-sized calcareous plateletscoccolithsproduced by unicellular haptophyte algae and on calcite particles of incertae sedis organisms) provides numerous and diverse examples of relative abundance data-based studies. Such kind of data together with paleoclimatic interpretations have been abundantly produced over the last decades (e.g., [START_REF] Mcintyre | Coccoliths as paleoclimatic indicators of Pleistocene glaciation[END_REF]. This type of data is analysed from the Late Triassicwhen these organisms first appear in the geological record [START_REF] Gardin | Where and when the earliest coccolithophores?[END_REF])up to Holocene microfossil studies. Among the micropaleontological studies dealing with relative abundances, only few works display the statistical CIs associated with empirical (i.e., sampled) proportions or percentages e.g., pollens (Beaudoin et al., 2007); calcareous nannofossils [START_REF] Beaufort | Past dynamics of the Australian monsoon: precession, phase and links to the global monsoon concept[END_REF][START_REF] Suchéras-Marx | Paleoenvironmental and paleobiological origins of coccolithophorid genus Watznaueria emergence during the Late Aalenian-Early Bajocian[END_REF]; phytoliths [START_REF] Strömberg | The Neogene transition from C3 to C4 grasslands in North America: assemblage analysis of fossil phytoliths[END_REF]. This recurrent lack of CIs undermines the statistical reliability of the interpretations and conclusions, a serious issue which has been often dealt with in the literature (e.g., [START_REF] Mosimann | Statistical methods for the pollen analyst: multinomial and negative multinomial techniques[END_REF][START_REF] Maher | Statistics for microfossil concentration measurements employing samples spiked with marker grains[END_REF][START_REF] Fatela | Confidence limits of species proportions in microfossil assemblages[END_REF][START_REF] Strömberg | Methodological concerns for analysis of phytolith assemblages: Does count size matter?[END_REF][START_REF] Heslop | Diagnosing the uncertainty of taxa relative abundances derived from count data[END_REF], although so far with little practical application within the micropaleontological community. In this note, we present a simple and rather classic practical solution on this issue, using the 'exact' Clopper-Pearson confidence interval method for relative abundance. Notwithstanding some software being already available (e.g., PALYHELP [http://www.ncdc.noaa.gov/paleo/softlib/palyhelp.html] or Dorai-Raj's (2014) R Package 'binom' [https://cran.r-project.org/web/packages/binom/index.html]), such calculation is now implemented in the free and user-friendly statistical software PAST [START_REF] Hammer | PAST: paleontological statistics software package for education and data analysis[END_REF]. In addition, we introduce two solutions (in the form of an excel sheet and also implemented in PAST) for the computation of the CI related to any ratio contrasting the abundances of two or more taxa as customarily defined and used over the last decades in many paleoecological, paleoclimatic or paleoenvironmental micropaleontological studies. In order to increase the statistical reliability, we hope that such user-friendly computational solutions will encourage the micropaleontological community to provide the statistical confidence intervals associated with relative abundance and abundance-based ratio results.

Methods

The 'exact' Clopper-Pearson confidence interval for a relative abundance

Binomial proportion CI calculation in the case of asymptotic normal approximation (also called Wald CI) is the simplest way to calculate a binomial proportion CI; it was recommended by some previous authors (e.g., [START_REF] Fatela | Confidence limits of species proportions in microfossil assemblages[END_REF]. If p is the empirical (i.e., sampled) proportion and n the sample size, the 95% CI related to p is:

𝐶𝐼 = [𝑝 -1.96√( 𝑝(1-𝑝) 𝑛 ) ; 𝑝 + 1.96√( 𝑝(1-𝑝) 𝑛 )],
where 1.96 is the approximate value of the 97.5 percentile point of the standard Normal distribution. The binomial proportion CI calculated with PAST since v. 3.06 is based on [START_REF] Clopper | The use of confidence or fiducial limits illustrated in the case of the binomial[END_REF] method, which is a slightly more complex technique than Wald's normal approximation. For  = 0.05 (i.e., for a 1 - = 95% CI), x being the number of successes and n the number of trials, i.e., 𝑝 = 𝑥 𝑛 is the empirical proportion of a taxon of interest (x counted specimens) within a given sample made of n specimens, we have:

[1 + 𝑛-𝑥+1 𝑥𝐹[1- 𝛼 2 ;2𝑥,2(𝑛-𝑥+1)] ] -1 < 𝑝 < [1 + 𝑛-𝑥+1 (𝑥+1)𝐹[ 𝛼 2 ;2(𝑥+1),2(𝑛-𝑥)] ] -1
, where F[c; d 1 , d 2 ] denotes the 1c quantile of the Fisher distribution with d 1 and d 2 degrees of freedom. The Clopper-Pearson method -referred as 'exact' since there is no approximationoutperforms the normal approximation methods for three main reasons: (i) it renders much more accurate confidence intervals for small-size samples; (ii) it cannot calculate CI boundaries beyond 0% and 100%; and (iii) it yields increasingly asymmetrical results as p departs from 50% whereas Wald CI is symmetric by definition. Other proportion CI computation techniques exist (e.g., [START_REF] Brown | Interval estimation for a binomial proportion[END_REF]. For instance, the Agresti-Coull, Wilson or Jeffreys CIs can be compared online to the Wald and Clopper-Pearson CI using the EpiTools web page (http://epitools.ausvet.com.au/content.php?page=CIProportion). However, these three alternate (and computationally more complex) methods usually return CI boundaries very close to the Clopper-Pearson method, which is slightly more conservative in most casesi.e., defining a slightly larger CI [START_REF] Brown | Interval estimation for a binomial proportion[END_REF].

Calculation procedure for a proportion CI using PAST v. 3.06 and higher

The free statistical software PAST [START_REF] Hammer | PAST: paleontological statistics software package for education and data analysis[END_REF]; https://folk.uio.no/ohammer/past/) is a steadily evolving and updated software for scientific data analysis. The latest released version is PAST 3.23 (February 2019) and the first version proposing computation of Clopper-Pearson 95% CI was PAST 3.06. The procedure for calculation is very simple and straightforwardsee the 'Reference manual' provided by Ø. Hammer and updated every new version (http://folk.uio.no/ohammer/past/past3manual.pdf). Input data should be presented as following (see also Supplementary Fig. 1 and Supplementary Material 1):

-Rows: each row corresponds to a sample; -Columns: except for the last one, each column corresponds to the empirical (i.e., sampled) percentage of a given taxon (ranging from 0 to 100%); -Last column corresponds to the total number of counted specimens per sample. Even if several taxa (columns) can be recorded in a single sample ordered by taxon in a spreadsheet, Clopper-Pearson 95% CI for one or more samples must be calculated for each taxon separately by selecting only two columns: the taxon of interest and the total number of specimens counted per sample (Supplementary Fig. 1). Users have to select the module 'Multiple Proportion CIs' from the menu 'Univariate' (Supplementary Fig. 1). The output result is displayed in a plot where empirical proportions are represented by dots (connected or not) and their associated 95% CIs are represented by whiskers ('Plot' tab; Figure 1, Supplementary Fig. 2). Samples can be plotted horizontally from left to right or vertically from bottom to top ('flip axes' option). As usual with PAST, a 'Graph settings' window offers several options to customize graphic preferences (e.g., graphic window size, font, font size, symbol size, min-max values and labels of the two axes, etc.). The resulting graph can finally be saved (various formats available: SVG, PDF, JPG, TIFF…) and/or copied and pasted as a bitmap image. Alternatively, users can display results in a five-column numerical table format (including sample name, empirical percentage value, 95% CI lower and upper limits, and total number of specimens per sample) which can be directly copied and pasted in any text editor or spreadsheet program of their choice ('Numbers' tab; Supplementary Figs. 3, Supplementary Material 2).

Confidence intervals for taxonomic abundance-based ratios

Over the last two decades, numerous paleoenvironmental indices have been proposed in the micropaleontological literature, corresponding to abundance ratios of selected taxa (e.g., [START_REF] Aguado | Late Barremian-early Aptian dark facies of the Subbetic (Betic Cordillera, southern Spain): Calcareous nannofossil quantitative analyses, chemostratigraphy and palaeoceanographic reconstructions[END_REF][START_REF] Aizawa | Late Quaternary paleoceanography deduced from coccolith assemblages in a piston core recovered off the central Japan coast[END_REF][START_REF] Bornemann | Reconstruction of short-term palaeoceanographic changes during the formation of the Late Albian 'Niveau Breistroffer' black shales (Oceanic Anoxic Event 1d, SE France)[END_REF][START_REF] Flores | Ocean-surface and wind dynamics in the Atlantic Ocean off Northwest Africa during the last 140 000 years[END_REF][START_REF] Gale | Marine biodiversity through the Late Cenomanian-Early Turonian: palaeoceanographic controls and sequence stratigraphic biases[END_REF][START_REF] Herrle | Reconstructing nutricline dynamics of mid-Cretaceous oceans: evidence from calcareous nannofossils from the Niveau Paquier black shale (SE France)[END_REF], Herrle et al., 2003;[START_REF] Linnert | Calcareous nannofossils from Eastbourne (Southeastern England) and the paleoceanography of the Cenomanian-Turonian boundary interval[END_REF][START_REF] Marino | Response of calcareous nannofossil assemblages to paleoenvironmental changes through the mid-Pleistocene revolution at Site 1090 (Southern Ocean)[END_REF][START_REF] Mutterlose | High sea-surface temperatures during the early Aptian Oceanic Anoxic Event 1a in the Boreal Realm[END_REF][START_REF] Tiraboschi | Origin of rhythmic Albian black shales (Piobbico core, central Italy): Calcareous nannofossil quantitative and statistical analyses and paleoceanographic reconstructions[END_REF][START_REF] Tremolada | Late Barremian to early Aptian calcareous nannofossil paleoceanography and paleoecology from the Ocean Drilling Program Hole 641C (Galicia Margin)[END_REF][START_REF] Villa | Middle Eocene-late Oligocene climate variability: Calcareous nannofossil response at Kerguelen Plateau, Site 748[END_REF][START_REF] Watkins | Calcareous nannofossil evidence for the existence of the Gulf Stream during the late Maastrichtian[END_REF]; see Supplementary Material 3). Such indices were based on the estimated ecological preferences of different taxa of interest. Even if they correspond to different taxon associations, all these paleoenvironmental indices can be expressed in a single general formula including three categories of taxa (see Supplementary Material 3):

-One or more taxa involved only in the numerator of the ratio; -One or more taxa involved only in the denominator of the ratio; -One or more taxa involved simultaneously in the numerator and denominator of the ratio. Whatever the actual number of taxa included in each category, let a, b, and c be the numbers of counted specimens (i.e., raw numbers of counted specimen) for each of these three categories, respectively; this gives the simple general form of the ratio 𝑟 = 𝑎+𝑐 𝑏+𝑐 . In some cases, r involves all three categories of taxa (e.g., [START_REF] Herrle | Forcing mechanisms for mid-Cretaceous black shale formation: evidence from the Upper Aptian and Lower Albian of the Vocontian Basin (SE France)[END_REF][START_REF] Aguado | Late Barremian-early Aptian dark facies of the Subbetic (Betic Cordillera, southern Spain): Calcareous nannofossil quantitative analyses, chemostratigraphy and palaeoceanographic reconstructions[END_REF]; in this case, r ranges between 0 (when 𝑎 + 𝑐 = 0) and a + 1 (when b = 0 and c = 1), and is undefined (infinite) when 𝑏 + 𝑐 = 0. Alternatively, there can be no taxa shared by the numerator and denominator (i.e., c = 0), leading to 𝑟 = 𝑎 𝑏 (e.g., [START_REF] Gale | Marine biodiversity through the Late Cenomanian-Early Turonian: palaeoceanographic controls and sequence stratigraphic biases[END_REF][START_REF] Linnert | Calcareous nannofossils from Eastbourne (Southeastern England) and the paleoceanography of the Cenomanian-Turonian boundary interval[END_REF]; in this case, r ranges between 0 (when a = 0) and a (when b = 1), and is undefined (infinite) when b =0. In most paleoenvironmental indices proposed so far, no taxa are involved only in the numerator (i.e., a = 0), giving 𝑟 = 𝑐 𝑏+𝑐 (see examples in Supplementary Material 3). In this case, r ranges between 0 (when c = 0) and 1 (when b = 0), and can be expressed either as a proportion or a percentage (100). Excluding the trivial case where 𝑏 + 𝑐 = 0, making 𝑟 = ∞, and considering the non-trivial situation where the numerator is non-zero, and thus r > 0, the statistical confidence interval of r results from the propagation of the binomial errors related to each sampled taxon ), r is not a proportion and the analytical combination of such asymmetric and dependent errors is not straightforward [START_REF] Barlow | Asymmetric Statistical Errors[END_REF]). Here we propose two solutions, one based on a Monte-Carlo procedure and using normal approximation of arcsine-transformed proportions and a second one based on a bootstrapping procedure. percentiles define the lower and upper limits of the x% CI, respectively (e.g., the 2.5 and 97.5 percentiles for the 95% CI). We provide an Excel sheet designed to calculate the 90%, 95% and 99% CIs associated to r, Log 10 (r), and √𝑟 (Supplementary Material 5). This method is also available in PAST since v. 3.23. Alternatively, confidence intervals for r can be estimated by bootstrapping (e.g., [START_REF] Davison | Bootstrap methods and their application[END_REF]. Random replicates of size T are constructed by resampling from the original sample with replacement. In practice this is done by selecting random a', b' and c' with probabilities proportional to a, b and c (i.e., from a trinomial distribution, or binomial if one of a, b, c is zero). A large number (e.g., 10,000) of replicates are computed, and r' is computed for each replicate. A 95% CI is then estimated from the distribution of r'. We have found that the confidence intervals estimated with bootstrapping and with the method described above generally converge for T>50, but for smaller sample sizes the bootstrapped CI is considerably smaller (Supplementary Material 6). Bootstrapped confidence intervals for r can be computed with PAST, v. 3.23 onwards.

Beyond counting strategy

Micropaleontologists use dedicated sample preparation techniques and counting procedures adapted to their goal (e.g., biostratigraphy, paleoecology, paleoclimatology) and to the different types of fossils under study. Since counting procedures are sampling procedures, all these techniques involve a minimum number of specimens to be recorded in order to achieve reliable results (e.g., [START_REF] Fatela | Confidence limits of species proportions in microfossil assemblages[END_REF][START_REF] Haidar | Sample size for concurrent species detection in a species-rich assemblage[END_REF]. For instance, [START_REF] Fatela | Confidence limits of species proportions in microfossil assemblages[END_REF] suggested that a minimum of n = 100 counted specimens is needed to sample species representing at least 5% of the original population. This easy-to-use conclusion has been highly cited so far (189 citations based on ScienceDirect and 283 citations based on Google Scholar in January 2019), although [START_REF] Fatela | Confidence limits of species proportions in microfossil assemblages[END_REF] main recommendation was "Generally, we suggest that percent abundance given in micropaleontological studies should include the binomial error estimate" (op. cit., p. 169), which was not followed by most of the studies citing them. Indeed, the CI related to any empirical proportion p directly depends on p and n, as illustrated in Table 1 for some selected p-n couples as well as in Fig. 2 for p = 50% and p = 5%: for a given sample size (n), the CI increases as p becomes closer to 50%, whereas for a given empirical proportion (p), the CI decreases as n increasesi.e., as the sampling effort increases, the empirical proportion p more and more accurately estimates the 'real' proportion value. 95% CI (see Table 1 for detailed values). (A) p = 50%. (B) p = 5%.

For the very same reason, it is exactly the same for a taxonomic abundance-based ratio, whose CI directly depends on the empirical ratio estimate r and sample size T, as well as on the relative number of specimens shared by the numerator and denominator of r (i.e., the c/T-value) (Table 2):

-For a given sample size (T), the CI increases as r increases; -For a given empirical ratio (r), the CI decreases as T increases; -For given T and r values, the CI decreases as c/T increases. It is worth noting that the smallest CIs are expected for r values close to the unity and large c/T-values corresponding to a large percentage of counted specimens involved simultaneously in the numerator and denominator of r. This could naively argue against 𝑎 𝑏 type ratios and favor ratios where a large number of abundant taxa are shared by the numerator and denominator, leading to large c-values with respect to a and b. However, in this case, r variations from one sample to another are expected to be very small, strongly reducing the between-sample discrimination power of such ratio. From a purely statistical point of view, a compromise must therefore be sought when defining a taxonomic abundance-based ratio by contrasting taxa with opposite characteristics at the numerator and denominator (increasing the between-sample discrimination power of r), and by also including taxa with 'neutral' characteristics shared by the numerator and denominator to control the extent of the CI. Table 2: 95% confidence intervals for selected values of r (taxonomic abundance-based ratio), T (sample size = number of counted specimens participating in the computation of r), and c/T (percentage of counted specimens involved both in the numerator and denominator of r). Note that the CI converges towards [1; 1] as T and c/T increases and r becomes closer to 1.

r = 0.2 r = 0.5 r = 1 r = 2 r = 5 T c

Example: Rødryggen section from Pauly et al. (2012) published in Marine Micropaleontology

In order to highlight the benefits of CI computation, we have calculated the Clopper-Pearson 95% CI corresponding to [START_REF] Pauly | Early Cretaceous palaeoceanography of the Greenland-Norwegian Seaway evidenced by calcareous nannofossils[END_REF] Rødryggen section dataset for nine coccolithophorid taxa of interest as well as for their Nutrient Index (Supplementary Material 4). These calculations are for illustrative purposes; we do not intend to discuss here [START_REF] Pauly | Early Cretaceous palaeoceanography of the Greenland-Norwegian Seaway evidenced by calcareous nannofossils[END_REF] paleoenvironmental or paleoecological conclusions. [START_REF] Pauly | Early Cretaceous palaeoceanography of the Greenland-Norwegian Seaway evidenced by calcareous nannofossils[END_REF] applied a common counting procedure by counting at least 300 specimens per sample. Figure 3 displays a plot similar to Pauly et al. (2012: fig. 3), but where the 95% CIs are added for each taxon as a gray interval, allowing for the direct observation of the statistical uncertainty related to each empirical relative abundance. For some taxa (e.g., Crucibiscutum spp.), there are indeed significant variations and long-term trends in relative abundance (i.e., percentage changes beyond the 95% CI of sample values). However, most of the variations recorded for less abundant taxa (e.g., Nannoconus spp. and Cretarhabdus spp.) fall within the 95% CI of the empirical percentages and therefore cannot be distinguished from random sampling noise of neither paleoenvironmental nor paleoecological significance. The 95% CI related to the Nutrient Index (NI; [START_REF] Pauly | Early Cretaceous palaeoceanography of the Greenland-Norwegian Seaway evidenced by calcareous nannofossils[END_REF]; Supplementary Material 3) also illustrates two contrasted situations (Fig. 3): in the Ryazanian to Lower Valanginian, a two-sample test between a single sample in the Ryazanian and a single sample in the Valanginian may remain within the range of the 95% CI estimates; on the contrary, the increase between the Lower Valanginian and the Hauterivian-Barremian falls beyond the 95% CI, and therefore appears statistically consistent.

By illustrating examples of statistically significant and non-significant changes, this case study highlights the usefulness of calculating confidence intervals and integrating them into graphs representing variations of taxonomic relative abundances or abundance-based ratios. There is not, and there cannot be, any rule-of-thumb to avoid these calculations because the statistical significance of the difference between two samples is not only a matter of absolute difference between them, but also directly depends on the compared values as well as on the sample sizes. On one hand, with n = 300, a 10% absolute difference between 10% (95% CI: [6.85%; 14.0%]) and 20% (95% CI: [15.6%; 25.0%]) is significant, whereas another between 40% (95% CI: [34.4%; 45.8%]) and 50% (95% CI: [44.2%; 55.8%]) is not; on the other hand, a 15% absolute difference between 35% and 50% is significant for n = 200 (95% CI: [28.4%; 42.1%] vs. [42.9%; 57.1%]) whereas it is not for n = 100 (95% CI: [25.7%; 45.2%] vs. [39.8%; 60.2%]). Therefore CIs should always be made available by authors as they are the only actual way for both authors and readers to critically assess the statistical reliability of the analyzed data and results, as well as the proposed interpretations.

Conclusions

This study highlights the importance of calculating and graphically representing statistical confidence intervals related to taxonomic relative abundances and abundance-based ratios as customarily studied by micropaleontologists since several decades. For this purpose, easy-to-use computational solutions now exist, including the free and user-friendly PAST software which offers a module for calculating relative abundance CIs for one or more samples. In addition, we provide here an excel sheet and added in PAST two taxonomic abundance-based ratio CI calculation solutions. Since a statistical uncertainty necessarily results from any sampling procedure, a confidence interval should always be associated with any empirical proportion or ratio estimate to allow authors and readers to objectively describe, compare, interpret and discuss the meaning and implications of the sample data at hand. We hope that the scientific benefits combined with the easy access to practical and user-friendly computational solutions such as PAST will help and encourage micropaleontologists, reviewers and publishers to produce, require and ultimately publish statistical confidence intervals in future studies.

Figure 1 :

 1 Figure 1: Representation of empirical (i.e., sampled) proportions (%) and associated binomial (Clopper-Pearson) 95% CIs in PAST v. 3.06 and higher (data: Biscutum constans from the Rødryggen section dataset published by Pauly et al., 2012).

  , and thus related to a, b and c with respect to 𝑇 = 𝑎 + 𝑏 + 𝑐. In the particular (but most frequent) case where = 𝑐 𝑏+𝑐 , its CI is actually the binomial CI related to 𝑝 = 𝑐 𝑛 where 𝑛 = 𝑏 + 𝑐, making it possible to use the Clopper-Pearson method as detailed above and implemented in PAST v. 3.06 and higher. In all other cases (including 𝑟

For

  Let e and g be the arcsine-transformed values related to the empirical proportions a/T and c/T, respectively, i.e.: 𝑒 = 𝑠𝑖𝑛 -1 By definition, e and g are normally distributed variables with a sample standard deviation (s) of √1 (4𝑇) ⁄[START_REF] Sokal | Biometry (4 th Edition)[END_REF]. Let e* and g* be normally distributed random variates with mean e and g, and standard deviation s, respectively. Then, e* and g* are back-transformed to proportion values of times (say, 10,000 times), ultimately leading to a Monte-Carlo distribution of the statistical uncertainty related to r, from which the

Figure 2 :

 2 Figure 2: Effect of the empirical proportion p and sample size n values on the binomial (Clopper-Pearson) 95% CI (see Table 1 for detailed values). (A) p = 50%. (B) p = 5%.

Figure 3 :

 3 Figure 3: Pauly et al.'s (2012) fig. 3 with the 95% CIs for nine coccolithophorid taxa of interest and for the Nutrient Index (gray envelopes).

Table 1 :

 1 Clopper

	n \ p	1%	5%	25%	50%	75%	95%	99%
	10	[~0%; 30.9%]	[0.25%;
			44.5%]

-Pearson 95% confidence intervals for selected p (sample proportion) and n (sample size = total number of counted specimens) values. Note that the CI is symmetric around p only for p = 50%, and becomes more and more asymmetric as p approaches 0% or 100%.

Acknowledgments

BSM thanks Fabienne Giraud and Emanuela Mattioli for support and discussions about the importance of confidence interval in micropaleontology; Luc Beaufort for discussions and references on the topic, Pierre Henry for his help, and Martin Tetard for comments on a previous version of the manuscript. This manuscript is a contribution to the team 'Climat' at CEREGE (BSM).

Supplementary Materials

Suppl. Fig. 1 to 3: How to calculate binomial confidence intervals using PAST v. 3.06 and higher. Suppl. Material 1: PAST input data example (.txt). Suppl. Material 2: PAST output data example (.txt). Suppl. Material 3: Examples of Temperature and Fertility/Nutrient/Productivity indices published over the last two decades and their corresponding taxonomic abundance-based ratio type (.pdf). Suppl. Material 4: Binomial Clopper-Pearson CIs calculated for the Rødryggen section dataset published by [START_REF] Pauly | Early Cretaceous palaeoceanography of the Greenland-Norwegian Seaway evidenced by calcareous nannofossils[END_REF] (.csv). Suppl. Material 5: Excel sheet for the calculation of confidence intervals (90%, 95% and 99% CIs) related to any taxonomic abundance-based ratio using Monte Carlo (.xlsx). Suppl. Material 6: Comparison of Monte Carlo and bootstrapped confidence intervals for r, with a=b=c (r = 1) and for increasing sample sizes T = 3a (.csv).

Supplementary figures

Supplementary Fig. 1


J Moritz

E Faudot

S Devaux

S Heuraux

The plasma-wall transition layers in the presence of collisions with a magnetic field parallel to the wall

Keywords: numbers: 52.40.Kh, 52.65.Rr, 52.40.Fd Plasma sheath, PIC simulations, Magnetic confinement a)

The plasma-wall transition is studied by mean of a particle-in-cell (PIC) simulations in the configuration of a parallel to the wall magnetic field (B), with collisions between charged particles vs. neutral atoms taken into account. The investigated system consists in a plasma bounded by two absorbing walls separated by 200 electron Debye lengths (λ d ). The strength of the magnetic field is chosen such as the ratio λ d /r l , with r l the electron Larmor radius, is smaller or larger than the unity. Collisions are modelled with a simple operator that reorients randomly ion or electron velocity, keeping constant the total kinetic energy of both the neutral atom (target) and the incident charged particle. The PIC simulations show that the plasma-wall transition consists in a quasi-neutral region (pre-sheath), from the center of the plasma towards the walls, where the electric potential or electric field profiles are well described by an ambipolar diffusion model, and in a second region at the vicinity of the walls, called the sheath, where the quasi-neutrality breaks down. In this peculiar geometry of B and for a certain range of the mean-free-path, the sheath is found to be composed by two charged layers, a first, positive, close to the walls, and a second one, negative, towards the plasma and before the neutral pre-sheath. Depending on the amplitude of B, the spatial variation of the electric potential can be non-monotonic and presents a maximum within the sheath region. More generally, the sheath extent as well as the potential drop within the sheath and the pre-sheath are studied with respect to B, the mean-free-path and the ion and electron temperature.

I. INTRODUCTION

Sheaths are encountered when a plasma is in contact with a surface, that can be an electrode, the wall of a laboratory reactor, or the wall of a spacecraft. Sheaths are one of the most well-known features of plasma physics, studied for decades, not only because of their theoretical interest, but because they are critical for many industrial and technological applications. In low temperature plasma for instance, they are of great importance in the context of directional etching, a key technology that has allowed the miniaturisation of semiconductor devices and the recent progresses of nano-technologies [START_REF] Chabert | Physics of radio-frequency plasmas[END_REF] . In fusion plasma physics, sheaths are also investigated because of their role in the heating and erosion issues of the plasma facing components.

Plasma sheaths forms because of the large difference between ions and electrons thermal velocities. Electron's one being usually larger, during a transient regime, the surface contacting the plasma becomes negatively charged and a depletion in the electron density appears at its vicinity. Consequently, a space charge builds up in front of the wall, in a region called "sheath", whose width is of several Debye length (λ d ). A large electric field, due to the space charge, accelerates ions towards the wall and slows down electrons. In the steady state, this space charge field allows the equality of both ion and electron current Γ i and Γ e at the wall, such as Γ i = Γ i = Γ w . It has also been demonstrated that the stability of such a charged layer depends on a criterion on the ion velocity: ions have to enter the sheath with a supersonic velocity (the Bohm criterion [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields[END_REF][3][4] ).

In the region connecting the sheath to the plasma, where the density can drop by a factor 2 or more [START_REF] Lieberman | Principles of plasma discharges and material processing[END_REF] , the quasi-neutrality does not breakdown (n i n e with n i and n e the ion and electron density respectively). In this region, called the pre-sheath, an ambipolar electric field maintains the current equality Γ i = Γ e at any position as long as the velocity towards the wall is subsonic (otherwise quasi-neutrality breaks down and the particles enter the sheath). Models calculating the potential drop in the pre-sheath has been derived, in the limiting cases of collisionless and high density plasma, taking into account ionization (a source term) [START_REF] Tonks | [END_REF]7 . Note that in this paper, and in its state of the art, we assume floating conditions for the walls, thus ambipolarity, but it is not mandatory for the existence of sheaths and a similar analysis could be carried out with non floating conditions.

In the presence of a magnetic field B tilted by an angle ψ with respect to the wall, another region appears between the Debye sheath and the pre-sheath. Ions are first accelerated along the magnetic line in the pre-sheath and then enter the new layer called magnetic pre-sheath or Chodura layer 8 . There, their velocity is reoriented towards the wall in order to satisfy the Bohm criterion at the Debye sheath entrance. The role of the magnetic field angle and strength on the different regions have been investigated by several authors since the original work of Chodura [9][10][11][12] .

When B is not connected to the collecting surface (ψ = 0), without collision or turbulence, the particle flux at the wall should vanish. Some authors have studied this limiting case of perfect confinement and given scaling variation of the sheath size s with the ion Larmor radius R l and showed that the potential drop within the sheath ∆φ s is positive (∆φ s = φ wφ s where φ w and φ s are the electric potential at the wall and at the sheath entrance respectively) and of the order of T i /e with T i the ion temperature expressed in energy units and e the electron charge [13][14][15] . In a previous study, we derived a simple law connecting both s and ∆φ s with the plasma parameters (both ion and electron temperature T i and T e , the plasma density n 0 and the ion cyclotron to plasma pulsation ratio ω ci /ω pi ), that has been checked against PIC simulations 16 . In this specific case of perfect alignment of the magnetic field with the wall, the space charge at the vicinity of the collecting surface is negative because the ion Larmor radius is larger than the electron one ; it leads to a positive electric field that pushes back ions into the plasma while accelerating electrons towards the wall.

This geometry eventually leads to a quasi-static situation where the flux at the wall in the steady state vanishes; it can be restored by turbulence or collisions.

Theilhaber and Birdsall showed in the context of bi-dimensional PIC simulations that Kelvin-Helmholtz instability due to the velocity shear within the sheath allows a particles drift across the magnetic field lines. They evidenced that the potential drop within the sheath keeps a positive value and that the sheath extent reaches 5R l [17][18][19] . When collisions are taken into account, both ions and electrons diffuse through the field line at a different rate. For high collisional plasma, when the ordering λ d < λ e ≤ r l < R l (λ e and r l are the electron mean-free-path and Larmor radius respectively) is respected, both ions and electrons are demagnetized by collisions: this almost leads to the classical situation where B = 0, ie. a negative potential drop within the sheath as well as in the pre-sheath 20 .

However, this scheme does not describe all possible situations that can arise from collision events and from the different orderings of the characteristic lengths.

It exists a mean-free-path threshold above which ions migrate faster across the field line than electrons, assuming that they are always magnetized (λ e ≥ r l ). In such a case, focussing on charged particle vs. neutral atoms collisions only, the ion motion towards the wall is ruled by a random walk whose step can be either the ion mean-free-path λ i (λ i does not necessarily equal to λ e ) or the ion Larmor radius R l , whereas the electron random step is always r l (see Fig. 1 for a geometric description of the problem). When λ i ≤ R l , neglecting the electric field effect, the characteristic time needed by a thermal ion to travel a distance d

is (d/λ i ) 2 × λ i /v ti while it is of about (d/r l ) 2
× λ e /v te for a thermal electron, with v ti and v te the ion and electron thermal velocity respectively. We can derive the following expression :

λ e λ i ≥ v te v ti r 2 l , (1)
that gives the mean-free-path threshold above which ions are faster than electrons. If ions are magnetized (λ i ≥ R l ), their characteristic time is (d/R l ) 2 ×λ i /v ti and the mean-free-path threshold becomes:

λ e λ i ≥ v te v ti r l R l 2 . (2)
Assuming both mean-free-paths of the same order of magnitude (λ e λ i ), Eq. 1 becomes:

λ e r l ≥ T e m i T i m e 0.25 , (3)
while Eq. 2 is always true as long as T e T i for λ e λ i . Then, above the critical meanfree-path given by Eq. 3, ions are expected to drift faster than electrons towards the wall.

However, at a characteristic distance of r l from the collecting surface, although the condition given by Eq. 3 is satisfied in the plasma, collisionless electrons travel faster than ions across the magnetic field line (see Fig. 1). This situation can eventuality lead to a double charged layer close to the wall, from a positive to a negative polarity towards the plasma.

In this paper we investigate thanks to PIC simulations the transition expected between the two different regimes, where ions are faster than electron, and reciprocally, with respect to the mean-free-path. A charged particle vs. neutral atoms collision operator has been developed and integrated to the PIC code allowing for the determination of the electric potential and field, as well as for the charged layer(s) at the vicinity of the walls, with respect to the mean-free-path. In the first part of the paper, the complete details of the simulations and the studied system are given; the results are presented in the second part.

We focus on the evolution of both potential polarity in the sheath and the pre-sheath, that can be opposite under a specific ordering of the characteristic lengths. An ambipolar diffusion model is also used to predict the potential drop amplitude in the pre-sheath and to determine the sheath extent in an original manner. Finally, our simulations and results give an overview of the transition between the two regimes, both independently studied in the literature for extreme values of the mean-free-path. It is important to notice also that the Debye sheath, defined as the sonic point location, where the plasma solution n i n e breaks away, is expected to disappear below a critical angle ψ * calculated by Stangeby 21 , which is of (m e /m i ) 0.5 . Indeed, according to Stangeby hypothesis, for such a shallow angle, the Debye sheath is no longer required to achieve an ambipolar plasma flow to the wall, the quasi-neutrality holds everywhere in the plasma, and the flow perpendicular-to-the wall is sub-sonic. In our specific geometry of the magnetic field, the Larmor radii define the area where the particles, ions and electrons, intercept the surface. In such a condition, a kinetic treatment, instead of a fluid one, is necessary to describe precisely the mechanisms arising in the sheath, extending over some Larmor radii. As it will be shown in this paper, thanks to the PIC simulations, the Bohm condition may not be verified, and the quasi-neutrality still breaks away. Then we will only refer to the sheath as the non-neutral area in front of the wall (positive, negative or both), that has to be distinguished to the "standard" Debye sheath whose edge is located at the sonic point. The different potential drops derived in this paper, in the sheath and in the pre-sheath, are then dependent on the method used to determine the sheath edge.

II. PIC SIMULATIONS

The simulations have been carried out using a 1D3V PIC code developed in the laboratory by the authors 16 . The studied system is a one-dimensional plasma bounded by two grounded absorbing walls separated of L = 200λ d . The geometry of the studied plasma as well as the orientation of the magnetic field parallel to the wall are shown in Fig. 1.

The simulation cell size is chosen such as 0.1 × min(r L , λ d ); such a small cell size allows an accurate description of both electron and ion motions in the range of magnetic field investigated. Initially, the superparticles are uniformly distributed on the grid (∼ 10 to 100 per cell) and their velocity randomly chosen from a Maxwellian distribution whose temperature are set by the user. Both T i and T e are chosen within the range 1 to 4 eV in

y x λ e E s , φ s E = 0 Electron Ion Collision x = L/2 x = s r l 1 4 λ e < d i < R l B Sheath n i = n e n i ≃ n e E E A λ i FIG.
1. Sketch of half of the studied plasma (left wall only). Schematic ion and electron trajectories are shown in the figure, assuming some collision events. Close to the wall, electrons do not suffer collision over a distance r l , while ions are collisionless over a distance d i varying between the ion mean-free-path (λ e /4) and the ion Larmor radius (R l ).

the different simulation runs, which is representative of the electron temperature in plasma discharges [START_REF] Lieberman | Principles of plasma discharges and material processing[END_REF] .

At each time step, the particles are moved under the action of both magnetic and electric field. While the former is set and kept constant in space and time, the latter depends on the ion and electron density and has to be calculated solving the Poisson equation by the finite difference method. The electric field on the particles position is then interpolated at the first order by using its value on the closest grid nodes. The Boris pusher method [START_REF] Boris | Proceedings of the Fourth Conference on Numerical Simulation of Plasmas[END_REF] is used next to update the position of all the particles and then time can be incremented. In order to reach a steady state, the number of ions in the simulation is kept constant ; at each time step the same exact amount of ions expelled out the plasma via both walls during the previous iteration are injected into the system at a random position. In order to preserve the neutrality of the plasma, an electron is also injected at the same random position. The velocity of the source particles is randomly set in the original Maxwellian distribution for both ions and electrons.

At the thermodynamic equilibrium, the plasma density and velocity distribution function are in a steady state, which means ionization and excitation processes are perfectly balanced, so that they can be replaced by simple elastic collisions. We are focussing here on elastic collisions with neutrals, as expected in plasma discharges. They are taken into account by using a simple hard sphere model, where the cross-section σ is independent of the velocity of the charged species and is only related to the spherical section of both incident and target particles. Within this description, assuming ions and neutrals of identical diameter, the mean-free-path defined as λ e = 1 σnn , where n n is the neutral density, is four times smaller for ions than for electrons. We express in our study the electron mean-free-path λ e in terms of electron larmor radius number in order to be able to compare the different characteristic lengths more easily and we assume the ion mean-free-path λ i = 1/4 × λ e . The collision operator in the code picks randomly within a Maxwellian distribution the neutral velocity components V nx , V ny and V nz . It is well known, because of the conversation of the total momentum, that the velocity of the center of mass -→ V g is conserved:

- → V g = -→ V * g = m n -→ V n + m p - → V p m n + m p , (4)
where the * denotes the different quantities after the collision, -→ V p is the velocity of the incident particle (ion or electron), m p its mass, and m n the mass of the neutral atom. It is

possible to express - → V p (resp. -→ V * p ) as a function of - →
V g and of the relative velocity between the incident and the target particles

- → V r = - → V p - -→ V n (resp. -→ V * r ): - → V p = - → V g + m n m n + m p - → V r . (5)
One can show, using the previous equation and the similar one connecting

-→ V n to - → V r ,
that the conversation of the total energy leads to a conversation of the magnitude of the

relative velocity | - → V r | = | -→ V * r |.
Moreover, for hard sphere collisions, all directions for -→ V * p are equiprobable as they are for -→ V * r . Picking randomly a direction for -→ V * r can be achieved by choosing both polar and azimuthal angles θ and φ such as:

cos θ = 1 -2 φ = 2π sin θ = √ 1 -cos 2 θ , (6)
with randomly chosen in [0 : 1]. The relative velocity after the collision is now:

-→ V * r = | - → V r | × (cos θ, sin θ cos φ, sin θ sin φ) (7)
Eq. 5, 6 and 7 allows the complete determination of the velocity of the charged particle after the collision -→ V * p . In case of electron vs. neutral collision, the velocity of the target -→ V n can be assumed as negligible with respect to the incident particle one [START_REF] Birdsall | [END_REF] ; in that case, the electron velocity magnitude is conserved after the collision and its direction is changed using Eq. 6. Finally a Monte Carlo sequence is used to determine which particle will collide with a neutral atom as detailed in reference 24 .

It is important to notice that such a hard sphere model for the charged particles vs neutrals collision is a crude approximation of the real collisional processes. We point out though that the aim of our study is not to give a complete physical description of these processes but to understand how the restored fluxes (collision enables particles drift across the magnetic field lines) change the potential drop in both the pre-sheath and the sheath and impact the sheath extent as well as the polarity of the space charge.

As the motion of the particles is expected to be driven by collisions through the magnetic field line, the time needed by ions to cross the system can be quite large, especially when λ e >> R l . That is why, in order to run the simulations in an acceptable time, we used reduced ion-to-electron mass ratio m i /m e from 50 to 100. Finally, the approximation used in the collision model making fastest particles have a higher probability to collide and then to migrate towards the walls (the cross-section being independent of the velocity 25 ), and the injection method we used (the effect of the particle injection on the velocity distribution functions is discussed for instance in reference 26,[START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF] ) can induce a relative cooling of the plasma with respect to the nominal loaded one. This cooling may occur together with a distortion of the velocity distribution function as shown in Fig. 7 of ref. [START_REF] Procassini | [END_REF] . It can be avoided by heating the plasma, reinitializing the velocity just as simulating Coulomb collisions 29 , or by changing the injected velocity distribution function as explained in ref. 26 . We did not observe in our simulations a significant difference between the ideal Gaussian distribution and the calculated one for electrons, and the difference for ions was relatively small, the root-mean-square deviation between the Gaussian estimator and the PIC calculated velocity distribution being of the order of 0.01. Then, actual ions and electrons temperatures, estimated by a numerical Gaussian fit from the PIC simulated velocity distribution functions,

have been used to normalized potential drops and velocities.

III. RESULTS

A. Space charge and electric field vs. λ e

The different diagnostics (electric potential, electric field, densities, velocity distribution functions...) are averaged on several tens of ion cyclotronic periods after the system has reached its steady state, ie. when the particles fluxes at the wall as well as the electric potential do not drift any more around their averaged values. As mentioned, the time needed to reach the steady state can be large, depending on the characteristic lengths of the system. In Fig. 2a is shown the particle flux at the wall vs. the simulation time for different mean-free-path to electron larmor radius ratios (λ e /r l ).

The particle flux at the wall Γ w decreases slowly with t, until it stabilizes around 100 ion gyroperiods, which is the waiting time before beginning the averaging of the different desired quantities. As a matter of fact, the time needed to reach the steady state increases with B:

the number of steps statistically required to travel across the studied plasma increases with (L/r l ) 2 or (L/R l ) 2 , depending on the slowest charged particle (L = 200λ d does not vary with B). The amplitude of Γ w depends on the ratio λ e /r l : it first increases up to λ e /r l = 10 and then decreases constantly; this trend will be explained in the following sections.

In Fig. 3 are depicted the space charge at the vicinity of the left wall (x = 0) as well as the electric field for 3 different λ e /r l ratios and a magnetic field such as λ d /R l = 0.1×λ d /r l = 0.03 (both electrons and ions are weakly magnetized, m i = 100m e and T i = T e = 2 eV). The situation evolves from a classical situation where the space charge is typically positive on 5 to 10λ d (highly collisional plasma with λ e /r l = 1 in Fig. 3a) to a situation where a bilayer, whose charge varies from positive (close to the wall) to negative (towards the plasma), takes place for a weakly collisional system (λ e /r l = 250 in Fig. 3e). For that collisionality, the spatial variation of the ion and electron density are shown in Fig. 2b, with a zoom close to the wall, exhibiting the charged bilayer. The electric field presents the same characteristics: no negative space charge is clearly visible in Fig. 3c, pushing back ions into the plasma.

There is obviously in this case a positive ambipolar field in the pre-sheath region where the simulation and ambipolar field plotted using Eq. 9 for the same λ e /r l ratios.

quasi-neutrality stands, while an opposite space charge field exists within the sheath.

The same evolution of the electric potential can be seen in Fig. 4a for λ e /r l = 10: the potential drop within the sheath is negative, while it is positive in the pre-sheath, leading to a non monotonic variation of φ(x). Below this threshold value, the total potential drop is negative for highly collisional plasma (as in Ref. 20 ) ; it becomes positive when the system is weakly collisional (larger λ e /r l ratio). In this situation, the particle flux at the wall is 10

times smaller than the one for intermediate value of the mean-free-path (see Fig. 2a) and would vanish for infinite mean-free-path. It would be close to the quasi-static situation that we previously studied assuming a null particle flux at the wall 16 . When the magnetic field is increased to B = 0.25, leading to strongly magnetized electrons (λ d /r l = 7.8) and an in-between situation for ions (λ d /R l = 0.78), the total potential drop presents the same variation than previously with an increasing λ e /r l ratio (see Fig. 4b), although we did not evidence a non monotonic variation of φ(x) for intermediate mean-free-path value. Indeed, electrons are always faster than ions over a characteristic distance r l from the wall, that induces, for a specific range of λ e /r l , the double charged layer. When the magnetic field increases, the positive space charge is confined closer to the collecting surface because of a less efficient screening (electron become less sensitive to the electric field as λ d < r l ). Added to the fact that r l is reduced when B increases, the positive space charge thickness decreases in absolute value (see the next section and Fig. 5 eventually).

The situation depicted in Fig. 3d where E changes from a negative polarity to a positive one towards the center of the plasma does not arise. As E does not take a null value in the sheath region, the situation of non monotonic variation of φ(x) disappears for the strongly magnetized electron case.

In order to study more accurately both potential drops within the sheath and the presheath, one has to determine the sheath size. Basically, two possibilities can be foreseen.

The first is to use the space charge calculated by the PIC code and to apply smoothing procedures for finally using an arbitrary criterion (eg. a percentage of the maximum space charge value) that would give the sheath extent. The second is to derive a suitable model for the electric field or potential in the pre-sheath and to quantify (with an arbitrary criterion as well), the spread between the PIC results and the expected variation (that will give the sheath size when the quasi-neutrality breaks down). We discuss the use of an ambipolar diffusion model applied to our simulations and our assumptions in the next section.

B. Modelling the pre-sheath

In floating wall conditions, an electric field must arise in the plasma in order to equalize electron and ion current and maintain the quasi-neutrality at any position. Close to the walls, a space charge field appears in the sheath because of the difference of velocity between the charged particles. From the center of the plasma (x = L/2) to the sheath edge x = s (resp. x = Ls), where n i n e , an ambipolar electric field allows the equality of ion and electron fluxes. Since both particles undergo collisions (L is always larger than r l or R l )

and since collisions are at the origin of the current perpendicular to the magnetic field line, the ambipolar diffusion model can be used to analyse our numerical results. Note that in this particular geometry of the magnetic field, because the magnetic lines are infinite, the particles can only diffuse across B. In most practical situations though, the diffusion in the perpendicular direction is not ambipolar, ions and electrons diffusing at their intrinsic classical rate. Quasi-neutrality is then usually maintained by electron flows along the field line (the short-circuit effect) 30 .

One can write in this specific case, for each specie, the continuity and momentum equations in the steady state. For ions, it comes (the derivatives are related to the spatial coordinate x only in the following equations):

∂Γ i ∂x = S n i m i V ix V ix = n i eE + n i eV iy B -T i n i -n i ν i m i V ix n i m i V ix V iy = -n i eV ix B -n i ν i m i V iy , (8)
where Γ i stands for the ion flux per unit area and S is a source term that has to be expressed relatively to our numerical model, -→ V i is the ion fluid velocity and ν i the ion-neutral collision frequency. The system of coordinate and the orientation of the E and B fields are shown in Fig. 1. Neglecting the inertia terms, one can extract the electric field E from Eq. 8 and from the identical set of equations for electrons. It comes:

E = C e D i n i -C i D e n e C e n i µ i + C i n e µ e , (9)
where µ = |e| νm is the mobility, D = T νm the diffusion coefficient and C = 1 + ω 2 c ν 2 , the subscript i or e stands for ion or electron with their own mass and temperature. In the ambipolar field region, where E = E A , we have n i n e n, so that:

E A = C e D i -C i D e C e µ i + C i µ e n n , (10)
Replacing E A in Eq. 8, assuming the ion and electron current such as Γ i = Γ e = Γ, it comes:

Γ = - µ i D e + µ e D i µ e C i + µ i C e n = -D ⊥ a n (11)
with D ⊥ a the ambipolar diffusion coefficient perpendicular to the magnetic field line. Now the continuity equation can be solved in the ambipolar region taking into account the geometry and the specificity of our system. As explained previously, the ions number is kept constant during the simulations: each time one reaches the walls, a couple (ion , electron) is injected randomly between 0 and L. The source term S is then |Γ w |/L with Γ w the ion (electron) current collected by both walls. Then using Eq. 8 and Eq. 11, we can write n = -Γw LD ⊥ a , that can be solved by assuming n(0) n(L) 0 as:

(x) = x|Γ w | 2D ⊥ a 1 - x L , (12)
that is coherent with the quadratic dependence of the density in the central region of the plasma shown in Fig. 2b. The ions number being constant in the simulation so is the particle density at the center ie. n(L/2) = n 0 . This steady state criterion allows the calculation of Γ w using Eq. 12 as:

|Γ w | = 8D ⊥ a n 0 L ( 13
)
The sheath size s is expected to be smaller than L, which implies that the source term within the sheath is almost zero. Hence the particle flux coming from the plasma at the sheath edge Γ s equals the particle flux leaving the plasma at the two walls Γ w . Then we have

Γ s = n s |V s | = |Γ w |/2,
with n s and V s the density and the velocity at the sheath entrance respectively, which gives using Eq. 13:

n s n 0 = 4D ⊥ a |V s |L (14)
The electric potential φ(x) in the ambipolar region can be calculated integrating Eq. 10:

φ(x) -φ(L/2) = - C e D i -C i D e C e µ i + C i µ e ln n(x) n 0 = -Λ ln n(x) n 0 ( 15
)
Using Eq. 14 and Eq. 15, the potential drop ∆φ p between the sheath edge and the center of the plasma is:

∆φ p = -Λ ln 4D ⊥ a V s L (16)
The sign of the potential drop in the pre-sheath region depends on the parameter Λ.

If Λ < 0, ∆φ s < 0 and the ambipolar field accelerates the positive charges along the magnetic field lines as in the classical ambipolar diffusion case. On the other hand, if Λ > 0, the positive electric field within the pre-sheath slows down the ions. Solving Λ = 0, ie.

C e D i -C i D e > 0, assuming ω ci /ν i << 1 and ω ce /ν e >> 1, it comes:

λ e r l > 2 T e m i T i m e 0.25 = θ c , (17)
which is the same expression that we inferred (Eq. 3) assuming a simple random walk for both ions and electrons (the factor 2 in Eq. 17 comes from the fact that the ion meanfree-path is four times smaller than the electron one in the model considered here). With the parameters used in the simulations, for an ion-to-electron mass ratio of 100 and T i = T e = 2 eV, the change in the ambipolar regime is expected to be of λe r l 6.32.

We can also verify our assumptions concerning the fact that the inertia term in Eq. 8 can be neglected with respect to the other terms by using the form of the velocity derived in Eq. 11, ie.

V x = -D ⊥ a n
n , and injecting it in Eq. 8 normalized to nm i . In the case of high collisionality where λ e << R l , it comes that D ⊥ a 2D i , the electric field term (from Eq. 10

) neE A nm i -v 2 ti Te T i n n , the pressure term -v 2 ti n
n and the collisional one 2v 2 ti n n . They have to be compared with

V x V x 4D 2 i n n n n -( n n ) 2 .
Normalizing the distance x to the mean-freepath λ e , it comes , normalizing the space coordinate x to R l with the same conclusion. It is finally important to notice that in this particular geometry of the magnetic field, collisions allow particles to drift perpendicularly to the wall over distances larger than the Larmor radius and so, participate as a matter of fact to the screening of electric fields. When λ e >> R l , even if λ e > L, collisions still occur in the plasma as the particles motion takes place parallel to the wall and the mean-free-path projection perpendicular to collecting surface is almost zero.

V x V x v 2 ti 4 n n n n -( n n ) 2 . It

C. Sheath Size and Electric Potential

The procedure to determine the sheath extents with Eq. 10 is as follow: the n /n ratio is calculated using the ion density that the PIC code has averaged over tens of ion cyclotronic periods; then the best agreement is searched between the electric field E P IC averaged by the PIC code and E A from Eq. 10, i.e the coefficient Λ has to be determined in the center of the plasma by minimizing the difference ∆ = |E P IC -E A |. In Fig. 3b, d andf are shown the electric field fits realized for different λ e /r l ratio using Eq. 10. A very good agreement is found in each case from the center of the plasma down to 10 or 20 λ d of the walls, depending on λ e /r l . The standard deviation σ ∆ is then calculated at the center of the plasma on 100 × λ d and the sheath extent determined on both sides (close to 0 and L) when ∆ > 5σ ∆ . The same procedure was applied to the electric potential profiles from the PIC simulations using Eq. 15 and gives comparable results. Another method based on the space charge n in e calculated by the PIC code has also been developed: the curve is first smoothed by using a locally weighted scatterplot smoother (LOESS) method 31 ; then the maximum h of the n in e is searched as shown in Fig. 3e and the sheath extent s arbitrary determined when n i (s)n e (s) = 0.5 × h.

λ d /R l = 0.03 λ d /R l = 0.78 λ d /R l = 0.03 λ d /R l = 0.78 Model FIG. 6.
a) Variation of the potential drop in the sheath ∆φ s with respect to the mean-free-path to the electron Larmor radius ratio λ e /r l for two different values of the magnetic field. b) Variation of the potential drop in the pre-sheath ∆φ p vs. λ e /r l . The curve labelled 'model' refers to Eq. 16.

The potential drops have been normalized to the real temperature measured during the simulations.

In Fig. 5a is shown the sheath extent s vs. λ e /r l extracted from both methods. A similar variation of s is found, but with a different amplitude. This is due to the arbitrary criteria used in the determination methods. The space charge method leads to a jump in the sheath extent between λ e /r l = 10 and 25, while the field method gives a smoother and continuous result. For instance, for λ e /r l = 10, the space charge is just positive and extends over 7.5λ d (see Fig. 3c): the space charge method gives logically a value for s that is of about 2.5λ d . Now the electric field in Fig. 3d varies from a negative to a positive polarity and keeps a significant value over 10 to 20λ d , whereas no space charge is visible in Fig. 3c at such an abscissa; the fit of the electric field gives for the sheath extent a larger value than the charge method (see Fig. 3d). When the ratio λ e /r l reaches 25, a negative charge is detected by the charge method between the plasma and the positively charged sheath (similarly to Fig. 3e), which gives a rapid increase of s. Hence the method based on the pre-sheath model gives a smoother result for the sheath extent, and is, in our opinion, a more reliable way of determining s. This method is used in the rest of this paper in Fig. 5b andc.

For all the considered values of the magnetic field, the sheath extent increases with λ e /r l , until it reaches a plateau when λ e /r l is above 10 to 50. This threshold range corresponds to the transition between demagnetized / magnetized ions due to collisions (the ion meanfree-path to the ion Larmor radius ratio is 1/40 of λ e /r l for m i /m e = 100 ). This increasing evolution of the sheath extent is due to the apparition of a second negatively charged layer between the positively charged sheath and the plasma (see Fig. 1). Moreover when the magnetic field increases, the sheath size decreases as shown in Fig. 5b. For high magnetic fields, the flux at the wall Γ w tends to vanish because of the higher confinement of the particles as seen in Fig. 7a, which in turn explains the reduction of s with B. Note that the evolution of the particle flux at the walls in the permanent regime Γ w vs. λ e /r l in Fig. 7a is proportional to D ⊥ a (see Eq. 13). Finally, the sheath extent on the plateau (large λ e /r l ratio) scales with λ d for weakly magnetized ions (ω ci /ω pi = 0.03) and with R l for magnetized ones (ω ci /ω pi = 0.78).

Using the previous value of s, the potential at the sheath edge φ(s) and consequently the potential drops within both the sheath and the pre-sheath (respectively ∆φ s and ∆φ p ) have been determined from the PIC simulations. In Fig. 6a is shown the evolution of ∆φ s with λ e /r l for two values of the magnetic field. As expected, the potential drop within the sheath vs. the mean-free-path to the electron Larmor radius ratio λ e /r l for several magnetic field strength.

The expected particle flux at the wall using the ambipolar diffusion coefficient D ⊥ a (Eq. 13) is plotted for a magnetic field such as λ d /Rl = 0.03. b) Variation of the velocity at the sheath entrance V s extracted from the PIC simulation using the sheath extent s given in Fig. 5 and normalized to the ion acoustic velocity C s . Inset : distribution function of the electron velocity on the x axis in the beginning (continuous black) and at the end (dashed blue) of the simulation for λ e /r l = 1 and λ d /Rl = 0.03. c) Spatial variation of the ion velocity close to the wall at x = 0 for 3 different collisionalities.

changes from a negative to a positive value with the mean-free-path. It is also reduced

T i = 2 eV.
in absolute value with the magnetic field amplitude in both regimes because B limits the diffusion towards the wall.

Note that for λ e /r l = 1, ions are demagnetized by collisions and electrons experience in average one collision per rotation around the magnetic field. It can be interesting to compare the potential drop within the sheath to the non-magnetized plasma case which is given by 21 :

∆φ s = T e 2e ln 2π m e m i 1 + T i T e (18)
In order to properly evaluate Eq. ( 18) with our plasma parameters, it is important to take into account the actual temperature of both ions and electrons in the simulations. Indeed, as explained previously, our injection and collisional models induce a numerical cooling of the plasma associated with a possible distortion of the velocity distributions. That is why T i

and T e may be different from the initial 2 eV. For instance, in the inset of Fig. 7b are shown the electron velocity distributions at the beginning and at the end of the PIC simulations for λ e /r l = 1. A clear cooling of the electrons is visible as the temperature decreases from 2 eV to 0.74 eV. Taking also into account the cooling of the ions (T i = 1.02 eV), Eq. ( 18)

gives a potential drop of e∆φ s /T e = -0.95. This has to be compared to e∆φ s /T e = -1.65

for λ d /R l = 0.03 and λ e /r l = 1 (see Fig. 6a). The ordering of the characteristic distances in this case is such as λ d < r l = λ e . In Fig. 5c, it can be seen that the sheath size extrapolated from the PIC results, in the same conditions of collisionality and magnetic field, is of s = 0.245 × R l = 2.45 × r l . Electrons experience then a few collisions within the sheath (λ e < s), and much more for ions because λ i = 0.25 × λ e . We infer that the difference between the expected potential drop given by Eq. 18 and our simulations arises for the most part of it from the collisions within the sheath. Indeed, it has be shown that potential drops increase with collisionality in previous studies 10,32 . The residual effects of the magnetic field onto electrons motion (λ e /r l = 1) can not be neglected neither, since the sheath (of 2.45×r l ) is still magnetized. Finally one particularity of this geometry of the magnetic field is the collisionality of the sheath: when the field line intercepts the surface, if the projection of the mean-free-path on the perpendicular-to-the wall direction is larger than the sheath extent (ie. some λ d ), the collisionless approximation holds. In our specific conditions, the particles are drifting parallel to the wall between collisions and, even if λ e,i > s, the collisionality has to be taken into account when modelling the sheath. Note that the potential drops in Fig. 6 (and in following figures) have been normalized to the electron temperature when they are negative and to the ion temperature otherwise, both temperatures being determined for each simulations independently.

In Fig. 6b is shown the potential drop in the pre-sheath for the same two values of the magnetic field than previously. The transition between the two regimes (ions faster than electrons) arises for λ e /r l = 5, which is coherent with our theoretical assumptions (Eq. 17).

For weakly magnetized ions and electrons (λ d /R l = 0.03) at λ e /r l = 1, e∆φ s /T e = -1.11

which decreases to 3.02 when electrons are strongly magnetized. The same tendency to an increase of the absolute value of the potential drop in the pre-sheath can be observed for the highest λ e /r l ratio in "the ion regime" with the magnetic field. Using Eq. 12, one can evaluate n (s)/n(s) 1/s, assuming the sheath extent s << L. With Eq. 11 and Eq. 13, and using the fact that |Γ w |/2 = |Γ s |, it comes that n (s) 4n 0 /s. Then, ∆φ p is found to scale with s as:

∆φ p -Λ ln 4s L (19)
As the sheath extent decreases with an increasing B for the whole range of λ e /r l investigated (see Fig. 5), the potential drop in the pre-sheath is also expected to increase in absolute value (Eq. 19 ), just as we evidenced it in our PIC simulations (Fig. 6b).

The spatial variation of the velocity perpendicular-to-the wall has been averaged over several ion gyroperiods during the simulations and shows a dependence on the collisionality as depicted in Fig. 7c. We have then extracted the velocity V s at the sheath entrance from the velocity profiles; it is shown in Fig. 7b and its evolution is very similar to Γ w vs. λ e /r l (Fig. 7a). Both the ambipolar flow (∝ D ⊥ a ) and the velocity are ruled by the slowest specie. When λe r l < θ c , D ⊥ a 2D i , and ions dominate the diffusion towards the walls. In this regime, an increase of λ e leads to an increase of the average velocity (note that this trend is possible because λ e < R l ), until λe r l = θ c . In this particular situation, D e /C e = D i /C i , the ambipolar flow (∝ D ⊥ a ) is maximum and the potential drop in the pre-sheath vanishes because both particle types have the same diffusion coefficient. Above the critical λ e /r l ratio,

D ⊥ a 2De
Ce and the velocity towards the wall decreases continuously with λ e /r l because the ambipolar regime is dominated by electrons (when the mean-free-path increases, electrons stick longer to their field line before undergoing a collision so the average electron velocity decreases). Also, using the extracted values of V s and Eq. 16, we found a very good agreement between the expected ∆φ p as given by the collision model and the simulated one in Fig. 6b. Finally, an equivalent ionization collision length λ iz can also be deduced from our model and compared to the mean-free-path. The considered source term in our simulations being Γ w /L, it can be compared to an equivalent ionization given by n 0 v te /λ iz . Using Eq. 13, it comes that λ iz = v te L 2 /8D ⊥ a . Taking the maximum of the ambipolar diffusion coefficient in our calculations (see Fig. 7a), for λ d /R l = 0.03 and λ e /r l = 25, the equivalent ionization collision length is of about 100 λ e .

D. Sheath properties vs. T e /T i ratio

In Fig. 8 is shown the evolution of the sheath size and of the potential drops in the sheath and the pre-sheath, ∆φ s and ∆φ p respectively, for two different initial electron temperature, T e = 1 eV and 4 eV, keeping T i = 2 eV, and a mass ratio m i /m e = 50.

The sheath size slightly increases with T e in the electron regime (λ e /r l < 6) because of the augmentation of the relative velocity between ions and electrons (electrons diffuse with a higher rate towards the walls). The opposite tendency is logically observed in the ion regime. The potential drop within the sheath follows the same evolution than the sheath extent for the same reasons, ie. in the electron regime, ∆φ s appears to increase with T e /T i and with T i /T e otherwise for larger λ e /r l ratios (see Fig. 8b). Concerning the potential drop in the quasi-neutral pre-sheath (Fig. 8c), it increases in both regimes when the ratio T i /T e decreases. The size L of the studied system depends on the electron temperature (L = 200λ d ) in this case and ∆φ p is expected to increase with L (see Eq. 16). Finally, as evidenced in Fig. 4 and 6 for another T e and mass ratio, both ∆φ s and ∆φ p can have opposite signs for a given λ e /r l , which leads to a non monotonic profile of the electric potential and opposite effects on ions velocity within the sheath (compare potential drops in Fig. 8b andc for λ e /r l = 1 to 10).

IV. CONCLUSION

In this paper, we studied by mean of PIC simulations how evolve the different layers existing in a one-dimensional magnetized plasma with respect to the collision rate, when the magnetic field is applied parallel to the walls. We showed that the space charge close to the wall is strictly positive for high collision rates. The sheath splits into two charged regions when the mean-free-path λ e increases, positively charged at the vicinity of the wall and negatively charged towards the plasma, this last sheath region expanding with λ e . This evolution of the sheath structure leads to a potential drop in the sheath that is first negative in order to accelerate ions towards the wall and that becomes positive for larger λ e . For all values of the collision rate or of the magnetic field strength, a pre-sheath exists, adjacent to the sheath, and where quasi-neutrality stands. The potential drop in the pre-sheath, which is well described by a diffusion ambipolar model, can be opposite, for a certain range of the λ e /r l ratio, to the potential drop in the sheath. This particular situation can lead to a nonmonotonic spatial variation of the electric potential within the sheath for weakly magnetized electrons (λ d < r l ). The sheath extent, ie. the total spatial extension of the non-neutral area in front of the wall, is found to decrease with the magnetic field strength B, scaling with λ d for weakly magnetized electrons and ions, and with R l for strongly magnetized electrons.

The sheath size also increases with the mean-free-path, for any value of the magnetic field, until it almost stabilizes. This behaviour is due to the apparition of the negative charged layer towards the plasma for the lowest collision rates. We have not investigated a regime of strongly magnetized ions (R l < λ d ) because of calculation time issues. However, the flux at the wall, following the variation of the ambipolar diffusion coefficient D ⊥ a , decreases dramatically with B. That would eventually lead to a quasi-static situation, where the flux at the wall vanishes, that we have studied in a previous work. Finally, this particular situation of a magnetic field strictly parallel to the wall is, of course, very restrictive and unlikely to happen with experimental devices, although it presents a theoretical appeal that deserves to be addressed. Moreover, as we pointed out in the third part of the paper, the infinite character of the magnetic lines leads to the impossibility for electrons to leave the plasma along them. Plasma neutrality is usually conserved by this electron short-circuit effect in most experiments, although this question is still being debated (see discussions and references inside 33 ). But the perfect alignment case of the magnetic field with the wall is the limiting case of the grazing incidence one, met in the scrape-off layers (SOL) of tokamaks 34 , or in linear machines such as ALINE 35 . Double sheaths structure and similar effects due to the finite size of both ion and electron Larmor radius, and evidenced in our study, may then arise in specific collision rate conditions, more specially when ions migrate faster towards the walls with collision and the associated random walk than following the tilted field line 36 .

In a grazing incidence of the magnetic field, the critical mean-free-path to electron Larmor radius ratio θ c derived in this paper may be valid (Eq. 17

λ

  FIG. 2. a) Particle flux at the wall Γ w normalized to the Bohm flux (Γ B = 0.5n 0 C s ) for different mean-free-path to electron larmor radius ratio (λ e /r l ) with respect to the simulation time. The parameters of the PIC simulations are m i /m e = 100, T i = T e = 2 eV, λ d /R l = 0.1 × λ d /r l = 0.03 (corresponding to an applied magnetic field of B = 0.01 T). b) Spatial variation of the particles density for λ d /R l = 0.03 and λ e /r l = 250. Inset: zoom close to the wall showing the double charged layer relatively to Fig. 3e.

λ e /r l = 1 λ e /r l = 1

 11 FIG. 3. a), c) and e) Space charge (n in e ) vs x/λ d at the vicinity of the left wall for λ e /r l = 1, 10 and 250. In e) is also shown the smoothed curve calculated via the LOESS method ; the dashed area correspond a confidence interval of the smoothed data. The sheath extent s is determined arbitrary at %50 of the maximum signal h. b), d) and f) Electric field extracted from the PIC

FIG. 4 .

 4 FIG. 4. Spatial variation of the electric potential φ for different values of the mean-free-path normalized to the electron Larmor radius λ e /r l . In a) λ d /R l = 0.1 × λ d /r l = 0.03 while in b) λ d /R l = 0.1 × λ d /r l = 0.78 ; the mass ratio is m i /m e = 100 and T i = T e = 2 eV.

FIG. 5 .

 5 FIG. 5. Variation of the sheath extent s normalized to the Debye length λ d in a) and b) and tothe ion Larmor radius R l in c) for different magnetic field amplitude. In a), two different methods are used to determine the size. The Field method is related to the use of Eq. 10 while the charge method is related to the use of the PIC space charge calculations.

  n

  FIG. 7. a) Variation of the particle flux at the wall Γ w normalized to the Bohm flux (Γ B = 0.5n 0 C s )

  FIG. 8. a) Sheath extent normalized to the ion Larmor radius s/R l for two different electron temperatures T e . b) and c) Potential drops in the sheath ∆φ s and in the pre-sheath ∆φ p for the same electron temperatures. The mass ratio is m i /m e = 50 and the nominal ion temperature

  ), and a change in the ambipolar regime occur. Applied to the case of the SOL of tokamaks (T e = 20 eV and T i = 40 eV), it comes, for a DT gas θ c = 2 Tem i T i me 0.25 13, and in linear machines such as ALINE with a gas of He at T e = 2 eV and T i = 0.026 eV, θ c 55, which are reasonable and accessible values.

  appears that n n is necessarily smaller than 1 because the density can not drop by n 0 over λ e and so 1 << n n < 1. The same comparison can be made in the low collisionality regime where λ e >> R l , with D ⊥

			4	n n -( n n ) 2 a	2De Ce and
	neE A nm i	v 2 ti	n n