V N Salomatov
email: sav@irgups.ru

Keywords: Relativistic particle, Helmholtz equation, Principle of correspondence, Discreteness of the rest mass spectrum

published or not. The documents may come

Klein-Gordon and Schrödinger equations for a free particle in the rest frame

Introduction

The main successes of quantum mechanics in the quantitative description of nonrelativistic systems are connected with the Schrödinger equation. Often this equation is regarded as the nonrelativistic limit of the Klein-Gordon equation [1][2][3][START_REF]Dong X[END_REF][5][START_REF] Pasquali | Annali di Matematica Pura ed Applicata[END_REF]. When considering this limit, as a rule, in the solutions of the Klein-Gordon equation, the speed of light tends to infinity. In [START_REF] Wharton | [END_REF], the relationship of the Klein-Gordon and Schrödinger equations is analyzed in detail and it is concluded that the Klein-Gordon equation cannot be reduced to the Schrödinger equation in any limit. We can agree with the author of [START_REF] Wharton | [END_REF], if we keep in mind an arbitrary reference frame. But it must be noted that the Klein-Gordon equation must strictly go over into the Schrödinger equation only at a particle velocity of zero, that is, in the particle's rest frame.

The objective of this paper is to find a solution joint to the Klein-Gordon and Schrödinger equations in the rest frame for a neutral spinless particle.

The rest frame in relativistic quantum mechanics

In classical relativistic mechanics, the particle's rest frame is a frame of reference where the particle velocity v is equal to zero. In relativistic quantum mechanics, its own reference frame was used in [8][START_REF] Edmonds | Quantum Uncertainties[END_REF][10][11][12]. In quantum mechanics, as in quantum field theory, there are problems with determining the eigenvalues of the velocity operator [13][14][15]. In various formulations of quantum theory, it is believed that the particle velocity has the same direction as the wave vector k. In any case, strictly speaking, the situation with the speed of a single particle does not satisfy the principle of correspondence between classical and quantum mechanics. There are formulations of quantum mechanics (primarily the De Broglie-Bohm approaches [16][17][18][19]), where the particle velocity is interpreted as the group velocity of the wave packet

𝐯 = 1  • 𝜕E 𝜕𝐤 . ( 1
)
When performing the relativistic dispersion relation

E = 𝑐  ( 𝑚 2 𝑐 2  2 + 𝐤 2 ) 1/2 (2)
particle velocity in the relativistic version of the De Broglie-Bohm approach is written as

𝐯 = 𝑐k( 𝑚 2 𝑐 2  2 + 𝐤 2 ) -1/2 . ( 3
)
Here,  is the Planck constant, c is the speed of light in vacuum, m is the rest mass of the particle, E = E (k) is the dispertion relation.

That is, v → 0 in the reference frame, where k → 0. Therefore, when passing from an arbitrary reference frame to its own reference frame, all the characteristics of the particle, including the wave function, do not contain the wave vector k in their record. If we find the wave function in the rest frame (not containing the vector k), then in other reference systems this function contains the vector k in its notation, which allows this function to be transformed according to irreducible representations of the Lorentz group. Such a situation should be fulfilled in other approaches when fulfilling the principle of correspondence between classical and quantum mechanics.

Joint solution of the Klein-Gordon and Schrödinger equations in the rest frame

Let's look at a system of two equations

𝑖  𝜕ᴪ 𝜕𝑡 = -  2 2𝑚 ∇ 2 ᴪ , (4)
𝜕 2 ᴪ 𝜕𝑡 2 + 2𝑖𝜔 𝜕ᴪ 𝜕𝑡 + 𝜔 2 = 0 . ( 5
)
Here 𝜔=𝑚c 2 /  , ∇ 2 is the Laplace operator.

Equation ( 4) is the Schrödinger equation for a free particle. By a free particle we mean a particle that does not interact with other particles, but interacts with a vacuum. Equation ( 5) is chosen so that when solving the system of equations (4, 5) by the substitution method, the Klein-Gordon equation

 2 𝑐 2 𝜕 2 ᴪ 𝜕𝑡 2 =  2 ∇ 2 ᴪ -𝑚 2 𝑐 2 ᴪ (6)
is fulfilled. Indeed, from ( 4)

𝜕𝑡 = 𝑖  2𝑚 ∇ 2 ᴪ . ( 7
)
Substituting ( 7) into (5), we obtain equation [START_REF] Pasquali | Annali di Matematica Pura ed Applicata[END_REF]. Thus, the solutions of the system of equations (4.5) are also solutions of the equation ( 6) (but not vice versa).

The solution of the system of equations (4.5) we begin with the solution of equation (5). Characteristic equation

𝑝 2 + 2𝑖𝜔𝑝 + 𝑝 2 = 0 (8)
has two roots

𝑝 1 =-(1-√2)iω 𝑝 2 =-(1+√2)iω . ( 9
)
Therefore, the general solution of equation ( 5) can be written as

ᴪ= 𝐶 1 (𝑟 ⃗)𝑒 𝑃 1 𝑡 + 𝐶 2 (𝑟 ⃗)𝑒 𝑃 2 𝑡 = 𝐶 1 (𝑟 ⃗)𝑒 -(1-√2)𝑖𝜔𝑡 + 𝐶 2 (𝑟 ⃗)𝑒 -(1+√2)𝑖𝜔𝑡 . ( 10
)
The coefficients C 1 (r) and C 2 (r) are generally complex and depend on spatial coordinates.

Solution (10) can also be written as

ᴪ=𝑒 -𝑖𝜔𝑡 [(𝐶 1 (𝑟 ⃗)𝑒 √2𝑖𝜔𝑡 + 𝐶 2 (𝑟 ⃗)𝑒 -√2𝑖𝜔𝑡 )] . ( 11
)
To determine C 1 (r) and C 2 (r) we substitute (10,11) into the Schrödinger equation ( 4). Then we get

𝐶 1 (𝑟 ⃗)(1 -√2)𝜔𝑒 -(1-√2)𝑖𝜔𝑡 + 𝐶 2 (𝑟 ⃗)(1 + √2)𝜔𝑒 -(1+√2)𝑖𝜔𝑡 = = -  2𝑚 ∇ 2 𝐶 1 (𝑟 ⃗)𝑒 -(1-√2)𝑖𝜔𝑡 -  2𝑚 ∇ 2 𝐶 2 (𝑟 ⃗)𝑒 -(1+√2)𝑖𝜔𝑡 . ( 12
)
This shows that (10,11) are solutions of equation ( 4) under the conditions

∇ 2 𝐶 1 (𝑟 ⃗) = 𝑄 1 2 • 𝐶 1 (𝑟 ⃗) (13)
∇ 2 𝐶 2 (𝑟 ⃗) = 𝑄 2 2 • 𝐶 2 (𝑟 ⃗) . ( 14
)
Here

𝑄 1 2 = 2(√2 -1) • 𝑚 2 𝑐 2  2 (15)
𝑄 2 2 = -2(√2 + 1) • 𝑚 2 𝑐 2  2 . ( 16
)

Discreteness of the rest masses of relativistic spinless particles as a consequence of the joint solution of the Klein-Gordon and Schrödinger equations in the rest frame

Thus, the solutions (10,11) are the same for the system of two equations, which includes the Schrödinger equation, and for the Klein-Gordon equation in the rest frame when (13,14) is satisfied. Conditions (13,14) are the Helmholtz equations for C 1 (r) and C 2 (r). Moreover, equation ( 14) has discrete solutions. These solutions are determined by the symmetry of the function C 2 (r) with an appropriate choice of the coordinate system (for example, spherical or cylindrical) and the shape of the boundary conditions [START_REF] Abramowitz | Handbook of Mathematical Physics[END_REF][START_REF] Korenev | Bessel Functions and Their Applications[END_REF][START_REF] Bayn | Essentials of Mathematical Methods in Science and Engineering[END_REF][START_REF] Salomatov | [END_REF][24]. In [START_REF] Salomatov | [END_REF]24], the simplest polar model of a particle in the rest frame was considered. The Helmholtz equation similar to (14) was obtained under the assumption that a particle that does not interact with other particles can be described by the simplest wave equation. In this case, a part depending on spatial coordinates was allocated in the solution of the equation. In the polar model, it was divided into radial and azimuthal components

          R q . ( 17
)
After the transformation q𝜌=x the equation for R(𝜌) was reduced to the Bessel equation for the function

F(x)=R(x/q). ( 18
)
When boundary conditions

  , 0  a R (19)      R under 0   (
         . x J B x N B x F n B nl n N nl nl   (21)
Here, J n (x), N n (x) are n-order Bessel and Neumann functions, B nl (N) , B nl (B) are expansion coefficients. In this case, the allowed discrete values of the rest masses were given by the expression

nl nl X c a m  1   . (22)
Here, X nl are the values of the roots of the  

x F nl functions in [START_REF] Korenev | Bessel Functions and Their Applications[END_REF].

In the present paper, the Helmholtz equation ( 14) is a consequence of the joint solution of the system of two equations and the Klein-Gordon equation and differs from the Helmholtz equation in [START_REF] Salomatov | [END_REF]24] by the factor 2(1+√2). In this regard, the discrete values of the rest masses allowed by equation ( 14) in the polar model have the form

𝑚 𝑛𝑒 = ± 1 √2(1+√2) • 1 𝑎 •  𝑐 𝑋 𝑛𝑒 . ( 23
)
The value of 𝑎 in (19,[START_REF] Bayn | Essentials of Mathematical Methods in Science and Engineering[END_REF][START_REF] Salomatov | [END_REF], generally speaking, is arbitrary. In [START_REF] Salomatov | [END_REF]24], it was chosen so that formula [START_REF] Bayn | Essentials of Mathematical Methods in Science and Engineering[END_REF] gives the value of the electron mass at X nl =X 01 (N)

. If now in formula [START_REF] Salomatov | [END_REF] the quantity 𝑎 √ 2(1 + √2) is chosen so that [START_REF] Salomatov | [END_REF] gives the value of the rest mass of an electron at X nl =X 01

(N)
, then all the conclusions and results of [START_REF] Salomatov | [END_REF]24] are valid and for formula [START_REF] Salomatov | [END_REF] in the polar model. In particular, the peculiarities of the values of the rest masses and lifetimes of the electron, proton, and neutron [START_REF] Salomatov | [END_REF], as well as the results for the current-quark masses of the first generation quarks [24] are valid. At the same time, it is of interest to study the features of the mass spectrum of elementary particles with the direct use of formula [START_REF] Salomatov | [END_REF], especially for other particle models, including the spherical model.

Discussions and conclusion

The system of equations (4, 5) determines the structure of the solution (10,11) of the Klein-Gordon equation ( 6), as well as additional conditions (13)(14)(15)(16) for the spatial components C 1 (r) and C 2 (r). Concretization of the solution of equation ( 6)

when (13)(14)(15)(16) are fulfilled should consist in the choice of a particle model (only the polar model was used in [START_REF] Salomatov | [END_REF]24]), as well as in the specification of boundary conditions of type (19,[START_REF] Abramowitz | Handbook of Mathematical Physics[END_REF] when solving the Helmholtz equations (13,14). While equation ( 14) can be considered partially investigated (it leads to discreteness of rest masses), equation ( 13) provides additional opportunities for studying the nature of microparticles. Thus, equations (4,5) can be considered as additional conditions for solving the Klein-Gordon equation in the rest frame and allow extracting additional information about microparticles, in particular, explaining the discreteness of the rest mass spectrum in the framework of relativistic quantum mechanics.

Note that equation ( 5) along with solutions of type (10,11) has "oscillating" solutions that differ from solutions of type (10.11). These solutions are not compatible with the Schrödinger equation ( 4). They can be of independent interest for the study of microparticles in the rest frame. It is also of interest to use the solutions (10,11) when creating operators in the framework of quantum field theory.

  20) are met the solutions of the Bessel equation for a   x F nl function are linear combinations of the Bessel and Neumann functions


Filippo Geraci
email: filippo.geraci@iit.cnr.it

Riccardo Massidda
email: r.massidda@studenti.unipi.it

Nadia Pisanti
email: pisanti@di.unipi.it

Technology and Species independent Simulation of Sequencing data and Genomic Variants

Keywords: NGS sequencing, simulation, genomic variants

Highly accurate genotyping is essential for genomic projects aimed at understanding the etiology of diseases as well as for routinary screening of patients. For this reason, genotyping software packages are subject to a strict validation process that requires a large amount of sequencing data endowed with accurate genotype information. In-vitro assessment of genotyping is a long, complex and expensive activity that also depends on the specific variation and locus, and thus it cannot really be used for validation of in-silico genotyping algorithms. In this scenario, sequencing simulation has emerged as a practical alternative. Simulators must be able to keep up with the continuous improvement of different sequencing technologies producing datasets as much indistinguishable from real ones as possible. Moreover, they must be able to mimic as many types of genomic variant as possible.

In this paper we describe OmniSim: a simulator whose ultimate goal is that of being suitable in all the possible applicative scenarios. In order to fulfill this goal, OmniSim uses an abstract model where variations are read from a .vcf file and mapped into edit operations (insertion, deletion, substitution) on the reference genome. Technological parameters (e.g. error distributions, read length and per-base quality) are learned from real data. As a result of the combination of our abstract model and parameter learning module, OmniSim is able to output data in all aspects similar to that produced in a real sequencing experiment.

I. INTRODUCTION

DNA sequencing is the task of determinating the ordered sequence of nucleotides of a biological sample. In spite of this very simple definition, the last fifty years have witnessed one of the biggest collaborative efforts to improve and refine sequencing technologies [START_REF] Heather | The sequence of sequencers: The history of sequencing dna[END_REF]. Increasing output quality and throughput has inspired researchers eager for understanding the mechanisms of life, and has fostered the development of astonishing clinical applications. However, to the fast enhancement of technology (in particular the so-called next generation sequencing [START_REF] Goodwin | Coming of age: ten years of next-generation sequencing technologies[END_REF]) did not always follow an equally fast development of algorithms and methods for genomic data analysis. In order to assess variant detection algorithms, designers need large collections of accurate genotype annotations, as well as the actual genomic sequences. Producing such essential supplementary information is still slow and expensive. In this scenario, simulation comes into play, promising a cheap and easy way to carry experiments on a large scale by producing genomic sequences endowed with genotypes. The drawback of using simulators is the need to be confident that results on real data would be comparable with those achieved on simulated ones.

Three generations [START_REF] Magi | The source of the data flood: Sequencing technologies[END_REF] of continuously evolving sequencing technologies (Sanger, high throughput short reads and long reads), several types of genomic variations (from single nucleotide polymorphisms to structural variants), and different error profiles, actually make the development and validation of accurate simulation pipelines a challenging task.

First generation (Sanger) sequencing produces reads slightly smaller than 1kbp that are merged in larger contigs by means of the shotgun sequencing technique [START_REF] Anderson | Shotgun dna sequencing using cloned dnase i-generated fragments[END_REF]. The low throughput of this technology makes it not suitable for large scale projects but enables very accurate genotypization because of its low error rate.

The market of the second generation of sequencing (sequencing by synthesis) has been dominated by Illumina (with many notable and valuable exceptions). In most applications, reads are produced in pairs, sequencing the two endpoints of long fragments. The read length is much smaller than that of the first generation, and started from about 35bps of the early stage of this technology to one or two hundred bps of its maturity. For this technology, substitutions are the most common type of error and tend to become more frequent in the last reading cycles and the higher error rate is compensated with coverage exploiting the massively higher throughput and (consequent) much lower costs.

The third generation of sequencing technologies has further revolutionized the scenario producing impressively long reads which span a large length range from few kilobases up to more than a hundred of kilobases. Errors are typically more frequent than those of the previous generations, and can more likely include insertions and deletions, making read mapping particularly challenging.

Besides technology-dependent constraints, there are other factors that influence the development of simulators. In particular, handling the largest possible set of variant types would be necessary to make the output more realistic and usable. However, as shown in table 4 of [START_REF] Escalona | A comparison of tools for the simulation of genomic next-generation sequencing data[END_REF], the majority of simulators is limited to few common types or handle only some specific one. There are two alternative approaches for variant generation. The most common method is that of embedding variations at random positions of the reference sequence controlling, feeding the tool with just a mutation rate parameter. Often, the spatial distribution across the genome and the type of variant is learned from a real dataset. The other complementary approach generates variants according to a list of types and genomic locations taken as input.

Secondary features of simulators include: the ability to generate datasets for metagenomics, the possibility of introducing PCR-based noise, and the generation of already aligned data.

In this paper we describe OmniSim: a novel simulator thatsimultaneously -fulfills the promise (i) to be independent from the sequencing technology, and (ii) to handle all the possible types of genomic variant. The technological parameters (e.g. errors and lengths distributions) for our tool are learned by a specialized module that analyzes a sequencing archive in .bam format and extracts information about: read length distribution, error types and probability, quality scores, and PCR noise. Genomic variants are generated according to a user-provided annotation .vcf file. We map each variant into a series of edit operations on the reference genome, so that our simulator can handle every possible type of variant. We also pursue secondary features such as (i) simulating the effect of PCRbased stutter noise in repeated sequences, (ii) enabling variants on different reference genomes to be mixed so as to simulate metagenomics data, and (iii) returning aligned data so as to speed up subsequent experiments.

Experiments, using the paired-end 100bp Illumina platinum genome [START_REF] Eberle | A reference data set of 5.4 million phased human variants validated by genetic inheritance from sequencing a three-generation 17-member pedigree[END_REF] of the NA12877 human sample and the Oxford Nanopore sequencing [START_REF] Jain | Nanopore sequencing and assembly of a human genome with ultra-long reads[END_REF] of the NA12878 human sample, both from the CEPH pedigree 1463, have shown that OmniSim is reasonably fast and completes the entire pipeline (i.e. generation of two human alleles including all the variations reported in dbSNP [START_REF] Sherry | dbsnp: the ncbi database of genetic variation[END_REF], learning of the error model, and generation of 100 million reads) in about 2/3 hours. Experiments have also shown that the nucleotide composition and the distribution of per-base quality scores of the simulated reads perfectly reflects that of a real sequencing. The source code of OmniSim is freely available at the URL: https://gitlab.com/geraci/omnisim II. RELATED WORK Introducing their work, in [START_REF] Johnson | A better sequence-read simulator program for metagenomics[END_REF] the authors complain that simulators are usually designed for specific needs, and are not general enough to be employed in different applicative scenarios. This complaint is somehow confirmed in [START_REF] Escalona | A comparison of tools for the simulation of genomic next-generation sequencing data[END_REF] where the authors provide a decision tree with the purpose of guiding users to an aware choice of the simulator that better fits their needs.

The main aspects on which simulators differ are: the sequencing technologies they can reproduce, the type of variants introduced, the adopted error model, and the richness of the output.

Short reads sequencing is the most represented technology with probably dozens of specialized tools. Among them, we cite GemSIM [START_REF] Mcelroy | Gemsim: general, errormodel based simulator of next-generation sequencing data[END_REF], IntSIM [START_REF] Yuan | Intsim: An integrated simulator of next-generation sequencing data[END_REF], SInC [START_REF] Pattnaik | Sinc: an accurate and fast error-model based simulator for snps, indels and cnvs coupled with a read generator for short-read sequence data[END_REF], pIRS [START_REF] Hu | pirs: Profile-based illumina pair-end reads simulator[END_REF], and ART [START_REF] Huang | Art: a next-generation sequencing read simulator[END_REF]. Short read simulation tools are usually vendor independent and differ on other features like the error model or the type of variants. In this category, Grinder [START_REF] Angly | Grinder: a versatile amplicon and shotgun sequence simulator[END_REF] has the unique feature to be able to accurately simulate shotgun and amplicon sequencing. The tools able to simulate long reads, instead, are usually specialized in either Pacific Bioscience SMRT sequencing or Oxford Nanopore sequencing with the notable exception of SiLiCO [START_REF] Baker | Silico: a simulator of long read sequencing in pacbio and oxford nanopore[END_REF] that can simulate both of them. Simulators tailored on Pacific Bioscience SMRT sequencing are more common. In particular, the SimLoRD [START_REF] Stöcker | Simlord: simulation of long read data[END_REF] can only simulate circular consensus sequence reads (CCS) while PaSS [START_REF] Zhang | Pass: a sequencing simulator for pacbio sequencing[END_REF] and NPBSS [START_REF] Wei | Npbss: a new pacbio sequencing simulator for generating the continuous long reads with an empirical model[END_REF] can also produce continuous long reads (CLR). Due to its complexity and rapid evolution, the Oxford Nanopore Sequencing simulation is more challenging and hence less common. In this case, simulators do not produce the sequence of nucleotides but the raw signal emitted by the sequencer. The interesting rationale of this choice is that of opening to the development and testing of base callers. Among the most valuable methods we mention [START_REF] Li | Deepsimulator: a deep simulator for nanopore sequencing[END_REF] and [START_REF] Rohrandt | Nanopore simulation-a raw data simulator for nanopore sequencing[END_REF]. It worth citing FASTQSim [START_REF] Shcherbina | Fastqsim: platform-independent data characterization and in silico read generation for ngs datasets[END_REF] which is the only attempt to simulate NGS data both for short and long reads (at least PacBio SMRT). Like for our method, this is made possible by learning the parameters from real datasets. However, as shown in [START_REF] Stöcker | Simlord: simulation of long read data[END_REF] figure 1, FASTQSim the error rate does not agree well with the raw data.

Handling metagenomics has also emerged as a useful feature of simulators. This task is complicated by the fact that reads come from several organisms, and hence, the generation needs to be able to deal with large data. Specific tools have been implemented for this task: MetaSim [START_REF] Richter | Metasima sequencing simulator for genomics and metagenomics[END_REF] is probably the first attempt in this sense; it can simulate either Sanger sequencing or sequencing by synthesis. The tool BEAR [START_REF] Johnson | A better sequence-read simulator program for metagenomics[END_REF] focuses on simulating a realistic abundance of samples without the need of learning it from real data; GemSIM [START_REF] Mcelroy | Gemsim: general, errormodel based simulator of next-generation sequencing data[END_REF] implements an accurate statistical model to handle sequencing errors but it only simulates single nucleotide polymorphisms (SNPs) variants. In our case, as we will explain in detail in section III, metagenomics simulation is achieved allowing alleles generated from different organisms to be mixed before reads simulation.

Mostly, simulators produce their output in fasta or fastq format, thus losing the alignment information. This has the disadvantage of leaving to the user the burden of alignment. An interesting exception is [START_REF] Earl | Assemblathon 1: a competitive assessment of de novo short read assembly methods[END_REF] that optionally returns also the aligned reads using the standard .bam format. Moreover, producing fastq does not necessarily mean producing accurate per-base qualities. Grinder [START_REF] Angly | Grinder: a versatile amplicon and shotgun sequence simulator[END_REF], for example, produce qualities with a fixed value. OmniSim has an option to keep track of the genomic position of the generated reads so as to produce the alignment in .bam format. Moreover, the per-base quality is first learnt from real data and then produced in accordance with the simulated sequencing errors.

III. METHOD

Our simulator consists of the three modules depicted in Figure 1: 1) A Genotype generator that takes in input a reference genome and a list of variants, and produces a customizable number of alleles; 2) A Learning module that extracts from a .bam file the technological parameters (e.g. distributions of lengths, errors, etc.), and 3) The Simulation module that actually generates the sequences from a set of alleles (e.g. coming from the first module) according to a given profile (which can be, for example, the outcome of the learning module).

The three modules are designed to be as independent from each other as possible, in order to be possibly employed for diverse usage scenarios. For example, developing de-novo genome assemblers, it is possible to simulate the sequencing of the same individual with different technologies by generating the reads from the same set of alleles but changing the technological parameters.

A further possible advantage of having the reads production independent from the genotype generator, is that the genotypes of different species can be mixed for metagenomics simulations purposes.

A. Genotype generation

Current simulators have two main drawbacks: they do not support all the possible types of genomic mutation at once (see table 4 in [START_REF] Escalona | A comparison of tools for the simulation of genomic next-generation sequencing data[END_REF]), and variations are spread across the genome according to predetermined statistical distributions that not necessarily reflect the real distribution for the species of interest. The first issue causes both the generation of unrealistic reads (where only a subset of variations is present) and the need to perform separate validations using a suitable simulator per each type of variation the user is interested in.

The second issue can become relevant when testing algorithms for which the distribution of the variations can affect the output. For example, as explained in [START_REF] Edge | Hapcut2: robust and accurate haplotype assembly for diverse sequencing technologies[END_REF], in the single individual haplotyping problem (SIH) the distance between consecutive SNPs can let some algorithms fail to reconstruct the correct haplotype.

The goal of the genotype generation module is that of overcoming both limitations. We started from the practical evidence that, for the large majority of the most commonly studied organisms, a reasonably complete map of the genomic variations with their type, estimated frequency, and exact position, is available in public databases like dbSNP [START_REF] Sherry | dbsnp: the ncbi database of genetic variation[END_REF]. For these organisms, maps can be used, rather than leveraging on a statistical model to decide the position, type and content of variations.

We further observed that all genomic variants are alterations of the DNA sequence that can be modeled as a series of unary edit operations (insertion, deletion, replacement). Therefore, variants types can be mapped into a suitable sequence of edit operations whose final outcome corresponds to them, and hence their application can be simulated by applying those.

The application of a specific variant on an allele is not always guaranteed because of uncontrollable aspects such as the allelic frequency and the previous application of another conflicting variant. Conflicts arise when two variants insist on the same locus. A typical example of this situation is that of a SNP taking place in the same DNA portion of -saya deletion: the latter let the nucleotide involved in the SNP disappear. Another case is when the same locus undergoes both an insertion and a deletion: although these variations could theoretically co-exist, their concurrent application would probably not make sense, and thus only one variation should be considered. According to different applicative scenarios, users may -or may not -want to leave to the faith the decision of which of the alternative variations to apply (either forcing a specific one to be selected regardless of the others or making a random choice). Forcing a specific variant is the case, for example, of [START_REF] Tang | Profiling of short-tandem-repeat disease alleles in 12,632 human whole genomes[END_REF] where the same tandem repeat is increasingly extended in order to evaluate the longest genotypable region.

Our genotype generation module takes as input a reference genome, a map (as a .vcf file) of all the annotated variants, and a list of user-defined variants, and returns a (customizable) number of allele sequences with the corresponding alignment to the reference genome. The allele' generation procedure iterates over the chromosomes copying the reference sequence until the position of a genomic variant is reached. The variation is applied to a randomly selected allele according to its observed frequency. Possible subsequent mutations on the same locus are applied to other alleles (again, according to their observed frequencies) until when all the alleles have been generated. If no variants are selected for an allele, then the reference sequence is just copied. Before applying a mutation to an allele, the procedure checks whether there is a conflict with previously introduced variants (this check is done by verifying that the allele's sequence matches the reference in the positions involved in the new mutation). In order to force user-defined variants, these mutations are set to omnipresent in the population (namely, observed with probability 1), and are evaluated before the others (so that they cannot be withdrawn because of conflicts). Once a mutation is introduced into an allele, the corresponding alignment sequence is updated.

Although fast, simulating a population with different genotypes, the above procedure could need to be run a large number of times. In this case, exploiting the CPU parallelism to speed up the alleles' generation would be useful. We achieve this goal by allowing the parallel generation of distinct chromosomes.

B. Parameter learning

Sequencing platforms produce reads with very different characteristics in terms of read length distribution, error types, and coverage. All these parameters need to be learned in order to accurately simulate the desired sequencing procedure. Interestingly, none of them depends on the organism under study, and thus it is possible to simulate sequencing of any sample regardless of the organism used to learn parameters.

In order to collect statistics useful to mimic the behavior of a given technology, the learning module leverages on a reference genome and a set of aligned reads in .bam format. Both single-end and paired-end are supported. In the latter case, the error distribution and types are learned separately from the two mate reads, since it is a common experience that the two quality profiles can differ.

Whole genome mapping software packages such as BWA [START_REF] Li | Fast and accurate short read alignment with burrows-wheeler transform[END_REF] are designed to provide suboptimal alignments sacrifying accuracy for speed. For applications where the optimal alignment is required, specialized tools like the Genome Analysis Toolkit (GATK) [START_REF] Mckenna | The genome analysis toolkit: a mapreduce framework for analyzing next-generation dna sequencing data[END_REF] can be used to post-process .bam files and refine the reads alignment. Following a similar philosophy, our parameter learning module re-aligns reads onto the reference by means of the ultrafast bit-vector-based semi-global alignment algorithm described in [START_REF] Myers | A fast bit-vector algorithm for approximate string matching based on dynamic programming[END_REF] and implemented in [START_REF] Šošić | Edlib: a c/c++ library for fast, exact sequence alignment using edit distance[END_REF]. We call non-matching characters in this alignment as sequencing errors. At first sight, this could be considered as a rough simplification since insertions, deletions and mismatches could be due to genomic variations rather than to errors. In practice, we observed that the higher probability of errors (compared to that of variations) makes this approximation acceptable. However, in order to reduce the bias due to misinterpretation of variants, in the genotype generation module we implemented a function to return a (multiploid) reference with the most probable alleles. This reference can be used as input for the parameter learning step, in order to reduce the possible wrong contribution of variants. When a read is processed, the module aligns it with all the alleles of the reference, and then decides for the one that maximizes the alignment score.

Due to the technology-specific distribution and type of errors, the learning module collects separate statistics for each position of the alignment string, as well as for each type of error. This approach allows to obtain a more accurate and realistic simulation. For example, in short-reads sequencing by synthesis, the probability of erroneously calling T and C as G increases as the chemistry degrades with the succession of cycles. In particular, we keep and update a transition table that stores the conditional probability of encountering a certain alignment symbol (match, mismatch, insert, delete, end of sequence), given the preceding position symbol. Moreover, we maintain a histogram for each alignment symbol/position with the distribution of the per-base quality score.

Mismatches require to decide which character to use for replacement. As mentioned above, substitutions are not equiprobable, and depend on several technological factors like the sequencing machinery and chemistry. Although a random choice of the replacement character would still be appropriate for most applications, it might not produce realistic results deceiving software for systematic bias removal (see [START_REF] Laehnemann | Denoising dna deep sequencing datahigh-throughput sequencing errors and their correction[END_REF] for an in-depth discussion on these tools). In order to drive the correct replacement, we maintain a (positional) substitution table with the probability of each base to be replaced with another.

Unlike deletions that do not require any additional information, insertions need a mechanism to produce the sequence to embed in the read at the simulation stage. Keeping track of sequences involved in insertion events would not be enough for such a mechanism since they depend on several contextual variables including the neighborhood. In this case, however, almost all the applications just need that the embedded string causes the introduction of a gap in the alignment of the reads. This can easily be achieved ensuring to insert immediately after the insertion point a string that completely mismatches the reference.

Although the above statistics would be enough to generate technology-specific reads, they do not capture parameters that are dependent on the specific biases introduced with library preparation. The two most important of these parameters are: stutter noise (introduced with PCR amplification), and fragment size distribution for paired-end sequencing.

Stutter noise introduces/removes spurious units on microsatellites (aka short tandem repeats). The entity of this phenomenon is more dependent on the unit length than on its content; thus, in order to quantify it, for each possible unit length (microsatellites are defined in the range from 1bp to 6bp), we keep track of the discrepancy in the number of units between pure tandem repeats located in the reads and the corresponding sequence in the reference.

When dealing with paired-end sequencing, the mate reads are derived from the two endpoints of a (possibly longer) DNA fragment. According to the experimental setup, some part of the sequence might either be read twice or not be sequenced. For a more accurate mimic of the DNA fragmentation, we collect statistics on the distribution of the insert sizes, as well as we store the information about the orientation the paired reads.

C. Reads generation

The read generation module takes as input the alleles (produced by the genotype generator) and a profile of sequencer, and returns a set of reads either in fastq or bam format. Alleles derived from different organisms and different executions of the genotype generator can be mixed in order to simulate metagenomics.

Conceptually, the reads simulation mimics the process of library preparation and sequencing. Given a chromosome, the module performs four steps

1) It selects one allele at random, 2) It partitions it into fragments, 3) It introduces PCR artifacts, 4) It extracts from each fragment the corresponding read(s).

Per each chromosome, these steps are repeated until the desired user-defind average coverage depth is reached.

Unfortunately, in real sequencing experiments, allelic frequencies of variations are neither equiprobable nor constant. This is due to two factors: first, the amount of DNA extracted for each allele cannot be tightly controlled and, secondly, fragments size selection causes coverage fluctuations. We simulate the differences of DNA amount modeling the selection of the alleles as a Bernoulli process of probability 1/p where p is the number of alleles. Coverage fluctuations are then obtained simulating fragmentation.

Fragmentation depends on the type of sequencing. For single-end sequencing, only the size of the sequenced portion can be known, while it is impossible to determine the real length of segments, and therefore, in this case we do not actually perform the fragmentation, but we rather constrain fragments to be as long as the corresponding read. On the other hand, for paired-end sequencing, reads are extracted from the endpoints of a long fragment whose length can be measured after the alignment with the reference. Hence, we can learn the distribution of fragment lengths in the learning phase, and then perform the fragmentation accordingly.

Before producing the read(s), fragments are preprocessed to introduce a bias similar to that of stutter noise due to PCR amplification [START_REF] Aponte | Sequence-based analysis of stutter at str loci: Characterization and utility[END_REF]. This technique, used during library preparation, tends to add or delete copies of the repeated unit of short tandem repeats deceiving tandem repeat genotyping algorithms. Pure tandem repeats are pinpointed through a simple regular-expression-based procedure.

Reads generation has the twofold goal of extracting the sequence of the nucleotides of the read and injecting sequencing errors. As mentioned in section III-B, these errors are modeled as non-matching characters in the alignment of the read with the corresponding allele. We thus generate a plausible alignment to drive the extraction of correct nucleotides from the reference, and errors from the parameters learned in the previous phase. The alignment symbols are generated as the sequence of nodes belonging to a random walk on a multi-partite graph, where the probability of passing from a node of the k-th independent set to a node of the k + 1th independent set is extracted from the k-th transition table learned in the previous step. The process stops when an end of sequence symbol is generated. Figure 2 shows a small multi-partite graph where an example of random walk and the corresponding alignment are highlighted.

The process of read creation consists of looping over the symbols of an allignment, copying the sequence of nucleotides from the endpoint of a fragment, or inserting an error. There are three possible types of errors: deletions, mismatches, and insertions. Deletions are the easiest type of mutation to cope with, since they do not require any supplementary information: when a deletion symbol is found, the corresponding character in the reference is just skipped. Mismatches require to select a replacement character: this selection is In order to model specific phenomena (e.g. polyG tails), we maintain a different substitution table for each possible position in the read. Insertions are complicated by the difficulty of generating the "most appropriate" sequence to be inserted in the read. An approach similar to that we used for mismatches would not work properly, as the inserted sequence depends on too many parameters (i.e. the length of the insertion, its position on the read, the neighborhood, etc.), and keeping track of all of them would require a huge training set. We observed that a "good" insertion sequence should force an alignment software to open a gap. This can easily be achieved by inserting a random sequence provided this does not match any neighborhood of the insertion in the reference.

The read generation process is iterated over all the chromosomes until the global average coverage is reached. As it is the case for real sequencing, the random length of the fragments/reads causes fluctuations of the local coverage, and it leaves open the possibility for certain loci to be under/over represented.

IV. RESULTS

In this section, we describe our experiments aimed at showing that OmniSim is able to produce realistic reads for different sequencing technologies while being fast enough to be used in practice. We used a workstation equipped with an Intel(R) Core(TM) i5-3570 CPU @ 3.40GHz and 8Gb of RAM, running Linux Fedora 29 (Workstation Edition). In order to provide an accurate estimation of the overall running time, all the tests have been run using only a single thread.

Our tests consisted of generating two realistic alleles for a human sample, learning a model for an Illumina paired-end 100bp reads sequencing and a model for an Oxford Nanopore sequencing, and finally producing 100 million of synthetic reads.

We Downloaded the latest snapshot of dbSNP [START_REF] Sherry | dbsnp: the ncbi database of genetic variation[END_REF] (build 152) consisting in a 14Gb gzip compressed .vcf file annotating 690M variants (including alternative variants of the same locus) and used it to generate two alleles from the hg38.p12 reference (accession GCF 000001405.38).

As for the model we used:

• a uniform subsample of 164.09 million of paired-end reads of length 100bp (14Gb) from the Illumina platinum genomes [START_REF] Eberle | A reference data set of 5.4 million phased human variants validated by genetic inheritance from sequencing a three-generation 17-member pedigree[END_REF] of the CEPH1463 NA12877 individual (sample accession SAMEA1573614) • a uniform subsample of 7.95 million of single-end reads (17Gb) of variable length (up to 1kbp) from an Oxford Nanopore sequencing [START_REF] Jain | Nanopore sequencing and assembly of a human genome with ultra-long reads[END_REF] of the CEPH1463 NA12878/GM12878 individual. Both the NA12877 and NA12878 samples have been aligned to the hg38 (GCA 000001405.15) primary build with BWA-MEM.

TABLE I RUNNING TIME OF THE THREE MODULES

Module

Illumina pe Nanopore se Alleles generation 40m.50s Error model learning 1h:07m.05s 8m.31s Production of 100k reads 16m.16s 1h:25m.05s

In Table I we report the running time of the three modules depicted in Figure 1. Alleles generation is independent of the simulated sequencing technology, but may need to be repeated several times when simulating a population. For Illumina reads, this could cause the first step of the pipeline to take a non-negligible fraction of the overall running time. However, it worth noting that the generation of the alleles has only short-range dependencies and, thus, this module can be easily accelerated by running parallel instances on different chromosomes. The running time of the learning module does not depend on the size (in bp) of the sequences, but rather on their number. In fact, the two samples we used for learning have roughly the same size (although a very different number of sequences) and the NA12878 has longer reads but smaller running time. This result is not surprising, given the linear time cost of the bit-vector-based semi-global alignment algorithm described in [START_REF] Myers | A fast bit-vector algorithm for approximate string matching based on dynamic programming[END_REF]. When generating the Illumina profile, most of the time is spent by edlib [START_REF] Šošić | Edlib: a c/c++ library for fast, exact sequence alignment using edit distance[END_REF] for the initialization of the data structure, and this cost is not amortized because of the little length of the reads. Reads generation running time is deeply affected by the length of the reads. In fact, most of the time is spent in the generation of the alignment sequence shown in Figure 2. However, even in this case the module can easily be accelerated by a parallel computation of multiple alignments. Summing up, as shown in table I, the entire pipeline of OmniSim can be run in a reasonable time independently from the simulated sequencing technology.

We have investigated the composition of the generated reads and compared it with that of a real sequencing. To this end, we computed the positional distribution of nucleotides and reported the average and standard deviation. As shown in Figure 3, OmniSim has the same average distribution of the corresponding sequencing but with a slightly lower standard deviation (notice that the lower standard deviation is due to the smaller number of generated reads compared to the real ones). Considering the uneven distribution of nucleotides on the human reference and the biases of sequencing technologies that tend to under-represent CG-rich regions, this result suggests that OmniSim was able to learn these biases and generate reads whose coverage distribution reflects that of a real sequencing. Comparing Illumina sequencing and Oxford Nanopore, we observed that the natural higher variability of the latter is preserved by OmniSim. In fact, the standard deviation of the simulated Nanopore sequencing is higher than that of Illumina with the same proportions of the real ones.

Finally, Figure 4 reports the positional average value of the quality scores (and the first and third quartile on the error bar) for both the real and the simulated reads. As for Figure 3 (b), in order to enhance visualization, we grouped the measures in buckets of 10bps. Figure 4 (a) shows that OmniSim was able to learn the typical error distribution of Illumina reads: few initial bases with increasing quality, a long plateau of high quality, and a 20bp-long slope with decreasing quality values. Figure 4 (b), instead, shows a comparison of a real Oxford Nanopore sequencing and OmniSim. Interestingly, in this case our simulator was not only able to learn the distribution of the average quality values, but also to mimic the variability generating exactly the same standard deviation.

V. CONCLUSIONS

In-silico simulation of sequencing has emerged as a promising practical tool to produce genomic data provided with accurate genotype annotations, in order to be used as a validation benchmark. This data is thus essential to drive bioinformaticians during the development and testing of new analysis methods. In the daily practice, however, developers have to face the fact that simulators are usually the result of specific needs, and they are hence far from general. Moreover, it is not uncommon for older simulators not to be maintained or updated, with the practical result of producing reads that no longer reflect the current technology. Selecting the simulator that better fits the user needs, has become so complex that the authors of the survey in [START_REF] Escalona | A comparison of tools for the simulation of genomic next-generation sequencing data[END_REF] ended up providing a decision tree to guide developers.

In this paper we presented OmniSim: a new simulation tool that is specifically designed to be independent on the sequencing technology and type of genomic variant. OmniSim uses a dedicated module to learn the technological parameters of the sequencing machinery. This allows to keep our simulator up-to-date by simply re-running the learning module on the sequencing data of the latest technology. Genomic modifications are mapped into sequences of edit operations on a reference, enabling OmniSim to handle every possible type of structural variant. Experiments on two real datasets (one paired-end 100bp Illumina sequencing, and one Oxford Nanopore sequencing) have demonstrated that OmniSim is able to accurately simulate very different sequencing technologies producing reads with the same nucleotide composition and error profile of the real ones. OmniSim has not demanding hardware requirements and can run on a relatively small work-station accomplishing the entire pipeline (alleles generation, error profile learing, reads generation) in a few hours.

Fig. 1 .

 1 Fig. 1. Structure of the simulation pipeline: the three modules of our tool are in rectangular boxes

Fig. 2 .

 2 Fig.2. Example of multi-partite graph able to generate an alignment string up to 6 characters. M is for matches, ! for mismatches, I for insertions, D for deletions and $ for the end of sequence. In red the example alignment that generate the sequence MD!I.

Fig. 3 .

 3 Fig. 3. Composition of reads

Fig. 4 .

 4 Fig. 4. Comparison of the error profile of a real sequencing and a simulated one


Charlotte Saverna

Yves Aurégan

Vincent Pagneux

Non-Reciprocal Scattering in Shear Flow

published or not. The documents may come

I. INTRODUCTION

Adding flow to an acoustic system is one of the simplest ways to make it non-reciprocal. Recently, it has been used to obtain sound isolation [START_REF] Fleury | Sound isolation and giant linear nonreciprocity in a compact acoustic circulator[END_REF] or asymmetric propagation [START_REF] Naghdi | Acoustic nonreciprocity in Coriolis mean flow systems[END_REF] , in the broad context of research on unidirectional transmission devices for acoustic waves [START_REF] Fleury | Nonreciprocal acoustics[END_REF][START_REF] Maznev | Reciprocity in reflection and transmission: What is a 'phonon diode' ?[END_REF] . It appears that sound-flow interactions have been studied for a very long time [START_REF] Ostashev | Acoustics in Moving Inhomogeneous Media, Second Edition[END_REF] , early works being motivated by the study of the deflection of sound in the wind [START_REF] Delaroche | Sur l'influence que le vent exerce dans la propagation du son, sous le rapport de son intensité[END_REF] . For instance, in the Theory of Sound, Rayleigh [START_REF] Rayleigh | Theory of Sound[END_REF] uses geometric acoustics to show in which direction the sound is deflected according to the direction of the wind. Subsequently, varieties of studies were then led on the effect of shear flows on acoustic propagation, from a simple flow velocity discontinuity [START_REF] Miles | On the reflection of sound at an interface of relative motion[END_REF][START_REF] Miles | On the disturbed motion of a plane vortex sheet[END_REF][START_REF] Ribner | Reflection, transmission, and amplification of sound by a moving medium[END_REF] , to piece-wise constant profiles [START_REF] Amiet | Refraction of sound by a shear layer[END_REF][START_REF] Yeh | A Further Note on the Reflection and Transmission of Sound Waves by a Moving Fluid Layer[END_REF] or continuous parallel shear flows [START_REF] Graham | Effect of a sheared layer on plane waves in a fluid[END_REF][START_REF] Jones | The Scattering of Sound by a Simple Shear Layer[END_REF][START_REF] Koutsoyannis | Acoustic resonances and sound scattering by a shear layer[END_REF] . Many other papers deal with the propagation of sound in parallel shear flow in the presence of walls, either for ducts [START_REF] Pridmore-Brown | Sound propagation in a fluid flowing through an attenuating duct[END_REF] or for reflection problems [START_REF] Brand | Reflection of sound by a boundary layer[END_REF][START_REF] Cantrell | Interaction between sound and flow in acoustic cavities : mass, momentum and energy considerations[END_REF][START_REF] Goldstein | Effect of shear on duct wall impedance[END_REF]25,[START_REF] Mariano | Effect of wall shear layers on the sound attenuation in acoustically lined rectangular ducts[END_REF][START_REF] Mungur | Acoustic wave propagation in a sheared fluid contained in a duct[END_REF][START_REF] Nagel | Boundary layer effects on sound in a circular duct[END_REF][START_REF] Nayfeh | Effect of the acoustic boundary layer on the wave propagation in ducts[END_REF] . For impedance walls, it can be noted that analytic expressions exist for very small boundary layers [START_REF] Brambley | A well-posed boundary condition for acoustic liners in straight ducts with flow[END_REF][START_REF] Gabard | A comparison of impedance boundary conditions for flow acoustics[END_REF] , which tend to the one given by Ingard 25 and Myers [START_REF] Myers | On the acoustic boundary condition in the presence of flow[END_REF] for vanishing boundary layers.

Though rarely mentioned as such, non-reciprocity is obviously an existing feature in most of the previously cited works. In this paper, we will inspect the effect of basic parallel shear flows on the non-reciprocity properties. In the model we use, both the mean Mach number and the height of the shear layer can be varied independently. Though not so realistic, this will allow us to highlight some key aspects of the influence of shear. We consider two specific scattering problems: i) reflection and transmission, when the fluid is moving at two different velocities above and under a layer of shear flow with a linear profile, ii) reflection, when the fluid is flowing above a lined wall with a shear layer. The first problem we consider is similar to the ones studied in jet-like configurations by Amiet 1 and Jones [START_REF] Jones | The Scattering of Sound by a Simple Shear Layer[END_REF] , except we consider it as a scattering problem. This enables us to get a unia) charlotte.saverna@univ-lemans.fr

M 2 M(y) O x y 𝜃 1 M 1 A 1 B 1 =R 1 A 1 A 2 B 2 =T 1 A 1 𝜹 𝜃 2 =-𝜃 1 B 2 =R 2 A 2 B 1 =T 2 A 2 a) A 2 =0 b) A 1 =0
FIG. 1. Scheme of the problem of scattering by a layer of shear flow a) in the case when the incident wave comes from below the shear layer, b) when it comes from above the shear layer.

tary relation on the full scattering matrix. In the second problem, we highlight other non-reciprocal effects such as unilateral perfect absorption or asymmetric lateral displacement of Gaussian beams. The latter phenomenon is known as Goos-Hänchen shift [START_REF] Artmann | Berechnung der Seitenversetzung des totalreflektierten Strahles (Calculation of the lateral displacement of the totally reflected beam)[END_REF][START_REF] Goos | Ein neuer und fundamentaler Versuch zur Totalreflexion (a new and fundamental attempt at total reflection)[END_REF] in optics [START_REF] Breazeale | Interaction of ultrasonic waves incident at the Rayleigh angle onto a liquidsolid interface[END_REF][START_REF] Declercq | Theory of the backward beam displacement on periodically corrugated surfaces and its relation to leaky Scholte-Stoneley waves[END_REF] or Schoch effect in ultrasounds [START_REF] Emile | Measurement of the nonlinear Goos-Hänchen effect for Gaussian optical beams[END_REF][START_REF] Herbison | Ultrasonic diffraction effects on periodic surfaces[END_REF][START_REF] Tamir | Lateral displacement of optical beams at multilayered and periodic structures[END_REF] . In the presence of shear flow, we show it is possible to reach large lateral displacement due to non-reciprocity.

II. GENERAL SCATTERING PROPERTIES OF A SHEAR FLOW

We consider the problem of the scattering of a plane wave by a layer of shear flow between two fluids moving at constant Mach number M 1 and M 2 , as described in Fig. 1. We study this problem using a scattering matrix S linking the incident amplitudes A 1,2 to the scattered amplitudes B 1,2 .

A. Energy conservation

The propagation equation for a wave in a parallel flow is given by 44 :

D t [(∂ 2 x + ∂ 2 y )p - 1 c 2 0 D 2 t p] -2∂ y U (y)∂ 2 x,y p = 0, (1)
where D t is a convective derivative term, D t = ∂ t + U ∂ x , and U (y) is the mean velocity of the flow, such that M (y) = U (y)/c 0 , with c 0 the constant sound speed. To study the energy conservation, non-trivial in shear flow, we consider the pressure p to be a plane harmonic wave, such that p(x, y, t) = [P (y)e -i(ωt-αx) ], where ω = 2πf is the circular frequency, α = K sin(θ i )/(1 + M i sin(θ i )) is the horizontal component of the wave wave number, K = ω/c 0 is the natural wave number, M i = M 1,2 the mean Mach in the medium of incidence and θ i = θ 1,2 , the angle of incidence. Then, we obtain the Pridmore-Brown 44 equation for P :

P + 2αM (y) K -M (y)α P + [(K -M (y)α) 2 -α 2 ]P = 0, ( 2
)
where the prime notation stands for a derivation with respect to y. The time-averaged energy equation, as given by Myers 38 is:

∇• < I >= -< D >, (3)
with I the energy flux given by I = (ρ 0 v + ρv 0 )(c 2 0 ρ/ρ 0 + v 0 .v), and D the dissipation, given by D

= ρ 0 v 0 • ((∇ × v) × v) -ρv • ((∇ × v 0 ) × v 0 )
. In these expressions, ρ 0 is the mean density, ρ is the density variation which in adiabatic conditions is equal to p/c 2 0 , v 0 is the mean velocity field and v is the acoustic particle velocity which we can evaluate using linearized Euler equations with flow as:

v = V x (y)e -i(ωt-αx) , V y (y)e -i(ωt-αx) , (4)
with:

ρ 0 c 0 V x (y) = [ αP K-M α -M P (K-M α) 2 ]. ρ 0 c 0 V y (y) = -iP K-M α . (5)
Substituting these values in the definition of D, we find:

< D >= M M 2(K -M α) 3 (P + 2αM K -M α P + [(K -M α) 2 -α 2 ]P )P * . (6)
Using Eq. 2, < D >= 0 and Eq. 3 can now be evaluated to:

(P P * ) (K -M (y)α) 2 y2 y1 = 0, (7)
meaning that the energy flux of plane waves is conserved across the shear layer.

B. Scattering matrix

In the domains of constant velocity, the pressure is given by: P (y) = A 1 e iβ1y + B 1 e -iβ1y for y < 0, P (y) = A 2 e -iβ2y + B 2 e iβ2y for y > δ,

where

β 2 j = (K -M j α) 2 -α 2 (j = 1, 2
). Using Eq. 7 and the formalism shown in Fig. 1 we obtain:

β 1 (K -M 1 α) 2 |A 1 | 2 + β 2 (K -M 2 α) 2 |A 2 | 2 = β 1 (K -M 1 α) 2 |B 1 | 2 + β 2 (K -M 2 α) 2 |B 2 | 2 . ( 9
)
We define a scattering matrix S by:

B 1 B 2 = S A 1 A 2 = R 1 T 2 T 1 R 2 A 1 A 2 (10)
where R 1 (respectively R 2 ) is the reflection coefficient for a wave coming from medium 1 (resp. 2) ; T 1 (respectively T 2 ) is the transmission coefficient for a wave coming from medium 1 (resp. 2). With the correct normalization on the amplitudes, one can re-write equation 9 as S †S = I, where the † exponent stands for Hermitian transposition, S is the normalized scattering matrix and I is the 4-by-4 identity matrix. This last relation is exactly the same relation as in the no-flow case. The scattering coefficients then satisfy the four relations:

         |R 1 | 2 + β2 β1 ( K-M1α K-M2α ) 2 |T 1 | 2 = 1, |R 2 | 2 + β1 β2 ( K-M2α K-M1α ) 2 |T 2 | 2 = 1, |R 1 | = |R 2 |, |T 1 | = β1 β2 ( K-M2α K-M1α ) 2 |T 2 |. (11)
Equations 11 show that the relation between the scattering coefficients is independent of the shear layer thickness. However, their values do depend on it (see appendix A). In the case when on of the incident waves has zeroamplitude, Eq. 11 verifies the relation found by Jones [START_REF] Jones | The Scattering of Sound by a Simple Shear Layer[END_REF] .

C. Dispersion relation and anomalous refraction properties

A useful representation of scattering between two media of different characteristics is to represent dispersion relations in both media in the (a, b) plane (a and b being respectively the horizontal and vertical component of the wave number, normalized by K). In our case, we have the following dispersion relation:

a 2 j + b 2 j = (1 -M j a j ) 2 , (12)
with j the index of the considered medium. In this formulation, we recognize a conic equation, with an eccentricity directly given by |M j |. Therefore, we will have a circle of radius 1 in the (a, b) plane if

M j = 0 ; an ellipse if |M j | < 1 ; a parabola if |M j | = 1 ; an hyperbola if |M j | > 1.
This property is illustrated in Fig. 2. Plotting such a relation gives access to both the group and the phase velocity directions for the refracted wave (see Fig. 3). As we study Gaussian beams, the direction of the group velocity, given by the direction of the normal to the ellipse at the point (a(θ), b(θ)) will have a particular interest as it will give the apparent direction of the beam. The beams are built by summing weighed plane waves as follows:

p(x, y) = +∞ -∞ e -(α-α0) 2 /w 2 p(α, x, y)dα, (13)
with α 0 the horizontal wave number in the (imposed) main direction of the beam, and w the waist of the beam. In Fig. 3, these beams are represented for the propagation from a medium 1 at rest to a medium 2 where M 2 = 0.4. By looking at the ellipses in the (a, b) plane, and since α is constant between the two media, we can see that for some incidence, θ 1 and θ 2 can be of opposite sign. This result can be interpreted as an equivalent negative refractive index, analog to what is observed for electromagnetic waves in left handed materials 49 and in hyperbolic meta-materials [START_REF] Christensen | Anisotropic metamaterials for full control of Acoustic Waves[END_REF][START_REF] Garcìa-Chocano | Negative refraction and energy funneling by hyperbolic materials: An experimental demonstration in acoustics[END_REF][START_REF] Ju | Asymmetric acoustic transmission with a lossy gradient-index metasurface[END_REF] . The shear layer thickness will influence the value of the moduli of R and T but not the direction of the transmitted beam which only depends on the Mach numbers.

III. REFLECTION ON AN IMPEDANCE WALL

We now consider an acoustically treated wall with uniform admittance Y 0 placed in y = 0 with a shear flow (see Fig. 4). Y 0 is such that its imaginary part stands for the reactance of the wall, and its real part stands for the resistance of the wall ( (Y 0 ) < 0 means that there are losses at the wall, (Y 0 ) > 0 means that there is gain at the wall).

A. Energy conservation

We know that Eq. 7 is still valid, and that the pressure field for y > δ is such that:

P (y) = Ae -iβy + Be +iβy , ( 14
)
where

β 2 = (K -M 0 α) 2 -α 2
. Now, due to the impedance wall, there is a mixed boundary condition involving the admittance Y 0 at y = 0:

P (0) = iKY 0 P (0). ( 15
)
The scattering problem that we consider is only with reflection and the energy flux conservation 7 applied between y 1 = 0 and y 2 = δ leads to:

|R| 2 = 1 + (α 2 + β 2 ) Kβ |P (0)| 2 (Y ). ( 16
)
As long as the flow is subsonic, one can check that β is positive. Therefore, the sign of |R| 2 -1 will only depend on the sign of (Y ), whatever the mean flow velocity or the shear layer thickness. In particular, as noted by Campos 7 , if the wall is purely reactive ( (Y ) = 0) then the reflection will always be perfect i.e. |R| = 1.

B. Non-reciprocal perfect absorption

The computation of R will be achieved numerically using the method detailed in appendix B. The addition of flow leads to non-reciprocity testified in the value of the absorption coefficient, defined as A = 1 -|R| 2 . Indeed, an incoming plane wave will not be absorbed to the same extend at the wall depending on whether it comes in the direction of the flow or against the flow. This nonreciprocal property is highlighted in Fig. 5, for a liner of admittance Y = 0.1i -0.1. This figure shows clearly that varying Kδ allows to reach high absorption for one incidence. This is a typical behaviour as long as the shear layer thickness is large enough. One can also remark on this figure that the higher the Mach number and the height of the shear layer, the more non-reciprocal the system becomes. For a clearer view of this asymmetry, we plot in Fig. 6 the value of the absorption coefficient as a function of both θ and Kδ and θ and M 0 . A clear difference can be noted whether incident waves come from the left (in the direction of the flow) or from the right (against the flow). In particular, what is striking on Fig. 6 is that, above a certain thickness of shear layer, waves coming against the flow with a wide angle are fully reflected. This could be explained by a strong refraction of the rays, such that the incident wave does not see the effects of the liner before being pushed back towards the "sky". This corresponds to the "turning points", evoqued by Rienstra [START_REF] Rienstra | Solutions and Properties of the Pridmore-Brown Equation[END_REF] . On the other side, waves coming with the flow (θ > 0) with an almost grazing incidence can either be perfectly reflected or perfectly absorbed depending on the shear layer thickness. The quasi-periodic phenomenon which appears is related to resonance phenomena inside the shear layer. For the previous Mach number, for small incidence angles, the absorption does not vary a lot. Therefore, we look at the variation of A with respect to both θ and M 0 to see if the peak absorption can be displaced towards normal incidence for higher Mach numbers (see Fig. 6). The position of the near-toperfect absorption for negative θ then depends linearly on M 0 . Using numerical simulations, we want to design a system which would enhance the absorption of a given liner by adding a specific flow above it, so that it reaches perfect absorption, or zero reflection. The admittance of the liner is fixed at Y 0 = -0.1 + 0.1i, and the incident angle such that |θ| = π 6 . Without flow, in such a situation A ≈ 0.36 (see Fig. 8). In order to get a physically consistent simulation, it is not the direction of the phase velocity which is held constant, but the direction of the group velocity θ g . A relation between θ and θ g can be found using the dispersion relation (as shown in Fig. 2):

θ g = arctan[tan(θ) + M 0 cos(θ) ], (17)
or θ = arcsin[sin(θ g ) cos(θ g )

1 cos 2 (θ) -M 2 0 - M 0 cos 2 (θ g )].
We now look at the variations of |R| as both M 0 and Kδ vary for θ g = π/6. In Fig. 7, an optimum is found, leading to no reflection when θ g = -π/6 (counter-flow propagation). We precisely read it as being linked to the values M 0 = 0.355 and Kδ = 69.68. Non-reciprocity can also be pointed out here, as there is no way to reach perfect absorption if the opposite incidence is considered. We can check that theses values are also associated with an optimal absorption in terms of incidence by plotting A as a function of θ without flow or with a flow corresponding to the optimum previously identified. Figure 8 gives us this confirmation: the absorption (the absorption coefficient is 0.36 without flow) can be enhanced all the way to perfect absorption with the right flow parameters.

Let us now illustrate the non-reciprocal perfect absorption with Gaussian beams following either upstream or downstream propagation. The results are plotted in Fig. 9. The reflection becomes highly non-reciprocal when flow is added: perfect absorption here is only achieved for waves coming counter-flow. Waves traveling in the direction of the flow actually encounter more reflection than they did in a medium at rest on the same liner (as pointed in Fig. 8). This behaviour is in agreement with the refraction of acoustic rays in shear flow: rays are pushed towards the regions of low local wave speed, consequently they are trapped near the wall when propagating upstream and pushed towards the sky when propagating downstream.

C. Non-reciprocal Goos-Hänchen shift

When the reflection implies an incidence-depending phase shift, Gaussian beams are shifted along the wall. This phenomenon is called Goos-Hänchen shift [START_REF] Goos | Ein neuer und fundamentaler Versuch zur Totalreflexion (a new and fundamental attempt at total reflection)[END_REF] , and can be quantified using Artmann's formula 2 :

d = - ∂φ R ∂α = - (1 + M sin(θ)) 2 K cos(θ) ∂φ R ∂θ , (18)
where φ R is the angle of the reflection coefficient, and d is the lateral shift, as displayed in Fig. 10. We also define the displacement of the beam perpendicularly to its direction of propagation: d = cos(θ)d.

Without flow

We know that, without flow, the reflection coefficient is given by:

R = cos(θ) + Y 0 cos(θ) -Y 0 . (19)
In the following, Y 0 = Y 1 + iY 2 . Then its angle φ R is given by:

φ R = arctan( Y 2 cos(θ) + Y 1 ) + arctan( Y 2 cos(θ) -Y 1
). (20) In the case of a lossless wall, Y 1 = 0, Eq. 20 reduces to

φ R = 2 arctan Y 2 cos(θ) . ( 21
)
We can then compute the value of d , and we get:

Kd = - 2Y 2 sin(θ) cos 2 (θ) + Y 2 2 . ( 22
)
By differentiating Eq. 22 with respect to Y 2 , we find a maximum displacement of the beam Kd max = ± tan(θ) for Y 2 = ± cos(θ). Thus, it appears that, for lossless walls, the value of Kd max is small for small incident angles. Nevertheless, it is possible to make this limitation disappear by adding losses to the liner. Indeed, in the case when Y 1 = 0, we find that Kd can take much larger values. With losses (Y 1 = 0), Eq. 22 takes the form:

Kd = -Y 2 sin(θ) cos 2 (θ) + |Y | 2 [(cos(θ) + Y 1 ) 2 + Y 2 2 ][(cos(θ) -Y 1 ) 2 + Y 2 2 ] (23)
where, for Y 1 = ± cos(θ) and Y 2 is close to zero, the denominator can go to zero leading to large Kd . This result is also the one leading to perfect absorption. Therefore, large lateral shift is associated with large absorption coefficient. This phenomenon is easily understandable if one is to look at the phase portrait of the reflection coefficient when it varies with θ (see Fig. 11): phase variation of R is enhanced near zero reflection.

With Flow

When flow is added, the Goos-Hänchen displacement has to be computed numerically. It is then interesting to study the impact of both the flow and the wall admittance on the beam displacement. Results shown in Fig. 12 illustrate that the effect of the flow is predominant. Indeed, the two maximal values of displacement visible without flow for θ ≈ ±π/6 are still visible when the flow is added but are negligible compared to the effect of the flow with a much larger displacement at θ ≈ -π/4. Remark that this displacement is comparable to the large displacement of beams observed in the case of reflection on multi-layered materials in optics [START_REF] Declercq | The Schoch effect to distinguish between different liquids in closed containers[END_REF][START_REF] Farmani | Graphene Plasmonic : Switching Applications[END_REF][START_REF] Tamir | Lateral displacement of optical beams at multilayered and periodic structures[END_REF] .

As in the previous section, where we illustrated nonreciprocal perfect absorption, Gaussian beams can nicely illustrate the Goos-Hänchen effect (Fig. 13). By choosing the angle for which the flow has the most visible effect (θ = -0.89), we see a large lateral shift. Besides, a strong non-reciprocal effect is demonstrated by taking the opposite angle.

IV. CONCLUSION

In the global context of acoustics with flow we have shown that the presence of a shear layer leads to nonreciprocal effects whose origin can be found in the known effects due to convection. We have illustrated some of these phenomena by considering diffraction by a shear layer that can lead to negative refraction, as well as absorption by an acoustic material and lateral displacement of a beam along this material which both depend on the direction of propagation relatively to the flow. In this paper, we have only considered homogeneous problems along the direction of flow. Interesting new effects appear as soon as inhomogeneities are considered along the flow direction.

FIG. 2 .

 2 FIG. 2. (color online) Dispersion relation curve in the (a; b) plane (a and b being the dimensionless horizontal and vertical wave numbers), when M0 = 0 (blue circle) ; M0 = 0.5 (red ellipse) and M0 = 1.2 (yellow hyperbola).

FIG. 3

 3 FIG. 3. (color online) Beams coming from medium 1 with an incident angle a) θ = -0.2 or b) θ = 0.2. c) represents the dispersion relations in the plane (ai, bi), with Mach numbers M1 = 0 and M2 = 0.4.Conservation of a at the interfaces allows us to find the direction of the phase velocity by finding the point of the ellipse with abscissa a. The normal to the ellipse in this point then gives the direction of the group velocity (i.e. the visible direction of the beam). Green items are linked to medium 1 and red items to medium 2. Solid arrows stand for group velocity direction and dashed arrows for phase velocity. In medium 1, only one arrow is represented as phase and group velocity are co-linear in a medium at rest.

FIG. 4 .

 4 FIG. 4. Geometry of the problem for the reflection of a plane on an admittance wall.

FIG. 5 .

 5 FIG. 5. (color online) Value of the absorption coefficient A as a function of θ for several flows, described by their mean Mach number M0 and the height of their shear layer δ for a liner of admittance Y = 0.1i -0.1.

FIG. 7 .

 7 FIG. 7. (color online) Value of A as a function of M0 and Kδ for waves coming either in the direction of the flow (case a) or counter-flow (case b). Black lines show the values of M0 and Kδ for which (R) = 0 (dashed lines) or (R) = 0 (solid lines). It appears that we can find one crossing point, corresponding to |R| = 0 (i.e. A = 1) for waves coming counter-flow only, for Kδ = 69.68 and M0 = 0.355.

FIG. 8 .

 8 FIG. 8. Value of A as a function of θg with and without flow (M0 = 0.355), for Y = -0.1 + 0.1i and Kδ = 69.68.

FIG. 10 .

 10 FIG.10. Scheme of the lateral shift phenomenon for an incident beam. d is the shift along the wall, d is the shift perpendicularly to the direction of the beam.

FIG. 11 .

 11 FIG. 11. Phase portrait of R as θ varies from -π/2 to +π/2, with M0 = 0.23, Kδ = 15.04 and a) Y0 = 0.1i, b) Y0 = 0.1i -0.2. When losses are added on the wall and in a situation where perfect absorption can be achieved, the reflection coefficient can cross the zero point, and therefore experience a sudden phase jump of amplitude π, resulting in a high phase gradient and therefore a large lateral displacement of the beam.

J. Acoust. Soc. Am. / 29 October 2019

J. Acoust. Soc. Am. / 29 October 2019

APPENDIX A: NUMERICAL DETERMINATION OF THE SCATTERING COEFFICIENTS IN A TRANSMISSION PROBLEM

To find the scattering coefficients, we need to solve Eq. 2. For this, we define Q:

Then Eq. 2 can be re-written as a first-order differential equation for Q, and becomes a Ricatti equation. Moreover, to make its integration easier, a dimensionless variable s is introduced, such that s = y δ . The integration will then be on [0, 1]. Finally,the equation to solve is: • if the incident wave comes from medium 1, A 2 is null, and Eq. 8 gives us the following boundary condition in s = 1:

where b 2 = β 2 /K. Equation A2 could have be solved with a common numerical method (such as a Runge-Kutta algorithm), but large heights of shear layer resulted in divergences in the result. Thus, we implemented a Magnus-Moebius scheme (as described in [START_REF] Pagneux | Multimodal admittance method in waveguides and singularity behavior at high frequencies[END_REF] ) for s going from 1 to 0. In order to do so, it is necessary to go back to Eq. 2 which, when written on s, becomes:

When put under a matrix form, Eq. A4 becomes:

The matrix involved in Eq. A5 will later be denoted M(s) Discretizing the domain in N +1 points (intervals of size h = 1/N ) allows us to apply the following scheme :

where

Note that the exponential in Eq. A6 is a matrix exponential. At each iteration, it is possible to compute Q(nh) from Q((n + 1)h):

where the E i,j are the coefficients of e Hn . After N iterations, we can access R 1 via Eq. 8, and we find:

where b 1 = β1 K .

• if the incident wave comes from medium 2, A 1 is null and we can now write a boundary condition in s = 0:

Once again, Eq. A5 is solved iteratively for s going from 0 to 1. This time, we have:

Thus:

and we can access R 2 using:

The determination of transmission coefficients T 1 and T 2 , now that Q is known in [0, 1], requires to integrate Eq. A4 one more time, but in the opposite direction. Here also, we will need to study the two cases of upward and downward traveling incident waves.

• if the incident wave comes from medium 1, we know from Eq. 8 that the boundary condition in s = 0 is:

) where P is the y-dependant component of the pressure field normalized by the incident amplitude A 1 . Then, by solving for P on [0, 1] using the scheme in Eq. A10, we can access T 1 using:

(A14)

• if the incident wave comes from medium 2, then this time the boundary condition is:

where P (s) is normalized by A 2 . Then, by solving for P on [1, 0] using the scheme in Eq. A6, we can access T 2 using:

The whole matrix can then be reconstituted. It turns out from Eq. A2 that three parameters will influence the values of the scattering coefficients: the dimensionless width of the shear layer Kδ, the relative Mach number M 2 -M 1 and the direction of the incident wave given by the horizontal wave number a.

APPENDIX B: NUMERICAL DETERMINATION OF R IN A REFLECTION PROBLEM

In order to determine R, we apply the same method as for the transmission problem : solving Eq. A2 for s in [0; 1], but this time the boundary condition will be given by the wall admittance:

Here also, the resolution is led to its end using a Magnus-Moebius scheme 30 of the second order as described in [START_REF] Pagneux | Multimodal admittance method in waveguides and singularity behavior at high frequencies[END_REF] .

In order to compute the solution, we go back to Eq. A5, and after writing P and P in discrete forms, we get: