The main objective of this paper is to investigate a discrete version of the eigenvalue problem on the Koch snowflake domain, which we denote by . To this end, we follow  and introduce a Dirichlet form (with a suitable domain) on ,
where is the usual Lebesgue measure on and denotes the Kusuoka-Kigami Dirichlet form on the Koch snowflake boundary . The novelty of this approach by comparison with past work [19, 23, 42, 37, 41, 53] is that both the Euclidean interior and the fractal boundary carry non-trivial Dirichlet forms. We approximate the Dirichlet form by a sequence of discrete energies (Theorem 4.1). This is done by inductively constructing trianguations of , the edges of which are then treated as a sequence of finite planar graphs equipped with a discrete energy and measure , and hence an inner product . We then define a discrete Laplacian such that
Note that takes into account both the interior and fractal boundary . Our approach is related to several recent works on diffusion problems involving fractal membranes [33, 13, 36]. Moreover, it can be viewed as a generalization of work of Lapidus et. al.  on the eigenstructure of the Dirichlet Laplacian on . Indeed, our numerical results on the spectra of are closely related to those for a discretization of the eigenvalue problem with Dirichlet boundary conditions due to a localization phenomenon that will be explored in Section 5.
One of the main goals of this paper is to investigate the impact of the fractal boundary on the eigenmodes of . This problem is of general interest, particularly in physics where questions like “How do ocean waves depend on the topography of the coastlines?” and “How do trees and wind interact?” have already been studied. For these and more examples the reader is referred to [49, 51] and references therein. We note in particular that Sapoval  conducted an experimental investigation of acoustic vibration modes of with soap bubbles placed on fractal drums, and made the striking observation that the fractal boundary causes some low-frequency wave modes to localize. Placing a soap bubble on the fractal boundary of a drum is mathematically equivalent to imposing Dirichlet boundary conditions; our situation is somewhat different, in that our model allows for non-trivial boundary energy, but nevertheless we find significant localization effects, which may be summarized as follows:
The eigenvalue counting function has two regimes with different scaling. The threshold for the change in regimes is approximately the largest eigenvalue of the Dirichlet Laplacian on .
Eigenfunctions with eigenvalues of below the threshold are not localized.
Eigenfunctions with eigenvalues of above the threshold begin to show localization near . As the eigenvalue increases, so does this localization.
Eigenfunctions corresponding to eigenvalues of that are significantly above the threshold are strongly localized on .
A class of boundary-localized high-frequency eigenmodes known as whispering-gallery modes are well-understood for convex domains, though an analytic explanation for their appearance in more general domains is not yet known . However, our boundary modes do not seem to be whispering gallery modes, or any other known class of localized high-frequency modes. Rather, they appear to arise because the boundary part of the Dirichlet form and the interior part have different scalings and hence interact only weakly. An elementary way to see that such modes should appear in our model is given using a variant of the landscape map of Filoche and Mayboroda  in Section 6.
The results given here are part of a long term study that aims to provide robust computational tools to address, in a fractal setting, a number of linear and nonlinear problems arising from physics [22, 10, 5, 7, 16, 6]
. The physics of magnetic fields, and of vector equations more generally, are particularly challenging on fractal spaces. Moreover, a discretization of the type considered here is expected to be essential in studying quantum walks[3, 47, 9, 4, 29]. On the mathematical side, our work is related to [1, 2, 24, 12, 15, 11, 48, 26, 27, 52, 8, 20].
The paper is organized as follows. Section 2, follows the treatment in  to introduce a Dirichlet form on the Koch snowflake. In Section 3 we give the Dirichlet form on the snowflake domain and discuss some of its properties, such as the Hölder continuity and uniform approximation of eigenfunctions. In Section 4 we construct a triangular grid to approximate the Koch snowflake domain and introduce the discrete Laplacian . Section 5 is concerned with our algorithm and numerical results, including the existence of boundary-localized eigenfunctions and their effect on the eigenvalue counting function. Finally, in Section 6 we show that the localized eigenfunctions can be predicted numerically using a variant of an argument from .
2. Dirichlet form on the Koch snowflake
The Koch snowflake and the associated snowflake domain are well-known. Here we introduce some notation and foundational results for our analysis, following . Let be the iterated function system defined on by
The Koch curve is the unique nonempty compact subset of such that . It can be approximated by a sequence of finite graphs for which the following notation is convenient.
Let . We define and call an infinite word. Similarly, a finite word of length is ; we write for its length.
We write , where and introduce finite graphs approximating the Koch curve as follows.
Let and for a finite word . Then define
We consider the points of as vertices of a graph in which adjacency, denoted by , means that there is a word of length such that .
On each of these graphs we define a graph energy for by
Following the general treatment in , to which we refer for all omitted details, we see that is nondecreasing, so is well-defined; setting its domain to be one obtains a resistance form. This non-negative definite, symmetric quadratic form extends to , where is the space of continuous functions on . There is a resistance metric on defined from and with the property that points with have comparable to , and thus is bi-Hölder to the Euclidean metric on with
. There is a resistance estimate for
so in particular these functions are -Hölder in the Euclidean metric, see also [21, Corollary 4.8].
Continuing our use of results from 
, we see that if equip the Koch curve with the standard Bernoulli probability measure, i.e. the self-similar measure with weights and for then gives a strongly local regular Dirichlet form on .
A particular collection of functions in of will be useful in what follows. A function on is called harmonic if it minimizes the graph energies for all . It is called piecewise harmonic at scale if it is harmonic on the complement of , or equivalently if is harmonic for each word of length . Piecewise harmonic functions are uniform-norm dense in and dense in with respect to the norm .
As in , we transfer the above definitions for the Koch curve to the boundary of the snowflake domain in a the obvious manner. Write the boundary as a union of three congruent copies of . Specifically, let where is the Euclidean translation and rotation such that and .
The boundary energy and its domain are defined by:
We then let
so that is a strongly local Dirichlet form on . Since we have .
3. Dirichlet form on the snowflake domain
We wish to consider a Dirichlet form on that incorporates our form as well as the classical Dirichlet energy on . The latter is, of course, simply , where is Lebesgue measure. The domain is the Sobolev space . The fact that a nice Dirichlet form of this type exists depends on results of Wallin  and Lancia , which we briefly summarize.
The trace of to is well defined and can be identified with a Besov space , the details of which will not be needed here . Moreover, the kernel of the trace map is , the -closure of functions with compact support in . The domain can be identified with a closed subspace of , and there is a bounded linear extension operator from to , see . Writing for the trace of we see that the following defines a Hilbert space and inner product, see [36, Proposition 3.2].
One consequence of the preceding is that if we let where is Lebesgue measure on and is from (2.3) then for any the quadratic form
with domain in is a Dirichlet form.
Another consequence that will be significant in the next section is that may be written as the sum of and the subspace of obtained by the extension operator from . In this context it is useful to describe an explicit extension operator given in  as the “second proof” of their Theorem 6.1, and to extract some further features of the extension and their consequences. We refer to  for a more detailed exposition of the construction, including diagrams of the hexagons and the triangulation.
The extension operator for a function on is defined as follows. There is an induction over that defines an exhaustion of by regular hexagons of sidelength which meet only on edges. As the induction proceeds, each new hexagon is subdivided into
equilateral triangles meeting at the center and a piecewise linear function is defined on each triangle by linear interpolation of the values at the vertices. The value at the center of the hexagon is defined to be the average of its vertices, and the induction is such that the hexagon vertices either lie on edges of the triangles from the previous stage of the construction, which provides the values at these points, or lie on. In the latter case the values of are used at these vertices, and we note that those vertices of a hexagon with sidelength that come from are points from the vertex set , see Definition 2.2. The resulting piecewise linear function is the extension of to , which we call . It is a special case of the result in [25, Theorem 6.1] that if is a piecewise harmonic function then is Lipschitz in the Euclidean metric on .
With regard to the following theorem, we note that the existence of a bounded linear extension in this setting is known , and we could develop the corollaries from this, but it will be useful for us to know the result for this specific extension.
The extension operator described above is a bounded linear extension operator from with norm to with its usual norm, and satisfies the seminorm bound
We need the following elementary lemma.
If is a linear function on an equilateral triangle then over the triangle is a constant multiple of the sum of the squared differences between vertices, independent of the size of the triangle. Thus over a hexagon in the extension construction is bounded by a constant multiple of the squared vertex differences summed over all pairs of vertices of the hexagon.
The first statement follows from a direct computation for a triangle of sidelength that is left to the reader, combined with the observation that when rescaling length the scaling of cancels with that for the area of the triangle. For the second, write the integral as a sum over the triangles, use the first statement on each triangle and apply the Cauchy-Schwarz inequality. ∎
Proof of Theorem 3.1.
We begin by observing that all operations in the definition of the extension are linear in , and consequently so is the extension operator. In the following argument is a constant that may change value from step to step, even within an inequality.
Observe that in the construction, the vertex values of a hexagon of sidelength come either from values at at points of , or from linear interpolation of values of at vertices from on the side of a hexagon of the previous scale, or from linear interpolation between a value at a point in and a value obtained as an interpolant between values at vertices from . The observation that drives the Lipschitz bound in  is that in any of these cases we can bound the pairwise difference between values at vertices by where the sum ranges over neighbor pairs in that are within a bounded distance from the hexagon. It follows that the number of terms in the sum has a uniform bound, and thus the squared pairwise differences are bounded by . In view of Lemma 3.2, the left side of (3.2) is bounded by a constant multiple of the sum of squares over edges in the triangulation. The preceding argument gives a bound for the sum over the edges that are in a hexagon of size , and summing over the hexagons yields (3.2).
For the norm bound we must control . It is immediate from the construction that . However the latter can be bounded in a standard manner from the resistance bound (2.2), because it implies that on an interval of size controlled by . Direct computation then gives . Since we obtain . Together with the seminorm bound (3.2) this proves the extension operator is bounded. ∎
The domain of the Dirichlet form may be written as a sum of compact subspaces of
where is the image of under the extension map.
The decomposition as a sum of closed subspaces follows from Wallin’s result  that the kernel of the trace map is , and from Lancia’s bounded linear extension  (or from Theorem 3.1 above). Compactness of is from the classical Rellich–Kondrachov theorem. Compactness of is an immediate consequence of the density of the finite dimensional space of harmonic functions in and the boundedness of the extension map in Theorem 3.1 because it exhibits as the completion of a sequence of finite dimensional spaces in . ∎
The following consequence is standard.
The non-negative self-adjoint Laplacian associated to the Dirichlet form by has compact resolvent and thus its spectrum is a sequence of non-negative eigenvalues accumulating only at .
We also note that the extension construction gives us an explicit Hölder estimate.
Functions in are -Hölder in the Euclidean metric on .
The resistance estimate (2.2) says that a function is -Hölder in the resistance metric. A pair of neighboring points in are separated by resistance distance , so for such points. Moreover these points are separated by Euclidean distance , so is -Hölder in the Euclidean metric.
The values of on are those used as boundary data on hexagons of sidelength in the extension construction, and they obviously contribute a term with gradient bounded by . Since the extension on such a hexagon involves terms from the extensions to hexagons of scale for , we may sum these to see that on this hexagon. In particular, . This is true for hexagons of all scales, hence for all points in , and was already known for points on , so the proof is complete. ∎
There are discontinuous functions in because it contains . However, eigenfunctions of , which exist by Corollary 3.4, can be shown to be Hölder continuous.
Sketch of the proof.
Suppose is an eigenfunction with eigenvalue , so for all , so by Corollary 3.3 it is in particular true for all . We will write for the usual seminorm in .
Let be the harmonic function on with boundary data , where we recall that the latter denotes the trace. Since is the extension of that minimizes and we know is an extension for which , we have and it follows immediately that , and in fact . Then for we compute as follows, using because is harmonic and that various boundary terms vanish because on and , the latter because is the kernel of the trace map.
This shows that satisfies , where is the Dirichlet Laplacian. This reduces the problem to determining the regularity of the harmonic function and , where is the Dirichlet Green’s operator on .
It is proved in  that the kernel of the Green’s operator is in a Sobolev space with , and it follows that is Hölder continuous. Hence so is . The Hölder continuity depends on the regularity of the boundary and can, in principle, be estimated using the quantity on page 337 of .
The proof of the Hölder continuity of discussed above admits a generalization which does not involve the Riemann mapping and is applicable in any dimension. We refer to [30, Chapter 1 Section 2] for the following reasoning; all references are to the definitions and results stated there. Our domain is class as in Definition 1.1.20, so a Lipschitz function on has harmonic extension in a Hölder class by Lemma 1.2.4. Taking a sequence of Lipschitz approximations to out boundary data and examining the constants in the proof of Lemma 1.2.4 we discover that these blow-up in manner reflecting the Hölder continuity of , and in particular that is Hölder continuous with exponent the worse of the Hölder exponent of the boundary data and the exponent for a Lipschitz function. A particular case is Remark 1.2.5. We are indebted to Tatiana Toro for pointing out this reference and sketching the argument.
A final consequence of Theorem 3.1 which will be significant later is that when high frequency oscillations on are extended to they do not penetrate the domain very far and the energy of the extended function is very close to the energy of
If the best piecewise harmonic approximations to at scale is the zero function then the extension is supported within distance of the boundary and satisfies a seminorm bound of the form
We saw in the construction that the values from were the only interpolation data for hexagons of sidelength , . If is zero at these points then on these hexagons, so its support is within distance at most from . Moreover, the terms corresponding to , in the middle term of the seminorm bound (3.2) are zero, in which case comparing this to the definition (2.1) of the boundary energy we see that the energy scaling factor provides an additional factor of . ∎
4. Inductive mesh construction and discrete energy forms
We inductively construct the usual approximations to the closed snowflake domain , along with a triangulating mesh, in the following manner. The scale approximation consists of a collection of equilateral triangles with side length , the mesh is their edges, and those edges which belong to only one triangle are called boundary edges. When there is exactly one equilateral triangle. The scale approximation is obtained as shown in Figure 1: each scale triangle is subdivided into nine triangles as on the left, and triangles are appended to the centers of boundary edges as on the right. (The figure shows the case ; at future steps at most two edges of a triangle can be boundary edges.) The approximation is (f) in Figure 2.
The mesh is treated as a graph with vertices . Note that the the vertices in the graphs in Section 2 are a subset of and indeed the graphs constructed there are simply the boundary subgraphs of the . Evidently, two vertices are connected by an edge of the mesh if and only if , where is the Euclidean norm, in which case we write . The vertices contained in boundary edges of are called boundary vertices, all other vertices are interior vertices.
Let and define the graph energy of by
where is the conductance between points , and is given by
Note that for boundary vertices both (4.2) and (4.3) agree with those in Section 2 and therefore the restriction of (4.1) to edges in the boundary of coincides with the energy defined in (2.1). Moreover, the terms corresponding to edges in have weight , and it follows from Lemma 3.2 that for a function which is piecewise linear on equilateral triangles the sum of squared differences of the function over edges in the triangulation is a constant multiple of the Dirichlet energy . After computing this constant and verifying that our triangulation by is obtained from the triangulation in the extension of Theorem 3.1 simply by dividing all triangles of sidelength larger than into subtriangles of sidelength the following consequence is immediate.
If is the extension of by the extension procedure of Theorem 3.1 then .
should be thought of as the graph energy of the restriction of to .
We now introduce a measure on by
This equips with an inner product
and we can then define a discrete Laplacian so that as follows.
5. Numerical Results
In order to investigate the spectral properties of the Dirichlet form we constructed the discrete Laplacian matrix from Definition 4.3 and solved
numerically. We compared the resulting eigenvalues and eigenvectors with those of the Dirichlet Laplacian, which we denote. In this section we almost exclusively present data for , and denote the eigenvector and eigenvalue of by and , where for all and eigenvalues are repeated according to their multiplicity. We label the eigenvectors and in the same manner.
In order to justify the above numerical approach to investigating the spectrum of we should like to prove a result of the following kind. We believe this ought to be possible by methods analogous to those in [35, 14, 34].
Our discrete and finite element approximations to eigenfunctions of converge uniformly, and hence the spectrum of converges to that of .
5.1. Algorithm and Implementation
We constructed by an iterative procedure written in the programming language Python and then computed and graphed the eigenvalues and eigenfunctions of using Mathematica. We also computed the matrix of the Dirichlet Laplacian so that comparison of and could be used to identify effects of the boundary Dirichlet form.
The construction of involved constructing the mesh graphs . The vertices of , which is shown in Figure 2(a), were hard-coded. Construction of the vertices of was achieved by the steps (b)–(f) of Figure 2: specifically, was scaled by , translated to cover the region in the first quadrant and reflected to the other quadrants. Duplicate vertices were deleted and the process repeated to construct the vertices of for and .
Construction of is equivalent to finding the vertex adjacencies of ; an elementary way to do this is to find pairs of points separated by distance . However it was not efficient to identify boundary vertices in this manner, and we used KD trees to improve run times for this aspect of the problem. With the adjacencies and boundary edges identified it is easy to weight these and construct the matrix . The matrix of the Dirichlet Laplacian is obtained by deleting the rows and columns corresponding to boundary vertices.
5.2. The eigenvalue counting function
The eigenvalue counting function for a general symmetric matrix is defined by
where eigenvalues are counted with multiplicity. This generalizes to operators with discrete spectrum in . It is well-known that if the operator in question is the (positive) Laplacian on a Euclidean domain then encodes a great deal of geometric information: for example, the Weyl law says that it grows as , where is the dimension and encodes the volume of the domain; more precise asymptotics can be obtained for certain classes of domains.
|Left:||Eigenvalue counting functions (blue) and (orange)|
|Right:||Log-Log plots of (blue) and (orange)|
The eigenvalue counting functions of and the Dirichlet Laplacian for the approximation are displayed in Figure 3, along with Log-Log plots that identify their growth behavior. It is immediately apparent that has two regimes, a low eigenvalue regime in which it fairly closely tracks the behavior of , though with a lower initial rate of growth, and a high eigenvalue regime in which its growth is entirely different. Moreover, there is an additional feature which is not apparent from the graphs. The largest eigenvalue of is , and the regime-change point on the graph occurs at the almost identical value . Graphs of the corresponding eigenfunctions are also very similar; they are both highly oscillatory, so their graphs look multi-valued, but the similarity in the macroscopic profile of the oscillations is readily apparent in Figure 4.
5.3. Eigenvectors and eigenvalues in the low eigenvalue regime
The eigenstructure of in the low eigenvalue regime holds few surprises. Comparing it to that for we found that the boundary Dirichlet form permitted many additional low-frequency configurations, some of which are shown in Figure 5. It is apparent from the counting function plots in Figure 3 that this difference in the growth of and decreases as increases.
The eigenvalues less than for each operator are in the table below. It is evident that for both operators some eigenvalues occur with multiplicity and some with multiplicity . This is true throughout the spectrum (not just in the low frequency regime) and is explained by symmetries of ; this may be verified by the argument given in . At higher frequencies these symmetries are difficult to see from the graphs but readily apparent in contour plots, as illustrated in Figure 6.
|Eigenvalues of||Eigenvalues of|
One striking observation about the spectra of and was the existence of eigenvectors and eigenvalues that were very similar in both. An example is shown in Figure 7, where it is apparent that with eigenvalue is very similar to with eigenvalue . Such “pairs” of eigenvectors, one from and one from , both with the same symmetries and similar eigenvalues, were found throughout the low eigenvalue regime. One possible explanation is that there might be Dirichlet eigenfunctions with symmetries that imply they are close to being eigenfunctions of .
To see this, let us try to compute how far a Dirichlet eigenfunction might be from being an eigenfunction of . Making a computation like that in Theorem 3.6 we see that a Dirichlet eigenfunction is such that for all . A general element of can be written as , where denotes the boundary values and is harmonic. Since has zero boundary trace we have
so would also be an eigenfunction of if for all . As importantly, we would expect to find an eigenfunction and eigenvalue near if was small for all . In consideration of Corollary 3.8 we might expect highly oscillatory to produce with support close to the boundary, where
is small because it has Dirichlet boundary values. This heuristic suggests that we need only consider slowly varying, and, moreover, that if the symmetries of are such that it is nearly orthogonal to a large enough subspace of harmonic extensions of slowly varying boundary value functions, then should be very close to being an eigenfunction of as well as . We do not know if this heuristic is the correct explanation for the observed phenomenon, or have estimates that would make the explanation precise.
5.4. Localization in the high eigenvalue regime
The high eigenvalue regime contains the eigenvalues of that are larger than the largest Dirichlet eigenvalue. In this regime we see a dramatic transition from eigenfunctions that are supported throughout the domain to eigenfunctions localized near the boundary. Figure 8 illustrates this change using contour plots for eigenfunctions of : the first row of images are of eigenfunctions just inside the regime, while the other three shows the progression of localization with increasing frequency. The last image is the highest eigenvalue eigenfunction of . Note that the onset of the localization phenomenon is rapid, occurring across just a few eigenfunctions, and that continues to sharpen as the frequency increases, but at a decreasing rate.
Several phenomena in physics involve the localization of eigenmodes (e.g. Anderson localization due to a random potential in the Schrödinger equation). A nice survey from a mathematical perspective is in . Some involve low frequency modes (e.g. weak localization due to irregular or complex geometry of the domain) and some involve high frequency modes (e.g. whispering-gallery modes, bouncing ball modes and focusing modes). The type that looks most similar to those observed for are the whispering-gallery modes. These are well understood for convex domains , but despite the fact that for many years similar modes have been observed in non-convex domains, including domains with pre-fractal boundaries, our understanding of why these occur is incomplete. One approach to high-frequency localized eigenmodes uses quantum billiards, but studying quantum billiards on sets with fractal boundaries seems to be a difficult problem . The only literature we are aware of for the Koch snowflake domain is [38, 39], and it does not include modes of the type we see for .
We believe that in our setting the high frequency boundary localizations are not due to a whispering gallery effect, but rather to the relatively weak coupling between the boundary energy and the domain energy. We can see this, at least heuristically, by repeating the sort of calculation done at the end of subsection 5.3 above, but this time considering how close a high frequency eigenfunction of the boundary energy might be to an eigenfunction of . Suppose satisfies for all . If we extend and harmonically to , obtaining and respectively, and add an arbitrary to so that ranges over , we obtain
Now if is of high frequency it is highly oscillatory and it is not unreasonable to expect that its piecewise harmonic approximation at large scales is small. In consideration of Corollary 3.8 one might expect that then is small, and therefore . Moreover, we have
and again we see from Corollary 3.8 that one should expect to be small if is high frequency, leading to . Combining this with the corresponding statement for the energies and the assumption that was an eigenfunction of we have
which, if exactly true, would say that was an eigenfunction of with eigenvalue . This argument is consistent with what we observe, especially the fact that oscillations of the most localized of the high frequency eigenvectors of appear to have wavelength approximately the length of an edge of and are indeed localized very closely to the boundary curve, as seen in Figure 10. However we do not have precise arguments or estimates to justify this heuristic reasoning. A less regular looking eigenfunction localized on the boundary is shown in Figure 9.
6. A landscape approach to high frequency localization
Although many aspects of wave localization problems remain open, pioneering work of Filoche and Mayboroda [17, 18] has led to a great deal of progress on low-frequency localization problems. For a suitable elliptic differential operator on a bounded open set , they introduce a “landscape” function by , where is the Dirichlet Green’s function, and show that an normalized Dirichlet eigenfunction with satisfies . In consequence, if there is a region where is small (a valley) that is enclosed by a set where is large (a set of peaks), then a low frequency eigenfunction that is non-zero in the valley must be confined to that valley. This reduces the problem of low frequency eigenfunction localization to studying the structure of the landscape function, and especially its level sets. The latter involves a great deal of sophisticated mathematics, but is numerically very simple to implement.
Our purpose here is to write a high frequency variant of the argument of Filoche and Mayboroda, and to see that its numerical implementation predicts the high frequency localization seen in our data. We do so by introducing a high frequency landscape function associated to a linear map .
For is a linear map expressed as a matrix with respect to the standard basis, we let denote the linear map with matrix and define the high frequency landscape vector to be . Thus with .
We remark that the high frequency landscape depends on whereas the low frequency landscape function depended on . The following elementary theorem is the analogue of the Filoche-Mayboroda bound using the high frequency landscape function.
Let be an eigenvector of corresponding to an eigenvalue normalized such that . Then for each ,
Since we may compute using the standard basis