\begingroup

%% --- Glazek macros ---
\newcommand{\bal}{\begin{align}}
\newcommand{\eal}{\end{align}}
\newcommand{\beq}{\begin{eqnarray}}
\newcommand{\eeq}{\end{eqnarray}}
\newcommand{\nneeq}{\nonumber \end{eqnarray}}
\newcommand{\nnn}{\nonumber}
\newcommand{\nn}{\nonumber \\}
\newcommand{\es}{& = &}
\newcommand{\eqv}{& \equiv &}
\newcommand{\rs}{\, = \,}
\newcommand{\ps}{& + &}
\newcommand{\ms}{& - &}
\newcommand{\ts}{& \times &}
\newcommand{\sims}{& \sim &}
\newcommand{\nt}{\nn \ts}
\newcommand{\np}{\nn \ps}
\newcommand{\nm}{\nn \ms}
\def\n{\nonumber}
\newcommand{\tl}{ & \to & }

\newcommand{\cA}{ {\cal A} }
\newcommand{\cB}{ {\cal B} }
\newcommand{\cAB}{ {\cA \2space{-3.7pt} \cB} }
\newcommand{\cC}{ {\cal C} }
\newcommand{\cD}{ {\cal D} }
\newcommand{\cM}{ {\cal M} }
\newcommand{\cH}{ {\cal H} }
\newcommand{\cE}{ {\cal E} }
\newcommand{\cF}{ {\cal F} }
\newcommand{\cG}{ {\cal G} }
\newcommand{\cT}{ {\cal T} }
\newcommand{\cV}{ {\cal V} }
\newcommand{\cQ}{ {\cal Q} }
\newcommand{\cX}{ {\cal X} }
\newcommand{\cJ}{ {\cal J} }
\newcommand{\cU}{ {\cal U} }
\newcommand{\cL}{ {\cal L} }
\newcommand{\cP}{ {\cal P} }
\newcommand{\cN}{ {\cal N} }
\newcommand{\cW}{ {\cal W} }
\newcommand{\cZ}{ {\cal Z} }
\newcommand{\cS}{ {\cal S} }

\newcommand{\f}{\hspace{-3pt} \not\!}
\newcommand{\fD}{D\hspace{-13pt} \not\!\hspace{9pt}}
\newcommand{\fd}{\partial\hspace{-11.25pt} \not\!\hspace{9pt}}

\newcommand{\2}{ \, { 1 \over 2} \, }
\newcommand{\4}{ \, { 1 \over 4} \, }
\newcommand{\8}{ \, { 1 \over 8} \, }

\newcommand{\green}{\color[rgb]{0,0.6,0}}
\newcommand{\blue}{\color[rgb]{0,0,0.7}}
\newcommand{\black}{\color[rgb]{0,0,0}}
\newcommand{\orange}{\color[rgb]{1,0.5,0}}
\newcommand{\red}{\color[rgb]{1,0,0}}

%% --- Body (figure in root directory, no graphicspath needed) ---
Renormalization has become a familiar procedure of regulating divergent expressions and absorbing effects of
the regularization in the bare parameters of the theory, ensuring that the computed physical quantities are finite
and independent of the regularization. Analysis of universality of renormalization-group limit-cycle behavior of
quantum Hamiltonians requires some additional insight, focusing on the transformation of Hamiltonians, called
the renormalization-group transformation (RGT). We first explain how the RGT works in the context of a simplified
example that exhibits asymptotic freedom but lacks a cycle. Then we modify the example to include a cycle in a
generically quantum way.

Consider the Hamiltonian operator $H=H_0+H_I$. Eigenstates of $H_0$ are assumed known. They are denoted
by $|i\rangle$, $H_0|i\rangle = E_i |i\rangle$, where $i$ is a natural number. The states $|i\rangle$ form a basis in
the space of states. To most easily describe logarithmic renormalization effects, we assume that $E_i = e^i$,
$\ln e =1$. Our model Hamiltonian is assumed to have the matrix elements
\beq
\label{Hij}
H_{ij} \es \langle i| H |j\rangle \rs e^i \delta_{ij} - g e^{i+j \over 2} \ .
\eeq
The unit of energy is set to 1. The number $g$ is called a coupling constant.  The eigenvalue problem for $H$,
$H| \psi \rangle=E| \psi \rangle$, where $|\psi\rangle = \sum_{j =1}^\infty \psi_j |j\rangle$, projected on the eigenstates of
$H_0$, reads
\beq
\sum_{j=1}^\infty H_{ij} \psi_j = E \psi_i ~~~~{\rm for ~ all ~ i} \ .
\eeq
Thus, for all $i$,
\beq
\label{EP}
e^i \psi_i - g \, e^{i/2} \, C \es E \, \psi_i\ ,
\eeq
where the number $C = \sum_{j=1}^\infty e^{j/2} \psi_j$ does not depend on $i$. Therefore, $\psi_i = g e^{i/2} C /(e^i -E)$
and $C$ drops out of the eigenvalue condition,
\beq
\label{div}
1 \es g \sum_{i=1}^\infty {e^i \over e^i - E} \ .
\eeq
This condition cannot be solved for $E$ because the sum is infinite. If one imposes a cutoff on the energy of
the basis states, $E_i \leq e^N$, and demands that the resulting finite Hamiltonian matrix has a negative eigenvalue,
$E = - \cE$ with $\cE \ll 1$, which would correspond to some observable binding energy in a realistic theory, then
the eigenvalue condition takes the regulated form
\beq
\label{ren}
1 \es g \sum_{i=1}^N {1 \over 1 + \cE e^{-i}} \ .
\eeq
The regulated condition can be satisfied by adjusting the magnitude of the coupling constant, which leads to demanding that
\beq
g \tl g_N \sim 1/N \ .
\eeq
The cutoff-dependent coupling constant falls off as inverse of the number of dynamically coupled degrees of freedom,
which are the coefficients $\psi_i$ with $1 \leq i \leq N$. In terms of the energy cutoff, $E_i \leq \Lambda$ with $\Lambda =
e^N$ in the adopted units, one obtains
\beq
\label{renormalizedg}
g_\Lambda & \sim & { 1 \over \ln \Lambda} \ .
\eeq
One can thus rescue the model by ``renormalizing'' the coupling constant - adjusting it to the data for $\cE$ and the
cutoff $\Lambda$ used in computations, {\it cf.}~\cite{Jackiw:1991je}. One can also adopt the approach of matching
the computed eigenvalue $\cE$ with the eigenvalue of some other model for which this model is meant to be an
effective theory. Scattering results can be used in the matching as well. In general, if the renormalization procedure
exemplified above leads to finite results for observables in the limit $\Lambda \to \infty$, the theory is called
``renormalizable.'' However, it would be helpful~\cite{Dirac:1965tvn} to explain why fixing one observable by making
one parameter a function of the cutoff leads to cutoff independence of other observables.

We now proceed to the Wilsonian approach to renormalization which qualitatively explains the ``renormalizability''.
One begins with a study of how the matrix $H_{ij}$ of size $n \times n$, with $n \ll N$, ought to be defined in order
to have the same spectrum of eigenvalues $E \ll e^n$ as the cut off matrix $N \times N$ has when $N \to \infty$.
This is accomplished using the Gaussian elimination, which is the simplest form of the RGT for Hamiltonian matrices.
One begins from using Eq.~(\ref{EP}) to express $\psi_N$ in terms of $\psi_i$ with $i \leq N- 1$.
\beq
\label{EP1}
\psi_N = g e^{N/2} C_{N-1} \left( e^N - g e^N - E \right)^{-1} ,
\eeq
where $C_K = \sum_{k=1}^K e^{k/2} \psi_k$. The result for $\psi_N$ is inserted in the remaining equations with $i \leq N-1$,
which implies
\beq
\label{EP2}
e^i \psi_i - g \, e^{i/2} \, C_{N-1} - g^2 e^{i/2} e^N C_{N-1}\left( e^N - g e^N - E \right)^{-1}  \es E \, \psi_i\ .
\eeq
Equations with $i \leq N-1$ can now be written in the form similar to Eq.~(\ref{EP}),
\beq
\label{EP3}
e^i \psi_i - g_{N-1} \, e^{i/2} \, C_{N-1} \es E \, \psi_i\ ,
\eeq
with
\beq
g_{N-1} \es g + g^2 \left( 1 - g - Ee^{-N} \right)^{-1}  \ .
\eeq
For eigenvalues negligible in comparison with $e^N$, one obtains the set of equations with $i \leq N-1$ in which the coupling
constant $g$ is replaced by
\beq
g_{N-1} = RGT(g) = g/(1-g) \ .
\eeq
The Gaussian elimination can be applied many times, each time $g_{k-1} = RGT(g_k) = g_k/(1-g_k)$.
After $k = N-n$ such steps, one obtains
\beq
g_n \es {g_N \over 1 - g_N (N - n) } \ .
\eeq
By construction, the matrix of size $n \times n$ and coupling constant $g_n$ has the same spectrum of eigenvalues
$E \ll e^n(1-g_n)$ as the $N \times N$ matrix with the
coupling constant $g_N$. The recursion can be inverted,
\beq
g_N \es {g_n \over 1 + g_n (N-n) } \ .
\eeq
When one uses some cutoff $n_0$ to fit $g_0 = g_{n_0}$
to the observed low-energy spectrum, the inverted relation says that
the coupling constant in the matrix with the cutoff $N \to
\infty$ should be
\beq
\label{gNn0}
g_N \es {g_0 \over 1 + g_0 (N - n_0) } \ .
\eeq
In terms of the energy cutoffs $\Lambda = e^N$ and $\lambda_0
= e^{n_0}$, the coupling constant corresponding to $\Lambda$
ought to be
\beq
\label{gLambda}
g_\Lambda \es {g_0 \over 1 + g_0 \ln (\Lambda/\lambda_0)  } \ .
\eeq
Writing $1/g_0 = \ln (\lambda_0/\mu_0)$, which is always possible
for a real coupling constant, one arrives at
\beq
\label{gLambdamu}
g_\Lambda \es {1 \over \ln (\Lambda/\mu_0)  } \ ,
\eeq
which explains Eq.~(\ref{renormalizedg}). It also provides some understanding of the concept of dimensional
transmutation~\cite{Coleman:1973jx} in the case of asymptotic freedom~\cite{Politzer:1973fx,Gross:1973id}.

In the approach to renormalization we discussed first, one uses solutions of a theory to adjust the coupling
constant to the cutoff. In the second approach we discussed, one uses the RGT and discovers that the cutoff
dependent coupling constant is the mathematically required modification of the Hamiltonian, when one
introduces a cutoff on the free energy of the basis states. Moreover, the RGT was developed without knowing any
solutions to the theory.  This is an important distinction because realistic theories requiring renormalization
are complex and it is not straightforward to obtain their exact solutions.

Our description of the RGT exhibiting the limit-cycle behavior starts here. The RGT is applied to the
Hamiltonian $H$ whose matrix elements in the basis of eigenstates of $H_0$ apparently differ only
a little from the ones given in Eq.~(\ref{Hij}), see App. B 2. in \cite{Glazek:1993qs}.
 \beq
\label{Hijh}
H_{kl} \es \langle k| H |l\rangle \rs e^k \delta_{kl} -  ( g + i h s_{kl})  e^{k+l \over 2} \ ,
\eeq
Here $i = \sqrt{-1}$ and $s_{kl}$ is the sign of $k-l$ for $k \neq l$ and 0 for $k=l$. The added coupling
constant $h$ is arbitrarily small and in the limit $h \to 0$ the Hamiltonian of Eq.~(\ref{Hijh}) seems to
approach the one of Eq.~(\ref{Hij}). However, the skew-symmetric imaginary addition produces the RGT
effects that can be large no matter how small $h$ is. Both $g$ and $h$ contribute to the model ultraviolet
behavior. The Gaussian elimination, carried out similarly to the asymptotically free model discussed earlier,
produces the recursion~\cite{Glazek:2002hq,Glazek:2003ljf}
\beq
g_{n-1} \es {g_n + h^2 \over 1 - g_n - E/E_n} \ .
\eeq
The coupling $h$ does not vary with $n$. For as long as the ratio $E/E_n$ is so small that its presence can
be neglected, or when the eigenvalue $E$ is 0, the recursion implies
\beq
g_n \es h \tan \left[ \, \arctan(g_N/h) + (N-n) \arctan h \, \right] \ .
\eeq
In that case, for $h = \tan(\pi/p)$ with a natural $p > 2$, the coupling constant $g_{n- kp}$ with natural $k$
is the same as $g_n$. This means that for any $n$ the set of coupling constants $g_{n-1}$, $g_{n-2}$,...,$g_{n-p}$
repeats itself as $g_{n-p-1}$, $g_{n-p-2}$,...,$g_{n-2p}$, and so on. We call the set of $p$ consecutive coupling
constants a cycle. In each cycle, there is a maximal coupling constant, of positive value. The subsequent coupling
constant in the recursion that lowers the running cutoff, is negative, of the largest negative value in the cycle.
There is only one such jump in every cycle. Every such jump is associated with a bound state, one per cycle.
The corresponding binding energies form a geometric sequence with the common ratio $e^p$. Positive
eigenvalues in a single cycle appear with a common ratio very close to $e^{p/(p-1)}$. But there is an unsolved
puzzle concerning the consecutive ratios. Their product is much closer to $e^{2p/(p-1)}$ than each one is to
$e^{p/(p-1)}$, see Sec. VII in~\cite{Glazek:2003ljf}.

Analysis of universality in the limit-cycle behavior is carried out similarly to the analysis in case of a fixed point behavior
except that there are $p$ fixed points to consider in a cycle. They are fixed with respect to the RGT applied $p$
times instead of one time. The coupling constants $g_n$ depend on the ratio $x=E/E_n$. Their dependence on
$x$ can be expanded in a Taylor series. The coefficients of powers of $x$, called Wegner's operators~\cite{Wegner:1972my},
are extracted from an eigenvalue condition for the linearized form of the RGT carried out $p$ times, keeping some
number of terms in the Taylor expansion. By inspection of how these terms vary with $n$ for every coupling constant
in a cycle, one finds that there are no operators that increase as $n$ is lowered. This is described by saying that there
are no Wegner's relevant operators, or operators that have positive anomalous dimensions and grow. One operator
per $g$ in the cycle stays constant, corresponding to the change of the coupling constant. It has anomalous dimension
0 and is called marginal. There are operators that decrease as inverse powers of $e^p$. These are called irrelevant.
Their anomalous dimensions are negative integers, independently of their position on the cycle.

Having understood the running of the model with the cutoff, one can design the initial Hamiltonian to cancel the
largest irrelevant operators. In consequence, the spectrum explicitly approaches the geometric series even for the
cutoffs small enough to not exceed the magnitudes manageable numerically. The model can thus be helpful in
developing a theory of the experimental methods for tuning physical systems to the cycle, especially the atomic
systems. It should be mentioned at this point that the limit cycle takes a convoluted form for rational $p$. When $p$
is irrational, the RGT leads to the chaotic dependence of coupling constants $g_n$ on the running cutoff $n$.

In the inverted form, analogous to Eq.~(\ref{gNn0}) in the asymptotically free model, using $\Lambda = e^N$, one obtains
\beq
\label{UVLC}
g_\Lambda \es h \tan \left[ \, \arctan(g_0/h) - \ln(\Lambda/\lambda_0) \arctan h \, \right] \ .
\eeq
Thus the coupling constant $g_\Lambda$ does not have a limit  when $\Lambda \to \infty$. There is no ultimate theory
with one coupling constant. Instead, one has a repeating sequence of theories with varying coupling constants and the same
low energy spectrum. We also observe that the space of real symmetric Hamiltonian matrices is of measure 0 in the space of
complex Hermitian matrices. Since quantum mechanics unavoidably involves complex numbers, to preserve probability, one
may expect that the cyclic renormalization-group behavior may actually be commonplace for quantum Hamiltonians requiring renormalization. Additional observation concerning Eq.~(\ref{UVLC}) is that it is obtained from an asymptotically free model
in which the coupling constant $h$ equals 0. But no matter how small $h$ is, even orders of magnitude smaller than any
{\it gedanken} experiment could detect, there is a number of the RGT steps within just one cycle that bring the model coupling
constant $g$ to the values comparable with the magnitudes known for the coupling constant $g_{\rm QCD}$ in the observed
$\alpha_{\rm QCD}$. At that stage of the renormalization-group evolution with the running cutoff, the coupling constant $g$
dominates over $h$. The presence of an extremely small $h$ is not noticeable. Hence, one can entertain the idea that the observed asymptotically free behavior of QCD may actually be a part of a larger, cyclic pattern of a renormalization-group behavior of
some more fundamental theory~\cite{Glazek:2006mg}.

We ought to stress, however, that the limit-cycle behavior of quantum Hamiltonians can occur without imaginary interaction
terms. Major example is provided by the theory of three-body systems with short range interactions and large scattering
length~\cite{Bedaque:1998kg,Braaten:2004rn}, describing the Efimov states~\cite{Efimov:1971zz} in nuclear physics and
more recently in atomic physics~\cite{Braaten:2006vd,Kunitski:2015qth}. Also, it has been suggested~\cite{Braaten:2003eu}
that QCD may involve an infrared limit cycle when the masses of quarks up and down are slightly increased; the pion mass
would increase, the deuteron binding energy would approach zero, and triton would develop infinitely many excited states
with binding energies forming a geometric sequence converging on 0. If one derived the nucleon-meson effective theory
from QCD, the three-body eigenvalue problem could be solved using numerical methods with precise control on the cutoff dependence~\cite{Mohr:2005pv}. Closely related but simpler example than the three-body model, is provided by the two-body Schr\"odinger Hamiltonian with an attractive potential inversely proportional to the square of the distance between the bodies.
One can solve for the spectrum of the regulated eigenvalue problem and subsequently adjust the interaction parameter to the
data and the regularization method one uses~\cite{Braaten:2004pg}, to make the spectrum not depend on the regularization.
One can also  compute the RGT~\cite{Dawid:2017ahd} and find out how the Hamiltonian ought to depend on the chosen
regularization to yield the regularization-independent results.

The question is then what are the prospects of using the renormalization-group methods for studying the renormalization-group
behavior, potentially of the limit-cycle type, in theories of complexity comparable to that of QCD. To contribute to such studies,
the RGT we discussed so far
is not sufficient. Cutoff on energy is not compatible with the special-relativity symmetry of QCD, because boosts can change
energy by the arbitrarily large amounts. The QCD basis states made of quarks and gluons have a highly degenerate, continuous
spectrum, as eigenstates of the free part of the canonical QCD Hamiltonian. The Gaussian elimination does not easily apply.
Expansion of effective Hamiltonians in a series of powers of even arbitrarily small coupling constant $g$, involves large terms
because of the small free-energy differences in denominators. Introduction of artificially large energy gaps between the eliminated
and retained states, used in~\cite{Wilson:1965zzb,Wilson:1970tp}, needs to be replaced with a more robust concept for discussion
of the relativistic bound states of quarks and gluons~\cite{Wilson:1994fk}.

Such concept is found in the new renormalization procedure for Hamiltonians~\cite{Glazek:1993rc,Glazek:1994qc}, which is
based on the idea of limiting the changes of energy an effective interaction can cause, instead of limiting the energy itself.
Mathematically, it means transforming the Hamiltonian matrix of a theory from its canonical form to the band diagonal form.
The width of the band, say $\lambda$, plays the role of a renormalization-group parameter. Transformations of the banded
Hamiltonian $H_\lambda = H_0 + H_{I\lambda}$ from one value of $\lambda$ to another form a group, called the similarity
renormalization group (SRG). Perturbative computation of $H_\lambda$ does not involve small energy denominators,
irrespective of the degeneracies, because the smallest difference of energies involved in the transformation is $\lambda$.
Since the largest jump in energy that $H_{I\lambda}$ can make, as measured in terms of the eigenvalues of $H_0$, is
$\lambda$, one needs a product of $n$ interactions, $H_{I\lambda}^n$, to reach energy $E_i + n \lambda$ starting from
$E_i$, where $E_i$ is a finite eigenvalue of $H_0$. To return to the energies of order $E_i$, say $E_j$, one needs another
$H_{I\lambda}^n$. If the cutoff on the canonical Hamiltonian matrix is $\Lambda$, one needs the power $n = 2\Lambda/
\lambda$ of $H_{I\lambda}$ to reach $\Lambda$ in a perturbative computation of the effective Hamiltonian matrix elements
$H_{\lambda ij}$. In the limit $\Lambda \to \infty$, the power $n$ grows to infinity, which means that the only dependence of
the effective theory perturbative predictions on $\Lambda$ comes from the explicit dependence of the matrix elements
$H_{\lambda ij}$ on $\Lambda$. Thus, it is sufficient to introduce the counterterms in the canonical $H$ with cutoff $\Lambda
\to \infty$ that counter the explicitly calculable divergences in $H_{\lambda ij}$, to obtain the $H_{\lambda}$ whose perturbative
predictions are free from regularization dependence to all orders of the weak-coupling expansion. Since the banded
Hamiltonian $H_\lambda$ limits the energy changes and does not limit the energy itself, it {\it in principle} allows the
dynamics of states of one energy scale to be related to the dynamics of states of a Lorentz-transformed energy scale,
both contained within the band but at different energies.

For the purposes of the SRG, one can employ Wegner's double-commutator flow of Hamiltonian matrices~\cite{Wegner:1994fdg,
Mielke:1998ag,Wegner_2006,Kehrein:2006ti}. The name ``double-commutator'' originates from the form of the equation for
$H = D + V$,
\beq
{dH \over ds} \es \left[ \left[ D, V \right], H \right] \ ,
\eeq
where $D$ and $V$ denote the diagonal and off-diagonal parts of the Hamiltonian matrix, respectively. The parameter
$s$ ranges {\it a priori} from 0 to $\infty$, corresponding to the inverse of the band-width squared, $s \sim 1/\lambda^2$.
The input matrix $H$ provides the initial condition at $s=0$, which corresponds to the band of  an infinite width. As $s$
increases, one obtains a matrix with a decreasing energy band-width $\sim 1/\sqrt{s}$~\cite{Wegner:1994fdg}. When the
Hamiltonian dependence on the parameter $s$ is expanded in a series of powers of the interaction strength, the band
formation can be seen already in the first-order solution for the matrix $V$, $V_{ij}(s) = V_{ij}(0) \exp[-s (E_i - E_j)^2]$,
where $E_k = D_{kk}(s=0)$. One can also observe that the double-commutator structure of the right-hand side of the
Wegner equation guarantees satisfaction of the cluster decomposition principle~\cite{Weinberg:1995mt} by the computed
effective Hamiltonian matrices. Similarly, the double-commutator equations used in the SRG procedure for renormalization
of Hamiltonians in quantum field theory guarantee that no disconnected interactions are generated.

Wegner's equation and similar double-commutator equations have been extensively used in nuclear physics computations
of softened NN interactions~\cite{Bogner:2006pc,Bogner:2007qb,Jurgenson:2007td,Glazek:2008pg,Anderson:2010aq,
Bogner:2009bt}, including the discussion of universality~\cite{Dainton:2013axa}. The double-commutator equation is used
to outline the SRG computational strategy for bound states in theories with asymptotic freedom~\cite{Glazek:1997gt}.
Videos illustrating the SRG limit-cycle behavior of the discussed earlier model Hamiltonian with coupling constant
$h$~\cite{bib9lc}, were generated using Wegner's equation, too. Alteration of double-commutator flow equations, needed
to improve convergence in the perturbative computations of $H_\lambda$, are described in~\cite{Glazek:2002iw}. Recent
examples of application of the SRG, in its operator form called the renormalization group procedure for effective particles
(RGPEP), directly to QCD, concern masses of mesons and baryons made of quarks $c$ and $b$~\cite{Serafin:2018aih}. The
quark masses are fitted to the quarkonium spectra, using optimized width $\lambda$. As a result, the formulas for masses
of baryons contain no free parameters. The early estimates for the masses of $\Omega_{ccc}$ obtained using the RGPEP
are recently confirmed by the lattice calculations claimed to be currently most precise~\cite{Dhindsa:2024erk}. The agreement
may be accidental, but it provides a point of focus for building a bridge between these methods, both of which are based on
the first principles. Discussion of how the RGPEP may be used to describe confinement is offered in~\cite{Serafin:2023pkf}.

Given the possibility that QCD, nuclear theory, and theory of atoms may involve universal cyclic behaviors, and taking into
account the expectation that new generations of computers will achieve greatly increased computational power, it seems
conceptually useful to outline here a mathematical scheme for obtaining effective-theory Hamiltonian matrices from the
underlying quantum dynamics. The issue is that the effective matrix applies to a different kind of degrees of freedom from
the ones in the underlying theory while one has to rely on quantum mechanics rather than classical intuition in approaching
the issue.

In the context of strong interactions, the RGT one needs should provide the Hamiltonian matrices corresponding
to nucleons and mesons~\cite{Erkelenz:1974uj,Machleidt:1995km} starting from the matrix describing the dynamics of
effective quarks and gluons. Figure  \ref{figQCDandNP} illustrates the procedure we have in mind~\cite{Glazek:2002jg}.
\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{QCDtoNP}
\caption{Sequence of renormalization group computations leading to the Hamiltonian
              matrix of an effective theory, square 4, starting from the Hamiltonian matrix of
              the underlying theory, square 1, {\it cf.}~\cite{Glazek:2002jg}.}\label{figQCDandNP}
\end{figure}
The banded matrix of $H_\lambda$ in QCD is represented by the large square labeled by 1 in Fig.~\ref{figQCDandNP}. The
figure is adapted from~\cite{Glazek:2002jg}. The matrix 1 is divided into square window matrices that can be diagonalized
separately. Changing the initial basis states to the basis states formed by the resulting eigenstates of the window matrices,
which are the precursors of the nucleon and meson states, one obtains the matrix elements shown in matrix 2. The varying
width of the band in matrix 2 indicates that the middle eigenstates of windows of matrix 1 approximate eigenstates of the
whole matrix 1 well. The width in the matrix 2 widens near the corners of the windows because the diagonalization of the
window matrices misses matrix elements outside the windows in matrix 1. New division of the matrix 2 into window matrices
is shown in the square 3. Corners of the new windows are made of the matrix elements between the middle eigenstates
of the windows in matrix 1. Diagonalization of the new window matrices shown in matrix 3 yields the band in the matrix 4.
The lowest window matrix in the square 4 stands for the theory of nucleons interacting through effective potentials, below
the threshold for producing pions. Coupling to the next higher window matrix in the square 4 includes emission and
absorption of a pion, etc.

\clearpage
\endgroup
