# Characterization of random features of chaotic eigenfunctions in unperturbed
basis

Jiaozi Wang 111 Email address<EMAIL_ADDRESS>and Wen-ge Wang 222 Email
address<EMAIL_ADDRESS>Department of Modern Physics, University of
Science and Technology of China, Hefei, 230026, China

###### Abstract

In this paper, we study random features manifested in components of energy
eigenfunctions of quantum chaotic systems, given in the basis of unperturbed,
integrable systems. Based on semiclassical analysis, particularly on Berry’s
conjecture, it is shown that the components in classically allowed regions can
be regarded as Gaussian random numbers in certain sense, when appropriately
rescaled with respect to the average shape of the eigenfunctions. This
suggests that, when a perturbed system changes from integrable to chaotic,
deviation of the distribution of rescaled components in classically allowed
regions from the Gaussian distribution may be employed as a measure for the
“distance” to quantum chaos. Numerical simulations performed in the LMG model
and the Dicke model show that this deviation coincides with the deviation of
the nearest-level-spacing distribution from the prediction of random-matrix
theory. Similar numerical results are also obtained in two models without
classical counterpart.

##  I Introduction

A commonsense in the field of quantum chaos is that energy eigenfunctions
(EFs) of chaotic systems should show certain random feature Haake ; CC94book ;
Berry77 ; Berry91 , though their Hamiltonian matrices are deterministic and
some of them even may show a sparse structure. This property has vast
applications in various fields Stockmann ; Meredith98 ; Iz96 ; Prosen03 ;
pre-98 ; Gnutzmann10 ; Kaplan09 ; pre02-LMG ; qd2 ; opt1 ; nuc3 . In
particular, it is of relevance to thermalization Izlev11 ; Haake12 ; Deutsch91
; Sr94 ; Rigol08 ; Rigol12 ; Gemmer15 ; ETH-review ; ETH-Chaos , a topic which
has attracted renewed interest in recent years.

According to Berry’s conjecture, for EFs of chaotic systems expressed in the
configuration space, their components in classically allowed regions can be
regarded as being given from certain Gaussian random numbers Berry77 . Based
on this conjecture, it is natural to expect that, when expanded in the bases
of unperturbed integrable systems, the EFs should show certain random feature
as well within appropriate regions. Indeed, numerical simulations revealed
such a feature for main bodies of EFs (see, e.g., Ref.Buch82 ). However, more
detailed study showed that the distribution of components of EFs usually
exhibit notable deviation from the Gaussian distribution, which is predicted
by the random-matrix theory (RMT) (see, e.g., Ref.Meredith98 ).

More recently, numerical simulations show that, if EFs are rescaled with
respect to their average shape, the above-discussed deviation can be
considerably reduced EFchaos-Benet03 ; EFchaos-WW . This gives a clue to a
solution to a long-standing problem in the field of quantum chaos, that is, in
which way statistical properties of EFs may be employed to give a quantitative
measure for the “distance” to chaos. The above-mentioned numerical simulations
suggest that deviation of the distribution of rescaled components of EFs from
the Gaussian distribution should be a candidate for such a measure of
“distance”. However, presently, the situation is not completely clear, because
in some cases this measure shows notable deviations from results obtained from
statistical properties of spectra, e.g., from deviation of the nearest-level-
spacing distribution from the prediction of the RMT EFchaos-Benet03 ; EFchaos-
WW .

In this paper, based on semiclassical analysis, particularly on the Berry’s
conjecture, we study random features manifested in components of EFs of
chaotic systems in integrable bases. Our analysis shows that the distribution
of the components in classically-allowed regions indeed should have a Gaussian
form, under a rescaling procedure which is more appropriate than that adopted
in Refs.EFchaos-Benet03 ; EFchaos-WW . Our numerical simulations performed in
the Lipkin-Meshkov-Glick (LMG) model and the Dicke model show that, adopting
this new rescaling procedure, deviation of the distribution of components from
the Gaussian distribution coincides quite well with that obtained from the
statistics of spectra. We also study some models without any classical
counterpart and find similar results.

The paper is organised as follows. In Sec.II, a detailed semiclassical
analysis is carried out for random features manifested in components of EFs of
chaotic systems in integrable bases. Numerical simulations in two models with
classical counterparts are discussed in Sec.III. Then, in Sec.IV, we discuss
numerical simulations performed in two models without classical counterpart.
Finally, conclusions and discussions are given in Sec.V.

##  II Random features of chaotic EFs

In Sec.II.1, based on semiclassical analysis we discuss random features of
chaotic EFs. Then, making use of results obtained, we discuss a quantitative
characterization of the random feature in Sec.II.2.

###  II.1 Semiclassical analysis of chaotic EFs

Consider a quantum system, which has an $f$-dimensional classical counterpart,
with a Hamiltonian

$H=H_{0}+\lambda V,$ (1)

where $H_{0}$ indicates the Hamiltonian of an integrable system and $V$ is a
perturbation. Within certain regime of the parameter $\lambda$, the classical
counterpart of the system $H$ undergoes a chaotic motion. In this section, we
consider a chaotic system $H$. In terms of action-angle variables, $H_{0}$ is
written as

$H_{0}=\boldsymbol{d}\cdot\boldsymbol{I}+c_{0},$ (2)

where $\boldsymbol{I}=(I_{1},I_{2},\cdots,I_{f})$ is the action variable,
$\boldsymbol{d}$ is a parameter vector,
$\boldsymbol{d}=(d_{1},d_{2},\cdots,d_{f})$, and $c_{0}$ is a constant
parameter.

In the quantum case, we use $|\boldsymbol{n}\rangle$ to denote the eigenbasis
of $\boldsymbol{I}$, with
$\boldsymbol{I}|\boldsymbol{n}\rangle=\boldsymbol{I_{n}}|\boldsymbol{n}\rangle$,
where $\boldsymbol{n}=(n_{1},n_{2},\cdots,n_{f})$ is an integer vector and
$\boldsymbol{I_{n}}=\boldsymbol{n}\hbar$. The Hamiltonian $H_{0}$ has a
diagonal form in this basis with eigenvalues denoted by
$E_{\boldsymbol{n}}^{0}$,

$H_{0}|\boldsymbol{n}\rangle=E_{\boldsymbol{n}}^{0}|\boldsymbol{n}\rangle.$
(3)

We use $|E_{\alpha}\rangle$ to denote eigenstates of $H$ with eigenvalues
$E_{\alpha}$ in energy order,

${H}|E_{\alpha}\rangle=E_{\alpha}|E_{\alpha}\rangle.$ (4)

The expansion of $|E_{\alpha}\rangle$ in the basis $|\boldsymbol{n}\rangle$ is
written as

$|E_{\alpha}\rangle=\sum_{\boldsymbol{n}}C_{\alpha\boldsymbol{n}}|\boldsymbol{n}\rangle,$
(5)

with $C_{\alpha\boldsymbol{n}}=\langle\boldsymbol{n}|E_{\alpha}\rangle$. Below
in this section, we discuss random features manifested in the components
$C_{\alpha\boldsymbol{n}}$ and their statistical properties in chaotic
systems.

In terms of the wave functions of $|E_{\alpha}\rangle$ and
$|\boldsymbol{n}\rangle$ in the momentum space, denoted by
$\psi_{\alpha}(\boldsymbol{p})$ and
$\psi_{\boldsymbol{n}}^{0}(\boldsymbol{p})$, respectively, the components
$C_{\alpha\boldsymbol{n}}$ are written as

$C_{\alpha\boldsymbol{n}}=\int(\psi_{\boldsymbol{n}}^{0}(\boldsymbol{p}))^{*}\psi_{\alpha}(\boldsymbol{p})d\boldsymbol{p}.$
(6)

Generically, a wave function $\psi_{\alpha}(\boldsymbol{p})$ can be written in
the following form,

$\psi_{\alpha}(\boldsymbol{p})=A_{\alpha}(\boldsymbol{p})\sqrt{\Pi_{\alpha}(\boldsymbol{p})},$
(7)

where $\Pi_{\alpha}(\boldsymbol{p})$ indicates local average of
$|\psi_{\alpha}(\boldsymbol{p})|^{2}$. Then, $C_{\alpha\boldsymbol{n}}$ is
written as

$\displaystyle C_{\alpha\boldsymbol{n}}=\int
A_{\alpha}(\boldsymbol{p})(\psi_{\boldsymbol{n}}^{0}(\boldsymbol{p}))^{*}\sqrt{\Pi_{\alpha}(\boldsymbol{p})}d\boldsymbol{p}.$
(8)

According to the Berry’s conjecture Berry77 , in a chaotic system the quantity
$A_{\alpha}(\boldsymbol{p})$ should have random phases. This implies that the
components $C_{\alpha\boldsymbol{n}}$ can be effectively regarded as some
random numbers.

In realistic physical models, the average shape of
$|C_{\alpha\boldsymbol{n}}|^{2}$ is usually not uniform with respect to the
perturbed and unperturbed energies. Due to this nonuniformity, the statistical
distribution of the components $C_{\alpha\boldsymbol{n}}$ can not have a
Gaussian shape Meredith98 . But, if they are rescaled such that the effect of
average shape of EFs is appropriately taken into account, it should be
possible for their distribution to have a Gaussian form. Below, we derive a
semiclassical expression for the average shape of
$|C_{\alpha\boldsymbol{n}}|^{2}$ that is suitable for this purpose.

The Wigner function supplies a useful tool in semiclassical analysis of
eigenstates. We use $\psi_{\alpha}(\boldsymbol{r})$ and
$\psi_{\boldsymbol{n}}^{0}(\boldsymbol{r})$ to indicate the wave functions of
$|E_{\alpha}\rangle$ and $|\boldsymbol{n}\rangle$ in the coordinate space,
respectively. The Wigner function corresponding to
$\psi_{\alpha}(\boldsymbol{r})$, denoted by
$W_{\alpha}(\boldsymbol{p},\boldsymbol{q})$, is written as

$W_{\alpha}(\boldsymbol{p,q})=\frac{1}{(2\pi\hbar)^{f}}\int_{-\infty}^{\infty}\psi_{\alpha}^{*}(\boldsymbol{q}+\frac{\boldsymbol{r}}{2})\psi_{\alpha}(\boldsymbol{q}-\frac{\boldsymbol{r}}{2})e^{i\boldsymbol{p\cdot
r}/\hbar}d\boldsymbol{r},$ (9)

and similar for the Wigner function corresponding to
$\psi_{\boldsymbol{n}}^{0}(\boldsymbol{r})$, denoted by
$W_{\boldsymbol{n}}^{0}(\boldsymbol{p},\boldsymbol{q})$. As is known, in a
chaotic system, the averaged Wigner function, with average taken within
certain small regions of the phase space, has the following expression,
Berry77 ; Berry91 ; Voros76 ; Voros77

$\overline{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})=\frac{\delta(H(\boldsymbol{p},\boldsymbol{q})-E_{\alpha})}{S({E_{\alpha}})},$
(10)

where $S(E)$ represents the area of an energy surface with
$H(\boldsymbol{p},\boldsymbol{q})=E$,

$S({E})=\int
d\boldsymbol{p}d\boldsymbol{q}\delta(E-H(\boldsymbol{p},\boldsymbol{q})).$
(11)

Equation (10) gives that

$\Pi_{\alpha}(\boldsymbol{p})=\frac{1}{S({E_{\alpha}})}\int\delta(E_{\alpha}-H(\boldsymbol{p},\boldsymbol{q}))d\boldsymbol{q}.$
(12)

Equation (10) implies that most eigenstates within a narrow energy window in a
chaotic system should have close shapes. Therefore, when computing the average
shape of $|C_{\alpha\boldsymbol{n}}|^{2}$ for the purpose discussed above, one
may perform an average within such a narrow energy window. For the convenience
in discussion, we write a coarse-grained $\delta$-function as
$\delta_{\epsilon}(E)$,

$\delta_{\epsilon}(E)=\begin{cases}\frac{1}{\epsilon}&E\in[-\frac{\epsilon}{2},\frac{\epsilon}{2}],\\\
0&{\rm otherwise},\end{cases}$ (13)

where $\epsilon$ is a small parameter. The choice of energy window $\epsilon$
should satisfy the following requirements: It is small in the classical case
such that the energy surface almost does not change within the window, while,
it is sufficiently large in the quantum case such that many energy levels are
included within the window.

Then, the average shape of EFs, denoted by
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$, is computed by

$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle=\frac{1}{N_{E_{\alpha}}}\sum_{\alpha^{\prime}}|C_{\alpha^{\prime}\boldsymbol{n}}|^{2}\delta_{\epsilon}(E_{\alpha^{\prime}}-E_{\alpha}),$
(14)

where

$N_{E_{\alpha}}=\sum_{\alpha{{}^{\prime}}}\delta_{\epsilon}(E_{\alpha{{}^{\prime}}}-E_{\alpha}).$
(15)

In order to derive an explicit expression for
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$, we make use of the following
well-known expression of $|C_{\alpha\boldsymbol{n}}|^{2}$,

$|C_{\alpha\boldsymbol{n}}|^{2}=(2\pi\hbar)^{f}\int
d\boldsymbol{p}d\boldsymbol{q}W_{\alpha}(\boldsymbol{p},\boldsymbol{q})W_{\boldsymbol{n}}^{0}(\boldsymbol{p},\boldsymbol{q}).$
(16)

Let us divide the phase space into small cells, denoted by $c_{\sigma}$ with a
label $\sigma$, each having a volume $\delta\Omega$, meanwhile, keep the ratio
$\delta\Omega/\hbar^{f}$ large such that there are many quantum states “lying”
within each small cell. Then, $|C_{\alpha\boldsymbol{n}}|^{2}$ is written as

$|C_{\alpha\boldsymbol{n}}|^{2}=(2\pi\hbar)^{f}\sum_{\sigma}\int_{c_{\sigma}}d\boldsymbol{p}d\boldsymbol{q}W_{\alpha}(\boldsymbol{p},\boldsymbol{q})W_{\boldsymbol{n}}^{0}(\boldsymbol{p},\boldsymbol{q}).$
(17)

Substituting Eq.(17) into Eq.(14) and performing the summation over the
perturbed states $|E_{\alpha^{\prime}}\rangle$, one gets that

$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle=(2\pi\hbar)^{f}\sum_{\sigma}\int_{c_{\sigma}}d\boldsymbol{p}d\boldsymbol{q}\langle{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})\rangle
W_{\boldsymbol{n}}^{0}(\boldsymbol{p},\boldsymbol{q}),$ (18)

where $\langle{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})\rangle$ indicates
the average of ${W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})$ over perturbed
states within a small energy window $\epsilon$. As many energy levels are
included within the window $\epsilon$,
$\langle{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})\rangle$ should vary slowly
within each small cell $c_{\sigma}$, such that

$\langle{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})\rangle\simeq\langle{W}_{\alpha}(\boldsymbol{p}_{\sigma},\boldsymbol{q}_{\sigma})\rangle\quad\text{
for }(\boldsymbol{p},\boldsymbol{q})\in c_{\sigma},$ (19)

where $(\boldsymbol{q}_{\sigma},\boldsymbol{p}_{\sigma})$ indicates the center
of $c_{\sigma}$. Then, Eq.(18) gives that

$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle\simeq(2\pi\hbar)^{f}\sum_{\sigma}\langle{W}_{\alpha}(\boldsymbol{p}_{\sigma},\boldsymbol{q}_{\sigma})\rangle\overline{W}^{0}_{\boldsymbol{n}}(\boldsymbol{p}_{\sigma},\boldsymbol{q}_{\sigma})\delta\Omega,$
(20)

where
$\overline{W}^{0}_{\boldsymbol{n}}(\boldsymbol{p}_{\sigma},\boldsymbol{q}_{\sigma})$
represents the average of the Wigner function of the integrable system within
the cell $c_{\sigma}$,

$\overline{W}^{0}_{\boldsymbol{n}}(\boldsymbol{p}_{\sigma},\boldsymbol{q}_{\sigma})=\frac{1}{\delta\Omega}\int_{c_{\sigma}}W_{\boldsymbol{n}}^{0}(\boldsymbol{p},\boldsymbol{q})d\boldsymbol{p}d\boldsymbol{q}.$
(21)

Due to the classical smallness and quantum mechanical largeness of the energy
windows $\epsilon$ discussed previously,
$\langle{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q})\rangle$ in Eq.(18) obeys
Eq.(10) in an approximate way, with

$\langle
W_{\alpha}(\boldsymbol{p},\boldsymbol{q})\rangle\simeq\overline{W}_{\alpha}(\boldsymbol{p},\boldsymbol{q}),$
(22)

and its dependence on $\epsilon$ can be neglected. It is known that Berry77

$\overline{W}_{\boldsymbol{n}}^{0}(\boldsymbol{p},\boldsymbol{q})=\frac{\delta(\boldsymbol{I}(\boldsymbol{p},\boldsymbol{q})-\boldsymbol{I_{n}})}{(2\pi)^{f}}.$
(23)

Substituting Eqs.(22) and (23) into Eq.(20) and noting that the smallness of
the cells $c_{\sigma}$ enables one to change the summation over $\sigma$ back
to the integration over phase space, one gets the following semiclassical
expression,

$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle\simeq\hbar^{f}\
\Pi(E_{\alpha},\boldsymbol{I_{n}}),$ (24)

where

$\displaystyle\Pi(E,\boldsymbol{I})=\frac{S(E,\boldsymbol{I})}{S({E})},$ (25)
$\displaystyle S(E,\boldsymbol{I})=\int
d\boldsymbol{p}d\boldsymbol{q}\delta(E-H(\boldsymbol{p},\boldsymbol{q}))\delta(\boldsymbol{I}-\boldsymbol{I}(\boldsymbol{p},\boldsymbol{q})).$
(26)

Here, $S(E,\boldsymbol{I})$ indicates the overlap of the energy surface of
$H(\boldsymbol{p},\boldsymbol{q})=E$ and the torus of
$\boldsymbol{I}(\boldsymbol{p},\boldsymbol{q})=\boldsymbol{I}$. Since Eq.(10)
works for classically allowed regions only, so does Eq.(24). Sometimes,
quantities like $\Pi(E,\boldsymbol{I})$ are called classical analog of
averaged EFs EFchaos-Benet03 ; EFchaos-Benet00 ; EFchaos-Borgonovi98 .

Finally, we consider rescaled components denoted by
$R_{\alpha\boldsymbol{n}}$, defined by

$R_{\alpha\boldsymbol{n}}=\frac{C_{\alpha\boldsymbol{n}}}{\sqrt{\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle}}.$
(27)

Discussions given above show that this quantity $R_{\alpha\boldsymbol{n}}$ can
be regarded as a Gaussian random number with mean zero. Note that
$\langle|R_{\alpha\boldsymbol{n}}|^{2}\rangle=1$ according to Eq.(27).

###  II.2 A measure for “distance” to quantum chaos

Let us use $f(R)$ to denote the distribution of $R_{\alpha\boldsymbol{n}}$.
According to results given in the above section, for a chaotic system, $f(R)$
should have a Gaussian form, i.e.,

$f(R)=f_{G}(R),$ (28)

where $f_{G}(R)$ is the Gaussian distribution,

$f_{G}(R)=\frac{1}{\sqrt{2\pi}}\exp(-R^{2}/2),$ (29)

In the RMT, the Gaussian distribution is predicted directly for components of
EFs Haake . But, for Hamiltonians in realistic models with chaotic classical
counterparts, as discussed above, it is the distribution of the rescaled
components $R_{\alpha\boldsymbol{n}}$ that should have a Gaussian form. On the
other hand, in a nearly integrable system, the quantity
$A_{\alpha}(\boldsymbol{p})$ on the rhs of Eq.(8) does not have random phases
and EFs with close energies may have quite different shapes. As a result, the
distribution $f(R)$ in nearly integrable systems should usually show notable
deviation from $f_{G}(R)$.

The above discussions suggest that deviation of $f(R)$ from $f_{G}(R)$ may be
employed as an measure for the “distance” to quantum chaos. In order to
quantitatively characterize the deviation, one may consider a quantity
$\Delta_{EF}$ defined by

$\Delta_{EF}=\int|I_{f}(R)-I_{f_{G}}(R)|dR,$ (30)

where $I_{f}(R)$ and $I_{f_{G}}(R)$ indicate the cumulative distributions of
$f(R)$ and $f_{G}(R)$, respectively, e.g.,
$I_{f}(R)=\int_{-\infty}^{R}drf(r)$. As is well known, cumulative
distributions usually exhibit less fluctuations compared with the origin
distributions.

Figure 1:  (a): The average shape of EFs,
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$, in the chaotic regime of the
LMG model with $\Omega=500$ and $\lambda=1$. The average is taken over 500 EFs
of $|E_{\alpha}\rangle$ with $\epsilon=0.4175$. (b):
$\Pi_{N}({E_{\alpha},\boldsymbol{I_{n}}})$, which is the normalized
$\Pi(E_{\alpha},\boldsymbol{I_{n}})$, as a classical analog of
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$ [see Eq.(24)].

In the field of quantum chaos, the most-often used criterion for quantum chaos
is given by statistical properties of spectra, e.g., by closeness of the
nearest-level-spacing distribution $P(s)$ to the prediction of RMT Haake . It
is known that the following distribution $P_{W}(s)$, which is obtained from
Wigner’s surmise,

$P_{W}(s)=\frac{\pi}{2}s\exp(-\frac{\pi}{4}s^{2}),$ (31)

gives a good approximation to the nearest-level-spacing distribution of the
Gaussian orthogonal ensemble (GOE) in the large size limit. Quantitatively,
the above-discussed closeness can be characterized by the following quantity
$\Delta_{E}$,

$\Delta_{E}=\int|I_{P}(s)-I_{P_{W}}(s)|ds,$ (32)

where $I_{P}(s)$ and $I_{P_{W}}(s)$ are cumulative distributions of $P(s)$ and
$P_{W}(s)$, respectively.

Figure 2: Similar to Fig.1, but for the Dicke model with $N=500$, $\lambda=1$,
and $\epsilon=0.007$.

In previous numerical simulations, deviation of the distribution of another
rescaled components, denoted by $R_{\alpha\boldsymbol{n}}^{\prime}$, from the
Gaussian distribution was studied as a measure for the “distance” to chaos,
where
$R_{\alpha\boldsymbol{n}}^{\prime}={C_{\alpha\boldsymbol{n}}}/{\sqrt{\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle^{\prime}}}$
EFchaos-Benet03 ; EFchaos-WW ; foot-region . Here, in the computation of
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle^{\prime}$, in addition to an
average over perturbed states with energies close to $E_{\alpha}$, a further
average is taken over unperturbed states $|\boldsymbol{n}^{\prime}\rangle$
with unperturbed energies close to $E^{0}_{\boldsymbol{n}}$ by a small
quantity $\epsilon_{0}$. Specifically,

$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle^{\prime}=\sum_{\alpha^{\prime},\boldsymbol{n}^{\prime}}\frac{|C_{\alpha^{\prime}\boldsymbol{n}^{\prime}}|^{2}}{N_{E_{\alpha}}N_{E_{\boldsymbol{n}}}}\delta_{\epsilon}(E_{\alpha^{\prime}}-E_{\alpha})\delta_{\epsilon_{0}}(E_{\boldsymbol{n}^{\prime}}-E_{\boldsymbol{n}}),$
(33)

where $N_{E_{\alpha}}$ is defined in Eq.(15) and

$N_{E_{\boldsymbol{n}}}=\sum_{\boldsymbol{n}^{\prime}}\delta_{\epsilon_{0}}(E_{\boldsymbol{n}^{\prime}}-E_{\boldsymbol{n}}).$
(34)

It was found that, in some cases (not rare) in which the classical
counterparts undergo chaotic motion and the distributions $P(s)$ are quite
close to $P_{W}(s)$, the distributions of $R_{\alpha\boldsymbol{n}}^{\prime}$,
denoted by $g(R^{\prime})$, deviate notably from the Gaussian distribution.

In view of the semiclassical analysis given in the previous section, it is
understandable that deviation of $g(R^{\prime})$ from $f_{G}(R^{\prime})$ may
be larger than that of $f(R)$ from $f_{G}(R)$. In fact, unperturbed basis
states $|\boldsymbol{n}\rangle$ with close energies $E^{0}_{\boldsymbol{n}}$
may correspond to quite different values of $\boldsymbol{I}_{\boldsymbol{n}}$,
meanwhile, according to Eq.(24), the values of
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$ of those $\boldsymbol{n}$, for
which $\boldsymbol{I}_{\boldsymbol{n}}$ are far from each other, are usually
quite different. As a result, taking average over unperturbed basis states
with close $E^{0}_{\boldsymbol{n}}$ may drive the distribution of rescaled
components away from the Gaussian distribution.

Therefore, in order to obtain rescaled components that have a Gaussian
distribution, no average should be taken over the unperturbed energies
$E^{0}_{\boldsymbol{n}}$. We would remark that, when the torus of
$\boldsymbol{I}=\boldsymbol{I}_{\boldsymbol{n}}$ does not change rapidly with
$\boldsymbol{n}$, an average over a neighborhood of $\boldsymbol{n}$ is
allowed. We did not mention this averaging procedure in the above discussions,
because it is unnecessary in the derivation of Eq.(24).

##  III Numerical simulations in two models with classical counterparts

In order to test the above analytical results, numerical simulations have been
performed in two models possessing classical counterparts, the LMG model and
the Dicke model. In this section, we first briefly discuss the two models,
then, present numerical results about the distribution $f(R)$ and about the
suggested “distance” to chaos, namely, $\Delta_{EF}$, in comparison with other
“distances”.

Figure 3: The distribution $f(R)$ (open circles) and $g(R^{\prime})$ (solid
blocks with dashed lines) for $\lambda=1$ in the LMG model with (a)
$\Omega=80$ and (b) $\Omega=1000$. The solid curves indicate the Gaussian
distribution $f_{G}(R)$.

###  III.1 Models

Figure 4: Similar to Fig.3, but for the Dicke model with $\lambda=1$, (a)
$N=80$, (b)$N=1000$.

The first model we employ is a three-orbital LMG model LMG . This model is
composed of $\Omega$ particles, occupying three energy levels labeled by
$r=0,1,2$, each with $\Omega$-degeneracy. Here, we are interested in the
collective motion of this model, for which the dimension of the Hilbert space
is $\frac{1}{2}(\Omega+1)(\Omega+2)$. We use $\epsilon_{r}$ to denote the
energy of the $r$-th level and, for brevity, we set $\epsilon_{0}=0$. The
Hamiltonian of the model, in the form in Eq.(1), has pre-98

$\displaystyle H_{0}=\epsilon_{1}K_{11}+\epsilon_{2}K_{22},$ (35)
$\displaystyle V=\sum_{t=1}^{4}\mu_{t}V^{(t)},$ (36)

where

$\displaystyle V^{(1)}=K_{10}K_{10}+K_{01}K_{01},\
V^{(2)}=K_{20}K_{20}+K_{02}K_{02},$ $\displaystyle
V^{(3)}=K_{21}K_{20}+K_{02}K_{12},\ V^{(4)}=K_{12}K_{10}+K_{01}K_{21}.\ \ $
(37)

Here, the operators $K_{rs}$ are defined by

$K_{rs}=\sum_{\gamma=1}^{\Omega}a_{r\gamma}^{\dagger}a_{s\gamma},\quad
r,s=0,1,2,$ (38)

where $a^{\dagger}_{r\gamma}$ and $a_{r\gamma}$ are fermionic creation and
annihilation operators obeying the usual anti-commutation relations.

Figure 5: “Distance” to chaos in the LMG model and the Dicke model. The
measures $\Delta_{EF}$ (solid squares) in Eq.(30) and $\Delta^{\prime}_{EF}$
(open circles) in Eq.(59) are computed from statistical properties of EFs and
the measure $\Delta_{E}$ (open triangles) in Eq.(32) is computed from the
statistics of spectra. (a) The LMG model with $\Omega=80$, (b) The LMG model
with $\Omega=1000$, (c) the Dicke model with $N=80$, and (d) the Dicke model
with $N=1000$. The effective Planck constants are given by $1/\Omega$ and
$1/N$, respectively, in the two modes. The two measures $\Delta_{EF}$ and
$\Delta_{E}$ give almost the same results for the “distance” to chaos, when
the systems are not far from chaos.

For symmetric states, the operators $K_{rs}$ can be written in terms of
bosonic creation and annihilation operators $b^{\dagger}_{r}$ and $b_{r}$
chaos_Xu ,

$K_{rs}=b_{r}^{\dagger}b_{s},\quad
K_{r0}=K_{0r}^{\dagger}=b_{r}^{\dagger}\sqrt{\Omega-
b_{1}^{\dagger}b_{1}-b_{2}^{\dagger}b_{2}}$ (39)

for $r,s=1,2$. Under the transformation,

$b_{r}^{\dagger}=\sqrt{\frac{\Omega}{2}}(q_{r}-ip_{r}),\ \ \
b_{r}=\sqrt{\frac{\Omega}{2}}(q_{r}+ip_{r})$ (40)

for $r=1,2$, it is easy to verify that $q_{r}$ and $p_{s}$ obey the following
commutation relation,

$[q_{r},p_{s}]=\frac{i}{\Omega}\delta_{rs}.$ (41)

Hence, $1/\Omega$ plays the role of an effective Planck constant,

$\hbar_{\rm eff}=\frac{1}{\Omega}.$ (42)

It is straightforward to find that the classical counterpart of model has the
following Hamiltonian pre-98 ; Meredith98 ,

$H(\boldsymbol{p},\boldsymbol{q})=H_{0}(\boldsymbol{p},\boldsymbol{q})+\lambda
V(\boldsymbol{p},\boldsymbol{q}),$ (43)

where

$\displaystyle
H_{0}(\boldsymbol{p},\boldsymbol{q})=\frac{\epsilon^{\prime}_{1}}{2}(p_{1}^{2}+q_{1}^{2})+\frac{\epsilon^{\prime}_{2}}{2}(p_{2}^{2}+q_{2}^{2}),$
$\displaystyle
V(\boldsymbol{p},\boldsymbol{q})=\mu^{\prime}_{1}(q_{1}^{2}-p_{1}^{2})(1-G/2)+\mu^{\prime}_{2}(q_{2}^{2}-p_{2}^{2})(1-G/2)$
$\displaystyle+\frac{\mu^{\prime}_{3}}{\sqrt{2}}[(q_{2}^{2}-p_{2}^{2})q_{1}-2q_{2}p_{1}p_{2}]\sqrt{1-G/2}$
$\displaystyle+\frac{\mu^{\prime}_{4}}{\sqrt{2}}[(q_{1}^{2}-p_{1}^{2})q_{2}-2q_{1}p_{1}p_{2}]\sqrt{1-G/2},$
(44)

with $G=q_{1}^{2}+p_{1}^{2}+q_{2}^{2}+p_{2}^{2}\leq 2$. Here,
$\epsilon^{\prime}_{1}=\epsilon_{1}\Omega,\epsilon^{\prime}_{2}=\epsilon_{2}\Omega,\mu^{\prime}_{1}=\mu_{1}\Omega^{2},\mu^{\prime}_{2}=\mu_{2}\Omega^{2},\mu^{\prime}_{3}=\mu_{3}\Omega^{2}$,
and $\mu_{4}^{\prime}=\mu_{4}\Omega^{2}$. In our numerical simulations, we set
$\epsilon^{\prime}_{1}=44.1,\epsilon^{\prime}_{2}=64.5,\mu^{\prime}_{1}=62.1,\mu^{\prime}_{2}=70.2,\mu^{\prime}_{3}=76.5$,
and $\mu^{\prime}_{4}=65.7$. Under this choice of the parameters, for a fixed
value of $\lambda$, different values of $\Omega$ correspond to a same
classical counterpart.

The second model is a single-mode Dicke modelDicke ; Emary03 , which describes
the interaction between a single bosonic mode and a collection of $N$ two-
level atoms. The system can be described in terms of the collective operator
${\bf{J}}$ for the $N$ atoms, with

${J}_{z}=\sum_{i=1}^{N}{s}_{z}^{(i)},\ \
{J}_{\pm}=\sum_{i=1}^{N}{s}_{\pm}^{(i)},$ (45)

where ${s}_{\pm}^{(i)}={s}_{x}^{(i)}\pm i{s}_{y}^{(i)}$ and
${s}_{x(y,z)}^{(i)}$ are Pauli matrices divided by $2$ for the $i$-th atom.
The Dicke Hamiltonian is written as Emary03

$H=\omega_{0}J_{z}+\omega
a^{\dagger}a+\frac{\lambda}{\sqrt{N}}\mu(a^{\dagger}+a)(J_{+}+J_{-}),$ (46)

which can also be written in the form of $H=H_{0}+\lambda V$. The operator
${J}_{z}$ and ${J}_{\pm}$ obey the usual commutation rules for the angular
momentum,

$[J_{z},J_{\pm}]=\pm J_{\pm},\ \ [J_{+},J_{-}]=2J_{z}.$ (47)

Figure 6: “Distances” to classical and quantum chaos in the LMG model and the
Dicke model. The measure $\Delta_{EF}$ (solid squares) is the same as that in
Fig.5. The measure $\Delta_{cl}$ (open diamonds) is defined in Eq.(60) and was
computed from properties of the corresponding classical phase spaces. (a) The
LMG model with $\Omega=80$, (b) the LMG model with $\Omega=1000$, (c) the
Dicke model with $N=80$, and (d) the Dicke model with $N=1000$.

The Hilbert space of this model is spanned by vectors $|j,m\rangle$ with
$m=-j,-j+1,\cdots,j-1,j$, known as Dicke states. They are eigenstates of
$\boldsymbol{J}^{2}$ and $J_{z}$, with $J_{z}|j,m\rangle=m|j,m\rangle$ and
$\boldsymbol{J}^{2}|j,m\rangle=j(j+1)|j,m\rangle$. Below, we take $j$ as its
maximal value, namely, $j=N/2$; it is a constant of motion, since
$[\boldsymbol{J}^{2},H]=0$. Another conserved observable in the Dicke model is
the parity $\Pi$, given by $\Pi=\exp(i\pi\hat{N})$, where
$\hat{N}=a^{\dagger}a+J_{z}+j$ is an operator for the “excitation number”,
counting the total number of excitation quanta in the system. In our numerical
study, we consider the subspace with $\Pi=+1$.

Figure 7: Poincaré surfaces of section in the LMG model for $E=12+\sqrt{2}$
and $q_{1}=0$. (a) $\lambda=0.1$; (b) $\lambda=0.3$; (c) $\lambda=0.4$; (d)
$\lambda=1.0$.

Making use of the Holstein-Primakoff representation of the angular momentum
operators,

$\displaystyle J_{+}=b^{\dagger}\sqrt{2j-b^{\dagger}b},\ \ \
J_{-}=\sqrt{2j-b^{\dagger}b}\ b,$ $\displaystyle J_{z}=(b^{\dagger}b-j),$ (48)

where $b$ and $b^{\dagger}$ are bosonic operators satisfying
$[b,b^{\dagger}]=1$, the Hamiltonian can be further written in the following
form,

$\displaystyle H=\omega_{0}(b^{\dagger}b-j)+\omega a^{\dagger}a$
$\displaystyle+\lambda\mu(a^{\dagger}+a)\left(b^{\dagger}\sqrt{1-\frac{b^{\dagger}b}{2j}}+\sqrt{1-\frac{b^{\dagger}b}{2j}}b\right).$
(49)

We write Fock states related to $a^{\dagger}$ and $b^{\dagger}$ as
$|n_{a}\rangle$ and $|n_{b}\rangle$, respectively, for which

$a^{\dagger}a|n_{a}\rangle=n_{a}|n_{a}\rangle,\ \ \
b^{\dagger}b|n_{b}\rangle=n_{b}|n_{b}\rangle.$ (50)

According to Eq.(48), $n_{b}$ should be truncated at $(n_{b})_{\rm max}=N$. In
numerical simulations, we set $(n_{a})_{\rm max}=N$. Other parameters are
$\omega_{0}=\omega=1/N$ and $\mu=1/N$. Under the transformation

$\displaystyle\begin{cases}b^{\dagger}=\sqrt{\frac{N}{2}}(q_{1}-ip_{1}),\ \ \
b=\sqrt{\frac{N}{2}}(q_{1}+ip_{1}),\\\
a^{\dagger}=\sqrt{\frac{N}{2}}(q_{2}-ip_{2}),\ \ \
a=\sqrt{\frac{N}{2}}(q_{2}+ip_{2}),\end{cases}$ (51)

one finds that

$[q_{r},p_{s}]=\frac{i}{N}\delta_{rs}$ (52)

for $r=1,2$, and, hence, gets an effective Planck constant

$\hbar_{\rm eff}=\frac{1}{N}.$ (53)

The Hamiltonian of the classical counterpart of the model is written as

$H(\boldsymbol{p},\boldsymbol{q})=H_{0}(\boldsymbol{p},\boldsymbol{q})+\lambda
V(\boldsymbol{p},\boldsymbol{q}),$ (54)

where

$\displaystyle
H_{0}(\boldsymbol{p},\boldsymbol{q})=\frac{1}{2}(q_{1}^{2}+p_{1}^{2}+q_{2}^{2}+p_{2}^{2}-1),$
$\displaystyle
V(\boldsymbol{p},\boldsymbol{q})=2q_{1}q_{2}\sqrt{1-\frac{q_{1}^{2}+p_{1}^{2}}{2}}.$
(55)

Figure 8: Poincaré surfaces of section in the Dicke model for
$E=\frac{\sqrt{3}-1}{2}$ and $q_{1}=0$. (a) $\lambda=0.2$; (b) $\lambda=0.4$;
(c) $\lambda=0.6$; (d) $\lambda=1.0$.

###  III.2 Numerical results

In this section, we discuss results of numerical simulations performed in the
LMG model and the Dicke model. We first test validity of the semiclassical
result given in Eq.(24). As seen in Fig.1 and Fig.2, the average shape
$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$ and its classical analog
$\Pi(E_{\alpha},\boldsymbol{I_{n}})$ indeed show similar features in these two
models. We have computed the difference between the two shapes, given by

$d_{c}=\sum_{\boldsymbol{n}}|\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle-\Pi_{N}(E_{\alpha},\boldsymbol{I_{n}})|,$
(56)

where $\Pi_{N}(E_{\alpha},\boldsymbol{I_{n}})$ is the normalised
$\Pi(E_{\alpha},\boldsymbol{I_{n}})$. We found that $d_{c}=0.065$ in the LMG
model and $d_{c}=0.08$ in the Dicke model.

Then, we discuss properties of the distribution $f(R)$ for components
$R_{\alpha\boldsymbol{n}}$ in classically allowed regions with
$\Pi(E_{\alpha},\boldsymbol{I_{n}})\neq 0$. We found that this distribution is
indeed quite close to the Gaussian form $f_{G}(R)$, when the underlying
classical dynamics is chaotic, as illustrated in Fig.3 and Fig.4 with
$\lambda=1$. In the computation of $f(R)$, $100$ EFs in the middle energy
region were used. The energy windows $\epsilon$ are as follows: In the LMG
mode, $\epsilon\approx 3$ for $\Omega=80$ and $\epsilon\approx 0.02$ for
$\Omega=1000$, in contrast to the total energy domain $\Delta E=64.5$ in the
unperturbed system; in the Dicke mode, $\epsilon\approx 0.2$ for $N=80$ and
$\epsilon\approx 0.002$ for $N=1000$, in contrast to the total energy domain
$\Delta E=2$.

Figure 9: The distribution $f(R)$ (open circles) and $g(R^{\prime})$ (solid
blocks with dashed lines) for $\lambda=0.5$ in the defect Ising model and the
defect XXZ model. The solid curve indicates the Gaussian distribution
$f_{G}(R)$. (a) The defect Ising model with $N=10$; (b) the defect Ising model
with $N=15$; (c) the defect XXZ model with $N=12$ and $S_{z}=-1$; (d) the
defect XXZ model with $N=19$ and $S_{z}=-3.5$.

For comparison, we have also computed the distribution $g(R^{\prime})$ given
by another rescaling procedure, in which an average over unperturbed states is
also performed (see discussions given in Sec.II.2). In this rescaling
procedure, as discussed in Refs.EFchaos-Benet03 ; EFchaos-Benet00 , the
average shape of EFs is expected to have the following semiclassical
approximation,

$\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle^{\prime}\simeq\frac{S(E,E_{0})}{(2\pi\hbar)^{f}\rho(E)\rho_{0}(E_{0})},$
(57)

where $\rho_{0}(E_{0})$ and $\rho(E)$ are the density of states of the two
systems $H_{0}$ and $H$, respectively, and $S(E,E_{0})$ indicates the overlap
of the perturbed energy surface of $H=E$ and the unperturbed energy surface of
$H_{0}=E_{0}$,

$S(E,E_{0})=\int
d\boldsymbol{q}d\boldsymbol{p}\delta(E-H(\boldsymbol{p},\boldsymbol{q}))\delta(E_{0}-H_{0}(\boldsymbol{p},\boldsymbol{q})).$
(58)

The difference between $\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle^{\prime}$
in Eq.(57) and $\langle|C_{\alpha\boldsymbol{n}}|^{2}\rangle$ in Eq.(24) is
quite clear.

In the computation of $g(R^{\prime})$, only those rescaled components
$R_{\alpha\boldsymbol{n}}^{\prime}$ in the region with nonzero
$S(E_{\alpha},E_{\boldsymbol{n}})$ were used. We found that, unlike the case
of $f(R)$ discussed above, the distribution $g(R^{\prime})$ usually shows
obvious deviation from $f_{G}(R^{\prime})$ when the classical system is in the
chaotic regime (Fig.3 and Fig.4). Here, in the additional average for the
unperturbed system, $100$ EFs in the middle energy region were used, with
$\epsilon_{0}\approx 0.645$ in the LMG model and $\epsilon_{0}\approx 0.02$ in
the Dicke model.

Variation of the measure $\Delta_{EF}$ in Eq.(30) with the controlling
parameter $\lambda$ is given in Fig.5, together with the often-used measure
given by $\Delta_{E}$ of the statistics of spectra. In order to improve the
statistics, for each value of $\lambda$ , we used data obtained from several
values of $\lambda^{\prime}$ in a neighborhood of $\lambda$,
$\lambda^{\prime}\in[\lambda-0.05,\lambda+0.05]$.

The agreement of the two measures $\Delta_{EF}$ and $\Delta_{E}$ is already
good in the case of not quite large $\Omega$ in the LMG model [Fig.5(a) with
$\Omega=80$]. The agreement becomes better, when the value of $\Omega$ is
increased such that the system becomes closer to its classical limit
[Fig.5(b)]. Similar results were also found in the Dicke model [Fig.5(c) and
(d)]. Therefore, in these two models, the difference $\Delta_{EF}$ can be
regarded as a good measure for the “distance” to chaos.

For comparison, we have also computed the difference $\Delta^{\prime}_{EF}$
given by the other rescaling procedure,

$\Delta_{EF}^{\prime}=\int|I_{g}(R^{\prime})-I_{f_{G}}(R^{\prime})|dR^{\prime},$
(59)

where $I_{g}(R^{\prime})$ denotes the cumulative distribution for
$g(R^{\prime})$. Due to the obvious difference between the distribution
$g(R^{\prime})$ and the Gaussian distribution shown in Fig.3 and Fig.4, one
expects a notable difference between $\Delta_{EF}^{\prime}$ and $\Delta_{E}$.
Indeed, as shown in Fig.5, unlike the case with $\Delta_{EF}$ discussed above,
the agreement between $\Delta_{EF}^{\prime}$ and $\Delta_{E}$ is not good.

We have also computed a “distance” to chaos in the classical counterparts,
denoted by $\Delta_{cl}$, which measures the proportion of regular region in
energy surface. The measure is defined by

$\Delta_{cl}=\lim_{N_{T}\to\infty}\frac{N_{R}}{N_{T}},$ (60)

where $N_{T}$ is a total number of points taken randomly in an energy surface
of interest and $N_{R}$ is the number of the points for which
$\lambda_{L}<\lambda_{m}$. Here, $\lambda_{m}$ is some small quantity and
$\lambda_{L}$ is the Lyapunov exponent, defined as follows in the long time
limit,

$\lambda_{L}=\lim_{t\rightarrow\infty}\lim_{d_{0}\rightarrow
0}\frac{1}{t}\ln\frac{|d_{t}|}{|d_{0}|},$ (61)

where $d_{0}$ denotes the initial phase-space distance and $d_{t}$ denotes the
distance at a time $t$. In our numerical simulation, we took $t=1000$,
$N_{T}=5000$, $\lambda_{m}=0.02$ in the LMG model and $t=50000$, $N_{T}=5000$,
$\lambda_{m}=0.001$ in the Dicke model. In Fig.6, it is seen that the
agreement between the distances to quantum and classical chaos, characterized
by $\Delta_{EF}$ and $\Delta_{cl}$, respectively, is quite good. Some examples
of Poincaré surfaces of sections in the two models are shown in Fig.7 and
Fig.8.

##  IV Numerical simulations in models without classical counterparts

In this section, we study the distribution of rescaled components of EFs in
models without any classical counterpart. It seems reasonable to expect that
the final result of Sec.II.1, that is, that the distribution of appropriately
rescaled components should have a Gaussian form, may be valid to some extent
in this type of models as well.

Figure 10: Similar to Fig.5, but for the defect Ising model with (a) $N=10$
and (b)$N=15$.

Here, a major problem met is the determination of the region of components
that should be taken into account. As discussed previously, in a system with a
classical counterpart, this region corresponds to the classically
energetically allowed region. For a system without any classical counterpart,
this is a highly nontrivial problem. In this paper, we do not intend to solve
this problem, but, to circumvent it by restricting ourselves to models, whose
EFs occupy almost the whole unperturbed energy region. In this type of models,
one can simply use all components of the EFs when computing $f(R)$.
Specifically, we study a defect XXZ model and a defect Ising model, adopting a
periodic boundary condition in numerical simulations.

The defect XXZ model d-XXZ  is a modified XXZ model, in which an external
magnetic field is applied on two sites of the $N$ spins. The unperturbed
Hamiltonian and the perturbation have the following forms,

$\displaystyle
H_{0}=\sum_{i=1}^{N}s_{x}^{i}s_{x}^{i+1}+s_{y}^{i}s_{y}^{i+1}+\mu_{z}\sum_{i=1}^{N}s_{z}^{i}s_{z}^{i+1}$
(62) $\displaystyle V=\mu_{1}s_{z}^{1}+\mu_{4}s_{z}^{4},$ (63)

where the periodic boundary condition implies that $s^{N+1}_{a}=s^{1}_{a}$ for
$a=x,y,z$. The system is a quantum chaotic system for $\lambda$ within an
appropriate regime, while, it exhibits the so-called many-body localisation
for $\lambda$ sufficiently large. The Hamiltonian $H$ is commutable with
$S_{z}$, the $z$-component of the total spin, and we consider a subspace with
a definite value of $S_{z}$ in our numerical study. Other parameters used in
this model are $\mu_{1}=1.11$, $\mu_{4}=1.61$, and $\mu_{z}=1$.

The defect Ising model is a transverse Ising model, in which an additional
magnetic field is applied on two sites of the $N$ spins, with

$\displaystyle
H_{0}=\sum_{i}^{N}s_{z}^{i}s_{z}^{i+1}+\mu_{x}\sum_{i=1}^{N}s_{x}^{i},$ (64)
$\displaystyle V=\mu_{1}s_{z}^{1}+\mu_{4}s_{z}^{4}.$ (65)

Similarly, it is a quantum chaotic system for $\lambda$ in an appropriate
regime and exhibits many-body localisation for $\lambda$ sufficiently large.
The parameters used are $\mu_{1}=1.11$, $\mu_{4}=1.61$, and $\mu_{x}=0.6$.

Figure 11: Similar to Fig.5, but for the defect XXZ model with (a) $N=12$,
$S_{z}=-1$, and (b)$N=19$, $S_{z}=-3.5$.

Our numerical simulations reveal that, in these two models, the distributions
$f(R)$ are also quite close to the Gaussian form $f_{G}(R)$, when the
statistics of the spectra is close to the prediction of RMT as illustrated in
Fig.9 with $\lambda=0.5$. Unlike the two models discussed in the previous
section, the distributions $g(R^{\prime})$ are also close to the Gaussian form
at $\lambda=0.5$.

In the computation of $f(R)$, $50$ EFs in the middle energy region were used.
The energy windows $\epsilon$ are as follows: In the defect Ising model,
$\epsilon\approx 0.2$ and $\epsilon_{0}=0.02$ for $N=10$ in contrast to the
total energy domain $\Delta E=7.11$, and $\epsilon\approx 0.07$ for $N=15$ in
contrast to $\Delta E=10.64$; in the defect XXZ model, $\epsilon\approx 0.3$
and $\epsilon_{0}=0.02$ for $N=12$, $S_{z}=-1$ in contrast to $\Delta E=8.03$,
and $\epsilon\approx 0.01$ for $N=19$, $S_{z}=-3.5$ in contrast to the total
energy domain $\Delta E=10.66$.

Figure 12: (a) The distributions of $f(R)$ (open circles) and $g(R^{\prime})$
(solid blocks with dashed lines) in the defect Ising model with $N=15$ and
$\lambda=0.06$. (b) The corresponding cumulative distribution of the nearest-
level-spacing distribution (open circles). The solid curve indicates the
cumulative distribution given by the Wigner surmise.

The two measures $\Delta_{EF}$ and $\Delta_{E}$ exhibit similar behaviors,
like the cases discussed in the previous section for the LMG and Dicke models
(Fig.10 and Fig.11). Thus, at least in these two models, the difference
$\Delta_{EF}$ can be regarded as a good measure for the “distance” to chaos.

In consistence with the behaviors of the distribution $g(R^{\prime})$
illustrated in Fig.9, the two quantities $\Delta^{\prime}_{EF}$ and
$\Delta_{EF}$ are close in most regions where the systems are chaotic systems
according to their spectra statistics. That is, in most cases, an average over
the unperturbed energy does not bring much difference in the defect Ising and
defect XXZ models. This may be partially related to the fact that EFs in these
two models occupy almost the whole energy region for $\lambda$ not small.

There are still some regions of $\lambda$ in Fig.10(b) and Fig.11(b) with
relatively large Hilbert spaces, in which $\Delta^{\prime}_{EF}$ shows some
notable deviation from $\Delta_{EF}$ and $\Delta_{E}$. Some examples of the
distributions $f(R)$ and $g(R^{\prime})$ in this case are shown in Fig.12 and
Fig.13, together with the corresponding distributions of $I_{P}(s)$ and
$I_{P_{W}}(s)$.

Figure 13: Similar to Fig.12, but for the defect XXZ model with $N=19$,
$S_{z}=-7$ and $\lambda=0.12$.

##  V Conclusions

In this paper, based on semiclassical analysis, it has been shown that those
components of EFs of quantum chaotic systems, which lie in classically-allowed
regions of integrable bases, can be regarded as random numbers in a sense
similar to that stated in the Berry’s conjecture. For the distribution $f(R)$
of these components to have a Gaussian form, which is predicted by the RMT, an
appropriated rescaling procedure with respect to the average shape of EFs is
needed, where the average should be taken over perturbed states with
neighbouring energies. It is found that an additional average over unperturbed
basis states with neighbouring unperturbed energies may cause deviation of the
distribution of rescaled components of EFs from the Gaussian form.

The above results suggest that deviation of the distribution $f(R)$ from the
Gaussian distribution may be used as a measure for the “distance” to quantum
chaos. In two models possessing classical counterparts, when the perturbed
system goes from integrable to chaotic with the increase of perturbation
strength, our numerical simulations show that this deviation coincides with
the deviation of the nearest-level-spacing distribution from the prediction of
RMT.

It is known that specific dynamics of the underlying classical systems may
induce certain modifications to the Berry’s conjecture KpHl ; Heller87 ;
Bies01 ; Sr96 ; Sr98 ; Backer02 ; Urb03 ; Kp05 . Since the main result of this
paper is based on this conjecture, specific underlying classical dynamics may
have some influence in results of this paper as well. In particular, it may
induce some deviation of the distribution $f(R)$ for some EFs from the
Gaussian distribution. However, if sufficiently many EFs are used in the
computation of $f(R)$, it is reasonable to expect that the induced deviation
should be small.

In two models without simple classical counterpart, we have found similar
numerical results about the distribution of $f(R)$. Analytical explanation of
this point is still lacking. It seems that the following feature of these two
models may be of relevance. That is, in both models the matrices of the
perturbations $V$ in the unperturbed bases do not have a clear band structure;
in other words, the perturbation couples basis vectors far separated in
energy. We hope that these numerical results may stimulate more
investigations.

###### Acknowledgements.

This work was partially supported by the Natural Science Foundation of China
under Grant Nos. 11275179, 11535011, and 11775210.

## References

  * (1) M. V. Berry, J. Phys. A: Math. Gen. 10, 2083 (1977); Proc. R. Soc. London, Ser. A 400, 229 (1985).
  * (2) M. V. Berry, in Les Houches LII, Chaos and Quantum Physics, edited by M. -J. Giannoni, A. Voros, and J. Zinn-Justin (North-Holland, Amsterdam, 1991).
  * (3) Quantum Chaos: Between Order and Disorder, edited by G. Casati and B.V. Chirikov (Cambridge University Press, Cambridge, England, 1994).
  * (4) F. Haake, Quantum Signatures of Chaos, 3rd ed. (Springer-Verlag, Berlin, 2010).
  * (5) D. C. Meredith, S. E. Koonin and M. R. Zirnbauer, Phys. Rev. A 37, 3499 (1988).
  * (6) V. V. Flambaum, G. F. Gribakin, and F. M. Izrailev, Phys. Rev. E 53, 5729 (1996); J. Phys. A: Math. Gen. 29, 5817 (1996).
  * (7) W.-G. Wang, F.M. Izrailev, and G. Casati, Phys. Rev. E 57, 323 (1998).
  * (8) H.-J. Stöckmann, Quantum Chaos: An Introduction (Cambridge University Press, Cambridge, England, 1999).
  * (9) Y. Alhassid, Rev. Mod. Phys. 72, 895, (2000).
  * (10) G. Hackenbroich, C. Viviescas, B. Elattari, and F. Haake, Phys.Rev.Lett. 86, 5262 (2001).
  * (11) W.-g. Wang, Phys. Rev.E 65, 036219 (2002).
  * (12) M. Horvat and T. Prosen, J. Phys. A 36, 4015 (2003).
  * (13) H. Olofsson, S. Aberg, O. Bohigas, and P. Leboeuf, Phys. Rev. Lett. 96, 042502 (2006).
  * (14) S. Gnutzmann, J. P. Keating, and F. Piotet, Ann. Phys. 325, 2595 (2010).
  * (15) A.M. Smith and L. Kaplan, Phys. Rev. E 80, 035205 (2009); ibid. 82, 016214 (2010).
  * (16) J.M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * (17) M.Srednicki, Phys. Rev. E 50, 888 (1994).
  * (18) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * (19) L. F. Santos, F. Borgonovi, F. M. Izrailev, Phys. Rev. Lett. 108, 094102 (2012); Phys. Rev. E 85, 036209 (2012).
  * (20) A. Altland, F. Haake, Phys. Rev. Lett. 108, 073601 (2012).
  * (21) M. Rigol, M. Srednicki, Phys. Rev. Lett. 108, 110601 (2012).
  * (22) A. Khodja, R. Steinigeweg, J. Gemmer, Phys. Rev. E 91, 012120 (2015).
  * (23) L. D’Alessio, Y. Kafri, A. Polkovnikov, M. Rigol, Advances in Physics 65(3), 239 (2016).
  * (24) F. Borgonovi, F. M. Izrailev, L. F. Santos, V. G. Zelevinsky, Phys. Rep. 626, 1-58 (2016)
  * (25) V. Buch, R. B. Gerber, and M. A. Ratner, J. Chem. Phys. 76, 5397 (1982).
  * (26) L. Benet, J. Flores, H. Hernandez-Saldana, F. M. Izrailev, F. Leyvraz and T. H. Seligman, J. Phys. A: Math. Gen. 36, 1289 (2003).
  * (27) J. Wang, W.-g. Wang, Chaos, Solitons & Fractals, 91, 291 (2016).
  * (28) A. Voros, Ann. Inst. Henri Poincaré Sect. A 24, 31 (1976).
  * (29) A. Voros, Ann. Inst. Henri Poincaré Sect. A 26, 343 (1977).
  * (30) F. Borgonovi, I. Guarneri, and F. M. Izrailev, Phys. Rev. E 57, 5291 (1998).
  * (31) L. Benet, F.M. Izrailev, T.H. Seligman, A. Suárez-Moreno, Phys. Lett. A 277, 87 (2000).
  * (32) H. J. Lipkin, N. Meshkov, and A. J. Glick, Nucl. Phys. 62, 188 (1965).
  * (33) R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * (34) C. Emary, T. Brandes, Phys. Rev. E 67, 066203 (2003).
  * (35) E. J. Torres-Herrera, L. F. Santos, Phys. Rev. E 89, 062110 (2014).
  * (36) The distribution of rescaled components computed in Ref.EFchaos-WW  is for the so-called nonperturbative regions of EFs. The relation between these regions and classically-allowed regions will be discussed in a different paper.
  * (37) Xu Gong-ou, Gong Jiang-bin, Wang Wen-ge, Yang Ya-tian, and Fu De-ji, Phys. Rev. E 51, 1770 (1995).
  * (38) E. J. Heller, Phys. Rev. Lett. 53, 1515 (1984).
  * (39) P. O’Connor, J. Gehlen and E. J. Heller, Phys. Rev. Lett. 58, 1296 (1987).
  * (40) M. Srednicki, Phys. Rev. E 54, 954 (1996).
  * (41) S. Hortikar and M. Srednicki, Phys. Rev. E 57, 7313 (1998).
  * (42) W. E. Bies, L. Kaplan, M.R. Haggerty, E.J. Heller, Phys. Rev. E 63, 066214 (2001).
  * (43) A. Bäcker and R. Schubert, J. Phys. A: Math. Gen. 35, 539 (2002).
  * (44) J. D. Urbina and K. Richter, J. Phys. A: Math. Gen. 36, L495 (2003); Phys. Rev. E 70, 015201 (2004); Phys. Rev. Lett. 97, 214101 (2006).
  * (45) L. Kaplan, Phys. Rev. E 71, 056212 (2005).


# Linear Nearest Neighbor Flocks with All Distinct Agents

R. Lyons , J. J. P. Veerman Digimarc Corporation, 9405 SW Gemini Drive,
Beaverton, OR, USA 97008-7192; Fariborz Maseeh Dept. of Math. and Stat.,
Portland State Univ.; e-mail<EMAIL_ADDRESS>Maseeh Dept. of Math. and
Stat., Portland State Univ.; e-mail<EMAIL_ADDRESS>
###### Abstract

This paper analyzes the global dynamics of 1-dimensional agent arrays with
nearest neighbor linear couplings. The equations of motion are second order
linear ODE’s with constant coefficients. The novel part of this research is
that the couplings are different for each distinct agent. We allow the forces
to depend on the positions and velocity (damping terms) but the magnitudes of
both the position and velocity couplings are different for each agent. We,
also, do not assume that the forces are “Newtonian” (i.e. the force due to A
on B equals the minus the force of B on A) as this assumption does not apply
to certain situations, such as traffic modeling. For example, driver A
reacting to driver B does not imply the opposite reaction in driver B.

There are no known analytical means to solve these systems, even though they
are linear, and so relatively little is known about them. This paper is a
generalization of previous work that computed the global dynamics of
1-dimensional sequences of identical agents [3] assuming periodic boundary
conditions. In this paper, we push that method further, similar to [2], and
use an extended periodic boundary condition to to gain quantitative insights
to the systems under consideration. We find that we can approximate the global
dynamics of such a system by carefully analyzing the low-frequency behavior of
the system with (generalized) periodic boundary conditions.

##  1 Introduction

A one-dimensional lattice of coupled agents is a model for many physical
systems. If all the agents are identical then connecting nearest neighbors
with a Hooke’s Law is a simple model of one-dimensional crystals [1]. If one
assumes periodic boundary conditions, then the eigenvectors of the system are
the Discrete Fourier Transform basis functions. In the 1950’s this nearest
neighbor crystal model was extended to include agents of different mass [6].

In the 1950’s simplified “microscopic” traffic models appeared with agents
coupled with a force dependent on spatial differences and an added force that
is a function of the difference of agent velocities [5] (see [9] for a survey
of traffic models). The velocity dependent force term originated as an
empirical law, observed when an automobile attempts to follow a leader. Since
then the subject of cooperative control has matured considerably [12] [13].

Recent technological advances make automated traffic platoons possible so
there has been renewed interest in one-dimensional lattice dynamics. There are
several works on both single [10] [11] and double integrator systems [14] [7]
[10]. The results for both single and double integrator systems with nearest
neighbor interactions are summarized in [11].

In the one-dimensional traffic platoon, one would like to know whether is it
possible, or even reasonable, to have a long platoon consisting of $N$ agents?
If we form a caravan of trucks, do we need to break the caravan into separate
small chunks or can we form a single caravan of, perhaps, over a thousand
trucks. There is a large literature on the topic, but almost all the
literature we are aware of addresses this question in the unrealistic case
where all cars are identical or distributed in some other highly improbable
way. Each agent is distinct and may have a unique mass. We can force the
forward and backward couplings to have a specific ratio but it is difficult
and certainly impracticable to insist that the force magnitudes are identical
for all agents. In section 3 we assign specific ratios for forward and
backward couplings as in equations (13) and (14), but the weights
$g_{x}^{(\alpha)},g_{v}^{(\alpha)}$ are chosen randomly from a distribution.

More specifically, in this paper we shall analyze a one-dimensional lattice of
agents with linear nearest neighbor couplings determined by the distance and
the velocity difference between neighbors. We will assume equations of motion
where the force on an agent is linear in position and velocity differences
(double integrator system). We shall not assume that the forces are
“Newtonian” (i.e. the force due to A on B equals the minus the force of B on
A).

In previous work ([3], [4]) it was shown that if the agents are identical and
the system has periodic boundary conditions then the equations of motion are
solvable. This system has equations of motion given by the ODE,

$\frac{d}{dt}\begin{pmatrix}Z\\\
\dot{Z}\end{pmatrix}=\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\
L_{x}&L_{v}\end{pmatrix}\begin{pmatrix}Z\\\ \dot{Z}\end{pmatrix}\;,$ (1)

where $Z$ is a vector positions and $L_{x},L_{v}$ are row-sum zero circulant
tri-diagonal matrices. Since the matrices are circulant, $L_{x}$ and $L_{v}$
commute. This is instrumental in the finding solutions and the conditions of
stability. If the forces are extended to include next-nearest neighbor terms
then, assuming periodic boundary conditions, the equations are motion are,
also, given by equation (1). As in the nearest neighbor case, $L_{x}$ and
$L_{v}$ are row-sum zero circulant matrices but, this time, they have 5 non-
zero diagonals. The solutions are more complicated as are the conditions of
stability, [8]. For both systems the matrices $L_{x}$ and $L_{v}$ commute. In
both these cases the characteristic polynomial has a double root at $0$, which
corresponds to the stable configuration where all agents are moving at a
constant velocity. To find the asymptotic behavior of the system you can
expand the zero locus around this point. On stable systems, roots near the
origin dominate the long-term behavior of the system as other roots have
larger negative real components and so decay faster. Expanding the
characteristic equation near the origin yields an approximation to the signal
velocity and a dispersion term.

However, in this work we do not assume the agents are identical. Instead we
introduce a repeating sequence of $p$ distinct agents. Duplicate this string
of $p$ agents and use an extension of periodic boundary conditions which was
first described in [2]. In particular, let $A_{0},\cdots,A_{p-1}$ be $p$ agent
types organized in a one-dimensional lattice,

$A_{p-1}\leftrightarrow A_{p-2}\leftrightarrow\cdots\leftrightarrow
A_{1}\leftrightarrow A_{0}.$

Then repeat this $p-$sequence $q$ times to get a total of $N=pq$ agents. The
general form for this system the matrices $L_{x}$ and $L_{v}$ do not commute.
The case for $p=2$ and $p=3$ is analyzed in [2] for both nearest neighbor and
next nearest neighbor interactions. Since $L_{x}$ and $L_{v}$ do not commute
the system is considerably more difficult to analyze, but some conditions
necessary for stability are derived.

In this paper we present a variety of tools to analyze this general system.
The goal is to first understand the periodic case and then use these results
to shed light on the general system of $N$ agents traveling on the real line.
In Section 2 Theorem 2.7 we prove a generalization of the stability condition
in [2]. To pursue the dynamics of a general system we start with the system in
equation (1), where $L_{x}$ and $L_{v}$ are the circulant matrices with $-1$
on the diagonal and $1/2$ on the sub and super diagonals. This is the system
in [4] with

$\rho_{x,1}=\rho_{x,-1}=-1/2\text{ and }\rho_{v,1}=\rho_{v,-1}=-1/2.$

This system is stable, [4]. We extend this by scaling each row by a distinct
value, which is the same as taking distinct weights $g_{x}^{(\alpha)}$ and
$g_{v}^{(\alpha)}$. In this case equation (1) becomes,

$\frac{d}{dt}\begin{pmatrix}Z\\\
\dot{Z}\end{pmatrix}=\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\
G_{x}L_{x}&G_{v}L_{v}\end{pmatrix}\begin{pmatrix}Z\\\ \dot{Z}\end{pmatrix}\;,$
(2)

where $G_{x}$ and $G_{v}$ are diagonal matrices with positive real values.
Again, we note that $G_{x}L_{x}$ and $G_{v}L_{v}$ do not commute, so this
system is more difficult to analyze. In Section 3 we prove stability for a
special case of this general problem. In Section 5 we analyze the dynamics by
expanding the characteristic polynomial root locus around the double root at
$0$. As in the problems above the asymptotic behavior of the system is given
by roots near the double root at $0$. Expansion around this point results in
an expression of the signal velocity and a dispersion term, given in Theorem
(5.1).

The expansion, used in Section 5, requires an extension of the periodic
boundary condition first found in [3] and [2]. We take the $p$ distinct agents
and repeat them $q$ times. This guarantees that the locus is well approximated
by a continuous curve as $q$ gets large, and this allows us to use a Taylor
expansion. In the simulations in section 6 we will set $q=1$ and show that the
results apply well to the general case of $N=p$ distinct agents.

##  2 Linear Nearest Neighbor Systems

In this section we describe the equations of motion for one-dimensional
lattice with linear equations and nearest neighbor interactions. In Theorem
2.7 below, we prove a necessary condition for stability. Later, in Section 3
we restrict to the case where the Laplacians are symmetric matrices that have
rows scaled by independent positive weights. In this case we can derive some
properties of the systems dynamics.

We consider $N$ agents consisting of $q$ groups of $p$ agents. The first
cluster of $p$ are unique and they are followed by $q$ identical clusters of
$p$ so that the total number of agents is $N=pq$. We follow the conventions in
[2], except we change the signs so that $g_{x},g_{v}\geq 0$. After adjusting
for the spacing, by subtracting $k\Delta$ from the $k$’th agent (see [4]), we
can write the nearest neighbor coupling as,

$\frac{d^{2}z_{k}^{(\alpha)}}{dt^{2}}=-g_{x}^{(\alpha)}\left(z_{k}^{(\alpha)}+\rho_{x,1}^{(\alpha)}z_{k+1}^{(\alpha+1)}+\rho_{x,-1}^{(\alpha)}z_{k-1}^{(\alpha-1)}\right)-g_{v}^{(\alpha)}\left(\dot{z}_{k}^{(\alpha)}+\rho_{v,1}^{(\alpha)}\dot{z}_{k+1}^{(\alpha+1)}+\rho_{v,-1}^{(\alpha)}\dot{z}_{k-1}^{(\alpha-1)}\right)$
(3)

where the arithmetic in $\alpha$ is $\mod(p)$. The coefficients satisfy
$\rho_{x,-1}^{(\alpha)}+1+\rho_{x,1}^{(\alpha)}=0$ and
$\rho_{v,-1}^{(\alpha)}+1+\rho_{v,1}^{(\alpha)}=0$ for all
$\alpha=0,\cdots,p-1$.

At this point we assume periodic boundary conditions. and re-arrange vector
components so the first $q$ coordinates are agents of type $\alpha=0$, the
next $q$ agents are type $\alpha=1$, etc. This is a generalization of
technique in [2].

To write the system in matrix form we start by writing an $N\times N$ matrix,

$L_{\eta}=\begin{pmatrix}-\mathbf{I}&-\rho_{\eta,1}^{(0)}\mathbf{I}&\mathbf{0}&\cdots&-\rho_{\eta,-1}^{(0)}P_{-}\\\
-\rho_{\eta,-1}^{(1)}\mathbf{I}&-\mathbf{I}&-\rho_{\eta,1}^{(1)}\mathbf{I}&\cdots&\mathbf{0}\\\
\mathbf{0}&-\rho_{\eta,-1}^{(2)}\mathbf{I}&-\mathbf{I}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
-\rho_{\eta,-1}^{(p-1)}P_{+}&\mathbf{0}&\mathbf{0}&\cdots&-\mathbf{I}\end{pmatrix}\text{
for }\eta=x,v,$ (4)

where $\mathbf{I}$ and $\mathbf{0}$ are $q\times q$ matrices and the matrix
$P_{+}$ and its inverse $P_{-}$ are $q\times q$ cyclic permutations matrices,

$\displaystyle P_{+}=\begin{pmatrix}0&1&0&\cdots&0\\\ 0&0&1&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&0\\\ 0&\ddots&\ddots&0&1\\\
1&0&\cdots&0&0\end{pmatrix}~{},~{}~{}~{}P_{-}=\begin{pmatrix}0&0&\cdots&0&1\\\
1&0&0&\ddots&0\\\ 0&1&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&0&0\\\
0&\cdots&0&1&0\end{pmatrix}.$ (5)

We write the scales $g_{x}^{(\alpha)},g_{v}^{(\alpha)}$ as an $N\times N$
matrix,

$G_{\eta}=\begin{pmatrix}g_{\eta}^{(0)}\mathbf{I}&\mathbf{0}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&g_{\eta}^{(1)}\mathbf{I}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&g_{\eta}^{(2)}\mathbf{I}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\cdots&g_{\eta}^{(p-1)}\mathbf{I}\end{pmatrix}\text{
for }\eta=x,v.$ (6)

The equations are motion are,

$\begin{pmatrix}\mathbf{\dot{z}}^{(0)}\\\ \mathbf{\dot{z}}^{(1)}\\\ \vdots\\\
\mathbf{\dot{z}}^{(p-1)}\\\ {\mathbf{\ddot{z}}}^{(0)}\\\
{\mathbf{\ddot{z}}}^{(1)}\\\ \vdots\\\
{\mathbf{\ddot{z}}}^{(p-1)}\end{pmatrix}=\left(\begin{array}[]{ccc|ccc}&&&&&\\\
&\mathbf{0}&&&\mathbf{I}&\\\ &&&&&\\\ \hline\cr&&&&&\\\
&G_{x}L_{x}&&&G_{v}L_{v}&\\\
&&&&&\end{array}\right)\begin{pmatrix}\mathbf{z}^{(0)}\\\ \mathbf{z}^{(1)}\\\
\vdots\\\ \mathbf{z}^{(p-1)}\\\ \dot{\mathbf{z}}^{(0)}\\\
\dot{\mathbf{z}}^{(1)}\\\ \vdots\\\ \dot{\mathbf{z}}^{(p-1)}\end{pmatrix}\;.$
(7)

The matrices $\mathbf{P}_{+}$ and $\mathbf{P}_{-}$ are $q\times q$ circulant
matrices and all circulant matrices have a common set of $q$ eigenvectors
given by,

$\mathbf{v}_{m}=\left[1,\exp\left(\frac{2\pi
i}{q}m\right),\exp\left(\frac{2\pi i}{q}2m\right),\cdots,\exp\left(\frac{2\pi
i}{q}(q-1)m\right)\right]^{T}$

where $m=0,1,\cdots,q-1$.

To simplify the notation, we introduce the following 1st degree polynomials,

$\displaystyle\psi_{0}^{(\alpha)}(\nu)$
$\displaystyle=\left(g_{v}^{(\alpha)}\nu+g_{x}^{(\alpha)}\right)$
$\displaystyle\psi_{1}^{(\alpha)}(\nu)$
$\displaystyle=g_{v}^{(\alpha)}\rho_{v,1}^{(\alpha)}\nu+g_{x}^{(\alpha)}\rho_{x,1}^{(\alpha)}$
$\displaystyle\psi_{-1}^{(\alpha)}(\nu)$
$\displaystyle=g_{v}^{(\alpha)}\rho_{v,-1}^{(\alpha)}\nu+g_{x}^{(\alpha)}\rho_{x,-1}^{(\alpha)}$
(8)

For all $\alpha,\nu$, we have,

$\psi_{-1}^{(\alpha)}(\nu)+\psi_{0}^{(\alpha)}(\nu)+\psi_{1}^{(\alpha)}(\nu)=0.$
(9)

The eigenvalue equation of the general system is simplified by using the
following,

######  Proposition 2.1.

The eigenvectors of the matrix in equation (10) are given by the $2N=2pq$
vectors,

$\mathbf{u}_{\nu}(m,\phi)=\begin{pmatrix}\epsilon_{0}\mathbf{v}_{m}\\\
\epsilon_{1}\mathbf{v}_{m}\\\ \vdots\\\ \epsilon_{p-1}\mathbf{v}_{m}\\\
\nu\epsilon_{0}\mathbf{v}_{m}\\\ \nu\epsilon_{1}\mathbf{v}_{m}\\\ \vdots\\\
\nu\epsilon_{p-1}\mathbf{v}_{m}\end{pmatrix}=\begin{pmatrix}\mathbf{v}_{m}\otimes\mathbf{\epsilon}\\\
\nu\left(\mathbf{v}_{m}\otimes\mathbf{\epsilon}\right)\end{pmatrix}\;.$ (10)

which have eigenvalue $\nu$ and where

$\mathbf{\epsilon}=\begin{pmatrix}\epsilon_{0}&\epsilon_{2}&\cdots&\epsilon_{p-1}\end{pmatrix}^{T}.$

For each $m\in\\{0,\cdots q-1\\}$ there are $2p$ eigenvalues given by the
determinant of the matrix in equation (11). For each of these roots the
eigenvector uses the values $\epsilon_{j}$ that satisfy,

$\begin{pmatrix}\nu^{2}+\psi_{0}^{(0)}(\nu)&\psi_{1}^{(0)}(\nu)&0&\cdots&\psi_{-1}^{(0)}(\nu)e^{-i\phi}\\\
\psi_{-1}^{(1)}(\nu)&\nu^{2}+\psi_{0}^{(1)}(\nu)&\psi_{1}^{(1)}(\nu)&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}^{(p-1)}(\nu)e^{i\phi}&0&0&\cdots&\nu^{2}+\psi_{0}^{(p-1)}(\nu)\end{pmatrix}\begin{pmatrix}\epsilon_{0}\\\
\epsilon_{1}\\\ \vdots\\\ \epsilon_{p-1}\end{pmatrix}=\begin{pmatrix}0\\\ 0\\\
\vdots\\\ 0\end{pmatrix}\;.$ (11)

###### Proof.

The proof of this proposition is a generalization of the proof of proposition
2.1 from [2]. Apply the matrix in equation (7) to the vector in equation (10).
Use,

$\begin{array}[]{cc}\mathbf{P}_{-}\mathbf{v}_{m}=e^{-i\phi}\mathbf{v}_{m}\;,&\mathbf{P}_{+}\mathbf{v}_{m}=e^{i\phi}\mathbf{v}_{m}.\end{array}$

where we define $\phi=\frac{2\pi}{q}m$. The top $N$ coordinates follow
immediately and the bottom $N$ coordinates yield $q$ copies of the equation
(11). ∎

######  Corollary 2.2.

The eigenvalues of the system are the roots to the polynomial,

$P_{\phi}(\nu)=\det\begin{pmatrix}\nu^{2}+\psi_{0}^{(0)}(\nu)&\psi_{1}^{(0)}(\nu)&0&\cdots&\psi_{-1}^{(0)}(\nu)e^{-i\phi}\\\
\psi_{-1}^{(1)}(\nu)&\nu^{2}+\psi_{0}^{(1)}(\nu)&\psi_{1}^{(1)}(\nu)&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}^{(p-1)}(\nu)e^{i\phi}&0&0&\cdots&\nu^{2}+\psi_{0}^{(p-1)}(\nu)\end{pmatrix}.$
(12)

For each value of $\phi=\frac{2\pi i}{q}m$, $m=0,\cdots q-1$, there are $2p$
roots.

######  Proposition 2.3.

When $\phi=0$ (e.g. $m=0$) the constant and linear terms for the polynomial
$P_{0}(\nu)$ both vanish.

###### Proof.

Neither the linear nor constant terms of the polynomial cannot have $\nu^{2}$
as a factor. Set $\phi=0$ and remove the terms with $\nu^{2}$ as a factor. The
resulting polynomial is,

$\det\begin{pmatrix}\psi_{0}^{(0)}(\nu)&\psi_{1}^{(0)}(\nu)&0&\cdots&\psi_{-1}^{(0)}(\nu)\\\
\psi_{-1}^{(1)}(\nu)&\psi_{0}^{(1)}(\nu)&\psi_{1}^{(1)}(\nu)&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}^{(p-1)}(\nu)&0&0&\cdots&\psi_{0}^{(p-1)}(\nu)\end{pmatrix}$

Every row in this matrix sums to zero, by equation (9). This means the vector,
consisting of all $1$’s is an eigenvector with eigenvalue $0$ and so the
determinant vanishes for all $\nu$. The constant and linear terms are
contained in this reduced polynomial and so must vanish. ∎

######  Remark 2.4.

The proof of Proposition 2.3 says the second order term of the polynomial
$P_{0}(\nu)$ must contain exactly one of the diagonal $\nu^{2}$ terms. The
third order term of $P_{0}(\nu)$ must, also, contain exactly one of $\nu^{2}$
terms. These facts will be used below.

######  Proposition 2.5.

$P_{\phi}(\nu)=s(\nu)+(-1)^{p}r_{\phi}(\nu),$

where all the $\phi$ dependence is in the polynomial,

$r_{\phi}(\nu)=(1-e^{i\phi})\left(\psi_{1}^{(0)}(\nu)\psi_{1}^{(1)}(\nu)\cdots\psi_{1}^{(p-1)}(\nu)\right)+(1-e^{-i\phi})\left(\psi_{-1}^{(0)}(\nu)\psi_{-1}^{(1)}(\nu)\cdots\psi_{-1}^{(p-1)}(\nu)\right)$

and $s(\nu)$ has zero constant and linear terms.

###### Proof.

The only terms of the expansion of equation (12) that depend on $\phi$ contain
either $\psi_{-1}^{(0)}(\nu)e^{-i\phi}$ or $\psi_{1}^{(p-1)}(\nu)e^{i\phi}$
but not both. The determinant is a sum over permutations $\sigma$ of terms
$\operatorname{sgn}(\sigma)M_{0\sigma 0}\cdots M_{p-1\sigma p-1}$. The only
non-zero permutations have $\sigma(k)\in\\{k-1,k,k+1\\}MOD(p)$. The
permutations that contain the term $\psi_{-1}^{(0)}(\nu)e^{-i\phi}$ but not
$\psi_{1}^{(p-1)}(\nu)e^{i\phi}$ must have $\sigma(0)=p-1$. The matrix is tri-
diagonal so the value $\sigma(p-1)\in\\{0,p-1,p-2\\}$ for all $\sigma$. But,
in this case we know $\sigma(p-1)\neq p-1$ and, by assumption,
$\sigma(p-1)\neq 0$ (or the term would contain the $e^{i\phi}$ term).
Therefore, $\sigma(p-1)=p-2$. By a similar logic,
$\sigma(p-2)\in\\{p-3,p-2,p-1\\}$ but $\sigma(p-2)\neq p-2$ and
$\sigma(p-2)\neq p-1$. So we get $\sigma(p-2)=p-3$. Proceed in this way to get
the permutation, $\sigma=\left(0,p-1,p-2,\cdots 1\right)$ which has
$\operatorname{sgn}(\sigma)=(-1)^{p-1}$. This corresponds to the term,

$-(-1)^{p}e^{-i\phi}\psi_{-1}^{(0)}(\nu)\cdots\psi_{-1}^{(p-1)}(\nu)$

The term that contains $\psi_{1}^{(p-1)}(\nu)e^{i\phi}$ but not
$\psi_{-1}^{(0)}(\nu)e^{-i\phi}$ is computed in a similar way, and seen to be,

$-(-1)^{p}e^{i\phi}\psi_{1}^{(0)}(\nu)\cdots\psi_{1}^{(p-1)}(\nu)$

We define $r_{\phi}$ so that $r_{0}(\nu)=0$.

From Proposition 2.3 it follows immediately that the constant and linear terms
of $s(\nu)$ both vanish. ∎

A consequence of this Proposition is the following.

######  Corollary 2.6.

$\left.\frac{d^{k}P_{\phi}}{d\phi^{k}}\right|_{\phi=0}=(-1)^{p+1}i^{k}\left(\psi_{1}^{(0)}(\nu)\psi_{1}^{(1)}(\nu)\cdots\psi_{1}^{(p-1)}(\nu)\right)+(-1)^{p+1}(-i)^{k}\left(\psi_{-1}^{(0)}(\nu)\psi_{-1}^{(1)}(\nu)\cdots\psi_{-1}^{(p-1)}(\nu)\right)$

From this we can prove a necessary condition for stability.

######  Theorem 2.7.

If, for a general (linear) nearest neighbor system,

$\displaystyle\prod_{i}\;\rho_{x,1}^{(i)}-\prod_{i}\;\rho_{x,-1}^{(i)}\neq 0$

the system is unstable in one sense or another.

###### Proof.

By [2] (specifically, see Appendix in [2]), the constant term of
$\left.\frac{dP_{\phi}}{d\phi}\right|_{\phi=0}$ must vanish. By Corollary 2.6
the derivative of the constant term is,

$(-1)^{p+1}i\left(\rho_{x,1}^{(0)}\rho_{x,1}^{(1)}(\nu)\cdots\rho_{x,1}^{(p-1)}-\rho_{x,-1}^{(0)}\rho_{x,-1}^{(1)}\cdots\rho_{x,-1}^{(p-1)}\right)$

The $\rho_{x,1}^{(\alpha)},\rho_{x,-1}^{(\alpha)}$ are all real so the theorem
follows. ∎

In this general case it is difficult to come up with sufficient conditions for
stability. If we simplify the problem a bit there is more that can be said.

##  3 $L_{x}$ and $L_{v}$ are Symmetric Laplacians

In [4] we showed that when $p=1$ then stable systems must have
$\rho_{x,1}=\rho_{x,-1}=-1/2$. Theorem 2.7 indicates that this is a reasonable
assumption and so, henceforth, we shall assume that,

$\rho_{x,1}^{(\alpha)}=\rho_{x,-1}^{(\alpha)}=-\frac{1}{2}\text{ for all
}\alpha=0,\cdots p-1.$ (13)

To make the problem tractable we shall also assume that,

$\rho_{v,1}^{(\alpha)}=\rho_{v,-1}^{(\alpha)}=-\frac{1}{2},\text{ for all
}\alpha=0,\cdots p-1.$ (14)

With these two assumptions $L_{x}$ and $L_{v}$ are symmetric row-sum zero
matrices. With these assumptions we, also, have,

$\psi_{1}^{(\alpha)}(\nu)=\psi_{-1}^{(\alpha)}(\nu)\text{ for all
}\alpha=0,\cdots p-1\text{ and for all }\nu.$ (15)

One nice feature of the symmetric case is that one can prove stability in a
restrictive sense using the following fact.

######  Proposition 3.1.

If $M$ is a diagonalizable matrix with eigenvalues in the left half complex
plane and $G$ is a positive definite matrix then $GM$ has eigenvalues in the
left half complex plane.

###### Proof.

Variants of this Proposition are known. We include the proof for completeness.
If $G$ is a positive definite matrix then there is a non-singular square root
$G^{1/2}$. For any $x$ there is a $y$ with $x=G^{1/2}y$. So, for any $x$ we
have,

$\displaystyle\langle{GMy,y}\rangle$
$\displaystyle=\langle{G^{1/2}G^{1/2}My,y}\rangle=\langle{G^{1/2}My,G^{1/2}y}\rangle=\langle{G^{1/2}MG^{-1/2}G^{1/2}y,G^{1/2}y}\rangle$
$\displaystyle=\langle{G^{1/2}MG^{-1/2}x,x}\rangle.$

The eigenvalues of $G^{1/2}MG^{-1/2}$ are the same as $M$. So, for any $y$ we
have,

$\Re\left(\langle{GMy,y}\rangle\right)\leq 0.$

∎

######  Theorem 3.2.

Let $L=L_{x}=L_{v}$ be a diagonalizable Laplacian with eigenvalues in the
closed left half plane. Let $G_{x}=G_{v}=G$ be a diagonal positive matrix and
let $\alpha>0$ be a real number. Then the roots of the characteristic
polynomial of,

$\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\ GL_{x}&\alpha GL_{v}\end{pmatrix},$
(16)

all lie in the closed left half plane. There is a double root at $0$ and all
other roots have real part that is strictly negative.

###### Proof.

In [4] the stability conditions are derived for Linear systems with a matrix,

$\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\ L_{x}&\alpha L_{v}\end{pmatrix}.$
(17)

Because $L_{x}$ and $\alpha L_{v}$ are commuting symmetric matrices there is a
complete set of eigenvectors with eigenvalues $\lambda_{x,m}$ and
$\lambda_{v,m}$ with $\Re(\lambda_{x,m})\leq 0$ and $\Re(\lambda_{v,m})\leq
0$. The eigenvalues of the matrix in equation (17) are given by,

$\nu_{m\pm}=\frac{\lambda_{v,m}}{2}\pm\sqrt{\frac{\lambda_{v,m}^{2}}{4}+\lambda_{x,m}}$

These roots are all in the closed left half plane.

If $G$ is a diagonal positive matrix then

$\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\
GL_{x}&GL_{v}\end{pmatrix}=\begin{pmatrix}\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&G\end{pmatrix}\begin{pmatrix}\mathbf{0}&\mathbf{I}\\\
L_{x}&L_{v}\end{pmatrix}$

Now use Proposition 3.1. ∎

######  Remark 3.3.

With the assumptions in Theorem 3.2 we have,

$\left[GL_{x},\alpha GL_{v}\right]=0.$

##  4 Characteristic Polynomial Expansion

The roots of the characteristic polynomial of the general system, described in
Section 2, are given by the roots of a series of $p$ degree polynomials
$P_{\phi}(\nu)$, as described in Corollary 2.2. The characteristic polynomial
$P_{\phi}(\nu)$ has a root of multiplicity $2$ at $z=0$. In this section we
will assume equations (13) and (14), so that
$\rho_{x,1}^{(\alpha)}=\rho_{x,-1}^{(\alpha)}$ and
$\rho_{v,1}^{(\alpha)}=\rho_{v,-1}^{(\alpha)}$. We will not assume
$G_{x}=\alpha G_{y}$, as in Theorem 3.2. Instead we will let $g_{x}^{\alpha}$
and $g_{v}^{\alpha}$ be independent random variables.

If the system is stable, then roots of the characteristic polynomial all have
non-positive real parts. Stable roots with large negative real part decay
quickly so a stable system is dominated by roots that lie near the imaginary
axis. In our system roots near the double zero will dominate the dynamics so
we will expand around this zero to approximate the system dynamics. The
details of the expansion are outlined in this section.

When $\phi=0$ there is a double root of $P_{0}(\nu)$ at $\nu=0$. We would like
to compute what happens to this double zero when $\phi=\frac{2\pi}{q}m$ is
small but non-zero. If $q$ is large enough we can approximate the system by
using a continuous variable for $\phi$. With this approximation, each of the
roots at $\phi=0,\nu=0$ form continuous zero loci as $\phi$ varies. We get two
zero loci, which are functions,

$\gamma:I\to\mathbb{C},$ (18)

where $I=(-\epsilon,+\epsilon)$ is some neighborhood of $0$. These curves
satisfy,

$P_{\phi}(\gamma(\phi))=0.$ (19)

The coefficients of the characteristic polynomial are analytic functions of
$\phi$ so we will expand everything in Taylor series and use this equation to
deduce conditions on the coefficients. Assume that $\gamma(0)=0$ and write the
expansion,

$\gamma(\phi)=\gamma^{\prime}(0)\phi+\frac{1}{2}\gamma^{\prime\prime}(0)\phi^{2}+\cdots.$
(20)

The coefficients of the polynomial $p_{\phi}(\nu)$ are real analytic functions
of $\phi$. This means we can expand each of them in a Taylor series. The
result is an expansion of the form,

$P_{\phi}(\nu)=\left(a_{00}+a_{01}\phi+\cdots\right)+\left(a_{10}+a_{11}\phi+\cdots\right)\nu+\left(a_{20}+a_{21}\phi+\cdots\right)\nu^{2}+\cdots$
(21)

where the coefficients $a_{0k},a_{1k},\cdots,$ arise as the coefficients of
the $k$th derivative of $p_{\phi}(\nu)$ with respect to $\phi$,

$\frac{1}{k!}\left.\frac{d^{k}P_{\phi}}{d\phi^{k}}\right|_{\phi=0}=a_{0k}+a_{1k}\nu+a_{2k}\nu^{2}+\cdots.$
(22)

######  Proposition 4.1.

The expansion of the zero locus to second order gives the coefficients,

$\displaystyle\gamma^{\prime}(0)$
$\displaystyle=\pm\sqrt{-\frac{a_{02}}{a_{20}}}$ (23)
$\displaystyle\gamma^{\prime\prime}(0)$
$\displaystyle=-\frac{a_{30}(\gamma^{\prime}(0))^{3}+a_{21}(\gamma^{\prime}(0))^{2}+a_{12}\gamma^{\prime}(0)+a_{03}}{a_{20}\gamma^{\prime}(0)}$
(24)

###### Proof.

The equation $P_{\phi}(\gamma(\phi))=0$ expands to a power series in $\phi$.
Set $\nu=\gamma(\phi)$ and expand using equation (20). Plug this value of
$\nu$ into the polynomial of equation (21). Condition 19 says this expansion
in $\phi$ vanishes identically. When you solve for the derivatives
$\gamma^{(m)}(0)$ in terms of $a_{jk}$ you get equations (23) and (24). This
is a straightforward calculation and you can check the results using algebraic
manipulation software, like SAGE. ∎

Notice that there are two solutions to the first order coefficient $c_{1}$.
There is a double root that splits into two distinct curves, so there are two
distinct values for $\gamma^{(1)}(0)$ and $\gamma^{(2)}(0)$.

At this point we find the coefficients $a_{jk}$ that are required for our
expansion. The reader may find it more digestible to jump to Section 5 and
refer to the remaining portion of this section as needed. We start with the
following Proposition.

######  Proposition 4.2.

If we assume equation (13) and (14) then we have,

$\displaystyle a_{mk}$ $\displaystyle=0,\text{ for }k\text{ odd.}$ (25)
$\displaystyle a_{02}$
$\displaystyle=(-1)^{p}(\rho_{x,1})^{p}g_{x}^{(0)}\cdots
g_{x}^{(p-1)}=2^{-p}\prod_{j=0}^{p-1}{g_{x}^{(j)}}$ (26) $\displaystyle
a_{12}$
$\displaystyle=2^{-p}\prod_{j=0}^{p-1}{g_{x}^{(j)}}\sum\limits_{k=0}^{p-1}{\frac{g_{v}^{(k)}}{g_{x}^{(k)}}}$
(27)

###### Proof.

Since $\psi_{1}(\nu)=\psi_{-1}(\nu)$ then, by Corollary 2.6, we have

$\frac{1}{k!}\left.\frac{d^{k}P_{\phi}}{d\phi^{k}}\right|_{\phi=0}=\frac{1}{k!}(-1)^{p+1}(i^{k}+(-i)^{k})\left(\psi_{1}^{(0)}(\nu)\psi_{1}^{(1)}(\nu)\cdots\psi_{1}^{(p-1)}(\nu)\right).$

The coefficients are just the coefficients of this polynomial in $\nu$. The
coefficient $a_{mk}$ is the $m$’th polynomial coefficient of the polynomial of
the $k$’th derivative. The coefficients $a_{mk}=0$ whenever $k$ is odd. When
$k=2$ we have the polynomial,

$\displaystyle\frac{1}{2}\left.\frac{d^{2}P_{\phi}}{d\phi^{2}}\right|_{\phi=0}$
$\displaystyle=(-1)^{p}\left(\psi_{1}^{(0)}(\nu)\psi_{1}^{(1)}(\nu)\cdots\psi_{1}^{(p-1)}(\nu)\right)$
$\displaystyle=a_{02}+a_{12}\nu+\cdots.$

To form the linear term in $\nu$ we can take the linear term in each
$\psi_{1}^{(\alpha)}(\nu)$ where all the other terms contribute a constant
term. The result is that the linear term is

$a_{12}=(-1)^{p}\sum\limits_{j=0}^{p-1}{g_{v}^{j}\rho_{v,1}^{j}\frac{g_{x}^{0}\cdots
g_{x}^{(p-1)}}{g_{x}^{j}}\left(\rho_{x,1}\right)^{p-1}}=\frac{1}{2^{p}}\prod_{j=0}^{p-1}{g_{x}^{(j)}}\sum\limits_{j=0}^{p-1}{\frac{g_{v}^{(j)}}{g_{x}^{(k)}}}.$

The term $a_{02}$ comes from the constant term which is computed in a similar
way. ∎

######  Proposition 4.3.

Assume equation (13) and (14). The degree $2$ term of the characteristic
polynomial, when $\phi=0$, is given by

$a_{20}=\frac{p}{2^{p-1}}\prod_{j=0}^{p-1}{g_{x}^{(j)}}\left(\sum\limits_{k=0}^{p-1}{\frac{1}{g_{x}^{(k)}}}\right)=\frac{p^{2}}{2^{p-1}}\prod_{j=0}^{p-1}{g_{x}^{(j)}}\operatorname{Avg}\left(\frac{1}{g_{x}^{(k)}}\right)$
(28)

The degree $3$ term of the characteristic polynomial, when $\phi=0$, has the
value

$\displaystyle a_{30}$
$\displaystyle=\frac{(2p)}{2^{p}}\left(g_{x}^{(0)}\cdots
g_{x}^{(p-1)}\right)\sum\limits_{k=0}^{p-1}\sum\limits_{\begin{subarray}{c}j=0\\\
j\neq k\end{subarray}}^{p-1}{\frac{g_{v}^{(j)}}{g_{x}^{(k)}g_{x}^{(j)}}}$
$\displaystyle=\frac{(2p)}{2^{p}}\prod_{j=0}^{p-1}{g_{x}^{(j)}}\left(\sum\limits_{k=0}^{p-1}{\frac{1}{g_{x}^{(k)}}}\sum\limits_{j=0}^{p-1}{\frac{g_{v}^{(j)}}{g_{x}^{(j)}}}-\sum\limits_{k=0}^{p-1}{\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}}\right)$
(29)

where the last sum is only over $(i,j)$ where $0=i<j\leq(p-1)$.

###### Proof.

To compute the second order term in $s(\nu)$ we see, from Proposition 2.3 that
all $2$nd order terms must contain a $\nu^{2}$ from the diagonal and all the
remaining terms are constants. Similarly, to compute the $3$rd order terms in
$s(\nu)$, by Proposition 2.3, all the cubic terms that do not include a
$\nu^{2}$ diagonal factor must sum to zero. So the $3$rd order terms include a
$\nu^{2}$ term from the diagonal and a $\nu$ factor from one of the $\psi$’s.

The co-factor of the diagonal term have the form assumed in Proposition 4.4.
The $\nu^{2}$ in the $(k,k)$ position has a co-factor of the form,

$\displaystyle(\operatorname{sgn}(\eta))^{2}\det\left|\begin{matrix}\psi_{0}^{(k+1)}(\nu)&\psi_{1}^{(k+1)}(\nu)&0&0&\cdots&0&0\\\
\psi_{-1}^{(k+2)}(\nu)&\psi_{0}^{(k+2)}(\nu)&\psi_{1}^{(k+2)}(\nu)&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\cdots&\vdots&\vdots\\\
0&0&0&0&\cdots&\psi_{0}^{(k-1)}(\nu)&\psi_{1}^{(k-1)}(\nu)\end{matrix}\right]$
$\displaystyle=$
$\displaystyle(-1)^{p-1}\left(\psi_{1}^{(k+1)}(\nu)\cdots\psi_{1}^{(k+p-1)}(\nu)+\cdots+\psi_{-1}^{(k+1)}(\nu)\cdots\psi_{-1}^{(k+p-1)}(\nu)\right)$

where all the coefficient arithmetic is assumed $\pmod{p}$ and $\eta\in S_{p}$
is the permutation that rotates the indices so that $k$ is rotated to position
$0$. Each term in the last sum has $p-1$ factors. The constant term is easily
isolated and all $p$ terms in the sum are equal to

$(-1)^{p-1}(\rho_{x,1})^{p-1}p\frac{g_{x}^{(0)}g_{x}^{(1)}\cdots
g_{x}^{(p-1)}}{g_{x}^{(k)}}=\frac{p}{2^{p-1}}\frac{g_{x}^{(0)}g_{x}^{(1)}\cdots
g_{x}^{(p-1)}}{g_{x}^{(k)}}.$

Equation (28) follows.

For the cubic term, exactly one of the $\psi$ terms will contribute a $\nu$,
and so a $\rho_{v,\pm 1}$, which we label the $j$th term. The contribution
from the $k$’th $\nu^{2}$ diagonal will have $p-2$ factors of $\rho_{x,\pm 1}$
and has the form,

$\displaystyle(-1)^{p-1}(\rho_{x,1})^{p-2}$
$\displaystyle\sum\limits_{j=k+1}^{k-1+p}{\frac{g_{x}^{(0)}\cdots
g_{x}^{(p-1)}}{g_{x}^{(k)}g_{x}^{(j)}}g_{v}^{(j)}\left(\rho_{v,1}(p-1-j)+\rho_{v,-1}(j+1)\right)}$
$\displaystyle=\frac{p}{2^{p-1}}\sum\limits_{\begin{subarray}{c}j=0\\\ j\neq
k\end{subarray}}^{p-1}{\frac{g_{x}^{(0)}\cdots
g_{x}^{(p-1)}}{g_{x}^{k}g_{x}^{j}}g_{v}^{j}}$

recall that all index arithmetic is $\mod(p)$. The formula for $a_{30}$
follows. ∎

######  Proposition 4.4.

Let $D_{n}$ be the determinant,

$D_{n}=\left|\begin{matrix}(c_{1}+d_{1})&-c_{1}&0&0&\cdots&0&0\\\
-d_{2}&(d_{2}+c_{1})&-c_{2}&0&\cdots&0&0\\\
0&-d_{3}&(d_{3}+c_{3})&-c_{3}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&(d_{n-1}+c_{n-1})&-c_{n-1}\\\
0&0&0&0&\cdots&-d_{n}&(d_{n}+c_{n})\end{matrix}\right|$ (30)

then the determinant is given by,

$D_{n}=(d_{1}\cdots d_{n})+(d_{1}\cdots d_{n-1}c_{n})+\cdots+(d_{1}c_{2}\cdots
c_{n-1}c_{n})+(c_{1}\cdots c_{n}).$

###### Proof.

We proceed by induction and use the general form for determinants of tri-
diagonal matrices. The formula is easily derived from the definition of the
determinant and has the form,

$D_{n}(1,n)=(d_{1}+d_{1})D_{n-1}(2,n)-c_{1}d_{2}D_{n-2}(3,n),$ (31)

where $D_{n-1}(k_{1},k_{2})$ is the determinant of the $k_{2}-k_{1}+1$ square
matrix with row and column indices between $k_{1}$ and $k_{2}$ inclusive. The
cases $n=3$ and $n=4$ are easily computed directly. The case for general $n$
can be proved by induction using equation (31). ∎

##  5 Characteristic Polynomial Near 0

The characteristic equation of the system is given in Corollary 2.2. The
Taylor expansion of this equation in the variable $\phi$, described in Section
4, results in an approximation for the roots that lie near the double root at
$\nu=0$. The main result is given in Theorem 5.1. See Section 4 for a detailed
derivation of some of the coefficients needed in the expansion.

######  Theorem 5.1.

With the assumptions described at the start of Section 3 the characterisitic
polynomial $P_{\phi}(\nu)$ has a double zero when $\phi=0$. The coefficients
$c_{1},c_{2}$, in the expansion of equation (20) are given by,

$\displaystyle\gamma^{\prime}(0)$
$\displaystyle=\pm\frac{i}{p}\sqrt{\frac{1}{2\operatorname{Avg}\left(1/g_{x}^{(k)}\right)}}$
(32) $\displaystyle\gamma^{\prime\prime}(0)$
$\displaystyle=-\frac{\operatorname{Avg}\left(\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}\right)}{2\left(p\operatorname{Avg}\left(\frac{1}{g_{x}^{(k)}}\right)\right)^{2}}$
(33)

This means that the zero locus, near the double root, is approximated by (see
equation (20)),

$\gamma(\phi)=\pm
i\sqrt{\frac{1}{2\operatorname{Avg}\left(1/g_{x}^{(k)}\right)}}\left(\frac{\phi}{p}\right)-\frac{\operatorname{Avg}\left(\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}\right)}{\left(2\operatorname{Avg}\left(1/g_{x}^{(k)}\right)\right)^{2}}\left(\frac{\phi}{p}\right)^{2}+\cdots,$
(34)

###### Proof.

To find $\gamma^{\prime}(0)$ we start with Proposition 4.1. Then use
Proposition 4.3 equation (28) and Proposition 4.2 equation (26). We get,

$\displaystyle\gamma^{\prime}(0)$
$\displaystyle=\pm\sqrt{-\frac{a_{02}}{a_{20}}}=\pm\sqrt{-\frac{1}{(2p)\sum\limits_{k=0}^{p-1}{1/g_{x}^{(k)}}}}$

Equation (32) follows.

To find $\gamma^{\prime\prime}(0)$ use Propositions 4.1 and 4.1.

$\displaystyle\gamma^{\prime\prime}(0)$
$\displaystyle=-\frac{a_{30}(\gamma^{\prime}(0))^{3}+a_{21}(\gamma^{\prime}(0))^{2}+a_{12}\gamma^{\prime}(0)+a_{03}}{a_{20}\gamma^{\prime}(0)}=-\frac{a_{30}(\gamma^{\prime}(0))^{2}+a_{12}}{a_{20}}=\frac{a_{30}a_{02}-a_{12}a_{20}}{a_{20}^{2}}$

The coefficients $a_{20}$ and $a_{30}$ are computed in Proposition 4.3. The
coefficients $a_{12}$ and $a_{02}$ are computed in Proposition 4.2. Assemble
the coefficients, cancel the terms $\prod_{j=0}^{p-1}{g_{x}^{(j)}}$ and
$2^{p}$ and we get

$\displaystyle\gamma^{\prime\prime}(0)$
$\displaystyle=\frac{a_{30}a_{02}-a_{12}a_{20}}{a_{20}^{2}}$
$\displaystyle=\frac{1}{4p^{4}\operatorname{Avg}\left(1/g_{x}^{(k)}\right)^{2}}\left((2p)\left(\sum\limits_{k=0}^{p-1}{\frac{1}{g_{x}^{(k)}}}\sum\limits_{j=0}^{p-1}{\frac{g_{v}^{(j)}}{g_{x}^{(j)}}}-\sum\limits_{k=0}^{p-1}{\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}}\right)-\left((2p)\sum\limits_{j=0}^{p-1}{\frac{g_{v}^{(j)}}{g_{x}^{(j)}}}\sum\limits_{k=0}^{p-1}{\frac{1}{g_{x}^{(k)}}}\right)\right)$
$\displaystyle=-\frac{(2p)\sum\limits_{k=0}^{p-1}{\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}}}{4p^{4}\operatorname{Avg}\left(1/g_{x}^{(k)}\right)^{2}}=-\frac{\operatorname{Avg}\left(\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}\right)}{2p^{2}\operatorname{Avg}\left(1/g_{x}^{(k)}\right)^{2}}$

∎

######  Remark 5.2.

Notice that the derivative $\gamma^{\prime}(0)$ is imaginary and the second
order $\gamma^{\prime\prime}(0)$ is real and negative. This is consistent with
the stability of Theorem 3.2.

We define two positive values $c_{1},c_{2}$ by,

$\displaystyle
c_{1}=p|\gamma^{\prime}(0)|=\sqrt{\frac{1}{2\operatorname{Avg}\left(1/g_{x}^{(k)}\right)}}$
(35) $\displaystyle
c_{2}=p^{2}|\gamma^{\prime\prime}(0)|=\frac{\operatorname{Avg}\left(\frac{g_{v}^{(k)}}{g_{x}^{(k)}g_{x}^{(k)}}\right)}{\left(2\operatorname{Avg}\left(\frac{1}{g_{x}^{(k)}}\right)\right)^{2}}$
(36)

Recall that $\phi=\frac{2\pi}{q}m$. We re-write equation (34) as,

$\displaystyle\gamma_{m}$ $\displaystyle=\pm
ic_{1}\left(\frac{2\pi}{pq}m\right)-c_{2}\left(\frac{2\pi}{pq}m\right)^{2}+\cdots=\pm
ic_{1}\left(\frac{2\pi}{N}m\right)-c_{2}\left(\frac{2\pi}{N}m\right)^{2}+\cdots$
(37)

Recall that we shuffled coordinates to get the eigenvector in equation 10.
This eigenvector, in the original coordinates, is

$\begin{pmatrix}\mathbf{\epsilon}\otimes\mathbf{v}_{m}(\phi)\\\
\nu\mathbf{\epsilon}\otimes\mathbf{v}_{m}(\phi)\end{pmatrix}$

The solutions have time dependence through the factor,

$\exp\left(\gamma(\phi)t\right)\exp\left(\frac{2\pi i}{q}\mu
m\right)=\exp\left(\frac{2\pi i}{pq}m\left(\pm
c_{1}t+p\mu\right)\right)\exp\left(-c_{2}\left(\frac{2\pi}{pq}m\right)^{2}t\right)$

The signal velocity is determined by the factor $\pm c_{1}t+p\mu$. There are
two waves going in opposite directions with equal velocities. We see that the
time for the signal to travel $p$ agents is given by $c_{1}T_{p}=p$ so that,

$T_{1}=\frac{p}{c_{1}}=p\sqrt{2\operatorname{Avg}\left(1/g_{x}^{(k)}\right)}.$
(38)

The factor, $\exp\left(-c_{2}\left(\frac{2\pi}{pq}m\right)^{2}t\right)$ is a
damping term that introduces a dispersion relation. The larger $c_{2}$ is the
larger the damping factor.

In our original problem we assumed the $N=pq$ agents were assembled as $q$
groups of $p$ randomly weighted agents. The values $c_{1}$ and $c_{2}$ are
averages over the $p$ agents which equals the average over the $N$ agents.
What if we have $N$ randomly weighted agents. This amounts to setting $q=1$.
Does equation (37) still apply? In the following section we shall run a
variety of simulations on $N$ uniquely weighted agents and compare to our
solutions.

##  6 Simulations

The computation in sections 4 and 5 assumed a condition that is a bit stronger
than periodic boundary conditions. It assumed that $p$ independently weighted
agents were duplicated $q$ times for a total sequence of $N=pq$ agents and
that this $N$-string satisfies periodic boundary conditions. In this
simulation section we shall forgo this assumption and shall set $q=1$ and set
the sequence on the real line without periodic boundary conditions. The
scales, $g_{x}^{(\alpha)},g_{v}^{(\alpha)}$, are chosen independently. We will
demonstrate that the periodicity of the $p$-sequence is not required in
practice. In these simulations we shall assume the conditions in equations
(13) and (14) so that $\rho_{x,1}=\rho_{x,-1}$ and $\rho_{v,1}=\rho_{v,-1}$.

(a) Flock with times indicated

(b) Flock of 1000 agents

Figure 1: Flock Examples

(a) Estimate for $T_{1}$ when $g_{x}$ is random from uniform distribution

(b) Estimate for $T_{1}$ when $g_{x}$ is a ramp

Figure 2: Prediction of the first response time $T_{1}$ compared to
computational predictions

In Figure 1(a) we show a linear flock with $N=p=200$ agents. The plot shows
the simulation estimate for $T_{1}$, which we derived in equation (38). The
signal velocity actually depends on the constituent waves in the wave packet
that moving from the leader to the tail. The simple model used to estimate
$T_{1}$ in the simulations is to pick the time where the distance to the
leader is the greatest. Also shown in the plot are the distance used to
estimate the period and the peaks used to evaluate amplitude ratios. We
estimate the amplitude by choosing the point with the largest distance to the
leader, just as we did to estimate $T_{1}$. The values $g_{x}^{(\alpha)}$ are
independently chosen from a uniform distribution on $[0.50,10.50]$ and the
$g_{v}^{(\alpha)}$ from $[0.50,1.50]$. For this particular snapshot of random
variables the value $T_{1}$ in equation (38) gives $T_{1}=153.0$ which is
extremely close to the observed value. In this graph each agent trajectory is
in a different color so that flock dynamics are evident.

In Figure 1(b) we have a basic simulation of $N=1000$. The weights for
$g_{x}^{(\alpha)}$ are independently chosen from a uniform distribution on
$[1,12]$ and $g_{v}^{(\alpha)}$ from a uniform distribution on $[1,12]$. In
this simulation the estimated $T_{1}=720$ and the computed value is
$T_{1}=721$, which is, again, quite close.

Equation (38) gives a prediction for the first response time. In figure 2(a)
we compare the computed value of $T_{1}$ to the simulation estimate, as we
vary the mean of $g_{x}^{(\alpha)}$. For each test we select $N=p$ values
$g_{x}^{(\alpha)}$ independently from a uniform distribution with a mean shown
in the x-axis. In figure 2(b) we run simulations where $g_{x}^{(\alpha)}$ form
a ramp as $\alpha=0,\cdots p-1$. The ramp starts with $g_{x}^{(0)}$ small and
increases so that $g_{x}^{(p-1)}$ has a maximum value according to the
formula,

$g_{x}^{(\alpha)}=g_{low}+\alpha\frac{(g_{hi}-g_{low})}{p-1},\;\;\alpha=0,1,\cdots,p-1.$

The $T_{1}$ estimates for the random distribution are very close to the
predicted value. The ramp “distribution” is not random and it has a large
discontinuity at $\alpha=0$. Although these are clear differences from the
uniform distribution, we do not yet understand why the ramp distribution makes
the $T_{1}$ prediction less accurate.

(a) Mean $g_{v}$ is smaller so $c_{2}$ is smaller

(b) Mean $g_{v}$ is larger so $c_{2}$ is larger

Figure 3: Changing $c_{2}$ changes dispersion relation

The value $c_{1}$ in equation (35) determines the first response as it appears
in an imaginary expential. The value $c_{2}$ in equation (36) appears in a
real exponential and is negative which indicates a stable solution. This is a
dispersion term and tends to disperse the pure frequencies so that the
amplitude edges are muted. Figure 3(a) shows a plot with small $c_{2}$ and
figure 3(b) is the same system except the mean value of $g_{v}^{(\alpha)}$ is
larger. This will increase the value of $c_{2}$ and the resulting peaks are
rounder.

(a) One agent (in green) has negative weight

(b) Set $\rho_{x,2}^{(j)}=0.25$ for agents $j=500,\cdots,500+(L-1)$ where
$L=0,\cdots 11$. The more rows with $\rho_{x,1}=0.25$ the faster the onset of
instability.

Figure 4: Simulations demonstrating instability

The numerics indicate that the system is stable only when $g_{x}^{(\alpha)}$
are non-negative for all $\alpha$. Setting $g_{x}^{(\alpha)}$ negative, for a
single agent has an immediate effect on system stability. Figure 4(a) shows a
simulation where a single agent has negative weight. The plot shows a subset
of the flock. The agent with negative weight is shown in green. When the
signal to move reaches this agent, the agent moves in the opposite direction
and the system starts on an unstable trajectory. If the weight of this single
agent was zero then the agents to the left of it would never react to the
motion of the leading agent since we are only including nearest neighbor
interactions.

Theorem 2.7 states a necessary condition for stability. The difficulty with
the theorem is that both $\prod_{i}\;\rho_{x,1}^{(i)}$ and
$\prod_{i}\;\rho_{x,-1}^{(i)}$ are very small and go to zero as $N\to\infty$.
If we have $\rho_{x,1}=\rho_{x,-1}$ then the condition is satisfied but
changing a single value of $\rho_{x,1}^{(\alpha_{0})}$ to $-0.25$ should
result in an unstable system. This is difficult to show numerically. In Figure
4, we plot the middle agent for a series of simulations. Each simulation has
$\rho_{x,1}=\rho_{x,-1}=-1/2$ except for $L$ agents centered around agent
$500$. For the given simulations, the system is unstable when $L\geq 5$. The
greater $L$ the faster the onset of the instability.

(a) Plot showing $\alpha$ as $N\to\infty$

(b) Simulation of truck convoy with 400 agents.

A previous work [8] derives a condition on the ratio of successive flock
amplitudes. There it is shown that as $N\to\infty$ then the ratio goes to
$1.0$. In Figure 5(a) we plot amplitudes $A_{4}/A_{3}$ for various
distributions of $g_{v}$. We do not use the ratio $A_{2}/A_{1}$ since the
first amplitude is close to the initial condition and behaves differently than
subsequent peaks. The curves depend on the mean of $g_{v}^{(\alpha)}$ as
altering $g_{v}^{(\alpha)}$ alters $c_{2}$ (see equation (36)). We cannot
conclusively conclude that this is exactly correct, and more simulation work
is required.

We conclude with the realistic system discussed in the introduction. With this
example we demonstrate that the tools presented in this paper can be used to
analyze more complicated and realistic problems. We make some rough estimates
in this next section. A automotive engineer could refine these numbers with
more realistic estimates. We model a convoy of $N$ trucks traveling on the
highway. The convey attempts to keep a fixed spacing between trucks and the
trucks are all different. As the convoy travels, lighter cars might,
inadvertently, enter the convoy creating a 1-dimensional convoy with very
different agents. To use our model, we must estimate the agent weights
$g_{x}^{k}$ and $g_{v}^{k}$. The weights for agent $k$ are force coupling
divided by the mass of the agent. The mass of a 18 wheel truck is somewhere
between 14 and 40 thousand kilograms and the coupling force is determined by
the torque of the engine. To make things simpler we shall assume the force
divided by the mass produces a given acceleration and we can estimate this
acceleration. For example, a truck can accelerate from $0$ to $60$ mph =
$26.8$ m/sec in 1 to 5 minutes. So we take our truck weights $g_{x}^{k}$ in
the range,

$g_{x}^{k}\in[26.8/60,26.8/300],\;\;k\text{ is a truck}.$

We insert cars into the convoy by randomly replacing $10\%$ of the agents with
lighter cars. Cars, typically, have higher power to mass and so have larger
weights. We take a collection of cars that accelerate from $0$ to $60$ in a
range of $6$ to $20$ seconds, so that for car agents,

$g_{x}^{k}\in[26.8/6,26.8/20],\;\;k\text{ is a car}.$

To guarantee stability we take $G_{v}=\alpha G_{x}$ where $\alpha=10.0$.
Increasing $\alpha$, as we’ve seen, will increase the damping. U.S. 18-wheel
trucks are typically around $23$ meters long. The convoy attempts to keep a
distance of $3\times 23=69$ meters between the agents.

The simulation results are shown in figure 5(b). The convey of $400$ starts
with a spacing of $69$ meters so is almost $28km$ long. The first track
suddenly increases its speed $10$ meters/second and it takes $1095$ seconds
for the signal to reach the last truck. The time duration is long because the
weights are small (e.g the trucks accelerate slowly). The distance from the
leader to the tail lengthens to $37.5km=94$m/truck before the tail starts to
catch up. The damping is not critical, and the tail overshoots the optimal
distance and the entire convey shrinks to $18km=45$ m/truck before expanding
again. This simulation assumed $\rho_{v,1}=\rho_{v,-1}$. There is some
evidence that system is more responsive if this condition is removed. We will
explore that in a subsequent paper.

## References

  * [1] N.W. Ashcroft and N.D. Mermin.  Solid State Physics.  Saunders College, Philadelphia, 1976.
  * [2] Pablo E. Baldivieso and J.J.P. Veerman.  Stability conditions for coupled autonomous vehicles formations.  IEEE Transactions on Control of Network Systems, in press.
  * [3] C.E. Cantos, D.K Hammond, and J.J.P. Veerman.  Transients in the synchronization of assymmetrically coupled oscillator arrays.  The European Physical Journal Special Topics, 225:1115–1125, 2016\.
  * [4] C.E. Cantos, J.J.P. Veerman, and D.K. Hammond.  Signal velocity in oscillator arrays.  The European Physical Journal Special Topics, 225:1199–1210, 2016\.
  * [5] Robert E. Chandler, Robert Herman, and Elliott W. Montroll.  Traffic dynamics: Studies in car following.  Operations Research, 6:165–184, 1957.
  * [6] Freeman J. Dyson.  The dynamics of a disordered linear chain.  Phys. Rev., 92:1331–1338, Dec 1953.
  * [7] He Hao and Prabir Barooah.  Stability and robustness of large platoons of vehicle with double-integrator models and nearest neighbor interaction.  International Journal of Robust and Nonlinear Control, 23, 12 2013\.
  * [8] J. Herbrych, A. Chazirakis, N. Christakis, and J.J.P. Veerman.  Dynamics of locally coupled agents with next nearest neighbor interaction.  Differential Equations and Dynamical Systems, 2017.
  * [9] B.S. Kerner.  The Physics of Traffic.  Springer-Verlag, 2004.
  * [10] S. E. Li, Y. Zheng, K. Li, Y. Wu, J. K. Hedrick, F. Gao, and H. Zhang.  Dynamical modeling and distributed control of connected and automated vehicles: Challenges and opportunities.  IEEE Intelligent Transportation Systems Magazine, 9(3):46–58, 2017\.
  * [11] Fu Lin, Makan Fardad, and Mihailo R. Jovanovic.  Optimal control of vehicular formations with nearest neighbor interactions.  IEEE Transactions on Automatic Control, 57(9):2203–2218, 2012.
  * [12] Uwe Mackenroth.  Robust Control Systems.  Springer-Verlag, 2004.
  * [13] Wei Ren and Randal Beard.  Distributed Consensus in Multi-vehicle Cooperative Control.  Springer-Verlag, 2008.
  * [14] S. Studli, M. M. Seron, and R. H. Middleton.  Vehicular platoons in cyclic interconnections.  Automatica, 94:283–293, 2018.