Multivariate Fractal Functions in Some Complete Function Spaces and Fractional Integral of Continuous Fractal Functions

There has been a considerable evolution of the theory of fractal interpolation function (FIF) over the last three decades. Recently, we introduced a multivariate analogue of a special class of FIFs, which is referred to as α-fractal functions, from the viewpoint of approximation theory. In the current note, we continue our study on multivariate α-fractal functions, but in the context of a few complete function spaces. For a class of fractal functions defined on a hyperrectangle Ω in the Euclidean space Rn, we derive conditions on the defining parameters so that the fractal functions are elements of some standard function spaces such as the Lebesgue spaces Lp(Ω), Sobolev spaces Wm,p(Ω), and Hölder spaces Cm,σ(Ω), which are Banach spaces. As a simple consequence, for some special choices of the parameters, we provide bounds for the Hausdorff dimension of the graph of the corresponding multivariate α-fractal function. We shall also hint at an associated notion of fractal operator that maps each multivariate function in one of these function spaces to its fractal counterpart. The latter part of this note establishes that the Riemann–Liouville fractional integral of a continuous multivariate α-fractal function is a fractal function of similar kind.

R n , we derive conditions on the defining parameters so that the fractal functions are elements of some standard function spaces such as the Lebesgue spaces L p (Ω), Sobolev spaces W m,p (Ω), and Hölder spaces C m,σ (Ω), which are Banach spaces. As a simple consequence, for some special choices of the parameters, we provide bounds for the Hausdorff dimension of the graph of the corresponding multivariate α-fractal function. We shall also hint at an associated notion of fractal operator that maps each multivariate function in one of these function spaces to its fractal counterpart. The latter part of this note establishes that the Riemann-Liouville fractional integral of a continuous multivariate

Preamble
This note aims to offer a modest contribution to the field of fractal interpolation. In particular, we consider a special class of fractal interpolation functions referred to as the α-fractal function, which has played a considerable role in the theory of univariate fractal approximation. Our work in the current note seeks to show that a few results on the construction of univariate α-fractal functions in various function spaces and associated fractal operator (see, for instance, [1]) carry over to higher dimensions.
For a prescribed data set D = {(x 0 , y 0 ), (x 1 , y 1 ), . . . , (x k , y k )} in R 2 with increasing abscissae, there are multitude of methods to construct a continuous function that maps each x i to y i -generally known as interpolation methods-available in the field of classical numerical analysis and approximation theory. Roughly speaking, the fractal interpolation function (FIF for short), as introduced by Barnsley in the original version [2], is a continuous function g : [x 0 , x k ] → R that interpolates D such that the graph of g, denoted by Gr(g), is a self-referential set (fractal set). Here the word fractal or self-referential is used to indicate that Gr(g) is the attractor of an iterated function system [3]. That is, roughly, Gr(g) is a finite union of tranformed copies of itself. For a compendium of the theory of FIF and its applications in interpolation and approximation, the reader is referred to the book and monograph [4,5]; the recent articles [6][7][8] may also be of interest.
In her research works on fractal interpolation, Navascués emphasized a special class of univariate FIFs, named α-fractal functions, (see, for instance, [9,10]) which garnered a significant amount of research attention in fractal approximation theory. It is our opinion that the notion of α-fractal functions assisted the field of fractal interpolation to find connections and consequences in other branches of mathematics such as approximation theory, harmonic analysis, functional analysis and the theory of bases and frames; see, for instance, [11,12]. In the research works reported in [13,14], authors utilized α-fractal functions to demonstrate that FIFs can be applied in various constrained approximation problems.
Several extensions of FIF to higher dimensions, in particular, bivariate FIFs or fractal surfaces, have been studied in the literature; see, for example, [4,[15][16][17][18][19]. Despite that the α-fractal function facilitated the theory of univariate FIF to merge seamlessly with various fields in mathematics, a similar approach to multivariate FIFs was not attempted except for a few research works on bivariate α-fractal functions reported lately in [20][21][22]. The aforementioned works on bivariate α-fractal functions find their origin, perhaps implicitly, in the general framework for the construction of fractal surfaces introduced in [23].
While an increasing amount of literature is being published in the field of univariate FIFs and fractal surfaces, the research in multivariate FIFs are still inadequate, especially in the framework of α-fractal functions. In the context of multivariate FIFs, the ingenious constructions appeared in [24,25], though worth mentioning, do not seem to be suitable for the implementation of the α-fractal function formalism. On the other hand, our acquaintance with the univariate and bivariate α-fractal functions revealed that the development of multivariate analogue of α-fractal function could be highly beneficial for the expansion of multivariate fractal approximation theory. Stimulated by the construction of fractal surface in [23], recently we put forward a satisfactory extension of the Barnsley's theory of univariate FIF to the multivariate case [26].
In this note, we continue to explore the notion of multivariate α-fractal functions. In the first part, we define multivariate α-fractal functions in various function spaces such as the Lebesgue spaces L p , Sobolev spaces W m,p , and Hölder spaces C m,σ . We also hint at some elementary properties of the fractal operator associated with the notion of multivariate α-fractal functions.
Fractal dimension is an important parameter of fractal geometry providing information about the geometric structure of the objects that it deals with. There are different notions of fractal dimension, the two most commonly used being the Hausdorff dimension and box dimension [27]. In particular, the Hausdorff dimension and box dimension of the graphs of fractal interpolation functions have been investigated; see, for instance, [2,4,15,28]. Since the aforementioned fractal dimensions are scale-independent, they may not be useful for describing scale-dependent laws and more complicated phenomena in nature. To this end, a new definition of fractal dimension, referred to as the two-scale dimension, is broached in [29], and it is perhaps more akin to physics than mathematics. However, we are forced to settle for less in the framework of multivariate α-fractal functions considered herein, because the analysis for the fractal dimension of the general nonaffine case is subtle. We shall just mention bounds for the Hausdorff dimension of the graph of the multivariate α-fractal function as an immediate consequence of its Hölder continuity for suitable choice of parameters.
On the other hand, fractional calculus, which broadly deals with derivatives and integrals of fractional order, is rather an old subject. During the last decades, fractional calculus has opened its wings wider to cover several real world applications in science and engineering. Despite being an old subject, fractional calculus continues to be a hot topic of research, resulting in a substantial body of literature; we refer the reader to the informative surveys [30,31]. Some recent developments made in the direction of fractional PDEs and their applications deserve a special mention; see, for instance, [32][33][34][35]. Studies on the interconnection between fractional calculus and fractal geometry have gained significant attention in recent years. For some links between the two-scale problem mentioned previously and fractional calculus, the reader may consult [36]. In the second part of this note, our modest aim is to show that the fractional integral of the multivariate fractal function considered herein is again a fractal function of a similar kind.
Overall, this note discusses how some results in univariate fractal interpolation, to be specific α-fractal functions, fractal operator and fractional calculus of fractal functions, carry over to higher dimensions. We strongly believe that these research findings may assist efforts to find interesting interconnections between multivariate FIFs and the theory of PDEs.

Preparatory Facts
To begin with, we list pertinent definitions and notation for use throughout the remainder of this note.
The set of first n natural numbers shall be denoted by Σ n . For l = (l 1 , . . . , l n ) ∈ (N ∪ {0}) n , called a multi-index, let |l| = ∑ n k=1 l k . Given two multi-indices l = (l 1 , . . . , l n ) and l = (l 1 , . . . , l n ), we say that l ≤ l if l k ≤ l k for all k ∈ Σ n . For l ≤ l, we define: Let n ≥ 2 and Ω = ∏ n k=1 I k be an n-dimensional hyperrectangle, where each I k = [a k , b k ] is a closed and bounded interval in R. For a function g : Ω → R and X 0 = (x 0 1 , . . . , x 0 n ) ∈ Ω, denoted by provided the right-hand side exists.

Function Spaces
The purpose here is to provide a short presentation of various function spaces that are used in this note. We refer to Triebel [37] for more information.
Let C(Ω) denote the Banach space of all real-valued continuous functions defined on Ω, endowed with the sup-norm . ∞ . For a positive integer m, we consider the linear space C m (Ω) defined by C m (Ω) = g ∈ C(Ω) : D l (g) exists and it is continuous for each multi-index l with |l| ≤ m .
For any g ∈ C m (Ω), we define It is well-known that C m (Ω) equipped with . C m (Ω) is a Banach space. Next, we recall the Lebesgue spaces. For 0 < p ≤ ∞, let L p (Ω) = g : Ω → R such that g is measurable and g p < ∞ , where g p is defined as It is a standard result in functional analysis that L p (Ω), . p is a Banach space for 1 ≤ p ≤ ∞. For 0 < p < 1, . p is a quasi-norm, that is, in place of the triangle inequality one has g + h p ≤ 2 1 p g p + h p , and L p (Ω) is a quasi-Banach space.
Let g ∈ L p (Ω). For a multi-index l, a function h : Ω → R is called the l th -weak derivative of g if it satisfies for all infinitely differentiable functions φ with compact support contained in Ω. By a slight abuse of notation, we write l th -weak derivative of g as D l (g) = h.
For 1 ≤ p ≤ ∞ and a non-negative integer m, W m,p (Ω) denotes the Sobolev space with smoothness m and integrability p defined by W m,p (Ω) = g : Ω → R such that D l (g) ∈ L p (Ω) for all multi-index l with |l| ≤ m .
The linear space W m,p (Ω) endowed with the norm is a Banach space. For p = 2, it is a Hilbert space, which shall be denoted by H m (Ω). A function g : Ω → R is Hölder continuous with exponent σ ∈ (0, 1] (or σ-Hölder for all X, Y ∈ Ω and some C g > 0, called a Hölder constant of g. Given a Hölder continuous g : Ω → R with exponent σ, the σ-Hölder semi-norm of g is defined by If m is a positive integer, then the Hölder space C m,σ (Ω) is defined as The space C m,σ (Ω) equipped with the norm is a Banach space. Note that C 0,σ (Ω) coincides with the space of all Hölder continuous functions with exponent σ.

Towards Multivariate FIF
Here we shall equip ourselves with a few rudiments needed for multivariate fractal functions that concern us. As mentioned previously, let I k = [a k , b k ], k = 1, 2, . . . , n be compact intervals in R and Ω = ∏ n k=1 I k be an n-dimensional hyperrectangle. Let n ≥ 2 be an integer and It is worth to note that I k = ∪ N k i k =1 I k,i k and each knot point in the partition of I k is exactly in one of the subintervals I k,i k , i k = 1, 2 . . . , N k mentioned above. We call such a set ∆ as a partition of Ω for an obvious reason.
For convenience, let us introduce the following notation. For a positive integer N, For each i k ∈ Σ N k , let u k,i k : I k → I k,i k be an affine map of the form When the interval I k,i k involved in the definition of the affine map is half-open, the above equation needs to be interpreted in terms of the one-sided limit. For instance, when Using the definition of the map u k,i k , one can verify that (4) Using the above notation, we see that It is easy to observe that the boundary of Ω in the usual metric of R n is

Multivariate α-Fractal Functions in Some Complete Function Spaces
This section targets to construct fractal functions (self-referential functions) in the complete function spaces C m (Ω), L p (Ω), W m,p (Ω) and C m,σ (Ω), which we recalled in the previous section. To this end, let X be any of the function space from the list C m (Ω), L p (Ω), W m,p (Ω), C m,σ (Ω) and f ∈ X, be a fixed function, which we shall refer to as the germ function. Let b ∈ X be a fixed function, called the base function.
For each g ∈ X, and X = ( where The ∏ n k=1 N k -tuple comprised of the real numbers α i 1 ...i n is called the scaling vector and it is denoted by α. We define The main objective in this section is to choose the scale vector α and base function b in (5) so that the Read-Bajraktarević (RB) operator T f is a well-defined map, and, in fact, T f is a contraction map on the function X or a suitable subspace of X. It is worth to emphasize that throughout the current note, a partition ∆ of the hyperrectangle Ω is chosen as mentioned in the previous section.
Suppose that the scaling vector α is so chosen that and b ∈ C m f (Ω). Then the following hold.

1.
The map T f given in (5) is well-defined on C m f (Ω).

2.
In fact, As a consequence, by the Banach fixed point theorem, there exists a unique function f α for all (i 1 , . . . , i n ) ∈ ∏ n k=1 Σ N k ,0 and multi-index l with |l| ≤ m. Moreover, the function f α ∆,b and its derivatives satisfy the self-referential equations given by for all X ∈ ∏ n k=1 I k,i k , (i 1 , . . . , i n ) ∈ ∏ n k=1 Σ N k and multi-index l with |l| ≤ m.
f (Ω) and X = (x 1 , . . . , x r , . . . , x n ) ∈ ∏ n k=1 I k,i k be such that x r ∈ I r,i r ∩ I r,i r +1 for some r ∈ Σ n and i r ∈ intΣ N r ,0 . Note that this is possible only when x r = x r,i r and in that case, by (3), we have for all multi-index l with |l| ≤ m. Thus, That is, D l T f (g) (X) = D l ( f )(X) irrespective of whether x r is considered as a point in I r,i r or as a point in I r,i r +1 . The above observation also yields the following: In particular, T f (g) ∈ C m f (Ω). Next, let g, h ∈ C m f (Ω) and l be a multi-index with |l| ≤ m. Then Taking sum over all |l| ≤ m, we get Since, max (i 1 ,...,i n )∈∏ n k=1 Σ N k Let us consider two base functions as follows: . Figure 1b is the graph of fractal perturbation f α ∆,b with base function b = b 1 and uniform scale vector α, where α i 1 i 2 = 0.7 for all i k ∈ Σ 4 for k = 1, 2. Figure 1c depicts the graph of f α ∆,b with base function b = b 2 and uniform scale vector α as taken previously. Finally, Figure 1d displays the graph of f α ∆,b with base function b = b 2 and uniform scale vector α, where α i 1 i 2 = 0.01. In this case, the parameters satisfy the conditions prescribed in Theorem 1, for m = 1. Thus, Figure 1a,b corroborate the technique demonstrated for the construction of smoothness preserving fractal functions in Theorem 1. Theorem 2. Let f ∈ C m,σ (Ω) and define C m,σ f (Ω) := g ∈ C m,σ (Ω) : g(X) = f (X) ∀ X ∈ ∂Ω . Choose the scale vector satisfying and the base function b ∈ C m,σ f (Ω). Then the RB operator T f defined in (5) is a contraction on C m,σ f (Ω), and its unique fixed point f α ∆,b satisfies the self-referential Equation (6).
Proof. Using Theorem 1 we see that T f (g) ∈ C m f (Ω) for all g ∈ C m,σ f (Ω). We shall show that for all multi-index l with |l| = m, D l (T f (g)) is Hölder continuous with exponent σ. Towards this, let X, Y ∈ ∏ n k=1 I k,i k be two points in the same rectangular mesh. We have where C D l (g) and C D l (b) denote the Hölder constants of D l (g) and D l (b), respectively. If X and Y lie in two distinct but in adjacent meshes, then by taking point on their common boundary and repeating the above steps we get Since the total number of rectangular meshes is ∏ n k=1 N k , for any X, Y ∈ Ω, we have which shows that T f C m,σ f (Ω) ⊆ C m,σ f (Ω). A similar computation reveals also that the map T f : C m,σ f (Ω) → C m,σ f (Ω) is a contraction map, completing the proof.
Corollary 1. Let f : Ω → R be a Hölder continuous function with exponent σ. Assume that a scaling vector α is so chosen that max (i 1 ,...,i n )∈∏ n k=1 Σ N k |α i 1 ...i n | ∏ n k=1 |a k,i k | σ < 1 and the parameter map b is a Hölder continuous function with exponent σ and b(X) = f (X) for all X ∈ ∂Ω. Then the Hausdorff dimension of the graph of the corresponding self-referential function Proof. With the stated hypotheses on α and b, it follows from the previous theorem (with m = 0) that the self-referential counterpart f α ∆,b of f is a Hölder continuous function with exponent σ. Define a map A : where we endow Gr( f α ∆,b ) ⊆ R n+1 and Ω ⊆ R n with the usual Euclidean norm. It is plain to see that A is a surjective Lipschitz map. From fundamental properties of the Hausdorff dimension given in ( [38], Theorem 2, Items (5), (8)) we have For the desired upper bound, let us recall that the Hausdorff dimension of the graph of a Hölder continuous function with Hölder exponent s ∈ (0, 1] whose domain is a compact subset of R n with the Hausdorff dimension equal to d is less than or equal to min{d + 1 − s, d s }( [39], Chapter 10). Therefore, completing the proof.
Theorem 3. Let f ∈ L p (Ω) for 0 < p ≤ ∞. Suppose the scaling vector α is so chosen that Then T f defined in (5) maps L p (Ω) to L p (Ω). Further, T f is a contraction map and hence by the Banach fixed point theorem, there exists a unique f α ∆,b ∈ L p (Ω) such that for X ∈ ∏ n k=1 I k,i k , and (i 1 , . . . , i n ) ∈ ∏ n k=1 Σ N k .
Proof. Using the stated hypotheses, it is easy to verify that the operator T f : L p (Ω) → L p (Ω) is well-defined. What remains is to show that T f is a contraction map. To this end, let g, h ∈ L p (Ω), 1 ≤ p < ∞. We have . . .
proving the claim for the case 1 ≤ p < ∞. The other cases can be dealt similarly.
Next, let us construct self-referential functions associated with a function f ∈ W m,p (Ω). First, let us recall the following result, popularly known as the Leibniz theorem.
If f ∈ W m,p (Ω) and φ is infinitely differentiable on Ω, then φ f ∈ W m,p (Ω) and Theorem 4. Let f ∈ W m,p (Ω) for 1 ≤ p ≤ ∞. Suppose that the base function b ∈ W m,p (Ω) and the scaling vector is chosen so that Then the RB operator T f given in (5) is a contraction map on W m,p (Ω). Consequently, T f has a unique fixed point f α ∆,b .

Proof.
A routine computation yields that the RB operator is well-defined and it maps the space W m,p (Ω) into itself. We shall just show that it is a contraction on W m,p (Ω). To this end, let 1 ≤ p < ∞ and l be a multi-index with |l| ≤ m. We note that . . .
The rest of the theorem follows from the Banach fixed point theorem and the assumption on the scale vector. The case p = ∞ can be worked out similarly.

Fractal Operator on Function Spaces
Let f ∈ X, where X is a fixed function space from the list C m (Ω), L p (Ω), W m,p (Ω), C m,σ (Ω) .
The results established in the previous section provide a self-referential counterpart to each f ∈ X, and consequently provide an operator. That is, for a prescribed set of parameters such as the partition, scale vector and the base function, there exists a fractal This section intends to record a few elementary properties of the operator F α ∆,b : X → X, what we call a multivariate selfreferential operator (fractal operator); see also [1]. We shall provide the details only for X = W m,p (Ω), as the other spaces can be similarly dealt with. For future reference, we introduce the notation Proposition 1. (Perturbation Error) Let f ∈ W m,p (Ω). Suppose that a partition ∆ of the hyperrectangle Ω, base function b ∈ W m,p (Ω), and scale vector α be chosen as in Theorem 4. Then Proof. Let us recall the self-referential equations satisfied by the fractal counterpart f α ∆,b ∈ W m,p (Ω) and its derivatives for all X ∈ ∏ n k=1 I k,i k , (i 1 , . . . , i n ) ∈ ∏ n k=1 Σ N k and multi-index l with |l| ≤ m. Assume that 1 ≤ p < ∞. By simple calculations Therefore, . . .
Similar analysis for p = ∞.
Now, let us take the multivariate base function b : Ω → R used in the construction of the self-referential function f α ∆,b through a suitable operator L : X → X. That is, we take b = L( f ) so that the conditions required for b are satisfied. In this case, the multivariate fractal operator will be denoted by F α ∆,L : X → X. In what follows, we intend to record some elementary properties of the multivariate fractal operator F α ∆,L . The following proposition provides a counterpart to the linearity property of the fractal operator well explored in the setting of univariate α-fractal functions on various function spaces; see, for instance, [11]. The proof follows almost verbatim, and hence omitted. Proposition 2. Let X = W m.p (Ω) and L : X → X be a linear operator. Choose the base function b in the construction of fractal function f α ∆,b via this operator L so that b = L( f ). Then the corresponding fractal operator, which shall be denoted by F α ∆,L : X → X, defined by Let X be a Banach space and A : X → X be a bounded linear operator such that I − A < 1, where I is the identity operator on X. Then, it is well-known that A is bijective and A −1 is bounded; see, for instance, [40]. The following result available in [41] is a generalization of the aforementioned Neumann's lemma.

Lemma 1. ([41]
, Lemma 1) Let A : X → X be a linear operator on a Banach space X such that for some λ 1 and λ 2 ∈ [0, 1). Then A is a topological automorphism (a bounded, invertible map that possesses a bounded inverse). Furthermore, The assertion is now immediate from the previous proposition.
The existence of Schauder bases consisting of appropriate functions for the Sobolev spaces is quite desirable in analysis of PDEs, for instance, for demonstrating the existence of solutions of various non-linear boundary value problems. We have the following result giving a Schauder basis consisting of self-referential functions for the Sobolev space W m,p (Ω). The heart of the matter is an elementary result in the theory of bases, which states that a topological isomorphism preserves Schauder bases; see, for instance, [37].

Corollary 2.
The Banach space W m,p (Ω) has a Schauder basis consisting of multivariate selfreferential functions.
Proof. Let { f m } m∈N be a Schauder basis of W m,p (Ω) whose existence is established and reported, for instance, in [42][43][44]. Choose the scale function α and operator L as in the previous proposition so that the fractal operator F α ∆,L is a topological automorphism. As an isomorphism, in particular, an automorphism, preserves Schauder bases, we conclude that is a Schauder basis consisting of self-referential functions for the Banach space W m,p (Ω).

Fractional Integral of Continuous Multivariate α-Fractal Function
As mentioned in the introductory section, exploration of interconnection between fractional calculus and fractal geometry has always been of interest. Our purpose in this section is limited; we shall observe that the Riemann-Liouville fractional integral of the continuous multivariate α-fractal function is also a fractal function. A similar result regarding univariate FIF can be found in [28].

Definition 1.
[45] Let f be a continuous function on the closed and bounded hyperrectangle Ω in R n . The left-hand-sided mixed Riemann-Liouville fractional integral of f of order γ is defined as . . where a = (a 1 , . . . , a n ) is a fixed point, X = (x 1 , x 2 . . . , x n ) and γ = (γ 1 , . . . , γ n ) with a k ≤ x k , γ k > 0 for each k ∈ Σ n .
Let f ∈ C(Ω). We write f (X) = f (x 1 , x 2 , . . . , x n ). From Theorem 1 it follows that by choosing b ∈ C f (Ω) and scaling vector α such that the fractal counterpart f α ∆,b of f belongs to C(Ω). Furthermore, since f α ∆,b is the fixed point of the RB operator T f : C(Ω) → C(Ω) defined by for all X ∈ ∏ n k=1 I k,i k , (i 1 , . . . , i n ) ∈ ∏ n k=1 Σ N k . Consequently, f α ∆,b satisfies the functional equation so that the self-referential equation for f α ∆,b becomes Since the multivariate fractal function f α ∆,b is continuous, we can talk about its Riemann-Liouville fractional integral. In what follows, we establish that the Riemann-Liouville fractional integral of f α ∆,b is again a fractal function. For the sake of convenience, we shall deal with the uniform scaling factor, that is, α i 1 ...i n = α for all (i 1 , . . . , i n ) ∈ ∏ n k=1 Σ N k . Then, with a slight abuse of notation, the above equation reduces to Theorem 5. Let ∆ be a partition of the hyperrectangle Ω in R n and f ∈ C(Ω). Assume that b : Ω → R is continuous and b(X) = f (X) for all X ∈ ∂Ω, the boundary of Ω. Choose a scaling vector α such that α ∞ < 1. Then I γ a f α ∆,b , the left-hand-sided mixed Riemann-Liouville fractional integral of order γ of the self-referential function f α ∆,b , satisfies the following equation: a n−k x n−k+1 a n−k+1 . . .
x n a n so that I γ a f α ∆,b u i 1 ...i n (X) = E 0 + J 1 . Turning our attention to J 1 , let us change the variable s n using the transformation s n = u n,i n (t n ). We have . .
Applying similar process to the variable s n−1 we get . . u n−1,i n−1 (a n−1 ) a n−1 x n a n u n−1,i n−1 (a n−1 ) x n a n . . , u n,i n (s n ) ds 1 . . . ds n =: E 1 + J 2 .
Proceeding in the same fashion, at the n th step we get . .  . .

Conclusions
The α-fractal formalism of fractal interpolation function is proved to be beneficial in expanding the applications of univariate fractal approximation theory. Through the construction of multivariate α-fractal functions on a few complete function spaces which are ubiquitous in the theory of partial differential equations and harmonic analysis, the present work intends to be a step forward in the theory of multivariate fractal approximation. The construction of self-referential analogue for each germ function in a complete function space under consideration leads naturally to an operator, referred to as the multivariate fractal operator. We have studied a few elementary properties of the fractal operator. The multivariate fractal operator introduced and studied here enabled us, in particular, to construct Schauder bases consisting of self-referential functions for the function spaces. Further, taking a slight detour from the main theme, it is shown that the Riemann-Liouville fractional integral of a self-referential counterpart of the given multivariate germ function will also be a self-referential function under some suitable conditions.