Quantum walks, the discrete wave equation
and Chebyshev polynomials
Simon ApersUniversité Paris Cité, CNRS, IRIF, Paris, France. Email<EMAIL_ADDRESS> and Laurent MicloToulouse School of Economics,
Institut de Mathématiques de Toulouse, CNRS and University of Toulouse. Email<EMAIL_ADDRESS>
A quantum walk is the quantum analogue of a random walk.
While it is relatively well understood how quantum walks can speed up random walk hitting times, it is a long-standing open question to what extent quantum walks can speed up the spreading or mixing rate of random walks on graphs.
In this expository paper, inspired by a blog post by Terence Tao, we describe a particular perspective on this question that derives quantum walks from the discrete wave equation on graphs.
This yields a description of the quantum walk dynamics as simply applying a Chebyshev polynomial to the random walk transition matrix.
This perspective decouples the problem from its quantum origin, and highlights connections to earlier (non-quantum) work and the use of Chebyshev polynomials in random walk theory as in the Varopoulos-Carne bound.
We illustrate the approach by proving a weak limit of the quantum walk dynamics on the lattice.
This gives a different proof of the quadratically improved spreading behavior of quantum walks on lattices.

§ INTRODUCTION

Random walks (and Markov chains, more generally) show up in many areas of science.
A particularly relevant question is how fast random walks spread on graphs.
This is related to topics such as the thermodynamic equilibration of physical systems, the time complexity of Markov chain Monte Carlo algorithms, or the diffusion speed of diseases and rumours in networks.
It is typically quantified by the random walk mixing time <cit.>, which captures the time that it takes for a random walk to converge to its stationary distribution.

Quantum walks are the quantum mechanical analogues of random walks.
They are a key tool in the design of quantum algorithms <cit.>, and they have been a source of interesting mathematics.
Building on a long line of work, it was recently established <cit.> that quantum walks can quadratically speed up the random walk hitting time, i.e., the expected time to hit a subset of vertices when starting from the stationary distribution.
This can be thought of as a graph-theoretic version of the quadratic speedup in Grover's quantum search algorithm <cit.>.
The random walk mixing time is somehow dual to the hitting time, requiring instead to converge to the stationary distribution starting from a single state.
A significantly smaller body of work has been devoted to proving a quadratic quantum walk speedup in the mixing time, paralleling the quadratic speedup in the hitting time.
Such a speedup was proven in a number of relevant but very restricted settings, including settings where the underlying graph is highly symmetric <cit.>, very poorly connected <cit.>, or when we have access to a simulated annealing schedule <cit.>.

Despite advances in these restricted settings, it seems that progress on a general quantum walk speedup in the random walk spreading rate or mixing time has largely stalled.
We revisit this question using an alternative perspective on the quantum walk dynamics, based on Chebyshev polynomials.[The connection between quantum walks and Chebyshev polynomials has been used before (see e.g. <cit.>) but not in the context of quantum walk mixing.]
Rather than following the standard exposition that is usually found in the quantum algorithms literature (see e.g. <cit.>), we present an alternative perspective based on the discrete wave equation on graphs.
This approach is inspired and largely follows a blog post on the discrete wave equation by Terence Tao <cit.> (which, on its turn, was inspired by discussions with Yuval Peres).
It also connects to earlier work on the discrete wave equation on graphs by Friedman and Tillich <cit.>.
Surprisingly, the connection with quantum walks seems to not have been made earlier, and formalizing this is one of the motivations for writing this paper.
As the main technical contribution and proof-of-concept, we use the analysis of Chebyshev polynomials to prove a weak limit on the quantum walk spreading behavior.
In particular, this implies that quantum walks have a ballistic spreading rate on lattices, in correspondence with earlier findings <cit.> that used very different techniques.

Throughout the manuscript, we put care in giving a self-contained exposition, with the hope of rejuvenating interest in the problem and attracting a broader (and potentially non-quantum) audience.

§ RANDOM WALKS FROM THE DIFFUSION EQUATION

As a warm-up, consider the random walk dynamics

\begin{equation} \label{eq:RW}
\fo t\in\ZZ_+,\qquad u(t+1)
= P u(t)
\end{equation}

where $P$ is the random walk transition matrix on a graph.
We will argue how these dynamics correspond to a discrete diffusion equation on graphs.

To this end, consider the diffusion equation
=   12d Δu
over $\real^d$ with $\Delta$ the Laplacian operator.
A discrete space and time analogue of (<ref>) is given by (<ref>) with $P$ the random walk transition matrix on the lattice $\integer^d$.
Indeed, one gets an approximation of (<ref>) by introducing a small parameter $\epsilon>0$ and replacing the time set $\ZZ$ by $\epsilon\ZZ$ and the lattice $\ZZ^d$ by $\sqrt{\epsilon}\ZZ^d$. Consider then the evolution
t∈ϵ_+,    u_ϵ(t+ϵ)   =    P_√(ϵ)[u_ϵ(t)]
where $P_{\sqrt{\epsilon}}$ is the random walk transition matrix on $\sqrt{\epsilon}\integer^d$.
This evolution coincides with (<ref>) for $\epsilon=1$.
Now assume e.g. that the initial condition $u(0)$ of (<ref>) is continuous and has a compact support and that the initial condition of (<ref>) is just the restriction of $u(0)$ to $\sqrt{\epsilon}\ZZ^d$.
Then the following well-known approximation result holds, as a consequence of Donsker's invariance principle <cit.>: for any given $(t,x)\in\RR_+\times\RR^d$, we have
lim_ϵ0_+ u_ϵ(t_ϵ,x_ϵ)   =   u(t,x)where for any $\epsilon>0$, $t_\epsilon\in\epsilon\ZZ_+$ and $x_\epsilon\in \sqrt{\epsilon}\ZZ^d$ are such that $\lim_{\epsilon\ri0_+}t_\epsilon = t$ and $\lim_{\epsilon\ri0_+}x_\epsilon = x$.

This motivates identifying a random walk on the lattice $\integer^d$ with a discrete diffusion equation.
Now, through a small leap of faith, we can generalize this by letting $P$ be the random walk transition matrix on a general graph.
We identify the corresponding dynamics in (<ref>) with the discrete wave equation on that graph.

§ QUANTUM WALKS FROM THE WAVE EQUATION

Now we consider the wave equation instead of the diffusion equation, given by

\[
\ddot{u}
= \f{1}{4d} \Delta u
\]

over $\real^d$.

\begin{align} \label{eq:wave}
\begin{split}
\dot{u} = v, \qquad
\dot{v} = \f{1}{4d} \Delta u.
\end{split}
\end{align}

We will consider a particular[There are alternative options, such as $u(n+1) = u(n) + v(n)$, $v(n+1) = v(n) - (I-P^2) u(n)$.] discrete space and time approximation of the diffusion equation, given by

\begin{align} \label{eq:discrete-wave}
\begin{split}
\fo t\in\ZZ_+,\qquad
u(n+1) &= P u(n) + v(n) \\
v(n+1) &= P v(n) - (I - P^2) u(n)
\end{split}
\end{align}

for $u(n),v(n) \in \ell_2(\ZZ^d)$.
The following exposition largely follows <cit.>.

Discrete approximation.
We can again motivate this approximation by considering the case with $P$ the random walk transition matrix on $\integer^d$.
We introduce a small parameter $\epsilon>0$ and replace the time set $\ZZ$ by $\epsilon\ZZ$ and the lattice $\ZZ^d$ by $\epsilon \ZZ^d$.
Consider then the evolution

\begin{align} \label{eq:discrete-wave-eps}
\begin{split}
\fo t\in\epsilon\ZZ_+,\qquad
u_\epsilon(n+1) &= P_\epsilon u_\epsilon(n) + \epsilon v_\epsilon(n) \\
v_\epsilon(n+1) &= P_\epsilon v_\epsilon(n) - \frac{1}{\epsilon} (I - P_\epsilon^2) u_\epsilon(n)
\end{split}
\end{align}

where $P_{\epsilon}$ is the random walk transition matrix on $\epsilon \integer^d$.
This evolution coincides with (<ref>) for $\epsilon=1$.
Now assume e.g. that the initial condition $(u(0),v(0))$ of (<ref>) is continuous and has a compact support and that the initial condition of (<ref>) is just the restriction of $(u(0),v(0))$ to $\epsilon\ZZ^d$.
Then the following again holds: for any given $(t,x)\in\RR_+\times\RR^d$, we have

\[
\lim_{\epsilon\ri0_+} u_\epsilon(t_\epsilon,x_\epsilon) = u(t,x), \qquad
\lim_{\epsilon\ri0_+} v_\epsilon(t_\epsilon,x_\epsilon) = v(t,x)
\]

where for any $\epsilon>0$, $t_\epsilon\in\epsilon\ZZ_+$ and $x_\epsilon\in \epsilon\ZZ^d$ are such that $\lim_{\epsilon\ri0_+}t_\epsilon = t$ and $\lim_{\epsilon\ri0_+}x_\epsilon = x$.

This argument motivates identifying the dynamics (<ref>) on $\integer^d$ with the discrete wave equation on $\integer^d$.
More generally, we identify (<ref>) with $P$ the random walk transition matrix on a general graph with the discrete wave equation on that graph.

*Unitary dynamics.
It is easy to check that the following “energy” quantity

\[
\| \sqrt{I-P^2} \, u(n) \|^2 + \| v(n) \|^2
\]

is conserved.
Indeed, if we define $w(n) = \sqrt{I-P^2} u(n)$, then we can equivalently write (<ref>) as

\begin{align*}
w(n+1) &= P w(n) + \sqrt{I-P^2} v(n) \\
v(n+1) &= P v(n) - \sqrt{I - P^2} w(n).
\end{align*}

The corresponding propagator

\begin{equation} \label{eq:unitary}
= \begin{bmatrix}
P & \sqrt{I - P^2} \\
- \sqrt{I - P^2} & P
\end{bmatrix}
\end{equation}

is then unitary, and hence conserves the norm or energy $\| w(n) \|^2 + \| v(n) \|^2$.
Such a unitary that encodes an operator $P$ as one of its blocks is called a “unitary dilation” in the mathematics literature.

The $n$-step propagator $U^n$ can be neatly expressed as

\begin{equation} \label{eq:Un}
= \begin{bmatrix}
T_n(P) & \sqrt{I-P^2} \, U_{n-1}(P) \\
- \sqrt{I-P^2} \, U_{n-1}(P) & T_n(P)
\end{bmatrix}.
\end{equation}

where $(T_n)_{n\in\ZZ_+}$ and $(U_n)_{n\in\ZZ_+}$ are respectively the families of Chebyshev polynomials of the first and second kind.

§.§ Quantum walks and block encodings

Arguably, such unitary dynamics are a natural map between quantum states, in which case we propose to call the resulting dynamics a quantum walk on the underlying graph.
In fact, the unitary dilation in (<ref>) is bread and butter in the quantum computing literature – there it is usually called a block encoding <cit.> of the operator $P$.
Such block encodings are central to our current understanding of quantum algorithms – see e.g. the recent “grand unification of quantum algorithms” in the form of the quantum singular value transformation <cit.>.
One of the first (if not the first) demonstrations of such block encoding was in the seminal work on quantum walks by Szegedy <cit.>.
Szegedy constructed a unitary quantum walk operator $W$, different from the propagator $U$ in (<ref>), but resembling it in that

* $W$ is a unitary dilation of the random walk operator $P$, and

* $W^n$ is a unitary dilation of $T_n(P)$.

Szegedy used some additional properties of their particular construction of $W$.
However, we now understand that any operator satisfying the properties 1. and 2. can be used as a quantum walk operator, allowing one to reproduce the algorithmic results on quantum walks referenced earlier <cit.>.
This motivates our interpretation of the discrete wave dynamics as a quantum walk – in fact, it is surprising that this connection seems to not have been made earlier.

We note that, while the $u$- and $v$-components (or $w$- and $v$ components) could be interpreted as position and momentum in the original wave equation, they lose that meaning in the quantum walk interpretation.
Rather, the two components should be interpreted as the components of a state with two chiralities.

§.§ Varopoulos-Carne bound and quantum fast-forwarding

A nice demonstration of this connection is the use of a particular identity for Chebyshev polynomials in random walk theory on the one hand (also discussed in <cit.>), and in quantum algorithms design on the other hand.
The identity is the following: for any self-adjoint contraction $P$ it holds that

\begin{equation} \label{eq:cheb-binom}
= \sum_{k \in \integer} q_n(k) T_{|k|}(P),
\end{equation}

where $q_n(k)$ is the probability that a random walk on the line, starting from the origin, is at position $k$ after $n$ steps.

In random walk theory, this identity is used to prove the Varopoulos-Carne bound (see <cit.> and a nice discussion in <cit.>).
The bound states that for any pair of vertices $x,y$ it holds that

\[
\leq 2 \exp(-d(x,y)^2/(2n)),
\]

where $d(x,y)$ is the shortest path distance between vertices $x$ and $y$.
The main idea of the bound is to apply Hoeffding's bound to the coefficients $q_n(k)$, showing that

\begin{equation} \label{eq:cheb-approx}
\approx \sum_{|k| \leq O(\sqrt{n})} q_n(k) T_{|k|}(P).
\end{equation}

Now notice that $[T_{|k|}(P)](x,y) = 0$ if $d(x,y) > k$, and so $P^n(x,y) \approx 0$ when $d(x,y) > O(\sqrt{n})$.

In quantum walk algorithms, the identities (<ref>) and (<ref>) were used recently in a technique called quantum fast-forwarding <cit.>.
Combining the expression (<ref>) for the quantum walk propagator with the approximation (<ref>), we see that

\begin{equation} \label{eq:trunc-cheb}
\sum_{|k| \leq O(\sqrt{n})} q_n(k) U^{|k|}
\approx \begin{bmatrix} P^n & \star \\ \star & P^n \end{bmatrix},
\end{equation}

where the $\star$ represent non-relevant entries.
Through an additional quantum algorithmic technique called “linear combination of unitaries” <cit.>, it is possible to algorithmically implement a unitary dilation of (<ref>) (and hence of $P^n$) while effectively only implementing $O(\sqrt{n})$ quantum walk steps.
This explains the terminology “quantum fast-forwarding”.

§ QUANTUM WALKS ON LATTICES

The preceding discussion shows that we can understand quantum walks through the lense of Chebyshev polynomials.
We demonstrate this connection by studying such Chebyshev polynomials for the particular case of lattices.
Specifically, we are interested in the dynamics

\[
U^n \begin{bmatrix} e_0 \\ 0 \end{bmatrix}
= \begin{bmatrix} T_n(P) \, e_0 \\ - \sqrt{I-P^2} \, U_{n-1}(P) \, e_0 \end{bmatrix},
\]

where $e_0$ is the indicator on some vertex of the lattice.
This corresponds to the discrete wave equation with initial conditions $u(0) = e_0$ and $v(0) = 0$.
Equivalently, it corresponds to an $n$-step quantum walk starting from $e_0$.
We will be interested in the spreading behavior of the quantum walk, and will focus on the $u$-component of the resulting state, i.e., the vector $T_n(P) e_0$.
This vector induces a (subnormalized) measure $\nu$ on the lattice by setting $\nu_x = \left| [T_n(P)]_{0,x} \right|^2$, and this measure will be the object of our study.
Operationally, it corresponds to the probability of measuring the quantum walk in a certain chirality at vertex $x$ (see discussion at the end of <ref>).

§.§ Quantum walk on the line

On the infinite line, the $T_n(P)$ component of the quantum walk takes a particularly simple form.
Specifically, consider the line with set of vertices $\integer$ and commuting operators $Q_1$ and $Q_1^{-1}$ defined by $Q_1 e_x = e_{x+1}$ and $Q_1 ^{-1} e_{x} = e_{x-1}$, where for any $x\in\ZZ$, $e_x$ is the function on $\ZZ$ taking the value 1 at $x$ and 0 elsewhere.
The random walk transition matrix is $P = (Q_1 + Q_1^{-1})/2$, and the random walk distribution $[P^n]_{0,\cdot}$ after $n$ steps essentially corresponds to a binomial distribution around $0$ with standard deviation $O(\sqrt{n})$.
In contrast, if we use the formula $T_n((x+x^{-1})/2) = (x^n + x^{-n})/2$ we can rewrite the Chebyshev polynomial $T_n(P)$ as

\begin{equation} \label{eq:ch-cycle}
= T_n\left(\frac{Q_1+Q_1^{-1}}{2}\right)
= \frac{Q_1^n + Q_1^{-n}}{2}.
\end{equation}

This implies that $[T_n(P)]_{0,x} = e_n(x)/2 + e_{-n}(x)/2$, and so the quantum walk shows two peaks precisely at a distance $n$ from the origin.
These dynamics represent the ballistic nature of the quantum walk dynamics, as opposed to the diffusive nature of the random walk dynamics.
The hope is that this improved spreading rate persists in more general graphs as well.

If instead of the simple random walk we consider the lazy random walk on the cycle, $P_L = (I + P)/2$, there is no longer a simple expansion of $T_n(P_L)$.
Instead, we retrieve the fickle (but still ballistic) behavior usually associated to the quantum walk on the line.
The behaviors of $P_L^n$ and $T_n(P_L)$ are illustrated in <ref>.

[black, sharp plot, line width = 1.5pt] table [col sep=comma, x=X, y=Y] RW-Z.csv;

[black, sharp plot, line width = 1.5pt] table [col sep=comma, x=X, y=Y] chebyshev-Z.csv;

(left) Probability measure induced by $[P_L^n]_{0,\cdot}$ for $n = 50$, corresponding to an $n$-step lazy random walk on the line. (right) Probability measure induced by $[T_n(P_L)]_{0,\cdot}^2$ for $n = 50$, associated to an $n$-step quantum walk on the line. (left and right) Lines drawn continuously for clarity.

§.§ Quantum walk on the lattice

As a more involved example, consider a random walk on the $2$-dimensional lattice with vertices $V = \{(x,y) \in \integer^2\}$ and commuting generators $Q_1,Q_1^{-1},Q_2,Q_2^{-1}$ defined by $Q_1 e_{x,y} = e_{x+1,y}$ and $Q_2 e_{x,y} = e_{x,y+1}$, and similarly for the inverses.
If we set $P_1 = (Q_1+Q^{-1}_1)/2$ and $P_2 = (Q_2+Q^{-1}_2)/2$, then the random walk transition matrix is given by

\[
= \frac{P_1 + P_2}{2}.
\]

In <ref> we compare the measure induced by the random walk transition matrix $P^n$ and that of the Chebyshev polynomial $T_n(P)$.

colormap=whitebluecolor=(white) color=(black) color=(black)]
3[surf, mesh/rows=71, mesh/cols=71, opacity=0.01,fill opacity=0.5] table RW-Z2.csv;

colormap=whitebluecolor=(white) color=(black) color=(black)]
3[surf, mesh/rows=71, mesh/cols=71, opacity=0.01,fill opacity=0.5] table cheb-Z2.csv;

(left) Probability measure induced by $[P^n]_{0,\cdot}$ for $n = 50$, corresponding to an $n$-step random walk on the lattice. (right) Probability measure induced by $[T_n(P)]_{0,\cdot}^2$ for $n = 50$, associated to an $n$-step quantum walk on the lattice.

To get a handle on the operator $T_n(P)$, we will use the following expansion of a degree-$n$ bivariate Chebyshev polynomial into univariate degree-$n$ Chebyshev polynomials:

\begin{equation} \label{eq:bivar-exp}
T_n\left(\frac{X + Y}{2}\right)
= \sum_{(p,q) \in \lin 0,n\rin} a_{n,p,q} T_p(X) T_q(Y),
\end{equation}

for a set of real coefficients $(a_{n,p,q})_{p,q\in\lin 0,n\rin}$.
Since $T_n(P) = T_n((P_1+P_2)/2)$ is a degree-$n$ polynomial in the commuting operators $P_1$ and $P_2$, we can use the same expansion to rewrite

\begin{align*}
T_n\left(\frac{P_1 + P_2}{2}\right)
&= \sum_{(p,q) \in \lin 0,n\rin} a_{n,p,q} T_p(P_1) T_q(P_2) \\
&= \sum_{(p,q) \in \lin 0,n\rin} a_{n,p,q} \left(\frac{Q^p_1 + Q_1^{-p}}{2}\right)
\left(\frac{Q^{q}_2 + Q_2^{-q}}{2}\right),
\end{align*}

where we reused the fact (<ref>) that $T_p(P_1) = (Q_1^p+Q_1^{-p})/2$.
Now notice that

\[
= a_{n,|p|,|q|}/2^{k(|p|,|q|)},
\]

where $k(|p|,|q|) \in \{0,1,2\}$ denotes the number of nonzero entries among $|p|,|q|$.
We can hence lower bound the probability for a quantum walk, starting from the origin $(0,0)$, to reach a vertex $(p,q)$ by

\begin{equation} \label{eq:qt-lattice}
\left| [T_n(P)]_{(0,0),(p,q)} \right|^2
= a_{n,|p|,|q|}^2/2^{2k(|p|,|q|)}
\geq a_{n,|p|,|q|}^2/16
\end{equation}

(since $k(p,q)\leq 2$).
These coefficients are fully determined by the Chebyshev expansion in (<ref>), and so the problem reduces to bounding these coefficients.
More specifically, we will consider the induced measure $\mu_n$ on $[-1,1]^2$ given by
μ_n       ∑_p,q a_n,|p|,|q|^2δ_(p/n,q/n) ∑_p,q a_n,|p|,|q|^2
where $\delta_{(p/n,q/n)}$ is the Dirac mass at $(p/n,q/n)\in[0,1]^2$.
As our main technical contribution, we prove a weak limit theorem on $\mu_n$.

There exists a continuous measure $\mu$ on $[-1,1]^2$ such that $\mu_n \rightarrow \mu$.

As the most importance consequence of this weak limit, it proves that the quantum walk spreads ballistically on the lattice – after $n$ steps, the expected distance from the origin is $\Omega(n)$.
Using very different machinery, similar conclusions were found in earlier work <cit.>.

§ PROOF OF WEAK CONVERGENCE ON THE LATTICE

Our goal is to investigate the weak convergence for $n$ large of $\mu_n$ toward a law $\mu$ on $[0,1]^2$.
To do so, it will be more convenient to consider the measure
on $[0,1]^2$ given by
γ_n      ∑_p,q∈0, n a_n,p,q^2δ_(p/n,q/n)
We will show there exists a probability measure $\gamma\neq 0$ on $[0,1]^2$ such that  for large $n$, $\gamma_n$ weakly converges toward $\gamma$.
Since for any $n\in\ZZ_+$,
μ_n   =   γ_nγ_n([0,1]^2)we deduce the weak convergence of $\mu_n$ toward $ \mu\df\gamma/\gamma([0,1]^2)=\gamma$.The method will consist in showing the convergence of the (bivariate) even moments of $\gamma_n$ toward those of $\gamma$.
More precisely, for any $K,L\in\ZZ_+$, consider the polynomial function
φ_K,L[0,1]^2∋(x,y)   ↦    x^2Ky^2LWe are to find a measure $\gamma\neq 0$ on $[0,1]^2$ such that for any $K,L\in\ZZ_+$,
lim_n γ_n[φ_K,L]   =   γ[φ_K,L]Taking into account the Stone-Weierstrass theorem (applied on the compact set $[0,1]^2$ with respect to the algebra generated by the $\varphi_{K,L}$, for $K,L\in\ZZ_+$), this convergence is sufficient to deduce the weak convergence of $\gamma_n$ toward $\gamma$.
The proof is split in three main parts corresponding to the following subsections.

§.§ Second moment as an integral

We start by approximating the second moments of $\mu_n$ by some simpler integral.
Consider the function

\begin{equation} \label{eq:def-h}
\coloneqq T_n\lt(\f{\cos(x)+\cos(y)}2\rt)
= \sum_{p,q\in\lin 0,n\rin} a_{n,p,q} \cos(px)\cos(qy),
\end{equation}

where the last equality follows from <ref> and the definition of the Chebyshev polynomials.
We show the following approximation.

We have
lim_n →∞ |
- 1/π^2 n^2(K+L) ∫_[0,2π]^2
( ∂^K_x ∂^L_y h(x,y) )^2 dxdy
|    =    0.

The proof consist of the two following results.

We have
1/π^2 n^2(K+L) ∫( ∂^K_x ∂^L_y h(x,y) )^2 dxdy   =
1/ n^2(K+L) ∑_p,q∈0,n a^2_n,p,qp^2Kq^2L w_K(p)w_L(q)

=   γ_n[φ_K,Lw_K(n·)⊗w_L(n·)]
k∈_+, r∈_+,    w_k(r)      {
2   , if $k$ is even and $r=0$

0   , if $k$ is odd and $r=0$

1   , otherwise.
We rewrite
∂^K_x ∂^L_y h(x,y)
= ∑_p,q∈0,n a_n,p,qp^Kq^L ξ_K(px)ξ_L(qy)where $\xi_r(z)$ denotes the $r$-th derivative $\cos^{(r)}(z)$.
Now recall the orthogonality relations:
p,p'∈_+,    12π∫_0^2π ξ_r(px)ξ_r(p'x)  dx
=    {

1   , if $p=p'=0$ and $r\equiv 0[2]$

1/2   , if $p=p'\geq 1$

0   , otherwise
.Taking the square of (<ref>) and integrating, we get
1(2π)^2∫(∂^K_x ∂^L_y h(x,y))^2dxdy

∑_p,q,p',q'∈0,n a_n,p,qp^Kq^L a_n,p',q'(p')^K(q')^L 1(2π)^2∫ξ_K(px)ξ_L(qy)ξ_K(p'x)ξ_L(q'y)  dxdy

=   ∑_p,q,p',q'∈0,n a_n,p,q a_n,p',q'p^K(p')^Kq^L(q')^L (12π∫ξ_K(px)ξ_K(p'x)  dx)(12π∫ξ_L(qy)ξ_L(q'y)  dy)

=   ∑_p,q∈0,n a^2_n,p,qp^2Kq^2L  (12π∫ξ_K^2(px)  dx)(12π∫ξ_L^2(qy)  dy)

=   14∑_p,q∈0,n a^2_n,p,qp^2Kq^2L w_K(p)w_L(q)
= 14n^2(K+L)γ_n[φ_K,Lw_K(n·)⊗w_L(n·)].
We deduce:

We have for any $n\in\ZZ_+$,
∑_p,q∈0,na^2_n,p,q   ≤   4

Considering the case $K=L=0$ in Lemma <ref>, we get
1π^2∫_[0,2π]^2 T_n^2(cos(x)+cos(y)2)  dxdy
=   ∑_p,q∈0,nw_0(p)w_0(q)a_n,p,q^2

≥    ∑_p,q∈0,na_n,p,q^2
To conclude to the announced bound, note that the l.h.s. is bounded above by $4$, since $T_n(z)$ takes its values in $[-1,1]$ for any $z\in[-1,1]$.The next step enables to end the proof of Proposition <ref>.

We have
lim_n |γ_n[φ_K,Lw_k(n·)⊗w_L(n·)]-γ_n[φ_K,L]|    =    0
To prove the above convergence, observe that for any given $K,L\in\ZZ_+$,
the quantity
|γ_n[φ_K,Lw_k(n·)⊗w_L(n·)]-γ_n[φ_K,L]|   =   |∑_p,q∈_+a^2_n,p,q(pn)^2k(qn)^2L(w_K(p)w_L(q)-1)|converges to zero for large $n\in\NN$.
Let us consider the term
p^2kq^2L(w_K(p)w_L(q)-1)It vanishes if either $p\geq 1$ and $q\geq 1$ or $p=0$ and $K\geq 1$ or $q=0$ and $L\geq 1$.
Thus it can be non zero only in the two following (non-disjoint) alternatives:
$\bullet$ If $p=0$ and $K=0$, then its value is
q^2L(2w_L(q)-1)   =   {
q^2L   , if $q\geq 1$

0   , if $q=0$ and $L\geq 1$

3   , if $q=0$ and $L=0$.$\bullet$ If $q=0$ and $L=0$, then its value is
p^2K(2w_K(p)-1)   =   {
p^2K   , if $p\geq 1$

0   , if $p=0$ and $K\geq 1$

3   , if $p=0$ and $K=0$.It follows that for any $n\in\NN$,

3a^2_n,0,0_{(0,0)}(K,L)+∑_q∈ a^2_n,0,q (qn)^2L_{0}(K)+∑_p∈ a^2_n,p,0 (pn)^2K_{0}(L)

=   3a^2_n,0,0_{(0,0)}(K,L)+∑_q∈ a^2_n,p,0 ((pn)^2K_{0}(L)+(pn)^2L_{0}(K))

≤   3a^2_n,0,0+2∑_q∈ a^2_n,p,0

≤    3∑_p∈_+ a^2_n,p,0
(where the the symmetry in $(p,q)$ of $a_{n,p,q}$ was taken in account in the second equality).
Thus it is sufficient to show that
lim_n ∑_p∈_+ a^2_n,p,0   =   0From the definition of the coefficients $(a_{n,p,q})_{n\in\ZZ_+, p,q\in\lin 0,n\rin}$ we get for any given $n\in \ZZ_+$ and $x\in[0,\pi]$,
12π∫_[0,2π] T_n(cos(x)+cos(y)2)  dy    =   ∑_p,q∈0,n a_n,p,qcos(px)12π∫_[0,2π] cos(qy)  dy

=   12∑_p∈0,na_n,p,0cos(px)Taking the square and integrating, we get
12π∫(12π∫_[0,2π] T_n(cos(x)+cos(y)2)  dy)^2dx    =   a_n,0,0^28+18∑_p∈0,n a_n,p,0^2so that (<ref>) amounts to
lim_n 12π∫_[0,2π](12π∫_[0,2π] T_n(cos(x)+cos(y)2)  dy)^2dx   =   0
or, by symmetry, to
lim_n∫_[0,π](∫_[0,π] T_n(cos(x)+cos(y)2)  dy)^2dx   =   0
Fix $x\in [0,\pi]$. To investigate the integral $\int_{[0,\pi]} T_n\lt(\f{\cos(x)+\cos(y)}{2}\rt)\, dy$ we consider the
change of variable $y\mapsto \theta_x(y)$ where $\theta_x(y)$ is the unique solution belonging to $[0,\pi]$ of the equation
cos(θ_x(y))   =   cos(x)+cos(y)2To simplify notations, from now on we remove the index $x$.
The mapping $\theta$ is a bijection from $[0,\pi]$ to
I       [arccos(cos(x)+12),arccos(cos(x)-12)]and
we have y∈(0,π),    |θ'(y)|   =   12|sin(y)sin(θ(y))|The change of variable formula implies that for any measurable function $f\in \LL^1(\ell_I)$ (where $\ell_I$ is the restriction of the Lebesgue measure on $I$),
12∫_[0,π]f(θ(y))|θ'(y)|  dy    =    ∫_I f(z)  dzTaking into account that for $y\in[0,\pi]$, we have
|sin(y)|   =   √(1-cos^2(y))

=   √(1-(2cos(θ(y))-cos(x))^2)we can write (<ref>) as
12∫_[0,π]f(θ(y))√(1-(2cos(θ(y))-cos(x))^2)|sin(θ(y))|  dy    =    ∫_I f(z)  dzDefine the function $g$ on $I$
z∈I,    g(z)       f(z) √(1-(2cos(z)-cos(x))^2)|sin(z)|to obtain
12∫_[0,π]g(θ(y))  dy    =    ∫_I g(z)|sin(z)|√(1-(2cos(z)-cos(x))^2)  dzequality which is valid as soon as $g\circ \theta\in\LL^1(\ell_{[0,\pi]})$.In particular, taking $g=T_n$, which satisfies this condition, we get
∫_[0,π]T_n(θ(y))  dy    =   2∫_I T_n(z)|sin(z)|√(1-(2cos(z)-cos(x))^2)  dznamely
∫_[0,π]T_n(cos(x)+cos(y)2)  dy    =   2∫_I cos(nz)|sin(z)|√(1-(2cos(z)-cos(x))^2)  dzConsidering this identity with $n=0$, we get that
the mapping
I∋z   ↦   |sin(z)|√(1-(2cos(z)-cos(x))^2)belongs to $\LL^1(\ell_I)$.
Thus we can apply Riemann-Lebesgue theorem to deduce
lim_n ∫_I cos(nz)|sin(z)|√(1-(2cos(z)-cos(x))^2)  dz   =   0and by consequence
lim_n ∫_[0,π]T_n(cos(x)+cos(y)2)  dy    =   0In addition, we have
n∈_+,x∈[0,π],    |∫_[0,π]T_n(cos(x)+cos(y)2)  dy|   ≤    πIt follows that (<ref>) is a consequence of the dominated convergence theorem.

§.§ Negligible boundary region

For technical reasons that appear later, we wish to limit the integral in <ref> to a subdomain $A_\epsilon \subset [0,2\pi]^2$, which for any $\epsilon\in (0,1)$ is defined by
A_ϵ       {(x,y)∈[0,2π]^2cos(x)+cos(y)2≤1-ϵ}.
This is justified since we have
lim_ϵ0_+lim_n |γ_n[φ_K,L]-1/π^2 n^2(K+L) ∫_A_ϵ
( ∂^K_x ∂^L_y h(x,y) )^2 dxdy
|   =   0.This is a consequence of the following lemma, which shows that the contribution of the boundary region $[0,2\pi]^2 \backslash A_\epsilon$ is negligible.

For $\epsilon\in (0,1)$, denote $B_\epsilon\df [0,2\pi]^2\setminus A_\epsilon$.
We have
lim_ϵ0_+sup_n≥1 1/n^2(K+L) ∫_B_ϵ (∑_p,q∈0,n a_n,p,qp^Kq^L ξ_K(px)ξ_L(qy))^2dxdy   =   0

Note that $B_\epsilon \subset B^-_\epsilon\cup B^+_\epsilon$ with
B_ϵ^-       {(x,y)∈[0,2π]^2cos(x)≤-1+2ϵ,  cos(y)≤-1+2ϵ}

B_ϵ^+       {(x,y)∈[0,2π]^2cos(x)≥1-2ϵ,  cos(y)≥1-2ϵ}We will prove the statement for the $B_\epsilon^+$ region, since the proof for $B_\epsilon^-$ is analogous.
Define $\eta$ as the unique solution in $(0,\pi)$ of $\cos(\eta)=1-2\epsilon$, we have that $\eta$ goes to $0_+$ when $\epsilon$ goes to $0_+$, and we can rewrite
∫_B_ϵ^+ (∑_p,q∈0,n a_n,p,qp^Kq^L ξ_K(px)ξ_L(qy))^2dxdy
∫_[-η,η]^2 (∑_p,q∈0,n a_n,p,qp^Kq^L ξ_K(px)ξ_L(qy))^2dxdy

=   ∑_p,q,p',q'∈0,na_n,p,qa_n,p',q'p^Kq^L (p')^K(q')^L ∫_-η^η
ξ_K(px)ξ_K(p'x)  dx ∫_-η^η
ξ_L(qy))ξ_L(q'y)  dy
Let us compute the integrals of the r.h.s.
When $K$ is even, we have
∫_-η^η ξ_K(px)ξ_K(p'x)  dx
=   2∫_0^η cos(px)cos(p'x)  dx

=   ∫_0^η cos((p+p')x)+cos((p-p')x)  dx

=   -sin((p+p')η)p+p'-sin((p-p')η)p-p'
with the usual convention that $\f{\sin((p-p')\eta)}{p-p'}=\eta$ when $p-p'=0$.
When $K$ is odd, a similar computation shows that
∫_-η^η ξ_K(px)ξ_K(p'x)  dx   =   sin((p+p')η)p+p'-sin((p-p')η)p-p'.
It follows that the main integral is a sum of terms of the form
∑_p,q,p',q'∈0,n±a_n,p,qa_n,p',q'p^Kq^L (p')^K(q')^L sin((p±p')η)p±p'sin((q±q')η)q±q',
for different choices of the signs.
We consider the case with all minus signs, as the other cases follow by essentially the same argument.
By a triangle inequality and the Cauchy-Schwarz inequality we can bound the sum by
∑_p,q,p',q'a_n,p,q   a_n,p',q'   |sin((p-p')η)p-p'sin((q-q')η)q-q'|

≤   √( ∑_p,q,p',q'a^2_n,p,q(sin((p-p')η)p-p'sin((q-q')η)q-q')^2)
√( ∑_p,q,p',q'a^2_n,p',q'(sin((p-p')η)p-p'sin((q-q')η)q-q')^2)

=   ∑_p,q,p',q'a^2_n,p,q(sin((p-p')η)p-p'sin((q-q')η)q-q')^2

∑_p,qa^2_n,p,q ∑_r,s∈(sin(rη)rsin(sη)s)^2

=   ∑_p,qa^2_n,p,q (∑_r∈(sin(rη)r)^2)^2

=   ∑_p,qa^2_n,p,q (η^2+2∑_r∈(sin(rη)r)^2)^2

≤    4(η^2+2∑_r∈(sin(rη)r)^2)^2
where we used Corollary <ref>.Dominated convergence now implies that
lim_η0_+∑_r∈(sin(rη)r)^2   =   0,which proves the claimed convergence.

§.§ Simplifying the integral

We will now approximate the integral $\int_{A_\epsilon} \lt(\partial^K_x \partial^L_y h(x,y)\rt)^2dxdy$ by a (arguably) simpler integral.

§.§.§ First simplification

We use Faà di Bruno's formula to rewrite $\partial_x^K h(x,y) = \partial_x^K T_n((\cos(x)+\cos(y))/2)$ as
(cos(x)+cos(y)2)∏_k∈K(ξ_k(x)2k!)^m_k,where the sum is over all nonnegative integers $(m_1, ..., m_K)$ satisfying the constraint $m_1+2m_2+\cdots +Km_K=K$.
Applying Faà di Bruno's formula again we can rewrite $\partial_x^K \partial_y^L h(x,y)$ as

T_n^(∑m_ℓ+ ∑n_ℓ)
(cos(x)+cos(y)2)∏_k∈K(ξ_k(x)2k!)^m_k∏_l∈L(ξ_l(y)2l!)^n_lHow terrible this expression may seem, as $n \to \infty$ the final integral will be dominated by a single term in the summand.
To see this, we use the following bound.

For any $\epsilon\in (0,1)$ and $r\in\NN$, there exists a constant $ C(\epsilon,r)>0$ such that
n∈_+, x∈[-1+ϵ,1-ϵ],    |T_n^(r)(x)|   ≤    C(ϵ,r)n^r
Recall that $T_n(x) = \cos(n\arccos(x))$ for all $x \in [-1,1]$.
The desired result follows from the fact that  any given  derivative of $\arccos$ is bounded on $[-1+\epsilon, 1-\epsilon]$.It follows that generically in (<ref>), as $n \rightarrow \infty$, the main term corresponds to $m_1=K$, $n_1=L$, and all other coefficients equal to zero.
We can prove the following lemma, where we define the integral
1π^2n^2(K+L)∫_A_ϵ  (T_n^(K+L)(cos(x)+cos(y)2)(-sin(x)2)^K(-sin(y)2)^L)^2  dxdy.

For any $\epsilon>0$ we have
lim_n 1π^2n^2(K+L)∫_A_ϵ ( ∂_x^K ∂_y^L h(x,y) )^2  dxdy
- I_K,L(ϵ,n)
=   0.

Using the expression (<ref>) and expanding the square, up to some factor not depending on $n$, we end up with products of the form
T_n^(∑m_ℓ+ ∑n_ℓ)
(cos(x)+cos(y)2)T_n^(∑m'_ℓ+ ∑n'_ℓ)
(cos(x)+cos(y)2)where $\sum_{\ell=1}^K \ell m_\ell = \sum_{\ell=1}^K \ell m'_\ell = K$ and $\sum_{\ell=1}^L \ell n_\ell = \sum_{\ell=1}^L\ell n'_\ell = L$.
According to Lemma <ref>, up to a constant depending on $\epsilon$, $K$ and $L$, the quantity (<ref>) is bounded independently for $(x,y)\in A_\epsilon$ by
n^∑m_ℓ+ m'_ℓ+ ∑n_ℓ+ n'_ℓ
which is negligible with respect to $n^{2(K+L)}$ unless the only nonzero coefficients are $m_1=m'_1=K$ and $n_1=n'_1=L$.
The only potentially non-vanishing term in the integral is hence
( T_n^(K+L)(cos(x)+cos(y)2)(-sin(x)2)^K(-sin(y)2)^L )^2
and this leads to the announced convergence.

§.§.§ Second simplification

We do a second simplification by approximating the integral $I_{K,L}(\epsilon,n)$.
Applying Faà di Bruno's formula to the composition of $T_n$ with cosine,
we get for any $\theta\in[0,\pi]$,
∑_m_1+2m_2+⋯+(K+L)m_K+L=K+L (K+L)!m_1!m_2!⋯m_K+L!T_n^(m_1+m_2+⋯+ m_K+L)
(cos(θ))∏_k∈K+L(ξ_k(θ)k!)^ m_kBy a similar argument as in <ref>, we have that when $\theta$ in $(0,\pi)$ and $K,L\in\ZZ_+$ are fixed then
lim_n|n^-(K+L)^K+L_θT_n(cos(θ))- T_n^(K+L)(cos(θ)) (-sin(θ))^K+L|   =   0Alternatively, taking into account that $T_n(\cos(\theta))=\cos(n\theta)$, we have that
= ξ_K+L(nθ)

Using some uniformity in the behavior (<ref>) with respect to $\theta$ away from 0 and $\pi$, as in the proof of <ref>, we deduce
that for any fixed $\epsilon\in(0,1)$,
∫_A_ϵ|(ξ_K+L(nθ(x,y))(-sin(θ(x,y)))^K+L)^2-n^-2(K+L)(T_n^(K+L)(cos(θ(x,y))))^2|sin^2K(x)sin^2L(y)  dxdy   =   0
with $\theta(x,y)$ the unique solution in $[0,\pi]$ of
cos(θ(x,y))   =    cos(x)+cos(y)2.In particular
we have
lim_ϵ0_+lim sup_n|I_K,L(ϵ,n)- J_K,L(ϵ,n)|   =   0with
1π^22^2(K+L)∫_A_ϵ  ξ_K+L^2(nθ(x,y))(sin(x)sin(θ(x,y)))^2K(sin(y)sin(θ(x,y)))^2L  dxdy

=   4π^22^2(K+L)∫_A_ϵ∩[0,π]^2  ξ_K+L^2(nθ(x,y))(sin(x)sin(θ(x,y)))^2K(sin(y)sin(θ(x,y)))^2L  dxdywhere symmetries $x\mapsto -x$ and $y\mapsto -y$ were taken into account.

Combining all of the above claims, we get that
lim_ϵ0_+lim_n |γ_n[φ_K,L]-J_K,L(ϵ,n)
|   =   0.

§.§ Integral as a second moment

In the final part of the proof we show that, as $n \to \infty$, the quantity $J_{K,L}(\epsilon,n)$ approximates the second moment of some continuous measure $\gamma$.
More specifically, we prove the following lemma.

It holds that

\[
\lim_{\epsilon\ri 0_+}\limsup_{n\ri\iy} J_{K,L}(\epsilon,n)
= \f1{\pi^2}\int_{[0,\pi]^2}  \lt(\f{\sin(x)}{\sin(\theta(x,y))}\rt)^{2K}\lt(\f{\sin(y)}{\sin(\theta(x,y))}\rt)^{2L}\, dxdy.
\]

It follows that we can define the measure $\gamma$ as the image of the measure $\f1{\pi^2} dxdy$ on $[0,\pi]^2$ by the mapping $\Psi$ defined by
Ψ : [0,π]^2∋(x,y)   ↦    (sin(x)sin(θ(x,y)),sin(y)sin(θ(x,y))).In particular it appears that $\gamma$ is a probability measure, so that $\mu=\gamma$, as announced previously.Combining the previous arguments, we see that $\lim_{n\to\infty} \gamma_n[\phi_{K,L}] = \gamma[\phi_{K,L}]$, as wanted.
Proof of Lemma <ref>Let us rewrite the integral $J_{K,L}(\epsilon,n)$ using the change of variables

\[
\st [0,\pi]^2\ni (x,y)\,\mapsto\, (x,\theta(x,y)).
\]

This is a bijection on its image, which is the set
Ω      {(x,θ)∈[0,π]^22cos(θ)-cos(x)∈[-1,1]}.
From the change of variables formula (see Lemma <ref> at the end of this proof), we have that for any measurable $f\st \Omega\ri \RR_+$,
chva∫_[0,π]^2 f(F(x,y))   dxdy    =   ∫_Ω f(x,θ)2 sin(θ)√(1-(2cos(θ)-cos(x))^2)  dxdθ.We can write the integral $J_{K,L}(\epsilon,n)$ in this form by noting that

=    ∫_|cos(θ(x,y))|≤1-ϵ ξ_K+L^2(nθ(x,y))(sin(x)sin(θ(x,y)))^2K(1-cos^2(y))^Lsin^2L(θ(x,y))  dxdy

=   ∫_|cos(θ(x,y))|≤1-ϵ ξ_K+L^2(nθ(x,y))(sin(x)sin(θ(x,y)))^2K(1-(2cos(θ(x,y))-cos(x))^2)^Lsin^2L(θ(x,y))  dxdy

=   ∫_[0,π]^2f_n(F(x,y))  dxdy
(x,θ)∈Ω,     f_n(x,θ)      _[0,1-ϵ](|cos(θ)|)ξ_K+L^2(nθ)(sin(x)sin(θ))^2K(1-(2cos(θ)-cos(x))^2)^Lsin^2L(θ)The change of variables (<ref>)  implies
π^22^2(K+L)J_K,L(ϵ,n)   =   ∫_Ω f_n(x,θ)2 sin(θ)√(1-(2cos(θ)-cos(x))^2)  dxdθNote that
r∈_+, z∈[0,π],    ξ^2_r(z)
1+cos(2z)2   , if $r\equiv 0[2]$

1-cos(2z)2   , if $r\equiv 1[2]$
.Consider the function $g$ defined on $\Omega$ via
(x,θ)∈Ω,     g(x,θ)
(sin(x)sin(θ))^2K(1-(2cos(θ)-cos(x))^2)^Lsin^2L(θ)2 sin(θ)√(1-(2cos(θ)-cos(x))^2)so that
2π^22^2(K+L)J_K,L(ϵ,n)   =   ∫_Ω (1+(-1)^K+Lcos(2nθ))g(x,θ)_[0,1-ϵ](cos(θ))  dxdθFor $\epsilon\in(0,1)$ still fixed, the mapping $\Omega\ni(x,\theta)\mapsto g(x,\theta)\un_{[0,1-\epsilon]}(\cos(\theta))$ is integrable with respect to the Lebesgue measure on $\Omega$. We can apply the Riemann-Lebesgue theorem to get
2π^22^2(K+L)lim_nJ_K,L(ϵ,n)   =   ∫_Ωg(x,θ)_[0,1-ϵ](cos(θ))  dxdθ,and by monotone convergence we have
∫_Ωg(x,θ)_[0,1-ϵ](cos(θ))  dxdθ   =   ∫_Ωg(x,θ)  dxdθ.Finally, note that the previous change of variable computations can be reversed, and yields
12π^22^2(K+L)∫_Ωg(x,θ)  dxdθ   =   1π^2∫_[0,π]^2  (sin(x)sin(θ(x,y)))^2K(sin(y)sin(θ(x,y)))^2L  dxdy.This proves  <ref>.
The following result gives the details for the change of variables formula used in Formula (<ref>).

We have for any measurable $f\st \Omega\ri \RR_+$,
∫_[0,π]^2 f(F(x,y))   dxdy    =   ∫_Ω f(x,θ)2 sin(θ)√(1-(2cos(θ)-cos(x))^2)  dxdθ.
In addition to the set $\Omega$ defined in (<ref>), for given $\theta\in[0,\pi]$, denote
Ω_θ      {x∈[0,π](x,θ)∈Ω}Thus for any measurable and non-negative mapping $g\st \Omega\ri\RR_+$, we have by the change of variables formula
∫_[0,π]^2 g(F(x,y)) |(( F))|(x,y)  dxdy    =   ∫_Ω g(x,θ)  dxdθwhere $\Jac( F)$ is the Jacobian matrix associated to $F$.
We have
( F(x,y))   =   (
1   _xθ(x,y)

0   _yθ(x,y)
Differentiating (<ref>) with respect to $y$, we get, when $\theta(x,y)\not\in\{0,\pi\}$ (i.e. $(x,y)\not\in \{(0,0), (\pi,\pi)\}$),
_yθ(x,y)   =   sin(y)2sin(θ(x,y))so that
|(( F))|(x,y)   =   |_yθ(x,y)|

=   √(1-cos^2(y))2 sin(θ(x,y))

=   √(1-(2cos(θ(x,y))-cos(x))^2)2sin(θ(x,y)) Formula (<ref>) rewrites
∫_[0,π]^2 g(F(x,y))√(1-(2cos(θ(x,y))-cos(x))^2)2 sin(θ(x,y))   dxdy    =   ∫_Ω g(x,θ)  dxdθor
∫_[0,π]^2 f(F(x,y))   dxdy    =   ∫_Ω f(x,θ)2 sin(θ)√(1-(2cos(θ)-cos(x))^2)  dxdθwhere $g$ was replaced by the non-negative and measurable mapping $f$ given by
(x,θ)∈Ω∖{(0,0), (π,π)},    f(x,θ)       g(x,θ)√(1-(2cos(θ)-cos(x))^2)2  sin(θ)
The fact that $\gamma$ is a probability implies (in fact, is equivalent to) the convergence
lim_n∑_p,q∈0,n a_n,p,q^2   =   γ() = 1
(compare with Corollary <ref>), which also amounts to
lim_n∑_p,q∈1,n a_n,p,q^2   =   1
according to (<ref>).

§ DISCUSSION AND FUTURE WORK

We discussed how quantum walks on graphs can be interpreted as solutions of a discrete wave equation on the corresponding graph.
This naturally gives rise to a description of the quantum walk dynamics using Chebyshev polynomials.
We then used the analysis of Chebyshev polynomials to prove a weak limit on the spreading behavior of quantum walks on lattices.

The main open direction of this work is to use Chebyshev polynomials for studying quantum walk spreading and mixing on more general graphs.
The analysis that we did for the lattice seems to naturally extend to lattices of higher dimensions, or Cayley graphs of Abelian groups more generally.
An interesting example out of this category, while still bearing some similarity with the lattice, is the Cayley graph of the Heisenberg group (as studied e.g. in <cit.>).
Preliminary simulations indicate a faster spreading rate of quantum walks, yet we are not aware of a proof of this.

Exploring the connection with the wave equation on graphs by Friedman and Tillich <cit.> is another open direction.
These authors associate to a Markov  generator $\tr_G$ (so that $I-\tr_G$ is a transition matrix, where $I$ is the identity matrix) on a discrete graph $G$ a diffusion generator $\tr_{\cG}$ on the geometric realization $\cG$ of $G$, where the edges have been replaced by segments of length 1 on which usual one-dimensional Laplacian operators are acting.
The domain of $\tr_{\cG}$ consists of functions satisfying some conditions at the vertices.
Friedman and Tillich then relate the eigenvalues of $\tr_G$ and $\tr_{\cG}$, a feature which enables them to interpret $T_k(I-\tr_G)$ as an operator related to the evolution of the wave equation associated to $\tr_\cG$ (see their Theorem 4.9).
We leave it as an open question whether this interpretation of the action of the Chebyshev polynomials on a transition matrix can be used to deduce our weak convergence result.

§ ACKNOWLEDGEMENTS

Simon Apers is partially supported by French projects EPIQ (ANR-22-PETQ-0007), QUDATA (ANR-18-CE47-0010) and QUOPS (ANR22-CE47-0003-01), and EU project QOPT (QuantERA ERA-NET Cofund 2022-25).


††thanks<EMAIL_ADDRESS>
# Neuronal heterogeneity modulates phase-synchronization between
unidirectionally coupled populations with excitation-inhibition balance

Katiele V. P. Brito  Fernanda S. Matias  Instituto de Física, Universidade
Federal de Alagoas, Maceió, Alagoas 57072-970 Brazil

###### Abstract

Several experiments and models have highlighted the importance of neuronal
heterogeneity in brain dynamics and function. However, how such a cell-to-cell
diversity can affect cortical computation, synchronization, and neuronal
communication is still under debate. Previous studies have focused on the
effect of neuronal heterogeneity in one neuronal population. Here we are
specifically interested in the effect of neuronal variability on the phase
relations between two populations, which can be related to different cortical
communication hypotheses. It has been recently shown that two spiking neuron
populations unidirectionally connected in a sender-receiver configuration can
exhibit anticipated synchronization (AS), which is characterized by a negative
phase-lag. This phenomenon has been reported in electrophysiological data of
non-human primates and human EEG during a visual discrimination cognitive
task. In experiments, the unidirectional coupling could be accessed by Granger
causality and can be accompanied by both positive or negative phase difference
between cortical areas. Here we propose a model of two coupled populations in
which the neuronal heterogeneity can determine the dynamical relation between
the sender and the receiver and can reproduce phase relations reported in
experiments. Depending on the distribution of parameters characterizing the
neuronal firing patterns, the system can exhibit both AS and the usual delayed
synchronization regime (DS, with positive phase) as well as a zero-lag
synchronization regime and phase bistability between AS and DS. Furthermore,
we show that our network can present diversity in their phase relations
maintaining the excitation-inhibition balance.

##  I Introduction

The coherent activity of different cortical areas has been considered related
to plenty of cognitive functions, such as, object recognition, visual-motor
integration, and working memory Varela  _et al._ (2001); Buzsaki (2006). Many
hypotheses about neuronal communication take into account the role of brain
oscillation and phase synchronization, for example: the binding by synchrony
Singer (1999) the communication through coherence Fries (2005); Bastos  _et
al._ (2015) the gating by inhibition Jensen and Mazaheri (2010), and nested
oscillations Bonnefond  _et al._ (2017). Although the mechanisms involved in
large-scale integration are still unknown, they have been extensively studied
with biologically inspired neuronal population models Wang (2010).

Spiking neurons are typically considered the building blocks of the population
dynamics and neuronal diversity is ubiquitous across the nervous system.
However, the functional significance of neuronal heterogeneity is still under
investigation Marsat and Maler (2010); Padmanabhan and Urban (2010); Savard
_et al._ (2011); Angelo  _et al._ (2012); Tripathy  _et al._ (2013). From the
experimental point of view, empirical observations show considerable
variability in the response properties of different neurons Beaulieu (1993);
Markram  _et al._ (2004); Angelo  _et al._ (2012), and its relation with
efficient neural coding Marsat and Maler (2010). From the modeling side, a
variety of computational network models of one neuronal population have been
employed to study the effect of heterogeneity in synchronization and coding
capabilities Golomb and Rinzel (1993); Shamir and Sompolinsky (2006); Pérez
_et al._ (2010); Mejias and Longtin (2012, 2014); Rossi  _et al._ (2020), in
self-sustained activity Borges  _et al._ (2020), and in persistent activity,
which could underlie the cognitive function of working memory Renart  _et al._
(2003); Hass  _et al._ (2019).

The previous studies mentioned above have investigated the effect of neuronal
heterogeneity in one neuronal population. Here we are interested in the role
of heterogeneity in the phase synchronization of a physiologically plausible
model of two neuronal populations unidirectionally coupled. This could be
potentially useful to understand the influence of neuronal variability on the
communication between distant brain areas Singer (1999); Fries (2005); Bastos
_et al._ (2015); Jensen and Mazaheri (2010); Bonnefond  _et al._ (2017). The
neuronal population model studied here have been employed before, in the light
of anticipated synchronization ideas Voss (2000); Matias  _et al._ (2014), to
explain electrophysiological results in non-human primate showing that
unidirectional Granger causality can be accompanied by both positive or
negative phase difference between cortical areas Matias  _et al._ (2014);
Montani  _et al._ (2015); Brovelli  _et al._ (2004); Salazar  _et al._ (2012).
However, the effect of neuronal variability on those phase relations have not
been explored.

The typical synchronized regime between two unidirectionally coupled systems
exhibits a positive phase-lag in which the sender is also the leader. This
regime is usually called delayed synchronization (DS) or simply lagged
synchronization. However, it has been shown that the sender-receiver
configuration can also synchronize with a negative phase-lag if the system can
be described by the following equations in which the receiver is subjected to
a delayed self-feedback Voss (2000):

$\displaystyle\dot{\bf{S}}$ $\displaystyle=$ $\displaystyle{\bf f}({\bf
S}(t)),$ (1) $\displaystyle\dot{\bf{R}}$ $\displaystyle=$ $\displaystyle{\bf
f}({\bf R}(t))+{\bf K}[{\bf S}(t)-{\bf R}(t-t_{d})].$

The stable solution ${\bf R}(t)={\bf S}(t-t_{d})$, characterizes the counter-
intuitive situation in which the receiver leads the sender and it is called
anticipated synchronization (AS). This means that the activity of the receiver
predicts the activity of the sender by an amount of time $t_{d}$. In the last
20 years, this solution have been extensively studied in physical systems both
theoretically Voss (2000, 2001a, 2001b); Ciszak  _et al._ (2003); Masoller and
Zanette (2001); Hernández-García  _et al._ (2002); Sausedo-Solorio and
Pisarchik (2014) and experimentally Sivaprakasam  _et al._ (2001); Ciszak  _et
al._ (2009); Tang and Liu (2003). In particular, AS in unidirectional circuits
have been employed for control and prediction of undesirable events Ciszak
_et al._ (2009).

It has been shown that AS can also occur if the delayed self-feedback is
replaced by different parameter mismatches at the receiver Kostur  _et al._
(2005); Pyragienè and Pyragas (2013); Simonov  _et al._ (2014), a faster
internal dynamics of the receiver Hayashi  _et al._ (2016); Dima  _et al._
(2018); Pinto  _et al._ (2019); Dalla Porta  _et al._ (2019), as well as by an
inhibitory loop mediated by chemical synapses Matias  _et al._ (2011, 2016,
2017); Mirasso  _et al._ (2017). Moreover, when the feedback is not hard-wired
in the equation but emerges from system dynamics, two neuronal population can
present phase diversity and a transition from AS to DS through zero-lag
synchronization, induced by synaptic properties Matias  _et al._ (2014); Dalla
Porta  _et al._ (2019). Different patterns of phase synchronization can also
emerge due to the time-delays in heterogeneous networks Petkoski  _et al._
(2018). More recently, it has been shown that a sender-receiver network can
also present phase bistability between DS and AS regimes Machado and Matias
(2020). It is worth to mentioning that in symmetric bidirectional coupled
systems it is not possible to separately classify AS or DS, but only a lagged
synchronization. In case of asymmetric bidirectional influence we can define
the stronger direction of influence as the sender-receiver direction. For
example, the effect of a bidirectional connection in the AS regime have been
studied in neuronal motifs coupled by chemical synapses Matias  _et al._
(2016). Since chemical synapses are unidirectional by their own nature, in
such circuits each direction of influence is independent of each other.

Here we show that the neuronal heterogeneity can promote diversity of phase
relations between two neuronal networks with excitation-inhibition balance. As
far as we know, this is the first verification that the neuronal spiking
properties of the population can promote AS and phase-bistability. In Sec. II,
we describe the neuronal population model as well as the parameters that we
use to change neuronal heterogeneity. In Sec. III, we report our results,
showing that the motif can exhibit phase-locking regimes: with positive,
negative, and zero phases and a bistable regime which alternates from AS to
DS. We also show that the excitatory and inhibitory conductances remain with a
fixed relationship during oscillations. Concluding remarks and a brief
discussion of the significance of our findings for neuroscience are presented
in Sec. IV.

##  II Network model with neuronal heterogeneity

Figure 1:  (a) Schematic representation of two cortical areas, composed of
hundreds of spiking neurons and chemical synapses, coupled in a sender-
receiver configuration. Both populations have neuronal heterogeneity but for
the Receiver (R) population we vary the distribution of parameters determining
the neuronal firing patterns. The inhibitory feedback is controlled by the
synaptic conductance $g_{I}$ at R, whereas the sender-receiver coupling is
determined by excitatory synaptic conductances with $g_{E}$.

### Each node is a neuron model with a specific firing pattern

Our neuronal motif is composed of two unidirectionally coupled cortical-like
neuronal populations: a sender (S) and a receiver (R), see Fig. 1(b). Each one
is composed of 400 excitatory and 100 inhibitory neurons Matias  _et al._
(2014) which is the typical employed proportion of $80\%$ excitatory neurons
and $20\%$ inhibitory neurons, based on anatomical estimates for neocortex
Brunel (2000). Each neuron is described by the Izhikevich model Izhikevich
(2003):

$\displaystyle\frac{dv}{dt}$ $\displaystyle=$ $\displaystyle
0.04v^{2}+5v+140-u+\sum_{x}I_{x},$ (2) $\displaystyle\frac{du}{dt}$
$\displaystyle=$ $\displaystyle a(bv-u).$ (3)

In Eqs. 2 and 3 $v$ is the membrane potential and $u$ the recovery variable
which accounts for activation of K+ and inactivation of Na+ ionic currents.
$I_{x}$ are the synaptic currents provided by the interaction with other
neurons and external inputs. If $v\geq 30$ mV, $v$ is reset to $c$ and $u$ to
$u+d$.

To study the effects of neuronal variability at the receiver population, we
use different values of (c,d) for the excitatory neurons in R:

$\displaystyle c$ $\displaystyle=$
$\displaystyle-55-X+[(5+X)\sigma^{2}]-[(10-X)\sigma^{2}],$ (4) $\displaystyle
d$ $\displaystyle=$ $\displaystyle 4+Y-[(2+Y)\sigma^{2}]+[(4-Y)\sigma^{2}].$
(5)

We use $Y=2X/5$ to guarantee that neuronal type distribution is along the line
changing from regular spikes (RS, $c=-65$ and $d=8$) to chattering (CH,
$c=-50$ and $d=2$) via intrinsically bursting (IB, $c=-55$ and $d=4$). As we
vary $X$ from $-5$ to $10$ the amount of neuron of each type in the R
population changes as shown in Fig. 2. We vary $X$ throughout the paper to
show how heterogeneity can affect phase relations between the two populations.
Illustrative examples of neuronal firings embedding in the neuronal population
but with different parameters are shown in Fig. 3. Equations were integrated
with the Euler method and a time step of $0.05$ ms.

### Each link is an excitatory or inhibitory chemical synapse

Figure 2:  Probability distribution of finding an excitatory neuron at R
population with model parameter $d$, which determines its spiking firing
pattern (see Eqs. 4 and 5). (a) For $X=-5$ it is very probable to find $d=2$
(and consequently $c=-50$) which ensures that there are more chattering (CH)
neuron in the population than other neuron types. (b) For $X=2$ the
distribution of $d$ has a peak around $d=5$ and there are more intrinsically
bursting neurons (IB, $d=4$ and $c=-55$) than other kinds. (c) For $X=10$ the
majority of neurons are regular spiking (RS, $d=8$ and $c=-65$).

The connections between neurons in each population are assumed to be fast
unidirectional excitatory and inhibitory chemical synapses mediated by AMPA
and GABA${}_{\text{A}}$. The synaptic currents are given by:

$I_{x}=g_{x}r_{x}(v-V_{x}),$ (6)

where $x=E,I$ (excitatory and inhibitory mediated by AMPA and
GABA${}_{\text{A}}$, respectively), $V_{E}=0$ mV, $V_{I}=-65$ mV, $g_{x}$ is
the synaptic conductance maximal strengths and $r_{x}$ is the fraction of
bound synaptic receptors whose dynamics is given by:

$\tau_{x}\frac{dr_{x}}{dt}=-r_{x}+D\sum_{k}\delta(t-t_{k}),\\\ $ (7)

where the summation over $k$ stands for pre-synaptic spikes at times $t_{k}$.
D is taken, without loss of generality, equal to $0.05$. The time decays are
$\tau_{E}=5.26$ ms $\tau_{I}=5.6$ ms. Each neuron is subject to an independent
noisy spike train described by a Poisson distribution with rate $R$. The input
mimics excitatory synapses (with conductances $g_{P}$) from $n$ pre-synaptic
neurons external to the population, each one spiking with a Poisson rate $R/n$
which, together with a constant external current $I_{c}$, determine the main
frequency of mean membrane potential of each population. Unless otherwise
stated, we have employed $R=2400$ Hz and $I_{c}=0$.

Figure 3:  Three exemplar neuronal time series to illustrate neuronal
heterogeneity embedded in the network. Firing pattern during the oscillatory
activity of (a) a chattering neuron (CH, $c=-50.4$), (b) an intrinsically
bursting neuron (IB, $c=-56.9$), and (c) a regular spiking neuron (RS,
$c=-64.9$). The purple line shows the mean activity of all neurons in the same
population (for fixed $X=2$ and $g_{I}=2.0$ nS). Figure 4:  Effect of neuronal
heterogeneity in the oscillatory period of one population. The S and R
populations are uncoupled ($g_{E}=0$ nS), and by changing $X$, we only change
the neuronal parameters $c$ and $d$ of neurons in the R population (see Eqs. 4
and 5). The oscillatory period of R depends on both inhibition and the
neuronal variability parameter $X$.  Figure 5:  Characterizing the three
different dynamical regimes exhibited by the two neuronal populations
depending on the parameter $X$ controlling the neuronal spiking variability.
There are two possible phase-locking regimes: delayed synchronization (DS,
$X=-5$) and anticipated synchronization (AS, $X=10$), as well as a bistable
regime in respect to the phase differences (BI, $X=2$). (a)-(c) Time series
for sender and receiver populations. (d)-(f) Time delay $\tau_{i}$ per cycle.
(g)-(i) Probability distribution of time delays. (g) The unimodal distribution
with a positive average characterizes DS, whereas (i) the negative average
represents AS. (j) The bimodal distribution together with an alternation
between cycles of DS randomly followed by cycles of AS characterizes the phase
bistability. Synaptic conductances are kept fixed: $g_{E}=0.5$ nS and
$g_{I}=2.0$ nS.

Connectivity within each population randomly targets 10% of the neurons, with
excitatory conductances set at $g^{S}_{E}=g^{R}_{E}=0.5$ nS. Inhibitory
conductances are fixed at the sender population $g^{S}_{I}=4.0$ nS and $g_{I}$
at the receiver population is varied throughout the study (see Fig. 1). Each
neuron at the R population receives 20 fast synapses (with conductance
$g_{E}$) from random excitatory neurons of the S population. It is worth
emphasizing that to analyze the excitation-inhibition balance at each neuron,
we use $G_{x}=g_{x}r_{x}$ as the effective synaptic conductance at each time
step, which is a fraction of the maximal possible value of $g_{x}$.

### Characterizing phase relations between sender and receiver populations

Since the mean membrane potential $V_{x}$ ($x=$S, R) of each population (which
we assume as a crude approximation of the measured LFP) is noisy, we average
within a sliding window of width 5-8 ms to obtain a smoothed signal, from
which we can extract the peak times $\\{t^{x}_{i}\\}$ (where $i$ indexes the
peak). The period of a given population in each cycle is thus $T^{x}_{i}\equiv
t^{x}_{i+1}-t^{x}_{i}$. For a sufficiently long time series, we compute the
mean period $T_{x}$ and its variance.

In a similar way we calculate the time delay in each cycle
$\tau_{i}=t^{R}_{i}-t^{S}_{i}$. Then, if $\tau_{i}$ obeys a unimodal
distribution, we calculate $\tau$ as the mean value of $\tau_{i}$ and
$\sigma_{\tau}$ as its variance. If $T_{S}\approx T_{R}$ and $g_{E}$ is
independent of the initial conditions, the populations exhibit oscillatory
synchronization with a phase-locking regime. We indistinguishably use the term
phase difference or time delay since it is always possible to associate both:
$\phi_{i}=2\pi\tau_{i}/T_{S}$. In all those calculations we discard the
transient time.

It is also possible to estimate the time-delay by using the peak of the
delayed cross-correlation function $C(V_{S},V_{R},\Delta t)$ between the the
mean membrane potential of the S and R populations. This function can be
calculated as:

$C(V_{S},V_{R},\Delta
t)=\frac{(\sum{V_{S}^{i}-\overline{V_{S}}})(\sum{V_{R}^{i+\Delta
t}-\overline{V_{R}}})}{\sqrt{\sum{(V_{S}^{i}-\overline{V_{S}})^{2}}}\sqrt{\sum{(V_{R}^{i}-\overline{V_{R}})^{2}}}}\;.$
(8)

The values of $C(V_{S},V_{R},\Delta t)$ also indicate the level of
synchronization between the two populations.

##  III Results

### The effect of neuronal heterogeneity in one population

To study the effect of neuronal heterogeneity in the oscillatory properties of
only one population, we analyze the uncoupled case ($g_{E}=0$). For
comparison, the Sender population, which has a fixed distribution of neuronal
spiking types, oscillates with a mean period of $T_{S}=130$ ms. As we change
the neuronal types distribution of the Receiver population by varying $X$ in
Eq. 4, the oscillatory period of R varies from more than $150$ ms to less than
$120$ ms (see Fig. 4 for two different values of inhibition: $g_{I}=2.0$ nS
and $g_{I}=4.0$ nS ). A decrease in the number of chattering neurons and an
increase in the number of intrinsically bursting neurons facilitate faster
oscillations.

This means that we can control the internal dynamics of the population by
changing the local properties of the neurons. In particular, we can turn the
receiver faster maintaining the excitation-inhibition balance by changing $X$.
It is worth mentioning that the E/I balance has been extensively related to
network dynamics, information processing in the nervous system, and social
dysfunction Brunel and Wang (2003); Landau  _et al._ (2016); Yizhar  _et al._
(2011). Moreover, it was not always possible to keep the E/I balance in
previous studies in which the internal dynamics of the receiver were
determined by the relatioship between inhibitory and excitatory synaptic
currents at the receiver Dalla Porta  _et al._ (2019); Machado and Matias
(2020).

### Local properties at the receiver population modulates global phase
relations between two population

The free-running properties of each population can influence the
synchronization patterns between them, when we turn the sender-receiver
coupling on ($g_{E}>0$) Dalla Porta  _et al._ (2019); Pinto  _et al._ (2019);
Dima  _et al._ (2018). In fact, for $g_{E}=0.5$ nS and $g_{I}=2.0$ nS the
motif can exhibit three different regimes depending on $X$. Illustrative
examples of these dynamics are shown in Fig. 5 for $X=-5$, $X=2$ and $X=10$.
The distribution of neuronal spiking parameters shown in Fig. 2 ensures that
the receiver population exhibits more chattering neurons for $X=-5$, more
intrinsically bursting neurons for $X=2$, and more regular spiking neuron for
$X=10$. For $X=-5$ the system presents the usual phase-locking regime with
positive mean time delay $\tau=13$ ms. This means that a peak of the mean
membrane potential of the sender population is followed by a peak of the
receiver. Therefore, the time delay in each cycle $\tau_{i}$ is positive and
fluctuates around a well defined mean value (see the left column in Fig. 5).
For this example, the delayed cross-correlation function $C(V_{S},V_{R},\Delta
t)$ calculated by Eq. 8 has a local peak of $0.92$ for $\Delta t=15$ ms.

For $X=10$, a peak of the Sender is, for the majority of the cycles, preceded
by a peak of the receiver ($\tau_{i}<0$) and, consequently, the mean time
delay is negative (see the right column in Fig. 5). This characterizes the
anticipated synchronization regime (AS, with mean time delay $\tau=-39$ ms for
this example). This counter intuitive regime explains the observed
unidirectional influence with negative phase difference verified in LFP monkey
data Brovelli  _et al._ (2004); Salazar  _et al._ (2012); Matias  _et al._
(2014) as well as in human EEG Carlos  _et al._ (2020). Moreover, AS could be
possibly related to commonly reported short latency in visual systems Orban
_et al._ (1985); Nowak  _et al._ (1995); Kerzel and Gegenfurtner (2003);
Jancke  _et al._ (2004); Puccini  _et al._ (2007); Martinez  _et al._ (2014),
olfactory circuits Rospars  _et al._ (2014), songbirds brain Dima  _et al._
(2018) and human perception Stepp and Turvey (2010, 2017). For this example,
the delayed cross-correlation function $C(V_{S},V_{R},\Delta t)$ has a local
peak of $0.84$ for $\Delta t=-39$ ms.

Figure 6:  Illustrative examples of the raster plots for all neurons in each
dynamical regime. Black dots are neurons from the sender population and purple
(light gray) dots are neurons at the receiver. All parameters as in Fig. 5 and
we only change the distribution of CH, IB, and RS spiking neurons in the
network: $X=-5$ generates DS, $X=2$ promotes the phase bistability, and for
$X=10$ the motif exhibits AS.

For intermediate values of $X$ the system exhibits a phase-bistability between
these two possible phase-locking regimes: DS and AS (see Fig. 5(b), (e) and
(h) with $X=2$). The activity alternates from a few periods of DS eventually
followed by a few periods of AS. In other words, the time delay $\tau_{i}$ is
positive for a few cycles, with a well-defined mean value and standard
deviation (close to $\tau=4.5$ ms in this case), which is similar to a DS
regime for a certain amount of periods. Then, the system randomly switches to
a different dynamics in which $\tau_{i}$ is negative during a few other cycles
(the mean time delay is close to $\tau=-36$ ms for this example). Therefore,
in this regime, the system cannot be simply characterized by the mean time
delay ($\tau$).

A phase-bistability between AS and DS has been previously reported for very
small values of inhibition Machado and Matias (2020). It has been proposed
that this phase bistability could be a plausible model for differences in
phase synchronization of MEG data during bistable perception Kösem  _et al._
(2014). Kosem et al. Kösem  _et al._ (2014) have shown that when participants
listening to bistable (or ambiguous) speech sequences that could be perceived
as two distinct word sequences repeated over time, their MEG recordings
present phase-differences related to which sequence they are perceiving.
Therefore, we propose that this phase-bistability can also be promoted by
neuronal variability in cortical regions.

Figure 7:  Mean time delay $\tau$ between neuronal population as a function of
the parameter controlling the neuronal heterogeneity $X$. The DS-AS transition
can occur via a bistable regime o through zero-lag synchronization for
sufficient large inhibition ($g_{I}=5$ nS and $g_{I}=6$ nS).

The three different regimes can also be visualized in a raster plot as in Fig.
6. In these plots, each dot represents a spike of the i-th neuron indexed by
the vertical axis at the time indicated by the horizontal axis. It is worth
mentioning a few points. Interneurons (indices 400 to 499) are spiking with
higher frequency than excitatory neurons (indices 0 to 399) in both
populations. The spiking pattern of the excitatory neurons at the receiver
population in the DS regime is more concentrated close to the peak of the mean
activity. This means that the neurons are more synchronized with each other.
On the other hand, the variability along the entire period is larger in the AS
regime. Bursting neurons and regular spiking neurons can more easily fire in
the middle of the period than chattering neurons (see an example in Fig. 3).
Since we are changing the proportion of CH, IB, and RS neurons to obtain the
transition from one regime to another, we could speculate that CH neurons
facilitate the synchronization of the entire population and the usual DS
regime, whereas RS neurons allow the observed larger diversity and the AS
regime.

By incremental changes in the neuronal heterogeneity, the system can undergo a
transition from DS to AS through the bistable regime or via zero lag-
synchronization depending on the amount of inhibition at the receiver. Fig. 7
shows the mean time delay as a function of $X$, the parameter controlling the
neuron type distribution, for different values of inhibitory conductance
$g_{I}$. Each curve corresponds to a horizontal line in Fig. 8 which displays
a two-dimensional projection of the parameter space of our model
($X$,$g_{I}$).

For $g_{I}>4.5$nS the transition from DS to AS is continuous and we can find
virtually any value of mean time delay between the two populations, including
a zero-lag synchronization, in which the peak of activity of both regions
occurs very close on average (see $g_{I}=5$nS and $g_{I}=6$nS in Fig. 7). In
Fig. 8 we can see that the end of the bistable regime between DS and AS, gives
place to a DS-AS transition via zero-lag, and coincides with a change in the
slope of the boundary between DS and AS. This re-entrant behavior allows the
system to have, for example for fixed $X=1$, a first DS-AS transition via
bistability followed by an AS-DS transition via zero-lag as a function of
$g_{I}$ (this would corresponds to a vertical line $X=1$ in Fig. 7).

The zero-lag regime has been extensively studied as a non-intuitive regime
that can overcome the synaptic delays between distant areas Vicente  _et al._
(2008); Gollo  _et al._ (2014). The first experimental results about the total
synchronization of distant neurons originated many theories about neuronal
communication as binding by synchrony Singer (1999) and communication through
coherence Fries (2005). With new experimental results about phase relation
diversity, the latter hypothesis has been adapted to include the effect of
non-zero phase relations Bastos  _et al._ (2015); Maris  _et al._ (2013,
2016). Furthermore, it has been shown that in unidirectionally coupled
cortical populations the DS-AS transition via zero-lag can be mediated by
synaptic conductances Matias  _et al._ (2014); Dalla Porta  _et al._ (2019)
and that spiking-timing dependent plasticity can promote auto-organized zero-
lag synchronization Matias  _et al._ (2015).

Figure 8:  Scanning parameter space to visualize the phase relation as a
function of inhibitory conductance $g_{I}$ and the parameter controlling
neuronal variability $X$. The three regimes: DS, AS and BI can be verified for
a reasonable set of parameters. Horizontal lines for integer vales of $g_{I}$
are shown in Fig. 7.

It is worth to emphasizing that previous experimental studies with brain
signals have shown that unidirectional influence can be accompanied by
positive, zero or negative phase differences in different frequency bands
Brovelli  _et al._ (2004); Salazar  _et al._ (2012); Matias  _et al._ (2014);
Carlos  _et al._ (2020). In our model, the populations oscillate with a period
of $130$ ms which is equivalent to a frequency close to $7.7$ Hz. This means
that the time delay of $13$ ms in the DS regime shown in Fig. 5(a) is
equivalent to a phase delay of $0.2\pi$ rad, whereas the anticipation in Fig.
5(c) has a phase delay of $-0.6\pi$ rad. Therefore our model is able to
reproduce unidirectional influence in the alpha band with both positive, zero
and negative phase relations. It has been recently shown that human EEG can
present unidirectional influence between signals with both positive and
negative phase relation Carlos2020]. All the 11 analyzed volunteers presented
pairs of electrodes synchronized in the alpha band (from $7$ to $13$ Hz) with
unidirectional influence and a diversity of phase relations including
positive, negative and zero phase lags. Two other studies with monkey LFP have
shown that unidirectional influence can be accopanied by the couterintuitive
negative phase difference reported here. In Brovelli et al. Brovelli  _et al._
(2004) they have shown that somatosensory motor cortex can influence motor
cortex during oscillatory activity with main frequency around $24$ Hz and a
negative time difference of $-8.7$ ms. In a different experiment Salazar  _et
al._ (2012), it has been shown the posterior parietal cortex can unidirectinal
influence the prefrontal cortex in a frequency of $17$ Hz and present a time
difference that can range from $-2.45$ ms to $-6.53$ ms.

### Excitation-Inhibition balance

In previous studies Matias  _et al._ (2014); Dalla Porta  _et al._ (2019) with
sender-receiver populations, the different regimes were achieved by changing
the excitation-inhibition relation via modifications in the internal
inhibition, the sender-receiver coupling, or the amount of external noise
simulated as excitatory synapses . In particular, the bistable regime reported
in Machado and Matias (2020) required very small values of inhibitory
conductance. Here we can obtain different phase relations ensuring that the
excitation-inhibition balance is maintained by only changing the neuronal
heterogeneity.

We can see the excitatory-inhibitory balance in Fig. 9 where the two effective
synaptic conductances are plotted against each other for four different
neurons of the receiver population: one chattering, one intrinsically
bursting, one regular spiking, and an inhibitory neuron. These plots yield a
linear relationship, which means that the excitatory and inhibitory synaptic
conductance increase and decrease in time in such a way that the ratio of the
two $G_{E}/G_{I}$ remains almost fixed. This ratio slightly varies from neuron
to neuron. The E/I balance has been related to a variety of dynamical features
such as oscillatory activity Brunel and Wang (2003); Compte  _et al._ (2003),
information processing and social dysfunction Yizhar  _et al._ (2011) as well
as cortical complexity Barbero-Castillo  _et al._ (2020).

Figure 9:  Excitatory-inhibitory balance for four illustrative neurons (three
excitatory: CH, IB, RS and one inhibitory) at the receiver population during
the bistable regime ($X=2$ and $g_{I}=2$ nS). Excitatory and inhibitory
effective synaptic conductances preserve an approximate proportional relation
to each other in the course of the oscillations.

##  IV Concluding remarks

We have shown that neuronal heterogeneity can promote different phase
relations between spiking neuron networks. We have verified that by changing
the proportion of different types of neuronal firing patterns the system can
present phase-locking regimes with positive (DS), negative (AS), and zero
phase differences as well as a bistable regime. The DS-AS transition could
possibly explain commonly reported short latency in visual systems Orban  _et
al._ (1985); Nowak  _et al._ (1995); Kerzel and Gegenfurtner (2003); Jancke
_et al._ (2004); Puccini  _et al._ (2007); Martinez  _et al._ (2014),
olfactory circuits Rospars  _et al._ (2014), songbirds brain Dima  _et al._
(2018) and human perception Stepp and Turvey (2010, 2017), whereas the
bistable regime can be associated with bistable perception during ambiguous
stimulation Kösem  _et al._ (2014); Machado and Matias (2020). Previous
studies on AS Matias  _et al._ (2014); Dalla Porta  _et al._ (2019) and phase-
bistability Machado and Matias (2020) in the neuronal population have not
explored the effects of neuronal properties on the network dynamics. Moreover,
the possibility of change the neuronal variability allows us to remain in the
excitation inhibition balance which has been considered a fundamental property
of the health brain Brunel and Wang (2003); Yizhar  _et al._ (2011); Compte
_et al._ (2003); Barbero-Castillo  _et al._ (2020).

Although examples of neuronal variability exist throughout the brain, their
importance for the information process remains a source of debate. The
functional roles of intrinsic neuronal heterogeneity on the network dynamics
have been more extensively studied in the last years Marsat and Maler (2010);
Padmanabhan and Urban (2010); Savard  _et al._ (2011); Angelo  _et al._
(2012); Tripathy  _et al._ (2013). Theoretically, the effect of heterogeneity
on synchronization properties of only one neuronal network has been studied in
different models Golomb and Rinzel (1993); Pérez  _et al._ (2010); Mejias and
Longtin (2012, 2014); Rossi  _et al._ (2020). In particular, Rossi et al.
Rossi  _et al._ (2020) argued that their results can be important in the light
of communication through coherence ideas Fries (2005); Bastos  _et al._
(2015). Here we give a step further in this direction by showing the specific
effect of heterogeneity in the phase synchronization patterns of two coupled
networks. In such a case, we could investigate the effect of neuronal
heterogeneity in phase relations and, consequently, in communication between
distant cortical areas.

Differently from the first papers about AS Voss (2000); Ciszak  _et al._
(2003), here the anticipation time is not hard-wired in the dynamical
equations, neither associated with the excitation-inhibition relation but
rather emerges from the neuronal dynamics and heterogeneity. This opens new
possibilities to study how neuronal variability modulates the phase relation
diversity between cortical areas which has been the object of more intense
investigation in the last years Maris  _et al._ (2013, 2016). In particular,
including effects from homeostasis and the variability of inhibitory neurons
are natural next steps that we are currently pursuing.

###### Acknowledgements.

The authors thank FAPEAL, UFAL, CNPq (grant 432429/2016-6), and CAPES (grant
88881.120309/2016-01) for financial support.

## References

  * Varela  _et al._ (2001) F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie, Nature Reviews Neuroscience 2, 229 (2001).
  * Buzsaki (2006) G. Buzsaki,  _Rhythms of the Brain_ (Oxford University Press (New York), 2006).
  * Singer (1999) W. Singer, Neuron Review 24, 49 (1999).
  * Fries (2005) P. Fries, Trends in Cognitive Sciences 9, 474 (2005).
  * Bastos  _et al._ (2015) A. M. Bastos, J. Vezoli, and P. Fries, Current Opinion in Neurobiology 31, 173 (2015).
  * Jensen and Mazaheri (2010) O. Jensen and A. Mazaheri, Frontiers in human neuroscience 4, 186 (2010).
  * Bonnefond  _et al._ (2017) M. Bonnefond, S. Kastner, and O. Jensen, eneuro 4 (2017).
  * Wang (2010) X. J. Wang, Physiological Reviews 90, 1195 (2010).
  * Marsat and Maler (2010) G. Marsat and L. Maler, Journal of neurophysiology 104, 2543 (2010).
  * Padmanabhan and Urban (2010) K. Padmanabhan and N. N. Urban, Nature neuroscience 13, 1276 (2010).
  * Savard  _et al._ (2011) M. Savard, R. Krahe, and M. Chacron, Neuroscience 172, 270 (2011).
  * Angelo  _et al._ (2012) K. Angelo, E. A. Rancz, D. Pimentel, C. Hundahl, J. Hannibal, A. Fleischmann, B. Pichler, and T. W. Margrie, Nature 488, 375 (2012).
  * Tripathy  _et al._ (2013) S. J. Tripathy, K. Padmanabhan, R. C. Gerkin, and N. N. Urban, Proceedings of the National Academy of Sciences 110, 8248 (2013).
  * Beaulieu (1993) C. Beaulieu, Brain research 609, 284 (1993).
  * Markram  _et al._ (2004) H. Markram, M. Toledo-Rodriguez, Y. Wang, A. Gupta, G. Silberberg, and C. Wu, Nature reviews neuroscience 5, 793 (2004).
  * Golomb and Rinzel (1993) D. Golomb and J. Rinzel, Physical review E 48, 4810 (1993).
  * Shamir and Sompolinsky (2006) M. Shamir and H. Sompolinsky, Neural computation 18, 1951 (2006).
  * Pérez  _et al._ (2010) T. Pérez, C. R. Mirasso, R. Toral, and J. D. Gunton, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368, 5619 (2010).
  * Mejias and Longtin (2012) J. Mejias and A. Longtin, Physical Review Letters 108, 228102 (2012).
  * Mejias and Longtin (2014) J. F. Mejias and A. Longtin, Frontiers in computational neuroscience 8, 107 (2014).
  * Rossi  _et al._ (2020) K. L. Rossi, R. C. Budzisnki, J. A. P. Silveira, B. R. R. Boaretto, T. L. Prado, S. R. Lopes, and U. Feudel, arXiv preprint arXiv:2003.03289 (2020).
  * Borges  _et al._ (2020) F. Borges, P. Protachevicz, R. Pena, E. Lameu, G. Higa, A. Kihara, F. Matias, C. Antonopoulos, R. de Pasquale, A. Roque, _et al._ , Physica A: Statistical Mechanics and its Applications 537, 122671 (2020).
  * Renart  _et al._ (2003) A. Renart, P. Song, and X.-J. Wang, Neuron 38, 473 (2003).
  * Hass  _et al._ (2019) J. Hass, S. Ardid, J. Sherfey, and N. Kopell, bioRxiv , 645663 (2019).
  * Voss (2000) H. U. Voss, Phys. Rev. E 61, 5115 (2000).
  * Matias  _et al._ (2014) F. S. Matias, L. L. Gollo, P. V. Carelli, S. L. Bressler, M. Copelli, and C. R. Mirasso, NeuroImage 99, 411 (2014).
  * Montani  _et al._ (2015) F. Montani, O. A. Rosso, F. S. Matias, S. L. Bressler, and C. R. Mirasso, Phil. Trans. R. Soc. A 373, 20150110 (2015).
  * Brovelli  _et al._ (2004) A. Brovelli, M. Ding, A. Ledberg, Y. Chen, R. Nakamura, and S. L. Bressler, Proc. Natl. Acad. Sci. USA 101, 9849 (2004).
  * Salazar  _et al._ (2012) R. F. Salazar, N. M. Dotson, S. L. Bressler, and C. M. Gray, Science 338, 1097 (2012).
  * Voss (2001a) H. U. Voss, Phys. Rev. E 64, 039904 (2001a).
  * Voss (2001b) H. U. Voss, Phys. Rev. Lett. 87, 014102 (2001b).
  * Ciszak  _et al._ (2003) M. Ciszak, O. Calvo, C. Masoller, C. R. Mirasso, and R. Toral, Phys. Rev. Lett. 90, 204102 (2003).
  * Masoller and Zanette (2001) C. Masoller and D. H. Zanette, Physica A 300, 359 (2001).
  * Hernández-García  _et al._ (2002) E. Hernández-García, C. Masoller, and C. Mirasso, Phys. Lett. A 295, 39 (2002).
  * Sausedo-Solorio and Pisarchik (2014) J. Sausedo-Solorio and A. Pisarchik, Physics Letters A 378, 2108 (2014).
  * Sivaprakasam  _et al._ (2001) S. Sivaprakasam, E. M. Shahverdiev, P. S. Spencer, and K. A. Shore, Phys. Rev. Lett. 87, 154101 (2001).
  * Ciszak  _et al._ (2009) M. Ciszak, C. R. Mirasso, R. Toral, and O. Calvo, Phys. Rev. E 79, 046203 (2009).
  * Tang and Liu (2003) S. Tang and J. M. Liu, Phys. Rev. Lett. 90, 194101 (2003).
  * Kostur  _et al._ (2005) M. Kostur, P. Hänggi, P. Talkner, and J. L. Mateos, Phys. Rev. E 72, 036210 (2005).
  * Pyragienè and Pyragas (2013) T. Pyragienè and K. Pyragas, Nonlinear Dynamics 74, 297 (2013).
  * Simonov  _et al._ (2014) A. Y. Simonov, S. Y. Gordleeva, A. Pisarchik, and V. Kazantsev, JETP Letters 98, 632 (2014).
  * Hayashi  _et al._ (2016) Y. Hayashi, S. J. Nasuto, and H. Eberle, Physical Review E 93, 052229 (2016).
  * Dima  _et al._ (2018) G. C. Dima, M. Copelli, and G. B. Mindlin, International Journal of Bifurcation and Chaos 28, 1830025 (2018).
  * Pinto  _et al._ (2019) M. A. Pinto, O. A. Rosso, and F. S. Matias, Physical Review E 99, 062411 (2019).
  * Dalla Porta  _et al._ (2019) L. Dalla Porta, F. S. Matias, A. J. dos Santos, A. Alonso, P. V. Carelli, M. Copelli, and C. R. Mirasso, Frontiers in Systems Neuroscience 13 (2019).
  * Matias  _et al._ (2011) F. S. Matias, P. V. Carelli, C. R. Mirasso, and M. Copelli, Phys. Rev. E 84, 021922 (2011).
  * Matias  _et al._ (2016) F. S. Matias, L. L. Gollo, P. V. Carelli, C. R. Mirasso, and M. Copelli, Phys. Rev. E 94, 042411 (2016).
  * Matias  _et al._ (2017) F. S. Matias, P. V. Carelli, C. R. Mirasso, and M. Copelli, Physical Review E 95, 052410 (2017).
  * Mirasso  _et al._ (2017) C. R. Mirasso, P. V. Carelli, T. Pereira, F. S. Matias, and M. Copelli, Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 114305 (2017).
  * Petkoski  _et al._ (2018) S. Petkoski, J. M. Palva, and V. K. Jirsa, PLOS Computational Biology 14, e1006160 (2018).
  * Machado and Matias (2020) J. N. Machado and F. S. Matias, Physical Review E 102, 032412 (2020).
  * Brunel (2000) N. Brunel, Journal of computational neuroscience 8, 183 (2000).
  * Izhikevich (2003) E. Izhikevich, IEEE Transaction on Neural Networks 14, 1569 (2003).
  * Brunel and Wang (2003) N. Brunel and X.-J. Wang, Journal of neurophysiology 90, 415 (2003).
  * Landau  _et al._ (2016) I. D. Landau, R. Egger, V. J. Dercksen, M. Oberlaender, and H. Sompolinsky, Neuron 92, 1106 (2016).
  * Yizhar  _et al._ (2011) O. Yizhar, L. E. Fenno, M. Prigge, F. Schneider, T. J. Davidson, D. J. O’shea, V. S. Sohal, I. Goshen, J. Finkelstein, J. T. Paz, _et al._ , Nature 477, 171 (2011).
  * Carlos  _et al._ (2020) F.-L. P. Carlos, M.-M. Ubirakitan, M. C. A. Rodrigues, M. Aguilar-Domingo, E. Herrera-Gutiérrez, J. Gómez-Amor, M. Copelli, P. V. Carelli, and F. S. Matias, Physical Review E 102, 032216 (2020).
  * Orban  _et al._ (1985) G. A. Orban, K. Hoffmann, and J. Duysens, Journal of Neurophysiology 54, 1026 (1985).
  * Nowak  _et al._ (1995) L. Nowak, M. Munk, P. Girard, and J. Bullier, Visual neuroscience 12, 371 (1995).
  * Kerzel and Gegenfurtner (2003) D. Kerzel and K. R. Gegenfurtner, Current Biology 13, 1975 (2003).
  * Jancke  _et al._ (2004) D. Jancke, W. Erlhagen, G. Schöner, and H. R. Dinse, The Journal of physiology 556, 971 (2004).
  * Puccini  _et al._ (2007) G. D. Puccini, M. V. Sanchez-Vives, and A. Compte, PLoS computational biology 3, e82 (2007).
  * Martinez  _et al._ (2014) L. M. Martinez, M. Molano-Mazo, X. Wang, F. T. Sommer, and J. Hirsch, Neuron 81, 943 (2014).
  * Rospars  _et al._ (2014) J.-P. Rospars, A. Grémiaux, D. Jarriault, A. Chaffiol, C. Monsempes, N. Deisig, S. Anton, P. Lucas, and D. Martinez, PLoS Computational Biology 10, e1003975 (2014).
  * Stepp and Turvey (2010) N. Stepp and M. T. Turvey, Cognitive Systems Research 11, 148 (2010).
  * Stepp and Turvey (2017) N. Stepp and M. T. Turvey, Journal of Experimental Psychology: Human Perception and Performance 43, 914 (2017).
  * Kösem  _et al._ (2014) A. Kösem, A. Gramfort, and V. van Wassenhove, Neuroimage 92, 274 (2014).
  * Vicente  _et al._ (2008) R. Vicente, L. L. Gollo, C. R. Mirasso, I. Fischer, and G. Pipa, Proc. Natl. Acad. Sci. USA 105, 17157 (2008).
  * Gollo  _et al._ (2014) L. L. Gollo, C. Mirasso, O. Sporns, and M. Breakspear, PLoS Computational Biology 10, e1003548 (2014).
  * Maris  _et al._ (2013) E. Maris, T. Womelsdorf, R. Desimone, and P. Fries, Neuroimage 74, 99 (2013).
  * Maris  _et al._ (2016) E. Maris, P. Fries, and F. van Ede, Trends in Neurosciences (2016).
  * Matias  _et al._ (2015) F. S. Matias, P. V. Carelli, C. R. Mirasso, and M. Copelli, PloS one 10, e0140504 (2015).
  * Compte  _et al._ (2003) A. Compte, M. V. Sanchez-Vives, D. A. McCormick, and X.-J. Wang, Journal of neurophysiology 89, 2707 (2003).
  * Barbero-Castillo  _et al._ (2020) A. Barbero-Castillo, P. Mateos-Aparicio, L. Dalla Porta, A. Camassa, L. Perez Mendez, M. V. Sanchez-Vives, and M. V. Sanchez-Vives, “Impact of gabaa and gabab inhibition on cortical dynamics and perturbational complexity during synchronous and asynchronous activity,” (2020).