# Semidefinite programming bounds for distance distribution of spherical codes

Oleg R. Musin

###### Abstract

We present an extension of known semidefinite and linear programming upper
bounds for spherical codes. We apply the main result for the distance
distribution of a spherical code and show that this method can work
effectively. In particular, we get a shorter solution to the kissing number
problem in dimension 4.

Mathematics Subject Classification (2010) 90C22, 52C17

##  1 Introduction

Let $G_{k}^{(n)}(t)$ (with $G_{k}^{(n)}(1)=1$ and $\deg(G_{k}^{(n)})=k$ ) be
Gegenbauer polynomials that are orthogonal on the interval $[-1,1]$ with
respect to the weight function $(1-t^{2})^{(n-3)/2}$.

Let $C$ be an $N$–element subset of the unit sphere
${\mathbb{S}}^{n-1}\subset{\mathbb{R}}^{n}$. We define the $k$-th moment of
$C$ as

$M_{k}(C):=\sum\limits_{(c,c^{\prime})\in C^{2}}{G_{k}^{(n)}(c\cdot
c^{\prime})}$

The positive semidefinite property of Gegenabauer polynomials yields that

$M_{k}(C)\geq 0\;\mbox{ for all }\;k=1,2,...$ $None$

Let $f$ be a non-negative linear combination of Gegenabauer polynomials:

$f(t)=\sum\limits_{k=0}^{d}{f_{k}G_{k}^{(n)}(t)},\,\,f_{k}\geq 0\,\mbox{ for
all }\;k=1,,..,d.$

Then

$S_{f}(C):=\sum\limits_{(c,c^{\prime})\in C^{2}}{f(c\cdot
c^{\prime})}=f_{0}N^{2}+\sum\limits_{k=1}^{d}{f_{k}M_{k}(C)}\geq f_{0}N^{2}.$
$None$

The distance distribution of $C$ with respect to $u\in C$ is the system of
numbers $\\{A_{t}(u):-1\leq t\leq 1\\}$, where

$A_{t}(u):=|\\{v\in C:v\cdot u=t\\}|,$

and the distance distribution of $C$ is the system of numbers $\\{A_{t}:-1\leq
t\leq 1\\}$, where

$A_{t}(C)=A_{t}:=\frac{1}{N}\sum\limits_{u\in C}{A_{t}(u)}.$

Denote

$R_{f}(C):=\sum\limits_{t\in[-1,1)}{f(t)A_{t}(C)}.$

Then

$S_{f}(C)=NR_{f}(C)+Nf(1)$

and (1.2) implies

$R_{f}(C)\geq f_{0}N-f(1).$ $None$

The semidefinite programming (SDP) method for spherical codes was proposed by
Bachoc and Vallentin [1] with further applications and extensions in [2, 3, 9,
10, 12, 15, 16, 18, 20].

The positive–semidefinite property of Gegenbauer polynomials yields the
positive–semidefinite property of matrices $S_{k}^{n}$. Now consider
polynomials $F$ that were defined by Bachoc and Vallentin. Let $F(t,u,v)$ be a
symmetric polynomial with expansion

$F(t,u,v)=\sum\limits_{k=0}^{d}{\langle H_{k},S^{n}_{k}(t,u,v)\rangle}$

in terms of the matrices $S^{n}_{k}$.

Suppose that all matrices $H_{k}$ with $k>0$ are positive semidefinite and for
a given $F_{0}\in\mathbb{R}$, $H_{0}-F_{0}E_{0}$ is also positive
semidefinite. (Here $E_{0}$ denote a matrix whose only nonzero entry is the
top left corner which contains 1). Then

$S_{F}(C):=\sum\limits_{(x,y,z)\in C^{3}}{F(x\cdot y,x\cdot z,y\cdot z)}\geq
F_{0}N^{3}.$ $None$

In [20] we consider an extension of known semidefinite and linear programming
upper bounds for spherical codes and a version of this bound for distance
graphs. This paper is a continuation and extension of [20]. We show how the
bounds from [20] can be applied to the distance distribution of spherical
codes.

In Section 2 for $(N,n,T)$ spherical codes $C$ we present a general 3-point
bound, see Theorem 2.2. Actually, this theorem can be considered as a lower
bound for $E_{g}(C)$ (or equivalently for $R_{g}(C)$, see Theorem 3.1), where
$g:T\to\mathbb{R}$. We obtain functions $g$ on $T\subset[-1,1)$ that can play
the same role as a nonnegative linear combination of Gegenbauer polynomials on
$[-1,1]$.

In Section 4 we show that this method can work effectively. For $T=[-1,0.5]$
we consider two polynomials $g_{1}$ and $g_{2}$ that found using the SDP.

The expansion of the first polynomial in Gegenbauer polynomials has negative
coefficients. However, the resulting boundary for it is almost exact and with
its help we get a shorter solution to the kissing number problem in dimension
4.

The second polynomial has only positive coefficients, but the SDP bound for it
is much stronger than (1.3).

Section 5 contains several possible applications of Theorems 2.2 and 3.1 and
their generalizations.

##  2 General bounds for spherical codes

Let $C$ be an $N$–element subset of the unit sphere
${\mathbb{S}}^{n-1}\subset{\mathbb{R}}^{n}$. Denote

$I(C):=\\{t=x\cdot y\,|\,x,y\in C\,\&\,x\neq y\\}.$

Let $T\subset[-1,1)$. We say that $C$ is an $(N,n,T)$ spherical code if
$I(C)\subset T$.

Let $g$ be a real function on $I(C)$. Define

$E_{g}(C):=\sum\limits_{(x,y)\in C^{2},x\neq y}{g(x\cdot y)}$

###  2.1 General 2-point bound

######  Theorem 2.1.

Let $C$ be an $(N,n,T)$ spherical code. Suppose $g:T\to\mathbb{R}$,
$f:[-1,1]\to\mathbb{R}$ and $f_{0}\in\mathbb{R}$ are such that

  1. 1.

$f(t)\leq g(t)$ for all $t\in T$.

  2. 2.

$S_{f}(C)\geq f_{0}N^{2}$.

Then

$Nf(1)+E_{g}(C)\geq f_{0}N^{2}.$

###### Proof.

Note that for all $x\in C\subset\mathbb{S}^{n-1}$ we have $x^{2}=x\cdot x=1$.
Then

$f_{0}N^{2}\leq S_{f}(C)=Nf(1)+E_{f}(C)\leq Nf(1)+E_{g}(C).$

∎

Note that if $f$ is a non-negative linear combination of the Gegenbauer
polynomials then (1.2) implies that $f$ satisfies assumption 2 in Theorem 2.1.

Suppose $f_{0}>0$ and $g(t)=0$ for all $t\in T$. Since $E_{g}(C)=0$, the
theorem yields that

$N\leq\frac{f(1)}{f_{0}}.$

This bound is called the linear programming (LP) or Delsarte’s bound for
spherical codes.

Let $q:(0,4]\to\mathbb{R}$ be any function. Then for positive semidefinite $f$
and $g(t)=q(2-2t)$ Theorem 2.1 implies that
Every set of $N$ points on ${\mathbb{S}}^{n-1}$ has potential energy
$E_{q}(C):=E_{g}(C)$ at least

$f_{0}N^{2}-Nf(1).$

This fact first proved by Yudin [24] and has a lot of applications.

###  2.2 General 3-point bound

The Gram matrix of a set of vectors $v_{1},…,v_{n}$ in $\mathbb{R}^{d}$ is the
matrix of inner products, whose entries are given by the inner product
$G_{ij}=v_{i}\cdot v_{j}$. The Gram matrix is symmetric and positive
semidefinite. Moreover, a symmetric matrix $M$ is positive semidefinite if and
only if it is the Gram matrix of some vectors $v_{1},…,v_{n}$.

Let

$M_{3}=\left(\begin{array}[]{ccccc}1&u&v\\\ u&1&t\\\ v&t&1\end{array}\right).$

If this matrix is positive semidefinite ($M_{3}\succeq 0$), then there are
three distinct points $x,y,z$ in $\mathbb{S}^{2}$ such that $t=x\cdot
y,u=x\cdot z$, and $v=y\cdot z$. It is easy to see that $M_{3}\succeq 0$ if
and only if $t,u,v\in[-1,1]$ and $\det(M_{3})=1+2tuv-t^{2}-u^{2}-v^{2}\geq 0$.

This fact explains the following definition.

$D_{3}(T):=\left\\{(t,u,v):t,u,v\in T\,\&\,1+2tuv-t^{2}-u^{2}-v^{2}\geq
0\right\\},\;\;T\subset[-1,1).$

######  Theorem 2.2.

Let $C$ be an $(N,n,T)$ spherical code and $F:[-1,1]^{3}\to\mathbb{R}$ be a
symmetric function. Suppose $f:T\to\mathbb{R}$ and $g:T\to\mathbb{R}$, are
such that

  1. 1.

$F(1,t,t)\leq f(t)$ for all $t\in T$,

  2. 2.

$F(t,u,v)\leq g(t)+g(u)+g(v)$ for all $(t,u,v)\in D_{3}(T)$.

If $S_{F}(C)\geq F_{0}N^{3},$ where $F_{0}\in\mathbb{R}$, then

$NF(1,1,1)+3E_{f}(C)+(3N-6)E_{g}(C)\geq F_{0}N^{3}.$

###### Proof.

We have $S_{F}(C)=S_{1}+S_{2}+S_{3}$, where

$S_{1}=\sum\limits_{x=y=z}{F(x\cdot y,x\cdot z,y\cdot z)},$
$S_{2}=\sum\limits_{(x,y,z)\in H}{F(x\cdot y,x\cdot z,y\cdot
z)},\;H=\\{(x,y,z)\in C^{3}\,|\,x=y\neq z\,or\,x=z\neq y\,or\,x\neq y=z\\},$
$S_{3}=\sum\limits_{x\neq y\neq z\neq x}{F(x\cdot y,x\cdot z,y\cdot z)}.$

Since $x^{2}=1$ for all $x\in C$, we have

$S_{1}=\sum\limits_{x\in C}{F(x^{2},x^{2},x^{2})}=NF(1,1,1),$
$S_{2}=3\sum\limits_{(x,y)\in C^{2},x\neq y}{F(1,x\cdot y,x\cdot y)}\leq
3\sum\limits_{x\neq y}{f(x\cdot y)}=3E_{f}(C)).$

By assumption 2

$S_{3}\leq\sum\limits_{x\neq y\neq z\neq x}{(g(x\cdot y)+g(x\cdot z)+g(y\cdot
z))}=3(N-2)E_{g}(C).$

Thus

$S_{F}(C)\leq NF(1,1,1)+3E_{f}(C)+3(N-2)E_{g}(C).$

∎

######  Corollary 2.1.

Under the assumptions of Theorem 2.2 let $f(t)=p(t)-q(t)$ with
$p:T\to\mathbb{R}$ and $q:[-1,1]\to\mathbb{R}$. If $S_{q}(C)\geq 0$, then

$NF(1,1,1)+3Nq(1)+3E_{p}(C)+(3N-6)E_{g}(C)\geq F_{0}N^{3}.$

###### Proof.

$S_{2}\leq 3E_{f}(C)=3E_{p}(C)-3E_{q}(C)).$

Since $q(C)=Nq(1)+E_{q}(C)\geq 0$, we have

$S_{2}\leq 3Nq(1)+3E_{p}(C).$

Thus

$S_{F}(C)\leq NF(1,1,1)+3Nq(1)+3E_{p}(C)+3(N-2)E_{g}(C).$

∎

######  Corollary 2.2.

Under the assumptions of Theorem 2.2 let $f(t)=B+2g(t)-q(t)$ with
$q:[-1,1]\to\mathbb{R}$. If $S_{q}(C)\geq 0$, then

$F(1,1,1)+3q(1)+3(N-1)B+3E_{g}(C)\geq F_{0}N^{2}.$

###### Proof.

$S_{2}\leq 3Nq(1)+3E_{p}(C)=3Nq(1)+3\sum\limits_{x\neq y}{(B+2g(x\cdot
y))}=3Nq(1)+3N(N-1)B+6E_{g}(C).$

Then

$F_{0}N^{3}\leq S_{F}(C)\leq
NF(1,1,1)+3Nq(1)+3N(N-1)B+6E_{g}(C)+3(N-2)E_{g}(C)$
$=NF(1,1,1)+3Nq(1)+3N(N-1)B+3NE_{g}(C).$

Thus

$F(1,1,1)+3q(1)+3(N-1)B+3E_{g}(C)\geq F_{0}N^{2}.$

∎

Note that if $(F-F_{0})$ and $q$ are positive semidefinite then $S_{F}(C)\geq
F_{0}N^{3}$ and $q(C)\geq 0$. This makes it possible to find new bounds for
$N$ and $E_{g}(C)$ using the SDP. We will look at these methods in more detail
in later sections.

Suppose $g(t)\equiv 0$ and $F_{0}>0$. Then Corollary 2.2 yields

$N^{2}\leq\frac{F(1,1,1)+3q(1)+3(N-1)B}{F_{0}}.$

This inequality as well as Corollary 2.1 with $g=0$ first were proposed by
Bachoc and Vallentin [1, 2] with further applications and extensions in [3,
10, 15, 16, 18]. In particular, Cohn and Woo [10] got three-point bounds for
potential energy minimization.

###  2.3 General $k$-point bound

Theorems 2.1 and 2.2 can be extensed for all $k$: $2\leq k\leq n-2$. Theorem
5.4 from our paper [18] is a particular case of this general theorem.

It is clear how to derive a generalization of Theorem 2.2, see some details in
Section 5 [18] and [12]. However, the resulting formulas for $k>3$ are quite
cumbersome. We decided not to present the general theorem here even for the
case $k=4$.

##  3 SDP bounds for the distance distribution

Let $C$ be an $(N,n,T)$ spherical code. In the Introduction we defined
$A_{t}(C)$ and $R_{f}(C)$. It is clear that $A_{1}=1$, $A_{t}=0$ for all
$t\neq 1$ and $t\notin T$,

$\sum\limits_{t\in T}{A_{t}}=N-1\quad\mbox{and}\quad E_{f}(C)=NR_{f}(C).$

The following theorem is a restatement of Theorem 2.2 for the distance
distribution of spherical codes.

######  Theorem 3.1.

Let $C$ be an $(N,n,T)$ spherical code and $F:[-1,1]^{3}\to\mathbb{R}$ be a
symmetric function. Suppose $h:T\to\mathbb{R}$ and $g:T\to\mathbb{R}$, are
such that

  1. 1.

$h(t)+h_{0}+F(1,t,t)\leq 2g(t)$ for all $t\in T$,

  2. 2.

$F(t,u,v)\leq g(t)+g(u)+g(v)$ for all $(t,u,v)\in D_{3}(T)$.

If $S_{F}(C)\geq F_{0}N^{3}$, then

$R_{g}(C)=\sum\limits_{t\in
T}{A_{t}\,g(t)}\geq\frac{1}{3}F_{0}N+\frac{1}{3}h_{0}-\frac{1}{3N}F(1,1,1)+\frac{1}{N^{2}}E_{h}$

######  Corollary 3.1.

Under the assumptions of Theorem 3.1 let $F_{0}=0$ and $h_{0}=1$. If
$S_{h}(C)\geq 0$ then

$R_{g}(C)=\sum\limits_{t\in T}{A_{t}\,g(t)}\geq B(N):=\frac{N-M}{3N},\quad
M=F(1,1,1)+3h(1)$

###### Proof.

Since $S_{h}(C)=Nh(1)+E_{h}(C)\geq 0$, we have $E_{h}(C)\geq-Nh(1)$. Thus,
Theorem 3.1 yields the inequality. ∎

##  4 Some applications. A shorter proof of the kissing number problem in four
dimensions.

In this section we show that the method discussed in Sections 2 and 3 works
effectively. We consider polynomials $g_{1}$ and $g_{2}$ found using SDP for
the corresponding inequalities on the distribution of distances of spherical
codes cannot be found using LP bounds. These polynomial satisfy assumptions of
Corollary 3.1 with $n=4$ and $T=[-1,0.5]$. Using $g_{1}$ and the same method
from our paper [17], discarding the most difficult case of five points we get
a shorter proof that the kissing number in four dimensions $k(4)=24$.

Remark. I would like to note that these polynomials were found by Maria
Dostert during our work on the uniqueness problem of the maximum kissing
arrangement in dimension 4 [13]. Maria also found several tight polynomials
for other intervals. To do this, she had to overcome a large number of
technical obstacles and, along the way, develop and improve an algorithm for
finding SDP bounds for the distance distribution.

###  4.1 Two examples.

Let $T=[-1,0.5]$. In this case we denote a spherical code $(N,n,T)$ by
$(N,n,\pi/3)$. Here we consider the bound given by Corollary 3.1 for
$(N,4,\pi/3)$ spherical codes.

1. Let the expansion of $g_{1}$ in Gegenbauer polynomials $G_{k}^{(4)}$ have the following coefficients:
$[c_{0},...,c_{22}]$= [-0.5438, -2.0024, -3.8887, -5.6414, -6.7025, -6.8508,
-6.0698, -4.6566, -3.0047,
-1.4686, -0.3226, 0.3704, 0.6521, 0.6486, 0.5104, 0.3361, 0.1911, 0.0963, 0.0411, 0.0157, 0.0056, 0.001, 0.0004].

We observe that $g_{1}(-1)=0.02$ and $g_{1}(t)\leq 0$ for all
$t\in[-\sqrt{2}/2,1/2]$. Fig. 1 shows the graph of $g_{1}$ with normalization
$\tilde{g}_{1}(-1)=100$. Since there are negative coefficients $c_{k}$, we
cannot use the LP (Delsarte) bound. The SDP bound in Corollary 3.1 gives
$M=M_{1}:=22.5689$.

Let $B_{1}(N):=(N-M_{1})/N$. Then $B_{1}(25)=0.0324$ and $B_{1}(24)=0.0199$.

2. The coefficients $[c_{0},...,c_{22}]$ in the expansion of $g_{2}$ in $G_{k}^{(4)}$
= [0.222, 0.8648, 1.8875, 3.1425, 4.5059, 5.7052, 6.5739, 6.9286, 6.7119,
6.0157, 4.9575, 3.7767, 2.6446, 1.6914, 0.9947, 0.5249, 0.2524, 0.1097,
0.0409, 0.0153, 0.0042, 0.001, 0.0002].

Figure 1: Graphs of $\tilde{g}_{1}$, $\tilde{g}_{2}$ and $\tilde{p}_{4}$

We have $g_{2}(-1)=0.02$ and $g_{2}(t)\leq 0$ for all $t\in[-0.73,0.5]$. Fig.
1 shows the graph of $g_{2}$ with normalization $\tilde{g}_{2}(-1)=100$.

This case is very interesting. Note that all coefficients of $g_{2}$ are
positive and therefore we can apply bound (1.3):

$R_{g_{2}}(C)\geq LP(N):=c_{0}N-g_{2}(1)=0.222N-57.5714.$

On the other side the SDP bound in Corollary 3.1 gives

$R_{g_{2}}(C)\geq B_{2}(N)=\frac{N-M}{3N}=\frac{N-22.6452}{3N}$

Then for the most interesting cases $N=24$ and $N=25$ we have

$B_{2}(24)=0.0188>-52.2431=LP(24),\quad B_{2}(25)=0.0314>-52.0211=LP(25).$

###  4.2 A shorter proof of the kissing number problem in four dimensions.

In [17] we proved that $k(4)=24$. Let $t_{0}\in(-1,-0.5)$ and $f(t)$ be is a
nonnegative combination of $G_{k}^{(4)}$ with coefficients $c_{k}\geq 0$ such
that $f(t)\leq 0$ for all $t\in[t_{0},0.5)$. From [17, Theorem 1] follows that

$k(4)\leq\frac{1}{c_{0}}\max\\{h_{0},h_{1},...,h_{\mu},\\}$

where $\mu=A(4,\pi/3,\psi_{0}),\,\psi_{0}=\arccos|t_{0}|$, and $h_{m}$ is the
maximum of $H_{f}(Y)=f(1)+f(e_{1},y_{1})+...+f(e_{1},y_{m})$ over all
configurations of $m$ unit vectors $y_{j}$ in the spherical cap in
$\mathbb{S}^{3}$ given by $e_{1}\cdot y_{j}\leq t_{0}$ whose pairwise scalar
products are at most $\frac{1}{2}$.

We considered a polynomial $p_{4}$ (see Fig. 1) with $t_{0}=-0.6058$. In this
case $\mu=6$.

Technically, the most difficult case turned out to be m=5. This case takes up
a significant part of the proof. We tried to exclude this case and reduce
$t_{0}$. However, numerous experiments with an extended LP bound did not lead
to success.

######  Theorem 4.1.

$k(4)=24$

###### Proof.

The 24–cell is an example of a kissing arrangement. Then $k(4)\geq 24$. It
remains to prove $k(4)<25$. Assume the converse: $N\geq 25$.

Let $C$ be an $(25,4,\pi/3)$ spherical code. Since $g_{1}(t)\leq 0$ for all
$t\in[-\sqrt{2}/2,1/2]$, we have $t_{0}=-\sqrt{2}/2$ and $\mu=4$. Using the
same method as in [17] we consider the cases $\mu=0,1,2,3,4$ to find the
maximum of $R_{g_{1}}(C)=\sum_{t}A_{t}g_{1}(t)$. This maximum is achieved at
$\mu=2$ and is $0.0266$, i.e. $R_{g_{1}}(C)<0.0266$. On the other side we have
$R_{g_{1}}(C)>B(25)=0.0324$, a contradiction. ∎

###  4.3 New bounds for the distance distribution on a $(24,4,\pi/3)$
spherical code

Now we consider a maximal kissing arrangement in dimension 4 that is a
$(24,4,\pi/3)$ spherical code. The long-standing open uniqueness conjecture on
this code states that this arrangement is isometric to the 24–cell. In fact,
the uniqueness conjecture is equivalent to the following

Conjecture. Let $C$ be a $(24,4,\pi/3)$ spherical code. Then

$A_{-1}=1,\quad A_{-1/2}=8,\quad A_{0}=6,\quad A_{1/2}=8,\ \quad
A_{t}=0,\,t\neq\pm 1,\pm 1/2,0.$ $None$

Note that in this dimension the equality $A_{-1}=1$ yields (3.1), see [7].
Moreover, in [13] we show that every equality in (3.1) implies all other
equalities.

Recently, using Corollaries 1 and 2 from [20] we found several bounds on the
distance distribution of a kissing arrangements in four dimensions [13]. For
several intervals these bounds are sharp.

Let $S\subset[-1,0.5]$. Denote

$A(S):=\sum\limits_{t\in S:A_{t}>0}{A_{t}}.$

######  Theorem 4.2.

Let $C$ be a $(24,4,\pi/3)$ – spherical code. Then

$A([-1,-0.45])\leq 9;\quad A([-1,0.05])\leq 15,\quad A([-0.55,0.05])\leq
14,\quad A([-0.05,0.5]\leq 14,$ $A([-1,-0.73])\geq 1,\quad A([0.35,0.5])\geq
8$

##  5 Concluding Remarks

In conclusion, we outline some applications of Theorems 2.2 and 3.1 and their
generalizations.

###  5.1 Towards a proof of the uniqueness conjecture

We know that $k(4)=24$ [17]. However, in dimension 4 the uniqueness of the
maximal kissing arrangement is conjectured to be the 24–cell but not yet
proven. Equivalently, the uniqueness conjecture is the following:

Denote by $s_{d}(n)$ the optimal SDP bound on $k(n)$ given by (3) with
$\deg(F)=d$ (see [16]). In the following table it is shown that this
minimization problem is a semidefinite program and that every upper bound on
$s_{d}(4)$ provides an upper bound for the kissing number in dimension 4.

  * •

$s_{7}(4)<24.5797$ – Bachoc & Vallentin [1];

  * •

$s_{11}(4)<24.10550859$ – Mittelmann & Vallentin [16];

  * •

$s_{12}(4)<24.09098111$ [16];

  * •

$s_{13}(4)<24.07519774$ [16];

  * •

$s_{14}(4)<24.06628391$ [16];

  * •

$s_{15}(4)<24.062758$ – Machado & de Oliveira Filho [15];

  * •

$s_{16}(4)<24.056903$ [15].

This table shows that $s_{d}$ with $d>12$ is relatively close to 24,
$s_{d}-24<2/N=1/12$. We think that our approach which is based on Theorems 2.2
and 3.1 can help to prove the uniqueness conjecture.

###  5.2 Towards a proof of the 24-cell conjecture

The sphere packing problem asks for the densest packing of ${\mathbb{R}}^{n}$
with unit balls. In four dimensions, the old conjecture states that a sphere
packing is densest when spheres are centered at the points of lattice $D_{4}$,
i.e. the highest density $\Delta_{4}$ is $\pi^{2}/16$, or equivalently the
highest center density is $\delta_{4}=\Delta_{4}/B_{4}=1/8$. For lattice
packings, this conjecture was proved by Korkin and Zolatarev in 1872.
Currently, for general sphere packings the best known upper bound for
$\delta_{4}$ is $0.130587$, a slight improvement on the Cohn–Elkies bound of
$\delta_{4}<0.13126$, but still nowhere near sharp.

In [19] we considered the following conjecture:

The 24–cell conjecture. Consider the Voronoi decomposition of any given
packing $P$ of unit spheres in ${\mathbb{R}}^{4}$. The minimal volume of any
cell in the resulting Voronoi decomposition of $P$ is at least as large as the
volume of a regular 24–cell circumscribed to a unit sphere.

Note that a proof of the 24-cell conjecture also proves that $D_{4}$ is the
densest sphere packing in 4 dimensions.

In [18, Sect. 4] and [19, 3.3] we considered polynomials $H_{k}$ that are
positive–definite in ${\mathbb{R}}^{n}$. Actually, $H_{k}$ are polynomials
that extend the Bachoc–Vallentin polynomials $S_{k}$. It is an interesting
problem to find generalizations of Theorems 2.2 and 3.1 for sphere packings in
${\mathbb{R}}^{n}$. Perhaps, these bounds for $n=4$ can help to prove the
24–cell conjecture.

###  5.3 SDP bounds for Thomson’s and related problems

The objective of the Thomson problem is to determine the minimum electrostatic
potential energy configuration of $N$ electrons constrained to the surface of
a unit sphere. More generally, a configuration is called $h$-optimal for a
potential interaction $h:[-1,1]\to\mathbb{R}$, if it minimizes the $h$-energy.

The logarithmic potential, Thomson, Newton potential, and more generally the
Riesz potential as well as as well as the Gaussian potential have been well
studied in the literature [6]. All of these potentials are absolutely monotone
potentials. Cohn and Kumar invented universally optimal configurations, namely
they minimize the energy for all absolutely monotone potentials $h$ [8].

The regular simplex, the cross polytope and the so-called isotropic spaces are
the only known classes of are universally optimal configurations, (see [8,
Table 1]). All other known optimal configurations in the literature, even when
the interacting potential $h$ is fixed, have particular values of the
dimension $d$ and the cardinality $N$, see the fundamental monograph on this
topic [6]. For instance, the Thomson problem is solved only for $N\leq 6$ and
$N=12$.

Note that one of the most powerful tool for lower bounding $h$-energy is the
Delsarte–Yudin linear programming method. This method was applied to most of
the configurations in [8, Table 1] to prove that they are universally optimal.

Let a potential $h$ is given. Let $g(t)\leq h(t)$ on some interval $T$. Then
Theorem 2.2 can be considered as a Yudin type bound the minimum energy. It is
a very interesting task to consider various cases of optimal configurations
for known potentials and see what kind of bounds can be obtained using this
method.

###  5.4 SDP bound for contact graphs and the Tammes problem

The following problem was first asked by the Dutch botanist Tammes in 1930:
Find the largest angular separation $\theta$ of a spherical code $C$ in
${\mathbb{S}}^{2}$ of cardinality $N$.
In other words,
How are $N$ congruent, non-overlapping circles distributed on the sphere when
the common radius of the circles is as large as possible?

The Tammes problem is presently solved for only a few values of $N$: for
$N=3,4,6,12$ by L. Fejes Tóth; for $N=5,7,8,9$ by Schütte and van der Waerden;
for $N=10,11$ by Danzer; for $N=24$ by Robinson.; and for $N=13,14$ by Musin &
Tarasov [21, 22].

The computer-assisted solution of Tammes’ problem for $N=13$ and $N=14$
consists of three parts: (i) creating the list $L_{N}$ of all planar graphs
with $N$ vertices that satisfy the conditions of [22, Proposition 3.1]; (ii)
using linear approximations and linear programming to remove from the list
$L_{N}$ all graphs that do not satisfy the known geometric properties of the
maximal contact graphs [22, Proposition 3.2]; (iii) proving that among the
remaining graphs in $L_{N}$ only one is maximal.

In fact, the list $L_{N}$ consists of a huge number of graphs. (For $N=13$ it
is about $10^{8}$ graphs.) We think that this paper can help to reduce the
number of graphs in $L_{N}$.

###  5.5 Generalization of the $k$–point SDP bound for spherical codes

In [18] we invented the $k$–point SDP bound for spherical codes. Note that for
$k=2$ that is the classical Delsarte bound. The 3–point SDP bound was first
considered by Bachoc and Vallentin [1]. Recently, this method with $k=4,5,6$
was apply for upper bounds of the maximum number of equiangular lines in $n$
dimensions [12]. It is an interesting problem to find generalizations of
results in this paper using the $k$–point SDP bounds and apply these bounds
for $s$–distance sets and equiangular lines.

Acknowledgments. I am very grateful to Maria Dostert, Alexander Kolpakov and
Philippe Moustrou for helpful discussions and useful comments. We spent a lot
of time together discussing ideas, algorithms and current results on the
uniqueness problem of kissing arrangements. We have encouraging results and I
hope that we can complete our paper [13].

## References

  * [1] C. Bachoc and F. Vallentin.  New upper bounds for kissing numbers from semidefinite programming.  J. Amer. Math. Soc., 21(3): 909–924, 2008.
  * [2] C. Bachoc and F. Vallentin.  Optimality and uniqueness of the $(4,10,1/6)$ spherical code.  J. Comb. Theory, Ser. A, 116(1): 195–204, 2009.
  * [3] C. Bachoc and F. Vallentin.  Semidefinite programming, multivariate orthogonal polynomials, and codes in spherical caps.  European J. Combin., 30(3): 625–637, 2009.
  * [4] E. Bannai and N. J. A. Sloane. Uniqueness of certain spherical codes, Canadian J. Math., 33: 437–449, 1981.
  * [5] A. Barg and O. R. Musin.  Codes in spherical caps.  Adv. Math. Comm., 1(1): 131–149, 2007.
  * [6] S. Borodachov, D. Hardin, E. Saff, _Discrete Energy on Rectifiable Sets_ , Springer Monographs in Mathematics, Springer, 2019.
  * [7] P.G. Boyvalenkov. Nonexistence of certain symmetric spherical codes, Designs, Codes and Cryptography, 3, 69–74, 1993.
  * [8] H. Cohn and A. Kumar, Universally optimal distribution of points on spheres, J. of AMS, 20 (2007), 99–148.
  * [9] H. Cohn, D. de Laat, A. Salmon, Three–point bounds for sphere packing, preprint, arXiv:2206.15373
  * [10] H. Cohn and J. Woo.  Three–point bounds for energy minimization.  J. Amer. Math. Soc., 25: 929–958, 2012.
  * [11] J.H. Conway and N.J.A. Sloane, Sphere Packings, Lattices, and Groups, New York, Springer-Verlag, 1999 (Third Edition).
  * [12] D. de Laat, F.C. Machado, F. M. de Oliveira Filho and F. Vallentin. k–point semidefinite programming bound for equiangular lines, Math. Program. 194, 533–567, 2022.
  * [13] M. Dostert, A. Kolpakov, P. Moustrou, and O. R. Musin. Uniqueness of spherical codes via SDP, in preparation.
  * [14] V.I. Levenshtein, On bounds for packing in $n$-dimensional Euclidean space, Sov. Math. Dokl. 20 (2): 417–421, 1979.
  * [15] F. C. Machado and F. M. de Oliveira Filho.  Improving the semidefinite programming bound for the kissing number by exploiting polynomial symmetry,  Experimental Math., 27: 362–369, 2018.
  * [16] H. D. Mittelmann and F. Vallentin,  High-accuracy semidefinite programming bounds for kissing numbers,  Experimental Math., 9: 175–179, 2010.
  * [17] O. R. Musin.  The kissing number in four dimensions.  Ann. of Math., 168: 1–32, 2008.
  * [18] O. R. Musin, Multivariate positive definite functions on spheres, in: Discrete Geometry and Algebraic Combinatorics, AMS Series: Contemporary Mathematics, 625: 177–190, 2014.
  * [19] O. R. Musin. Towards a proof of the 24-cell conjecture, Acta Math. Hungar., 155(1): 184–199, 2018.
  * [20] O. R. Musin.  An extension the semidefinite programming bound for spherical codes.  preprint, arXiv:1903.05767
  * [21] O. R. Musin and A. S. Tarasov.  The strong thirteen spheres problem.  Discrete Comput. Geom., 48(1): 128–141, 2012.
  * [22] O. R. Musin and A. S. Tarasov, The Tammes problem for N=14. Experimental Math., 24: 460–468, 2015.
  * [23] A.M. Odlyzko and N.J.A. Sloane, New bounds on the number of unit spheres that can touch a unit sphere in $n$ dimensions, J. of Combinatorial Theory A 26: 210–214, 1979.
  * [24] V. A. Yudin, Minimal potential energy of a point system of charges. Discret. Mat. 4, 115–121 (1992) (in Russian). English translation: Discr. Math. Appl., 3, 75–81 (1993)

O. R. Musin, School of Mathematical and Statistical Sciences, University of
Texas Rio Grande Valley, One West University Boulevard, Brownsville, TX,
78520.

E-mail address<EMAIL_ADDRESS>


# Oscillating behavior within the social force model

Mohcine Chraibi  Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH.
52425 Jülich, Germany<EMAIL_ADDRESS>
##  1\. Introduction

The social force model (SFM) [4] belongs to a class of microscopic force-based
pedestrian model for which the interaction with the neighbors solely depends
on the distance spacings. Yet, distance based models can lead to unrealistic
oscillating behaviors with collision and negative speed (especially in 1D, see
for instance [12] and the harmonic oscillator equation). With oscillation we
describe a situation where pedestrians perform, in the direction of the
intended movement, forward and backward movement. This behavior is to be
distinguished from the oscillations in the passing direction in bottlenecks
[5] and from lateral oscillations along the shoulder direction [6, 9].

The problem of oscillations in the movement direction was investigated in the
literature and shown to rise in general in simulations with force-based models
[1, 7, 3]. In [2] the model-induced oscillations were quantified by
introducing the “oscillation-proportion” introduced. Furthermore, it was
suggested that this phenomena is the dual-problem of another model-induced
model, namely the overlapping problem.

In this article we show on basis of analytical investigation of a simplified
one-dimensional scenario, that only particular values of the parameters allow
collision-free and non-oscillating simulation. However, ‘realistic” values for
the (physical) parameters of the SFM do not produce such simulations.

##  2\. Local stability analysis

The social force model in one dimension and with a coupled interaction with
the predecessor writes:

(1) $\ddot{x}_{n}=-a\exp\Big{(}\frac{-\Delta
x_{n}}{b}\Big{)}+\frac{v_{0}-\dot{x}_{n}}{\tau},$

with $v_{0},a,b,\tau>0$. Hereby pedestrians are presented as points and
$\Delta x_{n}=x_{n+1}-x_{n}$.

Equation (1) can be linearly approximated for $\Delta x_{n}/b\approx 0$ as
follows

(2) $\ddot{x}_{n}\approx-a\Big{(}1+\frac{-\Delta
x_{n}}{b}\Big{)}+\frac{-\dot{x}_{n}}{\tau}+\tilde{v}_{0},\;\;\;\tilde{v}_{0}=\frac{v_{0}}{\tau},$

or

(3) $\ddot{x}_{n}+\frac{\dot{x}_{n}}{\tau}-\frac{a}{b}\Delta
x_{n}+a-\tilde{v}_{0}\approx 0.$

Using the following variable substitution $\xi=-\frac{a}{b}\Delta
x_{n}+a-\tilde{v}_{0},$ and assuming that $x_{n+1}$ is constant we get

(4) $\frac{b}{a}\ddot{\xi}+\frac{b}{a\tau}\dot{\xi}+\xi=0,$

which can be brought in the form of a harmonic oscillator

(5) $\ddot{\xi}+r\dot{\xi}+\omega_{0}^{2}\xi=0,$

with $r=\frac{1}{\tau}$ and $\omega_{0}^{2}=\frac{a}{b}$.

Since the coefficients $r$ and $\omega_{0}$ are constant, (5) will have a
solution of the form $\xi=e^{\lambda t}$ which yields the following
characteristic polynomial:

(6) $\lambda^{2}+r\lambda+\omega_{0}^{2}=0.$

The roots of (6) are

(7) $\lambda_{1,2}=\frac{-r\pm\sqrt{r^{2}-4\omega_{0}^{2}}}{2},$

which gives a general solution of the form

$\xi=c_{1}e^{\lambda_{1}}+c_{2}e^{\lambda_{2}}.$

The system described by (1) does not oscillate if the imaginary part of the
solutions are nil, i.e. if $r^{2}-4\omega_{0}^{2}>0$ or

(8) $\frac{a}{b}\tau^{2}<\frac{1}{4}.$

Note that the solution does not depend $explicitly$ on $v_{0}$.

##  3\. Interpretation

In [13] is has been reported that evaluation of empirical data yields
$\tau=0.61$ s. A slightly different value of $\tau$ was measured in [10]
($\tau=0.54\pm 0.05$ s). It simulation $\tau=0.5$ s is often used as a
realistic value.

In the following examples, we always assume that $v_{0}=1.2$ m/s and
$\tau=0.5$ s. The parameter values in [4, 11] $a=2000$ N and $b=0.08$ m lead
to oscillations. See Fig. 1

Figure 1. Simulation with $a=2000$ N and $b=0.08$ m.

The model [8] with $a=300$ N and $b=0.5$ m also oscillates. See Fig. 2.

Figure 2. Simulation with $a=300$ N and $b=0.5$ m.

In [10] an experimental evaluation leads to $a=4.5\pm 0.3$, $b=1.25$ m (mean
value) and $\tau=0.54\pm 0.05$ s. That means that the system according to Eq.
(8) is oscillating.

##  4\. Conclusion

The presented results show that the social force model has to be extended by
using velocity terms in order to produce realistic collision-free behaviors
with reasonable values for the parameters.

In its original form without velocity dependency (reasonable) values of
parameters lead inevitably to erroneous oscillating movement.

## References

  * [1] M. Chraibi, S. A. Kemloh, U., and A. Schadschneider.  Force-based models of pedestrian dynamics.  Networks and Heterogeneous Media, 6(3):425–442, 2011.
  * [2] M. Chraibi, A. Seyfried, and A. Schadschneider.  Generalized centrifugal force model for pedestrian dynamics.  Physical Review E, 82:046111, 2010.
  * [3] F. Dietrich, G. Köster, M. Seitz, and I. Sivers.  Bridging the gap: From cellular automata to differential equation models for pedestrian dynamics.  In R. Wyrzykowski, J. Dongarra, K. Karczewski, and J. WaÅ›niewski, editors, Parallel Processing and Applied Mathematics, Lecture Notes in Computer Science, pages 659–668. Springer Berlin Heidelberg, 2014.
  * [4] D. Helbing, I. Farkas, and T. Vicsek.  Simulating dynamical features of escape panic.  Nature, 407:487–490, 2000.
  * [5] D. Helbing, P. Molnar, I. J. Farkas, and K. Bolay.  Self-organizing pedestrian movement.  Environment and Planning B, 28:361–383, 2001.
  * [6] S. P. Hoogendoorn and W. Daamen.  Pedestrian behavior at bottlenecks.  Transport. Sci., 39(2):147–159, 2005.
  * [7] G. Köster, F. Treml, and M. Gödel.  Avoiding numerical pitfalls in social force models.  Physical Review E, 87, 2013.
  * [8] T. I. Lakoba, D. J. Kaup, and N. M. Finkelstein.  Modifications of the helbing-molnar-farkas-vicsek social force model for pedestrian evolution.  Simulation, 81:339–352, 2005.
  * [9] X. Liu, W. Song, and J. Zhang.  Extraction and quantitative analysis of microscopic evacuation characteristics based on digital image processing.  Physica A: Statistical Mechanics and its Applications, 388(13):2717–2726, jul 2009.
  * [10] M. Moussaïd, D. Helbing, S. Garnier, A. Johansson, M. Combe, and G. Theraulaz.  Experimental study of the behavioural mechanisms underlying self-organization in human crowds.  Proc. R. Soc. B., 276(1668):2755–2762, 2009.
  * [11] D. R. Parisi, M. Gilman, and H. Moldovan.  A modification of the social force model can reproduce experimental data of pedestrian flows in normal conditions.  Physica A: Statistical Mechanics and its Application, 388(17):3600–3608, sep 2009.
  * [12] M. Treiber and A. Kesting.  Traffic Flow Dynamics.  Springer, Berlin, 2013.  ISBN: 978-3-642-32459-8.
  * [13] T. Werner and D. Helbing.  The social force pedestrian model applied to real life scenarios.  In E. R. Galea, editor, Pedestrian and Evacuation Dynamics 2003, pages 17–26. CMS Press, London, 2003.


# Machine Learning Assisted Characterization of Labyrinthine Pattern
Transitions

Kotaro Shimizu  Department of Physics, University of Virginia,
Charlottesville, Virginia 22904, USA  Department of Applied Physics, The
University of Tokyo, Tokyo 113-8656, Japan  Vinicius Yu Okubo  Dept.
Electronic Systems Engineering, Polytechnic School, University of São Paulo,
Brazil  Rose Knight  Department of Physics, University of Virginia,
Charlottesville, Virginia 22904, USA  Ziyuan Wang  Department of Physics,
University of Virginia, Charlottesville, Virginia 22904, USA  Joseph Burton
Department of Physics, University of Virginia, Charlottesville, Virginia
22904, USA  Hae Yong Kim  Dept. Electronic Systems Engineering, Polytechnic
School, University of São Paulo, Brazil  Gia-Wei Chern  Department of Physics,
University of Virginia, Charlottesville, Virginia 22904, USA  B. S. Shivaram
Department of Physics, University of Virginia, Charlottesville, Virginia
22904, USA

###### Abstract

We present a comprehensive approach to characterizing labyrinthine structures
that often emerge as a final steady state in pattern forming systems. We
employ machine learning based pattern recognition techniques to identify the
types and locations of topological defects of the local stripe ordering to
augment conventional Fourier analysis. A pair distribution function analysis
of the topological defects reveals subtle differences between labyrinthine
structures which are beyond the conventional characterization methods. We
utilize our approach to highlight a clear morphological transition between two
zero-field labyrinthine structures in single crystal Bi substituted Yttrium
Iron Garnet films. An energy landscape picture is proposed to understand the
athermal dynamics that governs the observed morphological transition. Our work
demonstrates that machine learning based recognition techniques enable novel
studies of rich and complex labyrinthine type structures universal to many
pattern formation systems.

Labyrinthine structures are ubiquitous in out-of-equilibrium nonlinear systems
ranging from biological and chemical reactions to fluid convection, crystal
growth, and magnetic ordering [1, 2, 3, 4, 5]. In such pattern forming
systems, the complex structures emerge as a result of energy injection and
dissipation in a highly nonlinear process. The labyrinthine patterns are
generally characterized by stripe domains of different orientations, sizes,
and grain-boundary structures. The predominance of periodic stripes indicates
breaking of translational symmetry locally. Yet, contrary to long-range
ordered states in an equilibrium phase transition, labyrinthine patterns are
essentially disordered and cannot be described by a well-defined order
parameter. Indeed, labyrinthine structures can be viewed as an intermediate
between featureless short-range correlated glassy state and long-range ordered
stripe or crystalline phases [6].

Despite their prevalence in pattern forming systems, a complete
characterization of labyrinthine structures is still lacking [7]. A defining
characteristic of labyrinthine patterns is the ring-like feature in its
structure factor obtained from conventional Fourier analysis [8]. The radius
and width of the ring correspond to the wavelength of local stripes and
characteristic size of stripe domains, respectively [8, 9, 10, 11]. While such
global Fourier analysis provides a basic characterization of labyrinths, it
fails to capture subtle differences of labyrinthine patterns which have
important structural or dynamical implications. Other useful measures, such as
the disorder functions [12, 13, 14], have been introduced to quantify
deviations from a perfect stripe order. Another important characterization
often employed is the density of topological defects of labyrinthine
structures [15, 16, 17, 18]. Indeed, the distribution and correlation between
topological defects, such as disclinations and dislocations, of the stripe
order encode important information about the labyrinths [17, 18]. However,
efficient and accurate identifications of such point-like defects in large-
scale experimental or simulation data remain a challenging task.

In this paper, we present a comprehensive framework for the characterization
of labyrinthine structures by combining Fourier analysis with machine learning
(ML) based template recognition methods. A two-step algorithm that consists of
rotation-invariant template matching followed by convolutional neural network
analysis was developed to precisely identify topological defects (junctions
and terminals) and their coordinates in the images. The high-precision real-
space configuration data allow us to compute pair distribution functions of
the topological defects, which provide valuable information about labyrinthine
patterns that are complementary to those extracted from structure factors.

We apply our approach to studying an intriguing nonequilibrium morphological
phase transition in the magnetic medium Ytrium Iron Garnet (YIG), a well known
technologically important ferromagnetic material. Magnetic interactions of YIG
are dominated by a short-range ferromagnetic exchange and long-range dipolar
interactions [19]. YIG films doped with bismuth are known to carry a strong
perpendicular magnetic anisotropy and form complex labyrinthine patterns [16,
17, 18]. Bismuth doping introduces a number of advantages desirable for
magneto-optic isolators and sensors [20, 21], enhances the Faraday rotation,
improves the perpendicular anisotropy and also lowers the saturation
magnetization [19, 20, 21, 22]. For films grown under appropriate conditions,
the magnetization can be saturated in fields less than 100 Oe. Thus the
labyrinthine stripe domain patterns are easily observed with a small
electromagnet in polarized light under a microscope [23].

Figure 1:  (a) Schematic figure of the external magnetic field stepping
protocol. The sequence starting from the positive (negative) magnetic field is
denoted by the red (blue) line. Domain images of (b) the quenched state at
step 0 and (c) the annealed state at step 36. The left and right insets in (c)
display the structures of a terminal and a junction, respectively.

Specifically, here we report an intriguing pattern transition in "annealing"
experiments performed on YIG films where a perpendicular magnetic field is
stepped down to zero starting from a fully saturated state. The final zero-
field state is reached with each subsequent step reduced from the previous
following an exponential decay pattern. These domain annealing experiments are
carried out with a specific protocol. Starting from the sample in the fully
magnetized state, we instantly drop the field to zero where it is held for 10
seconds during which an image is acquired. Further sequential de-magnetization
is carried out but with the height of each step reduced exponentially and with
the reversed field directions as indicated in Fig. 1(a). All measurements
reported here were performed at room temperature. The images obtained covered
an area of 2 mm $\times$ 1.8 mm. Two samples grown under similar conditions
were studied.

Representative images obtained during the beginning of the protocol and at the
end of the protocol, are shown in Fig. 1(b) and Fig. 1(c), respectively. Both
of them exhibit the labyrinthine stripe patterns and consist of a plethora of
defects depicted in Fig. 1(c). These defects where the dark domain ends and
three domains meet are termed terminals ($+1/2$ disclinations) and junctions
($-1/2$ disclinations), respectively, and are topological defects associated
with rotational symmetry breaking [16]. When introducing one of these defects
into perfectly ordered stripes, it inevitably gives rise to the appearance of
the other. In other words, the disorder manifests as a disclination pair.

Figure 2:  (a) Structure factor of the image $\mathcal{S}({\bf q})$ for (left)
the quenched state with step $0$ and (right) the annealed state with step
$36$. (b) The angle-averaged structure factor, $I(q)$, for different steps
represented by different colors. (c) Evolution of the peak width $\sigma$ and
the peak wave number $q_{0}$ for $I(q)$ represented by the red and the blue
circles, respectively.  Figure 3:  The schematic illustration of the algorithm
to detect the defects. First, we detect candidates for the defects (red and
blue dots) by using the rotation-invariant template matching algorithm. The
identified candidates in this process, while many are terminals and junctions,
also encompass structures that are not defects. Then, these candidates are
classified into terminal, junction, and false detection through the deep
learning filter.

The two (both in zero field) images show labyrinthine structures of similar
domains, but appear visually different at the same time, making it challenging
to quantify their morphology. For further analysis and discussion, we
distinguish the domain patterns in Figs. 1(b) and 1(c) as the quenched and
annealed states, respectively. In the quenched state, the border of bright and
dark domains exhibit a sinuous nature and do not appear as parallel. In
contrast, the annealed state consists of significant regions of very parallel
domains. This state exhibits roughly equal widths of dark and bright domains,
and the areas occupied by them are also roughly equal for any given sample
region. This net approximate equality can be achieved by a variety of
different patterns which more or less appear the same.

To characterize these morphologically distinct states, we first perform the
conventional Fourier analysis by calculating the structure factor

$\displaystyle\mathcal{S}({\bf q})=\Braket{\left
missing}{w_{\epsilon,n,m}({\bf q})}{{}^{2}}_{\epsilon,n},$ (1)

where $w_{\epsilon,n,m}({\bf q})$ is the Fourier transform of the image at
step $m$ ($m=0,1,\ldots,36$) of the $n$-th trial ($n=1,2,\ldots,6$) under the
field stepping protocol starting from the positive ($\epsilon=1$) and negative
($\epsilon=-1$) field, and $\braket{\cdots}_{\epsilon,n}$ represents the
average over $\epsilon$ and $n$. In Fig. 2(a), we show the structure factors,
manifesting circular patterns due to the isotropic nature of the system. The
circular pattern is smaller and thicker in the quenched state than in the
annealed state. We also computed the angle-averaged structure factor $I(q)$,
shown in Fig. 2(b). As discussed above, the demagnetization brings about the
larger radius of the circular pattern, which becomes apparent as the peak $q$
moves to the right. To capture the intricate evolution of the peak structure,
we calculated the peak wave number $q_{0}$ and peak width $\sigma$ obtained by
fitting $I(q)$ with a Gaussian given by
$A\exp\left(-(q-q_{0})^{2}/\left(2\sigma^{2}\right)\right)$.

The step dependence of $\sigma$ and $q_{0}$ is shown in Fig. 2(c), exhibiting
the decrease of $\sigma$ and the increase of $q_{0}$ with steps. Notably,
$\sigma$ and $q_{0}$ exhibit sharp transitions at distinct steps,
specifically, in the range of 6 to 12 steps represented by the yellow area and
10 to 20 steps represented by the green area, respectively. It is also clearly
seen that when the step height falls below 12 gauss (between steps 10 and 11),
the domains settle down to the annealed state. Considering that $1/\sigma$ and
$2\pi/q_{0}$ correspond to the correlation length and the magnetic modulation
period, respectively, we can infer that through the demagnetization process,
the labyrinthine structure, which is initially less compact in the quenched
state, aligns itself to increase the correlation length and subsequently
shortens the magnetic period, transitioning into a more compact annealed
state. Application of further smaller steps reorganizes the domain pattern
locally at various points but leaves the overall features intact.

The conventional analysis conducted thus far provides a certain degree of a
global picture in the process of domain compacting. However, it fails to yield
a comprehensive understanding of the detailed response of topological defects.
As the evolution of the labyrinth patterns is mediated by these defects
introduced in Fig. 1, the precise identification of their number and
coordinates is inevitable for further clarification.

Nowadays, we have high-resolution digital cameras and much greater
computational power compared to what was used in the pioneering work of
decades ago [16, 17, 18]. Modern popular algorithms for object detection, such
as Faster R-CNN [24] and YOLO [25], are based on convolutional neural networks
[26]. However, they are not well-suited for defect detection for two reasons.
First, they are designed to detect a few large objects, whereas in our
context, the goal is to detect thousands of small, closely clustered objects.
Second, they require a large number of manually annotated training images,
which is extremely labor-intensive to produce due to the high quantity of
defects in each image.

We developed a two-step detection algorithm capable of finding thousands of
small objects with little manual annotation, which is shown in Fig. 3. First,
we employ rotation-invariant template matching to identify candidate points
for junctions and terminals [27, 28]. This step is executed with a threshold
set to avoid any false negatives, even at the cost of generating several false
positives. This approach significantly reduces the annotation workload, as we
now only need to review the candidate points rather than locating and
classifying every defect in the image. We selected some images containing
candidate points, manually corrected the false positives, and used them to
train a convolutional neural network classifier. Therefore, in the second
step, the candidate points from the first stage are processed through the
classifier, which differentiates between terminals, junctions, and false
detections. Manual verification of numerous images indicates that the
algorithm’s detection accuracy is nearly 100%. More details about the proposed
algorithm will be provided elsewhere.

Figure 4:  (a) The radial distribution function for the positions of terminals
$g_{\rm t}(r)$ for step 0 and step 36, denoted by the dashed and solid lines,
respectively. (b) Evolution of (top) the effective coordination number for
terminals $n^{\rm eff}_{\rm t}$ and (bottom) the normalized density of
terminals $\tilde{\rho}_{\rm t}$.

From the positions of topological defects, we calculate the radial
distribution function for terminals $g_{\rm t}(r)$ given by

$\displaystyle g_{\rm t}(r)=\frac{1}{2\pi rN_{\rm t}\rho_{\rm
t}}\sum_{ij}f(r-|{\bf r}_{i}-{\bf r}_{j}|),$ (2)

where $N_{\rm t}$ is the number of terminals, $\rho_{\rm t}$ is the density of
terminals, $f(r)=\exp\left(-r^{2}/2\right)/\sqrt{2\pi}$, and ${\bf r}_{i}$ is
the position of the $i$-th terminal in the image. In Fig. 4(a), $g_{\rm t}(r)$
is presented for the quenched state at step 0 and the annealed state at step
36. Both instances of $g_{\rm t}(r)$ exhibit peaks around $r\sim 60$,
corresponding to the magnetic modulation period $\frac{2\pi}{q_{0}}$, which is
slightly larger in the quenched state. For $r\gtrsim 60$, $g_{\rm t}(r)$ for
the quenched state gradually approaches unity, whereas, in the annealed state,
peaks are also observed at $r\sim 95$ and $150$. In Fig. 4(b), in conjunction
with the computed $g_{\rm t}(r)$, we illustrate the effective coordination
number defined by

$\displaystyle n_{\rm t}^{\rm eff}=\Braket{\int_{0}^{r_{0}}2\pi r\rho_{\rm
t}g_{\rm t}(r)dr}_{n},$ (3)

where $r_{0}$ is the first bottom in $g_{\rm t}(r)$, along with the normalized
density of terminals given by 111 The characteristic length scale determined
by the stripe period decreases with increasing steps and hence, the number of
defects in the same area increases in general. To quantify the number of
defects in this length scale, we rescale $\rho_{\rm t}$ by the factor of
$\left(\frac{2\pi}{q_{0}}\right)^{2}$

$\displaystyle\tilde{\rho}_{\rm t}=\rho_{\rm
t}\left(\frac{2\pi}{q_{0}}\right)^{2}.$ (4)

Here, these results are shown after taking the average over all the trials
with $\epsilon=1$.

The effective coordination number signifies the number of nearest neighboring
terminals and $\tilde{\rho}_{\rm t}$ denotes the density of terminals
concerning the stripe length scale. Through the demagnetization process, both
$n_{\rm t}^{\rm eff}$ and $\tilde{\rho}_{\rm t}$ monotonically decrease,
whereas $n_{\rm t}^{\rm eff}$ is more significantly diminished. It is
noteworthy that $n_{\rm t}^{\rm eff}$ exhibits conspicuous variations in the
range of 6 to 12 steps denoted by the yellow-shaded area, similar to $\sigma$,
whereas $\tilde{\rho}_{\rm t}$ displays significant changes in the range of 10
to 20 steps represented by the green-shaded area, akin to $q_{0}$ [see also
Fig. 2(c)]. These results imply a two-stage evolution of the magnetic
structure: In the initial stage, the magnetic structure undergoes changes by
moving neighboring terminals apart while hardly altering the number of
terminals. The defects reduce their nearest-neighbor coordination number while
increasing the further-neighbor coordination number. This can be interpreted
as the topological defects working to restore a positional order, leading to
an increase in the correlation length. Meanwhile, in the later stage, the
density of defects decreases while maintaining the coordination number. By
shortening the period of the stripes and the density of the defects, the
system can gain the magnetostatic energy.

Figure 5:  Schematic illustration of the demagnetization process. The black
solid curve represents the energy landscape of the system as a function of an
abstract one-dimensional coordinate representing the spin configuration. The
magnetic structure is represented at points within a complex potential
landscape, exhibiting the transitions between different valleys through the
demagnetization process, denoted by the red and blue arrows.

YIG is renowned for exhibiting a perfect stripe phase in zero magnetic fields
(e.g., Fig. 6 of Ref. [18]). However, a labyrinthine structure is observed
within the time scale observed in the present study, and even after undergoing
a demagnetization process with alternating magnetic fields, the system falls
into a more compact labyrinthine pattern without achieving uniform stripe
alignment. This behavior implies the presence of a plethora of local minima in
the energy landscape, as depicted in Fig. 5. The demagnetization process, by
promoting transitions between these local minima, effectively bridges
qualitatively distinct quenched and annealed states. In systems governed by
such a complex energy landscape, nonlinear responses are expected, and in
recent years, they have much attention as a stage for physical reservoir
computing. The labyrinthine pattern in magnets, whose complexity can be
modulated by an external magnetic field, is expected to be one of the
promising venues for physical reservoir computing.

To summarize, we have identified a morphological transition of magnetic
domains in a technologically important system and characterized their
evolution in new ways by using modern machine learning and pattern recognition
analysis. Although similar morphological transitions have been noted before,
the present study involves a very subtle transition observed during a
demagnetization process. This points to the richness of the magnetic domains
in YIG (and other systems) but also is a testimony to the accuracy with which
its features can be accurately and quantitatively defined. This has
facilitated a detailed analyses about the positional correlation of the
defects thus resolving the complex evolution dominating the transition. This
is possible due to the high-resolution images and the fidelity of feature
recognition and quantification from machine learning tools. These versatile
tools will undoubtedly find applications in studies of diverse labyrinthine
structures.

###### Acknowledgements.

We thank Dr. Vincent Fratello and Dr. Shanthi Subramanian of Lucent-Bell Labs
for the supply of the LPE grown epitaxial thin films used in this work. We are
also grateful to David Huse of Princeton for many helpful comments. The
experimental work was partially funded by the NSF under grant DMR and by the
3Cavalier grant from the University of Virginia. This work was also supported
by JSPS KAKENHI Grant Number No. JP21J20812. Gia-Wei Chern was partially
supported by the US Department of Energy Basic Energy Sciences under Contract
No. DE-SC0020330. K.S. was supported by the Program for Leading Graduate
Schools (MERIT-WINGS).

## References

  * Nicolis and Prigogine [1977] G. Nicolis and I. Prigogine,  _Self Organization in Non Equilibrium Systems_ (Wiley, 1977).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * Koch and Meinhardt [1994] A. J. Koch and H. Meinhardt, Biological pattern formation: from basic mechanisms to complex structures, Rev. Mod. Phys. 66, 1481 (1994).
  * Pismen [2006] L. Pismen,  _Patterns and Interfaces in Dissipative Dynamics_ (Springer, 2006).
  * Cross and Greenside [2009] M. Cross and H. Greenside,  _Pattern Formation and Dynamics in Non-Equilibrium Systems_ (Cambridge University Press, 2009).
  * Le Berre  _et al._ [2002] M. Le Berre, E. Ressayre, A. Tallet, Y. Pomeau, and L. Di Menza, Example of a chaotic crystal: The labyrinth, Phys. Rev. E 66, 026203 (2002).
  * Echeverría-Alar and Clerc [2020] S. Echeverría-Alar and M. G. Clerc, Labyrinthine patterns transitions, Phys. Rev. Res. 2, 042036 (2020).
  * Elder  _et al._ [1992] K. R. Elder, J. Viñals, and M. Grant, Dynamic scaling and quasiordered states in the two-dimensional swift-hohenberg equation, Phys. Rev. A 46, 7618 (1992).
  * Ouyang and Swinney [1991] Q. Ouyang and H. L. Swinney, Transition from a uniform state to hexagonal and striped turing patterns, Nature 352, 610 (1991).
  * Cross and Meiron [1995] M. C. Cross and D. I. Meiron, Domain Coarsening in Systems Far from Equilibrium, Phys. Rev. Lett. 75, 2152 (1995).
  * Christensen and Bray [1998] J. J. Christensen and A. J. Bray, Pattern dynamics of Rayleigh-Bénard convective rolls and weakly segregated diblock copolymers, Phys. Rev. E 58, 5364 (1998).
  * Gunaratne  _et al._ [1995] G. H. Gunaratne, R. E. Jones, Q. Ouyang, and H. L. Swinney, An Invariant Measure of Disorder in Patterns, Phys. Rev. Lett. 75, 3281 (1995).
  * Hu  _et al._ [2004] S. Hu, D. I. Goldman, D. J. Kouri, D. K. Hoffman, H. L. Swinney, and G. H. Gunaratne, Stages of relaxation of patterns and the role of stochasticity in the final stage, Nonlinearity 17, 1535 (2004).
  * Hu  _et al._ [2005] S. Hu, G. Nathan, D. J. Kouri, D. K. Hoffman, and G. H. Gunaratne, Statistical characterizations of spatiotemporal patterns generated in the Swift–Hohenberg model, Chaos 15, 043701 (2005).
  * Hou  _et al._ [1997] Q. Hou, S. Sasa, and N. Goldenfeld, Dynamical scaling behavior of the swift-hohenberg equation following a quench to the modulated state, Phys. A: Stat. Mech. Appl. 239, 219 (1997).
  * Seul and Wolfe [1992a] M. Seul and R. Wolfe, Evolution of disorder in two-dimensional stripe patterns: “smectic” instabilities and disclination unbinding, Phys. Rev. Lett. 68, 2460 (1992a).
  * Seul and Wolfe [1992b] M. Seul and R. Wolfe, Evolution of disorder in magnetic stripe domains. ii. hairpins and labyrinth patterns versus branches and comb patterns formed by growing minority component, Phys. Rev. A 46, 7534 (1992b).
  * Seul  _et al._ [1992] M. Seul, L. R. Monar, and L. O’Gorman, Pattern analysis of magnetic stripe domains morphology and topological defects in the disordered state, Philosophical Magazine B 66, 471 (1992).
  * Hansen and Krumme [1984] P. Hansen and J.-P. Krumme, Magnetic and magneto-optical properties of garnet films, Thin Solid Films 114, 69 (1984), special Issue on Magnetic Garnet Films.
  * Fratello  _et al._ [1996] V. Fratello, S. Licht, and C. Brandle, Innovative improvements in bismuth-doped rare-earth iron garnet Faraday rotators, IEEE Trans. Mag. 32, 4102 (1996).
  * Fratello  _et al._ [1986] V. J. Fratello, S. E. G. Slusky, C. D. Brandle, and M. P. Norelli, Growth-induced anisotropy in bismuth: Rare-earth iron garnets, J. Appl. Phys. 60, 2488 (1986).
  * Fakhrul  _et al._ [2019] T. Fakhrul, S. Tazlaru, L. Beran, Y. Zhang, M. Veis, and C. A. Ross, Magneto-Optical Bi:YIG Films with High Figure of Merit for Nonreciprocal Photonics, Adv. Opt. Mater. 7, 1900056 (2019).
  * Puchalska  _et al._ [1978] I. B. Puchalska, G. A. Jones, and H. Jouve, A new aspect on the observation of domain structure in garnet epilayers, J. Phys. D: Appl. Phys. 11, L175 (1978).
  * Ren  _et al._ [2016] S. Ren, K. He, R. Girshick, and J. Sun, Faster r-cnn: Towards real-time object detection with region proposal networks (2016), arXiv:1506.01497 [cs.CV] .
  * Redmon  _et al._ [2016] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, You only look once: Unified, real-time object detection (2016), arXiv:1506.02640 [cs.CV] .
  * LeCun and Bengio [1998] Y. LeCun and Y. Bengio, Convolutional Networks for Images, Speech, and Time Series, in  _The Handbook of Brain Theory and Neural Networks_ (MIT Press, Cambridge, MA, USA, 1998) p. 255–258.
  * Kim and de Araújo [2007] H. Y. Kim and S. A. de Araújo, Grayscale Template-Matching Invariant to Rotation, Scale, Translation, Brightness and Contrast, in  _Advances in Image and Video Technology: Second Pacific Rim Symposium, PSIVT 2007 Santiago, Chile, December 17-19, 2007 Proceedings_ (Springer-Verlag, Berlin, Heidelberg, 2007) p. 100–113.
  * Kim  _et al._ [2013] H. Y. Kim, R. H. Maruta, D. R. Huanca, and W. J. Salcedo, Correlation-based multi-shape granulometry with application in porous silicon nanomaterial characterization, J. Porous Mater. 20, 375 (2013).
  * Note [1] The characteristic length scale determined by the stripe period decreases with increasing steps and hence, the number of defects in the same area increases in general. To quantify the number of defects in this length scale, we rescale $\rho_{\rm t}$ by the factor of $\left(\frac{2\pi}{q_{0}}\right)^{2}$.