Assyr Abdulle
email: assyr.abdulle@epfl.ch

Gilles Vilmart
email: gilles.vilmart@bretagne.ens-cachan.fr

PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise

Keywords: ROCK method, stabilized second-order integration method, partitioned Runge-Kutta methods, stiff problems, advection-diffusion-reaction problems, stochastic problems AMS subject classification (2010): 65L20, 65M12, 65M20

A partitioned implicit-explicit orthogonal Runge-Kutta method (PIROCK) is proposed for the time integration of diffusion-advection-reaction problems with possibly severely stiff reaction terms and stiff stochastic terms. The diffusion terms are solved by the explicit second order orthogonal Chebyshev method (ROCK2), while the stiff reaction terms (solved implicitly) and the advection and noise terms (solved explicitly) are integrated in the algorithm as finishing procedures. It is shown that the various coupling (between diffusion, reaction, advection and noise) can be stabilized in the PIROCK method. The method, implemented in a single black-box code that is fully adaptive, provides error estimators for the various terms present in the problem, and requires from the user solely the right-hand side of the differential equation. Numerical experiments and comparisons with existing Chebyshev methods, IMEX methods and partitioned methods show the efficiency and flexibility of our new algorithm.

Introduction

We consider systems of ordinary differential equations (ODEs) representing space discretizations of PDEs of the form 1 ẏ = F (y) = F D (y) + F A (y) + F R (y), y(0) = y 0 .

(

) 1
where F D (y), F A (y), F R (y) ∈ R n represent diffusion terms with eigenvalues close to the negative real axis, advection terms with eigenvalues close to the imaginary axis, and stiff reaction terms, respectively. Assuming that (1) comes from a spatial discretization of an advectiondiffusion-reaction problem in spatial dimension N and considering a spatial grid size ∆x that leads to a system of ODEs of size n = O((∆x) -N ), then the eigenvalues of the Jacobian of F D are typically distributed along the negative real axis in an interval that grows as [-O((∆x) -2 ), 0] (for a symmetric diffusion operator), the eigenvalues of F A are typically distributed on the imaginary axis in an interval that grows as [-iO((∆x) -1 ), iO((∆x) -1 )], i = √ -1, while the eigenvalues of the Jacobian of the reaction term F R are usually not related to ∆x but have a ratio (sometime called the stiffness ratio) max j |ℜλ j |/ min j |ℜλ j | that can be very large or vary over several orders of magnitude. The question of stability is central when applying an ODE solver to [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF] and it is well known that explicit solvers usually face step size restriction for such problems (usually called stiff), essentially due to the terms F D and F R [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF]. Implicit methods have usually no such restriction but come at the cost of solving a large nonlinear system at each time step. For such nonlinear systems, iterative methods that are preferably used can become quite involved, particularly when the system has a complicated nonlinear structure [START_REF] Hundsdorfer | Numerical solution of time-dependent advection-diffusion-reaction equations[END_REF].

Alternatives to fully implicit solvers have been proposed in the past. First, for problems involving only F D , explicit stabilized Runge-Kutta (RK) methods have been developed [START_REF] Lebedev | Explicit difference schemes with time-variable steps for solving stiff systems of equations[END_REF][START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF][START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF][START_REF] Abdulle | Fourth order Chebyshev methods with recurrence relation[END_REF][START_REF] Verwer | RKC time-stepping for advection-diffusionreaction problems[END_REF]. These methods considerably reduce the cost of standard explicit methods by allowing for much larger time steps, thanks to a special stabilization procedure (involving shifted Chebyshev-like polynomials), obtained by increasing the internal stages s of the RK method. As the stability domains grow quadratically with s (also a measure of the numerical work in terms of functions evaluation) along the negative real axis, this procedure is much more efficient than taking s steps of a standard explicit method (for which the stability domain would only grow linearly with the number of functions evaluation). These methods are much easier to implement and to use than implicit solvers and can usually just be substituted into an explicit Euler procedure. They have also been shown to be competitive with implicit solvers for some classes of problems. This is illustrated in [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF][START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF][START_REF] Abdulle | Fourth order Chebyshev methods with recurrence relation[END_REF] with reaction-diffusion problems and in [START_REF] Hundsdorfer | Numerical solution of time-dependent advection-diffusion-reaction equations[END_REF] for a combustion problems. In general, when the spatial discretization leads to problem that prevents the use of direct solver for implicit methods (e.g., (sparse) LU), the efficiency of implicit solvers will be problem dependent and rely on the efficiency of appropriate iterative methods, good pre conditioners. Furthermore, efficient parallelization of the linear algebra arising from implicit solvers is a non trivial task, while it is straightforward for explicit solvers. Stabilized methods can also accommodate an advection term F A , but remain efficient provided that the diffusion is dominant (small Peclet number regime). Moderate reaction terms F R can also be integrated with such methods.

However, when the reaction term becomes very stiff or in the large Peclet number regime, classical explicit stabilized methods are no longer attractive. They can still be used to integrate part of the equation (1) in a splitting strategy. In presence of very stiff reactions, integrating the diffusion term with a stabilized method and the reaction term with an implicit method can be efficient [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF]. Indeed, as the reaction term alone has no spatial connectivity, the nonlinear system for the reaction term can be decomposed in small independent systems that can be solved efficiently (e.g. with an LU decomposition). However, time step control and error estimates can be delicate with splitting methods. Due to the splitting error, the method proposed in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF] has second order even though each individual integrator used in the splitting strategy has a higher order. Another strategy based on partitioned RK methods (using the stabilized RKC method for the diffusion terms) has become popular. We mention the development of the implicit-explicit solver IRKC [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF], a modification [START_REF] Sommeijer | On stabilized integration for time-dependent PDEs[END_REF] as a two-step method with improved imaginary stability, and the partitioned method PRKC [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF]. We notice that the PRKC method [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF] is of second order, while the IRKC method [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF] is formally of first order (due to the implicit Euler solver used for the reaction term), but with a global error constant that scales as O(h 2 + h/s 2 ) where s is the number of internal stages. The IRKC method is intended for problems of the type [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF] with possibly stiff reactions treated implicitly, while the PRKC method is intended for problems of the type (1) with a non-stiff reaction term F R which is potentially expensive to evaluate (the F D and F R terms are then integrated by different explicit methods). Finally, generalizations of stabilized methods for problems [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF] that also include white noise have recently been proposed [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations[END_REF].

In this paper we introduce a new partitioned implicit-explicit integrator, called PIROCK, based on the explicit second order orthogonal Runge-Kutta Chebyshev method (ROCK2) introduced in [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF] and combining ideas from [START_REF] Verwer | RKC time-stepping for advection-diffusionreaction problems[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations[END_REF][START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF]. We derive a single algorithm that can combine the diffusion term F D with any combination of the term F A , F R and that also treats Itô stochastic systems of the form

ẏ = F (y) = F D (y) + F A (y) + F R (y) + m j=1 F j G (y) ξj , y(0) = y 0 , (2)
where ξ j , j = 1, . . . , m are independent one-dimensional Wiener processes. The main idea of the new method is to modify the finishing procedure of the standard ROCK2 method, i.e. the final stages of ROCK2 used to achieve the order two of accuracy. We introduce a partitioned RK method, where the diffusion terms F D and advection terms F A are treated explicitly, while the reaction terms F R are treated implicitly. Compared to similar existing stabilized methods, the PIROCK method has the following features:

• for problems with stiff reactions, the number of function evaluations of the reaction terms F R (solved implicitly) is independent of the stage number s used to handle the stiffness of the diffusion terms F D (in contrast, the number of implicit stages in each step of the IRKC method is equal to s);

• for advection dominated problems, the PIROCK method is more efficient than the RKC or ROCK2 solvers as it has better stability in the imaginary direction and requires a number of evaluations of the advection terms that is independent of the stage number of the method; compared to the PRKC method [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF], the PIROCK method has larger stability domains on both the real and the imaginary parts;

• for problems with expensive evaluation of (non-stiff) reaction terms PIROCK is more efficient than RKC [START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF] or ROCK2 [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF] as the number of evaluation of the reaction terms is independent of the stage number of the method; for such problems, it is comparable to the PRKC method [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF] but has larger stability along the negative real axis;

• for problems involving white noise, it is more efficient than previously constructed S-ROCK methods [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF], as PIROCK has a larger mean-square stability domain.

We emphasize that the methods RKC, ROCK2, PRKC are not efficient for problems with severely stiff reactions, while the RKC, IRKC, and ROCK2 methods are not efficient for problems for which the reaction is not stiff but is very costly to evaluate. In addition, notice that a common assumption for all known explicit stabilized integrators is that the differential operator is (nearly) symmetric and has eigenvalues close to the negative real axis. In the case of a non-symmetric differential operator the eigenvalues of the Jacobian of F D are typically located in a sector

S θ = {-ρe iτ ; ρ ≥ 0, -θ ≤ τ ≤ θ} (3)
of the left half complex plane, where θ ≤ π/2 is the angle of this sector. If θ is close to zero then the differential operator is nearly symmetric and the standard stabilized integrators like RKC and ROCK2 can be applied (using more damping if needed, to enlarge the stability domain in the imaginary direction). We shall show that the PIROCK method can also be extended to non-symmetric diffusion operators in sectors (3) for large angles up to θ = π/4. The proposed PIROCK algorithm is versatile and efficient (hence the " swiss-knife") in handling problems such as (1) for various regimes with a single code. It is fully adaptive and requires no tuning from the user.Appropriate error estimators take care of the stiff and nonstiff components of the problems as to deliver a variable step size aiming at an integration error of the size of a tolerance given by the user. While efficient stabilized integrators for special regimes of (1) are available, none has existed until now for the various potential regimes of (1). We also emphasize that PIROCK is more than a simple combination of integrators developed in [START_REF] Verwer | RKC time-stepping for advection-diffusionreaction problems[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations[END_REF][START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF], as the coupling of the different regimes requires new ideas to stabilize the various possible combinations of the dynamics in [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF].

The rest of the paper is organized as follows. In Section 2 we briefly recall the concept of partitioned RK methods, explicit stabilized RK methods and linear stability analysis for numerical integrators. In Section 3 we derive step by step the PIROCK method. Finally a number of numerical experiments taken from benchmark problems for ROCK2, SROCK, RKC, IRKC, PRKC, presented in Section 4 illustrate the flexibility and the performance of PIROCK.

Preliminaries

In this section we present the background concepts needed for the PIROCK integrators, namely a short description of the ROCK2 integrator, linear stability analysis and partitioned Runge-Kutta methods.

ROCK2 methods

We first describe an efficient explicit stabilized method for ODEs (2) for the problem

ẏ = F (y), y(0) = y 0 , (4)
assuming that the Jacobian of F (y) has its eigenvalues close to the real negative axis (with possibly large modulus). We recall that linear stability analysis for stiff ODEs is usually studied on the linear scalar test equation y ′ = λy, y(0) = 1, that can be obtained from linearization and diagonalization (or transformation to Jordan form) of a system of ODEs. The value of λ represents a typical eigenvalue of the Jacobian of the linearized system and a stable mode (i.e. for ℜλ ≤ 0) should be damped by the numerical solver for stability. A RK method applied to the linear test problem leads tp leads to the recursion

y n+1 = R(z)y n , (5)
where z = hλ, and R(z) is a rational function called the stability function of the numerical method (notice that it is a polynomial for an explicit method). This allows to define the stability domain of the numerical integrator as

S := {z ∈ C; |R(z)| ≤ 1}, (6)
and the numerical solution {y n } in ( 5) is bounded for n → ∞ if and only if z ∈ S. If {z ∈ C; ℜz ≤ 0} ⊂ S, then the method is called A-stable, and this is a desirable property for stiff problems. If in addition lim z→∞ R(z) = 0, the method is called L-stable. This latter property is also desirable for stiff problems because it permits a numerical damping of the rapidly decaying eigen modes. Classical explicit methods have a restricted stability domain [START_REF] Abdulle | Coupling heterogeneous multiscale FEM with Runge-Kutta methods for parabolic homogenization problems: a fully discrete space-time analysis[END_REF] along the negative real axis as illustrated in Figure 1 (left picture) for the second order Heun method, where R(z) = 1 + z + z 2 /2, resulting in a severe step size restriction for stiff problems. Stabilized methods allow to increase the stability domains along the negative real axis in an adaptive way, to accommodate stiff dissipative problems without facing step size restriction. An early first order stabilized method [START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF] is based on the following recursion

K 1 = y 0 + h ω 1 ω 0 F (y 0 ), K j = 2h T j-1 (ω 0 ) T j (ω 0 ) F (K j-1 ) + 2ω 0 T j-1 (ω 0 ) T j (ω 0 ) K j-1 - T j-2 (ω 0 ) T j (ω 0 ) K j-2 , j = 2, . . . , s y 1 = K s , (7)
where K 0 = y 0 , T s (cos x) = cos(sx) are the classical Chebyshev polynomials, s ≥ 1 is the number of stages, η ≥ 0 is a damping parameter (discussed below) and ω 0 = 1 + η s 2 , ω 1 = Ts(ω 0 ) T ′ s (ω 0 ) . Applied to the above linear test problem the method [START_REF] Abdulle | Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations[END_REF] gives y 1 = R s (z)y 0 , where z = λh and R s (z) is given by R s (z) = T s (ω 0 + ω 1 z)/T s (ω 0 ) and satisfies |R s (z)| ≤ 1 for all z ∈ (-d s , 0), with d s ≃ C • s 2 , for s large enough, where C depends on the damping parameter η (for η = 0, C = 2). Notice that (7) represents a family of numerical methods (indexed by the stage number s) and the stability domain size d s along the negative real axis grows quadratically with s.

Damping It has been realized in the early development of the stabilized method that it is desirable to introduce a damping of the higher frequencies and include an ellipse around the negative real axis in the stability domain, i.e., to consider

|R s (z)| ≤ r < 1 for all z ∈ (-d s , -ε), (8)
where ε is a small positive number (observe that R s (0) = 1). For instance, for the method (7), we have that r = T s (ω 0 ) -1 < 1 in (8) is a decreasing function of η > 0. This reduces the constant C in d s ≃ C • s 2 , but the boundary of the stability domains does not intersect with the negative real axis (except at the two endpoints). Constructing higher order stabilized methods is a non-trivial task and various strategies have been proposed. We mention the RKC method [START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF] and the DUMKA method [START_REF] Lebedev | Explicit difference schemes with time-variable steps for solving stiff systems of equations[END_REF]. We describe here the ROCK2 method introduced in [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF] and generalized to fourth order in [START_REF] Abdulle | Fourth order Chebyshev methods with recurrence relation[END_REF]. It will be used in what follows. It combines the second order optimal stability domains of [START_REF] Lebedev | Explicit difference schemes with time-variable steps for solving stiff systems of equations[END_REF] with a realization based on a three term recursion formula (similarly as in [START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF]). The stability function of the ROCK2 method reads

R s (z) = P s-2 (z)w 2 (z), (9)
where P s-2 (z) belongs to the family of polynomials {P j } j≥0 (depending on s) orthogonal with respect to the weight function w 2 (x) 2 √ 1-x 2 . The polynomial w 2 , which depends on s, is positive on R and has degree 2. It is chosen such that R s satisfies [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF]

R s (z) = 1 + z + z 2 2 + O(z 3 ), z → 0, (10)
together with a large stability interval [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]. The recurrence relation of the orthogonal polynomials {P j } j≥0 enables the construction of a RK method y 0 → y 1 of order two for (4) based on the following recursion for s ≥ 3,

K 1 = y 0 + µ 1 hF (y 0 ), K j = µ j hF (K j-1 ) -ν j K j-1 -κ j K j-2 , j = 2, • • • , s -2, K * s-1 = K s-2 + σhF (K s-2 ), K * s = K * s-1 + σhF (K * s-1 ), y 1 = K * s -σ(1 -τ /σ 2 )(hF (K * s-1 ) -hF (K s-2 )) (11)
where K 0 = y 0 . The parameters µ j , κ j , ν j (depending on s) are obtained from the three-term recurrence relation [5, eq. ( 24)-( 25)] of the orthogonal polynomials {P j } j≥0 , while σ, τ (which also depend on s) satisfy w 2 (z) = 1+2σz +τ z 2 and are chosen such that (10) holds. We notice that the polynomials P j (z) are the stability functions of the internal stages K j , j = 1, . . . , s-2.

We have for 3 ≤ s ≤ 200, σ ∈ (0.367, 0.410) and τ ∈ (0.2, 0.4). The ROCK2 method satisfies (8) with d s ≃ 0.81 • s 2 for a damping r = 0.95, where r is defined in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]. We have thus now a family of second order methods, whose stability domains increases quadratically with the stage number s (see Figure 1 where d 13 ≃ 135.1 ≃ 0.81 • 13 2 ). As shown in [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF], due to the aforementioned stability behavior this method is competitive with implicit solvers for diffusion problems, while remaining explicit. As the parameters µ j , ν j , κ j , τ j , σ j depend on r and are computed once and imported into the computer code, changing the value of r requires a priori to recompute these parameters.

We explain next a procedure that allows to change the damping with the already precomputed parameters included in the ROCK2 code http://anmc.epfl.ch. This idea, first used in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF] for weak second order stochastic stabilized methods, is to consider for (4) the following scheme for a fixed scalar parameter α:

K 1 = y 0 + αµ 1 hF (y 0 ), K j = αµ j hF (K j-1 ) -ν j K j-1 -κ j K j-2 , j = 2, • • • , s -2, K * s-1 = K s-2 + σ α hF (K s-2 ), K * s = K * s-1 + σ α hF (K * s-1 ), y 1 = K * s -σ α (1 -τ α /σ 2 α )(hF (K * s-1 ) -hF (K s-2 )), (12)
where K 0 = y 0 . Notice that for α = 1, we recover the original ROCK2 method [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF]. Applied to the linear test problem this method yields

y 1 = P s-2 (αz)(1 + 2σ α z + τ α z 2 ) =: R s,α (z), (13)
and it can be easily verified (see [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]Lemma 3.2]) that the method ( 12) has second order for the system of ODEs (2) for any α, provided that The oscillating polynomial is R s,α (z), the dotted lines correspond to the stability polynomials of the internal stages P j (αz), j = 1, . . . , 11. The stability functions P s-1 (αz), P s (αz) of the two supplementary stages K s-1 , K s are depicted in solid (red) lines oscillating with a small amplitude.

σ α = 1 -α 2 + ασ, τ α = (α -1) 2 2 + 2α(1 -α)σ + α 2 τ. (14)
For the PIROCK method we will need two additional stages associated to the family of orthogonal polynomials {P j } j≥0 ,

K j = αµ j hF (K j-1 ) -ν j K j-1 -κ j K j-2 , j = s -1, s. (15)
Applied to the linear test problem we obtain

K s-1 = P s-1 (αz), K s = P s (αz). ( 16
)
These polynomials decay faster to zero than P s-2 (αz) for z close to zero because P ′ s-2 (0) < P ′ s-1 (0) < P ′ s (0), as illustrated in Figure 2. It should be noted that the stability function [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF] of the ROCK2 methods oscillates for real p ≤ 0 around 0 (see Figure 2) while the stability functions of other stabilized methods such as [START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF][START_REF] Verwer | RKC time-stepping for advection-diffusionreaction problems[END_REF]) oscillates around the value a s > 0. In turn, these methods have less favorable damping, and the ellipse that can be included in the stability domain is smaller than the corresponding ellipse for the ROCK2 methods.

the RKC that reads R s (z) = a s + b s T s (ω 0 + ω 1 z), for appropriate ω 0 , ω 1 and a s = 1 -b s T s (ω 0 ), b s = T ′′ s (ω 0 )/(T ′ s (ω 0 )) 2 (see
Implementation of the ROCK2 method We emphasize that an efficient implementation relies on both stage and time step adaptivity as summarized in the following algorithm. Given the value y 0 , the step size h, and the required accuracy tol (prescribed by the user), 1. perform an integration step y 0 → y 1 ;

2. estimate the local error err := y 1 -ŷ1 ;

3. determine a new step size h new such that err ≃ tol based on a step size control strategy; 4. estimate the spectral radius ρ of the Jacobian of F and choose a stage number such that

h new • ρ ≃ C • s 2
new (where C depends on α);

back to 1.

A few comments are in order. First, ŷ1 is the numerical solution obtained from a so-called embedded method that allows to estimate the local error. This numerical solution given by ŷ1 = K * s in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF] does not involve any significant computational overhead and is obtained from the function evaluations needed for y 1 [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]II.4]). Second, sophisticated procedures are available to compute h new that also take into account previously computed step sizes (memory step size selection) [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF]IV.8]. These procedures have proven successful for ROCK2 methods [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF]. Third, a precise estimation of the spectral radius (which would be expensive) is not needed. It is sufficient to consider an upper bound that can be provided by the user or obtained from a power methods, using the function evaluations needed for the internal stages [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF]. Hence this procedure comes with a negligible cost.

Applicability and limit of the stabilized methods Stabilized methods as the ROCK, RKC or DUMKA methods have been primarily introduced for the solution of large dissipative systems (i.e., with eigenvalues of the Jacobian close to the negative real axis). They can be applied to systems such as (1) but become inefficient for problems with a large advection and for problems with a very stiff reaction. Increasing the damping (an idea also used for stochastic methods [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF]) allows to include larger ellipses in the stability domain of the method [START_REF] Verwer | RKC time-stepping for advection-diffusionreaction problems[END_REF]. But this comes at the cost of reducing the length of the stability domain along the negative real axis and even so, the size in the imaginary direction of the ellipse that can be included in the stability domain is limited (see Section 3.2). Furthermore, another drawback is that the number of function evaluations is proportional to the number of stages.

Partitioned Runge-Kutta methods

For the time integration of ( 1), an m-stage partitioned RK method2 (which maps y 0 to y 1 ) is given by

K i = y 0 + h m j=1 a ij F D (K j ) + h m j=1 âij F A (K j ) + h m j=1 āij F R (K j ), i = 1, . . . , m, (17)
y 1 = y 0 + h m i=1 b i F D (K i ) + h m i=1 bi F A (K i ) + h m i=1 bi F R (K i ). ( 18
)
Remark 2.1 For a non-autonomous problem, we need to evaluate the stages at discrete time t 0 + c i h, t 0 + ĉi h, t 0 + ci h, i = 1, . . . , m and we assume as usual that the conditions

c i = m j=1 a ij , ĉi = m j=1 âij , ci = m j=1 āij , i = 1, .
. . , m are satisfied [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]. The method is said to have order p if for sufficiently smooth problems (1), we have the local error bound

y(t 0 + h) -y 1 ≤ Ch p+1 . ( 19
)
Order conditions for a partitioned RK method ( 18) are algebraic conditions on the coefficients such that ( 19) is satisfied for a given p. As we will construct second order methods, we recall the corresponding order conditions, namely

m i=1 b i = m i=1 bi = m i=1 bi = 1, m i,j=1 b i a ij = m i,j=1 bi âij = m i,j=1 bi āij = 1 2 , ( 20
) m i,j=1 b i âij = m i,j=1 b i āij = m i,j=1 bi a ij = m i,j=1 bi âij = m i,j=1 bi a ij = m i,j=1 bi āij = 1 2 . ( 21
)
The conditions [START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF] ensure that each individual method (for F D , F A , F R ) has second order, while conditions [START_REF] Van Der Houwen | On the internal stage Runge-Kutta methods for large m-values[END_REF] ensure that the coupling of the different terms has the right accuracy.

It is customary to write the coefficients of a RK method in a so-called Butcher tableau

Applicability and limit of existing stabilized partitioned methods In the PRKC method [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF], a partitioned procedure allows to make the number of evaluations of the advection (or a non-stiff reaction) independent of the stage number. These methods are efficient for diffusion dominated advection-diffusion reaction problems with non-stiff reaction, particularly if the reaction or advection term is costly to evaluate. For dominant advection, the performance of the PRKC method deteriorates (as for the ROCK2 or RKC methods) because of the limited height of the ellipse that can be included in the stability domains of the method. On the other hand the implicit-explicit IRKC method [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF] can accommodate stiff reactions (done implicitly) but the number of nonlinear systems to be solved depends on the stage number used to treat the stiffness arising from the diffusion. As for the PRKC method, the performance of the IRKC also deteriorates for problems with dominant advection.

The PIROCK method

In this section we derive our new PIROCK method and analyze its stability behavior. We also discuss a posteriori error estimates via embedded RK methods and step size control.

Derivation of the PIROCK method

We explain step by step the construction of our new methods for the solution of (1).

Step1: choosing the F A and F R method We consider the following classical basic methods F A -method F R -method 0

1 3 1 3 2 3 2 3 1 4 0 3 4 γ γ 1 -γ 1 -2γ γ 1 2 1 2 (23)
where γ = 1 -√ 2/2. A 3-stage third order explicit method is taken for the advection (so that a non-empty portion (-i √ 3, i √ 3) of the imaginary axis is included in the stability domain of the F A method) and a 2-stage second order singly diagonally implicit RK method for the reaction. This latter method is L-stable and can be efficiently implemented: due to the advection F A 2.5 2.0 1.5 1.0 0.5 0.0 0.5 3 diagonal structure of RK coefficients, a single LU factorization needs to be done only once per step if using a quasi-Newton method [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF].

Step 2: Diffusion step and coupling In the PIROCK method, we perform the diffusion step first and introduce the advection and reaction steps as a"finishing procedure". This diffusion step needs then to be coupled with the F A and F R methods (see the coupling order conditions ( 21)). We explain the coupling by choosing the (explicit) Heun method for the diffusion steps. This choice is only done here temporarily to simplify the presentation and we emphasize that the Heun method will later be replaced by the ROCK2 method. For the time being we thus consider 0 1 1

1 2 1 2 (24)
The partitioned method with the Heun method for the diffusion step and the F A , F R methods as finishing procedure together with suitable coupling conditions reads

F D -method 0 1 1 0 1 1 0 0 2 3 2 3 1 1 1 2 1 2 -1 2(1-2γ) 0 1 2(1-2γ) 0 0 0 ( 25
)
F A -method F R -method 0 0 0 1 1 1 -2γ 1 -2γ 1 3 1 3 2 3 2 3 0 0 0 1 4 0 0 0 3 4 0 0 0 γ γ 1 -γ 1 -2γ γ 1 -γ 1 -γ 0 2-γ 3 2 3 -γ 2γ 3 0 0 0 1 2 1 2 0 0 0 0 (26)
where γ = 1 -√ 2/2. Of course this choice of coupling is by no mean unique. We make the above choice of coupling because it has the same cost as each of the 3 non-partitioned RK methods [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF], [START_REF] Verwer | Convergence properties of the Runge-Kutta-Chebyshev method[END_REF], with the exception that two extra evaluations of the diffusion function F D are needed at each step. The coefficients 2 3γ, 2γ 3 in the F R -method are chosen to make the corresponding internal stage L-stable with respect to the reaction F R .

Step 3: Introducing a perturbation of the initial value for F A and F R The above partitioned method ( 25)-( 26) takes y 0 as initial value for the advection and the reaction part. The idea of the PIROCK algorithm is to use internal stages of the ROCK2 method for the stabilization of the coupling with the reaction and advection terms. It is thus of interest to start with a suitably perturbed initial value, say y 0 + δhF D (y 0 ), where δ is a fixed real. The partitioned method is then given by the following F D -method together with the unchanged F A and F R tableau in [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF],

F D -method 0 δ δ δ δ δ + β δ β δ δ δ δ δ + 2 3 β δ 2 3 β 1 1 1 2 0 -1 2(1-2γ) 0 1 2(1-2γ) 0 0 1 2 (27)
A new parameter β has been introduced to maintain the second order accuracy for the coupling of the diffusion-advection and diffusion-reaction terms with an arbitrary value of δ. It can be again easily verified that the above method has second order (i.e., [START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF] and ( 21) hold) for arbitrary δ provided

β = 1 -2δ. ( 28
)
Efficient resolution of the nonlinear systems We notice that for each step of the above method ( 26)-( 27), we have to solve two nonlinear systems of the type

y = c j + γhF R (y), (29)
where c j (j = 1, 2) is independent of y. The resolution can be done efficiently because the nonlinear systems for each spatial component are independent. This means that (29) consists of decoupled systems in low dimension n P DEs ≪ n, where n P DEs is usually given by the number of scalar PDEs in a system of PDEs, and n is the dimension of the system after spatial discretization. It can be performed using a quasi-Newton method J R (y k+1y k ) = -y k + c j + γhF R (y k ) where the LU decomposition of the matrix

J R = I -γh ∂F R ∂y (c 1 ) (30)
has to be computed only once per step size. Notice that in contrast, the IRKC method requires at each time step the resolution of s nonlinear systems of the form (29) where s is the number of stages. In addition, the LU decomposition of the matrix J R has to be computed s times per step of IRKC because the systems (29) use different values of the coefficient γ in each stage.

Using J -1 R for stabilization It has been noticed by Shampine (as presented in [13, IV.8]) that J -1 R can also be used for stabilization purpose (originally used for the error control). In the PIROCK method, J -1 R will also be used to further stabilize the coupling of the diffusion step with the F A and F R methods. We emphasize that computing J -1 R v, for a given vector v represents a negligible cost as the LU factorization of J R is already available from the computation of the implicit stages (29). Likewise computing a higher order power such as J -2 R comes also with a negligible cost.

Step 4: Stabilizing the partitioned method and definition of the PIROCK method Although based on good methods for the advection and reaction part of the ODE (1), the method ( 26)-( 27) is useless because of its poor stability for the diffusion part. The PIROCK method is based upon the following modification of the method ( 26)-( 27):

• replacement of the Heun method [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF] with the ROCK2 method [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF] in which we compute the two additional stages K s-1 , K s in (15);

• choice of the starting value for the advection and reaction terms at K s (or K s-1 );

• use of J -ℓ R = (I -γh ∂F R ∂y ) -ℓ , ℓ = 1, 2 to stabilize the couplings F D -F R , F A -F R .
PIROCK integrator The partitioned integrator for ( 1) is defined by the following algorithm for s ≥ 3:

Diffusion stabilization procedure K 1 = y 0 + αµ 1 hF D (y 0 ), K 0 = y 0 , K j = αµ j hF D (K j-1 ) -ν j K j-1 -κ j K j-2 , j = 2, • • • , s -2 + ℓ (ℓ = 1 or 2)
Finishing procedure for diffusion

K * s-1 = K s-2 + σ α hF D (K s-2 ), K * s = K * s-1 + σ α hF D (K * s-1 ), Starting value for advection-reaction K = K s-2+ℓ ,
Finishing procedure for advection-reaction and coupling

K s+1 = K + γhF R (K s+1 ), K s+2 = K + βhF D (K s+1 ) + hF A (K s+1 ) + (1 -2γ)hF R (K s+1 ) + γhF R (K s+2 ), K s+3 = K + (1 -2γ)hF A (K s+1 ) + (1 -γ)hF R (K s+1 ), K s+4 = K + 1 3 hF A (K s+1 ), K s+5 = K + 2β 3 hF D (K s+1 ) + 2 3 hJ -1 R F A (K s+4 ) + ( 2 3 -γ)hF R (K s+1 ) + 2γ 3 hF R (K s+2 )
Computation of the integration step y 0 → y 1

y 1 = K * s -σ α (1 -τ α /σ 2 α )(hF D (K * s-1 ) -hF D (K s-2 )) + 1 4 hF A (K s+1 ) + 3 4 hF A (K s+5 ) + 1 2 hF R (K s+1 ) + 1 2 hF R (K s+2 ) + J -ℓ R 2 -4γ (hF D (K s+3 ) -hF D (K s+1 )) , (31)
where µ j , ν j , κ j are the same coefficients as for the standard ROCK2 method [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF], the coefficients σ α , τ α are defined in [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF],

γ = 1 - √ 2/2, β = 1 -2αP ′ s-2+ℓ (
0) and J R = Iγh ∂F R ∂y (K s ). We notice that the starting value for the advection and the reaction term requires the stage K of the ROCK2 method. We shall consider two choices of parameters α ≥ 1 (see [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]) and with corresponding values ℓ = 1 or 2 in (31). The first choice is

α = 1, ℓ = 2, ( 32
)
which permits to recover the standard ROCK2 method in the absence of advection and reaction terms (F A = F R = 0), with a (close to optimal) stability domain along the negative real axis of size d s ≃ 0.81 • s 2 . The second choice that we shall consider is

α = 1/(2P ′ s-1 (0)), ℓ = 1. ( 33
)
This is a regime with larger damping (α > 1) than in the standard ROCK2 method and allows to includes a larger ellipse in the stability domain of the method for the coupling F D -F A . It is thus suited for advection dominated problems. Observe also that for this choice of damping, β = 0.

Remark 3.1 The choice of values ℓ = 1 or 2, used for the definition of the internal stage K in (31), is made to provide a better damping in the stability domains close to the origin. Since P ′ s-2+ℓ (0) > P ′ s-2 (0), the stability polynomial P s-2+ℓ associated to K decays faster close to the origin than the one associated to K s-2 and involved in the ROCK2 method (see Figure 2). This permits to avoid a gap in the stability domains close to the origin for the F D -F A and F D -F R couplings. Also, the choice of the negative power J -ℓ R in (31) is made to avoid a gap in the stability domains involving the F D -F R coupling, while the term

J -1 R in K s+5 is used to stabilize the F A -F R coupling.
Complexity Compared to the standard ROCK2 method [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF], the PIROCK method requires at each time step 2 + ℓ additional evaluations of F D (with ℓ = 1 or 2), 3 evaluations of F A , and 2 resolutions of a nonlinear system of the form (29) that are usually computationally cheap (see (30)). In contrast, the IRKC method requires the resolution of s nonlinear system of the form (29) (with different values of γ at each internal stage) where s is the number of stages of the IRKC method.

Accuracy We next verify that the PIROCK method is indeed second order accurate. Theorem 3.2 The method (31) has second order of accuracy for (1).

Proof. We notice that (31) is a perturbation of the second order method ( 26)-( 27), now involving J -2 R , with the second order Heun method replaced by the second order ROCK2 method [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF] and with K 2 replaced with K. A simple expansion using [START_REF] Lebedev | Explicit difference schemes with time-variable steps for solving stiff systems of equations[END_REF] shows that

K = y 0 + αhP ′ s-2+ℓ (0)F D (y 0 ) + O(h 2 ),
and thus, in view of (28) where δ = αP ′ s-2+ℓ (0), we need to define β = 1 -2αP ′ s-2+ℓ (0). Finally, involving the matrix J -1 R only introduces a perturbation of size O(h 3 ) in y 1 because J -1 R is a O(h)-perturbation of the identity and the proof is complete. Remark 3.3 Theorem 3.2 states the second order of accuracy of the PIROCK method for systems of ODEs. Notice however that the analysis of the dependency of the error constant on the stiffness of the problem is out of the scope of the paper. We mention that the second order error constant of the F R -method itself, defined in [START_REF] Verwer | Convergence properties of the Runge-Kutta-Chebyshev method[END_REF], can be shown to be independent of the stiffness using the theory of B-convergence (see [10, p. 31]). The error constant of the RKC method is shown to be independent of the stiffness for a class of linear problems in [START_REF] Verwer | Convergence properties of the Runge-Kutta-Chebyshev method[END_REF]. This is extended to classes of explicit stabilized methods (including ROCK2) in [START_REF] Abdulle | Coupling heterogeneous multiscale FEM with Runge-Kutta methods for parabolic homogenization problems: a fully discrete space-time analysis[END_REF] in the context of multiscale homogenization problems.

Extension to non-symmetric diffusion operators. We explain now a simple modification of the PIROCK algorithm to treat the case of a non-symmetric diffusion operator. Assume that the diffusion operator can be decomposed as

F D (y) = F D S (y) + F D A (y). ( 34
)
where the eigenvalues of the Jacobians of F D S and F D A are located respectively on the negative real axis and the imaginary axis. For instance, in the case of the differential operator div(a∇•), where a is a non-symmetric N × N tensor, the corresponding decomposition reads div(

a + a T 2 ∇•) + div( a -a T 2 ∇•).
Then, we introduce the following modifications in the definition of the PIROCK algorithm (31) for the damping (33). First, we replace F D by F D S in all stages. Next, we replace F A by F A + F D A in the formulas for K s+2 , K s+3 , and y 1 . Finally, in the stage K s+5 , we substitute F A (K s+4 ) by F A (K s+4 ) + F D A (K s+1 ). Notice that this modification remains explicit. The stability analysis of this modification for non-symmetric tensors is addressed in the next Section 3.2.

Extension to stochastic problems. The PIROCK method can also be used for stiff meansquare stable Itô stochastic SDEs of the form [START_REF] Abdulle | Fourth order Chebyshev methods with recurrence relation[END_REF]. The stochastic version of PIROCK, a map from y 0 → y * 1 , is obtained by modifying the last line of the PIROCK integrator (31) as follows

y * 1 = y 1 + m j=1 F j G (K * s+1 )∆W j , (35)
where y 1 is defined in (31), and using the supporting value

K * s+1 = J -1 R (K s+1 + βhF D (K s+1 )),
where ∆W i ∼ N (0, h) are independent Wiener increments, and J R is the block diagonal matrix in (30) whose LU decomposition has been already computed. The term βhF D (K s+1 ), already computed before, is used to avoid a small gap in the mean-square stability domain close to the origin in the F D -F G coupling. The multiplication with J -1 R is used to stabilize the F R -F G coupling and make it mean-square A-stable with respect to this coupling (see below). It is easily seen that the integrator (35) has weak order 1 and strong order 1/2 for general Itô stochastic SDEs of the form (2). We refer to [START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF] for accuracy concepts and more details on stabilized stochastic methods. Notice that an explicit stabilized method of weak order 2 with extended mean-square stability domains was recently constructed in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF].

Stability analysis

We analyze here the stability of the various couplings in the method (31). To analyse the coupling of the diffusion term F D respectively with F A ,F R ,F G , we shall consider scalar linear test problems, where F D (y) = λy is linear and λ ∈ R -corresponds to the eigenvalue of a symmetric diffusion operator, and F A (y) = iµy with µ ∈ R, F R (y) = ρy with ρ < 0, F G = σy with σ ∈ R. We emphasize that such test problems only give insight on the stability of the PIROCK method, because the operators involved in problem [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF], cannot in general be put simultaneously into diagonal form (even in the case where they are linear). where λ ∈ R -, µ ∈ R. This test equation is relevant for advection-diffusion equations

∂ t u(x, t) + ν∂ x u(x, t) = ∂ xx u(x, t),
x ∈ (0, 1), t > 0, (36)

as the linear systems of ODEs arising from the (Fourier) spatial discretization on a uniform grid of size ∆x has eigenvalues that belong to the ellipse with half-height a = ν∆x -1 , and width d = 4∆x -2 . To study the stability of a given advection-diffusion integrator, a natural criteria is to search for the largest ellipse included in the stability domain of the method for the scalar test problem ẏ = λy + iµy, where (p, q) belongs to the ellipse (37), with p = hλ, q = hµ. Precisely, we search for the ellipse first with the largest width d, and then with the largest half-height a.

2p d + 1 2 + q a 2 = 1 (
The PIROCK method (31) then yields y n+1 = R(p, q)y n , where p = hλ, q = hµ, R(p, q) = R s,α (p) + P s-2+ℓ (αp) -

q 2 2 + i(q - q 3 6 + (β + 1 2 )pq) , (38)
and R s,α is the stability function ( 13) of the ROCK2 method. We plot in Figure 4 the corresponding stability domains in the (p, q)-plane, for s = 13 and the two choices of damping (32),(33). We also plot the largest stability ellipses (37) included in the stability domains.

In Figure 5, we plot the ellipse parameters a s , d s as functions of s for the choices (32) (solid lines) and (33) (dashed lines). We observe that these ellipse parameters as functions of s grow quadratically for the width d s and linearly for the half-height a s as

d s ≃ 0.
) 40
Notice that in [START_REF] Verwer | RKC time-stepping for advection-diffusionreaction problems[END_REF], the RKC method is considered with a large value η = 10 of the damping parameter. However, it reduces down to d s ≃ 0.34 • s 2 the length of the stability domain along the negative real axis, and only a growth of size O( √ s) is obtained for the ellipse half-height a s . In contrast, in the PRKC approach [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF], the coefficients of the PRKC method are chosen such that the stability domain of the diffusion-advection coupling includes in the (p, q)-plan

a rectangle [-d s , 0] × [-b, b] of half-height b = 1.7273 (close to √ 3
), independently of the number of stages s, where d s ≃ 0.65 • s 2 and while keeping the standard damping η = 2/13 of the standard RKC method. The stability constraint for advection of the PRKC method then reads h ≤ b/ρ A .

Remark 3.4

The linear growths (39),(40) of the ellipse half-height a s is a feature of PIROCK that takes advantage of the diffusion terms to stabilize the possibly large advection terms. Indeed, consider in particular the damping (33). Substituting hρ D = 0.43 • s 2 into hρ A ≤ max(0.53 • s, √ 3) yields the stability sufficient condition

hρ A ≤ max 0.53 2 0.43 ρ D /ρ A , √ 3 .
Considering again the advection-diffusion test problem (36) on a uniform grid of size ∆x, using ρ A = ν∆x -1 and ρ D = 4∆x -2 , we have that the PIROCK method is stable if hν ≤ max(2.61 • ν -1 , 1.73 • ∆x) and this stability constraint on the stepsize h is independent of the spatial grid size small enough (∆x ≤ ν -1 ), in contrast to the other stabilized explicit integrators RKC, PRKC and ROCK2.

Diffusion-reaction coupling For the F D -F R coupling, we consider the test problem

y ′ (t) = λy(t) + ρy(t)
where λ, ρ ≤ 0. Applied to this test problem, the PIROCK method yields the relation y n+1 = R(p, r)y n where p = hλ, r = hρ, and

R(p, r) = R s,α (p) + P s-2+ℓ (αp) r -γ 2 r 2 (1 -γr) 2 + βrp 2(1 -γr) 2 + pr 2(1 -γr) 1+ℓ
It can be checked that the stability domain in the (p, r) plane contains the subdomain (-d s , 0) × R -for both choices of damping (32),(33) (see Figure 6 for s = 13). This means that the PIROCK method is unconditionally stable with respect to ρ ≤ 0. Notice that the stability domain also includes points where r > 0 is large enough, this due to the L-stability of the F R -method. where λ ∈ R -, µ ∈ R. The values λ, iµ represent real and imaginary eigenvalues of the Jacobian of F D S , F D A , respectively. Applying the modified PIROCK method to the above test equation yields, for p = hλ, q = hµ, the stability function

R(p, q) = R s,α (p) + P s-1 (αp) - q 2 2 + i(q + 1 2 pq) ,
which is identical to (38) with damping (33), with the exception that the q 3 term no longer appears (thanks to an appropriate modification made in the stage K s+5 ). We have computed numerically for 3 ≤ s ≤ 200 the largest sector inscribed in the stability domain, given by the angle (see left picture of Figure 7)

θ s := sup{θ ; S θ ∩ A s ⊂ S}, where S = {p + iq ∈ C; |R(p, q)| ≤ 1}, A s = {p + iq ∈ C; -d s ≤ p ≤ 0}
and S θ is the sector defined in (3). We observe that θ s ≥ π/4 for all s ≥ 10. The corresponding stability domain is plotted in Figure 7 (right picture) for s = 13, where θ 13 ≃ 0.26π. Stochastic stability For SDEs such as (2), a widely used concept of stability is the meansquare stability considering the scalar test equation [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF][START_REF] Higham | A-stability and stochastic stability mean-square stability[END_REF][START_REF] Burrage | Numerical methods for strong solutions of stochastic differential equations: an overview[END_REF][START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF] ẏ = λy + σy ξ,

y(0) = 1, (41)
with fixed complex scalar parameters λ, σ. The exact solution of (41) is said mean-square stable if and only if lim t→∞ E |y(t)| 2 = 0 and for the problem (41) it can be characterized as the set

S M S = {(λ, σ) ∈ C 2 ; ℜ(λ) + 1 2 |σ| 2 < 0}. ( 42
)
Applied to the SDE (41), a numerical integrator yields the difference equation [START_REF] Higham | A-stability and stochastic stability mean-square stability[END_REF]

y n+1 = R(p, q, ξ n )y n , (43)
where p = λh, q = σ √ h, ξ n ∼ N (0, 1) are independent random variables, and the mean-square stability domain of a numerical method can be characterized as

lim n→∞ E(|y n | 2 ) = 0 ⇐⇒ S M S num := (p, q) ∈ C 2 ; E(|R(p, q, ξ)| 2 ) < 1 . (44)
If we restrict (p, q) ∈ R 2 then we consider the portion of the true mean-square stability domain (42), namely

S M S x = {(p, q) ∈ (-x, 0) × R ; p + 1 2 |q| 2 < 0}, (45)
and define for a given method

ℓ = sup{x > 0 ; S M S x ⊂ S M S num }, d = sup{x > 0 ; (-x, 0) × {0} ⊂ S M S num }, ( 46
)
where d is the size of the stability domain along the deterministic p-axis (observe that d ≥ ℓ). For the diffusion-noise coupling F D -F G , the PIROCK method applied to the scalar test problem (41) yields for the mean-square stability function For s = 13, this yields ℓ 13 ≃ d 13 ≃ 135.15, as illustrated in Figure 8 (left picture) where it can be observed that the portion S M S ℓ 13 below the dotted line is included in the stability domain. We have checked numerically that the mean-square stability domain S M S num in (44) contains a portion (45) of the true mean-square stability region of size

E(|R(p, q, ξ)| 2 ) = |R s,α (p)| 2 + |P s-2+ℓ (αp)| 2 |1 + βp| 2 |q| 2 . (47)
ℓ s ≃ d s ≃ 0.81 • s 2 using (32), ℓ s ≃ d s ≃ 0.43 • s 2 using (33),
as illustrated in Figure 9. Notice that the S-ROCK method [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF], with weak order one and strong order 1/2, and based on damped first order Chebyshev methods [START_REF] Abdulle | Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations[END_REF], has a shorter mean-square stability domain of size ℓ s ≃ 0.33 • s 2 . Considering the linear test problem (41) for the reaction-noise coupling F R -F G , we obtain the mean-square stability function

E(|R(r, q, ξ)| 2 ) = |1 + (1 -2γ)r| 2 |1 -γr| 4 + |q| 2 |1 -γr| 4 . (48)
The corresponding mean-square stability domain S M S num in (44) can be shown (adapting the proof of [START_REF] Abdulle | Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations[END_REF]Theorem 3.2]) to contain as a subset the exact mean-square stability domain (42), as illustrated in Figure 8 (right picture). We say that the reaction-noise coupling F R -F G is mean-square A-stable.

Numerical experiments and comparison with other explicit stabilized integrators

We compare our PIROCK integrator with the known time integrators that make an explicit stabilized treatment of the diffusion: RKC, IRKC, PRKC, ROCK2. All these integrators are implemented in FORTRAN and run on the same machine for CPU time comparisons. Unless specified, we use default parameters for all integrators.

4.1 The 2D Brusselator problem with a highly stiff reaction We consider a modification of the Brusselator problem [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF] in 2D with a stiff reaction,

∂u ∂t = ν∆u + A + u 2 v -(B + 1)u, ∂v ∂t = ν∆v + Bu -u 2 v, (51)
with x ∈ (0, 1) 2 , t ∈ (0, 2), A = 1.3, B = 2 • 10 7 , with periodic boundary conditions u(x 1 + 1, x 2 , t) = u(x 1 , x 2 , t) = u(x 1 , x 2 + 1, t) and the initial condition

u(x, 0) = 22x 2 (1 -x 2 ) 3/2 , v(x, 0) = 27x 1 (1 -x 1 ) 3/2 . ( 52
)
For the diffusion, we consider the parameter ν = 10 -1 . We discretize u, v in space with two n × n uniform meshes, where n = 200. Here, the vector field F D corresponds to the discretized Laplacians and the vector field F R consists of n 2 decoupled reaction ODEs in dimension n P DEs = 2. In Figure 10 and Table 1, we compare the integrators PIROCK and IRKC for the tolerances tol = 10 -r , r = 1, . . . , 5. We take the initial step size h = 10 -3 for PIROCK, while the step size is automatically selected for IRKC. Notice that the standard integrators RKC,ROCK2 and PRKC cannot be reasonably used because the reaction is too stiff (B ≫ 1). Indeed, the spectral radius of the Jacobians are for the diffusion ρ D = 8n 2 = 3.2 • 10 5 , and for the reaction ρ R ≃ 10 7 . We point out that large reaction terms, as chosen in the above example, arise in practical applications when very stiff reactions are modeled (see for example [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF]). We observe that the number of evaluations of the diffusion function F D (column "F D evals") is reduced by two orders of magnitude for IRKC and PIROCK compared to the standard RKC method. In addition, the PIROCK method compared to IRKC has a considerably reduced number of evaluations of the reaction function F R and its Jacobian (column "∂F R evals"), by one to two orders of magnitude. This is because the PIROCK method has only two implicit internal stages per step, whereas the IRKC method has all its s internal stages implicit. To check the accuracy of the integrators, we computed independently two reference solutions with tolerance tol = 10 -7 using PIROCK (95 sec. of CPU time) and IRKC (1662 sec. of CPU time), and we checked that the observed errors are identical for both reference solutions.

A radiation-diffusion problem with a highly stiff reaction

We consider the combustion problem taken from [START_REF] Mousseau | Physics-based preconditioning and the newton-krylov method for non-equilibrium radiation diffusion[END_REF], with spatial discretization taken from [START_REF] Hundsdorfer | Numerical solution of time-dependent advection-diffusion-reaction equations[END_REF]Chap. 5] and considered for numerical illustrations in [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF]Sect. 4]. This 2D model used in laser fusion applications describes the radiation energy E and the material temperature T , defined on the unit square domain (0, 1) 2 , given by two non-linear diffusion equations with a highly stiff reaction, for t > 0,

∂E ∂t = ∇ • (D 1 ∇E) + σ(T 4 -E), ∂T ∂t = ∇ • (D 2 ∇T ) -σ(T 4 -E), (53)
where σ = Z 2 /T 3 , D 1 = 1/(3σ + |∇E|/E), and D 2 = kT 5/2 with k = 0.005. Here, Z = Z(x) represents the atomic mass number which has the spatial inhomogeneity Z(x) = Z 0 = 10 if max j |x j -1/2| ≤ 1/6, and Z(x) = 0 else. For the boundary conditions, we consider homogeneous Neumann boundary conditions for T at all boundaries and for E on the top and bottom boundaries x 2 = 0, 1. For the left and right boundaries we consider the nonlinear conditions

1 4 E -1 6σ ∂E ∂x = 1 -x 1 , at x 1 = 0, 1.
The initial values are the constant functions E(x, 0) = 10 -5 , T (x, 0) = E(x, 0) 1/4 , which is not an equilibrium of the PDEs because the nonlinear boundary condition at x 1 = 1 is not satisfied for E.

We consider for the space discretization two n × n grids with n = 100. The spectral radius are taken from [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF]Sect. 4] and given by ρ D = 8.6 • 10 4 , ρ R = 6 • 10 6 . Similarly to the Brusselator problem considered in Section 4.1, the reaction is too stiff for a reasonable application of the standard methods ROCK2, RKC, and the partitioned method PRKC, and an implicit treatment of the reaction term is needed to avoid a severe step size restriction. We thus compare the PIROCK method only with the IRKC method. For the reference solution and plot of the numerical solution at final T = 3 in Figure 12, we used PIROCK with the tolerance tol = 10 -8 . We again compared this reference solution with the a reference solution given by IRKC (using the same tolerance) and checked that the significant digits where identical. The solution consists of a front moving from the left boundary (where the boundary condition is not satisfied) to the right boundary of the domain. It can be seen in Figure 12 that the inhomogeneity located at the center of the domain, where the reaction is very stiff while the diffusion is smaller, greatly affects the shape of the solution (both E and T ).

We take the initial step size h = 10 -5 for PIROCK, and compare the integrators PIROCK and IRKC in Figure 11. Notice that for PIROCK with tol = 10 -2 , the step size selection produced too large steps and the numerical solution became negative (the scalar D 2 = kT 5/2 is then no longer properly defined). We thus choose the largest tolerance for PIROCK as tol = 10 -2.5 and consider in Figure 11 the tolerances tol = 10 -r/2 , r = 5 . . . 9 (PIROCK) and tol = 10 -r/2 , r = 2 . . . 6 (IRKC). It can be seen in Figure 11 comparing the L 2 error versus CPU time, that the PIROCK integrator is more efficient than RKC for the tolerances tol = 10 -r , r = 2.5, 3, 3.5. Notice that the spatial L 2 -error for n = 100 is estimated as 3 • 10 -2 [24, Sect. 4.3.2] and that a uniform discretization in space is certainly not optimal (given the sharp front of the solution).

The 2D Brusselator problem with large advection

We consider the Brusselator problem in 2D where we add advection terms,

∂u ∂t = ν∆u + µU • ∇u + A + u 2 v -(B + 1)u,
∂v ∂t = ν∆v + µV • ∇v + Bu -u 2 v, (54)
with x ∈ (0, 1) 2 , t ∈ (0, 1), A = 1.3, B = 1 ,U = (-0.5, 1) T , V = (0.4, 0.7) T , with periodic boundary conditions and the initial condition (52). For the diffusion, we consider the parameter ν = 10 -2 . For the advection terms, we consider the values µ = 1 and µ = 10 -1 , respectively. We discretize u, v in space with two n × n uniform meshes, where n = 400. For the integrators PRKC and PIROCK, we consider the partitioning where F D corresponds to the discretized Laplacian ∆u, ∆v, and F A corresponds to the other terms (advection and reaction). For µ = 1, the solutions at point x = (0, 0) of PIROCK are plotted in Figure 13 (left pictures) together with the corresponding variable stepsizes along the integration interval (right picture) for the tolerances tol = 10 -1 , 10 -2 . As predicted in Remark 3.4 (for the test PDE (36)), we observe that large steps of sizes up to h ≃ 4 • 10 -2 are taken by the PIROCK method, whereas the PRKC method uses (nearly) constant stepsizes h ≃ 1.6 • 10 -3 because of its stability constraint hρ A ≤ 1.7273 (with ρ A ≃ 1040). For the tolerances tol = 10 -r , r = 2, . . . , 5, and the advection parameters µ = 1, 10 -1 , we compare in Figure 14 the integrators RKC, ROCK2, PRKC, PIROCK. For the advection parameter µ = 10 Figure 13: Problem (54) with advection parameter µ = 1. Left pictures: solutions with PIROCK at x = (0, 0) for the tolerances tol = 10 -1 , 10 -2 as functions of time (the steps are plotted in bullet points). Right picture: corresponding variable stepsizes for PIROCK and PRKC (Initial stepsize is h = 10 -5 ).

both PRKC and PIROCK have a significantly reduced number of evaluations of F A , corresponding to the advection and reaction, compared to RKC and ROCK2 (see the bottom right picture in Fig. 14). When the advection becomes large (µ = 1), we observe that this is still the case for PIROCK but not for PRKC where the advection ellipse half-height is limited to ≃ √ 3 in the stability domain of the F D -F A coupling (notice that in this case, PRKC yields the same results for all considered tolerances and we obtain a single point in Fig. 14). We took the initial step size h = 10 -5 for all integrators, and used for the reference solution DOPRI5 [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF] with tolerance tol = 10 -8 . Notice that the integrator IRKC is less efficient than RKC for this problem as the reaction terms are non-stiff. Thus, we did not include the results with this integrator in our comparisons.

A 1D integro-differential equation

We consider the 1D integro-differential problem from [START_REF] Vasudeva Murthy | Solving parabolic integro-differential equations by an explicit integration method[END_REF] and considered in [START_REF] Zbinden | Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems[END_REF]Sect. 6.3], with an integral source term that is computationally expensive. It models the temperature profile of air near the ground,

∂u ∂t (x, t) = ∂u ∂x 2 (x, t) -σ 1 0 u(s, t) 4 (1 + |x -s|) 2 ds, 0 ≤ x, t ≤ 1, u(x, 0) = cos 2 (πx/2), 0 ≤ x ≤ 1, u(0, t) = 1 - √ t/2, 0 ≤ t ≤ 1, u x (1, t) = 0, 0 ≤ t ≤ 1, (55)
where σ = 10 -2 . We discretize the space with a uniform grid with n = 100 mesh intervals. The integral is approximated by the second order trapeziodal rule this experiments is very short and we rather present our comparisons in terms of functions evaluations of the different components of the problem. In Table 2, we compare the integrators RKC, PRKC, ROCK2, PIROCK. For the integrators PRKC and PIROCK, we consider the partitioning where F D corresponds to the discrete Laplacian, and F A corresponds to the discretized integral source term. For RKC and ROCK2, we consider the single vector field F = F D +F A . We take the initial step size h = 10 -3 and we consider the tolerances tol = 10 -r , r = 1, 2, 3, 4, respectively. We report the L 2 and L ∞ errors compared to a reference solution computed with DOPRI5 with tolerance tol = 10 -8 . Notice that for large n, the evaluation of F D has complexity O(n), while the evaluation of F A has complexity O(n 2 ) and thus dominates the cost. Compared to RKC and ROCK2, we observe that for all considered tolerances, the partitioning F D -F A in PRKC and PIROCK permits to reduce by one order of magnitude the number of function evaluations of the most costly integral term F A , as illustrated in Table 2. We again emphasize that IRKC presents no advantage compared to RKC for this problem.

For this problem we thus have to open all the blades of the "swiss-knife". For the diffusion and advection parameters, we take ν = 0.1, µ = 0.1, U = (-0.5, 1) T , V = (0.4, 0.7) T . We also consider a stiff reaction with parameters A = 1.3, B = 10 7 , and with stiff noise parameters σ 11 = 3, σ 12 = 4.4 • 10 3 , σ 21 = 0.5, σ 22 = 1. Notice that since -B + σ 2 21 /2 < 0, the reactionnoise system can be shown to be mean-square stable. We also consider an inhomogenity defined as f (x) = 5 if (x 1 -0.3) 2 + (x 2 -0.6) 2 ≤ 0.3 2 , and f (x) = 0 else. We consider a space discretization with two 200 × 200 meshes and consider the constant time step size h = 10 -2 on the time interval (0, 1). The number of stages used at each step to treat the diffusion is s max = 28. We plot in Figure 15 one realisation of the problem (56). In picture 15(a), we plot the solutions u(x, t), v(x, t) as a function of time t for x = (0, 0) fixed, while in picture 15(b), we plot the solution v(x, t) at final time t = 1 as a function of the spatial variable x = (x 1 , x 2 ). It can be seen that the solution oscillates stochastically in time, while it remains smooth in space. Notice that for the standard Euler-Maruyama method, the step size restriction for mean-square stability can be estimated as h ≤ 0.64 • 10 -8 , which makes this method of no practical use for this problem.

Source code for PIROCK and some examples are available at http://anmc.epfl.ch Any feedback for this code is welcome.

Figure 1 :

 1 Figure 1: Complex stability domains (gray regions) of a standard second order explicit method (Heun method, left picture) versus an explicit stabilized method (ROCK2 with damping r = 0.95 (α = 1) for s = 13. Right picture: the dashed line indicate the stability domain boundary of the embedded method used for step size control (see Section 3.3)).

Figure 2 :

 2 Figure 2: Stability polynomials for the ROCK2 method with parameters s = 13, α = 1.2.The oscillating polynomial is R s,α (z), the dotted lines correspond to the stability polynomials of the internal stages P j (αz), j = 1, . . . , 11. The stability functions P s-1 (αz), P s (αz) of the two supplementary stages K s-1 , K s are depicted in solid (red) lines oscillating with a small amplitude.

c 1 a

 1 11 . . . a 1m . . . . . . . . . c m a m1 . . . a mm b 1 . . . b m

Figure 3 :

 3 Figure 3: Complex stability regions for the F A and the F R when applied to the linear test problem The dashed lines indicate the boundaries of the error estimator stability domains (see Section 3.3).

F 10 FFigure 4 :

 104 Figure 4: Diffusion-advection (F D -F A ) coupling in PIROCK with s = 13. Stability domain (dark gray) in the pq-axis. The dotted lines indicate the largest inscribed ellipses. Left picture: damping (32). Right picture: damping (33) where α ≃ 1.363. Diffusion-advection coupling For the F D -F A coupling, we consider the test problem ẏ(t) = λy(t) + iµy(t)

Figure 5 :

 5 Figure 5: Diffusion-advection (F D -F A ) coupling in PIROCK. Half-height and width of the largest ellipse inscribed in the stability domain as a function of the stage parameter s. Solid lines: damping (32). Dashed lines: damping (33).

FFFigure 6 :

 6 Figure 6: Diffusion-reaction (F D -F R ) coupling with s = 13. Stability domain in the pr-axis. Left picture: damping (32). Right picture: damping (33).

Figure 7 :

 7 Figure 7: Case of a non-symmetric differential operator with damping (33) (F D S -F D A coupling (34)). Left picture: angle θ s of the inscribed sector as a function of s. Right picture: stability domain in the pq-axis for s = 13.

FFigure 8 :Figure 9 :

 89 Figure 8: Mean-square stability domains (gray regions) (44). Left picture: diffusion-noise coupling (stability function (47) with (32) and s = 13). Right picture: reaction-noise coupling which is mean-square A-stable (stability function (48)).

Figure 10 :

 10 Figure 10: Comparison of IRKC and PIROCK for the 2D Brusselator problem with a stiff reaction (51). CPU time (seconds) versus L 2 and L ∞ errors for the tolerances tol = 10 -r , r = 1, . . . 5.

Figure 11 :

 11 Figure 11: Comparison of IRKC and PIROCK for the combustion problem with a nonlinear diffusion and a stiff reaction. CPU time (seconds) versus L 2 and L ∞ errors for the tolerances tol = 10 -r/2 , r = 2 . . . 6 (IRKC) and tol = 10 -r/2 , r = 5 . . . 9 (PIROCK).

  Material temperature T at t = 3. (c) Radiation energy E at t = 3. (d) Material temperature T at t = 3.

Figure 12 :

 12 Figure 12: Combusion problem with a nonlinear diffusion and a stiff reaction. Reference solution at t = 3. Spatial discretization: 100 × 100.

4 F 1 diffusion ρ = 10 - 2 L 2 FFigure 14 :

 41102214 Figure 14: Comparison of RKC (dashed-dotted lines), ROCK2 (dotted lines), PRKC (dashed lines), PIROCK (solid lines) for the 2D Brusselator problem with advection (54). Number of function evaluations versus L 2 errors for the tolerances tol = 10 -r , r = 2, . . . 5.

  Solution at the spatial point x = (0, 0). (b) Solution v(x, 1) at final time t = 1.

Figure 15 :

 15 Figure 15: Non-symmetric diffusion-advection-reaction-noise problem (56). Space discretization: two 200 × 200 meshes. Constant step size h = 10 -2 .4.5 A 2D brusselator with non-symmetric diffusion, advection, a highly stiff reaction, and stiff Itô stochastic noiseTo illustrate the versatility of the proposed PIROCK integrator, we consider the Brusselator problem with simultaneously all the difficulties of a non-symmetric diffusion operator, a stiff reaction, advection, and a two-dimensional stiff Itô stochastic noise, defined as∂u ∂t = ν∆u + ν/2 ∆v + µU • ∇u + A + u 2 v -(B + 1)u + (σ 11 + σ 12 u) Ẇ1 , ∂v ∂t = -ν/2 ∆u + ν∆v + µV • ∇v + f + Buu 2 v + (σ 21 + σ 22 uv) Ẇ2 .(56)

Table 1 :

 1 Comparison of IRKC and PIROCK for the 2D Brusselator problem with a stiff reaction (51).

	method	tol	F D evals	F R evals	F R Jacobian evals	steps	smax	L 2 /L ∞ errors	CPU
	IRKC PIROCK 10 -1 10 -1 IRKC 10 -2 PIROCK 10 -2 IRKC 10 -3 PIROCK 10 -3 IRKC 10 -4 PIROCK 10 -4 IRKC 10 -5 PIROCK 10 -5	1045 749 1175 912 1952 1400 3865 2845 8477 5889	2134 55 2443 75 4295 160 10402 913 23804 2363	1044 10 1174 14 1951 31 3864 159 8476 456	27(8) 10(0) 56(6) 14(0) 152(5) 34(3) 462(4) 161(6) 1450(4) 458(6)	206 182 171 150 152 137 106 114 71 74	6.2 • 10 -2 4.2 • 10 -2 3.8 • 10 -3 5.4 • 10 -3 8.0 • 10 -4 9.3 • 10 -4 6.6 • 10 -5 1.6 • 10 -4 7.5 • 10 -6 1.7 • 10 -5	2.8 • 10 -1 1.3 • 10 -1 1.2 • 10 -2 1.8 • 10 -2 2.2 • 10 -3 2.6 • 10 -3 1.9 • 10 -4 4.5 • 10 -4 2.1 • 10 -5 4.9 • 10 -5	29.3 0.8 33.6 1.0 57.5 1.7 129.4 5.0 282.4 11.8

Table 2 :

 2 Comparison of RKC,ROCK,PRKC,PIROCK for the 1D integro-differential problem (55).

In the literature, integrators of the type (17)-(18) are sometimes called additive RK methods, while partitioned RK methods refer to integrators for systems of the form ṗ = f (p, q, . . . ), q = g(p, q, . . .),. . . Notice however that the order conditions are equivalent[START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF].

0 f (s, u(s))ds ≃ n -1 ( n-1 i=1 f (x i , u i ) + f (x 0 , u 0 )/2 + f (x n , u n )/2). We use for the Laplacian the standard second order central finite difference formula. The CPU time for all the integrators used in

Acknowledgements. The authors are grateful to Ernst Hairer and Christophe Zbinden for helpful comments on an earlier version of the manuscript, and for providing the FORTRAN codes and drivers of IRKC and PRKC. The research of A. A. is partially supported by the Swiss National Foundation Grant 200021 140692.

Variable step size control

We discuss here variable step size control related to the application of PIROCK to the ODE [START_REF] Abdulle | Chebyshev methods based on orthogonal polynomials[END_REF]. Adaptive step size for stochastic problems is still to be developed for the PIROCK method and will not be discussed here. To obtain an a posteriori error depending on the computed solution, we use the idea of embedded methods [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]. Considering [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF], the idea is to define another set of coefficients b i , bi , bi say b e,i , be,i , be,i such that the second order non-partitioned order conditions [START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF] are not fulfilled, but only the first order conditions (except for the advection method where we consider a second order embedded method). We then consider three embedded methods y e,D , y e,A , y e,R where b i , bi , bi are changed respectively to b e,i , be,i , be,i . We obtain for the error estimators err D = y 1y e,D , err A = y 1y e,A , err R = y 1y e,R the values

The step size adjustments for adaptive integration are based on the standard strategy (50). Since the order is p = 2 for the F D and F R methods and p = 3 for the F A method, the step size is selected using

where ζ ∈ (0.1, 0.8) is a safety factor and p = 2 is the order of ROCK2 and using the error estimator err = max( err D , err A 2/3 , err R ).

Given the step size h (the initial step size or the adapted step size from (50)), the following selection procedure of PIROCK parameters is implemented in the code, based on the sizes d s , a s of the largest ellipses (39) and (40) inscribed in the stability domain of the F D -F A coupling:

1. Consider the parameters (32). Select the stage number s such that d s ≃ hρ D (with d s ≃ 0.81 • s 2 ), where ρ D is an estimate of the spectral radius of the Jacobian ∂F D /∂y.

2.

If F A = 0 and hρ A > a s where ρ A is an estimate of the spectral radius of ∂F A /∂y, then consider the damping (33) and select s such that

The choice of estimators (49) is motivated as follows. For the F D -method we use the embedded method K * s in (31) that is identical to the one of the standard ROCK2 method. For the estimator of the F A -method, we ask that the embedded method has the stability polynomial 1 + z + z 2 /2 + z 3 /5. It has order 2 with a similar stability domain as the F A method (see Figure 3). For the estimator of the F R -method, we impose that the embedded method has order one and is A-stable. Since the embedded method is not L-stable, we use Shampine's idea [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF]IV.8] and damp it by applying J -1 R , where the LU decomposition of the matrix J R in (30) (recall that this decomposition is already available, see (30)).


Anthony Martin

Florian Kaiser

Aline Vernier

Alexios Beveratos

Valerio Scarani

Sébastien Tanzilli
email: sebastien.tanzilli@unice.fr

Cross time-bin photonic entanglement for quantum key distribution

Keywords: numbers: 03.67.Bg, 03.67.Dd, 03.67.Hk, 42.50.Dv, 42.65.Lm, 42.65.Wi

Introduction. -Entanglement stands as an essential resource for quantum key distribution (QKD) [1,2] and quantum communication protocols [3], such as teleportation [4], entanglement swapping [5], relays [6], and repeaters [7]. Long-distance distribution of entangled photons has been demonstrated using various observables, namely polarization [8,9], energy-time [10], as well as time-bin [11]. Fully operational "out of the lab" quantum communication imposes stringent constraints on the source in terms of compactness, stability, and brightness, as well as quality of entanglement. Promising results were obtained using entangled photon sources built around nonlinear integrated waveguide generators [12]. In this framework, operating wavelengths around 1550 nm (telecom C-band) offer the possibility of long-distance entanglement distribution, taking advantage of low loss standard optical fibers, as well as highperformance telecommunications components [13][14][15][16][17].

In this paper, we demonstrate a novel time-bin entangled photon-pair source encompassing all the necessary criteria for real-world quantum applications. Unlike previous time-bin source realizations [10,14,18], our approach exploits cross-polarized paired photons generated from a type-II periodically poled lithium niobate waveguide (PPLN/W) [17]. A polarization to time-bin observable transcriber then enables preparing the photons in the time-bin Bell state

|Ψ - A,B = 1 √ 2 (|0 A , 1 B -|1 A , 0 B )
, where 0 and 1 denote short and long time-bins, and A and B the parties Alice and Bob, respectively. To our knowledge, the analysis of such a state has never been reported. This scheme, implemented in the continuous wave (CW) regime, is completely insensitive to phase fluctuations in the time-bin preparation stage (the transcriber), and allows 75% of the detected pairs to be exploited, as opposed to 50% for standard pulsed regime time-bin schemes [10,14,18]. By simply trusting the randomness of beam-splitters, it provides correlations that cannot be achieved with classical means even in a passive implementation.

The setup. -The setup of the source is shown in [17]. Note that various phase-matching conditions obtained in periodically poled materials can be found summarized in [12]. The brightness of the PPLN/W was measured to be ∼ 2 • 10 4 generated pairs per second, mW of pump power, and GHz of emission bandwidth, coupled into a single mode fiber [17], while multiple-pair emission probability is kept below 1% per time detection window. We then select only pairs of wavelength degenerate photons so as to prevent any polarization discernibility as a function of the wavelength [17]. To this end, a fiber Bragg grating (FBG) filter is used to reduce the natural SPDC emission bandwidth from ∼850 pm down to 200 pm, therefore avoiding undesirable spectral responses associated with alternative phase-matchings. This results in a single-photon coherence time of ∼17 ps (↔ coherence length of ∼5 mm).

After the filtering stage, the two photons are sent to a transcriber, arranged as a Michelson interferometer, so as to introduce a time delay between the |H and |V polarization modes. A fiber polarization controller (PC 1 ) allows optimizing the separation of the two polarization components at a fiber-pigtailed polarizing beam-splitter (PBS 1 ). In each arm, we use a fiber Faraday mirror (FM), rotating the associated polarization state by 90 • after a round trip. This ensures that both photons leave the transcriber through the output port of PBS 1 . For this experiment, the transcriber's optical path length difference is set to ∼60 cm (↔ ∼2 ns), which is much greater than both the coherence time of the single-photons and timing jitter of the employed detectors (∼0.4 ns). This prevents the photons to overlap temporally and the state at the output of the transcriber therefore reads |H, 0 |V, 1 , where 0 and 1 denote short and long time-bins, respectively. Next, the photons are sent to an additional polarizing beam-splitter (PBS 2 ) which is oriented at 45 • with respect to the {H; V } basis in which the photons are created. This way, the polarization modes are no longer associated with the timebins, therefore reducing the two-photon state to

|ψ A,B = 1 2 (|0 A |1 A + |0 A |1 B -|1 A |0 B -|1 B |0 B ).
Subsequently, the maximally entangled Bell state

|Ψ - A,B = 1 √ 2 (|0 A |1 B -|1 A |0 B
) can be post-selected among all other events using a coincidence detection electronics between the two parties. This state is insensitive to the phase accumulated by the two photons in the transcriber, unlike in standard time-bin schemes where entanglement is prepared using a stabilized unbalanced interferometer placed on the path of a pulsed pump laser [13,14,18].

The entanglement analysis is performed using a Franson setup with two equally unbalanced interferometers, one at each side [19]. Each interferometer is made of a 50/50 fiber coupler and two fiber-pigtailed FMs to compensate for the birefringence within the analyzer [10,13,14,18]. As the interferometers' path length difference are adjusted to suitably match the transcriber time-bin separation, such analyzers transform their respective incoming single-photon state in the following way, |τ j → 1 √ 2 |τ j + e iφj |τ + 1 j , where τ = {0, 1} represents the considered time-bin, j Alice or Bob, and φ j the phase in Alice's or Bob's interferometer, respectively. After the complete analysis apparatus including the post-selection procedure, the two-photon state reads:

|ψ = 1 √ 8 |0 A |1 B + e i φA |1 A |1 B +e i φB |0 A |2 B + e i (φA+φB ) |1 A |2 B -|1 A |0 B -e i φA |2 A |0 B -e i φB |1 A |1 B -e i (φA+φB ) |2 A |1 B . (1)
As can be understood from Eq. 1, one needs to consider five relative arrival times between Alice and Bob detectors. The time-dependent second order intensity correlation function is measured using a free running InGaAs avalanche photodiode (APD, IDQ-210) on Alice's side which is used as a trigger for Bob's gated InGaAs APD (IDQ-201). FIG. 2 . Remarkably, the peaks T ±1 on one hand, and T 0 on the other hand, naturally define two complementary basis useful for entanglementbased QKD protocols [1,2,20,21].

To quantify the quality of our analysis interferometers, we first perform a two-photon interference experiment using a genuine energy-time Franson setup [19]. In this case, the transcriber is bypassed (see FIG. 1), such that the cross-polarized pairs of photons are directly sent to PBS 2 . The polarization state of the two photons is adjusted using PC 3 for the pairs to be deterministically separated at this PBS. The relative phase between the analysers is varied by tuning the temperature of Alice's interferometer, while that of Bob's is kept constant. We achieve net and raw visibilities of 98±2% and 96±2%, respectively (curve not represented). The net visibility is obtained after subtraction of accidental coincidences due to the detectors dark-counts. We ascribe ∼1% net contrast reduction to the length mismatch between the two analysers which was measured to be 0.3±0.1 mm. Note that multiple-pair generation probabilities in this experiment are very small (<0.1%).

Eventually, we test the cross time-bin state |Ψ -

A,B

prepared with the transcriber. For the three coincidence peaks T -1 , T +1 , and T 0 , two-photon interference is also recorded by tuning the phase φ A , while keeping φ B ≈ 0.

As shown in FIG. 3 (a), the three coincidence rates are modulated as a function of the global phase (φ A +φ B ) and lead to net (raw) visibilities exceeding 97±2% (95±2%).

For this test, we ascribe an additional 1% visibility reduction to the path length mismatch between the transcriber and the analysers. In addition, the visibility associated with events detected in the T 0 time-bin is strongly dependent on the polarization state adjustments in front of the second PBS where the |Ψ -state is prepared. Also, the raw visibilities could be improved by employing detectors showing better dark-count figures, such as those offered by super-conducting devices. Furthermore, to show that the three contributions follow the expected phase dependency, we change the phase of Bob's interferometer to φ B = π/2. As expected, FIG. 3 (b) confirms that the phase relation of the interference patterns between T ±1 and T 0 is shifted by π. Note that the analysers are only temperature stabilized which explains the drifts between the patterns shown in FIG. 3 (a) and (b). This could be avoided by employing active interferometer stabilization schemes enabling accurate phase control and setting, as was demonstrated using a frequency stabilized laser combined to a phase measurement [22,23].

Theory. -For the characterization of the setup, we have plotted interference patterns, out of which one may check the violation of some Bell inequality [20] under some assumptions proper to time-bin implementations [24,25]. However, we want to operate this setup for QKD in a passive implementation, that is φ A and φ B are going to be fixed. The detections of any passive setup can always be reproduced by a local variable model if one does not make additional assumptions: indeed here, if Eve can choose the time of the detection, she can obviously fulfill all the requirements. So we make the following assumption on Alice's and Bob's measurements separately: if some signal comes in the analyzer at time τ , the detection can happen either at τ or at τ + 1 (up to the propagation time in the interferometer) and this time is not for Eve to control. In other words, we trust the randomness introduced by the beam-splitter. Under this assumption, we are going to provide a semi-deviceindependent criterion for quantumness, analog to a Bell inequality [START_REF]The specific local variable models developed[END_REF].

Consider Eve sending some signal in Alice's lab at time τ A and some signal in Bob's lab at times τ B . As per our assumption, detection at Alice's side can happen at τ A or at τ A + 1. Apart from this uncertainty in the time of firing, we assume the most favorable situation for Eve, namely that for both times she knows with certainty which detector will fire; we also assume that only one detector will fire (notice that there may be no state of light such that this can actually happen: we give more power to Eve by assuming only her impossibility of choosing the time of the firing). We do similar assumptions for Bob. Thus, Eve knows the four possible outcomes a(τ A ),

a(τ A + 1), b(τ B ), b(τ B + 1), with (a, b) = ±1.
For definiteness, we assume φ A = φ B = 0: for these values, ideally, the quantum source produces identical outcomes at T 0 and T ±1 , and completely uncorrelated outcomes at T ±2 . The QKD protocol we want to propose uses the detections at T 0 and T ±1 as the raw key, while the lateral peaks could be used to detect the presence of an eavesdropper. We need to see how close Eve can come to simulating the ideal quantum source (we assume that the beam-splitters in the Michelson interferometer have 50/50 coupling ratio, although this additional assumption can easily be relaxed):

(a) τ A = τ B = τ : in this case, T 0 happens with probability 1 2 , and T ±1 with probability 1 4 each. Eve can ensure that Alice and Bob get the same outcomes by fixing a(τ ) = a(τ + 1) = b(τ ) = b(τ + 1). But, if Eve does only this, there will never be any event T ±2 .

(b) τ A = τ B -1 = τ : in this case, T +1 happens with probability 1 2 , and T 0,+2 with probability 1 4 for each. Now, Eve cannot guarantee all the conditions: Eve can set a(τ ) = b(τ + 1) = a(τ + 1) without loss of generality, but she has to choose whether b(τ + 2) is also equal to those, so as to get the case (τ + 1, τ + 2) right; or random, to get the case (τ, τ + 2) right.

(c) τ A = τ B -2 = τ : in this case, T +2 happens with probability 1 2 , and T +1,+3 with probability 1 4 for each. Eve can set a(τ + 1) = b(τ + 2), while a(τ ) and b(τ + 3) are uncorrelated with anything. This way, she gets all the expected detections, but she may generate an unexpected event at T 3 .

As we see, one has to check for high correlations in T 0,±1 , for the presence of uncorrelated coincidences in T ±2 , and for the absence of coincidences in T ±3 . There are several ways to do that; among the simplest, one can form a single linear function like

S = x=0,±1 P (T x )E Tx + 2 x=±2 P (T x )(1 -E Tx ) - x=±3 P (T x ) (2)
where P (T x ) is the probability of a coincidence T x and E is the correlation coefficient. All the classical strategies above reach at most S = 1, while perfect experimental results would lead to S = 1.25. Our experimental raw data from FIG. 3 (a) lead to S exp raw = 1.20 ± 0.02, corresponding to 10 standard deviations above the classical boundary. To infer this value, we take into account the above mentioned visibilities for the coincidence peaks T ±1 and T 0 , as well as an average visibility of 0±2% for the coincidence peaks T ±2 (curves not shown). This result can be improved by employing lower dark-count detectors and by adjusting actively the path length difference of the analysis interferometers to match that of the transcriber. Note that a similar theoretical/experimental analysis with φ A = -π 2 and φ B = π 2 (case of FIG. 3 (b)) leads to the same results and conclusion. In this case, Alice and Bob expect the observation of coincidences in T 0 and simultaneously anti-coincidences in T ±1 . As Eve cannot know in which coincidence peak the photons are detected, she cannot mimic the correct coincidence/anticoincidence behavior.

Conclusion. -We have demonstrated, for the first time, a scheme that allows analysing the cross timebin entangled state |Ψ -. From the experimental side, the simplicity of the scheme associated with the phase-insensitivity of the prepared state to the transcriber fluctuations, make this source a promising candidate for time-bin entanglement based long-distance QKD. In this perspective, such a novel strategy offers naturally the two necessary complementary basis. This makes it possible to operate entanglement based QKD protocols [21,[START_REF] Ling | [END_REF] in a full passive fashion, requiring only one stabilized analyser plus two detectors at each user's location, and allowing the exploitation of 75% of the detected pairs of photons [28]. Moreover, the source can be operated in both CW and pulsed regimes, without any further experimental modification but the laser. The latter regime is particularly interesting when synchronization between different locations is required, as is the case in quantum relay schemes [6]. Using exclusively standard fiber telecom components allowed us reaching an overall loss figure of ∼5 dB from the output of the PPLN/W to the detectors. A simple theoretical model trusting only the randomness of the beam-splitters also confirms the quantum nature of the source, and therefore its potential for passive implementations of entanglement-based QKD.

FIG. 1 .

 1 FIG.1. Experimental setup for generating and analysing time-bin entangled states starting with cross-polarized pairedphotons. Two different Bell states can be prepared: on one hand the |Φ + state, using the bypass, and on the other hand the |Ψ -state, using the transcriber. The analysis is done using two equally unbalanced Mach-Zehnder interferometers in the Franson configuration[19]. APD: avalanche photodiode; FM: Faraday mirror.

FIG. 2 .

 2 FIG. 2. Coincidence histogram for the |Ψ -state showing the five expected peaks (see text for details).

2 .

 2 shows the coincidence histogram as a function of the time delay between the photon detection on Alice's and Bob's sides. The central peak, labelled as T 0 , contains the two contributions associated with photon-pairs in the state |1 A |1 B . The indistinguishablitiy of these two contributions results in usual timebin entanglement for which the coincidence rate follows R T0 c ∝ cos 2 φA-φB In our CW experimental setup, the pump coherence time (τ p c ≃ 400 ns) is significantly larger than the time-bin separation (2 ns), such that the emission times of the paired photons remain unknown within τ p c . As a result, the paths |0 A |1 B and |1 A |2 B , as well as |1 A |0 B and |2 A |1 B , are also indistinguishable. These contributions are respectively labelled as T -1 and T +1 in FIG. 2. In these cases, the coincidence rate is given by R T±1 c ∝ cos 2 φA+φB 2

FIG. 3 .

 3 FIG. 3. Bell tests performed for the T±1 and T0 conincidence peaks of FIG. 2 as a function of φA + φB, with φA continuously tuned and φB set to 0 (a) and to π/2 (b).

  FIG.1. Pairs of cross-polarized, i.e. horizontally (|H ) and vertically (|V ), photons are generated at the wavelength of 1540 nm from spontaneous parametric downconversion (SPDC) of a 770 nm CW, |H polarized, 2.5 mW, pump laser in a type-II PPLN waveguide. This degenerate phase-matching condition is reached in our case for a 3.6 cm long, 9.0 µm periodically poled, sample, heated at the temperature of 110 • C

We thank W. Sohler and H. Herrmann for providing the PPLN waveguide and for fruitful discussions. Financial support is acknowledged from the CNRS, the Conseil Régional PACA, the French Ministère de l'Enseignement Supérieur et de la Recherche (MESR), the European program ERA-SPOT "WASPS", the iXCore Research Foundation, as well as from the Ministry of Education and the National Research Fund of Singapore.

Note added: during the preparation of this manuscript, we became aware of a preprint paper in which cross timebin states are employed for studying, in a convincing manner, 'time-bin entangled photon holes' [29].