# Leading twist distribution amplitudes of $P$-wave nonrelativistic mesons.

V.V. Braguta<EMAIL_ADDRESS>A.K. Likhoded<EMAIL_ADDRESS>A.V. Luchinsky
<EMAIL_ADDRESS>Institute for High Energy Physics, Protvino, Russia

###### Abstract

This paper is devoted to the study of the leading twist distribution
amplitudes of $P$-wave nonrelativistic mesons. It is shown that at the leading
order approximation in relative velocity of quark-antiquark pair inside the
mesons these distribution amplitudes can be expressed through one universal
function. As an example, the distribution amplitudes of $P$-wave charmonia
mesons are considered. Within QCD sum rules the model for the universal
function of $P$-wave charmonia mesons is built. In addition, it is found the
relations between the moments of the universal function and the
nonrelativistic QCD matrix elements that control relativistic corrections to
any amplitude involving $P$-wave charmonia. Our calculation shows that
characteristic size of these corrections is of order of $\sim 30\%$.

###### pacs:

12.38.-t, 12.38.Bx, 13.66.Bc, 13.25.Gv

##  I Introduction

Hard exclusive processes are very interesting both from theoretical and
experimental points of view. Commonly, theoretical approach to the description
of such processes is based on the factorization theorem Lepage:1980fj ;
Chernyak:1983ej . Within this theorem the amplitude of hard exclusive process
can be separated into two parts. The first part is partons production at very
small distances, which can be treated within perturbative QCD. The second part
is hardronization of the partons at larger distances. This part contains
information about nonperturbative dynamic of strong interaction. For hard
exclusive processes it can be parameterized by process independent
distribution amplitudes (DA), which can be considered as hadrons’ wave
functions at light like separation between the partons in the hadron. It
should stressed that DAs are very important for the calculation of the
amplitude of any hard exclusive process.

Recently, the leading twist DAs of $S$-wave nonrelativistic mesons have become
the object of intensive study Bodwin:2006dm ; Ma:2006hc ; Braguta:2006wr ;
Braguta:2007fh ; Braguta:2007tq ; Choi:2007ze ; Feldmann:2007id ; Bell:2008er
. Knowledge about these DAs allowed one to build some models for $S$-wave
charmonia DAs, that can be used in practical calculations. In this paper
general properties of the leading twist DAs of $P$-wave nonrelativistc mesons
will be studied. The results of this study will be used to build the model for
the DAs of $P$-wave charmonia mesons, that can be used in calculations.

This paper is organized as follows. In the next section all definitions of the
DAs of $P$-wave mesons will be given. These DAs will be studied in section III
at the leading order approximation in relative velocity of quark-antiquark
pair inside $P$-wave meson. In section IV QCD sum rules will be applied to the
calculation of the moments of the $P$-wave charmonia DAs. Using the results of
this study the model of the $P$-wave charmonia DAs will be build in section V.
In the last section the results of this paper will be summarized.

##  II Definitions of the distribution amplitudes.

In this section the definitions of the leading twist distribution amplitudes
(DA) of $P$-wave nonrelativistic mesons will be given. In the conventional
quark model nonrelativistic mesons are quark-antiquark($Q\bar{Q}$) bound
states. In these mesons the quark-antiquark pair can be in the spin singlet or
spin triplet states. Since orbital momentum of the quark-antiquark pair is
unity one can conclude that there are four $P$-wave mesons:
$\chi_{0}(^{3}P_{0}),\chi_{1}(^{3}P_{1}),\chi_{2}(^{3}P_{2}),h(^{1}P_{1})$.
The leading twist DAs of these mesons can be defined as follows.
for the $\chi_{0}$ meson:

$\displaystyle{\langle\chi_{0}(p)|{\bar{Q}}(z)\gamma_{\mu}[z,-z]Q(-z)|0\rangle}$
$\displaystyle=$ $\displaystyle
f_{\chi_{0}}p_{\mu}\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\phi_{\chi_{0}}(\xi,\mu),$
(1)

for the $\chi_{1}$ meson:

$\displaystyle{\langle\chi_{1}(p,\epsilon_{\lambda=0})|{\bar{Q}}(z)\gamma_{\mu}\gamma_{5}[z,-z]Q(-z)|0\rangle}$
$\displaystyle=$ $\displaystyle
f^{\prime}_{\chi_{1}}p_{\mu}\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\phi^{\prime}_{\chi_{1}}(\xi,\mu),$
$\displaystyle{\langle\chi_{1}(p,\epsilon_{\lambda=\pm
1})|{\bar{Q}}(z)\sigma_{\mu\nu}[z,-z]Q(-z)|0\rangle}$ $\displaystyle=$
$\displaystyle
f_{\chi_{1}}e_{\mu\nu\alpha\beta}\epsilon^{\alpha}p^{\beta}\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\phi_{\chi_{1}}(\xi,\mu),$
(2)

for the $h$ meson:

$\displaystyle{\langle
h(p,\epsilon_{\lambda=0})|{\bar{Q}}(z)\gamma_{\mu}\gamma_{5}[z,-z]Q(-z)|0\rangle}$
$\displaystyle=$ $\displaystyle
f_{h}p_{\mu}\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\phi_{h}(\xi,\mu),$
$\displaystyle{\langle h(p,\epsilon_{\lambda=\pm
1})|{\bar{Q}}(z)\sigma_{\mu\nu}[z,-z]Q(-z)|0\rangle}$ $\displaystyle=$
$\displaystyle
f^{\prime}_{h}e_{\mu\nu\alpha\beta}\epsilon^{\alpha}p^{\beta}\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\phi^{\prime}_{h}(\xi,\mu),$
(3)

for the $\chi_{2}$ meson:

$\displaystyle{\langle\chi_{2}(p,\epsilon_{\lambda=0})|{\bar{Q}}(z)\gamma_{\mu}[z,-z]Q(-z)|0\rangle}$
$\displaystyle=$ $\displaystyle
f_{\chi_{2}}p_{\mu}\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\phi_{\chi_{2}}(\xi,\mu),$
$\displaystyle{\langle\chi_{2}(p,\epsilon_{\lambda=\pm
1})|{\bar{Q}}(z)\sigma_{\mu\nu}[z,-z]Q(-z)|0\rangle}$ $\displaystyle=$
$\displaystyle\tilde{f}_{\chi_{2}}M_{\chi_{2}}(\rho_{\mu}p_{\nu}-\rho_{\nu}p_{\mu})\int^{1}_{-1}d\xi\,e^{i(pz)\xi}\tilde{\phi}_{\chi_{2}}(\xi,\mu),~{}~{}~{}\rho_{\mu}=\frac{\epsilon_{\mu\sigma}z^{\sigma}}{pz},$
(4)

where the following designations are used: $x_{1},x_{2}$ are the fractions of
momentum of meson carried by quark and antiquark correspondingly,
$\xi=x_{1}-x_{2}$, $p,~{}\epsilon$ are the momentum and polarizations of
$P$-wave mesons. For the mesons $\chi_{1},h$ the polarization $\epsilon$ is
described by the four vector $\epsilon_{\mu}$, for the $\chi_{2}$ meson the
polarization $\epsilon$ is described by the tensor $\epsilon_{\mu\nu}$. The
factor $[z,-z]$, that makes matrix elements (1)-(4) gauge invariant, is
defined as

$\displaystyle[z,-z]=P\exp[ig\int_{-z}^{z}dx^{\mu}A_{\mu}(x)].$ (5)

In applications it is useful to rewrite the four-vector $\rho_{\mu}$ in the
following way. Evidently, one can write the polarization of the $\chi_{2}$
meson in terms of the polarization of two vector mesons. Thus for the
transverse polarization of the meson $\chi_{2}$ one has
$\epsilon^{\mu\nu}_{\lambda=\pm 1}=(\epsilon_{\lambda=\pm
1}^{\mu}\cdot\epsilon_{\lambda=0}^{\nu}+\epsilon_{\lambda=0}^{\nu}\cdot\epsilon_{\lambda=\pm
1}^{\mu})/\sqrt{2}$ ($\epsilon^{+}_{\mu\nu}\epsilon^{\mu\nu}=1$). If we
further contract the polarization tensor $\epsilon_{\mu\nu}$ with lightlike
four-vector $z$, to the leading twist accuracy we will get
$\epsilon^{\mu\nu}z_{\nu}=\epsilon_{\lambda=\pm
1}^{\mu}(pz)/(\sqrt{2}M_{\chi_{2}})$ or $\rho^{\mu}=\epsilon_{\lambda=\pm
1}^{\mu}/(\sqrt{2}M_{\chi_{2}})$. This form of the vector $\rho$ can be used
in the calculation with the leading twist accuracy. It should be noted that
the states of the $\chi_{2}$ meson with the polarizations $\lambda=\pm 2$ give
contribution only to higher twist DAs. Since, this paper is devoted to the
study of the leading twist DAs we don’t consider these states.

The functions without primes $\phi_{\chi 0}(\xi),\phi_{\chi
1}(\xi),\phi_{h}(\xi),\phi_{\chi 2}(\xi),\tilde{\phi}_{\chi 2}(\xi)$ are $\xi$
odd and they are normalized as

$\displaystyle\int_{-1}^{1}d\xi~{}\xi\phi(\xi)=1.$ (6)

The functions with primes $\phi^{\prime}_{\chi 1}(\xi),\phi^{\prime}_{h}(\xi)$
are $\xi$ even and they are normalized as

$\displaystyle\int_{-1}^{1}d\xi~{}\phi^{\prime}(\xi)=1.$ (7)

In this paper all DAs will be parameterized by their moments

$\displaystyle\langle\xi^{n}\rangle=\int_{-1}^{1}d\xi~{}\xi^{n}\phi(\xi).$ (8)

Evidently, for the DAs with primes all odd moments are zero. For the DAs
without primes all even moments are zero. To separate the moments of different
DAs of the $\chi_{2}$ meson, below we are going to use the following
designations: $\langle\xi^{n}\rangle$ for the moments of the function
$\phi_{\chi_{2}}(\xi)$, $\langle\tilde{\xi}^{n}\rangle$ for the moments of the
function $\tilde{\phi}_{\chi_{2}}(\xi)$.

The distribution amplitudes $\phi(\xi)$ and the constants $f$ that
parameterizes corresponding currents (1)-(4) are scale dependent objects. For
applications it is useful to write how they depend on the scale
Chernyak:1983ej . To do this we expand DA in the series

$\displaystyle\phi(\xi,\mu)=\frac{3}{4}(1-\xi^{2})\biggl{[}\sum_{n=0}^{\infty}a_{n}(\mu)C_{n}^{3/2}(\xi)\biggr{]},$
(9)

where $C_{n}^{3/2}(\xi)$ are Gegenbauer polynomials. At the leading
logarithmic accuracy the coefficients $a_{n}$ are renormalized
multiplicatively

$\displaystyle
a_{n}(\mu)=\biggl{(}\frac{\alpha_{s}(\mu)}{\alpha_{s}(\mu_{0})}\biggr{)}^{\gamma_{n}/{b_{0}}}a_{n}^{L,T}(\mu_{0}),$
(10)

where $\gamma_{n}$ are the anomalous dimensions. For the current
$\bar{Q}\sigma_{\mu\nu}[z,-z]Q$ the anomalous dimentions are

$\displaystyle\gamma_{n}$ $\displaystyle=$ $\displaystyle
C_{f}\biggl{(}1+4\sum_{j=2}^{n+1}\frac{1}{j}\biggr{)},~{}~{}b_{0}=11-\frac{2}{3}n_{\rm
fl},~{}~{}C_{f}=\frac{4}{3},$ (11)

for the other currents

$\displaystyle\gamma_{n}$ $\displaystyle=$ $\displaystyle
C_{f}\biggl{(}1-\frac{2}{(n+1)(n+2)}+4\sum_{j=2}^{n+1}\frac{1}{j}\biggr{)}.$

It is clear that the DAs without primes contain only n-odd terms in series
(9). DAs with primes contain n-even terms in the series.

The constants defined in equations (1)-(4) are multiplicatively
renormalizable. Using formulas (9)-(LABEL:an_dimVA) one can determine the
evolution of these constants

$\displaystyle
f(\mu)=\biggl{(}\frac{\alpha_{s}(\mu)}{\alpha_{s}(\mu_{0})}\biggr{)}^{\gamma/{b_{0}}}f(\mu_{0}).$
(13)

For the constants $f_{\chi_{0}},f_{h},f_{\chi_{2}}$ the anomalous dimentions
$\gamma$ is equal to $8/3C_{f}$, for the constants
$f_{\chi_{1}},\tilde{f}_{\chi_{2}}$ $\gamma=3C_{f}$, for the constant
$f^{\prime}_{h}$ $\gamma=C_{f}$, for the constant $f^{\prime}_{\chi_{1}}$
$\gamma=0$. It should be noted that if the anomalous dimensions $\gamma$ of
the constants $f$ are factored from sum (9), the anomalous dimensions of the
remaining terms equal to the difference $\gamma_{n}-\gamma$.

##  III The distribution amplitudes at the leading order approximation in
relative velocity.

In this section the DAs under study will be considered at the leading order
approximation in relative velocity of quark-antiquark pairs inside the mesons.
First let us consider the DA of the $\chi_{0}$ meson. The moments of this DA
can be represented as follows

$\displaystyle
f_{\chi_{0}}(pz)^{n+1}<\xi^{n}>_{\chi_{0}}=\langle\chi_{0}|\bar{Q}\hat{z}(-iz{\overset{\leftrightarrow}{D}})^{n}Q|0\rangle.$
(14)

To get the expressions for the constant $f_{\chi_{0}}$ and the moment
$<\xi^{n}>_{\chi_{0}}$ one needs to calculate the matrix element in the right
hand side. At the leading order approximation this calculation can be done
using projector Bodwin:2002hg ; Braaten:2002fi

$\displaystyle\bar{Q}(\bar{p})Q(p)\to\int
dq\frac{\varphi(-q^{2})}{\sqrt{3m_{Q}}}\frac{1}{4\sqrt{2}E(E+m_{Q})}(\hat{\bar{p}}-m_{Q})\Gamma(\hat{P}+2E)(\hat{p}+m_{Q}),$
(15)

where $P,q$ are the total and relative momentum of the $Q\bar{Q}$ pair,
$m_{Q}$ is the mass of the quark $Q$, $p=P/2+q,\bar{p}=P/2-q$,
$E^{2}=P^{2}/4=m_{Q}^{2}-q^{2}$. The matrix $\Gamma=\gamma_{5},\hat{e}_{S}$
for the spin singlet and spin triplet quark-antiquark pair correspondingly,
where $e_{S}$ is the spin polarization of this pair. The scalar products
$P\cdot e_{S}=0,P\cdot q=0$. In the center mass frame the $dq$ is reduced to
the $d^{3}{\bf q}/(2\pi)^{3}$ and the function $\varphi(-q^{2})$ is reduced to
the usual nonrelativistic wave function $\phi(\bf{q^{2}})$. For the $\chi_{0}$
meson the wave function $\varphi(\bf{q^{2}})$ can be written in the form

$\displaystyle\varphi({\bf{q^{2}}})=\frac{{\bf e_{S}\cdot
q}}{\sqrt{3}}\psi({\bf q}).$ (16)

At the leading order approximation in relative velocity the function $\psi(\bf
q)$ is universal function for all $P$-wave mesons. It is normalized as

$\displaystyle\int\frac{d^{3}{\bf q}}{(2\pi)^{3}}{\bf q^{i}}{\bf
q^{j}}|\psi({\bf q})|^{2}=\delta^{ij}.$ (17)

With this normalization of the function $\psi({\bf q})$, the function
$\varphi({\bf q})$ is normalized as

$\displaystyle\int\frac{d^{3}{\bf q}}{(2\pi)^{3}}|\varphi(\bf{q^{2}})|^{2}=1.$
(18)

The same normalization condition will be used for the wave functions of all
mesons under consideration. Using equations (14), (15) and (16) one gets the
result

$\displaystyle
f_{\chi_{0}}<\xi^{n+1}>_{\chi_{0}}=-2^{n+1}\sqrt{\frac{2}{m}_{Q}}\frac{A_{n+2}}{M^{n+1}_{\chi_{0}}}\frac{1}{n+3},$
(19)

where $M_{\chi_{0}}$ is the mass of the $\chi_{0}$ meson, $A_{n}$ equals to

$\displaystyle A_{n}=\int\frac{d^{3}{\bf q}}{(2\pi)^{3}}{|\bf q|}^{n}\psi({\bf
q}).$ (20)

It should be noted that because of the Coulombic part of the nonrelativistic
QCD potential the right hand side of equation (19) is ultraviolet divergent
Bodwin:1994jh ; Bodwin:2006dn . The moments of the DA in the left hand side
are QCD operators (see equation (14) ). In full QCD the Coulombic part of the
nonrelativistic potential corresponds to the rescattering of the quark-
antiquark pair of the QCD operators which is also ultraviolet divergent
Chernyak:1983ej . To control the divergences in the right and left hand sides
of equation (19) it is assumed that both sides are regularized within
dimensional regularization.

It is interesting to note that relation (19) is closely connected with
Brodsky-Huang-Lepage (BHL) Brodsky:1981jv  procedure. This fact can be seen as
follows. Let us rewrite this relation as follows:

$\displaystyle\int
d\xi\xi^{n+1}\phi_{\chi_{0}}(\xi)\sim\frac{1}{n+3}\frac{A_{n+2}}{m_{Q}^{n+1}}\sim\frac{1}{n+3}\int{\bf
q^{2}}d{\bf q}\frac{{\bf q}^{n+2}}{m_{Q}^{n+1}}\psi({\bf q})\sim\int d^{3}{\bf
q}\biggl{(}\frac{q_{z}}{m_{Q}}\biggr{)}^{n+1}q_{z}\psi({\bf q}).$ (21)

Note that $q_{z}\psi({\bf q})$ is the $L_{z}=0$ component of the wave function
$\varphi({\bf q^{2}})$ which is the only component important for the leading
twist DA. So, the last relation can be written as follows

$\displaystyle\int d\xi\xi^{n+1}\phi_{\chi_{0}}(\xi)\sim\int
dq_{z}\biggl{(}\frac{q_{z}}{m_{Q}}\biggr{)}^{n+1}\int d^{2}{\bf
q}_{\perp}\varphi_{L_{z}=0}({\bf q^{2}}).$ (22)

Further, we change the variables in the right side of this equation

$\displaystyle{\bf q}_{\perp}\to{\bf q}_{\perp},\quad q_{z}\to\xi{M_{0}},\quad
M_{0}^{2}=\frac{M_{Q}^{2}+{\bf q}_{\perp}^{2}}{1-\xi^{2}}.$ (23)

Note also that $\xi\ll 1,{\bf q}_{\perp}\ll M_{Q}$ and at the leading order
approximation in relative velocity of the quark-antiquark pair in the meson
relation (22) can be written as follows

$\displaystyle\int d\xi\xi^{n+1}\times\phi_{\chi_{0}}(\xi)\sim\int
d\xi\xi^{n+1}\times\int d^{2}{\bf
q}_{\perp}~{}\varphi_{L_{z}=0}\biggl{(}\frac{M_{Q}^{2}\xi^{2}+{\bf
q}_{\perp}^{2}}{1-\xi^{2}}\biggr{)}.$ (24)

So, at this level of accuracy the DA is just

$\displaystyle\phi_{\chi_{0}}(\xi)\sim\int d^{2}{\bf
q}_{\perp}~{}\varphi_{L_{z}=0}\biggl{(}\frac{M_{Q}^{2}\xi^{2}+{\bf
q}_{\perp}^{2}}{1-\xi^{2}}\biggr{)},$ (25)

what coincides with BHL procedure.

It is not difficult to get the relations for the other DAs and mesons

$\displaystyle{\chi_{1}~{}\mbox{meson:}}~{}~{}~{}~{}\phi({\bf{q^{2}}})$
$\displaystyle=$ $\displaystyle\frac{\epsilon_{ijk}{{\bf e}^{i}{\bf
e}_{S}^{j}{\bf q}^{k}}}{\sqrt{2}}\psi({\bf q}),$
$\displaystyle{h~{}\mbox{meson:}}~{}~{}~{}~{}\phi({\bf{q^{2}}})$
$\displaystyle=$ $\displaystyle\frac{{\bf e\cdot q}}{\sqrt{2}}\psi({\bf q}),$
$\displaystyle{\chi_{2}~{}\mbox{meson:}}~{}~{}~{}~{}\phi({\bf{q^{2}}})$
$\displaystyle=$ $\displaystyle{\bf e}_{ij}\bf{e_{S}}^{i}{\bf q}^{j}\psi({\bf
q}),$ (26)

where $\bf e$ is the polarization three vector of the $\chi_{1}$ and $h$
mesons, ${\bf e}_{ij}$ is the polarization tensor of the $\chi_{2}$ meson.
Using these expressions one can get

$\displaystyle
f_{\chi_{1}}<\xi^{n+1}>_{\chi_{1}}=-2^{n+1}\sqrt{\frac{3}{m_{Q}}}\frac{A_{n+2}}{M^{n+1}_{\chi_{1}}}\frac{1}{n+3},$
$\displaystyle
f^{\prime}_{\chi_{1}}<\xi^{n}>^{\prime}_{\chi_{1}}=2^{n+1}i\sqrt{\frac{12}{m_{Q}}}\frac{A_{n+2}}{M^{n+1}_{\chi_{1}}}\frac{1}{(n+1)(n+3)},$
$\displaystyle
f_{h}<\xi^{n+1}>_{h}=-2^{n+1}\sqrt{\frac{6}{m_{Q}}}\frac{A_{n+2}}{M^{n+1}_{h}}\frac{1}{n+3},$
$\displaystyle
f^{\prime}_{h}<\xi^{n}>^{\prime}_{h}=2^{n+1}i\sqrt{\frac{6}{m_{Q}}}\frac{A_{n+2}}{M^{n+1}_{h}}\frac{1}{(n+1)(n+3)},$
$\displaystyle
f_{\chi_{2}}<\xi^{n+1}>_{\chi_{2}}=-2^{n+1}\sqrt{\frac{4}{m_{Q}}}\frac{A_{n+2}}{M^{n+1}_{\chi_{2}}}\frac{1}{n+3},$
$\displaystyle\tilde{f}_{\chi_{2}}<\tilde{\xi}^{n+1}>_{\chi_{2}}=-2^{n+1}\sqrt{\frac{6}{m_{Q}}}\frac{A_{n+2}}{M^{n+1}_{\chi_{2}}}\frac{1}{n+3}.$
(27)

The moments with primes correspond to the DAs with primes. As it was noted
above, since all the DAs without primes are $\xi$-odd only odd moments of
these DAs are nonzero. Similarly, for the DAs with primes only even moments
different from zero. For these reason we took $(n+1)$th moments of DAs without
primes and $n$th moments with primes, where n is assumed to be even. It should
be noted that since the function $\psi({\bf q})$ is universal for the $P$-wave
mesons, the $A_{n}$ is the same for all mesons under study. From relation (27)
one can conclude that all the constants $f_{i}$ can be expressed through the
only constant which will be designated as $F$

$\displaystyle
F={\sqrt{3}}{f_{\chi_{0}}}={\sqrt{2}}{f_{\chi_{1}}}=i\frac{f^{\prime}_{\chi_{1}}}{\sqrt{2}}=f_{h}=if^{\prime}_{h}=\sqrt{\frac{3}{2}}f_{\chi_{2}}={\tilde{f}_{\chi_{2}}}.$
(28)

From equations (27) one can also find that at the leading order approximation
the $n$th moments of all functions without primes coincide. This means that
these DAs are equal to one universal DA which will be designated as
$\Phi(\xi)$

$\displaystyle\Phi(\xi)=\phi_{\chi_{0}}(\xi)=\phi_{\chi_{1}}(\xi)=\phi_{h}(\xi)=\phi_{\chi_{2}}(\xi)=\tilde{\phi}_{\chi_{2}}(\xi).$
(29)

The same is true for all DAs with primes, which are equal to one function
which will be designated below as $\Psi(\xi)$

$\displaystyle\Psi(\xi)=\phi^{\prime}_{\chi_{1}}(\xi)=\phi^{\prime}_{h}(\xi).$
(30)

In addition, one can relate the moments of the $\Phi(\xi)$ to the moments of
the $\Psi(\xi)$ as follows

$\displaystyle<\xi^{n}>_{\Psi}=\frac{<\xi^{n+1}>_{\Phi}}{n+1}.$ (31)

It should be noted that the constants and DAs in relations (28)-(31) depend on
scale in a different way. This means that relations (28)-(31), which are valid
at not to large scale, will be violated at sufficiently large scale.

Recursive relation (31) determines the function $\Psi(\xi)$ through the
function $\Phi(\xi)$. One can guess the solution of this relation:

$\displaystyle\Psi(\xi)=-\int_{-1}^{\xi}dt~{}\Phi(t).$ (32)

To prove that (32) is the solution of relations (31) one should put this
function to the definition of the n-th moment (8) and integrate the resulting
expression by parts. It is seen from equations (29), (30) and (32) that all
DAs of the $P$-wave mesons are defined through the universal function
$\Psi(\xi)$. It should be noted that this fact results from the
nonrelativistic spin-symmetry, which holds at leading order in the heavy-quark
velocity. Below equations (29) and (32) will be used to build the models for
the function $\Phi(\xi)$ and $\Psi(\xi)$ of $P$-wave charmonia.

At the end of this section it is interesting to discuss the question about
relativistic corrections to the matrix elements involving $P$-wave quarkonia.
If we ignore the contribution coming from the higher Fock states, relativistic
corrections to the matrix involving, for instance, the $\chi_{0}$ meson is
given by the matrix element

$\displaystyle\langle
v^{n}\rangle_{P}=\frac{1}{(m_{Q}^{*})^{2}}\frac{\langle\chi_{0}|\chi^{+}\bigl{(}{\bf\sigma
D}\bigr{)}{\bf D^{2}}\psi|0\rangle}{\langle\chi_{0}|\chi^{+}\bigl{(}{\bf\sigma
D}\bigr{)}\psi|0\rangle},$ (33)

where $\psi$ and $\chi^{+}$ are Pauli spinor fields that annihilate a quark
and an antiquark respectively, ${\bf\sigma}$ are Pauli matrixes, $m_{Q}^{*}$
is the quark pole mass. Using equation (19) it is not difficult to obtain
general formula that connects the moment $\langle\xi^{n+1}\rangle$ with the
matrix element $\langle v^{n}\rangle_{P}$ at the leading order approximation
in relative velocity

$\displaystyle\langle
v^{n}\rangle_{P}=\frac{A_{n+2}}{A_{2}}+O(v^{n+2})=\frac{n+3}{3}\langle\xi^{n+1}\rangle+O(v^{n+2}).$
(34)

Although the derivation was done for the $\chi_{0}$ meson, at the leading
order approximation in relative velocity the matrix element $\langle
v^{n}\rangle_{P}$ is universal for all $P$-wave mesons.

##  IV Charmonia distribution amplitudes within QCD sum rules.

In this section QCD sum rules approach Shifman:1978bx ; Shifman:1978by  will
be applied to the calculation of the moments of the $P$-wave charmonia DAs.
The problem which can be met if one tries to apply QCD sum rules to the
calculation of these moments is that it is not possible to write QCD sum rules
for one DA. Commonly, the contribution of different mesons and different DAs
mix in QCD sum rules, what does not allow us to calculate the constants and
the moments of DAs (1)-(4) separately. Only some combinations of different
constants and moments can be extracted from QCD sum rules. For instance, two
point QCD sum rules for the currents
$\bar{Q}(x)\gamma_{\mu}(z{\overset{\leftrightarrow}{D}})Q(x)\cdot\bar{Q}(0)\gamma_{\mu}(z{\overset{\leftrightarrow}{D}})^{n+1}Q(0)$
can be written as follows

$\displaystyle\frac{f_{\chi_{c0}}^{2}\langle\xi^{n+1}\rangle_{\chi_{c0}}}{(m_{\chi_{c0}}^{2}+Q^{2})^{m+1}}+\frac{f_{\chi_{c2}}^{2}\langle\xi^{n+1}\rangle_{\chi_{c2}}}{(m_{\chi_{c2}}^{2}+Q^{2})^{m+1}}=\frac{1}{\pi}\int_{4m_{c}^{2}}^{s_{0}}ds~{}\frac{\mbox{Im}~{}\Pi_{\rm
pert}(s,n)}{(s+Q^{2})^{m+1}}+\Pi^{(m)}_{\rm npert}(Q^{2},n).$ (35)

The expressions for the $\mbox{Im}~{}\Pi_{\rm pert}(s,n)$, $\Pi^{(m)}_{\rm
npert}(Q^{2},n)$ will be given below (equations (37),(38)). It is seen from
this example that it is not possible to extract the constants
$f_{\chi_{c0}},f_{\chi_{c2}}$ or the moments
$\langle\xi^{n+1}\rangle_{\chi_{c0}},\langle\xi^{n+1}\rangle_{\chi_{c2}}$ from
(35) separately. Evidently, this strongly restricts the accuracy of the
calculation. The only QCD sum rules which are free from this problem are two
point sum rules for the currents
$\bar{Q}(x)\gamma_{\mu}\gamma_{5}(z{\overset{\leftrightarrow}{D}})Q(x)\cdot\bar{Q}(0)\gamma_{\mu}\gamma_{5}(z{\overset{\leftrightarrow}{D}})^{n+1}Q(0)$

$\displaystyle\frac{f_{h_{c}}^{2}\langle\xi^{n+1}\rangle_{h_{c}}}{(m_{h_{c}}^{2}+Q^{2})^{m+1}}=\frac{1}{\pi}\int_{4m_{c}^{2}}^{s_{0}}ds~{}\frac{\mbox{Im}~{}\Pi_{\rm
pert}(s,n)}{(s+Q^{2})^{m+1}}+\Pi^{(m)}_{\rm npert}(Q^{2},n),$ (36)

where $\mbox{Im}~{}\Pi_{\rm pert}(s,n)$, $\Pi^{(m)}_{\rm npert}(Q^{2},n)$ can
be written as

$\displaystyle\mbox{Im}~{}\Pi_{\rm
pert}(s,n)=\frac{3}{8\pi}v^{n+3}(\frac{1}{n+3}-\frac{v^{2}}{n+5}),~{}~{}~{}~{}v^{2}=1-\frac{4m_{c}^{2}}{s},$
(37) $\displaystyle\Pi^{(m)}_{\rm npert}(Q^{2},n)$ $\displaystyle=$
$\displaystyle\Pi_{1}^{(m)}(Q^{2},n)+\Pi_{2}^{(m)}(Q^{2},n)+\Pi_{3}^{(m)}(Q^{2},n),$
(38) $\displaystyle\Pi_{1}^{(m)}(Q^{2},n)$ $\displaystyle=$
$\displaystyle\frac{\langle\alpha_{s}G^{2}\rangle}{24\pi}(m+1)\int_{-1}^{1}d\xi~{}\biggl{(}\xi^{n+2}+\frac{n(n+1)}{4}\xi^{n}(1-\xi^{2})\biggr{)}\frac{(1-\xi^{2})^{m+2}}{\bigl{(}4m_{c}^{2}+Q^{2}(1-\xi^{2})\bigr{)}^{m+2}},$
$\displaystyle\Pi_{2}^{(m)}(Q^{2},n)$ $\displaystyle=$
$\displaystyle-\frac{\langle\alpha_{s}G^{2}\rangle}{6\pi}m_{c}^{2}(m^{2}+3m+2)\int_{-1}^{1}d\xi~{}\xi^{n+2}\bigl{(}1+3\xi^{2}\bigr{)}\frac{(1-\xi^{2})^{m+1}}{\bigl{(}4m_{c}^{2}+Q^{2}(1-\xi^{2})\bigr{)}^{m+3}},$
$\displaystyle\Pi_{3}^{(m)}(Q^{2},n)$ $\displaystyle=$
$\displaystyle\frac{\langle\alpha_{s}G^{2}\rangle}{384\pi}n(n+1)(m+1)\int_{-1}^{1}d\xi~{}\xi^{n}\frac{(1-\xi^{2})^{m+3}}{\bigl{(}4m_{c}^{2}+Q^{2}(1-\xi^{2})\bigr{)}^{m+2}}.$

In the calculation we take $Q^{2}=4m_{c}^{2}$ Reinders:1984sr . In the
numerical analysis of QCD sum rules the values of parameters $m_{c}$ and
$\langle\alpha_{s}G^{2}/\pi\rangle$ will be taken from paper Reinders:1984sr :

$\displaystyle m_{c}=1.24\pm
0.02~{}\mbox{GeV},~{}~{}\langle\frac{\alpha_{s}}{\pi}G^{2}\rangle=0.012\pm
30\%~{}\mbox{GeV}^{4}.$ (39)

First sum rules (36) will be applied to the calculation of the constant
$f_{h_{c}}^{2}$. It is not difficult to express the constant $f_{h_{c}}^{2}$
from equation (36) at $n=0$ as a function of $m$. For too small values of $m$
($m<m_{1}$) there is large contributions from higher resonances and continuum
which spoil sum rules (36). Although for $m\gg m_{1}$ these contributions are
strongly suppressed, it is not possible to apply sum rules for too large $m$
($m>m_{2}$) since the contribution arising from higher dimensional vacuum
condensates rapidly grows with $m$ what invalidates our approximation. If
$m_{1}<m_{2}$ there is some region of applicability of sum rules (36)
$[m_{1},m_{2}]$ where the resonance and the higher dimensional vacuum
condensates contributions are not too large. Within this region
$f_{h_{c}}^{2}$ as a function of $m$ varies slowly and one can determine the
value of this constant. The value of the continuum threshold $s_{0}$ must be
taken so that to appear stability region $[m_{1},m_{2}]$. Our calculation
shows that for the central values of parameters (39) there exists stability
region for $s_{0}>(4.3$ GeV$)^{2}$. If the value of the continuum threshold
$s_{0}$ is varied in the region $s_{0}\in(4.3^{2},\infty)$ GeV2, the value of
the constant $f_{h_{c}}^{2}$ can be written as $f_{h_{c}}^{2}=(0.037\pm
0.005)$ GeV2. In addition to the uncertainty due to the variation of the value
of $s_{0}$, there are uncertainties due to the variation of the values of the
$m_{c}$ (which is $\pm 0.004$) and $\langle{\alpha_{s}}/{\pi}G^{2}\rangle$
(which is $\pm 0.001$). The last source of uncertainty is the radiative
corrections to the perturbative density $\mbox{Im}\Pi_{\rm pert}(s,n)$, which
will be estimated as $\alpha_{s}(m_{c})/\pi\sim 13\%$. Adding these
uncertainties in quadrature, one gets

$\displaystyle f_{h_{c}}^{2}(\mu\sim m_{c})=(0.037\pm 0.007)\mbox{GeV}^{2}.$
(40)

As it was noted above, the value of the constant $f_{h_{c}}^{2}$ is scale
dependent quantity. The characteristic scale of QCD sum rules is $\sim m_{c}$.
This means that the value of the constant $f_{h_{c}}^{2}$ is determined at the
scale $\sim m_{c}$, as it is shown in (40).

Next let us consider the moments
$\langle\xi^{3}\rangle_{h_{c}},\langle\xi^{5}\rangle_{h_{c}},\langle\xi^{7}\rangle_{h_{c}}$.
However, instead of considering QCD sum rules for $n=2,4,6$ we will consider
the ratios of sum rules at $n=2,4,6$ and the sum rules at $n=0$. Such approach
improves the accuracy of the calculation (see paper Braguta:2006wr  for
details). The analisys similar to that for the constant $f_{h_{c}}^{2}$ gives

$\displaystyle\langle\xi^{3}\rangle$ $\displaystyle=$ $\displaystyle 0.18\pm
0.03,$ $\displaystyle\langle\xi^{5}\rangle$ $\displaystyle=$ $\displaystyle
0.050\pm 0.010,$ $\displaystyle\langle\xi^{7}\rangle$ $\displaystyle=$
$\displaystyle 0.017\pm 0.004.$ (41)

It should be noted that these moments are defined at the scale $\sim m_{c}$.
It should be also noted that the values of the moments (41) are in good
agreement with potential model estimation (see paper Braguta:2006wr ). For
instance, within Buchmuller-Tye potential model Buchmuller:1980su
$\langle\xi^{3}\rangle=0.18,\langle\xi^{5}\rangle=0.047,\langle\xi^{7}\rangle=0.016$;
within Cornell potential model Eichten:1978tg
$\langle\xi^{3}\rangle=0.16,\langle\xi^{5}\rangle=0.040,\langle\xi^{7}\rangle=0.013$.

Using values (41) one can find the matrix elements that control relativisitic
corrections to any process with $P$-wave charmonia in the initial or final
state. The relationships between these matrix elements and the moments are
given in equation (34). These relations are valid up to the higher order
relativistic corrections, which can be estimated as $\langle v^{2}\rangle$.
Taking into account this additional source of uncertainty one gets

$\displaystyle\langle v^{2}\rangle_{P}$ $\displaystyle=$
$\displaystyle\frac{5}{3}{\langle\xi^{3}\rangle_{P}}=0.30\pm 0.10,$
$\displaystyle\langle v^{4}\rangle_{P}$ $\displaystyle=$
$\displaystyle\frac{7}{3}{\langle\xi^{5}\rangle_{P}}=0.12\pm 0.04,$
$\displaystyle\langle v^{6}\rangle_{P}$ $\displaystyle=$
$\displaystyle\frac{9}{3}{\langle\xi^{7}\rangle_{P}}=0.051\pm 0.018.$ (42)

In the next sections results (41) will be used to build the model for the DAs
$\Phi(\xi)$ and $\Psi(\xi)$.

##  V Model for charmonia distribution amplitudes.

To build the model of the function $\Phi(\xi,\mu\sim m_{c})$ we use Borel
version Shifman:1978bx ; Shifman:1978by  of sum rules (36) but without
continuum contribution and power corrections

$\displaystyle f_{h_{c}}^{2}\langle\xi^{n+1}\rangle
e^{-{m_{h_{c}}^{2}}/{M^{2}}}=\frac{M^{2}}{4\pi^{2}}\int_{-1}^{1}d\xi~{}\xi^{n+2}~{}\frac{3}{4}(1-\xi^{2})\mbox{exp}\biggl{(}-\frac{4m_{c}^{2}}{M^{2}}\frac{1}{1-\xi^{2}}\biggr{)}.$
(43)

Evidently, within this approximation the function $\Phi(\xi,\mu\sim m_{c})$
can be written in the form

$\displaystyle\Phi(\xi,\mu\sim
m_{c})=c(\beta_{P})(1-\xi^{2})~{}\xi~{}\mbox{exp}\biggl{(}-\frac{\beta_{P}}{1-\xi^{2}}\biggr{)},$
(44)

where $c(\beta_{P})$ is a normalization constant and $\beta_{P}$ is some
constant. We propose function (44) as the model for DAs $\Phi(\xi,\mu\sim
m_{c})$. To fix the constant $\beta_{P}$ the value of the moment
$\langle\xi^{2}\rangle$ (41) will be used. Thus we get
$\beta_{P}=3.4^{+1.5}_{-0.9}$. The constant $c(\beta)$ can be determined from
normalization condition (6). The moments of the function (44) are

$\displaystyle\langle\xi^{3}\rangle$ $\displaystyle=$ $\displaystyle~{}0.18\pm
0.03,$ (45) $\displaystyle\langle\xi^{5}\rangle$ $\displaystyle=$
$\displaystyle~{}0.047\pm 0.014,$ $\displaystyle\langle\xi^{7}\rangle$
$\displaystyle=$ $\displaystyle 0.015\pm 0.006.$

Using the model for the function $\Phi(\xi,\mu\sim m_{c})$ and equation (32)
one can get the model of the DA $\Psi(\xi,\mu\sim m_{c})$

$\displaystyle\Psi(\xi,\mu\sim m_{c})=-\int_{-1}^{\xi}dt\Phi(t,\mu\sim
m_{c})=\frac{c(\beta_{P})}{2}(1-\xi^{2})^{2}~{}E_{3}\biggl{(}\frac{\beta_{P}}{1-\xi^{2}}\biggr{)},$
(46)

where the function $E_{3}(z)$ is the exponential integral function

$\displaystyle E_{3}(z)=\int_{1}^{\infty}dt\frac{e^{-zt}}{t^{3}}.$ (47)

To determine DAs (1)-(4) at a scale different from $m_{c}$ one should use
relations (29), (30) and than evolution equations for the corresponding DAs.

##  VI Conclusion.

In this paper we have considered the leading twist distribution amplitudes
(DA) of $P$-wave nonrelativistic mesons. At the leading order approximation in
relative velocity of quark-antiquark pair inside the mesons these functions
can be expressed through one universal DA. We have derived the relations
between the DAs of $P$-wave nonrelativistic mesons and the universal DA.

As an example, we have considered the DAs of the $P$-wave charmonia mesons.
Within QCD sum rules we found the moments of the leading twist DA of $h_{c}$
meson, what allowed us to build the model of the universal DA for the $P$-wave
charmonia mesons.

In addition, we have found the relations between the moments and the
nonrelativistic QCD matrix elements that control relativistic corrections to
any amplitude involving $P$-wave charmonia. The calculation shows that
characteristic size of these corrections is $\sim 30\%$.

This work was partially supported by Russian Foundation of Basic Research
under grant 07-02-00417. The work of V. Braguta was partially supported by
CRDF grant Y3-P-11-05 and president grant MK-2996.2007.2. The work of A.
Luchinsky was partially supported by president grant MK-110.2008.2 and Russian
Science Support Foundation.

## References

  * (1) G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * (2) V. L. Chernyak and A. R. Zhitnitsky, Phys. Rept. 112, 173 (1984).
  * (3) G. T. Bodwin, D. Kang and J. Lee, Phys. Rev. D 74, 114028 (2006) [arXiv:hep-ph/0603185].
  * (4) J. P. Ma and Z. G. Si, Phys. Lett. B 647, 419 (2007) [arXiv:hep-ph/0608221].
  * (5) V. V. Braguta, A. K. Likhoded and A. V. Luchinsky, Phys. Lett. B 646, 80 (2007) [arXiv:hep-ph/0611021].
  * (6) V. V. Braguta, Phys. Rev. D 75, 094016 (2007) [arXiv:hep-ph/0701234].
  * (7) V. V. Braguta, arXiv:0709.3885 [hep-ph].
  * (8) H. M. Choi and C. R. Ji, Phys. Rev. D 76, 094010 (2007) [arXiv:0707.1173 [hep-ph]].
  * (9) T. Feldmann and G. Bell, arXiv:0711.4014 [hep-ph].
  * (10) G. Bell and T. Feldmann, JHEP 0804, 061 (2008) [arXiv:0802.2221 [hep-ph]].
  * (11) G. T. Bodwin and A. Petrelli, Phys. Rev. D 66, 094011 (2002) [arXiv:hep-ph/0205210].
  * (12) E. Braaten and J. Lee, Phys. Rev. D 67, 054007 (2003) [arXiv:hep-ph/0211085];
  * (13) G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125 (1995) [Erratum-ibid. D 55, 5853 (1997)] [arXiv:hep-ph/9407339].
  * (14) G. T. Bodwin, D. Kang and J. Lee, Phys. Rev. D 74, 014014 (2006) [arXiv:hep-ph/0603186].
  * (15) S. J. Brodsky, T. Huang and G. P. Lepage, In *Banff 1981, Proceedings, Particles and Fields 2*, 143-199.
  * (16) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 385 (1979).
  * (17) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 448 (1979).
  * (18) L. J. Reinders, H. Rubinstein and S. Yazaki, Phys. Rept. 127, 1 (1985).
  * (19) W. Buchmuller and S. H. H. Tye, Phys. Rev. D 24, 132 (1981).
  * (20) E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane and T. M. Yan, Phys. Rev. D 17, 3090 (1978) [Erratum-ibid. D 21, 313 (1980)].


wavelet bases $\Psi_{j}$ for the subspaces $\mathsf{W}_{j}$ by selecting a
subset of the columns of the band pass operator
$[\mathbf{I}_{\mathsf{V}_{j}}-\Phi_{j+1}\Phi_{j+1}^{\dagger}]_{\Phi_{j}}^{\Phi_{j}}$
which is the orthogonal projector on the complement of $\mathsf{V}_{j+1}$ into
$\mathsf{V}_{j}$. The wavelets capture the detail lost from going from
$\mathsf{V}_{j}$ to $\mathsf{V}_{j+1}$. As our framework falls into bi-
orthogonal scope, we need to build the dual wavelet bases. For a scale $j$, we
have a wavelet base $\Psi_{j}$ and we need to construct the associated dual
base $\widehat{\Psi}_{j}$ obtained such that we have the following relation

$\widehat{\Psi}_{j}^{\dagger}\Psi_{j}=\mathbf{I}.$

where $\widehat{\Psi}_{j}^{\dagger}$ is the pseudo-inverse of $\Psi_{j}$. We
mention that in [22], biorthogonal wavelet transform is proposed but only the
scaling functions are actually defined. Our construction is identical but we
also propose a definition for wavelets.

####  8.2.4. A generalization using more arbitrary scaling operators

We propose a generalization of the diffusion wavelet framework that enables it
to be combined for instance with the wavelet construction presented in Section
8.2. The idea is merely to replace the powers of $\mathbf{T}$ as
multiresolution scaling operators by more arbitrary low-pass filters. More
precisely, where the operator $\mathbf{T}^{2^{j}}$ is used to define the
approximation space $\mathsf{V}_{j+1}$ in section 8.2.3, we propose to use
instead a low-pass graph filter $\mathbf{H}_{j}$. If the graph filters
$\mathbf{H}_{j}$ correspond to the scaling operators defined in section 8.1
with appropriately increasing scales, this approach provides a way to reduce
the redundancy of the sets of scaling functions. Similarly, the same approach
could be used to reduce the redundancy of the wavelet functions defined in
Section 8.2.

###  8.3. Applications

This section is devoted to illustrate our multiresolution analysis framework
through the visualization of wavelets and scaling functions on regular-type
structures and semi-supervised learning on directed graphs. The application
framework will be the same as in section 7.

####  8.3.1. Multiresolution analysis on the directed cycle graph

We show an example of multi-resolution analysis on the directed cycle graph
$\mathcal{C}_{N}$ with $N=256$. We use the same assumptions made in section
6.4.1. We construct both orthogonal and biorthogonal multi-resolution analyses
on $\mathcal{C}_{N}$ through the framework of the diffusion wavelets applied
on the dyadic powers of $\mathbf{T}$, i.e.
$\\{\mathbf{T}^{2^{j}}\\}_{j=1}^{J}$. We set the number of scales at $J=6$.

Figure 10. Orthogonal and biorthogonal scaling functions on the directed cycle
graph $\mathcal{C}_{256}$. Figure 11. Orthogonal and biorthogonal wavelet
functions on the directed cycle graph $\mathcal{C}_{256}$.

Figure 10 shows the orthogonal and biorthogonal scaling functions at scales
1,3 and 5. At each scale, we represent 3 or 4 scaling functions.

Figure 11 shows the orthogonal and biorthogonal wavelet functions also at
scales 1,3 and 5.

We note at each scale that the support of the orthogonal scaling functions is
larger than the support of the biorthogonal scaling functions. We also note at
scale 3 and 5 that biorthogonal wavelet functions have support slightly
smaller than the orthogonal wavelets functions. Furthermore, we note that
orthogonal scaling functions as well orthogonal wavelet functions have more
oscillations than their biorthogonal counterparts. That means that orthogonal
scaling and wavelet functions have a poorer frequency localization. The
conclusions about the spatial localization of the scaling functions are
identical to the conclusions in [24]. The novelty here is we work on the
directed circle graph and we show the orthogonal and biorthogonal wavelets
functions in addition to the orthogonal and biorthogonal scaling functions.
Finally, the scaling and wavelets functions are not centered around the vertex
where they have been selected. They propagate towards the parents’ node of the
selected vertex. We can also note the orthogonal wavelet transform is more
robust with respect to perturbations than the biorthogonal wavelet transform.
Let $\mathcal{W}_{ortho}=[\Psi_{1},\ldots,\Psi_{5},\Phi_{5}]$ be the inverse
wavelet transform and $\mathcal{W}_{biorth}$ its biorthogonal counterpart. We
have that the condition numbers $\kappa(\mathcal{W}_{ortho})\approx 33$ and
$\kappa(\mathcal{W}_{biortho})\approx 2\times 10^{4}$. Consequently, the
biorthogonal wavelet transform is much more sensitive to perturbations than
its orthogonal counterpart.

####  8.3.2. Multiresolution analysis on the directed Watts-Strogatz graph

In this section, we show examples of multi-resolution analyses on the directed
version of a graph from the Watts-Strogatz model. The Watts-Strogatz model
[85] is a undirected random graph model exhibiting small world properties
including short average path lengths and high clustering [86]. The
construction of a directed graph from the Watts-Strogatz model starts with a
directed cycle graph with $N$ vertices. Each node is connected to its $k$ next
nodes following the direction of the directed cycle graph. For the sake of
simplicity, we consider the $k$ next nodes connected to a given node $i$
following the direction of the directed cycle graph as its "closest"
neighbors. Starting from an arbitrary vertex, we apply the following procedure
to each vertex in a clockwise manner. At vertex $i$, the edge that connects
$i$ to each of its next nodes is randomly rewired with probability $p$ or
remains untouched with probability $1-p$. This procedure is repeated
cyclically for each successive vertex until the vertex $i$ is selected again.
We denote by $\mathcal{G}\sim\operatorname{DWS}(N,K,\beta)$ a graph
constructed following a directed Watts-Strogatz model with $N$ vertices, $K$
nearest neighbors and rewiring probability $\beta$.

In a first instance, we analyze the scaling functions built from low-pass
filters based on $\bar{\mathbf{T}}_{\alpha}$ at a given scale of
$\mathcal{G}_{1}\sim\operatorname{DWS}(64,2,0)$. That is a special case of the
directed Watts-Strogatz model with no rewiring. In a second instance, we
analyze some scaling functions built from low-pass filters based on
$\bar{\mathbf{T}}_{\alpha}$ at a given scale of
$\mathcal{G}_{2}\sim\operatorname{DWS}(64,2,0.02)$. Lastly, we show both
orthogonal and biorthogonal scaling and wavelets functions built following the
construction Section 8.2. We summarize the following cases in table 6.

Case Graph Reference operators Case 1 : Scaling functions built wrt. sect. 8.1
$\mathcal{G}_{1}\sim\operatorname{DWS}(64,2,0)$
$\bar{\mathbf{T}}_{\alpha}\in\bar{\mathcal{T}}$ Case 2: Scaling functions
built wrt. sect. 8.1 $\mathcal{G}_{2}\sim\operatorname{DWS}(64,2,0.02)$
$\bar{\mathbf{T}}_{\alpha}\in\bar{\mathcal{T}}$ Case 3: Scaling functions
built wrt. sect. 8.2 $\mathcal{G}_{2}\sim\operatorname{DWS}(64,2,0.02)$
$\bar{\mathbf{T}}$

Table 6. Table of the different cases for the directed Watts Strogatz graph.

##### Case 1

In this case, we analyze some scaling functions built from low-pass filters
based on $\bar{\mathbf{P}}_{\alpha}$ at scale $t=2^{4}$ of a directed Watts-
Strogatz graph $\mathcal{G}_{2}\sim\operatorname{DWS}(64,2,0)$. We consider
the following low-pass filters

(33)
$\bar{\mathbf{T}}_{\alpha}=\boldsymbol{\Pi}^{1/2}\bar{\mathbf{P}}_{\alpha}\boldsymbol{\Pi}^{-1/2},\quad\forall\alpha\in[0,1].$

and we construct a collection of low-pass graph filters at a given scale
following the construction presented in section 8.1. More precisely, we build
the filters $\mathbf{H}_{\alpha}$ as follows

$\mathbf{H}_{\alpha}=\sum_{\omega\in\boldsymbol{\omega}}h(t\omega)\mathbf{S}_{\omega,\alpha},$

with $t=2^{4}$ and $h(x)=\exp(-x)$.

Figure 12. (Case 1) $50^{th}$ scaling function at scale t=16 on the graph
$\mathcal{G}\sim\operatorname{DWS}(64,2,0)$, for all $\alpha\in[0,1]$.

Figure 12 represents the $50^{th}$ scaling function
$\mathfrak{h}_{\alpha,50}=\mathbf{H}_{\alpha}\boldsymbol{\delta}_{50}$ at
scale $t=2^{4}$ built from low-pass filters $\mathbf{T}_{\alpha}$, for all
$\alpha\in[0,1]$ according the construction section 8.1 on the directed Watts-
Strogatz graph $\mathcal{G}\sim\operatorname{DWS}(64,2,0)$. The particularity
of this case is that the following graph has no rewiring. Consequently, the
associated adjacency matrix is circulant. All the low-pass filters
$\mathbf{T}_{\alpha}$ admit the same discrete Fourier basis and the associated
scaling functions $\mathfrak{h}_{\alpha,50}$ are exactly the same. We observe
that the scaling function is centered around the node 50 and has a symmetric
shape. In this case, for non symmetric graph filters $\mathbf{H}_{\alpha}$,
the scaling functions have symmetric shape and are centered around its given
node.

##### Case 2

In this case, we analyze some scaling functions built from low-pass filters
based on $\bar{\mathbf{P}}_{\alpha}$ at scale 10 of
$\mathcal{G}_{2}\sim\operatorname{DWS}(64,2,0.02)$. We also consider the
following low-pass filters

(34)
$\bar{\mathbf{T}}_{\alpha}=\boldsymbol{\Pi}^{1/2}\bar{\mathbf{P}}_{\alpha}\boldsymbol{\Pi}^{-1/2},\quad\forall\alpha\in\\{0,0.5,1\\}.$

and we construct a collection of low-pass graph filters at a given scale
following the construction presented in section 8.1. More precisely, we build
the filters $\mathbf{H}_{\alpha}$ as follows

$\mathbf{H}_{\alpha}=\sum_{\omega\in\boldsymbol{\omega}}h(t\omega)\mathbf{S}_{\omega,\alpha}$

with $t=2^{4}$ and $h(x)=\exp(-x)$. We observe the $50^{th}$ scaling functions
at the scale $t=2^{4}$, that is
$\mathfrak{h}_{\alpha,50}=\mathbf{H}_{\alpha}\boldsymbol{\delta}_{50}$.

Figure 13. $50^{th}$ scaling function at scale 4 on a graph
$\mathcal{G}\sim\operatorname{DWS}(64,2,0.02)$, $\alpha\in\\{0,0.5,1\\}$, eq.
33.

Figure 13 represent the $50^{th}$ scaling function at the scale $t=2^{4}$, for
different $\alpha$, $\mathfrak{h}_{\alpha,50}$. The graph $\mathcal{G}_{2}$
admits a directed edge from the node $50$ to the node $11$ and a directed edge
from the node $30$ to the node $50$. As $\alpha$ increases towards $1$, we
observe that the scaling function $\mathfrak{h}_{\alpha,50}$ propagates around
the child node $11$. This means the more $\alpha$ increases, the more the
influence of the backward random walk $\mathbf{P}^{*}$ is important. Wee note
for $\alpha=0.5$, the scaling function $\mathfrak{h}_{\alpha,50}$ is centered
around the child node $11$ an the parent node $30$. This case seems to be the
equilibrium. Finally, as $\alpha$ goes to 0, we observe that
$\mathfrak{h}_{\alpha,50}$ propagates around the parent node $30$.

##### Case 3

Here, we consider orthogonal and biorthogonal scaling functions on a graph
$\mathcal{G}\sim\operatorname{DWS}(64,2,0.02)$ used on the Case 2. These
orthogonal and biorthogonal scaling functions are built from the diffusion
wavelets procedure. We start the procedure with the low-pass filter
$\mathbf{H}=\mathbf{T}^{2}$ and we look at the dyadic powers of $\mathbf{H}$
i.e. $\\{\mathbf{H}^{2^{j}}\\}_{j=1}^{J}$. We set the number of scales at
$J=5$. We observe the orthogonal an biorthogonal scaling functions at node
$49$ at scale $3$ obtained by the diffusion wavelets procedure and we compare
to the scaling function at node $49$ built from a graph filter based on
$\bar{\mathbf{T}}$ at the scale $t=16$.

Figure 14. Orthogonal and bi-orthogonal scaling functions built w.r.t the
diffusion wavelet framework versus scaling function built w.r.t spectral
wavelets framework.

Figure 14 shows the orthogonal and bi-orthogonal scaling functions at node
$49$ built w.r.t. the diffusion wavelets framework and the scaling function at
node $49$ built w.r.t. the spectral graph wavelets framework. We note that the
biorthogonal scaling function and the scaling function built by the spectral
graph wavelet framework are similar, i.e. well localized around nodes $11$ and
$30$ and localized around $49$ with a large size support. By contrast, the
orthogonal scaling function has a poor spatial localization around the nodes
$11$ and $30$.

####  8.3.3. Semi-supervised learning on directed graphs with
$\ell_{1}$-regularization

We discuss the semi-supervised learning approach in the case of functions over
directed graphs with $\ell_{1}$ regularization on the wavelet coefficients.
Our aim is to show that the performance of the semi-supervised learning
problem with $\ell_{1}$ regularization is competitive to the existing
approaches, e.g. the semi-supervised learning problem studied in section 7. We
discuss the following semi-supervised learning problem, again on the example
of the political blogs dataset [28]. We use the same notations described in
section 7 and work on the subgraph $\mathcal{G}^{\prime}$. The
$\ell_{1}$-regularized semi-supervised learning problem is

(35)
$\boldsymbol{w}^{*}=\underset{\boldsymbol{w}}{\operatorname{argmin}}\>\|\tilde{\boldsymbol{y}}-\mathbf{M}\mathbf{K}\boldsymbol{w}\|_{2}^{2}+\lambda\|\boldsymbol{w}\|_{1},\quad\lambda\in\mathbb{R}^{+}.$

The graph signal $\tilde{\boldsymbol{y}}=\mathbf{M}\boldsymbol{y}$ is the
partially labeled graph signal with $\mathbf{M}=\\{m_{ij}\\}_{\\{1\leq i,j\leq
N\\}}$ is the mask operator, i.e. the diagonal matrix where
$m_{ii}=\mathbb{1}_{v_{i}\in\mathcal{O}}$ where
$\mathcal{O}\subset\mathcal{V}$ is the subset of known labels.

The matrix $\mathbf{K}=(\mathbf{H}_{J},\mathbf{G}_{1},\cdots,\mathbf{G}_{J})$
is the synthesis filter bank. If we set $\mathbf{X}=\mathbf{MK}$, the equation
(35) can be rewritten as

$\boldsymbol{w}^{*}=\underset{\boldsymbol{w}}{\operatorname{argmin}}\>\|\tilde{\boldsymbol{y}}-\mathbf{X}\boldsymbol{w}\|_{2}^{2}+\lambda\|\boldsymbol{w}\|_{1},\quad\lambda\in\mathbb{R}.$

The formulation (35) is identical to the formulation of the problem of signal
restoration with redundant wavelet transforms in [68], except that
$\mathbf{K}$ is the synthesis wavelet transform for functions defined over
directed graphs. Furthermore, previous approaches of semi-supervised learning
on undirected graphs have been investigated using overcomplete graph wavelets
[69] or critically sampled spline graph wavelets [70]. The
$\ell_{1}$-regularized synthesis semi-supervised learning problem is convex
and can be solved efficiently using e.g. proximal splitting methods [71, 72,
73]. From the solution $\boldsymbol{w}^{*}$ of eq. 35, we define the restored
signal $\boldsymbol{f}^{*}$ as

$\boldsymbol{f}^{*}=\operatorname{sign}(\mathbf{K}\boldsymbol{w}^{*}).$

where $\mathcal{L}$ is the directed normalized Laplacian defined at (4) and
the graph bandpass filters $\mathbf{G}_{1},\mathbf{G}_{2}$ are

We compare the performances of formulations (35) and (16). The formulation
(35) requires a synthesis filter bank $\mathbf{K}$. We choose a two scales
synthesis graph filter $\mathbf{K}$. We choose a two scales synthesis graph
filter bank $\mathbf{K}=(\mathbf{H}_{2},\mathbf{G}_{1},\mathbf{G}_{2})$ based
on a heat kernel construction [83]

$\mathbf{H}_{j}=e^{-\mathcal{L}t_{j}},\quad t_{j}=2^{j},\quad j=1,2.$

with $\mathcal{L}$ is the directed Laplacian defined at (4) and the graph
bandpass filters $\mathbf{G}_{1},\mathbf{G}_{2}$ as

$\mathbf{G}_{2}=\mathbf{H}_{1}-\mathbf{H}_{2},\quad\mathbf{G}_{1}=\mathbf{I}-\mathbf{H}_{1}.$

Figure 15. SSL on the largest connected graph of US political blogs.

Figure 15 shows the performance of semi-supervised learning between the
formulation (35) and (16) on the largest strongly connected graph obtained
from US political blogs. The performance is obtained by averaging 200
realizations and determining the value of the parameters $\lambda$ and
$\gamma$ associated respectively to the approaches (35) and (16) giving the
best performance rates. On fig. 15, the performance based on the $\ell_{1}$
penalization of redundant graph wavelets coefficients on directed graphs is
competitive against the problem with $\ell_{2}$ Dirichlet regularization term
for all percentage of known labels. Normally, one should expect better
performances in the $\ell_{1}$ case but it seems not the case because our
graph wavelets are not good in terms of spatial localization. Still, in order
to improve the performance of the method (35), it could be wise to build new
wavelets with a better space-frequency characterization, and choose the best
suitable number of scales with respect to the characteristics of the data set.

##  9\. Conclusion

We introduced a novel harmonic analysis on directed graphs. First, we proposed
a frequency analysis for functions defined on directed graphs based on the
eigenvectors of the random walk operator on a directed graph. From this
Fourier-type frequency interpretation, we showed how to construct redundant
wavelets on directed graphs as well as critically sampled wavelets by
generalizing the diffusion wavelets framework. Finally, we illustrated our
harmonic analysis through examples of semi-supervised learning and graph
signal modeling on directed graphs and showed the relevance of our framework
to existing approaches.

##  10\. Acknowledgements

The authors would like to thank Romain Cosentino, Randall Balestriero, Yanis
Bahroun, Benjamin Girault and Antonio Ortega for the useful discussions.

## References

  * [1] Steve Butler.  Interlacing for weighted graphs using the normalized laplacian.  Electronic Journal of Linear Algebra, 16(1):8, 2007.
  * [2] Fan Chung.  Laplacians and the cheeger inequality for directed graphs.  Annals of Combinatorics, 9(1):1–19, 2005.
  * [3] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288, 1996.
  * [4] Franziska Göbel, Gilles Blanchard, and Ulrike von Luxburg.  Construction of tight frames on graphs and application to denoising.  In Handbook of Big Data Analytics, pages 503–522. Springer, 2018\.
  * [5] Thierry Coulhon, Gerard Kerkyacharian, and Pencho Petrushev.  Heat kernel generated frames in the setting of dirichlet spaces.  Journal of Fourier Analysis and Applications, 18(5):995–1066, 2012\.
  * [6] D. B. H. Tay and A. Ortega.  Bipartite graph filter banks: Polyphase analysis and generalization.  IEEE Transactions on Signal Processing, 65(18):4833–4846, Sept 2017\.
  * [7] Antonio Ortega, Pascal Frossard, Jelena Kovačević, José MF Moura, and Pierre Vandergheynst.  Graph signal processing.  arXiv preprint arXiv:1712.00468, 2017.
  * [8] David I Shuman, Sunil K Narang, Pascal Frossard, Antonio Ortega, and Pierre Vandergheynst.  The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains.  IEEE Signal Processing Magazine, 30(3):83–98, 2013.
  * [9] Aliaksei Sandryhaila and José MF Moura.  Discrete signal processing on graphs.  IEEE transactions on signal processing, 61(7):1644–1656, 2013.
  * [10] Pierre Brémaud.  Markov chains: Gibbs fields, Monte Carlo simulation, and queues, volume 31.  Springer Science & Business Media, 2013.
  * [11] David Aldous and Jim Fill.  Reversible markov chains and random walks on graphs, 2002.
  * [12] László Lovász.  Random walks on graphs.  Combinatorics, Paul erdos is eighty, 2(1-46):4, 1993.
  * [13] James Allen Fill.  Eigenvalue bounds on convergence to stationarity for nonreversible markov chains, with an application to the exclusion process.  The annals of applied probability, pages 62–87, 1991.
  * [14] Thierry Coulhon and Alexander Grigoryan.  Random walks on graphs with regular volume growth.  Geometric and Functional Analysis, 8(4):656–701, 1998.
  * [15] Audrey Terras.  Fourier analysis on finite groups and applications, volume 43.  Cambridge University Press, 1999.
  * [16] MI Ostrovskii.  Sobolev spaces on graphs.  Quaestiones Mathematicae, 28(4):501–523, 2005.
  * [17] A Kaveh and H Rahami.  A unified method for eigendecomposition of graph products.  International Journal for Numerical Methods in Biomedical Engineering, 21(7):377–388, 2005.
  * [18] Ali Kaveh.  Optimal analysis of structures by concepts of symmetry and regularity.  Springer, 2013.
  * [19] M Maggioni and HN Mhaskar.  Diffusion polynomial frames on metric measure spaces.  Applied and Computational Harmonic Analysis, 24(3):329–353, 2008\.
  * [20] David K Hammond, Pierre Vandergheynst, and Rémi Gribonval.  Wavelets on graphs via spectral graph theory.  Applied and Computational Harmonic Analysis, 30(2):129–150, 2011\.
  * [21] Ronald R Coifman and Mauro Maggioni.  Diffusion wavelets.  Applied and Computational Harmonic Analysis, 21(1):53–94, 2006\.
  * [22] Mauro Maggioni, James C Bremer, Ronald R Coifman, and Arthur D Szlam.  Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs.  In Wavelets XI, volume 5914, page 59141M. International Society for Optics and Photonics, 2005.
  * [23] Ole Christensen et al.  An introduction to frames and Riesz bases.  Springer, 2016.
  * [24] James C Bremer, Ronald R Coifman, Mauro Maggioni, and Arthur D Szlam.  Diffusion wavelet packets.  Applied and Computational Harmonic Analysis, 21(1):95–112, 2006\.
  * [25] Ali Civril.  Column Subset Selection for approximating data matrices.  PhD thesis, Rensselaer Polytechnic Institute, 2009.
  * [26] Aliaksei Sandryhaila and Jose MF Moura.  Discrete signal processing on graphs: Frequency analysis.  IEEE Trans. Signal Processing, 62(12):3042–3054, 2014.
  * [27] Dengyong Zhou, Jiayuan Huang, and Bernhard Schölkopf.  Learning from labeled and unlabeled data on a directed graph.  In Proceedings of the 22nd international conference on Machine learning, pages 1036–1043. ACM, 2005.
  * [28] Lada A Adamic and Natalie Glance.  The political blogosphere and the 2004 us election: divided they blog.  In Proceedings of the 3rd international workshop on Link discovery, pages 36–43. ACM, 2005.
  * [29] Amy N Langville and Carl D Meyer.  Google’s PageRank and beyond: The science of search engine rankings.  Princeton University Press, 2011.
  * [30] Stefania Sardellitti, Sergio Barbarossa, and Paolo Di Lorenzo.  On the graph fourier transform for directed graphs.  IEEE Journal of Selected Topics in Signal Processing, 11(6):796–811, 2017.
  * [31] Rasoul Shafipour, Ali Khodabakhsh, Gonzalo Mateos, and Evdokia Nikolova.  A digraph fourier transform with spread frequency components.  arXiv preprint arXiv:1705.10821, 2017.
  * [32] Frank Bauer.  Normalized graph laplacians for directed graphs.  Linear Algebra and its Applications, 436(11):4193–4222, 2012.
  * [33] Delio Mugnolo.  Semigroup methods for evolution equations on networks.  Springer.
  * [34] Fan RK Chung.  Spectral graph theory.  Number 92. American Mathematical Soc., 1997.
  * [35] David A Levin and Yuval Peres.  Markov chains and mixing times, volume 107.  American Mathematical Soc., 2017.
  * [36] Stephane G Mallat.  A theory for multiresolution signal decomposition: the wavelet representation.  IEEE transactions on pattern analysis and machine intelligence, 11(7):674–693, 1989.
  * [37] Matthias Hein, Jean-Yves Audibert, and Ulrike von Luxburg.  Graph laplacians and their convergence on random neighborhood graphs.  Journal of Machine Learning Research, 8(Jun):1325–1368, 2007.
  * [38] Mikhail Belkin and Partha Niyogi.  Convergence of laplacian eigenmaps.  In Advances in Neural Information Processing Systems, pages 129–136, 2007.
  * [39] Mikhail Belkin, Irina Matveeva, and Partha Niyogi.  Regularization and semi-supervised learning on large graphs.  In International Conference on Computational Learning Theory, pages 624–638. Springer, 2004.
  * [40] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural computation, 15(6):1373–1396, 2003.
  * [41] Yonathan Aflalo and Ron Kimmel.  Spectral multidimensional scaling.  Proceedings of the National Academy of Sciences, 110(45):18052–18057, 2013.
  * [42] Ronald R Coifman and Stéphane Lafon.  Diffusion maps.  Applied and computational harmonic analysis, 21(1):5–30, 2006.
  * [43] Peter Schröder and Wim Sweldens.  Spherical wavelets: Efficiently representing functions on the sphere.  In Proceedings of the 22nd annual conference on Computer graphics and interactive techniques, pages 161–172. ACM, 1995.
  * [44] J-P Antoine and Pierre Vandergheynst.  Wavelets on the 2-sphere: A group-theoretical approach.  Applied and Computational Harmonic Analysis, 7(3):262–291, 1999\.
  * [45] Y Wiaux, JD McEwen, P Vandergheynst, and O Blanc.  Exact reconstruction with directional wavelets on the sphere.  Monthly Notices of the Royal Astronomical Society, 388(2):770–788, 2008.
  * [46] Nick Kingsbury.  Complex wavelets for shift invariant analysis and filtering of signals.  Applied and computational harmonic analysis, 10(3):234–253, 2001\.
  * [47] Emmanuel J Candès and David L Donoho.  New tight frames of curvelets and optimal representations of objects with piecewise c2 singularities.  Communications on pure and applied mathematics, 57(2):219–266, 2004\.
  * [48] Gabriel Peyré and Stéphane Mallat.  Orthogonal bandelet bases for geometric images approximation.  Communications on Pure and Applied Mathematics, 61(9):1173–1212, 2008.
  * [49] Wim Sweldens.  The lifting scheme: A custom-design construction of biorthogonal wavelets.  Applied and computational harmonic analysis, 3(2):186–200, 1996\.
  * [50] Albert Cohen, Ingrid Daubechies, and J-C Feauveau.  Biorthogonal bases of compactly supported wavelets.  Communications on pure and applied mathematics, 45(5):485–560, 1992\.
  * [51] Daryl Geller and Azita Mayeli.  Continuous wavelets on compact manifolds.  Mathematische Zeitschrift, 262(4):895, 2009.
  * [52] Bin Dong.  Sparse representation on graphs by tight wavelet frames and applications.  Applied and Computational Harmonic Analysis, 42(3):452–479, 2017\.
  * [53] Mark Crovella and Eric Kolaczyk.  Graph wavelets for spatial traffic analysis.  In INFOCOM 2003. Twenty-Second Annual Joint Conference of the IEEE Computer and Communications. IEEE Societies, volume 3, pages 1848–1857. IEEE, 2003.
  * [54] Maarten Jansen, Guy P Nason, and Bernard W Silverman.  Multiscale methods for data on graphs and irregular multidimensional situations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(1):97–125, 2009.
  * [55] Raif Rustamov and Leonidas J Guibas.  Wavelets on graphs via deep learning.  In Advances in neural information processing systems, pages 998–1006, 2013.
  * [56] Arthur D Szlam, Mauro Maggioni, Ronald R Coifman, and James C Bremer.  Diffusion-driven multiscale analysis on manifolds and graphs: top-down and bottom-up constructions.  In Wavelets XI, volume 5914, page 59141D. International Society for Optics and Photonics, 2005.
  * [57] Matan Gavish, Boaz Nadler, and Ronald R Coifman.  Multiscale wavelets on trees, graphs and high dimensional data: Theory and applications to semi supervised learning.  In ICML, pages 367–374, 2010.
  * [58] Raif M Rustamov.  Average interpolating wavelets on point clouds and graphs.  arXiv preprint arXiv:1110.2227, 2011.
  * [59] Franziska Göbel, Gilles Blanchard, and Ulrike von Luxburg.  Construction of tight frames on graphs and application to denoising.  arXiv preprint arXiv:1408.4012, 2014.
  * [60] Nicolas Tremblay and Pierre Borgnat.  Graph wavelets for multiscale community mining.  IEEE Transactions on Signal Processing, 62(20):5227–5239, 2014\.
  * [61] Nicolas Tremblay and Pierre Borgnat.  Subgraph-based filterbanks for graph signals.  IEEE Transactions on Signal Processing, 64(15):3827–3840, 2016\.
  * [62] Sunil K Narang and Antonio Ortega.  Perfect reconstruction two-channel wavelet filter banks for graph structured data.  IEEE Transactions on Signal Processing, 60(6):2786–2799, 2012.
  * [63] David I Shuman, Mohammad Javad Faraji, and Pierre Vandergheynst.  A multiscale pyramid transform for graph signals.  IEEE Transactions on Signal Processing, 64(8):2119–2134, 2016.
  * [64] Luca Avena, Fabienne Castell, Alexandre Gaudillière, and Clothilde Mélot.  Intertwining wavelets or multiresolution analysis on graphs through random forests.  arXiv preprint arXiv:1707.04616, 2017.
  * [65] Hrushikesh N Mhaskar.  A unified framework for harmonic analysis of functions on directed graphs and changing data.  Applied and Computational Harmonic Analysis, 2016.
  * [66] Jeff Irion and Naoki Saito.  Hierarchical graph laplacian eigen transforms.  JSIAM Letters, 6:21–24, 2014.
  * [67] Evarist Giné, Vladimir Koltchinskii, et al.  Empirical graph laplacian approximation of laplace–beltrami operators: Large sample results.  In High dimensional probability, pages 238–259. Institute of Mathematical Statistics, 2006.
  * [68] Ivan W Selesnick and Mario AT Figueiredo.  Signal restoration with overcomplete wavelet transforms: Comparison of analysis and synthesis priors.  In Wavelets XIII, volume 7446, page 74460D. International Society for Optics and Photonics, 2009.
  * [69] David I Shuman, Mohammadjavad Faraji, and Pierre Vandergheynst.  Semi-supervised learning with spectral graph wavelets.  In Proceedings of the International Conference on Sampling Theory and Applications (SampTA), number EPFL-CONF-164765, 2011.
  * [70] Venkatesan N Ekambaram, Giulia Fanti, Babak Ayazifar, and Kannan Ramchandran.  Wavelet-regularized graph semi-supervised learning.  In Global Conference on Signal and Information Processing (GlobalSIP), 2013 IEEE, pages 423–426. IEEE, 2013.
  * [71] Patrick L Combettes and Jean-Christophe Pesquet.  Proximal splitting methods in signal processing.  In Fixed-point algorithms for inverse problems in science and engineering, pages 185–212. Springer, 2011.
  * [72] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences, 2(1):183–202, 2009.
  * [73] Stephen Becker, Jérôme Bobin, and Emmanuel J Candès.  Nesta: A fast and accurate first-order method for sparse recovery.  SIAM Journal on Imaging Sciences, 4(1):1–39, 2011.
  * [74] Inam Ur Rahman, Iddo Drori, Victoria C Stodden, David L Donoho, and Peter Schröder.  Multiscale representations for manifold-valued data.  Multiscale Modeling & Simulation, 4(4):1201–1232, 2005.
  * [75] Benjamin Girault, Antonio Ortega, and Shrikanth Narayanan.  Irregularity-aware graph fourier transforms.  arXiv preprint arXiv:1802.10220, 2018.
  * [76] Aamir Anis, Akshay Gadde, and Antonio Ortega.  Efficient sampling set selection for bandlimited graph signals using graph spectral proxies.  IEEE Transactions on Signal Processing, 64(14):3775–3789, 2016\.
  * [77] Ravi Montenegro, Prasad Tetali, et al.  Mathematical aspects of mixing times in markov chains.  Foundations and Trends® in Theoretical Computer Science, 1(3):237–354, 2006.
  * [78] Marc Barthélemy.  Spatial networks.  Physics Reports, 499(1-3):1–101, 2011.
  * [79] Alain Barrat, Marc Barthelemy, and Alessandro Vespignani.  Dynamical processes on complex networks.  Cambridge university press, 2008.
  * [80] Steven Rosenberg.  The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds, volume 31.  Cambridge University Press, 1997.
  * [81] Olivier Lablée.  Spectral Theory in Riemannian Geometry.  2015\.
  * [82] Ulrike Von Luxburg.  A tutorial on spectral clustering.  Statistics and computing, 17(4):395–416, 2007.
  * [83] Martin T Barlow.  Random walks and heat kernels on graphs, volume 438.  Cambridge University Press, 2017.
  * [84] H Francis Song and Xiao-Jing Wang.  Simple, distance-dependent formulation of the watts-strogatz model for directed and undirected small-world networks.  Physical Review E, 90(6):062801, 2014.
  * [85] Duncan J Watts and Steven H Strogatz.  Collective dynamics of "small-world" networks.  nature, 393(6684):440, 1998.
  * [86] MEJ Newman.  Mej newman, siam rev. 45, 167 (2003).  Siam Rev., 45:167, 2003.
  * [87] Harry Sevi, Gabriel Rilling, and Pierre Borgnat.  Multiresolution analysis of functions on directed networks.  In Wavelets and Sparsity XVII, volume 10394, page 103941Q. International Society for Optics and Photonics, 2017.
  * [88] Carlos S Kubrusly.  Hilbert space operators.  In Hilbert Space Operators, pages 13–22. Springer, 2003.
  * [89] Gene H Golub and Charles F Van Loan.  Matrix computations, volume 3.  JHU Press, 2012.
  * [90] Yves Atchadé, Gersende Fort, Eric Moulines, and Pierre Priouret.  Adaptive markov chain monte carlo: theory and methods.  Preprint, 2009.
  * [91] Xiaojin Zhu, John Lafferty, and Ronald Rosenfeld.  Semi-supervised learning with graphs.  PhD thesis, Carnegie Mellon University, language technologies institute, school of computer science Pittsburgh, PA, 2005.