Numerical continuation for fractional PDEs: sharp teeth and bloated snakes

Partial differential equations (PDEs) involving fractional Laplace operators have been increasingly used to model non-local diffusion processes and are actively investigated using both analytical and numerical approaches. The purpose of this work is to study the effects of the spectral fractional Laplacian on the bifurcation structure of reaction-diffusion systems on bounded domains. In order to do this we use advanced numerical continuation techniques to compute the solution branches. Since current available continuation packages only support systems involving the standard Laplacian, we first extend the pde2path software to treat fractional PDEs. The new capabilities are then applied to the study of the Allen-Cahn equation, the Swift-Hohenberg equation and the Schnakenberg system (in which the standard Laplacian is each replaced by the spectral fractional Laplacian). Our study reveals some common effects, which contributes to a better understanding of fractional diffusion in generic reaction-diffusion systems. In particular, we investigate the changes in snaking bifurcation diagrams and also study the spatial structure of non-trivial steady states upon variation of the order of the fractional Laplacian. Our results show that the fractional order can induce very significant qualitative and quantitative changes in global bifurcation structures.


Introduction
In recent years increasing research attention has been devoted to the study of partial differential equations (PDEs) involving so called fractional Laplace operators, which generalize the notion of derivative of non-integer orders. Their study is not only interesting from a purely mathematical point of view, but has proven extremely useful for modeling super-diffusion processes, which naturally appear in many applications in physics, probability, biology, ecology, medicine and economics [39]. Examples include models for complex phase transitions [50], for chemical transport in heterogeneous aquifers [1], for pattern formation in coral reefs [51], for larval growth and recruitment in a turbulent environment [14], for foraging of marine predators [56], for cardiac electrical propagation [12] and for mechanisms in options pricing [34].
From a mathematical point of view, the fractional Laplacian does not have a unique definition. On unbounded domains there are several equivalent definitions (e.g. Fourier, singular integral, Balakrishnan formula, harmonic extension) [33]. On bounded domains the definitions are usually not equivalent (e.g. spectral, Riesz). In fact, there are several different associated microscopic stochastic processes from which on may derive a fractional Laplacian. On bounded domains, boundary conditions must be incorporated in these characterizations in mathematically distinct ways and there is currently no consensus in the literature. For a detailed review we refer to [35,19]. It is widely accepted that the fractional Laplacian constitutes the counterpart in the theory of non-local operators of the standard Laplacian in the classical theory of partial differential equations. In particular, the fractional Laplacian ∆ s is characterized by the parameter s ∈ (0, 1). Letting s → 0 + one obtains the identity, and letting s → 1 − one gets the classical Laplacian. Despite the difference of the two operators, most of the relevant questions associated to the Laplacian (e.g. comparison principles, linear and nonlinear boundary value problems, and regularity results) have an equivalent for the fractional Laplacian, and many techniques for solving partial differential equations also apply to fractional PDEs.
From a dynamical systems perspective, one of the most widely studied class of PDEs are standard reaction-diffusion problems where t ∈ R, x ∈ Ω ⊆ R n are the time and space variables, u ∈ R N is the vector of unknowns, the matrix D ∈ R N × R N contains the diffusion coefficients, ∆ is the classical Laplacian operator, F : R N → R N corresponds to non-linear reaction terms, and µ ∈ R is a parameter. In many cases a natural starting point in the study of such equations is to consider the non-linear elliptic problem which yields steady states. Of particular interest in applications is the question, how the steady states change under parameter variation. To this end, several methods (Lyapunov-Schmidt reduction, center manifolds, amplitude equations, Turing instability) can help to address this question analytically near trivial branches of homogeneous steady states [30]. In addition to theoretical methods and criteria, the global bifurcation structure of steady states reveals the behavior of solutions far from homogeneous solutions or the critical points. For classical reaction-diffusion systems on bounded domains (and n = 1, 2, 3), the global structure can be numerically computed with, for instance, pde2path [54]. This packacke is an advanced continuation/bifurcation software based on the FEM discretization of the stationary elliptic problem exploiting the package OOPDE [43] for the FEM discretization. Since pde2path allows branch point continuation and Hopf point continuation, continuation of relative equilibria (e.g., traveling waves and rotating waves), branch switching from periodic orbits (Hopf pitchfork/transcritical bifurcation, and period doubling), it is quite flexible. Recently, it has also been used to treat cross-diffusion systems [31,11]. Yet, non-standard diffusion problems are not within the classical realm of pde2path applications.
In this framework the aim of the paper is two-fold. First, we want to extend the continuation software pde2path which has been extensively used for classical non-linear reaction-diffusion equations, to treat fractional equations, i.e. equations involving the fractional Laplacian. Then, we test the new capabilities of the software onto three important benchmark problems: namely the Allen-Cahn equation [4], the Swift-Hohenberg equation [17] and the Schnakenberg system [46] on bounded domains in a non-local setting. All three equations have theoretical interest, concrete applications, and are currently widely used in various communities working on PDEs. The purpose is to qualitatively and quantitatively explore the features of these models in the fractional setting and to better understand the effect of fractional diffusion on the steady state bifurcation structure of generic fractional reactiondiffusion systems.
Let Ω ⊂ R n (n = 1, 2, 3) be a bounded domain with boundary ∂Ω and s ∈ (0, 1) be the fractional order. A generic fractional reaction-diffusion systems can be formulated as ∂ t u = −D(−∆) s u + F (u; µ) =: D∆ s u + F (u; µ), endowed with suitable boundary conditions and where we have introduced the shorthand notation ∆ s := −(−∆) s to be used throughout this work. In this paper we consider the spectral definition of the fractional Laplacian ∆ s [15], that is, for any ω ∈ C ∞ (Ω), where λ j and φ j the eigenvalues and the eigenfunctions, respectively, of the standard Laplacian ∆ on Ω, i.e., the solutions of the eigenvalue problem ∆φ j = λ j φ j with specific boundary conditions. In order to numerically treat fractional PDEs and to embed it naturally into the numerical continuation approaches based upon FEM, we have to provide a FEM discretization of the fractional operator. The approach we employ to represent (2) without a spectral Galerkin decomposition is to use the Balakrishnan representation formula [8], namely for s ∈ (0, 1) and u ∈ H −s (Ω) The numerical approximation of (3) is based on a sinc quadrature approximation of the additionally involved integral with respect to ν and a discretization of the operator νI − ∆ using the finite element method in space [10]. This strategy has been successfully applied recently in other contexts such as optimal control problems [22]. Here we integrate this discretization within the numerical continuation. Then, using the new capabilities of the continuation software, we study the effect of fractional derivative order on the steady state bifurcation structure of three well known PDEs in which the standard Laplacian is replaced by a spectral fractional Laplacian: -The fractional Allen-Cahn equation with cubic-quintic nonlinearities will be considered on a bounded interval Ω ⊂ R taken with homogenous Dirichlet boundary conditions. The model with the standard Laplacian, proposed in [4], is a well-studied equation for various polynomial nonlinearities [3,2,44,57]. In particular, the classical Allen-Cahn PDE (s = 1) with a cubic nonlinearity is also known as the (real) Ginzburg-Landau PDE [6], an amplitude equation or normal form for bifurcations from homogeneous states [40,48]. Moreover, it is also common to consider an additional quintic term added to the Ginzburg-Landau PDE [26,29].
-The fractional Swift-Hohenberg equation with competing cubic-quintic nonlinearities will also be studied on a bounded interval Ω ⊂ R, again taken with homogeneous Dirichlet boundary conditions. The classical Swift-Hohenberg equation (s = 1) is a widely studied model in pattern dynamics [17,47,52]. It is also a standard test case for deriving reduced amplitude/modulation equations [30,49,27,16,55]. Furthermore, it was found that global bifurcation diagrams of the Swift-Hohenberg equation exhibit a process called (homoclinic) snaking [13,7,25]. Recently, also the Swift-Hohenberg equation with nonlocal reaction terms has gained particular attention [32,41] but results for bifurcations in the space-fractional Swift-Hohenberg equation do not seem to be available up to now.
-The fractional Schnakenberg system is given by where the reaction part is As before, we use a bounded interval Ω ∈ R taken with homogeneous Neumann boundary conditions. Note that with σ = 0 the reaction part (7) reduces to the classical formulation of the Schnakenberg system [46], which is one of the prototype reaction-diffusion systems exhibiting Turing patterns [42]. A nonzero σ modifies the primary bifurcation from super-to subcritical, leading to to the appearance of snaking branches between periodic branches [53,54].
As indicated above, the classical Swift-Hohenberg equation and the Schnakenberg system both exhibit snaking behaviour [13,7,25,9,28,37,38,53]. Upon variation of the main bifurcation parameter, a non-trivial branch of solutions emerges from a branch point and later on undergoes several fold bifurcations leading to a snaking-type structure of the bifurcation diagram and to multistability. Since solutions on these branches are often important localized patterns, a detailed analysis is required. Generally it is important to understand, how the snaking scenario changes under nonlocal influence as already carried out for nonlocal reaction terms in [41]. Here we show that snaking structure in the Swift-Hohenberg equation gets "bloated" upon decreasing the fractional power s ∈ (0, 1). Furthermore, we study the deformation of the localized solutions on the nontrivial bifurcation branches. In the Schnakenberg system, we observe a similar bloating effect and further classes of deformations of the branches. The solutions themselves tend to sharpen in the transition layers, which is an effect already found for non-trivial branches of the fractional Allen-Cahn equation leading to "sharpening of the teeth" of spatially oscillating profiles.
The paper is organized as follows. In Section 2, we briefly recall the standard version of the benchmark problems (Allen-Cahn and Swift-Hohenberg equations, and the Schnakenberg system), including standard steady state bifurcation diagrams as well as typical solutions. The software pde2path as well as the discretization and the implementation of the fractional operator are then detailed in Section 3, followed by some preliminary results in order to validate the discretization method. In Section 4 we present the results obtained exploiting the new continuation capabilities, highlighting the effects of fractional diffusion on the bifurcation structure of the considered problems. Finally, in Section 5 some concluding remarks can be found. The Matlab code is reported in the Appendix, together with some explanations on the Matlab implementation. Some supplementary videos are available at [23].

Bifurcation diagrams with standard Laplacian
In this section we briefly recall the standard version of the benchmark problems which will be investigated in the fractional version. Standard steady state bifurcation diagrams are presented as well as typical solutions for a better comparison with the fractional versions. They are obtained using the continuation software pde2path, and they were already shown in [53,54,45].
We use several conventions for the figures. Hereafter, thicker lines in the bifurcation diagrams denote stable solutions, while we use thinner lines for unstable ones. Circles, crosses and diamonds indicate the presence of branch points, fold points and Hopf points, respectively.

The Allen-Cahn equation
The Allen-Cahn equation was proposed to model phase separation in alloy systems with multiple components. With standard diffusion and cubic-quintic non-linearity it reads where u = u(x, t) ∈ R is the unknown and µ, γ ∈ R are parameters. On a one-dimensional bounded domain Ω of length L, considering homogeneous Dirichlet boundary conditions,ū = 0 is a homogeneous stationary solution of equation (8) for all choices of µ, γ. Further, linearizing equation (8) aroundū, it can be shown that this state is stable as long as µ < (π/L) 2 , and first becomes unstable to spatially periodic modes with wavenumber k c = π/L when the parameter reaches the critical value µ c = (π/L) 2 . Subsequent modes become unstable, at µ j = (jπ/L) 2 , j ≥ 2.
The steady state bifurcation structure of equation (8) can be computed numerically to obtain a more global picture far from the homogeneous branch. We consider Ω = [−5, 5], i.e. L = 10, and set γ = 1. The branches of spatially periodic patterns bifurcating from the homogeneous branch are illustrated in Figure 1. They bifurcate via transcritical bifurcations and experience a fold at larger amplitude. The first bifurcating branch is unstable until the fold and becomes stable after the fold. All other branches are unstable for µ < 2. In addition, Figures 1a-1b show solution profiles along the second branch (originating at the second bifurcation point B 2 ) for increasing vales of µ: spatial subregions appear occupied by one of two phases between which there is an interface which sharpens as µ increases.

The Swift-Hohenberg equation
The Swift-Hohenberg equation is one of the simplest phenomenological model for pattern formation. With cubic-quintic non-linearity it reads with parameter λ ∈ R and ν > 0. This equation presents the following characteristics: it is reversible in space, i.e. it is equivalent under the transformation x → −x, u → u, and, due to the choice on nonlinearity, is also symmetric under the transformation x → x, u → −u. We consider the one-dimensional bounded domain Ω of length L = mπ, m ∈ N >0 with homogeneous Dirichlet boundary conditions. In this case the trivial stateū = 0 is the only spatially homogeneous steady state of equation (9). Linearizing equation (9) aroundū, it can be shown that the trivial state is stable as long as µ < 0 and first becomes unstable to spatially periodic modes with wavenumber k c = 1 at µ c = 0. On the bounded domain Ω further bifurcations occur at the parameter values µ j = (1 − (jπ/L) 2 ) 2 , j ≥ 1.
As in the previous section, the bifurcation structure can be obtained numerically and is shown in Figure 2 for Ω = [−5π, 5π] and ν = 2. The (spatially) periodic patterns arise on the blue branch bifurcating subcritically from the homogeneous (black) branch (as shown in [13, equation 9] for ν > 0). The subcriticality of the bifurcation is of fundamental importance for the appearance of the "snaking" (red) branch, which emerges on bounded domains from the branch of spatially periodic states close to the origin at µ < 0. This branch corresponds to front solutions shown in Figures 2a-2d. At each turn in the snaking, one additional oscillation is added to the pattern until the domain is filled with oscillations. Then, the snaking branch reconnects close to the fold of the periodic branch.
Note that additional bifurcations further up the periodic branch lead to other snaking branches corresponding to front or localized patterns. These branches are not illustrated here for the sake of clarity.

The Schnakenberg system
Finally, the Schnakenberg system has been used to model the spatial distribution of a morphogen for tissue patterning, and it presents the activator-inhibitor structure crucial in Turing instability. The system reads  (7), on a bounded domain with length L = 10π/k c , with homogeneous Neumann boundary conditions. For sake of clarity, only the first three branches bifurcating from the homogeneous states (black) are shown.
where the reaction part F is given in (7), and also depends on the parameter σ. The homogeneous stateū = (µ, 1/µ) is a steady state solution of system (10), independent of the value of σ in the reaction part. Further, following [20], we find thatū is stable for parameter values µ > µ c and that at µ c spatially periodic states bifurcate from the homogeneous branch with critical wavenumber k c , where Note that both critical values are independent of σ. On a finite domain, the first bifurcation is followed by a discrete number of bifurcations in µ < µ c as further modes become unstable. The bifurcation diagram for the steady states of (10) with σ = 0 is shown in Figure 3, considering the domain Ω = [−5π/k c , 5π/k c ] (chosen to accommodate the basic periodic patterns bifurcating from the homogeneous branch). The spatially periodic (blue) branch bifurcates supercritically at µ c from the homogeneous (black) branch. Two further branches are shown in red and magenta.
As mentioned above, µ c and k c are independent of σ. This is because σ does not affect the linearization. However, it affects the nonlinear terms in (10) and hence the restabilization of spatially periodic modes at higher amplitudes [53]. Thus, we can use σ to tune the primary bifurcation from super-to subcritical. This enables us to "create" snaking branches between periodic branches, as is illustrated in Figure 4. As for the Swift-Hohenberg equation, the snaking branch corresponds to front solutions illustrated in Figure 4a

Numerical continuation of fractional PDEs in pde2path
As mentioned in the introduction, we want to adapt numerical continuation in the context of the Matlab bifurcation package pde2path to treat fractional reaction-diffusion equations. In this section, we briefly present the software in its standard setting, that is for PDEs with standard Laplacian. Then we explain the FEM discretization of the fractional Laplacian and modifications made to the FEM discretization in pde2path in order to study fractional systems. We finally present numerical tests in order to validate our method. The purpose of this section is not to give a complete overview of the software for beginner users, hence we do not explain in detail the basic setup; see [21,45,54] for complete guides and for the notation adopted in the following.

Standard setting
The continuation software pde2path can deal with the general class of PDEs of the form and Ω a bounded domain, µ ∈ R p is a parameter vector, and the diffusion, c, advection, b, and linear tensors, a, as well as the nonlinearity f can depend on x, u, ∇u and parameters. The problem is endowed with generalized Neumann boundary conditions where n is the outer normal to the boundary, q ∈ R N ×N and g ∈ R N . This formulation covers Neumann and Robin boundary conditions, as well as Dirichlet boundary conditions via large stiff spring pre-factors in q and g. The software spatially discretizes the problem via the finite element method (FEM) exploiting the OOPDE toolbox [43], leading to a high-dimensional ODE problem. Note that a code extension is also available for periodic boundary conditions [21] and cross-diffusion terms [31].
In order to study the stationary problem related to equation (1), we restrict this introduction to the case a = 0, b = 0, c ∈ R N ×N ×n×n such that c ijkh = d i δ ij δ kh , i, j = 1, . . . , N, k, h = 1, . . . , n, and µ ∈ R, leading to where D ∈ R N ×N is the diagonal matrix of the diffusion coefficients d i . In pde2path the problem (12), together with boundary condition (11), is converted into the algebraic system where u ∈ R npN contains the nodal values of u with n p mesh points. The matrix K tot ∈ R npN ×npN corresponds to the finite element discretization of the D and q terms from equations (12) and (11), while F tot ∈ R npN corresponds to the discretization of f and g in (12) and (11). Precisely, where K is the stiffness matrix and corresponds to D, s BC is the stiff spring pre-factor use to approximate Dirichelet boundary conditions (BC) and F, Q BC , G BC corresponds to f, q, g respectively. Details on the discretization for different types of boundary conditions can be found in [45].

Fractional PDEs in pde2path
In order to adapt the continuation software pde2path to treat fractional reaction-diffusion equations (1) endowed with homogeneous Neumann or Dirichlet boundary conditions, we need to discretize the spectral fractional Laplacian (2) keeping the structure of equation (13). We exploit the Balakrishnan formula representation [8,33] in the version appearing in [58, formula (4), p.260], reported in equation (3). Note that "Balakrishnan formula" refers to both the direct fractional Laplacian representation, ∆ s u as above, and the inverse fractional Laplacian representation, ∆ −s u see [22,10].
In [22,10] the inverse formula has been discretized via a sinc quadrature for the integral as well as a finite element discretization in space. Here we adapt this method to the direct formula (3). First, the substitution ν = e η (dν = e η dη) in (3) leads to Then using sinc quadrature [36] one obtains which, choosing suitable values of n + and n − , can be approximated by where v l ∈ U(τ ) 1 is the solution of the Galerkin variational problem approximated via FEM as the solution of the linear system Algorithm 1: Algorithm that allows to build the matrix K s , which corresponds to a FEM discretization of the fractional Laplacian. A detailed description of the implementation of this algorithm in pde2path is given in Appendix A.
where M and K are the classical FEM mass and stiffness matrices respectively, and v l denotes the node vector for v l . Hence, we obtained an approximation of ∆ s u, which requires to solve n − + n + + 1 independent linear systems. However, to fit the pde2path setup (which necessitates the gradient expression of G in (13)) given by K tot − ∇F tot ), we need a matrix approximation of the operator ∆ s itself. Such matrix, called K s , can be built column by column by testing the proposed approximation method on the standard basis of R np of vectorsê i , 1 ≤ i ≤ n p 2 . In detail, the i-th column of K s is obtained computing −∆ sê i with the proposed scheme. Once the matrix K s has been built, we end up with an algebraic system of the form (13) by replacing the standard stiffness matrix K with M K s in the expression of K tot in formula (14).
The complete method is presented in Algorithm 1. A description of its implementation in pde2path can be found in Appendix A. In particular, the Matlab code is given in Listing 1 and the choice of the quantities k, n + , n − as well as further implementation issues are discussed.

Numerical tests
In order to validate the algorithm presented in the previous section, we first present some numerical results on the eigenvalues of the fractional Laplacian. The algorithm is then applied to a simple Poisson problem to exemplify the effect of fractional diffusion.

Spectrum of the discretized fractional Laplacian
A first simple test consists in checking the convergence of the eigenvalues of the discretized fractional operator K s towards those of the continuous operator as the mesh-size is decreased. In particular, the eigenvalues of the (continuous) spectral fractional Laplacian (−∆) s on a one-dimensional interval Ω of length L are known, both considering homogeneous Neumann and Dirichelet boundary conditions. For instance, with homogeneous Neumann boundary conditions they are given by λ j = (jπ/L) 2s for j ≥ 0. The first n e eigenvalues λ (h) j , j = 1, . . . , n e of the discretized fractional Laplacian K s can be numerically computed, and the maximum relative error can be used as a measure of the convergence We choose a domain Ω = [0, 1] of length L = 1 and we vary the number of discretization points n p between 50 and 250. For each n p , we compute K s according to the method presented in Section 3 (see Listing 1 for the Matlab code) and use the Matlab command eig to obtain the eigenvalues form which we select the first n e = 40 (in ascending order according to their absolute value). Figure 5 (left panel) shows the convergence of the maximum relative error (17) for fractional order s = 0.9, s = 0.5, s = 0.3 and s = 0.1. The convergence of eigenvalues for the standard Laplacian discretized via finite elements is shown as reference (right panel). Even if the a priori convergence rate of our method is not known, we observe that, as for the standard Laplacian, it decreases in O(h 2 ), which is the expected convergence rate for the finite element method with P1 Lagrange elements.

A simple Poisson problem
Another simple test can be done in order to exemplify the effect of different fractional orders. We consider the one-dimensional Poisson problem on the interval Ω = [0, 1] where f (x) = 6x + 2, the fractional operator (−∆) s is understood in the spectral definition (2) and with homogeneous Dirichlet boundary conditions u| ∂Ω = 0. We want to study the convergence of the solution computed applying the method presented in Section 3 as the mesh size decreases. While the solution to problem (18) for the standard Laplacian (s = 1) reads u D (x) = −x(x − 1)(x + 2), the analytical solution for the fractional Laplacian problem is not known. Therefore, we use the numerical solution on a mesh with 500 nodes,ū, as reference solution. We compute solutions for fractional orders ranging between s = 0.1 and s = 0.9, increasing the number of nodes from n p = 10 to n p = 250. The convergence towards the numerical solution as the mesh size decreases is shown in Figure 6 (left panel) for s = 0.75, s = 0.5, and s = 0.25. We observe that the convergence rate decreases from O(h 2 ) for s = 0.75 to O(h) for s = 0.25. Note that a similar effect has been observed in [22].
The deformation of the numerical solutions as s decreases can be seen in Figure 6 (right panel). As s → 0, the solution tends to 6x+2 which does not satisfy Dirichlet boundary conditions. Therefore this creates an abrupt change in the solution at the boundary values which could affect the convergence.

Bifurcation diagrams with fractional Laplacian
In this section we present the results obtained exploiting the new capabilities of the continuation software to the study of the bifurcation structure of three benchmark problems, namely the Allen-Cahn (4) and Swift-Hohenberg (5) equations, and the Schnakenberg system (6)- (7). The main purpose is to highlight the common effects of fractional diffusion on the bifurcation structure and on the solutions shape.

Fractional Allen-Cahn equation
We start with the fractional Allen-Cahn equation with cubic-quintic nonlinearity (4) on a 1D domain of length L with homogeneous Dirichlet boundary conditions. As for the standard Allen-Cahn equation (8), the homogeneous solutionū = 0 is a stationary steady state of equation (4). Furthermore, using the Fourier property of the fractional Laplacian, one can easily find the bifurcations from the homogeneous states to be located at from which we expect the bifurcation points to move towards µ = 1 as s tends to zero. To obtain a more global picture of the bifurcation diagram far from the homogeneous solution, numerical continuation  Figure 7 shows that, indeed, the shifting of the first and second bifurcation points as s decreases agrees with the values computed using formula (19). The maximum relative error for the first bifurcation point is 4 · 10 −3 and for the second bifurcation point it is 2 · 10 −3 . Figure 8 shows the evolution of the bifurcation diagram (with respect to the L 2 -norm of u) of equation (4) as the fractional order of the Laplacian decreases. We observe the gathering of branches which is a consequence of the shifting of bifurcation points. In addition, one can see that a stable region appears on the second and third bifurcating branches at s = 0.2 (Figure 8c). These stable states exist in fact for all fractional orders, however they appear at much smaller µ values when s is decreased. This is better understood by looking at the solution profiles. We have seen in Section 2.1 that the steady state solutions of the Allen-Cahn equation correspond to situations where the space is divided into subregions occupied by one of two phases. We also know that the interface between these phases sharpens as µ increases. Similar observations can be made for the steady states of the fractional Allen-Cahn equation (4). However, the "sharpening" of solutions occurs at much lower µ, as s decreases. Thus, it seems reasonable for the stability of solutions to be shifted to lower µ values too. Figure 9 illustrates the change in solution profiles as s decreases, for µ = 2 fixed, on the first, second and third branch. Note that for each solution profile u(x) in Figure 9, −u(x), is also a solution to (4). These two solutions lie on the same branch in the bifurcation diagram due to symmetry.

Fractional Swift-Hohenberg equation
We now study the steady state fractional Swift-Hohenberg equation with cubic-quintic non-linearity (5), on a 1D domain of length L with homogeneous Dirichlet boundary conditions. The homogeneous statē u = 0 is a stationary solution of the fractional equation (5) and it can be shown that the bifurcations from the homogeneous states occur at Thus, as in the standard problem, the homogeneous stateū is stable for µ < 0 and first becomes unstable at µ c = 0. In addition, note that as s tends to zero, µ j tends to µ c = 0. Substituting (u 1 , u 2 ) = (u, ∆ s u) in the stationary equation (5), we convert the Swift-Hohenberg equation into a system of second order PDEs which fit the pde2path framework based on finite elements. To obtain the following numerical result, we fix ν = 2 and the domain [−5π, 5π] of length L = 10π with boundary conditions: u(−5π) = u(5π) = 0 and ∆ s u(−5π) = 0 = ∆ s u(5π) = 0. We choose the meshsize h = 0.04, corresponding to np = 786 mesh points.
As mentioned above, the bifurcation points on the homogeneous branch accumulate at µ c = 0 as the fractional order s tends to 0. Figure 10 shows, as an example, the shifting of the second and third bifurcation points. The values obtained numerically are marked with dots and the analytical values, computed using formula (20) with L = 10π, are showed as lines on the figure. We find that the maximum relative error is on the order of 10 −2 for both points and observe that the numerical values indeed coincide with the ones expected from the analysis. In addition, we observe that the positions of the two bifurcation points seem to be interchanged close to s = 0.5.
In order to get a general impression of the effects of fractional diffusion on the bifurcation structure of problem (5), Figure 11 shows the bifurcation diagram for four different fractional orders (namely  The homogeneous branch is shown in black, the first branch bifurcating at µ = 0 is shown in blue and the first snaking branch in red. As we already said, secondary bifurcation points shift towards the primary bifurcation. Looking at the computed bifurcation structure, we notice two other important changes as s decreases. Firstly, the width of the snaking branch is significantly increasing. Secondly, for s = 0.5 and s = 0.3 the snaking branch does not reconnect to the blue branch anymore, but to different secondary branches. The increasing width of the snaking branch is better seen in Figure 12, where only the first snaking branch is shown for s = 0.9, s = 0.7, s = 0.5 and s = 0.3. In addition to the increase in width, one can see in Figure 12 that the back and forth oscillations forming the snaking tend to tilt left for smaller s. This phenomenon is particularly visible for s = 0.3. The solutions along the tilted branches are illustrated in Figure 13. We see that the oscillations significantly enlarge as we move up the snaking branch. They would then narrow again after the next turn in an "accordion" effect. Further, counting the number of turns on the snaking branches in Figure 12, we see that whereas for s = 0.9 we have 8 pair of turns, for s = 0.7 and s = 0.5 we have 9 pairs and for s = 0.3 only 6. As shown in Figure 13, the first snaking branch corresponds to front solutions. At each snaking turn, one oscillation is added to the front. Thus, the number of turns is directly related to the final number of oscillations in the solution. This is linked to the next observation which regards the reconnection of the snaking to a branch of spatially periodic solutions. In fact, as s decreases the snaking branch does not reconnect to the first periodic branch anymore but to periodic branches originating from subsequent bifurcations on the homogeneous branch. For instance, for s = 0.5 and s = 0.3 it connects to the branches originating from the the third and ninth bifurcation points, respectively, on the homogeneous (black) branch. Magnified views of the bifurcation diagram for s = 0.5 close to the reconnection and the starting points are reported in Figure 14.  Finally, the last comment has to be made on the shape of the solutions. Looking at Figure 11, one can notice that, as the fractional order s decreases, the spatially periodic states (blue-branch) have slightly lager L 2 -norm values. The solutions indeed tend to "sharpen" for smaller values of s closer resembling a square-profile. This is illustrated in Figure 15, where the solution profiles on the first (blue) branch at µ = 4 are shown for different values of the fractional order s. We see that, whereas the solution at s = 0.9 and s = 0.7 do not differ much, at s = 0.5 the profile significantly "flattens" towards the top of oscillation. Note that this pattern has also been observed in the previous section for the periodic solutions of the fractional Allen-Cahn equation.

Fractional Schnakenberg system
We finally study the fractional Schnakenberg system (6). The problem is considered on a one dimensional domain of length L with homogeneous Neumann boundary conditions. We study both classical (σ = 0) and modified (σ = 0) reaction terms (7) as presented in the Introduction.
As for the Schnakenberg system with standard Laplacian (10), the spatially homogeneous statē u = (µ, 1/µ) is a stationary solution to (6) . The bifurcation points on the homogeneous branch turn out to be located at Further, the critical parameter µ c at which the first instability occurs is not affected by the fractional diffusion, being The critical wavenumber, however, changes as a function of the fractional order s Thus, we expect the wavelength of spatially periodic solutions from the first bifurcation to increase as s decreases. In the following, we fix the diffusion coefficient d = 60. Thus, equation (22) leads to the critical parameter µ c ≈ 3.21. Further, we choose the domain i.e. the domain is "tuned" to the primary bifurcation for each fractional order s. Precisely, according to equation (23) we know that k c decreases with s, thus the size of Ω, L = 2mπ/k c , is increased with s in such a way that the first unstable mode always has m periods in the domain. This choice guarantees that the primary bifurcation always occur at µ c to modes with known wavenumber given by formula (23). Note that for the Allen-Cahn equation and for the Swift-Hohenberg equation, the domain size was always fixed since in these systems k c did not depend on s. Furthermore, because of this choice and to be able to compare results for different fractional orders, it is important to use a measure for the bifurcation diagram that does not depend on the domain size. The normalized L 2 -norm which was used in previous sections satisfies this criteria. In addition, we also use the normalized L 8 -norm [53] ||u|| , since in some situations it has been shown to be more suitable for plotting than the standard normalized L 2 -norm as it produces larger differences between states. Finally, note that since the domain size increases as s decreases, we must accordingly increase the number of discretization points to guarantee the accuracy of the computation (the values used in the following are listed in Table 1). This induces a significant increase in computational cost. For this reason we studied the problem with fractional order not smaller than s = 0.5.

Numerical results for the non-modified system
In this section we consider the fractional Schnakenberg system (6) with classical reaction terms, i.e. σ = 0. We choose a relatively small domain, namely m = 2 in (24), to limit the increase in computational cost as the fractional order s gets smaller.  Figure 16: Spatially periodic solutions of the Schnakenberg system (6) with classical reaction terms (σ = 0) close to the first bifurcation from the homogeneous states. Left: decrease of the critical wavenumber k c . Right: increase in the wavelength of the solutions for µ ≈ µ c for different fractional orders s. Remember that to each fractional order corresponds a different domain size according to (24); in this figure we show only part of the solutions for a comparison.
Before showing the bifurcation diagram of system (6), let us first illustrate the increase in the wavelength of solutions as s decreases. The critical wavenumber, computed via equation (23), for different values of s is plotted in Figure 16a. Figure 16b shows the corresponding change in wavelength for solutions close to the first bifurcation at µ ≈ 3.2. Further, one can see that the amplitude of the solution decreases with s. In fact, we notice that the L 2 -norm of these solutions is conserved for all s from 0.5 and above (as mentioned, solutions below s = 0.5 have not been computed).
In Figure 17 we report the bifurcation diagrams of system (6) for three different fractional orders (namely s = 0.9, s = 0.7 and s = 0.5). We observe several changes as s decreases. First, the bifurcation points on the homogeneous branch (in black) seem to move towards the first bifurcation at µ c . In addition, the second and third bifurcation points seem to exchange their position between s = 0.8 and s = 0.7. Correspondingly, the second (red) and third (magenta) bifurcating branches move closer and towards the first branch (blue) when s becomes smaller. The numerical values of the bifurcation parameter µ for the location of the second and third bifurcation points are shown in Figure 18. These values can easily be verified using formula (21) with k j = jπ/L for the domain size L = 4π/k c and the corresponding curves are plotted as reference in Figure 18. We observe that the numerical values agree with the analytical ones (with a maximum error of the order of 10 −3 ). In addition, the bifurcation points come closer as the fractional order s decreases and intersect close to s = 0.7. Note that the bifurcation points eventually converge towards µ c as s → 0. Second, we observe in Figure 17 that the L 2 -norm of the branches of solutions is increasing. This is better understood by looking at the solution profiles. Figure 19 shows the solution on the first bifurcating branch at µ = 3 for s = 0.9 and s = 0.5. We see once more the wavelength increasing. In addition, for lower s the solutions "valleys" tend to flatten. This is now a known effect, which we had already observed in the solutions of Allen-Cahn and Swift-Hohenberg equations. In contrast to previous systems, however, the amplitude of the solution increases and the "peak" sharpens.
The last observation arising from Figure 17 regards the stability of solutions. As for the standard problem, with fractional diffusion the first branch bifurcating from the homogeneous states (blue) is entirely stable and the second one (red) becomes stable at some µ < µ c . But, in contrast to the standard case, the third branch is unstable for s = 0.9 and s = 0.7 and a stable zone appears at s = 0.5. This can be a consequence of the position exchange of the second and third bifurcation points on the homogeneous branch (see Figure 18).

Numerical results for the modified system
We now turn to the fractional modified Schnakenberg system, i.e. equation (6) with σ = 0 in the reaction term. Recall from Section 2.3 that in this case we expect to observe a snaking branch in the bifurcation diagram. Since in small domains snaking branches have very few turns, we choose here a domain larger than in the previous section in order to have more turns and facilitate our observations, namely m = 5 in (24). Since we consider larger domains as the fractional order decreases, the computational costs increase, as explained in the previous section. Therefore we only consider fractional orders of the Laplacian greater than s = 0.7. The numerical and analytical values for the location of the bifurcation points for fractional order s between 0.7 and 0.9 are shown in Figure 20. The analytical values are computed using formula (21) with m = 5 and and agree with the numerical ones up to 3 decimal digits. However this accuracy could not be maintained for fractional orders bellow s = 0.7. Note that, from the analytical curves, we see that the bifurcation points cross each other close to s = 0.7.
In order to get a first impression of the effects of the fractional Laplacian on the bifurcation structure of system (6) with modified reaction terms (σ = 0), Figure 21 shows the bifurcation diagram for four different fractional orders, namely s = 0.95, s = 0.9, s = 0.8 and s = 0.73. In the bifurcation diagrams, the blue branch, called P 5, corresponds to spatially periodic solutions with five oscillations in the domain. The light blue branch, called P 4, corresponds to periodic solutions with four oscillations in the domain. Analogously, the red (P 4.5) and the green (P 5.5) branches correspond to periodic solutions with 4.5 and 5.5 oscillations in the domain, respectively. Further, the snaking branch generated at the first bifurcation point on P 5 is shown in purple. We observe three important changes in the bifurcation structure as s decreases. First, as expected from the previous section, the bifurcation points are accumulating towards the critical point µ c as s gets smaller. Second, looking closely at the bifurcation diagrams, we see that the stable solutions on the blue branch seem to move towards the fold when s decreases. Third, we observe that the snaking is widening as s becomes smaller. In addition, at s = 0.73, it neither reconnect to P 4 nor "climb" to other branches. In detail, Figure 22 illustrate how the stable solutions shift towards the fold, that is towards higher µ values, as s decreases. Recall that this effect was also observed in the fractional Allen-Cahn system (Section 4.1). Moreover, the fractional order appears to have a significant influence on the behavior of the snaking branch. As mentioned above, we have seen in Figure 21, that the number of turns and height of the snaking branch decreases progressively, and it does not "climb" anymore for s = 0.73   ( Figure 21d). In addition, we can observe an increase in the width of the snaking as s gets smaller. Figure 23 shows the evolution of the branch in more details for fractional orders s = 0.78, s = 0.73 and s = 0.7. Furthermore, Figure 24 shows solution profiles along the snaking branch at s = 0.73. As we had seen in Section 2.3, they initially correspond to front solutions. Going up the snake, we expect the front to move "forward", that is gaining more and more oscillations until the domain is filled, to match the number of oscillations of the reconnecting branch. However, here, we see that the oscillations does not gain the fourth complete bump. In fact, at the upper fold point the front stops moving forward. In addition, unstable localized solutions appear after the fourth fold, which progressively present oscillations in the central part of the domain (Figure 24f). This change of behavior can be due to an interaction between the first snaking branch of front solutions with the snaking branch of homoclinic solutions originating (in the standard Schnakenberg system) from the second bifurcation point on P 5 [53].

Conclusion and Outlook
In this paper we have successfully extended for the first time numerical continuation methods to nonlinear space-fractional PDEs. In particular, in the context of pde2path, we have interfaced FEM with the fractional Laplacian via a discretization of the Balakrishnan formula. This enabled us to investigate some of the effects of super-diffusion on the steady state bifurcation structure and solution profiles of the Allen-Cahn equation (4), the Swift-Hohenberg equation (5) and the Schnakenberg system (6). We highlight here the common effects we have observed, which contribute to a better understanding of the effects of the fractional Laplacian on the steady solutions of generic reactiondiffusion systems on bounded domains.
-Firstly, in all systems we have observed that the bifurcation points on the branch of homogeneous solutions accumulate at a precise value of the bifurcation parameter as the fractional order of the Laplacian decreases.
-Secondly, we have seen that, for a fixed value of the bifurcation parameter, the spatially periodic solutions to the three equations tend to have flatter peaks or valleys as the fractional order of the Laplacian decreases. This effect can be memorized as "sharpening the teeth" of the spatial solution profiles.
-Finally, in both systems featuring a snaking branch of localized solutions, that is the Swift-Hohenberg equation and the Schnakenberg system, we have observed that decreasing the fractional order of the Laplacian leads to a significant widening of the snaking. This can be memorized as "bloated snaking". Furthermore, the re-connection properties of the snaking branch change as the fractional order is decreased.
Hence, our results also pave the way for further analytical investigations of nonlinear reactiondiffusion systems involving fractional operators. Our investigation focused on one-dimensional fractional problems, where matrices associated with the discretized problem can be kept sufficiently small. Yet, these problems can be treated numerically on today's standard desktop computers. To discretize the fractional operator in higher dimensions poses additional numerical challenges. To obtain accurate results the Balakrishnan formula entails that one has to solve of a large number of linear systems of equations of significant size. In order to make the computation more efficient, a method based on Krylov subspaces has been implemented recently [22], which could be helpful in future work. Of course, it could be useful to use other discretizations of the spectral fractional Laplacian on bounded domains based on the integral formulation of the operator via the heat-semigroup formalism as proposed in [18]. Another way could be the very weak FEM [5]. Furthermore, from a numerical perspective, the choice of parameters k, N + , N − , appearing in the FEM discretization of the fractional Laplacian could be investigated. Similarly, the detailed study of the numerical accuracy is an issue that could be studied in the future; since we have benchmarked the method already against analytical results, which gave excellent agreement with the numerics, we expect that the numerical analysis is also going to yield favorable results.
Beyond perspectives of computational optimizations, the preliminary results obtained within this work bring up fascinating issues, which will be addressed in future works. For instance, the deformations and changes of the snaking branch in the Schnakenberg system for fractional order smaller than 0.7 is an open issue. A natural extension at this point is also the numerical investigation of these models in higher-dimensional domains. Regarding other models, it could be interesting to investigate if the fractional Laplacian can lead to the appearance of unexpected patterns (which are not shown with standard diffusion). Finally, implementing other FEM discretizations for different definitions of the fractional operator [10,24], we will be able to contrast and compare them by looking at the bifurcation diagrams and stationary solutions.
In conclusion, the paper constitutes the meeting point between theoretical results, bifurcation theory, numerical analysis and scientific computing. Extending advanced continuation methods, such as implemented in pde2path, allowed us to treat a much larger class of systems. In particular, we gained a new tool to investigate bifurcations and pattern formation in fractional PDEs. Thus, we established a novel interaction between continuation techniques and fractional PDEs. This point of view will shed new light on the field as it allows for global exploration of nonlinear systems involving the fractional Laplacian. We can now systematically apply our approach to investigate a large class of problems in many fields (biology, epidemiology, population dynamics). Furthermore, the paper can be viewed as one key contribution towards a dynamical system approach to fractional reaction-diffusion systems.

A The MatLab code
We provide here the Matlab code which implements the FEM discretization of the fractional Laplacian presented in Algorithm 1 in Section 3. We do not explain in detail the basic setup of pde2path; see [53,54,45] for a complete overview of the continuation software for beginner users.
In order to numerically study fractional reaction-diffusion equation of the form (1) in pde2path, we must adapt the standard setup providing a new function FractionalLaplacian which will be called by the function oosetfemops. We recall that the routine oosetfemops generates the matrices (14) appearing in the algebraic system (13). In the fractional case, this file differs from the standard setting only for the line by which the FEM stiffness matrix K, stored in the Matlab variable p.mat.K, is replaced by M K s , the matrix approximation of the fractional Laplacian multiplied with the mass matrix. Else the usage of pde2path stays the same for the basic continuation calls.
The new function FractionalLaplacian implementing Algorithm 1 is presented in Listing 1. The required input arguments are the standard stiffness and mass matrices K and M, the number of mesh points np, the (uniform) mesh size h and the fractional order s. The output of the function is the matrix approximation of the fractional Laplacian Ks.
Few remarks need to be made concerning the Matlab implementation. The first comment regards the choice of parameters k, n + and n − in the quadrature approximation (15). Given a spacial discretization with mesh size h, the value of k must be chosen such that the quadrature error balances with the finite element error. Following [22,Corollary 2] where a similar discretization method was used, we choose k = 1 |log(h)| , corresponding to Line 9 in Listing 1. Then, again following [22], the values n + and n − are chosen proportional to 1/k 2 : n + := π 2 4(1 − s)k 2 and n − := π 2 4sk 2 , and they are assigned in Lines 10, 11 in Listing 1. The second comment regards the efficiency and correctness of the computation. For each column of K s we need to solve n + + n − + 1 systems of equations (see equation (16)) of the type (e kl M + K)y = Kê i , where −n − ≤ l ≤ n + . Since all columns as well as all systems are independent, ideally, we would like to parallelize both loops in Lines 16 and 23. However, nested parallelism is not allowed in Matlab. Therefore, we only parallelize the outer loop (using parfor in Line 16).
Finally, when e kl becomes sufficiently small, we have that (e kl M + K) ≈ K. Since K is singular, there exist infinitely many solutions to systems of the form Kv = z. In the Matlab implementation, the typical backslash operator would produce one solution but we would not know which solution it produces. In order to have a more robust routine, we use lsqminnorm which computes the minimum norm least-squares solution to the system using the complete orthogonal decomposition of (e kl M + K) (see Line 26 in Listing 1).