Two-scale, non-local diffusion in homogenised heterogeneous media

We study how and to what extent the existence of non-local diffusion affects the transport of chemical species in a composite medium. For our purposes, we prescribe the mass flux to obey a two-scale, non-local constitutive law featuring derivatives of fractional order, and we employ the asymptotic homogenisation technique to obtain an overall description of the species’ evolution. As a result, the non-local effects at the micro-scale are ciphered in the effective diffusivity, while at the macro-scale the homogenised problem features an integro-differential equation of fractional type. In particular, we prove that in the limit case in which the non-local interactions are neglected, classical results of asymptotic homogenisation theory are re-obtained. Finally, we perform numerical simulations to show the impact of the fractional approach on the overall diffusion of species in a composite medium. To this end, we consider two simplified benchmark problems, and report some details of the numerical schemes based on finite element methods.

homogenisation methods permit to decouple the structural characteristic length scales [2,10,40,42,55,57], and, in particular, the asymptotic homogenisation technique [15,17,27,74] makes use of multiple scale expansions of the unknown fields to acquire an effective description of the medium at its coarser scales.
Usually, local constitutive laws are adopted in the description of the constituents of composite materials or heterogeneous media, thereby leading, in the majority of cases, to a homogenised local response of the composite. However, in certain circumstances, this has conducted to discrepancies with the experimental studies on heterogeneous media whose macro-scale properties could be more appropriately modelled by constitutive laws of non-local type. For instance, according to the experiments performed in [18,36,45,48], a medium with an involved structure may develop spontaneous restrictions on the way in which the transport processes occurring inside it are to take place. In turn, these restrictions, dictated by the medium's internal geometry, apply to the various length scales characteristing the transport processes and may result into nonlocal diffusion, thereby yielding non-Fickian diffusion [18,21,29,48]. Furthermore, it is worth noticing that, although the response of the constituents of a composite is often taken to be of local type at the lowest scale, in some cases, non-locality in time or space may arise as a result of homogenisation processes [9,26,37], or even by the adoption of standard concepts of solid mechanics [78], without having recourse to homogenisation techniques. For instance, as shown in [19], viscoelasticity can be obtained from suitable upscaling of a fluidstructure interaction problem between an elastic medium and a Newtonian fluid.

Scopes and novelties
To the best of our knowledge, there exist few works in which the constitutive laws of the constituents of composite media are assumed to be non-local already at the lower scales [20,79]. For instance, in [20], the homogenised properties of thermoelastic composites are studied by considering non-local integral operators for the characterisation of the stress-strain constitutive laws. In [20], the Author motivates the need for this constitutive choice by relating it to the complex geometrical and physical connections among the spatial length-scale of real materials, as is the case of hierarchical composite media [46].
In this work, motivated by the interest towards diffusion in highly heterogeneous biological media, we study a problem in which the mass flux of a chemical species in each constituent is related to the species' concentration gradient by means of a spatially non-local constitutive law. We do this by admitting that, in principle, two types of non-locality coexist: one pertains to the micro-scale and is thus associated with each constituent of the composite medium, while the other one is introduced a priori to allow for a non-local behaviour at the macro-scale. Note that these two types of non-locality are, in general, independent of each other. In the following, we will address each of these two non-localities through dedicated benchmark problems (see Sect. 5).
Several types of non-locality can be accounted for. For instance, as mentioned in [67], one can introduce higher-order gradients or integro-differential relations in the constitutive laws [4,8,14,31,34,38,47]. In the present work, continuing the research lines initiated in [67], we exploit Fractional Calculus [7,23,64] to describe diffusion processes that may deviate from Fick's law because of possible non-local behaviours in space. This modelling choice is motivated by the "success" of Fractional Calculus in addressing such phenomena [24,25,52,54]. In doing this, we have taken inspiration from the works [22,67,75].
In an attempt to realise how and to what extent non-local diffusion may affect the overall evolution of a given chemical substance in a composite medium, we prescribe a two-scale, non-local constitutive law for the mass flux of the considered substance. With this aim, we consider the asymptotic homogenisation technique to determine the effective diffusivity and the macroscopic evolution of the species. In doing this, we end up with an effective characterisation of the composite that is subjected to the existence of non-local interactions at both the micro-and the macro-scale. This is one of the novelties of this work.
As remarked in [67], the numerics of fractional diffusion in bounded domains requires special care because of the way in which the integro-differential operators featuring in the constitutive laws are to be handled, e.g. within FE methods. These difficulties increase if the medium in which fractional diffusion takes place is heterogeneous, as is the case in this work. A standard way of addressing numerically fractional differential equations in bounded domains is to have recourse to finite differences, specifically in the form of Grünwald-Letnikov schemes (see, e.g. [50,53]), although we are aware of works in which FE procedures are adopted [35,43,67,71]. However, this is not done for fractional differential equations in a multi-scale context, at least to our knowledge.
Another novelty of our work is that we address the FE discretisation of benchmark problems that include the non-local nature of the mass flux and the role of the heterogeneous structure of the medium under study. In this respect, we report the most important aspects of the FE discretisation and discuss some of the properties of the stiffness matrix and nodal force vector that feature in the resulting algebraic equations.
In summary, the main novelties of this work are: (i) The establishment of a combined framework in which some constitutive laws involving fractional derivatives are studied in conjunction with asymptotic homogenisation, in order to solve problems characterised by non-local diffusion at different scales. (ii) The realisation of a numerical scheme capable of putting together FE techniques with the integrodifferential operators from Fractional Calculus.
Our main result is the quantification of the impact of the spatially non-local diffusion of fractional type of chemical substances, resolved at the macro-and at the micro-scale of a strongly heterogeneous composite medium, on the transport of such substances within the medium. By employing asymptotic homogenisation to determine the effective diffusivity coefficient of the considered medium, our computations predict that, at the macro-scale, the attainment of the stationary state of the diffusion process is appreciably hindered by the nonlocal interactions accounted for by the operators of fractional differentiation that define the diffusion fluxes. This retardation manifests itself for decreasing values of a real parameter (the fractional order of differentiation) that defines the strength of such non-local interactions, at the micro-and at the macro-scale.
Finally, we mention that, in this work, we adopt a formalism that can be easily adapted to a two-or threedimensional context. However, since it can be challenging to study non-local phenomena in composite media via asymptotic homogenisation, we prefer, for the time being, to contextualise our mathematical model in a one-dimensional framework. Moreover, we remark that the results presented hereafter can be adapted in a straightforward manner to the study of thermal diffusion. Still, in the sequel, we shall only discuss diffusion of chemical species because the main problems that we have in mind come from the transport of chemical species in biological tissues.

Organisation of our work
The manuscript is organised as follows: in Sect. 2, some aspects of the topology of the composite are discussed, and we introduce the multi-scale governing equations describing the non-local diffusion of the chemical species. In Sect. 3, we consider the separation of scales between the macro-and the micro-scale, we illustrate the topology of the micro-structure, and discuss some aspects regarding periodicity in a two-scale context. Additionally, we reformulate the original governing equations to account for the two-scale nature of the non-local phenomena. In Sect. 4, the main mathematical tools of the asymptotic homogenisation technique are introduced, and we derive the effective properties and the homogenised equations for the composite under study. Furthermore, in Sect. 5, we specialise the general theory by presenting two different benchmark problems and discuss the results produced by numerical simulations. Moreover, we provide some details on the numerical schemes based on FE methods. Finally, in Sect. 6, we highlight the main results and outline some future developments.
form the periodic structure of B and, for every i, each pair of intervals ]X 2i , X 2i+1 [ and ]X 2i+1 , X 2i+2 [ represents two different constituents of the composite B. Moreover, it holds that B = B 1 ∪ B 2 and B 1 ∩ B 2 = B 1 ∩ B 2 = ∅, where the bar symbol indicates the closure of the set. In addition, we use the notation I to specify the interface separating the constituents B 1 and B 2 , namely Fig. 1).

Diffusion of chemical species
The diffusion of a chemical species in the composite B is described by with {X j = X 2i+1 } N i=0 ⊂ I, together with suitable initial and boundary conditions. Note that for ease of exposition these conditions will be specified later, when the benchmark problems are presented.
Equations (1b) and (1c) describe the contact on I, which in this case is assumed to be ideal, and the operator Φ(X j , t) denotes the jump of Φ across the interface I, i.e.
Moreover, Q denotes the mass flux of the chemical species and, as done in [67], we propose to express it in terms of the following non-local constitutive law, where D(X,X ) is referred to as non-local diffusivity, and is written as the product of the scalar quantity F(X −X ) and the fractional diffusivity D(X,X ), both taken to be strictly positive. We emphasise that F is defined for X =X , and that both D and D have, in general, physical dimensions different from those of standard diffusivity, depending on the prescription on F. Additionally, C and Q are continuous in B, which means that they are prolonged at the interfaces. It is worth noticing that further generalisations to the study of transport processes, involving for instance Darcy's law, can be found, e.g. in [1]. This work, however, pursues goals different from ours, since it considers constitutive laws that relate the time fractional derivative of the mass flux with the time fractional derivative of the classical pressure gradient. On the other hand, a one-dimensional diffusion problem in a bounded homogeneous medium is studied in [76], wherein Darcy's equation is generalised with a fractional integral in space. Furthermore, in the context of hierarchical materials, such as bones and ligaments, a generalised viscoelastic approach has been proposed to describe their rheological properties by using fractional derivatives and integrals [3,30], while numerical methods have been developed for the case of hereditary-ageing materials in [16]. Additionally, we notice that in [81] the analytical and numerical solution of a generalised heat conduction equation was studied by considering a fractional time derivative instead of the first-order partial time derivative of the temperature. Moreover, in [5], the authors considered a model in which, in addition to the fractional derivative in time, the heat conduction equation in a homogeneous material is extended by replacing the classical gradient of the temperature with its symmetrised Caputo fractional derivative. Finally, we point out that, in the context of viscoelastic composites, the Rabotnov exponential kernel [66], which is employed to construct a type of fractional derivative, has been considered in [70].

Separation of scales
In the characterisation of the two-scale nature of the composite, we assume the existence of two characteristic length scales, associated with the composite as a whole and its internal structure. Specifically, for our purposes, we denote by L c and the characteristic length scales of the composite medium and of its internal structure, respectively. Moreover, we require that the considered length scales are well separated by enforcing that /L c 1. Therefore, we introduce the dimensionless, smallness parameter ε, referred to as the scaling parameter, which is defined as the ratio We notice that ε characterises the heterogeneity of the composite, and permits to explicitly specify the twoscale nature of a given physical quantity Φ : B×[0, t f [→ R. In fact, following the discussion given in [32,63], one can take into account the multi-scale character of Φ(X, t) by rewriting it as Φ(X, t) =Φ(X, t; , L c ). As a particular case of this writing, we can impose thatΦ(X, t; , L c ) =Φ(X/L c , X/ , t), so that the dependence on the characteristic length scales is explicit. In this way, we have that where the dimensionless variables x := X/L c and y := x/ε are referred to as the macroscopic, or slow, variable, and the microscopic, or fast, variable, respectively. Note that, within this non-dimensional setting, B becomes X := ]0, L/L c [ and, accordingly, the non-dimensional variables x and y vary in X and in X/ε = ]0, L/ [=]0, 1 ε L/L c [, respectively. As stated in Eq. (5), one is able to express Φ as a function of two formally independent variables, thereby distinguishing the two scales characterising its nature. This means that, for every time t, the newly introduced function φ is defined, in general, as φ(·, ·, t) : Finally, we note that, by using the representation (5) and employing the chain rule, we can write 3.  Fig. 1). Furthermore, we consider that Here, for the sake of simplicity, we adopt the assumption of macroscopic uniformity [41,60,61]. This choice allows to choose the elementary cell, Y, independently of the macroscopic variable x, so thatY is representative of the composite's micro-structure (see Fig. 1). Moreover, for the type of functions φ(·, ·, t) : D x × D y → R used in the forthcoming calculations, we assume Y \ {y I } ⊂ D y and the existence of the lateral limits lim y→1 ± φ(x, y, t) and lim y→0 + φ(x, y, t). In the sequel, this property will be used to formalise the periodicity of φ with respect to its microscopic variable (this will be referred to as Y-periodicity), especially in the case in which D y has the form where N is a sufficiently large natural number. These considerations imply that it is sufficient to reformulate the problem at hand in the reference cell Y = ]0, 1[, along with the lateral limits outlined above, although for some physical quantities y I does not belong to the set in which they can be evaluated. In addition, since Y is chosen independently of the macroscopic variable x, also the following relation holds In general, however, if the hypothesis of macroscopic uniformity is not valid, the topology and geometry of the reference cell, Y, could vary with respect to the macroscopic spatial variable x and, thus, the reference cell should be regarded as a function of x, Y(x). In this case, Reynolds' transport theorem prescribes to rewrite the derivative of the left-hand side of Eq. (9) as where w n (x, y) is the normal "velocity" with which the boundary of the cell varies (see, e.g. [33,63] and the references therein for more details).

Periodicity
From the point of view of the small characteristic length scale , the body B can be approximated as unbounded, so that one can assume B = R. Within this approximation, a function Φ is said to be -periodic in the sense that Φ(X, t) = Φ(X + p , t), for all p ∈ Z, provided X and X + p are points in which the function can be evaluated [27]. Within the context of asymptotic homogenisation, one rephrases the periodicity of Φ in terms of the periodicity of the corresponding function φ with respect to the microscopic variable y. To this end, and to account for the fact that φ may be undefined for some values of y, it is necessary to express the periodicity of φ in the weaker sense supplied by with φ(x, y ± * , t) = lim y→y ± * φ(x, y, t), for all y * for which both lateral limits exist. This picture is consistent with the case in which φ(x, ·, t) is defined in a set D y of the type specified in (8) and y * is either p or p + y I , with p = 1, . . . , N − 2. In particular, the case y * = p + y I is important for performing the continuous prolongation at the interface of those physical quantities that have to be continuous at this point (for instance, the fluxes).
As anticipated above, the macroscopic uniformity, along with the Y-periodicity of the functions of interest for the problem at hand, enable us to restrict a given physical quantity to a single cell. For this purpose, one may choose the reference cell Y = ]0, 1[ , and take the restriction φ(x, · , t) |Y . Furthermore, to account for the presence of the interface, which splits the cell in the disjoint union of two materials with different properties, we define φ(x, · , t) |Y as the piecewise function In particular, to describe the periodicity at y * = 0, we invoke Eq. (11), so that, Granted this result, we notice that, if φ is a function for which the continuity condition at the boundary of the periodic cell must be respected (for example, a concentration or mass flux), we also find and, thus, because of periodicity,

Multi-scale non-local diffusion
Upon adopting the above considerations, and recalling the identities (5) and (6), we rephrase the original problem (1a)-(1c) as follows (see Remark 1 for further details), with (x, y, t) ∈ X × Y k ×]0, t f [, and where the index k ∈ {1, 2} indicates in which sub-cell Eq. (16a) and the quantities c k and q α,β k are defined. Particularly, the two-scale, non-local flux q α,β k is given by where d α,β k (x,x, y,ỹ) represents a two-scale "version" of the non-local diffusivity coefficient D(X,X ). More precisely, by using Eq. (5), we have that D(X,X ) is replaced by d α,β k (x,x, y,ỹ), which means that the parameters α and β k (see below) are already present in D(X,X ). Analogously, f α,β k (x −x, y −ỹ) and d α,β k (x,x, y,ỹ) replace F(X −X ) and D(X,X ) in the decomposition (3b), but describe the non-locality and the fractional diffusivity resolved on the two different scales accounted for in this work. Particularly, α ∈ R + is referred to as the macro-scale non-locality parameter and characterises the non-local interactions in the region X. On the other hand, β k ∈ R + , with k = 1, 2, is the micro-scale non-locality parameter describing the non-locality within the sub-cell Y k . Note that q α,β k absorbs the factor 1/L c that stems from the chain rule (6) when one switches to the two-scale representation of the flux.

Remark 1
The representation of the two-scale, non-local mass flux in Eqs. (17a) and (17b) does not follow directly from (3a). This is because the double integral over X × Y k defining q α,β k cannot be obtained by only applying the two-scale representation prescribed by (5) and (6) to the integrand of (3a). Rather, to account for the two-scale resolution of the flux, a further step is needed, which requires to pass from a single integration in the variableX to a double integration in the two auxiliary variablesx andỹ. In this respect, it must be clearly stated that the flux q α,β k is not equal to Q, and it is introduced ad hoc as a mathematical tool with the purpose of resolving the two-scale dependence of the original flux. Hence, the definition of q α,β k must be regarded as a conjecture, which in the limit ε → 0, and within the asymptotic homogenisation approach, converges to an effective flux that represents the limit of Q [refer to Eq. (45b)]. Proving this rigorously is part of our current investigations, which involve, among others, the concept of two-scale convergence [27,58,80].
We remark that the introduction of the non-local parameters α and β k follows from the fact that we interpret non-local effects by using the notions of Fractional Calculus [7], in which derivatives and integrals of fractional order are considered. The parameter α accounts for the intensity of non-locality at the macro-scale, whereas we have intentionally introduced two different non-locality parameters, β 1 and β 2 , at the microscopic level to describe the existence of "long-range" interactions even at the scale of each sub-cell Y k . This is indeed the essence of the micro-scale non-locality.
We further notice that, if the concentration c k is dimensionless, the flux q α,β k must have the physical dimensions of the reciprocal of time, and it follows from Eq. (17b) that the same must be true for the dimensions of d α,β k . To guarantee the latter condition, and keeping in mind that d α,β k is a (fractional) diffusivity, we take the dimensions of d α,β k to be [d α,β k ] = length ξ(α,β k ) /time, where ξ(α, β k ) is a real number expressed as a function of α and β k , and, consequently, we take [f α,β k ] = length −ξ(α,β k ) . In the local case, the non-locality function may be taken dimensionless, which means that ξ(α, β k ) must tend towards zero, and the fractional diffusivity becomes a pure rate, that is [d α,β k ] = time −1 .
In the sequel, we assume that every field is periodic with respect to the micro-scale variable y. Moreover, the spatial fractional diffusivity d α,β k is considered to be independent of x,x and y, and with a slight abuse of notation, we simply write d α,β k (ỹ). This simplification, however, has not major repercussions in the results of the following sections.

Two-scale asymptotic homogenisation approach
In this work, we adopt the asymptotic homogenisation technique and prescribe a formal two-scale expansion for c k in power series of the smallness parameter ε > 0, namely where each c (n) k (x, ·, t), n = 0, 1, 2, . . ., is assumed to be periodic with respect to y (said more rigorously: periodic with respect to its second argument).
Before substituting the formal expansion (18) into (16a)-(16c), we find it convenient to rewrite Eqs. (16a)-(16c) as follows where the following notation has been adopted Specifically, in (20a) and (20b), the uppercase and lowercase symbols Q α,β k and q α,β k indicate the partial differentiation of c k inside the integral with respect tox andỹ, respectively. Moreover, it holds that q α,β k = Q α,β k + ε −1 q α,β k . After substituting (18), truncated to the order ε 2 , into (19a)-(19c), (20a) and (20b), the problem reduces to finding the leading order coefficients c (n) k of the power series (18), which solve the boundary problems resulting from equating all the terms in the same powers of ε. To this end, it is useful to write explicitly the generic coefficients of the expansion of the fluxes q α,β k and Q α,β k , i.e.
In the sequel, to avoid the proliferation of indices, we simplify the notation as follows (i) To the order ε 0 , Equation (24a) implies that q (0) k is independent of the microscopic variable, since its partial derivative with respect to y is zero. One possible way of ensuring this condition could be to drop the dependence of f k on the micro-scale variables. However, this assumption would eliminate the possibility of keeping track of the nonlocality at the micro-scale, which is clearly in contrast with our purposes. Instead, to guarantee the fulfilment of Eq. (24a) and to make sure we remain within a non-local setting, we require c (0) k to be independent of y. Hence, with a slight abuse of notation, we set thereby satisfying Eqs. (24a)-(24c), since q (0) k = 0, without having to rebut the dependence of f k on the micro-scale. The above consideration is a standard result of linear asymptotic homogenisation, whereas it is often assumed for nonlinear problems (see, e.g. [28,65,68]). A direct consequence of (25) is that c (ii) To the order ε 1 By taking into consideration Eq. (25), we have that with (x, y, t) ∈ X × Y k ×]0, t f [. The structure of (27a) implies that, in general, this equation must be solved in the spatial domain X × Y k , which in the context of asymptotic homogenisation means that one needs to consider a problem defined in Y k for each x ∈ X.
Remark 2 (A comment on the solution to (27a)-(27c)) The local counterpart of the problem (27a)-(27c) (see [15,27]) can be obtained by choosing the non-locality function where δ is Dirac's delta. Specifically, with (x, y, t) ∈ X × Y k ×]0, t f [. In Eqs. (29b) and (29c) the evaluation in y I of d k and ∂ y c (1) k are to be understood in the sense of lateral limits y → y ± I . In this particular case, the problem (29a)-(29c) admits a unique solution, which is defined up to a function depending solely on time, t, and on the slow variable, x [15,27]. This unique solution is usually expressed through the ansatz where ϑ k is the new unknown of the problem (29a)-(29c) and ϕ is a function of x and t that spans the family of all the solutions [11,63]. To the best of our knowledge, in the non-local case there is no theorem that guarantees the existence and uniqueness (even in the sense explained above) of the solution. Still, in the absence of a supporting theory, we guess that, similarly to the local case, the solution should have the form (30), with ϑ k suitably parametrised by α and β k .
By substituting (30) into (27a)-(27b), we require the auxiliary functions ϑ k to satisfy the non-local cell problem, We notice that the structure of the non-local problem (31a)-(31c) does not permit, in general, to factorise the macroscopic term ∂x c (0) (x, t). This implies that one should account for the macroscopic contributions at the micro-structural level and, thus, for the interchange of information between the two length scales.
(iii) To the order ε 2 , Before going further in our analysis, we introduce, for a given field φ, defined in the cell Y or in a subset of it having the same measure, the operators such that the sum φ 1 + φ 2 = φ is the average of φ over the cell Y. Then, by applying these operators to (33a), we have Because of Eq. (26), c (0) depends only on x and t and then, Moreover, the assumption of macroscopic uniformity (see Sect. 3.2) implies that the differential operator ∂ x and the integral operator · k commute, so that, the second term of (35) rewrites as, Therefore, summing up Eq. (35) over k and taking into account the relations (36) and (37), we obtain We notice that the third term of (38) can be computed as where we have employed Gauss' theorem and the continuity of the fluxes at the interface and at the boundaries of the cell. Specifically, because of the continuity of the fluxes at the interface y I , it holds true that which eliminates the first and the fourth summand on the far right-hand side of Eq. (39). Moreover, the flux computed at the right boundary of the cell, i.e. q (2) , must be equal to the flux entering or leaving the neighbouring cell, which can be written as q (2) 1 Therefore, by invoking the Y-periodicity of the flux, we can conclude that the second and the third term of Eq. (39) also cancel themselves, i.e.
Equations (40) and (41) explain in detail the reason why Eq. (39) holds true. Before going further, we emphasise that the considerations done so far hold true also for all the other orders of the asymptotic expansion of the flux, like, for instance, q (1) k + Q (0) k . Then, the substitution of (39) into (38) yields the homogenised problem for the leading order term c (0) , i.e. where is referred to as the non-local effective mass flux, while d eff is defined through the expression and represents the non-local effective diffusivity. We notice that the homogenised equation (42) has the same structure as (1a), but, in this case, the contributions of the micro-structure are resolved by means of the non-local effective coefficient d eff .
Finally, according to [63,68] we introduce the notation X h to denote the homogenised version of the composite medium, and we reformulate the homogenised problem (42) and (43) as follows which has to be supplemented with appropriate initial and boundary conditions for the unknown c (0) .

Remark 3
Note that, since f k and d k are short-hand notations for f α,β k and d α,β k , both the non-local effective diffusivity, d eff , and the non-local effective flux, q eff , depend on the collection of all the parameters that describe the non-locality of the problem, i.e. α, β 1 and β 2 . Hence, the effective quantities d eff and q eff keep track simultaneously of the non-locality occurring both at the scale of the sub-cells, through β 1 and β 2 , and at the scale of the medium, through α. In the following, with the purpose of leaving the notation at a minimum level of complexity, we shall keep the symbols d eff and q eff , although we mean d eff α,β 1 ,β 2 and q eff α,β 1 ,β 2 , respectively.

Remark 4
In summary, the equations to be solved are (31a) for ϑ k and (45a) for c (0) . We notice that, in general, these equations are coupled in contrast to the local case of the standard asymptotic homogenisation in which the cell and the homogenised problems are uncoupled. As anticipated above, this is due to the fact that ∂ x c (0) (x, t) cannot be factorised out the integral in Eq. (32a).

Remark 5 It is worth noticing that if the non-locality function is
, we end up with classical results of homogenisation theory [15,27,63,69]. That is, by substituting this expression for f k into Eq. (31a), the cell problem reads where d y ϑ k (y) denotes the total derivative of ϑ k . Furthermore, the non-local effective diffusivity, d eff , is withd eff st being entirely defined by the sum over k in (47). In this case,d eff st is a constant coefficient that coincides with the effective diffusivity of a standard diffusion problem in a composite medium [15,51,56]. Furthermore, after substitution of (47) into (45a), we obtain the standard homogenised equations where standard Fick's law is re-obtained for the flux.
Remark 6 We notice that in the present framework, we do not take into account the timescales associated with the problems (24a), (27a) and (33a), since we intend to focus on the spatial connections between heterogeneity and non-locality. Nevertheless, one can notice that the characteristic length scales L c and associated with the composite medium and with its internal structure, respectively, induce different timescales (see, e.g. [49]). In fact, by virtue of a reference diffusivity d R , these can be expressed as Since in the sequel we specialise Eq. (49) to the case of media with different diffusivities inside the sub-cells By employing (49), we deduce the following relationship between the characteristic time scales, Now, before proceeding further, we mention that in this multi-scale framework, a given physical quantity Φ(X, t) can be rewritten as with ζ := t/ζ c and η := ζ /ε 2 (compare Eq. (51) with Eq. (5) in which time was not rescaled). Therefore, which, after substituting the two-scale expansion (18) and equating in the same powers of ε, up to the order ε 2 , yields As Eqs. (53a)-(53c) prescribe, this approach calls for the solution of diffusion problems at each order of ε. Moreover, the consideration of the separation of the time scales conduces to leading-order problems that are characterised by the presence of the rapid time variable η. The analysis of these equations is part of our current research.

Description of the benchmark problems
In this section, we present two simplified models in which the non-local effects are only present at the macroscale or at the micro-scale, and we report some details of the numerical schemes based on FE methods. It is worth mentioning that the cell and the homogenised problems obtained in this work feature integro-differential equations of fractional type in bounded domains and, therefore, the classical solution techniques, such as Laplace and Fourier transforms, used in Fractional Calculus are not suitable. Consequently, either do we need to develop dedicated numerical algorithms or we resort to well-established numerical methods, and we adapt them to take into account the presence of fractional differential operators in the considered problems. Here, we follow the second path and, indeed, we write a numerical scheme based on a finite element (FE) discretisation of the original integro-differential problems. In doing this, we need to emphasise that, partly because of the very easy geometry of the problems (we deal, in fact, with one-dimensional benchmark studies), and partly because the focus of our work is not on the numerics, the presentation of the FE scheme is very elementary. Indeed, it can be obtained by appropriately rephrasing the one-dimensional formulation of the FE method as presented, e.g. in [44]. Moreover, we do not fuss over some technical aspects of the finite element procedures, such as the "element point of view" [44] et similia, since our scope is solely meant to highlight how the symmetrised Caputo fractional derivatives affect the stiffness matrix and the nodal force of the discretisation. Clearly, a more detailed numerical study is required, and this is part of our current investigations. We highlight that previous works in this direction are [35,43,71]. In particular, the work we took major inspiration from is [43].
In the remainder of this work, the following considerations are adopted. (i) Fractional diffusivity. We prescribe d 1 and d 2 to be constant in Y 1 and Y 2 , respectively. Then, by recalling the discussion made in Sect. 3.4, we express each d k as where d kR is the constant reference diffusivity of Y k and it has the dimensions of a standard diffusivity, i.e. length squared over time [67].
(ii) Initial and boundary conditions for the homogenised equation. We enforce an initial spatial distribution for c (0) of the form where k, r and x 0 are model parameters.
To contextualise our work, we mention that the initial condition c in (x) in (55) is sometimes employed to simulate the initial concentration of molecules after photobleaching in a Fluorescence recovery after photobleaching (FRAP) experiment [13]. In this way, following [39], the model is prepared to describe the fluorescence recovery pattern of molecules surrounding a certain region of a tissue (in particular, articular cartilage) after being photobleached, by using a high-intensity laser beam. Here, we do not go into the technical details pertaining to a FRAP experiment, since this is not the focus of our work, and the benchmark proposed hereafter is also markedly different from the one developed in [39]. Thus, we do not claim that our results are meant to simulate a FRAP experiment. Still, since the setting presented in [39] refers to a tissue with hierarchical internal structure, as is the case of articular cartilage, our work might bring some new insight into the interpretation of the experimental results. To this end, we adapt the framework described in [39] to the setting of the homogenised problem (45a)-(45b) and, specifically, we identify k, r and x 0 with the bleaching depth parameter, the dimension of the bleached area, and the centre of the bleached region, respectively.
We notice that, since this work is framed in a one-dimensional setting, r characterises the measure of a line-segment of X h , and we choose x 0 as the centre of the macroscopic domain X h , namely x 0 = 1 2 (L/L c ). In addition, we adapt the boundary conditions given in [39] to the geometry of our problem, and impose Dirichlet boundary conditions for c (0) at x = 0 and x = L/L c . Specifically, we set which implies (iii) Parameters. In Table 1, we provide the values of the parameters used in our numerical simulations. We notice that the value of r is meant to "cover" 100 reference cells.

Benchmark problem I: micro-scale non-locality
Let us consider the case in which the non-locality is accounted for only at the micro-scale. This can be achieved by prescribing that is, we accept the existence of "long-range" interactions in each sub-cell Y k . Note that we use quotation marks because the concept of long-range interactions has to be understood with respect to each sub-cell, which is microscopic, and, in this context, as a synonym of non-locality. We also notice that the index k in g k allows to characterise two different non-local frameworks occurring in each sub-cell Y k . For instance, as discussed above, we could enforce that the non-local interactions exist only in one of the two sub-cells, although here we consider non-locality acting in both sub-cells. Clearly, different forms for g k can be considered, each of which leading to diverse non-local models of diffusion. In this work, we adopt the decaying power-law [6,22,59,67,76] where Γ (·) denotes the Euler Gamma function and β k ∈ ]0, 1[. From here on, we set β 1 = β 2 = β, thereby obtaining g 1 = g 2 = g and f 1 = f 2 = f. We notice that g scales multiplicatively with L 1−β c because it is expressed as a function of dimensionless variables. Accordingly, the physical dimensions of the fractional diffusivities

The non-local cell problem
By considering Eqs. (58) and (59), the non-local cell problem is given by where (x, y, t) ∈ X × Y k ×]0, t f [ and As shown in (62a) and (62b), the assumption of non-locality only at the micro-scale permits to factorise ∂ x c (0) from the integrals expressing q (1) k and Q (0) k . Consequently, the auxiliary unknowns ϑ 1 and ϑ 2 can be reformulated as functions of y only, and for further use, the following notation is introduced We notice that the quantities q k (y) and Q k (y) defined in (63a) and (63b) represent, up to the sign, "dressed" diffusivities rather than fluxes. In fact, we may write We further mention that, in the proper limit, q k (y) and Q k (y) return the negative of the standard diffusivities (see Remark 7). Particularly, in Eqs. (63a) and (63b), we recognise the symmetrised Caputo fractional derivative of order β ∈ ]0, 1[ [7]. Then, it holds true that where κ(y) = y, and denotes the symmetrised Caputo fractional derivative of order β of a generic differentiable function φ. Furthermore, we notice that Q k (y) can be computed explicitly for each k = 1, 2, and reads where the functions A 1 and A 2 are given by Because of these results, the flux (64b) admits the explicit expression while, for the time being, no explicit expression can be given to q k , since the functions ϑ k are still unknown. Upon recalling the definition (12), in which a given function restricted to the elementary cell is assigned in a piecewise manner, the fluxes Q (0) 1 and Q (0) 2 in (65) can be rejoined in a unique flux, whose restriction to Y is given by Now, for the function Q (0) given in (70), we can employ the definition of periodicity specified in (13), so that it holds It follows from this result that Q (0) is Y-periodic and that such periodicity does not depend on the point y I in which the interface is placed within the elementary cell. By using the above results, the non-local cell problem (61a)-(61c) can be rewritten as where the right-hand side of (72c) is the result of the difference Q 2 (y I )−Q 1 (y I ). We recall that all the expressions at the interface are to be understood for the values of the limits of the corresponding physical quantities for y → y ± I . Indeed, for instance, Q 1 (y I ) means, with a slight abuse of notation, Q 1 (y I ) = lim y→y − I Q 1 (y).

Remark 7
We notice that for y ∈ Y k , it holds Then, it follows from Eqs. (73)-(74b) that the convergence of Q k for β → 1 − is not uniform in Y k . By comparing with classical results of the asymptotic homogenisation technique, the above computations suggest that, for β → 1 − , the solution of the non-local cell problem must approach the solution of the classical, local cell problem. We notice that the 2 in the denominator of (74b) does not appear in the formulation of the standard cell problem. Nevertheless, it compensates with the 2 in the denominator of the left-hand side of (72c), hidden in the fractional derivatives defining q k , which can be determined after finding ϑ k .

The homogenised equation
By taking into account Eqs. (58) and (59), the non-local effective coefficient can be rewritten as and, therefore, according to Eq. (43), the effective flux is given by where we have setd Hence, the effective fractional diffusivity,d eff (β), can be expressed in terms of the symmetrised Caputo fractional derivative of ϑ k . We notice that, in this particular case,d eff (β) does not depend on space and time, while it is parametrised by β. Furthermore, from its mathematical expression it is clear that it ciphers the information on the non-local interactions at the micro-scale.

Remark 8
The form of the effective fractional diffusivity (77) recalls the relation obtained in the standard case by means of asymptotic homogenisation [12,15,69]. Particularly, by Fubini's theorem and Eq. (67),d eff (β) can be equivalently rewritten aŝ Therefore, taking into account Eq. (73), we obtain that which coincides with the effective diffusivity of a standard diffusion problem with unitary reference cell (see, e.g. [27]). Furthermore, since Eqs. (77) and (78) provide equivalent writings ford eff (β), from (79) it follows that, in the limit β → 1 − , the symmetrised Caputo fractional derivative of ϑ k converges to the first derivative of ϑ k , namely, and therefore, we can conclude that q (1) k (see Eq. (64a)) recovers Fick's law in bounded domains.

Numerical solution
In this section, we solve numerically the mathematical model given by the non-local cell problem (72a)-(72c) and the homogenised problem (81a)-(81c). In particular, the homogenised problem is characterised by a partial differential equation, while the non-local cell problem features an integro-differential equation of fractional type. For solving Eqs. (72a)-(72c) in a bounded domain, the techniques based on Fourier and Laplace transforms are of little help and, consequently, we implement a numerical algorithm. Specifically, we solve the non-local cell problem by means of a FE scheme which accounts for fractional derivatives and interface conditions. Since the homogenised problem (81a)-(81c) involves classical FE techniques, in what follows we report only on the numerical scheme to be used for solving the non-local cell problem (72a)-(72c), and for the computation of the non-local effective diffusivity (77).
Before going further, we notice that, in the classical homogenisation literature, the uniqueness of the solution of the cell problem is guaranteed by imposing that ϑ k k is equal to zero. However, from a computational point of view, a more feasible condition is to fix the value of the auxiliary variables ϑ k at one point in the cell [62]. Accordingly, here, we impose that ϑ 1 is zero at y = 0, and by periodicity ϑ 2 is also zero at y = 1. Now, let us introduce the following spaces of test functions derivatives up to order one [73]. Then, by multiplying Eq. (72a) by v k , integrating over Y k , and adding over k = 1, 2, we obtain Hence, due to the continuity condition at the interface (72c), and the restrictions made for v k , Eq. (83) reads Equation (84) represents the weak formulation of the non-local cell problem (72a)-(72c), and is discretised by employing the FE technique. Therefore, we introduce for Y 1 and Y 2 , N 1 +1 and (N 2 +1)− N 1 discretisation points, respectively, and the function bases with y N 1 = y I . Then, the test functions v k are approximated by v k (y) = where υ k j , k = 1, 2, are nonzero, arbitrary constants. Analogously, the approximated trial functionsθ k (y) are written as where ω k i are unknown constant coefficients representing the nodal values ofθ k . We notice that, in this particular case, the coefficients ω k i do not depend on time, whereas in a more general setting they should be defined as functions of time.
Next, by substituting expressions (87) and (86) into (84), we obtain the following system of equations for ω k i , where represent, for each k = 1, 2, the components of the fractional stiffness matrix of the FE discretisation of the sub-cell Y k , and of the nodal fractional force associated with the j-th node of Y k , respectively. Remark 9 (Density and limit of L k (β) and F k (β)) It is worth noting that, whereas in the standard cell problem the stiffness matrix is tridiagonal for each k = 1, 2, in the present framework it is dense because the cross integrations between the derivatives of the basis functions lead to nonzero components of for each pair of j and i. This is due to the non-locality introduced by the fractional derivatives D β k [ψ k i ], as reported in the far right-hand side of Eq. (90). Specifically, even though two discretisation nodes are far away from each other, the entries of the matrix corresponding to those nodes are nonzero. This results in a dense stiffness matrix and, the stronger the non-locality, the denser the matrix will be. Nevertheless, the fractional stiffness matrix will converge to a tridiagonal matrix when β → 1 − . Indeed, as discussed in Remark 8, when Analogously, from the definition of F k j (β), i.e. we infer that the existence of the non-locality function implies that the entries of F k j (β) are nonzero. Furthermore, recalling that lim β→1 Equation (93) returns 0 for all j = N 1 and d kR L −2 c if j = N 1 . To exemplify the limit of the symmetrised Caputo fractional derivative of the bases functions, we report in Fig. 2 a comparison of the symmetrised Caputo fractional derivative of order β ∈ ]0, 1[ of the function which recalls a Lagrange polynomial of the first order, with the classical first derivative of the same function.
Next, to obtain the algebraic form of the FE procedure, we introduce the notation where υ I N 1 = υ 1 N 1 = υ 2 N 1 and ω I N 1 = ω 1 N 1 = ω 2 N 1 are the nodal values of the virtual concentration and of the unknown concentration at the interface, and we write the final forms of the fractional stiffness matrix and of the fractional nodal force as follows Note that in (96a) and (96b), we have omitted the dependence on β , although this dependence is understood.
Then, by using the notation introduced in Eqs. (95a)-(96b), Eq. (88) can be rewritten as which leads to the algebraic equation On the other hand, by using the expression forθ k given in (87), the approximation of the effective fractional diffusivity (77) can be numerically calculated aŝ which we call numerical effective fractional diffusivity. We provide details about the explicit form of Eqs. (98) and (99) in "Appendix A".

Results and discussion
In this section, we show the numerical results for the benchmark problem described above, and we discuss the influence of the micro-scale non-local interactions on the homogenised behaviour of the concentration.
To begin with, in Fig. 3, we report the profile of the solution of the non-local cell problem (72a)-(72c), i.e. ϑ k , and compare it with the solution of the standard, local cell problem. Specifically, the solid lines distinguish the solutions of the non-local cell problem for different values of the non-locality parameter β ∈ ]0, 1[, and the dashed line represents the solution of the standard, local cell problem. In particular, the space discretisation of the computational domain was done by fixing the grid size to h := y i − y i−1 = 1.3 × 10 −3 uniformly with respect to i. We notice that the results of the finite element analysis are not affected appreciably by subsequent mesh refinements.
In Fig. 3, we observe that the spatial distribution ofθ k varies with β and it converges to the solution of the local cell problem as β → 1 − (dashed line in Fig. 3). This outcome is coherent with the theoretical results previously obtained in this section. Furthermore, we notice that the non-local solutions fluctuate around the local one, and they intersect each other and the local solution in symmetric points. Nevertheless, the non-local solutions are not symmetric with respect to the line y = y I .
Before going further, few words should be spent about the issue of symmetry. To this end, let us assume just for this discussion that the interface between the sub-cells, y I , is not the midpoint of the cell Y = ]0, 1[, and let us start with the local case (recovered for β → 1 − ). In the local case, the solution of the cell problem, here denoted by ϑ local and defined as ϑ local (y) = ϑ local 1 (y) for y ∈ [0, y I ] and ϑ local (y) = ϑ local 2 (y) for y ∈ ]y I , 1], is in general not symmetric because the diffusivity coefficients are distributed within the cell in a non-symmetric way (clearly, this asymmetry would disappear if the diffusivities were equal to each other). Within the framework studied here, the solution will have, indeed, the shape of a non-symmetric "roof", with an increasing straight line on [0, y I ] and a decreasing straight line on ]y I , 1], whose slopes have different signs and different absolute values. In other words, the inequality d 1R = d 2R and the position of the interface break the symmetry that the solution would have in the homogeneous case (the solution of the problem at hand would trivially boil down to a constant in the homogeneous case). Symmetry, however, can be partially restored if the interface is assumed to be placed at the midpoint, as is the case in our simulations. This condition, indeed, places a geometric constraint that forces the solution to be symmetric, thereby acquiring the shape of a symmetric "roof", with the slope of the straight line on [0, y I ] being the opposite of the slope of the straight line on ]y I , 1]. This means that we have passed from the continuous symmetry of the homogeneous solution to the discrete symmetry of the heterogeneous solution with interface in the midpoint of the cell.
The picture just described changes considerably when the non-local case is studied. Indeed, for β ∈ ]0, 1[, the fractional derivatives featuring in the non-local cell problem are an additional source of symmetry breaking that, together with the heterogeneity of the diffusivity, make the solution even more non-symmetric. This remains true even though y I is the midpoint of the cell. More importantly, since the local solution is symmetric  for y I = 1/2 in spite of the heterogeneity of the diffusivity, this setting singles out the contribution of the fractional derivatives to the symmetry breaking of the problem. We believe that the asymmetry of the non-local solution, which decreases with β, may be ascribable to an interplay between non-locality and heterogeneity: in the sub-cell in which the medium is more diffusive, the solution is "stiffer" which results into moderate deviations from the standard solution; on the other hand, in the sub-cell in which the diffusivity is smaller, the solution is more "compliant", thereby producing large deviations from the local solution. Now, onceθ k is known, we can compute the effective fractional diffusivityd eff num as prescribed by formula (99). Particularly, in Fig. 4, we plot the values of this homogenised coefficient for varying β ∈ [0, 1] and compare them with the classical effective diffusivity, i.e. the one resulting from the local case. Specifically, a closer look at the data reported in Fig. 4 reveals that, for increasing β, the value of the effective fractional diffusivity resulting from a non-local setting is higher. In particular, as discussed in Remark 8, as β tends to 1 from below, the approximated effective fractional diffusivity converges to the standard effective one given by the local case.
We notice that for β = 0 the auxiliary problem is ill-posed and, thus,θ k cannot be determined. This is also reflected by the fact that the stiffness matrix of the problem, L(0), becomes singular for β = 0, andθ k becomes non-differentiable at y = y I and at the boundaries of the cell. On the other hand, for β > 0, the gradients of ϑ k exist at these points but their magnitude increases for β → 0 + . Nevertheless, it is worth remarking that, for very small values of β, the numerical solution almost does not change. Particularly, the L ∞ -norm of the error between the numerical solutions for β = 10 −8 and β = 10 −3 is of the order of 10 −4 . In addition to these considerations, we would like to point out that neitherθ k nor its gradient are observable physical quantities. Rather,θ k (y) is just an auxiliary quantity for determining the observables q (1) k (x, y, t), q k (y), and, more importantly, the effective fractional diffusivity and the homogenised solution. To this end, we notice that, in fact,d eff (β) and c (0) are well behaved for all values of β ∈ ]0, 1[ as shown in Figs. 4 and 5, and also for β = 0. Specifically, in spite of the technical difficulty for β = 0, which makes the employment of the FE method impossible, it is still possible to determine the variations which return the value of ϑ k at the interface. This calculation allows us to compute, even in this limit case, the effective diffusivity coefficientd eff (β), which, as shown in Eq. (77), we rephrase as a function of β and, for Since, as per Fig. 4, which is the plot of Eq. (77),d eff (β) is a continuous and monotonically increasing function of β ∈ [0, 1], it occurs that the valued eff 0 represents the absolute minimum of the effective diffusivity coefficient, i.e.d eff 0 = min β∈[0,1] {d eff (β)}, and the absolute maximum isd eff (1) = max β∈[0,1] {d eff (β)}. The above result describes, for each β ∈ [0, 1], the influence of the micro-scale non-locality on the macroscopic distribution of the concentration c (0) (see Fig. 5). In particular, for β tending towards zero, i.e. for increasing "strength" of the micro-scale non-locality, the macro-scale diffusion of the considered substance is hindered, and c (0) (x, t) consistently tends to vary rather slowly in time. On the contrary, in the limit case β → 1 − , c (0) (x, t) varies more rapidly in time, since the diffusion tends to acquire the "classical" behaviour predicted by Fick's law (see Fig. 5). In this respect, the consideration of non-local interactions at the microscale influences the way in which diffusion takes place in the composite medium. Returning to the FRAP experiment in the context of the benchmark problem, this theoretical behaviour implies that the recovery pattern of chemical species after being photobleached is slower for β near 0, whereas it is faster for β close to 1, thereby simulating a standard diffusion process.

Benchmark problem II: macro-scale non-locality
In this section, we assume that the non-local interactions are present at the macro-scale only. Thus, by specialising f α,β k in (17b) to the limit case β k → 1 − , we consider the following form for the non-locality function Hence, since h α depends only on the difference x −x, the non-local character of the diffusion process is accounted for at the macroscopic level only, and similarly to what was done in the previous section, we write In this particular case, the physical dimensions of the fractional diffusivities d k are [d k ] = L −1+α c t −1 c and hence, from (54), we have that

The cell problem
By considering the expressions (102)-(104), the non-local cell problem (31a)-(31c) rewrites where In (106b), D α [c (0) ] represents the symmetrised Caputo fractional derivative of order α ∈ ]0, 1[ of c (0) . Particularly, the computational complexity of the above cell problem is significantly reduced if the solution ϑ k is x-constant (which in the present framework also implies that it is constant in time). Then, with a slight abuse of notation we write ϑ k (x, y, t) = ϑ k (y), and q (1) k in (106a) becomes while Eq. (105a) rewrites We notice that [6]. Therefore, by excluding this case, the cell problem can be written in the more standard form In this specific case, the analytical solution of the cell problem (109a)-(109c) can be found by using standard techniques for differential equations. However, since our scope is to find the effective coefficient, this is not necessary. Indeed, from (109a) we can deduce that where a k , with k = 1, 2, are two constants to be determined. In particular, substituting (110) in (109c) yields a 1 = a 2 ≡ a, and the constant a can be computed by invoking periodicity and (109b). In fact, from (110), it follows that and, by applying the operators defined in (34) to the last equation, we have which implies that Therefore, after substitution of (102) and (103) into Eq. (44), and using (110) and (113), the non-local effective coefficient can be computed as It is worth mentioning that, even though, in this particular formulation, the cell and the homogenised problems have been decoupled, the non-local effective diffusivity (114) is still influenced by the non-local interactions occurring at the macroscopic level through the scalar function |x −x| −α . Note also that the last two factors of d eff (x,x) in (114) are the kernel of the operator defined in (103).

The homogenised equation
By using the previous results, the effective non-local mass flux can be recast in the form which is thus entirely determined by the symmetrised Caputo fractional derivative of order α of the leading order concentration c (0) and by the effective diffusivity coefficient We notice that definition (116) coincides (not surprisingly) with the constant a defined in Eq. (113), and with the standard effective diffusivity [15,27,63]. Besides, the physical dimensions ofd eff st are those of the reciprocal of time.
Finally, the homogenised equation (45a), with the boundary and initial conditions given in (55) and (57), reduces to

Numerical solution
In this section, we find the numerical solution of the non-local, homogenised problem (117a)-(117c) by means of FE methods. As we previously mentioned, in this context, the effective diffusivity can be found without recurring to solve the cell problem (compare Eqs. (116) and (113)).
To start with, we discretise the time interval [0, t f ] in M subintervals, which we assume of equal amplitude τ . Then, for simplicity of notation, we set c (0) (x, t m ) = u m (x) and we adopt an implicit Euler scheme for Eq. (117a), which is thus approximated as Next, by introducing the space of test functions [73] where H 1 (X h ) is defined analogously to H 1 (Y k ), we put Eq. (118) in weak form. To this end, we multiply Eq. (118) by the test function v ∈ V, and after integrating over X h , we find Next, we discretise the spatial domain X h in N finite elements, and introduce the function basis {ψ i } N i=0 , with ψ i (x j ) = δ i j and i, j = 0, . . . , N . Then, we approximate v(x), the initial condition u 0 (x), and u m (x), for all m, asv where ω m i , with m = 1, . . . , M + 1, are constant coefficients to be determined and t M+1 = t f . Thus, by substituting (121a) and (121c) into (120), and adopting a standard procedure in FE, we find where It is worth to remark that both the stiffness matrix L ji (α) and the nodal force F j (α) depend on the parameter α ∈ ]0, 1[. Then, Eq. (122) can be rewritten as which, by factorising {υ} T , leads to the linear system Details about the explicit form of Eq. (125) are provided in "Appendix B".

Results and discussion
We notice that, in the present framework, to find the numerical solution of the homogenised problem, we only need to knowd eff st as prescribed by Eq. (116). Particularly, by using the values reported in Table 1, we obtain Then, in Fig. 6, we show the distribution of c (0) for different values of α ∈ ]0, 1[ and for different instants of time. According to the plots, and similarly to the previous benchmark problem, the variation of the non-locality parameter influences the way in which diffusion takes place, and in which the stationary state is attained. That is, the progression of the solution towards the stationary states for α = 0.1 is much slower than in the case determined by α = 0.9. In particular, when α approaches 1 from below, the standard diffusion is recovered. We remark that, although we have imposed an initial concentration with very small spatial derivative at the boundary, once time initiates to increase, the tails of the concentration profile tend to raise. This behaviour can be explained by the production of concentration gradients that are needed for the chemical species to diffuse, in this case, towards the centre of the specimen. However, such gradients tend to "turn off" themselves in the course of time since the concentration has to move towards its stationary value. It is worth noticing that the way in which the non-local interactions are introduced influences the diffusion profile of the chemical species (see Fig. 7). Indeed, when considering the existence of non-local interactions at the micro-scale, these are ciphered into the effective coefficientd eff (β), which is parametrised by β, while the effective mass flux has the classical form given by Fick's law. On the other hand, the consideration of long-range interactions at the macro-scale leads, as prescribed by (114), to a non-local effective diffusivity that depends on the spatial points, and thus to a homogenised equation of fractional type for the leading order of concentration. In this case, as shown in Fig. 7, there is a strong memory of the initial concentration, that is, the fractional operators appearing in Eq. (117a) help to preserve the information of the initial distribution of chemical species as time passes. This phenomenon is less pronounced when the non-locality is considered only at the micro-scale, and indeed, for t ≤ 6 h the diffusion near the boundary of the composite is slower. We further notice that, by comparing the curves resulting from the benchmark problem I with the ones from the standard, local framework, the assumption of the non-locality at the micro-scale produces a slower diffusion  Fig. 7 Comparison of the numerical solutions resulting from the benchmark problems I and II with the ones from the local framework. The way in which non-locality is introduced influences the diffusion of the chemical species of the species. Thus, we can conclude that the consideration of non-local effects at the micro-structure impacts the evolution of the concentration at the macro-scale.
Finally, we remark that for α → 1 − and β → 1 − both benchmark tests restore the standard diffusion given by Fick's law and, thus, they become indistinguishable in the limit. For this reason, in Fig. 7, we report only the cases in which there is a strong non-locality, that is when α and β are near zero.

Conclusions
In this work, we study the two-scale, non-local diffusion of a chemical species in a composite medium. This is addressed by prescribing a two-scale constitutive law of fractional type for the mass flux of the chemical species and, with the aid of the asymptotic homogenisation technique, by solving the homogenised version of the original problem. Specifically, we obtain an effective characterisation of the composite, which is subjected to the existence of non-local interactions at both length scales. In fact, we prove that if non-locality is neglected, we recover the classical results of homogenisation theory.
To quantify the effects of the non-locality in our model, we present two simplified benchmark tests in which long-range interactions are considered at the micro-scale or at the macro-scale only. In both cases, we prove that in the limit in which the non-locality parameters β (in the benchmark test I) and α (in the benchmark test II) tend to 1 from below, the fractional cell and homogenised problems lead to the standard ones given in the classical homogenisation literature.
Furthermore, we perform numerical simulations and report the most important steps leading to the finite element discretisation. We remark that the formulation of the FE scheme is very basic and was adapted from the one developed in [44] for one-dimensional problems. Nevertheless, the simplicity of the numerics allows us to discuss some of the specific properties of the algebraic equations resulting from the discretisation process. In particular, in the benchmark problems, the presence of the symmetrised Caputo fractional derivative results in fractional stiffness matrices and fractional nodal forces. More specifically, we prove that the fractional stiffness matrices are symmetric and, although they are dense because of the presence of the fractional derivatives, in the limits α → 1 − and β → 1 − they become the standard stiffness matrices of tridiagonal form. The numerical simulations are in harmony with these theoretical predictions.
It is important to emphasise that the way in which the non-local interactions influence the macroscopic behaviour of the system depends on the scale at which these interactions are introduced. When the non-local interactions are considered at the micro-scale, they emerge also at the macro-scale through the effective fractional diffusivityd eff (β) and by slowing down the diffusion. However, the effect of the non-local interactions is more evident when those are accounted for at the macro-scale, and occurs through a deceleration of the diffusion process that is stronger than in the previous case. Hence, the information enclosed in the initial distribution of the concentration is kept for a "longer" time.
As a first step in our investigations, we conceived a model in one dimension. Clearly, this model can be generalised to higher dimensions. However, there are some issues that must be tackled. One of them is that the non-locality function and the normalisation factors should be conceived in a symmetry-and dimensiondependent way. Moreover, a more detailed numerical study would be required. These issues are part of our current research.

A FE discretisation of the non-local cell problem. Benchmark problem I
Here, we calculate the matrices appearing in Eq. (98) and the non-local effective diffusivity (99). Firstly, we recall that the system (98) reads with For the sake of simplicity, let us consider that the basis functions are Lagrange polynomials of the first order, for i = 1, . . . , N 1 − 1 if k = 1, and i = N 1 + 1, . . . , N 2 − 1 for k = 2, and at the interface y N 1 = y I , we prescribe Then, from the above expressions, the fractional stiffness matrices L k ji (β) are computed in two steps. First, we calculate the symmetrised Caputo fractional derivatives of (129)-(130b) and then, we substitute the results into (128a) and (128b). However, before doing this, we find it convenient to compute the following integrals, From here on, for the sake of a lighter notation, we omit the dependence of H p,q k,l on β. Moreover, we notice that A 1 (y; β) = C 0,N 1 (y; β)/(1 − β), (132a) A 2 (y; β) = C N 1 ,N 2 (y; β)/(1 − β).
A.2 Computation of F k By recalling the definition of F k (β) given in (128b), and using the expressions (132a)-(134b), the components of the nodal fractional force are given by We notice that in Eq. (148), the coefficients ω k i are the solutions of the algebraic equation (98) and represent the nodal concentrations. Therefore, the effective coefficient can be computed after the non-local cell problem has been solved.

B FE discretisation of the non-local homogenised problem. Benchmark problem II
Analogously to what has been done above, we consider the basis functions to be defined by Lagrange polynomials of the first order, i.e.
Then, by using (149b), the mass matrix is given by Now, by recalling the expressions (131a) and (131b), with a slight abuse of notation, we have that where x replaces y, and the parameter α replaces β. In the following discussion, we omit the dependence of H p,q k,l on α. Thus, by using expressions (151a) and (151b), and the symmetrised Caputo derivative of order α of ψ i , i = 1, . . . , N , i.e.
the fractional stiffness matrix L(α) can be computed as follows Moreover, by taking into account that the elements of the fractional nodal force F(α) are given by As previously discussed, also in this case both the fractional stiffness matrix and the fractional nodal force tend to their classical counterparts when α → 1 − .