\section{Introduction}

Synchronization in networks of coupled oscillators is attractive phenomenon to study for specialists in various fields of science \cite{ARE08}. Synchronization underlies many natural phenomena and is the cornerstone of many technical concepts and engineering approaches. Examples of synchronization, include, among others, numerous forms of collective behavior in complex biological and artificial systems, such as flocks of birds, swarming, and rendezvous \cite{HER16,SUM10}; ensembles of oscillators \cite{HON11} and a group of mobile robots \cite{REN08}. Special attention is paid to synchronization in neural network dynamics. Synchronization depends on various network parameters, and a time delay in a signal propagation between the nodes plays a crucial role in this phenomenon.

Time delays are always present in real physical systems, therefore, in order to develop adequate realistic models of dynamic networks, one should to take into account delays in signal propagation in order to properly analyze the design of their dynamics. In neural networks, time delays can induce various rhythmic spatiotemporal patterns \cite{COO09,SON09}, change the stability of existing patterns \cite{ERM09}, and play a crucial role in synchronization behavior \cite{DAH09,SCH09}. Various works are devoted to study synchronization in delay coupled networks, just to mention a few, \cite{STE12,PRO13b,SEL15,PLO18}.

This paper continuous the work started in \cite{PLO21}. Here synchronization in heterogeneous networks of linear systems with nonlinear delayed diffusive coupling is considered. The network of neural mass model (NMM) populations \cite{JAN95} is an example of the networks of this type. The problem of network synchronization is reduced to studying the stability of synchronization error system which can be obtained using a coordinate transformation proposed in \cite{PAN17}. The analogue of the circle criterion for time-delay systems (TDS) proposed in \cite{CHU95} and generalized in \cite{BRY19} can be used to study the system stability.

The rest of the paper is organized as follows. Section~\ref{sec:prel} reminds some important concepts related to synchronization and the circle criterion for TDS. In Sec.~\ref{sec:main} synchronization conditions of linear network with nonlinear delayed diffusive couplings are obtained, and NMM network is considered as an example Section~\ref{sec:sim} provides numerical results on synchronization. Finally, conclusions are given in Sec.~\ref{sec:conclusion}.

{\bf Notation.} Throughout the paper the superscript $^{\rm\scriptscriptstyle T}$ ($^*$) stands for matrix transposition (complex conjugate); $\mathbb{R}^n$ denotes the $n$ dimensional real Euclidean space with vector norm $|\cdot|$; $j=\sqrt{-1}$ is the imaginary unit; the notation $z=\col(x,y)$ means that $z$ is a vector of two components $x,~y$; the notation $D=\diag\{d_1,\dots,d_n\}$ means that $D$ is a $n\times n$ diagonal matrix, where $d_i$ is its $i$th diagonal element; $A\otimes B$ means the Kronecker product of matrices $A$ and $B$; $I_n$ is an identity $n\times n$ matrix, while $0_n$ is a $n\times n$ matrix of zeros.

\section{Preliminaries}\label{sec:prel}
\subsection{Circle Criterion for Time-Delay Systems}\label{sec:circle}
In this section the details about studying the stability of linear system with multiple nonlinearities with time-varying delays are given. The analogue of the circle criterion proposed in \cite{CHU95} and its generalization for multi input -- multi output (MIMO) systems proposed in \cite{BRY19} will be used with the purpose of studying the global stability of nonlinear delay-coupled networks.

Consider the system, which is described by the following equations
\begin{equation}\label{c1}
\begin{aligned}
\dot x&(t)=Ax(t)+B\phi(t,\sigma(t-\tau(t))), \\
\sigma&(t-\tau(t))=C^{\rm\scriptscriptstyle T} x(t-\tau(t)),
\end{aligned}
\end{equation}
where $A$ is a constant $n\times n$ matrix, $B$, $C$ are constant $n \times m$ matrices and $x\in\mathbb{R}^n$ is a state vector, \mbox{$\sigma =\col\{\sigma_1,\dots,\sigma_m\}\in\mathbb{R}^m$} is an input of the system, \mbox{$\tau_i=\col\{\tau_1,\dots,\tau_m\}\in\mathbb{R}^m$}, $\tau_i(t)\in[0,T]$, $i=1,\dots,m$, $\forall t$ is a bounded time-varying delay, $\phi=\col\{\phi_1,\dots,\phi_m\}$ is a vector function of sector-bounded nonlinearities, i.e. the following inequalities are fulfilled:
\begin{equation}\label{c2}
\begin{gathered}
\mu_{11}\le\phi_1(t,\sigma_{1})/\sigma_{1}\le \mu_{21}, \\
\vdots \\
\mu_{1m}\le\phi_m(t,\sigma_{m})/\sigma_{m}\le \mu_{2m},
\end{gathered}
\end{equation}
for $\sigma_i\ne0$, $i=1,\dots,m$.
\begin{thm}
\noindent For the system \eqref{c1} denote the transfer function of its linear part $W(p)=C^{\rm\scriptscriptstyle T}(A-p I_n)^{-1}B$ and the characteristic polynomial \mbox{$\Delta(p)=\det(p I_n-A)$}. Introduce the following diagonal matrices $\mu_1=\diag\{\mu_{11},\dots,\mu_{1m}\}$, $\mu_2=\diag\{\mu_{21},\dots \mu_{2m}\}$.
Let the following assumptions be fulfilled:
\begin{enumerate}
\item nonlinearity $\phi(\sigma)$ in system \eqref{c1} satisfies the inequalities \eqref{c2} for $\sigma_i\ne0$, \mbox{$i=1,\dots,m$};
\item There exists diagonal $m\times m$ matrix $\mu_0$ such that each element $\mu_{0i}$ lies between $\mu_{1i}$ and $\mu_{2i}$, $i=1,\dots,m$, and matrix $A+B\mu_0 C^{\rm\scriptscriptstyle T}$ is Hurwitz;
\item For some diagonal $m \times m $ matrix $\nu$ with positive diagonal elements $\nu_i$ such that \mbox{$1-4\nu_i\mu_{1i}\mu_{2i}>0$}, \mbox{$i=1,\dots,m$} the function
\begin{multline}\notag
\pi(\omega)=W(j\omega)^*(\mu_1\mu_2\nu-T^2\omega^2I_m/4)(I_m-4\nu\mu_1\mu_2)\\
\times W(j\omega)+\Real[W(j\omega)(\mu_1+\mu_2)(I_m-4\nu\mu_1\mu_2)\nu]\\
+\nu(I_m-\nu(\mu_1+\mu_2)^2),
\end{multline}
satisfies the following conditions
\begin{subequations}
\begin{eqnarray}\label{c4a}
\lim\limits_{\omega\to\infty}\pi(\omega)&>&0,\\
|\Delta(j\omega)|^2\pi(\omega)&>&0, \quad \forall\omega\ge 0, \label{c4b}
\end{eqnarray}
\end{subequations}
(for $\omega$ such that $|\Delta(j\omega)|=0$ holds the inequality \eqref{c4b} is understood as a limiting one).
\end{enumerate}
Then there exist positive constants $C_1$, $C_2$, $\epsilon$ depending only on the coefficients of the linear part of the system \eqref{c1} $A$, $B$, $C$ and matrices $\mu_1$, $\mu_2$, $TI_n$, such that for all solutions $x(t)$ of the system \eqref{c1} with a continuous initial function $x_0(t)$ defined for $t\in[-\tau_{\max},0]$, where $\tau_{\max}=\max_{i=1,\dots,m}\tau_i(0)$, the following inequality holds
$$
\|x(t)\|\le(C_1\|x_0(0)\|+C_2\max\limits_{-\tau_{\max}\le t\le 0}\|C^{\rm\scriptscriptstyle T} x_0(t)\|)\e^{-\epsilon t},
$$
for all $t\ge0$.
\end{thm}

\subsection{Synchronization}\label{sec:synch}
Here the mathematical notion of synchronization will be introduced. The networks of linear systems with heterogeneous delayed nonlinearities are considered throughout this paper. Each system of this type can be presented in normal form as the following:
\begin{equation}\label{m1}
\begin{aligned}
\dot y_i(t)&=A^yx_i(t)+\phi_i(\sigma_i(t-\tau_i(t))), \\
\dot z_i(t)&=A^zx_i(t), \quad i=1,\dots,N,\\
\end{aligned}
\end{equation}
where $y_i\in\mathbb{R}^m$, $z_i\in\mathbb{R}^{n-m}$, $x_i=\col\{y_i,z_i\}\in\mathbb{R}^n$ and $\sigma_i\in\mathbb{R}^m$ denote the input, the zero-dynamics, the state and the output of the $i$th system, respectively. $A^y\in\mathbb{R}^{m\times n}$ and $A^z\in\mathbb{R}^{(n-m)\times n}$ are constant matrices. $\phi_i=\col\{\phi_{i1}\dots\phi_{im}\}:\mathbb{R}^m\to\mathbb{R}^m$ is a vector function, while $\tau_i=\col\{\tau_{i1},\dots,\tau_{iN}\}$ is a time-varying delay. As input diffusive coupling is considered, which is described by
\begin{equation}\label{contr}
\sigma_i(t-\tau_i(t))=\sum\limits_{k=1}^N \gamma_{ik}[y_i(t-\tau_{ik}(t))-y_k(t-\tau_{ik}(t))].
\end{equation}
Suppose that the graph of the network under consideration is connected and undirected, therefore $\gamma_{ik}=\gamma_{ki}$ \mbox{$\forall i\ne k$}, $i,k=1,\dots,N$. For this type of coupling both signals are time-delayed. Such type of coupling may be observed, for instance, when the systems are interconnected by a centralized control law.

Synchronization phenomenon is often defined as the asymptotically identical evolution
of the systems. One can easily introduce the notion of the asymptotic coordinate synchronization \cite{FRA07}:
\begin{equation}\label{goal1}
\lim\limits_{t\to\infty}(x_i(t)-x_k(t))=0,\quad i,k=1,\dots,N.
\end{equation}
The fulfilment of \eqref{goal1} means that the asymptotic behavior of all nodes of the network \eqref{m1}, \eqref{contr} is identical and can be described by the function $x_s\in\mathbb{R}^n$. Thus by defining the synchronization errors as $e_i=x_i-x_s$ one can study the problem of network synchronization as the problem of stability of synchronization-error system (SES).

To obtain the equations of SES the approach proposed in \cite{PAN17} can be used. The idea behind this approach is the following: the system state space is decomposed in two orthogonal subspaces, one on which is projected the behavior of the mean-field state and one in which lay the synchronization errors. To obtain the network equations in new coordinates suppose that the graph of considered network is connected and undirected. Following the steps described in \cite{PLO21} up to a time-delay one can present the equations of the network \eqref{m1}, \eqref{contr} in the following form:
\begin{subequations}
\begin{eqnarray}\notag
\dot{\tilde x}_1(t)&=& A\tilde x_1(t) + (\mathbf{1}_N\otimes E_m)\\\label{mf}
&\times&\Phi [(LU_1\otimes E_m^{\rm\scriptscriptstyle T})\tilde x_2(t-\tau(t))],\\\notag
\dot{\tilde x}_2(t)&=& (I_{N-1}\otimes A)\tilde x_2(t) + (U_1^\dagger\otimes E_m)\\
&\times&\Phi [(LU_1\otimes E_m^{\rm\scriptscriptstyle T})\tilde x_2(t-\tau(t))],\label{err}
\end{eqnarray}
\end{subequations}
where
\begin{equation}\label{main}
\begin{aligned}
A=\begin{bmatrix}A^y \\ A^z\end{bmatrix}\in\mathbb{R}^{n\times n},\quad \tau(t)=\begin{bmatrix}\tau_1(t) \\ \vdots \\ \tau_N(t),\end{bmatrix}\in\mathbb{R}^{N^2}\\
\Phi(\sigma)=\begin{bmatrix}\phi_1(\sigma_1) \\ \vdots \\ \phi_N(\sigma_N) \end{bmatrix}: \mathbb{R}^{mN}\to\mathbb{R}^{mN}
\end{aligned}
\end{equation}
are the matrix of the linear part of the individual system, the vector of delays and the nonlinear vector function, respectively;
\begin{equation}\label{lapl}
L=\begin{bmatrix}
\sum_{k=2}^{N}\gamma_{1k} &-\gamma_{12} &\cdots&-\gamma_{1N}\\
-\gamma_{21} &\sum_{k=1,k\ne2}^{N}\gamma_{2k}&\cdots&-\gamma_{2N} \\
\vdots &\vdots &\ddots&\vdots \\
-\gamma_{N1} &-\gamma_{N2} &\cdots&\sum_{k=1}^{N-1}\gamma_{Nk}
\end{bmatrix}\in\mathbb{R}^{N\times N},
\end{equation}
is the Laplace matrix defining the coupling links in the network, where $\gamma_{ik}=\gamma_{ki}$ $\forall i\ne k,~i,k=1,\dots,N$ by the assumption;
\begin{equation}\label{supp}
\begin{aligned}
\mathbf{1}_N=\begin{bmatrix}1 \\ \vdots \\1\end{bmatrix}\in\mathbb{R}^N, \quad E_m=\begin{bmatrix}I_m \\0_{(m-n)\times n}\end{bmatrix}\in\mathbb{R}^{n\times m}, \\
U=\begin{bmatrix}\mathbf{1}_N & U_1\end{bmatrix}\in\mathbb{R}^{N\times N}, \quad U^{-1}=\begin{bmatrix}\mathbf{1}_N^{\rm\scriptscriptstyle T} \\ U^\dagger_1 \end{bmatrix}\in\mathbb{R}^{N\times N} \\
\mbox{such that~}U^{-1}LU=\diag\{\lambda_1,\dots,\lambda_N\} \\
\end{aligned}
\end{equation}
are the auxiliary matrices and vectors;
$$
x=\begin{bmatrix}x_1 \\ \vdots \\x_N\end{bmatrix}\in\mathbb{R}^{nN}, ~{\tilde x}_1=\mathbf{1}_N^{\rm\scriptscriptstyle T}\otimes I_n x, ~{\tilde x}_2= U_1^\dagger\otimes I_n x
$$
 are the state space vector of the whole network, the vector proportional to the mean-field dynamics (MFD), and the vector proportional to the synchronization error.

\subsection{Neural Mass Model}
As an example of linear system with nonlinear delayed couplings a time-delay NMM will be considered. This model is based on the standard NMM proposed in \cite{JAN95} and incorporates a time delay \cite{GEN14}. The presence of a time delay in neuronal signal transmission could cause seizure-like activity in the brain. The NMM simulates the average firing activity of a population of pyramidal neurons that interacts with two populations of intercalary neurons and integrates inhibitory and excitatory signals from them \cite{JAN95,GRI06}. The dynamics of each neuronal population is described by two first-order nonlinear delay differential equations as follows:
$$
\begin{aligned}
\dot y_1(t) &= -2\alpha y_1(t) -\alpha ^2y_2(t) + \alpha\beta\phi[\sigma(t-\tau(t))], \\
\dot y_2(t) &= y_1(t).
\end{aligned}
$$
where $y_2$ refers to the post-synaptic potential, i.e. the deviation of the membrane from the resting potential, while $y_1$ is its derivative, $\sigma$ is a delayed input. The parameters $\alpha$ and $\beta$ are different in the excitatory and inhibitory cases. $\alpha$ is the reciprocal of the synaptic/membrane time constant; $\beta$ is the gain for the post-synaptic response kernel; $\tau(t)$ is a time-varying delay. The function $\phi(\sigma)$ is a nonlinear centered sigmoidal function relating the neuronal states
\begin{equation}\label{sigm}
\phi(\sigma) = \frac{g}{1+\e^{k_0-\sigma}}-\frac{g}{1+\e^{k_0}},
\end{equation}
where $k_0$ represents the ratio of average inhibitory synaptic gain and $g>0$ is the average excitatory synaptic gain.

\section{Main Result}\label{sec:main}
\subsection{Synchronization Conditions. General Case}
This paper considers a network of $N$ linear dynamical systems with heterogeneous delayed nonlinear couplings in normal form \eqref{m1}. The graph of considered network is supposed to be connected and undirected, and the connections between the nodes of the network are diffusive ones \eqref{contr}. As described in the Subsec.~\ref{sec:synch}, such a network can be represented in the form of MFD \eqref{mf} and SES \eqref{err}. To study the synchronization problem of such network one can apply the circle criterion for TDS (see Subsec.~\ref{sec:circle}) to the SES \eqref{err} to obtain the conditions of its stability. These conditions in turn will guarantee synchronization in the original network.

The transfer function of the linear part of the system \eqref{err} can be calculated just like in \cite{PLO21}:
\begin{equation}\label{tf}
W(p)=L\otimes[ E_m^{\rm\scriptscriptstyle T} (A-pI_{n})^{-1}E_m],
\end{equation}
where $L$ is the Laplace matrix \eqref{lapl}, $A$ is the matrix of linear part \eqref{main} and $E_m$ is the supplementary matrix \eqref{supp}.

The nonlinear part of the system \eqref{err} are described by the functions $\phi_i=\col\{\phi_{i1}\dots\phi_{im}\}$, $i=1,\dots, N$. Suppose that they belong to the two-cavity sector between two straight lines, i.e. the following inequalities are fulfilled:
\begin{equation}\label{bounds}
\begin{gathered}
\mu_{1i1}\le\phi_{i1}(\sigma_{i1})/\sigma_{i1}\le \mu_{2i1}, \\
\vdots \\
\mu_{1im}\le\phi_{im}(\sigma_{im})/\sigma_{im}\le \mu_{2im}, \\
i=1,\dots, N.
\end{gathered}
\end{equation}
Introduce diagonal matrices matrices by the following way
$$
\begin{gathered}
\mu_1=\diag\{\mu_{111},\dots,\mu_{1m1},\dots,\mu_{11N},\dots,\mu_{1mN}\}, \\
\mu_2=\diag\{\mu_{211},\dots,\mu_{2m1},\dots,\mu_{21N},\dots,\mu_{2mN}\}.
\end{gathered}
$$
Suppose that all delays $\tau_i=\col\{\tau_{i1},\dots,\tau_{iN}\}$, $i=1,\dots,N$ in system \eqref{err} are bounded functions, i.e. $\tau_i(t)\in[0,T]$, $i=1,\dots,m$, $\forall t$.

Thus the theorem about synchronization of linear networks with nonlinear delayed couplings can be formulated.
\begin{thm}\label{th2}
    If the following conditions are fulfilled
\begin{enumerate}
    \item The network \eqref{m1}, \eqref{contr} nonlinearities lie in the sector, i.e. inequalities \eqref{bounds} hold;
    \item The graph of the network \eqref{m1}, \eqref{contr} is connected and undirected.
    \item There exists matrix $\mu_0\in\mathbb{R}^{mN}\times\mathbb{R}^{mN}$ such that matrices $\mu_0-\mu_1$ and $\mu_2-\mu_0$ have only nonnegative elements, and matrix
$$
    \Psi=I_{N-1}\otimes A+[U_1^\dagger\otimes E_m]\mu_0 [LU_1\otimes E_m^{\rm\scriptscriptstyle T}]
$$
    is Hurwitz;
    \item
    For some diagonal $mN \times mN$ matrix $\nu$ with positive diagonal elements such that matrix $I_{mN}-4\nu\mu_{1}\mu_{2}$ has positive diagonal elements, the function
\begin{multline}\notag
\pi(\omega)=W(j\omega)^*(\mu_1\mu_2\nu-T^2\omega^2I_{mN}/4)\\
\times(I_{mN}-4\nu\mu_1\mu_2)W(j\omega)\\
+\Real[W(j\omega)(\mu_1+\mu_2)(I_{mN}-4\nu\mu_1\mu_2)\nu]\\
+\nu(I_{mN}-\nu(\mu_1+\mu_2)^2),
\end{multline}
satisfies the following conditions
\begin{subequations}
\begin{eqnarray}\label{m4a}
\lim\limits_{\omega\to\infty}\pi(\omega)&>&0,\\
|\Delta(j\omega)|^2\pi(\omega)&>&0, \quad \forall\omega\ge 0, \label{m4b}
\end{eqnarray}
\end{subequations}
(for $\omega$ such that $|\Delta(j\omega)|=0$ holds the inequality \eqref{m4b} is understood as a limiting one), where $W(p)$ is the transfer function \eqref{tf} and
\begin{multline}\notag
\Delta(p)=\det(pI_{n(N-1)}-I_{N-1}\otimes A)\\
=\det(I_{N-1}\otimes(pI_n-A))=(\det(pI_n-A))^{N-1}
\end{multline}
is the characteristic polynomial of the matrix of the linear part \eqref{err}.
\end{enumerate}
Then the systems in the network \eqref{m1}, \eqref{contr} are asymptotically synchronized.
\end{thm}

\subsection{Synchronization Conditions of Neural Mass Model Populations}

This section considers the heterogeneous network of nonlinearly delayed coupled NMM populations
\begin{equation} \label{net}
\begin{aligned}
\dot y_{1i}(t) &= -2\alpha y_{1i}(t) -\alpha ^2y_{2i}(t) \\
&+ \alpha\beta\phi_i\left\{\sum\limits_{k=1}^N \gamma_{ik}[y_{1i}(t-\tau_{ik}(t))-y_{1k}(t-\tau_{ik}(t))]\right\}, \\
\dot y_{2i}(t) &= y_{1i}(t), \quad i=1,\dots,N,
\end{aligned}
\end{equation}
where $x=\col\{y_{11},y_{21},\dots,y_{1N},y_{2N}\}$ is a state vector; $\alpha$, $\beta$ are system parameters; $\gamma_{ik}$ are coupling coefficients; $\tau_{ik}(t)\in[0,T]$, $\forall t$, $i,k=1,\dots,N$ are bounded time-varying delays (all delay functions have the same upper bound $T$). Functions $\phi_i$ are sigmoidal ones, which are descirbed by \eqref{sigm} with parameters $g_i$ and $k_{0i}$, $i=1,\dots,N$.

The network of NMMs can be presented in new coordinates \eqref{mf}, \eqref{err} with matrices
$$
A=\begin{bmatrix}-2\alpha & -\alpha ^2 \\1 & 0\end{bmatrix},\quad E_m=\begin{bmatrix} 1 \\ 0\end{bmatrix}.
$$

To study the network \eqref{net} synchronization one should check the conditions of Theorem~\ref{th2}:

Sigmoidal functions \eqref{sigm} are sector ones, which lie between two straight lines $0$ and $0.5g_i\sigma$ (see explanation in \cite{GOR17}), $i=1,\dots,N$. One can find $g_{\max}=\max_{i=1,\dots,N}g_i$. Suppose that $\alpha\beta>0$, then the matrices from the condition (1) of the Theorem~\ref{th2} can be expressed as $\mu_1=0_{N}$ and $\mu_2=0.25g_{\max}\alpha\beta I_{N}$.

Supposing that the graph of considered network \eqref{net} is connected and undirected guarantees the fulfillment of the condition (2) of the Theorem~\ref{th2}.

The condition (3) of Theorem~\ref{th2} is the same as for the case without delays: this fact follows from the conditions of Circle criterion. This condition was previously checked in \cite{PLO21}: If $\alpha>0$, then the matrix $\Psi$ in the condition (3) of the Theorem~\ref{th2} is Hurwitz.

To calculate the function $\pi(\omega)$ from the condition (4) of the Theorem~\ref{th2} the frequency transfer function $W(j\omega)$ should be found using formula \eqref{tf}:
\begin{multline}\label{tfex}
W(j\omega)=L\otimes\left[\begin{bmatrix} 1 & 0\end{bmatrix} \begin{bmatrix}-2\alpha-j\omega & -\alpha ^2 \\1 & -j\omega\end{bmatrix}^{-1} \begin{bmatrix} 1 \\ 0\end{bmatrix}\right]\\=
\frac{-j\omega }{\alpha^2-\omega^2+2j\alpha\omega }L=\frac{-2\alpha\omega^2-j\omega(\alpha^2-\omega^2) }{(\alpha^2+\omega^2)^2}L.
\end{multline}

Meaning that $\mu_1=0_{N}$ and $\mu_2=0.25g_{\max}\alpha\beta I_{N}$ and choosing matrix $\nu$ as $\nu_0 I_{N}$ one obtains:
\begin{multline}\notag
\pi(\omega)=\frac{-T^2\omega^2}{4}W(j\omega)^*W(j\omega)+\frac{g_{\max}\alpha\beta\nu_0}{4}\Real[W(j\omega)]\\
+\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)I_{N}.
\end{multline}

Since the graph of considered network is undirected the corresponding Laplace matrix $L$ is symmetric. Using this fact and \eqref{tfex} the function $\pi(\omega)$ equals:
\begin{multline}\notag
\pi(\omega)=\frac{-T^2\omega^4}{4(\alpha^2+\omega^2)^2}L^2-\frac{g_{\max}\alpha^2\beta\nu_0\omega^2}{2(\alpha^2+\omega^2)^2}L\\
+\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)I_{N}.
\end{multline}

Now check the matrix inequality \eqref{m4a}:
$$
\lim\limits_{\omega\to\infty}\pi(\omega)=\frac{-T^2}{4}L^2+\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)I_{N}>0
$$
Consider some symmetric matrix $P$: it is positive definite if the corresponding quadratic form $x^{\rm\scriptscriptstyle T} Px$ is positive $\forall x\ne 0$. $\exists$ matrix $U$: $P=UDU^{\rm\scriptscriptstyle T}$, where $D$ is a diagonal matrix. Then for $z=U^{\rm\scriptscriptstyle T} x$ one obtains $z^{\rm\scriptscriptstyle T} Dz>0$. Therefore the obtained inequality is fulfilled if and only if the following inequalities are fulfilled:
\begin{equation}\label{ref}
\frac{-T^2\lambda_i^2}{4}+\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)>0, ~ i=1,\dots,N,
\end{equation}
where $\lambda_i$ are the eigenvalues of Laplace matrix $L$. Then the maximal value of the delay $T$ can be estimated:
\begin{equation}\label{eqs}
T^2<\frac{16\nu_0-\nu_0^2g_{\max}^2\alpha^2\beta^2}{4\lambda_i^2}, \quad i=1,\dots,N.
\end{equation}
There is a quadratic equation depending on $\nu_0$ in the numerator of resulting fraction, which has the maximal value for $\nu_0=8/(g_{\max}^2\alpha^2\beta^2)$. All equations \eqref{eqs} hold if
\begin{equation}\label{result}
    T<\frac{4}{g_{\max}\lambda_{\max}\alpha\beta},
\end{equation}
where $\lambda_{\max}$ is the maximal eigenvalue of Laplace matrix $L$.

From \eqref{m4b} one obtains that inequality
\begin{multline}\notag
\frac{-T^2\omega^4}{4}L^2-\frac{g_{\max}\alpha^2\beta\nu_0\omega^2}{2}L\\
+\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)(\alpha^2+\omega^2)^2I_{N}>0
\end{multline}
should be fulfilled $\forall\omega\ge 0$. As before the set of the following inequalities can be consider instead of the obtained matrix inequality:
\begin{multline}\label{help}
\left[\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)-\frac{T^2\lambda_i^2}{4}\right]\omega^4\\
+\alpha^2\nu_0\left[2\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)-\frac{g_{\max}\beta\lambda_i}{2}\right]\omega^2\\
+\alpha^4\nu_0\left(1-\frac{\nu_0g_{\max}^2\alpha^2\beta^2}{16}\right)>0.
\end{multline}
The coefficient before $\omega^4$ is the same as \eqref{ref}, and it is positive if $\nu_0=8/(g_{\max}^2\alpha^2\beta^2)$ and the inequality \eqref{result} is fulfilled. Zero order term is also positive for chosen value of $\nu_0$. The coefficient before $\omega^2$ is positive if $g_{\max}\beta\lambda_{\max}<2$. Thus the inequality \eqref{help} is fulfilled for chosen parameters.

All conditions of the Theorem~\ref{th2} are fulfilled, therefore the network \eqref{net} is synchronized. The following theorem holds.
\begin{thm}\label{th3}
If the network \eqref{net} systems parameters $\alpha>0$ and $\beta>0$, $g_i>0$, $i=1,\dots,N$, the graph of the network is connected and undirected, the maximum eigenvalue of the Laplace matrix $L$ is less than $2/(\beta g_{\max})$, and the delays in the signal propagation are bounded \eqref{result}, then the network of NMMs is asymptotically synchronized.
\end{thm}
Note that while $T\to0$ and choosing $\nu_0\to0$ such that $T^2/\nu_0\to 0$, we get the similar conditions of network synchronization as in \cite{PLO21}.

\section{Simulation}\label{sec:sim}
This section presents the results of simulation. The network of NMM populations with $N=10$ is considered. The system parameters $\alpha$ and $\beta$ are equal to $1$ and $0.8$, respectively. The parameters defining the shape of sigmoidal function have uniform distribution: $g_i$, $i=1,\dots,N$ are distributed on the interval $[0;1]$, while $k_{0i}$, $i=1,\dots,N$ are distributed on the interval $[-1;1]$. The graph of considered network is weighted, connected and undirected, meaning that its adjacency matrix is a symmetric sparse matrix with density $0.7$, which means that it has approximately $0.7N^2$ nonzero entries. Let the delays $\tau_{ik}(t)\in[0;T]$, $i,k=1,\dots,N$ be time-varying functions
$$
\tau_{ik}(t)=h_{1ik}+h_{2ik}\sin(h_{3ik}t+h_{4ik}),
$$
which are uniformly distributed on the interval $[0;T/2]$ such that $h_{1ik}>h_{2ik}$, $i,k=1,\dots,N$. The maximum value of the delay $T$ will be defined later. The initial functions $y_{1i}(t)$, $y_{2i}(t)$, $t\in[-T;0]$, $i=1,\dots,N$ are constants, which are uniformly distributed on the interval $[-1;1]$.

\begin{figure}
\center{\includegraphics[width=1\linewidth]{fig1.eps}}
\caption{Synchronization of neural mass model (NMM) population network \eqref{net} with $N=10$ nodes. (a) and (b) dynamics of post-synaptic potential derivatives $y_1$ and post-synaptic potentials $y_2$ of all nodes, respectively. System parameters: $N=10$, $\alpha=1$, $\beta=0.8$, $T=2$, $g_i$, $i=1,\dots,N$ are uniformly distributed on the interval $[0;1]$, $k_{0i}$, $i=1,\dots,N$ are uniformly distributed on the interval $[-1;1]$, $\lambda_{\max}=2.5080$. Initial functions $y_{1i}(t)$, $y_{2i}(t)$, $t\in[-T;0]$, $i=1,\dots,N$ are constants, which are uniformly distributed on the interval $[-1;1]$.}
\label{fig1}
\end{figure}

First of all consider the case, when NMM network has parameters satisfying the conditions of the Theorem~\ref{th3}. The maximum eigenvalue $\lambda_{\max}$ of the Laplace matrix $L$ is equal to $2.5080$ is this case, while \mbox{$g_{\max}=0.9724$}. For these parameters of the network the inequality \mbox{$\lambda_{\max}<2/(\beta g_{\max})\approx 2.5709$} is fulfilled. Choosing $T=2$ one can ensure the fulfillment of the inequality \eqref{result}, thereby guarantee the fulfillment of all condition of the Theorem~\ref{th3}. This means that for these parameters the network of NMMs will synchronize. Figure~\ref{fig1} presents the results of simulation. As one can see, for the chosen parameter values, there is synchronization among the state variables of the network, and all system trajectories tend to equilibrium point, which confirms~Theorem~\ref{th3}.

\begin{figure}
\center{\includegraphics[width=1\linewidth]{fig2.eps}}
\caption{Desynchronization of neural mass model (NMM) population network \eqref{net} with $N=10$ nodes. (a) and (b) dynamics of post-synaptic potential derivatives $y_1$ and post-synaptic potentials $y_2$ of all nodes, respectively. System parameters: $T=10$, $\lambda_{\max}=1.7983$. Other parameters and initial functions are the same as in Fig.~\ref{fig1}.}
\label{fig2}
\end{figure}

Now consider the case, when the delays in signal propagation between the nodes are too large, that prevents the network synchronization. In this case $\lambda_{\max}=1.7983$ and $g_{\max}=0.9630$, which means that the inequality \mbox{$\lambda_{\max}<2/(\beta g_{\max})\approx 2.5960$} is also fulfilled. Choosing \mbox{$T=10$} one violates the Theorem~\ref{th3} condition \eqref{result}. Therefore for these parameters of the network, the Theorem~\ref{th3} doesn't guarantee the network synchronization. One can see the results of simulation in Fig~\ref{fig2}: there is no synchronization among the state variables of the network.

\section{Conclusion}\label{sec:conclusion}

In this paper, the problem of heterogeneous network synchronization of linear systems with delayed nonlinear diffusive couplings. The delays are supposed to be bounded time-varying functions. The nonlinear coupling functions can be different but their graphs should lie in two-cavity sector between two straight lines. As in \cite{PLO21} the synchronization problem is reduced to study the stability of the SES. The coordinate transformation approach proposed in \cite{PAN17} is used to present the network in coordinates "MFD - SES". To find the conditions for SES stability the circle criterion for TDS was applied. The theorem about network synchronization of this type was formulated and proven.

As an example, the dynamics of NMM populations connected via delayed nonlinear diffusive coupling was considered. Using the obtained theorem, the simple condition for network synchronization was derived. In the case of the delay absence this condition coincide with the result obtained in \cite{PLO21}. Also, the simulation of NMM network dynamics was performed. In the case, when the Theorem~\ref{th3} conditions are fulfilled, one can observe the synchronization between the network states. In the other case, when the delays are large enough, there is no synchronization between the network nodes.


\section{INTRODUCTION}

In medicine, there is a process called Clinical Diagnosis, which identifies a disease, condition, or injury based on the signs and symptoms a patient is having and the patient's health history and physical exam. Further testing, such as blood tests, imaging tests, and biopsies, may then follow \cite{2}. Our research targets those conditions that require no further tests after a physical examination.

Image analysis involves processing an image into fundamental components to extract meaningful information \cite{3} while Natural language processing (NLP) refers to the branch of computer science—and more specifically, the branch of artificial intelligence or AI—concerned with giving computers the ability to understand text and spoken words in much the same way human beings can \cite{4}. An indicator of the practicality of this is the level of accuracy that we will be able to achieve. Our goal is to get accuracy higher or equal to 50\%. Once we have achieved this, we shall combine it with natural language processing and use the two in a mobile application, which should perform the clinical diagnosis process on patients.

In the next section, we will briefly give an overview of the progress, advancements and gaps in image analysis, natural language processing and clinical diagnosis \ref{sec:litraturesurvey} Literature Survey. We will then proceed to \ref{sec:problem_definition} Problem Definition where we present the problem that we aim to solve then we shall proceed to \ref{sec:methodology_approach} Methodology/Approach where we indicate the steps and methods used to conduct our research. We will then go to \ref{sec:resuts_and_discussion} Results and Discussion where we show our findings as well as our thoughts on these findings. We will then \ref{sec:conclusion} Conclude our research and give our recommendations on what we think \ref{sec:future_scope} Future Scope may cover and finally provide all our \ref{sec:references} References.

\section{Literature Survey}
\label{sec:litraturesurvey}

Doctors use Skin Examination to identify skin disorders. The examination may include examination of the scalp, nails and mucus membranes while other times they may use a hand-held lens or a dermoscopy to see areas of concern. Characters that may give clues to the infection include size, shape, colour and the location of the abnormality. If looking at the skin does not provide enough information, doctors may then use other methods like a biopsy, scrapings, culture, wood light (blacklight), tzanck testing, diascopy and other skin tests \cite{5}.

\begin{figure}
  \centering
  \includegraphics[width=.6\textwidth]{images/dermascope.png}
  \caption{Analysing the Skin with a Dermascope.}
  \label{fig:dermascope}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=.6\textwidth]{images/biopsy.png}
  \caption{Punch Biopsy. To remove a small piece of skin for examination under a microscope, doctors may use a round cutting instrument. This procedure is what is known as a punch biopsy.}
  \label{fig:punch_biopsy}
\end{figure}

Image analysis has gained rapid adoption in modern medicine especially in computer-aided diagnostic approaches to provide possible disease mechanisms \cite{6} We have also seen its usage in images of MRI and CT scans as well as the detection of brain tumours. While the mentioned fields use Convolutional Neural Networks \cite{7}, ECG based predictions rely on recurrent neural networks, which are all advanced techniques in image analysis. According to Abhishek Sengupta and Priyanka Narad we expect to see more successful applications of AI-Based in the next 5-10 \cite{8} years. However, the research listed above did not target skin diagnosis but they did review a lot about image analysis in medicine.
It is also worth mentioning that the successful usage of image analysis requires huge volumes of related image data for training and evaluation which is almost impossible to obtain but thanks to modern emerging techniques like the use of synthetic data \cite{9}. A technique that we hope to use in this research.
Mobile applications have seen wide adoption in helping people manage their health and wellness to promote healthy living and gain access to useful medical information. According to industry estimates in 2017, 325,000-health care applications were available for download \cite{10}.
Some notable mobile applications include HnG a mobile app that displays and delivers a wide range of prescription medicine, medisearch, virtual doctor \cite{11}. IBM Watson for Oncology is an assistant that is being used to help provide individualized treatments for cancer patients in China. According to a study conducted by Fang-wen Zou, Babylon showed 78.2\% concordance with real-world clinical analysis. Although this is not good enough to replace humans, it does sufficed to work as a decision support system \cite{12}. Babylon is yet another tool that provides remote consultations with doctors and health care professionals via text and video messaging. Adam Baker also concluded that Babylon gives more accurate results as compared to human doctors \cite{13}. Medwhat provides medical diagnosis by referencing its large pool of expert diagnosis tips provided by doctors, engineers and researchers \cite{14}.

In a study by Petter Bae Brandtzaeg, we learn that people love Instant feedback, others feel more comfortable reviewing personal details to an AI than to a human. This shows great potential for a new paradigm that can change how Zambians receive medical services \cite{15}. In another study, Tara Qian Suna, and Rony Medagliab IBM Watson concludes that people from different cultures and societies have benefited differently from AI \cite{16}.
As indicated by the above review, we see that there has not been much research targeting skin infections in Zambia. The above It is also important to consider the local regulations provided by ZAMRA and HPCZ.

\section{Problem Definition}
\label{sec:problem_definition}

Zambia lacks adequate medical staff as reported by the World Health Organization (WHO) in 2017.  It has less than half the WHO recommended HRH workforce in all categories of the Health sector. 1:17000 Physicians, 1:13000 Nurses and Midwives; 80\% of which are in the public sector \cite{17}.
\begin{figure}[!ht]
  \centering
  \includegraphics[width=.4\textwidth]{images/ratios.jpg}
  \caption{Patient to physician ratios}
  \label{fig:patient_ratios}
\end{figure}

\newpage
\section{Methodology/Approach}
\label{sec:methodology_approach}
\subsection{Purposive Sampling}

The first step is to collect skin related information from specialist doctors/dermatologists. A process also known as purposive sampling. This data will review what kind of skin infections have distinct features, those that are very similar to each other and those that do not always require additional diagnostic physical/clinical examination. Below is a summary of the data that was collected using a survey:

\begin{table}[!ht]
 \caption{Common Diseases}
  \centering
  \begin{tabular}{lll}
    \toprule
    Name     & Requires a lab/biopsy (yes|no) & Similar to\\
    \midrule
    Mycosis & No & Mycoses\\
    Urticaria & No & Eczema,
Angioedema
\\
    Steven-johnson syndrome & No & \\
    Kaposi’s sarcoma & No & \\
    Acne vulgaris & No & Folliculitis\\
    Eczema & No & Acne, psoriasis\\
    Photodermatitis & No &\\
    Alopecia areata & No &\\
    Vitiligo & No \\
    Dermatophytes & No & Eczemas\\
    Psoriasis & No &\\
    \bottomrule
  \end{tabular}
  \label{tab:common_diseases}
\end{table}

This data will then act as a guide for what skin disease images to gather. A bulk image downloader is a good tool that can be used to achieve this.

\subsection{Machine Learning}

Once we have collected the data, the next step will be to prepare the data for the machine learning process. The process is depicted below:

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{images/ml.jpg}
  \caption{The Machine Learning Process}
  \label{fig:machine_learning}
\end{figure}

\newpage
\section{Results and Discussion}
\label{sec:resuts_and_discussion}

The machine learning approach showed very poor scores below 50\% accuracy as shown below:

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{images/mlresults.png}
  \caption{Machine Learning Evaluation Results}
  \label{fig:ml_evaluation_results}
\end{figure}

We also tried to reduce the number of diseases to train so we can understand the correlations between the different diseases. It appeared, that most images were indistinguishable.

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{images/ml2results.png}
  \caption{Machine Learning Evaluation Results for two classes}
  \label{fig:ml_2_classes}
\end{figure}

Further training did not also seem to show improvement as the running loss kept changing haphazardly

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{images/training_performance.png}
  \caption{Training Performance}
  \label{fig:training_performance}
\end{figure}

A confusion matrix of two diseases that seemed to have similar characteristics was as below

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{images/confusion_matrix.png}
  \caption{A Confusion Matrix for two classes}
  \label{fig:confusion_matrix}
\end{figure}

At this point, it was clear machine learning alone was not going to provide a solution so we tried a supervised learning approach.

\subsection{Object Detection}

In object detection, we are required to perform an additional data pre-processing step where we label the images and objects of interest on each image. A useful tool to do this is labellmg \cite{18}.

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{images/labeled_image.png}
  \caption{A labelled image of Kaposi sarcoma}
  \label{fig:labelled_kaposi}
\end{figure}

Once the labelling is done, the labelled images may be converted to different formats for processing. For this, roboflow provides a wide range of useful features. In our case, we need to convert to tf-records which is a TensorFlow data format that allows for the training with object detection.

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.4\textwidth]{images/roboflowmatrix.png}
  \caption{Average precision using roboflow}
  \label{fig:precision_roboflow}
\end{figure}

Once the images were labelled, we identified a suitable model to perform the training. Along the options were SSD, R-CNN and YOLO. Pulkit Sharma provides a good introduction to the various object detection models \cite{18}.

\subsubsection{R-CNN}

Region-based Convolutional Neural Networks are a family of neural networks that uses selective search to extract boxes from an image (these boxes are called regions). Again, Pulkit Sharma provides a more detailed overview of the working concepts of R-CNN.
R-CCN has a higher accuracy despite a slower training process. YOLO is faster with training but not as accurate as R-CNN.

\subsubsection{Colab}

Colab, or "Collaboratory", allows you to write and execute Python in your browser. It also provides access to a GPU. The training was then done using Colab by uploading a jupyter notebook, which is where we put all the preliminary python code ready for execution. With this, we were able to achieve accuracies higher than 60\%.

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.4\textwidth]{images/resultavitiligo.png}
  \caption{Results for Vitiligo }
  \label{fig:results_vitiligo}
\end{figure}

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.4\textwidth]{images/results_warts.png}
  \caption{Results for Warts Vulgaris }
  \label{fig:results_vulgaris}
\end{figure}

\newpage
\section{Conclusion}
\label{sec:conclusion}

The results show us that skin disease diagnosis is very possible with image analysis. We were able to detect the diseases from the images with very high accuracies: 63\%, 88\% and 57\% which are all reasonable scores. When these scores are backed with historical data, we can more confidently diagnose patients. However, we do not provide much information on the historical processing. This will follow in a future research paper. We also restricted the number of classes/diseases to five only so that we could operate within the limits of colab training time.

\section{Future Scope}
\label{sec:future_scope}

Colab did not allow us to train for longer periods so using a computer with a dedicated GPU would greatly improve these results and allow us to train for longer periods.
With more labelled data, we can achieve high scores. More data may have to be collected from health centres and taking into account skin race may radically improve these scores.
Once these two recommendations are met, more diseases can be trained to allow the model to detect several other diseases.
A mobile application, as well as a web application, will then be constructed to integrate these techniques and make them available to skin patients over the internet.

\bibliographystyle{unsrt}


\subsection{Cylinder eigenstates and CFT spectra}

The partition function of a torus CFT can be written as:
\begin{align}
     Z_{\text{CFT}}
  &= \Tr \Bigl[
       \exp \bigl( -2\pi\Im\tau \cdot H_{\text{CFT}} \bigr)
       \exp \bigl(  2\pi\Re\tau \cdot P \bigr)
     \Bigr] \notag \\
  &= \Tr \Bigl[
       \exp \Bigl( -2\pi\Im\tau \Bigl( L_0 + \bar{L}_0 - \frac{c}{12} \Bigr) \Bigr)
       \exp \Bigl(  2\pi\Re\tau \bigl( L_0 - \bar{L}_0 \bigr) \Bigr)
     \Bigr].
\end{align}
Here, operators $H_{\text{CFT}}$ and $P$ denote the Hamiltonian and the total momentum operators that generate translations along the time and space directions, while $c$ is the central charge. $\tau$ is a complex modular parameter that parametrizes the geometry of the torus.

For a local field $\phi_\alpha$ with conformal dimension $(h_\alpha,\bar{h}_\alpha)$, from Virasogo algebra we know that
\begin{align}
     Z_{\text{CFT}}
  &= \sum_\alpha \exp \Bigl[
       - 2\pi    \Im\tau \Bigl( h_0 + \bar{h}_0 - \frac{c}{12} \Bigr)
       + 2\pi\mathrm{i}\mkern1mu \Re\tau \bigl( h_0 - \bar{h}_0 \bigr)
     \Bigr] \notag \\
  &= \sum_\alpha \exp \Bigl[
       - 2\pi    \Im\tau \Bigl(\Delta_\alpha - \frac{c}{12} \Bigr)
       + 2\pi\mathrm{i}\mkern1mu \Re\tau \cdot s_\alpha
     \Bigr],
\end{align}
where $\Delta_\alpha=h_0+\bar{h}_0$ and $s_\alpha=h_0-\bar{h}_0$ are the scaling dimension and conformal spin of $\phi_\alpha$ respectively.

In the continuum limit, a critical classical lattice model with periodic boundary conditions can be described by the torus CFT, with modular parameter $\tau=\mathrm{i}\mkern1mu m/N$. The partition function is then
\begin{equation}
  Z = \sum_\alpha \exp \Bigl[
      - 2\pi \frac{m}{N} \Bigl(\Delta_\alpha - \frac{c}{12} \Bigr)
      + mNf + \mathcal{O} \Bigl( \frac{m}{N^\gamma} \Bigr)
    \Bigr],
\end{equation}
where $f$ is the free energy per each site at the thermodynamic limit, and $\mathcal{O}(m/N^\gamma)$ is the finite-size corrections. If the partition function can be written as $Z=\Tr M^m$ where $M$ is the transfer matrix, the eigenvalues of $M$ are then
\begin{equation}
  \lambda_\alpha = \exp \Bigl[
    - \frac{2\pi}{N} \Bigl( \Delta_\alpha-\frac{c}{12} \Bigr)
    + Nf + \mathcal{O} \Bigl( \frac{1}{N^\gamma} \Bigr)
  \Bigr].
\end{equation}
The conformal spin can be calculated by introducing the translation operator $\exp(2\pi\mathrm{i}\mkern1mu P/N)$, whose eigenvalues are $\exp(2\pi\mathrm{i}\mkern1mu s_\alpha/N)$. For a translationally invariant model, $\exp(2\pi\mathrm{i}\mkern1mu P/N)$ commutes with the transfer matrix $M$, so they can be diagonalized simultaneously to produce a set of eigenvectors with eigenvalues $ \{ \lambda, P_\alpha \}$. One way of estimate the scaling dimension $\Delta_\alpha$ is to fix the identity state $\ket{\phi_I}$ and state $\ket{\phi_T}$ corresponding to energy momentum tensor to have $\Delta_I=0$ and $\Delta_T=2$ so that
\begin{equation}
\Delta_\alpha = \frac{2}{\log \lambda_T - \log \lambda_I} \left( \log \lambda_\alpha - \log \lambda_I \right) .
\end{equation}
Conformal spin is obtained from the eigenvalues of the lattice translation operator $P$ as
\begin{equation}
s_\alpha = \frac{N}{2\pi} P_\alpha .
\end{equation}

\subsection{Numerical details for the Ising model}

\definecolor{blue1}{RGB}{10,120,180}
\definecolor{blue2}{RGB}{107,175,215}
\newcommand{\color{blue1}\bfseries}{\color{blue1}\bfseries}
\newcommand{\color{blue2}}{\color{blue2}}

\input{isingmel}
\clearpage

\subsection{Numerical details for the dimer model}

\input{dimermel}
\section{Introduction}

Critical lattice models provide a profound insight into the emergence of continuous field theory in a discrete world.
There are families of classical lattice partition functions that reproduce CFT data, most notably the series of minimal models.
These lattice models often are integrable, and their topological/categorical symmetries can be built in in many novel constructions developed recently based on fusion categories~\cite{Frank1, Aasen:2020jwb, frank_haagerup, Tachikawa_Haagerup}. These classical lattice partition functions generically
admit tensor network representations. The framework of tensor networks is a powerful tool for recovering the CFT. Particularly, it is found that the process of blocking and RG flow can be done in a controlled and efficient way with arbitrary accuracy depending on the bond dimension of the tensors kept. Tensor Renormalization Group (TRG)~\cite{Levin_2007} or Tensor Network Renormalization (TNR)~\cite{Evenbly_2015} are examples of these blocking algorithms. It is demonstrated that following a recursive blocking algorithm acting on these tensors, one can efficiently read off at least the low lying primaries and descendants from the ``fixed-point'' tensors. The set of tools thus promises an alternative route to obtaining and perhaps potentially classifying CFTs by constructing their lattice versions with various built-in (categorical symmetries) and solving them numerically to very high accuracy.

One very important property characterising 1+1\,d CFT is the infinite dimensional Virasoro symmetries and more generally the Kac--Moody symmetries. There have been discussions of how Virasoro~\cite{Koo:1993wz, Pasquier:1989kd, Milsted:2017csn, Zou:2017zce} and Kac--Moody operators~\cite{Wang} can be implemented in lattice models.
However, these constructions depend on the knowledge of quantum Hamiltonian retrieved from the transfer matrices of the classical lattice models by taking the highly inhomogenous limit of the parameters. In a generic tensor network representation, particularly one in which we obtained by extracting the fixed-point of the blocking procedure, the individual tensor is given by some numerical values --- there is no analytical parameter that one can expand in to obtain the quantum Hamiltonian. In cases as such, it is unclear how these Virasoro/Kac--Moody generators could be implemented.

In this paper, we develop a novel and efficient technique to construct these operators and insert them into the partition function in an arbitrary tensor network representation.
It has been observed that conformal eigenstates, primaries or descendants, are eigenstates of the partition functions made of tensors tiling a cylinder. It is known that to obtain accurate conformal data involving larger conformal dimensions, it is necessary that the total bond dimension of input/output legs of the cylinder is sufficiently large (i.e.\ by packing a large number of tensors into the cylinder, each with small bond dimensions, or a couple tensors with large bond dimensions).
In this paper, we demonstrate that one can obtain the stress tensor $T$ and conserved current $J$ as eigenstates of a very small cylinder, and use them to construct Virasoro/Kac--Moody descendants that are eigenstates of a very large cylinder, to surprising accuracy. The fact that it works suggests that there are interesting algebraic relations between eigenstates of cylinders of different sizes. As a numerical technique, this would open the door to many applications.

The organization of the paper is as follows. We will begin with a brief introduction of tensor network representations of partition functions.
Then we will dive into the precise methodology that has been tested on the Ising and dimer model.
We will end with some discussions. Some further technical details are relegated to the appendices.

\section{Tensor network and Partition functions}

Many 1+1\,d CFTs are continuous/critical limits of lattice models. These lattice partition functions admit tensor network realizations.
It is well known that the transfer matrices of these lattice models can be expressed as a tensor network, the most well-known being the Ising model, which will be the main focus of the current paper.

A generic tensor network representation of a partition function thus takes the following form, as depicted in Fig.~\ref{fig:tm}.

\section{Virasoro and Kac--Moody from Eigenstates of a small cylinder}
In a CFT, it is well known that there is a one-to-one correspondence between states and operators, i.e.\ consider performing a path-integral over a disk with an operator inserted at the origin. This is expected to produce the wavefunction of
a state on the boundary of the disk which corresponds to the inserted operator. If the inserted operator has a well-defined conformal dimension (i.e.\ either a primary or a descendant),
then the state produced would be an eigenstate of the dilation operator. The dilation operator is nothing but the path-integral on an annulus. An annulus with complex coordinates $z$ can be conformally mapped to a path-integral covering a cylinder with coordinates $\omega = \sigma + i \tau$ via the conformal map $z = \exp(\omega)$.
Therefore operator insertions in a planar path-integral can be achieved by cutting out a disk from the plane and replacing the disk by
eigenstates of the cylinder corresponding to the particular conformal operator.
In the case of a lattice partition function, the same local tensors $A_{ijkl}$ defining the lattice partition $Z$ on a plane can be used to cover a cylinder instead.
The cylinder contains one layer of tensors in the vertical direction, and with $N$ sites in the horizontal direction, see also Fig.~\ref{fig:tm}.
It is a matrix with in-coming and out-going indices corresponding to dangling legs of the local tensors at the two boundaries of the cylinder.
Up to normalizations, the eigenstates of this cylinder can be used to approximate operator insertion on a plane. To distinguish the spin of the operators, one can diagonalize
simultaneously with the translation operator $P$ to obtain the finite-size spectrum $\{ \phi_\alpha \}$ along with the conformal data~\cite{Levin_2007, Evenbly_2015, Hauru_2016, Evenbly_2016, Yang_2017, Frank1}, i.e.\ scaling dimensions $\Delta_\alpha$ and the conformal spins $s_\alpha$.

\begin{figure}[ht]
\centering
\includegraphics[width=0.8\linewidth]{transfermatrix.pdf}
\caption{(Left) Partition function constructed from four-leg tensors $A$ with an operator $\phi$ inserted. (Right) The tensors $\phi$ corresponding to a conformal operator are obtained by solving for eigenstates of a cylinder built from either $A$ tensors.}
\label{fig:tm}
\end{figure}

In analogy to the Virasoro generators in the continuum, the lattice version on an $N$-site cylinder can be written as
\begin{equation}
\label{eq:Vir}
L_n \sim \sum_{j=1}^N e^{\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} T(j), \quad \bar{L}_n \sim \sum_{j=1}^N e^{-\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} \bar{T}(j)
\end{equation}
where $T(j)$/$\bar{T}(j)$ are lattice representations of energy-momentum tensor acting on lattice site $j$.

However, for a general tensor network representation of the partition function, the energy-momentum tensor $T$ is often unknown analytically. To overcome this obstacle, we propose the following ansatz for constructing approximate Virasoro operators.

Consider a generic tensor network made up of local tensor $A_{ijkl}$ of bond dimension $\chi_A=d$. We apply exact diagonalization on a small cylinder of four copies of tensor $A$ and identify the eigenstate $\ket{\phi_T}$ that corresponds to the energy-momentum tensor. To reach satisfactory accuracy of the spectrum with a cylinder of just 4 sites, coarse-graining schemes such as direct blocking, TRG or TNR may be required to obtain this local tensor $A$. This state $\ket{\phi_T}$ of size $d^{4}$ is then reshaped into a four-leg tensor of bond dimension $\chi = d$ so it could fit back into the lattice as the approximate $T$ operator, see Fig.~\ref{fig:Vir}.

\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{construction.pdf}
\caption{Tensor network implementation of Eq.~\eqref{eq:Vir}. Eigenstate $\ket{\phi_T}$ or $\ket{\phi_J}$ that corresponds to energy-momentum tensor or Kac--Moody current is solved using a cylinder of 4 sites and reshaped into the four-leg blue tensor $T$ or $J$. The blue tensor $T$/$J$ is inserted at site $j$ and we sum over the contributions of different site $j$'s weighted by appropriate Fourier modes on the ring.}
\label{fig:Vir}
\end{figure}

Now using this approximate energy-momentum tensor $T$ obtained from cylinder of 4 sites, we can build the Virasoro generators for a larger cylinder of $N$ sites, as illustrated in Fig.~\ref{fig:Vir}. To represent $T(j)$, the energy-momentum tensor inserted at site $j$, we replace a local tensor $A$ at site $j$ in the transfer matrix with the approximate energy-momentum tensor $T$. The lattice Virasoro generator $L_n$ is then constructed as in Eq.~\eqref{eq:Vir}.

When the model possesses additional symmetries, the Virasoro symmetry in the underlying CFT theory can get extended to Kac--Moody symmetry. In those cases, the associated Kac--Moody generators on the lattice can be constructed as
\begin{equation}
\label{eq:KM}
J_n \sim \sum_{j=1}^N e^{\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} J(j), \quad \bar{J}_n \sim \sum_{j=1}^N e^{-\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} \bar{J}(j) .
\end{equation}
where $J$ and $\bar{J}$ are current operators with conformal dimensions $(1,0)$ and $(0,1)$.

Again we solve the eigenstate $\ket{\phi_J}$ that corresponds to the conserved current from a cylinder of 4 copies of $A$ by exact diagonlization. This vector $\ket{\phi_J}$ is subsequently reshaped into a four-leg operator and used to build the Kac--Moody generators $J_n$ on a larger cylinder, as described in Fig.~\ref{fig:Vir}.

\begin{figure}[ht]
\centering
\includegraphics[width=0.7\linewidth]{operator.pdf}
\caption{Lattice Virasoro and Kac--Moody generators acting on eigenstates from the cylinder.}
\label{fig:operator}
\end{figure}

We check the validity of our proposal by applying the so constructed Virasoro/Kac--Moody generators to the low energy eigenstates of the transfer matrix of the Ising and the Dimer models, as to be explained in further details in the next section.

\section{Examples}

\subsection{Ising model \& Ising CFT}

As our first test, we study the Ising model on the square lattice $Z = \sum_{\langle ij \rangle} \exp{(-\beta \sigma_i \sigma_j)}$, where $\sigma \in \{-1, +1\}$ is an Ising spin. This partition function can be encoded into a square tensor network consisting of copies of local tensor $A^{(0)}_{ijkl}$ with components
\begin{equation}
\label{eq:isingtensor}
A^{(0)}_{ijkl} = \exp{ \left[ \beta (\sigma_i\sigma_j+\sigma_j\sigma_k+\sigma_k\sigma_l+\sigma_l\sigma_i) \right] }
\end{equation}
where $\beta=\beta_c=\log{(1+\sqrt{2})}/2$ at criticality.

We start by blocking the tensor $A^{(0)}$ into $2\times2$ square block to form a new tensor $A$ of bond dimension $\chi=4$, which is used as the unit tensor in the following calculation. The idea is that we can use the eigenstates of small cylinder with $N=4$ copies of $A$ to construct the Virasoro generators $L_n$ and $\bar{L}_n$ for eigenstates of a larger cylinder (in this example, $N=8$). For both cylinders, low energy eigenstates together with conformal spins and scaling dimensions are obtained by exact diagonalization, see Fig.~\ref{fig:A4spectrum}.

In the spectrum of the small cylinder, we identify the state with conformal spin $s=2$ and scaling dimension $\Delta \approx 2$ as $\ket{\phi_T}$ that corresponds to the energy-momentum tensor $T$ and reshape it into a four-leg operator with bond dimension $\chi=4$. The approximate operator for the anti-holomorphic partner $\bar{T}$ is obtained similarly. (Technical subtlety involved in separating $\ket{\phi_T}$ and $\ket{\phi_{\bar{T}}}$ as degenerate states on a cylinder of 4 sites will be explained in the appendices.)

Using the approximate $T$ and $\bar{T}$ from cylinder of 4 sites, $L_n$ and $\bar{L}_n$ for the larger $N=8$ cylinder can be built via the construction as described in Fig.~\ref{fig:Vir} and Fig.~\ref{fig:operator}. The behaviors of $L_n$ and $L_{-n}$ are explored by examining the matrix elements $\mel{\phi_\beta}{L_n}{\phi_\alpha}$ in the low energy subspace of the $N=8$ transfer matrix.

\begin{figure}[ht]
\centering
\includegraphics[width=0.95\linewidth, left]{ising_spectrum.pdf}
\caption{Low energy spectra of Ising model in terms of scaling dimensions $\Delta$ and conformal spins $s$. Primary states are labeled with their corresponding primary operators in the CFT. We use the approximate lattice operators $T$ and $\bar{T}$ solved from a cylinder of 4 copies of $A$ (left) to construct lattice Virasoro generators $L_n$ and $\bar{L}_n$ for a cylinder of $N=8$ sites (right). As an example, actions of lattice $L_{-1}$ are illustrated on the right. Actions of other $L_n$ and $\bar{L}_{n}$ and the supporting matrix elements can be found in the appendices.}
\label{fig:A4spectrum}
\end{figure}

In the CFT, Virasoro generators $L_n$ act as ladder operators on eigenstates, $L_n \ket{\Delta_\alpha, s_\alpha} \propto \ket{\Delta_\alpha - n, s_\alpha - n}$, raising $\Delta$ for $n<0$ and lowering $\Delta$ for $n>0$. Our numerical results show that, overall $L_n$ and $\bar{L}_{n}$ successfully map the states to their descendants as predicted by CFT with high fidelity $\frac{\lVert \mel{\phi_\beta}{L_n}{\phi_\alpha} \rVert}{\lVert \ket{\phi_\beta}\rVert \cdot \lVert {L_n}\ket{\phi_\alpha}\rVert} \gtrsim 0.9$ (tested with $n=\pm1, \pm2$ on cylinder of $N=8$ sites and bond dimension $\chi=4$). We do observe few erroneous matrix elements (one or two orders of magnitude less than the correct matrix elements) within the same Virasoro towers due to the finite size effect of lattice model and that $T$ and $\bar{T}$ are approximated from a very small cylinder.

\subsection{Dimer model \& free boson CFT}

The second example is the dimer model~\cite{kasteleyn1961, kasteleyn1963, temperley1961} which can be mapped to a bosonic theory via the height mapping~\cite{Allegra:2014mgc, henley1997, Ioffe:1989kd}. The partition function of dimer model can be represented using a square lattice consisting of a local tensor $B^{(0)}_{ijkl}$ which has non-zero components $B^{(0)}_{1111}=B^{(0)}_{2211}=B^{(0)}_{2121}=B^{(0)}_{1212}=B^{(0)}_{2222}=1$ and $B^{(0)}_{1122}=2$
\footnote{This tensor can be obtained by applying one step of TRG on the standard dimer tensor $\tilde{B}_{ijkl}$ which has non-zero components $\tilde{B}_{1112}= \tilde{B}_{1121}=\tilde{B}_{1211} = \tilde{B}_{2111} = 1$.}.
In the continuum limit, this model corresponds to $c=1$ free boson CFT with vertex operators of scaling dimensions $\Delta_{e,m} = e^2 + \frac{1}{4} m^2$ and conformal spin $s_{e,m} = em$, where $e$ and $m$ are the electromagnetic charges in the Coulomb gas description of the dimer model. In addition to Virasoro symmetry, the continuum limit of dimer model exhibits an extended symmetry, Kac--Moody symmetry.

Since the construction of Virasoro generators has already been tested in the previous case, here we will just focus on the applications of our ansatz in constructing the Kac--Moody generators. Again, we take a $2\times2$ block of tensors $B^{(0)}$ to form a new tensor $B$ with bond dimension $\chi = 4$. From the spectrum of a cylinder of 4 copies of $B$, we identify the lowest state with conformal spin $s=1$ as $\ket{\phi_J}$ and the lowest state with $s=-1$ as $\ket{\phi_{\bar{J}}}$. They are then reshaped into four-leg tensors as approximate lattice current operators $J$ and $\bar{J}$. We construct lattice Kac-Moody generators using approximate $J$ and $\bar{J}$ as in Fig.~\ref{fig:Vir} and check their behaviors as in Fig.~\ref{fig:operator} on low energy eigenstates from cylinder of $N=8$ copies of $B$.
\begin{figure}[ht]
\centering
\includegraphics[width=0.95\linewidth, left]{dimer_spectrum.pdf}
\caption{Low energy spectra of dimer model in terms of scaling dimensions $\Delta$ and conformal spins $s$. Using approximate current operators $J$ and $\bar{J}$ obtained from a 4-site cylinder (left) we construct lattice Kac--Moody generators $J_n$ and $\bar{J}_n$ for a cylinder of $N=8$ sites (right) and examine their behaviors in the low energy subspace. Actions of lattice $J_{-1}$ that match with CFT predictions are illustrated using solid black arrows and the erroneous actions are marked with dashed gray arrows. See the appendices more examples.}
\label{fig:B4spectrum}
\end{figure}

In the CFT the Virasoro generators $L_n$ and Kac--Moody generators $J_n$ satisfy the commutation relation that $[L_m, J_n] = -n J_{m+n}$. Therefore $J_n$ act on eigenstates as ladder operators as well i.e. $J_n \ket{\Delta_\alpha, s_\alpha} \propto \ket{\Delta_\alpha - n, s_\alpha - n}$, raising $\Delta$ for $n<0$ and lowering $\Delta$ for $n>0$. We check the property of our lattice Kac--Moody generator by computing the matrix elements $\mel{\phi_\beta}{J_n}{\phi_\alpha}$. We confirm numerically that the correct Kac--Moody descendants $\ket{\phi_\beta} = \ket{\Delta_\beta, s_\beta}$ can be reached by acting $J_{n}$ and $\bar{J}_n$ on an initial state $\ket{\phi_\alpha} = \ket{\Delta_\alpha, s_\alpha}$. However, apart from those correct connections, errors in the matrix elements are also observed. In the case of dimer model, the errors are noticeably larger than those in the Ising case, as the dimer model suffers greater finite size effect. The errors fall in two classes: (i) the mixing of different Kac--Moody towers (mainly the tower of identity $\ket{s=0, \Delta=0}$ and that of $\ket{s=0, \Delta=1}$), which likely results from the fact that the lowest eigenstate of the cylinder is not the true vacuum but instead polluted by other primary states; (ii) the mixing of lowering and raising actions within the same Kac--Moody towers, suggesting that the holomorphic part of the lattice current and the anti-holomorphic part are not fully separated.
Having said that, one can refer to the appendix for the numerical results of the matrix element to see that the mixing is a consistent pattern of mixing -- meaning that the ground state is a definite mixing of the vacuum and several low lying degenerate excited states, and the mixing between $J$ and $\bar{J}$ is also consistent. Apart from this mixing, there is a 5-6 orders of magnitude separation in the matrix elements with states not supposed to be generated by the algebra from a given primary, suggesting that the method is working and that one can in principle isolate the eigenstates.

To enhance the precision of this method, one could increase both the bond dimension $\chi$ of the local tensor in the small cylinder that is used for obtaining the lattice energy-momentum tensor or current, and lattice size $N$ of the cylinder on which the Virasoro/Kac--Moody generators are constructed. Tensor network techniques such as TNR can also be employed for approaching the fixed-point tensor of the model and therefore achieving a more accurate representation of the local energy-momentum tensor/current. It would also be interesting to combine the idea here with the scaling map approach as introduced in~\cite{Evenbly_2016}.

\section{Conclusion}

In this paper, we proposed a method to insert descendant operators in a tensor network representation of the partition function of a critical lattice model. The method was inspired by ``discrete holomorphicity'' --- observations made in 1980s that many CFT symmetries are actually preserved all the way up to the lattice scale (see~\cite{Cardy_2009} for a review). In this paper, we demonstrate that the stress tensor and the Kac--Moody current constructed from eigenstates of a tiny cylinder can be used to generate descendents of states defined on a big cylinder. In fact, in the examples considered, namely the Ising model and also the dimer model corresponding to a $c=1$ free boson CFT, we only need to use eigenstates of a cylinder with 4 sites to construct descendent states in a bigger cylinder. These conserved quantities seem to be defined at the smallest scale, mirroring the situation when these quantum operators were constructed to act on lattice states on a ring~\cite{Koo:1993wz, Milsted:2017csn, Zou:2017zce, Wang}. In the case of discrete holomorphic operators, their conservation follows from co-product structures of quantum algebra which are symmetries of the integrable lattice models involved. It would be very interesting to explore if a new algebraic structure is involved in the realization of the Virasoro and Kac--Moody symmetries.

Aside of its theoretical appeal, the method proposed here is, to our knowledge, the first method available to date to deal with the Virasoro and Kac--Moody symmetries when the quantum Hamiltonian is not readily accessible, an issue particularly pronounced in numerical computations. This should open the door to many applications.

\section{Acknowledgements}

LYH acknowledges the support of NSFC (Grant No. 11922502, 11875111) and the Shanghai Municipal Science and Technology Major Project
(Shanghai Grant No. 2019SHZDZX01), and Perimeter Institute for hospitality as a part of the Emmy Noether Fellowship programme.
We thank Lin Chen, Jiaqi Lou and Bingxin Lao for useful discussions and comments.

\subsection{Cylinder eigenstates and CFT spectra}

The partition function of a torus CFT can be written as:
\begin{align}
     Z_{\text{CFT}}
  &= \Tr \Bigl[
       \exp \bigl( -2\pi\Im\tau \cdot H_{\text{CFT}} \bigr)
       \exp \bigl(  2\pi\Re\tau \cdot P \bigr)
     \Bigr] \notag \\
  &= \Tr \Bigl[
       \exp \Bigl( -2\pi\Im\tau \Bigl( L_0 + \bar{L}_0 - \frac{c}{12} \Bigr) \Bigr)
       \exp \Bigl(  2\pi\Re\tau \bigl( L_0 - \bar{L}_0 \bigr) \Bigr)
     \Bigr].
\end{align}
Here, operators $H_{\text{CFT}}$ and $P$ denote the Hamiltonian and the total momentum operators that generate translations along the time and space directions, while $c$ is the central charge. $\tau$ is a complex modular parameter that parametrizes the geometry of the torus.

For a local field $\phi_\alpha$ with conformal dimension $(h_\alpha,\bar{h}_\alpha)$, from Virasogo algebra we know that
\begin{align}
     Z_{\text{CFT}}
  &= \sum_\alpha \exp \Bigl[
       - 2\pi    \Im\tau \Bigl( h_0 + \bar{h}_0 - \frac{c}{12} \Bigr)
       + 2\pi\mathrm{i}\mkern1mu \Re\tau \bigl( h_0 - \bar{h}_0 \bigr)
     \Bigr] \notag \\
  &= \sum_\alpha \exp \Bigl[
       - 2\pi    \Im\tau \Bigl(\Delta_\alpha - \frac{c}{12} \Bigr)
       + 2\pi\mathrm{i}\mkern1mu \Re\tau \cdot s_\alpha
     \Bigr],
\end{align}
where $\Delta_\alpha=h_0+\bar{h}_0$ and $s_\alpha=h_0-\bar{h}_0$ are the scaling dimension and conformal spin of $\phi_\alpha$ respectively.

In the continuum limit, a critical classical lattice model with periodic boundary conditions can be described by the torus CFT, with modular parameter $\tau=\mathrm{i}\mkern1mu m/N$. The partition function is then
\begin{equation}
  Z = \sum_\alpha \exp \Bigl[
      - 2\pi \frac{m}{N} \Bigl(\Delta_\alpha - \frac{c}{12} \Bigr)
      + mNf + \mathcal{O} \Bigl( \frac{m}{N^\gamma} \Bigr)
    \Bigr],
\end{equation}
where $f$ is the free energy per each site at the thermodynamic limit, and $\mathcal{O}(m/N^\gamma)$ is the finite-size corrections. If the partition function can be written as $Z=\Tr M^m$ where $M$ is the transfer matrix, the eigenvalues of $M$ are then
\begin{equation}
  \lambda_\alpha = \exp \Bigl[
    - \frac{2\pi}{N} \Bigl( \Delta_\alpha-\frac{c}{12} \Bigr)
    + Nf + \mathcal{O} \Bigl( \frac{1}{N^\gamma} \Bigr)
  \Bigr].
\end{equation}
The conformal spin can be calculated by introducing the translation operator $\exp(2\pi\mathrm{i}\mkern1mu P/N)$, whose eigenvalues are $\exp(2\pi\mathrm{i}\mkern1mu s_\alpha/N)$. For a translationally invariant model, $\exp(2\pi\mathrm{i}\mkern1mu P/N)$ commutes with the transfer matrix $M$, so they can be diagonalized simultaneously to produce a set of eigenvectors with eigenvalues $ \{ \lambda, P_\alpha \}$. One way of estimate the scaling dimension $\Delta_\alpha$ is to fix the identity state $\ket{\phi_I}$ and state $\ket{\phi_T}$ corresponding to energy momentum tensor to have $\Delta_I=0$ and $\Delta_T=2$ so that
\begin{equation}
\Delta_\alpha = \frac{2}{\log \lambda_T - \log \lambda_I} \left( \log \lambda_\alpha - \log \lambda_I \right) .
\end{equation}
Conformal spin is obtained from the eigenvalues of the lattice translation operator $P$ as
\begin{equation}
s_\alpha = \frac{N}{2\pi} P_\alpha .
\end{equation}

\subsection{Numerical details for the Ising model}

\definecolor{blue1}{RGB}{10,120,180}
\definecolor{blue2}{RGB}{107,175,215}
\newcommand{\color{blue1}\bfseries}{\color{blue1}\bfseries}
\newcommand{\color{blue2}}{\color{blue2}}

\input{isingmel}
\clearpage

\subsection{Numerical details for the dimer model}

\input{dimermel}
\section{Introduction}

Critical lattice models provide a profound insight into the emergence of continuous field theory in a discrete world.
There are families of classical lattice partition functions that reproduce CFT data, most notably the series of minimal models.
These lattice models often are integrable, and their topological/categorical symmetries can be built in in many novel constructions developed recently based on fusion categories~\cite{Frank1, Aasen:2020jwb, frank_haagerup, Tachikawa_Haagerup}. These classical lattice partition functions generically
admit tensor network representations. The framework of tensor networks is a powerful tool for recovering the CFT. Particularly, it is found that the process of blocking and RG flow can be done in a controlled and efficient way with arbitrary accuracy depending on the bond dimension of the tensors kept. Tensor Renormalization Group (TRG)~\cite{Levin_2007} or Tensor Network Renormalization (TNR)~\cite{Evenbly_2015} are examples of these blocking algorithms. It is demonstrated that following a recursive blocking algorithm acting on these tensors, one can efficiently read off at least the low lying primaries and descendants from the ``fixed-point'' tensors. The set of tools thus promises an alternative route to obtaining and perhaps potentially classifying CFTs by constructing their lattice versions with various built-in (categorical symmetries) and solving them numerically to very high accuracy.

One very important property characterising 1+1\,d CFT is the infinite dimensional Virasoro symmetries and more generally the Kac--Moody symmetries. There have been discussions of how Virasoro~\cite{Koo:1993wz, Pasquier:1989kd, Milsted:2017csn, Zou:2017zce} and Kac--Moody operators~\cite{Wang} can be implemented in lattice models.
However, these constructions depend on the knowledge of quantum Hamiltonian retrieved from the transfer matrices of the classical lattice models by taking the highly inhomogenous limit of the parameters. In a generic tensor network representation, particularly one in which we obtained by extracting the fixed-point of the blocking procedure, the individual tensor is given by some numerical values --- there is no analytical parameter that one can expand in to obtain the quantum Hamiltonian. In cases as such, it is unclear how these Virasoro/Kac--Moody generators could be implemented.

In this paper, we develop a novel and efficient technique to construct these operators and insert them into the partition function in an arbitrary tensor network representation.
It has been observed that conformal eigenstates, primaries or descendants, are eigenstates of the partition functions made of tensors tiling a cylinder. It is known that to obtain accurate conformal data involving larger conformal dimensions, it is necessary that the total bond dimension of input/output legs of the cylinder is sufficiently large (i.e.\ by packing a large number of tensors into the cylinder, each with small bond dimensions, or a couple tensors with large bond dimensions).
In this paper, we demonstrate that one can obtain the stress tensor $T$ and conserved current $J$ as eigenstates of a very small cylinder, and use them to construct Virasoro/Kac--Moody descendants that are eigenstates of a very large cylinder, to surprising accuracy. The fact that it works suggests that there are interesting algebraic relations between eigenstates of cylinders of different sizes. As a numerical technique, this would open the door to many applications.

The organization of the paper is as follows. We will begin with a brief introduction of tensor network representations of partition functions.
Then we will dive into the precise methodology that has been tested on the Ising and dimer model.
We will end with some discussions. Some further technical details are relegated to the appendices.

\section{Tensor network and Partition functions}

Many 1+1\,d CFTs are continuous/critical limits of lattice models. These lattice partition functions admit tensor network realizations.
It is well known that the transfer matrices of these lattice models can be expressed as a tensor network, the most well-known being the Ising model, which will be the main focus of the current paper.

A generic tensor network representation of a partition function thus takes the following form, as depicted in Fig.~\ref{fig:tm}.

\section{Virasoro and Kac--Moody from Eigenstates of a small cylinder}
In a CFT, it is well known that there is a one-to-one correspondence between states and operators, i.e.\ consider performing a path-integral over a disk with an operator inserted at the origin. This is expected to produce the wavefunction of
a state on the boundary of the disk which corresponds to the inserted operator. If the inserted operator has a well-defined conformal dimension (i.e.\ either a primary or a descendant),
then the state produced would be an eigenstate of the dilation operator. The dilation operator is nothing but the path-integral on an annulus. An annulus with complex coordinates $z$ can be conformally mapped to a path-integral covering a cylinder with coordinates $\omega = \sigma + i \tau$ via the conformal map $z = \exp(\omega)$.
Therefore operator insertions in a planar path-integral can be achieved by cutting out a disk from the plane and replacing the disk by
eigenstates of the cylinder corresponding to the particular conformal operator.
In the case of a lattice partition function, the same local tensors $A_{ijkl}$ defining the lattice partition $Z$ on a plane can be used to cover a cylinder instead.
The cylinder contains one layer of tensors in the vertical direction, and with $N$ sites in the horizontal direction, see also Fig.~\ref{fig:tm}.
It is a matrix with in-coming and out-going indices corresponding to dangling legs of the local tensors at the two boundaries of the cylinder.
Up to normalizations, the eigenstates of this cylinder can be used to approximate operator insertion on a plane. To distinguish the spin of the operators, one can diagonalize
simultaneously with the translation operator $P$ to obtain the finite-size spectrum $\{ \phi_\alpha \}$ along with the conformal data~\cite{Levin_2007, Evenbly_2015, Hauru_2016, Evenbly_2016, Yang_2017, Frank1}, i.e.\ scaling dimensions $\Delta_\alpha$ and the conformal spins $s_\alpha$.

\begin{figure}[ht]
\centering
\includegraphics[width=0.8\linewidth]{transfermatrix.pdf}
\caption{(Left) Partition function constructed from four-leg tensors $A$ with an operator $\phi$ inserted. (Right) The tensors $\phi$ corresponding to a conformal operator are obtained by solving for eigenstates of a cylinder built from either $A$ tensors.}
\label{fig:tm}
\end{figure}

In analogy to the Virasoro generators in the continuum, the lattice version on an $N$-site cylinder can be written as
\begin{equation}
\label{eq:Vir}
L_n \sim \sum_{j=1}^N e^{\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} T(j), \quad \bar{L}_n \sim \sum_{j=1}^N e^{-\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} \bar{T}(j)
\end{equation}
where $T(j)$/$\bar{T}(j)$ are lattice representations of energy-momentum tensor acting on lattice site $j$.

However, for a general tensor network representation of the partition function, the energy-momentum tensor $T$ is often unknown analytically. To overcome this obstacle, we propose the following ansatz for constructing approximate Virasoro operators.

Consider a generic tensor network made up of local tensor $A_{ijkl}$ of bond dimension $\chi_A=d$. We apply exact diagonalization on a small cylinder of four copies of tensor $A$ and identify the eigenstate $\ket{\phi_T}$ that corresponds to the energy-momentum tensor. To reach satisfactory accuracy of the spectrum with a cylinder of just 4 sites, coarse-graining schemes such as direct blocking, TRG or TNR may be required to obtain this local tensor $A$. This state $\ket{\phi_T}$ of size $d^{4}$ is then reshaped into a four-leg tensor of bond dimension $\chi = d$ so it could fit back into the lattice as the approximate $T$ operator, see Fig.~\ref{fig:Vir}.

\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{construction.pdf}
\caption{Tensor network implementation of Eq.~\eqref{eq:Vir}. Eigenstate $\ket{\phi_T}$ or $\ket{\phi_J}$ that corresponds to energy-momentum tensor or Kac--Moody current is solved using a cylinder of 4 sites and reshaped into the four-leg blue tensor $T$ or $J$. The blue tensor $T$/$J$ is inserted at site $j$ and we sum over the contributions of different site $j$'s weighted by appropriate Fourier modes on the ring.}
\label{fig:Vir}
\end{figure}

Now using this approximate energy-momentum tensor $T$ obtained from cylinder of 4 sites, we can build the Virasoro generators for a larger cylinder of $N$ sites, as illustrated in Fig.~\ref{fig:Vir}. To represent $T(j)$, the energy-momentum tensor inserted at site $j$, we replace a local tensor $A$ at site $j$ in the transfer matrix with the approximate energy-momentum tensor $T$. The lattice Virasoro generator $L_n$ is then constructed as in Eq.~\eqref{eq:Vir}.

When the model possesses additional symmetries, the Virasoro symmetry in the underlying CFT theory can get extended to Kac--Moody symmetry. In those cases, the associated Kac--Moody generators on the lattice can be constructed as
\begin{equation}
\label{eq:KM}
J_n \sim \sum_{j=1}^N e^{\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} J(j), \quad \bar{J}_n \sim \sum_{j=1}^N e^{-\mathrm{i}\mkern1mu j n \frac{2\pi}{N}} \bar{J}(j) .
\end{equation}
where $J$ and $\bar{J}$ are current operators with conformal dimensions $(1,0)$ and $(0,1)$.

Again we solve the eigenstate $\ket{\phi_J}$ that corresponds to the conserved current from a cylinder of 4 copies of $A$ by exact diagonlization. This vector $\ket{\phi_J}$ is subsequently reshaped into a four-leg operator and used to build the Kac--Moody generators $J_n$ on a larger cylinder, as described in Fig.~\ref{fig:Vir}.

\begin{figure}[ht]
\centering
\includegraphics[width=0.7\linewidth]{operator.pdf}
\caption{Lattice Virasoro and Kac--Moody generators acting on eigenstates from the cylinder.}
\label{fig:operator}
\end{figure}

We check the validity of our proposal by applying the so constructed Virasoro/Kac--Moody generators to the low energy eigenstates of the transfer matrix of the Ising and the Dimer models, as to be explained in further details in the next section.

\section{Examples}

\subsection{Ising model \& Ising CFT}

As our first test, we study the Ising model on the square lattice $Z = \sum_{\langle ij \rangle} \exp{(-\beta \sigma_i \sigma_j)}$, where $\sigma \in \{-1, +1\}$ is an Ising spin. This partition function can be encoded into a square tensor network consisting of copies of local tensor $A^{(0)}_{ijkl}$ with components
\begin{equation}
\label{eq:isingtensor}
A^{(0)}_{ijkl} = \exp{ \left[ \beta (\sigma_i\sigma_j+\sigma_j\sigma_k+\sigma_k\sigma_l+\sigma_l\sigma_i) \right] }
\end{equation}
where $\beta=\beta_c=\log{(1+\sqrt{2})}/2$ at criticality.

We start by blocking the tensor $A^{(0)}$ into $2\times2$ square block to form a new tensor $A$ of bond dimension $\chi=4$, which is used as the unit tensor in the following calculation. The idea is that we can use the eigenstates of small cylinder with $N=4$ copies of $A$ to construct the Virasoro generators $L_n$ and $\bar{L}_n$ for eigenstates of a larger cylinder (in this example, $N=8$). For both cylinders, low energy eigenstates together with conformal spins and scaling dimensions are obtained by exact diagonalization, see Fig.~\ref{fig:A4spectrum}.

In the spectrum of the small cylinder, we identify the state with conformal spin $s=2$ and scaling dimension $\Delta \approx 2$ as $\ket{\phi_T}$ that corresponds to the energy-momentum tensor $T$ and reshape it into a four-leg operator with bond dimension $\chi=4$. The approximate operator for the anti-holomorphic partner $\bar{T}$ is obtained similarly. (Technical subtlety involved in separating $\ket{\phi_T}$ and $\ket{\phi_{\bar{T}}}$ as degenerate states on a cylinder of 4 sites will be explained in the appendices.)

Using the approximate $T$ and $\bar{T}$ from cylinder of 4 sites, $L_n$ and $\bar{L}_n$ for the larger $N=8$ cylinder can be built via the construction as described in Fig.~\ref{fig:Vir} and Fig.~\ref{fig:operator}. The behaviors of $L_n$ and $L_{-n}$ are explored by examining the matrix elements $\mel{\phi_\beta}{L_n}{\phi_\alpha}$ in the low energy subspace of the $N=8$ transfer matrix.

\begin{figure}[ht]
\centering
\includegraphics[width=0.95\linewidth, left]{ising_spectrum.pdf}
\caption{Low energy spectra of Ising model in terms of scaling dimensions $\Delta$ and conformal spins $s$. Primary states are labeled with their corresponding primary operators in the CFT. We use the approximate lattice operators $T$ and $\bar{T}$ solved from a cylinder of 4 copies of $A$ (left) to construct lattice Virasoro generators $L_n$ and $\bar{L}_n$ for a cylinder of $N=8$ sites (right). As an example, actions of lattice $L_{-1}$ are illustrated on the right. Actions of other $L_n$ and $\bar{L}_{n}$ and the supporting matrix elements can be found in the appendices.}
\label{fig:A4spectrum}
\end{figure}

In the CFT, Virasoro generators $L_n$ act as ladder operators on eigenstates, $L_n \ket{\Delta_\alpha, s_\alpha} \propto \ket{\Delta_\alpha - n, s_\alpha - n}$, raising $\Delta$ for $n<0$ and lowering $\Delta$ for $n>0$. Our numerical results show that, overall $L_n$ and $\bar{L}_{n}$ successfully map the states to their descendants as predicted by CFT with high fidelity $\frac{\lVert \mel{\phi_\beta}{L_n}{\phi_\alpha} \rVert}{\lVert \ket{\phi_\beta}\rVert \cdot \lVert {L_n}\ket{\phi_\alpha}\rVert} \gtrsim 0.9$ (tested with $n=\pm1, \pm2$ on cylinder of $N=8$ sites and bond dimension $\chi=4$). We do observe few erroneous matrix elements (one or two orders of magnitude less than the correct matrix elements) within the same Virasoro towers due to the finite size effect of lattice model and that $T$ and $\bar{T}$ are approximated from a very small cylinder.

\subsection{Dimer model \& free boson CFT}

The second example is the dimer model~\cite{kasteleyn1961, kasteleyn1963, temperley1961} which can be mapped to a bosonic theory via the height mapping~\cite{Allegra:2014mgc, henley1997, Ioffe:1989kd}. The partition function of dimer model can be represented using a square lattice consisting of a local tensor $B^{(0)}_{ijkl}$ which has non-zero components $B^{(0)}_{1111}=B^{(0)}_{2211}=B^{(0)}_{2121}=B^{(0)}_{1212}=B^{(0)}_{2222}=1$ and $B^{(0)}_{1122}=2$
\footnote{This tensor can be obtained by applying one step of TRG on the standard dimer tensor $\tilde{B}_{ijkl}$ which has non-zero components $\tilde{B}_{1112}= \tilde{B}_{1121}=\tilde{B}_{1211} = \tilde{B}_{2111} = 1$.}.
In the continuum limit, this model corresponds to $c=1$ free boson CFT with vertex operators of scaling dimensions $\Delta_{e,m} = e^2 + \frac{1}{4} m^2$ and conformal spin $s_{e,m} = em$, where $e$ and $m$ are the electromagnetic charges in the Coulomb gas description of the dimer model. In addition to Virasoro symmetry, the continuum limit of dimer model exhibits an extended symmetry, Kac--Moody symmetry.

Since the construction of Virasoro generators has already been tested in the previous case, here we will just focus on the applications of our ansatz in constructing the Kac--Moody generators. Again, we take a $2\times2$ block of tensors $B^{(0)}$ to form a new tensor $B$ with bond dimension $\chi = 4$. From the spectrum of a cylinder of 4 copies of $B$, we identify the lowest state with conformal spin $s=1$ as $\ket{\phi_J}$ and the lowest state with $s=-1$ as $\ket{\phi_{\bar{J}}}$. They are then reshaped into four-leg tensors as approximate lattice current operators $J$ and $\bar{J}$. We construct lattice Kac-Moody generators using approximate $J$ and $\bar{J}$ as in Fig.~\ref{fig:Vir} and check their behaviors as in Fig.~\ref{fig:operator} on low energy eigenstates from cylinder of $N=8$ copies of $B$.
\begin{figure}[ht]
\centering
\includegraphics[width=0.95\linewidth, left]{dimer_spectrum.pdf}
\caption{Low energy spectra of dimer model in terms of scaling dimensions $\Delta$ and conformal spins $s$. Using approximate current operators $J$ and $\bar{J}$ obtained from a 4-site cylinder (left) we construct lattice Kac--Moody generators $J_n$ and $\bar{J}_n$ for a cylinder of $N=8$ sites (right) and examine their behaviors in the low energy subspace. Actions of lattice $J_{-1}$ that match with CFT predictions are illustrated using solid black arrows and the erroneous actions are marked with dashed gray arrows. See the appendices more examples.}
\label{fig:B4spectrum}
\end{figure}

In the CFT the Virasoro generators $L_n$ and Kac--Moody generators $J_n$ satisfy the commutation relation that $[L_m, J_n] = -n J_{m+n}$. Therefore $J_n$ act on eigenstates as ladder operators as well i.e. $J_n \ket{\Delta_\alpha, s_\alpha} \propto \ket{\Delta_\alpha - n, s_\alpha - n}$, raising $\Delta$ for $n<0$ and lowering $\Delta$ for $n>0$. We check the property of our lattice Kac--Moody generator by computing the matrix elements $\mel{\phi_\beta}{J_n}{\phi_\alpha}$. We confirm numerically that the correct Kac--Moody descendants $\ket{\phi_\beta} = \ket{\Delta_\beta, s_\beta}$ can be reached by acting $J_{n}$ and $\bar{J}_n$ on an initial state $\ket{\phi_\alpha} = \ket{\Delta_\alpha, s_\alpha}$. However, apart from those correct connections, errors in the matrix elements are also observed. In the case of dimer model, the errors are noticeably larger than those in the Ising case, as the dimer model suffers greater finite size effect. The errors fall in two classes: (i) the mixing of different Kac--Moody towers (mainly the tower of identity $\ket{s=0, \Delta=0}$ and that of $\ket{s=0, \Delta=1}$), which likely results from the fact that the lowest eigenstate of the cylinder is not the true vacuum but instead polluted by other primary states; (ii) the mixing of lowering and raising actions within the same Kac--Moody towers, suggesting that the holomorphic part of the lattice current and the anti-holomorphic part are not fully separated.
Having said that, one can refer to the appendix for the numerical results of the matrix element to see that the mixing is a consistent pattern of mixing -- meaning that the ground state is a definite mixing of the vacuum and several low lying degenerate excited states, and the mixing between $J$ and $\bar{J}$ is also consistent. Apart from this mixing, there is a 5-6 orders of magnitude separation in the matrix elements with states not supposed to be generated by the algebra from a given primary, suggesting that the method is working and that one can in principle isolate the eigenstates.

To enhance the precision of this method, one could increase both the bond dimension $\chi$ of the local tensor in the small cylinder that is used for obtaining the lattice energy-momentum tensor or current, and lattice size $N$ of the cylinder on which the Virasoro/Kac--Moody generators are constructed. Tensor network techniques such as TNR can also be employed for approaching the fixed-point tensor of the model and therefore achieving a more accurate representation of the local energy-momentum tensor/current. It would also be interesting to combine the idea here with the scaling map approach as introduced in~\cite{Evenbly_2016}.

\section{Conclusion}

In this paper, we proposed a method to insert descendant operators in a tensor network representation of the partition function of a critical lattice model. The method was inspired by ``discrete holomorphicity'' --- observations made in 1980s that many CFT symmetries are actually preserved all the way up to the lattice scale (see~\cite{Cardy_2009} for a review). In this paper, we demonstrate that the stress tensor and the Kac--Moody current constructed from eigenstates of a tiny cylinder can be used to generate descendents of states defined on a big cylinder. In fact, in the examples considered, namely the Ising model and also the dimer model corresponding to a $c=1$ free boson CFT, we only need to use eigenstates of a cylinder with 4 sites to construct descendent states in a bigger cylinder. These conserved quantities seem to be defined at the smallest scale, mirroring the situation when these quantum operators were constructed to act on lattice states on a ring~\cite{Koo:1993wz, Milsted:2017csn, Zou:2017zce, Wang}. In the case of discrete holomorphic operators, their conservation follows from co-product structures of quantum algebra which are symmetries of the integrable lattice models involved. It would be very interesting to explore if a new algebraic structure is involved in the realization of the Virasoro and Kac--Moody symmetries.

Aside of its theoretical appeal, the method proposed here is, to our knowledge, the first method available to date to deal with the Virasoro and Kac--Moody symmetries when the quantum Hamiltonian is not readily accessible, an issue particularly pronounced in numerical computations. This should open the door to many applications.

\section{Acknowledgements}

LYH acknowledges the support of NSFC (Grant No. 11922502, 11875111) and the Shanghai Municipal Science and Technology Major Project
(Shanghai Grant No. 2019SHZDZX01), and Perimeter Institute for hospitality as a part of the Emmy Noether Fellowship programme.
We thank Lin Chen, Jiaqi Lou and Bingxin Lao for useful discussions and comments.