The Complex-Order Fractional Laplacian: Stable Operators, Spectral Methods, and Applications to Heterogeneous Media

Fractional differential equations have become a fundamental modelling approach for understanding and simulating the many aspects of nonlocality and spatial heterogeneity of complex materials and systems. Yet, while real-order fractional operators are nowadays widely adopted, little progress has been made in extending such operators to complex-order counterparts. In this work, we introduce new definitions for the complex-order fractional Laplacian, fully consistent with the theory of averaging of smooth functions over fractal sets, and present tailored spectral methods for their numerical treatment. The proposed complex-order operators exhibit a dual particle-wave behaviour, with solutions manifesting wave-like features depending on the magnitude of the imaginary part of the fractional order. Reaction-diffusion systems driven by the complex-order fractional Laplacian exhibit unique spatio-temporal dynamics, such as equilibrium of diffusion in random materials by interference of scattered waves, conduction block and highly fractionated propagation, or the generation of completely novel Turing patterns. Taken together, our results support that the proposed complex-order operators hold unparalleled capabilities to advance the description of multiscale transport phenomena in physical and biological processes highly influenced by the heterogeneity of


Introduction
Diffusion and diffraction are ubiquitous mesoscopic transport phenomena in any complex media, including heterogeneous, random, porous, and fractal materials [1]. A particle-like diffusive behaviour arises when the free path of the material is much longer than the propagating wavelength, while for media encompassing obstacles at the wavelength scale, diffusion theory breaks down and wave-like diffractive behaviour emerges [2]. However, deriving mathematical descriptions capable of linking both types of transport processes remains an open challenge, as the diffusion equation does not support diffraction effects. Attaining such descriptions could therefore have important implications in multiple fields of physics where diffusion and diffraction often coexist, including acoustics, atmospheric sciences, seismology, astrophysics, nanophotonics, telecommunications, and medical physics, to cite a few (see [1][2][3] and references therein). Especially, they could provide a unifying framework of diffusion-diffraction dualities, such as the formation of diffracting (instead of diffusion-driven) Turing patterns [4,5], or the generation of diffracting features in reaction-diffusion systems of phase transition [6][7][8], which still remain poorly understood to date.
In the last decades, fractional diffusion has gained increased recognition for the mesoscale modelling of processes deviating from standard diffusion in heterogeneous me-dia [9,10]. In particular, the real-order fractional Laplacian has been extensively exploited in the context of superdiffusion [11]. Importantly, recent theoretical results have established a reciprocity between complex-order fractional operators and the averaging of smooth functions over fractal sets [12][13][14]. Such results could represent a powerful methodology to account for inhomogeneities at a fractal (micro-)scale, while still preserving a mesoscale description of heterogeneous media. To the best of our knowledge, only the recent work by Ugarte et al. has aimed at connecting such results with complex-order definitions of the fractional Laplacian [15][16][17]. However, as will be discussed, their proposed complexorder fractional operator is heavily ill-conditioned for its application to diffusion processes. Therefore, suitable definitions of the complex-order fractional Laplacian are still needed.
In addition to the above, and even in the real-order case, the numerical treatment of the fractional Laplacian is challenging due to its non-local nature [11]. By embracing such a non-local nature, we previously proposed Fourier spectral methods as highly accurate and efficient numerical stencils for problems described by the real-order fractional Laplacian, scaling to any spatial dimension in the form of fully diagonal discretisations [18]. These methods have been broadly adopted in applications as diverse as imageprocessing [19], non-linear optics [20], or biology [21]. In the complex-order case, the numerical complexity of the fractional Laplacian is further exacerbated due to the presence of complex arithmetic. Appropriate numerical methods are therefore also required for these new classes of complex-order fractional operators.
In this work, we introduce suitable definitions for the complex-order fractional Laplacian, fully consistent with the aforementioned theory of averaging over fractal sets [12][13][14], and propose tailored spectral methods for their numerical treatment. Numerical simulations of diffusion processes driven by the complex-order fractional Laplacian demonstrate that the proposed operator attains the intended connection between diffusion and diffraction phenomena, with solutions exhibiting a dual particle-wave behaviour where wave-like features arise depending on the magnitude of the complex part of the fractional order. We further investigate the resulting oscillatory patterns in reaction-diffusion systems characteristic of wave propagation, relevant to the formation of phase transitions, excitable media, and morphogenesis. These complex-order fractional reaction-diffusion systems exhibit thought-provoking properties, such as light-like Anderson localisation (equilibrium of diffusion in random materials due to the interference of the scattered waves) [22,23], conduction block and highly-fractionated propagation, capable of wavefront fragmentation (pertaining to the presence of microfibrosis clefts in cardiac tissue) [24,25], and completely novel temporal dynamics and subsequent Turing patterns than their respective standard diffusion or real-order fractional diffusion counterparts [18]. Taken together, our results indicate that the proposed approach may have important implications in unifying different fields of physics, in yielding new applications of dynamical systems, and in unravelling the many facets of transport phenomena in research fields highly influenced by the heterogeneity of complex media.

Consistent Definition of the Complex-Order Fractional Laplacian with Diffusion on Fractals
Let us denote the standard and the real-order fractional Laplacians by the symbols (−∆) and (−∆) α /2 , defined on a bounded domain Ω ∈ R n . If the standard Laplacian has a complete set of orthonormal eigenfunctions {ϕ j } satisfying standard boundary conditions, with corresponding eigenvalues {λ j }, i.e., (−∆)ϕ j = λ j ϕ j , then both operators are related by the so-called spectral definition of the fractional Laplacian with power-law real fractional order 0 < α/2 ≤ 1, so that for α = 2 this recovers the standard Laplacian case. Here, we will concentrate on homogeneous Dirichlet and homogeneous Neumann boundary conditions in rectangular domains of R n , for which λ j ≥ 0, so we can omit the absolute value in equation (1) for clarity of notation. Let us consider at this point the real-order fractional diffusion equation ∂ t u(x, t) = −D(−∆) α /2 u(x, t), with diffusion coefficient D and initial data u 0 (x) = u(x, 0). By using a series expansion of the form u(x, t) = ∑ ∞ j=0ûj (t)ϕ j (x), where the functions ϕ j (x) satisfy the appropriate boundary conditions, then the solution to this equation is given by (see also § 2.2) which exhibits a bounded behaviour for all times given the positiveness of its spectrum. However, the direct extrapolation of equation (1) to complex power-law fractional orders of the form ν = (α ± iβ)/2 results in an inappropriate definition of the complexorder fractional Laplacian as an effective modelling technique in physical scenarios. This is clearly seen from the equality First, any real-valued initial condition (such as a density or concentration) would unequivocally evolve into a nonphysical complex-valued magnitude. Secondly and perhaps even more importantly, for | β /2 log λ j | > π/2, the periodic-logarithmic modulation induced by the imaginary part β would shift signs of the real part of the eigenvalues, resulting in an exponential blowup of the solutions when substituted into equation (2). Also note that, for a fixed domain size, the eigenvalues {λ j } numerically increase with increasing mesh resolution (see § 2.2). This implies that, if equation (3) were directly used for describing the complex-order fractional Laplacian, no numerical scheme for equation (2) or derived reaction-diffusion counterparts would be unconditionally stable when performing mesh refinement, even for small values of β.
In order to derive a suitable spectral definition for the complex-order fractional Laplacian, here we recall recent results on fractal sets [13]. In particular, for complex power-law orders ν = (α ± iβ)/2, the operator defining the spatial averaging of smooth functions over fractal sets is given by where ω n = 2πn/ log ξ is the set of frequencies providing periodicity with scale log ξ. Here, ξ < 1 defines the scaling parameter of the fractal object at its microscale, and z = ik · x represents the Fourier parameter in space. Equivalently, with C n andC n being complex conjugate constants. As shown in [12], the expression above can be approximated by where ω is the averaged frequency referring to the leading term. Therefore, in the limit where the spatial fractal scale ξ → 0, the averaging process yields At this point, and remembering that z represents the Fourier eigenvalues in space, we note the direct equivalence of equation (7) with the following definition for the complexorder fractional Laplacian with corresponding spectrum where the choice of constants A 0 = C 0 and A 1 = C 1 /2 is only for improved compactness of the expressions derived below. More rigorously, the equivalence above can also be established by exploiting the correspondence of each of the terms of equation (7) with the Green's functions defining the fractional Laplacian [15]. Very importantly, the definition of the complex-order fractional Laplacian given by equation (8) is more general than that suggested by Ugarte et al. [15][16][17]. In their work, they apply a combination of two complex-conjugate fractional Riesz operators in the form given by equation (3), which in our case would be equivalent to imposing A 0 = 0 and A 1 =Ā 1 = 1. Whereas such a choice circumvents the issue of real solutions becoming complex-valued, the resulting eigenvalues λ α /2 j cos β /2 log λ j can still lead to negative real parts and ill-conditioned behaviour when increasing mesh resolution, as previously discussed. Alternatively, a given mesh resolution severely constrains the range of possible values of β to be considered. In addition, the use of fractional Riesz operators instead of the fractional Laplacian can yield curvature artefacts on reaction-diffusion models of wavefront conduction, as previously discussed [26] Further insights on the definition of the complex-order fractional Laplacian given by equation (8) can be obtained by inserting A 1 = Ae −iθ . Together with equation (3), and using A 0 = 1 for clarity of notation (as this constant can be assimilated into the corresponding diffusion coefficient), it yields We note that, for β = 0, the complex-order operator must recover its real-order counterpart. This is only achieved by θ = ±π/2, which implies the constants A 1 andĀ 1 must be purely imaginary-conjugate numbers, instead of real ones. Using as convention A 1 = Ae −iπ/2 andĀ 1 = Ae iπ/2 , A ≥ 0, the complex-order fractional Laplacian can be finally written as where 0 < α/2 ≤ 1 and β ≥ 0, with spectrum where 0 ≤ A < 1 is the only stability condition in order to preserve its positiveness. The upper modulus bound A = 1 constitutes a degenerate case as an arbitrary number of high-frequency wave numbers would remain unattenuated, which is inconsistent with the expected behaviour of a diffusion operator. An important observation regarding the complex-order fractional Laplacian is that, despite its more intricate definition in operator form as provided by equation (11), its spectral decomposition becomes a straightforward extension of the real-order counterpart, with eigenvalues modulated by the periodic-logarithmic term given by equation (12). This indicates the suitability of Fourier spectral methods for the numerical treatment of dynamical systems driven by the complex-order fractional Laplacian, as presented next.

Fourier Spectral Methods for the Complex-Order Fractional Laplacian
As an extension of the methods presented in [18], and in order to illustrate the basis of Fourier spectral methods applied to the complex-order fractional Laplacian, let us consider for simplicity the one-dimensional complex-order fractional diffusion equation subject to u(x, 0) = u 0 (x), and either homogeneous Dirichlet or Neumann boundary conditions. This setting can be easily extended to any rectangular domain in R n by means of tensor products, as illustrated in the accompanying codes (see Supplementary Materials). The solution to this problem can then be approximated in the form of a truncated Fourier series where N is the number of discretisation points. The eigenfunctions ϕ j (x) and corresponding eigenvalues λ j of the standard Laplacian are determined by the specified boundary for homogeneous Dirichlet, and for the homogeneous Neumann type, where L = b − a. Inserting this ansatz in equation (13), one gets whereλ j = λ α /2 j 1 + A sin β /2 log λ j are the eigenvalues of the complex-order fractional Laplacian, as given by equation (12). Taking the scalar product of both sides of the last expression by a test function ϕ i (x), and recalling the orthonormality of the basis functions, i.e., ϕ i (x), ϕ j (x) = δ i,j , where f (x), g(x) denotes the scalar product and δ i,j is the Kronecker delta function, then each coefficientû j (t) satisfies with exact solutionû whereû j (0) = u 0 (x), ϕ j (x) . For each type of boundary data, the coefficients {û j (0)} (as well as the inverse reconstruction of u(x, t) in physical space from the set {û j (t)}) can be efficiently computed by fast and robust discrete sine/cosine transforms (see [18]). For complex-order fractional reaction-diffusion equations of the form subject to equivalent initial and boundary data, the application of the same arguments as above yields where {f j } are the Fourier coefficients of f (u) at the start of the time interval [t, t + ∆t]. This admits the following exact exponential integration stencil to advance the numerical solution of equation (18) The time-integration stencil above assumes a small variation of f (u) in the time interval [t, t + ∆t], for sufficiently small ∆t → 0. Other time-integration schemes, such as fixed point iteration [18,27], (semi-)implicit discretisations [28], or splitting methods [29,30], also constitute feasible alternatives for the time evolution of equation (18). As an example, the codes provided for the Gray-Scott model (see Supplementary Materials) illustrate how to accommodate easily a second-order Strang splitting method [31,32] into the proposed Fourier spectral methods for the complex-order fractional Laplacian. In addition to the efficiency derived from the use of discrete Fourier transforms, the numerical schemes given by equations (17) and (18) hold a number of important properties. Firstly, the orthogonality of the basis functions implies that each Fourier coefficient evolves independently to the others. This translates into a fully diagonal system, thus requiring no preconditioning or the solution of expensive non-sparse linear systems per time step, as it is often the case of finite-difference or finite-element counterparts. This property also holds when extending the proposed schemes to higher dimensions in R n , as illustrated in the codes in the Supplementary Materials. They also avoid associated numerical challenges for the treatment of singular Laplacians (containing the eigenvalue λ j = 0), as in the case of homogeneous Neumann boundary conditions (see [27]).

Complex-Order Fractional Diffusion Unifies Diffusive and Wave-Like Behaviours
As a starting point in our presentation of complex-order fractional diffusion dynamics, let us consider the diffusion equation driven by the complex-order fractional Laplacian defined in [0, 1] 2 under homogeneous Dirichlet boundary conditions, diffusion coefficient is the Dirac delta function. In terms of the fractional complex order, we will concentrate on exploring the upper range of admissible values for its real part α (typically kept in 1.5 ≤ α ≤ 2) to facilitate comparisons with previous works, while varying the imaginary part β ≥ 0. Figure 1 illustrates how a dual particle-wave behaviour emerges naturally from the proposed definition of the complex-order fractional Laplacian, given by equation (11). Recall its eigenvalues are given byλ j = λ α /2 j 1 + A sin β /2 log λ j , which introduces 0 ≤ A ≤ 1 as an additional controlling parameter of the modulation exerted by the imaginary part β of the fractional complex order. Results in Figure 1a (α = 2, A = 0.99, varying β) depict the periodic-logarithmic modulation of the complex-order fractional diffusion eigenvalues (λ (α±iβ)/2 ) against those associated with standard diffusion (λ α=2 ). Note that the eigenvalues always remain positive regardless of β, therefore ensuring the stability of the complex-order operator. Increasing values of β translate into a progressive periodic attenuation of the diffusion eigenvalues linked to smaller wave numbers, and therefore to the appearance of wider spatial wave-like features. A superposition of diffusion with a larger number of wave-like harmonics is present in the solutions for increasing values of β. In fact, it is important to remark that, regardless of the strong wave-like features exhibited by the solutions, these dynamics arise from the simulation of a (complex-order) diffusive process, and therefore u → 0 as t → ∞, compared to a pure (non-dissipating) wave equation.
Regarding the role of decreasing the real part α of the complex order in determining the observed dynamics, results in Figure 1b (α = 1.5, A = 0.99, varying β) illustrate a reduction of the fractional diffusion eigenvalues (note the decreased slope of the grey-dashed lines in the eigenvalues plots), especially of those associated with high wave numbers. This translates into the well-known super-diffusive effects of the real-order fractional Laplacian (wider diffusion tails, but with slower decay of peak values; see [18]), which prompt an accentuation of the manifested wave-like patterns. As mentioned before, the observed dual particle-wave behaviour is, in addition, controlled by the amplitude parameter A (subject to 0 ≤ A ≤ 1), which modulates the contribution of the imaginary part β to the complex-order fractional diffusion eigenvalues. Figure 2a shows the outcome of varying the magnitude of this parameter in the solutions to equation (21), with a strengthening in the exhibited wave-like behaviour as A approaches its upper limit. Unless otherwise specified, for the remainder of this manuscript we will concentrate for this reason on the value A = 0.99 in order to investigate nearly-maximal effects of the complex-order fractional Laplacian when driving the dynamics of different reaction-diffusion systems, while avoiding the degenerate case A = 1.
We conclude this section by analysing the numerical convergence of the proposed Fourier spectral methods for the complex-order fractional Laplacian. Convergence results in the solution to equation (21), subject to identical initial data, and varying α and β for A = 0.99, are presented in Figure 2b. Reference solutions were calculated using 2 12 × 2 12 discretisation points. In all situations, the Fourier approach is able to achieve spectral convergence up to machine precision, irrespective of α and β. Results for α = 2 (Figure 2b  lution requirements are slightly increased for α = 1.5 (Figure 2b, right panel), reflecting the slower rate of diffusion for α < 2. Still, mesh resolutions as coarse as 64 × 64 discretisation points are sufficient to achieve global L 2 -norm errors smaller than 10 −10 , even in situations strongly influenced by a marked wave-like behaviour (β = 2.5). Of note is that spectral convergence is achieved regardless of non-smooth initial data (Dirac delta function), and under the most oscillatory regime considered (A = 0.99). Even better convergence is therefore expected under smooth initial data and smaller values of the amplitude parameter A.

Strong Localisation Phenomena in Phase Motion by Complex-Order Fractional Diffusion
We now shift our interest to the study of different reaction-diffusion systems of general interest, when driven by the complex-order fractional Laplacian. Due to their wide use in this type of systems, we will concentrate here on the use of homogeneous Neumann boundary conditions, ∂ n u = 0.
We start by considering the Allen-Cahn equation with a quartic double well potential, a fundamental nonlinear reaction-diffusion model to describe domain coarsening kinetics in alloys, and with wider applications in the study of formation and motion of phase boundaries. The complex-order fractional version of this model takes the form where D is a small positive constant. The steady states u = ±1 are attracting, while u = 0 is unstable. Under standard diffusion (α = 2, β = 0), solutions exhibit flat areas around the attracting states, separated by interfaces of increasing sharpness as D is reduced, while solutions around u = 0 vanish or coalesce over long time scales. Figure 3 investigates the role of a complex fractional order in the dynamics of the Allen-Cahn equation. In these simulations, we used as initial data a Heaviside transition between the two stable states, separated by the polar curve r = 0.65 + 0.35 cos(4θ). This initial condition was chosen in order to incorporate the effect of non-constant curvature in the initial data in the resulting phase motion dynamics. The domain is taken to be [−1.5, 1.5] 2 , with Neumann boundary conditions, and a diffusion coefficient D = 0.01. In the real-order fractional case (0 < α/2 ≤ 1, β = 0, Figure 3a), the interfacial properties of the Allen-Cahn equation have been previously analysed [18,27]: in brief, decreasing α yields steeper but thicker interfaces, reflecting the long-tailed character of the fractional operator, together with a prolonged lifetime of the unstable interface. In the complex-order fractional case (0 < α/2 ≤ 1, β = 0, Figure 3a), these dynamics are complemented with the appearance of periodic oscillations originating from the phase boundary. For larger magnitudes of the real part α, these oscillations can even influence the interior phase, A much more interesting scenario arises when considering larger imaginary parts (Figure 3b). For specific values of β, solutions to the complex-order fractional Allen-Cahn equation can exhibit a steady-state equilibrium of diffusion, by interference of the waves localised at the phase boundary. This closely resembles the so-called Anderson localisation in non-linear optics, where light propagation can be effectively halted as a result of the interference of the waves scattered by a random material [1,22,23]. Anderson localisation is also known as strong localisation in condensed matter physics, to refer to the absence of diffusion due to scattering processes in disordered materials [33,34]. Such a wave-diffusion equilibrium behaviour appears intrinsic to the proposed definition of the complex-order fractional Laplacian. Importantly, it can generate a variety of steady-state structures under same initial data and governing equations (Figure 3b), just depending on the imaginary part β characterising the fractal structure of the medium.

Conduction Block and Fragmentation in Complex-Order Fractional Impulse Propagation
Another ubiquitous application of reaction-diffusion systems is in the investigation of impulse propagation in excitable media. The FitzHugh-Nagumo model undoubtedly represents the most studied non-linear system for describing conduction in nerve axons and cardiac tissue, where the propagation of the transmembrane potential u is modelled via a diffusion equation with cubic reaction term, coupled to a slow linear recovery variable v. The complex-order fractional version of this model takes the form where we consider the following choice of model parameters, a = 0.1, = 0.01, b = 0.5, γ = 1, δ = 0, and a diffusion coefficient D = 10 −4 . In the real-order fractional case (0 < α/2 ≤ 1, β = 0), decreasing α yields a variety of conduction effects, including reduced propagation speed, decreased wavelength, and a wider take-off of activation wavefronts, reflecting the long-tailed nature of the fractional Laplacian [18,21,28].
Moving to the complex-order case, we start our characterisation of the complex-order fractional FitzHugh-Nagumo model by analysing its patterns of radial conduction ( Figure  4). For this, we consider a localised stimulus u(x) = exp(−1000||x|| 2 2 ) as initial condition, with the tissue fully excitable, v(x) = 0, in [−0.5, 0.5] 2 with Neumann boundary conditions. The periodic-logarithmic modulation of diffusion eigenvalues by the imaginary part β results in a periodic regulation of conduction speed and wavelength for increasing values of β (note different time horizons and width of travelling solutions in Figure 4a). Two of the analysed cases deserve special consideration. Firstly, specific ranges of β can yield the complete block of impulse propagation (Figure 4a, β = 0.75), due to a marked reduction of the main wave numbers supporting conduction. Secondly, solutions for larger values of β can portrait considerable oscillations (Figure 4a, β = 2), with conduction proceeding in a nearly discontinuous manner. The magnitude and frequency of such oscillations can be controlled by the other parameters of the complex-order fractional Laplacian, as illustrated in Figure 4b for different values of the amplitude parameter A. Solutions to the complex-order FitzHugh-Nagumo model can also manifest other types of atypical conduction phenomena, such as the propagation of small amplitude ridges preceding the main front (Figure 4a, β = 2.5). However, we will concentrate on the above-mentioned cases characteristic of discontinuous conduction, due to their much wider interest for excitable media.
An important property of these types of systems is that they allow for self-sustained solutions in the form of re-entrant spiral waves. Re-entrant dynamics in the complexorder FitzHugh-Nagumo model under highly oscillatory wavefronts (β = 2) are shown in Figure 5. For these simulations, the trivial state (u, v) = (0, 0) was perturbed by setting the lower-left side of the domain to u = 1 and its upper part to v = 0.1. This allows the initial condition to curve and rotate clockwise, generating the spiral pattern. The rest of model parameters were left unaltered. For large real parts α of the complex order (Figure 5a, α = 2), oscillations are principally located around the spiral tip. Decreasing α leads to smaller wavelengths and more localised oscillations along the entire activation wavefront (Figure 5a, α = 1.75). Greater reductions in α further accentuate the oscillatory nature of the solutions, which combined with the role of α in reducing conduction speed, can eventually lead to wavefront fragmentation (Figure 5a, α = 1.5). The fragmented wavefronts then propagate into multiple self-sustained, chaotic wavelets across the domain, as a trait of impulse propagation in biological tissues with high densities of non-excitable inhomogeneities, such as microfibrosis clefts [24,25]. Importantly, the discussed mechanisms of slow conduction, oscillatory wavefronts, wavefront fragmentation, and onset of self-sustained wavelets are preserved in complex-order fractional representations of three-dimensional excitable media (Figure 5b), independently of the stronger loading effects of diffusion in three-dimensional domains.

Novel Turing Patterns in Complex-Order Fractional Morphogenesis
We conclude our investigation of reaction-diffusion systems when driven by the complex-order fractional Laplacian by considering the Gray-Scott model, as an exemplar in pattern formation and morphogenesis. Its complex-order version is given by where D u , D v , F and κ are positive constants. Here we select D u = 2 × 10 −5 , D u /D v = 2, feed rate F = 0.03, and varying depletion rate κ, as values known to yield interesting dynamics [18,35]. The domain of interest is taken to be [0, 1] 2 , under homogeneous Neumann boundary conditions. As initial condition, the trivial state (u, v) = (1, 0) is perturbed within a square of side length 0.08 to (u, v) = (1/2, 1/4). This initial perturbation then propagates outwards from the central square until the entire domain is affected. Despite its apparent simplicity, the dynamics of the Gray-Scott model are very rich. In brief, under real-order fractional diffusion (0 < α/2 ≤ 1, β = 0), it has been shown that decreasing α significantly alters the transient dynamics of pattern formation, as well as to generate finer-coarse steady-state structures [18,36,37]. Under complex-order fractional diffusion, we concentrate our numerical analysis of the Gray-Scott model in two different values of the depletion rate κ. For κ = 0.055, the model with standard diffusion is known to organise in a steady state field of negative solitons (Figure 6a, α = 2, β = 0). The wave-like behaviour associated with increasing the imaginary part β of the fractional complex order can shift these solutions to the appearance of long-lasting geometrical structures around the trivial state (u, v) = (1, 0) (Figure 6a, α = 2, β = 1). Larger values of β increase the complexity of both the transient dynamics of pattern formation and the resulting structures (Figure 6a, α = 2, β = 2.5). As previously reported [18], a process of nucleation of structures arises at the boundaries for smaller values of α in real-order fractional diffusion (Figure 6b, α = 1.75, β = 0), which then grows inwards until the entire domain reaches its final steady state. However, the presence of an imaginary part β in the fractional order can result in either pronounced backward propagation effects of elegant transient behaviour (Figure 6b, α = 1.75, β = 1), or a modulation of the initial propagation dynamics into intricate patterns of coexisting solitons and  (Figure 6b, α = 1.75, β = 2). For even smaller values of α, the reduced propagation speed (Figure 6c, α = 1.5, β = 0), combined with the highly oscillatory wavefronts, can produce even more elegant motifs with rich dynamics of transient growth ( Figure  6c, α = 1.5, β = 1), or influence the curvature of the generated patterns into labyrinthine, circuit-like configurations of the solutions (Figure 6c, α = 1.5, β = 2).
For a depletion rate κ = 0.061, standard diffusion in the Gray-Scott model yields curvature-driven dynamics of wavefront propagation (Figure 7a, α = 2, β = 0). However, under complex-order fractional diffusion, much more interesting phenomena arise. The localised oscillatory behaviour at the wavefront can induce further instabilities on its curvature, with solutions able to develop into dendritic-like structures (Figure 7a, α = 2, β = 0.75). Moreover, larger imaginary parts β can radically shift the dynamics of pattern formation into a mitosis process (Figure 7a, α = 2, β = 2.5), a dynamic behaviour previously only observed for model regimes at higher depletion rates [18]. Smaller values of α in real-order fractional diffusion result in the accentuation of the curvature-driven mechanism of pattern formation (Figure 7b

Discussion
In this work, we have established new and stable definitions for the complex-order fractional Laplacian as a representation of heterogeneous media with fractal microstructure, exploited their corresponding spectra to derive tailored Fourier methods for their numerical treatment, and investigated the dynamics of different systems of interest when driven by these operators. The proposed complex-order fractional Laplacian has been shown to encapsulate a dual particle-wave nature, with solutions manifesting wave-like features depending on the magnitude of the imaginary part of the fractional order, and to give rise to unique spatio-temporal behaviours.
The proposed definitions of the complex-order fractional Laplacian are fully consistent with the theory of averaging of smooth functions over fractal sets [12][13][14], therefore providing a link between fractals and fractional integro-differential operators. Our analytical derivation of the discussed operators highlights that, when combining fractional Laplacians on complex-law order, their average must be posed in terms of purely imaginary-conjugate weights. This is the only means to guaranteeing that solutions to the underlying diffusion processes remain real, as well as the positiveness in spectra of the resulting operators. In contrast, former works where fractional Laplacians of complex-law order were combined using real numbers result in an ill-conditioned behaviour when either increasing their imaginary parts or increasing mesh resolution [15][16][17]. A secondary benefit of the definitions of the complex-order fractional Laplacian presented in this work is the appearance of an additional amplitude parameter A (0 ≤ A ≤ 1), which provides further control of the periodic-logarithmic modulation exerted by the imaginary part of the fractional complex order on the diffusion eigenvalues.
By directly leveraging the spectra of the complex-order fractional Laplacian, and extending previous work [18], we have been able to introduce accurate and efficient Fourier spectral methods for the numerical simulation of complex-order fractional reactiondiffusion systems. These numerical schemes have been shown to be unconditionally stable regarding mesh resolution, and to achieve spectral convergence up to machine precision even under the presence of non-smooth initial data and highly oscillatory solution regimes. A further advantage of the proposed numerical stencils is their straightforward and efficient extension to higher dimensions, supported by a fully diagonal representation of the fractional operator in Fourier space (which avoids solving expensive linear systems per time step, as well as the need of preconditioners), and the exploitation of well-established Simulation results of the complex-order fractional Allen-Cahn equation display much richer dynamics of phase transition than its standard diffusion or real-order fractional diffusion counterparts. Perhaps the most interesting of the observed behaviours is the equilibrium of diffusion by interference of the waves localised at the phase boundary. Such an interference can yield a variety of steady-state structures solely based on the imaginary part of the fractional complex order, which characterises the fractal structure of the medium. We establish an analogy of this behaviour with the phenomena of Anderson localisation in non-linear optics (immobilisation in the diffusion of light in disordered materials due to interference between multiple scattered waves [1,22,23,38]), or strong localisation in condensed matter physics (absence of diffusion due to scattering processes in disordered materials [33,34]). Acoustic wave localisation of phase transitions has also been studied in viscoelastic media with random microstructures [39][40][41]. Besides localisation, the appearance of interfacial oscillations or ridges (known as periodic precipitation or the Liesegang phenomenon of pattern formation) is of interest on its own in the manufacturing of microstructures and nanostructures via reaction-diffusion processes [8,[42][43][44]. The Liesegang pattern of periodic precipitation typically arises in the reaction of two or more water-soluble substances diffusing in a solid hydrogel [8,44], and the role of heterogeneities in the substrate has been long acknowledged to markedly affect the obtained precipitation patterns [45,46]. Furthermore, while existing mathematical models have been shown to capture specific components of the experimentally observed precipitation patterns (such as nucleation, aggregation, or phase separation), a unified general model that collectively considers all these aspects is still lacking [43,44]. We believe that the wave-like nature encoded in the proposed definitions of the complex-order fractional Laplacian, coupled to the simulation of more comprehensive reaction models of solute precipitation, holds the potential to further advance the link between these reaction-diffusion phenomena and their underlying chemical model systems.
Solutions to the complex-order fractional FitzHugh-Nagumo model exhibit a variety of dynamic properties relevant to impulse propagation in heterogeneous biological tissues. In particular, slow and discontinuous conduction, conduction block, wavefront fragmentation, and break-up into self-sustained re-entrant wavelets are all characteristics of electrical conduction under patterns of microfibrosis [25,47]. Importantly, they are all acknowledged as major mechanisms for the onset and maintenance of deadly arrhythmias in structural heart disease. Computational studies of cardiac microfibrosis almost ubiquitously represent this type of substrates as small non-conductive infiltrations [25,[48][49][50][51][52], which translates into extremely fine mesh discretisations in order to capture the heterogeneities, and an associated increase in computational cost. On the other hand, the complex-order fractional operators presented in this work are capable of embedding such effects directly into a fractal nature of the operator itself, while maintaining a regular discretisation of the tissue. In fact, a fractal nature of fibrosis has been reported, both in cardiac [53,54] and other tissue types [55]. Other fractal attributes of cardiac tissue include its coronary vasculature and microvasculature, fractal morphology of microtubules and cytoskeleton in cardiac myocytes, or abnormal tissue growth in inherited cardiomyopathies, among others [56,57]. Moreover, the use of anomalous diffusion has gained a significant momentum over recent years for the description of cardiac tissue [15][16][17][18]26,28,[58][59][60][61], further supported by experimental confirmation of its findings [62,63]. In addition, previous studies have hypothesised the use of fractional calculus for the modelling of fibrotic tissues [64]. Altogether, our proposed operators and findings, enabling fractional diffusion operators as descriptions of discontinuous conduction and fibrosis, open new opportunities to investigate these topics from the lenses of a new modelling perspective. They are equally generalisable to other biological tissues characterised by anomalous diffusion behaviour, such as impulse propagation in nerve cells [65][66][67] or diffusion of species in the gastrointestinal system [68].
Finally, the above-discussed wave-like nature of the complex-order fractional Laplacian has been shown to be capable of generating completely novel spatio-temporal dynamics and Turing patterns in the complex-order fractional Gray-Scott model when compared to its standard diffusion or real-order fractional diffusion counterparts. These include intricate backward propagation effects, dendritic growth, multi-directional mitosis, or even the formation of star-shaped solitary structures, to cite a few. To the best of our knowledge, none of the above has been previously described in numerical simulation studies of pattern formation using the Gray-Scott model, not even under more pronounced regimes of realorder fractional diffusion [18,[35][36][37]69]. We anticipate that many other intriguing dynamics remain to be discovered in the complex-order fractional version of this system, either for different combinations of feed and depletion rates, or especially under unequal complex fractional orders for each of its diffusive species. Similar to the Gray-Scott model, many other reaction-diffusion systems have been shown to develop distinct dynamics of pattern formation when driven by real-order fractional operators [69][70][71][72][73][74][75]. We expect the dynamics of these and other systems of morphogenesis to also become significantly modulated by the wave-like nature of the complex-order fractional Laplacian, opening new avenues for investigation in the field of dynamical systems.

Conclusions
In this work, we have established new and stable definitions for the complex-order fractional Laplacian as a representation of heterogeneous media with fractal microstructure. The imaginary part of the fractional complex order induces a periodic-logarithmic modulation of the diffusion eigenvalues, which translates into a dual particle-wave behaviour in transport phenomena. Numerical simulation studies of different reaction-diffusion systems of general interest, by means of tailored Fourier spectral methods, demonstrate that these systems can exhibit marked wave-like features and novel spatio-temporal dynamics when driven by the proposed complex-order fractional Laplacian. Taken together, our methods and results support the thesis that the proposed complex-order operators hold unparalleled capabilities to advance the description of multiscale transport phenomena in physical and biological processes that are highly influenced by the heterogeneity of complex media.