Chuanju Xu
email: cjxu@xmu.edu.cn

Minghui Li

Mejdi Azaiez

NEW EFFICIENT TIME-STEPPING SCHEMES FOR THE ANISOTROPIC PHASE-FIELD DENDRITIC CRYSTAL GROWTH MODEL

Keywords: 2010 Mathematics Subject Classification. Primary 74A50, 65M06, 65M12, 76T10, 65Z05 Phase-field, Dendritic crystal growth, Time-stepping schemes, Unconditional stability

whether they are published or not. The documents may come

Introduction

On one side, dendritic growth is a very common phenomenon in nature. We are all familiar with the way how trees grow by spreading branches and roots from the main trunk. This is where the name "dendritic" comes from, although the term "dendrite" itself is used to describe branched projections of neurons. On the other side, dendritic growth phenomena and the shapes of growing crystals are of fundamental interest to physicists and are of practical importance to engineers.

Crystal dendritic growth is one of the most extensively studied topics in the scientific literature.

Crystallization proceeds through the competition between thermodynamics -driven by the local undercooling of the liquid ahead of the solidification front -and the ability of the system to diffuse latent heat of fusion away from the solid-liquid interface. It usually forms natural fractal microstructure, so-called dendrites, which are the ubiquitous crystal form in freezing alloys and supercooled melts. When the molten material is supercooled below the freezing point of the solid, a spherical solid nucleus grows in the undercooled melt initially. Along with some preferred directions of growth, the solid form begins to express some protrusion accompanied by steeper concentration gradients at its end. Dendritic microstructures formed during solidification/freezing play a key role in properties of the final solid material. Understanding these microstructures is therefore considered essential for controlling basic solidification and crystal growth processes.

The first phase-field models were suggested for numerically simulating dendritic growth in 1980s; see, e.g., [START_REF] Caginalp | An analysis of a phase-field model of a free boundary[END_REF][START_REF] Collins | Diffuse interface model of diffusion-limited crystal growth[END_REF][START_REF] Fix | Phase-field methods for the free boundary problems[END_REF][START_REF] Langer | Models of pattern formation in first-order phase transitions[END_REF]. This concept has been validated by comparison with theoretical predictions and experimental measurements and is applied to a broad range of investigations in materials science [START_REF] Boettinger | Phase-field simulations of solidification[END_REF]. Nowadays, the phase-field method has emerged as a powerful tool for modelling and simulation of crystal dendritic growth. In contrast to sharp interface approaches with interfaces of zero thickness, the phase-field model introduces a smooth phase-field variable by a diffuse interface profile to distinguish between the solid and liquid phases. In this model the complicated topological changes of a solid-liquid interface can be handled in an easy way without the need of the explicit tracking of the interface. The phase field is considered as an order parameter which is introduced to describe the moving interfacial boundary between unstable and stable phases during phase transformation processes. By asymptotic expansions, it can be shown that the phase-field methods relate to classical sharp interface models such as Hele-Shaw type models and Stefan problems in the limit of zero interfacial thickness, see, e.g. [START_REF] Caginalp | Phase-field computation of single needle crystals, crystal growth, and motion by mean curvature[END_REF].

Another advantage of the phase-field approach is that the governing set of equations in the model can be naturally derived from an energy-based variational principle. The variational framework of phase-field formulations makes them thermodynamically consistent and physically attractive in modeling the general phase-field dendritic crystal growth model [START_REF] Dantzig | Dendritic growth morphologies in al-zn alloys-part ii: phase-field computations[END_REF][START_REF] Demange | A phase field model for snow crystal growth in three dimensions[END_REF][START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF][START_REF] Kim | Universal dynamics of phasefield models for dendritic growth[END_REF][START_REF] Marinozzi | Phase-field model for dendritic growth in a channel[END_REF][START_REF] Mullis | A study of kinetically limited dendritic growth at high undercooling using phasefield techniques[END_REF][START_REF] Ramirez | Examination of binary alloy free dendritic growth theories with a phase-field model[END_REF][START_REF] Warren | Prediction of dendritic growth and microsegregation patterns in a binary alloy using the phase-field method[END_REF].

In this paper, we will focus on the numerical approximations for the anisotropic phase-field dendritic crystal growth model proposed in [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF][START_REF]Phase-field model of dendritic sidebranching with thermal noise[END_REF]. The model is composed of two coupled nonlinear equations. One is the phase-field equation that governs the anisotropy of the crystal. The other is the heat equation that controls heat diffusion of the system. It is shown that this nonlinear coupled system satisfies a thermodynamically consistent energy dissipation law. The main aim of this paper is to design efficient numerical schemes for this nonlinear crystal growth model, which satisfies a discrete version of the energy dissipation law. In fact, constructing schemes that preserve the discrete energy dissipation law for similar models has been subject of many recent papers [1, 5-7, 12, 25, 26, 29, 39]. Although large amounts of works have been devoted to numerical approximation for phase-field dendritic crystal growth models; see, e.g., [START_REF] Li | Energy and entropy preserving numerical approximations of thermodynamically consistent crystal growth models[END_REF][START_REF] Li | Phase-field simulations of crystal growth with adaptive mesh refinement[END_REF][START_REF] Shah | Numerical simulation of two-dimensional dendritic growth using phase-field model[END_REF][START_REF] Zhao | Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach[END_REF] and the references therein, there is still a need for efforts on developing low-cost, stable, and high order schemes for such models.

The main difficulties in constructing highly efficient schemes for the dendritic crystal growth models come from: 1) the double-well energy potential and the stiffness associated with the interfacial width in the phase equation; 2) the anisotropic coefficient; 3) the nonlinear interaction terms in both the heat equation and the phase field equation. Let's briefly review recent progress in this direction. Firstly, to overcome the difficulty caused by the nonlinearity and the thin interface in the phase field equation, schemes based on the invariant energy quadratization (IEQ) [START_REF] Yang | Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends[END_REF] and scalar auxiliary variable approach (SAV) [START_REF] Shen | The scalar auxiliary variable (SAV) approach for gradient flows[END_REF][START_REF]The ieq and sav approaches and their extensions for a class of highly nonlinear gradient flow systems[END_REF] have been proposed: a decoupled stable but only first-order scheme in [START_REF] Zhang | A novel decoupled and stable scheme for an anisotropic phase-field dendritic crystal growth model[END_REF] and a second-order stable scheme but fully coupled scheme in [START_REF]Efficient linear, stabilized, second-order time marching schemes for an anisotropic phase field dendritic crystal growth model[END_REF].

Due to the presence of the nonlinear phase term in the heat equation and the interaction term in the phase equation, it seems not easy to design a fully-decoupled, second-order accurate and energy stable scheme. For example, the splitting method used in [START_REF] Zhang | A novel decoupled and stable scheme for an anisotropic phase-field dendritic crystal growth model[END_REF] is not directly extendable to a second-order discretization for the time derivative of the phase function in the heat equation. In the case of constant mobility parameter, [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF] proposes a decoupling, linear, second-order accurate, and unconditionally stable scheme by using multi-auxiliary variables. A similar technique was used in [START_REF]A new efficient fully-decoupled and second-order time-accurate scheme for cahn-hilliard phase-field model of three-phase incompressible flow[END_REF][START_REF]A novel fully decoupled scheme with second-order time accuracy and unconditional energy stability for the navier-stokes equations coupled with mass-conserved allen-cahn phase-field model of two-phase incompressible flow[END_REF][START_REF]Numerical approximations of the navier-stokes equation coupled with volume-conserved multi-phase-field vesicles system: Fully-decoupled, linear, unconditionally energy stable and second-order time-accurate numerical scheme[END_REF] to deal with some coupling models. However, compared with the traditional SAV approach, multi-auxiliary variable approach means extra computational cost since more equations are to be solved. In particular, the second-order scheme proposed in [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF] is a four-step scheme, thus is much more computationally expensive.

The main purpose of the present paper is to propose easy-to-implement, second-order accurate, and unconditionally stable schemes for the anisotropic phase-field dendritic crystal growth model.

First of all, we rewrite the time derivative term of the phase function in the heat equation into an equivalent form, which allows to design a three-step second-order scheme. The idea is to introduce a suitable auxiliary variable to the Allen-Cahn equation and a new technique to treat the coupling terms. Then some carefully chosen stabilization terms are added to result in totally decoupled schemes that satisfy a discrete energy law without losing the convergence order. A careful examination shows that the proposed schemes can be implemented by only solving four linear elliptic equations and a simple algebraic equation. As far as we know this is the first second-order scheme that is totally decoupled, linear, unconditionally stable for the dendritic crystal growth model variable mobility parameter. In the case of constant mobility parameter, compared with [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF] (a four-layer scheme that requires solve five linear elliptic equations with constant coefficients and some algebraic equations), our scheme is a three-layer scheme that only needs solve four linear elliptic equations with constant coefficients and a simple algebraic equation.

The rest of the paper is organized as follows. In Section 2, we describe the phase-field dendritic crystal growth model, and present the equivalent reformulation using auxiliary variables. In Section 3 we propose a first-order unconditionally stable time-stepping scheme, and prove the energy decay property of the proposed scheme. Section 4 is devoted to construct and analyze a second-order, linear, decoupled, and unconditionally stable scheme. The implementation detail is also presented to show that the scheme can be efficiently realized through solving a set of decoupled, linear elliptic equations. We give in Section 5 some numerical examples to verify the efficiency of the proposed methods. Finally, the paper ends with some concluding remarks.

Governing equations and auxiliary variable reformulation

2.1. Governing equations. We are interested in numerically solving the following anisotropic phase-field dendritic crystal growth model equations in the domain Ω ⊂ R 2 :

(φ)φ t = - δE δφ - λ ε h (φ)T, (2.1)
T t = ∇ • (D∇T ) + Kh (φ)φ t , (2.2)
where φ(x, t) is the phase function to label the liquid and solid phases, (φ) > 0 is the mobility parameter that can be chosen either as a constant [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF], or as a function of φ [START_REF]Efficient linear, stabilized, second-order time marching schemes for an anisotropic phase field dendritic crystal growth model[END_REF]. ε > 0 is a parameter used to control the interface width, λ is the linear kinetic coefficient. In Equation (2.2), T (x, t) is the scaled temperature, D is the constant diffusion rate of the temperature, and K is the latent heat parameter that controls the speed of heat transfer along with the interface. It is worth noting that the efficiency of the schemes we propose below covers the case D is a function of φ [START_REF] Suzuki | Phase-field model of dendritic growth[END_REF].

The function h(φ) is defined by

h(φ) := 1 5 φ 5 - 2 3 φ 3 + φ,
which represents a generation of latent heat. Following the phenomenological free energy used in [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF], we consider here

E(φ, T ) = Ω 1 2 κ 2 (∇φ)|∇φ| 2 + 1 ε 2 F (φ) + λ 2εK T 2 dx, (2.3)
where

F (φ) = 1 4 φ 2 -1 2 is the double-well type Ginzburg-Landau potential. κ(•) in (2.
3) is a function describing the anisotropic property, which takes the form [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF][START_REF]Phase-field model of dendritic sidebranching with thermal noise[END_REF]:

κ(∇φ) = 1 + σ cos(mθ), (2.4)
where m is a model number of anisotropy, σ is the parameter for the anisotropy strength, and θ = arctan φy φx . The variational derivative of E with respect to φ is:

δE δφ = -∇ • κ 2 (∇φ)∇φ + κ(∇φ)|∇φ| 2 H(φ) + f (φ) ε 2 ,
where H(φ) is the variational derivative of κ(∇φ), and f (φ) = F (φ). In the case m = 4, a direct calculation shows

H(φ) := δκ(∇φ) δφ = 4σ 4 |∇φ| 6 φ x φ 2 x φ 2 y -φ 4 y , φ y φ 2 x φ 2 y -φ 4 x . (2.5)
For convenience, we only consider that the equations (2.1) and (2.2) are subject to the Neumann boundary conditions

∂φ ∂n | ∂Ω = 0, ∂T ∂n | ∂Ω = 0, (2.6)
although other boundary conditions such as the periodic conditions are possible.

Now we briefly recall main property of (2.1)-(2.6).

A key property of the model is that it satisfies an energy law, which can be derived by taking the inner product of (2.1) with -φ t and using integration by parts:

d dt Ω 1 2 κ 2 (∇φ)|∇φ| 2 + 1 ε 2 F (φ) dx + Ω λ ε h (φ)T φ t dx = - Ω (φ)φ 2 t dx.
Then taking the inner product of (2.2) by -λ εK T gives:

d dt Ω λ 2εK T 2 dx - Ω λ ε h (φ)T φ t dx = - λD εK Ω ∇T • ∇T dx.
Combining the above two equalities gives the following energy law

d dt E(φ, T ) = - (φ)φ t 2 - λD εK ∇T 2 , (2.7)
where • denotes the standard L 2 (Ω) norm. This means that the energy E(φ, T ) decays in time during the crystal-growing process.

Auxiliary variable reformulation.

The main purpose of this paper is to develop novel efficient schemes for the anisotropic crystal growth model (2.1)-(2.6). We start with an auxiliary variable approach, which will be used later to construct time-stepping schemes for the phase field equation (2.1). We define the variable

R(t) = E 1 (φ), E 1 (φ) = Ω 1 2 κ 2 (∇φ) -S 1 |∇φ| 2 + 1 ε 2 F (φ) - S 2 2 φ 2 + B dx, (2.8)
where S 1 and S 2 are two positive constants, 0 < S 1 < (1 -σ) 2 , B is a positive constant used to make E 1 positive. Notice κ 2 (∇φ) ≥ (1 -σ) 2 and F (φ) is a quartic polynomial, one can verify that

Ω ( 1 2 κ 2 (∇φ) -S 1 |∇φ| 2 + 1 ε 2 (F (φ) -S 2 2 φ 2 )
)dx is bounded from below. Therefore such a constant B exists. The introduction of the constants S 1 and S 2 is inspired by the work [START_REF]Efficient linear, stabilized, second-order time marching schemes for an anisotropic phase field dendritic crystal growth model[END_REF][START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF][START_REF] Zhang | A fully decoupled, linear and unconditionally energy stable numerical scheme for a melt-convective phase-field dendritic solidification model[END_REF]. We will see that these constants help in ensuring the H 1 -stability of the phase function.

Using the auxiliary variable R(t), the total free energy (2.3) can be rewritten as

E(φ, R, T ) = Ω λ 2εK T 2 + S 1 2 |∇φ| 2 + S 2 2ε 2 φ 2 -B dx + R 2 , (2.9)
and the original equations (2.1)-(2.2) can be reformulated into the following equivalent form:

φ t = M (φ)µ, (2.10a) µ = - R(t) E 1 (φ) g(φ) + S 1 ∆φ - S 2 ε 2 φ - R(t) E 1 (φ) λ ε h (φ)T, (2.10b) R t = Ω g(φ) 2 E 1 (φ) φ t dx, (2.10c)
T t = D∆T + R(t) E 1 (φ) Kh (φ)M (φ)µ, (2.10d)
where M (φ) = 1 (φ) , and

g(φ) = -∇ • κ 2 (∇φ) -S 1 ∇φ + κ(∇φ)|∇φ| 2 H(φ) + 1 ε 2 (f (φ) -S 2 φ) .
Obviously, the equation (2.10c) can be obtained by taking the time derivative of the auxiliary variable R(t). The initial conditions for (2.10) take

φ| t=0 = φ 0 , T | t=0 = T 0 , R| t=0 = E 1 (φ 0 ). Remark 2.1.
To separate the computation of different unknown functions, one may think about the splitting method used in [START_REF] Zhang | A novel decoupled and stable scheme for an anisotropic phase-field dendritic crystal growth model[END_REF] or explicit treatment of φ t as in [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF]. However the former is unlikely to lead to a second-order scheme, while the latter may result in expensive four-layer computation [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF]. To construct more efficient second-order scheme, our idea here is to split (2. (2.10d) with λ εK T , then we perform integration by parts and sum up all equalities to get

d dt E(φ, R, T ) = - (φ)φ t 2 - λD εK ∇T 2 ≤ 0, (2.11)
where E(φ, R, T ) is defined in (2.9). Now we are in a position to construct and analyze our schemes for the anisotropic phase-field dendritic crystal growth model (2.1)-(2.6). To better follow the main argument, let's start with a first-order scheme.

A first order scheme and stability analysis

Let τ > 0 be the time step size, t n = nτ, 0 ≤ n ≤ N , T = N τ . We propose the following scheme:

assuming {φ n , T n , R n } are known, {φ n+1 , T n+1 , R n+1 } is computed by solving: φ n+1 -φ n τ = M (φ n ) µ n+1 - S 3 ε 2 (φ n+1 -φ n ) + S 4 (∆φ n+1 -∆φ n ) , (3.1a
)

µ n+1 = -ξ n+1 g(φ n ) + S 1 ∆φ n+1 - S 2 ε 2 φ n+1 -ξ n+1 λ ε h (φ n )T n , (3.1b) R n+1 -R n τ = 1 2 E 1 (φ n ) g(φ n ), φ n+1 -φ n τ - λ ε h (φ n )M (φ n ), µ n T n+1 -µ n+1 T n - λ ε h (φ n )M (φ n )T n , S 3 ε 2 (φ n+1 -φ n ) -S 4 ∆(φ n+1 -φ n ) , (3.1c)
T n+1 -T n τ = D∆T n+1 + ξ n+1 Kh (φ n )M (φ n )µ n , (3.1d
)

∂φ n+1 ∂n ∂Ω = 0, ∂T n+1 ∂n ∂Ω = 0, (3.1e)
where

ξ n+1 = R n+1 √ E 1 (φ n )
, S 3 and S 4 are two extra positive stabilization parameters.

Before carrying out the stability analysis, the scheme (3.1) is worthy of some explanation. First, we notice that the coupling terms in both the phase field equation and the temperature equation are treated explicitly. This is for ease of calculation. The implicit treatment of the coupling terms in the auxiliary variable equation, i.e., eq.(3.1c), may make the implementation difficult. However, as we will see in the next section, the extra terms added to the scheme, i.e., terms involving the parameters S i , play a dual role. On one side, some of extra terms are useful in decoupling the calculation of different unknowns. On the other side, some other extra terms help in enhancing the stability. For example, the term

S 3 ε 2 φ n+1 -φ n is used to balance the explicit treatment of 1 ε 2 f (φ)
in the phase field equation, and the term S 4 ∆ φ n+1 -φ n has purpose to balance the explicit treatment of the gradient term. The last point we want to emphasize is that the extra terms have the same order as the approximation to the time derivatives, thus do not affect the overall accuracy.

For example, we can check that the term

1 2 E 1 (φ n ) λ ε h (φ n )M (φ n ), µ n T n+1 -µ n+1 T n - S 3 ε 2 T n (φ n+1 -φ n ) + S 4 T n ∆(φ n+1 -φ n ) in (3.1c) is of order O(τ ).
Therefore, formally the convergence of the scheme (3.1) is first order.

3.1. Stability analysis. In the following theorem, we establish the stability result for the scheme (3.1). That is, we prove that a discrete "energy" decays in time, and consequently the discrete solution remains bounded during the time stepping.

Theorem 3.1. Let {φ n , T n , R n } be the solution of the discrete problem (3.1). Then the following discrete energy law holds:

E n+1 -E n = -Q n+1 -τ || (φ n ) φ n+1 -φ n τ || 2 + λD εK ||∇T n+1 || 2 , (3.2)
where E n is defined by

E n = S 1 2 ||∇φ n || 2 + S 2 2ε 2 ||φ n || 2 + λ 2εK ||T n || 2 + |R n | 2 ,
Q n+1 is given by

Q n+1 = S 1 + S 4 2 ||∇φ n+1 -∇φ n || 2 + S 2 + S 3 2ε 2 ||φ n+1 -φ n || 2 + |R n+1 -R n | 2 .
Proof By taking the inner product of (3.1a) with 2(φ n+1 -φ n ) M (φ n ) , and (3.1b) with 2(φ n+1 -φ n ), then summing up the resulting equations, we obtain

2 τ φ n+1 -φ n M (φ n ) 2 + 2S 3 ε 2 φ n+1 -φ n 2 + 2S 4 ∇(φ n+1 -φ n ) 2 + 2(ξ n+1 g(φ n ), φ n+1 -φ n ) +S 1 (||∇φ n+1 || 2 -||∇φ n || 2 + ||∇(φ n+1 -φ n )|| 2 ) + S 2 ε 2 (||φ n+1 || 2 -||φ n || 2 + ||(φ n+1 -φ n )|| 2 ) +2ξ n+1 λ ε h (φ n )T n , φ n+1 -φ n = 0. (3.3)
Using (3.1a), we rewrite (3.1c) as follows

R n+1 -R n τ = 1 2 E 1 (φ n ) g(φ n ), φ n+1 -φ n τ - λ ε h (φ n )M (φ n ), µ n T n+1 + λ ε h (φ n )T n , φ n+1 -φ n τ .
Multiplying (3.1c) with 4τ R n+1 , we deduce

2(|R n+1 | 2 -|R n | 2 + |R n+1 -R n | 2 ) -2ξ n+1 (g(φ n ), φ n+1 -φ n ) -τ λ ε h (φ n )M (φ n ), µ n T n+1 + λ ε h (φ n )T n , φ n+1 -φ n = 0. (3.4)
Furthermore, by taking the inner product of (3.1d) with 2τ λ εK T n+1 , we obtain

λ εK (||T n+1 || 2 -||T n || 2 + ||T n+1 -T n || 2 ) + 2τ λD εK ||∇T n+1 || 2 - 2τ λξ n+1 ε h (φ n )M (φ n )µ n , T n+1 = 0. ( 3

Implementation technique.

It is clear that the efficiency of the scheme depends on whether it can be implemented in an efficient way. Besides the provable unconditional stability, we will show in this subsection that the proposed scheme can be equivalently reformulated into a set of linear elliptic equations, which can be easily solved.

Noticing that ξ n+1 is only a scalar variable, we decompose the solution {φ n+1 , µ n+1 , T n+1 } into the linear combinations as follows:

       φ n+1 = φ n+1 1 + ξ n+1 φ n+1 2 , µ n+1 = µ n+1 1 + ξ n+1 µ n+1 2 , T n+1 = T n+1 1 + ξ n+1 T n+1 2 . (3.6)
We impose for the components φ n+1 i and T n+1 i , i = 1, 2, the same boundary condition as for φ n+1

and T n+1 , respectively. Then the equations (3.1a), (3.1b), and (3.1d) can be rewritten as

       φ n+1 1 -φ n τ = M (φ n ) µ n+1 1 - S 3 ε 2 (φ n+1 1 -φ n ) + S 4 (∆φ n+1 1 -∆φ n ) , µ n+1 1 = S 1 ∆φ n+1 1 - S 2 ε 2 φ n+1 1 . (3.7)        φ n+1 2 τ = M (φ n ) µ n+1 2 - S 3 ε 2 φ n+1 2 + S 4 ∆φ n+1 2 , µ n+1 2 = -g(φ n ) + S 1 ∆φ n+1 2 - S 2 ε 2 φ n+1 2 - λ ε h (φ n )T n . (3.8)        T n+1 1 -T n τ = D∆T n+1 1 , (3.9)
T n+1 2 τ = D∆T n+1 2 + Kh (φ n )M (φ n )µ n . (3.10)
Note that R n+1 = ξ n+1 E 1 (φ n ), it follows from (3.1c):

ξ n+1 A n+1 1 = A n+1 2 , (3.11)
where

A n+1 1 = 2E 1 (φ n ) -(g(φ n ), φ n+1 2 ) + τ λ ε h (φ n )M (φ n ), µ n T n+1 2 -T n µ n+1 2 + τ λ ε h (φ n )M (φ n )T n , S 3 ε 2 φ n+1 2 -S 4 ∆φ n+1 2 , (3.12
)

A n+1 2 = 2 E 1 (φ n )R n + (g(φ n ), φ n+1 1 -φ n ) - τ λ ε h (φ n )M (φ n ), µ n T n+1 1 -T n µ n+1 1 - τ λ ε h (φ n )M (φ n )T n , S 3 ε 2 (φ n+1 1 -φ n ) -S 4 ∆(φ n+1 1 -φ n ) . (3.13)
To simplify the terms above, we use (3.8) to obtain

-(g(φ n ), φ n+1 2 ) - τ λ ε h (φ n )M (φ n ), T n µ n+1 2 + τ λ ε h (φ n )M (φ n )T n , S 3 ε 2 φ n+1 2 -S 4 ∆φ n+1 2 = -(g(φ n ), φ n+1 2 ) - λ ε h (φ n )T n , φ n+1 2 = µ n+1 2 -S 1 ∆φ n+1 2 + S 2 ε 2 φ n+1 2 , φ n+1 2 = φ n+1 2 τ M (φ n ) + S 3 ε 2 φ n+1 2 -S 4 ∆φ n+1 2 -S 1 ∆φ n+1 2 + S 2 ε 2 φ n+1 2 , φ n+1 2 = ε 2 (φ n ) + τ (S 2 + S 3 ) τ ε 2 φ n+1 2 2 + S 1 + S 4 ∇φ n+1 2 2 . (3.14)
Then taking the inner product of (3.10) with τ λ εK T n+1

2

, we have

τ λ ε h (φ n )M (φ n ), µ n T n+1 2 = λ εK T n+1 2 2 + τ D ∇T n+1 2 2 . (3.15)
Summing up (3.14) and (3.15), we get

A n+1 1 = 2E 1 (φ n ) + ε 2 (φ n ) + τ (S 2 + S 3 ) τ ε 2 φ n+1 2 2 + S 1 + S 4 ∇φ n+1 2 2 + λ εK T n+1 2 2 + τ D ∇T n+1 2
2 ) > 0.

(3.16)

Similarly, we can simplify A n+1 2 as follows

A n+1 2 = 2 E 1 (φ n )R n -µ n+1 2 -S 1 ∆φ n+1 2 + S 2 ε 2 φ n+1 2 , φ n+1 1 -φ n + λ εK (-T n+1 2 + τ D∆T n+1 2 , T n+1 1
).

(3.17 iii) Compute ξ n+1 and R n+1 by using (3.16), (3.17), and (3.11).

iv) Compute φ n+1 , µ n+1 and T n+1 by (3.6).

To summarize, the algorithm involves the solution of two variable coefficient linear elliptic equations, two constant coefficient linear elliptic equations, and one algebraic equation. Furthermore, it is a three-layer scheme, compared to the four-layer scheme proposed in [START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF]. In the case of constant mobility, our scheme is reduced to solve four elliptic equations with constant coefficients and one algebraic equation.

4.

A second order scheme 4.1. Construction of the scheme. For ease of notation we will use φn+1 to denote 2ϕ n -ϕ n-1 .

The second order scheme we propose reads:

3φ n+1 -4φ n + φ n-1 2τ = M ( φn+1 ) µ n+1 - S 3 ε 2 (φ n+1 -2φ n + φ n-1 ) + S 4 ∆(φ n+1 -2φ n + φ n-1 ) , (4.1a
)

µ n+1 = -ξ n+1 g( φn+1 ) + S 1 ∆φ n+1 - S 2 ε 2 φ n+1 -ξ n+1 λ ε h ( φn+1 ) T n+1 , (4.1b) 3R n+1 -4R n + R n-1 2τ = 1 2 E 1 ( φn+1 ) g( φn+1 ), 3φ n+1 -4φ n + φ n-1 2τ - λ ε h ( φn+1 )M ( φn+1 ), μn+1 T n+1 -µ n+1 T n+1 - λ ε h ( φn+1 )M ( φn+1 ) T n+1 , S 3 ε 2 (φ n+1 -2φ n + φ n-1 ) -S 4 ∆(φ n+1 -2φ n + φ n-1 ) . (4.1c) 3T n+1 -4T n + T n-1 2τ = D∆T n+1 + ξ n+1 Kh ( φn+1 )M ( φn+1 )μ n+1 , (4.1d) ∂φ n+1 ∂n ∂Ω = 0, ∂T n+1 ∂n ∂Ω = 0, ( 4.1e)
where

ξ n+1 = R n+1 E 1 ( φn+1 ) . (4.2)
Obviously, to start the calculation the scheme (4.1) must be accompanied by a suitable one-step scheme to compute φ 1 , T 1 , R 1 . This can be done, for example, by employing the first step of the scheme (3.1).

Intuitively this is a second-order scheme since all involved terms are approximated with secondorder precision. Although rigorous proof is not available for the time being, the convergence order of the scheme will be confirmed through a series of numerical tests.

Stability analysis.

Below we prove the stability of the scheme (4.1). The stability analysis will make use of the following well-known identities:

2ϕ n+1 (3ϕ n+1 -4ϕ n + ϕ n-1 ) = |ϕ n+1 | 2 -|ϕ n | 2 + |2ϕ n+1 -ϕ n | 2 -|2ϕ n -ϕ n-1 | 2 +|ϕ n+1 -2ϕ n + ϕ n-1 | 2 , (4.3) (ϕ n+1 -2ϕ n + ϕ n-1 )(3ϕ n+1 -4ϕ n + ϕ n-1 ) = |ϕ n+1 -ϕ n | 2 -|ϕ n -ϕ n-1 | 2 +2|ϕ n+1 -2ϕ n + ϕ n-1 | 2 . (4.4)
Theorem 4.1. Let {φ n , T n , R n } be the solution of the discrete problem (4.1). Then for n ≥ 1 it satisfies the discrete energy law:

E n+1 -E n = -Q n+1 -τ ( φn+1 ) 3φ n+1 -4φ n + φ n-1 2τ 2 + λD εK ||∇T n+1 || 2 , (4.5)
where E n is defined by

E n = 1 4 S 1 (||∇φ n || 2 + ||2∇φ n -∇φ n-1 || 2 ) + S 2 ε 2 (||φ n || 2 + ||2φ n -φ n-1 || 2 ) + 2S 3 ε 2 ||φ n -φ n-1 || 2 +2S 4 ||∇(φ n -φ n-1 )|| 2 + λ εK (||T n || 2 + ||2T n -T n-1 || 2 ) + 2(|R n | 2 + |2R n -R n-1 | 2 ) , (4.6)
and Q n+1 is defined by

Q n+1 = S 1 + 2S 4 4 ||∇φ n+1 -2∇φ n + ∇φ n-1 || 2 + S 2 + 2S 3 4ε 2 ||φ n+1 -2φ n + φ n-1 || 2 + λ 4εK ||T n+1 -2T n + T n-1 || 2 + 1 2 |R n+1 -2R n + R n-1 | 2 .
Proof First, we take the inner product of (4.1a) with 2 M ( φn+1 ) (3φ n+1 -4φ n + φ n-1 ), and (4.1b) with 2(3φ n+1 -4φ n + φ n-1 ). Then we sum up the resulting equations and use the identities (4.3) and (4.4) to obtain

4τ 3φ n+1 -4φ n + φ n-1 2τ M ( φn+1 ) 2 + 2S 3 ε 2 ( φ n+1 -φ n 2 -φ n -φ n-1 2 + 2 φ n+1 -2φ n + φ n-1 2 ) +2S 4 ∇(φ n+1 -φ n ) 2 -∇(φ n -φ n-1 ) 2 + 2 ∇(φ n+1 -2φ n + φ n-1 ) 2 +S 1 (||∇φ n+1 || 2 + ||∇(2φ n+1 -φ n )|| 2 -||∇φ n || 2 -||∇(2φ n -φ n-1 )|| 2 + ||∇(φ n+1 -2φ n + φ n-1 )|| 2 ) + S 2 ε 2 (||φ n+1 || 2 + ||2φ n+1 -φ n || 2 -||φ n || 2 -||2φ n -φ n-1 || 2 + ||φ n+1 -2φ n + φ n-1 || 2 ) +2ξ n+1 g( φn+1 ), 3φ n+1 -4φ n + φ n-1 + 4τ λξ n+1 ε h ( φn+1 ) T n+1 , 3φ n+1 -4φ n + φ n-1 2τ = 0. (4.7)
By virtue of (4.1c) and (4.1a), we can rewrite (4.1c) as

3R n+1 -4R n + R n-1 2τ - 1 2 E 1 ( φn+1 ) g( φn+1 ), 3φ n+1 -4φ n + φ n-1 2τ - λ ε h ( φn+1 )M ( φn+1 ), μn+1 T n+1 + λ ε h ( φn+1 ) T n+1 , 3φ n+1 -4φ n + φ n-1 2τ = 0.
Multiplying the both sides by 8τ R n+1 , we have

2(|R n+1 | 2 + |2R n+1 -R n | 2 -|R n | 2 -|2R n -R n-1 | 2 + |R n+1 -2R n + R n-1 | 2 ) -4τ ξ n+1 g( φn+1 ), 3φ n+1 -4φ n + φ n-1 2τ - λ ε h ( φn+1 )M ( φn+1 ), μn+1 T n+1 + λ ε h ( φn+1 ) T n+1 , 3φ n+1 -4φ n + φ n-1 2τ = 0. (4.8)
Furthermore, it follows from taking the inner product of (4.1d) with 4λτ εK T n+1 , and using (4.3):

λ εK (||T n+1 || 2 + ||2T n+1 -T n || 2 -||T n || 2 -||2T n -T n-1 || 2 + ||T n+1 -2T n + T n-1 || 2 ) +4τ D||∇T n+1 || 2 -4τ ξ n+1 (Kh ( φn+1 )M ( φn+1 )μ n+1 , T n+1 ) = 0. (4.9)
Finally we sum up (4.7), (4.8), and (4.9) to conclude. This completes the proof.

The second order scheme (4.1) can also be implemented in an efficient way. We are not going to describe the implementation details of this scheme since it is very similar to the first order scheme (3.1) explained in the previous section. Nevertheless we give here the algebraic equation to compute the auxiliary variable ξ n+1 :

ξ n+1 = A n+1 2 A n+1 1 ,
where

A n+1 1 = 3E 1 ( φn+1 ) + 3 2 3ε 2 ( φn+1 ) + 2τ (S 2 + S 3 ) 2τ ε 2 φ n+1 2 2 + 3 2 S 4 + S 1 ∇φ n+1 2 2 + λ εK 3 2 T n+1 2 2 + τ D ∇T n+1 2 2 > 0, A n+1 2 = E 1 ( φn+1 )(4R n -R n-1 ) - 1 2 µ n+1 2 -S 1 ∆φ n+1 2 + S 2 ε 2 φ n+1 2 , 3φ n+1 1 -4φ n + φ n-1 + λ εK - 3 2 T n+1 2 + τ D∆T n+1 2 , T n+1 1 .

Numerical experiments

In order to illustrate the performance of the proposed numerical method and confirm our analysis results, several numerical tests are carried out and presented in this section. For the sake of convenience, we consider the numerical examples with the square domain Ω = (-1, 1) 2 and fourfold anisotropy, i.e., m = 4 in (2.4). We consider constant mobility (φ) in the calculation. The spatial discretization is a Legendre-Galerkin spectral method. The approximation space for the phase function φ and the temperature T is IP N (Ω), where IP N (Ω) denotes the space of polynomials of degree ≤ N at each space direction. Since the focus of our numerical tests is the verification of the stability and convergence order of the time-stepping schemes, we will take N = 128, which is large enough so that the spatial discretization errors are negligible compared with the the temporal one.

Example 5.1. (Accuracy test) As the first example, we test the convergence order of the schemes (3.1) and (4.1). In the case one, we fabricate two forcing functions in (2.1) and (2.2) so that the exact solution to (2.1)-(2.6) is

Case-I φ(x, y, t) = sin(t) cos(πx) cos(πy), T (x, y, t) = sin(t) cos 2 (πx) cos 2 (πy).

(5.1)

We set the following parameters:

= 4e3, ε = 0.1, σ = 0.05, λ = 0.1, D = 2.25e -2, K = 0.01, S 1 = 0.9, S 2 = 10, S 3 = S 4 = 0, B = 1e4.
(5.2)

In the case two, we choose the initial conditions:

Case-II      φ(x, y, 0) = tanh r 0 -((x -x 0 ) 2 + (y -y 0 ) 2 ) ε 0 , T (x, y, 0) = -0.5φ(x, y , 0), (5.3)
where r 0 = 0.25, x 0 = y 0 = 0, ε 0 = 0.1. The model parameters are set as follows:

= 1e3, ε = 0.1, σ = 0.05, λ = 1, D = 5e -2, K = 0.1, S 1 = 0.9, S 2 = 10, S 3 = S 4 = 0, B = 5e3.
(

In the latter case the exact solution is unavailable, we will use the numerical solution computed with τ = 3e -5 to serve as the exact solution. In Figure 5.1 we plot the L 2 errors in log-log scale of the computed phase and temperature solutions at t = 1 as functions of the time step size τ . As expected, in both cases the obtained numerical convergence rates are in a perfect agreement with the claimed orders; i.e., first order for the scheme (3.1) and second order for the scheme (4.1). It is worth to mention that the calculation with some positive parameters S 3 and S 4 has produced similar results (not shown here). We know from (2.9) that E n defined in (4.6) can be regarded as a discrete version of the original energy functional E(φ n , T n ) defined in (2. However, the dissipation of E(φ n , T n ) loses monotonicity during some time period(τ ≥ 10) for the time step sizes bigger than 10. This is most likely caused by imprecise calculation with large time step sizes. This guess is supported by the observation from Figure 5.2(c), in which the time evolution of ξ n is given. Remember that ξ n , defined in (4.2), should be close to 1 if the approximation is good enough. The results presented in Figure 5.2(c) demonstrate good accuracy on ξ n when τ is not very large, say less than 10. However when τ = 10, ξ n becomes oscillatory and error becomes visible. When τ increases to 100, the computed auxiliary variable ξ n is not any more close to 1. This implies that the numerical solution is not accurate enough for τ ≥ 10, leading to a violation of the monotonic dissipation of the original energy. Thus in practice, it is not recommended to use time step sizes too large, although the calculation can always be stable.

Next test concerns impact of the parameters S 1 , S Example 5.3. (Fourfold anisotropy crystal growth) In this example, we carry out a simulation of crystal growth with fourfold anisotropy, and investigate how the anisotropic coefficient and the latent heat coefficient K affect the shape of the dendritic crystal.

We consider a benchmark problem, which has been extensively studied; see, e.g., [START_REF]Phase-field model of dendritic sidebranching with thermal noise[END_REF][START_REF] Kobayashi | Modeling and numerical simulations of dendritic crystal growth[END_REF][START_REF]A novel fully decoupled scheme with second-order time accuracy and unconditional energy stability for the navier-stokes equations coupled with mass-conserved allen-cahn phase-field model of two-phase incompressible flow[END_REF][START_REF] Zhang | A novel decoupled and stable scheme for an anisotropic phase-field dendritic crystal growth model[END_REF].

We set the initial conditions as follows: starting with the same small circle. Moreover, it is seen that the parameter K affects the width of the branch: larger is K, thinner is the width of the branches, and sharper are the tips.

φ(x, y, 0) = tanh r 0 -((x -x 0 ) 2 + (y -y 0 ) 2 ) ε 0 , T (x, y, 0) = 0, φ(x, y, 0) > 0, - 0
The isocontours of the temperature T at the last moment of each simulation is plotted in the figures 5.4(e). It is observed that the contours of the temperature T take similar dendrite crystal shape as the phase field. This is due to the fact that the heat is propagating only at the interface.

The dissipation behavior of the modified energy E n in time is shown in Figure 5.5(a). The monotonic decay feature of E n for all tested K reflects good stability property of the scheme used in the calculation. Finally we give in Figure 5.5(b) evolution of the area of the crystal, defined by the quantity Ω 1+φ 2 dx, for several values of K. We see that the area of the crystal keeps increasing during the simulation. This is in a good agreement with the existing results; see, e.g., [START_REF]Phase-field model of dendritic sidebranching with thermal noise[END_REF][START_REF]Efficient linear, stabilized, second-order time marching schemes for an anisotropic phase field dendritic crystal growth model[END_REF][START_REF]On a novel full decoupling, linear, second-order accurate, and unconditionally energy stable numerical scheme for the anisotropic phase-field dendritic crystal growth model[END_REF].

Concluding remarks

We have proposed a class of new time-stepping schemes for the anisotropic phase-field dendritic crystal growth model. The proposed schemes were constructed based on an auxiliary variable approach for the Allen-Cahn equation and sophisticated treatment of the terms coupling the Allen-Cahn equation and temperature equation. In particular, the new reformulation of the model introduced in the paper plays a key role in developing efficient schemes. Thanks to the carefully chosen extra terms added to the time discretization, we were able to construct a second-order scheme, which is linear, decoupled, uniquely solvable, and unconditionally stable. A detailed comparison with existing schemes is given, and the advantage of the new schemes are emphasized. The stability property of the proposed schemes was rigorously established, while the convergence rate was carefully examined through a series of numerical tests. Our analysis and numerical experiments demonstrated the efficiency of the proposed method. It seems to us that the approach proposed in this paper is extendable to more complex models such as those studied in [START_REF]A new efficient fully-decoupled and second-order time-accurate scheme for cahn-hilliard phase-field model of three-phase incompressible flow[END_REF][START_REF]A novel fully decoupled scheme with second-order time accuracy and unconditional energy stability for the navier-stokes equations coupled with mass-conserved allen-cahn phase-field model of two-phase incompressible flow[END_REF][START_REF]Numerical approximations of the navier-stokes equation coupled with volume-conserved multi-phase-field vesicles system: Fully-decoupled, linear, unconditionally energy stable and second-order time-accurate numerical scheme[END_REF].

. 5 )Remark 3 . 2 .

 532 Finally, the desired result (3.2) follows from summing up (3.3),(3.4), and (3.5). It is seen from the proof of Theorem 3.1 that the S 3 and S 4 -terms introduced in (3.1a) plays no role in stabilizing the scheme. In fact, it follows from (3.2) that the discrete energy remains dissipative even if S 3 = S 4 = 0.

1 . 7 ) for φ n+1 1 , µ n+1 1 .Solve the equation ( 3 . 8 ) for φ n+1 2 , µ n+1 2 . 1 .

 171138221 )Therefore, ξ n+1 is uniquely determined by dividing the both sides of (3.11) by A n+1 Based on the above discussion, we arrive at the decoupled algorithm for solving the equation set (3.1a)-(3.1e) as follows: Given {φ n , T n , R n , µ n }, we update {φ n+1 , T n+1 , R n+1 , µ n+1 } through: i) Solve the equation(3.ii) Solve the equation (3.9) for T n+1 Solve the equation (3.10) for T n+1 2 .

Example 5 . 2 .

 52 (Stability test) To investigate the stability property of the proposed schemes, we consider the problem Case-II, which has the initial conditions given in (5.3). The parameters used in this test are the same as in(5.4).

  3). According to (2.7) and Theorem 4.1, both E(φ n , T n ) and E n should be monotonically decreasing with the time step n ≥ 1. Notice that E n is usually called as modified energy functional, which is not necessarily an approximation to the original energy functional E(φ n , T n ).The modified energy functional E n , original energy functional E(φ n , T n ), and ξ n computed by the scheme (4.1) with different time step sizes are presented in Figure5.2 as functions of time. It is observed in Figure5.2(a) that the modified energy functional E n is indeed strictly dissipative in time as predicted by Theorem 4.1, even for very large time step sizes τ = 10, 100. As shown in First order scheme (3.1) for the Case-I.

  First order scheme (3.1) for the Case-II.

  Second order scheme (4.1) for the Case-I.

  Second order scheme (4.1) for the Case-II.

Figure 5 . 1 .Figure 5 . 2 .

 5152 Figure 5.1. (Example 5.1) Convergence order of the time-stepping schemes: L 2 errors of the phase field function and the temperature as functions of the time step size τ .

Figure 5 .

 5 Figure 5.2(b), the original energy E(φ n , T n ) is also dissipative for all relatively small time step sizes.

  (a)-(b) that for S 1 = S 2 = S 3 = S 4 = 0 and S 1 = 0.1, S 2 = 4, S 3 = S 4 = 0, the computed ξ n is quite inaccurate, i.e., far from the exact value 1, even with small time step sizes (0.001 and 0.01 respectively). However taking positive S 3 and S 4 , say S 3 = S 4 = 5, allows recovering the accuracy, as shown in Figures 5.3(c)-(d). It is notable that the presence of the S 3 -and S 4 -terms allows stable and accurate calculation even with τ = 10. This test clearly indicates the benefit of the stabilization terms. It is also worth to mention that all the cases above produced monotonically decreasing energy E n , which is consistent with what we have proved in Theorem 4.1.

  .6, otherwise, where r 0 = 9e -4, x 0 = y 0 = 0, ε 0 = 1.8e -4. The model parameters used in the simulation are: = 1e3, ε = 0.015, σ = 0.1, λ = 4e2, D = 2.5e -3, S 1 = 0.6, S 2 = 10, S 3 = S 4 = 4, B = 4e5, τ = 0.01.We let the latent heat parameter K vary. The spatial spectral discretization uses the polynomial space of degree 512 at each direction.InFigures 5.4

5 Figure 5 . 3 .

 553 Figure 5.3. (Example 5.2) History of ξ n with different chioces of the parameters S 1 , S 2 , S 3 , S 4 .

  (a) φ at t = 0, 3, 6, 9. K = 0.6 (b) φ at t = 3, 6, 9, 11. K = 0.8 (c) φ at t = 6, 9, 11, 14. K = 1 (d) φ at t = 9, 11, 14, 17. K = 1.2 (e) Temperature field T at the last moment of above cases. From left to right: K = 0.6, K = 0.8, K = 1, K = 1.2.

Figure 5 . 4 .

 54 Figure 5.4. (Example 5.3) Dendritic crystal growth with fourfold anisotropy for different values of the latent heat parameter K. (a)-(d): snapshots of the phase field φ at different times; (e): the temperature field T at the last moment of (a)-(d).

Figure 5 . 5 .

 55 Figure 5.5. (Example 5.3) Time evolution of the modified energy E n and crystal area Ω 1+φ 2 dx for different latent heat parameter K.

  1) into (2.10a) and (2.10b). Another notable idea is to replace φ t in (2.2) by using the equation (2.10a), resulting in an equivalent equation, i.e., (2.10d). Although, and obviously, (2.1)-(2.2) and (2.10) is strictly equivalent to each other at the continuous level, we will see in what follows that

the reformulation (2.10) facilitates construction of decoupled, stable, higher order convergent, and cheaper time stepping schemes. Since (2.1)-(2.2) and (2.10) are equivalent, the latter obviously satisfies the same energy dissipative law as (2.7). However, to better understand the discrete energy dissipative law, it is desirable to derive an alternative form of the energy law involving the auxiliary variables for (2.10). To this end, we take the L 2 inner product of (2.10a) with -µ, (2.10b) with φ t , (2.10c) with -2R, and

  2 , S 3 , and S 4 on the quality of numerical solutions. Since the auxiliary variable is indicative of the quality of the computed solutions, we only report the computed values of ξ n . Figure 5.3 shows the computed ξ n versus the time for four parameter sets {S i }. It is observed from Figures 5.3

* This research is partially supported by NSFC grant 11971408, NNW2018-ZT4A06 project, and NSFC/ANR joint program 51661135011/ANR-16-CE40-0026-01.


Maria Bernard

Olivier Rué

Mahendra Mariadassou

Géraldine Pascal

Geraldine Pascal
email: geraldine.pascal@inrae.fr

FROGS: A POWERFUL TOOL TO ANALYSE THE DIVERSITY OF FUNGI WITH SPECIAL MANAGEMENT OF INTERNAL TRANSCRIBED SPACERS FROGS efficiently manages ITS amplicons

She is a member of a platform team conducting NGS sequence analysis and designing software. She specializes in workflow development in particular for metabarcoding analysis.

Olivier Rué is a bioinformatics engineer. He is in charge of data analysis at the Migale bioinformatics facility. He specializes in the analysis of metabarcoding and metagenomics data.

Mahendra Mariadassou has a Ph.D. in statistics. He is involved in the development of new statistical methods and tools for metabarcoding analysis.

Géraldine Pascal has a Ph. D. in bioinformatics and coordinates the FROGS project. She is currently involved in designing solutions for long read problems, workflow development and metagenomics analysis.

INTRODUCTION

Using amplicon sequencing to describe the microbial composition of an environment is a time saving and cost-effective strategy and can be used even for very large-scale surveys [1]. Most studies currently focus on the bacterial fraction of microbial communities but the fungal fraction is equally important, as fungi are ubiquitous and provide several ecosystem services [START_REF] Naranjo-Ortiz | Fungal evolution: diversity, taxonomy and phylogeny of the Fungi[END_REF]. Unfortunately, studying the fungal fraction using metabarcoding has its own challenges. Indeed, in fungi, there is no equivalent of the 16S rRNA gene, which is widely used and highly suitable for bacteria. The best candidates are internal transcribed spacers (ITS), but these are more difficult to manipulate. The main problem with ITS is size polymorphism, with a size range of 361 to 1475 bases in UNITE 7.1 [START_REF] Koljalg | Towards a unified paradigm for sequencebased identification of fungi[END_REF] (unlike 16S where 95% of the sequences have a length between 1205 and 1556 bases). Most studies describing ITS data analyses process either (i) paired-end reads but filter out non-overlapping, non-mergeable reads, thus systematically discarding taxa with longer ITS, or (ii) single-end reads, thus limiting taxonomic resolution and losing the benefit of information contained in longer sequences [4,[START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF]. In both cases, the tools used, which only support mergeable paired-end reads or single-end reads, are unable to process all the sequences produced. The new version of FROGS [START_REF] Escudie | FROGS: Find, Rapidly, OTUs with Galaxy Solution[END_REF] solves these challenges. This metabarcoding analysis pipeline runs not only on the command line but also on the Galaxy platform, meaning it is easy to use for non-experts. FROGS tools make it possible to generate and explore the abundance of OTUs (Operational Taxonomic Units) in the target environment along with their taxonomic affiliation, FROGS also includes 7 tools dedicated to descriptive statistics and two tools for differential abundance analysis (supplemental data section 1). FROGS is widely used by the research community worldwide to study microbial environments i.e., water, skin, soil, gut, animal, human, food microbiomes [START_REF] De Lorgeril | Immune-suppression by OsHV-1 viral infection causes fatal bacteraemia in Pacific oysters[END_REF][START_REF] Teyssier | Inside the guts of the city: Urban-induced alterations of the gut microbiota in a wild passerine[END_REF][START_REF] Frere | Microplastic bacterial communities in the Bay of Brest: Influence of polymer type and size[END_REF][START_REF] Hoyles | Molecular phenomics and metagenomics of hepatic steatosis in non-diabetic obese women[END_REF][START_REF] Duron | Tick-Bacteria Mutualism Depends on B Vitamin Synthesis Pathways[END_REF]. Enriched with new features every year, FROGS is a very powerful tool to treat ITS amplicons. Here, we propose a smart method and the tools needed to deal with both overlapping and nonoverlapping paired-end reads of, in particular, ITS.

AVAILABILITY AND IMPLEMENTATION

FROGS can be used both on the Galaxy platform and on a command line. To process ITS data correctly, we recommend a seven-step process shown in Figure 1 (for command lines, see supplemental data section 2). Here, we only highlight the adaptations made to process ITS data and the novelties introduced since version 3.0 [START_REF] Escudie | FROGS: Find, Rapidly, OTUs with Galaxy Solution[END_REF].

FROGS Preprocessing: After the paired-end amplification process, ITS reads may not overlap.

In this step, merging tools such as PEAR [START_REF] Zhang | PEAR: a fast and accurate Illumina Paired-End reAd mergeR[END_REF], VSEARCH [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF] or Flash [START_REF] Magoc | FLASH: fast length adjustment of short reads to improve genome assemblies[END_REF] cannot merge these reads and they may thus be lost by the system. The solution we found to work best to avoid losing ITS amplicons is to both conventionally process mergeable data and to create artificial sequences with non-mergeable sequences. The idea is to join the two reads together with a stretch of 100 N letters. These "un-merged" sequences are tagged as "combined" sequences.

This option is activated only if reads do not entirely cover the target sequence. Otherwise, noise is retained in analyses. FROGS "combined" sequences are artificial and present particular features, especially their size. For a MiSeq sequencing of 2x250pb, FROGS will combine nonoverlapping reads into a single sequence 600 bp long. So, for ITS analysis to keep both merged and unmerged reads, choose --keep-unmerged parameter. We would normally advise using PEAR [START_REF] Zhang | PEAR: a fast and accurate Illumina Paired-End reAd mergeR[END_REF] to merge sequences because it produces the best results, but because non academics require a commercial license, we chose VSEARCH [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF] as our default paired-end read merger.

FROGS Clustering: To cluster both classical amplicons and unmergeable amplicon as ITS, we advise to choose the parameters of distance equal to 1 and the fastidious option (--distance 1 --fastidious) as recommended in Swarm v2 [START_REF] Mahe | Swarm v2: highly-scalable and high-resolution amplicon clustering[END_REF]. We have changed the guidelines from those given in FROGS v.1., and these new parameters provide an equally well defined cluster but are much faster. Swarm handles "combined" sequences very well. We simply and temporarily change the stretch of 100 N into a stretch of 50 A and 50 C. Then, Swarm treats these combined sequences as mergeable sequences. The principle of Swarm is that it creates links between sequences two by two when they have only one difference between them (with distance parameter equals 1). Thus, Swarm does the same for combined sequences, the 50 A 50 C stretches not being discriminating between combined sequences since they are identical, these differences being elsewhere on the sequence.

FROGS Chimera removal: there is no impact on ITS analysis with this tool, use it without specificity.

FROGS OTU Filter: Users must remove too low abundance clusters that merge with clusters of artifact sequences. The problem is doubtless not a bioinformatics problem but rather the use of metabarcoding to detect rare microbes in a large population of microbes. The best way to remove these small noisy clusters is to apply a 0.005% (--minabundance) abundance threshold and/or, if applicable, to keep the clusters only if they are present in a certain number of samples or replicated samples (--min-samplepresence). The user now has the possibility to provide a FASTA file containing potential contaminant sequences e.g., phiX, host sequences, mitochondrial sequences, OTUs that align with a contaminant sequence, are deleted.

FROGS ITSx: If the objective is to analyse fungal ITS, we advise only keeping sequences detected by ITSx [START_REF] Bengtsson-Palme | Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data[END_REF] to have an ITS signature with no conserved flanking sequences truncation (--check-its-only). ITSx identifies the highly conserved neighbouring rRNA sequences SSU, 5.8S and/or LSU. If the targeted variable region (ITS1 or ITS2) is not detected, the sequence is discarded. The ITSx step is time consuming and has to be done after clustering and filtering steps.

FROGS Affiliation: users have to select one of available ITS databanks. The affiliation process is specific for ITS, adapted to un-merged data thanks to a mix of local and global alignments.

Because of the inclusion of the 100 N stretch in the "combined" sequences, the affiliation of the OTU seed sequences against the databases by BLASTN+ [START_REF] Camacho | BLAST+: architecture and applications[END_REF], as in FROGS v1, was not possible. BLASTN+ aligns the sequences locally and, by definition, only keeps the best local alignment, i.e., either the 5' sequence or the 3' sequence, but in no case the entire "combined" sequence, since the N stretch collapses the BLAST score. We thus chose a double strategy. All normally merged sequences are treated like in FROGS v1 by BLASTN+ against a databank (e.g. UNITE) while the combined sequences, which are automatically detected by the system, are treated with needleall [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF], a global sequence aligner. To align with needleall takes too long, to reduce this time, we chose to align combined sequences not against the entire databank (e.g. UNITE) but against a drastically reduced bank. To this end, first, the two sequence parts and for all: the absence of the full or partial taxon name. Two modes are possible to process these impacted OTUs. In the masking mode (--mask), the taxonomic information is replaced by NA. In this case, affiliations with insufficient characteristics are not retained. In delete mode (--delete), the OTUs are deleted according to the filter settings. Among other things, this makes it possible to delete chimeric OTUs (too low percentage coverage). Another feature enabled by the tool is being able to mask or delete taxonomic information by specifying particular key words (e.g. unknown species). This functionality is particularly useful to remove contaminants that can be detected using taxonomic information.

(

DATASETS AND TESTS TO ASSESS THE MANAGEMENT OF ITS:

To test our methodology, we generated increasingly complex simulated datasets, available on https://doi.org/10.15454/AOT7UL (supplemental data section 3, supplemental Figure 1) and used a mock biological community (supplemental data section 4). These data were processed by FROGS but also by USEARCH, QIIME2 and DADA2. Since QIIME2 and DADA2 filter out nonoverlapping reads, we used QIIME2 and DADA2 in single-end and paired-end modes. For USEARCH, we followed the original guidelines by taking into account merged sequences and 5' R1 reads of non-overlapping paired-end sequences. All benchmarking command lines are available on http://frogs.toulouse.inrae.fr/ITS/frogs-its.html#Command_lines. We used four metrics to assess our results: (i) the divergence rate: computed as the Bray-Curtis distance between expected and observed abundance profiles at a given taxonomic level (observed composition is the abundance resulting from each tool process), (ii) the number of falsenegative taxa (FN): the number of expected taxa that were not recovered by the method, (iii)

the number of false positive taxa (FP): the number of recovered taxa that were not expected, and (iv) the number of true positive taxa (TP): the number of recovered taxa that were expected. We computed the precision (TP/(TP+FP)) and the recall rate (TP/(TP+FN)) of all methods, on all datasets: 35, 115 and 515 species, ITS1, ITS2, power law or uniform distribution of abundances, of simulated and real data. These two metrics allow us to evaluate the performance of each tool based on FP, TP, and FN metrics.

RESULTS AND DISCUSSION

On simulated data, whatever the complexity of the microbial diversity of the samples (35 As expected, species with low abundance rates were the most difficult to recover. The reference reconstruction and the number of OTUs with the expected length was very good with FROGS. The results of comparison with other tools showed, as expected, that tools that process paired-end reads (DADA2-pe, QIIME-pe) missed unmergeable ITS and tools that only process R1 reads (DADA2-se and QIIME-se) lost too much sensitivity and had more difficulty correctly identifying the species in the samples. This calls into question the strategy of processing only R1 reads. Despite good precision, DADA2 combined a low recall rate (supplemental Figure 2), mis-reconstruction of the reference (Figure 2D) and a high divergence rate. USEARCH reconstructed the references quite well and correctly identified the expected species, but at the price of more FP resulting in very limited precision (Figure 2A, B and supplemental Figures 3 to 6). The USEARCH filter applied to the bootstrap value at the end of the process, as recommended, removed some TP ZOTUs (Zero-radius OTUs). The great strength of FROGS is its ability to handle both mergeable and unmergeable reads but still being easy to use. However, it should be noted that the quality of the unmergeable sequences is lower because they do not benefit from the correction enabled by the overlapping parts. We can therefore expect less good identification of fungi with unmergeable ITS amplicons than of mergeable ones. Moreover, the processing unmergeable sequences requires global sequence alignment, which increases FROGS computing time.

(
On real data (http://frogs.toulouse.inrae.fr/ITS/frogs-its-meat.html, supplemental Figure 7 and supplemental Table ), FROGS produced very good results. The divergence rates of all the tool were worse than on simulated data because they were calculated on expected data while it is difficult to match theoretical and real quantities when preparing samples. The greatest divergences could therefore be due to this. Concerning FP, TP and FN, FROGS showed high precision and a satisfactory recall rate (median 0.99 and 0.92 respectively), even if, as expected, FROGS was not able to differentiate between species that are too similar, e.g.

Penicillium. Nevertheless, considering the taxonomies found, FROGS was the tool that detected the most expected species. FROGS also reconstructed the expected sequences best (100% identity and 100% coverage of the entire sequence).

CONCLUSION

FROGS is an easy-to-use suite that is perfectly suited and efficient not only for the processing of internal transcribed spacer amplicon sequences, but also of all other unmergeable amplicons i.e. rpb2 [START_REF] Větrovský | The rpb2 gene represents a viable alternative molecular marker for the analysis of environmental fungal communities[END_REF] and D1-D2 [START_REF] Sonnenberg | An evaluation of LSU rDNA D1-D2 sequences for their use in species identification[END_REF]. The test of FROGS data processing demonstrated excellent recall rates and precision thanks to smart ITS data management while accurately reconstructing the sequences as they are expected. FROGS can be used both on the command line and on the Galaxy interface, making it easy for everyone to access.

KEY POINTS

 FROGS processes ITS amplicons (and any other non-overlapping reads) with a smart method.

 FROGS was compared to many popular metabarcoding tools and proved to be very sensitive with confident OTU sequence reconstruction.

 FROGS produces very few false positives and false negatives, has an excellent recall rate and high precision.

 FROGS can be run both on the command line and on the Galaxy interface.

KEYWORDS

Fungi, ITS, metabarcoding, workflow, amplicon, metagenomics FIGURES: SUPPLEMENTAL DATA

NINE FROGS STATISTICAL TOOLS

The analysis of OTU (Operational Taxonomic Unit) tables with respect to the associated metadata (pH, temperature, weight, treatment, etc.) requires the use of statistical methods dedicated to compositional data. The OTU abundance table produced FROGS contains the abundance information for each OTU in each sample as well as its taxonomic affiliation. Questions such as "How many species are there in each sample i.e., species richness?" and "Which OTUs are abundant in my samples, i.e. structure" or "How does the microbial community vary with the different metadata?" can only be based on the calculation of distances, indices and descriptive statistics. FROGS therefore has seven tools for the statistical analysis of samples using the phyloseq R package, extended R functions in Rmarkdown scripts. The first tool imports data and transforms data in the abundance table and metadata into an Rdata object. The second tool uses phyloseq and a custom R function to build two plots to visualize the composition of the sample. The third tool builds richness plots to visualize the alpha diversity of the sample. The fourth computes the beta distance matrixes and displays them on a heatmap. The fifth builds a heatmap plot and an ordination plot using ordination methods to visualize data structure. The sixth builds the clustering plot for cluster analysis with different linkage methods, and the seventh uses phyloseq and vegan functions in R to perform multivariate analysis of variance (MANOVA).

As a preliminary step, users can obtain a phylogenetic tree based on OTUs using the FROGS Tree tool. This tool creates a multiple alignment of OTUs with Mafft and reconstructs phylogenetic tree with FastTree and phangorn R package.

Finally, FROGS include also two statistical tools dedicated to differential abundance analysis. Like previous tools, is uses Rmarkdown scripts, but based on DESeq2 R package. The first tool produces a DESeq2 object including the statistical model, and the second tool exports the comparison results between chosen conditions by producing differentially abundant OTU table, pie chart, MA plot and volcano plot.

As for bioinformatics tools, almost all tools return HTML reports with a variety of graphics and tables. # Abundance filter process otu-filters.py --min-abundance 0.00005 --nb-cpus $CPU --input-fasta remove_chimera.fasta --input-biom remove_chimera.biom --output-biom filters.biom --output-fasta filters.fasta --log-file filters.log --summary filters.html # ITSX filter itsx.py --nb-cpus $CPU --region $its --check-its-only --input-fasta filters.fasta --input-biom filters.biom --out-fasta itsx.fasta --out-abundance itsx.biom --summary itsx.html --log-file itsx.log # Affiliation versus Unite affiliation_OTU.py --nb-cpus $CPU --reference $REFERENCE --input-biom itsx.biom -input-fasta itsx.fasta --output-biom frogs.biom --summary affiliation.html --logfile affiliation.log # Tabular abundance table biom_to_tsv.py --input-biom frogs.biom --input-fasta itsx.fasta --output-tsv frogs.tsv --output-multi-affi multiaff.tsv --log-file biom_to_tsv.log # Affiliation filter: parameters is very dependent on your dataset and questions ./affiliation_filters.py --input-biom frogs.biom --input-fasta itsx.fasta --output-frogs_affi_filter.biom --output-fasta frogs_affi_filter.fasta --summary frogs_affi_filter.html --impacted frogs_affi_filter_impacted_OTU.tsv --impactedmultihit frogs_affi_filter_impacted_OTU_multihit.tsv --log-file affi_filter.log \ --delete / --mask # choose behavior # add blast affiliation filter \ #--min-blast-identity 0.8 --min-blast-coverage 0.8 --ignore-blast-taxa "known_contaminant_terme or key_species"

FROGS COMMAND LINES TO PROCESS ITS READS

TEST DATASET FOR ASSESSMENT OF THE FROGS ITS METHODOLOGY

Our objective was to create datasets which make it possible, among other things, to test methods and tools to process data on fungi for metabarcoding. Consequently, we obtained computer-generated paired reads 2 x 250 nucleotides in length from artificial PCRs produced with Grinder software [1] on ITS1 or ITS2 from fungi which are as diverse as possible. Datasets are available on https://doi.org/10.15454/AOT7UL. We extracted the reference sequences from the UNITE v7.1 database centred on the eukaryotic nuclear ribosomal ITS region, we selected species without ambiguous key words i.e. without unknown and incertae key words in their name. UNITE 7.1 from general FASTA releases contains 54 568 sequences composed of 8 179 and 46 389 sequences from respectively reference and representative species. It is composed of 54 291 Fungi with 10 phyla, 52 classes and 185 orders. We produced simulated data reflecting this distribution. The sequences used to support artificial amplification by in silico PCR were taken from the UNITE v7.1 database thanks to the python script treeSampling.py (https://github.com/geraldinepascal/FROGS/tree/master/assessment/bin) and see supplemental Figure 1 that describe the following procedure. We picked 35, 115, or 515 species of fungi to create three different types of diversity. We used Grinder, which is, among other things, a programme designed to create random amplicon sequence libraries based on DNA reference sequences. We produced paired-end reads R1 and R2, type Illumina Miseq. We incorporated sequencing errors. Reads were simulated by increasing the number of substitutions and indels linearly along the reads (Grinder parameters: -mutation_dist linear 2.54e-1 2.79e-1 and -mutation_ratio 98.6 1.4). We generated samples with two types of abundance distribution that follow either a power law or a uniform law. For each type of diversity, we merged these first datasets with exclusive datasets of chimera (3 types of chimera -Grinder parameters: -chimera_dist 314 38 1) to obtain an overall chimera rate of 2% per dataset. We focused on including reads from ITS1 and ITS2 of long lengths, not overlaid by two paired 250-nucleotide reads (min:132 nucl. to max 763 nucl. c.f. sections #Amplicon length distribution and #List of non-mergeable amplicons available from http://frogs.toulouse.inrae.fr/ITS/frogs-its.html#35_species, #115_species, #515_species). Thus, the resulting paired-end reads become unmergeable using sequence merging tools such as FLASH, VSEARCH or PEAR and reflect the reality of the sequences obtained in vitro.

Figure 1 :Figure 2 :

 12 Figure 1: Standard operation procedure for long and unmergeable reads from metabarcoding sequencing i.e. ITS, rpb2, D1-D2

  illumina --input-archive $ARCHIVE_FILE --keep-unmerged --mergesoftware pear --min-amplicon-size 20 --max-amplicon-size 490 --without-primers --R1-size 250 --R2-size 250 --nb-cpus $CPU --output-dereplicated preprocess.fasta -output-count preprocess.tsv --summary preprocess.html --log-file preprocess.log # Clustering clustering.py --input-fasta preprocess.fasta --input-count preprocess.tsv -fastidious --distance 1 --nb-cpus $CPU --output-biom clustering.biom --output-fasta clustering.fasta --log-file clustering.log # Remove Chimera remove_chimera.py --input-fasta clustering.fasta --input-biom clustering.biom --nbcpus $CPU --out-abundance remove_chimera.biom --non-chimera remove_chimera.fasta --log-file remove_chimera.log --summary remove_chimera.html

FalseSupplemental figure 5 :

 5 Negative and True Positive rates) are calculated in relation to the expected. Each tool was compared to FROGS using a Wilcoxon signed-rank test, ****: p <= 0.0001; ***: p <= 0.001; **: p <= 0.01; *: p <= 0.05; ns: p > 0.05). Affiliations and associated abundances are taken into account. (N.B. all results on ITS2 datasets are comparable). See complete results on: http://frogs.toulouse.inrae.fr/ITS/frogsits.html#Metrics_calculated_in_relation_to_the_expected19 Boxplots of the results of the processing of the synthetic datasets (515 species, ITS1, uniform abundance distribution) of the 4 tools tested (DADA2 paired-end (pe), DADA2 single-end (se) (pale and dark purple), FROGS (bluish green), QIIME-pe, QIIME-se (pale and dark khaki) and USEARCH (cyan)). Four metrics (Abundance divergence, False Positive, False Negative and True Positive rates) are calculated in relation to the expected. Each tool was compared to FROGS using a Wilcoxon signed-rank test, ****: p <= 0.0001; ***: p <= 0.001; **: p <= 0.01; *: p <= 0.05; ns: p > 0.05). Affiliations and associated abundances are taken into account. (N.B. all results on ITS2 datasets are comparable). See complete results on: http://frogs.toulouse.inrae.fr/ITS/frogsits.html#Metrics_calculated_in_relation_to_the_expected28

  /github.com/geraldinepascal/FROGS-wrappers. FROGS supports ITS since version 3.0 and is fully documented on the companion website http://frogs.toulouse.inrae.fr/.

	https:/FROGS: 7-STEP GUIDELINES FOR THE MANAGEMENT OF INTERNAL
	TRANSCRIBED SPACERS (ITS)
				Installing and
	updating FROGS are easy thanks to the FROGS repositories in CONDA
	(https://anaconda.org/bioconda/frogs)	and	Galaxy	Toolshed
	(https://toolshed.g2.bx.psu.edu/view/frogs/frogs). Two github FROGS repositories are
	effective	https://github.com/geraldinepascal/FROGS	and

  FROGS accurately recovered the species originally present in the samples. FROGS missed very few fungal sequences, had very few or no FP and very few or no FN. Precision (min. 0.991 max. 1) and the recall rates (min. 0.976 max. 1) were excellent (supplemental Figure2and supplementary Table), and better than those of the other tools.

	Figure	2),	115	or	515	species	(see	companion	website
	http://frogs.toulouse.inrae.fr/ITS/frogs-its.html#35_species, #115_species #515_species and
	supplemental Figures 2 to 6)),

ACKNOWLEDGMENTS

The authors are grateful to the genotoul bioinformatics platform Occitanie (doi: 10.15454/1.5572369328961167E12) and Migale bioinformatics facility (doi: 10.15454/1.5572390655343293E12) for providing help, computing and storage resources. The authors thank Valentin Marcon, Patrick Durand, Laure Quintric and Olivier Inizan for the packaging of FROGS under CONDA. The authors thank Ta Thi Ngan for her involvement in the development of FROGSSTAT tools. The authors also thank the members of the METABARFOOD project for allowing us to use the data generated as part of the project, particularly Monika Coton, Jérôme Mounier and Audrey Pawtowski who were in charge of the fermented meat samples. The authors are grateful to Daphne Goodfellow for attention to the English-language version.

Supplemental figure 7: boxplot showing the overall precision (blue) and recall rate (orange) of each of the tools used to process biological datasets (ITS1 and ITS2 datasets with ADN or PCR products).

NOTES ON SUPPLEMENTAL TABLE

Supplemental table is joined to the document. This document contains two tables of data from the "simulated data" and the "biological data" used to test the tools and which made it possible to produce the supplemental figure 2 and7. Respectively, the captions of these tables are:

"Companion table of the supplemental figure 2, indicating the performance (divergence rate, true positive, false negative counts and precision and recall rates) of the different methods (DADA2-pe, DADA2-se, FROGS, QIIME-pe, QIIME-se, USEARCH) used to process test datasets (all datasets: 35, 115 and 515 species, ITS1 and ITS2, power law and uniform abundance distribution). " And "Companion table of the supplemental figure 7, indicating the performance (divergence rate, true positive, false negative counts and precision and recall rates) of the different methods (DADA2-pe, DADA2-se, FROGS, QIIME-pe, QIIME-se, USEARCH) used to process meat datasets (ITS1 and ITS2 with ADN or PCR products)."