text:
\section{INTRODUCTION}
The extension of the classical {\it arrival time} concept to quantum
mechanics is problematic because of the absence of trajectories in the
standard interpretation of the quantum formalism.
The question is important since detection of particles in time-of-flight
and coincidence experiments are
quite common, and quantum mechanics should
be able to predict and interpret the statistics of the arrivals.
A recent review discusses the main theoretical approaches
and open questions \cite{MSP}.

The customary way to specify a dynamical variable in quantum
mechanics is to define the corresponding operator,
but Pauli showed that a
self adjoint time operator
conjugate to the Hamiltonian $\hat{H}$ cannot be constructed if the spectrum
of $\hat{H}$ is bounded from below \cite{Pauli}, see also
\cite{All,Gia,DM}.
There are however several ways to circumvent this objection and quantize
the classical arrival time, at least in simple cases such as free motion.
Aharonov and Bohm introduced the operator \cite{AB}
\begin{equation}\label{top}
\widehat{T}=-m\bigg{(}\widehat{q}\frac{1}
{\widehat{p}}
-\frac{\hbar}{2i}\frac{1}{\widehat{p}^2}\bigg{)}=
-\frac{m}{2}\left(\widehat{q}
\frac{1}{\widehat{p}}+\frac{1}{\widehat{p}}
\widehat{q}\right)\,,
\label{to}
\end{equation}
by symmetrizing the classical expression for the time-of-arrival at
$X=0$ for a particle with
position $q$ and momentum $p$ at time $t=0$.
The same result is obtained by using
the quantization rules of Weyl, Rivier or Born-Jordan
\cite{SMP}, but this coincidence does not justify (\ref{to})
since, a priori, no
rule is more fundamental than any other rule. They simply
provide a number of quantum generalizations of a classical quantity;
one (or some)
of these can be selected
by requiring either agreement with a particular experimental procedure
or that certain desirable conditions be satisfied, such as
properties obeyed by the classical counterpart.
It is thus essential to carefully examine the behaviour of a given
operator and of its eigenstates to ascertain their ultimate
physical content.
Paul and others have examined the mathematical
properties of $\widehat{T}$ as a linear operator in Hilbert
space \cite{Paul,Raza,Werner,Gia}.
In section II we shall first review these properties. Then the
relations
between $\widehat{T}$ and positive operator valued measures, or
its self-adjoint variants are examined, and the important
question of the domains of applicability of the theories
associated with the different operators
is clarified.
In section III less abstract aspects are studied, such as
the functional
form of the eigenstates in the coordinate representation and their
time dependence. Finally, the
behaviour of normalized wave packets ``peaked''
around one of the eigenstates is discussed in section IV.

\section{TIME OPERATORS}
$\widehat{T}$ is not a self-adjoint operator but it is in a sense
the next best thing, namely, a maximal symmetric operator
\cite{Werner,Buschbook,AG}.
This means that it is hermitean (or ``symmetric'' in the mathematical
literature) but it does not admit
a self-adjoint extension. Indeed a first way to circumvent Pauli's
argument is to relax the standard requirement and admit
hermitean but non self-adjoint operators as representations of
physical properties. By imposing hermiticity and square integrability
of the states in its range, its
domain $D(\widehat{T})$ is given by the set of states with
momentum representation obeying \cite{Paul}
\begin{equation}\label{3/2}
\lim_{p\to 0}\frac{\phi(p)}{p^{3/2}}\to 0\,.
\end{equation}
The eigenvalue equation in momentum representation,
${\langle} p|\widehat{T}|T{\rangle}=T{\langle} p|T{\rangle}$,
\begin{equation}
-mi\hbar\left(\frac{1}{p}\frac{d}{dp}-\frac{1}{2p^2}\right)
{\langle} p|\psi{\rangle}=T{\langle} p|\psi{\rangle}\,,
\end{equation}
is solved (for $p\ne$ 0) by
\begin{equation}\label{eigen}
{\langle} p|T,\alpha{\rangle}=(|p|/mh)^{1/2}e^{ip^2T/(2m\hbar)}\Theta(\alpha p)\,,
\end{equation}
with $\alpha=\pm$.
These (improper) eigenfunctions have the appropriate transformation
behaviour for a state arriving at
time $T$,
\begin{equation}\label{trans}
e^{-i\widehat{H}t/\hbar}|T,\alpha{\rangle}=|T-t,\alpha{\rangle}\,,
\end{equation}
namely, if the original state had already evolved for a
time $t$ then it would arrive at $T-t$ instead;
but they are not square integrable functions so they do not
define physically realizable
state vectors in the Hilbert space
${\cal H}=L^2$.
The constant factor in (\ref{eigen}) is
chosen to satisfy the resolution of the identity
\begin{equation}\label{one}
1_{\rm op}=\sum_\alpha\int_{-\infty}^{\infty}
dT\, |T,\alpha{\rangle}{\langle} T,\alpha|\,.
\end{equation}
In the physical interpretation of this expression, one assumes that
all particles arrive sooner or later.
Notice
the inclusion of negative times. In an experimental context, all measured
arrival times are positive if the single particle state $\psi(t)$
of interest is prepared at the instant $t=0$. However,
as discussed in considerable detail by Grot et al. \cite{GRT},
the theoretical analysis is much simpler if one considers a different
problem in which it is imagined that the state is prepared at $t=-\infty$
in the state which (in the absence of interaction) evolves to the
actually prepared state $\psi(t=0)$, and the arrival time distribution
includes both negative as well as positive arrival times. This is
very useful from
the mathematical point of view (e.g. see Eq. (\ref{one}),
or Eqs. (\ref{symx0}) and (\ref{tau})
below.)

The eigenfunctions $|T,\alpha{\rangle}$ are complete,
but they are not orthogonal,
\begin{equation}
\label{overlap}
{\langle} T',\alpha'|T,\alpha{\rangle}=\frac{\delta_{\alpha,\alpha'}}{2}
\bigg[\delta(T-T')+\frac{i}{\pi}
{\cal P}\left(\frac{1}{T-T'}\right)\bigg]\,.
\end{equation}
The non orthogonality of the eigenstates has been
associated with the
intrinsically ``unsharp'' character of
this ``observable'' by the proponents of
the ``operational approach to
quantum mechanics'' \cite{Buschbook}.
These authors, and
Gianitrappani in particular \cite{Gia},  point out that the
states $|T,\alpha{\rangle}$  provide a ``positive operator valued
measure'' (POVM) for the arrival time. This means that for an interval of
time $[T_1,T_2]$ the positive,
bounded operator
\begin{equation}\label{POVM}
\widehat{B}(T_2,T_1)=\sum_\alpha\int_{T_1}^{T_2} |T,\alpha{\rangle}\,
dT\,{\langle} T,\alpha|
\end{equation}
can be constructed. This is not a projector, but one does not
necessarily require projectors in order to introduce probabilities
in quantum mechanics \cite{Buschbook,Toller}.
By taking the trace with the normalized physical density
operator $\widehat{\rho}$,
the function
\begin{equation}
P(T_2,T_1)=tr[\widehat{B}(T_2,T_1)\widehat{\rho}]\,
\end{equation}
fulfills in principle the conditions of a probability for arrival
between $T_1$ and $T_2$: it is positive,
additive for disjoint sets, and, using (\ref{one}),
tends to 1 as $T_1\to-\infty$ and $T_2\to\infty$.
In particular, an arrival time distribution  (for a pure state
$\psi$ at time $t$) is defined by
\begin{equation}
\label{distri}
\Pi(T;\psi(t))=\sum_{\alpha} |{\langle} T,\alpha|\psi(t){\rangle}|^2\,.
\end{equation}

According to (\ref{trans}), the POVM satisfies
the ``covariance condition''
\begin{equation}\label{cc}
e^{-i\widehat{H}t/\hbar}\widehat{B}(T_2,T_1)e^{i\widehat{H}t/\hbar}
=\widehat{B}(T_2-t,T_1-t)\,.
\end{equation}
and this means in particular that
\begin{equation}\label{cova}
\Pi(T;\psi(0))=\Pi(T-t;\psi(t))\,,
\end{equation}
In words,
the probability for arriving at $T$ is equal to the probability
for arriving at $T-t$
when the original state has evolved a time $t$.
This is a basic physical requirement that any
candidate for an arrival time distribution should obey.
This is of course also true in classical mechanics.

The time-of-arrival operator $\widehat{T}$ is the ``first moment''
of the POV measure. For any $\phi\in\cal H$ and $\psi\in D(\widehat{T})$,
\begin{equation}\label{tb}
{\langle} \phi|\widehat{B}^{(1)}|\psi{\rangle}
\equiv\int T {\langle}\phi|\widehat{B}(dT)|\psi{\rangle}
={\langle} \phi|\widehat{T}|\psi{\rangle}=
\sum_\alpha\int_{-\infty}^{\infty}
dT\,{\langle} \phi |T,\alpha{\rangle} T {\langle} T,\alpha|\psi{\rangle}\,.
\end{equation}
There are in principle other POV measures compatible with
(\ref{cc}) and (\ref{tb}), but (\ref{POVM}) has the
following unique feature:
Its second moment operator obeys, for states in $D(\widehat{T})$,
\begin{equation}\label{var}
\Delta(\psi)\equiv \int T^2 {\langle}\psi|\widehat{B}(dT)|\psi{\rangle}-
{\langle}\widehat{T}\psi|\widehat{T}\psi{\rangle}=0\,,
\end{equation}
which is by no means an obvious relation since the eigenvectors of
$\widehat{T}$ are not orthogonal (An equation like (\ref{var})
would be satisfied trivially by a self-adjoint operator since in that case
the moments of the distribution of the observable are obtained as the
expectation values of the powers of the operator.)
$\Delta(\psi)$ is called the
``variance form'' \cite{Werner}.
In classical mechanics the true arrival time
distribution makes the analogous quantity minimal so this
property has also been invoked to select
a proper quantum distribution \cite{Ki,Werner}.

Another condition satisfied by $\Pi(T)$ is that, for the state $\psi_1$
defined by ${\langle} p|\psi_1{\rangle}={\langle} p|\psi{\rangle}^*$,
\begin{equation}\label{symx0}
\Pi(T,\psi(0))=\Pi(-T,\psi_1(0))\,,
\end{equation}
which follows immediately from the
symmetry property
\begin{equation}\label{pt}
{\langle} p|T,\alpha{\rangle}^*={\langle} p|-T,\alpha{\rangle}\,.
\end{equation}
Eq. (\ref{symx0}) is also a classically motivated relation
that must hold when the arrival
point is $X=0$ \cite{Ki,Werner}. (If only covariance and
minimum variance were imposed the arrival point would not be
specified.)

The distribution of arrival times (\ref{distri}),
was first obtained by Kijowski for states with positive momenta
or negative momenta only ($\alpha=1$ or $\alpha=-1$ but not both)
\cite{Ki},
and has been later rederived, studied or generalized by several authors
\cite{Werner,Busch,GRT,DM,Leon,D}. Werner in particular \cite{Werner},
with  the same restriction imposed by Kijowski, justified the
uniqueness of (\ref{distri}) subject to the conditions
(\ref{cova}), (\ref{var}), and (\ref{symx0}). Grot, Rovelli and Tate
\cite{GRT}
introduced a regularized self-adjoint operator and considered
the full expression (\ref{distri}) for possible application
to more general states having both positive and negative momenta
but vanishing in the proximity
of $p=0$. Delgado and Muga \cite{DM},
with the restriction imposed by Kijowski and Werner (i.e., for states with
purely positive or negative momenta) arrive at the
distribution using a different self-adjoint operator.
We shall next discuss these two proposals
for self-adjoint arrival time operators.

\subsection{SELF-ADJOINT VARIANTS OF $\widehat{T}$}
Grot, Rovelli and Tate \cite{GRT} trace the
non-orthogonality of the eigenstates to the singularity at $p=0$, and
produce a ``regularized'' self adjoint time operator
$\widehat{T}_{\varepsilon}$
with eigenstates that differ from (\ref{eigen}) only in a small momentum
region around $p=0$,
\begin{equation}\label{GRTs}
{\langle} p|T\pm{\rangle}_\varepsilon=\Theta(\pm p)(h m f_\varepsilon(p))^{-1/2}
\exp\bigg(\frac{iT}{m\hbar}\int_{\pm\varepsilon}
^p\frac{dp'}{f_\varepsilon(p')}\bigg)\,,
\end{equation}
where ${\varepsilon}$ is a small positive number,
\begin{equation}
f_{\varepsilon}=\cases{
1/p\,, & $|p|>\varepsilon$\cr
\varepsilon^{-2}p\,, &$|p|<\varepsilon$\cr}\,,
\end{equation}
and $\Theta$ is the Heaviside ``step'' function.
The integral in (\ref{GRTs}), that plays the role of the energy
when divided by $m$, becomes
\begin{equation}
\int_{\pm\varepsilon}
^p\frac{dp'}{f_\varepsilon(p')}=\cases{
(p^2-\varepsilon^2)/2\,, & $|p|>\varepsilon$\cr
\varepsilon^{2}\ln(|p|/\varepsilon)\,, &$|p|<\varepsilon$\cr}\,,
\end{equation}
so that the
regularization amounts to changing the energy spectrum in the region
around $p=0$ by introducing negative energies.
In this manner Pauli's objection is avoided. While Grot, Rovelli and
Tate modify the eigenstates (and therefore the time operator) in
the proximity of $p=0$,
Paul discussed a similar idea \cite{Paul},
namely, to modify the physical state vector wave functions
infinitesimally around this point so that
$\widehat{T}$ can be applied to them.
He concluded however that such modifications were not physically
meaningful since these
infinitesimal changes may lead, for example, to arbitrary values
of the expectation value of the square of $\widehat{T}$,
${\langle} \widehat{T}^2 {\rangle}$.
Note that the small $p$ region is responsible for
the long time asymptotic behaviour of the arrival time
distribution (see the Appendix A) and therefore any
infinitesimal change there
affects drastically
quantities such as ${\langle} \widehat{T}^n{\rangle}$,
even though it will affect only infinitesimally other
expectation values, for example
${\langle} \widehat{q}^n{\rangle}$, or ${\langle} \widehat{p}^n{\rangle}$,
with $n=1,2,3...$ \cite{Paul}.
Paul's observation is correct, but we shall argue in the discussion that
the consequences are not necessarily as negative as he thought.

A second proposal that was already pointed out by Kijowski
(section 8 of \cite{Ki})
and has been developed  further by Delgado and Muga \cite{DM} is to
take the two pieces of $\widehat{T}$,
\begin{equation}
\widehat{T}=\widehat{T}\Theta(\widehat{p})+\widehat{T}\Theta(-\widehat{p})\,,
\end{equation}
the first  acting on the positive momentum subspace and the
second acting on the negative momentum subspace,
and combining them with a negative sign instead,
\begin{equation}
\widehat{T}_{-}=\widehat{T}\Theta(\widehat{p})-\widehat{T}\Theta(-\widehat{p})\,.
\end{equation}
This operator is self-adjoint, and has eigenstates $|T{\rangle}$ formed by
the combination
\begin{equation}\label{tau}
|T{\rangle}\equiv |T+{\rangle}+|-T-{\rangle}\,.
\end{equation}
This time operator avoids Pauli's argument in a different manner.
It is not conjugate to $\widehat{H}$ but to an operator
$\widehat{\cal H}={\rm sgn}(\widehat{p})\widehat{H}$ related to
$\widehat{H}$ by
the change of the sign for negative momenta.
Technically, because $\widehat{T}$  is a maximal symmetric operator,
$\widehat{T}_-$ cannot be its self-adjoint
extension, and only gives the same result as $\widehat{T}$ when acting
on the subspace of positive momenta (there is a change of sign
for the negative momentum subspace).

So one can in principle construct self-adjoint variants of $\widehat{T}$.
But what do they mean? and what new information can we extract from them?
Unfortunately,
neither the states $|T\pm{\rangle}_\varepsilon$
nor the states $|T{\rangle}$ transform according to
Eq. (\ref{trans}). Equivalently, none of the self-adjoint operators
discussed satisfies the covariance condition, and the
``arrival time distributions'' computed with them
do not in general satisfy the basic
physical requirement (\ref{cova}).
To avoid this problem, in both approaches the domains of {\it physical
applicability} of the self-adjoint operators have to be restricted
with respect to the mathematical domains.
For the approach by Grot, Rovelli and Tate
the domain should be restricted to states with momentum support
outside the regularization region around $p=0$,
$(-\varepsilon, \varepsilon)$,
so that in fact the arrival time distribution is
again given by (\ref{distri}).
Similarly, the approach associated with $\widehat{T}_-$
is only physically meaningful for states with purely
positive/negative momenta
and the corresponding distributions are therefore contained in
(\ref{distri}). In summary, the apparent advantage of these
two proposals is misleading, since in practice they are
only applicable when their results are equivalent to
the ones provided by the POVM
related to the Aharonov-Bohm operator.
In fact, the limitations on the domains of applicability of these two
approaches are quite severe. The arrivals at a screen or detector
will occur for all states in
$\cal{H}$ and not just for an especial set of states. A complete theory
should provide the arrival time distribution in all cases.
The important point is that the
distribution associated with the POV measure applies for all states in
$\cal{H}$ (irrespectively of their behaviour at $p=0$), so it is in this
sense a more complete approach. This is perfectly
compatible with $D(\widehat{T})\ne {\cal H}$ since the time operator
is only one of its moments.
Let us recall that a probability
distribution exists independently of the existence of its moments.
In fact in classical mechanics ensembles with non-zero probability at
$p=0$ have no finite average arrival time, but the distribution is
nevertheless well defined. In this respect only the POVM approach provides a
correct classical limit.

In summary, the distribution $\Pi(T)$ is satisfactory in many
ways. There is however an important point that
has not been considered yet. How do the eigenstates of $\widehat{T}$
behave? Do they
really represent states that arrive at a given time for a given
position?
The meaning of these eigenstates, although central,
has not been sufficiently discussed.
In particular, they
have always been studied in momentum representation without paying
attention to its coordinate representation and time dependence.
These aspects are examined in the next section.
Before doing so, it is noted that , because of the symmetry (\ref{pt}),
the coordinate representations of the time evolved states
$|T_t\pm{\rangle}\equiv\exp(-i\widehat{H}t/\hbar)|T\pm{\rangle}$
are related by
\begin{equation}\label{sym1}
{\langle} x|T_t-{\rangle}={\langle} -x|T_t+{\rangle}\,.
\end{equation}
so that it is enough to study one of the cases. Further simplification
comes from the fact that,
as a consequence of (\ref{trans}), studying the time dependence
of one of the states, say  $|T'_t+{\rangle}$, from $t=0$ to the nominal
arrival time $t=T'$,
is equivalent to considering the sequence of eigenstates $|T+{\rangle}$
from
$T=T'$ to $T=0$.
The space-time analysis that we shall carry out in the next section
should be taken with
some precaution because
these are not square integrable states.
As occurs with the continuum stationary states used in
scattering theory, their physical interpretation requires the
construction of normalizable wave packets peaked at one of
them.
\section{COORDINATE REPRESENTATION OF THE EIGENSTATES}
The coordinate representation of the wavefunction (\ref{eigen})
is given by the integral,
\begin{equation}\label{xinte}
{\langle} x|T+ {\rangle}=\int_0^\infty
{\langle} x|p{\rangle}\left(\frac{p}{mh}\right)^{1/2}
e^{ip^2T/(2m\hbar)} dp\,,
\end{equation}
where delta function normalization, ${\langle} x|p{\rangle}=h^{-1/2}
\exp(ixp/\hbar)$, is used for the plane waves. Let us assume for the
time being that $T>0$. The original path of integration is very
inefficient numerically because of the rapid oscillations of the
exponentials. The easiest way to calculate the integral is to deform
the contour in the complex $p$-plane along the imaginary axis from the
origin up to the intersection with the steepest descent path defined by
\begin{equation}
p_I=mx/T+p_R\,,
\end{equation}
($p_R$ and $p_I$ are, respectively, the real and imaginary parts of $p$)
and then to follow this steepest descent path rightwards to
infinity in the first quadrant.
The saddle is on the real axis at $-mx/T$.
Two cases have to be distinguished:
For $x>0$ the path does not cross the
saddle and the value of the integral is small. In this case the
origin is the only relevant critical point.
Use of Watson's lemma \cite{asym} for large $x$ and small $T$
provides the leading term
\begin{equation}\label{xpo}
{\langle} x|T+{\rangle}\sim\frac{h^{1/2}e^{3i\pi/4}}{x^{3/2}\pi 2^{5/2}
m^{1/2}}\,,\:\:\:x\to\infty\,,
\end{equation}
note that this asymptotic behaviour is independent of $T$.
But for $x<0$ the saddle becomes the dominant critical point. Retaining
the leading term \cite{asym},
\begin{equation}\label{xne}
{\langle} x|T+{\rangle}\sim\frac{(m|x|/h)^{1/2}}{T}e^{i\pi/4}
e^{-\frac{ix^2m}{2\hbar T}}\,,\:\:\: x\to-\infty\,.
\end{equation}
In the first case
${\langle} x|T+{\rangle}$
decreases as $x^{-3/2}$, whereas in the second one it increases as
$|x|^{1/2}$ asymptotically.
In fact the integral can be expressed exactly in terms of
parabolic cylinder
functions, but the approximate critical point treatment
just outlined is worthwhile since it
allows a simple rationalization of the exact results.
By means of the change of variable
\begin{equation}
s=\alpha p\,,
\end{equation}
where
\begin{equation}
\alpha=-\left[\frac{T}{m\hbar}\right]^{1/2}e^{-i\pi/4}\,,
\end{equation}
(\ref{xinte}) becomes
\begin{equation}\label{in}
\frac{1}{hm^{1/2}(\alpha^{1/2})^3}\int_C ds\, s^{1/2}
e^{-s^2/2+zs}\,,
\end{equation}
where all square roots are calculated with a branch cut on
the negative real axis, and
\begin{equation}
z=xe^{-i\pi/4}\left(\frac{m}{\hbar T}\right)^{1/2}.
\end{equation}
The integration path $C$ in (\ref{in}) goes now from $0$ to $\infty$
along the bisector of the second quadrant in the complex $s$ plane.
Note that $z$ is the saddle point of $\exp(w)\equiv\exp(-s^2/2+sz)$. It
lies on the bisector of the second/fourth quadrant for
positive/negative $x$. The steepest descent paths from it
are parallel to the real axis, and
$C$ lies on the border between ``hill'' and ``valley''
($w_R\equiv {\rm Real}(w)=0$). It can however be deformed
into a line just above the branch cut since there is no contribution
at $\infty$. As the integrand across the branch cut simply
changes sign, (\ref{in}) is one half of the loop integral
around the cut.
In this fashion one of the integral forms of the parabolic
cylinder function $D_{-3/2}(z)$ can be recognized \cite{AS},
\begin{equation}\label{xt}
{\langle} x|T+{\rangle}=\frac{\Gamma(3/2)}
{hm^{1/2}(\alpha^{1/2})^3}e^{z^2/4}D_{-3/2}(z)\,.
\end{equation}
The asymptotic behaviour
for large $|z|$, see \cite{Grand}, is in agreement with the expressions
(\ref{xpo}) and (\ref{xne}).
The corresponding results for the case $T<0$ are simply
obtained by using the symmetry
\begin{equation}\label{sym2}
{\langle} x|T +{\rangle}^*={\langle} -x|-T+{\rangle}\,.
\end{equation}
The
asymptotic behaviour for negative and positive $x$ (growing and decaying
respectively) changes abruptly
at $T=0$. Since $|T+{\rangle}$  does not represent a physical state
vector, the discontinuity at $T=0$ is not problematic, but indicates again
that a literal physical interpretation of these states is not allowed.
We shall see later that normalized wavepackets formed with these
states do not present this singularity.

There is a region close to $x=0$ that cannot be described by
the asymptotic formulae for large argument. In that region, however,
the parabolic cylinder
function can be expressed by means of a power series
\cite{AS},
\begin{equation}\label{x0}
D_{-3/2}(z)=\frac{\Gamma(-1/4)}{\pi^{1/2} 2^{5/4}}
-\frac{\Gamma(1/4)}{\pi^{1/2} 2^{3/4}} z + {\cal O}(z^2)\,.
\end{equation}
Combining (\ref{x0}) and (\ref{xt}) explicit expressions
for ${\langle} x=0|T+{\rangle}$ and for its ``flux'' $J(x=0)$ can be obtained
\cite{note}. In particular,
\begin{equation}
J(x=0)=\frac{[\Gamma(3/2)\Gamma(-1/4)]^2}
{(2\pi)^3 T^2 2^{3/2}}\,.
\end{equation}
For $T>0$, neither the wave function ${\langle} x|T+{\rangle}$ nor the flux
are zero at $x=0$.
Both the ``density'' $|{\langle} x=0|T+{\rangle}|^2$ and
$J(x=0)$ grow monotonically
as $T\to 0$.
Figures 1a and 1b illustrate all the dependences discussed.
$|{\langle} x|T+{\rangle}|^2$ is depicted for a
series of decreasing times $T>0$ as a function of $x$
for two different scales:
between $x=-2$ and $x=0.2$ (1a), and between $x=-0.2$ and $x=0.2$
(1b) (atomic units are used in all numerical examples).
In the larger scale the density of the eigenstate
is essentially a straight line, pivoting at $x=0$,
that approaches the vertical
as $T\to 0$. The finer scale however shows that
the arrival is not sharply defined.
Even though, in a loose sense, ``most of the wave'' ${\langle} x |T+{\rangle}$
passes from $x<0$ to $x>0$ at
$T=0$, there is a tail at $x>0$ present for an arbitrary $T$.
Since these states are not normalized to one it is not possible to
quantify the fraction of particles that can be found to the
right of $x=0$ before $T=0$.

One might think that the eigenstates of the self-adjoint operators, which are
orthogonal, could avoid this
type of unsharpness \cite{sharp}. But they don't.
For arbitrary values of $\epsilon$ we have numerically checked
(see the Appendix B) that there is a non vanishing
density  $|{\langle} x|T+{\rangle}_\varepsilon|^2$ at $x>0$
for fixed $T$.
Similarly, the states ${\langle} x|T{\rangle}$ are also non
zero on the right hand side. Using their defining equation (\ref{tau}),
(\ref{sym1}) and (\ref{sym2}), one finds
\begin{equation}
{\langle} x|T{\rangle}=2{\rm Re} {\langle} x|T +{\rangle}\,.
\end{equation}
(It is to be noted that the corresponding time evolved state does
not take this form.
Instead, ${\langle} x|T_t{\rangle}\equiv{\langle} x|e^{-iHt/\hbar}|T{\rangle}=
{\langle} x|(T-t)+{\rangle}+{\langle} x|(T+t)+{\rangle}^*$.)
Figure 2 represents the square of this quantity for a relatively large
spatial interval (the small, but  non-zero density for
$x>0$ is not seen in this scale).
The figure corresponding to the imaginary part is very similar.
Note the increasingly
rapid oscillation as $x\to -\infty$. Cancellation of these oscillations
occurs in linear combinations of the states
${\langle} x|T{\rangle}$ or ${\langle} x|T+{\rangle}$ over a
nonzero band leading to localization in a region
closer to the origin, see Figure 3 and the discussion in the following
section.
There is a helpful classical association to understand the oscillatory
pattern: high velocity particles (fast oscillations)
have to start at longer distances
and low velocity particles (slow oscillations) at shorter distances
if they all have to arrive
at $X=0$ simultaneously.
\section{COORDINATE AND TIME DEPENDENCE OF THE
NORMALIZED QUASI-EIGENSTATES}
We shall here construct normalized wavepackets
by using a Gaussian distribution of the states
${\langle} x|T'_t +{\rangle}$ peaked at $T'=T$ \cite{Paulfn}.
\begin{equation}\label{norma}
{\langle} x|\Psi(t;T,\Delta T){\rangle}=
N\int_{-\infty}^{\infty}
e^{-\frac{(T-T')^2}{2(\Delta T)^2}}{\langle} x|T'_t +{\rangle}\,dT'\,.
\end{equation}
Carrying out the Gaussian integral over $T'$ and then determining the
normalization constant $N$ by evaluation of
$\int_{-\infty}^{\infty} dx\,|{\langle} x|\Psi(t;T,\Delta T){\rangle}|^2$ gives the
normalized wave packet
\begin{equation}\label{coor}
{\langle} x|\Psi(t;T,\Delta T){\rangle}=
\frac{2\pi^{1/4}(\Delta T)^{1/2}}{hm^{1/2}}
\int_0^\infty p^{1/2} e^{ipx/\hbar} e^{i(T-t)p^2/2m\hbar}
e^{-(\Delta T)^2p^4/8m^2\hbar^2}dp\,.
\end{equation}
This integral is readily evaluated numerically because of the
exponential dependence on $-p^4$. The momentum representation is
given explicitly by
\begin{equation}\label{momre}
{\langle} p|\Psi(t;T,\Delta T){\rangle}=\frac{2\pi^{1/4}(\Delta T)^{1/2}}
{(hm)^{1/2}}p^{1/2}e^{i(T-t)p^2/(2m\hbar)}
e^{-(\Delta T)^2p^4/(8m^2\hbar^2)}\Theta(p)\,.
\end{equation}
As $\Delta T\to 0$ these states are orthogonal
and satisfy the eigenvalue equation to any desired degree
of accuracy:
\begin{eqnarray}\label{overl}
{\langle}\Psi(t;T',\Delta T)|\Psi(t;T,\Delta T){\rangle}&=&w\left(
\frac{T-T'}{2\Delta T}\right)=\cases{
1 &${\rm if}\; T=T'$\cr
\sim {\cal{O}}\left(\frac{\Delta T}{T-T'}\right)& ${\rm as}\;
\frac{\Delta T}{T-T'}\to 0$\cr}\\
{\langle} p|\widehat{T}|\Psi(0;T,\Delta T){\rangle}&=&T{\langle} p|\Psi(0;T,\Delta T){\rangle}
+{\cal O} (\Delta T)^{5/2}\,,
\end{eqnarray}
where $w(z)=\exp(-z^2){\rm erfc}(-iz)$ is the ``$w$-function'' \cite{AS}.
The moments ${\langle} \widehat{p}^n {\rangle}$ of the momentum
distribution can be obtained from (\ref{momre}),
\begin{equation}
{\langle} \widehat{p}^n{\rangle}=\Gamma\left(\frac{n+2}{4}\right) \pi^{-\frac{1+n}{2}}
\left(\frac{\Delta T}{hm}\right)^{-n/2}\,.
\end{equation}
Both the average energy ${\langle} E{\rangle}\equiv{\langle} \widehat{p}^2{\rangle}/2m=
h/(2\pi^{3/2}\Delta T)$ and the variance
$\Delta E=[2^{-1}-\pi^{-1}]^{1/2}
\hbar/\Delta T$ diverge as $\Delta T$ goes to zero.

Similarly, from (\ref{coor}) one finds
\begin{equation}
{\langle} \widehat{x}{\rangle}
=-\frac{\Gamma(3/4)}{\pi}
\left(\frac{h}{m\Delta T}\right)^{1/2}(T-t)\,,
\end{equation}
whereas higher moments diverge. (A finite ``width'' can however be
defined as the half width at half height.)
The average velocity ${\langle} \widehat{p}/m{\rangle}$, which is also the
velocity of the centroid
${\langle} \widehat{x}{\rangle}$ is given by $\Gamma(3/4)[h/(\Delta T m
\pi^2)]^{1/2}$. The behaviour of these states is arbitrarily close to what
one could desire for an ideal arrival time eigenvector: They are normalized,
and there is no discontinuity at $t=T$;
for a fixed $\Delta T$,
they still obey a transformation law of the form (\ref{trans});
the wave travels towards the origin with constant velocity, and the
average position crosses the origin at time $t=T$, which is also the time
when the spatial width of the wave packet attains its minimum value.
Of course a certain
unsharpness remains: From
(\ref{coor}) one finds
${\langle} -x|\Psi(2T-t;T,\Delta T){\rangle}={\langle}\Psi(t;T,\Delta T)|x{\rangle}$, so the
probability density has inversion symmetry in $(x,t)$ with respect to
the space time ``origin'' $(0,T)$.
In particular, for $t=T$, half the norm is
to the right of $x=0$ independently of $\Delta T$
(The question of the existence of quantum states
where the particle stays strictly
on one side of $X$ before $T$ and on the other side after $T$
is addressed in the Appendix C.)
However, the
wave packet density and flux are
more and more peaked at the space-time point ($x=0,\, t=T$)
as $\Delta T \to 0$, so that
the passage of probability from
left to right is sharp to any desired accuracy,
\begin{eqnarray}
|{\langle} x=0|\Psi(t=T;T,\Delta T){\rangle}|^2&=&
\left(\frac{m}{h\Delta T}\right)^{1/2}\frac{\Gamma(3/8)^2}
{\pi 2^{5/4}}\\
J(x=0,t=T)&=&\frac{\Gamma(3/8)\Gamma(5/8)}{\Delta T \pi^{3/2}
2^{3/4}}\,.
\end{eqnarray}
The normalized wavepackets constructed with the eigenstates (\ref{GRTs})
have been examined in \cite{Jonathan}, where it is reported that this
arbitraryly sharp accuracy is not found in that case.

\section{Discussion}
Given the importance of the timed detection of particles
at screens or in time-of-flight experiments, finding a theoretical
description of the arrival times seems imperative.
Apparatus dependent results are available, for example via modelling
the detection with
complex absorbing potentials or other measurement
models \cite{AP,Aetall,Halli}, but it
is reasonable to inquire
if an intrinsic, apparatus-independent distribution can be
naturally defined by means of the usual operator approach to
quantum mechanics. We have seen in sections II and III that the POVM
and the corresponding distribution $\Pi(T)$ associated with the time
operator $\widehat{T}$ provide a rather satisfactory answer from the
point of view of the properties satisfied: covariance with respect to
time translations, minimum variance, appropriate symmetries,
physically correct
domain of applicability, and sharp space-time behaviour of the
normalized quasi-eigenstates.

One of the objections by Paul to $\widehat{T}$,
which is also a shortcomming of the theories based on self-adjoint operators,
namely the restrictive domain of the time operator(s) (which does not include,
for example, states such as minimum uncertainty product Gaussians)
is overcome by the POVM theory.
An important point is to consider the primary object as the POVM
(equivalently the resolution of the identity, or the arrival time
distribution) rather than the operator.
In this manner, the domain of applicability of the theory
is $\cal{H}$, and the classical limit is correct.
Moreover in this light the other problem indicated by
Paul, namely the
extreme sensibility of the expectation values of powers of $\widehat{T}$
to small perturbations, is relatively unimportant.
It is a fact that some quantities are very sensitive to certain small
perturbations and we are simply  dealing with one of them. The moments
in the classical case would also suffer from such a sensitivity.
The moral is that one should not pay as much attention to the moments
of a time-of-arrival distribution (highly unstable with respect to small
perturbations or changes in the apparatus resolution, and possibly
divergent) but to general
features of the distribution (peaks, global form, or width at half
height for example). A consequence is that an uncertainty principle in
terms of ${\langle} \widehat{T}^2{\rangle}$ is not of much use, since this quantity
will generally diverge. If it does not, it will be too unstable with
respect to small perturbations. It is preferable to express the
uncertainty principle in terms of other measurements of width, such as
the the half width at half height.

Apart from the the positive features of $\Pi(T)$,
it is also fair to point out several unclear points or
open questions.
For example, it is necessary to generalize the treatment
to higher dimensions or to scattering
problems.
The connection between $\Pi(T)$ and measurements
(especially in non-classical cases) is also pending, see also
\cite{Aetall}.
We can however advance as a preliminary analysis two
important features of
a hypothetical measurement of $\Pi(T)$:
The resolution of identity (\ref{one}) implies
the following structure of $\Pi(T)$,
\begin{equation}\label{ddd}
\Pi(T)=|{\langle} T+|\psi_+{\rangle}|^2+|{\langle} T-|\psi_-{\rangle}|^2\,,
\end{equation}
where $|\psi_\pm{\rangle}=\Theta(\pm \widehat{p})|\psi{\rangle}$.
This implies ignoring the ``interference terms'' $|\psi_-{\rangle}{\langle}\psi_+|$
and $|\psi_+{\rangle}{\langle}\psi_-|$, so that in a hypothetical operational procedure to
measure $\Pi(T)$, only the diagonal terms of the density operator
contribute,
\begin{equation}
\widehat{\rho}\to |\psi_+{\rangle}{\langle}\psi_+|+|\psi_-{\rangle}{\langle}\psi_-|\,.
\end{equation}
Apart from the fact that the
practical implementation of this diagonalization may be cumbersome,
the neglect of interferences is not a desirable feature, since
many different quantum states would give the
same distribution.
A second problematic aspect is associated with the interpretation of
$|{\langle} T+|\psi_+{\rangle}|^2$ as the contribution to $\Pi(T)$ from particles
arriving from the left and of $|{\langle} T-|\psi_-{\rangle}|^2$ as the contribution
from those arriving from the right. This is particularly evident when
the {\it backflow effect},
\begin{eqnarray}
J\le 0\;\; &{\rm for}&\;\; |\psi{\rangle}=|\psi_+{\rangle}\,,\\
J\ge 0\;\; &{\rm for}&\;\; |\psi{\rangle}=|\psi_-{\rangle}\,,
\end{eqnarray}
occurs. Then
$\Pi(T)$ assigns zero probability to arrivals from the ``anomalous side'',
(e.g. from the right when $|\psi{\rangle}=|\psi_+{\rangle}$).
This implies that if $|\psi{\rangle}=|\psi+{\rangle}$ or  $|\psi{\rangle}=|\psi-{\rangle}$
then either particles should be found to arrive either only from
the left or only
from the right, respectively, even during the time interval when $J$ has the
``wrong sign'', or that the theory is appropriate when the ``screen''
is ``one-sided'', failing to detect any particles arriving from the
``anomalous
side''. Now consider the corresponding implication
for the general state $|\psi{\rangle}=|\psi+{\rangle}+|\psi-{\rangle}$ with both $|\psi+{\rangle}$
and $|\psi-{\rangle}$ nonzero. One possibility is that the interference terms
do not in fact contribute to the intrinsic arrival time distribution.
The other is that the distribution (\ref{ddd}) is only
appropriate when the apparatus measures the sign of the momentum of each
incident particle, thus collapsing the wavefunction of that particle either
to $|\psi+{\rangle}$ or $|\psi_-{\rangle}$, and then switches on the appropriate
one-sided detecting screen.

\acknowledgments{
One of us, J. G. M., acknowledges a very useful correspondence with
Marco Toller. The content of this paper was presented as a talk in
the workshop ``Time in quantum mechanics'', La Laguna, May 1998.
The work has been supported by Gobierno Aut\'onomo
de Canarias (Spain)
(Grant No. PB2/95), and by CERION.}

redpajama_set_name: RedPajamaArXiv


text:
\section*{\refname}}

\begin{document}
\title{Enriching Visual with Verbal Explanations for Relational Concepts -- Combining LIME with Aleph}
\titlerunning{Enriching Visual with Verbal Explanations for Relational Concepts}
\author{Johannes Rabold \and
Hannah Deininger \and
Michael Siebers \and
Ute Schmid}
\authorrunning{J. Rabold et al.}
\institute{Cognitive Systems, University of Bamberg, Germany}
\maketitle
\begin{abstract}
With the increasing number of deep learning applications, there is a growing demand for explanations. Visual explanations provide information about which parts of an image are relevant for a classifier's decision. However, highlighting of image parts (e.g., an eye) cannot capture the relevance of a specific feature value for a class (e.g., that the eye is wide open). Furthermore, highlighting cannot convey whether the classification depends on the mere presence of parts or on a specific spatial relation between them. Consequently, we present an approach that is capable of explaining a classifier's decision in terms of logic rules obtained by the Inductive Logic Programming system Aleph. The examples and the background knowledge needed for Aleph are based on the explanation generation method LIME. We demonstrate our approach with images of a blocksworld domain. First, we show that our approach is capable of identifying a single relation as important explanatory construct. Afterwards, we present the more complex relational concept of towers. Finally, we show how the generated relational rules can be explicitly related with the input image, resulting in richer explanations.

\keywords{XAI \and Deep Learning \and Inductive Logic Programming.}
\end{abstract}
\section{Introduction}

\input{introductionLime.tex}

\section{Explaining Relational Concepts with LIME-Aleph}

\subsection{Core Concepts}

\subsubsection{The ILP System Aleph.}

To find symbolic explanations for relational concepts we use Aleph \cite{srinivasan2001aleph}. Aleph infers a logic theory $T$ given a set of positive ($E^+$) and negative ($E^-$) examples. An example is represented by the target predicate (e.g. \texttt{stack(e1).} or \texttt{not stack(e2).}) together with additional predicates (e.g. \texttt{contains(b1, e1).}) as background knowledge (BK). Predicates in BK are used to build the preconditions for the target rules.
Aleph is based on specific-to-general refinement search. It finds rules covering as many positive examples as possible, avoiding covering negative ones. Search is guided by modes which impose a language bias. The general algorithm is \cite{srinivasan2001aleph}:
{\small
\begin{enumerate}
	\item As long as positive  exist, select one. Otherwise halt.
	\item Construct the most-specific clause that entails the selected example and is within the language constraints
	\item Find a more-general clause which is a subset of the current literals in the clause.
	\item Remove  covered by the current clause.
	\item Repeat from step 1.
\end{enumerate}}
An example of a rule from $T$ in Prolog is {\small \texttt{stack(Stack) :- contains(Block1, Stack), contains(Block2, Stack),} \texttt{on(Block2, Block1).}} \\denoting that a \texttt{stack} is defined by one block on top of another.

\subsubsection{LIME's Identification of Informative Super-Pixels.}
LIME (\textbf{L}ocal \textbf{I}nterpre\-table \textbf{M}odel-Agnostic \textbf{E}xplanations) is an approach to explain the decision result of any learned model \cite{ribeiro2016should}. Explanations state the parts of an instance that are positively or negatively correlated to the class. It works by creating a simpler, local surrogate model around the instance to be explained. In case of an image, the explanation is a set of connected pixel patches called \textit{super-pixels}.

Let $x$ be an image and $x'$ be the binary vector that states whether super-pixels $x'_i \in x'$ are switched on or off (see below). LIME finds a sparse linear model $g(x')$ that locally approximates the unknown decision function $f(x)$ represented by a black-box classifier. It effectively finds the coefficients $\vec{w}$ for the super-pixel representations being variables in a simplified linear model. This is done by generating a pool of perturbed examples $z'$ by taking the original super-pixel representation $x'$ and randomly selecting elements in a uniformly distributed fashion. That way, images $z$ are obtained with some super-pixels still original and some altered according to a transform function $h$ effectively removing the information they contained (Switching them off). Each sample $z'$ (The binary vector indicating if super-pixels are switched off in this sample) is stored in a sample pool $\cal{Z}$ along with the classifier result $f(z)$ and a distance measure $\pi_x(z)$ that expresses the distance of the perturbed example $z$ to the original image $x$. For images this can be the Mean Squared Error. The distance is needed for the linear model to be locally faithful to the original function $f(x)$ and thus has to be minimized. The ``un-faithfulness'' of the model $g$ to the black-box model $f$ with respect to the distance measure $\pi_x(z)$ is expressed with the following formula \cite{ribeiro2016should}:

\begin{center}
	$\mathcal{L}(f, g, \pi_x) = \sum_{z, z' \in \mathcal{Z}} \pi_x(z) (f(z) - g(z'))^2$.
\end{center}

\noindent The goal is to find the coefficients $\vec{w}$ for $g$ that minimize this un-faithfulness $\mathcal{L}$. The coefficients ultimately translate back to weights for the super-pixels. LIME uses K-Lasso to find the weights~\cite{efron2004hastie}.

The original LIME uses the algorithm Quick Shift \cite{vedaldi2008quick} to find super-pixel. It imposes an irregular pixel mask over the input image that segments it in terms of pixel similarity. The segmentation is performed in a 5D space consisting of the image space and the color space. Quick Shift is only one of several segmentation algorithms that are available for LIME. They all share the attribute of imposing an irregular mask over an image. In many domains, this irregularity is not wanted. For the domain used in this paper it is preferable to use a segmentation algorithm that divides an image into a regular grid with square cells.

\begin{figure}[t]
	\centering
	\subfigure[All positions in the image, where a block has to be located in order to be \textbf{l}eft, \textbf{r}ight, \textbf{t}op, or \textbf{b}ottom of block a.]{\label{fig:rel_diagram_a}\includegraphics[width=30mm]{left_right_top_bottom_diagram.png}}
	\quad\quad\quad
	\subfigure[ for the relations \texttt{on} and \texttt{under}, illustrated with \texttt{on(b,a)}, \texttt{on(c,d)}, resp.~\texttt{under(a,b)} and \texttt{under(d,c)}.]{\label{fig:rel_diagram_b}\includegraphics[width=30mm]{on_under_diagram.png}}
	\caption{\label{fig:relation_diagrams} Diagrams to show the different concepts of the relations used.}
\end{figure}

\subsection{Extraction of Image Parts and Their Relations}\label{relex}
Based on image segmentation into a grid of super-pixels $i$ with domain-specific cell size, a set of attributes $A_i$ for cells and spatial relations between cells can be automatically extracted. Attributes $A_i$ are taken from a pool of attributes $\cal{A}$. An example for an attribute in $\cal{A}$ is the mean color of $i$ in the RGB color space. To find a human-comprehensible name, the nearest color according to the Euclidean distance in a pool of commonly known color names is assigned. Other extractable attributes are the size or the general location in the image. The coordinates of the center point of $i$ are stored for spatial reasoning. Extracted attributes are converted into predicates for BK. The attribute that a given super-pixel \texttt{SP} is blue is represented as \texttt{has\_color(SP, blue).}

Spatial relations can be defined between pairs of super-pixels. To restrict the number of pairs, we need a pre-selection $S$ of super-pixels that might be relevant for the concept. LIME's $\vec{w}$ describe the magnitude of relevance for either the true classification (positive weight) or the counter-class (negative weight). By introducing a user-defined constant $k$, we restrict how many super-pixels the selection $S$ should contain, taking the $k$ super-pixels with the highest values in $\vec{w}$. Spatial relations $r: S \times S$ are drawn from a pre-defined pool $\cal{R}$. For this work, we use the relations \texttt{left\_of, right\_of, top\_of, bottom\_of} as well as relations that represent an immediate adjacency in the regular grid mentioned earlier, namely \texttt{on} and \texttt{under}. Relations are defined with respect to the center coordinates of the super-pixels in $S$. Figure~\ref{fig:relation_diagrams} sketches the underlying semantics of these relations. It is possible to include additional relations as long as they are automatically extractable and their inverses are defined in the image space. In domains with super-pixels that differ in size, a \texttt{larger} relation between super-pixels could be defined. Also, a \texttt{not\_equal} relation can be considered.

\begin{comment}
This is achieved by approximating the complex decision boundary function $f$ of the classifier by a sparse linear function $g$ that is locally faithful to the instance to be explained. The super-pixels in the segmented version of the instance act as features with corresponding weights for the linear function. The goal of LIME is to assign high positive weights to super-pixels that are highly relevant for a decision and high negative weights to ones that are highly relevant for the counter-class. We can express local unfaithfulness $\mathcal{L}(f, g, \pi_x)$, with $\pi_x$ being a kernel for a distance function to the original instance $x$, with the following formula:

\begin{center}
	$\mathcal{L}(f, g, \pi_x) = \sum_{z, z' \in \mathcal{Z}} \pi_x(z) (f(z) - g(z'))^2$
\end{center}

\noindent $\mathcal{Z}$ is a sample of tuples $\langle z, z' \rangle$ where $z'$ is a binary vector indicating that a subset of super-pixels is present in an image $z$. The goal is to minimize the unfaithfulness $\mathcal{L}$ by finding the coefficients $\vec{w}$ for $g$ which ultimately translate back to weights for the super-pixels.

For our selection $S$ of important super-pixels, we take the $k$ super-pixels with the highest absolute $w_i \in \vec{w}$. $k$ is a user-defined parameter which can be tuned to match the domain. $S$ contains the $k$ super-pixels with the highest absolute value for the LIME weights. [We also include negatively correlated super-pixels since ...]
\end{comment}

\subsection{Learning Rules for Relational Concepts via Aleph}

To generate symbolic relational explanations for visual domains, we combine LIME's super-pixel weighting with Aleph's theory generation. The input into LIME-Aleph is an image $x$ and a model $f$ returning class probability estimations for $x$. Currently our approach is only applicable for explaining one class in contrast to all other classes, effectively re-framing the original classification as a concept learning problem. The output of LIME-Aleph is a theory $T$ of logic rules describing the relations between the super-pixels that lead to the class decision.

\begin{algorithm}[t]
	\caption{\label{alg:lime-aleph} Explanation Generation with LIME-Aleph.}
	\begin{algorithmic}[1]
		\State \textbf{Require:} Instance $x \in X$
		\State \textbf{Require:} Classifier $f$,  Selection size $k$, Threshold $\theta$
		\State \textbf{Require:} Attribute pool $\cal{A}$, Relation pool $\cal{R}$

		\State ~~~$S \leftarrow LIME(f, x, k)$ \algorithmiccomment{Selection of $k$ most important super-pixels.}
		\State ~~~$A \leftarrow$ extract\_attribute\_values$(S, \cal{A})$ \algorithmiccomment{Find all attribute values $A_i$ for all $i \in S$.}
		\State ~~~$R \leftarrow$ extract\_relations$(S, \cal{R})$ \algorithmiccomment{Find all relations $r: S \times S$ between all $i \in S$.}
		\State ~~~$E^+ \leftarrow \{\langle A, R \rangle\}$
		\State ~~~$E^- \leftarrow \{\}$
		\State ~~~\textbf{for each} $r(i,j) \in R$ \textbf{do}
		\State ~~~~~~$z \leftarrow$ flip\_in\_image$(x,i,j)$ \algorithmiccomment{Flip the super-pixels in the image space.}
		\State ~~~~~~$r' \leftarrow r(j,i)$ \algorithmiccomment{Obtain new predicate for the BK by flipping parameters.}
		\State ~~~~~~$R' \leftarrow R \setminus \{r\} \cup \{r'\}$ \algorithmiccomment{All relations in the BK; also the altered one.}
		\State ~~~~~~$R' \leftarrow calculate\_side\_effects(R', r')$ \algorithmiccomment{Re-calculate relations that are affected by the flipped relation.}
		\State ~~~~~~$c' \leftarrow f(z)$ \algorithmiccomment{Obtain new estimator for the perturbed image.}
		\State ~~~~~~\textbf{if} $c' \geq \theta$ \textbf{do} \algorithmiccomment{If estimator reaches threshold, add new positive example.}
		\State ~~~~~~~~~$E^+ \leftarrow E^+ \cup \{\langle A, R' \rangle\}$
		\State ~~~~~~\textbf{else} \algorithmiccomment{Else, add negative example.}
		\State ~~~~~~~~~$E^- \leftarrow E^- \cup \{\langle A, R' \rangle\}$
		\State ~~~\textbf{end for}
		\State ~~~ $T \leftarrow$ Aleph$(E^+, E^-)$ \algorithmiccomment{Obtain theory $T$ with Aleph.}
		\State \textbf{return} $T$
	\end{algorithmic}
\end{algorithm}

LIME's explanation relies on that linear surrogate model which contains the set of super-pixels with the highest positive weights for the true class. When dealing with the question which relations contribute most to the classification, identifying the most informative super-pixels has to be replaced by identifying the most informative \emph{pairs} of super-pixels. Instead of turning super-pixels on and off, LIME-Aleph inverts extracted relations between super-pixels and observes the effects on the classification. Algorithm~\ref{alg:lime-aleph} shows our approach. Given the selection $S$ of super-pixels together with the extracted attributes, our approach first finds all relations $R \subseteq \cal{R}$ that hold between them. For every relation $r(i,j) \in R$, a new perturbed example $z$ from the image space is created by flipping the super-pixels $i$ and $j$ in the image space. To generate a new example for Aleph, the resulting perturbed image is first put through the classifier $f$. If the estimator $f(z)$ exceeds a threshold $\theta$ for the class we want to explain, a new positive example is declared. Otherwise, the example is declared negative. All relations holding for the perturbed image are written in the BK characterizing this example. The initial positive example for Aleph is always generated for the unaltered constellation.

\section{Experiments and Results}

We investigate the applicability of LIME-Aleph in a blocksworld domain consisting of differently colored squares that can be placed in a regular-grid world. For a first investigation, we decided to focus on artificially generated images rather that real world domains.

\subsection{An Artificial Dataset for Relational Concepts}
We implemented a generator to create a huge variety of positive and negative example images for different blocksworld concepts.
All generated images are of size $32 \times 32$ pixels consisting of a single-colored red background, containing constellations of colored squares of dimension $4 \times 4$ pixels. The squares are single-colored (excluding red) with color-channel values either being set to 0.0 or 0.8. The squares are placed into the image according to an $8 \times 8$ uniform grid.
Positive examples are generated by first randomly placing a reference square. Then, the other squares are placed randomly following the relation conventions shown in Figure~\ref{fig:relation_diagrams}. For the experiments we restricted $\cal{A} = $ \texttt{\{color\}} with attribute values in \texttt{\{cyan, green, blue\}} and $\cal{R} = $ \texttt{\{left\_of, right\_of, top\_of, bottom\_of, on, under\}}.

\subsection{Training a Black-Box Model}

To obtain a black-box model for image classification, we used a small convolutional neural network \cite{krizhevsky2012imagenet} which we trained from scratch with commonly known best-practice hyper-parameters. The network consists of two convolution layers with kernel size $2 \times 2$ and ReLU activations. Each layer learns 16 filters to be able to robustly recognize the colored squares. After flattening the output, the convolution layers are followed by 2 fully connected layers each with a ReLU activation. The first layer consists of 256 neurons, the second one of 128 neurons. A small amount of dropout is applied past each layer to cope with potential overfitting \cite{srivastava2014dropout}.
The network does not contain a pooling layer. That way, fewer location information is lost in aggregation during the learning process which we believe is crucial for preserving spatial relationships~(see \cite{goodfellow2016deep} p. 331).
For the experiments, we generated perfectly balanced datasets with 7.000 training-, 2000 validation- and 1000 test-images for both the concept and the counter-examples. We trained the networks for a maximum of 10 epochs with early stopping if the validation loss did not decrease after 5 epochs.

\subsection{Experiment 1: Single Relation Concept}

\begin{figure}[b]
	\centering
	\subfigure[]{\label{fig:first_ex_a_pos}\includegraphics[width=20mm]{first_a_pos.png}}
	\quad\quad\quad
	\subfigure[]{\label{fig:first_ex_b_pos}\includegraphics[width=20mm]{first_b_pos.png}}
	\quad\quad\quad
	\subfigure[]{\label{fig:first_ex_c_neg}\includegraphics[width=20mm]{first_c_neg.png}}
	\caption{\label{fig:first_example} Positive ($a$, $b$) and negative ($c$)  for the first experiment.}
\end{figure}

The concept for the first experiment can be described by the single relation that a green square is left of a blue square in an image $x$. Figure~\ref{fig:first_example} shows two positive examples ($a$, $b$) and one negative example ($c$).
After the full 10 epochs, the accuracy on the validation set reached \texttt{93.47\%}. For image \ref{fig:first_example}.a the classifier gave the estimator for the concept to be \texttt{89.83\%}. For image \ref{fig:first_example}.b the estimator was \texttt{94.18\%}. The estimator output for belonging to the concept for \ref{fig:first_example}.c was \texttt{0.28\%} showing that the network is able to discriminate the positive and negative examples.
To generate explanations for these three images, each image is separately fed into LIME. The number of kept super-pixels $k$ is set to 3. We choose this value for $k$ because we were aware that there are 2 squares in the image that are distinguishable from the background. One additional super-pixel was taken to generate a richer pool for selection $S$ containing also some background. In general, for many domains it is not that easy to estimate good values for $k$. So in most of the cases it is preferable to over-estimate the value to not lose information for the explanation.

Finally, a symbolic explanation is generated with LIME-Aleph. We describe the procedure for the positive example \ref{fig:first_example}.a. First, Algorithm~\ref{alg:lime-aleph} extracts the colors of the selected super-pixels and the relations between them. Then, all the relations get flipped one after the other to produce the example set and BK for Aleph.
The original example \ref{fig:first_example}.a is used as the seed for a set of perturbed versions of the image. Threshold $\theta$ indicates whether the perturbed example is classified positive or negative. Based on the final validation accuracy of the trained $f$ and from the estimator for the original image \ref{fig:first_example}.a, it was set to $\theta = 0.8$.
For example \ref{fig:first_example}.a 3 positive and 4 negative examples were created. From these 7 examples, Aleph induced a theory $T$ consisting of a single rule with accuracy of \texttt{100\%}:

\begin{center}
	\texttt{concept(A) :- contains(B,A), has\_color(B,green), contains(C,A), has\_color(C,blue), left\_of(B,C).}
\end{center}

\noindent The learned rule accurately resembles the construction regulation of the wanted concept; a green square has to be left of a blue square in an example A. Also, this explanation matches the input image.

For image \ref{fig:first_example}.b we used the same hyper-parameters ($k=3$, $\theta=0.8$). Again, Aleph came up with an accuracy of \texttt{100\%} and a rule structurally different, but conveying the same concept as the first rule:

\begin{center}
	\texttt{concept(A) :- contains(B,A), has\_color(B, blue), contains(C,A), has\_color(C, green), left\_of(C,B).}
\end{center}

\subsection{Experiment 2: Tower Concept}\label{complex_tower}

In the second experiment, we investigated a specific concept of towers. Positive examples consist of three differently colored blocks with a given restriction on their stacking order. An example belongs to the concept \texttt{tower}, if a blue square is present as a foundation. Directly \texttt{on} the foundation (one grid cell above) there has to be a square of either cyan or green color. Directly \texttt{on} that square has to be the remaining square (green or cyan). Figure~\ref{fig:tower_experiments} gives a positive and two negative examples.

\begin{figure}[t]
	\centering
	\subfigure[]{\label{fig:second_ex_a_pos}\includegraphics[width=20mm]{second_a_pos.png}}
	\quad\quad\quad
	\subfigure[]{\label{fig:second_ex_b_neg}\includegraphics[width=20mm]{second_b_neg.png}}
	\quad\quad\quad
	\subfigure[]{\label{fig:second_ex_c_neg}\includegraphics[width=20mm]{second_c_neg.png}}
	\caption{\label{fig:tower_experiments} Positive ($a$) and negative ($b$, $c$)  for the tower experiment.}
\end{figure}

We again trained the CNN for 10 epochs. The final validation accuracy was \texttt{98.70\%}. The original estimator for example \ref{fig:tower_experiments}.a gave $f(a) = 94.88\%$. We first set $k=3$ being the smallest selection of which we know can contain a tower. Again setting $\theta=0.8$, LIME-Aleph came up with 5 positive and 6 negative examples (accuracy \texttt{81.82\%}) and the following rule:

\begin{center}
	\texttt{concept(A) :- contains(B,A), has\_color(B,cyan), contains(C,A), on(B,C).}
\end{center}

\noindent This rule expresses the fact, that the cyan square can not be the foundation.
When setting the selection size $k=4$, we let an additional background super-pixel be part of $S$. The resulting rule is:

\begin{center}
	\texttt{concept(A) :- contains(B,A), has\_color(B,cyan), contains(C,A), has\_color(C,blue), top\_of(B,C).}
\end{center}

\noindent This rule captures the fact, that a cyan block has to be above a blue one. The generated explanations are only partial representations of the intended concept. The symbolic explanations capture relevant aspects, but are too general.

\section{Bringing Together Visual and Symbolic Explanations}\label{fusion}

The generated rules give explanations in symbolic form which can be re-written into verbal statements. We postulate that helpful explanations for images should relate highlighting of relevant parts of the image with explicit symbolic information of attributes and relations. In this section we give an example on how this fusion might look like. Let us take the tower example from Section~\ref{complex_tower}. In Figure~\ref{fig:fusion_explanation}, the output of standard LIME is given with the 3 most important super-pixels matching the expected region in the image. Additionally, the relation from the instantiated rule from the experiment for $k=3$ is given. Since cyan is the only square that is mentioned in the rule, we take it as a reference. The relation \texttt{on} links the cyan square to another unknown square below. This relation is shown explicitly in the image by connecting the two squares and writing the instantiated relation.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.4\textwidth]{fusion_explanation.png}
	\caption{An example for the combination of visual and verbal explanations. Here it is explained, why and where this particular image shows evidence for belonging to the concept tower.} \label{fig:fusion_explanation}
\end{figure}

\section{Conclusion and Further Work}\label{conclusion}

We proposed an approach to extract symbolic rules from images which can be used to explain classifier decisions in a more expressive way than visual highlighting alone. For a simple artificial domain we gave a proof of concept for our approach LIME-Aleph. The work presented here significantly extends \cite{RaboldSS18} by providing a method of automated extraction of visual attributes and spatial relations from images. As a next step we want to also let the explanative power be evaluated by humans.
Also we plan to cover real world image datasets like explaining differences between towers and houses as shown in Figure~\ref{fig:visualverbal}. The challenge here is to come up with arbitrarily placeable segmentations that are easily interchangeable. While our algorithm relies on a regular grid, in an image ``in the wild'', the semantic borders of sub-objects can be irregular in shape and not easily be flipped in order to test for different relations. One idea to cope with this problems is to use relevance information from inner layers in a CNN (e.g., with LRP,~\cite{samek2017explainable}) to first pinpoint small important regions and sub-objects, then super-imposing a standardized selection shape (square, circle, etc.) over the pixel values to find interchangeable super-pixels for filling selection $S$.

In general, it might be useful to consider a variety of explanation formats to accommodate specific personal preferences and situational contexts. For example, visual highlighting is a quick way to communicate what is important while verbal explanations convey more details. Likewise, examples which prototypically represent a class and near-miss counter-examples could be used to make system decisions more transparent~\cite{adadi2018peeking}.
Explanations might also not be a one-way street. In many domains, it is an illusion that the labeling of the training  is really a ground truth. For example, in medical diagnosis, there are many cases where not even experts agree. Therefore, for many practical applications, learning should be interactive \cite{fails2003interactive}. To constrain model adaption, the user could mark-up that parts of an explanation which are irrelevant or wrong. Such a cooperative approach might improve the joint performance of the human-machine-partnership.

\bibliographystyle{splncs04}


text:
\section{Introduction}

The Dedekind eta-function, defined as:
\begin{displaymath}
\eta(z):=q^{\frac{1}{24}}\prod_{n=1}^{\infty}(1-q^n)
\end{displaymath}
where $q=2\pi i z$, plays an important role in modular form theories. Moreover, as the generating function of the Partition Function $p(n)$:
\begin{displaymath}
\sum_{n=0}^{\infty}p(n)q^n=\prod_{n=1}^{\infty}(1-q^n)^{-1}
\end{displaymath}
$\eta(z)$ is also very useful in partition theory and combinatorial theory. So are the eta-quotients, defined as the products and quotients of some eta functions:
\begin{displaymath}
f(z):=\prod_{\delta}\eta(\delta z)^{r_{\delta}}=q^s\prod_{n=1}^{\infty}\prod_{\delta}(1-q^{\delta n})^{r_{\delta}}
\end{displaymath}
where $\delta$ are positive integers, and $r_{\delta}$ ,related to $\delta$, are non-zero integers. $s:=\frac{1}{24}\sum_{\delta}\delta r_{\delta}$. When all $r_{\delta}$ are positive, $f(z)$ is also called eta-products.

An interesting property of a series $\sum_{n=0}^{\infty}a(n)q^n$ is the "lacunarity". A series $\sum_{n=0}^{\infty}a(n)q^n$ is said to be $lacunary$ if "almost all" of its
coefficients are zero, that is:
\begin{displaymath}
\lim_{X\rightarrow\infty}\frac{\sharp\{0\leq n \leq X: a(n)=0\}}{X}=1
\end{displaymath}
By far, there have been many researches on the lacunarity of eta-quotients. For example
in [1], Serre proved that given a positive even integer r, $\eta(z)^r$ is lacunary if and only if $r\in\{2,4,6,8,10,14,26\}$. In [2], Clader discussed the eta-quotients in the
form $\frac{\eta(az)^b}{\eta(z)}$, and proved the only 19 cases when they are lacunary. In [3], Martin researched the lacunarity of eta-quotients which are Hecke eigenforms. In this
article, we discuss the eta-products in the form $\eta(z)^2\eta(kz)^2$, and prove the following theorem:

\textbf{Theorem 1}  Let $f_b(z):=\eta(z)^2\eta(bz)^2$, $b$ is a positive integer and has no square integer divisor except 2,3. $23\nmid b$. Then $f_b(z)$ is lacunary if and only if $b\in\{1,2,3,4,16\}$.

We will prove this theorem in section 3, using the theory of modular form and CM forms. After that, we will give the method of computing the coefficients of $f_b(z)$ when it is lacunary
 in section 4.
~\\

\section{Preliminaries}
~\\

The Dedekind eta-function $\eta(24z)$, as we know, is a weight $\frac{1}{2}$ modular form on $\Gamma_0(576)$, with Nebentypus character $\chi_{12}(n):=\left(\frac{12}{n}\right)$  . As consequence, the eta-quotient may be modular form on some modular group.
In fact, we have the following theorem from Ono[4]:

\textbf{Theorem 2.1} If $f(z)=\prod_{\delta |N}\eta(\delta z)^{r_{\delta}}$ is an eta-quotient with $k=\frac{1}{2}\sum_{\delta | N}r_{\delta}$, with the additional properties that
\begin{displaymath}
\sum_{\delta | N}\delta r_{\delta}\equiv 0~(mod~24)
\end{displaymath}
and
\begin{displaymath}
\sum_{\delta | N}\frac{N}{\delta} r_{\delta}\equiv 0~(mod~24)
\end{displaymath}
then $f(z)$ satisfies
\begin{displaymath}
f\left(\frac{az+b}{cz+d}\right)=\chi(d)(cz+d)^k f(z)
\end{displaymath}
for every $\left(\begin{array}{ccc} a&b\\c&d \end{array}\right)\in \Gamma_0(N)$. Here the character $\chi$ is defined by $\chi(d):=\left(\frac{(-1)^k\cdot s}{d}\right)$, where
$s:=\prod_{\delta | N}\delta^{r_{\delta}}$.
Moreover, if $f(z)$ is holomorphic (resp. vanishes) at all of the cusps of $\Gamma_0(N)$, then $f(z)\in M_k(\Gamma_0(N),\chi)$ (resp. $S_k(\Gamma_0(N),\chi)$).

After theorem 2.1, we can easily get the following lemma about $f_b(z)$:

\textbf{Lemma 2.2}  Let $f_b(12z):=\eta(12z)^2\eta(12bz)^2$ as above. Then $f_b(12z)\in S_2(\Gamma_0(144b))$.

We still need the theory of CM forms. The basic idea of CM form can be found in Ribet[5]. Let K be an imaginary quadratic field of discriminant D, and $\varepsilon_K$ be the quadratic character associated to K. Next, let $c$ be a Hecke character of $K$ with exponent $k-1$ and conductor $f_c$. That means, view $c$ as a homomorphism:
\begin{displaymath}
(fractional~ideals~of~K~prime~to~f_c) \rightarrow C^{\ast}
\end{displaymath}
we have
\begin{displaymath}
c((\alpha))=\alpha^{k-1}
\end{displaymath}
for all $\alpha\in F^{\ast}$ such that $\alpha \equiv 1~(mod^{\ast}~f_c)$.

Next, associated to $c$, define a Dirichlet character $\omega_c$ mod $N(f_c)$ as:
\begin{displaymath}
\omega_c(n)=c((n))/n^{k-1}
\end{displaymath}
for all $n\in \mathbb{Z}$ coprime to $f_c$. $N(f_c)$ denotes the norm of $f_c$.

For any $\delta\in \mathbb{Z}^+$, define series $\varphi_{k,c,\delta}$ as:
\begin{displaymath}
\varphi_{K,c,\delta}=\sum_{a}c(a)q^{\delta\cdot N(a)}
\end{displaymath}
where the sum runs through all integral ideals of $K$ coprime to $f_c$. Then, by theorem 3.4 in [5], $\varphi_{K,c,\delta}$ is a cusp form of weight $k$ and
character $\omega_c\cdot \varepsilon_K$ on $\Gamma_0(\delta\cdot|D|\cdot N(f_c))$, moreover, $\varphi_{K,c,\delta}$ is an eigenform for all Hecke operators
$T_p$ in which $p\nmid \delta\cdot|D|\cdot N(f_c)$. So, in order to make $\varphi_{K,c,\delta}$ to be an element in $S_k(\Gamma_0(N),\varepsilon)$, the following two
conditions are necessary and sufficient:
\begin{displaymath}
\delta\cdot|D|\cdot N(f_c)~|~N
\end{displaymath}
\begin{displaymath}
\omega_c\cdot \varepsilon_K=\varepsilon
\end{displaymath}
All possible $\varphi_{K,c,\delta}$ satisfying the above two conditions generate a subspace of $S_k(\Gamma_0(N),\varepsilon)$, which denoted as $S_k^{cm}(\Gamma_0(N),\varepsilon)$.
Then, by Serre[6], an element of $S_k(\Gamma_0(N),\varepsilon)$ is lacunary if and only if it is an element of $S_k^{cm}(\Gamma_0(N),\varepsilon)$.
~\\

\section{Proof of Theorem 1}
~\\

With the above theory, we know that $f_b(12z)$ is lacunary if and only if $f_b(12z)=\sum \varphi_i \in S_2^{cm}(\Gamma_0(144b),\varepsilon)$, where $\varphi_i$ are some suitable
$\varphi_{K,c,\delta}$. First we have the following lemmas:

\textbf{Lemma 3.1} $f_b(12z)=\sum \varphi_i$ as above. $b$ as in Theorem 1. Then as a CM form, $\varphi_i$ can only be associated to one of the following four fields: $\mathbb{Q}(\sqrt{-1})$,
$\mathbb{Q}(\sqrt{-2})$, $\mathbb{Q}(\sqrt{-3})$ and $\mathbb{Q}(\sqrt{-6})$.

\textbf{Proof} Assume $\varphi_i$ is generated by $\mathbb{Q}(\sqrt{D})$, $k$, $f_c$ and $\varepsilon_K$, $\omega_c$. So we have $|D|\cdot N(f_c)~|144b$, and $\omega_c\cdot \varepsilon_K=1_N$,
where $1_N$ is the trivial character on $\Gamma_0(144b)$. By $\omega_c\cdot \varepsilon_K=1_N$, if there is some prime $p\neq 2,3$ such that $p\mid|D|$, that is $p$ is a divisor of the module of character $\varepsilon_K$, so p must be a divisor of the module of character $\omega_c$, too. But we already know that $\omega_c$ is a Dirichlet character mod $N(f_c)$, so there must be
$p\mid N(f_c)$, and $p^2\mid |D|\cdot N(f_c)~|144b$. This is a contradiction to our assumption that b is a square-free integer. This completes the proof. $\Box$

\textbf{Lemma 3.2} $\varphi_i$ is a CM form associated to $\mathbb{Q}(\sqrt{-1})$, $\mathbb{Q}(\sqrt{-2})$, $\mathbb{Q}(\sqrt{-3})$ or $\mathbb{Q}(\sqrt{-6})$. $T_{23}$ is the usual Hecke operator $T_p$ with $p=23$. Then $\varphi_i|T_{23}=0$.

This lemma is obvious because 23 remains prime in each of those four fields, so the 23-rd coefficient of $\varphi_i$ must be zero. Combined with the fact that $\varphi_i$ is an eigenform of
$T_{23}$, and the computing formula of $T_p$: $a_n(T_p(f))=a_{np}(f)+\chi(p)p^{k-1}a_{n/p}(f)$, this lemma can be easily obtained. $\Box$

\textbf{Lemma 3.3}  $f_b(12z)=\sum \varphi_i$ as in Theorem 1. Then if $b\geq 176$, $f_b(12z)$ cannot be lacunary.

\textbf{Proof} By Lemma 3.2, if $f_b(12z)$ is lacunary, then $f_b(12z)|T_{23}=0$. We will show that if $b$ is large enough, $f_b(12z)|T_{23}$ cannot be zero. Denote
$f_b(12z)=q^{1+b}\prod_{n=1}^{\infty}(1-q^{12n})^2(1-q^{12bn})^2=\sum_{n=1}^{\infty}a(n)q^n$ and $\prod_{n=1}^{\infty}(1-q^{n})^2=\sum_{n=1}^{\infty}b(n)q^n$. If $n<12b$, then $a(n)$ has nothing to do with $(1-q^{12bn})^2$, so we have $a(n)=b(\frac{n-(1+b)}{12})$. Denote $T_{23}(f)=\sum_{n=1}^{\infty}c(n)q^n$, Using the formula of $T_p$, we have $c(n)=a(23n)+p\cdot a(n/23)$. When $23\nmid n$, $c(n)=a(23n)=b(\frac{23n-(1+b)}{12})$. In order to make $\frac{23n-(1+b)}{12}$ an integer, we need $n\equiv -(1+b)~(mod~12)$.
If there exists some $c(n)\neq 0$, then $f_b(12z)|T_{23}\neq 0$. In summary, after three conditions: (1) $1+b\leq n<12b$; (2) $23\nmid n$; (3) $n\equiv -(1+b)~(mod~12)$; If there exists such $n$ that $b(\frac{23n-(1+b)}{12})\neq 0$, then $f_b(12z)$ cannot be lacunary.

Let $n_0$ be the smallest $n$ that makes $\frac{23n-(1+b)}{12}$ an positive integer. So $\frac{23n-(1+b)}{12}\in \{1,2,3,\ldots,23\}$, and when $n_0\rightarrow n_0+12$,
$\frac{23n-(1+b)}{12}\rightarrow \frac{23n-(1+b)}{12}+23$. Obviously, in $n_0+12t$, $t=0,1,2,\ldots,22$, these 23 numbers, only one can be divided by 23. We list the first 1000 coefficients of $b(n)$ in appendix 1.

From the table of $b(n)$ above, we only need $t\geq 5$, then there will be at least two numbers not equal to zero in each column, except the 21st column, which we will discuss later. So there must be at least one number that suits $23\nmid n$. In order to suit $1+b\leq n<12b$, we need $b>\frac{12*7*23}{11}$, that is $b\geq 176$. As to the 21st column, in this
condition, $\frac{23n_0-(1+b)}{12}=21$, that is $21n_0-(1+b)=21*12$. Taking module of 23 on the both sides, we get $k\equiv 0~(mod~23)$, which contradicts our assumption in Theorem 1. In summary, we can say that when $b\geq 176$, $f_b(12z)$ cannot be lacunary. $\Box$

Now, in order to prove Theorem 1, we only need to check the lacunarity of each case of $b<176$ by computer. The method is to check whether $f_b(12z)|T_{23}$ (or $T_{47}$, when $23|b$) is zero. At last we find that $f_b(12z)$ cannot be lacunary except $b=1,2,3,4$ or $16$. We will prove those five cases in Proposition 4.1. Theorem 1 is proved.
~\\

\section{Computation when $b=1,2,3,4,16$}
~\\

When b=1,2,3,4,16, $f_b(z)$ is lacunary. So it can be expressed as a linear combination of CM forms. We list these specific expressions in Proposition 4.1:

\textbf{Proposition 4.1}

(1) $b=1$; $f_b(6z)=\eta(6z)^4=\varphi_{K,c}(z)\in S_2(\Gamma_0(36))$, where $K=\mathbb{Q}(\sqrt{-3})$, $f_c=(2\sqrt{-3})$. This can be found in [1].

(2) $b=2$; $f_b(4z)=\eta(4z)^2\eta(8z)^2=\varphi_{K,c}(z)\in S_2(\Gamma_0(32))$, where $K=\mathbb{Q}(\sqrt{-1})$, $f_c=(2(1+i))$.

(3) $b=3$; $f_b(3z)=\eta(3z)^2\eta(9z)^2=\varphi_{K,c}(z)\in S_2(\Gamma_0(27))$, where $K=\mathbb{Q}(\sqrt{-3})$, $f_c=(3)$.

(4) $b=4$; $f_b(12z)=\eta(12z)^2\eta(48z)^2=\frac{1}{8}(\varphi_{K,c_-}(z)-\varphi_{K,c_+}(z))\in S_2(\Gamma_0(576))$, where $K=\mathbb{Q}(\sqrt{-1})$, $f_c=(12)$.

(5) $b=16$; $f_b(12z)=\eta(12z)^2\eta(16\cdot12z)^2=\frac{1}{16}(\varphi_{603}-\varphi_{203})-\frac{\sqrt{1-2\sqrt{-6}}}{16\cdot(6-2\sqrt{-6})}(\varphi_{130}-\varphi_{130}'+\varphi_{310}-\varphi_{310}')
\in S_2(\Gamma_0(144\cdot 16))$, where $\varphi_{603}$ and $\varphi_{203}$ are associated to $K=\mathbb{Q}(\sqrt{-1})$, $f_c=(24)$, and $\varphi_{130}$, $\varphi_{130}'$, $\varphi_{310}$, $\varphi_{310}'$ are associated to $K=\mathbb{Q}(\sqrt{-6})$, $f_c=(4\sqrt{-6})$.

\textbf{Proof}  First we need to compute the specific Hecke character on each case. On case (2), the hecke character $c$ is defined as follows:
choose $a+bi$ to be primary, that is, if $b\equiv 0~(mod~4)$, then $a\equiv 1~(mod~4)$; if $b\equiv 2~(mod~4)$, then $a\equiv 3~(mod~4)$. Then $c((a+bi))=a+bi$. On case (3),
the hecke character $c$ is defined as follows: choose $a+b\omega$ to be primary, that is, $a\equiv 2~(mod~3)$, $b\equiv 0~(mod~3)$. Then $c((a+b\omega))=-(a+b\omega)$. On case (4),
the hecke character $c_+$ and $c_-$ is defined as follows: choose $a+bi$ to be primary, $a+bi\equiv (1-i)^u~(mod~3)$, $u=0,1,\ldots,7$, and $a+bi\equiv (-1+2i)^v~(mod~4)$, $v=0,1$. Then $c_+((a+bi))=i^u\cdot(-1)^v\cdot(a+bi)$, $c_-((a+bi))=(-i)^u\cdot(-1)^v\cdot(a+bi)$. Case (5) is a little bit complicated. To $\varphi_{603}$ and $\varphi_{203}$, the hecke character $c_{603}$ and $c_{203}$ are defined as: choose $a+bi$ to be "standard", which means $a+bi\equiv 1,2+i,3+4i,2+3i,4+i,-1-2i,3i,or~-3-2i~(mod~8)$, if $a+bi\equiv (1-i)^u~(mod~3)$, $u=0,1,\ldots,7$, and $a+bi\equiv (2+i)^v(4+i)^w~(mod~8)$, $v=0,1,2,3$, $w=0,1,2,3$. Then $c_{rst}(a+bi)=\zeta_8^{ru}\cdot i^{sv} \cdot i^{tw}\cdot(a+bi)$, $\zeta_8=e^{\frac{2\pi i}{8}}$. To $\varphi_{130}$, $\varphi_{130}'$, $\varphi_{310}$, $\varphi_{310}'$, the hecke character $c_{130}$, $c_{130}'$, $c_{310}$ and $c_{310}'$ are defined as:
choose $a+b\sqrt{-6}$ to be "standard", which means $a+b\sqrt{-6}\equiv 1,1+\sqrt{-6},1-\sqrt{-6},5,or~the~product~of~them~(mod~4\sqrt{-6})$,
if $a+b\sqrt{-6}\equiv (1+\sqrt{-6})^u\cdot(1-\sqrt{-6})^v\cdot 5^w$, $u,v=0,1,2,3$, $w=0,1$. Then $c_{rst}(a+b\sqrt{-6})=i^{ru}\cdot i^{sv}\cdot (-1)^{tw}\cdot (a+b\sqrt{-6})$. However,
$\mathbb{Q}(\sqrt{-6})$ is not a Principal Ideal Domain, that means not every ideal in $\mathbb{Q}(\sqrt{-6})$ is principal. Let $5=(5,2+\sqrt{-6})\cdot(5,2-\sqrt{-6})$, and denote
$(5,2+\sqrt{-6})$ as $\alpha$. Every non-principal ideal plus $\alpha$ should be a principal ideal because $\mathbb{Q}(\sqrt{-6})$ has Class Number 2. Let $c(\alpha)=\sqrt{c(\alpha^2)}=\sqrt{c(-1+2\sqrt{-6})}$, and $c'(\alpha)=-\sqrt{c(-1+2\sqrt{-6})}$.
That would make the definition explicit.

Next we will verify the above result by Sturm's Theorem[7]:

\textbf{Sturm's Theorem}    Suppose that N is a positive integer, $p$ is prime, and $f(z)$, $g(z)\in M_k(\Gamma_0(N),\chi)\cap \mathbb{Z}[[q]]$. If
\begin{displaymath}
ord_p(f(z)-g(z))>\frac{k}{12}[SL_2(\mathbb{Z}):\Gamma_0(N)]
\end{displaymath}
in which
\begin{displaymath}
[SL_2(\mathbb{Z}):\Gamma_0(N)]=N \prod_{\ell~prime:\ell |N}(1+\frac{1}{\ell})
\end{displaymath}
then $f(z)\equiv g(z)~(mod~p)$.

Case (1)(2)(3) is trivial, because they are in one-dimension linear spaces. To case (4) and (5), although it is very complicated to get the above result, it is relatively easy to check it. By Sturm's Theorem, we only need to verify the first 192 and 768 coefficients, which are listed in appendix 2 and appendix 3. $\Box$

Here we make an example to show how to compute the coefficients of the case $b=16$. Let
$n=29645=5*7^2*11^2$. Denote $\varphi_{uvw}=\sum a_{uvw}(n)q^n$, and $\varphi_{uvw}'=\sum a_{uvw}'(n)q^n$. $5=(2+i)(-1-2i)$ in $\mathbb{Q}(\sqrt{-1})$, so $a_{603}(5)=2$, $a_{203}(5)=-2$.
By the formula $a(p^r)=a(p)a(p^{r-1})-p\cdot a(p^{r-2})$, we can get $a_{603}(49)=a_{203}(49)=-7$, $a_{603}(121)=a_{203}(121)=-11$, so $a_{603}(n)=2*7*11$, $a_{203}(n)=-2*7*11$.
Next, $5=(5,2+\sqrt{-6})\cdot(5,2-\sqrt{-6})$ in $\mathbb{Q}(\sqrt{-6})$, so $a_{130}(5)=a_{310}(5)=\frac{6-2\sqrt{-6}}{\sqrt{1-2\sqrt{-6}}}$, $a_{130}'(5)=a_{310}'(5)=-\frac{6-2\sqrt{-6}}{\sqrt{1-2\sqrt{-6}}}$. There are three
ideals whose norm are 49: $(-5+2\sqrt{-6})$, $(-5-2\sqrt{-6})$ and (7), so $a_{130}(49)=a_{130}'(49)=a_{310}(49)=a_{310}'(49)=17$. There are also three ideals whose norm are 121: $(5+4\sqrt{-6})$, $(5-4\sqrt{-6})$ and (11), so $a_{130}(121)=a_{130}'(121)=a_{310}(121)=a_{310}'(121)=21$. So $a_{130}(29645)=a_{310}(29645)=17*21*\frac{6-2\sqrt{-6}}{\sqrt{1-2\sqrt{-6}}}$, $a_{130}'(29645)=a_{310}'(29645)=-17*21*\frac{6-2\sqrt{-6}}{\sqrt{1-2\sqrt{-6}}}$. Put all these data in the formula, we get $a(29645)=\frac{1}{4}(7*11-17*21)=-70$.
Verifying it by directly computing with computer, we can see this is the right result.
~\\

~\\

\textbf{Acknowledgement}
~\\

The author is grateful to my tutor Chunlei Liu for the guide throughout the article and the help from Haobo Dai.
~\\

\centerline{References}
~\\

[1] J.P.Serre, Sur la lacunarite des puissances de eta, Glasgow Math. J.27(1985). 203-221

[2] E.Clader, Y.Kemper, M.Wage, Lacunarity of certain partition theoretic generating functions, Proceedings of the American Mathematic Society. Volume 137(9),2009. 2959-2968

[3] Y.Martin, Multiplicative eta-quotients, Trans. Amer. Math. Soc. 348(1996). 4825-4856

[4] K.Ono, The Web of Modularity: Arithmetic of the Coefficients of Modular Forms and q-series, CMBS Regional Conference Series in Mathematics, American Mathematical Society,
Providence, RI, 2004

[5] K.A.Ribet, Galois representations attached to eigenforms with Nebentypus

[6] J.P.Serre, Quelques applications du theoreme de densite de Chebatorev, Publ. Math. I.H.E.S.54 (1981). 123-201

[7] John J. Webb, Arithmetic of the 13-regular partition function modulo 3, Ramanujan J(2011)25:49-56

\newpage

\centerline{\textbf{Appendix 1}}

This is the first 1000 coefficients of $\prod_{n=1}^{\infty}(1-q^{n})^2=\sum_{n=1}^{\infty}b(n)q^n$. The first row, from left to right, are $b(1)$, $b(2)$, $\ldots$, $b(23)$, the second row are  $b(24)$, $b(25)$, $\ldots$, $b(46)$, and so on. We can see that in the first six rows, except for the 21st column, there are at least two numbers not equal to zero in each column. That's just what we need in Lemma 3.3.

\begin{displaymath}
\begin{tabular}{*{23}{c}}
-2&-1&2&1&2&-2&0&-2&-2&1&0&0&2&3&-2&2&0&0&-2&-2&0&0&-2\\
-1&0&2&2&-2&2&1&2&0&2&-2&-2&2&0&-2&0&-4&0&0&0&1&-2&0\\
0&2&0&2&2&1&-2&0&2&2&0&0&-2&0&-2&0&-2&2&0&-4&0&0&-2\\
-1&2&0&2&0&0&0&-2&2&4&1&0&0&2&-2&2&-2&0&0&2&0&-2&0\\
-2&-2&0&-2&0&0&0&2&-2&-1&-2&-2&4&0&0&2&0&2&0&0&0&3&-2\\
0&0&4&2&0&-2&0&0&-2&0&0&-2&0&2&0&-2&0&-2&-2&-2&0&0&2\\
-2&-1&2&0&0&0&2&-2&0&-2&2&2&2&2&0&1&2&-2&0&0&0&0&2\\
0&0&0&0&2&0&-2&-2&-4&2&0&0&-2&0&-2&0&2&0&-2&0&0&-4&1\\
-2&0&0&4&2&-2&2&0&0&0&-2&4&0&-2&2&1&0&2&2&0&0&2&0\\
0&-4&2&0&0&-2&0&0&2&0&-2&0&0&0&0&-2&-2&-4&-2&2&0&2&0\\
0&0&-2&-1&0&2&0&-2&0&0&0&0&2&0&0&2&0&4&2&-2&0&-1&2\\
-2&2&0&0&0&2&-2&4&0&0&2&-2&0&0&-2&-2&0&-2&0&0&0&2&-2\\
0&0&0&-2&-2&0&0&0&0&-2&0&-2&-2&1&0&0&2&2&2&0&0&-2&-4\\
4&2&0&-2&0&0&2&0&2&2&3&-2&0&2&2&0&-2&0&0&0&0&-2&0\\
2&2&0&0&-2&0&0&0&0&0&0&0&-2&-4&0&2&-4&0&-2&0&0&2&0\\
2&0&-2&0&-2&0&-3&0&0&2&-2&0&2&0&0&0&0&2&0&4&0&0&0\\
0&2&-4&0&2&1&0&2&2&-4&2&2&0&0&2&0&2&0&0&-2&0&0&-2\\
0&0&-2&0&-2&0&0&2&2&2&-2&0&-4&-2&0&0&0&-2&0&-2&0&-2&2\\
0&-2&0&0&0&1&0&0&2&0&-2&2&0&0&0&-4&0&0&2&2&0&2&0\\
0&0&2&0&0&4&3&-2&0&0&0&0&0&-2&0&-2&2&0&4&0&0&0&2\\
0&0&-2&-2&2&0&0&0&-4&-2&2&0&0&-2&-2&0&2&-2&-2&0&0&0&0\\
4&-2&0&0&-2&-2&-2&0&0&-4&1&4&-2&0&0&0&0&0&2&2&0&0&2\\
0&0&0&-2&2&0&0&0&0&0&0&-2&2&1&6&0&2&0&-2&0&0&2&0\\
-2&2&0&2&0&-2&0&0&0&-2&-2&0&0&0&2&0&2&0&-2&0&0&0&0\\
2&-4&-2&-2&0&-4&2&0&-2&0&0&0&-2&2&0&0&2&-2&0&0&-2&1&0\\
0&2&0&-2&-2&0&0&-2&0&0&4&0&2&-2&4&0&0&0&-2&4&0&0&2\\
0&0&2&1&-2&0&0&0&0&-2&2&2&-2&0&0&2&2&2&0&2&0&0&0\\
0&-2&-2&-4&0&0&2&2&-2&-2&0&0&0&-2&0&0&0&-2&2&2&0&-2&0\\
-2&-2&0&0&2&0&-4&-2&0&0&0&-2&0&0&0&-1&2&0&4&-4&0&2&0\\
2&0&0&-2&0&0&2&2&0&-2&0&0&-2&0&0&0&2&2&2&0&0&2&3\\
2&0&2&-2&0&-2&0&-2&2&0&0&2&0&0&2&-4&0&0&0&2&0&0&0\\
0&2&-4&0&0&2&2&0&0&-2&-2&0&-2&-2&2&-4&-2&0&0&0&0&2&-2\\
0&0&-2&0&2&0&-4&2&2&0&0&0&0&-2&-2&-1&0&-2&0&0&0&0&2\\
0&-2&0&2&0&2&0&0&4&0&2&0&0&-2&0&0&0&-2&2&2&0&0&0\\
2&3&-2&-2&0&2&0&0&0&0&0&4&0&0&0&4&0&2&-2&0&0&-2&-2\\
0&0&-2&0&-2&2&0&0&-2&2&0&0&-2&-2&4&0&0&0&0&-2&0&-2&0\\
2&0&-2&0&-4&-2&0&0&0&0&-2&0&-2&0&0&0&2&0&0&2&0&-1&0\\
2&0&-4&0&-2&0&0&6&2&-2&0&-2&-2&0&0&0&0&2&0&2&0&-2&0\\
\end{tabular}
\end{displaymath}
\begin{displaymath}
\begin{tabular}{*{23}{c}}
2&2&2&0&0&0&4&0&2&1&0&0&0&-2&0&0&-2&0&0&-2&0&2&0\\
4&2&0&0&-2&2&0&-2&0&0&0&2&0&2&0&-2&0&0&-2&0&0&0&-4\\
0&0&0&-4&0&2&0&0&-2&-2&2&0&2&0&-2&0&0&-2&0&0&0&0&-2\\
2&0&-2&0&0&0&-4&0&2&-2&1&-4&0&0&2&0&0&0&0&0&0&4&2\\
0&2&0&0&0&0&-2&2&-2&0&0&-2&0&0&4&0&0&-2&2&2&0&0&0\\
3&0&0&0&0&2&2&0&0&0&0&
\end{tabular}
\end{displaymath}

\newpage

\centerline{\textbf{Appendix 2}}

Let $f_4(12z)=\sum a(n)q^n$, and $\varphi_{K,c_+}(z)=\sum b(n)q^n$, $\varphi_{K,c_-}(z)=\sum c(n)q^n$. By proposition 4.1(4), we only need to verify $a(n)=\frac{1}{8}(c(n)-b(n))$.

We omit those $n$ that $2|n$ or $3|n$, because those coefficients are all zero.
\begin{displaymath}
\begin{array}{cccc}
\quad n\quad&\quad a(n)\quad&\quad b(n)\quad&\quad c(n)\quad\\
1&0&1&1\\
5&1&-4&4\\
7&0&0&0\\
11&0&0&0\\
13&0&6&6\\
17&-2&8&-8\\
19&0&0&0\\
23&0&0&0\\
25&0&11&11\\
29&-1&4&-4\\
31&0&0&0\\
35&0&0&0\\
37&0&2&2\\
41&2&-8&8\\
43&0&0&0\\
47&0&0&0\\
49&0&-7&-7\\
53&-1&4&-4\\
55&0&0&0\\
59&0&0&0\\
61&0&10&10\\
65&6&-24&24\\
67&0&0&0\\
71&0&0&0\\
73&0&6&6\\
77&0&0&0\\
79&0&0&0\\
83&0&0&0\\
85&0&-32&-32\\
89&-4&16&-16\\
91&0&0&0\\
95&0&0&0\\
97&0&-18&-18\\
101&-5&20&-20\\
103&0&0&0\\
107&0&0&0\\
109&0&-18&-18\\
113&-4&16&-16\\
\end{array}
\end{displaymath}
\begin{displaymath}
\begin{array}{cccc}
\quad n\quad&\quad a(n)\quad&\quad b(n)\quad&\quad c(n)\quad\\
115&0&0&0\\
119&0&0&0\\
121&0&-11&-11\\
125&6&-24&24\\
127&0&0&0\\
131&0&0&0\\
133&0&0&0\\
137&-2&8&-8\\
139&0&0&0\\
143&0&0&0\\
145&0&-16&-16\\
149&5&-20&20\\
151&0&0&0\\
155&0&0&0\\
157&0&-22&-22\\
161&0&0&0\\
163&0&0&0\\
167&0&0&0\\
169&0&23&23\\
173&1&-4&4\\
175&0&0&0\\
179&0&0&0\\
181&0&-18&-18\\
185&2&-8&8\\
187&0&0&0\\
191&0&0&0\\
\end{array}
\end{displaymath}

\newpage

\centerline{\textbf{Appendix 3}}

Let $f_{16}(12z)=\sum a(n)q^n$, and $\varphi_{603}(z)=\sum b_1(n)q^n$, $\varphi_{203}(z)=\sum b_2(n)q^n$, $\frac{\sqrt{1-2\sqrt{-6}}}{6-2\sqrt{-6}}\varphi_{130}(z)=\sum c_1(n)q^n$,
$\frac{\sqrt{1-2\sqrt{-6}}}{6-2\sqrt{-6}}\varphi_{130}'(z)=\sum c_2(n)q^n$, $\frac{\sqrt{1-2\sqrt{-6}}}{6-2\sqrt{-6}}\varphi_{310}(z)=\sum c_3(n)q^n$,
$\frac{\sqrt{1-2\sqrt{-6}}}{6-2\sqrt{-6}}\varphi_{310}'(z)=\sum c_4(n)q^n$.
Denote $\frac{\sqrt{1-2\sqrt{-6}}}{6-2\sqrt{-6}}=t$. By proposition 4.1(5), we need to verify that $a(n)=\frac{1}{16}\left(b_1(n)-b_2(n)\right)-\frac{1}{16}\left(c_1(n)-c_2(n)+c_3(n)-c_4(n)\right)$.

We omit those $n$ that $2|n$ or $3|n$, because those coefficients are all zero.

"$\ast$" means that because of the "minus" in the formula, "$*-*=0$", so we don't need to compute these specific values. These case are:

(1) If $p=(a+bi)$ is an ideal in $\mathbb{Q}(\sqrt{-1})$, $3|a$ or $3|b$, then $c_+(p)=c_-(p)$.

(2) If $p$ is a principle ideal in $\mathbb{Q}(\sqrt{-6})$, then $c_{130}(p)=c_{130}'(p)$, $c_{310}(p)=c_{310}'(p)$.

\begin{displaymath}
\begin{array}{cccccccc}
\quad n \quad&\quad a(n)\quad & \quad b_1(n) \quad & \quad b_2(n)\quad & \quad c_1(n)\quad & \quad c_2(n) \quad & \quad c_3(n)\quad & \quad c_4(n)\quad\\
1&0&1&1&1&1&1&1\\
5&0&2&-2&1&-1&1&-1\\
7&0&0&0&-2\sqrt{6}t&-2\sqrt{6}t&2\sqrt{6}t&2\sqrt{6}t\\
11&0&0&0&\frac{4}{\sqrt{6}}&-\frac{4}{\sqrt{6}}&-\frac{4}{\sqrt{6}}&\frac{4}{\sqrt{6}}\\
13&0&4&4&0&0&0&0\\
17&1&8&-8&0&0&0&0\\
19&0&0&0&0&0&0&0\\
23&0&0&0&0&0&0&0\\
25&0&-1&-1&7t&7t&7t&7t\\
29&-2&-10&10&3&-3&3&-3\\
31&0&0&0&-2\sqrt{6}t&-2\sqrt{6}t&2\sqrt{6}t&2\sqrt{6}t\\
35&0&0&0&-2\sqrt{6}&2\sqrt{6}&2\sqrt{6}&-2\sqrt{6}\\
37&0&12&12&0&0&0&0\\
41&-1&-8&8&0&0&0&0\\
43&0&0&0&0&0&0&0\\
47&0&0&0&0&0&0&0\\
49&0&-7&-7&17t&17t&17t&17t\\
53&2&14&-14&-1&1&-1&1\\
55&0&0&0&\ast&\ast&\ast&\ast\\
59&0&0&0&\frac{8}{\sqrt{6}}&-\frac{8}{\sqrt{6}}&-\frac{8}{\sqrt{6}}&\frac{8}{\sqrt{6}}\\
61&0&12&12&0&0&0&0\\
65&1&8&-8&0&0&0&0\\
67&0&0&0&0&0&0&0\\
71&0&0&0&0&0&0&0\\
73&0&6&6&14t&14t&14t&14t\\
77&2&0&0&-8&8&-8&8\\
79&0&0&0&6\sqrt{6}t&6\sqrt{6}t&-6\sqrt{6}t&-6\sqrt{6}t\\
83&0&0&0&-\frac{4}{\sqrt{6}}&\frac{4}{\sqrt{6}}&\frac{4}{\sqrt{6}}&-\frac{4}{\sqrt{6}}\\
85&0&16&16&0&0&0&0\\
\end{array}
\end{displaymath}
\begin{displaymath}
\begin{array}{cccccccc}
\quad n \quad&\quad a(n)\quad & \quad b_1(n) \quad & \quad b_2(n)\quad & \quad c_1(n)\quad & \quad c_2(n) \quad & \quad c_3(n)\quad & \quad c_4(n)\quad\\
89&-2&-16&16&0&0&0&0\\
91&0&0&0&0&0&0&0\\
95&0&0&0&0&0&0&0\\
97&0&18&18&2t&2t&2t&2t\\
101&0&-2&2&-1&1&-1&1\\
103&0&0&0&-6\sqrt{6}t&-6\sqrt{6}t&6\sqrt{6}t&6\sqrt{6}t\\
107&0&0&0&\frac{8}{\sqrt{6}}&-\frac{8}{\sqrt{6}}&-\frac{8}{\sqrt{6}}&\frac{8}{\sqrt{6}}\\
109&0&20&20&0&0&0&0\\
113&-2&-16&16&0&0&0&0\\
115&0&0&0&0&0&0&0\\
119&0&0&0&0&0&0&0\\
121&0&-11&-11&21t&21t&21t&21t\\
125&-2&-12&12&2&-2&2&-2\\
127&0&0&0&2\sqrt{6}t&2\sqrt{6}t&-2\sqrt{6}t&-2\sqrt{6}t\\
131&0&0&0&\frac{16}{\sqrt{6}}&-\frac{16}{\sqrt{6}}&-\frac{16}{\sqrt{6}}&\frac{16}{\sqrt{6}}\\
133&0&0&0&0&0&0&0\\
137&1&8&-8&0&0&0&0\\
139&0&0&0&0&0&0&0\\
143&0&0&0&0&0&0&0\\
145&0&-20&-20&36t&36t&36t&36t\\
149&0&-14&14&-7&7&-7&7\\
151&0&0&0&-10\sqrt{6}t&-10\sqrt{6}t&10\sqrt{6}t&10\sqrt{6}t\\
155&0&0&0&-2\sqrt{6}&2\sqrt{6}&2\sqrt{6}&-2\sqrt{6}\\
157&0&12&12&0&0&0&0\\
161&0&0&0&0&0&0&0\\
163&0&0&0&0&0&0&0\\
167&0&0&0&0&0&0&0\\
169&0&3&3&\ast&\ast&\ast&\ast\\
173&2&26&-26&5&-5&5&-5\\
175&0&0&0&\ast&\ast&\ast&\ast\\
179&0&0&0&-\frac{8}{\sqrt{6}}&\frac{8}{\sqrt{6}}&\frac{8}{\sqrt{6}}&-\frac{8}{\sqrt{6}}\\
181&0&20&20&0&0&0&0\\
185&3&24&-24&0&0&0&0\\
187&0&0&0&0&0&0&0\\
191&0&0&0&0&0&0&0\\
193&0&14&14&\ast&\ast&\ast&\ast\\
197&-2&-2&2&7&-7&7&-7\\
199&0&0&0&\ast&\ast&\ast&\ast\\
203&0&0&0&-6\sqrt{6}&6\sqrt{6}&6\sqrt{6}&-6\sqrt{6}\\
205&0&-16&-16&0&0&0&0\\
209&0&0&0&0&0&0&0\\
211&0&0&0&0&0&0&0\\
215&0&0&0&0&0&0&0\\
217&0&0&0&\ast&\ast&\ast&\ast\\
221&4&32&-32&0&0&0&0\\
223&0&0&0&\ast&\ast&\ast&\ast\\
\end{array}
\end{displaymath}
\begin{displaymath}
\begin{array}{cccccccc}
\quad n \quad&\quad a(n)\quad & \quad b_1(n) \quad & \quad b_2(n)\quad & \quad c_1(n)\quad & \quad c_2(n) \quad & \quad c_3(n)\quad & \quad c_4(n)\quad\\
227&0&0&0&\frac{20}{\sqrt{6}}&-\frac{20}{\sqrt{6}}&-\frac{20}{\sqrt{6}}&\frac{20}{\sqrt{6}}\\
229&0&4&4&0&0&0&0\\
233&2&16&-16&0&0&0&0\\
235&0&0&0&0&0&0&0\\
239&0&0&0&0&0&0&0\\
241&0&-30&-30&\ast&\ast&\ast&\ast\\
245&-6&-14&14&17&-17&17&-17\\
247&0&0&0&0&0&0&0\\
251&0&0&0&\frac{4}{\sqrt{6}}&-\frac{4}{\sqrt{6}}&-\frac{4}{\sqrt{6}}&\frac{4}{\sqrt{6}}\\
253&0&0&0&0&0&0&0\\
257&-4&-32&32&0&0&0&0\\
259&0&0&0&0&0&0&0\\
263&0&0&0&0&0&0&0\\
265&0&-28&-28&\ast&\ast&\ast&\ast\\
269&-4&-26&26&-3&3&-3&3\\
271&0&0&0&\ast&\ast&\ast&\ast\\
275&0&0&0&\frac{28}{\sqrt{6}}&-\frac{28}{\sqrt{6}}&-\frac{28}{\sqrt{6}}&\frac{28}{\sqrt{6}}\\
277&0&\ast&\ast&0&0&0&0\\
281&4&32&-32&0&0&0&0\\
283&0&0&0&0&0&0&0\\
287&0&0&0&0&0&0&0\\
289&0&47&47&-17&-17&-17&-17\\
293&-2&-34&34&-9&9&-9&9\\
295&0&0&0&\ast&\ast&\ast&\ast\\
299&0&0&0&0&0&0&0\\
301&0&0&0&0&0&0&0\\
305&3&24&-24&0&0&0&0\\
307&0&0&0&0&0&0&0\\
311&0&0&0&0&0&0&0\\
313&0&\ast&\ast&\ast&\ast&\ast&\ast\\
317&4&22&-22&-5&5&-5&5\\
319&0&0&0&\ast&\ast&\ast&\ast\\
323&0&0&0&0&0&0&0\\
325&0&-4&-4&0&0&0&0\\
329&0&0&0&0&0&0&0\\
331&0&0&0&0&0&0&0\\
335&0&0&0&0&0&0&0\\
337&0&\ast&\ast&0&0&0&0\\
341&2&0&0&-8&8&-8&8\\
343&0&0&0&\ast&\ast&\ast&\ast\\
347&0&0&0&-\frac{20}{\sqrt{6}}&\frac{20}{\sqrt{6}}&\frac{20}{\sqrt{6}}&-\frac{20}{\sqrt{6}}\\
349&0&\ast&\ast&0&0&0&0\\
353&-2&-16&16&0&0&0&0\\
355&0&0&0&0&0&0&0\\
359&0&0&0&0&0&0&0\\
361&0&-19&-19&\ast&\ast&\ast&\ast\\
365&-2&12&-12&14&-14&14&-14\\
\end{array}
\end{displaymath}
\begin{displaymath}
\begin{array}{cccccccc}
\quad n \quad&\quad a(n)\quad & \quad b_1(n) \quad & \quad b_2(n)\quad & \quad c_1(n)\quad & \quad c_2(n) \quad & \quad c_3(n)\quad & \quad c_4(n)\quad\\
367&0&0&0&\ast&\ast&\ast&\ast\\
371&0&0&0&2\sqrt{6}&-2\sqrt{6}&-2\sqrt{6}&2\sqrt{6}\\
373&0&\ast&\ast&0&0&0&0\\
377&-5&-40&40&0&0&0&0\\
379&0&0&0&0&0&0&0\\
383&0&0&0&0&0&0&0\\
385&0&0&0&\ast&\ast&\ast&\ast\\
389&6&34&-34&-7&7&-7&7\\
391&0&0&0&0&0&0&0\\
395&0&0&0&6\sqrt{6}&-6\sqrt{6}&-6\sqrt{6}&6\sqrt{6}\\
397&0&\ast&\ast&0&0&0&0\\
401&-5&-40&40&0&0&0&0\\
403&0&0&0&0&0&0&0\\
407&0&0&0&0&0&0&0\\
409&0&\ast&\ast&\ast&\ast&\ast&\ast\\
413&4&0&0&-16&16&-16&16\\
415&0&0&0&\ast&\ast&\ast&\ast\\
419&0&0&0&-\frac{28}{\sqrt{6}}&\frac{28}{\sqrt{6}}&\frac{28}{\sqrt{6}}&-\frac{28}{\sqrt{6}}\\
421&0&\ast&\ast&0&0&0&0\\
425&-1&-8&8&0&0&0&0\\
427&0&0&0&0&0&0&0\\
431&0&0&0&0&0&0&0\\
433&0&\ast&\ast&\ast&\ast&\ast&\ast\\
437&0&0&0&0&0&0&0\\
439&0&0&0&\ast&\ast&\ast&\ast\\
443&0&0&0&-\frac{20}{\sqrt{6}}&\frac{20}{\sqrt{6}}&\frac{20}{\sqrt{6}}&-\frac{20}{\sqrt{6}}\\
445&0&-32&-32&0&0&0&0\\
449&5&40&-40&0&0&0&0\\
451&0&0&0&0&0&0&0\\
455&0&0&0&0&0&0&0\\
457&0&\ast&\ast&\ast&\ast&\ast&\ast\\
461&-2&-38&38&-11&11&-11&11\\
463&0&0&0&\ast&\ast&\ast&\ast\\
467&0&0&0&-\frac{28}{\sqrt{6}}&\frac{28}{\sqrt{6}}&\frac{28}{\sqrt{6}}&-\frac{28}{\sqrt{6}}\\
469&0&0&0&0&0&0&0\\
473&0&0&0&0&0&0&0\\
475&0&0&0&0&0&0&0\\
479&0&0&0&0&0&0&0\\
481&0&48&48&0&0&0&0\\
485&4&36&-36&2&-2&2&-2\\
487&0&0&0&\ast&\ast&\ast&\ast\\
491&0&0&0&-\frac{16}{\sqrt{6}}&\frac{16}{\sqrt{6}}&\frac{16}{\sqrt{6}}&-\frac{16}{\sqrt{6}}\\
493&0&-80&-80&0&0&0&0\\
497&0&0&0&0&0&0&0\\
499&0&0&0&0&0&0&0\\
503&0&0&0&0&0&0&0\\
\end{array}
\end{displaymath}
\begin{displaymath}
\begin{array}{cccccccc}
\quad n \quad&\quad a(n)\quad & \quad b_1(n) \quad & \quad b_2(n)\quad & \quad c_1(n)\quad & \quad c_2(n) \quad & \quad c_3(n)\quad & \quad c_4(n)\quad\\
505&0&-4&-4&\ast&\ast&\ast&\ast\\
509&2&-10&10&-13&13&-13&13\\
511&0&0&0&\ast&\ast&\ast&\ast\\
515&0&0&0&-6\sqrt{6}&6\sqrt{6}&6\sqrt{6}&-6\sqrt{6}\\
517&0&0&0&0&0&0&0\\
521&-5&-40&40&0&0&0&0\\
523&0&0&0&0&0&0&0\\
527&0&0&0&0&0&0&0\\
529&0&-23&-23&\ast&\ast&\ast&\ast\\
533&-4&-32&32&0&0&0&0\\
535&0&0&0&\ast&\ast&\ast&\ast\\
539&0&0&0&\frac{68}{\sqrt{6}}&-\frac{68}{\sqrt{6}}&-\frac{68}{\sqrt{6}}&\frac{68}{\sqrt{6}}\\
541&0&\ast&\ast&0&0&0&0\\
545&5&40&-40&0&0&0&0\\
547&0&0&0&0&0&0&0\\
551&0&0&0&0&0&0&0\\
553&0&0&0&\ast&\ast&\ast&\ast\\
557&-8&-38&38&13&-13&13&-13\\
559&0&0&0&0&0&0&0\\
563&0&0&0&\frac{20}{\sqrt{6}}&-\frac{20}{\sqrt{6}}&-\frac{20}{\sqrt{6}}&\frac{20}{\sqrt{6}}\\
565&0&-32&-32&0&0&0&0\\
569&-5&-40&40&0&0&0&0\\
571&0&0&0&0&0&0&0\\
575&0&0&0&0&0&0&0\\
577&0&\ast&\ast&\ast&\ast&\ast&\ast\\
581&-2&0&0&8&-8&8&-8\\
583&0&0&0&\ast&\ast&\ast&\ast\\
587&0&0&0&\frac{32}{\sqrt{6}}&-\frac{32}{\sqrt{6}}&-\frac{32}{\sqrt{6}}&\frac{32}{\sqrt{6}}\\
589&0&0&0&0&0&0&0\\
593&2&16&-16&0&0&0&0\\
595&0&0&0&0&0&0&0\\
599&0&0&0&0&0&0&0\\
601&0&\ast&\ast&\ast&\ast&\ast&\ast\\
605&-8&-22&22&21&-21&21&-21\\
607&0&0&0&\ast&\ast&\ast&\ast\\
611&0&0&0&0&0&0&0\\
613&0&\ast&\ast&0&0&0&0\\
617&4&32&-32&0&0&0&0\\
619&0&0&0&0&0&0&0\\
623&0&0&0&0&0&0&0\\
625&0&-19&-19&\ast&\ast&\ast&\ast\\
629&12&96&-96&0&0&0&0\\
631&0&0&0&\ast&\ast&\ast&\ast\\
635&0&0&0&2\sqrt{6}&-2\sqrt{6}&-2\sqrt{6}&2\sqrt{6}\\
641&1&8&-8&0&0&0&0\\
643&0&0&0&0&0&0&0\\
647&0&0&0&0&0&0&0\\
\end{array}
\end{displaymath}
\begin{displaymath}
\begin{array}{cccccccc}
\quad n \quad&\quad a(n)\quad & \quad b_1(n) \quad & \quad b_2(n)\quad & \quad c_1(n)\quad & \quad c_2(n) \quad & \quad c_3(n)\quad & \quad c_4(n)\quad\\
649&0&0&0&\ast&\ast&\ast&\ast\\
653&2&26&-26&5&-5&5&-5\\
655&0&0&0&\ast&\ast&\ast&\ast\\
659&0&0&0&-\frac{32}{\sqrt{6}}&\frac{32}{\sqrt{6}}&\frac{32}{\sqrt{6}}&-\frac{32}{\sqrt{6}}\\
661&0&\ast&\ast&0&0&0&0\\
665&0&0&0&0&0&0&0\\
667&0&0&0&0&0&0&0\\
671&0&0&0&0&0&0&0\\
673&0&\ast&\ast&\ast&\ast&\ast&\ast\\
677&4&2&-2&-15&15&-15&15\\
679&0&0&0&\ast&\ast&\ast&\ast\\
637&0&-28&-28&0&0&0&0\\
683&0&0&0&-15i&15i&15i&-15i\\
685&0&16&16&0&0&0&0\\
689&7&56&-56&0&0&0&0\\
691&0&0&0&0&0&0&0\\
695&0&0&0&0&0&0&0\\
697&0&-64&-64&0&0&0&0\\
701&-4&-10&10&11&-11&11&-11\\
703&0&0&0&0&0&0&0\\
707&0&0&0&2\sqrt{6}&-2\sqrt{6}&-2\sqrt{6}&2\sqrt{6}\\
709&0&\ast&\ast&0&0&0&0\\
713&0&0&0&0&0&0&0\\
715&0&0&0&0&0&0&0\\
719&0&0&0&0&0&0&0\\
721&0&0&0&\ast&\ast&\ast&\ast\\
725&-4&10&-10&21&-21&21&-21\\
727&0&0&0&\ast&\ast&\ast&\ast\\
731&0&0&0&0&0&0&0\\
733&0&\ast&\ast&0&0&0&0\\
737&0&0&0&0&0&0&0\\
739&0&0&0&0&0&0&0\\
743&0&0&0&0&0&0&0\\
745&0&-28&-28&\ast&\ast&\ast&\ast\\
749&4&0&0&\ast&\ast&\ast&\ast\\
751&0&0&0&\ast&\ast&\ast&\ast\\
755&0&0&0&-10\sqrt{6}&10\sqrt{6}&10\sqrt{6}&-10\sqrt{6}\\
757&0&\ast&\ast&0&0&0&0\\
761&5&40&-40&0&0&0&0\\
763&0&0&0&0&0&0&0\\
767&0&0&0&0&0&0&0\\
\end{array}
\end{displaymath}

\end{document}