1 Introduction
Nonlocal models are becoming a popular alternative to partial differential equations (PDEs) when the latter fail to capture effects such as multiscale and anomalous behavior. In fact, several scientific and engineering applications exhibit hierarchical features that cannot be described by classical models. As an example we mention applications in continuum mechanics
[ha2011characteristics, littlewood2010simulation, ouchi2015fully][Bates1999, Burkovska2020CH, Chen_nonlocalmodels], corrosion [Rokkam2019], turbulence [Bakunin2008, Pang2020, Schekochihin2008], and geoscience [Benson2000, Meerschaert2006, Schumer2001, Schumer2003, WeissBloemenEtAl2020_FractionalOperatorsAppliedToGeophysicalElectromagnetics]. In this work we consider nonlocal operators of fractional type, i.e. integral operators characterized by singular and nonintegrable kernels, that correspond to differential operators of fractional orders, as opposed to the integer order in the PDE case.Note that fractional differential operators are almost as old as their integer counterpart [Gorenflo1997]; however, their usability has increased during the last decades thanks to progress in computational capabilities and to a better understanding of their descriptive power. The simplest form of a fractional operator is the fractional Laplacian; in its integral form, its action on a scalar function is defined as [Gorenflo1997]
(1) 
where is the fractional order, the spatial dimension, a constant and where p.v. indicates the principal value. It follows from (1) that in fractional modeling the state of a system at a point depends on the value of the state at any other point in the space; in other words, fractional models are nonlocal. Specifically, fractional operators are special instances of more general nonlocal operators [DElia2020, DElia2013, Du2012, Pang2020] of the following form [DElia2017]:
(2) 
Here, interactions are limited to a Euclidean ball of radius , often referred to as horizon or interaction radius. The kernel is a modeling choice and determines regularity properties of the solution. Note that for and for the fractionaltype kernel the nonlocal operator in (2) is equivalent to the fractional Laplacian in (1). Also, it has been shown in [DElia2013] that for that choice of and solutions corresponding to the nonlocal operator (2) converge to the ones corresponding to the fractional operator (1) as (see [DElia2020] for more convergence results and for a detailed classification of these operators and relationships between them).
In recent years, with the purpose of increasing the descriptive power of fractional operators, new models characterized by a variable fractional order have been introduced for both space and timefractional differential operators [Razminia2012, AntilRautenberg2019_SobolevSpacesWithNon, Zheng2020] and several discretization methods have been designed [Chen2015, Zeng2015, Zheng2020optimalorder, Zhuang2009, SchneiderReichmannEtAl2010_WaveletSolutionVariableOrderPseudodifferentialEquations]. However, the analysis of variableorder models is still in its infancy, with [Felsinger2015] and [Schilling2015] being perhaps the only relevant works that address theoretical questions such as wellposedness for spacefractional differential operators with variable order
. The improved descriptive power of variableorder operators has been demonstrated in some works on parameter estimation
[Pang2020, Pang2019fPINNs, Pang2017discovery, WANG]; here, via machinelearning algorithms or more standard optimization techniques, the authors estimate the variable fractional power for operators of the form
(3) 
where is either infinite or finite.
Models such as the one in (3) are convenient to describe anomalous diffusion in case of heterogeneous materials or media, where different regions may be characterized by different diffusion rates, i.e. different values of , see Figure 1 for twodimensional illustrations where is a piecewise constant function defined over the domain. In other words, these models can describe physical interfaces by simply tuning the function . Modeling nonlocal interfaces is a nontrivial task due to the unknown nature of the nonlocal interactions across the interfaces. Only a few works in the literature have addressed this problem; among them, we mention [Alali2015] for nonlocal models in mechanics and [Capodaglio2020] for nonlocal diffusion models with integrable kernels.
In this work, we tackle the nonlocal interface problem for fractionaltype operators using a generalization of the nonlocal operator in (3). Specifically, we add variability to the fractional order and the kernel function itself. Our main contributions are:

The introduction of a new variableorder fractional operator characterized by and . This choice enables modeling of a much broader set of interface behaviors and, for symmetric , symmetrizes the kernel, making analysis and implementation a much easier task. In fact, note that in (3) the kernel is no longer symmetric, requiring more sophisticated quadrature rules for the integration (regardless of the discretization technique of choice).

The analysis of wellposedness of the Poisson problem associated with the new operator in the general nonsymmetric case.

The design of an efficient finiteelement matrix assembly technique for the practical case of piecewise constant definition of the fractional order that often appears in, e.g., subsurface applications.
Outline of the paper
In Section 2 we first recall some stateoftheart results; then, we introduce the new variableorder fractional operator and the corresponding strong and weak forms of the nonsymmetric diffusion problem. We also describe properties of the associated energy norm and space. In Section 3 we report the main result of the paper that proves via Fredholm alternative the wellposedness of the diffusion problem. In Section 4 we introduce a finite element discretization and propose efficient quadrature rules for the integration of the nonsymmetric kernel. In Section 5 we report several one and twodimensional numerical tests that illustrate both the improved variability of our model and the accuracy of the numerical scheme. Finally, in Section 6 we summarize our findings.
2 A new variableorder fractional model
In this section we first introduce the notation that will be used throughout the paper and recall stateoftheart results on nonsymmetric kernels and fractional kernels with variable order. Then, we introduce a new variableorder fractional kernel and discuss the properties of the associated norm and energy space.
2.1 Preliminaries
The purpose of this section is twofold. First, we recall the formulation of a symmetric nonlocal diffusion problem with finite interaction radius following the theory presented in [Du2012, Du2013]; then, we recall the theory presented in [Felsinger2015] on variableorder fractional operators (with infinite interaction radius). Our new model, presented in the next subsection, is a combination of the two.
Let be an open bounded domain. We define its interaction domain as the set of points outside that interact with points inside , i.e.
(4) 
where , referred to as interaction radius, determines the extent of the nonlocal interactions. Also, let be a nonnegative, symmetric, kernel function, not necessarily integrable. The diffusion operator introduced in [Du2013] is defined as
(5) 
From now on we remove the explicit reference to the principal value when the operator is used within an equation. The strong form of a truncated, symmetric, diffusion problem is then given by: for , and , find such that
(6)  
The wellposedness of problem (6) was studied in [Du2012]. Extensions of this model to the nonsymmetric case have been proposed in [Duadvection], for integrable kernels, and in [DElia2017] for nonintegrable kernels by using the operator (2) introduced in the previous section. This operator is often referred to as nonsymmetric nonlocal diffusion operator or, more appropriately, nonlocal convectiondiffusion operator as it introduces a nonlocal convection term. Since the nonsymmetry of our problem is not related to convection, but only to the presence of an interface that affects the way a quantity diffuses, in this work, we do not employ (2); instead, we use the operator defined in (5) with a nonsymmetric kernel .
Among the stateoftheart formulations of nonlocal diffusion problems with nonsymmetric kernels, we mention the one introduced in [Felsinger2015] for integrable and nonintegrable kernels with infinite support. Here the authors consider the nonlocal operator (5) where the kernel is a nonnegative measurable function, possibly nonsymmetric. Among the kernels studied in [Felsinger2015], we recall the following nonsymmetric variableorder fractional kernel, whose analysis is relevant for the theory presented in this paper:
(7) 
where the nonsymmetry with respect to is due to the terms and . The latter is a normalization functions with the same role as in the standard fractional Laplacian in (1) [Schilling2015]. Note that the choice of corresponds to the variableorder version of the standard fractional Laplacian operator introduced in (3). Furthermore, when and are constant functions, such operator is a multiple of the standard fractional Laplacian. The wellposedness of problem (6) for nonsymmetric kernels was analyzed in [Felsinger2015] in a very general setting and in [Schilling2015] for kernels of the form (7).
2.2 The variableorder model
We the purpose of increasing the descriptive power of the operator , we consider a fractionaltype kernel with added variability. We define the new variableorder nonsymmetric fractionaltype kernel as follows:
(8) 
with constants , and such that
(9)  
(10) 
Thus, the resulting nonlocal operator is given by the following expression
(11) 
Remark 1
As opposed to the kernel model presented in [Schilling2015], the function does not have a normalizing purpose but depends on the type of phenomenon that we are targeting. As an example, could describe a material property, such as diffusivity or permeability in the subsurface; in this context, the dependence on and describes changes in material properties across material interfaces.
Remark 2
For a specific choice of , the operator defined in (11) corresponds to a variableorder, tempered fractional Laplacian. Specifically, we let
(12) 
where is the socalled tempering parameter, whose effect is to fasten the decay of the fractional kernel. Since the overall effect of tempering is similar to the one of truncation[Olson2020], in our numerical tests we do not study the tempered case and we limit ourselves to piecewise constant definitions of for .
Remark 3
When the kernel is symmetric. This case is much easier to deal with both in terms of analysis and computation. In fact, in the symmetric case, the theory introduced in [Du2012] can be readily applied and guarantees wellposedness of problem (6). Furthermore, numerical integration of a symmetric function poses minor challenges compared with the nonsymmetric case, which requires the employment of more sophisticated quadrature rule. As we show in our numerical tests, a symmetric still guarantees improved model variability across interfaces, compared to the singlevariable model in (7).
We define the spaces
where the seminorm is defined as
Here, is the symmetric part of the kernel. We recall that the symmetric and antisymmetric parts of the kernel are defined as follows:
(13)  
(14) 
The following bounds hold trivially:
(15)  
(16) 
In what follows, we will also frequently make use of the following lower bound on the horizon:
(17) 
The following Lemma shows that the seminorm above satisfies a Poincaré inequality and, hence, equips with a norm.
Lemma 1
For all , the following Poincaré inequality holds:
Proof
The weak form. Let for simplicity; then, we can write the weak form of problem (6) as
(18) 
We denote the bilinear form in (18) by
(19) 
While this form is formally equivalent to the one introduced in [Felsinger2015], the kernel (8) does not belong to class of kernels studied in that paper. In the next section we show how to extend the wellposedness theory developed in [Felsinger2015] and [Schilling2015] to our new operator.
For discretization purposes, it is convenient to further rewrite as
(20) 
This form allows for simpler numerical integration and is used in our numerical implementation.
3 Wellposedness analysis
In this section we prove the wellposedness of problem (18) by adapting the theory developed in the works of Felsinger, Kassmann and Voigt [Felsinger2015] and of Schilling and Wang [Schilling2015]. We first report the building blocks of our main wellposedness result.
The following proposition establishes a relation between the symmetric and antisymmetric part of the kernel and is a generalization of Proposition 3.1 in the paper by Schilling and Wang [Schilling2015].
Proposition 1
Assume that
(21) 
satisfies
(22) 
Moreover, assume that
is Hölder continuous for arbitrarily small. Then there exists a constant such that
(23) 
Proof
First, by definition of symmetric and antisymmetric parts of the kernel in (13), we have that . This and assumptions (9), (10) and inequality (16) imply that there exists a positive constant such that
Moreover, due to the definition of ,
so that we obtain
(24) 
Next, for , we write:
Then, by (10) and the Hölder continuity of we have that there exists a positive constant such that
(25) 
Here, we have used (15). On the other hand, since for
then, for , we have that
Hence, by (10) and assumption (22), there exists a positive constant such that
(26) 
Here, we have again used (15). Hence, from (24), (25) and (26), we obtain:
and the result follows.
The next Gårding inequality follows directly from [Felsinger2015, Lemma 3.1]; we report its proof for completeness.
Lemma 2 (Gårding inequality)
There exist constants such that
(27) 
Proof
By using Young’s inequality for , we obtain
Hence, for small enough, we obtain
The Poincaré inequality and Proposition 1 yield the continuity of the bilinear form , as illustrated in the following lemma.
Lemma 3
For all
i.e., the bilinear form is continuous on .
Proof
Before proceeding with the wellposedness analysis, we prove a result for a perturbation of the weak problem (18). First, let
Then, there exists some such that
and hence for all ,
(28) 
Hence, is continuously embedded in . Since is compactly embedded in , and is compactly embedded in , forms a Gelfand triple, and the embedding of into is compact. This fact, and the Gårding inequality (27), directly imply the wellposedness of a perturbation of problem (18) and a bound on its solution, as shown in the following lemma. In what follows, we denote the duality pairing between and its dual in the usual manner, i.e. .
Lemma 4
Given a positive constant and a function , the problem
is wellposed and its solution is such that , for a positive constant .
The proof is simply based on application of (27) to the perturbed bilinear form, and is hence omitted. Note that (28) also allows us to apply [Felsinger2015, Theorem 4.1], which states that the bilinear form satisfies a weak maximum principle, which we report for completeness.
An immediate consequence of Lemma 4 is the fact that the operator mapping , is a compact operator. This fact allows us to use the Fredholm alternative theorem to prove the wellposedness of problem (18).
Theorem 1
Given a function , the problem
has a unique solution .
Proof
Consider the operator . We have the following identity:
By the Fredholm alternative, either
is an eigenvalue of
, or is a bounded linear operator. We show that the first option leads to a contradiction.Let and assume that is an eigenvalue of . This implies that for all and, by (29), that . Hence, is not an eigenvalue of . This ends the proof since we can conclude that , and subsequently , are bounded linear operators. In other words, the variational problem associated with the operator is wellposed.
4 Finite element discretization and efficient matrix assembly
In this section we briefly introduce the finite element (FE) discretization of problem (18) and provide details regarding numerical integration.
We choose a family of finitedimensional subspaces , parameterized by the mesh size . We define the Galerkin approximation as the solution of the following problem: find such that
(30) 
As in [Du2012], we consider finite element approximations for the case that both and are polyhedral domains. We partition into finite elements and denote by the diameter of the largest element in the partition. We assume that the partition is shaperegular and quasiuniform [Brenner2008] as the grid size and we choose the subspace to consist of piecewise polynomials of degree no more than defined with respect to the partition.
We perform numerical integration on pairs of elements, by rewriting the bilinear form (20) as the sum of integrals over the partition, i.e.
The design of quadrature rules for the integrals above depends on the regularity properties of the functions and . If and are constant on a pair of elements , we can use the stateoftheart quadrature rules described in [AinsworthGlusa2018, AinsworthGlusa2017] designed for the constantorder fractional Laplacian, i.e. for the case . However, when this is not the case, the use of more sophisticated and often more expensive quadrature rules, such as adaptive quadrature rules [PiessensDonckerKapengaEtAl2012_Quadpack], becomes mandatory to prevent the integration error to be dominant. Conditions for wellposedness derived in Proposition 1 are such that any nonsymmetric kernel falls into the latter category, hence requiring higher computational effort. In addition, note that, as for the constantorder fractional Laplacian, when , the integrand function features a singularity at . In this case, the domain of integration must be partitioned into subdomains in order to avoid the singularity.
In the specific case of piecewise constant fractional order, we speed up the assembly of the discrete operator and subsequent matrixvector multiplications by generalizing the panelclustering approach of
[AinsworthGlusa2018, AinsworthGlusa2017]developed for the constantorder fractional Laplacian. Informally speaking, in the standard version of this approach, the unknowns are recursively grouped into nodes of a tree structure, with the root node containing all unknowns. Pairs of clusters correspond to subblocks of the discrete operator. If certain conditions on their location in physical space are satisfied, subblocks are approximated using Chebyshev interpolation, otherwise they are assembled using quadrature
^{a}^{a}aSee the works by Ainsworth and Glusa[AinsworthGlusa2018, AinsworthGlusa2017] for a detailed description of this approach for constant parameters and infinite horizon.. We generalize this algorithm to finite horizon and piecewise constant coefficients by restricting the approximation of matrix subblocks to cluster pairs that are strictly contained within the horizon of each other and whose interaction reduces to a constant fractional kernel.5 Numerical tests
5.1 Comparison of different types of kernels
In order to highlight the different behaviors of available fractional kernels, we restrict ourselves to a onedimensional problem where model parameters are constant over two subdomains. This configuration can be considered a onedimensional counterpart of the twodimensional configurations displayed in Figure 1. Specifically, we let , , and consider kernels given by (8) that are piecewise constant with respect to the partition . In order for the assumptions of Proposition 1 to be satisfied, such kernels need to be symmetric, i.e. and need to take the form
Note that this excludes the interesting cases, that can be found in the literature, of piecewise constant or that only depend on . Hence, we also consider functions of the following form and compare them with the proposed model.
Here, determines the size of the transition region between the two states and , and the steepness of the transition. Since is Lipschitz continuous, we can use for both and and satisfy the conditions of Proposition 1.
We numerically solve the diffusion problem (6) for homogeneous Dirichlet data and forcing term . In all tests conducted in this section we use continuous piecewise linear FE on a uniform mesh of size .
With the purpose of highlighting how a variable fractional order affects the qualitative behavior of the solution across the interface, we first compare the family of symmetric kernels reported in Table 1 and display the corresponding solutions in Figure 2. We observe that the solutions for constant seem to be differentiable, whereas the cases with piecewise defined are only continuous at .
identifier 

Comments
There are no comments yet.