# Interacting Dark Side of Universe Through Generalized Uncertainty Principle

M<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of
Physics, Faculty of Science, University of Zabol, Zabol, Iran.  M. Fathi
Department of Physics, Shahid Beheshti University, G. C., Evin, Tehran, 19839,
Iran.  B. Mostaghel  Department of Physics, Shahid Beheshti University, G. C.,
Evin, Tehran, 19839, Iran.  S<EMAIL_ADDRESS>Federal University of Latin-American Integration, Technological Park of Itaipu
PO box 2123, Foz do Iguaçu-PR, 85867-670, Brazil,
Departamento de Física, Universidade Federal de Pernambuco, Recife, PE
52171-900, Brazil.

(Day Month Year; Day Month Year)

###### Abstract

We investigate the impact of the generalized uncertainty principle proposed by
some approaches to quantum gravity such as string theory and doubly special
relativity on the cosmology. Using generalized Poisson brackets, we obtain the
modified Friedmann and Raychaudhuri equations and suggest a dynamical dark
energy to explain the late time acceleration of the Universe. After
considering the interaction between dark matter and dark energy, originated
from the minimal length, we obtain the effective cosmological parameters and
equation of state parameter for dark matter and dark energy. Finally, we show
that the resulting model is equivalent to the Phantom and Tachyon fields.

###### keywords:

Generalized Uncertainty Principle; Dark Matter; Dark Energy

PACS numbers: 04.50.Kd, 95.30.Sf, 95.35.+d, 98.80.-k

##  1 Introduction

The existence of minimal length [1] has been predicted by various theoretical
models such as string theory [2] and Doubly Special Relativity (DSR) [3]. The
presence of this minimal length may arise from imposing the quantum gravity
effects in quantum mechanics through modification of the Heisenberg
uncertainty principle, which is well-known as the generalized uncertainty
principle (GUP) [4, 5]. From the perturbative string theory, this length is
due to the fact that the strings cannot influence through distances smaller
than their size. An interesting property of the existence of the minimal
length is the modification of the standard commutation relation between
canonical position and momentum in quantum mechanics and also Poisson brackets
algebra in classical mechanics.

Various authors have been used the deformed Poisson brackets (DPB) in
classical cosmological models [6]. Also, the quantum commutator of GUP is
applied for a variety of quantum cosmological models [7]. We should note that
GUP deformation of Poisson brackets, for a point particles, implies a clear
violation of the equivalence principle (EP) [8]. Consistency between DPB and
EP can be partially recovered only for composite systems [9] and at the price
of defining different $\beta$ parameters for different species of particles.
That is expected that the EP to be violated to some degree by any quantum
theory of gravity and by many alternative classical theories as well. The
violation of this principle in scalar-tensor theories of classical gravity is
linked to the existence of scalar gravitational fields in addition to the
usual tensor field. On the other hand, one can expect the violation of EP in
classical cosmological models at the presence of the deformed PBs. In fact,
the origin of this violation is in the usage of the minimal length of quantum
gravity in the context of the classical gravitational theory.

Astrophysical observations show that the Universe is in an accelerating
expansion phase [10] which is an evidence to the existence of dark energy.
There exist several proposals which explain the origin of dark energy such as
cosmological constant, the expectation value of the vacuum ground state, which
leads to the hierarchy problem [11]. Modified theories of gravity such as
$f(R,T)$ [12], DGP braneworld [13], the extrinsic cosmology [14] where dark
energy could actually be the manifestation of the local extrinsic shape of the
spacetime, scalar-tensor theories [15, 16], modified matter models such as
quintessence [17, 18], K-essence, coupled dark energy [19] and dynamical dark
energy models [20] are some examples. The simplest candidates to explain the
dynamical dark energy are the minimal coupling scalar fields like as the
quintessence and Phantom models [21].

In this article, using the classical GUP, we modify the Friedmann and
Raychaudhuri equations to propose a new way to explain the interaction of the
dark side of Universe. By using the reconstruction approach to scalar fields
[20], we obtain a Phantom potential and Tachyon field in the model and show
that this potential originated from the GUP effects. In Section II, we present
the model that assists in our investigation. In section III we consider an
interaction between dark energy and dark matter trough GUP, and finally in the
last two sections we reconstruct Phantom and Tachyon fields which are
consistent to the model.

##  2 Classical cosmology at the presence of GUP

Let us consider the isotropic and homogeneous Friedmann-Lemaître-Robertson-
Walker (FLRW) Universe with line element

$\displaystyle
ds^{2}=-N^{2}(t)d{t}^{2}+{a^{2}(t)}\left(\frac{dr^{2}}{{1-kr^{2}}}+r^{2}d\Omega^{2}\right),$
(1)

where $N(t)$ denotes the lapse function, $a(t)$ is the scale factor of the
Universe and $k=0,+1,-1$ represents the usual spatial curvature. The action
functional corresponding to the line element (1) displays in the gravitational
and matter sectors (with the latter as perfect fluid) [22]

$\displaystyle\begin{split}S&=\frac{1}{16\pi
G}\int{\sqrt{-g}{R}d^{4}x}-\int{\sqrt{-g}\rho d^{4}x}\\\
&=\frac{3\pi\sigma}{4G}\int(-\frac{a\dot{a}^{2}}{N}+{k}Na)dt-2\pi^{2}\sigma\int{Na^{3}}\rho
dt,\end{split}$ (2)

where $\rho=\sum\limits_{i}\rho_{i}$ is the total matter density of the
Universe, $2\pi^{2}\sigma$ is the spatial volume of the metric which remains
finite for the observable Universe with the radius about $H_{0}^{-1}$
$(\sigma=1\,\,\text{for}\,\,k=1)$ and overdot denotes differentiation respect
to cosmic time $t$. According to the energy conservation the matter density
content of the Universe takes the form
$\rho_{i}=\rho_{0i}(a/a_{0})^{-3(w_{i}+1)}$ where $w_{i}$ denotes the equation
of state (EoS), $w_{i}={p_{i}}/{\rho_{i}}$, for $i$-th component of fluid and
$\rho_{i0}$ is the present density. The action functional of the model will be

$S=\frac{3\pi\sigma}{4G}\int{\left[\left(\frac{-\dot{a}^{2}a}{N}+kNa\right)-Na^{3}\sum\limits_{i}\rho_{0i}(\frac{a}{a_{0}})^{-3(\omega_{i}+1)}\right]}dt\,.$
(3)

One can rewrite the present density of $i$-th component of density,
$\rho_{0i}$, in terms of density parameter, $\Omega_{i}$, as
$\rho_{0i}=\dfrac{3{H}^{2}_{0}\Omega_{i}}{8\pi G}$ . If we use the
dimensionless time defined by $d\eta=H_{0}dt$ and ignoring the overall
constant factor, Lagrangian of the model in terms of density parameters will
be

$\displaystyle\mathcal{L}=-\frac{M}{2}\left(\frac{x\dot{x}^{2}}{N}+N\Omega_{k}x+N\sum\limits_{i}\Omega_{i}x^{-3w_{i}}\right)\,,$
(4)

where $\Omega_{k}=\frac{-{k}}{a_{0}^{2}H_{0}^{2}}$, $x=\frac{a}{a_{0}}$
denotes the rescaled dimensionless scale factor

$\displaystyle M=\frac{3\pi a_{0}^{3}H_{0}}{2G},$ (5)

and after this overdot denotes differentiation respect to $\eta$. The
conjugate momentum to the scale factor $x$ and the primary constraint are
given by

$p=\frac{\partial\mathcal{L}}{\partial\dot{x}}=-\frac{Mx}{N}\dot{x},\hskip
11.38092ptp_{N}=\frac{\partial\mathcal{L}}{\partial\dot{N}}=0.$ (6)

Consequently, the Hamiltonian corresponding to the Lagrangian (4) will be

${\mathcal{H}}=N\left[-\frac{p^{2}}{2Mx}+\frac{M}{2}\Omega_{k}x+\frac{M}{2}\sum\limits_{i}\Omega_{i}x^{-3w_{i}}\right].$
(7)

In Eq. (7), $N$ is a Lagrange multiplier, hence it compel the Hamiltonian
constraint

$\displaystyle-\frac{p^{2}}{2Mx}+\frac{M}{2}\Omega_{k}x+\frac{M}{2}{\sum\limits_{i}}\Omega_{i}x^{-3w_{i}}=0.$
(8)

###  2.1 The Generalized Uncertainty Principle

The generalized uncertainty principle between canonical conjugates
$(X,\Pi_{X})$ which leads to a minimum length in one dimension is given by
[23]

$\displaystyle\triangle
X\triangle\Pi_{X}\geqslant\frac{1}{2}\left(1+\beta(\triangle\Pi_{X})^{2}+\gamma\right),$
(9)

where $\beta$ and $\gamma$ are two parameters, independent of $\Delta{X}$ and
$\Delta{\Pi_{X}}$, but may, in general, depend on the expectation value of the
position and the momentum operators. The above generalization of uncertainty
principle corresponds to the following commutation relation

$\displaystyle[X,\Pi_{X}]=i(1+\beta\Pi_{X}^{2}),$ (10)

where by comparing Eqs. (9) and (10), we obtain
$\gamma=\beta\langle\Pi_{X}\rangle^{2}$. A more general case of such
commutation relations are studied in Ref. [Hinvichsen1996]. The various
applications of the low energy effects of the modified Heisenberg uncertainty
principle relations have been extensively studied in Refs.
[Brau,Akhoury2003,Hofmann2003,Nozari2006]. The above equations show that the
smallest uncertainty in position has the value

$\displaystyle\Delta X_{min}=\sqrt{\beta}\,,$ (11)

wherein generally, the parameter of $\beta$ can’t be fixed by the theory.

One major feature of Eq. (11) is that the physics below $\sqrt{\beta}$ becomes
inaccessible and therefore it defines a natural cut-off which prevents from
the usual UV divergences. On the other hand, the second consequence of GUP is
the appearance of an intriguing UV/IR mixing, first noticed in the ADS/CFT
correspondence [29]. As we know, the UV/IR mixing means that one can probe
short distances physics (high energy physics) by long distances physics (low
energy physics), and therefore, it justifies the use of classical mechanics in
the presence of a minimal length. We point that the UV/IR mixing is also a
feature of noncommutative quantum field theory [30]. On the other hand, some
scenarios have been proposed where non zero minimal length is related to large
extra dimensions [31], to the running coupling constant [32] and to the
physics of black holes production [33].

Let us now investigate the effects of deformed Poisson algebra in the presence
of a minimal length. Since there is a UV/IR mixing embodied in the deformed
commutation relation [6], it is also important to study the effects of the
minimal length in a classical context. Note that to obtain the classical and
commutative limit, it is necessary to consider simultaneously
$\hbar\rightarrow 0$ and $\beta\rightarrow 0$. Specifically, if $\beta\neq 0$
there is no classical limit. Therefore, to embedding the effect of UV-IR
mixing of deformation in the classical cosmology at the presence of minimal
length, it is enough to consider the limit of $\hbar\rightarrow 0$. Then, the
modified Poisson brackets corresponding to the GUP will be [34]

$\displaystyle\\{X,X\\}=\\{\Pi_{X},\Pi_{X}\\}=0,$ (12)
$\displaystyle\\{X,\Pi_{X}\\}=1+\beta\Pi_{X}^{2}.$ (13)

Such deformed Poisson algebra is used in Ref. [Benczik] to investigate the
effects of the deformation on the classical orbits of particles in a central
force field and on the Kepler third law. Also, the cosmological constant
problem and removability of initial curvature singularity are investigated in
Ref. [Bina]. The effect of modified Poisson algebra defined in Eq. (12) is
also investigated in a dilatonic cosmological model at Ref. [Vakili1] and in
Bianchi I and II cosmological models in Ref. [Farag]. It is examined in de
Sitter and Anti-de Sitter classical cosmological models in Ref. [Vakili2]. The
effect of above algebra on inflation parameters is investigated in Ref.
[Nozari]. In addition, the authors of Ref. [Zeynali] are studied the effects
of the deformed Poisson algebra Eq. (12) in multidimensional cosmology. And as
the last example, we should mention, the signature changing classical
cosmology at the presence of minimal length is described in Ref. [Ghaneh]. The
result of the above investigations is that the modifications of physics at
small scales, via introducing the minimal length in deformed Poisson algebra,
have a rather profound effect on classical physics at large scales, something
similar to the UV/IR mixing.

One can easily show that the deformed Poisson bracket algebra given in Eq.
(12) is satisfied by the classical dynamical position $x$ and conjugate
momentum $p$ given by [42]

$X=x,\hskip 8.5359pt\Pi_{X}=p+\frac{\beta}{3}p^{3},$ (14)

in which $x$ and $p$ obey the usual Poisson bracket $\\{x,p\\}=1$. Note that
in the natural units that we used in this article all quantities are
dimensionless: $[S]=[\mathcal{L}]=[x]=[p]=[M]=[\beta]=1$.

By substituting relations (14) into Eq.(7) the deformed Hamiltonian will be

$\displaystyle{\mathcal{H}}_{def}=N\left(-\frac{1}{2Mx}\left[p^{2}+\frac{2\beta}{3}p^{4}\right]+\frac{M}{2}\Omega_{k}x+\frac{M}{2}\sum\limits_{i}\Omega_{i}x^{-3w_{i}}\right).$
(15)

Consequently, according to Eq.(15) the equations of motion in the gauge of
$N=1$ will be

$\displaystyle\dot{x}=\dfrac{\partial{{{\mathcal{H}}_{def}}}}{\partial{p}}=-\dfrac{p}{Mx}{\left(1+\frac{4\beta
p^{2}}{3}\right)}\,$ (16)
$\displaystyle\dot{p}=-\dfrac{\partial{{\mathcal{H}}_{def}}}{\partial{x}}=-\dfrac{p^{2}}{2Mx^{2}}\left(1+\dfrac{2\beta}{3}p^{2}\right)+\dfrac{M}{2}\Omega_{k}-\dfrac{3M}{2}\sum\limits_{i}\Omega_{i}w_{i}x^{-3w_{i}-1}.$
(17)

The first equation in (16) gives us the relation of momentum with $\dot{x}$.
By inserting the obtained momentum in Hamiltonian constraint (15) the deformed
Friedmann equation will be

$\displaystyle\left(\frac{H}{H_{0}}\right)^{2}=\Omega_{k}x^{-2}+\sum\limits_{i}\Omega_{i}x^{-3(w_{i}+1)}+2\beta
M^{2}x^{4}\left(\Omega_{k}x^{-2}+\sum\limits_{i}\Omega_{i}x^{-3(w_{i}+1)}\right)^{2},$
(18)

where $H=\frac{1}{a}\frac{da}{dt}$ denotes the Hubble parameter in the
comoving coordinates. The second equation of (16) gives the modified
Raychaudhuri equation

$\displaystyle\begin{split}\frac{1}{H_{0}^{2}}\frac{\ddot{x}}{x}=&-\frac{1}{2}\sum\limits_{i}\Omega_{i}(1+3w_{i})x^{-3(w_{i}+1)}+2\beta
M^{2}x^{4}\left(\frac{\Omega_{k}}{x^{2}}+\sum\limits_{i}\Omega_{i}x^{-3(w_{i}+1)}\right)\times\\\
&\left(\frac{\Omega_{k}}{x^{2}}-\sum\limits_{i}\Omega_{i}w_{i}x^{-3(w_{i}+1)}\right)\,.\end{split}$
(19)

Also, the age of Universe will be

$\displaystyle
t_{0}=\frac{1}{H_{0}}\int_{0}^{1}{dx[\Omega_{k}+\sum\limits_{i}\Omega_{i}x^{-3w_{i}-1}+2\beta
M^{2}x^{2}(\Omega_{k}+\sum\limits_{i}\Omega_{i}x^{-3w_{i}-1})^{2}]^{-\frac{1}{2}}}.$
(20)

Let us consider a spatially flat Universe dominated by vacuum energy with EoS
$w=-1$. According to Eq.(18) the modified Friedmann equation becomes

$\displaystyle\begin{array}[]{cc}{H}^{2}={H}_{0}^{2}\left(\Omega_{\Lambda}+2\beta
M^{2}\Omega^{2}_{\Lambda}x^{4}\right)\equiv{H}_{0}^{2}\Omega_{{DE}}{(x)},\end{array}$
(22)

where we defined the modified dark energy parameter as

$\displaystyle\Omega_{{DE}}(x)=\Omega_{\Lambda}+2\beta
M^{2}\Omega^{2}_{\Lambda}x^{4},$ (23)

and $\Omega_{\Lambda}$ denotes the density parameter of cosmological constant
at the present epoch. From Eq.(23) we define

$\displaystyle\Omega_{\beta}^{0}\equiv 2\beta
M^{2}\Omega^{2}_{\Lambda}=\Omega^{0}_{{DE}}-\Omega_{\Lambda}\,,$ (24)

which states the dependence of GUP parameter $\beta$ to the density parameters
of dark energy and the cosmological constant. As we see in our model the GUP
parameter could take the positive and negative value according to the
difference between $\Omega_{DE}^{0}$ and $\Omega_{\Lambda}$. With substituting
value of $M$ from Eq. (5) in Eq.(24) we obtain the order of magnitude of
$\beta$ as

$\displaystyle\beta=\frac{\Omega_{\beta}}{2M^{2}{\Omega_{\Lambda}}^{2}}\approx\frac{G^{2}}{a_{0}^{6}H_{0}^{2}}\approx\frac{G^{2}}{d_{(H_{0})}^{6}H_{0}^{2}}\approx
10^{-84}\,,$ (25)

where $d_{H_{0}}$ denotes the Hubble distance at the present epoch and
$\Omega_{\beta}$ is in order of $10^{-2}$. The order of the GUP parameter is
reasonable in our model because we apply the GUP approach to the evolution of
the Universe in the present epoch where the quantum effects are small. Note
that this limit is 49 orders of magnitude below the Planck scale.

Benczik and et.al in Ref. [Benczik] obtain the bound $\beta\lesssim 10^{-66}$
from investigation the precession perihelion of Mercury by GUP-deformed
Newtonian mechanics. While Scardigli and et.al in Ref. [Scardigli2] compute
upper bound for the GUP parameter in light deflection by deformation
Schwarzschild metric as $\beta\lesssim 10^{78}$, and for perihelion precession
by that deformed metric obtained the lower bound for GUP parameter as
$\beta\lesssim 10^{69}$. As well as, in Ref. [Farag Ali] the author
investigated derivation of the modified Hawking temperature from the GUP, and
find an upper bound on GUP parameter for the deflection of light
$\alpha_{0}<10^{41}$, which there is $\alpha_{0}$ in contrast to Ref.
[Scardigli2] namely $\alpha_{0}^{2}\sim\beta$, and upper bound for of the
perihelion precession of Mercury is $\alpha_{0}<10^{35}$. Recently Scardigli
et.al in Ref. [Scardigli3] compute the value of GUP parameter by comparing two
different low energy corrections. They, first involved $\beta$ through the GUP
procedure and then insert the corrections of Newtonian potential which is
leads to a quantum correction for Schwarzschild metric that it computed by
Donoghue[50]. From this approach, they obtain that the order of GUP parameter
is $\beta=\frac{82\pi}{5}$.

The dynamic of dark energy in our model originates from the parameter $\beta$.

We can write the dark energy density parameter in terms of red-shift as

$\displaystyle\Omega_{{DE}}(z)=\Omega_{\Lambda}+\Omega_{\beta}^{0}(1+z)^{-4}\,,$
(26)

We plot the evolution of $\Omega_{{DE}}(z)$ in the Fig. 1. To make more sense,
we compare the behavior of GUP model with standard $\Lambda$CDM model.
According to this figure, there is a meaningful difference between standard
$\Lambda$CDM model and GUP dark energy. The modified EoS for dark energy
becomes

$\displaystyle
w_{{DE}}=-1+\frac{1}{3}\frac{d\ln\Omega_{{{DE}}}(z)}{d\ln(1+z)}=-1-\frac{4}{3}\frac{\Omega^{0}_{\beta}(1+z)^{-4}}{\Omega^{0}_{\Lambda}+\Omega^{0}_{\beta}(1+z)^{-4}}.$
(27)

This relation shows that in the early Universe, i.e., $z=\infty$, we have
$\Omega_{{DE}}\approx\Omega^{0}_{\Lambda}$ and the effect of the second term
is important in the very late times. In Fig. 1, the equation of state of the
proposed model is plotted for two values of the $\beta$ parameter.

Figure 1:  panel a: $\Omega_{{DE}}(z)$ as a function of redshift. panel b:The
evolution of EoS of dark energy respect to red-shift. We change the value of
$\beta$ while other parameters have been fixed to $\Omega^{0}_{DE}=0.70$ and
$\Omega^{0}_{DM}=0.30$.

According to Eq. (27), at the present epoch, the modified EoS parameter of the
dark energy is

$\displaystyle
w^{0}_{{DE}}=-1-\frac{4}{3}\frac{\Omega^{0}_{{\beta}}}{\Omega^{0}_{{DE}}}\,.$
(28)

For the positive value of the parameter $\beta$, we get
${{\Omega^{0}_{\beta}\geq 0}}$ and consequently, ${w^{0}_{DE}<-1}$. It implies
that in our approach, Universe enters the Phantom phase at the late time.

##  3 Interacting Dark Energy and Dark Matter

Let us now consider the interaction between dark energy and the dark matter
content of the Universe (where we neglect the contribution of baryonic matter
and radiation) [43]. As we know the continuity equation for dark energy and
dark matter in standard cosmology are

$\displaystyle\dot{\rho}_{\Lambda}=0,$ (29)
$\displaystyle\dot{\rho}_{DM}+3H\rho_{DM}=0.$ (30)

Substituting the effective dark energy parameter from Eq. (23) into Eq. (29)
we obtain the modified continuity equations which can be interpreted as the
interaction between dark energy and dark matter. The simplest form of the
interaction equations are

$\displaystyle\dot{\rho}_{DE}=-3H(1+w_{DE})\rho_{DE}=Q,$ (31)
$\displaystyle\dot{\rho}_{DM}+3H\rho_{DM}=-3Hw_{DM}=-Q,$ (32)

where $w_{DE}$ and $w_{DM}$ denote the effective EoS parameters for dark
energy and dark matter respectively and $Q$ represents the interacting term.
Generally speaking, the function $Q$ represents the energy transition between
dark energy and dark matter. Using Eq. (31), we obtain the following two
equations for EoS and energy density parameter of dark matter content

$\displaystyle w_{DM}=-(1+w_{DE})\frac{\Omega_{DE}}{\Omega_{DM}},$ (33)
$\displaystyle\frac{d\Omega_{DM}}{dz}-\frac{3}{1+z}\Omega_{DM}=-\frac{d\Omega_{DE}}{dz}.$
(34)

Substituting $\Omega_{DE}$ from Eq.(26) into the second Equation of (34) gives

$\Omega_{DM}=\Omega_{DM}^{0}(1+z)^{3}+\frac{4}{7}\Omega_{\beta}^{0}\left((1+z)^{3}-(1+z)^{-4}\right),$
(35)

The evolution of interacting dark matter respect to the red-shift, shown in
the Fig. (2). The effective EoS parameter for dark matter takes the form

$\displaystyle\begin{split}w_{DM}&=-1+\frac{1}{3}\frac{d\ln(\Omega_{DM})}{d\ln(1+z)}\\\
&=\frac{28\Omega_{\beta}^{0}(1+z)^{-4}}{21\Omega_{DM}^{0}(1+z)^{3}+12\Omega_{\beta}^{0}((1+z)^{3}-(1+z)^{-4})}.\end{split}$
(36)

The evolution of $w_{DM}$ respect to red-shift is shown in Fig.(2). For the
positive value of the parameter $\Omega^{0}_{\beta}$, dark matter behaves like
as a stiff fluid with $w^{0}_{DM}\approx 1$. When the GUP parameter, $\beta$
(and consequently $\Omega^{0}_{\beta}$), tends to zero, the EoS parameter
$w_{DM}$ and density parameter of dark matter, $\Omega_{DM}$, take the values
zero and $\Omega_{DM}^{0}$ respectively which shows that the deviation of EoS
for ordinary dark matter is a consequence of the GUP and the special form of
interaction in the Eq. (34) for the proposed dark energy model. To make more
sense and in order to compare the role of the parameter $\Omega^{0}_{\beta}$
on the dynamics of Universe respect to the standard model, in the Fig. 3 we
plot the ratio of dark energy density parameter to the dark matter density
parameter, $\Omega_{{DM}}(z)/\Omega_{{DE}}(z)$. For the positive value of the
parameter $\Omega^{0}_{\beta}$, the contribution of the dark matter in the
Universe is lower than the standard $\Lambda$CDM model.

Figure 2: panel a: The evolution of the dark matter density parameter with
respect to red-shift. panel b: Equation of state of dark matter respect to the
red-shift. We change the value of $\Omega^{0}_{\beta}$ while other parameters
have been fixed to $\Omega^{0}_{DE}=0.70$ and $\Omega^{0}_{DM}=0.30$.  Figure
3: The ratio of dark energy density parameter over dark matter density
parameter as a function of red-shift. This figure shows that for
$\Omega^{0}_{\beta}>0$ the contribution of dark matter in the model is lower
than the standard $\Lambda$CDM model, while for $\Omega^{0}_{\beta}<0$ the
contribution of dark matter in the model is more than a standard $\Lambda$CDM
model. We change the value of $\Omega^{0}_{\beta}$ while other parameters have
been fixed to $\Omega^{0}_{DE}=0.70$ and $\Omega^{0}_{DM}=0.30$.

Deceleration parameter defined by

$\displaystyle
q(z)=\frac{(1+3w_{DM})\rho_{DM}+(1+3w_{DE})\rho_{DE}}{2(\rho_{DM}+\rho_{DE})}.$
(37)

In Fig. 4 we plot the deceleration parameter $q(z)$ of the proposed model
respect to the red-shift. According to this figure, Universe enters to the
late time acceleration expansion with a time delay, respect to the concordance
$\Lambda$CDM model. As we see in Fig. 4, for high red-shifts, the deceleration
parameter tends to $1$ and the dynamic of Universe determined by dark matter
component in our model.

Figure 4: The behavior of deceleration parameter with respect to the red-
shift. We change the value of $\Omega^{0}_{\beta}$ while other parameters have
been fixed to $\Omega^{0}_{DE}=0.70$ and $\Omega^{0}_{DM}=0.30$.

##  4 Reconstruction of Scalar Field

###  4.1 Reconstruction of Phantom Fields

Now, we try to reconstruct a scalar field which coincides with results of the
previous section. We use the reconstruction method which proposed in Refs.
Tsujikawa2005ju,Guo2005ata. Since for different values of the parameter
$\Omega^{0}_{\beta}$, there is no phantom crossing, we can consider scalar
fields with minimal coupling [46]. In minimal coupling with Einstein gravity,
energy density and pressure for scalar fields are

$\displaystyle\begin{array}[]{cc}\rho_{\phi}=\frac{1}{2}\epsilon\dot{\phi}^{2}+V(\phi),\\\
P_{\phi}=\frac{1}{2}\epsilon\dot{\phi}^{2}-V(\phi),\end{array}$ (40)

where $\epsilon=\pm 1$ for quintessence and phantom field respectively and
hereafter we will set $\epsilon=-1$. Using the EoS parameter,
$w_{\phi}={P_{\phi}}/{\rho_{\phi}}$, the kinetic term becomes

$\dot{\phi}^{2}=\epsilon\left(1+w_{\phi}\right)\rho_{\phi}.$ (41)

In the proposed dark energy model, energy density always is positive and the
EoS is smaller than $-1$, respect to the value of the parameter
$\Omega^{0}_{\beta}$. and, $\dot{\phi}$ is a real filed. If, we want to apply
such equation for our model, since $\omega_{{DE}}<-1$, it forces that
$\dot{\psi}^{2}<0$, which means that the field must be imaginary. To remove
this difficulty, we use the replacement
$\dot{\psi}^{2}\longrightarrow-\dot{\Psi}^{2}$, therefore, the Tachyonic
Lagrangian changes to

$\displaystyle\mathcal{L}=-V(\Psi)\sqrt{1+\dot{\Psi}^{2}}.$ (42)

On the other hand the potential term will be

$\displaystyle V(\phi)=\frac{1}{2}\left(1-w_{\phi}\right)\rho_{\phi}.$ (43)

To construct the scalar field, we assume that $\rho_{\phi}=\rho_{{DE}}$.
According to the Eqs. (26) and (35) we obtain the kinetic term as

$\displaystyle\phi^{\prime 2}=\frac{3}{8\pi
G_{N}}\frac{\epsilon(1+w_{\phi})\Omega_{\phi}}{(1+z)^{2}\left(\Omega_{{DM}}(z)+\Omega_{\phi}(z)\right)},$
(44)

which over-prime denotes differentiation respect to red-shift. Substituting
${\tilde{\phi}=\sqrt{\frac{8\pi G_{N}}{3}}\phi}$, in Eq. (44) changes to the
form

$\displaystyle\tilde{\phi}^{\prime
2}=\frac{\epsilon(1+w_{\phi})\Omega_{\phi}}{(1+z)^{2}\left(\Omega_{{DM}}(z)+\Omega_{\phi}(z)\right)}.$
(45)

Therefore, we have

$\displaystyle\frac{d\tilde{\phi}}{dz}=\pm\sqrt{\frac{\epsilon(1+w_{\phi})\Omega_{\phi}}{(1+z)^{2}\left(\Omega_{{DM}}(z)+\Omega_{\phi}(z)\right)}}.$
(46)

The sign of the above equation is arbitrary, from now on we choose the
positive sign. By numerical analysis of Eq. (46), we show that the scalar
field dynamics become important at the late times, as shown in Fig. 5. On the
other hand, using Eq. (43) we define a dimensionless scalar potential
according to

$\displaystyle U(\tilde{\phi})=\frac{8\pi G_{N}}{3H^{2}_{0}}V(\tilde{\phi}),$
(47)

which evolves with respect to the scalar field as shown in Fig.(6).

Figure 5: Evolution of scalar field with respect to the red-shift. We change
the value of $\Omega^{0}_{\beta}$ while other parameters have been fixed to
$\Omega^{0}_{DE}=0.70$ and $\Omega^{0}_{DM}=0.30$ and we set
$\tilde{\phi}(z=0)=0$.

Figure 6: panel a: The behavior of the dimensionless potential
$U(\tilde{\phi})$ respect to the red-shift. panel b: The behavior of the
dimensionless potential $U(\tilde{\phi})$ respect to the field $\tilde{\phi}$.
We change the value of $\Omega^{0}_{\beta}$ while other parameters have been
fixed to $\Omega^{0}_{DE}=0.70$ and $\Omega^{0}_{DM}=0.30$.

From the above discussion, we conclude that the kinetic and potential terms at
early times were constant which represent a cosmological constant, with the
passage of time kinetic and potential terms start to evaluate with respect to
time. So the Universe experience a smooth transition from dark matter
dominated epoch to the late time dark energy dominated epoch.

###  4.2 Reconstruction of Tachyon Field

Now, we use a similar approach in the previous section to construct another
kind of scalar filed, namely Tachyon fields. Tachyon fields appear in string
theory, using effective field theory [47]. The interaction between Tachyon
field and dark matter were studied in Ref. Sheykhi2011fe. If the field in this
theory was an ordinary field, it is shown that such filed can be considered as
a candidate for dark energy. It is shown that in the standard Tachyonic field,
the equation of state changes between $-1$ and $0$. Tachyon field Lagrangian
density is

$\displaystyle\mathcal{L}=-V(\psi)\sqrt{1-g^{\mu\nu}\partial_{\mu}\psi\partial_{\nu}\psi}.$
(48)

It is important to notice that mass dimension of the scalar field is $M^{-1}$.
If we consider that Tachyon field depends only on cosmic time, $t$, then the
Lagrangian density becomes

$\displaystyle\mathcal{L}=-V(\psi)\sqrt{1-\dot{\psi}^{2}}.$ (49)

Figure 7: panel a:The evolution of dimensionless Tachyon field,
$\tilde{\psi}$, as a function of red-shift. We assume that
$\tilde{\psi}(z=0)=0$. panel b: The evolution of the dimensionless potential
$V(\tilde{\psi})$ as a function of red-shift. panel c: The behavior of
dimensionless potential $V(\tilde{\psi})$ as a function of Tachyon filed
$\tilde{\psi}$.

Hence, the energy density and the pressure associated to the Tachyon field
will be

$\displaystyle\rho_{\psi}$
$\displaystyle=\frac{V(\psi)}{\sqrt{1-\dot{\psi}^{2}}},$ (50) $\displaystyle
P_{\psi}$ $\displaystyle=-V(\psi)\sqrt{1-\dot{\psi}^{2}}.$ (51)

Using above equations, the EoS parameter is

$\displaystyle w_{\psi}=\frac{P_{\psi}}{\rho_{\psi}}=-1+\dot{\psi}^{2}.$ (52)

If, we want to apply such equation for our model, since $w_{{DE}}<-1$, it
forces that $\dot{\psi}^{2}<0$, which means that the field must be imaginary.
To remove this difficulty, we use the replacement
$\dot{\psi}^{2}\longrightarrow-\dot{\Psi}^{2}$, therefore, the Tachyonic
Lagrangian changes to

$\displaystyle\mathcal{L}=-V(\Psi)\sqrt{1+\dot{\Psi}^{2}}.$ (53)

By considering Lagrangian density (53) we can rewrite energy density and
pressure as

$\displaystyle\rho_{\Psi}$
$\displaystyle=\frac{V(\Psi)}{\sqrt{1+\dot{\Psi}^{2}}},$ (54) $\displaystyle
P_{\Psi}$ $\displaystyle=-V(\Psi)\sqrt{1+\dot{\Psi}^{2}}.$ (55)

So the EoS parameter for Tachyonic field becomes

$\displaystyle w_{\Psi}=\frac{P_{\Psi}}{\rho_{\Psi}}=-1-\dot{\Psi}^{2},$ (56)

which satisfy the condition $w_{\Psi}<-1$. Combining above equations, we
obtain

$\displaystyle V(\Psi)=\pm\sqrt{-w_{\Psi}}\rho_{\Psi}.$ (57)

Similar to the Phantom case, we choose the plus sign.

Now to establish a correspondence between our model and Tachyon field, we use
the previous formula for EoS and energy density. If we define
$\tilde{\psi}=H_{0}\psi$, the kinetic term rewrite as

$\displaystyle\frac{d\tilde{\psi}}{dz}=\frac{1}{1+z}\sqrt{\frac{-(1+w_{\psi})}{\Omega^{0}_{{DM}}(z)+\Omega_{\psi}(z)}}\,.$
(58)

Again by numerical integration, we obtain the behavior of Tachyon field and
potential as shown in Figs. (7) and (7).

##  5 Conclusion

In this paper, we have applied the existence of minimal observable length to
modify the classical field equations in standard cosmology. Considering a
spatially flat Universe filled with the dark energy we obtain the effective
equation of state for the dark side of the Universe containing dark energy and
dark matter. We considered the interaction between dark energy and dark matter
and showed that this interaction is a consequence of the effective EoS
parameters which originated from the GUP parameter. Finally, we constructed
two types of scalar fields for the said dynamical dark energy model, namely
Phantom and Tachyon fields and compared their behavior during the evolution of
the Universe. Strictly speaking, we showed that the what is responsible for a
dynamical dark energy model and interaction between dark components of the
Universe could be the existence of minimal length in Universe.

## References

  * [1] C. A. Mead, Phys. Rev. D 135 (1964) 849.
  * [2] D. J. Gross, P. F. Mende, Phy. Lett. B 197 (1987) 1 ; K. Konishi, G. Paffuti, and P. Provero, Phys. Lett. B 234 (1990) 276.
  * [3] G. Amelino-Camelia, Int. J. Mod. Phys. D 11 (2002) 35, [arXiv:0012051 [gr-qc]]; G. Amelino-Camelia, Phys. Lett. B 510 (2001) 255, [arXiv:0012238 [hep-th]]; C. Rovelli and L. Smolin, Nucl. Phys. B 442 (1995) 593 ; C. Rovelli and L. Smolin, Nucl. Phys. B 456 (1995) 734.
  * [4] D. Amati, M. Ciafaloni and G. Veneziano, Phys. Lett. B 216 (1989) 41; M. Maggiore, Phys. Lett. B 304 (1993) 65 [arXiv:9301067 [hep-th]]; L. J. Garay, Int. J. Mod. Phys. A 10 (1995) 145 [arXiv:9403008 [gr-qc]]; A. Kempf, G. Mangano and R. B. Mann, Phys. Rev. D 52 (1995) 1108 [arXiv:9412167[ hep-th]], C. Kiefer, Quantum Gravity Clarendon, Oxford (2004).
  * [5] S. Hossenfelder, M. Bleicher, S. Hofmann, J. Ruppert, S. Scherer and H. Stoecker, Phys. Lett. B 575 (2003) 85 [arXiv:0305262 [ hep-th]]; C. Bambi and F.R. Urban, Class. Quantum Grav. 25 (2008) 095006 [arXiv: 0709.1965 [ gr-qc]]; A. Kempf, J. Phys. A 30 (1997) 2093 [arXiv:9604045[hep-th]]; F. Brau, J. Phys. A 32 (1999) 7691 [arXiv:9905033 [ quant-ph]]; F. Scardigli, Phys. Lett. B 452 (1999) 39 [arXiv:9904025 [hep-th]]; M. Maggiore, Phys. Lett. B 319 (1993) 83; M. Maggiore, Phys. Rev. D 49 (1994) 5182; A. Kempf, J. Math. Phys. 35 (1994) 4483; R. J. Adler and D. I. Santiago, Mod. Phys. Lett. A 14 (1999) 1371; A. Kempf, Phys. Rev. D 63 (2001) 083514; A. Kempf and J. C. Niemeyer, Phys. Rev. D 64 (2001) 103501; L. N. Chinvichsenhang, D. Minic, N. Okamura and T. Takeuchi, Phys. Rev. D 65, (2002) 125027;
  * [6] S. Benczik,L. N. Chang, D. Minic, N. Okamura, S. Rayyan and T. Takeuchi, Phys. Rev. D 66 (2002) 026003 [arXiv:0204049[hep-th]].
  * [7] M. Kober, Int. J. Mod. Phys. A 27 (2012) 1250106 [arXiv: 1109.4629 [gr-qc]]; P. Pedram, Phys. Rev. D 91 (2015) 063517 [arXiv: 1502.07320 [gr-qc]]; O. Obregon and I. Quiros, Phys. Rev. D 84 (2011) 044005 [arXiv: 1011.3896 [gr-qc]]; S. Jalalzadeh, M. A. Gorji and K. Nozari, Gen. Relativ. Gravit. 46 (2014) 1632 [arXiv: 1310.8065 [gr-qc]]; S. Jalalzadeh, S. M. M. Rasouli and P.V. Moniz, Phys. Rev. D 90 (2014) 023541 [arXiv:1403.1419 [gr-qc]]; S. M. M. Rasouli, A .H. Ziaie, J. Marto and P. V. Moniz, Phys. Rev. D 89 (2014) 044028 [arXiv:1309.6622 [gr-qc]].
  * [8] A.Farag Ali, 2011 Class. Quantum Grav. 28 065013 [arXiv:1101.4181]; F. Scardigli and R. Casadio, Eur.Phys.J.C 75 (2015) [arXiv:1407.0113].
  * [9] V. M. Tkachuk, Phys. Rev. A 86 (2012) 062112.
  * [10] A. Balbi et al., Astrophys. J. 545 (2000) L1; Astrophys. J. 558 (2001) L145 [arXiv:0005124 [astro-ph]]; S. Perlmutter et al., Nature 391 (1998) 51 [arXiv: astro-ph/9712212]; A. G. Riess et al. Astron. J. 116, (1998) 1009 [arXiv:9805201 [astro-ph]].
  * [11] S. Weinberg, Rev. Mod. Phys. 61 (1989) 1.
  * [12] H. Shabani and M. Farhoudi, Phys. Rev. D 90 (2014) 044031 [arXiv: 1407.6187 [gr-qc]].
  * [13] G. R. Dvali, G. Gabadadze and M. Porrati, Phys. Lett. B 485 (2000) 208 [arXiv:0005016 [hep-th]].
  * [14] T. Rostami and S. Jalalzadeh, Phys. Dark Univ. 9-10, (2015) 31 [arXiv:1510.02068 [gr-qc]]; S. Jalalzadeh and T. Rostami, Int. J. Mod. Phys. D 24 (2015) 1550027 [arXiv:1307.1913 [gr-qc]]; M. Heydari-Fard, M. Shirazi, S. Jalalzadeh and H. R. Sepangi, Phys. Lett. B 640, (2006) 1 [arXiv:0607067[gr-qc]].
  * [15] L. Amendola, Phys. Rev. D 60 (1999) 043501 [arXiv:9904120[astro-ph]], J. P. Uzan, Phys. Rev. D 59 (1999) 123510 [arXiv: 9903004[gr-qc]].
  * [16] B. Boisseau, G. Esposito-Farese, D. Polarski, A.A. Starobinsky, Phys.Rev.Lett. 85 (2000) 2236. [arXiv:gr-qc/0001066].
  * [17] Y. Fujii, Phys. Rev. D 26, (1982) 2580.
  * [18] A. A. Starobinsky, JETP Lett. 68, (1998), 757 [astro-ph/9810431]; D. Huterer and M. S. Turner, Phys. Rev. D 90, (1999), 081301 [astro-ph/9808133]; T. Nakamura and T. Chiba, MNRAS 306, (1999) , 696 [astro-ph/9810447].
  * [19] T. Chiba and T. Nakamura, Phys. Rev. D 62, (2000) 121301 [arXiv: 0008175[astro-ph]].
  * [20] U. Alam, V. Sahni and A. A. Starobinsky, JCAP 0406 (2004) 008 [arXiv: 0403687[astro-ph]].
  * [21] R. R. Caldwell, R. Dave and P. J. Steinhardt, Phys. Rev. Lett. 80 (1998) 1582 [arXiv: 9708069[astro-ph]].
  * [22] S. W. Hawking and G. F. R. Ellis,  The Large Scale Structure of Space-Time, Cambridge University Press, Cambridge, (1973).
  * [23] P. Pedram, Int. J. Mod. Phys. D 19 (2010) 2003 [arxiv: 1103.3805 [hep-th]]; P. Pedram, Phys. Lett. B 714 (2012) 317 [arXiv: 1110.2999 [hep-th]]; P. Pedram, Phys. Rev. D 91 (2015) 063517 [arXiv: 1502.07320 [gr-qc]]; B. Vakili and H. R. Sepangi, Phys.Lett.B 651 (2007) 79 [arXiv: 0706.0273 [gr-qc]].
  * [24] H. and A. Kempf, J. Math. Phys. 37, (1996) 2121.
  * [25] F. Brau, J.Phys. A 32 (1999) 7691 [arXiv: 9905033[quant-ph]].
  * [26] R. Akhoury and Y.P. Yao, Phys. Lett. B 572, (2003) 37-42 [arXiv:0302108[ hep-ph]].
  * [27] S. Hossenfelder, M. Bleicher, S. Hofmann, J. Ruppert, S. Scherer and H. St¨ocker, Phys. Lett. B 575 (2003) 85 [arXiv: 0305262[hep-th]].
  * [28] K. Nozari and T. Azizi, Gen. Rel. Grav. 38, (2006) 735.
  * [29] L. Susskind and E. Witten arXiv: 9805114[hep-th]; A. W. Pet and J. Polchinski, Phys. Rev. D 59 (1999) 065011.
  * [30] M. R. Douglas and N. A. Nekrasov, Rev. Mod. Phys. 73, (2001) 977; S. Minwalla, M. Van Raamsdonk and N. Seiberg, JHEP 0002 (2000) 020; R. J. Szabo, Phys. Rept. 378 (2003) 207; A. Micu, M. M. Sheikh-Jabbari, JHEP 0101 (2001) 025.
  * [31] S. Hossenfelder, Mod. Phys. Lett. A 19 (2004) 2727.
  * [32] S. Hossenfelder, Phys. Rev. D 70 (2004) 1054003.
  * [33] S. Hossenfelder, Phys. Lett. B598 (2004) 92.
  * [34] L. Nam Chang, D. Minic, N. Okamura and T. Takeuchi, Phys. Rev. D 65 (2002) 125028 [arXiv:0201017[ hep-th]].
  * [35] A. Bina, K. Atazadeh, S. Jalalzadeh, Int. J. Theor. Phys. 47 (2008) 1354 [arXiv: 0709.3623 [gr-qc]].
  * [36] B. Vakili, Phys. Rev. D 77 (2008) 044023 [arXiv: 0801.2438 [gr-qc]]; H. R. Sepangi, B. Shakerin and B. Vakili, Class. Quant. Grav. 26 (2009) 065003 [arXiv: 0901.3829 [gr-qc]]; B. Majumder, Phys. Rev. D 84 (2011) 064031 [arXiv: 1106.4494 [gr-qc]].
  * [37] A. Farag Ali and B. Majumder, Class. Quant. Grav. 31 (2014) 215007 [arXiv: 1402.5104 [gr-qc]].
  * [38] B. Vakili, Int. J. Mod. Phys. D 18 (2009) 1059 [arXiv: 0811.3481 [gr-qc]].
  * [39] K. Nozari and S. Akhshabi, Int. J. Mod. Phys. D 19 (2010) 513 [arXiv: 0910.2808 [gr-qc]].
  * [40] K. Zeynali, F. Darabi and H. Motavalli, JCAP 12, (2012) 033 [arXiv: 1206.0891 [gr-qc]].
  * [41] T. Ghaneh, F. Darabi and H. Motavalli, Int. J. Mod. Phys. D 22 (2013) 1350026 [arXiv: 1212.1577 [gr-qc]].
  * [42] R. C. S. Bernardo and J. P. H. Esguerra, Few-Body Systems 56 (2015) 219.
  * [43] M. S. Berger and H. Shojaei, Phys. Rev. D 73 (2006) 083528 [arXiv:0601086[gr-qc]]; A. Awad and A. F. Ali, JHEP 1406, (2014) 093 [arXiv: 1404.7825 [gr-qc]]; R. G. Cai, L. M. Cao and N. Ohta, Phys. Rev. D 81, (2010) 061501 [arXiv: 1001.3470 [hep-th]]; A. Sheykhi, Phys. Rev. D 81, (2010) 104011 [arXiv: 1004.0627 [gr-qc]]; S. H. Hendi and A. Sheykhi, Phys. Rev. D 83 (2011) 084012 [arXiv: 1012.0381 [hep-th]].
  * [44] S. Tsujikawa, Phys. Rev. D 72 (2005) 083512 [arXiv: 0508542[astro-ph]].
  * [45] Z.K. Guo, N. Ohta and Y.Z. Zhang, Phys. Rev. D 72, (2005) 023504 [arXiv:0505253 [astro-ph]].
  * [46] Mostaghel, Behrang and Moshafi, Hossein and Movahed, S. M. S., The European Physical Journal C, 77, (2017), [ arXiv:1611.08196v2 [astro-ph.CO]].
  * [47] G.W. Gibbons, Phys. Lett. B 537, (2002) 1 [arXiv:0204008 [hep-th]].
  * [48] F. Scardigli and R. Casadio, Eur. Phys. J. C 75, 425, (2015), [arXiv:1410.8062v2 [hep-th]]; G. Lambiase and F. Scardigli; Phys. Rev. D 97, (2018), 075003, [arXiv:1709.00637v1 [hep-th]].
  * [49] F. Scardigli, G. Lambiase and E. Vagenas, Phys.Lett. B 767, (2017) , 242-246 [arXiv:1611.01469].
  * [50] J. F. Donoghue, Phys. Rev. Lett. 72 ,(1994), 2996; J. F. Donoghue, Phys. Rev. D 50 (1994) 3874.
  * [51] Ahmed Farag Ali, Mohammed M. Khalil, Elias C. Vagenas, EPL (Europhysics Letters), 112, (2015), 2 [arXiv:1510.06365v2 [gr-qc]] .
  * [52] A. Sheykhi, M. Sadegh Movahed and E. Ebrahimi, Astrophys. Space Sci. 339 (2012) 93 [arXiv1105.0325 [hep-th]].


# A SAT+CAS Approach to Finding Good Matrices:
New Examples and Counterexamples

Curtis Bright
University of Waterloo
&Dragomir Ž. Đoković
University of Waterloo
&Ilias Kotsireas
Wilfrid Laurier University
&Vijay Ganesh
University of Waterloo

###### Abstract

We enumerate all circulant good matrices with odd orders divisible by $3$ up
to order $70$. As a consequence of this we find a previously overlooked set of
good matrices of order $27$ and a new set of good matrices of order $57$. We
also find that circulant good matrices do not exist in the orders $51$, $63$,
and $69$, thereby finding three new counterexamples to the conjecture that
such matrices exist in all odd orders. Additionally, we prove a new
relationship between the entries of good matrices and exploit this
relationship in our enumeration algorithm. Our method applies the SAT+CAS
paradigm of combining computer algebra functionality with modern SAT solvers
to efficiently search large spaces which are specified by both algebraic and
logical constraints.

##  1 Introduction

In 2002, circulant “good” matrices were searched for in all odd orders $n$
with $n<40$ (?) and were found to exist in all such orders. This gave evidence
to the conjecture that such matrices actually exist in _all_ odd orders. The
mathematician George Szekeres studied circulant good matrices in terms of an
equivalent type of object which he termed $E$-sequences and in his classic
paper (?) he mentioned this hypothesis as being worthwhile to study using
computers:

> …it is conceivable that $E$-sequences exist for all $n=2m+1$, $m\geq 1$ and
> it is worth testing this hypothesis at least for those orders which are
> accessible to present day computers…

Unfortunately, this conjecture was recently shown to be false in (?) where it
was shown that the orders $n=41$, $47$, and $49$ are counterexamples, i.e., no
circulant good matrices of orders $41$, $47$, and $49$ exist.

In this paper we find the additional three larger counterexamples $n=51$,
$63$, and $69$. In addition, we verify the claim of Đoković and Kotsireas that
there are exactly four inequivalent sets of circulant good matrices of order
$45$ and determine that there is a single set of inequivalent circulant good
matrices of order $57$. We also find a previously undiscovered set of good
matrices of order $27$ demonstrating that Szekeres’ claimed exhaustive search
(?) was incomplete.

The matrices that are now known as good matrices were first used in the PhD
thesis of Jennifer Seberry Wallis (?); we give their formal definition in
Section 2. She gave a construction using good matrices that enabled the
construction of _Hadamard matrices_ —square matrices with $\pm 1$ entries and
whose rows are pairwise orthogonal. In her thesis she gave examples of good
matrices in all odd orders up to $15$ as well as order $19$ and shortly later
gave an example in order $23$ (?). Subsequently, (?) ran a complete search in
the odd orders up to $21$ and gave examples in order $25$, (?) ran a complete
search up to order $31$, (?) gave examples in orders $33$ and $35$, (?) ran a
complete search up to order $39$ and most recently (?) ran a complete search
up to order $49$. The result of these searches have shown that good matrices
exist in all odd orders up to $39$ as well as $43$ and $45$. Additionally, (?)
gave a construction showing they exist in order $127$.

In this paper we extend the search to all odd orders divisible by $3$ up to
order $69$. Note that increasing $n$ by $2$ increases the size of the search
space (at least to a first approximation) by a factor of $16$, meaning that
the size of the search space in order $69$ is an enormous $2^{40}$ times
larger than the search space in order $49$, the largest order which has
previously been exhaustively searched. There are two primary reasons our
approach is able to scale to such orders:

  1. 1.

We use a number of theoretical mathematical results to cut down the search
space, including filtering theorems, recently proven compression theorems, and
a new product theorem given in this paper.

  2. 2.

We use state-of-the-art programmatic SAT solvers, computer algebra systems,
and mathematical libraries to very efficiently search spaces specified by both
logical and algebraic constraints.

This approach of coupling computer algebra systems and SAT solvers was first
proposed in 2015 at the conferences CADE and ISSAC (?; ?) and has since
appeared in papers at the conferences IJCAI, ISSAC, and the journal of
automated reasoning (?; ?; ?) and is the aim of the SC2 project (?). In
particular, the SAT+CAS method in this paper is an adaption of one presented
at AAAI (?).

Furthermore, we employ a SAT solver which can learn clauses _programmatically_
, through a piece of code compiled with the SAT solver and which makes calls
to an external mathematical library. This approach permits one to efficiently
encode constraints that would be much too cumbersome to encode with native SAT
clauses, as we will see in Section 3. For example, it is much faster and
easier to compute a discrete Fourier transform using a numerical library or
computer algebra system than it would be to encode a DFT circuit using logical
clauses, the format normally accepted by SAT solvers.

A detailed description of our enumeration algorithm will be given in Section
4, followed by our results in Section 5. In particular, we explicitly give two
new examples of good matrices, one in order $27$ which has gone undetected
since 1988 and one in order $57$, an order which had previously been out of
reach of exhaustive search algorithms; the largest order searched prior to our
work was $49$. We also provide a table enumerating the number of sets of good
matrices in all odd orders divisible by $3$ up to $69$. These counts are given
up to an equivalence which is described, along with the required background on
good matrices, in Section 2.

##  2 Background

In this section we define good matrices and review the properties and theorems
of good matrices that we use in our enumeration algorithm.

To begin, we recall the definitions of symmetric, skew, and circulant
matrices. Let $X$ be a square matrix of order $n$ with entries $x_{i,j}$ with
$0\leq i,j<n$. We say that $X$ is _symmetric_ if all entries satisfy
$x_{i,j}=x_{j,i}$, that $X$ is _skew_ if diagonal entries are $1$ and non-
diagonal entries satisfy $x_{i,j}=-x_{j,i}$, and that $X$ is _circulant_ if
all entries satisfy $x_{i,j}=x_{i-1,j-1}$ with indices taken mod $n$ as
necessary.

### Good matrices

######  Definition 1.

Four matrices $A$, $B$, $C$, $D\in\\{\pm 1\\}^{n\times n}$ are known as _good
matrices_ if they have the following properties.

  1. (a)

They are pairwise amicable ($XY^{T}$ is symmetric for $X,Y\in\\{A,B,C,D\\}$).

  2. (b)

$A$ is skew and $B$, $C$, $D$ are symmetric.

  3. (c)

$AA^{T}+B^{2}+C^{2}+D^{2}$ is the identity matrix scaled by $4n$.

The primary reason that good matrices have attracted sustained interest for
almost 50 years is because of the following theorem of (?) which provides a
construction for skew Hadamard matrices.

######  Theorem 1.

Let $A$, $B$, $C$, $D$ be good matrices of order $n$. Then

$\begin{bmatrix}A&B&C&D\\\ -B&A&D&-C\\\ -C&-D&A&B\\\ -D&C&-B&A\end{bmatrix}$

is a skew Hadamard matrix of order $4n$.

Skew Hadamard matrices are conjectured to exist in all orders $4n$ with $n\geq
1$ and much effort has gone into constructing them in as many orders as
possible; the current smallest unknown order is $n=69$ and the previous
smallest unknown order $n=47$ was solved ten years ago by (?).

From a computational perspective the issue with good matrices is that the
search space is too large to effectively search. The standard remedy for this
is to focus on _circulant_ matrices for which the search can be performed more
effectively. For example, in (?) good matrices are defined as in Definition 1
except that condition (a) is replaced with the strictly stronger condition

  1. (a′)

$A$, $B$, $C$, $D$ are circulant.

Technically matrices which satisfy conditions (a′), (b), and (c) are not good
matrices but one can recover good matrices from them simply by reversing the
order of the rows of $B$, $C$, and $D$ (thereby making them amicable with
$A$). Since using (a′) is convenient to state our results we will employ it
for the remainder of this paper.

We often consider good matrices in terms of their first rows
$(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$, $(c_{0},\dotsc,c_{n-1})$,
and $(d_{0},\dotsc,d_{n-1})$ which we call their _defining rows_. Since $A$ is
skew and $B$, $C$, and $D$ are symmetric, the defining rows of a set of good
matrices must satisfy $a_{0}=1$ and

$a_{i}=-a_{n-i},\quad b_{i}=b_{n-i},\quad c_{i}=c_{n-i},\quad d_{i}=d_{n-i}$

for $1\leq i<n/2$. Furthermore, we often use $A$, $B$, $C$, and $D$ to denote
the defining rows themselves, in which case we say that $A$ is a skew sequence
and that $B$, $C$, and $D$ are symmetric sequences.

### Equivalence operations

The following three invertible operations are well-known and can be applied to
a set of good matrices $A$, $B$, $C$, $D$ to produce another set of good
matrices we consider _equivalent_ to $A$, $B$, $C$, $D$.

  1. 1.

Reorder $B$, $C$, $D$ in any way.

  2. 2.

Negate any of $B$, $C$, $D$.

  3. 3.

Apply an automorphism of $\mathbb{Z}_{n}=\\{0,\dotsc,n-1\\}$ to the indices of
the defining rows of each of $A$, $B$, $C$, and $D$.

For example, the negation equivalence operation implies that any set of good
matrices is equivalent to one with $b_{0}=c_{0}=d_{0}=1$.

### Filtering theorems

The search space for good matrices of order $n$ is enormous. Even if one takes
into account the symmetry properties of the defining rows and fixes the first
entries to be positive, each defining row contains $\lfloor n/2\rfloor$
unspecified entries. Since each unspecified entry can be one of two values,
there are a total of $2^{4\lfloor n/2\rfloor}=4^{n-1}$ possibilities in the
search space.

In order to cut down the size of the search space so that we can run
exhaustive searches in larger orders we need filtering theorems. One of the
most powerful filtering theorems is based on the following which gives an
alternative characterization of good matrices (?).

######  Theorem 2.

Let $A$, $B$, $C$, $D$ be sequences of length $n$ with $A$ skew and $B$, $C$,
$D$ symmetric. Then $A$, $B$, $C$, $D$ are the defining rows of a set of good
matrices if and only if for all integers $k$ we have

$\operatorname{PSD}_{A}(k)+\operatorname{PSD}_{B}(k)+\operatorname{PSD}_{C}(k)+\operatorname{PSD}_{D}(k)=4n.$

Here
$\operatorname{PSD}_{X}(k)\coloneqq\big{\lvert}\sum_{j=0}^{n-1}x_{j}\omega^{jk}\big{\rvert}^{2}$
with $\omega\coloneqq\exp(2\pi i/n)$ is the power spectral density of
$X=[x_{0},\dotsc,x_{n-1}]$.

This is important because of the following well-known corollary:

######  Corollary 3.

Let $A$, $B$, $C$, $D$ be the defining rows of a set of good matrices of order
$n$. If $S$ is a subset of $\\{A,B,C,D\\}$ then

$\sum_{X\in S}\operatorname{PSD}_{X}(k)\leq 4n$

for all integers $k$.

For example, this corollary says that if $\operatorname{PSD}_{A}(k)>4n$ for
some $k$ then $A$ cannot be part of a set of good matrices. In other words, we
can filter $A$ and ignore it in our enumeration method. A second well-known
corollary of Theorem 2 comes from setting $k=0$ in the statement of the
theorem:

######  Corollary 4.

Let $A$, $B$, $C$, $D$ be the defining rows of a set of good matrices of odd
order $n$. Then

$\operatorname{sum}(B)^{2}+\operatorname{sum}(C)^{2}+\operatorname{sum}(D)^{2}=4n-1$

where $\operatorname{sum}(X)$ denotes the sum of entries in $X$.

Since $\operatorname{sum}(X)$ is an integer when $X$ has integer entries this
corollary tells us that every set of good matrices gives a decomposition of
$4n-1$ into a sum of three integer squares which is a rather restrictive
condition. For example, when $n$ is $69$ there are only three ways to write
$4n-1$ as a sum of exactly three squares of nonnegative integers, namely,

$1^{2}+7^{2}+15^{2},\quad 5^{2}+5^{2}+15^{2},\quad\text{and}\quad
5^{2}+9^{2}+13^{2}.$

One must also consider decompositions with squares of negative integers but
the negative integers can be found with the following simple (often unstated)
lemma. It gives a criterion to determine the signs of $\operatorname{sum}(B)$,
$\operatorname{sum}(C)$, and $\operatorname{sum}(D)$.

######  Lemma 5.

Let $(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$,
$(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$ be the defining rows of a
set of good matrices $A$, $B$, $C$, $D$ of odd order $n$ with
$b_{0}=c_{0}=d_{0}=1$. Then

$\operatorname{sum}(B)\equiv\operatorname{sum}(C)\equiv\operatorname{sum}(D)\equiv
n\pmod{4}.$

### Compression

Although the filtering theorems greatly decrease the number of possibilities
for the defining rows of good matrices there are still typically a large
number of possibilities which are not filtered by Corollaries 3 and 4. An
effective way of shrinking the search space even farther is to apply
properties that “compressions” of good matrices must satisfy. Such theorems
are only applicable in composite orders; we will assume that the order $n$ is
a multiple of $3$. To our knowledge this is the first time compression has
been used in the search for good matrices and it allows us to search larger
orders than other methods. For example, (?) use the aforementioned filtering
corollaries but do not apply compression and are only able to search orders up
to $49$.

The _3-compression_ of a sequence $X=(x_{0},\dotsc,x_{n-1})$ of length $n$ is
the sequence $X^{\prime}$ of length $m\coloneqq n/3$ whose $k$th entry is
$x_{k}+x_{k+m}+x_{k+2m}$. Since we are concerned with the case where the
entries of $X$ are $\pm 1$ the entries of $X^{\prime}$ are either $\pm 1$ or
$\pm 3$. Compression is useful because (?) showed that Theorem 2 holds for
compressions of the defining rows of good matrices.

######  Theorem 6.

Let $A$, $B$, $C$, $D$ be the defining rows of a set of good matrices of order
$n$ and let $A^{\prime}$, $B^{\prime}$, $C^{\prime}$, $D^{\prime}$ be
compressions of those rows. Then for all integers $k$ we have

$\operatorname{PSD}_{A^{\prime}}(k)+\operatorname{PSD}_{B^{\prime}}(k)+\operatorname{PSD}_{C^{\prime}}(k)+\operatorname{PSD}_{D^{\prime}}(k)=4n.$

### A new product theorem

Our study of good matrices uncovered the following relationship that the
entries of good matrices must satisfy. This is an analog of the product
theorem that was proven by John Williamson for Williamson matrices (?).

######  Theorem 7.

Let $(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$,
$(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$ be the defining rows of a
set of good matrices $A$, $B$, $C$, $D$ of odd order $n$. Then
$a_{k}b_{k}c_{k}d_{k}=-a_{2k\bmod n}b_{0}c_{0}d_{0}$ for all $1\leq k<n$.

Let $\bar{x}\coloneqq(1-x)/2$, i.e., the mapping $x\mapsto\bar{x}$ is the
group isomorphism between $\mathbb{Z}_{4}^{*}=\\{\pm 1\\}$ and
$\mathbb{Z}_{2}=\\{0,1\\}$. Our enumeration method uses Theorem 7 in the
following form.

######  Corollary 8.

Let $(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$,
$(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$ be the defining rows of a
set of good matrices $A$, $B$, $C$, $D$ of odd order $n$ with
$b_{0}=c_{0}=d_{0}=1$. Then for $1\leq k<n/2$ we have in $\mathbb{Z}_{2}$

$\bar{a}_{k}+\bar{a}_{2k}+\bar{b}_{k}+\bar{c}_{k}+\bar{d}_{k}=1$

and when $n=3m$ we have $\bar{b}_{m}+\bar{c}_{m}+\bar{d}_{m}=1$.

We give proofs of Lemma 5, Theorem 7, and Corollaries 3, 4, and 8 in the
appendix.

##  3 The SAT+CAS paradigm

The SAT+CAS paradigm is a new method of solving certain kinds of computational
problems which harnesses the power of two fields of computer science:
satisfiability checking and symbolic computation. The paradigm is particularly
useful for problems that have a significant search component and a rich
mathematical component.

Briefly, this is due to the fact that modern SAT solvers have extremely
efficient search procedures but lack the domain knowledge to effectively deal
with rich mathematical concepts. On the other hand, computer algebra systems
like Maple, Mathematica, and SageMath as well as special-purpose numerical
libraries can very effectively deal with rich mathematics.

The fact that satisfiability checking and symbolic computation had great
potential synergy was first pointed out by Erika Ábrahám in an invited talk at
the conference ISSAC in 2015 (?). At almost the same time this synergy was
demonstrated by the system MathCheck presented at the conference CADE (?). The
system MathCheck coupled a SAT solver with a computer algebra system and
solved open cases of two conjectures in graph theory and was later extended to
solve open cases in combinatorial conjectures (?). Since then, the
$\text{SC}^{2}$ project (?) has organized an annual workshop on this topic for
the last three years.

To illustrate the usefulness of the SAT+CAS paradigm in the search for good
matrices specifically, consider Corollaries 3, 4, and 8 from Section 2. These
corollaries all state various properties that good matrices must satisfy but
each property arises from a different branch of mathematics. Consequently,
each property is best dealt with using systems which have been optimized to
deal with that particular branch:

  * •

Corollary 3 is a statement in Fourier analysis and is best checked using a
system which has been designed to efficiently compute Fourier transforms.

  * •

Corollary 4 (and Lemma 5) are number theoretic statements and possibilities
for the values $\operatorname{sum}(B)$, $\operatorname{sum}(C)$, and
$\operatorname{sum}(D)$ can be found by using a system which has been designed
to solve quadratic Diophantine systems.

  * •

Corollary 8 is a statement using arithmetic in $\mathbb{Z}_{2}$ and can
efficiently be encoded directly in a SAT instance.

### The programmatic SAT paradigm

A programmatic SAT solver, introduced in (?), is a new kind of SAT solver
which allows solving satisfiability problems with more expressiveness than the
kinds of problems that can be solved with a standard SAT solver. Briefly, a
programmatic SAT solver will use some custom code to examine the current
partial assignment and if it cannot be extended into a full assignment (based
on some known theorem from the problem domain) a clause will be learned
encoding that fact. Because one can think of “calling” the programmatic SAT
solver with the custom code this code is known as a “callback” function.

Programmatic SAT solvers are conceptually similar to SMT (SAT modulo theories)
solvers and in some ways could be considered a “poor man’s SMT”. However, they
also differ from SMT solvers in some key ways:

#### Simplicity.

SMT solvers are necessarily more complicated than SAT solvers because they
solve problems in first-order logic instead of problems in propositional
logic. This increase in complexity is justifiable for many problems but some
problems contain mostly just propositional logic with a handful of “theory
lemmas” that are cumbersome to state in propositional logic but can be given
to a programmatic SAT solver as necessary to help guide its search.

#### Flexibility.

Typically SMT solvers only support theories which are fixed in advance; for
example, the SMT-LIB standard (?) specifies a fixed number of theories and a
common input and output language for those theories. However, not all problems
can naturally be expressed in those theories; e.g., the current SMT-LIB
standard does not include theories involving Fourier transforms (or
transcendental functions or complex numbers). In contrast, a programmatic SAT
solver can learn theory lemmas derived using Fourier transforms so long as an
appropriate CAS or library can be called.

#### Tailored solving.

Programmatic SAT solvers can be tailored to solve specific classes of problems
(in addition to using the same state-of-the-art techniques which make modern
SAT solving so efficient). For example, in this paper we develop a SAT solver
tailored to searching for good matrices. In principle SMT solvers could work
in the same way but we are not aware of any SMT solvers that offer a
programmatic interface like this.

### A programmatic SAT encoding of good matrices

We now describe how our programmatic SAT solver encoding of good matrices
works in detail. We encode the entries of the defining rows of a set of good
matrices $A$, $B$, $C$, $D$ of order $n$ using $4n$ Boolean variables where we
let true values denote $1$ and false values denote $-1$. By abuse of notation
we use the same names for the Boolean variables and the $\pm 1$ entries that
they represent but it will be clear from the context if we are referring to a
Boolean or integer value.

In fact, using the symmetry constraints from Section 2 we only need to define
the $2(n+1)$ variables $a_{i}$, $b_{i}$, $c_{i}$, $d_{i}$ that have indices
$i$ with $0\leq i<n/2$. In what follows we implicitly use this whenever
necessary; i.e., any variables with indices larger than $n/2$ are used for
clarity only.

A programmatic SAT solver includes a callback function that is run whenever
the SAT solver’s usual conflict analysis fails to find a conflict. The
callback function examines the current partial assignment and learns clauses
that block the current assignment (and extensions of the current assignment)
from the search in the future.

If $x$ is a variable that appears in the current partial assignment we let
$x^{\text{cur}}$ denote either the literal $x$ or $\lnot x$, whichever is true
under the current assignment. The callback function used in our encoding of
good matrices works by encoding the property given in Corollary 3 and its
steps are described in detail below.

  1. 1.

Let $S$ be a list containing the sequences $(a_{0},\dotsc,a_{n-1})$,
$(b_{0},\dotsc,b_{n-1})$, $(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$
that have all their variables currently assigned.

  2. 2.

If $S$ is empty then return control to the SAT solver.

  3. 3.

Let $X$ be the first sequence in $S$ and compute the power spectral density of
$X$. If $\operatorname{PSD}_{X}(k)>4n$ for some $k$ then learn a clause
blocking the sequence $X$ from the search in the future. Explicitly, learn the
clause

$\lnot x_{0}^{\text{cur}}\lor\lnot x_{1}^{\text{cur}}\lor\dotsb\lor\lnot
x_{n-1}^{\text{cur}}$

and return control to the SAT solver.

  4. 4.

If $S$ contains at least two sequences then let $Y$ be the second sequence in
$S$ and compute the power spectral density of $Y$. If
$\operatorname{PSD}_{X}(k)+\operatorname{PSD}_{Y}(k)>4n$ for some $k$ then
learn a clause blocking the sequences $X$ and $Y$ from occurring together
again in the future. Explicitly, learn the clause

$\lnot x_{0}^{\text{cur}}\lor\dotsb\lor\lnot x_{n-1}^{\text{cur}}\lor\lnot
y_{0}^{\text{cur}}\lor\dotsb\lor\lnot y_{n-1}^{\text{cur}}$

and return control to the SAT solver.

  5. 5.

If $S$ contains at least three sequences then let $Z$ be the third sequence in
$S$ and compute the power spectral density of $Z$. If
$\operatorname{PSD}_{X}(k)+\operatorname{PSD}_{Y}(k)+\operatorname{PSD}_{Z}(k)>4n$
for some $k$ then learn a clause blocking the sequences $X$, $Y$, and $Z$ from
occurring together again in the future and return control to the SAT solver.

  6. 6.

If $S$ contains four sequences then we know every entry of $A$, $B$, $C$, and
$D$. If the relationship in Theorem 2 holds then record $A$, $B$, $C$, and $D$
as rows that define good matrices. Whether or not the sequences define good
matrices learn a clause blocking these sequences from occurring together again
in the future.

If only a single instance of good matrices is desired then one can stop
searching in Step 6 if a satisfying $A$, $B$, $C$, $D$ is found. However, we
learn a blocking clause and continue the search because we want to provide an
exhaustive search.

The programmatic approach was essential to our algorithm because it allowed us
to easily and efficiently apply Corollary 3 which would otherwise be difficult
to apply in a SAT instance. Corollary 3 was also extremely effective; the SAT
solver was usually able to learn a clause in Step 3 when enough variables were
assigned. This enormously cut down the search space—the SAT solver could often
learn conflicts with just $n$ variables instead of conflicts with all $4n$
variables and in practice this appeared to give us an exponential speedup in
$n$.

However, using this method alone was not sufficient for us to derive results
beyond order $21$ and a more efficient enumeration method (using the
programmatic encoding as a subroutine) was required to scale to order $69$. We
describe this method in the next section.

##  4 Enumeration algorithm

Although SAT solvers are useful tools for searching large spaces they have
their limits and the main reason that we were not able to derive any new
results only using SAT solvers is because the search space is too large, even
using the programmatic filtering encoding.

One way that the performance of SAT solvers can be improved on instances with
large search spaces is by splitting the instance into many smaller instances
of approximately equal difficulty. The cube-and-conquer paradigm (?) does this
and has achieved some impressive successes such as solving the Boolean
Pythagorean triples problem (?) and determining the fifth Schur number (?).

Our enumeration algorithm will use a similar paradigm to cube-and-conquer
except that instead of splitting the search space into _cubes_ we split the
search space into _compressions_. Not only does splitting via compressions
increase the performance of the SAT solver it also allows for easy
parallelization because instances generated in this way are independent.
Additionally, we use Theorem 6 to apply filtering to the compressions which
also has the effect of further decreasing the size of the search space.

We now describe our algorithm in four steps followed by some postprocessing.
Given an odd order $n$ divisible by $3$ our algorithm enumerates all
inequivalent circulant good matrices of order $n$ with positive first entries.
Each step is designed to search through a larger subspace than the previous
step: step 1 finds every possible rowsum of a good matrix, step 2 finds every
possible $3$-compressed row of a good matrix, step 3 finds every possible
$3$-compressed quadruple of good matrices, and step 4 finds all good matrices
up to equivalence.

### Step 1: Find possible rowsums

If $A$, $B$, $C$, $D$ are good matrices then from Corollary 4 and Lemma 5 the
rowsums of the matrices $B$, $C$, and $D$ must satisfy the quadratic
Diophantine system

$\displaystyle x^{2}+y^{2}+z^{2}=4n-1$ $\displaystyle x\equiv y\equiv z\equiv
n\pmod{4}$

where $x$, $y$, and $z$ are the rowsums of $B$, $C$, and $D$.

Many computer algebra systems contain number theoretic functions that allow
this Diophantine system to be easily solved. For example, the Mathematica
function PowersRepresentations or the Maple script nsoks (?).

### Step 2: Find possible compressed rows

In this step we will generate sets $S_{\text{sk}}$ and $S_{\text{sy}}$ where
$S_{\text{sk}}$ will contain all possible $3$-compressions of skew defining
rows and $S_{\text{sy}}$ will contain all possible $3$-compressions of
symmetric defining rows. These sets will be generated by an essentially brute-
force enumeration through all skew or symmetric $\\{\pm 1\\}$-sequences of
length $n$, using Corollaries 3 and 4 to filter sequences which cannot be a
defining row.

Let $d\coloneqq\lfloor n/2\rfloor$. In detail, we do the following:

  1. 1.

For all $X\in\\{\pm 1\\}^{d}$:

    1. (a)

Let $A$ be the skew sequence whose first entry is $1$ and next $d$ entries are
given by $X$, i.e.,
$A\coloneqq(1,x_{0},\dotsc,x_{d-1},-x_{d-1},\dotsc,-x_{0})$.

    2. (b)

If $\operatorname{PSD}_{A}(k)\leq 4n$ for all $k$ then add the $3$-compression
of $A$ to $S_{\text{sk}}$.

    3. (c)

Let $B$ be the symmetric sequence whose first entry is $1$ and next $d$
entries are given by $X$, i.e.,
$B\coloneqq(1,x_{0},\dotsc,x_{d-1},x_{d-1},\dotsc,x_{0})$.

    4. (d)

If $\operatorname{PSD}_{B}(k)\leq 4n$ for all $k$ and $\operatorname{sum}(B)$
is one of the solutions $x$, $y$, $z$ from step 1 then add the $3$-compression
of $B$ to $S_{\text{sy}}$.

### Step 3: Find possible compressed quadruples

In step 2 we generated all possible 3-compressions of the defining rows of
good matrices. However, most quadruples of possible sequences found in step 2
will violate the relationship in Theorem 6 and can therefore be filtered. In
this step we identify all quadruples of 3-compressions which can’t be filtered
and store them in a set $S_{\text{q}}$. We do this by examining all
$(A^{\prime},B^{\prime},C^{\prime},D^{\prime})$ in $S_{\text{sk}}\times
S_{\text{sy}}^{3}$ and adding the quadruples that satisfy the relationship in
Theorem 6 to $S_{\text{q}}$.

### Step 4: Uncompress via programmatic SAT

In step 3 we generated all $3$-compressed quadruples
$(A^{\prime},B^{\prime},C^{\prime},D^{\prime})$ that satisfy the relationship
in Theorem 6. Given such a quadruple we need to determine if it is possible to
_uncompress_ the quadruple. In other words, we want to find all quadruples
$(A,B,C,D)$ of defining rows of good matrices (if any) whose $3$-compression
is $(A^{\prime},B^{\prime},C^{\prime},D^{\prime})$. To do this we formulate
the uncompression problem as a programmatic SAT problem.

We already discussed our programmatic encoding of the constraint saying that
$A$, $B$, $C$, $D$ define rows of good matrices of order $n=3m$. The remaining
constraints can easily be expressed in the standard format accepted by SAT
solvers (conjunctive normal form). We now discuss how to encode the constraint
that they are $3$-compressions of $A^{\prime}$, $B^{\prime}$, $C^{\prime}$,
and $D^{\prime}$. (For simplicity we only consider the compression constraints
of $A^{\prime}=(a^{\prime}_{0},\dotsc,a^{\prime}_{m-1})$ but the other
constraints are analogous.)

Consider the entry $a^{\prime}_{k}=a_{k}+a_{k+m}+a_{k+2m}$. There are four
possible cases:

  1. 1.

$a^{\prime}_{k}=3$. In this case $a_{k}=a_{k+m}=a_{k+2m}=1$ which is encoded
as the three unit clauses

$a_{k},\quad a_{k+m},\quad a_{k+2m}.$

  2. 2.

$a^{\prime}_{k}=1$. In this case two of $\\{a_{k},a_{k+m},a_{k+2m}\\}$ are
true and one is false. This is encoded as the four clauses

$\displaystyle\lnot a_{k}\lor\lnot a_{k+m}\lor\lnot a_{k+2m},$ $\displaystyle
a_{k}\lor a_{k+m},\quad a_{k}\lor a_{k+2m},\quad a_{k+m}\lor a_{k+2m}.$

  3. 3.

$a^{\prime}_{k}=-1$. In this case one of $\\{a_{k},a_{k+m},a_{k+2m}\\}$ is
true and two are false. This is encoded as the four clauses

$\displaystyle a_{k}\lor a_{k+m}\lor a_{k+2m},$ $\displaystyle\lnot
a_{k}\lor\lnot a_{k+m},\quad\\!\\!\lnot a_{k}\lor\lnot
a_{k+2m},\quad\\!\\!\lnot a_{k+m}\lor\lnot a_{k+2m}.$

  4. 4.

$a^{\prime}_{k}=-3$. In this case $a_{k}=a_{k+m}=a_{k+2m}=-1$ which is encoded
as the three unit clauses

$\lnot a_{k},\quad\\!\\!\lnot a_{k+m},\quad\\!\\!\lnot a_{k+2m}.$

Finally, we discuss how to encode Corollary 8. The simplest case of Corollary
8 tells us that an even number of $\\{b_{m},c_{m},d_{m}\\}$ are true. This is
encoded as the four clauses

$\displaystyle b_{m}\lor c_{m}\lor d_{m},\quad\lnot b_{m}\lor\lnot c_{m}\lor
d_{m},$ $\displaystyle\lnot b_{m}\lor c_{m}\lor\lnot d_{m},\quad
b_{m}\lor\lnot c_{m}\lor\lnot d_{m}.$

The general case of Corollary 8 is an XOR constraint of length $5$. While SAT
solvers are not good with arbitrary XOR constraints we found that these
constraints were short enough that they could be effectively encoded as the
sixteen clauses with an even number of negated literals, i.e.,

$a_{k}\lor a_{2k}\lor b_{k}\lor c_{k}\lor d_{k},\dotsc,a_{k}\lor\lnot
a_{2k}\lor\lnot b_{k}\lor\lnot c_{k}\lor\lnot d_{k}.$

### Postprocessing: Remove equivalent good matrices

After all the SAT instances generated in step 4 are solved we will have a list
of all circulant good matrices of order $n$ up to equivalence. However, some
of them may be equivalent to each other since we have not encoded all of the
equivalence operations of Section 2. In order to provide a list which contains
only inequivalent good matrices it is necessary to test the matrices for
equivalence and remove any that are found to be equivalent to one which has
previously been found.

To do this, we define a canonical representative of each equivalence class of
defining rows of good matrices. This way, to check if two good matrices are
equivalent we compute the canonical representative of each and check if they
are equal.

First, consider the first two equivalence operations (reorder and negate).
Given a quadruple of defining rows $(A,B,C,D)$ we apply the negation operation
to set the first entry of each defining row to be positive. Then we apply the
reorder operation to sort $B$, $C$, $D$ in lexicographically increasing order.
This gives us a canonical representative of each equivalence class formed by
the first two equivalence operations which we denote by $M_{A,B,C,D}$.

Next, consider the third equivalence operation (permute indices). Let $\sigma$
be an automorphism of $\mathbb{Z}_{n}$ and let $\sigma$ act on $(A,B,C,D)$ by
permuting the indices of the defining rows. This operation commutes with the
first two operations, so the lexicographic minimum of the set

$S\coloneqq\\{\,M_{\sigma(A,B,C,D)}:\sigma\in\operatorname{Aut}(\mathbb{Z}_{n})\,\\}$

is the lexicographic minimum of all quadruples equivalent to $(A,B,C,D)$ and
is therefore canonical. We compute this by explicitly generating all elements
of $S$ and selecting the lexicographic minimum.

### Optimizations

For the benefit of those who would like to implement our algorithm we now
discuss some optimizations that we found useful.

In the programmatic encoding and in step 2 of our enumeration algorithm we
need to check if there exists an integer $k$ such that
$\operatorname{PSD}_{X}(k)>4n$ for a given sequence $X$. We only need to check
$k$ with $0\leq k<n/2$ because the symmetries of the Fourier transform on real
input $X$ imply that $\operatorname{PSD}_{X}(k)=\operatorname{PSD}_{X}(\pm
k\bmod n)$. In fact, checking $k=0$ is also unnecessary because
$\operatorname{PSD}_{X}(0)=\operatorname{sum}(X)^{2}$ and step 1 will filter
any sequences $X$ with $\operatorname{sum}(X)^{2}>4n$. Also, when checking the
PSD criterion for each $k$ we sort the elements of $S$ so that the sequences
with the largest PSD values appear first, making it easier to learn short
clauses.

In step 3, the matching procedure can be done very efficiently by a string
sorting technique similar to the method used in (?). First we enumerate all
$(A^{\prime},B^{\prime})\in S_{\text{sk}}\times S_{\text{sy}}$ that satisfy
$\operatorname{PSD}_{A^{\prime}}(k)+\operatorname{PSD}_{B^{\prime}}(k)\leq 4n$
for all $k$ and all $(C^{\prime},D^{\prime})\in S_{\text{sy}}^{2}$ that
satisfy
$\operatorname{PSD}_{C^{\prime}}(k)+\operatorname{PSD}_{D^{\prime}}(k)\leq 4n$
for all $k$. At this point we form one list that contains strings that contain
the values of
$\operatorname{PSD}_{A^{\prime}}(k)+\operatorname{PSD}_{B^{\prime}}(k)$ and
another list that contains strings that contain the values of
$4n-(\operatorname{PSD}_{C^{\prime}}(k)+\operatorname{PSD}_{D^{\prime}}(k))$
and look for strings which occur in both lists. This can be done by sorting
the two lists and doing a linear scan through the lists to find duplicates.
The fact that PSDs are real numbers can make this tricky, but we can deal
solely with integers by applying the inverse Fourier transform to the PSD
values; this produces integer values known as PAF values and the relationship
given in Theorem 6 becomes (for $0<k<n$)

$\operatorname{PAF}_{A^{\prime}}(k)+\operatorname{PAF}_{B^{\prime}}(k)+\operatorname{PAF}_{C^{\prime}}(k)+\operatorname{PAF}_{D^{\prime}}(k)=0.$

In fact, since this step was not a bottleneck we computed the PAF values using
the slower but more straightforward relationship
$\operatorname{PAF}_{X}(k)=\sum_{j=0}^{n-1}x_{j}x_{j+k\bmod n}$.

In step 4, before calling the SAT solver we partition the compressed
quadruples into equivalence classes using the reorder and automorphism
equivalence operations from Section 2 and keep one quadruple from each
equivalence class. If any satisfiable SAT instances are removed by this
process they will have equivalent solutions (also with positive first entries)
in the instance that was kept.

##  5 Results

In this section we discuss our implementation of our enumeration algorithm,
the software and hardware we used to run it, and the results that we achieved.
Our code is available from our website uwaterloo.ca/mathcheck.

Step 1 of our algorithm was done using nsoks in Maple (?). Steps 2 and 3 were
done with custom C++ code and the string sorting was done with the GNU core
utility sort. Step 4 was done with the programmatic SAT solver MapleSAT (?)
and the PSD values were computed using the library FFTW (?). Because FFTW uses
floating point values which are inherently inaccurate all comparisons were
checked to a tolerance of $\epsilon=10^{-2}$ which is small but larger than
the precision of the FFT used. The computations were run on a cluster of
64-bit Opteron 2.2GHz and Xeon 2.6GHz processors running CentOS 6.9 and using
500MB of memory. By far the most computationally expensive step of our
algorithm was running the SAT solver. This accounted for 99.9% of the total
running time, underscoring the importance of using a state-of-the-art SAT
solver such as MapleSAT. This step was parallelized across 250 cores and
completed in about two weeks. Despite step 4 being the bottleneck the other
steps were found to be essential could not be removed without significantly
increasing the computation time.

The timings for running our implementation in each odd order divisible by $3$
up to order $70$ are given in Table 1. This table also contains the number of
SAT instances generated in each order and the number of inequivalent good
matrices found (denoted by $\\#G_{n}$).

$n$ | Time (h) | # inst. | $\\#G_{n}$
---|---|---|---
3 | 0.00 | 1 | 1
9 | 0.00 | 2 | 1
15 | 0.00 | 11 | 11
21 | 0.00 | 39 | 10
27 | 0.00 | 186 | 13
33 | 0.01 | 840 | 15
39 | 0.07 | 1934 | 5
45 | 1.91 | 19205 | 4
51 | 11.35 | 23611 | 0
57 | 233.56 | 102402 | 1
63 | 11115.70 | 808642 | 0
69 | 72366.85 | 918940 | 0
Table 1: The running time in hours, number of SAT instances used, and number
of inequivalent good matrices found ($\\#G_{n}$) in each odd order $n<70$
divisible by $3$.

Our exhaustive search in order $27$ produced 13 sets of inequivalent good
matrices, though Szekeres’ classic paper (?) contains only 12 sets of such
matrices and Seberry reports that this was checked by at least one other
researcher (?). Comparing the matrices generated using our method with those
in his paper shows that Szekeres missed a set of good matrices. Furthermore,
we found a new set of good matrices of order $57$; the two new sets of good
matrices produced by our method are given in Figure 1.

+++++++--+++-+-+---++------
+-+---++--+--+-+-+++-------++--+++++++---+-+-++-++--+++-+
+-++-+---+--++++--+---+-++-
+++-+--++---+---+--+-+-++++----++++-+-+--+---+---++--+-++
+-+++---+--+----+--+---+++-
++++---+--+--+--+-+-+----+++--+++----+-+-+--+--+--+---+++
+----++-+---+--+---+-++----
+++-++++++--+-+++-+-++----++--++----++-+-+++-+--++++++-++

Figure 1: The defining rows of new good matrices of order 27 and order 57
where + encodes $1$ and - encodes $-1$.

##  6 Conclusion

In this paper we have demonstrated the applicability of the SAT+CAS and
programmatic SAT paradigms to the problem of finding good matrices in
combinatorics. This problem has been well-studied for almost 50 years and the
algorithm we’ve developed using SAT solvers coupled with computer algebra
systems and mathematical libraries is a practical method of enumerating
circulant good matrices (at least those with orders divisible by $3$). This is
evidenced by the large gap between the prior state-of-the-art and the results
of this paper—prior to this year the largest order enumerated was $39$ and the
orders up to $49$ were only enumerated this year, while our algorithm was able
to enumerate circulant good matrices of order $69$. Since the search space
grows exponentially in the order (approximately like $4^{n}$) this is a much
larger search space than the order $49$ search space.

Additionally, we’ve demonstrated the effectiveness of our algorithm by
constructing two new sets of good matrices, including one that escaped
detection in a 1988 search (?) and a double-check by Koukouvinos in 1999 (?).
It is surprising that this set of good matrices was overlooked for so long
though the fact that the double-check also failed to produce it likely
dissuaded researchers from spending more computational resources to triple-
check the same search in order $27$.

Unfortunately, the fact that there are no known certificates for an exhaustive
search means that it is difficult to verify that a search completed
successfully. In our case, we relied on multiple pieces of software including
Maple, MapleSAT, FFTW, GNU sort, and custom-written C++ code. A bug in any one
of these programs could cause a good matrix to be overlooked. The fact that
our code confirms (and even corrects) the results of previous searches gives
some assurance that it is working as intended. However, verification of our
nonexistence results from multiple independent sources would be ideal.

## Acknowledgments

We thank the reviewers for their detailed comments which improved this paper’s
clarity. The computations were made possible by the high-performance computer
clusters administered by Compute Canada and SHARCNET.

## References

  * [Ábrahám et al. 2016] Ábrahám, E.; Abbott, J.; Becker, B.; Bigatti, A. M.; Brain, M.; Buchberger, B.; Cimatti, A.; Davenport, J. H.; England, M.; Fontaine, P.; Forrest, S.; Griggio, A.; Kroening, D.; Seiler, W. M.; and Sturm, T.  2016\.  $\text{SC}^{2}$: Satisfiability checking meets symbolic computation.  In Intelligent Computer Mathematics: 9th International Conference, CICM 2016, Bialystok, Poland, July 25–29, 2016, Proceedings, 28–43.
  * [Ábrahám 2015] Ábrahám, E.  2015\.  Building bridges between symbolic computation and satisfiability checking.  In Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation, 1–6.  New York: ACM.
  * [Awyzio and Seberry 2015] Awyzio, G., and Seberry, J.  2015\.  On good matrices and skew Hadamard matrices.  In Colbourn, C. J., ed., Algebraic Design Theory and Hadamard Matrices, 13–28.
  * [Barrett, Fontaine, and Tinelli 2016] Barrett, C.; Fontaine, P.; and Tinelli, C.  2016\.  The Satisfiability Modulo Theories Library (SMT-LIB).
  * [Bright et al. 2018] Bright, C.; Kotsireas, I.; Heinle, A.; and Ganesh, V.  2018\.  Enumeration of complex Golay pairs via programmatic SAT.  In Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Computation, 111–118.
  * [Bright, Kotsireas, and Ganesh 2018] Bright, C.; Kotsireas, I.; and Ganesh, V.  2018\.  A SAT+CAS method for enumerating Williamson matrices of even order.  In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, New Orleans, Louisiana, USA, February 2–7, 2018, 6573–6580.
  * [Đoković and Kotsireas 2015] Đoković, D. Ž., and Kotsireas, I. S.  2015\.  Compression of periodic complementary sequences and applications.  Designs, Codes and Cryptography 74(2):365–377.
  * [Đoković and Kotsireas 2018] Đoković, D. Ž., and Kotsireas, I. S.  2018\.  Goethals–Seidel difference families with symmetric or skew base blocks.  Mathematics in Computer Science 12(4):373–388.
  * [Đoković 1993] Đoković, D. Ž.  1993\.  Good matrices of order 33, 35 and 127 exist.  J. Combin. Math. Combin. Comput 14:145–152.
  * [Đoković 2008] Đoković, D. Ž.  2008\.  Skew-Hadamard matrices of orders 188 and 388 exist.  International Mathematical Forum 3(22):1063–1068.
  * [Frigo and Johnson 2005] Frigo, M., and Johnson, S. G.  2005\.  The design and implementation of FFTW3.  Proceedings of the IEEE 93(2):216–231.
  * [Ganesh et al. 2012] Ganesh, V.; O’Donnell, C. W.; Soos, M.; Devadas, S.; Rinard, M. C.; and Solar-Lezama, A.  2012\.  Lynx: A programmatic SAT solver for the RNA-folding problem.  In International Conference on Theory and Applications of Satisfiability Testing, 143–156.  Springer.
  * [Georgiou, Koukouvinos, and Stylianou 2002] Georgiou, S.; Koukouvinos, C.; and Stylianou, S.  2002\.  On good matrices, skew Hadamard matrices and optimal designs.  Comput. Statist. Data Anal. 41(1):171–184.
  * [Heule, Kullmann, and Biere 2018] Heule, M. J.; Kullmann, O.; and Biere, A.  2018\.  Cube-and-conquer for satisfiability.  In Handbook of Parallel Constraint Reasoning.  31–59.
  * [Heule, Kullmann, and Marek 2016] Heule, M. J.; Kullmann, O.; and Marek, V. W.  2016\.  Solving and verifying the Boolean Pythagorean triples problem via cube-and-conquer.  In International Conference on Theory and Applications of Satisfiability Testing, 228–245.
  * [Heule 2018] Heule, M. J.  2018\.  Schur number five.  In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, New Orleans, Louisiana, USA, February 2–7, 2018, 6598–6606.
  * [Hunt 1972] Hunt, D. C.  1972\.  Skew-Hadamard matrices of order less than 100.  In Proceedings of the First Australian Conference on Combinatorial Mathematics, 23–27.
  * [Kotsireas, Koukouvinos, and Pardalos 2010] Kotsireas, I. S.; Koukouvinos, C.; and Pardalos, P. M.  2010\.  An efficient string sorting algorithm for weighing matrices of small weight.  Optimization Letters 4(1):29–36.
  * [Liang et al. 2016] Liang, J. H.; Ganesh, V.; Poupart, P.; and Czarnecki, K.  2016\.  Exponential Recency Weighted Average Branching Heuristic for SAT Solvers.  In Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, 3434–3440.
  * [Riel 2006] Riel, J.  2006\.  nsoks: A Maple script for writing $n$ as a sum of $k$ squares.
  * [Seberry 1999] Seberry, J.  1999\.  Library of good matrices.  www.uow.edu.au/~jennie/good.html.
  * [Szekeres 1988] Szekeres, G.  1988\.  A note on skew type orthogonal $\pm 1$ matrices.  In Combinatorics, Colloquia Mathematica Societatis, Janos Bolyai, volume 52, 489–498.
  * [Wallis 1970] Wallis, J. S.  1970\.  Combinatorial matrices.  Ph.D. Dissertation, La Trobe University.
  * [Wallis 1971] Wallis, J. S.  1971\.  A skew-Hadamard matrix of order 92.  Bulletin of the Australian Mathematical Society 5:203–204.
  * [Williamson 1944] Williamson, J.  1944\.  Hadamard’s determinant theorem and the sum of four squares.  Duke Mathematical Journal 11(1):65–81.
  * [Zulkoski et al. 2017] Zulkoski, E.; Bright, C.; Heinle, A.; Kotsireas, I. S.; Czarnecki, K.; and Ganesh, V.  2017\.  Combining SAT solvers with computer algebra systems to verify combinatorial conjectures.  J. Autom. Reasoning 58(3):313–339.
  * [Zulkoski, Ganesh, and Czarnecki 2015] Zulkoski, E.; Ganesh, V.; and Czarnecki, K.  2015\.  MathCheck: A math assistant via a combination of computer algebra systems and SAT solvers.  In Felty, A. P., and Middeldorp, A., eds., Automated Deduction - CADE-25, volume 9195 of Lecture Notes in Computer Science.  607–622.
  * [Zulkoski, Ganesh, and Czarnecki 2016] Zulkoski, E.; Ganesh, V.; and Czarnecki, K.  2016\.  MathCheck: A math assistant via a combination of computer algebra systems and SAT solvers.  In Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, 4228–4233.

## Appendix: Proofs

######  Theorem 2.

Let $A$, $B$, $C$, $D$ be sequences of length $n$ with $A$ skew and $B$, $C$,
$D$ symmetric. Then $A$, $B$, $C$, $D$ are the defining rows of a set of good
matrices if and only if for all integers $k$ we have

$\operatorname{PSD}_{A}(k)+\operatorname{PSD}_{B}(k)+\operatorname{PSD}_{C}(k)+\operatorname{PSD}_{D}(k)=4n.$

Here
$\operatorname{PSD}_{X}(k)\coloneqq\big{\lvert}\sum_{j=0}^{n-1}x_{j}\omega^{jk}\big{\rvert}^{2}$
with $\omega\coloneqq\exp(2\pi i/n)$ is the power spectral density of
$X=[x_{0},\dotsc,x_{n-1}]$.

###### Proof.

Using the notation of (?) an alternative way of stating condition (c) in our
definition of good matrices is that

$\operatorname{PAF}_{A}+\operatorname{PAF}_{B}+\operatorname{PAF}_{C}+\operatorname{PAF}_{D}=[4n,0,\dotsc,0].$

Then (?, Theorem 2) gives the desired result. ∎

######  Corollary 3.

Let $A$, $B$, $C$, $D$ be the defining rows of a set of good matrices of order
$n$. If $S$ is a subset of $\\{A,B,C,D\\}$ then

$\sum_{X\in S}\operatorname{PSD}_{X}(k)\leq 4n$

for all integers $k$.

###### Proof.

By Theorem 2 we have

$\sum_{X=A,B,C,D}\operatorname{PSD}_{X}(k)=4n.$

The desired inequality follows from the fact that PSD values are nonnegative.
∎

######  Corollary 4.

Let $A$, $B$, $C$, $D$ be the defining rows of a set of good matrices of odd
order $n$. Then

$\operatorname{sum}(B)^{2}+\operatorname{sum}(C)^{2}+\operatorname{sum}(D)^{2}=4n-1$

where $\operatorname{sum}(X)$ denotes the sum of entries in $X$.

###### Proof.

Since $\operatorname{PSD}_{X}(0)=\operatorname{sum}(X)^{2}$ and
$\operatorname{sum}(A)=1$ (as $A$ is skew) the desired equality follows by
setting $k=0$ in Theorem 2. ∎

######  Lemma 5.

Let $(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$,
$(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$ be the defining rows of a
set of good matrices $A$, $B$, $C$, $D$ of odd order $n$ with
$b_{0}=c_{0}=d_{0}=1$. Then

$\operatorname{sum}(B)\equiv\operatorname{sum}(C)\equiv\operatorname{sum}(D)\equiv
n\pmod{4}.$

###### Proof.

When $X$ is symmetric and $x_{0}=1$ we have

$\operatorname{sum}(X)=1+2(x_{1}+\dotsb+x_{(n-1)/2}).$

Since each $x_{k}$ is $\pm 1$ and congruent to $1\pmod{2}$ the sum in
parenthesis above is congruent to $(n-1)/2\pmod{2}$ and twice the sum is
congruent to $n-1\pmod{4}$. Thus $\operatorname{sum}(X)\equiv n\pmod{4}$. ∎

######  Theorem 7.

Let $(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$,
$(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$ be the defining rows of a
set of good matrices $A$, $B$, $C$, $D$ of odd order $n$. Then
$a_{k}b_{k}c_{k}d_{k}=-a_{2k\bmod n}b_{0}c_{0}d_{0}$ for all $1\leq k<n$.

###### Proof.

We may assume that $b_{0}=c_{0}=d_{0}=1$; all sets of good matrices are
equivalent to one in which this holds and the once the theorem has been proven
for this class of good matrices the general theorem follows by applying the
negation equivalence operation appropriately.

Let $P_{X}$ be the matrix whose positive entries are identical to the entries
of $X$ and whose entries are $0$ otherwise. Since $X\in\\{A,B,C,D\\}$ has $\pm
1$ entries we have $P_{X}=(X+J)/2$ where $J$ is a matrix of all $1$s and
$P_{X}=\sum_{x_{i}=1}T^{i}$ where $T$ is the circulant matrix with first row
$(0,1,0,\dotsc,0)$. Let $p_{X}$ be the number of positive entries in the first
row of $X$. In (?, §3.4) it is shown that we have

$-P_{A}^{2}+P_{A}+P_{B}^{2}+P_{C}^{2}+P_{D}^{2}=-(n-p_{B}-p_{C}-p_{D})J+nI.$

Since $n$ is odd and $p_{B}$, $p_{C}$, and $p_{D}$ are odd by the symmetry of
$B$, $C$, and $D$, we have that the right side reduces to $I\pmod{2}$. Since
$P_{X}^{2}\equiv\sum_{x_{i}=1}T^{2i}\pmod{2}$ the left side mod $2$ reduces to

$\sum_{a_{i}=1}(T^{2i}+T^{i})+\sum_{b_{i}=1}T^{2i}+\sum_{c_{i}=1}T^{2i}+\sum_{d_{i}=1}T^{2i}.$

Let $i/2$ denote the unique integer in $\\{0,\dotsc,n-1\\}$ such that
$2(i/2)\equiv i\pmod{n}$ which is defined since $n$ is odd. Putting together
the above, we derive (mod $2$)

$\sum_{a_{i/2}=1}T^{i}+\sum_{a_{i}=1}T^{i}+\sum_{b_{i/2}=1}T^{i}+\sum_{c_{i/2}=1}T^{i}+\sum_{d_{i/2}=1}T^{i}\equiv
I.$

Consider the $i$th entry of the initial row of the matrix on the left where
$1\leq i<n$. This entry will be the sum of the $1$s in the set
$\\{a_{i/2},a_{i},b_{i/2},c_{i/2},d_{i/2}\\}$. The right side says that this
entry is $0\pmod{2}$ and therefore an even number of entries in this set are
$1$ (and an odd number must be $-1$). In other words,

$a_{i/2}a_{i}b_{i/2}c_{i/2}d_{i/2}=-1$

which a rewriting of the required result with the definition $k=i/2$. ∎

Let $\bar{x}\coloneqq(1-x)/2$, i.e., the mapping $x\mapsto\bar{x}$ is the
group isomorphism between $\mathbb{Z}_{4}^{*}=\\{\pm 1\\}$ and
$\mathbb{Z}_{2}=\\{0,1\\}$.

######  Corollary 8.

Let $(a_{0},\dotsc,a_{n-1})$, $(b_{0},\dotsc,b_{n-1})$,
$(c_{0},\dotsc,c_{n-1})$, $(d_{0},\dotsc,d_{n-1})$ be the defining rows of a
set of good matrices $A$, $B$, $C$, $D$ of odd order $n$ with
$b_{0}=c_{0}=d_{0}=1$. Then for $1\leq k<n/2$ we have in $\mathbb{Z}_{2}$

$\bar{a}_{k}+\bar{a}_{2k}+\bar{b}_{k}+\bar{c}_{k}+\bar{d}_{k}=1$

and when $n=3m$ we have $\bar{b}_{m}+\bar{c}_{m}+\bar{d}_{m}=1$.

###### Proof.

Applying the mapping $x\mapsto\bar{x}$ to $a_{k}b_{k}c_{k}d_{k}=-a_{2k}$ from
Theorem 7 gives

$\bar{a}_{k}+\bar{b}_{k}+\bar{c}_{k}+\bar{d}_{k}=1+\bar{a}_{2k}$

which implies the first result since $-1=1$ in $\mathbb{Z}_{2}$. When $k=n/3$
we have $\bar{a}_{2k}=-\bar{a}_{n-2k}=-\bar{a}_{k}$ and the second result
follows. ∎