Strong Stability Preserving IMEX Methods for Partitioned Systems of Differential Equations

We investigate strong stability preserving (SSP) implicit-explicit (IMEX) methods for partitioned systems of differential equations with stiff and nonstiff subsystems. Conditions for order p and stage order q=p\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=p$$\end{document} are derived, and characterization of SSP IMEX methods is provided following the recent work by Spijker. Stability properties of these methods with respect to the decoupled linear system with a complex parameter, and a coupled linear system with real parameters are also investigated. Examples of methods up to the order p=4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=4$$\end{document} and stage order q=p\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=p$$\end{document} are provided. Numerical examples on six partitioned test systems confirm that the derived methods achieve the expected order of convergence for large range of stepsizes of integration, and they are also suitable for preserving the accuracy in the stiff limit or preserving the positivity of the numerical solution for large stepsizes.


Introduction
Strong stability preserving (SSP) implicit-explicit (IMEX) Runge-Kutta methods for differential problems in additive form have been investigated in recent literature [9,26].
In this manuscript, we consider a partitioned system of ordinary differential equations (ODEs) t ∈ [t 0 , t end ] , where the solution vector is separated into two components x and z. We assume that the first subsystem of (1) corresponding to the function f ∶ ℝ m × ℝ n → ℝ m is nonstiff, and the second subsystem of (1) corresponding to the function g ∶ ℝ m × ℝ n → ℝ n is stiff. Introducing the notation the system (1) can be written in a more compact form t ∈ [t 0 , t end ] , corresponding to the function F ∶ ℝ m+n → ℝ m+n .
The problem (1) will be integrated by the class of IMEX methods, where the first subsystem of (1) is integrated by explicit general linear method (GLM), and the second subsystem of (1) is integrated by diagonally implicit GLM. These methods for subsystems of (1) and the overall IMEX scheme, are defined in Sect. 2. In Sect. 3, we discuss stage order and order conditions for IMEX methods, and prove that the IMEX scheme has order p and stage order q = p if and only if the explicit and implicit "components" of the scheme have order p and stage order q = p . Stability properties of IMEX GLMs are investigated in Sect. 4. In Sect. 5, we define the SSP property of IMEX methods. A large class of transformed IMEX GLMs is defined in Sect. 6. In Sect. 7, we discuss the characterization of SSP IMEX methods, which is based on the recent work by Spijker [44]. Construction of transformed SSP IMEX GLMs is described in Sect. 8, and starting and finishing procedures are discussed in Sect. 9. The results of some numerical experiments with these methods are reported in Sect. 10. Finally, some concluding remarks and plans for future research are presented in Sect. 11.

Partitioned GLMs and IMEX Schemes
For the numerical solution of (1), we consider the partitioned method, where the first and the second subsystems of (1) are integrated by the explicit and diagonally implicit GLMs, respectively. These methods on the uniform grid are defined by n = 0, 1, ⋯ , N − 1 , and (1) x � (t) = f x(t), z(t) , x(t 0 ) = x 0 , z � (t) = g x(t), z(t) , z(t 0 ) = z 0 , (2) y � (t) = F y(t) , y(t 0 ) = y 0 , These methods are also characterized by four integers p, q, r, and s, where p is the order, q the stage order, r the number of external approximations, and s the number of internal approximations or stages. Introducing the notation (4) q k = q 1,k q 2,k ⋯ q r,k T ∈ ℝ r , k = 1, 2, ⋯ , q.
the methods (3) and (4) can be written in the vector form n = 0, 1, ⋯ , N − 1 , and n = 0, 1, ⋯ , N − 1 . Here, I is the identity matrix of dimension m, and Î is the identity matrix of dimension n. The relations (5) and (6) take now the form and We will reformulate next the methods (7) and (8) as GLM for the system (2). To simplify this reformulation, as well as the theoretical analysis of order and stage order conditions which is undertaken in the next section, and SSP properties discussed in Sects. 5 and 7, we will assume that the number of equations m of the nonstiff subsystem of (1) is equal to the number of equations n of the stiff subsystem. This assumption can be made without loss of generality, since if this is not the case, the subsystem with smaller number of equations can be padded by adding the equations of the form y � (t) = 0 , y (t 0 ) = 0 , or repeating some of the equations of this subsystem. Putting and taking into account that for any matrices Q and Q , we have where I and stand for identity matrices of dimension m, leads to Here, the abscissa vector and the coefficient matrices , , , and are defined by Observe that for this method, the relations (9) and (10) assume the form and where k are now matrices defined by

Stage Order and Order Conditions
In the remainder of the paper, we will be concerned with GLMs (11) of order p and stage order q = p . As in [49] it will be always assumed that the partitioned GLMs (7) and (8)  i approximate the solutions x and z to (1) at the same time points t n + c i h , i = 1, 2, ⋯ , s.
To derive order conditions for such methods, we consider local discretization errors (t n , h) and (t n , h) of internal and external stages of (11), which are defined by the relations and Here, and x(t n + ch) and z(t n + ch) are defined in Sect. 2. Then, the method (11) has order p and stage order q = p if and only if We have the following theorem.

Theorem 1
The internally consistent GLM (11) defined by abscissa vector with c =ĉ , and coefficient matrices , , , given by (12), and the matrices k , k = 0, 1, ⋯ , p , given by (15), has order p and stage order q = p if and only if the explicit GLM defined by c, A, U, B, V, and q k , k = 0, 1, ⋯ , p , and implicit GLM defined by c, Â , Û , B , V , and q k , k = 0, 1, ⋯ , p , have order p and stage order q = p.
Proof Let us partition the vectors (t n , h) and (t n , h) as where E (t n , h) , E (t n , h) correspond to the first m non-stiff components of (t n , h) , (t n , h) , and I (t n , h) , I (t n , h) correspond to the last n stiff components of (t n , h) , (t n , h) . Then, taking into account that y � (t n + h) = F(y(t n + h)) , we obtain and These relations imply that

and
We recognize E (t n , h) and E (t n , h) as local discretization errors of internal and external stages of explicit GLM (7), and I (t n , h) and I (t n , h) as local discretization errors of internal and external stages of implicit GLM (8). Hence, the explicit and implicit methods (7) and (8) have order p and stage order q = p , i.e., and if and only if the IMEX scheme (11) has order p and stage order q = p , i.e., These arguments complete the proof.
We can derive stage order and order conditions for explicit and implicit GLMs (7) and (8) by expanding x(t n + ch) , x � (t n + ch) , x (k) (t n+1 ) , and z(t n + ch) , z � (t n + ch) , z (k) (t n+1 ) , appearing in E (t n , h) , E (t n , h) , and I (t n , h) , I (t n , h) , into Taylor series around t n , and comparing to zero the terms with powers h k , k = 0, 1, ⋯ , p . Such stage order and order conditions were derived in [3,4] and in the monograph [36]. To make this paper self-contained, these stage order and order conditions are also listed in Theorem 2 below. (7) and (8)

3
For k = 0 , we obtain the pre-consistency conditions where e = [1, ⋯ , 1] T ∈ ℝ r . For the methods investigated in this paper, the pre-consistency vectors q 0 and q 0 are given by It follows from Theorem 1 that the IMEX scheme (11) has order p and stage order q = p if the "component" GLMs (7) and (8) independently meet the conditions (18), (19), (20), and (21), for order p and stage order q = p . Contrary to partitioned Runge-Kutta methods, considered for example in [21], no coupling conditions are required, i.e., conditions involving parameters of both methods. This result was also proved, using somewhat different arguments, in [49], for internally consistent partitioned GLMs of order p and stage order q = p or q = p − 1.

Stability Properties of IMEX Schemes
We consider first the decoupled linear test system t ⩾ 0 , where is a complex parameter. Then applying IMEX scheme (11)  Hence, it follows that the absolute stability properties of IMEX method with respect to (22) are combinations of absolute stability properties of explicit and implicit parts of the scheme. In this paper we will search for IMEX schemes, where the explicit and implicit GLMs have so-called inherent Runge-Kutta stability (IRKS) [5,36,47,48], and the Uq 0 = e, Vq 0 = q 0 ,Ûq 0 = e,Vq 0 =q 0 , implicit GLM is A-and L-stable. Moreover, we will monitor the size of the region of absolute stability of the explicit GLM.
We consider next the coupled linear test system of the form t ⩾ 0 , where a and b are real parameters. The solutions (x(t), z(t)) tend to zero as t tends to the infinity if and only if and we will investigate whether this property is inherited by the numerical approximations (x [n] , z [n] ) obtained by application of IMEX method (11) to the system (23). Applying (11) to (23)  Observe that for = 0 we have where M E ( ) and M I ( ) are stability functions of explicit and implicit GLMs, respectively, with respect to (22) corresponding to Re(z) = . We also define stability polynomial IMEX (w; , ) with respect to the test system (23) by the relation and the stability region with respect to (23) is then given by We recall that IMEX (w; , ) is a Schur polynomial if all its roots w i , i = 1, 2, ⋯ , 2s , are inside of the unit circle. Stability regions A IMEX defined by (29), with respect to the test system (23), for SSP IMEX GLMs of order p = 1 , 2, 3, and 4, and stage order q = p , are displayed in Sect. 8.
The characterization of SSP IMEX methods is investigated in Sect. 7 for the large class of IMEX transformed GLMs (TGLMs) defined in Sect. 6

IMEX TGLMs
Similarly as in [34,35], a class of IMEX TGLMs is defined by multiplying the relation for y [n+1] in (11)  Observe that these matrices have the following block structure: Then the method (35) can be rewritten in the following form: Similarly as in [34,35], we can show that the IMEX TGLMs (36) preserve the order and stage order of the original methods (11), and have the same stability properties with respect to the basic test system (22).
To investigate stability properties of IMEX TGLMs with respect to the test system (23), we examine the stability matrix * IMEX ( , ) corresponding to this scheme. We have and taking into account the relation it follows that This relation implies that the stability properties of the transformed IMEX schemes with respect to the linear test system (23) are the same as that of the original IMEX methods. However, in general, SSP properties of transformed methods are different from SSP properties of original methods, and in Sect. 8, we will search for IMEX schemes with maximal SSP coefficients in a large class of IMEX TGLMs. * =

Characterization of SSP IMEX TGLMs
To investigate SSP properties of IMEX TGLMs (36), we will first reformulate these methods using the notation introduced by Spijker [44]. Define the coefficient matrices Then (11) can be written in the form In what follows, we will always assume that p = q and r = s . It follows from the preconsistency conditions that and the coefficient matrix * satisfies as required in [44]. We describe next the characterization of SSP coefficient C = C( * , * ) of GLMs (37) which was discovered by Spijker in his seminal paper [44]. Denote by the identity matrix of dimension 2(s + r) × 2(s + r) and by the matrix whose first 2r columns are equal to * and the last 2(s + r) columns are equal to * . Following [44] consider the condition where the inequality in (38) should be interpreted componentwise. Then it was demonstrated by Spijker [44] that the SSP coefficient of GLM (37) is given by * = * * ∈ ℝ 2(s+r)×2(s+r) , * = * * ∈ ℝ 2(s+r)×2r .

3
This condition can be reformulated in terms of coefficient matrices * , * , * , and * of TGLMs (37). Since where 2s and 2r are identity matrices of dimensions 2s × 2s and 2r × 2r , it follows from the structure of the coefficient matrix * that the condition det( + * ) ≠ 0 is automatically satisfied. It can be also verified, compare [31,32], that the second condition of (38) is equivalent to the following inequalities of lower dimensions: Hence, the SSP coefficient of GLM (11), which will be denoted now by C IMEX , can be expressed in terms of * , * , * , and * by the formula These inequalities, reformulated in terms of coefficient matrices of explicit and diagonally implicit TGLMs (36), take the form and where I s is the identity matrix of dimension s × s , and the SSP coefficients of explicit and diagonally implicit TGLMs are given by and It follows from the above discussion that the IMEX method (36) is SSP if and only if the explicit and diagonally implicit TGLMs are SSP, and we have

Construction of SSP IMEX GLMs
In this section, we describe the construction of internally consistent SSP IMEX TGLMs up to the order p = 4 and stage order q = p . We assume that the vectors q k and q k , where e 1 , e 2 , ⋯ , e p+1 are the canonical bases in ℝ p+1 . For such methods, the vectors y [n] approximate the Nordsieck vectors N y (t n , h) , n = 0, 1, ⋯ , N , defined by It follows from the discussion in Sect. 7, in particular from the formulas (45) and (44), that implicit parts of these schemes can be obtained by solving the minimization problem with non-linear inequality constrains (43), and the implicit parts of these schemes can be obtained by solving the minimization problem with non-linear inequality constrains (42). As already mentioned in Sect. 4, we will always assume that both implicit and explicit methods have IRKS property introduced in [5,47] (see also [36,48]), and that implicit methods are A-and L-stable. Such methods are obtained following the algorithm for construction of formulas with IRKS presented in [5,36,47].

SSP IMEX TGLMs with p = q = 1 and r = s = 2
If the abscissa vector is fixed a priori, the explicit and implicit methods can be constructed independently. Here, we show the results corresponding to the abscissa vector = [1∕2, 1] T . Consider first the implicit methods with r = s = 2 . Following the algorithm for construction of GLMs with IRKS, and assuming that the diagonal element of the coefficient matrix Â is equal to 1/4, we obtain an A-and L-stable family of methods with p = q = 1 depending on the parameter ̂ 1 , appearing in the algorithm for construction of such methods. Then solving the minimization problem (46) with inequality constrains (43) with respect to ̂ 1 and coefficients of the transformation matrix T defined in Sect. 6 we obtain the method with coefficients The transformation matrix T is For this method, C I = 2 and C I,eff = 1.
Consider next the explicit methods with r = s = 2 . Following again the algorithm for construction of GLMs with IRKS presented in [5,36,47], we obtain a family of methods with p = q = 1 depending on the parameter 1 appearing in the algorithm for construction of such methods. Then, assuming that the method is internally consistent with the implicit formula, and solving the minimization problem (47) with inequality constrains (42), with respect to 1 , and coefficients of the transformation matrix T defined in Sect. 6, we obtain the method with coefficients The transformation matrix T isT  Fig. 4 stability region A E of the numerical scheme, where both equations of (23) are treated by the explicit GLM. As before, we have also plotted by a thick line the boundary = − , ⩽ 0 , of stability region of the test system (23). We can see that this region is bounded, and the region A IMEX is unbounded for ⩽ 0 , and ≈ or ≈ − . Observe that stability intervals on Figs. 1, 2, 3, 4 are the same. This is a consequence of the fact that for b = 0 or = 0 the test system (23) reduces to (22) with Re( ) = a or Re(z) = . If the abscissa vector is fixed a priori, the explicit and implicit methods can be constructed independently. Here, we show the results corresponding to the linearly spaced abscissa vector = [0, 1∕2, 1] T . Consider first the implicit methods with r = s = 3 . Following the algorithm for construction of GLMs with IRKS presented in [5,36,47], and assuming that = 0 , we obtain a possibly A-and L-stable family of methods with p = q = 2 depending on , ̂ 1 , ̂ 2 , and l 21 . Then, solving the minimization problem (46) with inequality constrains (43) and coefficients of the transformation matrix T defined in Sect. 6, we obtain the method with coefficients The transformation matrix T is For this method C I = 2.131 and C I,eff = 0.710.
Consider next the explicit methods with r = s = 3 . Following again the algorithm for construction of GLMs with IRKS presented in [5,36,47], assuming that the diagonal element of the coefficient matrix Â is equal to 1/4 and that = 1∕6 , we obtain a family of methods with p = q = 2 depending on the parameters 1 and 2 , l 21 , appearing in the algorithm for construction of such methods. Then, solving the minimization problem (47) with inequality constrains (42), with respect to 1 , 2 , l 21 , and coefficients of the transformation matrix T defined in Sect. 6 we obtain the method with coefficientŝ  The SSP coefficients for the resulting IMEX scheme are C IMEX = 1.193 and C IMEX,eff = 0.398 . Stability region A IMEX of this IMEX scheme is plotted in Fig. 6. We have also plotted on this figure by a thick line the boundary = − , ⩽ 0 , of stability region of the test system (23). We have also plotted on

SSP IMEX GLMs with p = q = 3 and r = s = 4
Here, we show the construction of methods with linearly spaced abscissa vector = [0, 1∕3, 2∕3, 1] T . Consider first the implicit methods with r = s = 4 . Following the algorithm for construction of GLMs with IRKS, and assuming = 0 , we obtain a possibly A-and L-stable family of methods with p = q = 3 depending on , and additional parameters ̂ 1 , ̂ 2 , ̂ 3 , and l 21 , l 31 , l 32 . Then solving the minimization problem (46) with inequality constrains (43), additional The transformation matrix T is For this method C I = 1.51 and C I,eff = 0.38. Consider next the explicit methods with r = s = 4 . Following again the algorithm for construction of GLMs with IRKS, and assuming = 1∕24 , we obtain a family of methods with p = q = 3 depending on the parameters 1 and 2 , 3 , l 21 , l 31 , and l 32 appearing in the algorithm for construction of such methods. Then, assuming that the method is internally consistent with the implicit formula, and solving the minimization problem (47) with inequality constrains (42), with respect to 1 , 2 , 3 , l 21 , l 31 , l 32 , and coefficients of the transformation matrix T defined in Sect. 6, we obtain the method with coefficientŝ constrains ⩽ 3∕2 , with respect to , ̂ 1 , ̂ 2 , ̂ 3 , ̂ 4 , l 21 , l 31 , l 32 , l 41 , l 42 , l 43 , and coefficients of the transformation matrix T defined in Sect. 6 we obtain the method with coefficients The transformation matrix T is For this method C I = 1.50 and C I,eff = 0.30. Consider next the explicit methods with r = s = 5 . Following again the algorithm for construction of GLMs with IRKS, assuming = 1∕120 , we obtain a family of methods with p = q = 4 depending on the parameters 1 and 2 , 3 , 4 , l 21 , l 31 , l 32 , l 41 , l 42 , l 43 , appearing in the algorithm for construction of such methods. Then, assuming that the method is internally consistent with the implicit formula, and solving the minimization problem (47) with inequality constrains (42), with respect to 1 , 2 , 3 , 4 , l 21 , l 31 , l 32 , l 41 , l 42 , l 43 , and coefficients of the transformation matrix T defined in Sect. 6 we obtain the method with coefficientŝ          For this method C E = 0.63 , C E,eff = 0.13 , and the region of absolute stability is plotted in Fig 13. The SSP coefficients for the resulting IMEX scheme are C IMEX = 0.63 and C IMEX,eff = 0.13 . Stability region A IMEX of the resulting IMEX scheme is plotted in Fig. 14. This region is bounded and extends to about = ± = −6.00 along the lines = ± . We have also plotted on this figure by a thick line the boundary = − , ⩽ 0 , of stability region of the test system (23). We have also plotted on Fig. 15   The line corresponding to ( IMEX ) = 1 is plotted by a thin dash-dot line. For comparison we have also plotted on Fig. 16 stability region A E of the numerical scheme, where both equations of (23) are treated by the explicit GLM. As before, we have also plotted by a thick line the boundary = − , ⩽ 0 , of stability region of the test system (23).

Starting and Finishing Procedures
To start the integration with IMEX method (36), we have to compute the starting vectors where This can be done in many different ways and in what follows, we describe the strategy to compute x [0] and z [0] using initial values x 0 = x(t 0 ) and z 0 = z(t 0 ) and approximations x k and z k of order p to x(t 0 + kh) and z(t 0 + kh) , k = 1, 2, ⋯ , p . This strategy was first suggested in the context of diagonally implicit multistage integration methods by Butcher [3], and then described in more detail in [6,7] in the context of GLMs. Putting we have where N x (t, h) and N z (t, h) are Nordsieck vectors corresponding to x(t) and z(t), i.e., For methods with r = p + 1 , these Nordsieck vectors can be approximated by where the matrix S can be obtained by expanding x(t 0 + kh) and z(t 0 + kh) into Taylor series around t 0 and comparing terms of the same order on both sides of the resulting relations. The matrices S corresponding to p = 1 , 2, 3, and 4, are listed in [6,7] but to make this paper self-contained these matrices are also listed below. We have for p = 1 and p = 2 , and for p = 3 and p = 4.
We describe now the finishing procedure following the presentation in [7]. After completing the integration with IMEX GLMs (3) and (4) 1 are the required approximations to x(t end ) and z(t end ) , i.e.,   x � (t) = z(t),  t ∈ [0, t end ] , t end = 5 . The initial conditions are Example 4 The differential system from [27]:  where h is the water height with respect to the bottom and hv is the flux of the velocity field. We use periodic boundary conditions and initial conditions at t 0 = 0 with x ∈ [0, 1] . For this problem the space derivative has been discretized by a fifth order finite difference weighted essentially non-oscillatory (WENO) scheme following the implementation described in [42]. The evolution of the model is shown in Figs. 21 and 22.
To verify the experimental order of the proposed methods, we considered the case with N = 100 spatial points, that results in a system of first-order ODEs of dimension d = 200 . The first hundred components, corresponding to h t , are integrated by the explicit method, while the rest of the components, corresponding to (hv) t , are integrated by the implicit method.
The numerical results with = 10 −2 and = 10 −8 , reported in Figs. 23 and 24, confirm that the methods proposed in this paper have the asymptotic preserving property, but they are also asymptotically accurate in the stiff limit for → 0 (compare [40]). In other words, for small values of the expected order of convergence matches the theoretical one and no order reduction occurs (compare [33,35]).
We consider the evolution of a population density P according to the model where x ∈ [0, 1] and t ⩾ 0 . We assume periodic boundary conditions and zero initial condition, P(0, x) = 0 , x ∈ [0, 1]. This differential equation models the evolution of an ecological population using birth, death and migration. A description of the model can be found in [29]. The birth rate, b(x, P), has taken to decline nonlinearly with the population density: The death rate is set equal to one, i.e., r d = 1 . The birth rate r b varies by position according to The forcing term, f(t, x), is taken to be zero for all times t ≠ 0 . At time t = 0 , and for each grid point, f is assigned a random value in the interval [0.8, 1.2] according to a uniform distribution. Finally, note that we treat models with diffusion, d = 0.01, 0.04 , and without, d = 0 . Without diffusion, the population eventually must tend to zero in the first half of the domain. Conversely, the diffusive model evolves as shown in Fig. 25 and tends to a nonzero profile as t increases, see Fig. 26. In all cases, the analytical solution is non-negative for all times t ≠ 0. For these tests, we set Dx = 1∕100 and use second-order centered differences to discretize the diffusion term. The problem is reformulated as (1), where stiff diffusion term is treated implicitly, while the other terms are treated explicitly.
As stated in [29], the preservation of positivity is particularly desirable in this problem since it leads to more biologically meaningful population densities and avoids the possibility of having reaction terms grow unboundedly (which would happen for P tending to − ). Examples of IMEX methods which, for some choice of the stepsize h, return negative approximations that cause huge error can be found in [29].
We have experimentally determined the largest step size h for which the numerical solution maintains the positivity (at least) up to t = 10 . The results reported in Table 1 show that the proposed IMEX-GLMs methods preserve the positivity of the numerical solution for larger stepsizes with respect to methods analyzed in [29].

Concluding Remarks
We considered the integration of partitioned systems of ODEs, with stiff and nonstiff components, by IMEX general linear methods. For this class of methods we provide formulation, order and stage order conditions, as well as a characterization of SSP IMEX GLMs. We then analyzed the stability of these methods with respect to two, uncoupled and coupled, linear test systems, and we introduced a suitable transformation which does not modify the order, the stage order and the linear stability properties of the methods, but allow to obtain methods with larger SSP coefficients. Afterwards, we showed how to determine good methods inside the considered class, and provided examples of methods of order p = 1, 2, 3 and p = 4 . Finally, after giving information about the starting and the finishing procedure, we reported the results of numerical tests on a selection of stiff problems, that confirm the good performance of the methods here introduced.
Funding Open access funding provided by Università degli Studi di Napoli Federico II within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.