Pierre Peterlongo

Nicolas Schnel

Nadia Pisanti

Marie-France Sagot

Vincent Lacroix

Identifying SNPs without a reference genome by comparing raw reads

Next generation sequencing (NGS) technologies are being applied to many fields of biology, notably to survey the polymorphism across individuals of a species. However, while single nucleotide polymorphisms (SNPs) are almost routinely identified in model organisms, the detection of SNPs in non model species remains very challenging due to the fact that almost all methods rely on the use of a reference genome. We address here the problem of identifying SNPs without a reference genome. For this, we propose an approach which compares two sets of raw reads. We show that a SNP corresponds to a recognisable pattern in the de Bruijn graph built from the reads, and we propose algorithms to identify these patterns, that we call mouths. We outline the potential of our method on real data. The method is tailored to short reads (typically Illumina), and works well even when the coverage is low where it reports few but highly confident SNPs. Our program, called kisSnp, can be downloaded here: http://alcovna.genouest.org/kissnp/.

Introduction

Biology in general, and genomics more particularly, witnessed a revolution in the middle 1970s with the development of rapid DNA sequencing techniques, notably the Sanger method which remained the standard approach for sequencing including whole genomes until the early years of the twenty first century. We have since then been witnessing a second revolution, various orders of magnitude bigger than the first, with the advent of the so-called "next generation sequencers" (NGS for short) which enable to obtain up to several hundred million bases in one single run at increasingly lower costs. These include (not exclusively) the 454 Life Sciences, SOLiD Applied Biosystems and Illumina technologies, each with its own characteristics in terms of read length and error rate. Such characteristics are however evolving extremely fast, faster indeed than the algorithms developed to handle the data such technologies produce.

This incredible acceleration has two implications that motivate the work presented in this paper: first it is now possible to obtain data for various individuals of a same species and thus to investigate the genetic differences among such individuals, and second, increasingly more often this will concern species for which we have no genome of reference, that is no genome already fully sequenced and assembled that could guide the investigation.

The genetic markers that will be of interest in this paper are so-called Single Nucleotide Polymorphisms (SNP for short). These correspond to a DNA sequence variation that occurs when a single nucleotide -A, T, C, or G -in a genome differs among members of a species or between paired chromosomes in an individual. There are two types of SNPs: substitutions or insertions/deletions. We focus here on the first type, that is on substitutions of single nucleotides.

Identifying SNPs in a population may have a wide range of applications that goes from assessing the polymorphism of the population, linking this polymorphism to phenotype information, or selecting SNPs as markers of subpopulations. However, while SNPs are almost routinely identified in model organisms, the detection of SNPs in non model species remains very challenging due to the fact that almost all methods to identify SNPs rely on the use of a reference genome.

Our objective is, given high-throughput read data for a pair of individuals, to identify a set of SNPs with good confidence, without having to perform an assembly of the reads with all the possible mistakes this entails, in a context where we do not have a reference sequence to help the identification.

We are aware of only two publications, dating both from 2010, that deal with the same problem [START_REF] Cannon | Assembly free comparative genomics of short-read sequence data discovers the needles in the haystack[END_REF][START_REF] Ratan | Calling SNPs without a reference genome[END_REF]. Recognisably, the major difficulty one faces is due to the presence of errors in the reads, which may be mistaken for a SNP. Additionally, the presence of inexact repeats in the genomes of the studied individuals, may further harden the task. In this paper, we restrict to the case where there is only one genomic variant for each individual (we say that the individuals are homozygous). In this context, the issue of repeats is greatly reduced.

Ratan et al. [START_REF] Ratan | Calling SNPs without a reference genome[END_REF] first filter the reads in order to remove the repetitive sequences, then create clusters of overlapping reads which they assemble using a short read assembler. The SNPs are finally identified in the micro-assembled regions using a combination of filters, based on the number of reads supporting each variant or the distance of the SNP w.r.t. the end of the contig.

Unlike Ratan et al., we chose not to use an assembler, which we think can make undesired choices as to sequence variants to remove during the assembly. Indeed, the purpose of an assembler is not to identify SNPs but to propose one reference sequence compatible with the data. Similar to Canon et al. [START_REF] Cannon | Assembly free comparative genomics of short-read sequence data discovers the needles in the haystack[END_REF], we work with raw reads, but we go further than a statistical description of the reads and propose to locally reconstruct the de Bruijn graph in order to identify SNPs. The use of a de Bruijn graph in computational biology was introduced by Pevzner et al. in 2001 [START_REF] Pevzner | An Eulerian path approach to DNA fragment assembly[END_REF] and used since then as a first step by many short read assemblers.

The key point of our method is that a SNP corresponds to a recognisable pattern in the de Bruijn graph, which we call a mouth, each lip of the mouth representing an individual variant of a same genomic locus.

Our aim is to directly find the mouths that may be reliably associated to a SNP without making use of any preliminary filters that may eliminate repeats. This is important because, although not explicitely stated, such filters seem to strongly rely on the assumption of an approximately uniform coverage of each sequence position by the reads in the available data. This assumption is usually not true. Moreover, the biases in read coverage may even vary across two sequencing experiments of the same genomic sample. This means that filters may remove sequences which in fact do not belong to repeats.

We thus present in this paper an algorithm which takes as input two sets of short reads (Illumina or AB/SOLiD) and outputs candidate SNPs (i.e. mouths in the de Bruijn graph), without performing any filtering nor using a short read assembler. This is what we call a comparative micro-assembly. This method is new as, as far as we know, no other treats data coming from distinct sequencing experiments. This approach presents the interest of taking advantage of differences in the data directly into the heart of the algorithm and not in a post-treatment step. SNPs are thus detected on raw read data instead of on pre-assembled sequences. We applied our algorithm on data simulated using MetaSim [START_REF] Richter | MetaSim -A Sequencing Simulator for Genomics and Metagenomics[END_REF], where we show under which sets of parameters the method works best. We finally apply the method to real data for Escherichia coli, for which experimentally validated SNPs are available [START_REF] Barrick | Genome evolution and adaptation in a long-term experiment with escherichia coli[END_REF], which is very rare. We show that our method successfully identifies the previously known SNPs, but also predicts new SNPs missed by the conservative method used in the original publication [START_REF] Barrick | Genome evolution and adaptation in a long-term experiment with escherichia coli[END_REF].

Preliminaries

Sequence, k-mers, prefix, suffix. A sequence is composed by zero or more symbols from an alphabet Σ containing |Σ| distinct characters. A sequence s of length n on Σ is denoted also by s[0]s [START_REF] Lacroix | coli long-term experimental evolution project site[END_REF] . . . s[n -1], where s[i] ∈ Σ for 0 ≤ i < n. The length of s is denoted by |s|. Finally, we denote by s[i, j] the substring s[i]s[i + 1] . . . s[j] of s. In this case, we say that the substring s[i, j] occurs at position i in s. We call k-mer a substring of length k. If s = uv for u, v ∈ Σ * , we say that v is a suffix of s and that u is a prefix of s.

De Bruijn graph. Each node of a de Bruijn graph stores exactly one k-mer. An edge connects a node n 0 to a node n 1 if the suffix of length k -1 of the k-mer corresponding to node n 0 is equal to the prefix of length k -1 of the k-mer corresponding to node n 1 .

A category of de novo read-assembly methods such as SOAPdenovo [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF], Euler [START_REF] Pevzner | An Eulerian path approach to DNA fragment assembly[END_REF] and Velvet [START_REF] Zerbino | Velvet: Algorithms for de novo short read assembly using de Bruijn graphs[END_REF] (to mention a few) uses the de Bruijn graph as a fundamental data structure. In a few words, reads are first divided into overlapping k-mers, then the associated de Bruijn graph is created and finally Eulerian paths are found in the graph for reconstructing the initial genomic sequence, or fragments thereof that are as large as possible (contigs).

One of the main difficulties encountered by such methods comes from the sequencing errors that generate substitutions and insertions/deletions in the data that must then be assembled. Such errors lead to loops in the de Bruijn graph which may hinder the Eulerian path detection. A first step in such algorithms consists thus in "cleaning the data" by removing suspicious reads and substituting suspect nucleotides. This cleaning step may be problematic when looking for SNPs as it may remove a significant part of them that will be mistakenly considered as sequencing errors.

Comparative micro-assembly model

Our method compares reads generated by two distinct sequencing experiments, and creates parts of the de Bruijn graph potentially linked to a SNP between these two experiments, thereby detecting such SNPs.

The main idea is that the de Bruijn graph of k-mers stemming from two sequences that contain a SNP presents a mouth shape as shown in Fig. 1. The algorithm described in Section 4 detects and constructs such graph shapes directly from the non-assembled k-mers coming from the sets of reads of two distinct sequencing experiments. It is important to notice that the algorithm does not reconstruct the full de Bruijn graph but focuses only on putative SNPs by building mouths. Taking into account the k-mer counts. Let us first consider the case where the sequencing is perfect: uniform coverage C and no sequencing errors nor repeats. In such case, all k-mers from A covering a SNP have C occurrences more than the same k-mers from B, and vice-versa. We considered this theoretical perfect coverage C = 20 in the example of Fig. 1. In practice, the coverage is not uniform and the reads contain errors. The consequence is that the k-mer count difference between experiments will not be constant along the mouth. To account for this, we introduce a parameter called δ, which is meant to capture a deviation from the exact case. Below, we describe the mouth model with counts.

Mouth model integrating k-mer counts. For any k-mer ω, let mult A (ω) (resp. mult B (ω)) be its number of occurrences in the set of reads A (resp. B). We define dif f (ω) = mult A (ω)mult B (ω). The SNP mouth model integrates the k-mer number of occurrences as follows. A special k-mer a op called the opening k-mer is chosen as reference (see Section 4). If a op is in the upper (resp. lower) lip, for any k-mer a i contained in a node of the upper (resp. lower) lip, we have dif f (a i ) = dif f (a op ) ± δ while for any k-mer b i contained in a node of the lower (resp. upper) lip, we have dif f -Check that the found mouths are coherent with the reads (detailed in Section 4.3).

(b i ) = -dif f (a op ) ± δ.

Creating the mouths

Selection of the k-mers opening the mouth. The opening k-mer a op is selected such that max (mult A (a op ), mult B (a op )) < max (mult A (a i ), mult B (a i )) for any k-mer a i = a op and dif f (a op ) = 0. In other words, a op is the k-mer having the smallest number of occurrences in either set A or B (possibly zero occurrence in one of the two sets), and is such that it occurs more in a set than in the other, possibly due to a SNP. The rationale for choosing the k-mer with the smallest count is to avoid choosing a k-mer involved in a repeat for opening the mouth. The opposite opening k-mer b op is selected such that a op and b op are distant by exactly one substitution and dif f (b op ) = -dif f (a op ) ± δ. Notice that the substitution position between the k-mers a op and b op may be anywhere. We denote by p this substitution position (p ∈ [0, k -1]). It is worth noticing that, for a given opening k-mer a op , several (at most (|Σ| -1)k) distinct k-mers b op may satisfy these conditions. They are all iteratively tested as mouth openers.

Extending and closing the mouth. Once a pair of opening k-mers (a op , b op ) is selected, a recursive procedure extends them to the right and left with other k-mers fulfilling the following conditions (also shown in Fig. 2). The k-mers a i+1 and b i+1 may extend a i and b i iff:

p ≥ 0 (the closing k-mer c k has not yet been reached), and

-a i [1, k -1] = a i+1 [0, k -2] and b i [1, k -1] = b i+1 [0, k -2]
(the new k-mers overlap on k -1 characters with their predecessors), and

-a i+1 [k -1] = b i+1 [k -1]
(the extension is done with a same character), and

-dif f (a i+1 ) = dif f (a op ) ± δ and dif f (b i+1 ) = -dif f (a op ) ± δ (the counting model is fulfilled).
Similar conditions apply to a i-1 and b i-1 for extending a i and b i on the left. The two lips of a mouth have to be closed. A mouth can be right-closed (resp. left-closed) if there exists a k-mer c k (resp. c -1) ) whose prefix (resp. suffix) of length k -1 is equal to the suffix (resp. prefix) of length k -1 of a k-1 (resp. a 0 ), by definition itself equal to the suffix (resp. prefix) of length k -1 of b k-1 (resp. b 0 ). Once the mouth is right-and left-closed, the procedure stops. Disabling the use of a same opening k-mer more than once Once a k-mer a op was used for opening a mouth, it is flagged and never used anymore either as an opening or as an extending k-mer. The underlying idea is to avoid detecting twice the same mouth. Moreover, we can safely discard this k-mer because, since it was the one with smallest count, it should not belong to another mouth.

Complexity

The time complexity can then be divided into two main parts as follows.

-Indexation: if N is the total number of distinct k-mers, indexing them into a tree can be done in O(N log N ) time using heap sort. This then provides access in time O(k) to any k-mer information. -Mouth creation: Given one opening k-mer a op , at worst k × (|Σ| -1)

k-mers b op may fulfill the opening conditions. Any k-mer may be opening. Thus at worst, O(N × k × |Σ|) mouth opening pairs must be tried. Given an opening pair of k-mers a op and b op , in the worst case, for each of the k -1 steps of the extension, |Σ| k-mers muts be tested. Access to a k-mer information being in time O(k), the time complexity for one mouth extension is thus O(k|Σ| k ). Thus, at worst, the time complexity for mouth finding is

O N × k 2 × |Σ| 2k .
Each k-mer being stored in O(k), the space complexity is O(N k).

Checking for read coherence

Two k-mers are linked in the de Bruijn graph if they overlap over k-1 characters, without checking whether the created k+1-mer indeed exists in the set of reads. This may lead to false-positive results. In order to remove those k-mers, in a post-treatment step, we check for read-coherency of the identified mouths. The upper (resp. lower) lip of a mouth is said to be read-coherent if it is 100% covered by reads from set A (resp. B) and, moreover, if in the upper (resp. lower) lip the SNP position is covered by at least two distinct reads from set A (resp. B). We keep only the mouths for which both lips are read-coherent. The rationale for restricting to mouths covered by at least 2 reads is to minimize the number of sequencing errors that are mistaken for SNPs. Indeed, it is unlikely that a sequencing error occurs at the same nucleotide for 2 distinct reads, as shown in [START_REF] Ratan | Calling SNPs without a reference genome[END_REF].

Applications to simulated read data

We developed the algorithm in a program called kisSnp, coded in Java, that was used for testing our approach. kisSnp is available for download at: http://alcovna.genouest.org/kissnp/.

To test our approach on controlled datasets, we applied the following procedure. Given an input sequence s ref , we generated a sequence s snp such that s ref and s snp differ by

|s ref | 1000
substitutions, each considered as a SNP. The average distance between two virtual SNPs is then 1000 nucleotides, in agreement with [START_REF] Cooper | An estimate of unique DNA sequence heterozygosity in the human genome[END_REF]. The substitutions are randomly distributed over s snp , and each substitution is introduced following the transition/transversion probabilistic model [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF].

The sequence s ref and s snp are then virtually sequenced into a set of reads r ref and r snp using the MetaSim [START_REF] Richter | MetaSim -A Sequencing Simulator for Genomics and Metagenomics[END_REF] program. Among other parameters, MetaSim enables to tune the sequencing errors model as well as the average coverage. In all our experiments, we generated reads of length 62, in agreement with the Illumina technology. kisSnp is then tested using sets r ref and r snp .

Finding the SNPs

Human chrX portion We extensively tested the parameters of kisSnp on a small portion of the human chromosome X of length 137897 bp, corresponding to a kinesin family member 4A (KIF4A). We applied MetaSim to the pair s ref , s snp distant by 137 simulated SNPs with i) no sequencing errors (Fig. 3(a)) and ii) errors following an empirical Illumina error model (Jean Marc Aury, Genoscope, personal communication -Fig. 3(b) and(c)).

One may start by observing that the quality of the results is relatively robust to the choice of parameters. With a good coverage (> 4x), large distinct sets of parameters thus enable to find almost all SNPs with no false positives. The main lessons learnt about such parameters from these results are the following:

-The δ value has not a strong influence for δ ≥ 20 (Fig. 3(b) vs Fig. 3(c)).

However, for smaller values of δ (data not shown), the specificity decreases rapidly due to sequencing errors and a non uniform coverage. -As expected, for small values of k (≤ 20), false positives are found. For larger values of k (say, ≥ 30), less SNPs are found as more k-mers involve sequencing errors and/or more positions are not covered by any k-mer.

Concerning the data, coverage is of major importance as a small coverage leads to lower sensitivity (in each case, the lower the coverage, the less sensitive is kisSnp). Illumina sequencing errors have a small influence on the results (Fig. 3(a) vs Fig. 3(b)). One important message is that, using k = 25 and δ ≥ 20, all experiments obtained 100% specificity (no false positives) for various values of sensitivity, the latter depending in particular on the coverage.

Neisseria meningitidis strain MC58 One main limitation of our mouth model could be the presence of repeats in the genomes considered. To measure the effect of repeats on the performance of kisSnp, we performed experiments on the bacterium Neisseria meningitidis (strain M C58) of length 2.27 Mbp. The size of the repeated elements in this genome range from 10 bases to more than 2000 bases, and their number may reach more than 200 copies. We performed tests on the original MC58 sequence introducing 2272 SNPs and simulating reads using MetaSim with an Illumina error model. Moreover, we performed exactly the same tests on a randomised sequence obtained from the MC58 sequence (same length and nucleotide frequencies, distribution of nucleotides following a Bernoulli model). Results for both MC58 and the randomised MC58 are presented in Fig. 4.

One may observe that even on a difficult dataset like MC58, a large part of the SNPs are identified (26%, 86% and 97% respectively with coverage 4x, 10x, and 20x with k = 20). Another important remark is that the difference between the results obtained on MC58 and the randomised MC58 is small, showing that the algorithm is robust to repeats.

Execution time

The kisSnp program, coded in Java, is a prototype and is not yet time optimised. The performance results below enable only to give an idea of the evolution of the running time with different parameters. All tests were done on a DELL laptop, quad-core 2.67GHz with 4Gb memory running under Fedora Core 12.

We started by testing on the human KIF4A dataset (simulating reads with an Illumina error model), the influence of the δ parameter on the execution time. We observed (data not shown) that this has no influence whatsoever for δ ≥ 20. On the same dataset, we fixed δ = 30 and checked the execution time for values of k varying from 0 to 40. The results, presented in Fig. 5, show that two phases may be distinguished. For k from 0 to 6, we observe an exponential time growth that is in agreement with the theoretical complexity. During this phase, the worst-case behaviour is reached, each mouth extension tests |Σ| k lips. The second phase is observed for bigger values of k, when a large number of possible k-mers are no longer present in the data, thus limiting the number of tested lips extensions. This greatly reduces the execution time, which starts decreasing for k ≥ 25 as less and less mouths are successfully created.

The execution time also highly depends on the number N of distinct k-mers we have to deal with. We thus performed experiments on random sequences of growing size. The results are presented in Fig. 6. The execution time grows linearly with N while N remains below a threshold of around 15 million reads. Above this threshold, one observes an exponential growth that could be explained by the fact that the kisSnp prototype uses a hash table instead of a tree for storing and accessing the k-mers information. With a large number of k-mers, the hash table load factor becomes higher than 0.75 increasing the lookup cost, because of an important number of collisions.

Applications to real read data

To test our approach on a real dataset now, we used raw reads from the Escherichia coli Long-term Experimental Evolution Project [START_REF] Lacroix | coli long-term experimental evolution project site[END_REF] whose purpose was to grow Escherichia coli during more than 20 years, conserving a sample each 500 generations. The SNPs found over these generations are listed in the Barrick et. al paper [START_REF] Barrick | Genome evolution and adaptation in a long-term experiment with escherichia coli[END_REF]. An Illumina 1G platform was used for sequencing the samples with reads of length 36 and a high coverage (50x). We focused our attention on the raw reads from the first generation sample, and those from the 20.000 th generation sample. The existence of experimentally validated SNPs is very rare which is the main reason that led us to work on this dataset, for which true positives are known.

Using a custom-made computational pipeline called BRESEQ, Barrick et. al identified 28 SNPs by mapping these two generations of reads against the reference genome CP000819. These 28 SNPs were then experimentally validated.

We used kisSnp for comparing these two sets of reads, forgetting the reference genome. Using as parameters k = 26 and δ = 20, 88 SNPs were found. Of these, 27 of the 28 SNPs found by Barrick et. al were also identified by us, giving a sensitivity of 96%. Our kisSnp method missed one SNP, located position 430835.

To evaluate the potential interest of the remaining 61 SNPs we identified, we mapped them against the reference genome using Maq [START_REF] Li | Mapping short DNA sequencing reads and calling variants using mapping quality scores[END_REF]. Among the 61, 43 correspond to a SNP structure not detected in the Barrick et. al project: the two lips map at the same position with one substitution in the 20.000 th generation sequence. The remaining 18 correspond to suspicious SNPs. Indeed, the two lips do not map to the same position in the genome. Without experimental validation, one however cannot conclude on those 61 detected putative SNPs.

The results obtained with kisSnp are very good on this real dataset, as without a reference genome it was able to find back 27 of the 28 experimentally verified SNPs, and 41 additional ones that correspond to real SNP structures.

Conclusion and future work

We proposed an algorithm for comparing the raw outputs of short reads experiments, typically Illumina ones, with the purpose of finding SNPs between individuals of a same species. This is of particular interest for quickly designing genomic tags without waiting and/or paying for a full genome assembly.

Preliminary results on both simulated and real experimental data are particularly promising. In both cases, kisSnp identifies the SNPs, while not being too sensitive to the parameters used. On a real dataset, kisSnp enabled to find 96% of the SNPs initially detected by mapping to the reference genome. In addition, we propose new SNPs, which could be tested experimentally.

There is clearly room for improvement. For now, the method does not handle heterozygous SNPs, does not take sequencing qualities into account and is limited to pairwise comparison while sets of more than two individuals may be compared. More generally, the three challenges of SNP identification are that the reads contain errors, the genome contains repeats and the read coverage is not uniform. This last item is usually disregarded whereas we notice that it has a significant impact on the results. We expect that several algorithms in the area of NGS bioinformatics could be improved by taking this observation into account.

  Mouth model definition. In a de Bruijn graph of k-mers coming from reads of two sequencing experiments (reads A and reads B), a mouth is composed by: -an upper path of k overlapping k-mers {a 0 ..a k-1 } resulting from the reads of at least set A. This path is called the upper lip of the mouth; -a lower path of k overlapping k-mers noted {b 0 ..b k-1 } resulting from the reads of at least set B. The a i 's and the b i 's differ by one substitution. This path is called the lower lip of the mouth; -a left (resp. right) node, noted c -1 (resp. c k ) that corresponds to a k-mer present in both sets A and B and is connected to both a 0 and b 0 (resp. a k-1 and b k-1 ). These k-mers are called the closing k-mers of the mouth.

Fig. 1 .

 1 Fig. 1. A SNP between two genome fragments (Seq A and Seq B) generates a mouth shape (rightmost frame) in the de Bruijn graph of the k-mers (here k = 4) extracted from Seq A and Seq B. It is assumed in this example that the coverage is exactly 20 (each position of each sequence is covered by 20 reads, thus each position gives rise to 20 k-mers). In the rightmost frame, the number above (resp. below) the nodes indicates the number of occurrences in Seq A (resp. Seq B).

  The left and right closing k-mers c -1 and c k have no counting properties. 4 Algorithm kisSnp for mouth detection Algorithm outline. The algorithm kisSnp takes as input two sets of reads (A and B) coming from two distinct sequencing experiments. The output is a set of pairs of micro-assembled sequences, each of length 2k -1, differing by exactly one substitution located at the central position. Those correspond to putative SNPs detected thanks to the mouth model. The algorithm is divided into three main steps:-For each set A and B, extract the k-mers and their reverse complement and store them in a tree together with their number of occurrences. -Create the mouths (detailed in Section 4.1):• For each possible opening k-mer a op , detect all possible opposite opening k-mers b op distant by one substitution from a op and fulfilling the counting model. • For each pair (a op , b op ), construct the mouth by extending the k-mers to the right and to the left with coherent k-mers (i.e. overlapping on k -1 characters and fulfilling the counting model). • Stop the right and left extensions once the mouth is closed or no extension can be found.

Fig. 2 .

 2 Fig.2. A mouth with k = 4. The symbol ′ -′ stands for a match between two k-mers positions while the symbol ′ X ′ stands for a mismatch. The k-mer a2 is the opening k-mer and the k-mer b2 is the opposite opening k-mer. Here, p = 1 as aop[START_REF] Lacroix | coli long-term experimental evolution project site[END_REF] = bop[START_REF] Lacroix | coli long-term experimental evolution project site[END_REF].

Fig. 3 .

 3 Fig. 3. Number of SNPs (read-coherent mouths) found by kisSnp. Fragment of the human chromosome X, 137 SNPs to find (symbolised by the horizontal line). The "n×" values indicate the coverage used while simulating the reads.

Fig. 4 .

 4 Fig. 4. Results on MC58 and the randomised MC58 sets while looking for 2272 SNPs (horizontal line) with δ = 20. The "n×" values indicate the coverage used while simulating the reads, "rand." stands for results on the randomised MC58 set.

Fig. 5 .

 5 Fig. 5. Influence of k on the execution time. Data: set KIF4A (described Section 5.1), δ = 30.

Fig. 6 .

 6 Fig. 6. Execution time with respect to the number N of distinct k-mers of length k = 25, with δ = 20, Bernoulli random sequences.

Acknowledgments

We are grateful to Jean Marc Aury, who provided the empirical Illumina error model, to Matteo Brilli who pointed out to us the E.coli experiment, to the GenOuest platform who supported extensive computations and to the INRIA ARC Alcovna for financial support.


Gaston Miño

Thomas E Mallouk

Thierry Darnige

Mauricio Hoyos

Jeremy Dauchet

Jocelyn Dunstan

Rodrigo Soto

Yang Wang

Annie Rousselet

Eric Clement

Gastón Miño

Enhanced diffusion due to active swimmers at a solid surface

Keywords: numbers: 87, 16, -b, 05, 65, +b

des établissements d'enseignement et de recherche

Since the pioneering work of Wu and Libchaber [1] considerable efforts have been made to understand hydrodynamic properties of active suspensions. Generally speaking, this is the name borne by fluids laden with selfswimming entities such as bacteria [2][3][4][5], algae [6,7] or collections of active artificial swimmers [8]. Assemblies of microscopic motors dispersed in a fluid display emergent properties that differ strongly from passive suspensions. The momentum and energy transfer balances as well the constitutive transport properties are deeply modified by the momentum sources distributed in the bulk [2,9]. Some of these anomalous properties have already been identified such as active diffusivity [1,2], anomalous viscous response [7,10], active transport and mixing [11] as well as the possibility to use fluctuations to extract work [12]. The presence of living and apparently gregarious entities also offers the possibility to move collectively and organize at the mesoscopic or macroscopic level in the form of flocks and herds [9,14]. Similar collective effects were also identified in suspensions of self-propeled inorganic particles [15]. In the bulk, swimming bacteria with flagella such as E. coli create in the far field-limit, a force-dipole velocity field and consequently, experience a hydrodynamic attraction toward surfaces [16]. Then, it has been observed that E. coli smooth out their runand-tumble movement and spend long times parallel to the surface undergoing circular motion as a consequence of the torque-free condition [17,18]. When the concentration becomes large, the E. coli population eventually associates collectively to form a bio-film. Even in the low concentration limit, the quantitative analysis of the near surface motion increase tremendously in complexity, a reason being the close field hydrodynamic forces that become prevalent and require a complex treatment of the lubrication hydrodynamic fields. However, even in this frame of description, it remains unclear whether the motion close to the surface is hydrodynamically stable and if the presence of thermal noise is essential to account for the bacterium dwelling time at a surface [19]. Beyond the hydrodynamic interactions, more complex ingredients may come into play such has the surface interaction potentials (electrotastic or van der Waals) [20] or more refined details of the bacterium physiology such as swimming speed variations and desynchronization during bacteria cell cycle [21,22]. From the perspective of providing a fully consistent treatment of active hydrodynamics, with important applications for understanding bacterial transfer in biological micro-vessels, microfluidic devices, or the formation of bio-films, a reliable description of fluid activity in the vicinity of a solid surface is strongly needed. In this letter, to tackle this open and timely question, we compared the behavior of two kinds of active micrometric swimmers: wild type E. coli K12 and artificial self-propelling rods [8], with completely different propulsion mechanisms. In both cases, we monitor the swimmers' motions and their ability to activate, beyond Brownian motion, passive tracers, hence characterizing the active momentum transfer to the fluid.

Following the experimental procedure described in Ref. [16], wild type E. coli K12 were grown overnight in rich medium (LB). After washing, they were transferred into MMAP, a motility medium supplemented with K-acetate (0.34 mM) and polyvinyl pyrolidone (PVP: 0.005%). They were incubated for at least an hour in that medium and, in some cases, so-called "baby cells" were selected by centrifugation and resuspended in MMAP. To avoid bacterial sedimentation (isodense conditions), Percoll was mixed with MMAP (1vol/1vol). We checked that under these conditions, the suspending fluid was still Newtonian (viscosity η = 1.28 × 10 -3 Pa s at 22 • C). The overall concentration of bacteria was controlled such as to prepare suspensions between 10 9 and 10 10 bact/ml. To study the The colormap goes from blue for vanishing probability to red for high probability. Two clusters are identified, centered at (0.9, 0.3) and (0.17, π/2), corresponding respectively to the "active" and "diffusive" swimmers.

effect of bacterial activity on the diffusivity of passive tracers, latex beads of 1 or 2 µm diameter (Beckman-Coulter, density ρ = 1.027 g/ml) were added to the suspensions. Experiments were performed in 110µm thick chambers, built with two horizontal microscope cover slips separated by a glass spacer. To avoid sticking, cover slips were coated with PVP. The biological sample consisted in a drop of liquid (20 µl) placed between the two slides. The suspension was observed under an inverted microscope (Zeiss-Obzerver, Z1-magnification 40X) connected to a digital camera. The observation field ∆V was 96×128µm 2 and 5µm in depth. Short video sequences of 20s duration at 20 frames/s (fps) were recorded to track the bacterial motion. In a first series of experiments, we measured the bacterial density profile through the entire height of the chamber. We obtained profiles similar to the ones published by Berke and coworkers [16], namely a flat density in the bulk and a strong density increase near the surfaces. However, the wild-type E. coli we used was significantly less attracted by the surfaces (2.5 times increase in density within 10µm of the surface) than a mutant E. coli strain that does not display tumbling motion [17,18]. Another series of experiments was performed with bimetallic Au-Pt self-propelled rods (length 1.2 µm and diameter 0.4 µm) that are very similar in size to the E. coli cell body (1 to 2 µm long, 0.8 µm diameter) but have no flagella (15 µm long for E. coli) (see inserts in Figs. 1c and1f). They also have a much higher density (ρ = 17 g/ml). In the presence of 1 to 10% hydrogen peroxide, these particles are propelled in the axial direction towards the platinum end by the catalytic decomposition of the peroxide fuel [8]. Recent experiments and simulations are consistent with self-electrophoresis as the dominant propulsion mechanism [23,24]. While other mechanisms have been proposed, it is clear that the mode of propulsion is intrinsically different from the flagellar propulsion of bacteria. The experiments were conducted in a similar fashion to those involving E. coli, but in an open chamber (without the upper wall), in order to allow the oxygen bubbles produced in the reaction to escape from the cell. The concentration of H 2 O 2 was varied, as well as the concentration of active rods (n = 3 -20 × 10 6 rods/ml). We also used passive tracers (1µm diameter beads, Dynal-My-one, density ρ = 1.8 g/ml) or 2 µm diameter latex beads (Beckman-Coulter, density ρ = 1.027 g/ml)) to follow the activation of the fluid by the Au-Pt rods. In all cases, all the particles in the suspension were localized at the bottom of the chamber as a consequence of sedimentation. Again, short videos (20s duration at 20f ps) were used to track the self-propelled rods.

In the following, we only focus on bacteria and rods moving close to the surface (less than 5µm). In both cases, we observed that not all swimmers display similar trajectories. This was expected as for wild type bacteria, the run or tumble dynamics may depend strongly on the microenvironment or on the position in the cell cycle. For Au-Pt rods, this is also consistent with previous observations that even within a single batch, electrochemically grown rods have a range of catalytic activity [8]. We developed a tracking program to analyze the short videos and obtained tracks for each swimmer present in the field. We identified two major types of motion: a ballistic and a diffusive one (see Fig. 1(a,b,d,e)), and the swimmers that follow these motions are called "active" swimmers and "diffusive" swimmers, respectively. To discriminate systematically all the tracks, two parameters were defined. The first parameter ⟨|θ|⟩ is the mean angle between two successive steps. For example, ⟨|θ|⟩ = 0 for straight trajectories and ⟨|θ|⟩ = π/2 for a purely random walk. The second discriminating parameter is based on the minimal circle diameter L that encompasses a given trajectory of duration T . For an acquisition time δt and a mean step size δr, the number N c = Lδt T δr is computed. When N c is close to 1, the trajectory is associated with a straight line, whereas when N c is small, its value points to diffusive motion. Therefore, each track is associated with these two numbers and in the (N c , ⟨|θ|⟩) parameter-space, we could identify two clusters that clearly differentiated the active and diffusive swimmers (see Fig. 1(g)). Nevertheless, for very small or interrupted trajectories, the separation procedure remained ambiguous, so we systematically discarded tracks shorter than 10 steps. In the case of artificial swimmers, we also managed to control a priori the fraction of active swimmers by adding inactive rods (made only of gold) and keeping the total number of rods at the surface constant. According to the trajectory classification, a fraction ϕ A of active swimmers was determined. Thus, for a mean number ⟨N ⟩ of swimmers, identified in the field of vision, we define a density of "active swimmers" as n A = ϕ A ⟨N ⟩/∆V . In Fig. 1, we display tracks during a time lag τ = 1.5 s, for two populations of swimmers (a-b, bacteria and d-e, Au-Pt rods), having different ϕ A (a,d small active-fraction and b,e high active fraction). In Figs. 1(c) and 1(f), we present the mean velocities of active swimmers V A as a function of their fraction ϕ A , for different experiments with bacteria and active rods. In the case of the bacterial suspension, we also tried several synchronization protocols to select bacteria, at different position in the cell cycle, showing different swimming characteristics. We were able to produce "baby-bacteria" populations (1N short cells: 1.12 µm long) which were found to have a ϕ A larger than the more mature bacteria populations (2N long cells: 2.5 µm long). We took advantage of this difference to look at the influence of ϕ A on V A in bacterial suspensions. On Fig. 1(c) we display 6 independent experiments performed with populations having a majority of 1N cells (brown,red and green), a mixture of 1N and 2N cells (black), a majority of 2N cells (pink, blue). For each sample, ϕ A was taken from suspensions showing an average of 100 or 200 bacteria in the observation field (represented by ( 1) and ( 2) on Fig. 1c, respectively). We observed that ϕ A has little influence on V A but that V A could be different according to bacteria position within the cell cycle. In the case of the artificial swimmers suspension, Fig. 1(f) shows a stronger independence of V A on ϕ A , which is changed by varying the proportion of inactive Au-rods.

In the following, we will relate the motion of the passive tracers to the number of active swimmers and their mean velocity. Passive tracer trajectories were analyzed (about 10 tracks per long video of 300 images at 1f ps for bacteria and 40s sequences at 8f ps for Au-Pt rods). From these tracks, the mean passive tracer diffusion coefficient D P was extracted consistently using two independent methods. In the first case, the mean square displacements were measured for individual particles and in the other case, for particles as pairs, in order to eliminate residual drift. In Fig. 2 the passive tracer diffusivities are displayed for all the experiments presented in Fig. 1(c,f). D P values are displayed as a function of the product J A = n A V A , that we call the "activity flux". For experiments performed with the same tracer size, we observe a collapse of all data onto a linear curve.

D P = D B P + βJ A (1)
where D B P is the Brownian diffusivity of the latex particles in the vicinity of the surface, in the absence of swimmers. Note that due to lubrication forces, this value is smaller than the Brownian diffusivity expected in the bulk (D B = k B T /3πηd) and, for the parallel motion, they are related by D B P = αD B , where α < 1 is the parallel drag correction factor [25][26][27]. The α factor depends on the bead distance to the surface, vanishing at contact and going asymptotically to one at large distances. The beads are not at a fixed distance to the surface, but they are distributed according to the Boltzmann's factor exp(-m * gz/k B T ), where m * is the buoyant mass. Therefore α must be averaged with this factor. In the active rods experiments we obtained α = 0.7 both for the the buoyant d = 1µm and d = 2µm latex spheres, value that agrees with the theoretical prediction given above (α = 0.64 for d = 1µm and α = 0.74 d = 2µm). In the experiments with bacteria, the suspended beads are almost isodense but they sediment anyway. The experimental fit gives α = 0.85. This value allows us to infer the density mismatch ∆ρ = 0.008 g/ml, which is consistent with the experimental preparation. The collapse holds also for bacterial populations at different maturation stages (1N, 2N or unsynchronized mixtures). In the expression (1), the prefactor β is a length to the fourth power. It varies from 5µm 4 = (1.5µm) 4 for bacteria, to 13µm 4 = (1.9µm) 4 for active rods, but seems to be almost independent of the passive tracer size. The enhanced diffusivity in (1), proportional to n A and V A , can be understood as a result of a series of elementary encounters between active swimmers and the tracers: the number of encounters per unit time is proportional to n A V A . On the other hand, low Reynolds dynamics points out that the tracer displacement at each encounter is independent of the swimmer velocity and depends only on geometrical factors: the impact parameter, the swimmer dimensions and weakly on the tracer size through the Faxén correction of passive transport [28,[START_REF] Happel | Low Reynolds Number Hydrodynamics: with special applications to particulate media[END_REF]. The β factor comes out from averaging the tracer's displacements but its computation is difficult because it requires a correct modeling of the near field interactions between the swimmer and the tracer, taking into account the detailled swimmer geometry and the effect of the surface.

In conclusion, we have characterized active momentum transfer in the vicinity of a solid surface for two active suspensions (wild type bacteria and artificial selfpropelled swimmers). The effect was measured using the diffusion enhancement of a passive tracer. In spite of the a priori complexity of the hydrodynamics and essential differences in the propulsion modes, we demonstrated that the effect emerges quantitatively in a similar way. The resulting diffusion coefficient is the sum of the Brownian contribution near the wall and an active part, proportional to the product of the density of active swimmers and their mean velocity at the surface. The proportionality factor, scaling as the 4 th power of a micron size length, encompasses the details of momentum transfer for each swimmer and is found to be weakly (if at all) sensitive to the probe diameter. Importantly, discriminating between so-called "active" and "diffusive" swimmer trajectories was crucial for predicting the induced transport phenomenon and we have developed a protocol to make such a distinction. The functional dependence of the enhanced diffusivity is explained in terms of successive interactions between a single swimmer and the tracer, each one producing a net displacement. Our result provides a justification to pursue a quantitative determination of such encounters based on simple hydrodynamic models [28].

FIG. 1 :

 1 FIG. 1: Identification of the swimmer populations by tracking "active" swimmers (red tracks) and "diffusive" swimmers (blue tracks), ϕA is the corresponding fraction of active swimmers. Figs. 1(a,b,g) correspond to E. coli (see inset in 1c) and Figs. 1(d,e) to Au-Pt rods (see inset in 1f). The round black circles in (a,b) are 2µm latex beads; the white small circles in (d,e) are 1µm Dynal beads. Figs. (c) and (f) display the relation between the active swimmers mean velocity VA and ϕA. On Fig. 1(c), 6 independent experiments with E. coli: 1N cells (brown, red and green), mixture of 1N and 2N cells (black) and 2N cells (pink, blue). Labels (1) and (2) are for N = 100 and N = 200 cells in the observation field respectively. Fig.1(f), for Au-Pt rods with varying proportions of inactive Au rods (2 independent experiments). Fig.1(g), density probability of the observed bacterial tracks in the (Nc, ⟨|θ|⟩) space. The colormap goes from blue for vanishing probability to red for high probability. Two clusters are identified, centered at (0.9, 0.3) and (0.17, π/2), corresponding respectively to the "active" and "diffusive" swimmers.

FIG. 2 :

 2 FIG. 2: Enhanced diffusivity DP of passive tracers measured as a function of JA the "activity flux". Squares and circles symbols are related repsectively to tracer diameters 2µm and 1µm. Fig.2(a) corresponds to the bacterial suspensions; 1N cells (red and green) and unsynchronized (black, blue and pink). Each color defines an experiment performed over a range of bacterial dilution. Fig.2(b) corresponds to suspensions of Au-Pt-rods. Mixture of active and inactive Au-rods (red circles, brown squares); Various H2O2 concentration from 2,5% to 20%, (green, blue, black circles). The dashed lines represent linear fits and the error bars one standard deviation.

Acknowledgment. We thank Pr. David Grier for discussions on the tracking programs, the financial support of the PGDG Foundation, the Alfa-SCAT program, Sesame Ile-de-France, Fondecyt Grants No. 1061112, No. 1100100, Fondap Grant No. 11980002, ECOS C07E08.


Prof Jean-Luc L Autran

Karim Nehari

Daniela Munteanu

Irphe Nehari

Compact Modeling of the Threshold Voltage in Silicon Nanowire MOSFET including 2D-Quantum Confinement Effects

Keywords: 2D-quantum confinement effects, threshold voltage, quantum nanowire MOSFET, compact modeling, 3D numerical simulation compact modeling, threshold voltage, quantum nanowire MOSFET, 2D-quantum confinement effects, 3D numerical simulation

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

Introduction

Multi-gate MOS transistors are widely recognized as one of the most promising solutions for meeting the roadmap requirements in the deca-nanometer scale of device integration [1]. Recently, these architectures are also presenting a growing interest for the integration of post-conventional Silicon devices and molecular-based devices, since the active region of the transistor, i.e. the conduction channel, can be alternatively played by a semiconductor nanowire, a Carbon nanotube, a molecule or an atomic linear chain [START_REF] Bescond | 3D Quantum Modeling and Simulation of Multiple-Gate Nanowire MOSFETs[END_REF][START_REF] Bescond | Atomic-scale modeling of Double-Gate MOSFETs using a tight-binding Green's function formalism[END_REF]. In all cases, these structures exhibit a superior control of short channel effects resulting from an exceptional electrostatic coupling between the conduction channel and the surrounding gate electrode.

In the case of Silicon-based CMOS technologies, a wide variety of multi-gate architectures, including Double-Gate (DG), Gate-All-Around (GAA), PiFET, FinFET, Rectangular or Cylindrical nanowire MOSFETs, linear atomic-chain transistors has been proposed in the recent literature [1][START_REF] Bescond | 3D Quantum Modeling and Simulation of Multiple-Gate Nanowire MOSFETs[END_REF]. For all these structures, one of the identified challenges remains the development of compact models taking into account the main physical phenomena governing the devices at this scale of integration.

In this work, a compact model for the threshold voltage (V T ) of long-channel quantum Si nanowire (QW) MOSFETs is developed. Schematic representations of different QW architectures are shown in Fig. 1. In such devices for which the lateral dimensions, i.e. the Si film thickness and its width, are ultimately downscaled to a few nanometers, quantummechanical carrier confinement effects are exacerbated in the two (y-and z-) directions perpendicular to the transport one (x-). As we will detail in the following, our present approach is thus based on the decoupled solution of the Poisson and Schrödinger equations in the (y,z) plane of the structure. Threshold voltage obtained from this new model will be

Threshold voltage modeling

V T definition and criterion

The band diagram in the vertical (y-) direction (middle of the structure) for a n-channel QW transistor is schematically represented in Fig. 2. The potential reference is assumed to be the electron Fermi level of the n+ source region. The potential, Ψ, is defined as the band bending with respect to the intrinsic Fermi level in the silicon film. It is important to note that the surface potential, Ψ s , as defined here, is different from the surface potential in conventional bulk MOSFET, Ψ sc , which is given by the band bending between the silicon surface and the substrate volume. The relation between the two potentials is Ψ sc =Ψ s +φ F where φ F =(kT/q)ln(N A /n i ) is the Fermi potential, N A is the film doping and n i is the intrinsic carrier concentration. In conventional MOSFET, the threshold voltage is usually defined as the gate voltage necessary for obtaining a band bending of Ψ sc =2φ F , which means Ψ s =φ F . This definition does no more apply to QW MOSFET: our extensive numerical simulation fully confirms this remark and shows that surface potential at threshold, Ψ s T , is different from φ F . Other definitions for V T exist, as derived from the constant-current method or the maximum of ∂ 2 I D /∂V G 2 . In a recent work [START_REF] Munteanu | Unified Analytical Model of Threshold Voltage in Symmetric and Asymmetric Double-Gate MOSFETs[END_REF], we used this latter criterion, but the resulting analytical expression of the threshold voltage is quite complicated and its demonstration demands heavy mathematical manipulations. For obtaining a more simple expression of V T , we consider here the well-known definition where V T is the gate voltage for which the inversion charge Q iT in the transistor channel (inversion charge at threshold) is given by [START_REF] Munteanu | Unified Analytical Model of Threshold Voltage in Symmetric and Asymmetric Double-Gate MOSFETs[END_REF][START_REF] Harrison | Electrical characterization and modelling of high-performance SON DG MOSFETs[END_REF]:

ox ox T G iT C q kT C ) V V ( Q = - ≅ (1)
where V G -V T is assumed to be ~kT/q around threshold. We verified that this value of Q iT is very close to the inversion charge obtained by the numerical simulation of multi-gate structures with different channel doping levels, channel thicknesses and gate oxide thicknesses.

Once the inversion charge at threshold has been defined, the determination of V T can be easily conducted in two steps: i) express this inversion charge (Q i ) in the QW channel as a function of V G and ii) numerically solve the equation Q i (V G )=Q iT . In the following, the calculation of Q i is obtained from the decoupled solution of the Poisson equation and Schrödinger equations. The first equation leads to the band bending profile in the Si film;

the second one gives energy levels in the QW (considered as a 2D quantum-well) used to evaluate the inversion charge.

Device electrostatics

Considering a multi-gate QW structure represented in Fig. 1, our approach is based on a parabolic approximation of the vertical (y-) potential distribution in the film at threshold, which has been fully confirmed by numerical simulation in previous studies [START_REF] Munteanu | Unified Analytical Model of Threshold Voltage in Symmetric and Asymmetric Double-Gate MOSFETs[END_REF][START_REF] Autran | Quantum-Mechanical Analytical Modeling of Threshold Voltage in Long-Channel Double-Gate MOSFET with Symmetric and Asymmetric Gates[END_REF]. This result can be, of course, generalized for the other direction (z-axis) perpendicular to the channel axis. We can write: With definitions in Fig. 1 and applying the Gauss law, the boundary conditions at the two oxide/silicon interfaces are identical and give:

2 Si s y y t β + β - Ψ = Ψ (2)
F s ox ox Si FB G E t V V φ + Ψ + ε ε = - (4)
where V FB is the flat band voltage, ε Si and ε ox are the silicon and oxide permitivities, respectively, t ox is the oxide thickness and E is the electric field at the interface :

Si 0 y t dy d E β = Ψ - = = (5)
Expressions ( 2) to (5) form a complete set of equations to analytically evaluate the band bending profile in the silicon film under threshold condition. As shown in Fig. 3 in comparison with numerical simulation, relation ( 2) is a very good approximation for the potential distribution in the film at threshold for various doping levels. The potential is symmetric with respect to the middle of the silicon film and has a quasi-parabolic dependence.

Quantum-mechanical evaluation of the inversion charge

The second step of our model consists in calculating the inversion charge in the silicon quantum-wire. In a first approach, this QW can be simply depicted as a 2D quantum-well with infinite walls, for example a rectangular well for DG or GAA structures, a cylindrical well for the cylindrical QW transistor (Fig. 1), etc. The calculations are developed here for the rectangular infinite well (i.e. the GAA structure) with carriers confined in the (y,z) plane. Assuming, in a first approximation, that the conduction band profile is flat inside the well, simple analytical expression for the quantum energy levels can be used. In the particular case of <100>-oriented Silicon, three different series of quantum energy levels must be considered in the evaluation of the inversion charge, to properly take into account the effective mass ellipsoid characteristics (see Fig. 4):

⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ π + π = π + π = π + π = 2 t 2 2 2 2 Si t 2 2 2 3 2 l 2 2 2 2 Si t 2 2 2 2 2 t 2 2 2 2 Si l 2 2 2 l W m j t m i ) j , i ( E W m j t m i ) j , i ( E W m j t m i ) j , i ( E h h h h h h (6)
where m l ≈0.916×m 0 and m t ≈0.19×m 0 are the longitudinal and transversal effective mass of electrons, respectively.

For each energy level E k (i,j), referenced with respect to the bottom of the quantum-well, the contribution to the total charge integrated over the Si film cross-section (W×t Si ) is:

∫ ρ × = +∞ k E k D 1 k dE ) E ( f ) E , E ( 2 ) j , i , E ( n (7)
where the factor 2 accounts for the number of equivalent valleys, ρ 1D (E,E k ) is the 1D density-of-states and f(E) is the Fermi-Dirac distribution function:

k 2 / 1 2 D k E E 1 m 2 1 ) E , E ( - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ π = ρ h (8) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - + = kT E E exp 1 1 ) E ( f F (9)
where m D is the 1D density-of-states effective mass (m D =m t for k=1,2 and m D =m l for k=3).

In Eq. ( 7), the calculated charge depends on V G via the position of the Fermi level with respect to the bottom of the conduction band (assimilated to the bottom of the quantum well) and thus with respect to the different energy levels. The inversion charge to be considered for V T determination using Q i (V G )=Q i T is finally:

F o r P e e r R e v i e w O n l y ∑ ∑ ∑ = k i j k G i ) j , i , E ( n q ) V ( Q (10)
In order to enhance the accuracy of the model, a first-order perturbation correction of the energy levels given by Eq. ( 6) has been conducted [START_REF] Autran | Quantum-Mechanical Analytical Modeling of Threshold Voltage in Long-Channel Double-Gate MOSFET with Symmetric and Asymmetric Gates[END_REF][START_REF] Baccarani | A Compact Double-Gate MOSFET Model Comprising Quantum-Mechanical and Nonstatic Effects[END_REF]. The Hamiltonian V of the perturbation is simply evaluated from the electrostatic potential profile (Eq. ( 2)), assuming that a similar profile is obtained in the second perpendicular direction:

) z W ( z ) y t ( y ) z , y ( V z Si y - β + - β = (11)
where β y and β z are two constants depending on β.

The first-order correction to apply to the energy levels in the well is given by:

j , i j , i j , i V E ϕ ϕ = ∆ ( 12
)
with ϕ i,j is the wave function related to levels defined with quantum numbers i and j:

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ π ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π = ϕ W z j sin W 2 t y i sin tSi 2 ) z , y ( Si j , i (13)
An analytical evaluation of

j , i E ∆ immediately gives: [ ] ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π + β + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ π + β = ∫∫ ϕ = ∆ × 2 2 2 z 2 2 2 Si y tSi W 2 j , i j , i j 3 1 6 W i 3 1 6 t dydz ) z , y ( V ) z , y ( E (14)
The corrected energy levels to be considered in Eq. ( 10) for the evaluation of Q i are finally:

j , i k k E ) j , i ( E ) j , i ( E ~∆ + = (15)

Model validation and discussion

The present model has been validated by comparison with numerical simulation using a 2D/3D quantum-mechanical simulation code coupling the Poisson-Schrödinger system with the Drift-Diffusion transport equation (Balmos3D, [START_REF]Balmos3D simulation code[END_REF]). First, considering a single direction for carrier confinement (t Si downscaling at fixed infinitely large W). Fig. 6 shows that a very good match is obtained between analytical and 2D numerical results in a rectangular QW structure. In particular, it is shown that the model is able to very well predict the dependence of V T with t Si and also with N A . In a second step, both y-and zdirections for carrier confinement have been considered and the compact model has been validated using 3D numerical simulation, as shown in Fig. 7. One again the fit between compact model and numerical data is very satisfactory. Fig. 7 and Fig. 8 show that the threshold voltage increase, due to carrier confinement in the z-direction, becomes significant for W lower than ~20nm. It is also interesting to remark that, when W is scaled down, the threshold voltage increase becomes more important as long as the film thickness decreases. Indeed, when switching from 1D to 2D confinement, the threshold voltage increase is more pronounced due not only to the increase of the distance between the energy levels and the Fermi level, but in the same time, to the change of the density-ofstates. In other words, the band bending must be accentuated to obtain the same inversion charge with higher energy levels and lower number of quantum states per level. Finally, Fig. 9 shows the variations of V T predicted by the compact model as a function of both W and t Si for a rectangular GAA QW.

Conclusion

In conclusion, a compact modeling of the threshold voltage for quantum-wire MOSFETs

Ψ s Ψ E c E v E i V G V G E F y Gate Gate

  directly calculated using a 2D/3D quantum-mechanical simulation code coupling the Poisson-Schrödinger system with the Drift-Diffusion transport equation (Balmos3D).

Figure 1 .

 1 FIGURE CAPTIONS

Figure 2 :

 2 Figure 2: Schematic band diagram for a QW transistor in the y-direction (middle part of the

Figure 3 :

 3 Figure 3: Potential distribution in the silicon film at threshold in a symmetric DG-QW

Figure 4 :

 4 Figure 4: Schematic illustration of electron confinement in a Silicon quantum-wire with a

Figure 5 :

 5 Figure 5: Illustration of the first-order perturbation method used to correct energy level

Figure 6 .

 6 Figure 6. Quantum V T predicted by the analytical model in an infinitely large DG-QW

Figure 7 .

 7 Figure 7. Quantum V T predicted by the analytical model as a function of W and t Si for a

Figure 8 .

 8 Figure 8. ∆V T predicted by the analytical model as a function of the channel width in thick

Figure 9 .

 9 Figure9. V T (W,t Si ) surface map for a rectangular GAA QW long-channel transistor (L=100nm, V DS =50mV, N A =10 16 cm -3 ) as calculated by the proposed model.

  Figure 1.

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

Figure 5 .

 5 Figure 5.

Figure 6 .

 6 Figure 6.

Figure 7 .

 7 Figure 7.

  Figure 8.

  Figure 9.

http://mc.manuscriptcentral.com/tandf/jenmol

ACKNOWLEDGMENTS

This work is supported by the European Commission in the framework of the Network of Excellence on Silicon-based nanodevices (SINANO, contract IST-506844).


A C Brańka
email: branka@ifmpan.poznan.pl

D M Heyes
email: d.heyes@surrey.ac.uk

The role of particle softness in determining the value of Poisson's ratio for soft sphere solids

Keywords: inverse-power potential, elastic constants, Poisson, auxetics, molecular dynamics simulations PACS numbers: 62.20.Dc, 02.70.Ns, 81.40.Jj, 07.05.Tp, 81.40.Wx inverse-power potential, elastic constants, Poisson's ratio, auxetics, molecular dynamics simulations

The influence of particle softness on the Poisson's ratio of model solids has been investigated.

We have used the repulsive inverse power potential (∼ r -n for particle separations, r) between the particles, which is conveniently characterised by one adjustable parameter, ε = 1/n. For large ε the interaction is soft whereas in the ε → 0 limit the particles approach hard spheres. The pressure and elastic constants of the solid phase have been calculated at various densities with constant temperature molecular dynamics simulation for a range of the softness parameter in the range, n > 12. Densitysoftness surfaces of these quantities were determined which revealed hitherto unrecorded trends in the behaviour of the elastic moduli and Poisson's ratio. It was found that the pressure and some elastic properties e.g. the C 12 elastic constant and the bulk modulus, manifest a maximum value or 'ridge' on this surface. The height of the maximum increases with density and interaction steepness (small ε). The Poisson's ratio varies essentially linearly with softness and is relatively insensitive to density. However, at higher densities and for larger steepness a considerable lowering of the Poisson's ratio is observed. In order to identify possible mechanisms for reducing the value of Poisson's ratio, ν, the fluctuation and Born-Green contributions were analyzed. Changes in the Poisson's ratio are mainly determined by the fluctuation contribution which can cause a considerable decrease as well as increase of its value.

How a body behaves under applied stress or tension is an important question for technical and industrial applications as well as having scientific interest. For relatively small deformations the elastic response is characterised by a set of elastic moduli [START_REF] Wallace | Statistical Mechanics of Elasticity[END_REF]. These quantities are also relevant in the study of phase stabilities as the symmetry elements of solid phases are well reflected in the number of independent elastic constants. In fact, elastic constants can be considered to be order parameters indicative of phase stability and are therefore relevant in the study of structural transformations [START_REF] Landau | Theory of Elasticity[END_REF]. Indeed the fundamental distinction between a solid and fluid is marked by the absence of a (zero frequency) shear modulus in the latter [START_REF] Hess | Shear modulus of fluids and solids[END_REF].

Knowledge of the elastic constants greatly facilitates the design and manufacture of materials with desirable physical properties. Most of the naturally occurring materials decrease their transverse dimensions under uniaxial load. This is such a commonly observed phenomenon that only relatively recently has it been realised that the inverse situation can be observed in nature and introduced in constructed materials [START_REF] Lakes | Foam Structures with a Negative Poisson's Ratio[END_REF], [START_REF] Caddock | Microporous materials with negative Poisson's ratios. I. Microstructure and mechanical properties[END_REF] [6], [START_REF] Alderson | How to make auxetic fibre reinforced composites[END_REF]. Materials that expand in all directions when pulled in one direction are called auxetics (from the Greek 'auxetos' which can be increased). Auxetic materials are not only interesting scientifically for their rare and counter-intuitive elastic behaviour but also for a number of potentially useful technological applications [START_REF] Evans | Auxetic Materials: Functional Materials and Structures from Lateral Thinking![END_REF]. This has focussed attention on what microstructurally is required to give rise to auxetic behaviour.

Even though auxetic behaviour originates on a wide range of length scales in very different naturally-occurring and man-made materials, its defining characteristic, in terms of the Pois-

ν = - η T η L , (1)
where η T and η L are the strains in orthogonal (transverse and longitudinal) directions.

Poisson's ratio, the Young's, shear and bulk moduli are the basic parameters defining the elastic response of isotropic materials. An auxetic response means a negative Poisson's ratio, which is a convenient situation as one can basically concentrate on one well defined quantity irrespective of the complexity of the considered system or the lengthscale on which the effect originates. In general, in many materials, physical properties depend on direction and also the Poisson's ratio depends on the direction of the applied stretching. In anisotropic materials, such as crystals, ν can be positive in one direction and negative in another direction.

Materials with negative Poisson's ratio are quite rare in nature so the impetus is for them to be manufactured. However, despite considerable progress, the designing and preparation of materials with ν < 0 still remains a nontrivial task and a real challenge. In particular little progress has been reported on the synthesis of auxetics that derive this behaviour at the molecular level [START_REF] He | Toward molecular auxetics: Main chain liquid crystalline polymers consisting of laterally attached para-quaterphenyls[END_REF]. This situation can be attributed to a lack of sufficient fundamental knowledge on mechanisms leading to ν < 0. Most of the modelling to date has concentrated on idealized microstructures which by cooperative motion of the members exhibit expansion when they are stretched. The re-entrant honeycomb structure is the most well known example that gives auxetic behaviour, and has almost become synonymous with auxeticity [11] [12].

Cellular structures composed of geometrical elements which are vital for a negative ν, such as the re-entrant honeycomb structure, can be achieved on different length scales and constitute today the primary source of auxeticity [START_REF] Spadoni | Global and local linear buckling behavior of a chiral cellular structure[END_REF]. Recently, however, a new approach has been discovered for generating auxetics which does not require re-entrant features but instead employs the relative rotation of rigid units [START_REF] Alderson | Molecular Origin of Auxetic Behavior in Tetrahedral Framework Silicates[END_REF]. These rigid units appear not to be restricted to just a few geometries. In tetrahedral frameworks, the auxetic effect is achieved through synchronization of the cooperative movements of the tetrahedral units [START_REF] Alderson | Molecular Origin of Auxetic Behavior in Tetrahedral Framework Silicates[END_REF]. Thus, it is likely that the auxetic effects may also occur if a particular form of translational-rotational coupling exists in the system. In contrast to several proposed models which are 'structural' or mechanical in nature there are only a few thermodynamic models exhibiting auxeticity (by 'thermodynamic' here we mean that the effect originates at the molecular level). The example in this category is the thermodynamically stable system of 2D hard cyclic hexamers [START_REF] Wojciechowski | Constant thermodynamic tension Monte Carlo studies of elastic properties of a two-dimensional system of hard cyclic hexamers[END_REF][START_REF] Wojciechowski | Negative Poisson's ratio in a two-dimensional "isotropic" solid[END_REF].

The connection between auxeticity and the form of interaction between constituent objects or particles still needs further clarification. In particular the relationship between an interaction, or its different parts, and ν is still not well understood. Establishing such a connection would help us understand mechanisms leading to auxetic behaviour and possibly help to explain why auxetic materials are so rare in nature and not easily manufactured. In fact, the role played by particle interactions in determining the elastic constants and Poisson's ratio behaviour in general, needs to be investigated in a more systematic manner as part of this exercise.

As a step in this direction, we investigate the influence of repulsive interaction softness on the elastic properties, and in particular on Poisson's ratio. The softness of the repulsive core is an important feature of any particle interaction. For any ensemble of interacting particles to be thermodynamically stable a repulsive core is required in the interparticle interaction to prevent the system from collapsing [START_REF] Ruelle | Statistical Mechanics Rigorous Results[END_REF]. The nature of this core depends on the chemical architecture of the particle, and can vary considerably from very soft to extremely hard. It is worth noting that particles of variable softness, usually defined through an effective repulsive potential, can now be manufactured. There is, for example growing interest in soft semi-solid systems made from colloidal and microgel particles with variable softness [17, 18][19]. So in these systems as well, establishing any general correlation between elastic properties and microscopic particle softness is desirable and can be helpful in developing better ν-controlled materials.

Among the different model interparticle potentials, the inverse-power or soft-sphere potential form is particularly suitable for investigating the role of particle softness. First of all it has the simple analytic form,

φ(r) = ( σ r ) n , ( 2
)
where r is the separation between two particles, σ is the particle diameter, sets the energy scale and n is a parameter determining the potential steepness (the softness is ε ∼ n -1 ).

Thus, the single parameter, n, changed continuously can cover a wide spectrum of practically important systems from the very soft to the extremely hard, i.e., from the long-ranged Coulomb interaction (n = 1) to the hard sphere system (n → ∞). Apart from its simplicity the inversepower potential has many features that make this potential attractive as a model system.

From its definition in (2), we see that the inverse power potential is a self-similar function, which mixes the energy and distance scales. It can be simply represented as r * -n , where r * = r/σ 1/n . An exceptional feature of the inverse power system is that the configurational properties do not depend upon the density and temperature separately but upon a particular dimensionless combination of them [START_REF] Hoover | Thermodynamic Properties of the Fluid and Solid Phases for Inverse Power Potentials[END_REF]. In effect the properties computed along one isotherm are sufficient to determine the entire phase behaviour. Depending on the softness, the r -n system freezes either into the fcc or bcc crystal structure, and this softness-driven transition has been extensively studied [START_REF] Hoover | Statistical Mechanics of Phase Diagrams. I. Inverse Power Potentials and the Close-Packed to Body-Centered Cubic Transition[END_REF][START_REF] Laird | Phase diagram for the inverse sixth power potential system from molecular dynamics computer simulation[END_REF][START_REF] Kofke | Solid-Fluid Coexistance for Inverse-Power Potentials[END_REF].

6 Elastic constants are macroscopic characteristics of materials which can be determined experimentally using various techniques including Brillouin scattering, ultrasonic wave propagation, and neutron scattering. From a microscopic point of view they are N-particle averages which in many cases can be obtained from molecular simulation. In the last two decades considerable progress has been made in developing reliable simulation methods for the calculation of elastic constants of many-particle systems [START_REF] Ray | Calculation of elastic constants using isothermal molecular dynamics[END_REF], [START_REF] Gusev | Fluctuation formula for elastic constants[END_REF], [START_REF] Wojciechowski | Elastic properties of twodimensional hard disks in the close-packing limit[END_REF], [START_REF] Van Workum | Improved simulation method for the calculation of the elastic constants of crystalline and amorphous systems using strain fluctuations[END_REF]. Both Molecular Dynamics (MD) and Monte Carlo (MC) calculations can be performed for specified thermodynamic conditions and in various statistical mechanical ensembles [START_REF] Ray | Statistical ensembles and molecular dynamics studies of anisotropic solids[END_REF][START_REF] Zhou | Fluctuations and thermodynamis properties of the constant shear strain ensemble[END_REF]. For example, the isoenthalpic-isotension ensemble was devised essentially to determine elastic constants by simulation [START_REF] Ray | Statistical ensembles and molecular dynamics studies of anisotropic solids[END_REF]. So the development of new simulation methods for calculating elastic properties has also had a wider impact on the development of statistical thermodynamics.

F
The various approaches for calculating elastic constants can be divided into those exploiting the defining stress-strain relations by deformation of the simulation cell, and those methods exploiting fluctuation formulas for different ensembles. For example, in the strain-strain fluctuation methods the elastic constants are obtained from the fluctuations in the shape of the simulation cell i.e., from fluctuations of the h tensor constructed from the three vectors forming the simulation cell [START_REF] Ray | Statistical ensembles and molecular dynamics studies of anisotropic solids[END_REF][START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF]. It was found, however that the strain-strain fluctuation schemes converge slowly, at least in the dynamical MD scheme used [START_REF] Sprik | Second-order elastic constants for the Lennard-Jones solid[END_REF][START_REF] Ray | Elastic constants and statistical ensembles in molecular dynamics[END_REF]. The elastic constants are more efficiently determined in the ensemble in which the h tensor and hence the strain are kept fixed and the stress tensor is allowed to fluctuate [START_REF] Ray | Calculation of elastic constants using isothermal molecular dynamics[END_REF][START_REF] Zhou | Stability criteria for homogeneously stressed materials and the calculation of elastic constants[END_REF].

For our purposes the most suitable approach is the method based on the fluctuation formulas in the canonical ensemble. This approach was invented by Squire, Holt and Hoover, [35] and has been developed and exploited in many later works (e.g., see refs., [START_REF] Hess | Shear modulus of fluids and solids[END_REF] and [START_REF] Zhou | Stability criteria for homogeneously stressed materials and the calculation of elastic constants[END_REF]). It yields isothermal elastic constants and has the advantage that all of the elastic constants can be obtained from a single simulation. Another advantage of this approach is that the reference state is not involved in the calculations. Also in the expressions for elastic constants the role of the interaction potential is more transparent than in the h-fluctuation methods, and allows for a resolution of the various contributions in the hamiltonian to the elastic properties. For a system of particles interacting via a pairwise potential, Φ(r) the elastic constants can be expressed as [START_REF] Zhou | Stability criteria for homogeneously stressed materials and the calculation of elastic constants[END_REF]:

C αβλτ = 1 V Σ i<j ∆r α ij ∆r β ij ∆r λ ij ∆r τ ij 1 r 2 (Φ -Φ /r) - 1 k B T V δ( ∆r α ij ∆r β ij Φ /r)δ( ∆r λ ij ∆r τ ij Φ /r) + N k B T V (δ αλ δ βτ + δ ατ δ βλ ), (3)
where ∆r α ij = r α i -r α j and Greek indices refer to cartesian components. The symbol δ indicates the deviation of the quantity from its mean i.e., δX = X-< X > and δ αβ is the Kronecker delta (unity when α = β, zero otherwise). V is the volume of the system of N particles, k B is Boltzmann's constant and T is the temperature. The first term in (3) is referred to as the Born-Green term and the second is the fluctuation term. The last term is the kinetic contribution, which is relatively small in solids. Thus, the elastic constants C αβλτ can be considered to be composed of the sum of two terms, C BG αβλτ + C F L αβλτ , where the kinetic contribution is included in the first, the Born-Green term.

The above formula is for an unstressed system where the components of the reference stress tensor S αβ are zero. This quantity can be computed according to

S αβ = 1 V Σ i<j ∆r α ij ∆r β ij Φ /r - N k B T V δ αβ . ( 4
)
For a system with isotropic reference stress, i.e., S αβ = -P δ αβ = 0 the appropriate expression

B αβλτ = C αβλτ + P (δ αβ δ λτ -δ αλ δ βτ -δ ατ δ βλ ). (5)
For not too soft interactions (i.e., n > 6) the inverse power solids form the fcc crystal structure, which has cubic symmetry. In what follows we restrict our discussion to a system with this symmetry. For a system of spherical particles in a cubic box and axes parallel to the coordinate axes, the following elastic constants, in the more condensed Voigt's notation (where the indices

C 11 = C 22 = C 33 , C 12 = C 13 = C 23 , C 44 = C 55 = C 66 . ( 6
)
Thus, in this symmetry only three elastic constants need to be considered. The B αβλτ have the same symmetry as the elastic constants C αβλτ , so the Voigt notation can be used and

B 11 = C 11 -P, B 12 = C 12 + P, B 44 = C 44 -P. ( 7
)
The elasticity of a system with cubic symmetry can be equivalently described by three elastic

G c = (B 11 -B 12 -2B 44 )/2. ( 8
)
The Poisson's ratio for a stretch along one of the cube axes is, [START_REF] Tokmakova | Stereographic projections of Poisson's ratio in auxetic crystals[END_REF]

ν = B 12 B 12 + B 11 = C 12 + P C 12 + C 11 . ( 9
)
Similarly as for the elastic constants, the Born-Green and the fluctuation contribution to the above moduli can be considered separately using the decomposition,

B αβ = B BG αβ + B F L αβ ,
where

B BG 11 = C BG 11 -P, B BG 12 = C BG 12 + P, B BG 44 = C BG 44 -P and B F L 11 = C F L 11 , B F L 12 = C F L 12 , B F L 44 = C F L 44 .
Knowing As a general remark, it is most convenient to keep formulae and calculated properties in a set of reduced units (uniquely appropriate for the r -n system) where we do not need to include the temperature in a formula. Temperature effects are removed by taking into account the unique scaling feature of the inverse power system. We have in these units, for length,

σ = σ(β ) -1 n , so r = rσ -1 (β ) -1 n = r/σ.
For density ρ = ρσ 3 and packing fraction ζ = π ρ/6.

For energy in k B T , we have ũ = βu and for the potential φ = βφ = r-n . For the force there In the simulations, the symmetry relations [START_REF] Baughman | Negative Poisson's ratios as common feature of cubic materials[END_REF] were exploited in the calculations of C 11 , C 12 and C 44 . It was verified that the calculated pressure P = (S xx + S yy + S zz )/3 obeys the exact relation (for the inverse power interaction), P/k B T ρ = 1 + n < u > /3, where < u > is the average potential energy per particle.

is F = F σ/k B T , and X = X -1 σ 3 (β ) 1+ 3 n = X σ3 /k B T ,
The results for the pressure, P , are shown in Fig. 1 in which P is plotted as a function of softness for nine packing fractions in the solid phase. The lowest (dashed) curve is the melting line, determined using the Gibbs-Duhem method by Agrawal and Kofke [START_REF] Kofke | Solid-Fluid Coexistance for Inverse-Power Potentials[END_REF] [38] for a number of values of the potential softness parameter. The points at ε = 0 are from the equation of state formula for the hard sphere solid [START_REF] Hall | Another hard sphere equation of state[END_REF]. From Fig. 1 it can be seen that P (ζ, ε) is not a monotonically evolving surface but has an increasingly strong ε or softness dependence at higher densities. The pressure surface contains, at any packing fraction, a maximum P max which grows in height with density. The maximum indicates that the surface separates into two regions, one where (∂P/∂ε) ζ is positive and the other (at higher ε) where it is negative.

We can characterise the two regions respectively by either a positive or negative softness 'compressibility', defined as χ = (∂P/∂ε) -1 ζ . In other words to increase pressure at a given packing fraction in the left region on the figure we have to make the potential softer, and to increase pressure in the right region it needs to be harder (i.e., by increasing n).

The maximum pressure curve, P max , marked on the figure as a bold solid line, starts on the melting line at n ≈ 19 or ε ≈ 0.053). It increases significantly for smaller ε < 0.02 (i.e., n > 50)

and goes towards ∞ for ε → 0, i.e. the hard sphere limit. Figure 2 shows P max as a function of the ε value where it occurs. This curve is well represented by the second order polynomial, In Fig. 4 examples of the elastic moduli defined in equation ( 8) are shown. The general softness-density dependence of these moduli is similar to that of the elastic constants, with a weak dependence at low densities and significant differences at higher packing fractions.

ζ(ε) = ζ 0 -4.2367ε + 12.
Furthermore there are also two regions in the ε dependence at high density. In fact the distinction is even more pronounced as now both B and G c display non-monotonic character, in having a maximum. Values of the bulk modulus can also be inferred from the ratio of the pressure and density differences. Thus, taking into account the definition B = ζ(∂P/∂ζ) and the results for the pressure surface in Fig. 1, one can expect that, for a given packing fraction, The above results suggest that there might be a non-trivial softness dependence in other elastic properties. The quantity of primary interest here, the Poisson's ratio, is shown in Fig. 5. In its softness dependence two regions can also be distinguished. For ε > 0.02, ν(ε) is nearly linear in ε and, apart from the lowest density on the figure, is almost independent of density.

Furthermore, Poisson's ratio is well represented by the static model in which the particles are fixed on the fcc lattice sites [START_REF] Tretiakov | Poisson's ratio of the fcc hard spheres at high densities. to apear[END_REF], ν static = (n + 6)/(3n + 6). For less soft systems Poissons's ratio loses its linear character and becomes strongly density dependent. For each packing fraction ν goes towards its limiting, hard sphere value, ν HS which decreases from about 1/3 close to the melting density to about 1/5 at the close packing [START_REF] Tretiakov | Poisson's ratio of the fcc hard spheres at high densities. to apear[END_REF].

ν BG = B BG 12 /(B BG 12 + B BG 11 ), (10)
which is shown in Fig. 6 for several packing fractions. From the figure it can be seen that ν BG is a fairly regular, almost linear function of ε with a slope practically equal to that of the static model line. Only for the soft systems (ε > 0.05) is a slight bend away from linearity observed. Also we note that the density dependence is considerably weaker than the softness

ν ν BG = (1 - B F L 12 B BG 12 )(1 - B F L 12 + B F L 11 B BG 12 + B BG 11 ) -1 . ( 11
)
From the above definition it follows that to lower the value of the Poisson's ratio , i.e., to get ν/ν BG < 1, in the inverse power system, the following condition has to be obeyed,

B F L 12 B BG 12 > B F L 11 B BG 11 , (12)
or

B F L 12 B F L 11 > B BG 12 B BG 11 . ( 13
)
Figure 7 shows a plot of the two ratios in equation ( 13) against ε, to explore the above criterion. The division of two fluctuation quantities unfortunately leads to some scatter in the data points. Nevertheless certain trends can be identified, and these are represented From Fig. 7 it can be seen that at lower densities, the fluctuation ratio bends toward the Born-Green ratio for ε → 0 (see the data for ζ = 0.58). This indicates that for some range of the softness and packing fractions close to the melting the Poisson's ratio is lower than that of the hard spheres. That ν < ν HS (and the possible presence of a small minimum in Fig. 5 for

ζ = 0.58
) is an unexpected result. The lowest value of ν for the lower packing fractions close to the melting value is therefore not the hard sphere value. 17 Poisson's ratio is predominately linear and hardly modified by density. In the second region of larger interaction steepness (smaller ε), considerable variation in ν with density and softness is evident. From the point of view of the auxetic behaviour this may be an interesting ν but also may causes an increase of its value.

To obtain the desired effect, (i.e., a decrease of Poisson's ratio value) augmentation of the fluctuation part should take place in an appropriate range of softness and state parameters. It is not sufficient to just make the system harder (i.e., decrease ε) to get a noticeable lowering of Poisson's ratio (see the data for ζ = 0.58 in Fig. 5), but it is also necessary to increase the packing fraction. For the inverse power solid the appropriate range can be considered to be in [START_REF] Agrawal | Phase diagram for the inverse sixth power potential system from moleculoar dynamics computer simulation[END_REF]. The edge (ε = 0) points are taken from the hard sphere equation of state by Hall [START_REF] Hall | Another hard sphere equation of state[END_REF]. The bold line is the demarcation or maximum pressure line. 3) and ( 6) against the interaction softness ε for the fcc solid of the inverse power system, and for two packing fractions (indicated on the figure).

  Elastic properties constitute one of the most fundamental physical properties of any material.

  ratio, ν, is always the same,

  moduli: the bulk modulus B, the shear modulus G and the cubic modulus G c[START_REF] Hess | Shear modulus of fluids and solids[END_REF]. They are related to the elastic constants in the following way,B = (B 11 + 2B 12 )/3, G = (B 11 -B 12 + 3B 44 )/5,

  where X stands for pressure, elastic constant or modulus. Time is in t = t(k B T /mσ 2 ) 1/2 = t(T * /mσ 2 ) 1/2 . In what follows we omit the tilde, so keep in mind that we are dealing with temperature-scaled quantities (also for k B T / = 1 the conventional reduced units are recovered).Molecular dynamics simulations were performed for particles interacting via the inverse power potential (2) for a range of densities in the solid phase, ζ ≥ 0.56 and several values of n ranging between 14 to 144. The MD simulations were carried out on N = 256 particle systems with periodic boundary conditions and a standard interaction cut-off of 2.5, and 2.0 for n > 70. The canonical Nosé-Hoover equations of motion[START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] were solved using the 4-th order Runge-Kutta algorithm, with a time step, ∆t = 0.001 or 0.0005 (for n > 70). The reduced temperature was T = 1. The data were calculated from averages over 0.5 -1 million time step runs after equilibration and the initial positions of the particles were on fcc lattice sites. The statistical convergence of the elastic moduli by the fluctuation method is quite slow compared with the usual thermodynamic averages (e.g., energy and pressure). These runs are long enough to obtain elastic properties with sufficient accuracy to observe trends on the ζ, ε surface.

1813ε 2 ,

 2 where ζ 0 = 0.7405 ∼ = π √ 2/6 is the value of the close-packed 13 the hard sphere fcc structure, ζ HS cp[START_REF] Hall | Another hard sphere equation of state[END_REF]. Thus, for any ζ, a maximum on the P (ζ, ε) surface exists at some value of the softness and its magnitude tends to infinity in the hard sphere limit, i.e., ζ = ζ HS cp at the boundary ε → 0.In Fig.3the softness dependence of the elastic constants is shown. At the lower density (ζ = 0.58), the ε-dependence of all three quantities is relatively weak and they are almost flat lines, with the relative order C 11 > C 44 > C 12 > 0, indicating a weak ε dependence. This behaviour changes significantly on increasing density to ζ = 0.66, where the ε-dependence become stronger and more importantly its form is different for the three elastic constants. Additionally, two regions become visible, depending on ε . In the first, representing the harder systems for which ε is less than 0.02, the values of C 11 and C 44 are relatively insensitive to ε while the C 12 increases with softness. For softer systems, with ε > 0.02, all three quantities decrease significantly and monotonically with increasing softness. In this second region the form of the softness dependence is similar for the three elastic constants and can be approximated well by an exponential. The emergence of a well-defined maximum in the softness dependence of the C 12 is a new observation.

  B and P have similar ε values.

  Figure 5 also shows that the static limiting case (bold line on the figure) is a relatively good description of the Poisson's ratio, at least for some range of the ζ, ε. The static model means that the fluctuation contributions are not taken into account, and that for the Poisson's ratio formula, only the Born-Green terms are retained. This approximation gives a 'Born-Green Poisson's ratio',

  dependence. A variation of ζ over the entire solid phase causes only a minor change (typically a few percent) in the ν BG value for any ε. As can be seen in the figure, on increasing ζ, the nearly colinear curves approach from below the static model line. Thus, the static result can 15 starting point for quantitative description of ν BG , and it can be seen that significant changes in Poisson's ratio are achieved by varying ε rather than ζ in the available parameter range. A comparison between Figs. 5 and 6 indicates that changes in ν(ζ, ε) are mainly determined by the fluctuation terms and more importantly they can cause a considerable decrease as well as increase of the Poisson's ratio. Their role becomes more transparent after writing Poisson's ratio in a form explicitly featuring ν BG ,

  figure. Hence, respectively ν > ν BG , ν = ν BG and ν < ν BG , can be satisfied at any given density.

  In this work the influence of interparticle interaction softness on the solid phase elastic properties has been investigated, particularly for Poisson's ratio, ν. The canonical ensemble MD simulation method of Nosé and Hoover was used and the stress fluctuation scheme was applied to obtain the isothermal elastic constants of the fcc phase. The elastic constants and elastic moduli were determined for a range of the softness parameter (ε ≡ 1/n) of the soft-sphere potential and covering most of the solid phase. This data was used to construct density-softness surfaces of these properties. The determination of the elastic constants from the simulations is time consuming and so the results to date only indicate trends. More quantitative conclusions will have to wait for more extensive simulations to be carried out. One characteristic feature which emerges from these studies is the existence of a maximum in the pressure and the elastic constant C 12 on a plot of these quantities against the particle softness (ε) for constant density (ζ) systems. Locating the elastic constant maxima in terms of the value of softness for each density (as we have done for the pressure) requires more extensive simulations.Two other moduli, B and G c defined in Section 2 also exhibit maxima. In all cases the height of the maximum increases with increasing density. In the case of the pressure the maximum line starts at the melting density for n ≈ 19 and terminates at close packing of the hard sphere fcc structure (ζ = 0.7405). For all calculated elastic properties it is possible to distinguish two regions in their density-softness dependence. One, roughly for ε > 0.02 where a weak ζ-dependence and fairly regular ε-dependence is present. In this region the

FIG. 2 :FIG. 3 :

 23 FIG.2:The softness and packing fractions for which the maximum pressure line in Fig.1is the locus.

FIG. 4 :FIG. 5 :[FIG. 6 : 11 FLFIG. 7 :

 456117 FIG.4: The bulk, shear and cubic moduli defined in equation (8) of the inverse power solid against

  In this work we concentrate on thermodynamic properties, elastic constants and moduli of the fcc solid phase of the inverse power potential. We investigate the basic thermodynamic and mechanical quantities to establish how the interaction softness, ε, influences the elastic moduli and Poisson's ratio.Methods for calculation of elastic constants are briefly discussed in Sec. II. Simulation results are presented in Sec. III and concluding remarks are made in Sec. IV.

	o r F o r
	P P
	e e
	e r e r
	R R
	e e
	v i e v i e
	w w
	O n l O n l
	y y
	7

II. CALCULATION OF ELASTIC CONSTANTS

  1, 2, 3, 4, 5, 6 are used for xx, yy, zz, xy, xz, yz, respectively) are not zero,

  B BG αβ and B F l αβ we also have B, G, G c as a sum of Born-Green and fluctuation terms, e.g., B = B BG + B F L where B F L = (B F L 11 + 2B F L 12 )/3, and B BG = (B BG 11 + 2B BG 12 )/3. Obviouslythe Poisson's ratio ν is not ν BG + ν F L , but it is nevertheless informative to define and inspect the quantities ν BG and ν F L , as we discuss in the next section.

	III. SIMULATION RESULTS
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	11

  region as considerable lowering of the value of the Poisson's ratio can be achieved by slight modification of the interaction softness, particularly at higher packing fraction. Therefore for steep interaction particle systems, small variations in the density and softness can have a significant affect on the Poisson's ratio. This possibility seems not to have been fully exploited yet.It was demonstrated that the Born-Green Poisson's ratio, ν BG , is well represented by the static case and is almost density independent, suggesting that the fluctuation part plays an important role in determining Poisson's ratio. It was found that it is the ratio of the fluctuation part of the elastic constants, B F L 12 /B F L 11 that is responsible for significant changes both in magnitude and form of ν(ζ, ε). We have found that for the inverse power potential fcc solid, ν > ν BG

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	18

for such ζ, ε for which the ratio B F L 12 /B F L 11 < B BG 12 /B BG 11 . Therefore an important, potentially useful outcome of this work is the observation that the fluctuation part can not only diminish

  the region ca. n > 50 and ζ ≥ 0.66 (it worth to add that in this region the relative contributions of the fluctuation part to the Poisson's ratio becomes greater than that of the Born-Green part or ν F L /ν BG > 1). It is perhaps not surprising that the fluctuation term is so important in potentially lowering Poisson's ratio, as it is the cooperative motion of the building components that has been found responsible for giving auxeticity in most the known examples.

	40 60 80 100 pressure P	F o r F o r F o r F o r F o r F o r	ζ 0.70 0.69 0.68 0.66 0.64 0.62 0.60 0.58 0.56
			P P P P P P
	20 0	0.02	e e r R e e r R e e r R e e r R e e r R e e r R 0.04 ε	0.06
			e e e e e e
			v i e v i e v i e v i e v i e v i e
				w w w w w w
				O n l O n l O n l O n l O n l O n l
				y y y y y y
			19

FIG. 1: The pressure, P (ζ, ε) -surface, of the inverse-power solid. The lowest line is the melting curve determined by Agrawal and Kofke

http://mc.manuscriptcentral.com/tandf/jenmol

Acknowledgements

The authors thank The Royal Society (London) and the Polish Academy of Sciences for partly funding this collaboration. We are grateful to K.V. Tretiakov and K.W. Wojciechowski for providing computer simulation data for the hard sphere Poisson's ratio prior to publication. The work has been partially supported by the Polish Committee for Scientific Research (KBN) grant No. 4T11F01023. 20