Elisabeth Chevreau
email: elisabeth.chevreau@inra.fr

Nicolas Dousset

Clément Joffrion

Andréa Richer

Aurélie Charrier

Emilie Vergne

Agroinfiltration is a key factor to improve the efficiency of apple and pear transformation

Keywords: Apple, pear, transformation, agroinfiltration, surfactant

Genetic transformation of apple and pear is mainly performed via Agrobacterium tumefaciens. The average efficiency of transformation remains low in most laboratories and the success of the transformation experiments is very variable. Improving the efficiency and the reproducibility of apple and pear transformation is thus highly desirable. Adventitious regeneration ability of the explants is not a limiting factor for apple or pear transformation. The present study focuses on improving the frequency of stably transformed cells in the explants following A. tumefaciens inoculation. We report here the results of 36 independent transformation experiments on 'Gala' apple, performed with 10 different binary plasmids, comparing three methods of A. tumefaciens inoculation. Agroinfiltration (Agrobacterium-mediated vacuum infiltration) of the explants in a

Introduction

Apple (Malus domestica Borkh) and European pear (Pyrus communis L.) are two closely related species belonging to the subtribe Malinae [START_REF] Potter | Phylogeny and classification of Rosaceae[END_REF]. Recent analysis of the European pear genome [START_REF] Chagné | The draft genome sequence of European pear (Pyrus communis L. 'Bartlett')[END_REF] confirmed its similarity to apple. Because of its high commercial importance, apple is one of the first woody plant species which was successfully transformed via Agrobacterium in 1989 [START_REF] James | Genetic transformation of apple (Malus pumila Mill) using a disarmed Ti-binary vector[END_REF]. Despite this early success and the high number of apple genotypes amenable to genetic transformation, the rate of transformation of apple has remained low. Most published data indicate a rate of transformation below 2% on a per explant basis and transformation rates above 10% are exceptional (Malnoy et al. 2008a). Furthermore, the repeatability of apple transformation experiments is low and it is still difficult for most laboratories to produce large number of transgenic lines in a short period of time. With the recent production of a high-quality de novo assembly of the apple genome [START_REF] Daccord | High-quality de novo assembly of the apple genome and methylome dynamics of early fruit development[END_REF], analysis of gene function has become a major research topic in apple. Thus, improving the efficiency of apple transformation rates is a key issue. The first successful gene transfer in European pear was in 1996 [START_REF] Mourgues | Efficient Agrobacterium-mediated transformation and recovery of transgenic plants from pear (Pyrus communis L.)[END_REF]). The transformation rates of pear are generally low (1 -4%) with a few exceptions reaching 40-50% (Malnoy et al 2008b). Currently, the number of clonal genotypes amenable to genetic transformation is much lower in Pyrus (< 20) than in Malus (> 50). However this field of research is very active in pear and the need for high transformation rates is important.

Most of the biotechnology protocols that have been first developed on apple were later adapted to pear with success. For both species, the efficiency of Agrobacterium-mediated transformation is strongly linked to the regeneration ability, which is highly genotype dependent. In apple as well as in pear, adventitious regeneration from young leaf explants is very efficient and several genotypes reach rates of regeneration of 100% with a high number of buds per regenerating leaves in absence of transformation [START_REF] Yepes | Factors that affect leaf regeneration efficiency in apple, and effect of antibiotics on morphogenesis[END_REF]Aldwinckle 1994, Chevreau et al. 1989). The second key parameter for the recovery of transgenic plants is the frequency of stably transformed cells in the explant following Agrobacterium tumefaciens inoculation. Despite a high natural susceptibility of apple and pear to wild A. tumefaciens strains in orchards, stable T-DNA insertion following A. tumefaciens inoculation is a rate-limiting factor for the production of apple and pear transgenic plants. The first inoculation method used for apple transformation was immersion of prewounded leaves in the inoculum for a few minutes (Malnoy et al 2008a). Crushing the leaf blade with non-traumatic forceps prior to immersion into the inoculum increased the efficiency of apple transformation [START_REF] Norelli | Leaf wounding increases efficiency of Agrobacteriummediated transformation of apple[END_REF], probably by increasing the number of infection sites for Agrobacterium in the leaves. For pear transformation, the main inoculation procedure is leaf wounding with a scalpel dipped into the inoculum (Malnoy et al. 2008b).

Agrobacterium-mediated vacuum infiltration, also called agroinfiltration, is a standard method for transient expression of foreign genes in plant tissues. It has been used as a simple and rapid method for assaying gene function in a variety of plant species such as tobacco, grapevine, potato, rose, rice and poplar [START_REF] Palanichelvam | Agroinfiltration of Cauliflower Mosaic Virus gene VI elicits hypersensitive response in Nicotiana species[END_REF][START_REF] Santos-Rosa | Development of a transient expression system in grapevine via agro-infiltration[END_REF][START_REF] Bashkar | Agrobacterium-mediated transient gene expression and silencing: a rapid tool for functional gene assay in potato[END_REF][START_REF] Yasmin | Transient gene expression in rose petals via Agrobacterium infiltration[END_REF][START_REF] Andrieu | Anin planta, Agrobacterium-mediated transient gene expression method for inducing gene silencing in rice (Oryza sativa L.) leaves[END_REF], Takata and Erikson 2012). In this case, non-integrated copies of T-DNA present in the nucleus of plant cells permit to reach a level of expression 1000 fold higher than in stable transformants (Jansen and Gardner 1989). By applying a selection pressure after the period of co-culture, agroinfiltration has also be used for the production of stable transgenic lines of recalcitrant plants species such as cotton, grapevine or common bean (Haq 2004[START_REF] Zottini | Agroinfiltration of grapevine leaves for fast transient assays of gene expression and for longterm production of stable transformed cells[END_REF][START_REF] Mukeshimana | Factors influencing regeneration and Agrobacterium tumefaciens-mediated transformation of common bean (Phaseolus vulgaris L.)[END_REF]). In the present report, we took advantage of the capacity of agroinfiltration to enhance gene delivery to increase both the efficiency and the regularity of transformation of apple and pear and we propose an optimized protocol based on agroinfiltration of leaf explants.

Material and Methods

Biological material

The experiments were performed on two genotypes: the apple 'Gala' and the pear 'Conference'. Both genotypes have a very high adventitious regeneration potential. In vitro proliferating shoot cultures of the apple 'Gala' were micropropagated on Murashige and Skoog (MS) (1962) medium supplemented with 2.22 µM 6-benzyladenine (BA) and 0.5µM 3indolebutyric acid (IBA). Cultures of the pear 'Conference' were micropropagated as described by Leblay et al. (1991) on a derivative of Lepoivre's medium supplemented with 2.22 µM 6-BA and 0.5 µM IBA. All cultures were grown at 22-24° C with a 16:8 h light:dark photoperiod (cool white fluorescent tubes, 40-60 µmol m -2 s -1 ) and transferred to fresh medium every 4 weeks. The youngest leaves of 4-week-old in vitro shoots were used for A. tumefaciens inoculation. The bacterial strain used in all transformation experiments was the A. tumefaciens EHA105 containing the ternary plasmid pBBR1MCS-5 with a constitutive virG gene (van der [START_REF] Van Der Fits | The ternary transformation system: constitutive virG on a compatible plasmid dramatically increases Agrobacterium-mediated plant transformation[END_REF]. The binary plasmids used in the different transformation experiments are described in Supplementary Table 1.

Transformation methods

For apple transformation, the bacteria were re-suspended in liquid apple regeneration medium consisting of MS medium containing 22.7 µM thidiazuron (TDZ), 2.15 µM IBA and 100 µM acetosyringone. Two methods of inoculation were tested: "Crushing" consisted in crushing the leaf blade with non-traumatic forceps prior to immersion into the inoculum at 10 8 bacteria/ml (Fig. 1b); "Agroinfiltration" consisted in vacuum-infiltrating the leaves in an inoculum at 10 8 bacteria/ml containing or not a wetting agent (Silwet L-77® (Lehle Seeds, USA) at 0.002 % v/v) under -0.9 bar during one minute. After inoculation, the leaves were wounded transversely with a scalpel and plated adaxial side down on co-cultivation medium for two days in the dark at 22-24 °C. The cocultivation medium was apple regeneration medium containing 100 µM acetosyringone.

For pear transformation, the bacteria were re-suspended in liquid pear micropropagation medium (described above) with 100 µM acetosyringone. Two methods of inoculation were tested: "Scalpel" consisted in wounding with a scalpel dipped into an inoculum at 10 8 bacteria/ml ; "Agroinfiltration" consisted in vacuum-infiltrating the leaves in an inoculum at 10 8 bacteria/ml without wetting agent. After inoculation, the leaves were wounded transversely with a scalpel in case of "Agroinfiltration" and plated abaxial side down on co-cultivation medium for two days in the dark at 22-24 °C. The co-cultivation medium was pear regeneration medium [START_REF] Mourgues | Efficient Agrobacterium-mediated transformation and recovery of transgenic plants from pear (Pyrus communis L.)[END_REF]) containing 9 µMl TDZ, 2.68 µM naphthalene acetic acid (NAA) and 100 µM acetosyringone.

For apple and pear, at the end of the co-culture, the leaves were plated (adaxial side down for apple and abaxial side down for pear) on their respective regeneration medium containing 300 mg/l cefotaxime, 150 mg/l timentin and 100 mg/l kanamycin. The explants were kept in the dark and transferred to fresh medium every month for six months. When the inoculated A. tumefaciens strain contained a GUS marker gene, histochemical GUS test was performed one month after inoculation.

Samples were incubated overnight at 37°C, in a solution containing 3 mM 5-bromo-4-chloro-3indolyl-b-D-glucuronic acid, 4 mM potassium ferricyanide, 0.05 mM potassium ferrocyanide, 10 mM EDTA and phosphate buffer (0.02 M, pH 7.2). Clearing was achieved using 70% (v/v) ethanol.

Appearance of adventitious buds was monitored for a period of six months. All regenerated buds were micropropagated on the same medium as their mother plants, with the addition of 300 mg/l cefotaxime, 150 mg/l timentin and 100 mg/l kanamycin. Because all the regenerated buds appeared on the wounds, it was possible to isolate very clearly independent regeneration events occurring on separate wounds. If multiple buds arose from the same wound, only one shoot was cloned and counted as a transformation event, to avoid overestimation of the transformation rate. Presence of transgenes and absence of contaminating agrobacteria were monitored by PCR. After PCR analysis of the lines which survived on the kanamycin containing medium, the final transformation rate was estimated as the number of PCR positive lines per 100 leaves.

DNA extraction and PCR

Presence of transgenes and absence of contaminating agrobacteria were monitored by PCR.

Genomic DNA of apple and pear leaves was extracted as described in [START_REF] Fulton | Microprep protocol for extraction of DNA from tomato and other herbaceous plants[END_REF] [START_REF] Karimi | GATEWAY vectors for Agrobacterium-mediated plant transformation[END_REF]) and (iii) elongation factor 1α (EF1α) coding gene as a marker of plant DNA suitability for PCR (CTCTTGGTGTCAGGCAAATG and TCAAGGTTGGTGGACCTCTC on AJ223969). Amplifications were performed using GoTaq®

Flexi DNA Polymerase (Promega, Madison, WI, USA) according to the manufacturer's recommendations. The PCR reaction conditions were identical for the three genes except the hybridization step which was at 55°C and not 58°C for A. tumefaciens detection primers: 95°C for 5 min, followed by 35 cycles at 95°C for 30 s, 58°C for 45 s, 72°C for 1 min, with a final extension at 72°C for 5 min. The PCR products were separated on a 1.5 % agarose gel. A typical result of this detection method applied during transformation experiment 249 (Supplementary Table 2) is given in Supplementary Fig. 1.

Statistical analysis

Statistical analyses were performed with XLSTAT by using the nonparametric Wilcoxon test (p < 0.05) for apple data and a Khi 2 independence test (α = 0.05) for pear data.

Results

We report here the synthesis of 36 independent transformation experiments on 'Gala' apple, performed with 10 different binary plasmids on a total of 12,971 explants. Results of each inoculation method were recorded on 12 independent experiments performed with 4 to 5 different binary plasmids. The rate of failure (experiments producing none transgenic lines) was 50% for "Crushing", 33.3% for "Agroinfiltration without Silwet" and 16.7% for "Agroinfiltration with Silwet". PCR analyses indicated that 90% of the kanamycin resistant lines contained a complete T-DNA and no Agrobacterium (Supplementary Table 2 and Fig. 1). They were counted as transgenic lines. Despite the great variability of transformation rates during these experiments, box-plots in Fig. 1a indicate a slight increase of transformation rate (non-significant using Wilcoxon test) with "Agroinfiltration without Silwet" compared to "Crushing". "Agroinfiltration with Silwet" was significantly more efficient than "Crushing" (Wilcoxon test, p = 0.049), with a mean rate of 5.8 % and a maximum rate of 30.2%. No correlation between the rate of transformation and the plasmid used was observed. The success of this method can be linked to a very high level of GUS expression one month after inoculation (Fig. 1 d ande) compared to the "Crushing" method (Fig. 1c).

The protocol described here for 'Gala' apple was successfully tested in one experiment for the transformation of other apple genotypes: 'Ariane', 'Golden Delicious', 'Greensleeves' and 'M26' with transformation rates between 1 and 9% (Table 1). In our hands, 'Conference' pear presents a very high ability for adventitious bud regeneration (Fig. 2a) and susceptibility to kanamycin selection (Fig. 2b). This genotype is amenable to transformation with rates of success between 2 and 10% (Fig. 2c andd). Here we tested the "Agroinfiltration" method without Silwet compared to the standard "Scalpel" method, in three independent experiments performed with the same binary plasmid (392p9N35s-Mr5FB1orf). PCR analyses indicated that 93% of the kanamycin resistant lines analyzed contained a complete T-DNA and no medium increases the transient transformation efficiency in rose [START_REF] Yasmin | Transient gene expression in rose petals via Agrobacterium infiltration[END_REF] or poplar (Takata and Erikson 2012). But in some cases, the addition of a surfactant at high concentrations can also be toxic to plant tissues [START_REF] Wu | Factors influencing successful Agrobacteriummediated genetic transformation of wheat[END_REF]). In our study, higher concentrations of Silwet (0.005% and 0.01%) caused severe necrosis of apple leaves in a preliminary regeneration experiment without agroinfiltration. Finally, we also showed that agroinfiltration of pear leaves led to very high rates of transformation, even in the absence of surfactant. Positive results were obtained by [START_REF] Habashi | Comparing vacuum Agroinoculation and common Agroinoculation in two pear (Pyrus communis L.) cultivars "Bartlett" and "Harrow Delight[END_REF] using agroinfiltration on the pear cultivars 'Bartlett' and 'Harrow Delight' with much lower transformation efficiencies (7 to 12%).In conclusion, the protocols of transformation described here are efficient for 'Gala' apple and 'Conference' pear. They permit to reduce the frequency of failed experiments and to increase the maximum rate of transformation of these genotypes. They should be efficient for other genotypes after minor optimization of some parameters such as hormonal concentrations.

Fig. 1

 1 Fig. 1 Main transformation results on apple. a Transformation rate (%) of 'Gala' apple according to the method of A. tumefaciens inoculation: "Crushing", "Agroinfiltration without Silwet" and "Agroinfiltration with Silwet". Each box-plot corresponds to 12 independent experiments performed on 150 to 850 explants. +: mean, : maximum, ×: outlier. b Numerous wounds visible on 'Gala' apple leaf just after crushing with non-traumatic forceps. c Low efficiency of GUS histochemical test performed on 'Gala' apple leaf one month after "Crushing" inoculation. d, eHigh efficiency of GUS histochemical test performed on 'Gala' apple leaf one month after "Agroinfiltration with Silwet" inoculation: adaxial side (d) and abaxial side (e).

Table 1Transformation rates of several apple varieties after agroinfiltration with Silwet

 1Transformation

	Genotype	# leaves inoculated	# regenerated buds	rate of regeneration	# transgenic lines	rate of transformation
	M26	140	23	16.4 %	12	8.6 %
	Greensleeves	141	8	5.7 %	3	2.1 %
	Golden Delicious	245	18	7.3 %	4	1.6 %
	Ariane	329	5	1.5 %	3	0.9 %

Acknowledgments

This project has received funding from the French Government managed by the Research National Agency (ANR) under the Investment for the Future program (GENIUS project ANR11-BTBR-0001).

Agrobacterium (Supplementary Table 3). They were counted as transgenic lines. Results in Fig. 2d indicate that "Agroinfiltration" was significantly more efficient than "Scalpel" in the three experiments (Khi 2 test, p < 0.0001) with rates of transformation between 50 and 80%.

Discussion

Agroinfiltration is frequently used as an efficient method for transient transformation assays because the high density of agrobacteria forced into the intercellular spaces increases the efficiency of T-DNA delivery. Here we demonstrated that this method also improved the stable integration of T-DNA into apple and pear cells. Penetration of a bacterial solution by a simple syringe infiltration is not efficient in apple or pear leaves. It is necessary to use a low-pressure vacuum to efficiently force the bacterial suspension into the apple or pear leaves. In addition we demonstrated that, for apple, the addition of a surfactant at a low concentration (0.002%) boosts the transformation

Disclosure of potential conflicts of interest

The authors declare that they have no conflict of interest.


Chang Yang

Fabrice Deluzet

Jacek Narski

On the numerical resolution of anisotropic equations with high order differential operators arising in plasma physics

Keywords: Plasma physic, Anisotropic equation, Asymptotic-Preserving schemes

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

This paper is devoted to the study of Asymptotic-Preserving numerical methods for the resolution of equations incorporating high order differential operators. This class of numerical methods have been introduced for the efficient resolution of singular perturbation problems, initially for diffusive limit of kinetic equations [START_REF] Klar | An Asymptotic-Induced Scheme for Nonstationary Transport Equations in the Diffusive Limit[END_REF][START_REF] Klar | An Asymptotic Preserving Numerical Scheme for Kinetic Equations in the Low Mach Number Limit[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF] and then extended to other frameworks, including plasma physics (see for instance [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review. Lecture Notes for Summer School on "Methods and Models of Kinetic Theory[END_REF][START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF][START_REF] Degond | Asymptotic-Preserving methods and multiscale models for plasma physics[END_REF] for reviews), the context addressed by the present work. The specific applications targeted here are indeed related to simulation of tokamak plasmas by means of fluid models [START_REF] Minjeaud | High order approximation of a tokamak edge plasma transport minimal model with Bohm boundary conditions[END_REF] and more specifically by reduced models [START_REF] Ricci | Simulation of plasma turbulence in scrape-off layer conditions: the GBS code, simulation results and code validation[END_REF][START_REF] Halpern | The GBS code for tokamak scrape-off layer simulations[END_REF][START_REF] Nespoli | Non-linear simulations of the TCV Scrape-Off Layer[END_REF][START_REF] Isoardi | Penalization modeling of a limiter in the Tokamak edge plasma[END_REF][START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code SolEdge2d-Eirene[END_REF][START_REF] Galassi | Flux expansion effect on turbulent transport in 3d global simulations[END_REF] derived for the investigation of the plasma edge dynamics in the so-called scrape-of-layer [START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF]. These equations are characterized by large anisotropies due to the intense magnetization of the plasma. The particles are indeed bounded to the magnetic field lines in the directions perpendicular to the magnetic field while their motion is almost free along this direction. On the macroscopic level, this is translated by operators with an anisotropy all the more severe than the magnetic field is intense. This feature originates difficulties in the derivation of efficient numerical methods. This is particularly evident when the boundary conditions imposed at each magnetic field line end translate either the periodicity of the torus or the Bohm criterion [START_REF] Rogers | Low-Frequency Turbulence in a Linear Magnetized Plasma[END_REF][START_REF] Tamain | TOKAM-3d: A 3d fluid code for transport and turbulence in the edge plasma of Tokamaks[END_REF][START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF]. In these specific frameworks, the operator associated with the parallel dynamics does not define, by itself, a well posed problem. Since, direct discretizations of these problems are only able to capture the parallel dynamics for large anisotropies, they give rise to a system matrix with a large condition number, hardly invertible (see for instance [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF][START_REF] Ratnani | Anisotropic Diffusion in Toroidal geometries[END_REF]).

Similar problems have been investigated in precedent realizations (see [START_REF] Degond | An asymptotic preserving scheme for strongly anisotropic elliptic problems[END_REF][START_REF] Besse | Efficient numerical methods for strongly anisotropic elliptic equations[END_REF][START_REF] Yang | Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths[END_REF] and [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF]). These works are devoted to the derivation of Asymptotic-Preserving methods for anisotropic elliptic or diffusion equations. The purpose is to derive system matrices with a condition number uniformly bounded with respect to the ansitropy strength. This is achieved thanks to the introduction of an auxiliary variable. By this means the stiffness of the problem is cancelled and the degeneracy of the equation in the limit of infinite anisotropy prevented. One way to implement these ideas is referred to as the Micro Macro decomposition introduced in [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF]. The purpose here is therefore to develop the concept of AP methods for the numerical resolution of the vorticity equation proposed in [START_REF] Tamain | TOKAM-3d: A 3d fluid code for transport and turbulence in the edge plasma of Tokamaks[END_REF]. In this model, the electric potential is solution to an equation with a fourth order differential operator along the transverse directions. The transposition of the Micro Macro decomposition into this framework was first proposed in [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF] and then resumed in [START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF]. However, this approach revealed to be inoperative, leading to a matrix condition number scaling badly with respect to the mesh refinement, precisely 1/h 6 if h is the typical mesh size. This trend is to be compared to the conditioning of standard discretizations proportional to both 1/h 4 and ε -1 /h 2 , ε -1 being the anisotropy strength. The issue raised here is therefore twofold. On the one hand, standard discretizations suffer from a bad condition number for large anisotropies. On the other hand, the proposed AP methods are only effective for meshes with a coarse resolution.

The aim of the present work is therefore to introduce numerical methods that offer the advantage of a condition number bounded irrespective of the anisotropy strength together with a dependency to the mesh resolution comparable to standard methods. These methods share the Asymptotic-Preserving property but they are derived with a different approach compared to [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF]. The first one consists in elaborating a Duality Based decomposition [START_REF] Degond | An asymptotic preserving scheme for strongly anisotropic elliptic problems[END_REF][START_REF] Besse | Efficient numerical methods for strongly anisotropic elliptic equations[END_REF] for this problem. This approach is specific to coordinates aligned with the magnetic field as proposed in some works [START_REF] Dudson | BOUT++: A framework for parallel plasma fluid simulations[END_REF]. It could also be transposed to coordinates aligned with the main (poloidal) component of the magnetic field which are routinely used (see [START_REF] Tamain | The TOKAM3x code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF]). The advantage of this approach is its high efficiency as outlined in [START_REF] Yang | Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths[END_REF]. The resolution of the anisotropic problems may therefore be much more computationally efficient than the resolution of the corresponding isotropic problem.

The second method proposed herein is a two field iterated method in the spirit of [START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF]. This method can address arbitrary field geometries with non adapted coordinates and meshes. These two methods are compared against the Micro Macro decomposition proposed in [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF] and a standard discretization of the problem. Beyond the computational efficiency which is assessed herin, in a two dimensional framework, we investigate the condition number of these different methods. This is indeed a crucial property to envision three dimensional computations. As emphasized in [START_REF] Yang | Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths[END_REF], the computational efficiency of sparse direct solvers (used for two dimensional computations either in [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF] or in the present work), is poor for three dimensional computations. The scaling of both the memory requirements and computational costs increasing extremely non linearly with the number of unknowns. It is therefore essential to develop numerical methods with controlled condition numbers. The simulation of plasma turbulence also requires to refine the mesh to capture the scale at which the turbulence develops. The scaling of the numerical method with respect to the mesh refinement is therefore an important property to investigate. This defines the main purpose of the present work. The last purpose of this paper is to provide comprehensive arguments to highlight the weaknesses of the methods proposed so far for this class of problems and how they are circumvented thanks to the ones introduced in the present paper.

The organization of this paper is the following. The definition of the problem at hand is introduced in Sec 2, together with a simplified model problem containing the difficulty that the numerical methods need to tackle. The degeneracy of the systems in the asymptotic of infinite anisotropy is highlighted. The Asymptotic-Preserving formulations are stated in Sec. 3. An emphasis is made on how the degeneracy of the system is prevented for each formulation. In particular, the Two Field Iterated method is introduced. It advances the two component in a fixed-point sequence for which a proof of convergence is provided. The material related to the discretization is specified in Sec 4 with analytic estimates of the condition number outlining their dependence to both the anisotropy strength and the mesh size. The properties of these methods are numerically investigated in Sec 5.

From the vorticity equation to a simplified model problem

The vorticity equation

The starting point is the vorticity equation derived in [START_REF] Tamain | TOKAM-3d: A 3d fluid code for transport and turbulence in the edge plasma of Tokamaks[END_REF]. The computational domain is two dimensional and square Ω = Ω x × Ω y = [0, 1] × [0, 1]. In a simplified magnetic field geometry, assuming a field straight and aligned with the y axis, the electric potential φ is the solution to the problem (MP)

                                     - ∂ ∂t ∂ 2 φ ε ∂x 2 + ∂ 4 φ ε ∂x 4 - 1 ε ∂ 2 φ ε ∂y 2 = f ε , (x, y) ∈ Ω , ∂φ ε ∂x (0, y) = 0 , ∂ 3 φ ε ∂x 3 (0, y) = 0 , y ∈ Ω y , ∂φ ε ∂x (1, y) = 0 , ∂ 3 φ ε ∂x 3 (1, y) = 0 , y ∈ Ω y , ∂φ ε ∂y (x, 0) = ε 1 -exp Λ -φ ε (x, 0) , x ∈ Ω x , ∂φ ε ∂y (x, 1) = -ε 1 -exp Λ -φ ε (x, 1) , x ∈ Ω x . (1a) (1b) (1c) (1d) (1e)
where, ε is an asymptotic parameter defining the anisotropy strength of the problem, Λ is given data to the problem, representing a reference potential. The purpose here is to extract a minimal problem containing the difficulty to be solved for an efficient resolution of this problem in the limit of vanishing ε. The targeted value for the applications are defined by the range ε ∈ [10 -8 , 10 -6 ]. The purpose of this paper is therefore to derive efficient numerical methods for the resolution of this problem for these severe anisotropies. We refer to [START_REF] Tamain | TOKAM-3d: A 3d fluid code for transport and turbulence in the edge plasma of Tokamaks[END_REF] for the physical background and to [START_REF] Negulescu | Existence and uniqueness of the electric potential profile in the edge of tokamak plasmas when constrained by the plasma-wall boundary physics[END_REF] for the analysis of the well posedness of this problem. The difficulty in handling this problem from a numerical point of view can easily be pointed out thanks to the analysis of the dominant operator in the limit ε → 0. Setting formally ε = 0 in (1) yields

       ∂ 2 φ 0 ∂y 2 = 0 (x, y) ∈ Ω x × Ω y , ∂φ 0 ∂y (x, 0) = 0 , ∂φ 0 ∂y (x, 1) = 0 , x ∈ Ω x .
(

) 2
where φ 0 is the limit of φ ε when ε → 0. This system admits an infinite number of solutions since all the functions χ with no aligned gradient, i.e. ∂χ/∂y = 0, satisfy this system. For large anisotropies (ε ≪ 1) the system matrix issued from the discretization of the system (1) is consistent with the degenerate problem (2). This explains a deterioration of the condition number of this matrix with vanishing ε.

The limit problem is defined as the well posed system providing the limit of the solution when ε → 0. First, from Eq. ( 2) we conclude that φ 0 does not depend on the aligned coordinate y. Then integrating Eq. (1a) with respect to y and owing to Eqs. (1d) and (1a), and assuming that φ ε (x, 1) and φ ε (x, 0) do not depend on time, the limit problem can be stated as

                 ∂ 4 φ 0 ∂x 4 (x) = 1 0 f 0 (x, y)dy + G 1 (x) -G 0 (x) , x ∈ Ω x , ∂φ 0 ∂x (0) = 0 , ∂ 3 φ 0 ∂x 3 (0) = 0 , ∂φ 0 ∂x (1) = 0 , ∂ 3 φ0 ∂x 3 (1) = 0 , (3a) (3b) (3c)
with

G 1 = -1 + exp Λ -φ 0 (x) , G 0 = 1 -exp Λ -φ 0 (x) . (3d)

A simplified model problem

The aim here is to simplify the problem (1) to isolate the difficulty raised by the degeneracy of the system for vanishing ε. First, the stationary problem will be considered. The reason for this choice is the following. The discretization of the time derivative in Eq. (1a) provides a contribution on the diagonal of the system matrix that may offset the anisotropy of the problem. It is therefore frequent, to adjust the value of the time step to preserve a tractable matrix conditioning rather than setting its value according to the physics of interest. Tackling directly the stationary system is the best way to develop numerical methods for which the time step can be chosen with respect to the dynamic of interest and arbitrarily with respect to the anisotropy strength.

Second, the boundary conditions Eqs. (1d) and (1e) are non linear in the unknown. A linearization is routinely implemented for the numerical resolution of the problem, with

φ(x, 0) ε - 1 ε ∂φ ε ∂y (x, 0) = g 0 (x) , φ(x, 1) ε + 1 ε ∂φ ε ∂y (x, 1) = g 1 (x) ,
where the following notations g 0 (x) = Λ(x, 0) and g 1 (x) = Λ(x, 1) are used. The simplified toy problem can thus be stated as

(SP)                                      ∂ 4 φ ε ∂x 4 - 1 ε ∂ 2 φ ε ∂y 2 = f ε , (x, y) ∈ Ω x × Ω y , ∂φ ε ∂x (0, y) = 0 , ∂ 3 φ ε ∂x 3 (0, y) = 0 , x ∈ Ω x , ∂φ ε ∂x (1, y) = 0 , ∂ 3 φ ε ∂x 3 (1, y) = 0 , x ∈ Ω x , φ ε (x, 0) - 1 ε ∂φ ε ∂y (x, 0) = g 0 (x) , y ∈ Ω y , φ ε (x, 1) + 1 ε ∂φ ε ∂y (x, 1) = g 1 (1) , y ∈ Ω y . (4a) (4b) (4c) (4d)
(4e) Note that the system (4) degenerates into a problem similar to (2), as the original problem (1) one does.

The limit problem related to (4) is defined as

(L)                  2φ 0 (x) + ∂ 4 φ 0 ∂x 4 (x) = 1 0 f 0 (x, y)dy + g 1 (x) + g 0 (x) , x ∈ Ω x , ∂φ 0 ∂x (0) = 0 , ∂ 3 φ 0 ∂x 3 (0) = 0 , ∂φ 0 ∂x (1) = 0 , ∂ 3 φ0 ∂x 3 (1) = 0 , (5a) (5b) (5c)
3. Asymptotic-Preserving formulations

Introduction

This section is devoted to the introduction of Asymptotic-Preserving methods and to a comparison with the Micro Macro method used in [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF]. The anisotropic problem defined by Eqs. ( 4) is a singular perturbation problem: setting ε = 0 in this system yields the ill posed problem [START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]. The aim of Asymptotic-Preserving methods is to introduce reformulated sets of equations equivalent to (4) for ε > 0. However, setting formally ε = 0 in these reformulated systems, the limit problem (5) is recovered. Therefore, the limit ε → 0 is regular in the reformulated systems.

The methodology operated in this work consists in cancelling the stiffness of the problem by introducing an auxiliary variable. This auxiliary variable does not always carry any physical meaning and can therefore be chosen arbitrarily. These choices define different numerical methods with properties making them more or less attractive regarding their computational efficiency.

In the sequel the dependence of the solution and right hand side with respect to ε will be drop in the notations.

The Duality Based decomposition

The Duality Based formulation consists in decomposing the solution of the problem into its mean and fluctuating parts φ(x, y) = φ(x) + φ ′ (x, y) ,

where φ = Π(φ) :=

1 0 φ(x, y)dy , φ ′ (x, y) := φ(x, y) -φ(x) . (6b)
The component φ is the mean of the solution along the magnetic field lines, while φ ′ represents the fluctuations of the solution around its mean value. Note that the fluctuation verifies

Π(φ ′ ) = 0 . ( 6c
)
The following duality property holds true for any function ψ

1 0 Π(ψ)ψ ′ dy = 0 . (7)
The fluctuating part ψ ′ and the mean part Π(ψ) of any function ψ are orthogonal for the scalar of L 2 (Ω y ). This property is used (more specifically in problems with heterogeneous diffusion coefficients) to derive the equation providing φ. It is derived by applying the mean operator Π onto Eq. (4a). This yields

∂ 4 φ ∂x 4 (x) = Π(f )(x) + g 0 (x) + g 1 (x) -φ ′ (x, 0) -φ ′ (x, 1) -2 φ(x) .
Substituting (6a) into Eq.(4a) provides the one giving φ ′ :

∂ 4 φ ′ ∂x 4 (x, y) - 1 ε ∂ 2 φ ′ ∂y 2 (x, y) = f (x, y) - ∂ 4 φ ∂x 4 (x, y) .
This finally provides the Duality Based reformulated problem consisting of a first system verified by the fluctuation (DB')

                                         ∂ 4 φ ′ ∂x 4 (x, y) - 1 ε ∂ 2 φ ′ ∂y 2 (x, y) = f (x, y) - ∂ 4 ∂x 4 φ(x, y) , (x, y) ∈ Ω x × Ω y , ∂φ ′ ∂x (0, y) = 0 , ∂ 3 φ ′ ∂x 3 (0, y) = 0 , y ∈ Ω y , ∂φ ′ ∂x (1, y) = 0 , ∂ 3 φ ′ ∂x 3 (1, y) = 0 , y ∈ Ω y , φ ′ (x, 0) - 1 ε ∂φ ′ ∂y (x, 0) = g 0 (x) -φ(x) , x ∈ Ω x , φ ′ (x, 1) + 1 ε ∂φ ′ ∂y (x, 1) = g 1 (x) -φ(x) , x ∈ Ω x , Π(φ ′ )(x) = 0 , x ∈ Ω x . (8a) (8b) (8c) (8d) (8e) (8f)
This set of equations is coupled to the system providing the mean part of the solution (DB)

                     2 φ(x) + ∂ 4 φ ∂x 4 (x) = Π(f )(x) + g 0 (x) + g 1 (x) , -φ ′ (x, 0) + φ ′ (x, 1) , x ∈ Ω x , ∂ φ ∂x (0) = 0 , ∂ 3 φ ∂x 3 (0) = 0 , ∂ φ ∂x (1) = 0 , ∂ 3 φ ∂x 3 (1) = 0 . (9a) (9b) (9c)
These two systems are coupled, on the one hand by the definition of the right hand side of Eqs. (8a), (8d) and (8e) and, on the other hand, by that of Eq. (9a). Note that the system providing φ does not incorporate stiff terms any more. It defines a well posed problem for the computation of the mean part provided that φ ′ is known. The system verified by the fluctuation is very similar to the original problem. However, setting ε = 0 into (8) yields

                                 ∂ 2 φ ′ ∂y 2 (x, y) = 0 , (x, y) ∈ Ω x × Ω y , ∂φ ′ ∂x (0, y) = 0 , ∂ 3 φ ′ ∂x 3 (0, y) = 0 , y ∈ Ω y , ∂φ ′ ∂x (1, y) = 0 , ∂ 3 φ ′ ∂x 3 (1, y) = 0 , y ∈ Ω y , ∂φ ′ ∂y (x, 0) = 0 , ∂φ ′ ∂y (x, 1) = 0 , x ∈ Ω x , Π(φ ′ )(x) = 0 , x ∈ Ω x . (10a) (10b) (10c) (10d) (10e)
which admits the null function as unique solution thanks to the zero mean constraints (10e) verified by the fluctuation. This proves that the system (8-9) defines a well posed set of equations for the computation of ( φ, φ ′ ) for ε ≥ 0. Note also that the limit problem is recovered from the "Duality Based" reformulated system when ε = 0. Indeed, from the problem [START_REF] Degond | Asymptotic-Preserving methods and multiscale models for plasma physics[END_REF] we conclude that the fluctuation vanishes with ε.

Inserting φ ′ = 0 into (9) yields the limit problem (5).

The Two Field Iterated formulation

The Two Field Iterated reformulation of the problem (4) is constructed by means of a different decomposition of the solution. Indeed, the zero mean value verified by the fluctuation is difficult to discretize when the coordinate system is unrelated to the magnetic field. To circumvent this difficulty, the solution is decomposed into φ(x, y) = p(x) + εq(x, y) .

In this decomposition, the solution components, namely p and q, are no more orthogonal in the sense of Eq. [START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code SolEdge2d-Eirene[END_REF]. The auxiliary variable q is introduced to rescale the parallel derivatives, with ∂q ∂y = 1 ε ∂φ ∂y cancelling thus the stiffness from the equation and yielding

∂ 4 φ ∂x 4 - ∂ 2 q ∂y 2 = f , (12a)
∂ 2 φ ∂y 2 = ε ∂ 2 q ∂y 2 . ( 12b
)
Note that Eq. (12b) enforces zero aligned derivatives for φ in the limit ε → 0. The next step consists in working these two equations to derive a well posed problem for both φ and q. To this end, each component of this two field system is advanced alternatively, until convergence, in a fixed-point iteration process defining (φ (k) , q (k) ) k>0 the sequence of iterates. The first system is intended to compute φ (k+1) from q (k) . It is derived from [START_REF] Degond | An asymptotic preserving scheme for strongly anisotropic elliptic problems[END_REF], multiplying Eq. (12a) by ε 0 > 0 and subtracting Eq. (12b) to obtain (TFI,φ)

                                       ε 0 ∂ 4 φ (k+1) ∂x 4 - ∂ 2 φ (k+1) ∂y 2 = ε 0 f -(ε -ε 0 ) ∂ 2 q (k) ∂y 2 , in Ω ∂φ (k+1) ∂x (0, y) = 0 , ∂ 3 φ (k+1) ∂x 3 (0, y) = 0 , ∂φ (k+1) ∂x (1, y) = 0 , ∂ 3 φ (k+1) ∂x 3 (1, y) = 0 , ∂φ (k+1) ∂y (x, 0) = ε 0 φ (k+1) (x, 0) -g 0 + (ε -ε 0 ) ∂q (k) ∂y (x, 0) , ∂φ (k+1) ∂y (x, 1) = ε 0 -φ (k+1) (x, 1) + g 1 + (ε -ε 0 ) ∂q (k) ∂y (x, 1) , (13a) (13b) (13c) (13d) (13e)
Note that ε 0 is a non vanishing parameter intended to verify ε ≪ ε 0 < 1. Therefore, this set of equations defines a well posed problem for φ (k+1) . The equation for q consists of Eq. (12b) in which ∂ 2 y φ is substituted by ε (f -∂ 4

x φ) upgraded with ∂ 4

x q (k+1) on the left hand side as well as its counterpart ∂ 4

x q (k) on the right hand side, to obtain

(TFI,q)                                        ε 0 ∂ 4 q (k+1) ∂x 4 - ∂ 2 q (k+1) ∂y 2 = f - ∂ 4 φ (k+1) ∂x 4 + ε 0 ∂ 4 q (k) ∂x 4 , in Ω ∂q (k+1) ∂x (0, y) = 0 , ∂ 3 q (k+1) ∂x 3 (0, y) = 0 , ∂q (k+1) ∂x (1, y) = 0 , ∂ 3 q (k+1) ∂x 3 (1, y) = 0 , ∂q (k+1) ∂y (x, 0) = φ (k+1) (x, 0) -g 0 , ∂q (k+1) ∂y (x, 1) = -φ (k+1) (x, 1) + g 1 . (14a) (14b) (14c) (14d) (14e)
Remark 3.1 (Uniqueness of q in the Two Field Iterated formulation). In the system (14), q (k+1) is determined up to a constant. Since only the derivatives of the auxiliary variable are used to compute (φ (k) ) k>0 , we can choose arbitrarily this constant. The value of q in one point of the domain (x 0 , y 0 ) may be set to 0 to restore the well posedness of the problem.

The asymptotic preserving property of the Two Field Iterated system is easily verified. First, we note that if the sequence (φ (k) , q (k) ) k>0 converges to (φ, q) then Eq. ( 12b) is recovered from Eq. (14a) and Eq. (12b) from both Eqs. (13a) and (14a). Now setting ε = 0, owing to Eq. (12b), the limit of the solution does not depend on y. Finally integrating Eq. (13a) with respect to y the limit problem ( 5) is recovered.

The convergence of the fixed point iterations can be proved regardless of the initial guess (arbitrarily) chosen to start the sequence and the anisotropy strength. Theorem 3.1. Let V be the functional space defined as

V := {v ∈ H 1 (Ω) | ∂ 2 x v ∈ L 2 (Ω)}.
For any (φ (0) , q (0) ) ∈ V × V, the sequence (φ (k) , q (k) ) k>0 defined by the iterative method [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF] converges to a fixed point (φ, q). The component φ of the stationary point solves uniquely the initial singular problem (4) for ε > 0 and the limit problem (5) when ε = 0.

The proof of this theorem is deferred to Appendix A.

A comparison with the Micro Macro formulation

The Micro Macro reformulation of the problem (4) proposed in [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF] is based on the decomposition (11) but with a different choice for q. It may be understood as a generalization of the Duality Based formulation in which

εq(x, y) = φ ′ (x, y) -φ ′ (x, 0) . ( 15
)
By this means, q verifies a homogeneous boundary condition, the so-called "inflow" condition [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF], rather than the zero mean condition which is more intricate to discretize with coordinates unrelated to the anisotropy direction. Note that q is re-scaled by a factor ε, hence the name "Micro Macro" decomposition. The sets of equations is inherited from ( 12) with

(MM)                                                                      ∂ 4 φ ∂x 4 (x, y) - ∂ 2 q ∂y 2 (x, y) = f (x, y) , (x, y) ∈ Ω x × Ω y , - ∂ 2 φ ∂y 2 (x, y) = -ε ∂ 2 q ∂y 2 , (x, y) ∈ Ω x × Ω y , ∂φ ∂x (0, y) = 0 , ∂ 3 φ ∂x 3 (0, y) = 0 , y ∈ Ω y , ∂φ ∂x (1, y) = 0 , ∂ 3 φ ∂x 3 (1, y) = 0 , y ∈ Ω y , ∂φ ∂y (x, 1) = ε -φ(x, 1) + g 1 (x) , x ∈ Ω x , ∂φ ∂y (x, 0) = ε φ(x, 0) -g 0 (x) , x ∈ Ω x , ∂q ∂y (x, 0) = φ(x, 0) -g 0 (x) , x ∈ Ω x , ∂q ∂y (x, 1) = -φ(x, 1) + g 1 (x) , x ∈ Ω x , q(x, 0) = 0 , x ∈ Ω x . (16a) (16b) (16c) (16d) (16e)
(16f) (16g) (16h) (16i)
Remark 3.2 (Uniqueness of q in the Micro Macro formulation). The auxiliary variable q is involved in the set of Eqs. (16a-16h) through its aligned derivatives. Considering this subset of equations, q is determined up to a function of the x-coordinate. The inflow condition (16i) prescribed at one end of each magnetic field line restores the unicity of the problem providing q and the well posedness of the Micro Macro system.

As stated before, the Micro Macro method may be regarded as a generalization of the Duality Based concepts to coordinate systems non adapted to the anisotropy direction. Note that this approach does not require a discretization of the limit problem (or a problem equivalent to it, as the system providing the mean part of the solution in the Duality Based approach). In the system [START_REF] Dimarco | Asymptotic preserving Implicit-Explicit Runge-Kutta methods for non linear kinetic equations[END_REF], the problem providing φ and q are strongly coupled. The Two Field Iterated formulation allows to decouple the computation of each field. They are advanced independently, thanks to mildly anisotropic problems. The conditioning of the linear system issued from the discretization of these systems is parametrized by ε ≪ ε 0 < 1 and is therefore bounded with respect to the anisotropy strength. These two systems are classical elliptic problems for which very efficient solvers (e.g. Krylov subspace, multi-grid methods) can be used for their efficient resolution. To date only direct sparse solver have been successful for the resolution of the linear system issued from the Micro Macro formulation. The inflow condition mandatory for the unicity of q is imposed on one end of each magnetic field lines. This prevents it from addressing closed magnetic field lines (that are not connected to the domain boundaries). In this respect, the Two Field Iterated method is more versatile, no restrictions on the magnetic field topology being required.

A finite difference space discretization

The Cartesian mesh is defined by the nodes at position (x i , y j ) with I and J the sets of indices for interior nodes and Ī and J the sets for all nodes, as defined by

I = {3, . . . , N x -2} , Ī = {1, . . . , N x } , J = {2, . . . , N y -1} , J = {1, . . . , N y } , (17)
the position of the nodes being defined, thanks to ∆x = L x /(N x -4), ∆y = L y /(N y -4), as

x i = (i -5 /2)∆x , y j = (j -5 /2)∆y , ∀(i, j) ∈ Ī × J , ,
Finally, we introduce the bijection I: (i, j) ∈ Ī × J → I(i, j) ∈ {1, . . . , N x • N y } that defines the ordering of the unknowns. With these notations, the following discrete differential operators can be introduced. For any vector v h ∈ R Nx×Ny , we define

(∂ 4,h x v h ) I(i,j) = 1 ∆x 4 v I(i-2,j) -4v I(i-1,j) + 6v I(i,j) -4v I(i+1,j) + v I(i+2,j) , (18a)
(∂ 3,h x,+ v h ) I(i,j) = 1 ∆x 3 -v I(i-1,j) -3v I(i,j) + 3v I(i+1,j) + v I(i+2,j) , (18b)
(∂ 3,h x,-v h ) I(i,j) = 1 ∆x 3 -v I(i-2,j) -3v I(i-1,j) + 3v I(i,j) + v I(i+1,j) , (18c)
(∂ h x,+ v h ) I(i,j) = 1 ∆x v I(i+2,j) -v I(i+1,j) , (∂ h x,-v h ) I(i,j) = 1 ∆x v I(i-1,j) -v I(i-2,j) , (18d)
(∂ 2,h y v h ) I(i,j) = 1 ∆y 2 v I(i,j-1) -2v I(i,j) + v I(i,j+1) , (18e)
(∂ h y,+ v h ) I(i,j) = 1 ∆y v I(i,j+1) -v I(i,j) , (∂ h y,-v h ) I(i,j) = 1 ∆y v I(i,j) -v I(i,j-1) ; (18f)
the same definitions being deduced for any vector wh ∈ R Nx with for instance

(∂ h,4 x wh ) i = 1 ∆x 4 wi-2 -4 wi-1 + 6 wi -4 wi+1 + wi+2 . ( 19a
)
At last, the mean value constraint is discretized as

Π h (v h ) = Ny-1 j=2 v I(i,j) . (19b)
Next, we define the approximation of the different system and state properties related to the conditioning of their associated matrices.

Definition 4.1 (Discrete Duality Based system). Let Φh = ( φh ) i denote the approximation of φ(x i ) for i ∈ I, and Φ ′ h = (Φ ′ h ) I(i,k) that of φ ′ (x i , y j ) for (i, j) ∈ I × J . The set of equations [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF] are discretized by the linear system

A DB Φ(k+1) h = F DB h + B DB ′ (Φ ′ h ) (k) , (20a)
A DB ′ (Φ ′ h ) (k + 1) = F DB ′ h + B DB Φ(k+1) h . ( 20b
)
Definition 4.2 (Discrete Two Field Iterated, Micro Macro and Singular Perturbation systems). Let Φ h = (φ h ) I(i,k) be the approximation of φ(x i , y j ) for (i, j) ∈ I × J , and

Q h = (Q h ) I(i,k
) the approximation of q(x i , y j ) for (i, j) ∈ I × J . The Two Field Iterated system [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF] is discretized thanks to the linear systems

A TFI,φ Φ (k+1) h = F TFI,φ h + B TFI,q Q (k) h , (21a)
A TFI,q Q (k+1) h = F TFI,q h + B TFI,φ Φ (k+1) h + C TFI,q Q (k) h . ( 21b
)
The set of equations associated to the Micro Macro formulation [START_REF] Dimarco | Asymptotic preserving Implicit-Explicit Runge-Kutta methods for non linear kinetic equations[END_REF] gives rise to the linear system

A MM Φ h Q h = F MM h 0 (22a)
Finally the Singular Perturbation problem (4) is discretized by the following linear system

A SP Φ h = F SP h . (23)
The detailed expressions of all the matrices introduced in these definitions will be omitted. They rely on the discrete partial differential operators defined by [START_REF] Filbet | Numerical study of a nonlinear heat equation for plasma physics[END_REF][START_REF] Galassi | Flux expansion effect on turbulent transport in 3d global simulations[END_REF] and refer to Appendix B for more details. We now state some properties of these matrices.

• M nn is semi-positive definite with (λ i ) i=3,...,Nx-2 ≥ 0 , and
   min i λ i = 0 , max i λ i ∼ 1 ∆x 4 ; (24a)
• M rb,ε is positive definite for ε > 0 and semi-positive definite for ε = 0 with (ℓ ε k ) k=2,...,Ny-1 ≥ 0 , and

   min k ℓ ε k ∼ ε , max k ℓ ε k ∼ 1 ∆y 2 ; (24b)
• M DB ′ is definite positive:

(ℓ DB ′ k ) k=2,...,Ny-1 > 0 ,    min k ℓ DB ′ k ∼ 1 , max k ℓ DB ′ k ∼ 1 ∆y 2 . (24c)
This makes possible to precise the condition number of the matrix issued from the Singular Pertubation problem discretization.

Proposition 4.1. The eigenvalues of the matrix A SP denoted (Λ SP i,k ) (i,k)∈{3,...,Nx-2}×{2,...,Ny-1} are defined thanks to (λ i ) i=3,...,Nx-2 and (ℓ ε k ) k=2,...,Ny-1 as

Λ SP i,k = λ i + 1 ε ℓ ε k , (i, k) ∈ {3, . . . , N x -2} × {2, . . . , N y -1}.
The matrix A SP is positive definite for ε > 0, and its corresponding condition number is stated as follows

Cond(A SP ) = max i λ i + 1 ε max k ℓ ε k 1 ε min k ℓ ε k ∼ 1 (∆x) 4 + 1 ε β SP (∆y) 2 , (25)
β SP being a constant uniformly bounded with respect to ε, ∆x and ∆y.

Note that the matrix A T F I,φ can be obtained by substituting ε by ε 0 . This gives rise to the following result.

Corollary 4.1. The eigenvalues of the matrix A T F I,φ denoted (Λ T F I,φ i,k

) (i,k)∈{3,...,Nx-2}×{2,...,Ny-1} Λ T F I,φ i,k = λ i + 1 ε 0 ℓ ε0 k , (i, k) ∈ {3, . . . , N x -2} × {2, . . . , N y -1}.
The condition number of the matrix A T F I,φ is

Cond(A T F I,φ ) = max i λ i + 1 ε 0 max k ℓ ε0 k 1 ε 0 min k ℓ ε0 k ∼ 1 (∆x) 4 + 1 ε 0 β T F I,φ (∆y) 2 , (26)
β T F I,φ being a constant uniformly bounded with respect to ε 0 , ∆x and ∆y.

Proposition 4.2. The eigenvalues of the matrix A DB ′ denoted Λ DB ′ i,k , (i, k) ∈ {3, . . . , N x -2}×{2, . . . , N y -1} can be expressed thanks to (λ i ) i=3,...,Nx-2 and (ℓ DB ′ k ) k=2,...,Ny-1 thanks to

Λ DB ′ i,k = 1 (∆x) 4 λ i + 1 ε 1 (∆y) 2 ℓ DB ′ k , (i, k) ∈ {3, . . . , N x -2} × {2, . . . , N y -1}.
The matrix A DB ′ is positive definite for ε > 0, and its corresponding condition number is stated as follows

Cond(A DB ′ ) = max i λ i + 1 ε max k ℓ DB ′ k 1 ε min k ℓ DB ′ k ∼ ε (∆x) 4 + β DB ′ (∆y) 2 , (27)
β DB ′ being a constant uniformly bounded with respect to ε, ∆x and ∆y.

The demonstration of these results is provided in Appendix B. This series of results is concluded by two conjectures. The first one originates from the analysis carried out in [START_REF] Lozinski | Numerical analysis of an asymptotic-preserving scheme for anisotropic elliptic equations[END_REF]. The inf-sup condition issued in this work demonstrate that the condition number of the matrix stemming from the Micro Macro system can be estimated thanks to the eigen systems of the perpendicular and parallel operators.

(A M M ) ∼ 1 (∆x) 4 1 (∆y) 2 . ( 28
)
The second conjecture is related to the matrix A T F I,q stating properties similar to that of A T F I,φ .

Conjecture 4.2. The matrix A T F I,q is positive definite with

Cond(A T F I,q ) ∼ 1 (∆x) 4 + 1 ε 0 β T F I,q (∆y) 2 , (29)
β T F I,q being a constant uniformly bounded with respect to ε 0 , ∆x and ∆y.

From [START_REF] Klar | An Asymptotic Preserving Numerical Scheme for Kinetic Equations in the Low Mach Number Limit[END_REF], it is clear that the condition number of A SP is not bounded with respect to ε. Conversely, from Corollary 4.1, we conclude that the condition number of the matrix A T F I,φ is uniformly bounded with respect to ε, for any given ε 0 . This property also holds true for the matrix A T F I,q . This will be assessed by numerical investigations in the next section. Finally, the Asymptotic-Preserving property of the Duality Based formulation is also manifest thanks to Prop. 4.2 and noting that the problem (9) providing φ does not depend on ε.

Numerical investigations of the schemes properties

Setup definition

An analytic solution is manufactured in order to validate the numerical methods by comparison with the exact solutions. The following solution

φ(x, y) = cos (2πx) (1 + ε sin (2πy)) , (30)
is used to compute analytically the right hand side of Eqs. (4a), (4d) and (4e). The precision of a numerical method is evaluated by computing the error between the numerical approximation Φ h = ((φ h ) i,j ) and the exact solution evaluated in (x i , y j ) in a norm specified in the sequel.

System matrices condition number

The conditioning of the system matrices issued from the discretization of the numerical methods introduced in Sec. 3 is analyzed numerically and compared to that of the original problem and to their analytic expression.

First the deterioration of the conditioning of the matrix A SP issued from the Singular Perturbation problem discretization is illustrated on Fig. 1. The first investigations illustrated on Fig. 1(a) are carried out with meshes equally refined in both directions and for strong anisotropies (ε ≤ 10 -6 ). The condition number of the matrix A SP is observed to be dominated by the aligned differential operator and scales as 1/(∆y 2 ). The meshes used for these computations are not refined enough (in the perpendicular direction) for the transverse operator to influence the matrix condition number. The outputs plotted on Fig. 1(b) are related to meshes with different resolutions along the parallel and perpendicular directions. The number of nodes in the perpendicular direction is denoted N x , N y being that of the parallel direction. The outputs of this figure are carried out by setting either N x or N y to a prescribed value and then increasing the resolution in the other direction. The behaviour of the condition number derived in Prop. 4.1 is recovered from these numerical investigations. For meshes with a coarse resolution along the aligned coordinate but a larger number of nodes in the perpendicular direction, the condition number is influenced by the fourth order operator and scales as 1/(∆x 4 ). The condition number of the matrix A SP increases linearly with ε -1 with a condition number hardly tractable for the targeted anisotropy strengths (ε < 10 -6 ) even with coarse meshes. This motivates the needs to develop numerical methods with condition numbers bounded with respect to ε.

The properties of the matrix A M M issued from the discretization of the Micro Macro formulation are ivestigated on Fig. 2. The Asymptotic-Preserving property of the methods is also verified thanks to the computations of Fig. 2(a), the condition number being unaffected by the anisotropy strength value (ε = 10 -10 and ε = 10 -20 ). The plots for intermediate anisotropies (ε ∼ 10 -6 ) are omitted since they are very similar to those actually presented. The bad scaling of the Micro Macro matrix condition number with respect to the mesh refinement is particularly manifest with the discretization of high order differential operators. The condition number is shown to deteriorate rapidly with the grid size. Indeed, if h denotes the typical mesh size, the conditioning of the matrix is proportional to 1/h 6 . The variations of the condition number with respect to a mesh refinement either in the direction perpendicular or parallel to the magnetic field is illustrated in Fig. 2(b). These results are in perfect agreement with Conjecture 4.1 stating a condition number proportional to both (1/∆x 4 ) and (1/∆y 2 ). For the Duality Based formulation, only the matrix providing the fluctuating part of the solution, namely A DB ′ is analysed here. Indeed, the mean part φ is a one dimensional unknown (to be compared to the two dimensional component φ ′ ), therefore the numerical cost needed for its computation is negligible. The plots relating the properties of the matrix A DB ′ are presented on Fig. 3. The condition number of the matrices at hand are computed for isotropic as well as anisotropic problems with the following trends emerging. For a given value of N y the condition number is proportional to both ε and 1/(∆x) 4 (see Fig. 3(b)). Therefore the deterioration of the matrix conditioning due to the fourth order discrete operator may be offset by the anisotropy: the condition number of the matrix improves with the anisotropy strength. For severe anisotropies, defined by ε = 10 -10 for these computations, the influence of the fourth order differential operator on the condition number is completely balanced by the anisotropy strength and only depends on 1/(∆y) 2 (see also Fig. 3(a)). Conversely, for vanishing anisotropies (ε ∼ 1), the variations of the conditioning are totally correlated with the discretized fourth differential operator. This is totally in line with the analytic expression stated in Prop. 4.2

The condition number of the matrices A T F I,φ and A T F I,q solved for the Two Field Iterated formulation is analysed on Figs. 4 and5. The asymptotic preserving property is manifest also for this formulation, with both matrices conditioning uniformly bounded with respect to the values of ε as depicted by the plots of Fig. 4. The condition numbers of both matrices scale as 1/(∆x) 4 when the mesh is refined in the direction perpendicular to the magnetic field and to 1/(∆y) 2 when the mesh is refined along the aligned direction. The condition number of these matrices deteriorates with the vanishing of ε 0 , the conditioning being proportional to 1/(ε 0 ∆y 2 ) as stated by Conjecture 4.2 and illustrated by the plots of Fig. 5. This latter property is similar to that of the matrix A SP issued from the discretization of the Singular Perturbation problem (with a trend equal to 1/(ε∆y 2 )). However, for the typical values of the anisotropy strength ε = 10 -8 -10 -6 and ε 0 (10 -3 ) the condition number of the TFI formulation is improved by three to five order of magnitudes.

A synthesis of the condition numbers for the different numerical methods is proposed in Tab. 1. The analytic estimates are gathered with two asymptotic regimes. The first one provides the condition number in the limit of infinite anisotropy. This provides a clear overview of the Asymptotic-Preserving property. The second regime is related to mild anisotropies but with refined meshes and is somehow more representative of the applications. From these outputs, it appears that either the discretized Singurlar-Perturbation problem or the Micro Macro formulation are not well suited to address severe anisotropies on refined meshes. The Duality Based as well as the Two Field Iterated formulations define better choices in this context.

Convergence of the fixed-point iterations

The properties of the different methods with respect to the system matrices condition number show clear advantages to use either the Duality Based or the Two Field Iterated formulations. Both methods Fig. 5. Condition number of the matrix A T F I,q providing the auxiliary variable q in the TFI scheme as a function of ε 0 for different refinement of the perpendicular and parallel mesh sizes.

define an iterative computation of the solution. For the Duality Based method a direct computation of both component may be proposed. However, solving each problem at once gains different advantages outlined in the next lines. The system providing φ is a one dimensional problem. Its numerical resolution is costless compared to that of the fluctuation φ ′ which is a two dimensional problem. This property has been harnessed in [START_REF] Yang | Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths[END_REF] to improve the computational efficiency of the Duality Based method and render its cost comparable to, or lower than, the resolution of an isotropic problem. Another advantage should be outlined here. The system matrix providing Φ ′ h has a condition number scaling as 1/h 2 . This means that the two dimensional problem resolution required for the iterative process offers a better conditioning scaling than a full system matrix providing ( Φh , Φ ′ h ) which is expected to be at least that of the system providing Φh , hence 1/h 4 . The only linear system with a condition number scaling as 1/h 4 is limited to the resolution of the one dimensional problem giving Φh in the iterative resolution proposed herein.

Table 1. Condition number of matrices issued form the discretization of the Singular problem (SP) as well as the Duality Based, Micro Macro and Two Field iterated formulations: General expression of the condition number, asympotics for vanishing mesh steps (h := ∆x = ∆y → 0, assuming mild anisotropies) as well as for vanishing anisotropies ε → 0.

SP

MM DB DB ′ TFI,φ TFI,q Cond. 1 (∆x) 4 + 1 ε β SP (∆y) 2 1 (∆x) 4 (∆y) 2 1 (∆x) 4 ε (∆x) 4 + βDB ′ (∆y) 2 1 (∆x) 4 + 1 ε0 β T F I,φ (∆y) 2 1 (∆x) 4 + 1 ε0 β T F I,q (∆y) 2 ε → 0 1 ε 1 h 2 → ∞ 1 h 6 1 h 4 1 h 2 1 h 4 1 h 4 h → 0 1 h 4 1 ε 1 h 2 1 h 6 1 h 4 1 h 2 ε h 4 1 h 4 1 h 4
The questions of interest now are therefore: how fast the iterations defined by those schemes are converging and how to easily control the convergence of the iterates ?

The control of the fixed point iterations is managed thanks to a stopping criterion. It consists in monitoring the error between two consecutive iterates of the solution related to the norm of the iterate and to consider that the convergence is reached, when

Φ (k+1) h -Φ (k) h ∞ Φ (k+1) h ∞ ≤ Tol , (31)
k denoting the iterate number and Tol ≪ 1 is a given value prescribing the desired precision. In Table 2, the influence of this tolerance parameter on the approximation precision is investigated for the Two Field Iterated method. This table relates the number of iterations required to meet the stopping criterion defined by Eq. ( 31) with the values of the parameter Tol specified in first column. The number of iterations hardly depends on ε (therefore only the computations carried out with ε = 10 -8 are presented), but increases with tighter precision requirements (smaller values of Tol). Note that, imposing a too restrictive convergence criterion induces a large number of iterations without improving the quality of the numerical approximation.

Over the full set of computations performed, the optimal value is Tol = 10 -7 for the most refined meshes. For this value the minimal error norm is obtained in a moderate number of iterations (less than 10) for ε 0 ∼ 10 -3 . The same conclusion, regarding the choice of Tol, can be drawn from the investigations carried out with the Duality Based formulation as reported on Tab. 3. The convergence of the fixed point iterations is observed to be very fast, within a couple of iterations to a tight convergence (Tol=10 -6 ) for intermediate mesh sizes. Furthtermore, the number of iterations to convergence does not depend on anisotropy strength. Computations carried out on more refined meshes (with 1000 × 1000 nodes) show that the convergence rate is hardly affected by the mesh size. To reach the tightest precision on those refined meshes the number of fixed point iterations is marginally increased to at most 4 to 5. This section is concluded by an analysis of the influence of the parameter ε 0 on the convergence rate of the TFI method. Indeed, this parameter plays a central role for the method efficiency. Therefore, the convergence of this formulation with respect to different sets of parameters (ε, ε 0 ) is investiagted in Table 4. It appears that small ε 0 -values increase the convergence rate of the fixed point iterations and deteriorates the condition number of the matrices. The choice of this parameter is thus a trade off between fast convergence of the fixed iterations and a tractable condition number for the inner problems resolution. Note that, for ε 0 = ε, the equation ( 13) is equivalent to the original problem (1) hence the convergence in one iteration. The value ε 0 = 10 -3 seems to define the best choice and is selected in the sequel. The convergence is secured in no more than 6 iterations and the matrices condition number is preserved whatever the anisotropy strength. This latter property is mandatory to extract the optimal precision from the scheme. Finally, we point out, that the properties obtained for this meshes remain quite the same for larger mesh size, the number of iterations required to reach convergence on a 1000 × 1000 mesh being increased to 8 iterations with ε 0 = 10 -3 .

A comparison of the precision and computational efficiency of the different methods

In this section, we shall focus on representative computations. To this end, the targeted anisotropy strength (ε = 10 -8 ) is selected and refined meshes will be used to proceed with the comparisons of the different methods. The Duality Based formulation [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF], the Micro Macro formulation [START_REF] Dimarco | Asymptotic preserving Implicit-Explicit Runge-Kutta methods for non linear kinetic equations[END_REF] and the Two Field Iterated formulation [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF] are compared against the discretized Singular Perturbation problem. The precision of the approximation carried out by the different methods is plotted on Fig. 6(c) for meshes refined equally in both directions, on Fig. 6(d) for meshes more refined in the perpendicular direction than in the parallel one. The condition number of the matrices issued from the different numerical methods is also plotted on Figs. 6(a) and 6(b) for the two mesh refinements. The linear systems are solved thanks to a direct sparse solver (MUMPS [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]) which also furnishes the estimate of the matrices condition number. The Duality Based as well as the Two Field Iterated methods offer a precision up to the second order accuracy of the space discretization. This is illustrated by a regular convergence rate on the plots of Fig. 6(c). Both the discretization of the Singular-Pertubation problem and the Micro Macro method can not preserve the accuracy of the numerical approximation. This is explained by the large condition number of the associated system matrices. The estimate provided on Figs. 6(a) and 6(b) are not accurate for these two methods. Indeed, the condition number of these matrices is too large to be computed accurately with the precision of the computer arithmetic. Therefore the scaling of the Micro Macro matrix conditioning (1/h 6 ) is not recovered on Fig. 6(a), contrary to the computations performed on coarsest meshes (see Fig. 2(a)). The same remark applies for the matrix issued from the distretization of the Singular Perturbation problem. The advantages of the Two Field Iterated and Duality Based schemes are manifest on these plots. The conditioning of the matrices issued from these two methods is manageable for very refined meshes. Note also the property of the Duality Based approach which gives rise to matrices with a very low condition number for the computation of the fluctuating (two dimensional) component. On the contrary, due to the bad conditioning of the matrices issued from the Micro Macro and the discretized Singular Perturbation problem, the precision of the computations carried out by either method can not be guaranteed as one observes on Fig. 6(d). This is even more sticking on Fig. 6(c) with a precision deteriorated for coarse meshes (from 400 × 400). Note that for meshes with more than 1000 nodes in the perpendicular direction, the condition number of the discrete fourth order differential operator is too large and deteriorates the precision of the numerical approximation for all the methods. This explains the plateau reached by the plots on Fig. 6(d) when refining the grid from 1000 × 200 to 1200 × 200 nodes.

The computational efficiency of these numerical methods are compared in Fig. 7 with the same mesh refinement strategies. The computational times are given relative to those of the Singular Perturbation problem. The Duality Based scheme is less than twice as expensive as the discretized Singular Perturbation scheme [START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF], the Micro Macro method (16) and the Two Field Iterated scheme [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF]. The computation time is relative to that of the SP problem resolution.

problem on the whole set of meshes. The Two Field Iterated schemes is roughly 2.5 times as expensive as the singular perturbation problem. Finally the Micro Macro method is the less efficient, more than 3 times as expensive. Note that, these computations are carried out with very conservative parameters and the efficiency of the iterative methods may be improved by more aggressive parameters. However the data collected in Fig. 7 give the right scales. Few concluding remarks can be stated at this point. First, when one coordinate is aligned with the magnetic field, it makes sense to dissociate the mesh step along the aligned and perpendicular directions. Since the gradients are expected to be small in the parallel direction it is likely that the mesh resolution remains coarse along the y-direction. Conversely, in the perpendicular direction, the variations of the solution call for a small mesh step to capture accurately the gradients in this direction. Therefore, the Micro Macro formulation may be competitive compared to standard methods. Indeed, if the refinement along the y-coordinate is coarse and the value of ε small enough (ε ≪ ∆x) , the Micro Macro formulation is likely to give rise to a system matrix with a better condition number than a discretization of the original system. However, in this context (adapted coordinates), it is much more efficient to use the Duality Based decomposition.

Second, for complex magnetic field geometries and coordinates unrelated to the anisotropy direction, the parallel and perpendicular directions are not oriented along the y or x axis. Therefore, there is no reason to assume that ∆x and ∆y have different values. The conditioning grows as h -6 (h = ∆x = ∆y) for the Micro Macro method (see Fig 2(a)) while standard discretizations provide a growth rate proportional to ∆x -4 + ε -1 β SP ∆y -2 . This feature explains why, though a direct discretization of the initial problem is not Asymptotic-Preserving, it may provide a better means of solving the problem numerically, for refined meshes and moderate anisotropy strengths. However, in this framework too, the Two Field Iterated method offers a more effective way to carry out precise numerical approximations for large anisotropies on refined meshes.

Third, to address more refined meshes a preconditioning strategy is mandatory to offset the bad scaling of the fourth order differential operator. This is routinely implemented for standard elliptic equations. The construction of preconditionner for the mildly anisotropic problems of the Two Field Iterated methods can be readily implemented using classical tools for elliptic problems. This is a straightforward task for the Duality Based method. We refer to [START_REF] Yang | Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths[END_REF] for conclusive demonstrations of this last two statements. In this regard, the Two Field Iterated method offers a decisive advantage over the Micro Macro approach. The discretization of the latter gives rise to an augmented system matrix, similar to those stemming from Saddle Point problems (see for instance [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]). However, the efficient construction of preconditionner for the linear system issued from the Micro Macro method remains an open question.

Towards plasma physics simulations

The model investigated in this work is simplified in order to analyze specifically the weaknesses of the AP methods proposed so far and outline their differences with respect to standard discretizations. It is important to note that the AP property makes possible to choose the discretization parameters independently of the anisotropy strength, hence the motivation for developing such numerical methods. Nonetheless, further developments are mandatory before addressing significant physics, they are examined in the following lines.

The first approximation of this study consists in considering a stationary problem. The derivation of time discretization compliant with the asymptotic-preserving property has been studied extensively (see for instance [START_REF] Pareschi | Implicit-Explicit Runge-Kutta Schemes and Applications to Hyperbolic Systems with Relaxation[END_REF][START_REF] Boscarino | Implicit-Explicit Runge-Kutta Schemes for Hyperbolic Systems and Kinetic Equations in the Diffusion Limit[END_REF][START_REF] Dimarco | Asymptotic preserving Implicit-Explicit Runge-Kutta methods for non linear kinetic equations[END_REF] or [START_REF] Lozinski | Highly anisotropic nonlinear temperature balance equation and its numerical solution using asymptotic-preserving schemes of second order in time[END_REF][START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF] for implementations in a context close to the one addressed herein). The class of time discretizations with L-stable properties guaranties an AP discretization. The present investigations can therefore be extended to time dependent problems thanks to these methods.

The topology of the magnetic field considered for these first investigations is also crude. Though an extension of the duality based approach to meshes unrelated to curved magnetic field lines is proposed in [START_REF] Degond | Duality-based asymptotic-preserving method for highly anisotropic diffusion equations[END_REF], this method is not efficient for non aligned coordinates (see the discussion in [START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF]). Actually, the use the Duality-Based approach for curved field lines should be considered together with adapted coordinates. Nevertheless, precedent works have also demonstrated the efficiency of AP methods to cope with complex magnetic field geometries with non adapted coordinates. The Micro Macro as well as the TFI methods have been introduced with this specific goal: dealing with meshes and coordinates unrelated to the magnetic field. In [START_REF] Narski | Asymptotic Preserving scheme for strongly anisotropic parabolic equations for arbitrary anisotropy direction[END_REF] [14] magnetic islands are for instance efficiently addressed thanks to Cartesian grids.

A last difficulty should be pointed out. In general, plasma physics is non linear and the formation of internal layers is a common feature in this context. Such an occurrence of local areas where the solution develops steep gradients is for instance dealt with in [START_REF] Deluzet | Numerical study of the plasma tearing instability on the resistive time scale[END_REF] using AP-methods. More sophisticated numerical methods remain yet to be developed along with AP schemes to address efficiently such problems: either local mesh refinements in the framework of unstructured grids, or patch techniques to perform local mesh refinement for Cartesian meshes [START_REF] Lozinski | Numerical zoom for advection diffusion problems with localized multiscales[END_REF].

Conclusion

In this document two new numerical methods have been introduced for the numerical approximation of the so-called vorticity equation. This is an anisotropic equation with unprecedented difficulties that are two-fold. First, the aligned differential operator has a kernel which is not reduced to the null function. Therefore, for severe anisotropies (ε ≪ 1), the system matrix issued from the discretization of this equation exhibits a large condition number preventing from computing an accurate approximation. Second, the perpendicular dynamic is accounted for by a fourth order differential operator. This causes the ineffectiveness of Asymptotic Preserving methods classically operated to overcome the difficulty originating from the aligned differential operator. The analysis conducted within this document demonstrates that the Micro Macro schemes proposed in [START_REF] Angot | Asymptotic-Preserving Methods for an Anisotropic Model of Electrical Potential in a Tokamak[END_REF][START_REF] Mentrelli | Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modeling[END_REF] give rise to system matrices with a condition number bounded with respect to ε, but scaling as 1/h 6 , h being the typical mesh step. This destroys the advantage of this approach, compared to a straight discretization of the original problem, as soon as the mesh is mildly refined. This calls for the derivation of new numerical methods, as the one introduced in this work. The first approach implements a decomposition of the two dimensional solution into a (one dimensional) mean component corrected by a (two dimensional) fluctuating part. This method, designed for anisotropy aligned to one coordinate, is very efficient. Indeed, it gives rise to a system matrix with a condition number scaling as 1/h 4 for the one dimensional problem and 1/h 2 for the two dimensional problem. A second method is proposed, operating a decomposition of the solution into two (two dimensional) components. This method can address arbitrary anisotropy, unrelated to both the coordinate system and the mesh. The condition number of the system matrices issued for each component scales as 1/h 4 and is bounded with respect to the anisotropy strength. These methods are proved to be effective for the computation of the vorticity equation solution for arbitrary anisotropy strength and refined meshes. These two methods are also well suited to address three dimensional problems with iterative resolutions of the associated linear systems. These extensions, together with the possibility to cope with complex magnetic field geometries, transient problems are mandatory to address significant physics. Therefore they define the next steps for this work.

Multiplying by ν i and integrating by parts yields

λ i = a (ν i , ν i ) a (ν i , ν i ) + ε 0 a ⊥ (ν i , ν i ) + ε 0 1 0 ν i (x, 1) 2 + ν i (x, 0) 2 dx .
This proves that, on the one hand λ i ∈ [0, 1] and, on the other hand, ν 0 the eigenfunctions associated to the eigenvalue 0 are in G.

Similarly, the eigenvalues of the operator A -1 ε0,q ∆ are given by

λ i = a (ν i , ν i ) a (ν i , ν i ) + ε 0 a ⊥ (ν i , ν i ) .
Here again λ i ∈ [0, 1] with λ i = 0 for any function belonging to G.

The elementary properties of the operators involved in the iteration operator allow to bound its spectrum. Let us first remark, that the eigenvalues of the operator

ε ε 0 - ε -ε 0 ε 0 A -1 ε0,φ ∆
are larger than ε/ε 0 > 0 and smaller than 1 since the spectrum of A -1 ε0,φ ∆ belongs to the interval [0, 1]. The operator

A -1 ε0,q ∆ ε ε 0 - ε -ε 0 ε 0 A -1 ε0,φ ∆
is therefore a product of two operators with the spectrum in the interval [0, 1] and as a consequence, the iteration matrix A I has all its eigenvalues in the same interval. Moreover, the eigenvalue of the iteration operator A I equal to one is associated to eigenfunctions in G, the kernel of ∆ . Following arguments similar to [START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF], we can prove that q n+1 -q n is orthogonal to the kernel of ∆ for any n ≥ 0 with respect to the semi-norm induced by the bilinear form a ε0,q . This proves the convergence of the iterative schemes.

Appendix B. Characterization of the discrete operators eigenvalues

In this section, the properties of the matrices issued from the discretization of the different formulations are investigated. This analysis is conducted thanks to the following result. Let (λ i ) i∈{1,...,n} be the eigenvalues of M 1 and (ℓ k ) k∈{1,...,m} that of the matrix M 2 , then the eigenvalues of the matrix A denoted by (Λ i,k ) (i,k)∈{1,...,n}×{1,...,m} verify Λ i,k = aλ i + bℓ k . The next steps consist in identifying the matrices M 1 and M 2 for the numerical methods introduced in Sec. 4 and characterizing their eigenvalues to highlight their properties. Let Φ h = (φ i,k ) (i,k)∈{3,...,Nx-2}×{2,...,Ny-1} and M nn a matrix of M Nx-4 (R), then the matrix vector product M nn Φ h defines a discretization of the differential operator ∂ 4

The condition number of the matrix A is given by

x φ, where The matrix M nn is semi-positive definite, as stated by the following lemma.

M nn := 1 ∆x 4                2
Lemma 5. The eigenvalues of the matrix M nn denoted (λ i ) i∈{3,...,Nx-2} are real and non negative:

λ i ≥ 0 , i ∈ {3, . . . , N x -2} . (B.2)
Proof of Lemma 5. Let λ i ∈ R and q i ∈ R Nx-4 be a couple of eigenvalue and eigenvector of the matrix M nn . The identity q T i M nn q i = λ i q T i q i gives λ i = 1 ∆x 4 1 q T i q i   (q i ) 3 -(q i ) 4 2

+

Nx-3 j=4 -(q i ) j-1 + 2(q i ) j -(q i ) j+1 ) 2 + (q i ) Nx-3 -

(q i ) Nx-2 2   ,
It is easy to find that the eigenvalue λ i is non-negative. Then taking constant eigenvector defined as q i = (c, • • • , c) T , for any c = 0, yields λ i = 0.

Considering now the differential operator along the y-coordinate, as the one involved in the definition of the Singular Perturbation problem (4), the Robin boundary conditions are discretized thanks to We can prove that the matrix M rb,ε is positive definite for any ε > 0, while M rb,0 is semi-positive definite. For severe anisotropies ε ≪ 1, the following estimation can be stated for the smallest eigenvalue

φ i,1 + φ i,2 2
min k ℓ ε k ∼ ε. (B.6)
Proof of Lemma 6. The property of the eigenvalues (ℓ ε k ) k∈{2,...,Ny-1} to be real and non negative for ε ≥ 0 is demonstrated using similar arguments to the ones developed for the proof of Lemma 5. The equations (B.4) and (B.5) are direct consequences of [START_REF] Yueh | Eigenvalues of several tridiagonal matrices[END_REF].

To derive the estimate of the smallest eigenvalue for ε ≪ 1, we assume θ ≪ 1. Expanding sin (N y -1)θ , sin (N y -2)θ and sin (N y -3)θ into Taylor series and truncating to second order, Eq. (B.5) yields θ 2 = 2β(ε) + β(ε) 2 (N y -3) N y -2 + β(ε)(N y -2) 2 -β(ε)(N y -2) + 1 3 β(ε) + 1 6 β(ε) 2 (N y -1) 3 . Now considering that β(ε) ≪ 1, we have

θ 2 ∼ 2β(ε) N y -2 ∼ 2ε∆y 2 .
Finally, the smallest eigenvalue of (B.4) is estimated as min k ℓ ε k = 4 ∆y 2 sin 2 θ 2 ∼ 1 ∆y 2 θ 2 ∼ ε.

The linear system associated with the discretization of the singular perturbation problem (4) can thus be recast as M nn Φ h + 1 ε Φ h (M rb,ε ) T = F SP h , The above lemmas imply directly the statements of Prop. 4.1.

The discretization system of Eq. ( 8) providing the fluctuation component in the Duality Based decomposition, can be written as Proof of Lemma 7. Owing to [START_REF] Yang | Iterative Solvers for Elliptic Problems with Arbitrary Anisotropy Strengths[END_REF], the matrix M DB ′ nn is positive definite, while D DB ′ is a non-negative diagonal matrix. Therefore, the matrix M DB ′ is positive definite. The eigenvalues defined by Eq. (B.8) are those of M DB ′ nn . This is an approximation all the more precise than β ≪ 1. This parameter is indeed very small, since ε ≪ 1 and ∆y ≪ 1.

M nn Φ h + 1 ε Φ h (M DB ′ ) T = F DB ′
The Proposition 4.2 is an outcome of Lemma 7.

Definition 4 . 3 (

 43 Discrete differential operators). We denote by • (λ i ) i=3,...,Nx-2 the eigenvalues of the matrix M nn discretizing the fourth order partial differential operator (with respect to x: ∂ 4 x ψ) equipped with Neumann boundary conditions (∂ x ψ = 0 , ∂ 3 x ψ = 0); • (-ℓ ε k ) k=2,...,Ny-1 the eigenvalues of -M rb,ε the finite differenced Laplacian (with respect to y: ∂ 2 y φ) equipped with the Robin boundary conditions ∂ n ψ+εψ = 0 (n being the outward normal to the domain); • (-ℓ DB ′ k ) k=2,...,Ny-1 the eigenvalues of -M DB ′ the finite differenced Laplacian (with respect to y: ∂ 2 y φ) equipped with Robin boundary conditions and implementing the zero mean constraint. These elementary matrices are defined in Appendix B by Eqs. (B.1), (B.3), (??) and (B.7). A first result related to these matrices can be stated. Lemma 1 (Eigenvalues of the elementary matrices). The elementary matrices introduced in Definition 4.3 satisfy the following properties:

Conjecture 4 . 1 .

 41 The condition number of the matrix A M M satisfies, for ε ≪ 1 Cond

2 Slope 4 (

 24 Condition Number as a function of N = Nx = Ny. y = 30, ε = 10 -6 ) SP (N x = 30, ε = 10 -6 ) SP (N x = 30, ε = 10 -8 ) Slope b) Condition Number as a function of N = Nx (Ny = 30) or N = Ny (Nx = 30).

Fig. 1 .Fig. 2 .

 12 Fig. 1. Condition number of the matrix A SP issued from the discretization of the Singular Perturbation problem as a function of the grid size (N ) for different anisotropy strengths ε -1 .

7 DB

 7 ′ (ε = 1) DB ′ (ε = 10 -2 ) DB ′ (ε = 10 -10 ) Slope 2 (a) Condition Number as a function of N = Ny (Nx = 10).

15 DB

 15 ′ (ε = 1) DB ′ (ε = 10 -2 ) DB ′ (ε = 10 -10 ) Slope 4 (b) Condition Number as a function of N = Nx (Ny = 10).

Fig. 3 .

 3 Fig. 3. Condition number of the matrix A DB ′ issued from the discretization of the fluctuation equation (8) of the Duality Based formulation as functions of the parallel (left) and perpendicular (right) mesh sizes (N ) for different anisotropy strengths (ε -1 ).

  Condition number as a function of N = Ny (Nx = 10) for ε = 10 -20 .

  Condition number as a function of N = Nx (Ny = 10) for ε = 10 -20 .

Fig. 4 .

 4 Fig. 4. Condition number of the Two Field Iterated formulation (ε 0 = 10 -3 ) as a function of the perpendicular and parallel mesh sizes.

  Condition number for different values of Nx (Ny = 20).

  Condition number for different values of Ny (Nx = 20).

Table 2 .× 10 - 5 10 - 5 5 10 - 6 5 10 - 7 5 10 - 8 5 10 - 9

 2105556575859 Control of the TFI fixed-point iterations by the criteria defined by Eq. 31: Its is the number of fixed point iterations, Error the error in ℓ∞-norm. The computations are carried out on a mesh with Nx = Ny = 300 with ε = 10 -8 . 25 1.31 × 10 -4 4 7.81 × 10 -5 2 7.51 × 10 -40 7.99 × 10 -5 6 7.52 × 10 -5 3 7.50 × 10 -54 7.55 × 10 -5 8 7.50 × 10 -5 3 7.50 × 10 -72 7.50 × 10 -5 10 7.50 × 10 -5 4 7.50 × 10 -91 7.50 × 10 -5 19 7.50 × 10 -5 18 7.50 × 10 -5

  Mesh size (N = N x = N y ).

  Mesh size (N x ). (b) Cond. numbers as functions of Nx for Ny = 200. N = N x = N y ).

( d )

 d Precision as a function of Nx for Ny = 200.

Fig. 6 .

 6 Fig. 6. Condition number of the system matrices issued from the Singular perturbation problem, Duality Based (DB and DB ′ ) formulation (8-9), Micro Macro (MM) formulation (16) and Two Field Iterated formulation (13-14) as a function of the mesh step as well as error in ℓ∞-norm. These computations are carried out with ε = 10 -8 for meshes with either N = Nx = Ny with N ∈ [200, 1000] or Nx ∈ [200, 1200] and Ny = 200.

8 .Fig. 7 .

 87 Fig. 7. Computational time for the resolution of the Singular Pertubation (SP) problem, Duality Based (DB) scheme[START_REF] Crouseilles | Comparison of Numerical Solvers for Anisotropic Diffusion Equations Arising in Plasma Physics[END_REF][START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF], the Micro Macro method (16) and the Two Field Iterated scheme[START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF]. The computation time is relative to that of the SP problem resolution.

Lemma 4 .

 4 Let M 1 ∈ M n (R) and M 2 ∈ M m (R) be two square matrices, with n and m two integers, a, b ∈ R and a, b > 0. Then we define Ψ := aM 1 Φ + bΦM T 2 , where Ψ, Φ are matrices of ∈ M n,m (R). Reshaping Ψ and Φ into vectors Ψ ′ and Φ ′ ∈ R n×m respectively, then there exists a unique matrix A ∈ M n×m (R) such that Ψ ′ := AΦ ′ .

Cond 2 (

 2 A) = a max i λ i + b max k

  giving rise to the following identities:φ 1,k = φ 4,k , φ 2,k = φ 3,k , φ Nx-3,k = φ Nx,k , φ Nx-2,k = φ Nx-1,k , k = 2, . . . , N y -1.

- 1 ε φ i, 2 - 2 .

 22 φ i,1 ∆y = (g 0 ) i , φ i,Ny-1 + φ i,Ny 2 + 1 ε φ i,Ny -φ i,Ny-1 ∆y = (g 1 ) i , i = 3, . . . , N x -2,or equivalently, with α(ε) := ε∆y 2φ i,1 + α(ε) -1 α(ε) + 1 φ i,2 = ε∆y α(ε) + 1 (g 0 ) i , φ i,Ny + α(ε) -1 α(ε) + 1 φ i,Ny-1 = ε∆y α(ε) + 1 (g 1 ) i , i = 3, . . . , N x -The matrix associated with the second order differential operator is denoted M rb,ε and defined as M rb,ε :

Lemma 6 . 4 )

 64 The eigenvalues (ℓ ε k ) k∈{2,...,Ny-1} of the matrix M rb,ε defined by Eq. (B.3) are real and non negative. Furthermore, they satisfy the following property with θ = mπ, m ∈ Z, such thatsin (N y -1)θ -2 1 -β(ε) sin (N y -2)θ + 1 -β(ε)2 sin (N y -3)θ = 0. (B.5)

h, 7 )Lemma 7 . 2

 772 where the matrixM DB ′ ∈ M Ny-2 (R) is defined thanks to M DB ′ nn ∈ M Ny-2 (R) and D DB ′ ∈ M Ny-2 (R) as M DB ′ nn := M DB ′ nn + D DB ′It is possible to characterize the eigenvalues of M DB ′ in the asymptotic β → 0. The matrix M DB ′ is positive definite, with eigenvalues satisfying when (ε∆y) ≪ 1 sin 2 π(k -1) 2(N y -2) , k = 2, . . . , N y -

Table 3 .

 3 Number of fixed-point iterations (Its) for the Duality Based method: Its is the number of fixed point iterations, Error the error in ℓ∞-norm. The computations are carried out on a mesh with Nx = Ny = 300 with ε = 10 -8 .

			Tol= 10 -5		Tol= 10 -6		Tol= 10 -7	Tol= 10 -8
	ε	Its Error	Its Error	Its Error	Its Error
	10 -4	2	7.50 × 10 -5 2	7.50 × 10 -5 2	7.50 × 10 -5 2	7.50 × 10 -5
	10 -8	2	7.50 × 10 -5 2	7.50 × 10 -5 2	7.50 × 10 -5 2	7.50 × 10 -5
	10 -15 2	7.50 × 10 -5 2	7.50 × 10 -5 2	7.50 × 10 -5 2	7.50 × 10 -5

Table 4 .

 4 Number of fixed point iterations and error in ℓ∞-norm for the Two Field iterated formulation on different meshes. The stopping criteria is set to 10 -6 (Divergence is reported by X). Test case with Nx = Ny = 300.

Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work has been supported by the french "Agence Nationale pour la Recherche (ANR)" in the frame of the contract ANR-11-MONU-009-01 "MOONRISE: MOdels, Oscillations and NumeRIcal SchEmes" (2015-2019) as well as the "labex CIMI" (International Centre for Mathematics and Computer Science in Toulouse) in the frame of the project "SCANISO: SCalable solvers for ANISOtropic equations arising in magnetized plasma simulations" (2017-2019). Support from the "Fédération de Fusion pour la Recherche par Confinement Magnétique" (FrFCM) in the frame of the project "NEMESIA: Numerical mEthods for Macroscopic models of magnEtized plaSmas and related anIsotropic equAtions" (2018-2019) is also acknowledged. CY is partially supported by National Natural Science Foundation of China (Grant No. 11401138) and FD acknowledges invitations as a scholar professor by Harbin Institute of Technology in 2016, 2017 and 2018. The authors acknowledge fruitful discussions with Pierre Jolivet and Luc Giraud.

Appendix A. Proof of convergence of the fixed point iterations

The convergence of the sequence (φ k , q k ) k≥0 is proved in this section. We only provide the main results for the proof of Thoerem 3.1 and refer to [START_REF] Deluzet | A two field iterated asymptotic-preserving scheme for highly anisotropic elliptic equations[END_REF] for the details.

For any function φ ∈ V,

(Ω)} we introduce the operators ∆ ⊥ and ∆ :

along with the bilinear forms defined for (u, v) ∈ V × V by

Moreover, for ε 0 ∈ R, ε 0 > 0, A ε0,φ = -ε 0 ∆ 2 ⊥ + ∆ denotes the fourth order operator corresponding to the system (13), supplemented with Robin boundary conditions and equivalent to the following bilinear form

where Γ denotes the part of the boundary with imposed Robin condition.

Finally, we introduce A ε0,q = -ε 0 ∆ 2 ⊥ + ∆ the fourth order operator corresponding to the system ( 14) supplemented with Neumann (instead of Robin) boundary conditions. This operator is invertible on the space of functions with a value fixed in one point of the domain:

The bilinear form associated to A ε0,q is

The proof of convergence of the iterative schemes relies on the following results.

Lemma 2 (Iteration operator). The iterative method defined by the equations (13-14) yields the following recurrence

where Tr(A -1 ε0,φ ) and g are the traces of the operator A -1 ε0,φ and the function φ on Γ, the operator A I being defined as

The eigenvalues of A I , denoted ℓ i , are real with ℓ i ∈ [0, 1]. The eigenfunctions associated to the largest eigenvalue ℓ i = 1 belong to G the kernel of the operator ∆ .

Lemma 3 (Elementary properties). The operators A ε0,φ and A ε0,q are invertible on V and V 1 respectively. The eigenvalues of the operators A -1 ε0,φ ∆ and A -1 ε0,q ∆ are real, non-negative and bounded by 1. The eigenfunctions ν 0 associated to the null eigenvalue of both operators belong to the kernel of the operator ∆ .

Proof. The eigenvalues of the operator A -1 ε0,φ ∆ are solution of the following problem: find λ i ∈ R and ν i ∈ V such that A -1 ε0,φ ∆ ν i = λ i ν i , that is to say