X_{k}\rangle\right)\right)^{2}\right\rangle-\left\langle 2\sum_{i\neq j\neq
k}X_{i}^{2}X_{j}X_{k}\right\rangle^{2}.$ (783)

The dominant terms of the first part are the cross terms, whose number is at
the same order of $N$ as the second part. So to evaluate the approximation we
can only consider the two terms below

$\displaystyle\langle\left(X_{i}^{2}-\phi_{i}\right)X_{j}X_{k}+\langle
X_{i}\rangle^{2}\langle X_{j}\rangle\langle X_{k}\rangle\rangle,\quad\langle
X_{i}^{2}X_{j}X_{k}\rangle,$ (784)

where the variables are identical and independent. Recall that we have the
equation (747) therefore the two expressions are equal under the average.
Explicitly the cross terms also contain many parts, we can also write down the
one with the largest number of terms,

$\displaystyle\left(X_{i}^{2}-\phi_{i}\right)X_{j}X_{k}\left(X_{l}^{2}-\phi_{l}\right)X_{m}X_{n},\quad
4{N\choose 2}(N-2){N-3\choose 2}(N-5),$ (785)

while the second part in (783) contains

$\displaystyle 4{N\choose 2}^{2}(N-2)^{2}.$ (786)

Therefore the dominant contribution in (783) will cancel out, which makes the
proposal hold.

We can also consider the problem that occurs in the case with $\mu=0$ when the
average of $X^{2}-\phi$ is not zero, and we can find that there’s no such
problem here. As the dominant term (783) contains the largest number of $X$
terms, therefore even when $X^{2}-\phi$ is a constant the total number of $X$
terms in (783) will not change.

Next we consider the case with general $q$. Following the previous procedure
for the error we can get an expression similar to (780)

$\displaystyle\langle\text{Error}^{2}\rangle=\left\langle\left(\sum_{s=0}^{q-1}(Y_{s}^{2}+\langle
Y\rangle^{2}_{s}-\Phi_{s})\right)^{2}\right\rangle-\langle Y^{2}\rangle^{2},$
(787)

and the dominant sum will also be the last non-zero one, i.e. the sum labelled
with $s=q-1$. To find out the behavior of the dominant terms the relation of
the two functions below is important,

$\displaystyle\langle\left(X_{i_{1}}^{2}-\phi_{i_{1}}\right)X_{i_{2}}\dots
X_{i_{2q-1}}+\langle X_{i_{1}}\rangle^{2}\langle X_{i_{2}}\rangle\dots\langle
X_{i_{2q-1}}\rangle\rangle,\quad\langle X_{i_{1}}^{2}X_{i_{2}}\dots
X_{i_{2q-1}}\rangle,$ (788)

we want to determine whether they are equal. And we can find that the relation
(747) still holds thus the approximation is also proper in general $q$.

### SYK with one time point with $\mu=0$

Like in the previous section given any distribution we can always construct a
variable with zero mean, such as

$\displaystyle\tilde{X}=X-\langle X\rangle,$ (789)

so that we can consider any distribution we want.

We first consider a simple case with $q=2$, the Hamiltonian has the form

$\displaystyle Y=\sum_{i\neq j}\text{sgn}(ij)X_{i}X_{j},$ (790)

where $X_{i}$’s are identical and independent variables and the sign term
mimics the sign function in the SYK model. Note that here we have no need to
compute the explicit form of the sign function, since it has no effect yet. It
is very similar to the previous case, but the difference is that here all
$X_{i}$’s only appear once in the sum. Therefore for the square we have

$\displaystyle Y^{2}=\sum_{i\neq j}X_{i}^{2}X_{j}^{2}+2\sum_{i\neq j\neq k\neq
l}\text{sgn}(ijkl)X_{i}X_{j}X_{k}X_{l},$ (791)

comparing to (734) it lacks the middle term. Our previous computation shows
that the approximation fails in the Gaussian case with $\mu=0$ due to the
middle term, thus the approximation is proper in this case no matter what the
distribution is.

We can have an explicit computation for the proposal, but first the convention
needs to be consistent with the SYK model. Now the total number $N$ counts the
whole Majorana fermions and the number $q$ is for the number of the fermions
in the interaction. Then for $N/q=2$ the expression (790) is still valid for
the partition function, except that the $N$ and $q$ have different meanings.
Then the number of $X_{i}$ in (790) becomes

$\displaystyle\frac{N!}{\left(\left(N/2\right)!\right)^{2}},$ (792)

and the sum contains the number of the terms

$\displaystyle\frac{N!}{2!\left(\left(N/2\right)!\right)^{2}}.$ (793)

Following the previous process we have

$\displaystyle Y^{2}-\Phi=\sum_{i\neq
j}\left(X_{i}^{2}X_{j}^{2}-\phi_{i}\phi_{j}\right)$ (794)

where $\phi$ is defined in (744) and for the error

$\displaystyle\langle\text{Error}^{2}\rangle=\left\langle\left(\sum_{i\neq
j}\left(X_{i}^{2}X_{j}^{2}-\phi_{i}\phi_{j}\right)\right)^{2}\right\rangle-\left\langle\sum_{i\neq
j}X_{i}^{2}X_{j}^{2}\right\rangle^{2}.$ (795)

Since there are only two parts in the above function and the sums are the
same, we can only compare the two expressions in the brackets. The average of
$\phi$ is zero due to the zero mean of $X_{i}$, so the proposal holds in the
case. Explicitly the dominant terms in (795) come from the cross term in the
first part and the whole second part, the number becomes

$\displaystyle\frac{N!}{2!\left(\left(N/2\right)!\right)^{2}}\left(\frac{N!}{2!\left(\left(N/2\right)!\right)^{2}}-1\right)-\left(\frac{N!}{2!\left(\left(N/2\right)!\right)^{2}}\right)^{2}.$
(796)

We can find the dominant terms cancel out so that the proposal is valid in
this simple case.

Then we can consider a more general Hamiltonian with arbitrary $N$ and $q$.
Here we take the more familiar convention in SYK and still define $p=N/q$,
then the partition function becomes

$\displaystyle z=\sum_{A}\text{sgn}(A)J_{A_{1}}\dots J_{A_{p}}.$ (797)

Before we check the approximation we can first have some computation about the
numbers of the different terms in this partition function. The total number of
the $J_{A}$’s is

$\displaystyle{N\choose q},$ (798)

and the number of the combinations of $J_{A}$’s is

$\displaystyle\frac{N!}{p!(q!)^{p}},$ (799)

while the times for a single $J_{A}$ appearing in the Hamiltonian is

$\displaystyle\frac{(N-q)!}{(p-1)!(q!)^{p-1}}.$ (800)

The square of the partition function can be given as

$\displaystyle z^{2}$
$\displaystyle=\sum_{A}\text{sgn}(A)^{2}J_{A_{1}}^{2}\dots
J_{A_{p}}^{2}+\sum_{A,B}\text{sgn}(A)\text{sgn}(B)J_{A_{1}}^{2}\dots
J_{A_{p-2}}^{2}J_{A_{p-1}}J_{A_{p}}J_{B_{p-1}}J_{B_{p}}$
$\displaystyle\quad+\dots+\sum_{A,B}\text{sgn}(A)\text{sgn}(B)J_{A_{1}}\dots
J_{A_{p}}J_{B_{1}}\dots J_{B_{p}}$ (801)

the first non-diagonal term appears with four different $J_{A}$’s as the
combination $A$ should be the permutation of $1,\dots,N$. And in the sums
there may be many identical terms, here we may not give the explicit forms. To
figure out the explicit expressions for the numbers in each sum we can count
the numbers of different combinations of $r$ different $J_{A}$’s in the
product of the combinations of $A$ and $B$. Obviously there’s only one
combination in the case $r=1$, actually $r=1$ is identical to $r=0$ since the
last $J_{A_{i}}$ in a combination $A$ can be determined by the other $p-1$
$J_{A_{i}}$’s. For the larger $r$ we have to subtract the combinations in
$r-i$ $(i<r)$ to get the correct ones, as it can give $r-i$ cases in the
product of two combinations. Therefore we have

$\displaystyle r=2,\quad\frac{(2q)!}{2!(q!)^{2}}-1,$ (802) $\displaystyle
r=3,\quad\frac{(3q)!}{3!(q!)^{3}}-{3\choose
2}\left(\frac{(2q)!}{2!(q!)^{2}}-1\right)-1,$ (803)

from the above we can introduce the function $n_{r}$ for this counting

$\displaystyle n_{r}=\frac{(rq)!}{r!(q!)^{r}}-\sum_{s=2}^{r-1}{r\choose
s}n_{s}-1.$ (804)

From the calculation we can take mathematical induction to derive a simpler
form

$\displaystyle n_{r}=\sum_{s=1}^{r}\frac{(sq)!}{s!(q!)^{s}}{r\choose
s}(-1)^{r-s}+(-1)^{r},\quad r>1,$ (805)

and

$\displaystyle n_{0}=1,$ (806)

where $n_{1}$ is absent according to the equation (801).

Then we calculate the numbers of the terms in the each sum in (801),

$\displaystyle\frac{N!}{p!(q!)^{p}}n_{0},\frac{N!}{p!(q!)^{p}}{p\choose
2}n_{2},\frac{N!}{p!(q!)^{p}}{p\choose
3}n_{3},\dots,\frac{N!}{p!(q!)^{p}}n_{p},$ (807)

and note that the first one is different to the others. And about the last one
the dominant term of the subtrahend should be at the order $(N-q)!$, therefore
the subtrahend is almost at the order $N(N-q)!$. Actually this number is far
larger than the accurate one, since in (805) there’s a factor $(-1)^{r-s}$
then the accurate number is approximately at $(N-q)!$.

Here we also have two conditions similar to (754),(756), and first we consider
the former one. As in this section we assume the mean of the $J_{A_{i}}$ is
zero, therefore we have

$\displaystyle\langle
z^{2}\rangle=\left\langle\sum_{A}\text{sgn}(A)^{2}J_{A_{1}}^{2}\dots
J_{A_{p}}^{2}\right\rangle=\frac{N!}{p!(q!)^{p}}\bar{J}^{2p}$ (808)

where

$\displaystyle\langle J_{A_{i}}^{2}\rangle=\bar{J}^{2}.$ (809)

For the error we evaluate it as

$\displaystyle\langle\text{Error}^{2}\rangle$
$\displaystyle=\langle\left(z^{2}-\Phi\right)^{2}\rangle-\langle
z^{2}\rangle^{2},$ (810)
$\displaystyle=\left\langle\sum_{s=0,2}^{p}\left(\sum_{i_{1}\dots
i_{q+s}}\left(z^{2}_{s}-\Phi_{s}\right)\right)^{2}\right\rangle-\langle
z^{2}\rangle^{2},$ (811)

where $s$ does not take value $1$ for clarity and the cross terms of different
sums are zero due to the zero mean of $J_{A_{i}}$. Then we find that the
computation is similar to before that for $s>0$ the nonzero contribution under
the average is the diagonal terms in the square, while for $s=0$ the whole
terms contribute to the computation of the error. Which means to find out the
dominant sum $s$ we should compare the square of the first term with the rest
terms in (807), and note that the last one cancels out in (811) so that the
non-trivial largest number should be the one labeled by $p-1$. And we give an
upper bound for it

$\displaystyle\frac{N!}{p!(q!)^{p}}{p\choose
p-1}n_{p-1}<p^{2}\frac{N!}{p!(q!)^{p}}\frac{(N-q)!}{(p-1)!(q!)^{p-1}}.$ (812)

Note that the number of the identical terms in this sum may have an effect but
we expect the upper bound is still valid. Then we can find that when $N$,$q$
are sufficiently large it is subordinate to the square of the first term in
(807). Therefore we write the dominant terms in (811)

$\displaystyle\frac{N!}{p!(q!)^{p}}n_{p}-\left(\frac{N!}{p!(q!)^{p}}\right)^{2}$
(813)

where the cross terms of the first part and the whole second part contribute.
Recall the expression (805) the dominant terms in the above equation cancel
out, which implies the proposal holds.

Then we consider the second condition (756) that for a single variable we can
view $J_{A_{i}}^{2}-\Phi_{i}$ as a constant. The expression (811) is still
valid so that the thing we need to consider is to find the correction to
(813). As an example we consider the behavior of the second sum labeled by
$s=2$, the diagonal contribution after the subtraction should be

$\displaystyle\alpha\sum J_{A_{p-1}}J_{A_{p}}J_{B_{p-1}}J_{B_{p}},$ (814)

where $\alpha$ is determined by the quotient of the two numbers of the two sum

$\displaystyle\alpha=\frac{N!}{p!(q!)^{p}}{p\choose 2}n_{2}/\left({N\choose
2q}\frac{(2q)!}{2!(q!)^{2}}n_{2}\right).$ (815)

Therefore after the square the total number becomes

$\displaystyle\left(\frac{N!}{p!(q!)^{p}}{p\choose
2}n_{2}\right)^{2}/\left({N\choose 2q}\frac{(2q)!}{2!(q!)^{2}}n_{2}\right),$
(816)

comparing to the two terms in (813) the above term is subordinate. Similarly
we can find that after the subtraction the largest number in the sums $s>0$
occurs in the case $s=p-1$, which gives

$\displaystyle\left(\frac{N!}{p!(q!)^{p}}{p\choose
p-1}n_{p-1}\right)^{2}/\left({N\choose
N-q}\frac{(N-q)!}{(p-1)!(q!)^{p-1}}n_{p-1}\right).$ (817)

The above contribution is still subordinate to the two terms in (813), which
means the proposal is proper in this case.

Thus for the SYK model we can conclude that when the mean is zero the proposal
works well no matter what the distribution is.

### SYK with one time point with $\mu\neq 0$

Like what we have discussed before the case with non-zero mean is different to
last section. The whole computation is similar to the case with $\mu=0$ except
the mean terms, we can take relevant expressions from the previous sections.
We consider the case with general $N/q$,

$\displaystyle z=\sum_{A}\text{sgn}(A)J_{A_{1}}\dots J_{A_{p}},$ (818)

the main task here is find the dominant terms. The sign functions may have
some effect which makes $s=q-1$ not dominant, such as when $N/q$=2 we have

$\displaystyle\left\langle\sum_{A}\text{sgn}(A)J_{A_{1}}J_{A_{2}}\right\rangle=\sum_{A}\text{sgn}(A)\langle
J_{A_{i}}\rangle^{2}={N/2-1\choose[N/4]}\langle J_{A_{i}}\rangle^{2},$ (819)

where the square bracket means the integer part. Actually here $N/4$ is $q/2$,
in the later computation we always have even $q$ therefore we’ll omit the
square bracket.

To derive this parameter explicitly we have two different ways, one is to find
the sequence functions over different $N$ which can be solved by Mathematica.
Another is to consider different combinations, we can define the positive and
the negative combinations then their difference gives the parameter. Given $N$
different numbers we can have a permutation group which can be divided into
the even and odd parts, and the two parts contain the same number of elements.
If the indices of $J$ have no order then the numbers of the even and odd parts
are equal, so that the parameter is zero. But in the SYK model the indices are
listed in a particular order, it makes the numbers of the two parts different.

To proceed the computation we define the sign of a $q$-sequence,

$\displaystyle\text{sgn}(J_{A_{i}})=(-1)^{a_{i_{1}}+\dots+a_{i_{q}}},\quad
q=0\text{mod}(4)$ (820)

where when $q=2\text{mod}(4)$ the sign will be opposite. Note that here we
take the $a_{i}$ from the indices of $J_{A_{i}}$, but what in the superscript
should be the initial sites of the indices. When the indices take the values
$1,\dots,N$, the values themselves naturally label their initial sites. But
when the indices take arbitrary numbers or symbols, we can define an initial
order and the sites with any permutation. After the assignment complete we can
take the definition (820) to define the sign of any given $q$-sequence.

Then for a combination $A$ containing $r$ negative q-sequences we have the
argument

$\displaystyle\text{sgn}(A)=1,\quad\text{if
}\sum_{i=1}^{r}s_{i}=0\text{mod}(2),$ (821)

where $s_{i}$ is the site of the $i$-th negative $q$-sequence. We give the
condition above as the negative $q$-sequences seems appearing in pairs. For
$N/q=2$ the positive combination $A$ contains two positive q-sequences, so our
mission is to find the first sequence with positive sign. As the combination
has an order the first index in $J_{A_{1}}$ is 1, we only need to consider the
others. The positive $q$-sequence has the number

$\displaystyle{N/2-1\choose 1}{N/2\choose q-2}+{N/2-1\choose 3}{N/2\choose
q-4}+\dots+{N/2-1\choose q-1}{N/2\choose 0},$ (822)

where $q=N/2$ the first chooses odd number of odd integers while the second
chooses even number of even integers. Then the negative $q$-sequence has

$\displaystyle{N/2-1\choose 0}{N/2\choose q-1}+{N/2-1\choose 2}{N/2\choose
q-3}+\dots+{N/2-1\choose q-2}{N/2\choose 1},$ (823)

the difference of the two numbers will give the parameter

$\displaystyle d(p,q)=\sum_{i=0}^{q-1}(-1)^{i+1}{pq/2-1\choose i}{pq/2\choose
q-1-i},$ (824)

which coincides with (819) when $N/q=2$. For general $p=N/q$ we can also find
an expression for finding the first positive $q$-sequence

$\displaystyle d(p,q)=\frac{1}{p-1}{pq/2-1\choose q/2},$ (825)

where $q$ is even and it can be verified by numerics.

And note that for different $q$’s there may be a difference with the minus
sign, which is explained in (820).

For general $N/q$ we can try to compute $\gamma$ in the expression below by
the recursion

$\displaystyle\left\langle\sum_{A}\text{sgn}(A)J_{A_{1}}\dots
J_{A_{p}}\right\rangle=\sum_{A}\text{sgn}(A)\langle
J_{A_{i}}\rangle^{p}=\gamma\langle J_{A_{i}}\rangle^{p}.$ (826)

In the previous computation we have calculated the case with $p=2$, we only
need to derive the $p+1$ case from the $p$ one. Given $N$ numbers we can give
the difference between the positive combinations and the negative ones which
is denoted as $D(p,q)$ and we have

$\displaystyle D(2,q)=d(2,q),\quad D(p+1)=f(D(p,q)).$ (827)

where $f$ is an unknown function. When we add extra $q$ numbers to the $N$
case which is already known, in the new permutation we can fix one number in
the extra $q$ ones at the $p+1$ site in case of repetition. For simplicity we
let the additional $q$ numbers be $1,\dots,q$ and put it at the first site,
and we always fix the number $1$ in this site.

We determine the first site through the procedure in the case of $p=2$, which
can be divided into the positive and negative parts. After choosing $q$
numbers from the $N+q$ ones, about the left $N$ numbers we assume we already
know the difference between the positive combinations and the negative ones.
To illustrate the product of the sequences, we define the positive (negative)
sequences as $P_{+}$ ($P_{-}$) so that we have

$\displaystyle P_{+}P_{+}=P_{+},\quad P_{-}P_{-}=P_{+},\quad
P_{+}P_{-}=P_{-}.$ (828)

The product of two positive or two negative sequences gives a positive one
while the product of one positive and one negative sequences gives a negative
one. Then we can derive the sequences in the case $N+q$ as

$\displaystyle P_{+}^{N+q}=P^{q}_{+}P^{N}_{+}+P^{q}_{-}P^{N}_{-},\quad
P_{-}^{N+q}=P^{q}_{+}P^{N}_{-}+P^{q}_{+}P^{N}_{-},$ (829)

therefore the difference of the positive and negative sequences can be
illustrated as

$\displaystyle
P_{+}^{N+q}-P_{-}^{N+q}=(P_{+}^{q}-P_{-}^{q})(P_{+}^{N}-P_{-}^{N}).$ (830)

Which means the recursion can be expressed as

$\displaystyle D(p,q)=d(p,q)D(p-1,q)=\prod_{s=2}^{p}d(s,q),$ (831)

inserting the equation (825) it becomes

$\displaystyle
D(p,q)=\frac{(pq/2-1)!}{(p-1)!(q/2-1)!((q/2)!)^{p-1}}=\frac{(pq/2)!}{p!((q/2)!)^{p}}.$
(832)

Then we can try to evaluate the cross terms in the square of the partition
function (801). Consider the terms with $p-s$ squares in the sequence,

$\displaystyle\sum_{A,B}\text{sgn}(A)\text{sgn}(B)J_{A_{1}}^{2}\dots
J_{A_{p-s}}^{2}J_{A_{p-s+1}}\dots J_{A_{p}}J_{B_{p-s+1}}\dots J_{B_{p}},$
(833)

after the average over $J$ the number in the sum becomes

$\displaystyle d_{s}=\sum_{r=2}^{s}{p-r\choose s-r}{N\choose pq-rq}\frac{(pq-
rq)!}{(p-r)!(q!)^{p-r}}D(r,q)^{2}(-1)^{s-r}+{p\choose
s}(s-1)(-1)^{s-1}\frac{(pq)!}{p!(q!)^{p}},$ (834)

where the derivation is similar to the equation (805).

We can find an upper bound for the number (834)

$\displaystyle{N\choose pq-sq}\frac{(pq-sq)!}{(p-s)!(q!)^{p-s}}D(s,q)^{2},$
(835)

Then for $s=2,\dots,p-1$ we can compare the above number to the first one in
(807), from this we may find out the dominant one.

By a few numerical comparisons it seems that the dominant terms do not locate
at a particular $s$ as the two parameters $p$,$q$ vary. In this computation we
can directly take the equation (834) rather than the one with approximation.

##  Appendix H Some details about the fermion integral

The relevant integral is

$\displaystyle I_{f}^{N}$ $\displaystyle=$
$\displaystyle\int\text{d}^{2N}\psi\exp\left(T[\sum_{i}[\Sigma_{LR}\psi_{i}^{L}\psi_{i}^{R}+\Sigma_{L}\theta_{i}\psi_{i}^{L}+\Sigma_{R}\theta_{i}\psi_{i}^{R}]\right)\,$
(836) $\displaystyle=$
$\displaystyle\left(\int\text{d}\psi^{L}\text{d}\psi^{R}\exp(T\Sigma_{LR}\psi^{L}\psi^{R})\exp(T\Sigma_{L}\theta\psi^{L})\exp(T\Sigma_{R}\theta\psi^{R})\right)^{N}$
(837)

Before evaluating this integral we need some specifications. Here
$\psi^{L(R)}$ are not Grassmann numbers but can be transferred into Dirac
fermions as

$\displaystyle c=\frac{\psi^{L}+\text{i}\psi^{R}}{2},\quad
c^{\dagger}=\frac{\psi^{L}-\text{i}\psi^{R}}{2},\quad\psi_{L}=c+c^{\dagger},\quad\psi_{R}=\text{i}(c^{\dagger}-c)\,$
(838)
$\displaystyle\\{c,c^{\dagger}\\}=1,\quad\\{c,c\\}=\\{c^{\dagger},c^{\dagger}\\}=0\,.$
(839)

Then the integral $I_{f}$ becomes

$\displaystyle I_{f}$ $\displaystyle=$ $\displaystyle
2\int\text{d}c\text{d}c^{\dagger}e^{a_{1}(cc^{\dagger}-c^{\dagger}c)+a_{-}c+a_{+}c^{\dagger}},\quad
a_{1}=\text{i}T\Sigma_{LR},\quad
a_{\pm}=T(\Sigma_{L}\pm\text{i}\Sigma_{R})\theta,$ (840) $\displaystyle=$
$\displaystyle
4\cosh(\sqrt{a_{1}^{2}+a_{+}a_{-}})=4\cosh(T\sqrt{\Sigma_{L}^{2}+\Sigma_{R}^{2}-\Sigma_{LR}^{2}})$
(841)

Next let us use the same method to compute (382):

$\displaystyle I_{4}^{N}$ $\displaystyle=$
$\displaystyle\int\text{d}^{4N}\psi\exp\left(\Sigma_{ab}\psi_{i}^{a}\psi_{i}^{b}\right)$
(843) $\displaystyle=$
$\displaystyle\\{\int\text{d}^{4}\psi\exp\left(\frac{1}{2}\Sigma_{12}\psi^{1}\psi^{2}+\Sigma_{13}\psi^{1}\psi^{3}+\Sigma_{14}\psi^{1}\psi^{4}\right.$
$\displaystyle\quad\left.\Sigma_{23}\psi^{2}\psi^{3}+\Sigma_{24}\psi^{2}\psi^{4}+\Sigma_{34}\psi^{3}\psi^{4}\right)\\}^{N}.$

Generally to compute $I_{4}$ let us introduce two Dirac fermions

$\displaystyle\psi_{1}=\sqrt{h}(c_{1}+c_{1}^{\dagger}),\quad\psi_{2}=\text{i}\sqrt{h}(c_{1}-c_{1}^{\dagger}),\quad\psi_{3}=\sqrt{h}(c_{2}+c_{2}^{\dagger}),\quad\psi_{4}=\text{i}\sqrt{h}(c_{2}-c_{2}^{\dagger}),$
(844)

such that $\Sigma_{ab}\psi_{i}^{a}\psi_{i}^{b}$ can be written as a $4\times
4$ matrix $H$ with non-vanishing elements

$\displaystyle H_{11}=-\text{i}h(\Sigma_{12}+\Sigma_{34}),\quad
H_{14}=-h(\Sigma_{13}-\Sigma_{24})-\text{i}h(\Sigma_{14}+\Sigma_{23})$ (845)
$\displaystyle H_{22}=\text{i}h(\Sigma_{34}-\Sigma_{12}),\quad
H_{23}=-h(\Sigma_{13}+\Sigma_{24})+\text{i}h(\Sigma_{14}-\Sigma_{23})$ (846)
$\displaystyle
H_{32}=h(\Sigma_{13}+\Sigma_{24})+\text{i}h(\Sigma_{14}-\Sigma_{23}),\quad
H_{33}=-\text{i}h(\Sigma_{34}-\Sigma_{12})$ (847) $\displaystyle
H_{41}=h(\Sigma_{13}-\Sigma_{24})-\text{i}h(\Sigma_{14}+\Sigma_{23}),\quad
H_{44}=\text{i}h(\Sigma_{12}+\Sigma_{34}).$ (848)

By diagonalizing this matrix we can compute $I_{4}$ as the trace of the
exponential of the matrix

$\displaystyle
I_{4}=2\left(\cos\left(h\sqrt{(\Sigma_{14}-\Sigma_{23})^{2}+(\Sigma_{13}+\Sigma_{24})^{2}+(\Sigma_{12}-\Sigma_{34})^{2}}\right)\right.$
$\displaystyle\qquad\left.+\cos\left(h\sqrt{(\Sigma_{14}+\Sigma_{23})^{2}+(\Sigma_{13}-\Sigma_{24})^{2}+(\Sigma_{12}+\Sigma_{34})^{2}}\right)\right).$
(849)

And note that the change of the variable from $\psi$ to $c$,$c^{\dagger}$ will
introduce an additional coefficient in the integral

$\displaystyle(2h)^{2},$ (850)

the final result is the product of the two terms.

## References

  * (1) J. M. Maldacena, “The Large N limit of superconformal field theories and supergravity,” Adv. Theor. Math. Phys. 2, 231-252 (1998) doi:10.1023/A:1026654312961 [arXiv:hep-th/9711200 [hep-th]].
  * (2) E. Witten, “Anti-de Sitter space and holography,” Adv. Theor. Math. Phys. 2, 253-291 (1998) doi:10.4310/ATMP.1998.v2.n2.a2 [arXiv:hep-th/9802150 [hep-th]].
  * (3) S. S. Gubser, I. R. Klebanov and A. M. Polyakov, “Gauge theory correlators from noncritical string theory,” Phys. Lett. B 428, 105-114 (1998) doi:10.1016/S0370-2693(98)00377-3 [arXiv:hep-th/9802109 [hep-th]].
  * (4) G. Penington, “Entanglement Wedge Reconstruction and the Information Paradox,” JHEP 09, 002 (2020) doi:10.1007/JHEP09(2020)002 [arXiv:1905.08255 [hep-th]].
  * (5) A. Almheiri, N. Engelhardt, D. Marolf and H. Maxfield, “The entropy of bulk quantum fields and the entanglement wedge of an evaporating black hole,” JHEP 12, 063 (2019) doi:10.1007/JHEP12(2019)063 [arXiv:1905.08762 [hep-th]].
  * (6) A. Almheiri, R. Mahajan, J. Maldacena and Y. Zhao, “The Page curve of Hawking radiation from semiclassical geometry,” JHEP 03, 149 (2020) doi:10.1007/JHEP03(2020)149 [arXiv:1908.10996 [hep-th]].
  * (7) G. Penington, S. H. Shenker, D. Stanford and Z. Yang, “Replica wormholes and the black hole interior,” [arXiv:1911.11977 [hep-th]].
  * (8) S. R. Coleman, “Black Holes as Red Herrings: Topological Fluctuations and the Loss of Quantum Coherence,” Nucl. Phys. B 307, 867-882 (1988) doi:10.1016/0550-3213(88)90110-1
  * (9) S. B. Giddings and A. Strominger, “Baby Universes, Third Quantization and the Cosmological Constant,” Nucl. Phys. B 321, 481-508 (1989) doi:10.1016/0550-3213(89)90353-2
  * (10) S. B. Giddings and A. Strominger, “Loss of Incoherence and Determination of Coupling Constants in Quantum Gravity,” Nucl. Phys. B 307, 854-866 (1988) doi:10.1016/0550-3213(88)90109-5
  * (11) J. Polchinski and A. Strominger, “A Possible resolution of the black hole information puzzle,” Phys. Rev. D 50, 7403-7409 (1994) doi:10.1103/PhysRevD.50.7403 [arXiv:hep-th/9407008 [hep-th]].
  * (12) P. Saad, S. H. Shenker and D. Stanford, “JT gravity as a matrix integral,” [arXiv:1903.11115 [hep-th]].
  * (13) P. Saad, S. H. Shenker and D. Stanford, “A semiclassical ramp in SYK and in gravity,” [arXiv:1806.06840 [hep-th]].
  * (14) P. Saad, “Late Time Correlation Functions, Baby Universes, and ETH in JT Gravity,” [arXiv:1910.10311 [hep-th]].
  * (15) C. Yan, “Crosscap Contribution to Late-Time Two-Point Correlators,” [arXiv:2203.14436 [hep-th]].
  * (16) R. Bousso, X. Dong, N. Engelhardt, T. Faulkner, T. Hartman, S. H. Shenker and D. Stanford, “Snowmass White Paper: Quantum Aspects of Black Holes and the Emergence of Spacetime,” [arXiv:2201.03096 [hep-th]].
  * (17) J. M. Maldacena and L. Maoz, “Wormholes in AdS,” JHEP 02, 053 (2004) doi:10.1088/1126-6708/2004/02/053 [arXiv:hep-th/0401024 [hep-th]].
  * (18) D. Stanford and E. Witten, “JT gravity and the ensembles of random matrix theory,” Adv. Theor. Math. Phys. 24, no.6, 1475-1680 (2020) doi:10.4310/ATMP.2020.v24.n6.a4 [arXiv:1907.03363 [hep-th]].
  * (19) L. V. Iliesiu, “On 2D gauge theories in Jackiw-Teitelboim gravity,” [arXiv:1909.05253 [hep-th]].
  * (20) D. Kapec, R. Mahajan and D. Stanford, “Matrix ensembles with global symmetries and ’t Hooft anomalies from 2d gauge theory,” JHEP 04, 186 (2020) doi:10.1007/JHEP04(2020)186 [arXiv:1912.12285 [hep-th]].
  * (21) H. Maxfield and G. J. Turiaci, “The path integral of 3D gravity near extremality; or, JT gravity with defects as a matrix integral,” JHEP 01, 118 (2021) doi:10.1007/JHEP01(2021)118 [arXiv:2006.11317 [hep-th]].
  * (22) E. Witten, “Matrix Models and Deformations of JT Gravity,” Proc. Roy. Soc. Lond. A 476, no.2244, 20200582 (2020) doi:10.1098/rspa.2020.0582 [arXiv:2006.13414 [hep-th]].
  * (23) E. Mefford and K. Suzuki, “Jackiw-Teitelboim quantum gravity with defects and the Aharonov-Bohm effect,” JHEP 05, 026 (2021) doi:10.1007/JHEP05(2021)026 [arXiv:2011.04695 [hep-th]].
  * (24) A. Altland and J. Sonner, “Late time physics of holographic quantum chaos,” SciPost Phys. 11, 034 (2021) doi:10.21468/SciPostPhys.11.2.034 [arXiv:2008.02271 [hep-th]].
  * (25) L. Eberhardt, “Summing over Geometries in String Theory,” JHEP 05, 233 (2021) doi:10.1007/JHEP05(2021)233 [arXiv:2102.12355 [hep-th]].
  * (26) D. Stanford, Z. Yang and S. Yao, “Subleading Weingartens,” JHEP 02, 200 (2022) doi:10.1007/JHEP02(2022)200 [arXiv:2107.10252 [hep-th]].
  * (27) I. Aref’eva and I. Volovich, “Gas of Baby Universes in JT Gravity and Matrix Models,” Symmetry 12, no.6, 975 (2020) doi:10.3390/sym12060975 [arXiv:1905.08207 [hep-th]].
  * (28) P. Betzios and O. Papadoulaki, “Liouville theory and Matrix models: A Wheeler DeWitt perspective,” JHEP 09, 125 (2020) doi:10.1007/JHEP09(2020)125 [arXiv:2004.00002 [hep-th]].
  * (29) D. Anninos and B. Mühlmann, “Notes on matrix models (matrix musings),” J. Stat. Mech. 2008, 083109 (2020) doi:10.1088/1742-5468/aba499 [arXiv:2004.01171 [hep-th]].
  * (30) M. Berkooz, V. Narovlansky and H. Raj, “Complex Sachdev-Ye-Kitaev model in the double scaling limit,” JHEP 02, 113 (2021) doi:10.1007/JHEP02(2021)113 [arXiv:2006.13983 [hep-th]].
  * (31) T. G. Mertens and G. J. Turiaci, “Liouville quantum gravity – holography, JT and matrices,” JHEP 01, 073 (2021) doi:10.1007/JHEP01(2021)073 [arXiv:2006.07072 [hep-th]].
  * (32) G. J. Turiaci, M. Usatyuk and W. W. Weng, “Dilaton-gravity, deformations of the minimal string, and matrix models,” doi:10.1088/1361-6382/ac25df [arXiv:2011.06038 [hep-th]].
  * (33) D. Anninos and B. Mühlmann, “Matrix integrals & finite holography,” JHEP 06, 120 (2021) doi:10.1007/JHEP06(2021)120 [arXiv:2012.05224 [hep-th]].
  * (34) P. Gao, D. L. Jafferis and D. K. Kolchmeyer, “An effective matrix model for dynamical end of the world branes in Jackiw-Teitelboim gravity,” [arXiv:2104.01184 [hep-th]].
  * (35) V. Godet and C. Marteau, “From black holes to baby universes in CGHS gravity,” JHEP 07, 138 (2021) doi:10.1007/JHEP07(2021)138 [arXiv:2103.13422 [hep-th]].
  * (36) C. V. Johnson, F. Rosso and A. Svesko, “Jackiw-Teitelboim supergravity as a double-cut matrix model,” Phys. Rev. D 104, no.8, 086019 (2021) doi:10.1103/PhysRevD.104.086019 [arXiv:2102.02227 [hep-th]].
  * (37) A. Blommaert and M. Usatyuk, “Microstructure in matrix elements,” [arXiv:2108.02210 [hep-th]].
  * (38) K. Okuyama and K. Sakai, “JT gravity, KdV equations and macroscopic loop operators,” JHEP 01, 156 (2020) doi:10.1007/JHEP01(2020)156 [arXiv:1911.01659 [hep-th]].
  * (39) S. Forste, H. Jockers, J. Kames-King and A. Kanargias, “Deformations of JT Gravity via Topological Gravity and Applications,” [arXiv:2107.02773 [hep-th]].
  * (40) A. Maloney and E. Witten, “Averaging over Narain moduli space,” JHEP 10, 187 (2020) doi:10.1007/JHEP10(2020)187 [arXiv:2006.04855 [hep-th]].
  * (41) N. Afkhami-Jeddi, H. Cohn, T. Hartman and A. Tajdini, “Free partition functions and an averaged holographic duality,” JHEP 01, 130 (2021) doi:10.1007/JHEP01(2021)130 [arXiv:2006.04839 [hep-th]].
  * (42) J. Cotler and K. Jensen, “AdS3 gravity and random CFT,” JHEP 04, 033 (2021) doi:10.1007/JHEP04(2021)033 [arXiv:2006.08648 [hep-th]].
  * (43) A. Pérez and R. Troncoso, “Gravitational dual of averaged free CFT’s over the Narain lattice,” JHEP 11, 015 (2020) doi:10.1007/JHEP11(2020)015 [arXiv:2006.08216 [hep-th]].
  * (44) N. Benjamin, C. A. Keller, H. Ooguri and I. G. Zadeh, “Narain to Narnia,” Commun. Math. Phys. 390, no.1, 425-470 (2022) doi:10.1007/s00220-021-04211-x [arXiv:2103.15826 [hep-th]].
  * (45) J. Cotler and K. Jensen, “AdS3 wormholes from a modular bootstrap,” JHEP 11, 058 (2020) doi:10.1007/JHEP11(2020)058 [arXiv:2007.15653 [hep-th]].
  * (46) M. Ashwinkumar, M. Dodelson, A. Kidambi, J. M. Leedom and M. Yamazaki, “Chern-Simons Invariants from Ensemble Averages,” doi:10.1007/JHEP08(2021)044 [arXiv:2104.14710 [hep-th]].
  * (47) N. Afkhami-Jeddi, A. Ashmore and C. Cordova, “Calabi-Yau CFTs and Random Matrices,” [arXiv:2107.11461 [hep-th]].
  * (48) S. Collier and A. Maloney, “Wormholes and Spectral Statistics in the Narain Ensemble,” [arXiv:2106.12760 [hep-th]].
  * (49) N. Benjamin, S. Collier, A. L. Fitzpatrick, A. Maloney and E. Perlmutter, “Harmonic analysis of 2d CFT partition functions,” doi:10.1007/JHEP09(2021)174 [arXiv:2107.10744 [hep-th]].
  * (50) J. Dong, T. Hartman and Y. Jiang, “Averaging over moduli in deformed WZW models,” doi:10.1007/JHEP09(2021)185 [arXiv:2105.12594 [hep-th]].
  * (51) A. Dymarsky and A. Shapere, “Comments on the holographic description of Narain theories,” JHEP 10, 197 (2021) doi:10.1007/JHEP10(2021)197 [arXiv:2012.15830 [hep-th]].
  * (52) V. Meruliya, S. Mukhi and P. Singh, “Poincaré Series, 3d Gravity and Averages of Rational CFT,” JHEP 04, 267 (2021) doi:10.1007/JHEP04(2021)267 [arXiv:2102.03136 [hep-th]].
  * (53) R. Bousso and E. Wildenhain, “Gravity/ensemble duality,” Phys. Rev. D 102, no.6, 066005 (2020) doi:10.1103/PhysRevD.102.066005 [arXiv:2006.16289 [hep-th]].
  * (54) O. Janssen, M. Mirbabayi and P. Zograf, “Gravity as an ensemble and the moment problem,” JHEP 06, 184 (2021) doi:10.1007/JHEP06(2021)184 [arXiv:2103.12078 [hep-th]].
  * (55) J. Cotler and K. Jensen, “Wormholes and black hole microstates in AdS/CFT,” JHEP 09, 001 (2021) doi:10.1007/JHEP09(2021)001 [arXiv:2104.00601 [hep-th]].
  * (56) D. Marolf and H. Maxfield, “Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information,” JHEP 08, 044 (2020) doi:10.1007/JHEP08(2020)044 [arXiv:2002.08950 [hep-th]].
  * (57) V. Balasubramanian, A. Kar, S. F. Ross and T. Ugajin, “Spin structures and baby universes,” JHEP 09, 192 (2020) doi:10.1007/JHEP09(2020)192 [arXiv:2007.04333 [hep-th]].
  * (58) J. G. Gardiner and S. Megas, “2d TQFTs and baby universes,” JHEP 10, 052 (2021) doi:10.1007/JHEP10(2021)052 [arXiv:2011.06137 [hep-th]].
  * (59) A. Belin and J. de Boer, “Random statistics of OPE coefficients and Euclidean wormholes,” Class. Quant. Grav. 38, no.16, 164001 (2021) doi:10.1088/1361-6382/ac1082 [arXiv:2006.05499 [hep-th]].
  * (60) A. Belin, J. De Boer, P. Nayak and J. Sonner, “Charged Eigenstate Thermalization, Euclidean Wormholes and Global Symmetries in Quantum Gravity,” [arXiv:2012.07875 [hep-th]].
  * (61) A. Altland, D. Bagrets, P. Nayak, J. Sonner and M. Vielma, “From operator statistics to wormholes,” Phys. Rev. Res. 3, no.3, 033259 (2021) doi:10.1103/PhysRevResearch.3.033259 [arXiv:2105.12129 [hep-th]].
  * (62) A. Belin, J. de Boer, P. Nayak and J. Sonner, “Generalized Spectral Form Factors and the Statistics of Heavy Operators,” [arXiv:2111.06373 [hep-th]].
  * (63) C. Peng, J. Tian and J. Yu, “Baby universes, ensemble averages and factorizations with matters,” [arXiv:2111.14856 [hep-th]].
  * (64) A. Banerjee and G. W. Moore, “Comments on Summing over bordisms in TQFT,” [arXiv:2201.00903 [hep-th]].
  * (65) J. J. Heckman, A. P. Turner and X. Yu, “Disorder averaging and its UV discontents,” Phys. Rev. D 105, no.8, 086021 (2022) doi:10.1103/PhysRevD.105.086021 [arXiv:2111.06404 [hep-th]].
  * (66) C. V. Johnson, “The Microstate Physics of JT Gravity and Supergravity,” [arXiv:2201.11942 [hep-th]].
  * (67) S. Collier and E. Perlmutter, “Harnessing S-Duality in $\mathcal{N}=4$ SYM & Supergravity as $SL(2,\mathbb{Z})$-Averaged Strings,” [arXiv:2201.05093 [hep-th]].
  * (68) J. Chandra, S. Collier, T. Hartman and A. Maloney, “Semiclassical 3D gravity as an average of large-c CFTs,” [arXiv:2203.06511 [hep-th]].
  * (69) J. M. Schlenker and E. Witten, “No Ensemble Averaging Below the Black Hole Threshold,” [arXiv:2202.01372 [hep-th]].
  * (70) R. Jackiw, “Lower Dimensional Gravity,” Nucl. Phys. B 252, 343-356 (1985) doi:10.1016/0550-3213(85)90448-1
  * (71) C. Teitelboim, “Gravitation and Hamiltonian Structure in Two Space-Time Dimensions,” Phys. Lett. B 126, 41-45 (1983) doi:10.1016/0370-2693(83)90012-6
  * (72) S. Sachdev and J. Ye, “Gapless spin fluid ground state in a random, quantum Heisenberg magnet,” Phys. Rev. Lett. 70, 3339 (1993) doi:10.1103/PhysRevLett.70.3339 [arXiv:cond-mat/9212030 [cond-mat]].
  * (73) A. Kitaev, ”A simple model of quantum holography.” Talks at KITP http://online.kitp.ucsb.edu/online/entangled15/kitaev/ and http://online.kitp.ucsb.edu/online/entangled15/kitaev2
A. Kitaev, ”Hidden correlations in the Hawking radiation and thermal noise.”
Talk at KITP http://online.kitp.ucsb.edu/online/joint98/kitaev/

  * (74) D. Stanford, “More quantum noise from wormholes,” [arXiv:2008.08570 [hep-th]].
  * (75) A. Almheiri and H. W. Lin, “The Entanglement Wedge of Unknown Couplings,” [arXiv:2111.06298 [hep-th]].
  * (76) P. Saad, S. Shenker and S. Yao, “Comments on wormholes and factorization,” [arXiv:2107.13130 [hep-th]].
  * (77) P. Saad, S. H. Shenker, D. Stanford and S. Yao, “Wormholes without averaging,” [arXiv:2103.16754 [hep-th]].
  * (78) B. Mukhametzhanov, “Half-wormholes in SYK with one time point,” [arXiv:2105.08207 [hep-th]].
  * (79) A. M. García-García and V. Godet, “Half-wormholes in nearly AdS2 holography,” [arXiv:2107.07720 [hep-th]].
  * (80) S. Choudhury and K. Shirish, “Wormhole calculus without averaging from $O(N)^{q-1}$ tensor model,” [arXiv:2106.14886 [hep-th]].
  * (81) B. Mukhametzhanov, “Factorization and complex couplings in SYK and in Matrix Models,” [arXiv:2110.06221 [hep-th]].
  * (82) K. Okuyama and K. Sakai, “FZZT branes in JT gravity and topological gravity,” JHEP 09, 191 (2021) doi:10.1007/JHEP09(2021)191 [arXiv:2108.03876 [hep-th]].
  * (83) K. Goto, Y. Kusuki, K. Tamaoka and T. Ugajin, “Product of random states and spatial (half-)wormholes,” JHEP 10, 205 (2021) doi:10.1007/JHEP10(2021)205 [arXiv:2108.08308 [hep-th]].
  * (84) A. Blommaert, L. V. Iliesiu and J. Kruthoff, “Gravity factorized,” [arXiv:2111.07863 [hep-th]].
  * (85) K. Goto, K. Suzuki and T. Ugajin, “Factorizing Wormholes in a Partially Disorder-Averaged SYK Model,” [arXiv:2111.11705 [hep-th]].
  * (86) A.I. Barvinok, _New algorithms for linear k-matroid intersection and matroid k-parity problems_ , _Mathematical Programming_ 69 (1995) 449A470.
  * (87) C. Peng, “Ensemble averages, Poisson processes, and microstates,” Phys. Rev. D 103, no.6, L061901 (2021) doi:10.1103/PhysRevD.103.L061901 [arXiv:2010.11192 [hep-th]].
  * (88) N. N. Bogolyubov, “On the theory of superfluidity,” J. Phys. (USSR) 11, 23-32 (1947) Valatin:1958ja
  * (89) J. G. Valatin, “Comments on the theory of superconductivity,” Nuovo Cim. 7, 843-857 (1958) doi:10.1007/BF02745589
  * (90) E. Witten, “Analytic Continuation Of Chern-Simons Theory,” AMS/IP Stud. Adv. Math. 50, 347-446 (2011) [arXiv:1001.2933 [hep-th]].
  * (91) A. Blommaert and J. Kruthoff, “Gravity without averaging,” [arXiv:2107.02178 [hep-th]].
  * (92) A. Blommaert, T. G. Mertens and H. Verschelde, “Eigenbranes in Jackiw-Teitelboim gravity,” JHEP 02, 168 (2021) doi:10.1007/JHEP02(2021)168 [arXiv:1911.11603 [hep-th]].
  * (93) A. Goel, L. V. Iliesiu, J. Kruthoff and Z. Yang, “Classifying boundary conditions in JT gravity: from energy-branes to $\alpha$-branes,” JHEP 04, 069 (2021) doi:10.1007/JHEP04(2021)069 [arXiv:2010.12592 [hep-th]].
  * (94) A. Blommaert, L. V. Iliesiu and J. Kruthoff, “Alpha states demystified: Towards microscopic models of AdS2 holography,” [arXiv:2203.07384 [hep-th]].
  * (95) D. Das, S. Pal and A. Sarkar, “(Half) Wormholes under Irrelevant Deformation,” [arXiv:2203.14988 [hep-th]].


# Variations on a result of Granville about polynomial invariance

Labib Haddad and Charles Helou  120 Rue de Charonne, 75011 Paris, France;
email: labib.haddadwanadoo.fr  Penn State Univ., 25 Yearsley Mill Rd, Media,
PA 19063, USA; e-mail: cxh22psu.edu

## Abstract

In 2004, Andrew Granville proved the following: A polynomial in variables
$x,y,z,$ invariant under the map $\varphi:z\to-(x+y+z)$ is a polynomial in
$b=z(x+y+z)$. We tell, curtly, its history and then give a variation on this
theme.

## History

In March 2004, the second author received an inquiry from a colleague about
the factorization of some symmetric polynomials, with power monomials, in two
or three variables. There was a statement about a factorization in
$\mathbb{Z}[x,y]$ of the polynomial $F_{p}(x,y)=(x+y)^{p}-x^{p}-y^{p}$, for a
prime number $p>3$, namely

$F_{p}(x,y)=pxy(x+y)\left(x^{2}+xy+y^{2}\right)^{e}C_{p}(x,y),$

where $e=1\ \text{or}\ 2$ according to whether $p\equiv-1\ \text{or}\
1\pmod{6}$, respectively, and $C_{p}(x,y)\in\mathbb{Z}[x,y]$ is a so-called
Cauchy polynomial, described as a well-known elementary result. In fact,
$F_{p}(x,y)$ is a homogeneous polynomial of degree $p$, which, upon division
by $y^{p}$, gives a polynomial $K_{p}(t)$, of degree $p-1$ in $t=\frac{x}{y}$,
namely

$\frac{F_{p}(x,y)}{y^{p}}=\left(\frac{x}{y}+1\right)^{p}-\left(\frac{x}{y}\right)^{p}-1=(t+1)^{p}-t^{p}-1=\sum_{k=1}^{p-1}\binom{p}{k}t^{k}=K_{p}(t).$

Clearly, $K_{p}(0)=K_{p}(-1)=0$ and $K^{\prime}_{p}(0)=p,\
K^{\prime}_{p}(-1)=-p\neq 0$, so $0$ and $1$ are simple roots of $K_{p}$, so
that $K_{p}$ is divisible by $t(t+1)$ in $\mathbb{Z}[t]$, and also by $p$
since the binomial coefficients $\binom{p}{k}$ are divisible by $p$ for $1\leq
k\leq p-1$. Moreover, letting $\zeta_{3}=\exp\left(\frac{2i\pi}{3}\right)$ be
a primitive cubic root of unity, we have

$K_{p}(\zeta_{3})=(\zeta_{3}+1)^{p}-\zeta_{3}^{p}-1=-(1+\zeta_{3}^{p}+\zeta_{3}^{2p})=0$

and

$K^{\prime}_{p}(\zeta_{3})=p((\zeta_{3}+1)^{p-1}-\zeta_{3}^{p-1})=p\zeta_{3}^{p-1}(\zeta_{3}^{p-1}-1),$

so that $K^{\prime}_{p}(\zeta_{3})=0$ if and only if $p\equiv 1\pmod{3}$,
i.e., $p\equiv 1\pmod{6}$. Thus $\zeta_{3}$ is a root of $K_{p}(t)$, and it’s
a double root if and only if $p\equiv 1\pmod{6}$, i.e., $K_{p}(t)$ is
divisible by $t^{2}+t+1$ if $p\equiv-1\pmod{6}$ and by $(t^{2}+t+1)^{2}$ if
$p\equiv 1\pmod{6}$, in $\mathbb{Z}[t]$. Thus

$K_{p}(t)=pt(t+1)(t^{2}+t+1)^{e}C_{p}(t),$

where

$e=\begin{cases}1,&\text{if $p\equiv-1\pmod{6}$}\\\ 2,&\text{if $p\equiv
1\pmod{6}$},\end{cases}$

and $C_{p}(t)\in\mathbb{Z}[t]$, of degree $p-3-2e$, is called a Cauchy
polynomial, because it was introduced by A. L. Cauchy in 1839 and reconsidered
in 1841 ([2,  3,  5]). The Cauchy polynomials were studied by D. Mirimanoff in
1903 ([4]), who conjectured their irreducibility over $\mathbb{Q}$, a still
open conjecture. For more details, one can consult the 1989 paper of G.
Terjanian ([6]). This factorization of $K_{p}(t)$ implies the factorization of
$F_{p}(x,y)$ cited above. To be accurate, $C_{p}(t)$ should be $C_{p}(t,1)$ in
the factorization of $K_{p}(t)$, but by “abuse of notation”, we identified the
two.

This leads to introduce the polynomial

$E_{2}=\frac{(x+y)^{p}-x^{p}-y^{p}}{pxy(x+y)}=\left(x^{2}+xy+y^{2}\right)^{e}C_{p}(x,y),$

which lies in $Z[x,y]$, and, by analogy, in the case of 3 variables $x,y,z$,
the polynomial

$E_{3}=\frac{(x+y+z)^{p}-x^{p}-y^{p}-z^{p}}{p(x+y)(y+z)(z+x)},$

which lies in $Z[x,y,z]$. Indeed, letting

$H_{p}(x,y,z)=(x+y+z)^{p}-x^{p}-y^{p}-z^{p}=\sum\Sb 0\leq j,k,l<p\\\
j+k+l=p\frac{p!}{j!k!l!}x^{j}y^{k}z^{l},$

it is easily seen that all the coefficients of $H_{p}$ are divisible by $p$,
and that $H_{p}(x,-x,z)=H_{p}(x,y,-y)=H_{p}(x,y,-x)=0$, so that $H_{p}$ is
divisible by $p$, $x+y$, $y+z$ and $z+x$, in $Z[x,y,z]$, which are pairwise
coprime. Therefore $H_{p}$ is divisible by their product $p(x+y)(y+z)(z+x)$,
i.e., $E_{3}$ has its coefficients in $\mathbb{Z}$. Moreover, it is noted that

$E_{3}=E_{2}+Q(x,y,z),$

where $Q\in z\mathbb{Z}[x,y,z]$. Indeed,

$\displaystyle
E_{3}-E_{2}=\frac{(x+y+z)^{p}-x^{p}-y^{p}-z^{p}}{p(x+y)(y+z)(z+x)}-\frac{(x+y)^{p}-x^{p}-y^{p}}{pxy(x+y)}=$
$\displaystyle=\frac{xy((x+y+z)^{p}-x^{p}-y^{p}-z^{p})-(y+z)(z+x)((x+y)^{p}-x^{p}-y^{p})}{pxy(x+y)(y+z)(z+x)}=$
$\displaystyle=\frac{xy\left((\sum_{k=1}^{p-1}\binom{p}{k}(x+y)^{k}z^{p-k}+(x+y)^{p}-x^{p}-y^{p}\right)-(xy+xz+yz+z^{2})((x+y)^{p}-x^{p}-y^{p})}{pxy(x+y)(y+z)(z+x)}=$
$\displaystyle=\frac{xy\sum_{k=1}^{p-1}\binom{p}{k}(x+y)^{k}z^{p-k}-z(x+y+z)((x+y)^{p}-x^{p}-y^{p})}{pxy(x+y)(y+z)(z+x)},$

which is a polynomial in $\mathbb{Z}[x,y,z]$ all of whose terms are multiples
of $z$.

Then, came the main question. Setting $b=z(x+y+z)$, my correspondent made the
observation that

$E_{3}=a_{n}+b(a_{n-1}+b(a_{n-2}+\cdots+b(a_{2}+b(a_{1}+b))\cdots))=b^{n}+a_{1}b^{n-1}+\cdots
a_{n-2}b^{2}+a_{n-1}b+a_{n},$

where $n=\frac{p-3}{2}$ and all the $a_{k}\in\mathbb{Z}[x,y]$, with
$a_{n}=E_{2}$. This amounts to stating that $E_{3}$ is a polynomial in $b$
with coefficients in $\mathbb{Z}[x,y]$, and asking for a proof.

In April 2004, a second message came, extending the first one. In it, he makes
the observation that $a_{n}+xya_{n-1}=(x+y)^{p-3}$, and he conjectures that,
for $0\leq m\leq n-1$, where $n=\frac{p-3}{2}$,

$a_{n-m}+xya_{n-m-1}=a_{n-m}^{(1)}(x+y)^{p-2m-3},$

where $a_{k}^{(j)}$ is the coefficient of the $j$-th term in the polynomial
$a_{k}$. Then it is claimed that a result of E. Catalan, cited in L. E.
Dickson’s History of the Theory of Numbers, vol. II, gives the coefficients of
$C_{p}(x,y)=\sum_{k=0}^{p-3-2e}c_{k}x^{p-3-2e-k}y^{k}$, namely (with a
correction)

$c_{k}=\frac{1}{p}\left(\binom{p-1}{k}-(-1)^{k}\right).$

He then claims to deduce that $a_{n-1}^{(1)}=a_{n}^{(2)}$ and
$a_{n}^{(2)}=\frac{p-3}{2}$. All computations were done, using Mathematica.

In May 2004, a third message came saying that he sent his query to Andrew
Granville, who replied that any polynomial $F(x,y,z)$ which is invariant under
the map $\varphi:z\mapsto-(x+y+z)$ (as is $E_{3}$) can be written as a
polynomial in $b=z(x+y+z)$, providing a simple proof by induction on the
degree of $F$. Indeed, setting $a_{0}=F(x,y,0)$, Granville notes that $a_{0}$
and $F(x,y,z)-a_{0}$ are invariant under $\varphi$. But $F(x,y,z)-a_{0}$ is
divisible by $z$, and, applying $\varphi$, it is also divisible by $-(x+y+z)$,
and therefore it is divisible by $b$. Noting that $b$ is also invariant under
$\varphi$, he deduces that $G(x,y,z)=\frac{F(x,y,z)-a_{0}}{b}$ is invariant
under $\phi$, and has lower degree than $F$. Thus, by the induction
hypothesis, $G$ can be written as a polynomial in $b$, and therefore so can
$F$.

It is this brilliantly simple argument of Granville that we attempt to
generalize below.

Remarks. The polynomials denoted by $E_{2}$ and $E_{3}$ will henceforth be
designated by $E_{2,p}$ and $E_{3,p}$ to emphasize their dependence on $p$.
Note that

$\displaystyle\frac{(x+y)^{p}-(x^{p}+y^{p})}{x+y}$
$\displaystyle=(x+y)^{p-1}-\frac{x^{p}+y^{p}}{x+y}=\sum_{k=0}^{p-1}\binom{p-1}{k}x^{k}y^{p-1-k}-\sum_{k=0}^{p-1}(-1)^{k}x^{k}y^{p-1-k}$
$\displaystyle=\sum_{k=0}^{p-1}\left(\binom{p-1}{k}-(-1)^{k}\right)x^{k}y^{p-1-k}=\sum_{k=1}^{p-2}\left(\binom{p-1}{k}-(-1)^{k}\right)x^{k}y^{p-1-k},$

since $\binom{p-1}{k}-(-1)^{k}=0$ for $k=0$ and $k=p-1$. Hence

$\displaystyle E_{2,p}(x,y)$
$\displaystyle=\frac{(x+y)^{p}-(x^{p}+y^{p})}{pxy(x+y)}=\frac{1}{pxy}\sum_{k=1}^{p-2}\left(\binom{p-1}{k}-(-1)^{k}\right)x^{k}y^{p-1-k}=$
$\displaystyle=\sum_{k=1}^{p-2}\frac{\binom{p-1}{k}-(-1)^{k}}{p}x^{k-1}y^{p-k-2},$

where the coefficients $\frac{\binom{p-1}{k}-(-1)^{k}}{p}$ are integers, for
$1\leq k\leq p-2$, since

$\binom{p-1}{k}=\frac{(p-1)(p-2)\cdots(p-k)}{k!}\equiv\frac{(-1)^{k}k!}{k!}\equiv(-1)^{k}\pmod{p}.$

On the other hand,

$\displaystyle H_{p}(x,y,z)$
$\displaystyle=(x+y+z)^{p}-(x^{p}+y^{p}+z^{p})=\sum_{m=0}^{p}\binom{p}{m}(x+y)^{m}z^{p-m}-(x^{p}+y^{p}+z^{p})=$
$\displaystyle=\sum_{m=1}^{p}\binom{p}{m}(x+y)^{m}z^{p-m}-(x^{p}+y^{p}),$

so that

$\displaystyle E_{3,p}(x,y,z)$
$\displaystyle=\frac{H_{p}(x,y,z)}{p(x+y)(y+z)(z+x)}=\frac{1}{p(y+z)(z+x)}\left(\sum_{m=1}^{p}\binom{p}{m}(x+y)^{m-1}z^{p-m}-\frac{x^{p}+y^{p}}{x+y}\right)$
$\displaystyle=\frac{1}{p(y+z)(z+x)}\left(\sum_{m=1}^{p}\binom{p}{m}(x+y)^{m-1}z^{p-m}-\sum_{k=0}^{p-1}(-1)^{k}x^{k}y^{p-1-k}\right)=$
$\displaystyle=\frac{1}{p(y+z)(z+x)}\sum_{k=0}^{p-1}\left(\binom{p}{k+1}(x+y)^{k}z^{p-k-1}-(-1)^{k}x^{k}y^{p-k-1}\right),$

the last equality being obtained by setting $k=m-1$ in the first of the two
sums.

Furthermore, multiplying and dividing by $(x+y)$, distributing in the
numerator, and setting $l=k+1$, we get

$\displaystyle E_{3,p}(x,y,z)$
$\displaystyle=\frac{1}{p(x+y)(y+z)(z+x)}\left(\sum_{l=1}^{p}\binom{p}{l}(x+y)^{l}z^{p-l}+\sum_{l=1}^{p}(-1)^{l}x^{l}y^{p-l}-\sum_{k=0}^{p-1}(-1)^{k}x^{k}y^{p-k}\right)$
$\displaystyle=\frac{1}{p(x+y)(y+z)(z+x)}\left(\sum_{l=0}^{p}\binom{p}{l}(x+y)^{l}z^{p-l}-x^{p}-y^{p}-z^{p})\right).$

Expanding the binomial $(x+y)^{l}$, and replacing all binomial coefficients by
their defining expressions, then reversing the order of summations, and
factoring out the common factors, we get

$E_{3,p}(x,y,z)=\frac{1}{p(x+y)(y+z)(z+x)}\left(p!z^{p}\sum_{k=0}^{p}\frac{1}{k!}\left(\frac{x}{y}\right)^{k}\sum_{l=k}^{p}\frac{1}{(p-l)!(l-k)!}\left(\frac{y}{z}\right)^{l}-x^{p}-y^{p}-z^{p}\right).$

Hence the following two expressions for $E_{3,p}$, established above,

$\displaystyle E_{3,p}(x,y,z)$
$\displaystyle=\frac{1}{p(y+z)(z+x)}\sum_{k=0}^{p-1}\left(\binom{p}{k+1}(x+y)^{k}z^{p-k-1}-(-1)^{k}x^{k}y^{p-k-1}\right)=$
$\displaystyle=\frac{1}{p(x+y)(y+z)(z+x)}\left(p!z^{p}\sum_{k=0}^{p}\frac{1}{k!}\left(\frac{x}{y}\right)^{k}\sum_{l=k}^{p}\frac{1}{(p-l)!(l-k)!}\left(\frac{y}{z}\right)^{l}-x^{p}-y^{p}-z^{p}\right).$

Similarly, exchanging $x+y$ and $y+z$, then $x+y$ and $x+z$, gives,
respectively,

$\displaystyle E_{3,p}(x,y,z)$
$\displaystyle=\frac{1}{p(x+y)(x+z)}\sum_{k=0}^{p-1}\left(\binom{p}{k+1}(y+z)^{k}x^{p-k-1}-(-1)^{k}y^{k}z^{p-k-1}\right)=$
$\displaystyle=\frac{1}{p(x+y)(y+z)(z+x)}\left(p!x^{p}\sum_{k=0}^{p}\frac{1}{k!}\left(\frac{y}{z}\right)^{k}\sum_{l=k}^{p}\frac{1}{(p-l)!(l-k)!}\left(\frac{z}{x}\right)^{l}-x^{p}-y^{p}-z^{p}\right),$

and

$\displaystyle E_{3,p}(x,y,z)$
$\displaystyle=\frac{1}{p(x+y)(y+z)}\sum_{k=0}^{p-1}\left(\binom{p}{k+1}(x+z)^{k}y^{p-k-1}-(-1)^{k}x^{k}z^{p-k-1}\right)=$
$\displaystyle=\frac{1}{p(x+y)(y+z)(z+x)}\left(p!y^{p}\sum_{k=0}^{p}\frac{1}{k!}\left(\frac{z}{x}\right)^{k}\sum_{l=k}^{p}\frac{1}{(p-l)!(l-k)!}\left(\frac{x}{y}\right)^{l}-x^{p}-y^{p}-z^{p}\right).$

Moreover, after a search for the result of E. Catalan mentioned above, we
located it in an article ([1]), originally published in 1861, then taken up
again, and complemented, in 1885, by Catalan, where he gives an explicit
expression for the polynomials

$E_{3,n}(x,y,z)=\frac{(x+y+z)^{n}-x^{n}-y^{n}-z^{n}}{(x+y)(y+z)(z+x)},$

for odd integers $n>3$. He defines the homogeneous polynomials of degree
$m\geq 0$, with all coefficients equal to 1

$H_{m}=H_{m}(x,y,z)=\sum_{j,k,l\in\mathbb{N}:j+k+l=m}x^{j}y^{k}z^{l}=\sum_{j=0}^{m}\sum_{k=0}^{m-j}x^{j}y^{k}z^{m-j-k},$

with $H_{0}=1$, he sets $P=x+y+z$, and he establishes that, for odd $n>3$,

$E_{3,n}(x,y,z)=\sum_{m=0}^{n-3}H_{m}(x,y,z)P^{n-3-m}+2H_{\frac{n-3}{2}}(x^{2},y^{2},z^{2}).$

## The variation

Let $K$ be a commutative field. Take $n+1$ variables $x=(x_{1},...,x_{n})$ and
$z$, then consider the $K[x]$-algebra of polynomials $A=K[x,z]$. Choose a
polynomial, $p(x,z)$, in $A$ and define the following transformation $T$ of
$A$ into $A$: To each polynomial $F(x,z)$, it matches the polynomial

$TF(x,z)=F(x,p(x,z)).$ $None$

This transformation $T$ is a $K[x]$-homomorphism of the algebra $A$ into
itself. Le $B$ be the set of polynomial invariant under $T$, i.e.,

$B=\\{F\in A:TF=F\\}.$ $None$

Of course, $B$ is a $K[x]$-subalgebra of $A$ and it contains $K[x]$. Then set

$C=B\setminus K[x].$ $None$

How can the elements of $B$ be characterized?

What is needed for $C$ not to be empty?

The case $p(x,z)=z$ is trivial since $T=$ idA and $B=A$.

We shall say that polynomial $p(x,z)$ is à la Granville when

$p(x,z)=-z+r(x).$ $None$

Notice that $T$ is an involution (i.e., $T^{2}=$ idA) in that case since we
will have $Tp(x,z)=z$.

The case studied by Granville is that in which $p(x,z)=-(x+y+z)$.

Now, on to the “variation”!

### 1 Preparation

Introduce the sequence $T^{k}$ ot iterates of the transformation $T$, and the
related sequence of polynomials $p_{k}(x,z)=T^{k}z$.

Of course, $T^{0}$ is idA, and $T^{1}=T$.

Similarily, $p_{0}(x,z)=z$ and $p_{1}(x,z)$ is $p(x,z)$ itself.

For short, we say that $p(x,z)$ is $1$-adequate when $p(x,z)=z$.

Let $m>1$ be an integer. Henceforth, we assume that $K$ contains the $m$-th
roots of unity. We say that a polynomial $p(x,z$) is $m$-adequate when it has
the following form:

$p(x,z)=qz+r(x),$ $None$

where $q$ is a primitive $m$-th root of unity, and $r(x)$ is in $K[x]$.

So, $p(x,z)$ is “à la Granville” if and only if it is $2$-adequate.

When $p(x,z)$ is $m$-adequate for some $m\geq 1$, just say it is adequate.

### 2 Remark

When $p(x,z)$ is $m$-adequate, we have $p_{m}(x,z)=z$ and therefore $T^{m}=$
idA. More precisely, we have

$p_{0}(x,z)=z$

$p_{1}(x,z)=qz+r(x)$

$p_{2}(x,z)=q^{2}z+(q+1).r(x)$

$\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots$

$p_{k}(x,z)=q^{k}z+(q^{k-1}+\dots+q+1).r(x)$

$\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots\ldots$

$p_{m}(x,z)=z$.

When $r(x)\neq 0$, polynomials $p_{0},p_{1},\dots,p_{m-1}$ are two by two
mutually prime.

Let $d(P)$ be the degree in $z$ of each polynomial $P(x,z)$.]

### 3 Theorem

Let $p(x,z)$ be $m$-adequat. Let $b$ be the product of polynomials
$\\{p_{k}:0\leq k<m\\}$. A polynomial $F(x,z)$ is in $B$ if and only if it has
the following form:

$F(x,z)=G(b),\ \text{with}\ G(y)\in K[x][y].$ $None$

### Proof

Observe that polynomial $b$ is in $B$, so $TG(b)=G(b)$, and every polynomial
which has the form (6) is then in $B$.

For the converse, distinguish two cases.

If $r(x)=0$, then $F$ is in $B$ means $F(x,z)=F(x,qz)$; since $q^{m}=1$, we
conclude that $F$ is a polynomial in $z^{m}$.

Otherwise, proceed by induction on the degree, $d(F)$, of $F$ in $z$.

The case $d(F)=0$ is trivial.

### In four steps

### Lemma 1

If $F(x,z)$ is in $B$, each of the $p_{k}(x,z)$, $k\geq 0$, divides the
polynomial $W(x,z):=F(x,z)-F(x,0)$.

Indeed, $W(x,z)$ also is in $B$ and is divisible by $z$. Therefore, we must
have $W(x,z)=z.U(x,z)$, and for $k\geq 1$, we have

$W(x,z)=T^{k}U(x,z)=p_{k}(x,z).T^{k}U(x,z)$.

If $F(x,z)$ is in $B$, polynomial $F(x,z)-F(x,0)$ is divisible by $b$, i.e.,
we have $F(x,z)-F(x,0)=b(x,z).V(x,z)$. So, $V(x,z)$ is in $B$, moreover we
have $d(V)<d(F)$. Whence the conclusion.∎

### Complement

The following three conditions are equivalent.

(i) The set $C$ is not empty.

(ii) There exists an integer $m\geq 1$ such that $T^{m}=$ idA.

(iii) Polynomial $p(x,z)$ is adequate.

### Proof

If $p(x,z)$ is $m$-adequate, we already know that $C$ is not empty, and that
$T^{m}=$ idA.∎

### Lemma 2

If $d(p)\neq 1$, then $B=K[x]$.

Indeed, $d(TF)=d(F).d(p)$. If $F$ is in $B$, i.e., $TF=F$, we have
$d(F)=d(F).d(p)$, and if, moreover, $d(p)\neq 1$, then $d(F)=0$. ∎

### Lemma 3

If $C$ is not empty, the polynomial $p(x,z)$ is adequate.

### Proof

By the previous lemma, we know that $p(x,z)$ has the following form:

$p(x,z)=q(x).z+r(x),\ \text{with}\ q(x)\neq 0.$ $None$

Take $F\in B$ such that $d(F)=m\geq 1$.

Ordering $F$ in decreasing powers of $z$, write

$F=u.z^{m}+v.z^{(}m-1)+\cdots$

whence

$TF=u.(qz+r)^{m}+v.(qz+r)^{m-1}+\cdots$

$=u.q^{m}z^{m}+(u.m.q^{m-1}.r+v.q^{m-1}z^{m-1}+\cdots$

then, identifying the first two coefficients, obtain both equalities

$q^{m}=1$ $None$

$q^{m-1}.(m.r+v)=v.$ $None$

For $q=1$, we therefore have $r=0$. ∎

### Lemma 4

If $T^{m}=$ id A for some $m\geq 1$, the product $b$ of polynomials

in $\\{p_{k}:0\leq k<m\\}$ is invariant by $T$, and $d(b)\geq 1$.

This ends the proof ∎∎

## References

  * 1 E. Catalan, Mélanges Mathématiques, vol. 12, 1885, pp. 179-187.
  * 2 A. L. Cauchy, Rapport sur un Mémoire de M. Lamé, relatif au dernier théorème de Fermat, C. R. Acad. Sci. Paris 9 (1839), 359-363.
  * 3 A. L. Cauchy, Note sur Quelques Théorèmes d’Algèbre, Exerc. d’Analyse et de Phys. Math., vol. 2, Bachelier, 1841, pp. 137-144.
  * 4 D. Mirimanoff, Sur l’équation $(x+1)^{l}-x^{l}-1=0$, Nouv. Ann. Math. 3 (1903), 385-397.
  * 5 P. Ribenboim, 13 Lectures on Fermat’s Last Theorem, Springer, 1979.
  * 6 G. Terjanian, Sur la loi de réciprocité des puissances $l$-èmes, Acta Arith. 54 (1989), 87-125.


# Concurrent Emergence of Clustering, Flocking and Synchronization in Systems
of Interacting Agents

Trenton Gerew  _Department of Applied Mathematics, Illinois Institute of
Technology_ Ming Zhong  _Department of Applied Mathematics, Illinois Institute
of Technology_

###### Abstract

We present some recent development in modeling concurrent emergence of
collective behaviors, namely, the emergence of clustering, flocking and
synchronization at the same time. We derive two new models, namely
Swarmalator-Vicsek and Swarmalator-Cucker-Smale, which can produce the
synchronization and flocking with interesting spatial patterns. We present
extensive numerical insights into how the synchronization of phases can affect
both the spatial patterns and the flocking behavior.

##  1 Introduction

Collective behaviors (also known as self organization), such as clustering [9,
2, 13], flocking [7, 8, 1, 16], milling [4], swarming [6], synchronization[10,
20, 18, 17], describes how global orders, can emerge from initially randomized
configuration through local interactions between pairs of agents. It has
gathered considerable amount of research interests lately, see the Section $9$
in [13] and references therein. Moreover it also has a great number of of
applications in Physics (such as self assemble nanoparticles,
superconductivity, crystal growth), Chemistry (molecular self-assembly, self-
assembled monolayers), Biology (morphogenesis, social behavior of insects,
chemotaxis/phototaxis), human society (herd behavior, self-referentiality,
self-organizing market economy), etc. With proper modeling tools, single kind
of behavior can be replicated through carefully designed interacting agent
systems. We present in this paper that a couple new models that can reproduce
complicated behaviors with the mixture of three collective behaviors, namely
synchronization, flocking and clustering, which is partially motivated by the
synchronized fireflies behavior observed in the nature [17, 21, 15].

One of the first models which combines two collective behaviors together is
the swarmalator model proposed in [15]. This swarmalator model combines
clustering and synchronization together. To be precise, it is about the how
synchronization of phases affects the spatial pattern of the agents. The
system is described by an ODE system, i.e., for a system of $N$ interacting
agents, their states are evolved according to the following ODE equations,

$\begin{cases}\dot{\mathbf{x}}_{i}&=\sum_{j=1,j\neq
i}^{N}\phi^{E}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\xi_{j}-\xi_{i})(\mathbf{x}_{j}-\mathbf{x}_{i}),\\\
\dot{\xi}_{i}&=\sum_{j=1,j\neq
i}^{N}\phi^{\xi}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\xi_{j}-\xi_{i})(\xi_{j}-\xi_{i}),\end{cases}\quad
i=1,\cdots,N.$

Here the pair $(\mathbf{x}_{i},\xi_{i})\in{\mathbb{R}}^{d}\times{\mathbb{R}}$
presents the position/phase for the $i^{th}$ agent respectively in the
systems. Furthermore, the functions,
$\phi^{E},\phi^{\xi}:{\mathbb{R}}^{+}\times{\mathbb{R}}\rightarrow{\mathbb{R}}$,
are called the interaction functions, which govern how the $j^{th}$ agent
influences the change of position and phase of the $i^{th}$ agent. Moreover,
they depend only on the pairwise distance between the two agents as well as
the relative phase difference. Although whole system has local (pairwise)
interactions between the agents, it can produce interesting spatial patterns
with various phase distribution among the agents (static or non-static with
either synchronized distribution of phases), see [15] for detailed examples.

However, individuals in nature tend to align and move as others do in their
neighborhood [14]. The synchronized flapping of migratory birds in the V-shape
[1], [16], synchronized firefights (most significant sighting in the smokey
mountain area) are intriguing examples of concurrent emergence of more than
two collective behaviors. Hence the combination of only two behaviors is not
enough to model natural phenomena. Therefore we propose two new models, where
the the states of the agents now include $(\mathbf{x}_{i},\xi_{i},\theta_{i})$
or $(\mathbf{x}_{i},\xi_{i},\mathbf{v}_{i})$, where $\theta_{i}$ is the
heading (the polar angle of velocity of the $i^{th}$ agent) and
$\mathbf{v}_{i}$ is the velocity of the $i^{th}$ agent respectively. The ODE
system is changed to incorporate the additional state. For example, the
Swarmalator-Viscek model is governed by the following ODE system

$\begin{cases}\dot{\mathbf{x}}_{i}&=\mathbf{v}_{i}+\begin{bmatrix}\cos(\theta_{i})\\\
\sin(\theta_{i})\end{bmatrix}+\frac{1}{\left\lvert\Lambda_{i}\right\rvert}\sum_{j\in\Lambda_{i}}\phi^{E}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\theta_{j}-\theta_{i},\xi_{j}-\xi_{i})(\mathbf{x}_{j}-\mathbf{x}_{i}),\\\
\dot{\theta}_{i}&=\frac{1}{\left\lvert\Lambda_{i}\right\rvert}\sum_{j\in\Lambda_{i}}\phi^{\theta}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\theta_{j}-\theta_{i},\xi_{j}-\xi_{i})(\theta_{j}-\theta_{i}),\\\
\dot{\xi}_{i}&=\omega_{i}+\frac{1}{\left\lvert\Lambda_{i}\right\rvert}\sum_{j\in\Lambda_{i}}\phi^{\xi}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\theta_{j}-\theta_{i},\xi_{j}-\xi_{i})(\xi_{j}-\xi_{i}),\end{cases}\quad
i=1,\cdots,N.$

This new system will have the states, namely position/heading angle/phase, for
the $i^{th}$ agent, denoted as $(\mathbf{x}_{i},\theta_{i},\xi_{i})$.
Furthermore, the interaction of these internal states together with states
from their neighboring agents can produce complex behaviors, especially the
concurrent emergence of clustering, flocking and synchronization. With the
inclusion of three different states for each individual agent, these two new
models can now combine the flocking behavior described in [19, 3] to
synchronization of phases and spatial patterns described in [15]. We are able
to show through extensive numerical experiments that when phases also enter
into the interaction with position and velocity, intriguing spatial patterns
can emerge. For detailed discussion as well as the model equations for the
Swarmalator-Cucker-Smale model, see section 2.

The remaining of this paper is arranged as follows: Section 2 describes in
details how the two new models are derived; Section 3 shows detailed
experiments of the different patterns which our new models can produce; and we
conclude our paper in Section 4.

##  2 Model Description

We describe in this section in details about the two new models, namely
Swarmalator-Viscek and Swarmalator-Cucker-Smale. The two models govern how the
state tuples, i.e. $(\mathbf{x}_{i},\xi_{i},\theta_{i})$ or
$(\mathbf{x}_{i},\xi_{i},\mathbf{v}_{i})$, would interact with each other and
with those from neighboring agents, and such intertwined relationship for
these state variables would lead to interesting flocking and synchronization
with complicated spatial patterns.

###  2.1 Swarmalator-Viscek Model

First, we consider a second-order swarmalator model with headings, i.e.
$\theta_{i}\in\mathcal{S}^{2}$ (the $2D$ unit sphere), for each agent. The
governing equations for a population of $N$ swarmalators are

$\begin{cases}\dot{\mathbf{x}}_{i}&=\mathbf{v}_{0}+\begin{bmatrix}\cos(\theta_{i})\\\
\sin(\theta_{i})\end{bmatrix}+\frac{1}{\left\lvert\Lambda_{i}\right\rvert}\sum_{j\in\Lambda_{i}}\phi^{E}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\theta_{j}-\theta_{i},\xi_{j}-\xi_{i})(\mathbf{x}_{j}-\mathbf{x}_{i}),\\\
\dot{\theta}_{i}&=\frac{1}{\left\lvert\Lambda_{i}\right\rvert}\sum_{j\in\Lambda_{i}}\phi^{\theta}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\theta_{j}-\theta_{i},\xi_{j}-\xi_{i})(\theta_{j}-\theta_{i}),\\\
\dot{\xi}_{i}&=\omega_{0}+\frac{1}{\left\lvert\Lambda_{i}\right\rvert}\sum_{j\in\Lambda_{i}}\phi^{\xi}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\theta_{j}-\theta_{i},\xi_{j}-\xi_{i})(\xi_{j}-\xi_{i}),\end{cases}\quad
i=1,\cdots,N.$ (1)

Here the tuple $(\mathbf{x}_{i},\xi_{i},\theta_{i})$ describes the
position/phase/heading (in terms of polar angle of the velocity vector) of the
$i^{th}$ agent respectively. Moreover, $\mathbf{v}_{0},\omega_{0}$ are
constants describing its self-propulsion velocity and natural frequency,
respectively. The functions $(\phi^{E},\phi^{\xi},\phi^{\theta})$ define the
interaction between agent $j$ and agent $i$. $\Lambda_{i}$ is the neighborhood
around agent $i$, and for the purposes of this paper, we define the
neighborhood of interaction as in [19]:

######  Definition 2.1 (Neighborhood of Interaction).

Given an interaction radius $r>0$, the neighborhood of interaction for
swarmalator $i$ is

$\Lambda_{i}=\left\\{j\in[N]\mid\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert\leq
r\right\\},\quad[N]=\\{1,\cdots,N\\}.$ (2)

######  Remark 2.1.

Note that the choice of the neighborhood parameter $r>0$ will affect the
spatial patterns significantly, see section 3 for details.

When compared to the original second-order swarmalator model [14], our mode
has the local dependence on neighborhood to both the position and phase
equations. It is a more natural modeling of the swarming oscillators, since
the decision of neighborhood for each agent will have effect on all three
states instead of only on the local heading. In fact, a more precision
description of the swarming oscillators should include three different kinds
of neighborhood for each state equation. However, we will then introduce too
many changing parameters in this model. We will leave such discussion in the
future work. Moreover, the dependence of the interaction functions on pairwise
state variables can produce more complex patterns. Difference choices of
$(\phi^{E},\phi^{\xi},\phi^{\theta})$ can also lead to interesting patterns,
see [15]. Moreover, these interaction functions can be also data-driven i.e.
learned from observation data [11, 12, 5].

###  2.2 Swarmalator-Cucker-Smale Model

An important step in modeling self organization for local models is in
choosing a proper neighborhood. For the Swarmalator-Viscek model in
particular, it is crucial on how to choose the radius of interaction, i.e.
$r>0$ in (1). As demonstrated in section 3, the choice of $r$ has a
significant influence on the behavior of the whole system. We address the
issue by switching the emphasis of local neighborhood from selection of a
subset of agents to proper choices of interaction functions, i.e. the
interaction function on the alignment of velocities of agents (similar to that
in [3]). Therefore, we present the following Swarmalator-Cucker-Smale model,
for $i=1,\cdots,N$,

$\begin{cases}\dot{\mathbf{x}}_{i}&=\mathbf{v}_{i},\\\
\dot{\mathbf{v}}_{i}&=\frac{1}{N}\sum_{j=1,j\neq
i}^{N}\Big{[}\phi^{E}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\left\lvert\mathbf{v}_{j}-\mathbf{v}_{i}\right\rvert,\xi_{j}-\xi_{i})(\mathbf{x}_{j}-\mathbf{x}_{i})\\\
&\quad+\phi^{\mathbf{v}}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\left\lvert\mathbf{v}_{j}-\mathbf{v}_{i}\right\rvert,\xi_{j}-\xi_{i})(\mathbf{v}_{j}-\mathbf{v}_{i})\Big{]}\\\
\dot{\xi}_{i}&=\omega_{i}+\frac{1}{N}\sum_{j=1,j\neq
i}^{N}\phi^{\xi}(\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert,\left\lvert\mathbf{v}_{j}-\mathbf{v}_{i}\right\rvert,\xi_{j}-\xi_{i})(\xi_{j}-\xi_{i}),\\\
\end{cases}$ (3)

Here $\mathbf{v}_{i}$ is the velocity for the $i^{th}$ agent, which can change
in both direction and magnitude (as compared to the Swarmalator-Viscek model,
where the velocity has fixed magnitude). The functions
$(\phi^{E},\phi^{\xi},\phi^{\theta})$ can be compactly supported to reflect
the local neighborhood.

(Add more)

##  3 Numerical Results

We conduct numerical tests of the two models we propose in the previous
sections.

###  3.1 Swarmalator-Vicsek

In order to compare our results fairly to the results in [15], we choose the
pair, $(\phi^{E},\phi^{\xi})$, to be the following

$\phi^{E}(r,s)=\frac{A+J\cos(\xi_{j}-\xi_{i})}{r}-\frac{B}{r^{2}}\quad\text{and}\quad\phi^{\xi}(r,s)=\frac{K\sin(s)}{r}.$

Here $r_{ij}=\left\lvert\mathbf{x}_{j}-\mathbf{x}_{i}\right\rvert$ and
$s_{ij}=\xi_{j}-\xi_{i}$ represents the pairwise distance and pairwise
difference of phases variables. The interaction function on the heading
variable $\theta$, i.e. $\phi^{\theta}$, is taken as

$\phi^{\theta}(z)=1,\quad\text{for all $z_{ij}=\theta_{j}-\theta_{i}$}.$

This is a continuation of the setup in [19]. One can use other interaction
functions. We then set the natural frequency $\omega_{i}=0$ for each agent.
Following the setup in [15], we set $A=B=1$. Therefore, our model has a set of
three parameters, i.e. $(r,J,K)$, to choose for different types of system
behaviors. The phase coupling strength is represented by the parameter $K$.
For $K>0$, the phase coupling tends to minimize the phase difference between
swarmalators, and for $K<0$ the phase difference is maximized. When $K=0$, the
phases of the swarmalators do not change. The strength at which phase
similarity enhances spatial attraction is measured by parameter $J$. If $J>0$,
swarmalators in the same phase are attracted to each other. For $J<0$,
swarmalators move towards others with the opposite phase. When $J=0$, the
spatial attraction is independent of phase. We constrain $J\in[-1,1]$ to keep
$\frac{A+J\cos(\xi_{j}-\xi_{i})}{r}-\frac{B}{r^{2}}\geq 0$. We performed
numerical experiments to identify the behavior of our model. All the
simulations were run using MATLAB’s standard ODE solver `ode45`. Initially,
the swarmalators were positioned in a box of length $1$ with phases in
$[-\pi,\pi]$ and orientations in $[0,2\pi]$, all drawn uniformly at random. We
choose $\mathbf{v}_{i}=0.003$ (the self propelling velocity) to be within the
range of velocities explored in [19]. We found the system’s behavior is highly
dependent on the chosen radius of interaction, and new types of grouping
behavior emerge. To compare our model to the swarmalator model, we use the
parameters that describe the five unique states found in [15].

  1. 1.

Static asynchronous state. In the original swarmalator model the parameters
$K>0$ and all $J$ form a static synchronous phase in which the swarmalators
form a circularly symmetric distribution and are fully synchronized in phase
and orientation. In our model, this combination of $J$ and $K$ form what is
essentially the static asynchronous state from [15], as shown in Figure 1. The
swarmalators separate faster than they can synchronize, and end up forming a
crystal-like distribution in which all phases $\theta$ can occur. At small
radii $r$, no extremely little interaction occurs at all and the orientations
$\beta$ are random as well. However, after about $r=0.65$ the interactions are
strong enough that the orientations converge rapidly. Note that in Figure 6(a)
the swarmalators are distributed uniformly, so that every phase occurs
everywhere.

(a) $r=0.20$

(b) $r=0.65$

(c) $r=1.40$

Figure 1: Static asynchronous state: $J=0.1,\ K=1$ with $N=300$ swarmalators
after $T=50$ time units.

  2. 2.

Static gradient state. Swarmalators can also form a gradient state, as seen in
Figure 2. This state occurs for $J<0,\ K<0$ and also for $J>0$ if
$J<\left\lvert K_{c}\right\rvert$ as shown in [15]. In the original
swarmalator model, this parameter set forms the completely asynchronous state.

(a) $r=0.20$

(b) $r=0.65$

(c) $r=1.40$

Figure 2: Static gradient state: $J=0.1,\ K=-1$ with $N=300$ swarmalators
after $T=50$ time units.

  3. 3.

Static wave state. The special case with $K=0$ and $J>0$ forms a curvy wave-
like structure. The swarmalators are frozen in their initial phases, and $J>0$
causes the swarmalators to be attracted towards others with similar phases. In
the infinite-radius case, this interaction forms an annular structure called
the static phase wave. With the limited radius, however, swarmalators can’t
“see” all others of the same phase and thus settle into a repeating phase-
wave-like pattern. The new structure is not annular, but the spatial angle
$\phi$ and the phase $\theta$ of each swarmalator are still moderately
correlated, as seen in Figure 3 and 6(c).

(a) $r=0.20$

(b) $r=0.65$

(c) $r=1.40$

Figure 3: Static wave state: $J=1,\ K=0$ with $N=300$ swarmalators after
$T=50$ time units.

  4. 4.

Clustered state. Moving into $K<0$, we encounter the clustered state. As seen
in Figures 4 and 7(a), the population splinters into groups of distinct
phases. In the movies included in the supplemental materials, it can be seen
that some nearby groups will collide and merge into one phase. It is unclear
what determines the number of clusters, but in general fewer are found when a
larger radius is used. This state corresponds with the splintered phase wave
state in [15].

(a) $r=0.20$

(b) $r=0.65$

(c) $r=1.40$

Figure 4: Clustered state: $J=1,\ K=-0.1$ with $N=300$ swarmalators after
$T=50$ time units.

  5. 5.

Static synchronous state. The final state is the synchronous state which
occurs after $K$ is further increased. The parameter set for this state
corresponds to the active phase wave for the unaltered swarmalator model. In
contrast, here the swarmalators rapidly adopt the same phase and condense into
a tight group as seen in Figures 5 and 7(b).

(a) $r=0.20$

(b) $r=0.65$

(c) $r=1.40$

Figure 5: Static synchronous state: $J=1,\ K=-0.75$ with $N=300$ swarmalators
after $T=50$ time units.

Conclusion: in all cases, pattern formation does not begin to occur until
about $r=0.6$. The reason for this is not clear, and the critical radius
appears to occur for all starting densities of swarmalators. Figure 6 shows
the static spatial patterns with different phases distribution for the
relationship between the distribution of the polar angles of position vs the
phase variable.

(a) Static asynchronous state
$(r,J,K)=(1.25,0.1,1)$

(b) Static gradient state
$(r,J,K)=(1.10,0.1,-1)$

(c) Static wave state
$(r,J,K)=(1.25,1,0)$

Figure 6: Scatter plots in $(\psi,\xi)$ space for $N=300$ swarmalators and
$T=50$ time steps, where $\psi=\tan^{-1}(y/x)$. The color of each point
represents the orientation $\theta$.

Figure 7 shows the astatic spatial patterns with different phases distribution
for the relationship between the distribution of the polar angles of position
vs the phase variable.

(a) Clustered state
$(r,J,K)=(1.10,1,-0.1)$

(b) Static synchronous state
$(r,J,K)=(1.25,1,-0.75)$

Figure 7: Scatter plots in $(\psi,\xi)$ space for $N=300$ swarmalators and
$T=50$ time steps, where $\psi=\tan^{-1}(y/x)$. The color of each point
represents the orientation $\theta$.

###  3.2 Swarmalator-Cucker-Smale

The model in (3) behaves essentially the same as a global swarmalator-Vicsek
model posed in [15], where the parameter controlling the size of the
neighborhood, i.e. $r$, is going to infinity. We performed simulations of the
model with $N=500$ swarmalators. The swarmalators were initially positioned in
a box of length 2 with phases in $[-\pi,\pi]$ and velocity vectors in
$[-0.1,0.1]\times[-0.1,0.1]$, all drawn uniformly at random. The same three
stationary states and two active states from [15] are clearly identifiable.
First, the three static states are present.

(a) Static sync state
$(J,K)=(0.1,1)$

(b) Static async state
$(J,K)=(0.1,-1)$

(c) Static phase wave state
$(J,K)=(1,0)$

Figure 8: Scatter plots of the three stationary states in the $(x,y)$ plane
with $N=500$ swarmalators for $T=300$ time units.

Next, the moving states are also present.

(a) Splintered phase wave state
$(J,K)=(1,-0.1)$

(b) Active phase wave state
$(J,K)=(1,-0.75)$

Figure 9: Scatter plots of the two active states in the $(x,y)$ plane with
$N=500$ swarmalators for $T=300$ time units.

##  4 Conclusion

We have examined the dynamics of two types of systems which give concurrent
emergence of clustering, flocking and synchronization. The two systems
describe orientable swarmalators with a local or global neighborhood. These
agents have spatial, phase, and orientation degrees of freedom which grants
them the ability to flock (common velocity) and cluster (spatial pattern), as
well as synchronize (common phase). The position/velocity (or heading)/phase
are coupled through a set of interaction functions. These models would be the
first step in studying the behavior of real-world swarming oscillators (for
example synchronized firefights). We also look to biological systems for
motivation to choose the radius of interaction, based on sensory limitations.
There are also a number of theoretical questions left to open, to be answered
in future work.

## References

  * [1] I. L. Bajec and F. H. Heppner.  Organized flight in birds.  Animal Behaviour, 78(4):777–789, 2009.
  * [2] V. D. Blondel, J. M. Hendricks, and J. N. Tsitsiklis.  On krause’s multi-agent consensus model with state-dependent connectivity.  IEEE Trans Automat. Control, 54:2586 – 2597, 2009.
  * [3] F. Cucker and S. Smale.  Emergent behavior in flocks.  Automatic Control, IEEE Transactions on, 52:852 – 862, 06 2007\.
  * [4] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes.  Self-propelled particles with soft-core interactions: Patterns, stability, and collapse.  Phys. Rev. Lett., 96:104302, Mar 2006.
  * [5] J. Feng, M. Maggioni, P. Martin, and M. Zhong.  Learning interaction variables and kernels from observations of agent-based systems.  IFAC-PapersOnLine, 55(30):162–167, 2022.  25th International Symposium on Mathematical Theory of Networks and Systems MTNS 2022.
  * [6] A. Filella, F. m. c. Nadal, C. Sire, E. Kanso, and C. Eloy.  Model of collective fish behavior with hydrodynamic interactions.  Phys. Rev. Lett., 120:198101, May 2018.
  * [7] H. Hildenbrandt, C. Carere, and C. Hemelrijk.  Self-organized aerial displays of thousands of starlings: a model.  Behavioral Ecology, 21(6):1349–1359, 10 2010.
  * [8] A. J. King and D. J. Sumpter.  Murmurations.  Current Biology, 22(4):R112–R114, 2 2012.
  * [9] U. Krause.  A discrete nonlinear and non-autonomous model of consensus formation.  In S. Elyadi, G. Ladas, J. Popenda, and J. Rakowski, editors, Communications in Difference Equations, pages 227–236. Gordon and Breach Pub., Amsterdam, 2000.
  * [10] Y. Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In H. Araki, editor, International Symposium on Mathematical Problems in Theoretical Physics, pages 420–422, Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.
  * [11] F. Lu, M. Zhong, S. Tang, and M. Maggioni.  Nonparametric inference of interaction laws in systems of agents from trajectory data.  Proceedings of the National Academy of Sciences, 116(29):14424–14433, 2019.
  * [12] M. Maggioni, J. J. Miller, H. Qiu, and M. Zhong.  Learning interaction kernels for agent systems on riemannian manifolds.  In M. Meila and T. Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 7290–7300. PMLR, 18–24 Jul 2021.
  * [13] S. Motsch and E. Tadmor.  Heterophilious dynamics enhances consensus.  SIAM Review, 56(4):577–621, 2014.
  * [14] D. O’Brien.  Analysis of the internal arrangement of individuals within crustacean aggregations (euphausiacea, mysidacea).  Journal of Experimental Marine Biology and Ecology, 128(1):1–30, 1989.
  * [15] K. O’Keeffe, H. Hong, and S. Strogatz.  Oscillators that sync and swarm.  Nature Communications, 8, 11 2017.
  * [16] S. J. Portugal, T. T. Hubel, J. Fritz, S. Heese, D. Trobe, B. Voelkl, S. Hailes, A. M. Wilson, and J. R. Usherwood.  Upwash exploitation and downwash avoidance by flap phasing in ibis formation flight.  Nature, 505:399 – 402, January 2014.
  * [17] S. Strogatz.  Sync: How Order Emerges From Chaos In the Universe, Nature, and Daily Life.  Hyperion, 2012.
  * [18] S. H. Strogatz.  From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 143(1):1–20, 2000.
  * [19] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet.  Novel type of phase transition in a system of self-driven particles.  Phys. Rev. Lett., 75:1226–1229, Aug 1995.
  * [20] A. T. Winfree.  Biological rhythms and the behavior of populations of coupled oscillators.  Journal of Theoretical Biology, 16(1):15–42, 1967.
  * [21] A. T. Winfree.  The Geometry of Biological Time.  Springer, 1980.