Qinmeng Zou

Guillaume Gbikpi-Benissan

Frédéric Magoulès
email: frederic.magoules@hotmail.com

Asynchronous Parareal Algorithm Applied to European Option Pricing

Keywords: asynchronous iterations, parareal method, parallel computing, European options, Black-Scholes equation, time-dependent problems

Asynchronous iterations arise naturally in parallel computing if one wants to solve large problems with a minimization of the idle times. This paper presents an original model of asynchronous iterations for a time-domain decomposition method, namely the parareal method. The asynchronous parareal algorithm is here applied to European option pricing, and numerical experiments performed on a parallel supercomputer, illustrate the performance and efficiency of this new method.

Introduction

An option is a contract that gives the right to buy or sell the underlying asset at a specific price within a period of time, which gives the buyer the right, but not the obligation to perform the transaction. Option pricing is a crucial target in financial decision-making. In the practical side, however, it is also one of the most complicated problems in the mathematical domain of financial engineering. The breakthrough came in 1973, when Black and Scholes proposed the original Black-Scholes equation [START_REF] Black | The pricing of options and corporate liabilities[END_REF], which is the first complete European option pricing model. With an alternative method, Merton derived the formula and generalized it in important ways [START_REF] Merton | Theory of rational option pricing[END_REF][START_REF] Merton | On the pricing of corporate debt: the risk structure of interest rates[END_REF]. The Black-Scholes model has a huge influence to the financial market and drives an unexpected prosperity in the trading of derivatives.

With the advent of parallel computers, many complicated problems can be solved efficiently by the improving processing power and new algorithms. An example is the asynchronous iterative scheme. The pioneering research of such approach dates back to 1969 by Chazan and Miranker [START_REF] Chazan | Chaotic relaxation[END_REF] that arises in the face of communication delays. Their works have been extended to the nonlinear cases in [START_REF] Miellou | Algorithmes de relaxation chaotique à retards[END_REF][START_REF] Baudet | Asynchronous iterative methods for multiprocessors[END_REF][START_REF] Tarazi | Some convergence results for asynchronous algorithms[END_REF] based on the notion of contraction. A well known systematic formalizations is contained in the book of Bertsekas and Tsitsiklis [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]. Furthermore, two-stage theories were proposed in the 1990s by Frommer and Szyld [START_REF] Frommer | Asynchronous two-stage iterative methods[END_REF][START_REF] Frommer | Asynchronous iterations with flexible communication for linear systems[END_REF]. We refer the reader to [START_REF] Frommer | On asynchronous iterations[END_REF][START_REF] Kollias | Asynchronous iterative algorithms[END_REF] as surveys on asynchronous iterative methods. Recently, asynchronous iterations have been extended to efficient domain decomposition methods in space [START_REF] Toselli | Domain Decomposition Methods: Algorithms and Theory[END_REF], including sub-structuring methods [START_REF] Magoulès | Asynchronous iterative sub-structuring methods[END_REF], optimized Schwarz method [START_REF] Magoulès | Asynchronous optimized Schwarz methods with and without overlap[END_REF]. Another issue is the extension to domain decomposition methods in time. In this paper, we are interested in the asynchronous parareal scheme, which will be explored in Section 2. In Section 3 we present the transformation from the original Black-Scholes equation to an appropriate form which applies to our asynchronous solver. Section 4 is devoted to the experimental results. Finally, concluding remarks and future directions are presented in Section 5.

Asynchronous Parareal Algorithm

Classical Parareal Scheme

The parareal algorithm is a method to solve the time-dependent partial differential equations in parallel, first presented by Lions, Maday and Turinici in 2001 [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF]. The ideas of parareal are relevant to both multiple shooting method and predictor-corrector scheme, see the historical remarks in [START_REF] Gander | Analysis of the parareal time-parallel time-integration method[END_REF] whereby a comprehensive analysis of parareal has been mathematically studied.

Given a second-order linear elliptic operator L, consider the following time-dependent problem

∂u ∂t + Lu = f, t ∈ [0, T ], u = u 0 , t = 0, (1)
where the boundary conditions are unnecessary to be mentioned. By introducing a function λ n , the above mathematical description can be decomposed into N sequential problems, according to the rule T n = n∆T with n = 0, . . . , N . We now rewrite our problem (1) as

∂un ∂t + Lu n = f n , t ∈ [T n , T n+1 ], u n = λ n , t = T n , (2)
where n = 0, . . . , N -1, together with the condition λ n+1 = lim →0 u n (T n+1 -). The parareal scheme is established by two propagators. Let δt be a fine time step corresponding to the problem [START_REF] Ahamed | Iterative methods for sparse linear systems on graphics processing unit[END_REF]. Suppose that the time-dependent problem is approximated by some appropriate classical discretization schemes. Let F be a precise integration method based on δt, often called a fine propagator. Similarly, let G be a coarse propagator with respect to the broad time step ∆T . It is seen that the discretization schemes used for such two levels might be different. Finally, parareal replaces the accurate sequential method by an iteration

λ k+1 n+1 = G(λ k+1 n ) + F (λ k n ) -G(λ k n ),
which is a predictor-corrector iterative scheme with n = 0, . . . , N -1 and λ 0 0 = u 0 . Additionally, we are obliged to respect two more conditions λ 0 n+1 = G(λ 0 n ) and λ k+1 0 = λ k 0 throughout the iterative process. We note that the time-dependent problems might benefit from such scheme only if the number of iterations is less than N . Following the mathematical model, we arrive at the iterative process below

Algorithm 1. (Parareal Iterative Algorithm) λ 0 0 = u 0 for n = 0 to N -1 do λ 0 n+1 = G(λ 0 n ) end for k = 0 while not convergence do Solve F (λ k n ) in parallel for n = 0, . . . , N -1 λ k+1 0 = λ k 0 for n = 0 to N -1 do Solve G(λ k+1 n ) λ k+1 n+1 = G(λ k+1 n ) + F (λ k n ) -G(λ k n ) end for k = k + 1 end while
Notice that F can be evaluated independently using parallel computers. Thus, parareal is often preferable than the serial setting for the time-dependent problems.

Asynchronous Parareal Scheme

To derive the corresponding asynchronous scheme, now we present the asynchronous iterative model. The fundamental asynchronous model (see, e.g., [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]) is

x k+1 i = f i (x µ i 1 (k) 1 , . . . , x µ i p (k) p ), i ∈ P k , x k i , i / ∈ P k , (3)
where p is the number of processor, P k ⊆ {1, . . . , p}, P k = ∅ and µ i j (k) are integers satisfying 0 ≤ µ i j (k) ≤ k. In addition, the model ( 6) is subject to the following conditions

lim k→+∞ µ i j (k) = +∞, i, j ∈ {1, . . . , p}, {k ∈ N | i ∈ P k } = +∞, i ∈ {1, . . . , p}, (4)
which indicates that each processor reads eventually the latest information for each component and no processor stops updating. Unfortunately, this model fails to be applied to the two-stage problems. Therefore, we explore an alternative method which is often called asynchronous iterations with flexible communication, see [START_REF] Frommer | Asynchronous iterations with flexible communication for linear systems[END_REF]. To the end we will consider the model

x k+1 i = f i (x µ i 1 (k) 1 , . . . , x µ i p (k) p , x ρ i (k) ), i ∈ P k , x k i , i / ∈ P k , with x ρ i (k) = [x ρ i 1 (k) 1 , . . . , x ρ i p (k) p ],
where ρ i j (k) are integers satisfying 0 ≤ ρ i j (k) ≤ k and similarly we have

lim k→+∞ ρ i j (k) = +∞, i, j ∈ {1, . . . , p}. (5)
In the context of the asynchronous two-stage iterative model, we can eventually derive the asynchronous parareal method. Consider the following iteration

λ k+1 n+1 = G(λ µn(k) n ) + F (λ ρn(k) n ) -G(λ ρn(k) n ), n ∈ P k , λ k n+1 , n / ∈ P k ,
where λ k+1 0 = λ k 0 , P k ⊆ {1, . . . , p} and P k = ∅, with the similar conditions

0 ≤ µ i j (k) ≤ k + 1, 0 ≤ ρ i j (k) ≤ k.
We supposed that the essential assumptions ( 7) and ( 8) are satisfied.

Asynchronous Parareal Scheme

To derive the corresponding asynchronous scheme, now we present the asynchronous iterative model. The fundamental asynchronous model (see, e.g., [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]) is

x k+1 i = f i (x µ i 1 (k) 1 , . . . , x µ i p (k) p ), i ∈ P k , x k i , i / ∈ P k , ( 6
)
where p is the number of processor, P k ⊆ {1, . . . , p}, P k = ∅ and µ i j (k) are integers satisfying 0 ≤ µ i j (k) ≤ k. In addition, the model ( 6) is subject to the following conditions

lim k→+∞ µ i j (k) = +∞, i, j ∈ {1, . . . , p}, {k ∈ N | i ∈ P k } = +∞, i ∈ {1, . . . , p}, (7)
which indicates that each processor reads eventually the latest information for each component and no processor stops updating. Unfortunately, this model fails to be applied to the two-stage problems. Therefore, we explore an alternative method which is often called asynchronous iterations with flexible communication, see [START_REF] Frommer | Asynchronous iterations with flexible communication for linear systems[END_REF]. To the end we will consider the model

x k+1 i = f i (x µ i 1 (k) 1 , . . . , x µ i p (k) p , x ρ i (k) ), i ∈ P k , x k i , i / ∈ P k , with x ρ i (k) = [x ρ i 1 (k) 1 , . . . , x ρ i p (k) p ],
where ρ i j (k) are integers satisfying 0 ≤ ρ i j (k) ≤ k and similarly we have

lim k→+∞ ρ i j (k) = +∞, i, j ∈ {1, . . . , p}. (8)
In the context of the asynchronous two-stage iterative model, we can eventually derive the asynchronous parareal method. Consider the following iteration

λ k+1 n+1 = G(λ µn(k) n ) + F (λ ρn(k) n ) -G(λ ρn(k) n ), n ∈ P k , λ k n+1 , n / ∈ P k ,
where λ k+1 0 = λ k 0 , P k ⊆ {1, . . . , p} and P k = ∅, with the similar conditions

0 ≤ µ i j (k) ≤ k + 1, 0 ≤ ρ i j (k) ≤ k.
We supposed that the essential assumptions ( 7) and ( 8) are satisfied.

Formalization of the Option Pricing Model

We now turn to the study of the option pricing problem that is expressed by appropriate evolution models. Accordingly, they apply to the above described scenario. To the end we will concentrate on the European call option, which can be exercised only on the expiration date, using the Black-Scholes model as following

∂V ∂t + rS ∂V ∂S + 1 2 σ 2 S 2 ∂ 2 V ∂S 2 = rV, (9)
where V is the option price as a function of underlying asset price S and time t. The parameters are volatility σ and risk-free interest rate r. Furthermore, the Black-Scholes equation is backward parabolic satisfying

     V (S, T ) = max(S -E, 0), S ≥ 0, V (0, t) = 0, t ≥ 0, V (S, t) ∼ S -Ee -r(T -t) as S → +∞, t ≥ 0,
where the first formula is a final condition, T is the time to maturity and E is the exercise price. We are particularly interested in the solution of its discretization. However, the irregular form poses some difficulties for the solution. Therefore, we now try to transform the formula (9) to a simplest case.

A well known special case is the heat equation, which leads to a really simple parabola. Thus, the transformation from the original to the latter can simplify the solving process. Therefore, a change of variables is

S = Ee x , t = T - 2τ σ 2 , V = Ev(x, τ ),
substituting into the Black-Scholes equation [START_REF] Frommer | Asynchronous iterations with flexible communication for linear systems[END_REF] gives

∂v ∂τ = (κ -1) ∂v ∂x + ∂ 2 v ∂x 2 -κv, where κ = 2r σ 2 . Setting α = 1 2 (κ -1), β = 1 2 (κ + 1), (10)
yields v = e -αx-β 2 τ u(x, τ ),
then gives the heat equation in an infinite interval

∂u ∂τ = ∂ 2 u ∂x 2 , x ∈ R, τ ∈ [0, T σ 2 2 ]. ( 11
)
heat equation is forward parabolic as opposed to the original model, with the corresponding conditions u(x, 0) = max(e βx -e αx , 0),

x ∈ R, u(x, τ ) ∼ 0 as x → ±∞, τ ∈ [0, T σ 2 2 ]
. We notice here that equation ( 11) corresponds to the initial problem (1) with f = 0 and L = -∂ 2 ∂x 2 . Hence, the parareal schemes can be used to approximate the solution. We assume that the equation ( 11) is discretized in time using a two-stage scheme. Generally, the same discretization approach as coarse solver or a method of higher order can be applied to the fine grid propagator. Nevertheless, the spatial discretization is not straightforward enough. We give the precise boundaries as following instead of limits to infinity

x -= min(x 0 , 0) -log(4), x + = max(x 0 , 0) + log [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF],

such that u(x -, τ ) = 0, τ ∈ [0, T σ 2 2 ], u(x + , τ ) = e βx + +β 2 τ -e αx + +α 2 τ , τ ∈ [0, T σ 2 2 ]
, where α and β are defined by [START_REF] Frommer | On asynchronous iterations[END_REF].

Experimental Results

In this section we present numerical results for the above asynchronous iterative scheme in order to solve the Black-Scholes equation. Let us give the volatility σ = 0.2 and the risk-free interest rate r = 0.05, which can be considered as constants since they have little relation with the asynchronous properties. Another intuitive configuration is to set δt equals one day, while taking into account T = 1 representing one year. However, such assignment seems inflexible and yields a somewhat complicated constant. Hence, we keep δt = 0.001 fixed and vary ∆T to simulate the different remaining life of the financial instrument. Finally, the spatial discretization is implemented by the finite difference method with suitable equal sub-intervals.

The asynchronous implementation requests some advanced capabilities. We conduct our tests based on JACK (an asynchronous communication kernel) library [START_REF] Magoulès | JACK: an asynchronous communication kernel library for iterative algorithms[END_REF], which provides high-level supports for both synchronous and asynchronous iterative algorithms, developed on the top of MPI library. Then, mathematical operations and linear systems solvers are implemented by Alinea (An Advanced Linear Algebra) library [START_REF] Magoulès | Alinea: An advanced linear algebra library for massively parallel computations on graphics processing units[END_REF] The Alinea library is implemented in C++, MPI, CUDA and OpenCL. For both central processing unit and graphic processing unit devices, there are different matrix storage formats, and real and complex arithmetics in single-and double-precision. It includes several linear algebra operations [START_REF] Ahamed | Fast sparse matrix-vector multiplication on graphics processing unit for finite element analysis[END_REF] and numerous algorithms for solving linear systems such as iterative methods [START_REF] Magoulès | Fast iterative solvers for large compressed-sparse row linear systems on graphics processing unit[END_REF], [START_REF] Ahamed | Iterative methods for sparse linear systems on graphics processing unit[END_REF], [START_REF] Magoulès | Auto-tuned Krylov methods on cluster of graphics processing unit[END_REF], together with some energy consumption optimization [START_REF] Magoulès | Green computing on graphics processing units[END_REF], and domain decomposition methods in space [START_REF] Magoulès | Optimized Schwarz method without overlap for the gravitational potential equation on cluster of graphics processing unit[END_REF]. We carry out two sets of experiments on a SGI ICE X cluster connected with InfiniBand (56 Gbit/s). Each computing node consists of two Intel Xeon E5-2670 v3 2.30 GHz CPUs. Furthermore, the version of the MPI library used is SGI-MPT 2.14. The first is performed by changing the length of ∆T on which we simulate the different time to maturity. The backward Euler method is chosen for both the coarse and the fine propagator. Given the initial stoke price S = 15 and the strike price E = 20, table 1 illustrates the some average results for 16 cores and 250 sub-intervals, where V a is the approximate option prices, V e is the exact price, a is the absolute error and r is the relative error with r = a Ve . Table 2 reports the similar results with S = 25 and E = 30. The absolute errors and the relative errors are small enough, which leads to the conclusion that the asynchronous iterations converge to the exact solution of the Black-Scholes equation. Therefore, the fact above illustrates that the asynchronous parareal scheme is suitable to the option pricing problem, formalized by the Black-Scholes model. Notice that T = 16 when ∆T = 1, which is beyond the scope of common requirements, even in a more general context.

The second is a comparison of synchronous and asynchronous parareal scheme. We choose S = 25, E = 30, ∆T = 0.1 and 150 sub-intervals. Table 3 reports the average experimental results, where N is the number of processor cores. The results illustrate that the number of iterations of the synchronous version remains the same throughout the test, whereas the asynchronous scheme requires further iterations with the increase of cores. Let us mention here that the context is not distributed, which rejects the possibility of large communication delays. We notice that the asynchronous version is faster.

Conclusions

In this paper we show the convergence performance of an original asynchronous extension of the parareal scheme applied to the Black-Scholes equation, and illustrate some cases that the asynchronous method takes on superiority than its synchronous counterpart. The asynchronous iterative methods performs well even if there exists large communication delays or the computation bottlenecks in the systems, so long as all the processors keep updating using eventually the latest values. We illustrate the excellent precision and the conditional efficiency for the target approach. A theoretical and comprehensive analysis is under investigation by the authors to explain the results.

Table 1 :

 1 Results of the Asynchronous Parareal Scheme for 16 Cores, with Approximate Option Prices V a , Exact Option Prices V e , Absolute Error a and Relative Error r (δt = 0.001, S = 15, E = 20, m = 250)

	∆T	V a	V e	a	r	Time
	0.1	0.4857 0.4853 0.0004 0.0008 0.958
	0.2	1.3950 1.3947 0.0003 0.0002 2.012
	0.3	2.3145 2.3140 0.0005 0.0002 3.014
	0.4	3.1932 3.1925 0.0007 0.0002 4.091
	0.5	4.0212 4.0203 0.0009 0.0002 5.041
	0.6	4.7972 4.7961 0.0011 0.0002 6.084
	0.7	5.5226 5.5213 0.0013 0.0002 6.873
	0.8	6.1995 6.1981 0.0014 0.0002 7.925
	0.9	6.8306 6.8291 0.0015 0.0002 8.526
	1.0	7.4184 7.4169 0.0015 0.0002 9.969

Table 2 :

 2 Results of the Asynchronous Parareal Scheme for 16 Cores (δt = 0.001, S = 25,

	E = 30, m = 250)
	∆T	V a	V e	a	r	Time
	0.1	1.5204	1.5179 0.0025 0.0016 0.969
	0.2	3.3163	3.3141 0.0022 0.0007 1.972
	0.3	4.9525	4.9504 0.0021 0.0004 3.022
	0.4	6.4497	6.4476 0.0021 0.0003 3.817
	0.5	7.8264	7.8241 0.0023 0.0003 5.221
	0.6	9.0963	9.0939 0.0024 0.0003 6.108
	0.7	10.2699 10.2673 0.0026 0.0003 6.953
	0.8	11.3557 11.3530 0.0027 0.0002 7.731
	0.9	12.3611 12.3584 0.0027 0.0002 8.871
	1.0	13.2923 13.2898 0.0025 0.0002 10.396

Table 3 :

 3 Comparison of Synchronous and Asynchronous Parareal Scheme (S = 25, E = 30, δt = 0.001, ∆T = 0.1, m = 150)

		Synchronous	Asynchronous
	N	Iter. Time Iter. (min, max, mean) Time
	16	11	0.620 22 30	26	0.490
	32	11	0.781 30 47	40	0.677
	64	11	0.971 44 77	60	0.947


Qinmeng Zou

Guillaume Gbikpi-Benissan

Frédéric Magoulès
email: frederic.magoules@hotmail.com

Asynchronous Communications Library for the Parallel-in-Time Solution of Black-Scholes Equation

Keywords: asynchronous iterative algorithms, parallel computing, parareal method, time domain decomposition, option pricing

The advent of asynchronous iterative scheme gives high efficiency to numerical computations. However, it is generally difficult to handle the problems of resource management and convergence detection. This paper uses JACK2, an asynchronous communication kernel library for iterative algorithms, to implement both classical and asynchronous parareal algorithms, especially the latter. We illustrate the measures whereby one can tackle the problems above elegantly for the time-dependent case. Finally, experiments are presented to prove the availability and efficiency of such application.

Introduction

The parareal algorithm is a time domain decomposition method aiming at solving PDEs over a splitting of the original domain into multiple subdomains, over which the equation could be solved independently. The algorithm is therefore suited for parallel programmin by distributing the computational cost over several processors. The parallelization with respect to the time variable was earlier attempted by Nievergelt [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF] and by Miranker and Liniger [START_REF] Miranker | Parallel methods for the numerical integration of ordinary differential equations[END_REF] for the ordinary differential equations. The parareal algorithm is another extension of decomposition methods in the time direction, first introduced in [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF] that has grown in popularity as it has been successfully applied to a wide variety of problems. The algorithm consists of two solvers with different time steps, the coarse and the fine, to produce a convergent iterative method for parallel computations. The coarse solver solves the equation sequentially on the coarse time step while the fine solver uses the results from the coarse solution to solve, in parallel, over the fine time steps. The solution can therefore be obtained from independent time subdomain, over which the algorithm will iteratively use a coarse and a fine solver to compute the solution on the entire domain. Several extensions have been proposed following the theoretical analysis. Bal and Maday in [START_REF] Bal | A "parareal" time discretization for non-linear pde's with application to the pricing of an american put[END_REF] slightly modified the algorithm, and proposed an analysis of the convergence properties of the algorithm for parabolic equations.

On the other hand, the asynchronous iterative scheme appeared on the horizon since Chazan and Miranker proposed their pioneering ideas applied to the linear systems [START_REF] Chazan | Chaotic relaxation[END_REF], leading to a prosperity to the high performance applications, such as sub-structuring methods [START_REF] Magoulès | Asynchronous iterative sub-structuring methods[END_REF] and optimized Schwarz methods [START_REF] Magoulès | Asynchronous optimized Schwarz methods with and without overlap[END_REF]. We are interested in the asynchronous iterations of the parareal method. Compared to the classical parareal algorithm, the asynchronous version gives further attractive properties including high efficiency and flexibility. Nevertheless, the implementation of the asynchronous iterative methods involves huge amount of work besides the theoretical analysis, which leads to a frustrating down time for the programming activities. Fortunately, some libraries have been developed to fill this gap. JACE Figure 1: Time domain of the parareal scheme [START_REF] Bahi | Jace: a Java environment for distributed asynchronous iterative computations[END_REF] is a multi-threaded library aiming at the execution of asynchronous iterative algorithms based on Java, followed by a centralized volatility tolerant extension named JaceV [START_REF] Bahi | JaceV: A programming and execution environment for asynchronous iterative computations on volatile nodes[END_REF]. Recently, they have developed some P2P and decentralized versions to enlarge the applicable scope [START_REF] Bahi | JaceP2P: an environment for asynchronous computations on peer-to-peer networks[END_REF][START_REF] Charr | JACEP2P-V2: A fully decentralized and fault tolerant environment for executing parallel iterative asynchronous applications on volatile distributed architectures[END_REF]. Additionally, CRAC [START_REF] Couturier | CRAC: a grid environment to solve scientific applications with asynchronous iterative algorithms[END_REF] is another library designed to build the asynchronous applications based on C++. Finally, JACK [START_REF] Magoulès | JACK: an asynchronous communication kernel library for iterative algorithms[END_REF] is an MPI-based C++ library, which provides many advanced properties designed for parallel iterative algorithms. A new version called JACK2 has been released. It gives an improved user-friendly API and several new convergence detection methods with global residual computation.

In this paper we concentrate on the implementation of asynchronous parareal solver based on JACK2 library. Section 2 recalls the mathematical and computational framework of the parareal algorithm. In Section 3, we presents the implementation details of such methods by JACK2. In Section 4, we choose the option pricing models to illustrate the experimental results of our asynchronous solver. Finally, our conclusions and remarks follow in Section 5.

Mathematical and Computational Framework

Consider a time-dependent problem

∂u ∂t + Lu = f, t ∈ [0, T ], u = u 0 , t = 0,
where L is a second-order linear elliptic operator. After a decomposition of the time domain, we arrive at

∂un ∂t + Lu n = f n , t ∈ [T n , T n+1 ], u n = λ n , t = T n , (1)
where the time domain is defined in a series of [T n , T n+1 ], with T n = n∆T , n = 0, . . . , N . It is easily seen that the problem ( 1) is also an exact model. Nevertheless, we often choose an accurate enough sub-interval δt to approximate the original problem as shown in Figure 1. We note here that the limit of each subproblem satisfies

λ n+1 = lim →0 u n (T n+1 -),
which enables us to execute a continuous computation. A well known model proposed in [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF] can be described as the following

λ k+1 n+1 = G(λ k+1 n ) + F (λ k n ) -G(λ k n ),
where G is a coarse propagator, and F is a fine propagator. Accordingly, ∆T represents a coarse time step, whereas δt is a fine time step. We notice that this is a synchronous parareal iterative scheme, over which one can come up with the algorithmic model forthwith

Algorithm 1. (Classical Parareal Algorithm) n := rank of current processor λ 0 = u 0 for i = 0 to n -1 do λ 0 = G(λ 0 ) end for w = G(λ 0 ) while not convergence do v = F (λ 0 )
wait for the update of λ 0 from processor n -1 w = G(λ 0 ) λ = w + v -w send λ to processor n + 1 as λ 0 w = w end while Evidently, the first processor has only successor, whereas the last one has none but predecessor.

We now turn to the asynchronous mode. Formally, we view parareal as a two-stage problem which can be formalized carefully using the asynchronous scheme, illustrated as below

λ k+1 n+1 = G(λ µn(k) n ) + F (λ ρn(k) n ) -G(λ ρn(k) n ), n ∈ P k , λ k n+1 , n / ∈ P k ,
where P k ⊆ {1, . . . , p} and P k = ∅, with the conditions 0 ≤ µ i j (k) ≤ k+1, 0 ≤ ρ i j (k) ≤ k. We assume that each processor keeps updating rather than resting permanently. Furthermore, the elements used by this processor will be renewed from time to time. Similarly, such computational model can be written as

Algorithm 2. (Asynchronous Parareal Algorithm) n := rank of current processor λ 0 = u 0 for i = 0 to n -1 do λ 0 = G(λ 0 ) end for w = G(λ 0 ) while not convergence do v = F (λ 0 ) if detect λ 0 from processor n -1 then update λ 0 end if w = G(λ 0 ) λ = w + v -w send λ to processor n + 1 as λ 0 w = w end while
It is seen that Algorithm 2 has no waiting process in the internal loop, whereas Algorithm 1 requires a latency time, which might cause a waste of time on coordination, but gain through the number of iterations.

Implementations

Preprocessing

The parareal algorithm is indeed a special case of domain decomposition methods, since each time frame only depends upon its predecessor, and essentially needed by its successor. We separate the neighbors into the outgoing links and the incoming links, and thus write the code as Listing 1. It is seen that sneighb_rank and rneighb_rank have respectively one element in the communication graph. Take note that the first processor has no predecessor, therefore numb_rneighb = 0; while the last one has no successor, obviously, with numb_sneighb = 0. The key idea behind this project is that there exist two levels over the computation process. The upper object is the instance of Parareal, which provides the interface of JACK2's power, whereas the lower level is the PDESolver, used to solve the equations. We notice that it is easy to enlarge the number of solvers without changing the code of Parareal. Accordingly, we need two instances of PDESolver in Parareal to obtain respectively the coarse results and the fine results, illustrated in Listing 2.

Listing 2: PDE solver in Parareal /* template <typename T, typename U> */ PDESolver<T,U> coarse_pde; PDESolver<T,U> fine_pde; Vector<T,U> coarse_vec_U; Vector<T,U> fine_vec_U; Vector<T,U> vec_U; // solution vector. Vector<T,U> vec_U0; // initial vector.

coarse_vec_U is the solution of coarse_pde, while fine_vec_U is the solution of fine_pde. vec_U is the final solution, which is also the output of Parareal solver. Finally, vec_U0 is the incoming value for both coarse_pde and fine_pde.

Overview of the Iterations

We now present the iterative process of parareal algorithm. Notice that the beginning parts of Algorithm 1 and Algorithm 2 are the same, we can therefore write in a unified scheme in Listing 3. Let us mention here that Listing 3 gives us the exact position of vec_U0 and the initial estimation of G(λ 0 n ). Nevertheless, we need to distinguish the synchronous and the asynchronous parareal algorithm afterwards.

Synchronous mode

We follow the scheme of Algorithm 1 except for the convergence detection mode, since there is a rule of thumb to lighten the communication. It is easily seen that the processor n will stop updating after (n + 1) th iteration. Hence, we can apply a supplementary condition to the judging area in the synchronous mode, with the final code shown as Listing 4.

Listing 4: Synchronous iterative process // --synchronous parareal iterations numb_iter++; }

Experimental Results

In this section, we give the experiments for the asynchronous parareal scheme, based on the Black-scholes equation, which is a well-known time-dependent equation in the domain of option pricing. Consider the following problem

∂V ∂t + rS ∂V ∂S + 1 2 σ 2 S 2 ∂ 2 V ∂S 2 = rV,
where V is the option price, depending on stock price S and time t. Volatility σ and riskfree interest rate r are the constant parameters. Then, we employ some trivial variable substitutions to obtain the heat equation

∂u ∂τ = ∂ 2 u ∂x 2 , x ∈ R, τ ∈ [0, T σ 2 2 ], (2)
where T is the time to maturity, with initial and boundary conditions u(x, 0) = max(e βx -e αx , 0),

x ∈ R, u(x, τ ) ∼ 0 as x → ±∞, τ ∈ [0, T σ 2 2 ].
where α = 1 2 (κ -1), β = 1 2 (κ + 1), κ = 2r σ 2 . Obviously, equation ( 2) can be applied to the classical and the asynchronous parareal scheme. We note here that the input variables are the current stock price S and the exercise price E. Finally, we broaden the spatial domain and let the dimension correspond to the precision of our problem.

The mathematical operations are supported by Alinea library [START_REF] Magoulès | Alinea: An advanced linear algebra library for massively parallel computations on graphics processing units[END_REF], which is implemented in C++ for both central processing unit and graphic processing unit devices. It includes several linear algebra operations [START_REF] Ahamed | Fast sparse matrix-vector multiplication on graphics processing unit for finite element analysis[END_REF] and numerical linear algebra solvers [START_REF] Magoulès | Fast iterative solvers for large compressed-sparse row linear systems on graphics processing unit[END_REF], [START_REF] Ahamed | Iterative methods for sparse linear systems on graphics processing unit[END_REF], [START_REF] Magoulès | Auto-tuned Krylov methods on cluster of graphics processing unit[END_REF].

The experiments are exercised on a SGI ICE X cluster connected with InfiniBand. Each node consists of two Intel Xeon E5-2670 v3 2.30 GHz CPUs with SGI-MPI 2.14 installed. We assume that σ = 0.2, r = 0.03, δt = 0.001 with 250 sub-intervals. Given S = 100, E = 80, N = 16, several average results are reported in Table 1 with Approximate Option Prices V a , Exact Option Prices V e , Absolute Error a and Relative Error r . It is seen that

Conclusions

In this paper, we illustrated the implementation of synchronous and asynchronous parareal algorithm using JACK2, an asynchronous communication kernel library for iterative algorithms. We discussed the particularity of such time-dependent problem, which leads to a rather different configuration than the general parallel context. As expected, the experimental results at the end verified the correctness of the asynchronous scheme.

Listing 1 :

 1 Communication graph /* template <typename T, typename U> */ // T: float, double, ... // U: int, long, ... U numb_sneighb = 1; U numb_rneighb = 1; U* sneighb_rank = new U[1]; U* rneighb_rank = new U[1];

Listing 3 :

 3 Initialization of parareal iterative process /* template <typename T, typename U> */ for (U i = 0; i < rank; i++) { coarse_pde.Integrate(); vec_U0 = coarse_vec_U; } coarse_pde.Integrate(); vec_U = coarse_vec_U;

Table 1 :

 1 Results of Asynchronous Parareal Scheme (σ = 0.2, r = 0.03, δt = 0.001, S = 100, E = 80, m = 250, N = 16) 23.9476 23.9426 0.0050 0.0002 0.547 0.15 31.5512 31.5477 0.0035 0.0001 1.663 0.25 37.7225 37.7192 0.0033 0.0001 2.881 0.35 42.9995 42.9960 0.0035 0.0001 3.830 0.45 47.6375 47.6339 0.0036 0.0001 4.814 the asynchronous parareal algorithm implemented by JACK2 is efficient and accurate.

We note that it is necessary to copy the values of recv_buf to vec_U0 and the values of vec_U to send_buf whenever needed. However, we can also ease the process by setting directly the two vectors as communication buffers. In practice, the choice relates to the mathematical library used throughout the project. Moreover, using the results of Listing 4, it is also shown in Listing 5 that we need to finalize the fine solution and wait for the global termination.

Listing 5: Synchronous finalized process

Asynchronous mode

The asynchronous parareal scheme can be implemented similarly with Listing 4 except for the matching conditions. There is no more facile code in practice, and we present the process as Listing 6. From the entry of lconv_flag, JACK2 can interact with the asynchronous process and give back results as res_norm.

Listing 6: Asynchronous iterative process


L Von Helden

L Bogula

P.-E Janolin

M Hanke

T Breuer

M Schmidbauer

S Ganschow

J Schwarzkopf

Huge impact of compressive strain on

We present a study in which ferroelectric phase transition temperatures in epitaxial K x Na 1-x NbO 3 films are altered systematically by choosing different (110)-oriented rare-earth scandate substrates and by variation of the potassium to sodium ratio. Our results prove the capability to continuously shift the ferroelectric-to-ferroelectric transition from monoclinic M C to orthorhombic c-pha se by about 400 °C via application of anisotropic compressive strain. The phase transition was investigated in detail by monitoring the temperature-dependence of ferroelectric domain patterns with piezoresponse force microscopy (PFM) and upon analyzing structural changes by means of high resolution X-ray diffraction (HR-XRD) including X-ray reciprocal space mapping. Moreover, the temperature evolution of the effective piezoelectric coefficient d 33,f was determined using double beam laser interferometry exhibiting a significant dependence on the particular ferroelectric phase.

The key to predict and subsequently deliberately tailor phase transition temperatures in epitaxial K x Na 1-x NbO 3 thin films is to understand the relationship between phase formation, transition temperature and incorporated lattice strain. However, no experimental investigation addressing this topic has been available for this material, yet. In this paper, we present a study in which the impact of different anisotropic epitaxial strains in K x Na 1-x NbO 3 films on phase transition temperatures is investigated. In this context, the choice of substrate material allows for a gradual modulation of strain while the change of potassium to sodium ratio enables a continuous and precise fine tuning. Via the combination of temperature-dependent piezoresponse force microscopy (PFM), temperature-dependent in-situ X-ray diffraction (XRD) and laser interferometry we show that we can indeed continuously shift the transition from the ferroelectric, monoclinic M C to the ferroelectric c-phase in the range of -15 °C to 400 °C. K x Na 1-x NbO 3 films with a thickness of 30 (±10) nm were grown on (110) DyScO 3 , (110) TbScO 3 , (110) GdScO 3 and (110) SmScO 3 substrates by means of liquid-delivery metalorganic vapor phase epitaxy (MOVPE). This technique allows for growth conditions close to the thermodynamic equilibrium providing smooth film surfaces and highly regular ferroelectric domain formation. The samples were prepared with a molar potassium content x ranging from 0.54 to 0.77, which was varied by adding different amounts of K and Na precursors to the MOVPE growth process, which was carried out at a growth temperature of 700 °C. Details are described in Ref. 4. The rare-earth scandate substrates have been grown by Czochralski technique as explained in Ref. 21.

The ferroelectric domain patterns of K x Na 1-x NbO 3 films on different substrates have been analyzed using the lateral piezoresponse force microscopy (LPFM) mode of an Asylum MFP3Dinstrument with dual AC resonance tracking (DART). 12 In Figs. 1(a)-1(d) representative piezoresponse (PR=amplitude ⋅ sin(phase)) micrographs recorded at room temperature are displayed. The K x Na 1-x NbO 3 film shown in Fig. 1(a) grown on (110) DyScO 3 does not exhibit a periodic domain structure at room temperature. In contrast the films presented in Figs. 1(b)-1(e) reveal periodic stripe domain patterns (indicated by white double arrows) forming up to four superdomain variants. A detailed structural analysis of the domain configuration as described in Ref. 11 proves monoclinic M C domains with Pm symmetry. This structure contains four types of superdomains, which are denoted as 0°, 90°, 180° and 270° variants, as sketched in Fig. 1(e). Each of them is formed by lamellar nanodomains whose in-plane polarization components alternate by 90° as indicated by small, red arrows. The resulting in-plane net polarization of the superdomains is schematically given by blue arrows. In order to determine the induced epitaxial strain in K x Na 1-x NbO 3 thin films it is most convenient to describe the film structure in pseudocubic (pc) notation. 13 The strain state of each specimen is determined as follows: first, the variation of pseudocubic (pc) bulk lattice parameters a pc , b pc , and c pc is calculated by applying Vegard's law 14 for 0.5 < x < 1.0 with the end members given in Table 1. 15,16 Second, we analyzed the vertical lattice parameter by X-ray diffraction, proving c pc -orientation for all tested films, which allows for an evaluation of the chemical composition and thus the experienced in-plane strains ε xx and ε yy . Thereby, a Poisson ratio of = 0.38 as reported for K 0.5 Na 0.5 NbO 3 by Matsubara et al. 17 is applied. This calculation method was found to yield almost identical values to the alternative evaluation based on the elastic compliance. However, both methods produce a systematic uncertainty since the Poisson ratio and the elastic compliance are only available for the orthorhombic phase of K 0.5 Na 0.5 NbO 3 . Nonetheless, they allow for a reliable comparison among the different K x Na 1-x NbO 3 films investigated in this study. For the (110)-oriented rare-earth scandate (REScO 3 ) substrates the average in-plane unit cells parameters are taken (Tab. 1), since the possible in-plane unit cell orientations of the film occur with equal likelihood in an M C domain structure. The potassium contents and in-plane strain values for the samples shown in Figs. 1(a Avg. in-pl. lattice parameter (Å) DyScO 3 21 3.949 TbScO 3 21 3.959 GdScO 3 21 3.968 SmScO 3 21 3.987 In order to depict the characteristics of the phase transition in these films upon heating, its investigation by means of different techniques will be presented in the following on the example of the K 0.66 Na 0.34 NbO 3 film on (110) TbScO 3 which was presented in Fig. 1

(Table 1: Pseudocubic (pc)

(b). In Figs. 2(a)-2(d)

LPFM measurements recorded within one particular superdomain of the named film reveal the temperature-dependence of the ferroelectric domain pattern between 80 °C and 120 °C. Although the noise in PFM measurements increases when heating the sample, the stripe domains can be identified up to 110 °C before they abruptly vanish at 120 °C. Upon cooling back down to 110 °C the stripe domain pattern in PFM reappears, Fig. 2(d). Due to thermal drift of the sample stage it was not viable to perform temperature-dependent PFM measurements at a fixed sample position, e.g., at a superdomain boundary recording different superdomains at once. However, individual investigations of the particular superdomains proved that all variants vanish at the same temperature as expected for an M C domain structure where all superdomains hold the same overall in-plane strain, which we define as (ε xx + ε yy ).

For a more detailed determination of the observed phase transition we have conducted temperature-dependent X-ray diffraction experiments. High-resolution X-ray diffraction (HRXRD) under elevated temperatures was performed on a diffractometer equipped with an 18-kW rotating-anode generator in Bragg-Brentano geometry and using synchrotron radiation at the PHARAO beamline of BESSY II. Additional measurements under cryogenic temperatures were performed on a Bruker Discover D8 diffractometer equipped with a Bruker MTC-Lowtemp stage. Furthermore, X-ray reciprocal space mappings (RSMs) were collected at the KMC2 beamline of BESSY II applying an Anton Paar DHS 1100 hot stage.

As known from several studies on various ferroelectric oxide materials like K x Na 1-x NbO 3 , BaTiO 3 or PbTiO 3 the thermal expansion coefficients in different crystalline phases of ferroelectrics significantly differ, which allows identifying phase transitions. 18,19,20 In Fig. 3(a) the evolution of the vertical lattice parameter of the K 0.66 Na 0.34 NbO 3 film presented in Fig. 1(a) and Fig. 2, and of the underlying (110) TbScO 3 substrate are shown in dependence of the temperature. While the substrate reveals a rather linear thermal expansion with an expansion coefficient of about 7 10 -6 K -1 well compatible with literature values, 21 a distinct change of the slope is observed which corresponds to the phase transition in the K 0.66 Na 0.34 NbO 3 film. The position of the kink at 115 °C is in good agreement with the phase transition temperature observed by means of LPFM. Below 115 °C the vertical lattice parameter strongly increases. For T > 115 °C the curve is clearly flattened which indicates a very small or even negative thermal expansion. This effect is known for ferroelectrics like PbTiO 3 22 or BaTiO 3 23 and is attributed to the spontaneous volume ferroelectrostriction. 22,24,25 While this effect is present in both phases, its impact is more pronounced when the polarization vector coincides with the measured unit cell axis. As here the vertical lattice parameter is measured, this indicates a transition into a c-phase where the polarization vector is aligned along the c-axis of the unit cell. In order to clarify the film structure below and above the phase transition, temperaturedependent X-ray reciprocal space mappings were performed in the vicinity of the asymmetric (444 ) TSO substrate Bragg reflection containing both in-plane and out-of-plane scattering vector components, as shown in Figs. 3(b) and 3(c). At room temperature, the intensity distribution consists of a complex diffraction pattern resulting from the periodic arrangement of the M C domains. Their monoclinic shearing can be identified from the different vertical positions of the two satellite branches (dashed, white lines), which do not intersect in their center position (red dot). 11 At 150 °C well above the transition temperature, this diffraction pattern has vanished confirming the absence of a periodic domain structure in agreement with the LPFM images shown in Fig. 2. Meanwhile, the center position of the X-ray intensity scattered by the film is still at the same in-plane q 001 position as the (444 ) TSO substrate Bragg reflection proving that the film does not undergo any plastic relaxation but perfectly preserves its fully strained relation to the substrate. However, no peak splitting can be identified in the intensity pattern anymore. This implies the absence of a monoclinic distortion, proving an orthorhombic symmetry of the K x Na 1-x NaNbO 3 unit cells in the high temperature phase. Note that the (110) TbScO 3 surface exhibits an anisotropic unit cell, so no tetragonal symmetry is possible in the clamped layer.

To evaluate the ferroelectric properties below and above the phase transition, macroscopic measurements of the effective, vertical piezoelectric coefficient d 33,f were performed with an aixACCT double beam laser interferometer in combination with a TF Analyzer 2000 and a Linkam HFS temperature stage. Thereby, so-called small-signal measurements were conducted by superimposing an AC electric field with 0.2 V amplitude and a frequency of 1 kHz to a DC field, which was gradually varied in the range of ±0.5 V. 26 This technique allows one to identify the material's response mainly caused by the intrinsic reverse piezoelectric effect. The measurements were conducted for a K 0.58 Na 0.42 NbO 3 film on a previously grown SrRuO 3 bottom electrode layer on (110) TbScO 3 with Pt top electrodes.

In Fig. 4, the evolution of the piezoelectric coefficient d 33,f is displayed in the temperature range between 20 °C and 180 °C. From room temperature up to 50 °C d 33,f is fairly constant with a value of 18 (±1) pm/V, which is comparable to values for clamped, lead-free BaTiO 3 thin films. 27 Above this temperature until 70 °C a distinct decrease is observed. For the high temperature cphase an effective piezoelectric coefficient of 7 (±2) pm/V is found. This behavior is in good agreement with the expectation to gain the highest piezoelectric response for the monoclinic phase. We further note that the low transition temperature in this particular sample might be attributed to screening charges in the additional metallic intermediate layer, which are known to stabilize tetragonal phases. [28][29][30] Accordingly, PFM measurements confirm the transition from a periodic stripe pattern to an irregular domain arrangement below 60 °C (insets in Fig. 4).

(Figure 4: Macroscopic effective vertical piezoelectric coefficient d 33,f as function of temperature measured for a K 0.58 Na 0.42 NbO 3 film on a SrRuO 3 bottom electrode layer grown on (110) TbScO 3 .) Summarizing the PFM, HR-XRD and interferometer results, we conclude a thermally induced phase transition from the ferroelectric, monoclinic M C phase into a ferroelectric, orthorhombic c-phase in K x Na 1-x NbO 3 films on (110) TbScO 3 . From several experimental and theoretical studies, it is well-known that the phase transition temperature in a ferroelectric material can be systematically shifted by the application of strain. [31][32][33] Therefore, we have grown several K x Na 1-x NbO 3 thin films with a molar potassium content x ranging from 0.54 and 0.77 on different rare-earth scandate substrates (without intermediate bottom electrode). In Fig. 5 the phase transition temperatures for K x Na 1-x NbO 3 films with an overall compressive in-plane strain between -1.8 % and 0.0 % recorded by means of different techniques are gathered. We remark that all provided strain values were calculated for room temperature. The strain change upon temperature variation due to stronger thermal expansion coefficients of the substrates compared to the layers is less than 0.1 % per 100 K. The evolution of the transition temperature with increasing compressive strain unveils a systematic and apparently linear dependence, whereby a shift over a wide temperature range from -15 °C to 400 °C was experimentally realized.

The lowest phase transition temperature is found for the K 0.71 Na 0.29 NbO 3 film on (110) DyScO 3 being significantly below room temperature. In contrast, the highest phase transition temperature, which is observed for the K 0.74 Na 0.26 NbO 3 film on SmScO 3 , is already in the vicinity of the Curie temperature of bulk K 0.5 Na 0.5 NbO 3 of T C = 415 °C. [34] However, we suppose to observe a ferroelectric-to-ferroelectric phase transition also in this case. Even though this sample is assigned a zero overall in-plane strain, it is fully epitaxially and highly anisotropically strained with ε xx = 0.8 % and ε yy = -0.8 %. It is widely accepted that any strain, which draws the unit cell symmetry away from being cubic, leads to an increase of the Curie temperature in perovskite-like ferroelectrics. For this reason, the ferroelectric-to-paraelectric phase transition in the layer is expected to have a significantly higher T C than the unstrained bulk material. In between the DyScO 3 and SmScO 3 samples the influence of the varying molar potassium amount can be nicely obtained by comparing the different K x Na 1-x NbO 3 films on, e.g., (110) GdScO 3 . Here, the variation of x between 0.54 and 0.77, which corresponds to a change of the overall in-plane strain from -1.0 % to -0.3 %, leads to a shift of the phase transition temperature from 105 °C to 300 °C. We further note that the slope of the phase transition temperature versus overall inplane strain (ca. 200 °C/%) is in quite good agreement with calculations recently published by Zhou et al.. 35 However, we underline the restricted comparability between both studies, since Zhou et al. conducted their theoretical investigation on isotropically strained monodomain K 0.5 Na 0.5 NbO 3 revealing an M A to c-phase transition. In contrast, the system of anisotropically strained K x Na 1-x NbO 3 films with polydomain configuration as present in our films holds a distinctly higher complexity.

In summary, we have studied a thermally induced phase transition in K x Na 1-x NbO 3 thin films with 0.54 < x < 0.77 grown by MOVPE. Compressive strains were varied between -1.8 % and 0.0 % by using different (110)-oriented REScO 3 substrates (with Re = Dy, Tb, Gd, Sm). By applying (Figure 5: Phase transition temperature of the M C to c-phase transition versus overall in-plane strain for K x Na 1-x NbO 3 films with x varying between 0.54 and 0.77 grown on different (110)-oriented rareearth scandate substrates. Blue, red, violet and green color indicate the use of DyScO 3 , TbScO 3 , GdScO 3 , and SmScO 3 substrates, respectively. The transition temperatures have been determined by LPFM (circles), via the kink in thermal expansion measurements (stars) and by X-ray RSMs (squares).) temperature-dependent PFM, different X-ray diffraction methods and laser interferometry we verified a ferroelectric-to-ferroelectric phase transition from monoclinic to orthorhombic symmetry in these films. Thereby, a linear correlation between the phase transition temperature and the applied overall in-plane lattice strain has been revealed. Our results underline the great potential to strain engineer K x Na 1-x NbO 3 layers for technological application. The detailed knowledge of the relation between lattice strain, phase formation and phase transition temperature will enable the customization of functional properties like the transmission coefficients to the particular temperature regime in which an SAW sensor is desired to operate.

(Figure 1 :

 1 Lateral piezoresponse force micrographs for K x Na 1-x NbO 3 films at room temperature presented in order of decreasing compressive strains on (a) (110) DyScO 3 , (b) (110) TbScO 3 , (c) (110) GdScO 3 , (d) (110) SmScO 3 . White double arrows indicate the orientation of stripe domains. (e) Schematic presentation of an M C domain configuration. The in-plane polarization within the stripe domains is represented by red arrows and the resulting in-plane net polarization of superdomains by blue arrows.)

  )-1(d) are as follows: (110) DyScO 3 with x = 0.71, ε xx = -0.1 % and ε yy = -1.7 %, (110) TbScO 3 with x = 0.67, ε xx = 0.2 % and ε yy = -1.3 %, (110) GdScO 3 with x = 0.54, ε xx = 0.6 % and ε yy = -0.9 %, (110) SmScO 3 with x = 0.74, ε xx = 0.8 % and ε yy = -0.8 %. a pc (Å) b pc = c pc (

(Figure 2 :

 2 lattice parameters of K x Na 1-x NbO 3 for x = 0.5 and x = 1.0 and averaged in-plane lattice parameters (1/4 • (d[ 11 0 ] REScO3 +d[ 001 ] REScO3 )) of different rare-earth scandate substrates.) Lateral piezoresponse force micrographs for the K 0.66 Na 0.34 NbO 3 film on (110) TbScO 3 presented in Fig. 1(b) measured within one particular superdomain at different temperatures (a)-(d). Upon heating above 110 °C the periodic M C domain structure vanishes until it reappears when cooling down again.)

(Figure 3 :

 3 (a) Vertical lattice parameter of the K 0.66 Na 0.34 NbO 3 film shown in Fig. 1(b) and Fig. 2 and of the TbScO 3 substrate versus temperature. (b)-(c) X-ray reciprocal space maps in the vicinity of the (444 ) TSO Bragg reflection at 25 °C and 150 °C, respectively.)

We thank DFG (project funding reference FE 1438|2-1) and EFRE (Project number 1.8/15) for funding the project. We further thank A. Kwasniewski, D. Többens, Z. Cheng, and T. Demuth for assistance with X-ray measurements and M. Klann for MOCVD sample preparation. We are grateful to BESSY (KMC-2 projects 17206091-ST, 18106401-ST, and in-house research at PHARAO) for providing beam time. Finally, we acknowledge J. Boschker and R. Wördenweber for fruitful discussion.


Charlène Puzin

Julien Pétillon
email: julien.petillon@univ-rennes1.fr

Contrasted responses of dominant ground-dwelling arthropods to landscape salt-marsh fragmentation

Keywords: fragmentation, Araneae, Carabidae, Amphipoda, pitfall traps, grazing Western Europe, France, Normandy,

Contrasted responses of dominant ground-dwelling arthropods to landscape salt-marsh fragmentation

Charlène Puzin, Julien

Introduction

Salt marshes play important roles in terms of environmental functionality and conservation.

Their primary productivity is one of the highest in the world (up to 3900 g C m -2 yr -1 : [START_REF] Laffaille | Composition of fish communities in a European macrotidal salt marsh (the Mont Saint Michel bay, France)[END_REF]. As other estuarine and coastal ecosystems such as mangroves or dunes, salt marshes provide many environmental and economical benefits. Indeed they act like a physical barrier against erosion, protecting lands against waves and storm events [START_REF] Adam | Saltmarshes in a time of change[END_REF]. They also provide a useful shelter for growth and survival of young fishes, act like natural filter to purify water and are used by humans for fisheries, harvesting or pasture [START_REF] Garbutt | European salt marshes: ecology and conservation in a changing world[END_REF]. They are yet submitted to several human threats: climate change, sea level rise, increasing sea and air temperature, pollution, vegetation disturbance, biological invasion or eutrophication. The degradation of 50% of worldwide salt marshes that has been reported over the last three decades has become a conservation issue, especially in Western Europe [START_REF] Dijkema | Salt marshes in Europe[END_REF]. Salt marshes indeed host several specialist species (i.e. only found there), some of them being also rare with low abundances and/or restricted distribution (e.g. [START_REF] Irmler | Zonation of ground beetles (Coleoptera: Carabidae) and spiders (Araneida) in salt marshes at the North and the Baltic M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Sea and the impact of the predicted sea level increase[END_REF]. They also constitute ideal systems to test for the effects of habitat fragmentation on biodiversity because they have a linear, highly patchy distribution along coasts [START_REF] Desender | Diversity and conservation of terrestrial arthropods in tidal marshes along the River Schelde: a gradient analysis[END_REF]. Their fragmented repartition severely limits exchanges between populations from different salt marshes. Locally, they are submitted to various human disturbances (grazing, mowing, invasive species) that created within-site heterogeneity, and therefore differences in habitat suitability (either at local or at landscape scales). In Western Europe, salt marshes are mainly covered by Atriplex portulacoides, but for the last thirty years this natural vegetation is more and more replaced by the invasive species Elymus athericus, which was at first limited to the upper part of the marsh [START_REF] Valéry | Impact of the invasive native species Elymus athericus on carbon pools in a salt marsh[END_REF], leading to a patchy landscape (Veeneklaas et al., 2012). This invasion does not seem to have strong negative impacts on diversity, as it favors in fact the settlement of continental

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
species, but it undoubtedly modifies the functioning of the salt-marsh ecosystem and fragments populations living in A. portulacoides [START_REF] Pétillon | Impact of an invasive plant (Elymus athericus) on the conservation value of tidal salt marshes in western France and implications for management: responses of spider populations[END_REF]. Sheep grazing also induces frequent changes in salt-marsh natural vegetation, A. portulacoides being replaced by the short grass Puccinellia maritima. If local effects of both invasions by Elymus athericus and grazing were quite studied in the past (e.g. [START_REF] Kiehl | Impact of different sheep grazing intensities on salt marsh vegetation in northern Germany[END_REF]Veeneklaas et al., 2012;[START_REF] Leroy | Complementarity of rarity, specialisation and functional diversity metrics to assess responses to environmental changes, using an example of spider communities in salt marshes[END_REF], no study investigated how the fragmentation of natural landscape they induced affects dominant species of arthropods so far.

In this study, we wanted to assess the consequences of within-site fragmentation (occurring at a scale inferior to 1km²) on salt-marsh biodiversity. The effects of patch size was investigated on the dominant species of ground-active arthropods. We first expect the effect of patch size to be different among species, weaker of salt-marsh predators with a high mobility [START_REF] Bonte | Aerial dispersal plasticity under different wind velocities in a salt marsh wolf spider[END_REF]Van Belleghem et al., 2015) and stronger for salt-marsh decomposers with a low mobility [START_REF] Mantzouki | Accuracy of pitfall traps for monitoring populations of the amphipod Orchestia gammarella (Pallas 1766) in a salt marsh[END_REF]. Because patches were distributed along transects set in two contrasted environments, we then expect an especially negative effect of the grazed environment on their abundances (see e.g. [START_REF] Rickert | αand β-diversity in moth communities in salt marshes is driven by grazing management[END_REF]Ford et al., 2012;van Klink et al., 2013).

Material and Methods

Sampling design

Sampling took place in the eastern part of the Mont Saint-Michel Bay (48°38′10″N; 1°30′41″W). To assess if there was an effect of patch size, three sizes of relictual patches of A. portulacoides were sampled: small (area<0.1ha), medium (0.1ha<area<1ha), and large (area>1ha). For each patch size, three replicates were sampling in both environments. Thus, there were nine patches of A. portulacoides distributed along each of two transects (see the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT map in ESM: Fig. A.1
). The impact of environment (i.e. surrounding vegetation of the patch) on abundances was also considered because transects were surrounded by either invaded or grazed vegetation (around 10 sheep per hectare, see [START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript Pétillon | Impact of cutting and sheep-grazing on ground-active spiders and ground beetles in some intertidal salt marshes (Western France)[END_REF] for a more detailed description of the study site). In each patch, four pitfall traps were placed at 10 m from each other (to avoid interference between them). Pitfall traps were polypropylene cups (10 cm diameter, 17 cm deep) filled of ethylene-glycol. All 64 traps were visited weekly from April 27 th to June 15 th of 2012, except during one week, when salt marshes were inundated due to strong tides, for a total of 5 sampling weeks. Because catches in pitfall are both function of local densities and mobilities of arthropods, abundances of the dominant species (based on the 5 sampling weeks pooled together) will be referred to "activity-densities" further.

Data analyses

A total of 66,299 arthropods were sorted and identified to family or species level. Three species, the wolf spider Pardosa purbeckensis (Araneae, Lycosidae), the ground beetle Pogonus chalceus (Coleoptera, Carabidae) and the beach-hopper Orchestia gammarellus (Amphipoda, Talitridrae) were the dominant ground-dwelling macro-arthropods, representing 96% of all individuals caught, with 30584, 2055 and 30989 individuals respectively).

In order to assess differences in the activity-density of the three main species according to patch size and environment, quasi-Poisson Generalized Linear Models were performed with activity-density of the species as dependent variable, patch size and environment and their interaction as fixed factors. Post-hoc Tukey tests, with Bonferroni correction, were performed in case of patch size significant effect. All data analyses were performed using R software packages (R Development Core Team, 2013).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Results

Patch size had no effect on the activity-density of P. purbeckensis (Table 1, Fig. 1). It had a significant effect on the activity-density of P. chalceus and O. gammarellus, with more individuals in large patches than in small ones independent of the transect (Figs. 2 & 3,respectively).

Environment had a significant effect on the activity-density of P. purbeckensis, P. chalceus, and O. gammarellus (Table 1), with systematically more individuals in patches surrounded by an invaded than by a grazed environment (Figs 1, 2 & 3, respectively).

Discussion

The two predator species had very contrasted responses to changes in patch size. The latter had indeed no effect on the wolf spider Pardosa purbeckensis, and a significant effect on the ground beetle Pogonus chalceus, with higher abundances in larger patches. This result should certainly be related to the mobility of the two species, P. purbeckensis having high dispersal capacities (for both short-and long-distances, see [START_REF] Richter | Differences in breeding and motility between Pardosa pullata (Clerk) and Pardosa prativaga (L. Koch), (Lycosidae, Araneae) in relation to habitat[END_REF]Puzin et al., 2018, respectively) and P. chalceus being mostly brachypterous in the Mont Saint-Michel Bay [START_REF] Desender | Age and size of European saltmarshes and the population genetic consequences for ground beetles[END_REF]. Salt-marsh fragmentation has consequently more drastic effects for low-mobility species. The short-distance dispersal and ground-level movements of P. chalceus are still to be investigated, yet.

Although transect interacted with patch size for Orchestia, results are quite similar for both environments, i.e. more individuals in the largest patches of natural vegetation. The fact that there was no difference in activity-density for both the amphipod and the carabid species between small-and medium-sized patches might indicate that a minimum patch-size is

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

needed for these two species. Dispersal of Orchestia gammarellus is probably water-mediated (see [START_REF] Fanini | Coastal talitrids and connectivity between beaches: a behavioural test[END_REF] on closely related species) and the low mobility of this species probably explains the strong effect of patch size on population density. This result has strong consequences for salt-marsh functioning, Orchestia being indeed an important prey item for several fish species, including some of high commercial value like the sea bass [START_REF] Laffaille | Does the invasive plant Elymus athericus modify fish diet in tidal salt marshes? Estuarine[END_REF]. Interestingly, it may help understanding why relictual areas of natural vegetation (e.g. like along the creeks) do not host amphipods enough to fulfil the nursery function, especially reduced in a grazed salt marsh [START_REF] Laffaille | Composition of fish communities in a European macrotidal salt marsh (the Mont Saint Michel bay, France)[END_REF].

Transect had an important effect on the activity-density of all dominant species in natural patches, populations being less dense in a grazed environment than in an invaded salt marsh.

Although our design did not allow for spatial replication and statistical inference, we assume that the differences found between transects can be related to the effects of invasion vs grazing. Differences in humidity and organic matter (both lowered in grazed salt marshes) can also explain the differences in activity-density between environments for the three model species. That would overall confirm at a broader (landscape) level what was previously reported at a local scale, i.e. over-grazing is actually more deleterious to salt-marsh specialist arthropods than invasions by Elymus athericus (with very few exceptions, like e.g. the aeronaut linyphiid spider Erigone longipalpis, found to be more abundant in heavily grazed salt marshes: [START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript Pétillon | Impact of cutting and sheep-grazing on ground-active spiders and ground beetles in some intertidal salt marshes (Western France)[END_REF]. Intensive grazing usually leads to extremely homogeneous habitats, and the low variance of data in the grazed environment for all three species also supports this statement. Our study finally shows that the impacts of both grazing and plant invasion are not only local, but also influence specialist biodiversity at a landscape scale though the fragmentation they induce, which was not previously reported in salt marshes.

The way fragmentation had such adverse effects on population size of salt-marsh specialist arthropods is still to be investigated, and can notably encompass changes in micro-climate

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
and/or strong edge effects. We limited our sampling to spring because it is the best season for sampling the target species, but patterns are likely similar in autumn (see [START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript Pétillon | Impact of cutting and sheep-grazing on ground-active spiders and ground beetles in some intertidal salt marshes (Western France)[END_REF]).

Yet, our results are based on one sampling during one year only, and thus need to be confirmed by a long-term monitoring. Further studies should also get interested in interspecific dynamics in such patchy environments, and especially changes in predator-prey densities. Metacommunity theory was indeed mostly tested on assembly rules, and more rarely on interspecific relationships [START_REF] Vanschoenwinkel | The role of metacommunity processes in shaping invertebrate rock pool communities along a dispersal gradient[END_REF]. Salt marshes are thus ideal model systems for such studies because they experience a high level of stress due regular tides, and have a highly patchy structure, both between and within sites. • The effects of landscape on arthropods has been little studied in salt marshes

• The response of populations to fragmentation was assessed for 3 dominant species

• Patch size of natural vegetation had an effect on low-mobility species only

• More individuals were found in patches surrounded by invaded than by grazed marshes

• Results suggest that disturbances impact arthropods also at a landscape scale

  S.M.,Roelofs, D., Hendrickx, F. 2015. Evolutionary history of a dispersal associated locus cross sympatric and allopatric divergent populations of a wingpolymorphic beetle across Atlantic Europe. Molecular Ecology 24, 890-908.van Klink, R., Rickert, C., Vermeulen, R., Vorst, O., WallisDeVries, M.F., Bakker, J.P., 2013.Grazed vegetation mosaics do not maximize arthropod diversity: Evidence from salt marshes. Biological Conservation 164, 150-157.Veeneklaas, R.M., Dijkema, K.S., Hecker, N., Bakker, J.P., Schaminée, J., 2013. Spatiotemporal dynamics of the invasive plant species Elytrigia atherica on natural salt marshes. Applied Vegetation Science 16, 205-216.

Fig. 1 .

 1 Figure captions

Fig. 2 .

 2 Fig. 2. Activity-density of Pogonus chalceus according to patch size and environment.

Fig. 3 .

 3 Fig. 3. Activity-density of Orchestia gammarellus according to according to environment and

  Figure 1. Puzin & Pétillon

Table 1 .

 1 GLMs of patch size (small, medium or large) and environment (surrounded by invaded or grazed vegetation) effects on the activity-density of the dominant ground-dwelling arthropods in salt marshes (significant effects at α=0.05 are in bold). F-ratio and P-values of 1 st order terms are those from the full model if the interaction was significant, and those from the additive model if the interaction was non-significant.

	Arthropod species	Factor

Acknowledgments

We are grateful to Etienne Camenen, Romain Calleau, El Aziz Djoudi, Anita Georges and Romain Lahaye for help in field and sorting, and Dries Bonte for useful comments. This study was partly founded by the Région Bretagne (ARED 7023) and by the Agence de l'eau Loire-Bretagne ("PEPS" project).

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at