# Assessment of Models for Pedestrian Dynamics with Functional Principal
Component Analysis

Mohcine Chraibi  Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH,
52428 Jülich, Germany<EMAIL_ADDRESS>Tim Ensslen  Department of
Mathematics and Computer Science, Bergische Universität Wuppertal, 42119
Wuppertal, Germany<EMAIL_ADDRESS>Hanno Gottschalk
<EMAIL_ADDRESS>Mohamed Saadi<EMAIL_ADDRESS>Armin
Seyfried<EMAIL_ADDRESS>Department of Civil Engineering, Bergische
Universität Wuppertal, Pauluskirche 7, 42285 Wuppertal, Germany

###### Abstract

Many agent based simulation approaches have been proposed for pedestrian flow.
As such models are applied e.g. in evacuation studies, the quality and
reliability of such models is of vital interest. Pedestrian trajectories are
functional data and thus functional principal component analysis is a natural
tool to asses the quality of pedestrian flow models beyond average properties.
In this article we conduct functional PCA for the trajectories of pedestrians
passing through a bottleneck. In this way it is possible to asses the quality
of the models not only on basis of average values but also by considering its
fluctuations. We benchmark two agent based models of pedestrian flow against
the experimental data using PCA average and stochastic features. Functional
PCA proves to be an efficient tool to detect deviation between simulation and
experiment and to asses quality of pedestrian models.

###### keywords:

pedestrian dynamics; statistical analysis; comparison with experiment;
functional PCA; model quality

###### PACS:

89.75-k Complex Systems, 50.40-a Stochastic Models

††journal: Physica A

##  1 Introduction

Most of force-based models qualitatively describe the movement of crowds of
pedestrians. Self-organization phenomena e.g., lane formations [7, 6, 29],
oscillations at bottlenecks [7, 6], clogging at exit doors [6, 29] etc., are
reproduced. From a physical point of view it is of interest how simple model
reproduce qualitatively self-organization phenomena of driven multi-particle
systems. That contributes to a better understanding of the investigated
systems and the essential interactions. In addition numerical simulations
basing of these models are used to address safety related issues, concerning
e.g. design and conception of escape routes in buildings [26, 28] or optimal
organization of mass events or public transport facilities (VISWalk [31],
Legion [32], $\ldots$). For such utilization a thorough quantitative
validation of the models is obligatory to ensure a reliable layout,
dimensioning or evaluation of pedestrian facilities. In most known cases this
is fulfilled by reproducing the fundamental diagram [24, 15, 2, 9] or
measuring the flow through bottlenecks [8, 2, 13]. An overview of quantitative
validation of models by means of the fundamental diagram is given in [21]. On
one hand, the common point between these quantitative methods is the fact that
they are based on calculating specific traffic quantities, e.g. density, flow
and velocity. On the other hand, these measurements are performed based on
locally averaged values over time or space. [30] and [22] provide examples how
the measurement methods could influence the resulting empirical relations of
such granular and heterogeneous systems of finite size. The differences
between the measurement methods suggest that important information on the
system may be lost during the measurement process. Moreover state of the art
models describe pedestrian dynamics on a more detailed level by simulating
trajectories of every single pedestrian allowing in principle a validation
method assessing average pedestrian or traffic flow behavior, but also
accounting for the amount and typical nature of fluctuation around this
average.

A first methodology based on exploiting information of individual trajectories
was introduced in [10] to calibrate the social force model. While one
pedestrian was moved according to the model the others were moved according to
real trajectories. By means of an evolutionary algorithm the deviations of the
resulting trajectories from the experimental ones was used to calibrate the
parameters of the model. But this approach doesn’t allow an assessment of the
quality of a model.

While an abundance of agent-based models in the field of pedestrian and
traffic dynamics were developed in the last years [3, 17] the question of
systematic comparison of experimental evidence and model generated results has
not caught the same attention. This would however be important for the ranking
of models into more or less adequate ones. As argued above methodology of the
evaluation should provide a comparison of model results and empirical data
corresponding to the level of detail of the model. It is desirable that such a
validation method should not only be able to asses average pedestrian or
traffic flow behavior, but also account for the amount and typical nature of
fluctuation around this average.

Among the difficulties in this validation process is the fact that in agent-
based pedestrian or traffic flow data is functional, i.e. to each individual
we associate data in the infinite dimensional space of trajectories $x(t)$.
The adequate statistical approach for the study of pedestrian or traffic flow
data is thus the well established method of functional data analysis [20]. In
this method, the variation in the trajectories of different agents is
interpreted as random fluctuations. Thus, the measured or simulated
trajectories are interpreted as realizations of some stochastic process
$X(t)\in\mathbb{R}^{2}$, where $t$ stands for a time parameter and
$X(t)=X(t,\omega)$ tacitly depends on some random parameter $\omega$ from a
probability space $(\Omega,{\cal A},P)$. For more details the reader is
conferred to [1]. Although there are infinitely many trajectories available
for an agent to move from point $A$ to point $B$, it often turns out that a
few typical modes of variation around the average movement are responsible for
the bulk of fluctuation of trajectories between different individuals. As a
classical method in the analysis of functional data, the functional principal
component analysis (PCA) is the standard method to find and analyze these
typical variations.

The scope of this article is to use functional PCA analysis to study the
performance of agent-based models of pedestrian motion with respect to
experimental data. In order to demonstrate the methodological approach, two
models – social force model (SFM) [14] and generalized centrifugal force model
(GCFM) [2] – are used to simulate pedestrian movement through a bottleneck of
the same dimensions. In the following we apply functional PCA using the open
source extension fda by Ramsey, Hooker and Graves [19] to conduct the
analysis. We present the results and give a detailed comparison of average
values for locations and velocities and their respective principal components.
For the latter we separately compare strength, distribution of total
variation, and morphology of principal components.

We show that functional PCA in fact can be used to make statistically
significant statements about model quality. Functional PCA reveals significant
deviation between both models and the experiment already on the level of
average values. While the morphology of principal components for locations is
more or less adequately represented by both models, there are significant
deviations in the strength of fluctuations around the mean behavior with the
GCFM model underestimating the experimentally observed fluctuations while the
SFM mostly overestimates fluctuation strength. These empirical observations
can be confirmed with statistical testing for significance using the PCA-
bootstrap methodology [4, 5].

In this article, for the fist time we combine functional PCA in the sense of
[20] with the bootstrapping of scores in order to calculate the fluctuations
of specific statistics that describe and distinguish characteristic features
of fluctuations of individual pedestrian behavior in a crowd. Also on the PCA-
side, benchmarking and testing with specific statistics evaluated in
functional PCA is a new strategy, to the best of our knowledge.

The article is organized as follows. In section 2 we review the pedestrian
flow experiment [23] as the benchmark case for this study. Section 3 gives a
brief account on the SFM and the GCFM model. In Section 4 reviews the
functional PCA and its numerical implementation. Section 5 is the main part of
this article. After some introductory remarks on data formatting and smoothing
(Subsection 5.1), we compare average data for $x$ and $y$ position data
(Subsection 5.2) and velocities directed in the main direction of motion,
which is the $x$-direction. We then compare fluctuations strength via PCA
eigenvalues (Subsection 5.3) and morphology for the first PCA harmonics for
$x$\- and $y$\- positions and $x$-velocities (Subsection 5.4). Section 6
presents the PCA-bootstrap approach in the context of spline-based PCA
(Subsection 6.1) and applies this to total variation and Gini index
(Subsection 6.2) as well as the $L^{2}$-distance of the average trajectories
and the Hilbert Schmid distance of the empirical correlation functions
(Subsection 6.3). In Section 7 we summarize or findings and give some
conclusions on model quality in the specific case and general applicability of
functional PCA in the given context.

##  2 Experiment

In this work we use as a reference the experimental data extracted from the
experiment [23], that was performed in 2006 in the wardroom of the “Bergische
Kaserne Düsseldorf”. See Fig. 1.

Figure 1: The simulation set-up: pedestrians start from the shaded area and
move through the bottleneck ($l=4$ m, $h=4.5$ m, $b=6$ m and $w=0.9$ m). An
adjacent area of length $e=2.5$ m is added to consider the backward effect of
leaving pedestrians on those still in the bottleneck.

A waiting area was used to distribute the attendees before the start of each
run of the experiment. For simulation purposes we enlarge the area of the set-
up by an extra room of length $e$. This is necessary to take into
consideration the effects of pedestrians that leaved the bottleneck on the
pedestrians still in the system.

The flow through the bottleneck is measured as follows:

$J=\frac{N}{\Delta t},$ (1)

with $N$ the number of pedestrians and $\Delta t=t_{\rm last}-t_{\rm first}$
the time gap between the first and the last pedestrian passing the bottleneck
at the measurement line.

##  3 Models

Force-based models describe the movement of pedestrians as a superposition of
forces. Given the state variables of pedestrian $i$ at time $t$
$(\overrightarrow{x_{i}}(t),\overrightarrow{v_{i}}(t))$ and considering
Newton’s second law of dynamics the state of each pedestrian $i$ is defined
by:

$m_{i}\overrightarrow{a_{i}}(t)=\sum_{i\neq
j}\overrightarrow{f_{ij}^{r}}+\sum_{W}\overrightarrow{f_{iw}^{r}}+\overrightarrow{f_{i}^{d}},$
(2)

and

$\overrightarrow{v_{i}}(t)=\frac{d\overrightarrow{x_{i}}(t)}{dt},$ (3)

where $\overrightarrow{f_{ij}^{r}}$ denotes a repulsive force acting from the
$j^{th}$-pedestrian on the $i^{th}$-pedestrian, $\overrightarrow{f_{iw}^{r}}$
is a repulsive force emerging from borders, walls etc. and
$\overrightarrow{f_{i}^{d}}$ is a driven force. $m_{i}$ is the mass of
pedestrian $i$.

The superposition of the forces reflect the fact that pedestrians move towards
a certain point in space (e.g. an exit) and meanwhile try to avoid collisions
with each other or with walls and objects.

The driving force $\overrightarrow{f_{ij}^{r}}$ models, at low densities, an
exponential acceleration towards a desired speed $v_{0}$: The following
expression [18] is used:

$\overrightarrow{f_{i}^{d}}=m_{i}\frac{v_{0}\overrightarrow{e_{i}^{0}}-\overrightarrow{v_{i}}}{\tau},$
(4)

with a relaxation time $\tau$ typically equal to 0.5 s, and a desired
direction $\overrightarrow{e_{i}^{0}}$ of pedestrian $i$.

The repulsive force between pedestrians $\overrightarrow{f_{ij}^{r}}$ is
defined differently from one model to another [7, 24, 29, 2, 25, 11].

In this work we study a variation of the SFM and the GCFM. Both models are
microscopic and continuous in space. In the GCFM the agents have an elliptical
form with velocity-dependent semi-axes, whereas the shape of agents in the SFM
is circular. In the general case, the distance
$\parallel\overrightarrow{d_{ij}}\parallel$ is defined as the distance between
the borders of the ellipses $i$ and $j$ along a line connecting their centers.
See Fig. 2. For the SFM the semi-axis orthogonal to the movement direction is
equal to the other semi-axis in the direction of movement.

Figure 2: The effective distance $d_{ij}$ of two pedestrians represented by
two ellipses.

For simplicity we write $d_{ij}$ to denote the norm of the vector
$\overrightarrow{d_{ij}}$.

###  3.1 The social force model (SFM)

The SFM as originally published by Molnár [14] describes the movement of
circular agents as superposition of different factors e.g. influence of
neighboring pedestrians, walls, attractions and groups. In this work we reduce
the complexity of the model to a minimum, considering only the influence of
pedestrians and walls and assuming only circular potentials.

The repulsive force in the SFM between agent $i$ and $j$ is defined as

$\overrightarrow{f_{ij}^{r}}=-m_{i}k_{ij}A\exp(\frac{d_{ij}}{B})\overrightarrow{e_{ij}}.$
(5)

with

$\overrightarrow{d_{ij}}=\parallel\overrightarrow{x_{j}}-\overrightarrow{x_{i}}\parallel-
r_{i}-r_{j},$ (6)

and

$\overrightarrow{e_{ij}}=\frac{\overrightarrow{x_{j}}-\overrightarrow{x_{i}}}{\parallel\overrightarrow{x_{j}}-\overrightarrow{x_{i}}\parallel},$
(7)

with the parameters $A$ and $B$ the strength and the range of the force are
adjusted. The limited vision of pedestrians ($180^{\circ}$) is modeled by the
coefficient $k_{ij}$:

$k_{ij}=\Theta\Big{(}\overrightarrow{v_{i}}.\overrightarrow{e_{ij}}\Big{)}.$
(8)

$\Theta(\cdot)$ is the Heaviside function. The repulsive force between
pedestrians and static objects is defined similarly to (5).

###  3.2 The generalized centrifugal force model (GCFM)

The repulsive force in the GCFM is inversely proportional to the distance of
two ellipses representing moving pedestrians $i$ and $j$ and depends on their
relative velocity:

$\overrightarrow{f_{ij}^{r}}=-m_{i}\frac{\Big{(}\alpha
v_{0}+v_{ij}\Big{)}^{2}}{d_{ij}},$ (9)

where
$v_{ij}=\Theta\Big{(}(\overrightarrow{v_{i}}-\overrightarrow{v_{j}})\cdot\overrightarrow{e_{ij}}\Big{)},\;$
is the relative velocity. The use of the Heaviside function $\Theta(\cdot)$
ensures, that faster pedestrians are not effected by slower pedestrians. By
means of the parameter $\alpha$ the strength of the force can be adjusted. As
mentioned earlier in the GCFM the space requirement in the direction of
movement is modeled by the semi-axis

$a=a_{\text{min}}+\tau_{a}v_{i},$ (10)

with two parameters $a_{\text{min}}$ and $\tau_{a}$, whereas the lateral
swaying of pedestrians is modeled by the semi-axis

$b=b_{\text{max}}-(b_{\text{max}}-b_{\text{min}})\frac{v_{i}}{v_{0}}.$ (11)

###  3.3 Model parameters

As mentioned earlier the original SFM includes several forces e.g. physical
contact forces and attractive forces. For our purpose we use a simplified
version of the SFM as presented in Sec. 3.1. We choose $A=5$ N for pedestrian-
pedestrian interactions (5) and $A=7$ N for pedestrian-wall interactions. The
range of the function defined by the parameter $B$ in (5) was chosen to be
$0.08$ m for pedestrian-pedestrian interactions and $0.05$ m for pedestrian-
wall interactions. The parameter $\alpha$ in (9) is set to 0.2 for pedestrian-
pedestrian interactions and 0.33 for pedestrian-wall interactions. The desired
speed $v_{0}$ is set to $\mu=1.1$ m/s. For simplicity we set for both models
$m_{i}=1$ Kg. Table 1 gives a resume of the parameters used.

Parameter | Equation | Value
---|---|---
$A_{\rm ped}$ | (5) | 5 N
$A_{\rm wall}$ | Similar to (5) | 7 N
$B_{\rm ped}$ | (5) | 0.08 m
$B_{\rm wall}$ | Similar to (5) | 0.05 m
$\alpha_{\rm ped}$ | (9) | 0.2
$\alpha_{\rm wall}$ | Similar to (9) | 0.33
$\tau$ | (4) | 0.5 s
$v_{0}$ | (4) | 1.1 m/s
$m$ | (4) | 1 Kg
$a_{\tau}$ | (10) | 0.12 s
$a_{\rm min}$ | (10) | 0.15 m
$b_{\rm min}$ | (11) | 0.15 m
$b_{\rm max}$ | (11) | 0.2 m
Table 1: Parameter values in simulations with both GCFM and SFM.

The values chosen in Tab. 1 differ from the values published in other works
[12, 2]. Our choice of the above mentioned values is supported by qualitative
reasons, ensuring minimal overlapping among pedestrians, as well by
quantitative consideration of the flow through the bottleneck. See Fig. 3.

Figure 3: The flow through the bottleneck measured at the middle of the
corridor after the entrance to the bottleneck. The empirical value is a
reference value for the calibration of the GCFM and the SFM.

For safety relevant simulations a careful calibration of the used models is
needed. Having calibrated two different models based on usual qualitative and
quantitative criteria, we strive to apply a new technique to assert the
goodness of the investigated models and verify whether the aforementioned
validation is sufficient to ensure a trustworthy and safe use of the produced
simulations.

##  4 Functional PCA: Foundations

###  4.1 What is functional PCA?

In this section we give some details of functional PCA following [20]. The
principal component analysis uses the principal axis transformation for
multivariate, correlated numerical data using the (empirical) covariance
information between the single random variables. Eigenvalues then sort the
importance of the single eigenvectors (also called harmonics or modes)
according to the variance.

This concept needs to be adapted to the case where the observed data from $n$
individuals are functions – as it is the case with the trajectories of
pedestrians. The variability of the data can still be described in with the
eigenvalues and eigenvectors of the covariance function seen as an operator on
the function space of square integrable functions. Given the stochastic
process of random trajectories $X(t)$, with $t\in(0,L)$, the covariance
function is defined as

$C(s,t)=\mathbb{E}\left[(X(s)-\mathbb{E}[X(s)])(X(t)-\mathbb{E}[X(t)])\right],$
(12)

with $\mathbb{E}$ the expected value with respect to the underlying
probability space. This covariance function needs to be estimated out of the
data $x_{j}(t)$

$\hat{C}(s,t)=\frac{1}{n-1}\sum_{j=1}^{n}(x_{j}(s)-\bar{x}(s))(x_{j}(t)-\bar{x}(t)),$
(13)

with $x_{j}(t)$ the $j$-th observation is one realization of the random
process $X(t)$ and $\bar{x}(t)=\frac{1}{n}\sum_{j=1}^{n}x_{j}(t)$.

In the following we assume that average values have already been removed from
the stochastic signal, i.e. we consider transformed random quantities $X(t)\to
X(t)-\mathbb{E}[X(t)]$ with estimated observations $x_{i}(t)-\bar{x}(t)$. The
eigenvalues of $C(s,t)$ can then be calculated after solving the follow
eigenvalue equation

$\int_{0}^{L}C(s,t)\xi(t)\,dt=\lambda\,\xi(s).$ (14)

This results in a set of eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq
0$ and corresponding eigenfunctions $\xi_{i}(s)$. These eigenfunctions are
orthonormal $\int_{0}^{L}\xi_{i}(t)\xi_{j}(t)\,dt=\delta_{i,j}$, where
$\delta_{i,j}=1$ for $i=j$ and zero otherwise. The eigenfunctions and
eigenvalues can now approximately be determined from the observations
$x_{j}(t)$ by replacing $C(s,t)$ by its empirical counterpart $\hat{C}(s,t)$.

###  4.2 Numerical approximation

The problem (14) is an infinite dimensional eigenvalue problem and its
empirical counter part is potentially very high dimensional (of dimension
$n$). A frequently used method to make this problem numerically tractable is
to project the covariance matrix on the space spanned by some finite basis, -
e.g. a sufficiently fine B-spline or Fourier basis. Then one solves for the
eigenvalues and functions in the given finite dimensional space of basis
functions. Therefore we approximate the observed functions $x_{i}(t)$ with a
suitable linear combination of basis functions

$x_{i}(t)=\sum_{k=1}^{K}c_{i,k}\Phi_{k}(t)\Leftrightarrow
x_{i}(t)=\bf{C}\Phi(t),$ (15)

with basis function vector $\Phi(t)=(\Phi_{1}(t),\ldots,\Phi_{K}(t))^{T}$ and
coefficient matrix ${\bf C}=(c_{j,k})_{j=1,\ldots,n\atop k=1,\ldots,K}$
obtained e.g. by orthogonal projection of $x_{i}(t)$ to the space spanned by
the basis functions and a subsequent basis decomposition, in which case ${\bf
C}={\bf W}^{-1}{\bf v}_{j}$ with
$W_{i,j}=\left\langle\Phi_{i},\Phi_{j}\right\rangle=\int_{0}^{L}\Phi_{i}(t)\Phi_{j}(t)\,ds$
and $c_{j,k}=\langle
x_{j},\Phi_{k}\rangle=\int_{0}^{L}x_{j}(t)\Phi_{k}(t)\,dt$. Here some
numerical quadrature may be employed for the integrals involved in the
definition of $v_{j,k}$, whereas in most cases analytic formulae are available
for $W_{i,j}$. This projection method implies that the covariance function can
be approximated by

$\hat{C}(s,t)=\frac{1}{n-1}\Phi(s)^{T}{\bf\bar{C}}^{T}{\bf\bar{C}}\Phi(t).$
(16)

Here we use $\bar{C}$ for the matrix of coefficients of $x_{j}(t)-\bar{x}(t)$
with respect to the basis $\Phi(t)$, namely
$\bar{c}_{j,k}=c_{j,k}-\frac{1}{n}\sum_{j=1}^{n}c_{j,k}$. Now we expand the
eigenfunction with the same basis functions to a good approximation:

$\xi(t)=\sum_{k=1}^{K}b_{k}\Phi_{k}(t)=\Phi(t)^{T}{\bf b}.$ (17)

The approximate eigen value equation can be written as

$\int_{0}^{L}\hat{C}(s,t)\xi(t)\,dt=\frac{1}{n}\Phi(s)^{T}{\bf\bar{C}}^{T}{\bf\bar{C}}{\bf
W}{\bf b}=\lambda\Phi(s)^{T}{\bf b}=\lambda\xi(s).$ (18)

leading to the eigenvalue equation ${\bf\bar{C}}^{T}{\bf\bar{C}}{\bf W}{\bf
b}=\lambda{\bf b}\Leftrightarrow{\bf W}^{1/2}{\bf\bar{C}}^{T}{\bf\bar{C}}{\bf
W}^{1/2}{\bf u}=\lambda{\bf u},~{}{\bf b}={\bf W}^{-1/2}{\bf u}$, which can be
solved numerically. The result is a number of eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$ and coefficient vectors ${\bf
b}_{i}$ for approximate principal components $\xi_{i}(t)=\Phi(t)^{T}{\bf
b}_{i}$, for $i=1,2,\ldots,n-1$, which is the maximal rank of ${\bf\bar{C}}$.

###  4.3 Statistics of the eigenvalues

In this subsection, we discuss how to reduce the information from the set of
eigenvalues to a few significant characteristics. In particular we will focus
on code figures that measure the strength of fluctuations and their
concentration to a few, active modes.

The eigenvalue $\lambda_{i}$ represents the strength of fluctuations in the
respective mode of characteristic shape $\xi_{i}(t)$. The relative strength
$\rho_{i}$ of the variation in the mode $\xi_{i}(t)$ and the cumulative
relative strength $L_{j}$ up to the $j$-th mode $\xi_{j}(t)$ is given by

$\rho_{i}=\frac{\lambda_{i}}{\sum_{j=1}^{n}\lambda_{j}}\mbox{ and
}L_{j}=\sum_{i=1}^{j}\rho_{i}.$ (19)

Two quantities that can be derived form the eigenvalues of the PCA are of
special interest: First, the total variation strength is simply the sum of all
eigenvalues $\Lambda=\sum_{j=1}^{n}\lambda_{j}$ whereas the Gini index is a
measure of concentration that is build from the Lorenz curve quantities
$L_{j}$ via $G=\frac{2}{n-1}\sum_{j=1}^{n}(L_{j}-j/n)$. Geometrically, the
Gini index measures the area between the diagonal and the Lorentz curve, cf.
e.g. Figure 8 in the right panel. It is normalized such that it takes the
value one if only one mode is active and takes the value zero when all modes
are equally activated $\lambda_{1}=\lambda_{2}=\cdots=\lambda_{n}$. Note that
the order of $\lambda_{j}$ is descending in contrast to the usual definition
of the Gini index, where the order is ascending. As an alternative, one could
also consider the entropy of the distribution of the total activity to the
single modes. The result of the observation however remain largely unchanged.

###  4.4 Deviation measures

In this section we derive some quantities that can be used to measure the
distance between one set of functional data and another such data set. In
particular we will utilize these distances for benchmarking models with
respect to their distance to the experiment. Two distance measures will be
employed in the following: First, the mean quadratic deviation between the
average trajectories of the model on the one hand and the experimental data on
the other. Secondly, we consider the Hilbert-Schmidt norm between the
respective empirical covariance functions as a measure of the distance of the
fluctuation behavior of the experiment and the simulation. In the following we
work with the data after projection to a finite spline basis $\Phi_{k}(t)$.

We start with the mean quadratic difference in the average behavior. The mean
of the observed function $x_{i}(t)$ is:

$\overline{x}(t)=\frac{1}{n}\sum_{j=1}^{n}x_{j}(t)=\sum_{k=1}^{K}\biggl{(}\underbrace{\frac{1}{n}\sum_{j=1}^{n}c_{jk}}_{=C_{k}}\biggr{)}\Phi_{k}(t)=\sum_{k=1}^{K}C_{k}\Phi_{k}(t).$
(20)

The mean quadratic distance, the squared $L^{2}$ norm, of the difference
between $\overline{x}(t)=\sum_{k=1}^{K}C_{k}\Phi_{k}(t)$ and
$\overline{y}(t)=\sum_{k=1}^{K}C_{k}^{{}^{\prime}}\Phi_{k}(t)$ is:

$\|\overline{x}(t)-\overline{y}(t)\|_{2}^{2}=\left\langle\sum_{k=1}^{K}(C_{k}-C_{k}^{{}^{\prime}})\Phi_{k}(t),\sum_{k=1}^{K}(C_{k}-C_{k}^{{}^{\prime}})\Phi_{k}(t)\right\rangle=({\bf
C}-{\bf C}^{\prime})^{T}{\bf W}({\bf C}-{\bf C}^{\prime}).$ (21)

We now derive formulae for measuring the distance between experiment and
simulation in the covariance structure. Let
$\hat{D}(s,t)=\hat{C}(s,t)-\hat{C}^{\prime}(s,t)=\sum_{k=1}^{K}\sum_{l=1}^{K}(\underbrace{\bar{c}_{k,l}-\bar{c}^{\prime}_{k,l}}_{=d_{k,l}})\Phi_{k}(t)\Phi_{l}(s)$
be the difference of covariance functions. The Hilbert-Schmidt norm of
$\hat{D}(s,t)$ is:

$\|\hat{D}\|_{HS}^{2}=\int_{0}^{L}\int_{0}^{L}\left(\hat{D}(s,t)\right)^{2}dsdt=\operatorname{Tr}\bigl{(}(\bf{D}{\bf
W})^{T}\bf{D}{\bf W}\bigr{)}.$ (22)

Here $\operatorname{Tr}(A)$ stands for the trace of the matrix $A$ and
$\bf{D}$ is the matrix with entries $d_{j,k}$.

##  5 PCA Results

In this section we show that functional PCA is a useful tool for detailed
validation of models for pedestrian dynamics. Ideally, the variability in the
data can be described with the help of the PCA with a few principal
components. These main components can be interpreted as the schemes for the
deviation of individual trajectories from the mean flow. This allows a
comparison of simulated and experimental data beyond averaged flow features.
Therefore, we apply the PCA to the experimental and simulated data from two
models SFM and GCFM and compare the results. Here we apply the PCA for $x$ and
$y$ coordinates over time separately, as this approach is somewhat more
accessible to the interpretation. For the alternative approach of jointly
analysing $x$ and $y$ trajectories and a discussion of the pros and cons of
both approaches, see [20, 19].

The analysis of the data is based on the R package fda developed by J.O.
Ramsay et al. [19] with some minor extensions by the authors.

###  5.1 Preparation of the Data

The pedestrian trajectory data in the experiment is recorded electronically
with video tracking at the rate of 25 frames per second. A total number of 149
trajectories has been recorded. Likewise, the SFM and GCFM models have been
simulated with 25 time steps per second using Euler integration with a time
step $\Delta t=0.01$ s. For both models, a total of 149 trajectories have been
simulated. See the trajectories in Figure 4.

Figure 4: XY plots of pedestrian trajectories generated by the SFM model
(left), experiment (middle) and GCFM model (right).

For the analysis, the pedestrian motion has to be stationary, i.e. all agents
move under the same conditions.

Furthermore, the considered time interval has to be the same for all agents.
Therefore we account only to the pedestrians, who need more than 12 seconds to
reach the exit. Furthermore, the trajectory is tracked only two seconds after
the passage of the bottleneck entrance. The analysis of the experimental data
and the models thus is based on trajectories in individual time intervals that
range from $-12$ s before passage through the bottleneck by the individual to
$+2$ s afterwards. This makes a total time of $14$ s. In order to avoid
negative time values, we start each pedestrian trajectory at time $t=0$ s such
that passage through the door for each individual occurs at $t=12$ s, exactly.
From now on we work with this time scale. Figure 5 visualizes the formatting
steps. After reformatting, a total of 118 pedestrian trajectories were
available for the experiment, 121 for the SFM and also 121 for the GCFM model,
respectively.

Figure 5: Plots of $x$-coordinates of pedestrians (experimental data): raw
data (left), stationary data (middle) and with individual time in the interval
[0,14] s (right). One Frame in the first two panels corresponds to 1/25 s.

The experimental data contain the swaying caused by the bipedal locomotion of
pedestrians in combination with the tracking of markers on the head. But the
SFM as well as the GCFM model the movement of the centre of mass neglecting
the bipedal locomotion, which produce swaying-free trajectories. Therefore, we
smooth the data before the analysis in order to filter out the lateral
swaying. It turns out that the regression with a B-spline basis containing 10
elements with nodes equally distributed over the underlying time interval
$[0,14]$ effectively removes swaying while properly reproducing the other
features of the individual’s trajectories, see Figure 6.

Figure 6: Smoothing of a experimental trajectory (dotted blue) with a B-spline
basis of dimension 10 (solid red).

###  5.2 Average trajectories

As PCA components describe variation around some mean value, it is essential
to analyse average functions $\bar{x}(t)$ and $\bar{y}(t)$.

Figure 7: Average curves for position vs time for the SFM (left), experiment
(middle) and GCFM (right): $x$-ccordinate (top) and $y$-coordinate (right).

Figure 7 shows the mean functions $\bar{x}(t)$ and $\bar{y}(t)$ of the $x$\-
and $y$-components for the experiment and the models.

When we examine the $x$-component, we identify clearly discrepancies between
the experiment and both models. The average pedestrian in the experiment shows
a nearly linear progress to reach the exit. The acceleration after passing the
bottleneck, i.e. the increase in the slope, is modest. In contrast, both SFM
and GCFM show a slower progress of the average pedestrian through the crowd
and a much more pronounced acceleration after the passage of the bottleneck.
While the latter deviation from experiment is about the same for both models,
the underestimation in the slope from the experimentally observed value is
bigger for the GCFM model. Thus both models overestimate the dwell time
indicating a missing anticipation and cooperation of the modelled pedestrians.
However the SFM produces this effect by a lesser amount.

The $y$-component is described by both models in a satisfactory manner, as
trend lines only move at a scale of a few centimeters from the center of the
bottleneck. At least for the simulated data this is due to the left hand -
right hand reflection symmetry of the agents in both models and the
(approximate) symmetry of initial positions with respect to the $y=0$ axis,
i.e. the center line through the bottleneck. In the experiment, a certain
asymmetric behavior is visible for the mean $y$-position of the trajectories
over time. In average, the pedestrians approach the bottleneck coming slightly
from the left seen from the direction of progress. Interestingly, this
asymmetry can not be traced back to the initial conditions, as these are the
same for the experimental and the simulated trajectories.

This behavior is absent in both models, which have left-right symmetries in
their respective constituting equations.

###  5.3 PCA eigenvalues

Having analyzed the average behavior, we now turn to the question, how well
the models describe fluctuations in pedestrian data around the averages. We
start with the absolute strength of PCA variability, which is represented by
the PCA eigenvalues $\lambda_{i}$, $i=1,\ldots,10$, as we are using a 10
dimensional spline basis. At the same time, we also consider the cumulative
relative strength $\varrho_{j}$ in order to measure the concentration or
dispersion of variability in experimental or the simulated data.

Figure 8: Left: Barplot of absolute fluctuation strength (eigenvalues) of the
10 PCA-modes (harmonics) for $x$-position over time for the SFM-model,
Experiment and GCFM-model. Right: Cumulative relative strength of PCA-modes
over all 10 harmonics.

We first consider PCA mode strength for $x$-position over time as given in
Figure 8. These modes describe typical patterns of pedestrians lagging behind
and being in front of the trajectory of the average pedestrian. The SFM
overestimates the total amount of statistical deviation in $x$-position from
the average $x$-position by approximately 12% as compared with the experiment.
Also, concentration of variability in the first mode is slightly lower than in
the experiment. In the GCFM, the total level of $x$-position variation is
underestimated by 37% of the total variation. The relative concentration in
the fist mode is higher as in the experimental data by an amount comparable to
the SFM, but in the opposite direction. For the values of total variation and
Gini indices confer Table 2.

Figure 9: Left: Barplot of absolute fluctuation strength (eigenvalues) of the
10 PCA-modes (harmonics) for $y$-position over time for the SFM-model,
Experiment and GCFM-model. Right: Cumulative relative strength of PCA-modes
over all 10 harmonics.

In Figure 9 the PCA-modes for the statistical $y$-fluctuation around the
average $y$-position (essentially $y=0$) is displayed. The experiment and both
simulations all show that basically only one mode is active representing the
axially symmetric shape of the jammed area in front of the bottleneck. The
size of this area is underestimated by both models. The GCFM predicts a
pronouncedly reduced area in the $y$-direction covered by trajectories passing
the bottleneck in the next 12 seconds, showing a total $y$-variation of
$\approx 1/3(33.6\%)$ of the experimental data. The same figure of
underestimation of $y$-variation for the SFM compared with experiment is 57%.
Total variations and Gini indices can again be found in Table 2.

Tot.Var | SFM | Experiment | GCFM
---|---|---|---
$x$-position | 10863 | 9550 | 6835
$y$-position | 35469 | 61619 | 20749

Gini | SFM | Experiment | GCFM
---|---|---|---
$x$-position | 0.83 | 0.85 | 0.87
$y$-position | 0.89 | 0.89 | 0.89

Table 2: Total variation (left) and Gini indices (right) for the eigenvalues
of the PCA.

###  5.4 PCA modes

Figure 10: PCA components in for $x$-position over time for the SFM model
(left), experiment (middle) and GCFM model (right). The first three harmonics
are displayed.

Figure 10 shows the principal fluctuation components of the $x$-position for
the experiment and the models. From the aforementioned eigenvalue analysis
(Figure 8) we can conclude that fluctuations can be mainly described by the
first three principal components.

Let us now give an interpretation to the morphology of the PCA components. The
first principal component describes the tolerance between the initial
positions of the pedestrians 12 seconds before passing the bottleneck. Two
pedestrians reach the bottleneck at the same reference time $t=12$. Due to the
fact that some pedestrian starts with a higher or lower $x$-distance to the
bottleneck than the the others, at $t=0$, a statistical variation in
$x$-positions occurs. This could be called “slipping-through-effect” because
faster pedestrians find more favorable configurations of fellow pedestrians
ahead which allows a faster passage through the crowd.

The second and third principal component describe an effect which we can
associate to long stop and go behavior in different lanes in a traffic jam:
One trajectory is temporary faster than the other, but afterwards it is the
other way round. In the case of the experiment and the GCFM, the third
principal component also shows different velocity patterns after the
bottleneck.

The morphological comparison of the experimental data with the SFM and GCFM
shows that the points of intersection of the first two principal components
are nearly at the same times. Thus both models reproduce the qualitative
behavior of statistical fluctuations in the pedestrians $x$-positions over
time quite well. The main difference thus lies in the different activation
strength of the “slipping-through” mode.

Figure 11: PCA components in the y direction for the SFM model (left),
experiment (middle) and GCFM model (right). The first three harmonics are
displayed.

Figure 11 shows the PCA of the $y$-components for the experiment and the
models. We observe that the variability of the data can be mainly described by
the first principal component which represents the shape of the crowd in front
of the door. The first principal component of both models describe the
experiential data acceptably well. Also the higher modes are of quite similar
shape, although they should be neglected since they hardly contribute to the
total variation.

###  5.5  Evaluation of deviation measures

Lastly in this section, we want to compare the deviation measures of the
respective simulation model with the experiment. The results are summarised in
3.

$L^{2}$-norm | EXP-GCFM | EXP-SFM
---|---|---
$x$-position | 169.8 | 166.0
$y$-position | 15.34 | 13.20

HS-norm | EXP-GCFM | EXP-SFM
---|---|---
$x$-position | 2291 | 1282
$y$-position | 4449 | 2880

Table 3: Deviations between experimental data and models data using the
$L^{2}$-norm (left) Hilbert-Schmidt norm (right).

The average squared distance of $x$\- and $y$\- coordinates of the
trajectories as function over time is of the same order of magnitude for both
models.

A slight advantage can however be attributed to the SFM-model. This effect is
even more pronounced in the Hilbert-Schmidt norm that measures the distance to
the experiment in the fluctuation structure of measured and simulated data.

##  6 Statistical inference based on the bootstrap

In the previous section, we evaluated the total variation and the Gini
coefficient or deviation measures for the average behavior and the fluctuation
structure ($L^{2}$\- and HS-norms, respectively) in order to compare simulated
and experimental data. This descriptive approach however leaves open the
question, to which extent these findings depend on the intrinsic stochastic
nature of pedestrian trajectories and to which extent they are due to
structural differences between simulated agents in the models and real
pedestrians observed in the experiment. In the present section, we describe
and apply a simulation-based test procedure in order to clarify, to what
extent the observed differences between models and experiment are
statistically significant.

###  6.1 Bootstrapping PCA scores

As the basis of our statistical testing procedure, we use the bootstrap over
the matrix of principal components from [4, 5]. We now shortly describe the
bootstrap approach. Given, e.g., the $i$-th $x$-value of the trajectory over
time, $x_{i}(t)$, the score of this trajectory with respect to the principal
component $\xi_{j}(t)$ is

$\displaystyle s_{ij}^{(x)}$ $\displaystyle=\langle
x_{i}(t)-\bar{x}(t),\xi_{j}(t)\rangle$
$\displaystyle=\left\langle\sum_{l=1}^{K}\underbrace{(c_{i,l}-C_{l})}_{{\bf\Delta}_{i,l}}\Phi_{l}(t),\sum_{k=1}^{K}{\bf
b}_{jk}\Phi_{k}(t)\right\rangle$ (23) $\displaystyle=({\bf\Delta}^{T}{\bf
W}{\bf b}^{T})_{i,j}~{}.$

By construction, the scores $s_{ij}^{(x)}$ and $s_{i,j^{\prime}}^{(x)}$ are
(linearly) uncorrelated for $j\not=j^{\prime}$. Neglecting potential higher
order correlations, we construct a virtual bootstrap sample from the scores of
the experimental data by drawing with replacement, for $i,j$ fixed,
$s_{i,j}^{(x,{\rm boot})}$ from the $N$ samples $s_{i,j}^{(x)}$ of the
original scores with respect to the $j$-the principal component $\xi_{j}(t)$.
Doing this independently for $i=1,\ldots,N$ and $j=1,\ldots,K$ (remember, in
our case $N=110$ is the number of experimental pedestrian trajectories and
$K=10$ the number of principal components) we obtain the $N\times K$ bootstrap
score matrix ${\bf s}^{(x,{\rm boot})}$. The corresponding bootstrapped
trajectories then are

$\displaystyle x_{i}^{({\rm boot})}(t)$ $\displaystyle=[{\bf s}^{(x,{\rm
boot})}\Phi(t)]_{i}+\bar{x}(t),~{}~{}i=1,\ldots N.$ (24)

Figure 12 shows the $x$-coordinates plots of pedestrian trajectories by
bootstrapped and experimental data.

Figure 12: Plots of $x$-coordinates of pedestrians by bootstrapped data (left)
and experimental data (right).

With this virtual data set, the PCA analysis is then repeated. In particular,
we obtain bootstrapped quantities for total variation and Gini index, as well
as distance measures for the average behavior of the actual experiment and its
virtual bootstrap replica. This entire process is then repeated a sufficiently
high number of times, such that $p$-values in the range of usual significance
levels $\approx 1-5\%$ can safely be determined. Here we generate $10^{4}$
bootstrap samples, each containing $N=118$ virtual trajectories, and thereby
obtain a simulated distribution for each of the aforementioned quantities.

The same iteration is repeated for the $y$-coordinate and the $x$-coordinate
of the velocity $v_{x}$.

For statistical testing, we generate two-sided confidence intervals for the
total variation and the Gini coefficient and left open confidence intervals
for the distance measures based on the empirical distributions of the
respective quantities. If the related quantities for the SFM and GCFM model
are not contained in these confidence regions, we consider this as a positive
test result for a deviation between experiment and model.

###  6.2 Testing Gini indices and total variations

One of the advantages of using this bootstrap technique is to have the
opportunity to examine the distributions of Gini indices and the total
variations. We compute for every bootstrap sample the Gini indices and total
variations by the experimental data for $x$-coordinates, $y$-coordinates and
$x$-velocities. Afterwards, we are able to compute their empirical cumulative
distribution functions (ECDF). Figure 13 shows the ECDF of Gini indices and
total variations of $x$-coordinates by bootstrapped experimental data. The
blue lines show the values of Gini index and total variation by original
experimental data. The corresponding $p$-values, i.e. the critical level of
statistical significance where the difference between model and experiment
becomes significant, are calculated on the basis of two-sided confidence
regions of the bootstrapped distribution. The $p$-values are summarized in
Table 4. Note that $p$-values below $10^{-3}$ become numerically unreliable
for $10^{4}$ bootstrap samples and are set to zero.

Figure 13: ECDF of Gini indices (left) and ECDF of total variation (right)
from bootstrapped experimental data. Green, grey and blue vertical lines mark
the values for the SFM, experiment and GCFM, respectively.

p-values for Gini index
---
| EXP/GCFM | EXP/SFM
$x$-position | 0 | 0.055
$y$-position | 0.610 | 0.007

p-values for TotalVar
---
| EXP/GCFM | EXP/SFM
$x$-position | 0.034 | 0.297
$y$-position | 0 | 0

Table 4: The p-value of Gini index and total variations by bootstrapped
experimental data.

The $p$-values in Table 4 show that statistical testing reveals significant
differences between experiment and model with respect to several Gini indices
and total variations. Only the Gini index of the GCFM and the total variation
of the SFM do not produce significant differences (significance level $5\%$).

###  6.3 Tests based on deviation measures

We are now interested to measure the deviations between bootstrapped and
experimental data. Firstly we compute for every bootstrap sample the $L^{2}$
norm of the difference between the mean trajectory of the bootstrap sample
based on the virtual experimental data and the mean trajectory of the
experiment. In this way we obtain an empirical distribution of the
$L^{2}$-norm distance due to natural fluctuation inside the experiment. This
is then compared with the $L^{2}$-norm distance of average trajectories
between the experiment and the SFM and GCFM model. The same procedure is also
carries through for the Hilbert-Schmidt distance between the estimated
correlation functions, see Figure 14.

Again, $p$-values are calculated as in the previous subsection, however this
time we have to use one-sided regions of confidence for statistical testing.
The $p$-values are displayed in Table 5.

p-values for $L^{2}$ norm
---
| EXP/GCFM | EXP/SFM
$x$-position | 0 | 0
$y$-position | 0.512 | 0.576

p-values for HS norm
---
| EXP/GCFM | EXP/SFM
$x$-position | 0.122 | 0.434
$y$-position | 0 | 0

Table 5: The p-values for $L^{2}$ norm and Hilbert-Schmidt norm .

Figure 14: The ECDF of the $L^{2}$-norm difference of average $x(t)$ values
(left) and the Hilbert Schmidt (HS) norm difference between covariance
functions of $x(t)$ fluctuations of experimental and bootstrapped data. Green
and blue vertical lines show $L^{2}$\- and HS-norm distance between the
experiment SFM and GCFM models respectively.

All statistical tests of the $L^{2}$-norm distance between the average
$x$-values and $v_{x}$-values of the experiments and the models are highly
significant. We find thus a clear indication that model and experiment are
statistically distinguished. Not unexpectedly, the situation is different for
the average of the $y$-values over time. Here, due to the axial symmetry of
the experimental set up, no major differences of the average $y$-trajectory
can be observed. This also shows that the slight asymmetry in the average
$y$-trajectory in the experimental data is not statistically significant.

With regard to the fluctuation structure, the $x$-position fluctuations
encoded by the empirical correlation function can not be easily distinguished
between the SFM model and the experiment. The difference between the
fluctuation of the $x$-trajectories of the GCFM and the experiment also shows
a very marginal $p$-value of $\approx 12\%$, which is not significant when
compared to the usual $5\%$-level of significance. All other fluctuation
structures significantly differ between the experiment and the models.

##  7 Summary

The functional PCA has been applied as a diagnostic tool to assess model
quality for agent-based simulations of pedestrian flows with respect to
average behavior and beyond. Here we applied it to experimentally measured
pedestrian trajectories passing through a bottleneck and agent trajectories
simulated by two different force-based models. Both models are a-priori
calibrated to satisfy qualitative and quantitative criteria.

Already in the analysis of mean flow behavior, the PCA reveals considerable
and statistically significant deviation of both models from the experiment. In
the $x$-direction, the SFM and the GCFM predict slower progress and lower
velocities leading to a longer dwell time before the bottle neck, as
experimentally observed, although total through flow is being produced
correctly. This effect is more pronounced for the GCFM average behavior so
that the SFM model reproduces the average behavior of the experiment – for the
given set of parameters – in a relatively better way.

Coming to the statistical variations simulated in the SFM for $x$-position, we
find a quite reasonable match in the qualitative behavior (the PCA mode
shapes) of the $x$\- and $y$-position over time. A certain deviation in the
quantitative variation strength is observed as well. Again the SFM predictions
are a little closer, when measured in terms of the deviation in the total
variation. Also, the concentration of variability to dominating modes is under
estimated by the SFM and over estimated by the GCFM for the $x$-variations,
while in the velocity variations both models show a higher degree of
concentration as compared with the experiment.

Summarizing the PCA gives none of the models a clear “pass”, as the SFM and
the GCFM both significantly differ from the experiment, although both models
were validated quantitatively with respect to the experimental flow through
the bottleneck. The SFM however performs relatively better than the GCFM,
which is mostly due to a gradually better prediction of the average
$x$-positions and $x$-velocities.

Caution is needed when applying these (and presumably also other) models to
evacuation studies, as the evidently do not capture all features of real life
pedestrian flows, as has been shown by our functional data analysis. The
overall picture of all the qualitative metrics derived from the PCA however
slightly favours the SFM as the more accurate model over the GCFM, given our
set of model parameters.

In this case study, applying functional PCA for the first time to pedestrian
flows, we have thus shown that it is in fact a useful tool to benchmark and
statistically test agent based pedestrian flow models. Given the amount of
deviation between experiment and model, it is certainly of interest to use
this methodology for the future refinement of pedestrian flow models and for
the critical assessment of models used by practitioners.

## Acknowledgments

We would like to thank Stefan Rosbach for his preparatory work on the
functional PCA analysis for pedestrian data.

## References

  * [1] Bauer, H., Probability Theory, de Gruyter 1995.
  * [2] Chraibi, M., Seyfried, A., and Schadschneider, A., The generalized centrifugal force model for pedestrian dynamics, _Phys. Rev. E_ 82 (2010) 046111\.
  * [3] Traffic and Granular Flow ’13, Eds. M. Chraibi, M. Boltes, A. Schadschneider and A. Seyfried, Springer (2015)
  * [4] Diaconis, P. and Efron, B., Computer-intensive methods in statistics. Scientific American, 248 No. 5 (1983) 116–130.
  * [5] Fisher, A., Caffo, B., Schwartz, B. and Zipunnikov, V., Fast, Exact Bootstrap Principal Component Analysis for $p>1$ million, arXiv:1405.0922v3.
  * [6] Helbing, D., Collective phenomena and states in traffic and self-driven many-particle systems, _Comp. Mater. Sci._ 30 (2004) 180–187.
  * [7] Helbing, D. and Molnár, P., Social force model for pedestrian dynamics, _Phys. Rev. E_ 51 (1995) 4282–4286.
  * [8] Hoogendoorn, S. P. and Daamen, W., Microscopic Parameter Identification of Pedestrian Models and its Implications to Pedestrian Flow Modeling (Transporttation Research Board Annual Meeting, 2006).
  * [9] Jian, M., Weiguo, S., Jun, Z., Siuming, L., and Guahgxuan, L., k-Nearest-Neighbor interaction induced self-organized pedestrian counter flow, _Physica A: Statistical Mechanics and its Applications_ 389 (2010) 2101–2117.
  * [10] Johansson, A., Helbing, D., and Shukla, P. K. Specification of the social force pedestrian model by evolutionary adjustment to video tracking data, _Advances in Complex Systems_ 10 (2007) 271–288.
  * [11] Karamouzas, I., Skinner, B., Guy, S. J, A universal power law governing pedestrian interactions, _Phys. Rev. Lett._ 113 (2014) 238701
  * [12] Lakoba, T. I., Kaup, D. J., and Finkelstein, N. M., Modifications of the Helbing-Molnár-Farkas-Vicsek social force model for pedestrian evolution, _Simulation_ 81 (2005) 339–352.
  * [13] Liao, W., Chraibi, M., Seyfried, A., Zhang, J., Zheng, X. and Zhao, Y. Validation of FDS+Evac for pedestrian simulations in wide bottlenecks, _Intelligent Transportation Systems (ITSC), 2014 IEEE 17th International Conference_ , (2014) 554–559.
  * [14] Molnár, P., Modellierung und Simulation der Dynamik von Fußgängerströmen, Ph.D. Thesis. University Stuttgart 1996.
  * [15] Nishinari, K., Sugawara, K., and Kaza, Modelling of self-driven particles: Foraging ants and pedestrians, _Physica A_ 372 (2006) 132–141.
  * [16] Parisi, D. R. and Dorso, C. O., Morphological and dynamical aspects of the room evacuation process, _Physica A_ 385 (2007) 343–355.
  * [17] Daamen, W., Duives, C. and Hoogendoorn, S. P. (eds.) _The Conference on Pedestrian and Evacuation Dynamics 2014 (PED 2014)_
  * [18] Pipes, L. A., An operational analysis of traffic dynamics, _J. Appl. Phys._ 24 (1953) 274 – 281.
  * [19] Ramsay, J., Hooker, G. and Graves, S. , Functional Data Analysis with R and Matlab, UseR! Series, Springer Dortrecht Heidelberg London New York, 2009.
  * [20] Ramsay, J. and Silverman, B., Functional Data Analysis, Springer Verlag New York, 2005.
  * [21] Seyfried, A. and Schadschneider, A., Fundamental Diagram and Validation of Crowd Models, in _Cellular Automata_ , eds. Umeo, H., Morishita, S., Nishinari, K., Komatsuzaki, T., and Bandini, S., _Lecture Notes in Computer Science_ , Vol. 5191/2008 (Springer, Berlin Heidelberg, 2008), pp. 563–566, doi:10.1007/978-3-540-79992-4.
  * [22] Seyfried, A.; Boltes, M.; Kähler, J.; Klingsch, W.; Portz, A.; Rupprecht, T.; Schadschneider, A.; Steffen, B. and Winkens,A. Klingsch, W. W. F.; Rogsch, C.; Schadschneider, A. and Schreckenberg, M. (Eds.)Enhanced empirical data for the fundamental diagram and the flow through bottlenecks _Pedestrian and Evacuation Dynamics_ Springer Berlin Heidelberg (2010) 145-156
  * [23] Seyfried, A. and Schadschneider, A., Empirical results for pedestrian dynamics at bottlenecks, in _Parallel Processing and Applied Mathematics_ , eds. Wyrzykowski, R., Dongarra, J., Karczewski, K., and Wasniewski, J., _Lecture Notes in Computer Science_ , Vol. 6068 (Springer, Berlin Heidelberg, 2010), pp. 575–584, http://dx.doi.org/10.1007/978-3-642-14403-5_62.
  * [24] Seyfried, A., Steffen, B., and Lippert, T., Basics of modelling the pedestrian flow, _Physica A_ 368 (2006) 232–238.
  * [25] Shiwakoti, N., Sarvi, M., Rose, G., and Burd, M., Animal dynamics based approach for modelling pedestrian crowd egress under panic conditions, _Transportation and Traffic Theory_ 17 (2011) 438–461.
  * [26] Schneider, V. and Könnecke, R. Simulating evacuation processes with ASERI, _Pedestrian and Evacuation Dynamics_ Springer (2002) 303–314
  * [27] Steffen, B. and Seyfried, A., Methods for measuring pedestrian density, flow, speed and direction with minimal scatter, _Physica A_ 389 (2010) 1902–1910.
  * [28] TraffGo HT GmbH, Handbuch PedGo 2, PedGo Editor 2, www.evacuation-simulation.com (2005)
  * [29] W. J. Yu, R. Chen, L. Y. Dong and S. Q. Dai, Centrifugal force model for pedestrian dynamics, _Phys. Rev. E_ , 72 (2005) p. 026112.
  * [30] Zhang, J., Klingsch, W., Schadschneider, A., and Seyfried, A., Transitions in pedestrian fundamental diagrams of straight corridors and t-junctions, _J. Stat. Mech._ (2011).
  * [31] http://vision-traffic.ptvgroup.com/de/produkte/ptv-viswalk/
  * [32] http://www.legion.com/


# Spectral heat content on a class of fractal sets for subordinate killed
Brownian motions

Hyunchul Park and Yimin Xiao

###### Abstract

We study the spectral heat content for a class of open sets with fractal
boundaries determined by similitudes in ${\mathbb{R}}^{d}$, $d\geq 1$, with
respect to subordinate killed Brownian motions via $\alpha/2$-stable
subordinators and establish the asymptotic behavior of the spectral heat
content as $t\to 0$ for the full range of $\alpha\in(0,2)$. Our main theorems
show that these asymptotic behaviors depend on whether the sequence of
logarithms of the coefficients of the similitudes is arithmetic when
$\alpha\in[d-\mathfrak{b},2)$, where $\mathfrak{b}$ is the interior Minkowski
dimension of the boundary of the open set. The main tools for proving the
theorems are the previous results on the spectral heat content for Brownian
motions and the renewal theorem.

††2020 Mathematics Subject Classification: 35K05, 28A80, 60G52††Key words:
spectral heat content, subordinate killed Brownian motions, renewal theorem

##  1 Introduction

Spectral heat content on an open set $D\subset{\mathbb{R}}^{d}$ measures the
total heat that remains on $D$ at time $t>0$ when the initial temperature is
one with Dirichlet boundary condition outside $D$. The spectral heat content
with respect to Brownian motions has been studied intensively, not only for
domains with smooth boundary (for example, see [5]) but also for certain
domains with fractal boundaries such as the $s$-adic von Koch snowflake (see
[3, 4, 6, 7]).

Recently there have been increasing interests in the spectral heat content for
more general Lévy processes (see [2, 8, 9, 13]). In [8] the authors studied
the spectral heat content for some class of Lévy processes on bounded open
sets in ${\mathbb{R}}$. In particular, it is proved in [8, Theorems 4.2 and
4.14] that when the underlying open set has infinitely many components (or
infinitely many non-adjacent components in the case of the Cauchy process) the
decay rate of the spectral heat content is strictly bigger than that of the
spectral heat content with respect to open sets with finitely many components.
Hence, a natural question is to determine the exact decay rate of the spectral
heat content for Lévy processes when there are infinitely many components in
$D\subset{\mathbb{R}}$.

Many Lévy processes can be realized as subordinate (time-changed) Brownian
motions, where the time change is given by an independent subordinator. Since
we need two operations to define the spectral heat content for subordinate
Brownian motions, time-change and killing, there are two objects that can be
called the spectral heat content for subordinate Brownian motions; One is
related to the killed subordinate Brownian motions (do time-change first then
kill the time-changed Brownian motions when they exit the domain under
consideration) and the other is related to subordinate killed Brownian motions
(kill Brownian motions when they exit the domain, then make time-change for
the killed Brownian motions). Even though the spectral heat content for killed
subordinate Brownian motions is a natural object to study as it covers a large
class of the spectral heat content for killed Lévy processes, the spectral
heat content for subordinate killed Brownian motions is also important as it
oftentimes gives useful information on the spectral heat content for the
killed subordinate Brownian motions. For example, in [14] the asymptotic
behavior of the spectral heat content for subordinate killed Brownian motions
with respect to stable subordinators provides crucial information on the
spectral heat content for killed stable processes.

In this paper, we study the spectral heat content for the subordinate killed
Brownian motions (see (2.10) below for the definition) when the underlying
subordinator is a stable subordinator
$S^{(\alpha/2)}=\\{S_{t}^{(\alpha/2)}\\}_{t\geq 0}$ whose Laplace transform is
given by

${\mathbb{E}}[e^{-\lambda
S_{t}^{(\alpha/2)}}]=e^{-t\lambda^{\alpha/2}},\quad\lambda>0,\
\alpha\in(0,2),$ (1.1)

and the underlying open sets have fractal boundaries which are determined by
similitudes in ${\mathbb{R}}^{d}$. Our main results answer the aforementioned
question and they show that the exact decay rates of the spectral heat content
depend on whether the sequence of logarithms of the coefficients of the
similitudes is arithmetic when $\alpha\in[d-\mathfrak{b},2)$, where
$\mathfrak{b}\in(d-1,d)$ is the interior Minkowski dimension of the boundary
of the underlying set (cf. (2.1) below for definition). It is noteworthy to
observe that when $\alpha\in(0,d-\mathfrak{b})$ the theorem is independent of
whether the sequence $\\{\ln(1/r_{j})\\}$ of the logarithms of the
coefficients of the similitudes is arithmetic or not.

The main technique for studying the asymptotic behavior of the spectral heat
content for open sets with fractal boundaries in this paper is the renewal
theorem in [12]. Two crucial properties that we need are the additivity
property of the spectral heat content under disjoint union (Lemma 2.5) and the
scaling property of the subordinate killed Brownian motions with respect to
stable subordinators (Lemma 2.7). However, in order to apply the renewal
theorem, one needs an exponential decay condition (2.9) and this is only valid
when $\alpha\in(d-\mathfrak{b},2)$. In order to establish the asymptotic
behavior of the spectral heat content as $t\to 0$ for the full range of
$\alpha\in(0,2)$, we will use the weak convergence of Lévy measure (see
Proposition 3.7) to establish the asymptotic behavior for the case when
$\alpha\in(0,d-\mathfrak{b})$. The result is given as Theorem 3.8. The
remaining case when $\alpha=d-\mathfrak{b}\in(0,1)$ is proved in Theorems 3.4
and 3.5 for the arithmetic and non-arithmetic cases, respectively. We observe
that there is an extra logarithm term $\ln(1/t)$ in the decay rate of the
spectral heat content when $\alpha=d-\mathfrak{b}$. This is due to the fact
that the heat loss $|G|-Q_{G}^{(2)}(u)$ for Brownian motions on the open set
$G$ with fractal boundaries, where
$|G|-Q_{G}^{(2)}(u)=\int_{G}{\mathbb{P}}_{x}(\tau_{G}^{(2)}\leq u)dx$ and
$\tau_{G}^{(2)}$ is the first exit time of the Brownian motions out of $G$, is
barely non-integrable with respect to the law of stable subordinator
$S_{t}^{(\alpha/2)}$. Note that this occurrence of extra logarithm term
$\ln(1/t)$ is observed for smooth open sets in [1, 13, 14], but it happens at
a different index $\alpha$. More specifically, for the spectral heat content
on smooth sets this phenomenon happens when $\alpha=1$, whereas for open sets
with fractal boundaries as in our case, this happens when
$\alpha=d-\mathfrak{b}$, which is strictly less than 1.

The organization of this paper is as follows. In Section 2 we set up
notations, define the class of open sets with fractal boundaries in (2.3), and
recall some facts which will be used for proving our main theorems. In
particular, we recall the renewal theorem in [12] and the result for the
spectral heat content for Brownian motions in Theorem 2.4. Section 3 is the
main part of this paper and here we study the spectral heat content for
subordinate killed Brownian motions. The main results are Theorems 3.2, 3.4,
3.5, and 3.8. This section is divided into three subsections for the cases
$\alpha\in(d-\mathfrak{b},2)$, $\alpha=d-\mathfrak{b}$, and
$\alpha\in(0,d-\mathfrak{b})$, respectively.

We use $c_{i}$ to denote constants whose values are unimportant and may change
from one appearance to another. The notations ${\mathbb{P}}_{x}$ and
${\mathbb{E}}_{x}$ mean probability and expectation of the underlying
processes started at $x\in{\mathbb{R}}^{d}$, and we use
${\mathbb{P}}={\mathbb{P}}_{0}$ and ${\mathbb{E}}={\mathbb{E}}_{0}$ to
simplify notations.

##  2 Preliminaries

In this section, we introduce notations and recall some facts that will be
used for proving the main theorems in Section 3.

###  2.1 Some geometric notions

We first recall some geometric notions and the definition of the class of open
sets with fractal boundaries from [12]. See also [10, 11] for more recent
developments.

For any bounded (open) set $G\subset{\mathbb{R}}^{d}$ with boundary $\partial
G$ and any $\varepsilon>0$, let

$G_{\varepsilon}^{\rm int}=\\{x\in G:{\rm dist}(x,\partial G)<\varepsilon\\}$

be the interior Minkowski sausage of radius $\varepsilon$ of the boundary
$\partial G$. We denote by $\mu(\varepsilon;G)$ the $d$-dimensional Lebesgue
measure of $G_{\varepsilon}^{\rm int}$. For any $s>0$, define

${\mathscr{M}}^{*}(s,\partial G)=\limsup_{\varepsilon\to
0}{\varepsilon}^{-(d-s)}\mu(\varepsilon;G)$

and

${\mathscr{M}}_{*}(s,\partial G)=\liminf_{\varepsilon\to
0}{\varepsilon}^{-(d-s)}\mu(\varepsilon;G).$

Following [12, Definition 1.2], the interior Minkowski dimension of $\partial
G$ (which is also called the Minkowski dimension of $\partial G$ relative to
$G$) is defined by

$\dim_{\rm M}^{\rm int}(\partial G)=\inf\\{s>0:\,{\mathscr{M}}^{*}(s,\partial
G)=0\\}=\sup\\{s>0:\,{\mathscr{M}}_{*}(s,\partial G)=\infty\\}.$ (2.1)

If, for $s=\dim_{\rm M}^{\rm int}(\partial G)$, we have
$0<{\mathscr{M}}^{*}(s,\partial G)={\mathscr{M}}_{*}(s,\partial G)<\infty$,
then, $\partial G$ is said to be Minkowski measurable relative to $G$.

###### Definition 2.1

A map $R:{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ is called a similitude with
coefficient $r>0$ if

$|Rx-Ry|=r|x-y|\quad\hbox{ for all }x,\,y\in{\mathbb{R}}^{d}.$

It is well known (cf. e.g., [12, p.191]) that any similitude is a composition
of a homothety with coefficient $r$, an orthogonal transform, and a
translation.

Now we define the class of open sets with fractal boundaries that we will
consider in this paper. Let $G_{0}\subset{\mathbb{R}}^{d}$ be a bounded open
set. When $d=1$ we assume that $G_{0}$ is a bounded open interval and when
$d\geq 2$ we assume that $G_{0}$ is a bounded $C^{1,1}$ open set. Let $R_{j}$
($1\leq j\leq N$) be similitudes with coefficients $r_{j}$, respectively.

For each $n\geq 1$, define
$\Upsilon_{n}=\\{{\boldsymbol{j}}=(j_{1},\ldots,j_{n}),\ 1\leq j_{i}\leq
N\\}$. We define the set $G$ by

$G=\bigg{(}\bigcup_{n=1}^{\infty}\bigcup_{{\boldsymbol{j}}\in\Upsilon_{n}}{\mathscr{R}}_{{\boldsymbol{j}}}G_{0}\bigg{)}\cup
G_{0},$ (2.2)

where, for every ${\boldsymbol{j}}=(j_{1},\ldots,j_{n})\in\Upsilon_{n}$,
${\mathscr{R}}_{{\boldsymbol{j}}}$ is the similitude defined by
${\mathscr{R}}_{{\boldsymbol{j}}}=R_{j_{1}}\circ\cdots\circ R_{j_{n}}$. It
follows from (2.2) that $G$ can be represented as

$G=\bigg{(}\bigcup_{j=1}^{N}R_{j}G\bigg{)}\cup G_{0}.$ (2.3)

We assume that all the sets $R_{j}G$, $1\leq j\leq N$, and $G_{0}$ in (2.3)
are pairwise disjoint. As in [12, Equation (1.7)] we also assume that
$\sum_{j=1}^{N}r_{j}^{d}<1<\sum_{j=1}^{N}r_{j}^{d-1}$. Since the expressions
in (2.3) are pairwise disjoint, we have

$|G|=\sum_{j=1}^{N}r_{j}^{d}|G|+|G_{0}|.$ (2.4)

Also, the condition $\sum_{j=1}^{N}r_{j}^{d}<1<\sum_{j=1}^{N}r_{j}^{d-1}$
ensures that $G$ has a finite volume and there exists a unique number
$\mathfrak{b}\in(d-1,d)$ such that

$\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}=1.$ (2.5)

It follows from [12, Theorem A] that the number $\mathfrak{b}$ is equal to the
interior Minkowski dimension of $\partial G$.

As an illustration, we consider the following examples. Let $R_{1}$ and
$R_{2}$ be two similitudes on ${\mathbb{R}}$ defined by

$R_{1}(x)=\frac{1}{3}x\quad\hbox{ and }\quad
R_{2}(x)=\frac{1}{3}x+\frac{2}{3}.$

We take $G_{0}=(\frac{1}{3},\frac{2}{3})$. Then, it is easy to observe that
the set defined in (2.2) is given by $G=(0,1)\setminus\mathfrak{C}$, where
$\mathfrak{C}$ is the standard ternary Cantor set and $G$ satisfies (2.3), and
$\mathfrak{b}=\log 2/\log 3$. Similarly, the open set $G\subset(0,\infty)^{2}$
with the Sierpínski gasket as its boundary can be obtained from $G_{0}$ being
the open triangle with vertices $(1/4,\sqrt{3}/4),\,(1/2,0)$ and
$(3/4,\,\sqrt{3}/4)$ (notice that the boundary $\partial G_{0}$ is not
$C^{1,1}$) and three similitudes on ${\mathbb{R}}^{2}$ defined by

$R_{1}(x)=\frac{1}{2}x,\quad
R_{2}(x)=\frac{1}{2}x+\Big{(}\frac{1}{4},\frac{\sqrt{3}}{4}\Big{)}\ \ \ \hbox{
and }\ \ R_{3}(x)=\frac{1}{2}x+\Big{(}\frac{1}{2},0\Big{)}.$

The (interior) Minkowski dimension of the boundary $\partial G$ is
$\mathfrak{b}=\log 3/\log 2$.

###  2.2 The renewal theorem

Now we state a version of the renewal theorem from [12]. Let
$f:{\mathbb{R}}\to{\mathbb{R}}$ be a map. For any $\gamma\in{\mathbb{R}}$
define

$L_{\gamma}f(z)=f(z-\gamma),$

and

$Lf(z)=\sum_{j=1}^{N}c_{j}L_{\gamma_{j}}f(z)=\sum_{j=1}^{N}c_{j}f(z-\gamma_{j}),$

where $c_{j}>0$, $\gamma_{j}$ are distinct points in ${\mathbb{R}}$, and
$\sum_{j=1}^{N}c_{j}=1$.

Consider the following renewal equation

$f=Lf+\phi.$ (2.6)

Intuitively, it is natural to expect that the solution of the renewal equation
is given by

$f(z)=\sum_{n=0}^{\infty}L^{n}\phi(z)=\phi(z)+\sum_{n=1}^{\infty}\sum_{c_{i_{1}},\cdots,c_{i_{n}}}c_{i_{1}}\cdots
c_{i_{n}}L_{\gamma_{i_{1}}}\cdots L_{\gamma_{i_{n}}}\phi(z).$ (2.7)

The following renewal theorem says it is indeed the case under certain
conditions. We say a set of finite real numbers
$\\{\gamma_{1},\cdots,\gamma_{N}\\}$ is arithmetic if
$\frac{\gamma_{i}}{\gamma_{j}}\in\mathbb{Q}$ for all indices. The maximal
number $\gamma$ such that $\frac{\gamma_{i}}{\gamma}\in\mathbb{Z}$ is called
the span of $\\{\gamma_{1},\cdots,\gamma_{N}\\}$. If the set is not
arithmetic, it is called non-arithmetic.

######  Theorem 2.2 (Renewal Theorem [12])

Suppose that a map $f:{\mathbb{R}}\to{\mathbb{R}}$ satisfies the renewal
equation (2.6) and it satisfies

$\lim_{z\to-\infty}f(z)=0,$ (2.8)

and

$|\phi(z)|\leq c_{1}e^{-c_{2}|z|},\quad z\in{\mathbb{R}},$ (2.9)

for some constants $c_{1},c_{2}>0$. Then, the solution of the renewal equation
(2.6) is given by (2.7). Furthermore, if $\\{\gamma_{j}\\}$ is non-arithmetic,
then

$f(z)=\frac{1}{\sum_{j=1}^{N}c_{j}\gamma_{j}}\int_{-\infty}^{\infty}\phi(x)dx+o(1),\text{
as }z\to\infty.$

If $\\{\gamma_{j}\\}$ is arithmetic with span $\gamma$, then

$f(z)=\frac{\gamma}{\sum_{j=1}^{N}c_{j}\gamma_{j}}\sum_{k=-\infty}^{\infty}\phi(z-k\gamma)+o(1),\text{
as }z\to\infty.$

###  2.3 The spectral heat content of subordinate killed Brownian motions

For an open set $G\subset{\mathbb{R}}^{d}$ we define the spectral heat content
for Brownian motion $W=\\{W_{t}\\}_{t\geq 0}$ on $G$ as

$Q_{G}^{(2)}(t)=\int_{G}{\mathbb{P}}_{x}\big{(}\tau_{G}^{(2)}>t\big{)}dx,\quad\tau_{G}^{(2)}=\inf\\{t>0:W_{t}\notin
G\\}.$

We record the following lemma from [12, Lemma 4.4]. Note that there was a typo
there and $r^{2}$ should be written as $r^{d}$.

######  Lemma 2.3 (Lemma 4.4 [12])

Let $G$ be an open set in ${\mathbb{R}}^{d}$ and $R$ be a similitude with
coefficient $r$. Then,

$Q_{RG}^{(2)}(t)=r^{d}Q_{G}^{(2)}(t/r^{2}).$

We recall the following theorem for the spectral heat content for Brownian
motion from [12, Theorem D].

######  Theorem 2.4 (Theorem D [12])

Let $G$ be a set defined as in (2.2) with $R_{j}$ being similitude with
coefficient $r_{j}$, and $G_{0}$ is either a bounded open interval when $d=1$,
or a bounded $C^{1,1}$ open set when $d\geq 2$.

  1. 1.

If $\\{\ln(\frac{1}{r_{j}})\\}_{j=1}^{N}$ is non-arithmetic, then

$Q_{G}^{(2)}(t)=|G|-Ct^{\frac{d-\mathfrak{b}}{2}}+o(t^{\frac{d-\mathfrak{b}}{2}}),\quad\text{
where
}C=\frac{\int_{0}^{\infty}\left(|G_{0}|-Q_{G_{0}}^{(2)}(u)\right)u^{-(1+\frac{d-\mathfrak{b}}{2})}du}{\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\ln(\frac{1}{r_{j}^{2}})}.$

  2. 2.

If $\\{\ln(\frac{1}{r_{j}})\\}_{j=1}^{N}$ is arithmetic with span $\rho$, then

$Q_{G}^{(2)}(t)=|G|-s(-\ln
t)t^{\frac{d-\mathfrak{b}}{2}}+o(t^{\frac{d-\mathfrak{b}}{2}}),\quad\text{ as
}t\to 0,$

where the function $s(\cdot)$ is defined by

$\displaystyle s(z)$ $\displaystyle:=$
$\displaystyle\frac{2\rho}{\sum_{j=1}^{N}(r_{j})^{\mathfrak{b}}\ln(\frac{1}{r_{j}^{2}})}\sum_{n=-\infty}^{\infty}\left(|G_{0}|-Q_{G_{0}}^{(2)}\big{(}e^{-(z-2n\rho)}\big{)}\right)e^{\frac{d-\mathfrak{b}}{2}(z-2n\rho)}.$

Now we introduce the spectral heat content for subordinate killed Brownian
motions. Let $W=\\{W_{t}\\}_{t\geq 0}$ be Brownian motion in
${\mathbb{R}}^{d}$ and let $S^{(\alpha/2)}=\\{S_{t}^{(\alpha/2)}\\}_{t\geq 0}$
be an $(\alpha/{2})$-stable subordinator with Laplace transform given by
(1.1), which is independent of $W$. Let $D$ be any open set in
${\mathbb{R}}^{d}$. Then, the spectral heat content
$\tilde{Q}_{D}^{(\alpha)}(t)$ for subordinate killed Brownian motions with
respect to stable subordinator $S^{(\alpha/2)}$ on $D$ is defined by

$\tilde{Q}_{D}^{(\alpha)}(t)=\int_{D}{\mathbb{P}}_{x}\Big{(}\tau_{D}^{(2)}>S_{t}^{(\alpha/2)}\Big{)}dx,$
(2.10)

where $\tau_{D}^{(2)}=\inf\\{t>0:W_{t}\notin D\\}$.

We will need the following important properties for
$\tilde{Q}_{D}^{(\alpha)}(t)$; One is the additivity under disjoint union and
another is the scaling property.

###### Lemma 2.5

Let $D_{1},D_{2}$ be open sets in ${\mathbb{R}}^{d}$ with $D_{1}\cap
D_{2}=\emptyset$. Then,

$\tilde{Q}^{(\alpha)}_{D_{1}\cup
D_{2}}(t)=\tilde{Q}^{(\alpha)}_{D_{1}}(t)+\tilde{Q}^{(\alpha)}_{D_{2}}(t).$

Proof. Note that

$\begin{split}\tilde{Q}^{(\alpha)}_{D_{1}\cup D_{2}}(t)&=\int_{D_{1}\cup
D_{2}}{\mathbb{P}}_{x}\Big{(}\tau^{(2)}_{D_{1}\cup
D_{2}}>S_{t}^{(\alpha/2)}\Big{)}dx\\\
&=\int_{D_{1}}{\mathbb{P}}_{x}\Big{(}\tau^{(2)}_{D_{1}\cup
D_{2}}>S_{t}^{(\alpha/2)}\Big{)}dx+\int_{D_{2}}{\mathbb{P}}_{x}\Big{(}\tau^{(2)}_{D_{1}\cup
D_{2}}>S_{t}^{(\alpha/2)}\Big{)}dx.\end{split}$

Note that under ${\mathbb{P}}_{x}$ with $x\in D_{1}$ we have

$\tau^{(2)}_{D_{1}\cup D_{2}}=\inf\\{t>0:W_{t}\notin D_{1}\cup
D_{2}\\}=\inf\\{t>0:W_{t}\notin D_{1}\\}=\tau_{D_{1}}^{(2)}.$

Hence, we have

$\int_{D_{1}}{\mathbb{P}}_{x}\Big{(}\tau^{(2)}_{D_{1}\cup
D_{2}}>S_{t}^{(\alpha/2)}\Big{)}dx=\int_{D_{1}}{\mathbb{P}}_{x}\Big{(}\tau^{(2)}_{D_{1}}>S_{t}^{(\alpha/2)}\Big{)}dx=\tilde{Q}_{D_{1}}^{(\alpha)}(t).$

It can be proved in the same way that the integral on $D_{2}$ gives
$\tilde{Q}^{(\alpha)}_{D_{2}}(t).$ $\Box$

###### Remark 2.6

Let
$Q_{D}^{(\alpha)}(t):=\int_{D}{\mathbb{P}}_{x}\Big{(}\tau_{D}^{(\alpha)}>t\Big{)}dx$
be the spectral heat content for killed stable processes, where
$\tau_{D}^{(\alpha)}$ is the first exit time of the $\alpha$-stable process
$W_{S^{(\alpha/2)}}=\\{W_{S_{t}^{(\alpha/2)}}\\}_{t\geq 0}$. This is the
spectral heat content related to the killed subordinate Brownian motions by
stable subordinators. Note that for disjoint sets $D_{1}$ and $D_{2}$ we have

$\displaystyle Q^{(\alpha)}_{D_{1}\cup D_{2}}(t)$
$\displaystyle=\int_{D_{1}\cup
D_{2}}{\mathbb{P}}\Big{(}\tau^{(\alpha)}_{D_{1}\cup
D_{2}}>t\Big{)}dx=\int_{D_{1}}{\mathbb{P}}\Big{(}\tau^{(\alpha)}_{D_{1}\cup
D_{2}}>t\Big{)}dx+\int_{D_{2}}{\mathbb{P}}\Big{(}\tau^{(\alpha)}_{D_{1}\cup
D_{2}}>t\Big{)}dx$
$\displaystyle\geq\int_{D_{1}}{\mathbb{P}}\Big{(}\tau^{(\alpha)}_{D_{1}}>t\Big{)}dx+\int_{D_{2}}{\mathbb{P}}\Big{(}\tau^{(\alpha)}_{D_{2}}>t\Big{)}dx=Q^{(\alpha)}_{D_{1}}(t)+Q^{(\alpha)}_{D_{2}}(t).$

Furthermore, the inequality can be strict as $\tau^{(\alpha)}_{D_{1}\cup
D_{2}}\neq\tau^{(\alpha)}_{D_{1}}$ when the process starts at $x\in D_{1}$
because the process starting at $x\in D_{1}$ can jump into $D_{2}$ without
visiting the complement of $D_{1}\cup D_{2}$. Hence, the spectral heat content
for killed subordinate Brownian motions does not satisfy the additivity
property under disjoint union.

###### Lemma 2.7

Let $R$ be a similitude with coefficient $r$ and $G$ is any open set in
${\mathbb{R}}^{d}$. Then, we have

$\tilde{Q}^{(\alpha)}_{RG}(t)=r^{d}\tilde{Q}^{(\alpha)}_{G}(t/r^{\alpha}),\quad
t>0.$

Proof. By the scaling property and rotational invariance of Brownian motions,
we observe that

$\tau_{RG}^{(2)}\text{ under }{\mathbb{P}}_{Rx}\text{ is equal in distribution
to }r^{2}\tau_{G}^{(2)}\text{ under }{\mathbb{P}}_{x}.$

By the change of variable $x=Ry$ and the scaling property of
$S_{t}^{(\alpha/2)}$ we have

$\begin{split}\tilde{Q}^{(\alpha)}_{RG}(t)&=\int_{RG}{\mathbb{P}}_{x}\Big{(}\tau_{RG}^{(2)}>S_{t}^{(\alpha/2)}\Big{)}dx=\int_{G}{\mathbb{P}}_{Ry}\Big{(}\tau_{RG}^{(2)}>S_{t}^{(\alpha/2)}\Big{)}r^{d}dy\\\
&=\int_{G}{\mathbb{P}}_{y}\Big{(}r^{2}\tau_{G}^{(2)}>S_{t}^{(\alpha/2)}\Big{)}r^{d}dy=\int_{G}{\mathbb{P}}_{y}\Big{(}\tau_{G}^{(2)}>r^{-2}S_{t}^{(\alpha/2)}\Big{)}r^{d}dy\\\
&=\int_{G}{\mathbb{P}}_{y}\Big{(}\tau_{G}^{(2)}>S_{tr^{-\alpha}}^{(\alpha/2)}\Big{)}r^{d}dy=r^{d}\tilde{Q}_{G}^{(\alpha)}(t/r^{\alpha}).\end{split}$

This proves the lemma. $\Box$

##  3 Asymptotic behavior of the spectral heat content

###  3.1 The case of $\alpha\in(d-\mathfrak{b},2)$

Analogous to [12, Theorem D], we will prove that the spectral heat content
$\tilde{Q}^{(\alpha)}_{G}(t)$ as defined in (2.10) has the form

$\tilde{Q}^{(\alpha)}_{G}(t)=|G|-f(-\ln
t)t^{\frac{d-\mathfrak{b}}{\alpha}}+o(t^{\frac{d-\mathfrak{b}}{\alpha}}),$

when $\alpha\in(d-\mathfrak{b},2)$, where $\mathfrak{b}$ is the constant in
(2.5). We start with the following lemma.

###### Lemma 3.1

Assume that $\alpha\in(d-\mathfrak{b},2)$. Suppose that $G_{0}$ is an open
interval when $d=1$ or a bounded $C^{1,1}$ open set when $d\geq 2$. Define

$\psi(z)=\left(|G_{0}|-\tilde{Q}_{G_{0}}^{(\alpha)}(e^{-z})\right)e^{\frac{d-\mathfrak{b}}{\alpha}z}.$

Then, there exists a constant $c=c(\alpha,\mathfrak{b},d)>0$ such that

$|\psi(z)|\leq ce^{-c|z|}\text{ for all }z\in{\mathbb{R}}.$

Proof. We define the function

$\phi(t)=\left(|G_{0}|-\tilde{Q}^{(\alpha)}_{G_{0}}(t)\right)t^{-\frac{d-\mathfrak{b}}{\alpha}},\
\ \ t>0,$

so that $\psi(z)=\phi(e^{-z})$.

The case when $z\to-\infty$ is easy since we have

$|\phi(t)|\leq|G_{0}|t^{-\frac{d-\mathfrak{b}}{\alpha}}$

and this gives

$\psi(z)=\phi(e^{-z})\leq|G_{0}|e^{\frac{d-\mathfrak{b}}{\alpha}z}=|G_{0}|e^{-\frac{d-\mathfrak{b}}{\alpha}|z|}$

for all $z\leq 0$ and $\alpha\in(0,2)$.

Now we handle the case when $z\to\infty$, or $t=e^{-z}\to 0$. First, assume
that $\alpha\in(1,2)$. Since $G_{0}$ is an interval when $d=1$ or a bounded
$C^{1,1}$ open set when $d\geq 2$, it follows from [13, Theorem 1.1] that
there exists a constant $c_{1}>0$ such that

$|G_{0}|-\tilde{Q}^{(\alpha)}_{G_{0}}(t)\leq c_{1}t^{1/\alpha}$

for all $0<t\leq 1$. Hence,

$\phi(t)\leq c_{1}t^{\frac{(\mathfrak{b}+1)-d}{\alpha}}\text{ for }0<t\leq 1.$

Since $\mathfrak{b}\in(d-1,d)$ we note that
$\frac{(\mathfrak{b}+1)-d}{\alpha}>0$. Let $z=-\ln t$ and we conclude that

$\phi(t)=\phi(e^{-z})=\psi(z)\leq c_{2}e^{-c_{3}|z|},\quad z\in{\mathbb{R}},$

where $c_{2}=\max(c_{1},|D_{0}|)$ and
$c_{3}=\min(\frac{\mathfrak{b}+1-d}{\alpha},\frac{d-\mathfrak{b}}{\alpha})>0$.

Second, when $\alpha=1$ we have from [13, Theorem 1.1]

$|G_{0}|-\tilde{Q}^{(\alpha)}_{G_{0}}(t)\leq c_{2}t\ln(1/t)$

for all $0<t\leq 1$. Hence,

$\phi(t)\leq
c_{2}t^{1-(d-\mathfrak{b})}\ln(1/t)=c_{2}t^{(\mathfrak{b}+1)-d}\ln(1/t)\quad\text{
for }0<t\leq 1,$

and this implies

$\psi(z)=\phi(e^{-z})\leq c_{2}ze^{-z(\mathfrak{b}+1-d)}\text{ for }z\geq 0.$

Since $\mathfrak{b}+1-d>0$ there exists $c_{4}$ and $\eta>0$ such that

$\psi(z)\leq c_{4}e^{-\eta z}\ \ \text{ for all }z\geq 0.$

Finally, we handle the case when $\alpha\in(d-\mathfrak{b},1)$. From [13,
Theorem 1.1] we have

$|G_{0}|-\tilde{Q}^{(\alpha)}_{G_{0}}(t)\leq c_{5}t$

for all $0<t\leq 1$, and this implies $\phi(t)\leq
c_{5}t^{\frac{\alpha-(d-\mathfrak{b})}{\alpha}}\text{ for }t\leq 1$, which in
turn implies $\psi(z)\leq
c_{5}e^{-\frac{\alpha-(d-\mathfrak{b})}{\alpha}z}\text{ for }z\geq 0$. $\Box$

Here is the main theorem for the case of $\alpha\in(d-\mathfrak{b},2)$.

###### Theorem 3.2

Let $\alpha\in(d-\mathfrak{b},2)$, where $\mathfrak{b}$ is the constant in
(2.5) and $G$ is a set given as (2.2) with $G_{0}$ being an open interval when
$d=1$ or a bounded $C^{1,1}$ open set when $d\geq 2$. If
$\\{\ln\frac{1}{r_{j}}\\}_{j=1}^{N}$ is non-arithmetic, then we have

$\tilde{Q}^{(\alpha)}_{G}(t)=|G|-C_{1}t^{\frac{d-\mathfrak{b}}{\alpha}}+o(t^{\frac{d-\mathfrak{b}}{\alpha}})\text{
as }t\to 0,$

where

$\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle\frac{\int_{-\infty}^{\infty}\left(|G_{0}|-\tilde{Q}_{G_{0}}^{(\alpha)}(e^{-z})\right)e^{\frac{(d-\mathfrak{b})z}{\alpha}}dz}{\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\ln(1/r_{j}^{\alpha})}=\frac{\int_{0}^{\infty}\left(|G_{0}|-\tilde{Q}_{G_{0}}^{(\alpha)}(t)\right)t^{-1-\frac{d-\mathfrak{b}}{\alpha}}dt}{\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\ln(1/r_{j}^{\alpha})}.$

If $\\{\ln\frac{1}{r_{j}}\\}_{j=1}^{N}$ is arithmetic with span $\rho$, then
we have

$\tilde{Q}^{(\alpha)}_{G}(t)=|G|-f(-\ln
t)t^{\frac{d-\mathfrak{b}}{\alpha}}+o(t^{\frac{d-\mathfrak{b}}{\alpha}})\text{
as }t\to 0,$

where

$f(z)=\frac{\alpha\rho}{\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\ln(1/r_{j}^{\alpha})}\sum_{n=-\infty}^{\infty}\left(|G_{0}|-\tilde{Q}_{G_{0}}^{(\alpha)}(e^{-(z-\alpha
n\rho)})\right)e^{\frac{d-\mathfrak{b}}{\alpha}(z-\alpha n\rho)}.$

Proof. We set

$\tilde{Q}_{G}^{(\alpha)}(t)=|G|-f(-\ln t)t^{\frac{d-\mathfrak{b}}{\alpha}},$
(3.1)

and will show that $f(\cdot)$ satisfies the renewal equation and conditions
for the renewal theorem.

Since $G=\bigcup_{j=1}^{N}R_{j}G\cup G_{0}$ and all expressions are disjoint,
it follows from Lemmas 2.5 and 2.7

$\displaystyle\tilde{Q}_{G}^{(\alpha)}(t)=\tilde{Q}_{\bigcup_{j=1}^{N}R_{j}G\cup
G_{0}}^{(\alpha)}(t)=\sum_{j=1}^{N}\tilde{Q}_{R_{j}G}^{(\alpha)}(t)+\tilde{Q}_{G_{0}}^{(\alpha)}(t)=\sum_{j=1}^{N}r_{j}^{d}\tilde{Q}_{G}^{(\alpha)}(t/r_{j}^{\alpha})+\tilde{Q}^{(\alpha)}_{G_{0}}(t).$

Using (3.1) we have

$\begin{split}&|G|-f(-\ln
t)t^{\frac{d-b}{\alpha}}=\sum_{j=1}^{N}r_{j}^{d}\left(|G|-f(-\ln\frac{t}{(r_{j})^{\alpha}})(\frac{t}{(r_{j})^{\alpha}})^{\frac{d-\mathfrak{b}}{\alpha}}\right)+\tilde{Q}^{(\alpha)}_{G_{0}}(t)\\\
&=\sum_{j=1}^{N}r_{j}^{d}|G|-\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\cdot
f\left(-\ln
t-\ln(\frac{1}{(r_{j})^{\alpha}})\right)t^{\frac{d-\mathfrak{b}}{\alpha}}+\left(|G_{0}|-(|G_{0}|-\tilde{Q}^{(\alpha)}_{G_{0}}(t))\right).\end{split}$

Using (2.4) we conclude that

$f(-\ln t)=\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\cdot f\left(-\ln
t-\ln(\frac{1}{(r_{j})^{\alpha}})\right)+\phi(t),\quad\text{where
}\phi(t)=\left(|G_{0}|-\tilde{Q}^{(\alpha)}_{G_{0}}(t)\right)t^{-\frac{d-\mathfrak{b}}{\alpha}}.$

By changing the variable $z=-\ln t$ we have

$f(z)=\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\cdot
f\left(z-\ln(\frac{1}{(r_{j})^{\alpha}})\right)+\phi(e^{-z}).$

Note that from (3.1) we have

$\lim_{z\to-\infty}f(z)=\lim_{t\to\infty}f(-\ln
t)=\lim_{t\to\infty}\left(|G|-\tilde{Q}^{(\alpha)}_{G}(t)\right)t^{-\frac{d-\mathfrak{b}}{\alpha}}=0,$

and this shows that the condition (2.8) holds. It follows from Lemma 3.1 that
for any $\alpha\in(d-\mathfrak{b},2)$ there exist two constants
$c_{1},c_{2}>0$ such that

$\psi(z)=\phi(e^{-z})\leq c_{1}e^{-c_{2}|z|}\ \ \text{ for all
}z\in{\mathbb{R}},$

and the condition (2.9) holds. Now the conclusions of the theorem follow
immediately from the Renewal Theorem 2.2. $\Box$

###  3.2 The case of $\alpha=d-\mathfrak{b}$

In this subsection, we study the case when $\alpha=d-\mathfrak{b}\in(0,1)$. We
need a simple lemma which is similar to [13, Lemma 3.2]. The proof is
essentially the same with obvious modifications and will be omitted.

###### Lemma 3.3

For any $\delta>0$ and $\alpha\in(0,2)$, we have

$\lim_{t\to
0}\frac{{\mathbb{E}}\left[\big{(}S_{1}^{(\alpha/2)}\big{)}^{\alpha/2},0<S_{1}^{(\alpha/2)}<\delta
t^{-2/\alpha}\right]}{\ln(1/t)}=\frac{1}{\Gamma(1-\frac{\alpha}{2})}.$

###### Theorem 3.4

Let $\alpha=d-\mathfrak{b}\in(0,1)$, where $\mathfrak{b}$ is the constant in
(2.5) and $G$ is a set given as (2.2) with $G_{0}$ being an open interval when
$d=1$ or a bounded $C^{1,1}$ open set when $d\geq 2$. Assume that
$\\{\ln(1/r_{j})\\}_{j=1}^{N}$ is arithmetic with span $\rho$. Define
$A=\sup_{z\in{\mathbb{R}}}s(z)$ and $B=\inf_{z\in{\mathbb{R}}}s(z)$, where
$s(z)$ is from Theorem 2.4.

  1. (1)

Let
$g(t):=\int_{0}^{t^{-2/\alpha}}s(-\ln(t^{2/\alpha}v))v^{\alpha/2}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)$. Then, we have

$|G|-\tilde{Q}_{G}^{(\alpha)}(t)=tg(t)+o(t\ln(1/t)).$ (3.2)

  2. (2)

We have

$\limsup_{t\to
0}\frac{g(t)}{\ln(1/t)}=\frac{A}{\Gamma(1-\frac{\alpha}{2})}\quad\text{ and
}\quad\liminf_{t\to
0}\frac{g(t)}{\ln(1/t)}=\frac{B}{\Gamma(1-\frac{\alpha}{2})}.$ (3.3)

Proof. Note that by the scaling property of $S_{t}^{(\alpha/2)}$ we have

$\begin{split}&|G|-\tilde{Q}_{G}^{(\alpha)}(t)=\int_{0}^{\infty}\left(|G|-Q_{G}^{(2)}(u)\right){\mathbb{P}}(S_{t}^{(\alpha/2)}\in
du)=\int_{0}^{\infty}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)\\\
&=\int_{0}^{t^{-2/\alpha}}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)+\int_{t^{-2/\alpha}}^{\infty}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)\\\
&=\int_{0}^{t^{-2/\alpha}}\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)+\int_{t^{-2/\alpha}}^{\infty}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv).\end{split}$

Hence, we have

$\begin{split}|G|-\tilde{Q}_{G}^{(\alpha)}(t)-tg(t)&=t\int_{0}^{t^{-2/\alpha}}\bigg{(}\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}-s(-\ln(t^{2/\alpha}v))\bigg{)}v^{\frac{d-\mathfrak{b}}{2}}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)\\\
&\qquad\qquad+\int_{t^{-2/\alpha}}^{\infty}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
du).\end{split}$

It follows from [13, Equation (2.8)] we have

$\int_{t^{-2/\alpha}}^{\infty}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
du)\leq
c\int_{t^{-2/\alpha}}^{\infty}|G|u^{-1-\frac{\alpha}{2}}du=o(t\ln(1/t)).$
(3.4)

In this case, by applying Theorem 2.4 we have

$\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}=s(-\ln(t^{2/\alpha}v))+o(1)\quad\text{
as }t\to 0.$ (3.5)

From (3.5) for any $\varepsilon>0$ there exists $t_{0}(\varepsilon)$ such that

$\left|\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}-s(-\ln(t^{2/\alpha}v))\right|<\varepsilon$

for all $t\leq t_{0}$. Hence it follows from Lemma 3.3

$\begin{split}&\limsup_{t\to
0}\frac{t\int_{0}^{t^{-2/\alpha}}\Big{(}\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}-s(-\ln(t^{2/\alpha}v))\Big{)}v^{\frac{d-\mathfrak{b}}{2}}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)}{t\ln(1/t)}\\\ &\leq\varepsilon\limsup_{t\to
0}\frac{\int_{0}^{t^{-2/\alpha}}v^{\alpha/2}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)}{\ln(1/t)}\leq\frac{\varepsilon}{\Gamma(1-\frac{\alpha}{2})}.\end{split}$

This establishes (3.2).

For (3.3), note that

$\begin{split}|G|-\tilde{Q}_{G}^{(\alpha)}(t)&=\int_{0}^{t^{-2/\alpha}}\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)\\\
&\qquad\qquad+\int_{t^{-2/\alpha}}^{\infty}\left(|G|-Q_{G}^{(2)}(t^{2/\alpha}v)\right){\mathbb{P}}(S_{1}^{(\alpha/2)}\in
du).\end{split}$

As in (3.4) the second expression above is $o(t\ln(1/t))$ as $t\to 0$. For any
$\varepsilon>0$ it follows from (3.5) we have
$\frac{|G|-Q_{G}^{(2)}(t^{2/\alpha}v)}{(t^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}<A+\varepsilon$
for all sufficiently small $t$. This fact, together with (3.2) and Lemma 3.3
gives

$\limsup_{t\to 0}\frac{g(t)}{\ln(1/t)}=\limsup_{t\to
0}\frac{|G|-\tilde{Q}_{G}^{(\alpha)}(t)}{t\ln(1/t)}\leq\frac{A+\varepsilon}{\Gamma(1-\frac{d-\mathfrak{b}}{2})}.$

Since $\varepsilon>0$ is arbitrary, we conclude that

$\limsup_{t\to 0}\frac{g(t)}{\ln(1/t)}=\limsup_{t\to
0}\frac{|G|-\tilde{Q}_{G}^{(\alpha)}(t)}{t\ln(1/t)}\leq\frac{A}{\Gamma(1-\frac{d-\mathfrak{b}}{2})}.$
(3.6)

For the lower bound, it follows from Theorem 2.4 that for any $\varepsilon>0$
there exists a sequence $t_{n}\to 0$ such that

$\frac{|G|-Q_{G}^{(2)}(t_{n}^{2/\alpha}v)}{(t_{n}^{2/\alpha}v)^{\frac{d-\mathfrak{b}}{2}}}\geq
A-\varepsilon.$

Hence, from Lemma 3.3, (3.4),and (3.2) we have

$\begin{split}\limsup_{n\to\infty}\frac{g(t_{n})}{\ln(1/t_{n})}&=\limsup_{n\to\infty}\frac{|G|-\tilde{Q}_{G}^{(\alpha)}(t_{n})}{t_{n}\ln(1/t_{n})}\\\
&\geq(A-\varepsilon)\limsup_{n\to\infty}\int_{0}^{t_{n}^{-2/\alpha}}v^{\frac{\alpha}{2}}{\mathbb{P}}(S_{1}^{(\alpha/2)}\in
dv)\geq\frac{A-\varepsilon}{\Gamma(1-\frac{d-\mathfrak{b}}{2})}.\end{split}$

Since $\varepsilon>0$ is arbitrary, we have

$\limsup_{t\to
0}\frac{g(t)}{\ln(1/t)}\geq\frac{A}{\Gamma(1-\frac{d-\mathfrak{b}}{2})}.$
(3.7)

Hence, the $\limsup$ version of (3.3) follows from (3.6) and (3.7), and the
$\liminf$ version can be proved in the same way. $\Box$

Here is the result for the non-arithmetic case. The proof is very similar to
the proof of Theorem 3.4, hence it will be omitted.

###### Theorem 3.5

Let $\alpha=d-\mathfrak{b}\in(0,1)$, where $\mathfrak{b}$ is the constant in
(2.5) and $G$ is a set given as (2.2) with $G_{0}$ being an open interval when
$d=1$ or a bounded $C^{1,1}$ open set when $d\geq 2$. Assume that
$\\{\ln(1/r_{j}\\}_{j=1}^{N}$ is non-arithmetic. Then, we have

$\lim_{t\to
0}\frac{|G|-\tilde{Q}_{G}^{(d-\mathfrak{b})}(t)}{t\ln(1/t)}=\frac{1}{2\Gamma(1-\frac{d-\mathfrak{b}}{2})\sum_{j=1}^{N}r_{j}^{\mathfrak{b}}\ln(1/r_{j})}\int_{0}^{\infty}\left(|G_{0}|-Q_{G_{0}}^{(2)}(u)\right)u^{-1-\frac{d-\mathfrak{b}}{2}}du.$

###  3.3 The case of $\alpha\in(0,d-\mathfrak{b})$

Now we handle the case when $\alpha\in(0,d-\mathfrak{b})$. We need the
following simple lemma for the continuity of the map $t\to\tilde{Q}_{D}(t)$
which is proved in [9, Lemma 3.11].

###### Lemma 3.6

For any open set $D$ with $|D|<\infty$, the map
$t\to\tilde{Q}^{(\alpha)}_{D}(t)$ is continuous.

The following proposition is proved in [9, Proposition 3.12].

###### Proposition 3.7

Let $f$ be a bounded continuous function on $(0,\infty)$ such that
$\displaystyle\lim_{x\downarrow 0}\frac{f(x)}{x^{\gamma}}$ exists as a finite
number for some constant $\gamma>\frac{\alpha}{2}$. Then, we have

$\lim_{t\downarrow
0}\int_{0}^{\infty}f(u)\frac{{\mathbb{P}}(S_{t}^{(\alpha/2)}\in
du)}{t}=\frac{\alpha}{2\Gamma(1-\frac{\alpha}{2})}\int_{0}^{\infty}f(u)\,u^{-1-\frac{\alpha}{2}}du.$

###### Theorem 3.8

Let $\alpha\in(0,d-\mathfrak{b})$, where $\mathfrak{b}$ is the constant in
(2.5) and $G$ is a set given as (2.2) with $G_{0}$ being an open interval when
$d=1$ or a bounded $C^{1,1}$ open set when $d\geq 2$. Then, we have

$\lim_{t\to
0}\frac{|G|-\tilde{Q}_{G}^{(\alpha)}(t)}{t}=\frac{\alpha}{2\Gamma(1-\frac{\alpha}{2})}\int_{0}^{\infty}\big{(}|G|-Q_{G}^{(2)}(u)\big{)}\,u^{-1-\frac{\alpha}{2}}du.$

Proof. Note that we have

$|G|-\tilde{Q}_{G}^{(\alpha)}(t)=\int_{0}^{\infty}\big{(}|G|-Q_{G}^{(2)}(u)\big{)}{\mathbb{P}}(S_{t}^{(\alpha/2)}\in
du).$

It follows from Theorem 2.4 that there exists constants $c_{1}$ such that

$|G|-Q_{G}^{(2)}(u)\leq c_{1}u^{\frac{d-\mathfrak{b}}{2}}\text{ for }u\leq 1.$

Since $\alpha\in(0,d-\mathfrak{b})$, we can take
$\gamma\in(\frac{\alpha}{2},\frac{d-\mathfrak{b}}{2})$ and this implies

$\lim_{u\to 0}\frac{|G|-Q_{G}^{(2)}(u)}{u^{\gamma}}=0.$

Now the conclusion of the theorem follows immediately from Proposition 3.7.
$\Box$

Acknowledgements The research of Y. Xiao is partially supported by NSF grant
DMS-1855185.

## References

  * [1] L. Acuña Valverde. On the one dimensional spectral heat content for stable processes. J. Math. Anal. Appl., 441 (2016), 11–24.
  * [2] L. Acuña Valverde. Heat content for stable processes in domains of ${\mathbb{R}}^{d}$. J. Geom. Anal., 27 (2017), 492–524.
  * [3] M. van den Berg. Heat content and Brownian motion for some regions with a fractal boundary. Probab. Th. Rel. Fields, 100 (1994), no. 4, 439–456.
  * [4] M. van den Berg. Heat equation on the arithmetic von Koch snowflake. Probab. Th. Rel. Fields, 118 (2000), no. 1, 17–36.
  * [5] M. van den Berg and E. B. Davies. Heat flow out of regions in ${\mathbb{R}}^{m}$. Math. Z., 202 (1989), 463–482.
  * [6] M. van den Berg and P. B. Gilkey. A comparison estimate for the heat equation with an application to the heat content of the $s$-adic von Koch snowflake. Bull. London Math. Soc., 30 (1998), no. 4, 404–412.
  * [7] J. Fleckinger, M. Levitin, and D. Vassiliev. Heat equation on the triadic von Koch snowflake: asymptotic and numerical analysis. Proc. London Math. Soc. (3), 71 (1995), no. 2, 372–396.
  * [8] T. Grzywny, H. Park, and R. Song. Spectral heat content for Lévy processes. Math. Nachr. 292 (2019), 805–825.
  * [9] K. Kobayashi and H. Park. Spectral heat contents for time-changed killed Brownian motions. Preprint. arxiv.org/abs/2007.05776.
  * [10] M. L. Lapidus and M. van Frankenhuijsen, Fractal Geometry, Complex Dimensions and Zeta Functions. Geometry and spectra of fractal strings. Second edition. Springer, New York, 2013.
  * [11] M. L. Lapidus, G. Radunović, and D. Z̆ubrinić, Fractal Zeta Functions and Fractal Drums. Higher-dimensional Theory of Complex Dimensions. Springer, Cham, 2017.
  * [12] M. Levitean and D. Vassiliev. Spectral asymptotics, renewal theorem, and the Berry conjecture for a class of fractals. Proc. London Math. Soc. 72 (3) (1996), no. 1, 188–214.
  * [13] H. Park and R. Song. Small time asymptotics of spectral heat contents for subordinate killed Brownian motions related to isotropic $\alpha$-stable processes. Bul. London Math. Soc., 51, Issue 2 (2019), 371–384.
  * [14] H. Park and R. Song. Spectral heat content for $\alpha$-stable processes in $C^{1,1}$ open sets. Preprint. arxiv.org/abs/2007.02815.

Hyunchul Park

Department of Mathematics, State University of New York at New Paltz, NY
12561, USA

E-mail<EMAIL_ADDRESS>

Yimin Xiao

Department of Statistics and Probability, Michigan State University, East
Lansing, MI 48824, USA

E-mail<EMAIL_ADDRESS>