Spectrum of MATLABs magic squares
This article looks at the eigenvalues of magic squares generated by the MATLAB's magic($n$) function. The magic($n$) function constructs doubly even ($n = 4k$) magic squares, singly even ($n = 4k+2$) magic squares and odd ($n = 2k+1$) magic squares using different algorithms. The doubly even magic squares are constructed by a criss-cross method that involves reflecting the entries of a simple square about the center. The odd magic squares are constructed using the Siamese method. The singly even magic squares are constructed using a lower-order odd magic square (Strachey method). We obtain approximations of eigenvalues of odd and singly even magic squares and prove error bounds on the approximation. For the sake of completeness, we also obtain the eigenpairs of doubly even magic squares generated by MATLAB. The approximation of the spectra involves some interesting connections with the spectrum of g-circulant matrices and the use of Bauer-Fike theorem.
\section{Introduction}
Magic squares are given historical~\cite{anderson2001mathematical}, cultural~\cite{cammann1969islamic} and recreational~\cite{benefiel2012magic} importance. Matrix properties of the magic square are also studied~\cite{loly2009magic}\cite{van1990magic} for various applications in image processing~\cite{liu2012study}, privacy \cite{ranjani2017data},  statistics \cite{hunter2010some}, and geophysics \cite{adetokunbo20163d}. Spectral properties of some special magic squares have also been studied in the past~\cite{loly2009magic},~\cite{mattingly2000even} and ~\cite{lee2006linear}.

In this article, we discuss the eigenvalues of a class of magic squares. More specifically, we obtain approximations of eigenvalues of magic squares generated by MATLAB's magic() command~\cite{moler2004numerical} and prove error bounds on the approximation. The magic() command uses the
\begin{enumerate}
	\item
	Siamese method~\cite{weisstein2002magic} for generating odd magic squares
	\item
	Strachey method~\cite{ball2016mathematical} for generating singly even magic squares
	\item
	A form of \emph{criss-cross} method for generating doubly even magic squares.
\end{enumerate}

\section{Notations and definitions}
The following notations and definitions are used throughout this article.
\begin{enumerate}
	\item
	$e_k$ is a unit column vector whose $k^{th}$ entry is $1$.
	\item
	$\mathbb{I}_{p \times q}$ is a $p \times q$ matrix with all entries as $1$.
	\item
	$\mathbb{I}_{p}$ is a $p \times 1$ column vector with all entries as $1$. Note that $\mathbb{I}_{p \times 1} = \mathbb{I}_p$.
	\item
	$\mathbb{O}_{r \times s}$ is a $r \times s$ matrix with all entries as $0$.
	\item
	$\mathbb{O}_{r}$ is a $r \times 1$ column vector with all entries as $0$. Note that $\mathbb{O}_{r \times 1} = \mathbb{O}_r$.
	\item
	$I_t$ is a $t \times t$ identity matrix and $J_t$ is a $t \times t$ reverse identity matrix, i.e.,
	\begin{align}
		I_t(i,j) & = \begin{cases}
		1 & \text{ if }i=j\\
		0 & \text{ otherwise}
		\end{cases}\\
		J_t(i,j) & = \begin{cases}
		1 & \text{ if }i+j=t+1\\
		0 & \text{ otherwise}
		\end{cases}
	\end{align}
	\item
	A set $P$ is said to be $\epsilon$-approximated by set $Q$ if for each $p \in P$, there exists a $q(p) \in Q$ such that
	\begin{align}
		\dfrac{\abs{p-q}}{\abs{q}} \leq \epsilon
	\end{align}
	\item
	A set $P$ is said to be $\epsilon$-well-approximated by set $Q$ if for each $p \in P$, there exists a $q(p) \in Q$ such that
	\begin{align}
		\dfrac{\abs{p-q}}{\abs{q}} \leq \epsilon
	\end{align}
	and the map $q: P \mapsto Q$ is a bijection.
    \item
    $\omega$ will denote the $n^{th}$ root of unity with the greatest real part and negative imaginary part, i.e., $\omega = e^{-2 \pi i/n}$.
    \item
    The vector $v_j = \begin{bmatrix}
    1 & \omega^j & \omega^{2j} & \cdots & \omega^{(n-1)j}
    \end{bmatrix}^T \in \mathbb{R}^n$ will be termed as the $j^{th}$ Fourier mode. Note that $v_{n-j} = v_{-j}$.
    \item
	$F_n \in \mathbb{C}^{n \times n}$ is the discrete Fourier matrix given by
    \begin{align}
        F_n & = \begin{bmatrix}
        v_0 & v_1 & \cdots & v_{n-1}
        \end{bmatrix}
    \end{align}
	\item
	$M_n$ is the $n \times n$ magic square obtained from MATLAB's magic($n$) command.
	\item
	$\mu_0 = \dfrac{n\bkt{n^2+1}}2$ is the magic sum.
	\item
	$\left\{\mu_i \right\}_{i=0}^{n-1}$ will denote the $n$ eigenvalues of $M_n$.
\end{enumerate}
\section{Main results}
We prove the following new results as part of this article. Note that the magic sum $\mu_0=\dfrac{n\bkt{n^2+1}}2$ is always an eigenvalue for the magic square and $\mathbb{I}_{n \times 1}$ is the associated eigenvector. We will be looking at the remaining $n-1$ non-trivial eigenvalues. The three main results are stated below.
\begin{tcolorbox}[title=Non-trivial eigenvalues and eigenvectors of odd magic squares]
\begin{theorem}
	\label{odd_theorem}
	Let $n=2m+1$ and let $$\Lambda = \left\{\lambda_j = \dfrac{n^2}{2\sin\bkt{j\pi/n}}: j \in\{\pm1,\pm2,\ldots,\pm m\} \right\}$$
	The set of all non-trivial eigenvalues of an odd magic square ($M_n$ with $n=2m+1$) generated using the \textbf{Siamese method} is $\frac1n$-approximated by $\Lambda$.
\end{theorem}
\end{tcolorbox}
\begin{tcolorbox}[title=Non-trivial eigenvalues and eigenvectors of singly even magic squares]
\begin{theorem}
	\label{singly_even_theorem}
	Let $M_n$ be a singly even magic square, where $n=2m=4k+2$, generated using the \textbf{Strachey method}. Let $$\Lambda = \left\{\lambda_j = \dfrac{n^2}{4\sin\bkt{j\pi/m}}: j \in\{\pm1,\pm2,\ldots,\pm k\} \right\}$$
	\begin{enumerate}
		\item
		Two of the non-trivial eigenpairs of $M_n$ are $\bkt{-3m^2,\begin{bmatrix}\mathbb{I}_m\\-\mathbb{I}_m\end{bmatrix}}$ and $\bkt{3m^2,\begin{bmatrix}\mathbb{I}_m - 3e_{k+1}\\-\bkt{\mathbb{I}_m - 3e_{k+1}}\end{bmatrix}}$.
		\item
		$2k-1$ of the non-trivial eigenvalues are zero.
		\item
		The set of remaining $2k$ eigenvalues are $\frac2n$-approximated by the set
		$$\left\{\pm\dfrac{n^2}{4\sin\bkt{\pi/m}},\pm\dfrac{n^2}{4\sin\bkt{2\pi/m}},\ldots,\pm\dfrac{n^2}{4\sin\bkt{k\pi/m}}\right\}$$
	\end{enumerate}
\end{theorem}
\end{tcolorbox}
\begin{tcolorbox}[title=Non-trivial eigenvalues and eigenvectors of doubly even magic squares]
\begin{theorem}
	\label{doubly_even_theorem}
	The rank of the doubly even magic square generated by MATLAB's \textbf{criss-cross} method is $3$ and the non-zero non-trivial eigenvalues are
	$$\pm \dfrac{n}2 \sqrt{\dfrac{n^3-n}3}$$
	The corresponding eigenvectors are $\sqrt{n^2-1}v \pm \sqrt{3n}u$, where
\setcounter{MaxMatrixCols}{12}
	$$u = \begin{bmatrix} w\\
	-\text{flip}\bkt{w} \end{bmatrix}_{4m \times 1}$$ and
	$$v^T = \begin{bmatrix} 1 & -1 & -1 & 1 & \cdots & 1 & -1 & -1 & 1 \end{bmatrix}_{1 \times 4m}$$
	with $w^T = \begin{bmatrix}(n-1) & -(n-3) & -(n-5) & (n-7) & \cdots & (-1)^m\end{bmatrix}$.
	\begin{note}
		This result on doubly even magic squares has already been obtained in~\cite{kirkland1995group}. We however rederive it for the sake of completeness.
	\end{note}
\end{theorem}
\end{tcolorbox}
\section{Odd magic squares}
The first two odd magic squares generated by MATLAB's magic() command are given in Equation~\eqref{odd_magic}.
\begin{align}
M_3 = \begin{bmatrix}
\teal{8} & \red{1} & \blue{6} \\
\red{3} & \blue{5} & \teal{7} \\
\blue{4}  & \teal{9} & \red{2}
\end{bmatrix},
M_5 =  \begin{bmatrix}
\teal{17} & \red{24} & \blue{1} & \brown{8} & \cyan{15} \\
\red{23} & \blue{5} & \brown{7} & \cyan{14} & \teal{16} \\
\blue{4} & \brown{6} & \cyan{13} & \teal{20} & \red{22} \\
\brown{10} & \cyan{12} & \teal{19} & \red{21} & \blue{3} \\
\cyan{11} & \teal{18} & \red{25} & \blue{2} & \brown{9} 
\end{bmatrix}
\label{odd_magic}
\end{align}
Two different MATLAB codes (one using for-loops and the other one completely vectorized) for obtaining odd magic square is presented here~\cite{clevemagic2}. We present the for-loop code below for the sake of clarity and continuity.
\begin{tcolorbox}
	\begin{verbatim}
    A = zeros(n,n);
    i = 1;
    j = (n+1)/2;
    for k = 1:n^2
       is = i;
       js = j;
       A(i,j) = k;
       i = n - rem(n+1-i,n);
       j = rem(j,n) + 1;
       if A(i,j) ~= 0
          i = rem(is,n) + 1;
          j = js;
       end
    end
	\end{verbatim}
\end{tcolorbox}

Before we get to the proof, we will play around with an example to get a better understanding of such odd magic squares. Note that $M_5$ can be written as
\begin{align}
	M_5 & = \begin{bmatrix}
	\teal{17} & 0 & 0 & 0 & 0 \\
	0 & 0 & 0 & 0 & \teal{16} \\
	0 & 0 & 0 & \teal{20} & 0 \\
	0 & 0 & \teal{19} & 0 & 0 \\
	0 & \teal{18} & 0 & 0 & 0
	\end{bmatrix} +
	\begin{bmatrix}
	0 & \red{24} & 0 & 0 & 0 \\
	\red{23} & 0 & 0 & 0 & 0 \\
	 0 & 0 & 0 & 0 & \red{22}\\
	 0 & 0 & 0 & \red{21} & 0\\
	 0 & 0 & \red{25} & 0 & 0
	\end{bmatrix} +
	\begin{bmatrix}
	0 & 0 & \blue{1} & 0 & 0 \\
	0 & \blue{5} & 0 & 0 & 0 \\
	\blue{4} & 0 & 0 & 0 & 0 \\
	0 & 0 & 0 & 0 & \blue{3} \\
	0 & 0 & 0 & \blue{2} & 0
	\end{bmatrix}\\
	& 
	+ \begin{bmatrix}
	0 & 0 & 0 & \brown{8} & 0 \\
	0 & 0 & \brown{7} & 0 & 0 \\
	0 & \brown{6} & 0 & 0 & 0 \\
	\brown{10} & 0 & 0 & 0 & 0\\
	0 & 0 & 0 & 0 & \brown{9} 
	\end{bmatrix} +
	\begin{bmatrix}
	0 & 0 & 0 & 0 & \cyan{15} \\
	0 & 0 & 0 & \cyan{14} & 0 \\
	0 & 0 & \cyan{13} & 0 & 0 \\
	0 & \cyan{12} & 0 & 0 & 0\\
	\cyan{11} & 0 & 0 & 0 & 0
	\end{bmatrix}
\end{align}
Let $J_5$ denotes the reversed $5 \times 5$ identity matrix, i.e.,
\begin{align}J_5 = \begin{bmatrix}
0 & 0 & 0 & 0 & 1\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 1 & 0 & 0\\
0 & 1 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0
\end{bmatrix}
\end{align}
We then have
\begin{align}
M_5 & = J_5\bkt{\begin{bmatrix}
0 & \teal{18} & 0 & 0 & 0\\
0 & 0 & \teal{19} & 0 & 0 \\
0 & 0 & 0 & \teal{20} & 0 \\
0 & 0 & 0 & 0 & \teal{16} \\
\teal{17} & 0 & 0 & 0 & 0 \\
\end{bmatrix} +
\begin{bmatrix}
0 & 0 & \red{25} & 0 & 0\\
0 & 0 & 0 & \red{21} & 0\\
0 & 0 & 0 & 0 & \red{22}\\
\red{23} & 0 & 0 & 0 & 0 \\ 
0 & \red{24} & 0 & 0 & 0 \\ 
\end{bmatrix} +
\begin{bmatrix}
0 & 0 & 0 & \blue{2} & 0\\
0 & 0 & 0 & 0 & \blue{3} \\
\blue{4} & 0 & 0 & 0 & 0 \\
0 & \blue{5} & 0 & 0 & 0 \\
0 & 0 & \blue{1} & 0 & 0 \\
\end{bmatrix}}\\
& +
J_5\bkt{\begin{bmatrix}
0 & 0 & 0 & 0 & \brown{9}\\ 
\brown{10} & 0 & 0 & 0 & 0\\
0 & \brown{6} & 0 & 0 & 0 \\
0 & 0 & \brown{7} & 0 & 0 \\
0 & 0 & 0 & \brown{8} & 0 \\
\end{bmatrix} +
\begin{bmatrix}
\cyan{11} & 0 & 0 & 0 & 0\\
0 & \cyan{12} & 0 & 0 & 0\\
0 & 0 & \cyan{13} & 0 & 0 \\
0 & 0 & 0 & \cyan{14} & 0 \\
0 & 0 & 0 & 0 & \cyan{15} \\
\end{bmatrix}}
\end{align}
Let $C_5$ denote the matrix corresponding to the cycle
$$1\mapsto 2\mapsto 3\mapsto 4\mapsto 5 \mapsto 1$$
i.e.,
\begin{align}
	C_5 = \begin{bmatrix}
0 & 0 & 0 & 0 & 1\\
1 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
\end{bmatrix}
\end{align}
Note that pre-multiplying by $C_5$ permutes the rows along the forward cycle, while post multiplying by $C_5$ permutes the columns along the reverse cycle.
Now note that we have
\begin{align}
M_5 & = J_5\bkt{C_5^3\begin{bmatrix}
0 & 0 & 0 & 0 & \teal{16} \\
\teal{17} & 0 & 0 & 0 & 0 \\
0 & \teal{18} & 0 & 0 & 0\\
0 & 0 & \teal{19} & 0 & 0 \\
0 & 0 & 0 & \teal{20} & 0 \\
\end{bmatrix} +
C_5\begin{bmatrix}
0 & 0 & 0 & \red{21} & 0\\
0 & 0 & 0 & 0 & \red{22}\\
\red{23} & 0 & 0 & 0 & 0 \\ 
0 & \red{24} & 0 & 0 & 0 \\ 
0 & 0 & \red{25} & 0 & 0\\
\end{bmatrix}}\\
& + J_5\bkt{ 
C_5^{-1}\begin{bmatrix}
0 & 0 & \blue{1} & 0 & 0 \\
0 & 0 & 0 & \blue{2} & 0\\
0 & 0 & 0 & 0 & \blue{3} \\
\blue{4} & 0 & 0 & 0 & 0 \\
0 & \blue{5} & 0 & 0 & 0 \\
\end{bmatrix} +
C_5^{-3}
\begin{bmatrix}
0 & \brown{6} & 0 & 0 & 0 \\
0 & 0 & \brown{7} & 0 & 0 \\
0 & 0 & 0 & \brown{8} & 0 \\
0 & 0 & 0 & 0 & \brown{9}\\ 
\brown{10} & 0 & 0 & 0 & 0\\
\end{bmatrix}}\\
& + J_5
\begin{bmatrix}
\cyan{11} & 0 & 0 & 0 & 0\\
0 & \cyan{12} & 0 & 0 & 0\\
0 & 0 & \cyan{13} & 0 & 0 \\
0 & 0 & 0 & \cyan{14} & 0 \\
0 & 0 & 0 & 0 & \cyan{15} \\
\end{bmatrix}
\end{align}
\begin{align}
M_5 & = J_5\bkt{C_5^{-2}\begin{bmatrix}
\teal{16} & 0 & 0 & 0 & 0\\
0 & \teal{17} & 0 & 0 & 0\\
0 & 0 & \teal{18} & 0 & 0\\
0 & 0 & 0 & \teal{19} & 0\\
0 & 0 & 0 & 0 & \teal{20}\\
\end{bmatrix}C_5 +
C_5^{-4}\begin{bmatrix}
\red{21} & 0 & 0 & 0 & 0\\
0 & \red{22} & 0 & 0 & 0\\
0 & 0 & \red{23} & 0 & 0\\
0 & 0 & 0 & \red{24} & 0\\ 
0 & 0 & 0 & 0 & \red{25}\\
\end{bmatrix}C_5^2}\\
& + J_5\bkt{
C_5^{4}\begin{bmatrix}
\blue{1} & 0 & 0 & 0 & 0\\
0 & \blue{2} & 0 & 0 & 0\\
0 & 0 & \blue{3} & 0 & 0\\
0 & 0 & 0 & \blue{4} & 0\\
0 & 0 & 0 & 0 & \blue{5}\\
\end{bmatrix}C_5^3
+
C_5^{2}
\begin{bmatrix}
\brown{6} & 0 & 0 & 0 & 0\\
0 & \brown{7} & 0 & 0 & 0\\
0 & 0 & \brown{8} & 0 & 0\\
0 & 0 & 0 & \brown{9} & 0\\ 
0 & 0 & 0 & 0 & \brown{10}\\
\end{bmatrix}C_5^4}\\
&+
J_5\begin{bmatrix}
\cyan{11} & 0 & 0 & 0 & 0\\
0 & \cyan{12} & 0 & 0 & 0\\
0 & 0 & \cyan{13} & 0 & 0 \\
0 & 0 & 0 & \cyan{14} & 0 \\
0 & 0 & 0 & 0 & \cyan{15} \\
\end{bmatrix}
\end{align}
From the example above and from the MATLAB code~\cite{clevemagic2}, we see that for any odd $n$ of the form $2m+1$, we have
\begin{align}
M_n & = J_n\bkt{\displaystyle\sum_{k=0}^{2m}C_n^{2m-2k}\bkt{D_n+nkI_n}C_n^{m+1+k}}\\
& = J_n\bkt{\displaystyle\sum_{k=0}^{2m}C_n^{2m-2k}D_nC_n^{m+1+k}} + nJ_n\bkt{\displaystyle\sum_{k=0}^{2m}kC_n^{3m+1-k}}
\end{align}
where $I_n \in \mathbb{R}^{n \times n}$ is the identity matrix, $J_n \in \mathbb{R}^{n \times n}$ is the reverse identity matrix, $C_n = \begin{bmatrix}0_{1 \times (n-1)} & 1\\
I_{n-1} & 0_{(n-1) \times 1}\end{bmatrix} \in \mathbb{R}^{n \times n}$ is the matrix corresponding to the cycle
$1 \mapsto 2 \mapsto 3 \mapsto \cdots \mapsto n \mapsto 1$
and $D_n \in \mathbb{R}^{n \times n}$ is a diagonal matrix such that $D_n(i,i) = i$. Let $X_n=\displaystyle\sum_{k=0}^{2m}kC_n^{3m+1-k}$ and $Y_n=\displaystyle\sum_{k=0}^{2m}C_n^{2m-2k}D_nC_n^{m+1+k}$, i.e., we have
\begin{align}
    M_n & = nJ_nX_n + J_n Y_n
    \label{relation}
\end{align}
$X_n$ is a circulant matrix (since it's a sum of circulant matrices) whose first row is
\begin{align}
	\mathbf{a} = \begin{bmatrix}m & m+1 & m+2 & \cdots &2m & 0 & 1 & 2 & \cdots & m-1\end{bmatrix}
\end{align}
and $Y_n = \bkt{\displaystyle\sum_{k=0}^{2m}C_n^{2m-2k}D_nC_n^{m+1+k}}$ is a $\bkt{m+1}$-row circulant~\cite{andrade2020spectra} matrix with first row being
\begin{align}
	\mathbf{b} = \begin{bmatrix}
1 & 3 & \cdots & 2m-1 & 2m+1 & 2 & 4 & \cdots & 2m
\end{bmatrix}
\end{align}

For the $5\times5$ case, we have
\begin{align}
    M_5 & = \begin{bmatrix}
\teal{17} & \red{24} & \blue{1} & \brown{8} & \cyan{15} \\
\red{23} & \blue{5} & \brown{7} & \cyan{14} & \teal{16} \\
\blue{4} & \brown{6} & \cyan{13} & \teal{20} & \red{22} \\
\brown{10} & \cyan{12} & \teal{19} & \red{21} & \blue{3} \\
\cyan{11} & \teal{18} & \red{25} & \blue{2} & \brown{9} 
\end{bmatrix} =
\begin{bmatrix}
\teal{15} & \red{20} & \blue{0} & \brown{5} & \cyan{10}\\
\red{20} & \blue{0} & \brown{5} & \cyan{10} & \teal{15}\\
\blue{0} & \brown{5} & \cyan{10} & \teal{15} & \red{20}\\
\brown{5} & \cyan{10} & \teal{15} & \red{20} & \blue{0}\\
\cyan{10} & \teal{15} & \red{20} & \blue{0} & \brown{5}\\
\end{bmatrix}
+
\begin{bmatrix}
\teal{2} & \red{4} & \blue{1} & \brown{3} & \cyan{5}\\
\brown{3} & \cyan{5} & \teal{2} & \red{4} & \blue{1}\\
\red{4} & \blue{1} & \brown{3} & \cyan{5} & \teal{2}\\
\cyan{5} & \teal{2} & \red{4} & \blue{1} & \brown{3}\\
\blue{1} & \brown{3} & \cyan{5} & \teal{2} & \red{4}\\
\end{bmatrix}
\\
& = \begin{bmatrix}
0 & 0 & 0 & 0 & \bf{1}\\
0 & 0 & 0 & \bf{1} & 0\\
0 & 0 & \bf{1} & 0 & 0\\
0 & \bf{1} & 0 & 0 & 0\\
\bf{1} & 0 & 0 & 0 & 0\\
\end{bmatrix}
\bkt{
5 \begin{bmatrix}
\cyan{2} & \teal{3} & \red{4} & \blue{0} & \brown{1}\\
\brown{1} & \cyan{2} & \teal{3} & \red{4} & \blue{0}\\
\blue{0} & \brown{1} & \cyan{2} & \teal{3} & \red{4}\\
\red{4} & \blue{0} & \brown{1} & \cyan{2} & \teal{3}\\
\teal{3} & \red{4} & \blue{0} & \brown{1} & \cyan{2}\\
\end{bmatrix}
+
\begin{bmatrix}
\blue{1} & \brown{3} & \cyan{5} & \teal{2} & \red{4}\\
\cyan{5} & \teal{2} & \red{4} & \blue{1} & \brown{3}\\
\red{4} & \blue{1} & \brown{3} & \cyan{5} & \teal{2}\\
\brown{3} & \cyan{5} & \teal{2} & \red{4} & \blue{1}\\
\teal{2} & \red{4} & \blue{1} & \brown{3} & \cyan{5}\\
\end{bmatrix}
}
\end{align}

We rewrite Equation~\eqref{relation} as
\begin{align}
    M_n & = J_n \bkt{\overbrace{nX_n + (m+1)\mathbb{I}_n\mathbb{I}_n^T}^{S_n} + \underbrace{Y_n - (m+1)\mathbb{I}_n\mathbb{I}_n^T}_{T_n}}
    \label{UVeqn}
\end{align}
where $\mathbb{I}_n$ is column vector of ones. Note that $S_n$ is a circulant matrix and $T_n$ is $(m+1)$-row circulant matrix.

For the $5 \times 5$ case, we have
\begin{align}
    = \begin{bmatrix}
0 & 0 & 0 & 0 & \bf{1}\\
0 & 0 & 0 & \bf{1} & 0\\
0 & 0 & \bf{1} & 0 & 0\\
0 & \bf{1} & 0 & 0 & 0\\
\bf{1} & 0 & 0 & 0 & 0\\
\end{bmatrix}
\bkt{
\begin{bmatrix}
\cyan{13} & \teal{18} & \red{23} & \blue{3} & \brown{8}\\
\brown{8} & \cyan{13} & \teal{18} & \red{23} & \blue{3}\\
\blue{3} & \brown{8} & \cyan{13} & \teal{18} & \red{23}\\
\red{23} & \blue{3} & \brown{8} & \cyan{13} & \teal{18}\\
\teal{18} & \red{23} & \blue{3} & \brown{8} & \cyan{13}\\
\end{bmatrix}
+
\begin{bmatrix}
\blue{-2} & \brown{0} & \cyan{2} & \teal{-1} & \red{1}\\
\cyan{2} & \teal{-1} & \red{1} & \blue{-2} & \brown{0}\\
\red{1} & \blue{-2} & \brown{0} & \cyan{2} & \teal{-1}\\
\brown{0} & \cyan{2} & \teal{-1} & \red{1} & \blue{-2}\\
\teal{-1} & \red{1} & \blue{-2} & \brown{0} & \cyan{2}\\
\end{bmatrix}
}
\end{align}
The main theorem proved in this section is below (Theorem~\ref{main_odd_theorem}).
\begin{tcolorbox}[title=Non-trivial eigenvalues of odd magic squares]
\begin{theorem}
	\label{main_odd_theorem}
	Let $n=2m+1$ and let $$\Lambda = \left\{\lambda_j = \dfrac{n^2}{2\sin\bkt{j\pi/n}}: j \in\{\pm1,\pm2,\ldots,\pm m\} \right\}$$
	The set of all non-trivial eigenvalues of an odd magic square ($M_n$ with $n=2m+1$) generated using the \textbf{Siamese method} is $\frac1n$-approximated by $\Lambda$.
\end{theorem}
\end{tcolorbox}
The reminder of this section discusses the proof of the above theorem (Theorem~\ref{main_odd_theorem}).
\begin{proof}
From Equation~\eqref{UVeqn}, we have $M_n=J_nS_n+J_nT_n$. \noindent\textbf{We now look at $J_nT_n$ as a perturbation to $J_nS_n$ to get $M_n$ and claim that the eigenvalues of $M_n = J_nS_n+J_nT_n$ are well approximated by the eigenvalues of $J_nS_n$. The eigenvalues and eigenvectors of $J_nS_n$ can be easily computed since $S_n$ is a circulant matrix. We now rely on Bauer-Fike theorem to make the above statement precise.}
\begin{tcolorbox}
\begin{theorem}
\label{Bauer-Fike}
\textbf{Bauer-Fike theorem}: Let $U_n \in \mathbb{C}^{n \times n}$ be a diagonalizable matrix. Let $\Lambda_n \in \mathbb{C}^{n \times n}$ be a diagonal matrix with eigenvalues of $U_n$ and $W_n \in \mathbb{C}^{n \times n}$ be the matrix containing the eigenvectors of $U_n$ such that $U_n = W_n \Lambda_n W_n^{-1}$. Let $\kappa_p\bkt{W_n} = \magn{W_n}_p\magn{W_n^{-1}}_p$ be the condition number of the matrix $W_n$ in $p$-norm. 

If $\mu$ is an eigenvalue of $U_n+V_n$, where $V_n \in \mathbb{C}^{n \times n}$, then there exists $\lambda$, an eigenvalue of $U_n$, such that
\begin{align}
	\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \kappa_p\bkt{W_n} \magn{U_n^{-1}V_n}_p
\end{align}
\end{theorem}
\end{tcolorbox}
We utilize Theorem~\ref{Bauer-Fike}, by choosing $U_n=J_nS_n$ and $V_n=J_nT_n$. The matrix $J_nS_n$ is a reverse circulant matrix and is diagonalizable. The eigenvalues and eigenvectors of $J_nS_n$ are derived in the Appendix (See Note~\ref{note12}). Further, the eigenvectors can be chosen to be orthonormal (See Note~\ref{note12}). Hence, the matrix $W_n$ in Bauer-Fike theorem (Theorem~\ref{Bauer-Fike}) is an orthonormal matrix. Hence from Theorem~\ref{Bauer-Fike}, and choosing $p=2$, we have that if $\mu$ is an eigenvalue of $M_n$, then there exists $\lambda$, an eigenvalue of $J_nS_n$, such that
\begin{align}
	\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \magn{\bkt{J_nS_n}^{-1}\bkt{J_nT_n}}_2 = \magn{S_n^{-1}T_n}_2
\end{align}
\begin{lemma}
	\label{lemma1}
	The eigenvalues of $U_n=J_nS_n$ are given
	\begin{align}
		\lambda_j = \begin{cases}
		\dfrac{n\bkt{n^2+1}}2 & j=0\\
		\pm\dfrac{n^2}{2\sin\bkt{\pi j/n}} & j \in \{1,2,\ldots,m\}
	\end{cases}
	\end{align}
	The eigenvectors of $U_n=J_nS_n$ are orthonormal.
\end{lemma}
\begin{proof}
	From Lemma~\ref{lem1}, Lemma~\ref{lem_ortho_eig} and Note~\ref{note12} in the Appendix, we immediately obtain the result.
\end{proof}
\begin{lemma}
	\label{lemma2}
	\begin{align}
		\magn{S_n^{-1}T_n}_2 = \dfrac1n
	\end{align}
\end{lemma}
\begin{proof}
First note that $X_n\mathbb{I} = \dfrac{n(n-1)}2\mathbb{I}$ and $Y_n\mathbb{I} = \dfrac{n(n+1)}2\mathbb{I}$. We have $S_n = nX_n+\bkt{m+1}\mathbb{I}\mathbb{I}^T$ and $T_n=Y_n-(m+1)\mathbb{I}\mathbb{I}^T$. Now using the fact that the inverse of a circulant matrix is circulant matrix, the product of circulant matrices is another circulant matrix, the Sherman–Morrison formula and some tedious algebra, we obtain
\begin{align}
	S_n^{-1}T_n = \dfrac1{n^2}\bkt{nG-\mathbb{I}\mathbb{I}^T}
\end{align}
where $G$ is a $\bkt{m+1}$-row circulant matrix and the first row of $G$ is given by
\begin{align}
	\begin{bmatrix}0&1&0&0&\cdots & 0\end{bmatrix}
\end{align}
Now we seek bound on $\magn{S_n^{-1}T_n}_2$. We have
\begin{align}
	\magn{S_n^{-1}T_n}_2^2 & = \max_{x \in \mathbb{R}^n \backslash \{0\}} \dfrac{\magn{S_n^{-1}T_nx}_2^2}{\magn{x}_2^2}
\end{align}
We have
\begin{align}
	S_n^{-1}T_nx & = \dfrac1{n^2}\bkt{nGx - \mathbb{I}\mathbb{I}^Tx} = \dfrac1{n^2}\bkt{nGx - \bkt{x_1+x_2+\cdots+x_n}\mathbb{I}}
\end{align}
Hence, we have
\begin{align}
	\magn{S_n^{-1}T_nx}_2^2 & = \dfrac1{n^4}\bkt{\displaystyle\sum_{j=1}^n \bkt{nx_j-\bkt{x_1+x_2+\cdots+x_n}}^2}\\
	& = \dfrac{\displaystyle\sum_{j=1}^n\bkt{\bkt{nx_j}^2 + \bkt{x_1+x_2+\cdots+x_n}^2-2nx_j\bkt{x_1+x_2+\cdots+x_n}}}{n^4}\\
	& = \dfrac{n^2 \displaystyle\sum_{j=1}^n x_j^2 - n\bkt{x_1+x_2+\cdots+x_n}^2}{n^4}\\
	& = \dfrac{\magn{x}_2^2}{n^2} - \dfrac{\bkt{x_1+x_2+\cdots+x_n}^2}{n^3}
\end{align}
Hence, we have
$$\dfrac{\magn{S_n^{-1}T_nx}_2^2}{\magn{x}_2^2} = \dfrac1{n^2} - \dfrac{\bkt{\mathbb{I}^Tx}^2}{n^3\magn{x}_2^2}$$
Hence, we obtain that $\magn{S_n^{-1}T_n}_2 = \dfrac1{n}$.
\end{proof}
\begin{tcolorbox}
Hence, from Theorem~\ref{Bauer-Fike}, Lemma~\ref{lemma1} and Lemma~\ref{lemma2}, we have that if $\mu$ is a non-trivial eigenvalue of $M_n$, then there exists
\begin{align}
	\lambda \in \left\{\pm\dfrac{n^2}{2\sin\bkt{\pi/n}},\pm\dfrac{n^2}{2\sin\bkt{2\pi/n}},\ldots,\pm\dfrac{n^2}{2\sin\bkt{m\pi/n}}\right\}
\end{align}
such that
\begin{align}
	\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac1n
\end{align}
\end{tcolorbox}
This proves Theorem~\ref{main_odd_theorem}.
\end{proof}
\noindent \begin{remark} \label{remark1} We, in fact, numerically observe (see Tables~\ref{magic_3},~\ref{magic_5},~\ref{magic_7},~\ref{magic_9} and~\ref{magic_11}) a stronger version of the above result (Theorem~\ref{main_odd_theorem}), namely, for each eigenvalue of $J_nS_n$, say $\lambda$, there is an eigenvalue of $M_n$, say $\mu$, such that
	\begin{align}
		\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac1n
	\end{align}
	or equivalently the set of non-trivial eigenvalues of an odd magic square ($M_n$ with $n=2m+1$) generated using the \textbf{Siamese method} is $\frac1n$-well-approximated by the set
	\begin{align}
		\left\{\pm\dfrac{n^2}{2\sin\bkt{\pi/n}},\pm\dfrac{n^2}{2\sin\bkt{2\pi/n}},\ldots,\pm\dfrac{n^2}{2\sin\bkt{m\pi/n}}\right\}
	\end{align}
\end{remark}
\noindent \begin{remark} \label{remark2} A slightly stronger form of Theorem~\ref{main_odd_theorem} can be proven by diagonalizing the matrix $J_nS_n$ and making use of Gershgorin circle theorem. However, the result still doesn't prove the observation in Remark~\ref{remark1}.
\end{remark}
Let 
\begin{align}
	e_n & = \displaystyle\max_{i=1,2,\ldots,n-1} \displaystyle\min_{j=\pm1,\pm2,\ldots,\pm m} \dfrac{\abs{\mu_i-\lambda_j}}{\abs{\lambda_j}}
\end{align}
We term $e_n$ as the maximum relative error in eigenvalues between $M_n$ and $J_nS_n$. Figure~\ref{fig_rel_error} plots $e_n$ versus $n$ for odd values of $n$. There is a lot to decipher from Figure~\ref{fig_rel_error}. Some observations are listed below.
\begin{enumerate}
	\item
	The bound of $1/n$ is very conservative.
	\item
	$e_n$ for $n \equiv 3\pmod6$ is much higher than the rest.
	\item
	$e_n$ when $n$ is a prime is close to machine precision.
\end{enumerate}
We hope to address some of these observations in the future.
\begin{table}
	\caption{Eigenvalue of $3 \times 3$ magic square}
	\small
	\begin{center}
	\begin{tabular}{|c|c|c|}
		\hline
		Eigenvalues of $M_n$ ($\mu$) & Eigenvalues of $J_nS_n$ ($\lambda$) & Relative error $\bkt{\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac13}$\\
		\hline
        -4.89897949 &      -5.196152422707 & 5.719096e-02\\ 
 \hline 
          4.89897949 &       5.196152422707 & 5.719096e-02\\ 
 \hline 
         15 &      15 & 0\\ 
 \hline 
	\end{tabular}
	\end{center}
	\label{magic_3}
\end{table}
\begin{table}
	\small
	\caption{Eigenvalue of $5 \times 5$ magic square}
	\begin{center}
	\begin{tabular}{|c|c|c|}
		\hline
		Eigenvalues of $M_n$ ($\mu$) & Eigenvalues of $J_nS_n$ ($\lambda$) & Relative error $\bkt{\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac15}$\\
		\hline
        -21.27676547 &     -21.266270208801 & 4.935168e-04\\ 
 \hline 
        -13.12628093 &     -13.143277802978 & 1.293199e-03\\ 
 \hline 
         13.12628093 &      13.143277802978 & 1.293199e-03\\ 
 \hline 
         21.27676547 &      21.266270208801 & 4.935168e-04\\ 
 \hline 
         65 &      65 & 0\\ 
 \hline 
	\end{tabular}
	\end{center}
	\label{magic_5}
\end{table}
\begin{table}
	\small
	\caption{Eigenvalue of $7 \times 7$ magic square}
	\begin{center}
	\begin{tabular}{|c|c|c|}
		\hline
		Eigenvalues of $M_n$ ($\mu$) & Eigenvalues of $J_nS_n$ ($\lambda$) & Relative error $\bkt{\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac17}$\\
		\hline
       -56.48482668 &     -56.466739338581 & 3.203186e-04\\ 
 \hline 
        -31.08815856 &     -31.336676188403 & 7.930568e-03\\ 
 \hline 
        -25.39666812 &     -25.130063150178 & 1.060901e-02\\ 
 \hline 
         25.39666812 &      25.130063150178 & 1.060901e-02\\ 
 \hline 
         31.08815856 &      31.336676188403 & 7.930568e-03\\ 
 \hline 
         56.48482668 &      56.466739338581 & 3.203186e-04\\ 
 \hline 
        175 &     175 & 0\\
 \hline 
	\end{tabular}
	\end{center}
	\label{magic_7}
\end{table}
\begin{table}
	\small
	\caption{Eigenvalue of $9 \times 9$ magic square}
	\begin{center}
	\begin{tabular}{|c|c|c|}
		\hline
		Eigenvalues of $M_n$ ($\mu$) & Eigenvalues of $J_nS_n$ ($\lambda$) & Relative error $\bkt{\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac19}$\\
		\hline
       -118.41407217 &    -118.414078206605 & 5.096328e-08\\ 
 \hline 
        -63.00687636 &     -63.006814987847 & 9.740991e-07\\ 
 \hline 
        -46.47580015 &     -46.765371804360 & 6.192010e-03\\ 
 \hline 
        -41.12470113 &     -41.124777781373 & 1.863973e-06\\ 
 \hline 
         41.12470113 &      41.124777781373 & 1.863973e-06\\ 
 \hline 
         46.47580015 &      46.765371804360 & 6.192010e-03\\ 
 \hline 
         63.00687636 &      63.006814987847 & 9.740991e-07\\ 
 \hline 
        118.41407217 &     118.414078206605 & 5.096329e-08\\ 
 \hline 
        369 &     369 & 0\\
 \hline 
	\end{tabular}
	\end{center}
	\label{magic_9}
\end{table}
\begin{table}
	\small
	\caption{Eigenvalue of $11 \times 11$ magic square}
	\begin{center}
	\begin{tabular}{|c|c|c|}
		\hline
		Eigenvalues of $M_n$ ($\mu$) & Eigenvalues of $J_nS_n$ ($\lambda$) & Relative error $\bkt{\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac1{11}}$\\
		\hline
       -214.74266474 &    -214.742664739496 & 1.005878e-14\\ 
 \hline 
       -111.90424039 &    -111.904240387785 & 6.442010e-12\\ 
 \hline 
        -80.05297263 &     -80.052972642020 & 1.428606e-10\\ 
 \hline 
        -66.51041338 &     -66.510413341849 & 5.782562e-10\\ 
 \hline 
        -61.12213518 &     -61.122135205551 & 4.611025e-10\\ 
 \hline 
         61.12213518 &      61.122135205551 & 4.611032e-10\\ 
 \hline 
         66.51041338 &      66.510413341849 & 5.782560e-10\\ 
 \hline 
         80.05297263 &      80.052972642020 & 1.428577e-10\\ 
 \hline 
        111.90424039 &     111.904240387786 & 6.441375e-12\\ 
 \hline 
        214.74266474 &     214.742664739495 & 9.529374e-15\\ 
 \hline 
        671 &     671 & 0\\ 
 \hline 
	\end{tabular}
	\end{center}
	\label{magic_11}
\end{table}
\begin{table}
	\small
	\caption{Eigenvalue of $13 \times 13$ magic square}
	\begin{center}
	\begin{tabular}{|c|c|c|}
		\hline
		Eigenvalues of $M_n$ ($\mu$) & Eigenvalues of $J_nS_n$ ($\lambda$) & Relative error $\bkt{\dfrac{\abs{\mu-\lambda}}{\abs{\lambda}} \leq \dfrac1{11}}$\\
		\hline
        -353.09013412 &    -353.090134118702 & 1.609884e-16\\ 
  \hline 
        -181.82867798 &    -181.828677981311 & 1.406793e-15\\ 
  \hline 
        -127.42740570 &    -127.427405699228 & 7.449615e-14\\ 
  \hline 
        -102.67515987 &    -102.675159868710 & 7.413024e-13\\ 
  \hline 
         -90.37276161 &     -90.372761608196 & 2.140605e-12\\ 
  \hline 
         -85.12062520 &     -85.120625195150 & 1.493364e-12\\ 
  \hline 
          85.12062520 &      85.120625195150 & 1.493531e-12\\ 
  \hline 
          90.37276161 &      90.372761608196 & 2.139190e-12\\ 
  \hline 
         102.67515987 &     102.675159868710 & 7.421328e-13\\ 
  \hline 
         127.42740570 &     127.427405699227 & 7.159660e-14\\ 
  \hline 
         181.82867798 &     181.828677981311 & 2.032035e-15\\ 
  \hline 
         353.09013412 &     353.090134118702 & 1.931861e-15\\ 
  \hline 
        1105 &    1105 & 0\\ 
  \hline
	\end{tabular}
	\end{center}
	\label{magic_13}
\end{table}
\begin{figure}[!htbp]
	\caption{Maximum relative error in eigenvalues between $M_n$ and $J_nS_n$ for odd $n$}
	\includegraphics[width=0.8\textwidth]{./figure.pdf}
	\label{fig_rel_error}
\end{figure}
\FloatBarrier
\section{Singly even magic squares}
In this section, we look at the eigenvalues of singly even magic square. The singly even magic squares in MATLAB are constructed using the Strachey method~\cite{ball2016mathematical}. The MATLAB code for obtaining singly even magic square is presented here~\cite{clevemagic2}. We also present it below for the sake of clarity and continuity.
\begin{tcolorbox}
\begin{verbatim}
		function [A] = magic(n)
	
		m = n/2; 
	
		A = magic (m);   
	
		A = [A, A+2*m*m; A+3*m*m, A+m*m];  
		k = (m-1)/2; 
		if (k > 1) 
			I = 1:m; 
			J = [2:k, n-k+2:n]; 
		
			A([I,I+m],J) = A([I+m,I],J);  
		end 
		I = [1:k, k+2:m];
		A([I,I+m],1) = A([I+m,I],1);
		I = k + 1;
		A([I,I+m],I) = A([I+m,I],I);
\end{verbatim}
\end{tcolorbox}

Let $M_n$ denote the singly even magic square of size $n \times n$, where $n=2m=4k+2$, where $n,m,k \in \mathbb{Z}^+$. The first two singly even magic squares generated by MATLAB are in Equations~\eqref{M_6} and~\eqref{M_10}.

\begin{align}
	M_6 & =
	\begin{bmatrix}
	    35 &    1 &    6 &   26 &    19 &    24\\
	     3 &   32 &    7 &   21 &    23 &    25\\
	    31 &    9 &    2 &   22 &    27 &    20\\
	     8 &   28 &   33 &   17 &    10 &    15\\
	    30 &    5 &   34 &   12 &    14 &    16\\
	     4 &   36 &   29 &   13 &    18 &    11
	\end{bmatrix}
	\label{M_6}
\end{align}
\begin{align}
	M_{10} & = \begin{bmatrix}
	92 &    99 &     1 &     8   &  15  &   67  &   74  &   51  &   58  &   40 \\
    98 &    80  &    7  &   14  &   16  &   73  &   55  &   57  &   64   &  41 \\
     4 &    81  &   88   &  20  &   22  &   54  &   56  &   63  &   70  &   47 \\
    85 &    87  &   19  &   21  &    3  &   60  &   62  &   69  &   71  &   28 \\
    86 &    93  &   25 &     2  &    9  &   61  &   68  &   75  &   52  &   34 \\
    17 &    24  &   76  &   83  &   90  &   42  &   49  &   26  &   33  &   65 \\
    23 &     5  &  82  &   89  &   91   &  48  &   30  &   32  &   39  &   66 \\
    79 &     6  &  13  &   95  &   97   &  29  &   31  &   38  &   45  &   72 \\
    10 &    12  &  94  &   96   &  78   &  35   &  37   &  44   &  46   &  53 \\
    11 &    18  &  100  &   77  &   84  &   36  &   43  &   50  &   27  &   59 
\end{bmatrix}
\label{M_10}
\end{align}

The singly even magic square thus obtained can be written as follows.
\begin{align}
	M_n & = \begin{bmatrix}
	M_m & M_m+2m^2\\
	M_m+3m^2 & M_m+m^2
	\end{bmatrix}\\
	& +
	\begin{bmatrix}
		\mathbb{I}_{m \times 1}\\
		-\mathbb{I}_{m \times 1}
	\end{bmatrix}
	\begin{bmatrix}
		3m^2 \mathbb{I}_{1 \times k} & \mathbb{O}_{1 \times \bkt{2k+3}} & -m^2 \mathbb{I}_{1 \times \bkt{k-1}}
	\end{bmatrix}
	-3m^2uv^T
\end{align}
where $u,v \in \mathbb{R}^{n \times 1}$ and
\begin{align}
	u_j = \begin{cases}
1 & j=k+1\\
-1 & j=n-k\\
0 & \text{ otherwise}
\end{cases},
v_j = \begin{cases}
1 & j=1\\
-1 & j=k+1\\
0 & \text{ otherwise}
\end{cases}
\end{align}



\begin{align}
	M_n & = \begin{bmatrix}M_m+\dfrac{3m^2}2\mathbb{I}_{m \times m} & M_m+\dfrac{3m^2}2\mathbb{I}_{m \times m}\\M_m+\dfrac{3m^2}2\mathbb{I}_{m \times m} & M_m+\dfrac{3m^2}2\mathbb{I}_{m \times m}\end{bmatrix} + \begin{bmatrix}-\dfrac{3m^2}2\mathbb{I}_{m \times m} & \dfrac{m^2}2\mathbb{I}_{m \times m}\\\dfrac{3m^2}2\mathbb{I}_{m \times m} & -\dfrac{m^2}2\mathbb{I}_{m \times m}\end{bmatrix}\\
	& + \begin{bmatrix}
		\mathbb{I}_{m \times 1}\\
		-\mathbb{I}_{m \times 1}
	\end{bmatrix}
	\begin{bmatrix}
		3m^2 \mathbb{I}_{1 \times k} & \mathbb{O}_{1 \times \bkt{2k+3}} & -m^2 \mathbb{I}_{1 \times \bkt{k-1}}
	\end{bmatrix}
	- 3m^2uv^T\\
	& = \begin{bmatrix}
	I_m\\
	I_m
	\end{bmatrix}
	\bkt{M_m+\dfrac{3m^2}2\mathbb{I}_{m \times m}}\begin{bmatrix}
	I_m&
	I_m
	\end{bmatrix}
	+\dfrac{m^2}2\begin{bmatrix}
		\mathbb{I}_{m \times 1}\\
		-\mathbb{I}_{m \times 1}
	\end{bmatrix}
	\begin{bmatrix}
		-3\mathbb{I}_{1 \times m} & \mathbb{I}_{1 \times m}
	\end{bmatrix}\\
	& + \dfrac{m^2}2\begin{bmatrix}
		\mathbb{I}_{m \times 1}\\
		-\mathbb{I}_{m \times 1}
	\end{bmatrix}
	\begin{bmatrix}
		6 \mathbb{I}_{1 \times k} & \mathbb{O}_{1 \times \bkt{2k+3}} & -2 \mathbb{I}_{1 \times \bkt{k-1}}
	\end{bmatrix} - 3m^2uv^T
\end{align}

\begin{align}
	M_n & = \begin{bmatrix}
	I_m\\
	I_m
	\end{bmatrix}
	\bkt{M_m+\dfrac{3m^2}2\mathbb{I}_{m \times m}}\begin{bmatrix}
	I_m&
	I_m
	\end{bmatrix} - 3m^2uv^T\\
	& +\dfrac{m^2}2\begin{bmatrix}
		\mathbb{I}_{m \times 1}\\
		-\mathbb{I}_{m \times 1}
	\end{bmatrix}
	\begin{bmatrix}
		3 \mathbb{I}_{1 \times k} & -3 \mathbb{I}_{1 \times \bkt{k+1}} & \mathbb{I}_{1 \times \bkt{k+2}} & - \mathbb{I}_{1 \times \bkt{k-1}}
	\end{bmatrix}




\end{align}
Hence, we see that rank of $M_n$ is at most $m+2$ or equivalently at least $m-2$ of the eigenvalues are zeros. We will now prove the following result on the eigenvalues of $M_n$.

\begin{tcolorbox}[title=Non-trivial eigenvalues of singly even magic squares]
\begin{theorem}
	\label{main_singly_even_theorem}
	Let $M_n$ be a singly even magic square, where $n=2m=4k+2$, generated using the \textbf{Strachey method}. Let $$\Lambda = \left\{\lambda_j = \dfrac{n^2}{4\sin\bkt{j\pi/m}}: j \in\{\pm1,\pm2,\ldots,\pm k\} \right\}$$
	\begin{enumerate}
		\item
		Two of the non-trivial eigenpairs of $M_n$ are $\bkt{-3m^2,\begin{bmatrix}\mathbb{I}_m\\-\mathbb{I}_m\end{bmatrix}}$ and $\bkt{3m^2,\begin{bmatrix}\mathbb{I}_m - 3e_{k+1}\\-\bkt{\mathbb{I}_m - 3e_{k+1}}\end{bmatrix}}$.
		\item
		$2k-1$ of the non-trivial eigenvalues are zero.
		\item
		The set of remaining $2k$ eigenvalues are $\frac2n$-approximated by the set
		$$\left\{\pm\dfrac{n^2}{4\sin\bkt{\pi/m}},\pm\dfrac{n^2}{4\sin\bkt{2\pi/m}},\ldots,\pm\dfrac{n^2}{4\sin\bkt{k\pi/m}}\right\}$$
	\end{enumerate}
\end{theorem}
\end{tcolorbox}


\begin{proof}
Consider the orthonormal matrix $Q_n = \dfrac1{\sqrt{2}}\begin{bmatrix}I_m & -I_m\\I_m & I_m\end{bmatrix}$. The eigenvalues of $M_n$ and $Q_n M_nQ_n^T$ are the same. Let us now evaluate $Q_nM_nQ_n^T$. Note that
\begin{align}
	Q_n\begin{bmatrix}I_m\\I_m\end{bmatrix} = \dfrac1{\sqrt2}\begin{bmatrix}\mathbb{O}_{m \times m}\\ 2I_m \end{bmatrix}
\end{align}
\begin{align}
	Q_n\begin{bmatrix}\mathbb{I}_{m\times 1}\\-\mathbb{I}_{m\times 1}\end{bmatrix} = \dfrac1{\sqrt2}\begin{bmatrix}2\mathbb{I}_{m \times 1} & \mathbb{O}_{m \times 1}\end{bmatrix}
\end{align}
\begin{align}
	Q_nu = \sqrt2 e_{k+1}
\end{align}
\begin{align}
	Q_nv = \dfrac1{\sqrt2} \tilde{v}
\end{align}
where $\tilde{v} \in \mathbb{R}^{n \times 1}$ with $\tilde{v}_j = \begin{cases}
1 & j=1,m+1\\
-1 & j=k+1,m+k+1\\
0 & \text{ otherwise}
\end{cases}$
\begin{align}
	Q_n\begin{bmatrix}
		\mathbb{I}_{m \times 1}\\
		-\mathbb{I}_{m \times 1}
	\end{bmatrix} & = \sqrt2 \begin{bmatrix}
	\mathbb{I}_{m\times 1}\\
	\mathbb{O}_{m\times 1}
	\end{bmatrix}
\end{align}
\begin{align}
	Q_n\begin{bmatrix}
		3 \mathbb{I}_{k \times 1} \\ -3 \mathbb{I}_{\bkt{k+1} \times 1} \\ \mathbb{I}_{\bkt{k+2} \times 1} \\ - \mathbb{I}_{\bkt{k-1} \times 1}
	\end{bmatrix} & =\sqrt2
	\begin{bmatrix}
		\mathbb{I}_k\\
		-2\\
		-2\\
		-\mathbb{I}_{k-1}\\
		2\mathbb{I}_k\\
		-1\\
		-1\\
		-2\mathbb{I}_{k-1}
	\end{bmatrix}
\end{align}

Hence,
























	\begin{align}
		\boxed{Q_n M_n Q_n^T = \begin{bmatrix}
	m^2A_{m} & m^2B_{m}\\
	\mathbb{O}_{m \times m} & 2M_m + 3m^2\mathbb{I}_{m \times m}
	\end{bmatrix}
	}
	\end{align}
where 
\begin{align}
	A_m = \begin{bmatrix}\mathbb{I}_{m \times k} & -2\mathbb{I}_{m \times 1} & -2\mathbb{I}_{m \times 1} & -\mathbb{I}_{m \times \bkt{k-1}}\end{bmatrix}-3e_{k+1}w^T
\end{align}
and
\begin{align}
	B_m=\begin{bmatrix}2\mathbb{I}_{m \times k} & - \mathbb{I}_{m \times 1} & - \mathbb{I}_{m \times 1} & -2\mathbb{I}_{m \times \bkt{k-1}}\end{bmatrix}-3e_{k+1}w^T
\end{align}
where $w,e_{k+1} \in \mathbb{R}^{m \times 1}$ with
\begin{align}
	w_j = \begin{cases}
1 & j=1\\
-1 & j=k+1\\
0 & \text{ otherwise}
\end{cases}
\end{align} and $e_{k+1}$ is a unit vector whose $\bkt{k+1}^{th}$ entry is $1$. Hence, the eigenvalues of $M_n$ are the eigenvalues of $Q_nM_nQ_n^T$, which in-turn are the eigenvalues of $m^2A_m$ and $2M_m + 3m^2\mathbb{I}_{m \times m}$.

Note that $A_m$ is a rank $2$ matrix, since
\begin{align}
	A_m = \mathbb{I}_{m \times 1}\begin{bmatrix}\mathbb{I}_{1 \times k} & -2 & -2 & -\mathbb{I}_{1 \times \bkt{k-1}}\end{bmatrix} -3e_{k+1}w^T
\end{align}

Further, note that
\begin{align}
	A_m \mathbb{I}_{m \times 1} & = \mathbb{I}_{m \times 1}\begin{bmatrix}\mathbb{I}_{1 \times k} & -2 & -2 & -\mathbb{I}_{1 \times \bkt{k-1}}\end{bmatrix} \mathbb{I}_{m \times 1} -3e_{k+1}w^T \mathbb{I}_{m \times 1}\\
	& = -3\mathbb{I}_{m \times 1}
\end{align}
since $w^T \mathbb{I} = 0$. Also, note that the trace of $A_m$ is zero. Hence, the non-zero eigenvalues of $A_m$ are $\pm3$. Hence, two of the non-zero eigenvalues of $M_n$ are $\pm3m^2$, i.e., $\pm \dfrac{3n^2}4$. The eigenvector corresponding to the eigenvalue $3$ is $\begin{bmatrix}\mathbb{I}_m - 3e_{k+1}\\-\bkt{\mathbb{I}_m - 3e_{k+1}}\end{bmatrix}$. The remaining non-zero eigenvalues of $M_n$ are the eigenvalues of $2M_m+3m^2\mathbb{I}_{m \times m}$. From the previous section, we have that the eigenvalues of $M_m$ to be its magic sum, $\dfrac{m\bkt{m^2+1}}2$ and the remaining are $\frac1m$-approximated by
\begin{align}
\left\{\pm\dfrac{m^2}{2\sin\bkt{j\pi/m}} \text{ where }j\in\{1,2,\ldots,k\} \right\}
\end{align}
Of course the magic sum of $M_n = \dfrac{n\bkt{n^2+1}}2$ is nothing but $2\bkt{\dfrac{m\bkt{m^2+1}}2} + 3m^3$. The remaining non-trivial non-zero eigenvalues of $M_n$ are $\frac2n$-approximated by
\begin{align}
	\left\{\pm \dfrac{n^2}{4\sin\bkt{\pi/m}},\pm \dfrac{n^2}{4\sin\bkt{2\pi/m}},\ldots,\pm \dfrac{n^2}{4\sin\bkt{k\pi/m}} \right\}
\end{align}
\end{proof}

\noindent \begin{remark} \label{remark3} If we could prove the stronger $\frac1n$-well-approximated remark as in Remark~\ref{remark1}, then Theorem~\ref{main_singly_even_theorem} could also be further strengthened appropriately.
\end{remark}
\section{Doubly even magic squares}
The eigenvalues of doubly even magic squares generated by MATLAB have already been obtained in~\cite{kirkland1995group}. However, for the sake of completeness we rederive these as part of this article. The code for generating doubly even magic squares is presented here~\cite{clevemagic2}. We present the code below for the sake of clarity and continuity.
\begin{tcolorbox}
\begin{verbatim}
	    M = reshape(1:n^2,n,n)';
	    [I,J] = ndgrid(1:n);
	    K = fix(mod(I,4)/2) == fix(mod(J,4)/2);
	    M(K) = n^2+1 - M(K);
\end{verbatim}
\end{tcolorbox}
The first two doubly even magic squares generated by MATLAB are in Equations~\eqref{M_4} and~\eqref{M_8}.
\begin{align}
	\label{M_4}
	M_4 & = \begin{bmatrix}
	16 & 2 & 3 & 13\\
     5 & 11 & 10 & 8\\
     9 & 7 & 6 & 12\\
     4 & 14 & 15 & 1
	 \end{bmatrix}\\
	 	\label{M_8}
	 M_8 & =\begin{bmatrix}
     64  &   2   &   3  &    61 &    60 &    6  &    7  &    57\\
      9  &   55  &   54 &    12 &    13 &    51 &    50 &   16\\
     17  &   47  &   46 &    20 &    21 &    43 &    42 &    24\\
     40  &   26  &   27 &    37 &    36 &    30 &    31 &    33\\
     32  &   34  &   35 &    29 &    28 &    38 &    39 &    25\\
     41  &   23  &   22 &    44 &    45 &    19 &    18 &    48\\
     49  &   15  &   14 &    52 &    53 &    11 &    10 &    56\\
      8  &   58  &   59 &     5 &     4 &    62 &    63 &     1
	 \end{bmatrix}
\end{align}
The main theorem proved in this section is below (Theorem~\ref{main_doubly_even_theorem}).
\begin{tcolorbox}[title=Non-trivial eigenvalues of doubly even magic squares]
\begin{theorem}
	\label{main_doubly_even_theorem}
	The rank of the doubly even magic square generated by MATLAB's \textbf{criss-cross} method is $3$ and the non-zero non-trivial eigenvalues are
	$$\pm \dfrac{n}2 \sqrt{\dfrac{n^3-n}3}$$
	The corresponding eigenvectors are $\sqrt{n^2-1}v \pm \sqrt{3n}u$, where
\setcounter{MaxMatrixCols}{12}
	$$u = \begin{bmatrix} w\\
	-\text{flip}\bkt{w} \end{bmatrix}_{4m \times 1}$$ and
	$$v^T = \begin{bmatrix} 1 & -1 & -1 & 1 & \cdots & 1 & -1 & -1 & 1 \end{bmatrix}_{1 \times 4m}$$
	with $w^T = \begin{bmatrix}(n-1) & -(n-3) & -(n-5) & (n-7) & \cdots & (-1)^m\end{bmatrix}$.
\end{theorem}
\end{tcolorbox}
The proof of the above theorem (Theorem~\ref{main_doubly_even_theorem}) constitutes this section.

Let $P_4 = I_4+J_4$, where $I_4 \in \mathbb{R}^{4 \times 4}$ is the identity matrix and $J_4 \in \mathbb{R}^{4 \times 4}$ is the reverse identity matrix. If $S_n$ denotes the $n \times n$ with numbers arranged row wise from $1$ to $n^2$, i.e., $S_n\bkt{i,j} = j+(i-1)n$, it is fairly easy to show that
\begin{align}
	M_n\bkt{i,j} & = S_n\bkt{i,j} \bkt{1-P_4\bkt{i\bmod4,j\bmod4}}\\
	 & + S_n\bkt{n+1-i,n+1-j} P_4\bkt{i\bmod4,j\bmod4}
\end{align}
It is a fairly elementary algebra to show that 
\begin{align}
	M_n = \dfrac{n^2+1}2 \mathbb{I} \mathbb{I}^T + \dfrac{n}2 uv^T + \dfrac{vu^T}2
\end{align}
where
\begin{align}u = \begin{bmatrix} w\\
	-\text{flip}\bkt{w} \end{bmatrix}_{4m \times 1}
\end{align}
where $w^T = \begin{bmatrix}(n-1) & -(n-3) & -(n-5) & (n-7) & \cdots & (-1)^m\end{bmatrix}$ and
\begin{align}
	v^T = \begin{bmatrix} 1 & -1 & -1 & 1 & \cdots & 1 & -1 & -1 & 1 \end{bmatrix}_{1 \times 4m}
\end{align}
This gives us
\begin{align}
	M_n & = \begin{bmatrix}
		\mathbb{I} & u & v
	\end{bmatrix}
	\begin{bmatrix}
		\dfrac{n^2+1}2 & 0 & 0\\
		0 & 0 & \dfrac{n}2\\
		0 & \dfrac12 & 0
	\end{bmatrix}
	\begin{bmatrix}
		\mathbb{I}^T \\ u^T \\ v^T
	\end{bmatrix}\\
	& = USU^T
\end{align}
where $U =  \begin{bmatrix} \mathbb{I} & u & v \end{bmatrix}$ and $S = \dfrac12 \begin{bmatrix}n^2+1 & 0 & 0\\ 0 & 0 & n\\ 0 & 1 & 0\end{bmatrix}$.

We have
\begin{align}
	D & = U^TU = \dfrac{n}3 \begin{bmatrix}
3 & 0 & 0\\
0 & n^2-1 & 0\\
0 & 0 & 3
\end{bmatrix}
\end{align}
Note that $M_n=USU^T$ is almost the SVD of the matrix $M_n$ except for ordering of the vectors and scaling. The singular values of $M_n$ are
\begin{align}
	\dfrac{n\bkt{n^2+1}}2, \dfrac{n^2}2\sqrt{\dfrac{n^2-1}3} \text{ and }\dfrac{n}2\sqrt{\dfrac{n^2-1}3}
\end{align} 
Since $M_n$ is rank $3$, we immediately obtain that $n-3$ eigenvalues are zeros and the corresponding eigenvectors are the vectors orthogonal to $\mathbb{I}$, $u$ and $v$. Now let's obtain the non-zero eigenvalues of $M_n$. Let $\bkt{\lambda,x}$ be a non-zero non-trivial eigenpair of $M_n$. We then need $M_nx = \lambda x$. This immediately gives us that $x=Uy$, where $y \in \mathbb{R}^{3 \times 1}$. Hence, we need
\begin{align}
	USU^TUy & = \lambda U y\\
	SDy = \lambda y
\end{align}
Hence, $\bkt{\lambda,y}$ is an eigenpair of $SD \in \mathbb{R}^{3 \times 3}$. We have
\begin{align}
	SD & = \dfrac{n}{6}\begin{bmatrix}
	3\bkt{n^2+1} & 0 & 0\\
	0 & 0 & 3n\\
	0 & n^2-1 & 0
	\end{bmatrix}
\end{align}
We hence obtain the eigenpairs corresponding to non-zero eigenvalues as
\begin{enumerate}
	\item
	$\bkt{\dfrac{n\bkt{n^2+1}}2,\mathbb{I}_n}$
	\item
	$\bkt{\dfrac{n}2\sqrt{\dfrac{n^3-n}3},\sqrt{3n}u+\sqrt{n^2-1}v}$
	\item
	$-\bkt{\dfrac{n}2\sqrt{\dfrac{n^3-n}3},-\sqrt{3n}u+\sqrt{n^2-1}v}$
\end{enumerate}
This gives us the non-zero eigenvalues to be $\dfrac{n\bkt{n^2+1}}2, \pm \dfrac{n}2\sqrt{\dfrac{n^3-n}3}$ and thereby proving Theorem~\ref{main_doubly_even_theorem}.
\section{Conclusion}
Eigenvalues of the magic squares obtained using MATLAB's magic() function are investigated.  We obtain approximations of eigenvalues of odd and singly even magic squares with error bounds. Further, the eigenvalues of doubly even magic squares are exactly obtained. The approximation of the spectra involves some interesting connections with the spectrum of g-circulant matrices and the use of Bauer-Fike theorem. We also highlight some observations which we hope to prove in future.
\section{Appendix}
Throughout the appendix, we will adopt the following.
\begin{enumerate}
    \item $n$ is an odd number and will be written as $n=2m+1$
	\item
	$g \in \{1,2,\ldots,n-1\}$ is a number relatively prime to $n$, i.e., $\gcd\bkt{g,n}=1$
	\item
	$C \in \mathbb{C}^{n\times n}$ is a circulant matrix
\end{enumerate}
\begin{lemma}
	\label{lem_magic_sum}
	For any $n \times n$ magic square, $M_n$, we have $M_n \mathbb{I}_{n \times 1} = \dfrac{n\bkt{n^2+1}}2 \mathbb{I}_{n \times 1}$ and $\mathbb{I}_{1 \times n} M_n = \dfrac{n\bkt{n^2+1}}2 \mathbb{I}_{1 \times n}$.
\end{lemma}
\begin{proof}
	Follows from the fact that row sum and column sum of a magic square is the magic sum.
\end{proof}
\begin{lemma}
	\label{lem_magic_sum_1}
	The magic sum is the largest eigenvalue with multiplicity $1$.
\end{lemma}
\begin{proof}
	Follows immediately from Perron-Frobenius theorem.
\end{proof}
\begin{lemma}
	\label{lem_ortho_eig}
	Every eigenvector of $M_n$ other than $\mathbb{I}_{n \times 1}$ is orthogonal to $\mathbb{I}_{n \times 1}$.
\end{lemma}
\begin{proof}
	Let $\bkt{\lambda,v}$ be an eigenpair of $M_n$ other than $\bkt{\dfrac{n\bkt{n^2+1}}2,\mathbb{I}_{n \times 1}}$. We then have
	$$M_nv = \lambda v$$
	Premultiplying by $\mathbb{I}_{1 \times n}$, we obtain
	$$\mathbb{I}_{1\times n}M_nv = \lambda \mathbb{I}_{1\times n} v \implies \dfrac{n\bkt{n^2+1}}2 \mathbb{I}_{1\times n} v = \lambda \mathbb{I}_{1\times n} v$$
	Since $\lambda \neq \dfrac{n\bkt{n^2+1}}2$, we have $\mathbb{I}_{1 \times n}v = 0$.
\end{proof}
\begin{note}
	\label{note1}
	\begin{align}
		F_n\bkt{p,q} & = \omega^{\bkt{p-1}\bkt{q-1}}
	\end{align}
\end{note}
\begin{note}
	\label{note2}
	\begin{align}
		v_i^*v_j = n\delta_{ij}
	\end{align}
\end{note}
\begin{note}
	\label{note3}
	\begin{align}
		F_n^*F_n = nI_{n} = F_nF_n^*
	\end{align}
\end{note}
\begin{note}
	\label{note4}
	The vectors $v_j$ are the eigenvectors of the circulant matrix $C$. The eigenvalue corresponding to the vector $v_j$ is given by
	\begin{align}
		\lambda_j = \dfrac{v_j^*Cv_j}{v_j^*v_j} = \dfrac{v_j^*Cv_j}n
	\end{align}
\end{note}
\begin{note}
	\label{note5}
	If $a$ is the first row of the circulant matrix $C$, we then have
	\begin{align}
		\lambda_j = av_j
	\end{align}
\end{note}
\begin{note}
	\label{note6}
	If $\Lambda_n \in \mathbb{C}^{n \times n}$ is the matrix of eigenvalues of $C$, i.e.,
	\begin{align}
		\Lambda_n\bkt{i,j} = \begin{cases}
	\lambda_j & \text{ if }i=j\\
	0 & \text{ otherwise}
	\end{cases}
	\end{align}
	we then have
	\begin{align}
		\Lambda_n = \dfrac{F_n^*CF_n}n\\
		C = \dfrac{F_n\Lambda_n F_n^*}n
	\end{align}
\end{note}
\begin{lemma}
\label{lem1}
Let $X_n \in \mathbb{R}^{n \times n}$ be a circulant matrix whose first row is
\begin{align}
    \bf{a} & =\begin{bmatrix}
m & m+1 & m+2 & \cdots & 2m & 0 & 1 & \cdots & m-1
    \end{bmatrix}
\end{align}
We then have the eigenvalues of $X_n$ to be given by
\begin{align}
    \lambda_j & = \begin{cases}
    \dfrac{n\bkt{n-1}}2 & \text{ if }j=0\\
    \dfrac{\bkt{-1}^{j} in}2\csc\bkt{\dfrac{\pi j}n} & \text{ if }j \in \{1,2,\ldots,n-1\}
    \end{cases}
\end{align}
\end{lemma}
\begin{proof}
The $j^{th}$ eigenvalue is given by
\begin{align}
    \lambda_j & = av_j
\end{align}
Hence, we have
\begin{align}
    \lambda_j & = \displaystyle\sum_{k=0}^{2m} \bkt{\bkt{m+k}\mod n}\omega^{kj} = \omega^{-mj}\bkt{\displaystyle\sum_{k=0}^{2m} \bkt{\bkt{m+k}\mod n}\omega^{(m+k)j}}\\
    & = \omega^{-mj}\bkt{\displaystyle\sum_{k=0}^{2m} k\omega^{kj}}
\end{align}
For $j=0$, we have $\lambda_n = 0+1+\cdots+2m = \dfrac{n\bkt{n-1}}2$. Let $j \in \{1,2,\ldots,n-1\}$. We then have
\begin{align}
    \lambda_j & =\omega^{-mj}\bkt{\displaystyle\sum_{k=0}^{2m} k\omega^{kj}} = \dfrac{n\omega^{-mj}}{\omega^{2j}-1} = \dfrac{n\exp\bkt{\dfrac{2\pi imj}n}}{\exp\bkt{-\dfrac{2\pi i j}n}-1} = \dfrac{n\exp\bkt{\dfrac{2\pi ij(n-1)}{2n}}}{\exp\bkt{-\dfrac{2\pi i j}n}-1}\\
    & = \dfrac{(-1)^{j+1} n}{2i\sin\bkt{\dfrac{\pi j}n}} = \bkt{-1}^{j} \dfrac{i n}{2\sin\bkt{\dfrac{\pi j}n}} = \boxed{\dfrac{(-1)^j i n}2 \csc\bkt{\dfrac{\pi j}n}}
\end{align}
\end{proof}
\begin{note}
	\label{note7}
	Let $R_g$ be a $g$-row circulant matrix, where $g\in\{1,2,\ldots,n-1\}$, if for all $i,j,k \in \{1,2,\ldots,n\}$, we have
	\begin{align}
		R_g\bkt{i,j} = R_g\bkt{\bkt{i+k}\bmod{n},\bkt{j+gk}\bmod{n}}
	\end{align}
	Note that $R_1$ is the usual circulant matrix.
\end{note}
\begin{note}
\label{note9}
	$R_g = P_g R_1$, where $P_g$ is a $g$-row circulant matrix with first row being $e_1^T$, $R_g$ is a $g$-row circulant matrix and $R_1$ is a circulant matrix having the same first row as $R_1$.
\end{note}
\begin{note}{\textbf{Spectrum of reverse circulant matrix}}
\label{note12}
	Let $C$ be the circulant matrix with eigenpairs given by $\bkt{\lambda_j,v_j}_{j=1}^n$, i.e.,
	\begin{align}
		Cv_j & = \lambda_j v_j
	\end{align}
	We now need the eigenpairs of the matrix $JC$, where $J$ is the flipped identity matrix. From above, we have
	\begin{align}
		JCv_j & = \lambda_j Jv_j
	\end{align}
	We have
	\begin{align}
		Jv_j & = \begin{bmatrix}\omega^{(n-1)j} & \omega^{(n-2)j} & \cdots & \omega^{2j} & \omega^{j} & 1\end{bmatrix}^T\\
		& = \omega^{-j} \begin{bmatrix}1 & \omega^{-j} & \cdots & \omega^{-(n-2)j} & \omega^{-(n-1)j} \end{bmatrix}^T\\
		& = \omega^{-j} v_{n-j}
	\end{align}
	
	Note that $v_{n-j} = \overline{v_j}$ and $\lambda_{n-j} = \overline{\lambda_j}$. Hence, we obtain
	\begin{align}
		JCv_j & = \lambda_j Jv_j = \omega^{-j} \lambda_j v_{n-j}
	\end{align}
	Replacing $j$ by $n-j$, we obtain
	\begin{align}
		JCv_{n-j} & = \lambda_{n-j} Jv_{n-j} = \omega^{-\bkt{n-j}} \lambda_{n-j} v_{j} = \omega^{j} \lambda_{n-j} v_{j} = \overline{\omega^{-j} \lambda_{j}} v_{j}
	\end{align}
	Hence, we see that
	\begin{align}
		JC \bkt{v_j + \dfrac{\omega^{-j} \lambda_j}{\abs{\omega^{-j} \lambda_j}} v_{n-j}} & = \omega^{-j}\lambda_j v_{n-j} + \abs{\omega^{-j} \lambda_j} v_j\\
		& = \abs{\omega^{-j} \lambda_j} \bkt{v_j + \dfrac{\omega^{-j} \lambda_j}{\abs{\omega^{-j} \lambda_j}} v_{n-j}}
	\end{align}
	Hence, we obtain the eigenpairs of $JC$ to be
	$$\bkt{\abs{\lambda_j},v_j + \dfrac{\omega^{-j} \lambda_j}{\abs{\lambda_j}} v_{n-j}}$$
	where $\bkt{\lambda_j,v_j}$ are the eigenpairs of $C$. Note that if we let $u_j = \dfrac{v_j + \dfrac{\omega^{-j} \lambda_j}{\abs{\lambda_j}} v_{n-j}}{\magn{v_j + \dfrac{\omega^{-j} \lambda_j}{\abs{\lambda_j}} v_{n-j}}_2}$, we have
	\begin{align}
		u_j^*u_k & = 1
	\end{align}
\end{note}







2022-06-01T02:10:44.000
2205
2205.15563
en
https://arxiv.org/abs/2205.15563
This article looks at the eigenvalues of magic squares generated by the MATLAB's magic($n$) function. The magic($n$) function constructs doubly even ($n = 4k$) magic squares, singly even ($n = 4k+2$) magic squares and odd ($n = 2k+1$) magic squares using different algorithms. The doubly even magic squares are constructed by a criss-cross method that involves reflecting the entries of a simple square about the center. The odd magic squares are constructed using the Siamese method. The singly even magic squares are constructed using a lower-order odd magic square (Strachey method). We obtain approximations of eigenvalues of odd and singly even magic squares and prove error bounds on the approximation. For the sake of completeness, we also obtain the eigenpairs of doubly even magic squares generated by MATLAB. The approximation of the spectra involves some interesting connections with the spectrum of g-circulant matrices and the use of Bauer-Fike theorem.
General Mathematics (math.GM)
Spectrum of MATLABs magic squares
Qwen/Qwen-72B
1. YES
2. YES