text:
\section{Introduction}
\label{sec:intro}

Though they were first discovered 50 years ago \citep{Hewish:1968_709}, neutron stars are still not well understood compared to other more common stellar objects like the sun.
A major challenge to our understanding of neutron stars is identifying the properties of the poorly understood high dense nuclear matter that exists in their cores. This uncertainty is reflected
by the fact that a large number of nuclear-matter equation of state (EOS) models with various different predicted properties of neutron stars have been proposed in the past 50 years or so.
The observed properties of neutron stars, such as their masses and radii, can thus be used to
put constraints on the theoretical EOS models \citep[see, e.g.,][]{Lattimer:2012_485}.

As the structure and properties (both static and dynamical) of a neutron star in general depend sensitively on the underlying EOS model, it is thus quite interesting that various approximately
EOS-insensitive universal relations for neutron stars have indeed been found. These relations generally connect different physical quantities of neutron stars, and they are said to be universal in the sense that they are insensitive to EOS models to the $O(1\%)$ level (see \citet{Yagi:2017_1} for a review).

For instance, \citet{Lau:2010_1234} found a pair of universal relations to connect the frequency and damping rate of the quadrupolar $f$-mode to the mass $M$ and moment of inertia $I$ of nonrotating
neutron stars. It is also shown in \citet{Lau:2010_1234} that the values of $M$, $I$, and the stellar radius $R$ of a neutron star can be inferred accurately from the $f$-mode gravitational wave signals (if detected) emitted from the star.
More recently, the so-called I-Love-Q relations discovered by \citet{Yagi:2013_023009,Yagi:2013_365} have gained a lot of interest. These relations connect the moment of inertia $I$, the tidal deformability (also called the Love number), and the spin-induced quadrupole moment $Q$ of slowly rotating neutron stars.
In particular, the {\it I}$-${\it Q} relation has been extended to include rapid rotation.
It was initially found by \citet{Doneva:2014_L6} that the {\it I}$-${\it Q} relation is broken and becomes more EOS-dependent when considering rapidly rotating stars with a fixed rotation frequency $f$.
However, it was then found by \citet{Pappas:2014_121101} that the {\it I}$-${\it Q} relation remains approximately EOS-insensitive if, instead of the dimensional quantity $f$, one uses the dimensionless spin parameter $j$ to characterize rotation. This conclusion is extended by \cite{Chakrabarti:2014_201102} to include other dimensionless parameters such as $Mf$ or $Rf$.
In effect, there still exists a universal {\it I}$-${\it Q} relation for each value of $j$ or $Mf$ for rapidly rotating neutron stars.

It is noted that previously known universal relations connect the intrinsic stellar properties of a neutron star itself (e.g., its mass) and the dynamical response of the star (e.g., its tidal deformability) to external perturbations.
In this work, we propose a pair of universal relations that connect the radius and orbital frequency of the innermost stable circular orbit (ISCO) to the mass and spin frequency of
rapidly rotating neutron stars.
The ISCO is an important prediction of general relativity concerning the strong field spacetime around a compact stellar object\footnote{It should be mentioned that an ISCO can also appear
in Newtonian gravity around highly oblate objects, such as rapidly rotating low-mass strange stars
\citep{Zdunik:2001_087501,Kluzniak:2013_2825}.}.
It is linked to the geodesic motion of a test particle that orbits close to the compact object.
As we shall discuss below, the ISCO universal relations are practical from an observational perspective since the ISCO may be closely related to the kilohertz quasi-periodic oscillations (kHz QPOs) well observed from neutron stars in low-mass X-ray binaries (LMXBs).

The plan of this paper is as follows. Section~\ref{sec:method} presents our numerical methods and
chosen EOS models used in this work. The ISCO around rotating compact objects is also reviewed
briefly.
Our main numerical results and the proposed ISCO universal relations are discussed in Section~\ref{sec:results}. In Section~\ref{sec:analy}, we compare the relation for the ISCO
radius to the corresponding relation for the Kerr black hole and also study its connection to
other known universal relations for rotating neutron stars.
Finally, Section~\ref{sec:discuss} discusses
the astrophysical relevance of the ISCO universal relations. Unless otherwise noted, we use geometric units where $G=c=1$.

\section{Rotating neutron stars and ISCO in general relativity}
\label{sec:method}

The computation of rotating stellar models in general relativity is a nontrivial task. Nevertheless, a few public codes using different formulations and numerical methods to construct rapidly rotating neutron stars in general relativity are readily available (see \citet{Paschalidis:2017_7} for a review). In this work, we use the numerical code {\tt rotstar} from the publicly available C++ LORENE library\footnote{http://www.lorene.obspm.fr} that solves the
Einstein equations in a stationary and axisymmetric spacetime assuming a perfect-fluid matter source
using a multidomain spectral method \citep{Bonazzola:1993_421,Bonazzola:1998_104020}.

One needs to provide an EOS for nuclear matter to compute neutron star models. While the EOS
in the high-density core of neutron stars is still not well understood, the observations
of neutron stars with masses $M \approx 2 M_\odot$ \citep{Demorest:2010_1081,Antoniadis:2013_6131} have already ruled out many soft EOS models in general relativity.
Note, however, that some EOS models that are ruled out by this $2 M_\odot$ constraint would be revived in other theories of gravity \citep[see, e.g.,][]{Pani:2012_084020,Sham:2012_064015}.
In this study, we assume general relativity is the correct theory of gravity and employ 12 different nuclear-matter EOS models based on various theoretical approaches, including nuclear many-body theory, nuclear energy density functional theory, and Skyrme mean-field models.
The first nine EOSs are APR \citep{Akmal:1998_1804}, AU (the AV14+UVII model in \citet{Wiringa:1988_1010}), BSk20, BSk21 \citep{Potekhin:2013_A48}, GM1 \citep{Glendenning:1991_2414}, SKa, SKI2 \citep{Gulminelli:2015_055803},
SLy4 \citep{Douchin:2000_107}, and UU (the UV14+UVII model in \citet{Wiringa:1988_1010}). The remaining three EOSs are the three representative models (soft, intermediate, and stiff) presented in \citet{Hebeler:2013_11}, which are based on nuclear interactions derived from chiral
effective field theory combined with observational constraints. All of our chosen EOS models can support a nonrotating neutron star with a maximum mass larger than $2 M_\odot$ in order to be consistent with the $2 M_\odot$ observational constraint \citep{Demorest:2010_1081,Antoniadis:2013_6131}.

In this work, we are interested in the ISCO of test particles in a prograde orbit on the equatorial plane around rotating neutron stars
\citep{Miller:1998_793,Zdunik:2000_612,Bhattacharyya:2011_3247,Torok:2014_L5,Cipolletta:2017_024046}.
The spacetime outside a nonrotating neutron star is described by the Schwarzschild metric and it is well known that the circumferential radius of the ISCO is given by $R_{\rm ISCO} = 6 M$ in this case. In contrast to coordinate radius, the circumferential radius is more meaningful physically and hence we shall adopt this as our definitions for both the ISCO and stellar radii in this work.
While the ISCO radius for a nonrotating star depends only on $M$, the situation for rotating stars
is more complicated. To first order in the dimensionless spin parameter $j\equiv J/M^2$, where $J$ is the angular momentum, the spacetime of a slowly rotating star is determined uniquely by $M$ and $j$ and the ISCO radius is given by \citep[e.g.,][]{Miller:1998_793}
\begin{equation}
R_{\rm ISCO} = 6 M  \left[ 1 - j \left({2\over 3} \right)^{3/2} \right] .
\label{eq:Risco_slow}
\end{equation}
As we shall see below, the exact numerical results of $R_{\rm ISCO}$ for rotating neutron stars
start to deviate from the prediction of Equation~(\ref{eq:Risco_slow}) when $j\gtrsim 0.1$.
For a rapidly rotating neutron star, nonspherical deformation of the star is significant and
the contributions from higher-order multipole moments of the star to the exterior vacuum spacetime
become more important (see Section~\ref{sec:other_relations}).

For comparison, the ISCO radius around a Kerr black hole depends only on the two parameters $M$ and $j$, and it can be obtained analytically. The Boyer$-$Lindquist radial coordinate $r_{+}$ of the ISCO for a prograde orbit around a Kerr black hole is given by \citep{Bardeen:1972_347}
\begin{equation}
r_{+}  = M \left[ 3 + Z_2 - \sqrt{ (3 - Z_1)(3+Z_1 + 2Z_2 ) } \right] ,
\end{equation}
where $Z_1$ and $Z_2$ are defined by
\begin{equation}
Z_1 = 1 + (1 - j^2)^{1/3} \left[ (1+j)^{1/3} + (1-j)^{1/3} \right] ,
\end{equation}
\begin{equation}
Z_2 = \sqrt{3 j^2 + Z_1^2 } .
\end{equation}
The more physical circumferential radius $R_{\rm ISCO}^{\rm Kerr}$ of the ISCO is related to
$r_{+}$ by
\begin{equation}
 R_{\rm ISCO}^{\rm Kerr} = \sqrt{ r_{+}^2 + a^2 + { 2 M a^2 \over r_{+} } } ,
 \label{eq:Risco_kerr}
\end{equation}
where the so-called Kerr parameter $a = J/M=Mj$ is the angular momentum per unit mass.

\section{Numerical results}
\label{sec:results}

\subsection{EOS sensitive relations}
\label{sec:eos_sen_rel}

Before presenting our universal relations for the radius and orbital frequency of ISCO, we first study and review how the ISCO generally depends on the EOS and stellar spin frequency
\citep{Miller:1998_793,Bhattacharyya:2011_3247,Torok:2014_L5}.
Let us first consider the situation for nonrotating neutron stars where
the ISCO radius is given analytically by $R_{\rm ISCO} = 6 M$. We must require the stellar radius $R < R_{\rm ISCO}$ in order for the ISCO to exist outside the star. The compactness of a
star needs to be $M/R > 1/6$ in order for the ISCO to exist. Nonrotating neutron stars constructed
from typical EOS models can easily achieve this condition when $M$ is high enough.
The situation for rotating stars is more complicated since one also needs to consider the stellar spin frequency as a parameter. In general, for a given EOS model, there exists a minimum mass $M_{\rm min}$ at each value of the stellar spin frequency $f$ in order for the ISCO to exist (i.e., $R < R_{\rm ISCO}$). In Figure~\ref{fig:Min_M_f}, we plot $M_{\rm min}$ against $f$ to illustrate this fact using some of our EOS models.
It is seen that $M_{\rm min}$ rises with $f$ generally and it also depends sensitively on the EOS. In particular, $M_{\rm min}$ can change by as much as 100\% for different EOSs at high values
of $f$.

While the value of $M_{\rm min}$ depends sensitively on $f$ and the EOS, Figure~\ref{fig:Min_M_f} shows that an ISCO can exist generally only for high mass neutron stars.
For rapidly rotating models with $f > 1500$ Hz, an ISCO exists only for $M_{\rm min} > 2 M_\odot$ for the EOS models plotted in Figure~\ref{fig:Min_M_f}.
Nevertheless, it should be noted that the currently observed fastest rotating neutron star has a
frequency of about 700 Hz \citep{Hessels:2006_1901}. The corresponding smallest possible value
of $M_{\rm min}$ for the existence of ISCO at that frequency is about $1.4 M_\odot$ as shown in
Figure~\ref{fig:Min_M_f}.

\begin{figure}
  \centering
  \includegraphics*[width=9.0cm]{figure1.eps}
  \caption{Minimum required mass $M_{\rm min}$ for the appearance of ISCO is plotted against the neutron star spin frequency $f$ for some of our EOS models.}
  \label{fig:Min_M_f}
\end{figure}

Let us now study the dependence of the ISCO radius on the EOS and stellar properties. Figure~\ref{fig:Risco_j} plots $R_{\rm ISCO}$ (normalized by $6M$) against the spin parameter $j$
for three EOS models. We also consider three different masses 1.6, 1.8, and 2.0 $M_\odot$ for each
EOS. In the figure, we also plot the results given by the slowly rotating limit and the Kerr black hole for comparison.
It is seen that the numerical results for rotating stars agree quite well to the prediction of
Equation~(\ref{eq:Risco_slow}) for slow spin cases $j \lesssim 0.1$. The numerical results also
converge to the nonrotating limit $R_{\rm ISCO} / 6 M = 1$ as $j$ decreases to zero.
However, the results deviate significantly from the prediction of the slowly rotating limit for high values of $j$. It is also apparent that $R_{\rm ISCO}$ becomes very sensitive to the EOS
and the stellar mass as the value of $j$ increases.

It can also be seen from Figure~\ref{fig:Risco_j} that $R_{\rm ISCO}$ for the Kerr black hole (dashed line) behaves quite differently compared to rapidly rotating neutron stars. In particular, for a given value of $j$, the ISCO radius of a rapidly rotating neutron star (black hole) is larger (smaller) than that predicted by the slowly rotating limit (solid line). This difference may not be surprising since it is well known that the exterior spacetime of rapidly rotating neutron stars cannot be approximated by the Kerr metric.

The contributions from higher-order multipole moments of the rotating star become significant for high values of $j$.
Figure~\ref{fig:Risco_j} illustrates that the effects of higher-order multipole moments, presumably dominated by the quadrupole moment $Q$, generally increase the ISCO radius for rotating
neutron stars.
The properties of the ISCO in terms of the multipole moments of the rotating neutron star and the spacetime have been studied \citep{Shibata:1998_104011,Berti:2004_1416,Berti:2005_923}.
In particular, the dependence of $R_{\rm ISCO}$ on $j$ and $Q$ has been studied by \citet{Pappas:2015_4066}. We shall discuss the connection of this work to our proposed universal relations in Section~\ref{sec:other_relations}.

\begin{figure}
  \centering
  \includegraphics*[width=9.0cm]{figure2.eps}
  \caption{ISCO radius $R_{\rm ISCO}$ (normalized by $6M$) is plotted against the dimensionless
  spin parameter $j$ for different masses and EOS models. The results for the slowly rotating
  limit (Equation~(\ref{eq:Risco_slow})) and the Kerr black hole (Equation~(\ref{eq:Risco_kerr}))
  are also shown.   }
  \label{fig:Risco_j}
\end{figure}

\subsection{Universal relations}
\label{sec:uni_relation}

Having seen that the ISCO radius $R_{\rm ISCO}$ for rapidly rotating neutron stars depends sensitively on the EOS and stellar properties such as $M$ and $f$, it is thus somewhat surprising that approximate EOS-insensitive relations connecting the ISCO to the global stellar quantities
can indeed be found. To motivate our proposed universal relations, let us first note that previously known universal relations, as mentioned in Section~\ref{sec:intro}, are relations connecting suitable dimensionless quantities.

In order to search for possible universal relations concerning the ISCO, we construct dimensionless quantities out of the relevant physical variables such as $R_{\rm ISCO}$, $M$, and $f$ in this situation. In particular, we consider the
two quantities $R_{\rm ISCO} f$ and $M f$, which are dimensionless in geometric units.
In Figure~\ref{fig:Risco_uni}, we plot $R_{\rm ISCO} f$ against $M f$ using our 12 different EOS models. In contrast to Figure~\ref{fig:Risco_j}, we now see that the results are relatively
EOS-insensitive for a large range of $Mf$. To facilitate any potential use of this universal relation in astrophysical situations, we have restored physical units in the figure so that $R_{\rm ISCO}$, $M$, and $f$ are expressed in km, $M_\odot$, and Hz, respectively.
The lower panel of Figure~\ref{fig:Risco_uni} shows the relative error, $({\hat y} - y)/y$, between the numerical data $\hat y$ and a fitting curve (solid line) $y$ given by
\begin{equation}
y = a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4  ,
\label{eq:Risco_fit}
\end{equation}
where $y = R_{\rm ISCO} f$ and $x = M f$.
The fitting parameters are $a_1 = 8.809$, $a_2 = -9.166 \times 10^{-4}$, $a_3 = 8.787 \times 10^{-8}$, and $a_4 = - 6.019 \times 10^{-12}$.
The vertical dashed line in the figure corresponds to the current observational upper bound of
$Mf$ in the sense that its value is obtained by combining the maximum measured mass at about $2 M_\odot$ \citep{Demorest:2010_1081,Antoniadis:2013_6131} and the largest spin frequency $f = 716$ Hz for the fastest rotating neutron star PSR J1748-2446ad \citep{Hessels:2006_1901}. We see that the numerical data can generally be fitted by Equation~(\ref{eq:Risco_fit}) to about 2\% for a large range of $Mf$, except for a few outsiders between $Mf\approx 2000$ and 3000 $M_\odot {\rm Hz}$
that have relative errors close to 6\%.

\begin{figure}
  \centering
  \includegraphics*[width=9.0cm]{figure3.eps}
  \caption{Universal relation connecting the ISCO radius $R_{\rm ISCO}$ to the neutron star spin   frequency $f$ and mass $M$. The lower panel shows the relative error between the numerical data and the fitting curve given by Equation~(\ref{eq:Risco_fit}). }
  \label{fig:Risco_uni}
\end{figure}

As we shall see in Section~\ref{sec:other_relations}, to second order in the spin parameter $j$, the ISCO radius $R_{\rm ISCO}$ can be determined by $M$, $j$, and the quadrupole moment $Q$ of the rotating star (see Equation~(\ref{eq:Risco_shibata})). While $M$ is one of the common observables for neutron stars, the other two quantities $j$ and $Q$ cannot be obtained easily from observations. By comparison, our approximate EOS-insensitive relation (Equation~(\ref{eq:Risco_fit})) is more appealing and would be more useful since the relation is also valid for rapidly rotating stars and the relevant quantities $M$ and $f$ are more promising observables for neutron stars. In particular, these quantities can usually be measured to high accuracy if they are observed.

Having seen that there exists a universal relation connecting $R_{\rm ISCO}$ to the global
stellar quantities $M$ and $f$, it is thus not surprising that a similar universal relation should
also exist for the orbital frequency $f_{\rm ISCO}$ of the ISCO, since $f_{\rm ISCO}$ is determined
by the exterior spacetime metric and the position of ISCO around a rotating star.
As we shall discuss in Section~\ref{sec:discuss}, $f_{\rm ISCO}$ is indeed more interesting than $R_{\rm ISCO}$ since it may be directly related to the observed QPOs in LMXBs.
In Figure~\ref{fig:Fisco_uni}, we plot the ratio $f/f_{\rm ISCO}$ against $Mf$ to demonstrate the
universality for $f_{\rm ISCO}$. As before, we have restored physical units so that $M$ is
expressed in $M_\odot$ and the frequencies $f$ and $f_{\rm ISCO}$ are expressed in Hz. The vertical dashed line corresponds the current observational upper bound of $Mf$. The solid line is a fitting curve given by
\begin{equation}
y = b_1 x + b_2 x^2 + b_3 x^3 + b_4 x^4 ,
\label{eq:Fisco_fit}
\end{equation}
where $y = f / f_{\rm ISCO}$ and $x = Mf$.
The fitting parameters are given by $b_1 = 4.497 \times 10^{-4}$, $b_2 = -6.130 \times 10^{-8}$, $b_3 = 4.527 \times 10^{-12}$, and $b_4 = - 1.446 \times 10^{-16}$. The relative error between the numerical data and the fitting curve is plotted in the lower panel of Figure~\ref{fig:Fisco_uni}. Similar to Figure~\ref{fig:Risco_uni}, we see that the relation is approximately EOS-insensitve to
a few percent level in a large range of $Mf$.

\section{Analysis}
\label{sec:analy}

\subsection{Connection to Kerr black hole}
\label{sec:kerr_BH}

As we have seen in Figure~\ref{fig:Risco_j}, the dependence of $R_{\rm ISCO}$ on $M$ and $j$ for rotating neutron stars is very different from that for the Kerr black hole.
The difference is due to the fact that the multipolar structure of the spacetime of a rapidly rotating star \citep{Laarakkers:1999_282,Berti:2004_1416,Pappas:2012_231104,Yagi:2014_124013} is in general different from that of a (uncharged) rotating black hole, which depends only on two parameters, $M$ and $j$, because of the no-hair theorem.
It should be noted that approximate analytic formulas for $R_{\rm ISCO}$ and $f_{\rm ISCO}$ around
a rotating star in terms of the various moments of the star have been obtained by \citet{Shibata:1998_104011}.
For comparison, the approximately EOS-insensitive relations that we find for $R_{\rm ISCO}$ (Equation~(\ref{eq:Risco_fit})) and $f_{\rm ISCO}$ (Equation~(\ref{eq:Fisco_fit})) depend
only on two parameters, $M$ and $f$. It is thus interesting to see whether these universal relations for rotating neutron stars share greater similarities to the corresponding relations for the Kerr black hole or not.

\begin{figure}
  \centering
  \includegraphics*[width=9.0cm]{figure4.eps}
  \caption{Universal relation connecting the ISCO frequency $f_{\rm ISCO}$ to the neutron star
  spin frequency $f$ and mass $M$. The lower panel shows the relative error between the numerical
  data and the fitting curve given by Equation~(\ref{eq:Fisco_fit}). }
  \label{fig:Fisco_uni}
\end{figure}

The ISCO radius $R_{\rm ISCO}$ around a Kerr black hole is determined analytically by
Equation~(\ref{eq:Risco_kerr}) in terms of the mass $M$ and spin parameter $j$ of the black hole.
Since our universal relation (Equation~(\ref{eq:Risco_fit})) for $R_{\rm ISCO}$ around rotating stars involves the spin frequency $f$ instead of $j$, we replace $j$ by the spin frequency
$f_{\rm H} = \Omega_{\rm H} / 2 \pi$ of the black hole horizon \citep[e.g.,][]{Wald:book}
in order to obtain the corresponding equation for a rotating black hole. The angular frequency $\Omega_{\rm H}$ of the black hole horizon can be expressed in terms of $M$ and $j$:
\begin{equation}
2 M \Omega_{\rm H} = { j \over 1 + \sqrt{1 - j^2} } .
\label{eq:Omega_H}
\end{equation}
Solving Equation~(\ref{eq:Omega_H}) for $j$ and using the result in Equation~(\ref{eq:Risco_kerr})
then provide us with the desired relation connecting $R_{\rm ISCO}$, $M$, and $\Omega_{\rm H}$ (and
hence $f_{\rm H}$).
In Figure~\ref{fig:Risco_uni_Kerr}, we compare the universal relation for $R_{\rm ISCO}$ (Equation~(\ref{eq:Risco_fit})) for rotating neutron stars and the corresponding relation for the
Kerr black hole. We see that the two relations agree quite well over a large range of
$Mf$. The relative difference between the two relations is within 6\% up to $Mf = 5000 M_\odot {\rm Hz}$, which is well above the current observational upper bound for rotating
neutron stars.

In contrast to the plot of $R_{\rm ISCO}$ against $j$ in Figure~\ref{fig:Risco_j}, we now see that the relations between the dimensionless quantities $R_{\rm ISCO}f$ and $Mf$ for neutron stars and the Kerr black hole are very close to each other for the range of $Mf$ allowed theoretically
for neutron stars.
This similarity is quite intriguing if one recalls that the multipolar structure of the spacetime
of a rapidly rotating neutron star is generally different from that of the Kerr black hole.
Furthermore, the spin frequency $f_{\rm H}$ of the black hole horizon is conventionally defined
as the orbital frequency of a zero-angular-momentum observer at the horizon due to the effect of
frame dragging \citep{Wald:book}. The concept of $f_{\rm H}$ is thus also quite different from the spin frequency of a rotating star, which is associated with the fluid motion of the star itself and has a clear Newtonian limit.

\begin{figure}
  \centering
  \includegraphics*[width=9.0cm]{figure5.eps}
  \caption{Comparison between the universal relation for $R_{\rm ISCO}$ for neutron stars
  (Equation~(\ref{eq:Risco_fit})) and the corresponding relation for the Kerr black hole.}
  \label{fig:Risco_uni_Kerr}
\end{figure}

\subsection{Connection to other universal relations}
\label{sec:other_relations}

While there is no exact analytic solution for the spacetime of rapidly rotating compact stars
in general relativity, it was shown by \citet{Pappas:2013_3007} that the two-soliton analytic solution of \citet{Manko:1995_3063} can provide an accurate representation of the exterior spacetime of rotating neutron stars.
By further making use of the EOS-insensitive ($\sim 10$\% variability) universal relations between the multipole moments of rapidly rotating neutron stars \citep{Pappas:2014_121101,Yagi:2014_124013}, the formally four-parameter two-soliton solution can lead to an EOS-independent description of the spacetime in terms of the mass $M$, spin parameter $j$, and the reduced quadrupole moment $\alpha\equiv - Q/(j^2 M^3)$ of a rotating neutron star \citep{Pappas:2015_4066}. In retrospect, it may thus not be entirely surprising that there can exist universal relations for the ISCO radius and frequency, although the fact that these relations (Equations~(\ref{eq:Risco_fit})-(\ref{eq:Fisco_fit})) have relatively simple forms is very interesting and unexpected.

In \citet{Pappas:2015_4066}, the properties of ISCO are studied using the above mentioned analytic
spacetime. In particular, the region of the parameter space $(j, \alpha)$ for which the ISCO is outside the stellar surface ($R_{\rm ISCO} > R$) is investigated. It is found that $\alpha$ is
restricted to be smaller than about 4.5 for small spins at $j \approx 0.1$ and it is
further restricted to $\alpha \lesssim 2$ at $j \approx 0.7$ (see Figure 3 of \citet{Pappas:2015_4066}). For comparison, we note that $\alpha=1$ for the Kerr black hole.
As we shall discuss later, the restriction on the value of $\alpha$ might be responsible for the
similarity between Equation~(\ref{eq:Risco_fit}) and the corresponding relation for the Kerr
black hole.

A formula for the ISCO radius around rotating relativistic stars as an expansion in the mulitpole moments of the spacetime has been developed by \citet{Shibata:1998_104011}.
\citet{Pappas:2015_4066} found that, for the region of the parameter space $(j, \alpha)$ where $R_{\rm ISCO} > R$, the \cite{Shibata:1998_104011} formula agrees to the prediction of the two-soliton solution to within 1\%. In the following, we shall make use of the expansion formula and other universal relations for rotating neutron stars to provide a connection to our universal relation for $R_{\rm ISCO}$.

The \citet{Shibata:1998_104011} expansion formula for $R_{\rm ISCO}$, up to the third order, is
given by
\begin{eqnarray}
R_{\rm ISCO} &=& 6 M \left( 1 - 0.54433 j - 0.22619 j^2 + 0.17989 \alpha j^2 \right. \nonumber \\
\cr
&& \left. - 0.23002 j^3  + 0.26296 \alpha j^3 - 0.05317 j^3 \beta \right) ,
\label{eq:Risco_shibata}
\end{eqnarray}
where $\beta \equiv - J_3 /(j^3 M^4)$ and $J_3$ is the spin octupole moment. As discussed by
\citet{Pappas:2015_4066}, the parameter $\beta$ can be given in terms of $\alpha$ using the
universal relation between the quadrupole and spin octupole moments \citep{Pappas:2014_121101}:
\begin{equation}
\beta^{1/3} = -0.36 + 1.48 \alpha^{0.65/2} .
\label{eq:Q_J3}
\end{equation}
Equation~(\ref{eq:Risco_shibata}) provides us with a starting point to construct a universal relation for the ISCO radius. For the Kerr black hole, we use Equation~(\ref{eq:Omega_H}) to replace $j$ by $Mf$ in order to obtain the corresponding universal relation for the ISCO radius. For rotating neutron stars, we can make use of the relation connecting $Mf$, $j$, and $\alpha$ found numerically by \citet{Pappas:2015_4066}:
\begin{eqnarray}
{ {Mf} \over j} &=& B_0 + B_1 j + B_2 j^2 + \left( A_0 + A_1 j + A_2 j^2 \right)
\left( \sqrt{\alpha} \right)^{n_1}  \nonumber \\
\cr
&& + \left( C_0 + C_1 j + C_2 j^2 \right) \left( \sqrt{\alpha} \right)^{n_2} ,
\label{eq:Pappas_Mfj}
\end{eqnarray}
where the fitting parameters are $A_0=15.0297$, $A_1=-0.114154$, $A_2=-7.72439$, $B_0=-1.48338$,
$B_1=-1.07874$, $B_2=1.64592$, $C_0=-2.45303$, $C_1=3.40995$, $C_2=7.39354$, $n_1=-1.1698$,
and $n_2 = -4.50216$. This relation is EOS-insensitive to within the 1\% level.

In order to use Equations~(\ref{eq:Risco_shibata})$-$(\ref{eq:Pappas_Mfj}) to construct a relation
between $R_{\rm ISCO}f$ and $Mf$, we still need one more relation to connect $\alpha$ and $j$.
As discussed above, \citet{Pappas:2015_4066} studied the properties of ISCO based on the two-soliton
analytic solution and in particular he considered the contours of constant relative difference
between $R_{\rm ISCO}$ and $R$, defined by the parameter
$\delta R = (R_{\rm ISCO} - R)/R_{\rm ISCO}$, in the parameter space of $(j, \alpha)$. Note that the region with $\delta R \geq 0$ corresponds to the parameters where the ISCO can exist outside the
stellar surface.
We refer the reader to Figure 3 of \citet{Pappas:2015_4066} to see how different values of $\delta R$ define different EOS-independent contours in the parameter space.
Those contours with $\delta R \geq 0$ thus provide us with the desired relations between $\alpha$ and $j$ connected to the ISCO. We consider the contours defined by $\delta R = 0$, 0.11, 0.22, and 0.33 in Figure 3 of \citet{Pappas:2015_4066} and use them to express $\alpha$ in terms of $j$ numerically in Equation~(\ref{eq:Pappas_Mfj}), which is then solved together with Equations~(\ref{eq:Risco_shibata}) and (\ref{eq:Q_J3}) in order to construct a relation between $R_{\rm ISCO}f$ and $Mf$.

\begin{figure}
  \centering
  \includegraphics*[width=9.0cm]{figure6.eps}
  \caption{Universal relations constructed from
  Equations~(\ref{eq:Risco_shibata})$-$(\ref{eq:Pappas_Mfj}) for different constant
  $\delta R$ contours in the parameter space $(j,\alpha)$ extracted from Figure 3 of
  \citet{Pappas:2015_4066}. Equation~(\ref{eq:Risco_fit}) and the corresponding relation for the Kerr black hole are also plotted for comparison. The lower panel shows the relative errors, taking
  Equation~(\ref{eq:Risco_fit}) as a reference.   }
  \label{fig:Risco_uni_pappas}
\end{figure}

It should be pointed out that different contours of constant $\delta R$ lead to different relations between $Mf$ and $j$ using Equation~(\ref{eq:Pappas_Mfj}). There is thus no a priori reason to guarantee that our procedure would lead to a single universal relation connecting $R_{\rm ISCO}f$
and $Mf$.
In the upper panel of Figure~\ref{fig:Risco_uni_pappas}, we plot the resulting relations between $R_{\rm ISCO}f$ and $Mf$ obtained for the different contours.
Equation~(\ref{eq:Risco_fit}) and the corresponding relation for the Kerr black hole are also
plotted in the figure for comparison. The relative errors, taking Equation~(\ref{eq:Risco_fit}) as a reference, are shown in the lower panel of Figure~\ref{fig:Risco_uni_pappas}.
It is interesting to note that the relations corresponding to different contours of $\delta R$ turn out to match each other very well and establish approximately the existence of a single universal relation.
Deviations between the different $\delta R$ curves start to become noticeable at the high end of $Mf$. These $\delta R$ curves also match Equation~(\ref{eq:Risco_fit}) very well to within about
2\% for a large range of $Mf$, which is consistent with the accuracy level of Equation~(\ref{eq:Risco_fit}) as a fitting curve to the numerical data.

In Figure~\ref{fig:Risco_uni_pappas}, it can be seen that the different $\delta R$ curves generally tend toward the Kerr case as the value of $\delta R$ increases.
This is consistent with the fact that the value of $\alpha$, at a fixed value of $j$, on a contour of constant $\delta R$ decreases as the value of $\delta R$ increases as shown in Figure 3 of
\citet{Pappas:2015_4066}. In particular, the maximum value of $\alpha$ on a contour decreases from about 4.5 to 2 as $\delta R$ increases from 0 to 0.33, and thus the value of $\alpha$ indeed tends toward the Kerr case $\alpha = 1$.
The restriction on the value of $\alpha$ for the existence of an ISCO might be responsible for the
similarity between Equation~(\ref{eq:Risco_fit}) and the corresponding relation for the Kerr black
hole. However, further work is needed to study this connection in more detail.

\section{Discussion}
\label{sec:discuss}

As mentioned in Section~\ref{sec:eos_sen_rel}, an ISCO can exist only around high mass neutron stars. Depending on the underlying EOS model, the minimum mass for the existence of an ISCO can increase from about 1.4 to $2 M_\odot$ as the spin frequency of the star increases from 700 Hz to 1500 Hz.
In this paper, we have proposed two universal relations to connect the ISCO radius and orbital
frequency (when an ISCO exists) to the spin frequency and mass of rotating neutron stars. We have
also compared the universal relation for the ISCO radius to the corresponding relation for the Kerr black hole and analyze it from the perspective of other known universal relations for rotating neutron stars. In the following, we shall discuss how the ISCO universal relations can be applied practically.

As mentioned briefly in Section~\ref{sec:intro}, the ISCO may be relevant to the kHz QPOs well observed from neutron stars in LMXBs. These QPOs often come in pairs, one with a lower frequency $f_{\rm l}$ and another with a higher frequency $f_{\rm u}$. Many models have been proposed to explain the kHz QPOs. We refer the reader to \citet{Klis:2006_39}, \citet{Torok:2016_273}, and
references therein for reviews.

While there is still no general consensus on the mechanism for generating the kHz QPOs, it is believed that the QPOs are associated with the orbital motion and/or oscillations near the inner edge of an accreting disk surrounding the central neutron star.
The ISCO defines the smallest inner radius that an accreting disk around a neutron star can have. Indeed, there is also evidence to suggest that the ISCO is responsible for some of the
observed properties associated to the kHz QPOs. For instance, it has been observed
that the quality factor of the kHz QPOs drops sharply in some systems, which is consistent with
the signature that the orbits involved approach the ISCO \citep{Barret:2006_1140,Barret:2007_1139, Boutelier:2009_1901}. Our proposed universal relations for the ISCO are relevant and applicable to these systems.

As an illustration of potential applications of our proposed universal relations, we apply Equation~(\ref{eq:Fisco_fit}) to the system 4U 0614+09. The spin frequency $f=414$ Hz of the neutron star in the system has been inferred from the detection of burst oscillations during a thermonuclear X-ray burst \citep{Strohmayer:2008_L37}. The highest frequency for the upper kHz QPO is measured to be $f_{\rm u} \approx 1220$ Hz \citep{Boutelier:2009_1901}. Using these two observed frequencies $f$ and $f_{\rm u}$, together with the assumption that $f_{\rm u}$ is equal to the ISCO orbital frequency $f_{\rm ISCO}$, we apply Equation~(\ref{eq:Fisco_fit}) to determine that the mass of the neutron star is $2.0 M_\odot$. We can then determine the ISCO radius to be about 16 km by using Equation~(\ref{eq:Risco_fit}), and hence put a constraint on the radius of the star to be
$R < 16$ km.

If, in reality, the QPO frequency $f_{\rm u}$ corresponds to the orbital frequency at an orbit outside the ISCO, instead of being equal to $f_{\rm ISCO}$, then our value $2.0 M_\odot$ would
set an upper bound for the mass of the star.
Previously, \citet{Boutelier:2009_1901} also estimated the mass of 4U 0614+09 to be about $1.9 M_\odot$ by using the value $f_{\rm ISCO}=1250$ Hz, which is obtained by using the observed drop of the quality factor of the lower kHz QPO and extrapolating it to zero, together with the constant
frequency difference between the lower and upper kHz QPOs.
In their analysis, \citet{Boutelier:2009_1901} determined the mass by assuming slow rotation and included the first-order correction due to the dimensionless spin parameter $j$, the value of which generally depends on the mass and EOS for a given spin frequency of the star. For comparison, besides the minimal assumption that $f_{\rm u} = f_{\rm ISCO}$, our analysis applies generally to rapidly rotating stars and is insensitive to the underlying EOS model to within about the 2\% level.

If the existence of an ISCO around some neutron stars in LMXBs can be confirmed, together with
the measurement of the ISCO frequency (or radius), our universal relations can be used to yield a
precise determination of the masses of the neutron stars, bypassing the assumption of slow
rotation and the uncertainty related to the dimensionless spin parameter, which are commonly
required in the literature.

We thank the anonymous referee for very helpful comments that improved the quality of the manuscript. In particular, we thank the referee for bringing the work of \citet{Pappas:2015_4066} to our attention. Some of our chosen EOSs are obtained from the EOS repository CompOSE (compose.obspm.fr). We thank the team members of CompOSE and LORENE for their development.


text:
\section{Introduction} \label{sec.Intro}
A complex Hadamard matrix is a square matrix
$W$ of order $n$ which satisfies $W\overline{W}^\top= nI$ and
all of whose entries are complex numbers of absolute value $1$.
A complex Hadamard matrix is said to be Butson-type,
if all of its entries are roots of unity.
In our earlier work \cite{MI}, we proposed a method to classify symmetric
complex Hadamard matrices belonging to the Bose--Mesner algebra
of a symmetric association scheme.
In this paper, we propose an analogous method to classify nonsymmetric
hermitian complex Hadamard matrices belonging to the Bose--Mesner algebra
of a commutative nonsymmetric association scheme.
First we consider nonsymmetric hermitian complex Hadamard matrices
belonging to the Bose--Mesner algebra of
a commutative nonsymmetric association scheme of class $3$,
and then we show that the first eigenmatrix of such a scheme is characterized as
\begin{equation}\label{d3P1}
\begin{pmatrix}
1&a(2a-1)&a(2a-1)&2a-1\\
1&ai&-ai&-1\\
1&-ai&ai&-1\\
1&-a&-a&2a-1
\end{pmatrix},
\end{equation}
where $a$ is a positive integer and $i=\sqrt{-1}$.
Moreover, we show that such a complex Hadamard matrix is necessarily
a Butson-type complex Hadamard matrix
whose entries are $4$-th roots of unity.
An association scheme with the first eigenmatrix (\ref{d3P1}) is
a nonsymmetric amorphous association scheme
belonging to $L_{2a;1}$, according to \cite{IMY}.
An example of such an association scheme was
constructed from a Galois ring of characteristic $4$  in \cite[Theorem 9]{IMY}.
Galois rings have been used to construct certain association
schemes. Yamada \cite{Y} used Galois rings of characteristic $4$
to construct distance-regular digraphs, which are nonsymmetric
association schemes of class $3$. This construction was generalized in
\cite{IMY} to produce amorphous association schemes. Ma \cite{Ma}
gave association schemes of class $3$ on Galois rings of characteristic $4$,
which do not come from fusions in any amorphous association schemes.
Moreover, certain properties of association schemes
obtained from Galois rings of odd characteristic
have been investigated in \cite{EP}.

Our second result of this paper is a construction of
nonsymmetric association schemes $\mathfrak{X}$ of
class $6$ on Galois rings of characteristic 4,
whose first eigenmatrix is given by
\begin{equation} \label{P6}
(p_{i,j})_{\substack{0\leq i\leq 6 \\ 0\leq j\leq 6}}=\begin{pmatrix}
1&2b(b-1)&2b(b-1)&b&b&b-1&b\\
1&bi&-bi&0&0&-1&0\\
1&-bi&bi&0&0&-1&0\\
1&0&0&bi&-bi&b-1&-b \\
1&0&0&-bi&bi&b-1&-b \\
1&-2b&-2b&b&b&b-1&b \\
1&0&0&-b&-b&b-1&b
\end{pmatrix},
\end{equation}
where $b$ is a power of $4$.
We also classify hermitian complex Hadamard matrices belonging to
the Bose--Mesner algebra of $\mathfrak{X}$.
We show that such a matrix is necessarily a Butson-type matrix
whose entries are $4$-th roots of unity.

\section{Association schemes and complex Hadamard matrices} \label{sec.AS}
In this section,
we consider hermitian matrices
belonging to the Bose--Mesner algebra of a commutative association scheme.
Assuming that all
entries are complex numbers with absolute value $1$,
we find conditions under which such a matrix is a complex Hadamard matrix.
We refer the reader to \cite{BI,BCN} for undefined terminology in the
theory of association schemes.

Let $X$ be a finite set with $n$ elements, and let
$\mathfrak{X}=(X,\{R_i\}_{i=0}^d)$ be a commutative association scheme
with the first eigenmatrix
$P=(P_{i,j})$.
We let $\mathfrak{A}$ denote the Bose--Mesner algebra
spanned by the adjacency matrices $A_0,A_1,\ldots,A_d$ of $\mathfrak{X}$.
Then the adjacency matrices are expressed as
\begin{equation}\label{eq:pij}
A_j=\sum_{i=0}^d P_{i,j}E_i \quad (j=0,1,\ldots,d),
\end{equation}
where $E_0=\frac{1}{n}J,E_1,\ldots,E_d$ are the primitive idempotents of $\mathfrak{A}$.
Let $w_0=1$ and $w_1,\ldots,w_d$ be complex numbers with absolute value $1$.
Set
\begin{equation}\label{eq:W}
W=\sum_{j=0}^d w_jA_j\in\mathfrak{A}.
\end{equation}
Define
\begin{equation}\label{eq:gamma}
\gamma_k=\sum_{j=0}^d w_jP_{k,j} \quad (k=0,1,\ldots,d).
\end{equation}
By \eqref{eq:pij}, \eqref{eq:W} and \eqref{eq:gamma} we have
\begin{equation}\label{eq:WE}
W=\sum_{k=0}^d \gamma_kE_k.
\end{equation}
Let $X_j$ {\rm($1\leq j\leq d$)} be indeterminates. For $k=1,2,\dots,d$,
let $e_k$ be the polynomial defined by
\begin{equation}\label{eq:ek}
e_k=1+2\left(\sum_{j=1}^d P_{k,j}X_j+
\sum_{1\leq j_1<j_2\leq d}P_{k,j_1}P_{k,j_2}X_{j_1}X_{j_2}\right)
+\sum_{j=1}^dP_{k,j}^2X_j^2-n.
\end{equation}
Then we have the following.

\begin{lem}\label{lem:equiv}
Assume that the matrix $W$ given in \eqref{eq:W} is a hermitian matrix.
Then the following statements are equivalent:
\begin{itemize}
\item[\rm{(i)}] $W$ is a complex Hadamard matrix,
\item[\rm{(ii)}] $\gamma_k^2=n$ for $k=1,\ldots,d$,
\item[\rm{(iii)}] $(w_j)_{1\leq j\leq d}$ is a common zero of $e_k$ {\rm($k=1,\ldots,d$)}.
\end{itemize}
\end{lem}
\begin{proof}
By \eqref{eq:WE} we have $W\overline{W}^T=W^2=\sum_{k=0}^d\gamma_k^2E_k$.
Therefore, (i) implies (ii).
To prove the converse, it is enough to show that $\gamma_0^2=n$.
Since
\begin{align*}
(W^2)_{j,j}&=\sum_{k=1}^nW_{j,k}\overline{W_{j,k}}\\
&=n,
\end{align*}
the diagonal entries of $W^2$ are all $n$.
Thus,
\begin{align*}
n^2&=\tr W^2\\
&=\sum_{k=0}^d\gamma_k^2\tr E_k\\
&=\gamma_0^2+\sum_{k=1}^dn\tr E_k\\
&=\gamma_0^2+n\tr(I-E_0)\\
&=\gamma_0^2+n(n-1).
\end{align*}
Hence $\gamma_0^2=n$.
By \eqref{eq:gamma} we have
\[
\gamma_k^2=1+2\left(\sum_{j=1}^d P_{k,j}w_j+
\sum_{1\leq j_1<j_2\leq d}P_{k,j_1}P_{k,j_2}w_{j_1}w_{j_2}\right)
+\sum_{j=1}^dP_{k,j}^2w_j^2
\]
for $k=1,\dots,d$.
Therefore, the equivalence of (ii) and (iii) follows.
\end{proof}

\section{Hermitian complex Hadamard matrices
 and nonsymmetric association schemes of class $3$}\label{sec.3-class}
In this section, we classify nonsymmetric hermitian complex Hadamard matrices
belonging to the Bose--Mesner algebra of a nonsymmetric commutative association scheme
$\mathfrak{X}=(X,\{R_i\}_{i=0}^3)$ of class $3$,
where $R_1^\top=R_2$, $R_3$ symmetric.
We use the following result.

\begin{lem}[{\cite[Lemma in (5.3)]{S}}]\label{lem:song}
Let
\[
\begin{pmatrix}
1&\frac{k_1}{2}&\frac{k_1}{2}&k_2\\
1&\frac12(r+bi)&\frac12(r-bi)&-(r+1)\\
1&\frac12(r-bi)&\frac12(r+bi)&-(r+1)\\
1&\frac{s}{2}&\frac{s}{2}&-(s+1)
\end{pmatrix}
\]
be the first eigenmatrix of $\mathfrak{X}$, where $r,s$ are integers and $i^2=-1$.
Then one of the following holds.
\begin{itemize}
\item[\rm{(i)}] $(r,s,b^2)=(0,-(k_2+1),\frac{k_1(k_2+1)}{k_2})$,
\item[\rm{(ii)}] $(r,s,b^2)=(-(k_2+1),0,(1+k_2)(1+k_1+k_2))$,
\item[\rm{(iii)}] $(r,s,b^2)=(-1,k_1,k_1+1)$.
\end{itemize}
\end{lem}

Note that, if we put
\begin{equation}\label{eq:0930}
(k_1,k_2)=(2a(2a-1),2a-1),
\end{equation}
in part (i) of Lemma~\ref{lem:song},
we have the matrix (\ref{d3P1}).

Let $\mathfrak{A}$ be the Bose--Mesner algebra of $\mathfrak{X}$
which is the linear span of the adjacency matrices
$A_0,A_1,A_2,A_3$ of $\mathfrak{X}$, where $A_1^\top=A_2$, $A_3$ symmetric.
Let $w_0=1$ and $w_1,w_2,w_3$ be complex numbers of absolute value $1$.
Set
\begin{equation} \label{eq:w0W}
W=w_0A_0+w_1A_1+w_2A_2+w_3A_3\in\mathfrak{A}.
\end{equation}
Then our first main theorem is the following.

\begin{thm}\label{thm:main0925}
Assume that the matrix \eqref{eq:w0W} is a
hermitian complex Ha\-da\-mard matrix and not a real Ha\-da\-mard matrix.
Then $\mathfrak{X}$ is a nonsymmetric association scheme
whose unique nontrivial symmetric relation
consists of $2a$ cliques of size $2a$,
and the first eigenmatrix of $\mathfrak{X}$ is given by \eqref{d3P1}.
Moreover, $w_1=\pm i$ and $w_3=1$.
\end{thm}

\begin{proof}
Suppose that the matrix $W$ is a complex Hadamard matrix.
Without loss of generality, we may assume $w_0=1$ in (\ref{eq:w0W}).
Since $W$ is a hermitian matrix, we have $w_1w_2=1$ and $w_3=\pm1$.
Since $W$ is not a real Hadamard matrix, we have
$w_1-w_2\not=0$.
By Lemma~\ref{lem:equiv}, $(w_i)_{i=1}^3$ is a common zero of the polynomials $e_k$ ($k=1,2,3$) defined in (\ref{eq:ek}).
Since
\[ e_1-e_2=bi(X_1-X_2)(r(X_1+X_2-2X_3)-2X_3+2), \]
we have
\begin{equation}\label{eq:1-1}
r(w_1+w_2-2w_3)-2w_3+2=0.
\end{equation}

First assume that $w_3=1$.
Then by (\ref{eq:1-1}) we have $r(w_1+w_2-2)=0$.
Hence $r=0$.
Therefore we have case (i) of Lemma~\ref{lem:song}.
Then
\begin{equation}\label{eq:nk1}
(k_1,k_2)=\left(\frac{(s+1)b^2}{s},-(s+1)\right).
\end{equation}
After specializing $X_3=1$, we have
\[
e_3-e_1=\frac{1}{4}((s-bi)X_1+(s+bi)X_2-2s)((s+bi)X_1+(s-bi)X_2-2s).
\]
Hence $((s-bi)w_1+(s+bi)w_2-2s)((s+bi)w_1+(s-bi)w_2-2s)=0$.
Then by $w_2=w_1^{-1}$ we have
\begin{equation}\label{eq:1-2}
(w_1-1)^2((s-bi)w_1-(s+bi))((s+bi)w_1-(s-bi))=0.
\end{equation}
Put $w=(s+2bi)/(s-2bi)$.
By $w_1\not=1$ we have $w_1\in\{w,\overline{w}\}$ by (\ref{eq:1-2}).
We may assume $w_1=w$.
After specializing $X_1=w$, $X_2=\overline{w}$ and $X_3=1$, we have
\[
e_1=-(1+k_1+k_2)-\frac{4b^4s^2}{(b^2+s^2)^2}.
\]
Hence
\begin{equation}\label{eq:1-3}
k_1+k_2=-1-\frac{4b^4s^2}{(b^2+s^2)^2}.
\end{equation}
Then by (\ref{eq:nk1}), (\ref{eq:1-3}) we have
\[
(b-s)(b+s)((s+1)b^4-(s-2)s^2b^2+s^4)=0.
\]
Since $b>0$ and $s=-(k_2+1)<-1$, we have $b-s\not=0$.
First assume that $s=-b$.
Then, by putting $b=2a$, we have $(k_1,k_2)=(2a(2a-1),2a-1)$ and $w_1=-i$.
Therefore we have the first eigenmatrix (\ref{d3P1}).
Next assume that
\begin{equation}\label{eq:0921}
(s+1)b^4-(s-2)s^2b^2+s^4=0.
\end{equation}
Then the discriminant of (\ref{eq:0921}) as a quadratic equation in $b^2$ is $s^4((s-4)^2-16)$.
Since $b^2$ is rational, $(s-4)^2-16$ is a square. This implies $s\in\{0,-1\}$.
This contradicts $s=-(k_2+1)<-1$.

Next assume that $w_3=-1$.
Then by (\ref{eq:1-1}) we have $r(w_1+w_2+2)+4=0$.
Put $g(x)=rx^2+2(r+2)x+r$.
By $w_2=w_1^{-1}$ we have $g(w_1)=0$.
Since $w_1\not\in\mathbb{R}$, we have $r+1<0$.
Thus we have case (ii) of Lemma~\ref{lem:song}.
After specializing $X_1=w_1$, $X_2=w_1^{-1}$, $X_3=-1$ and $s=0$, we have
\[
 e_1-e_3=\frac{(w_1+1)f(w_1)}{4w_1^2},
\]
where $f(x)=((r+bi)x+r-bi)((r+bi)x^2+2(r+4)x+r-bi)$.
Since $w_1\not=-1$, we have $f(w_1)=0$.
Since
\begin{align*}
f(x)&=\frac{(r+bi)((r+bi)rx+r(r+4)-(3r+4)bi)}{r^2}g(x)\\
&\quad-\frac{4((r+1)b^2+r^2)}{r^2}((r+4)x-r)
\end{align*}
and $(r+4)w_1-r\not=0$,
we have
\begin{equation}\label{eq:0918}
(r+1)b^2+r^2=0.
\end{equation}
Substituting $(r,b^2)=(-k_2+1,(1+k_1+k_2)(k_2+1))$ in (\ref{eq:0918}),
we have $(k_2+1)(k_2^2+k_1k_2-1)=0$.
This is a contradiction.
\end{proof}

In the next section,
we will construct a nonsymmetric association scheme $\mathfrak{X}$ of class $6$
on a Galois ring of characteristic $4$, with the first eigenmatrix
\eqref{P6}.
An association scheme with the first eigenmatrix (\ref{d3P1}) can be obtained
by fusing some relations of $\mathfrak{X}$.

\section{Association schemes on Galois rings} \label{sec.Galoisring}

For the reminder of this section, we let $e\geq 3$ be an odd positive integer.
We refer the reader to \cite{Y} for basic theory of Galois rings.
Let $F=\text{GF}(2)$ be the prime field of characteristic $2$ and
$K=\text{GF}(2^e)$ be an extension of degree $e$.
We identify $K$ with $F[x]/(\varphi(x))$,
where $\varphi(x)$ is a primitive polynomial of degree $e$ over $F$,
and we denote by $\zeta$ a root of $\varphi(x)$ in $K$.

Let $\mathcal{A}=\mathbb{Z}/4\mathbb{Z}$.
There exists a monic polynomial $\Phi(x)$ over $\mathcal{A}$ such that $\Phi(x)\equiv \varphi(x)$ mod $2\mathcal{A}[x]$
and $\Phi(x)$ divides $x^{2^e-1}-1$ in $\mathcal{A}[x]$ (see \cite[Theorem 1]{CS}).
The ring $\mathcal{R}=\mathcal{A}[x]/(\Phi(x))$ is called a Galois ring,
and it is a local ring with maximal ideal $\mathcal{P}=2\mathcal{R}$ and
residue field $\mathcal{R}/\mathcal{P}\cong K$.

Let $\xi$ be the image of $x$ in $\mathcal{R}$, so that $\xi+\mathcal{P}$ is mapped to $\zeta$
under the isomorphism $\mathcal{R}/\mathcal{P}\cong K$.
Then $\mathcal{R}=\mathcal{A}[\xi]$ and $\xi$ has order $2^e-1$ in $\mathcal{R}$.
Let $\mathcal{T}$ be the cyclic group of order $2^e-1$ generated by $\xi$.
Since $\mathcal{T}$ is mapped bijectively onto $K\setminus\{0\}$ under the natural homomorphism $\mathcal{R}\rightarrow K$,
each element $\alpha\in \mathcal{R}$ is uniquely expressed as
\[
\alpha=\alpha_0+2\alpha_1, \quad \alpha_0,\alpha_1\in\mathcal{T}\cup\{0\}.
\]
The unit group $\mathcal{R}^{\ast}$ is the direct product of $\mathcal{T}$
and the principal unit group $\mathcal{E}=1+\mathcal{P}$.
Let $\mathcal{T}_\delta=\{\xi^j\mid 0\leq j\leq 2^e-2,\;\mathrm{Tr}(\zeta^j)=\delta\}$ ($\delta=0,1$).
We set $\mathcal{P}_0=2\mathcal{T}_0\cup\{0\}$ and $H=1+\mathcal{P}_0$.
Thus $\mathcal{P}_0$ and $H$ are subgroups of index $2$ in the additive group $\mathcal{P}$ and
the multiplicative group $\mathcal{E}$, respectively.
Setting $b=2^{e-1}$, we have $b=|\mathcal{P}_0|=|H|$.

The mapping $\xi\mapsto\xi^2$ can be extended to a ring automorphism $f$ of $\mathcal{R}$
which fixes $\mathcal{A}$ elementwise.
The ring automorphism $f$ is called the {\em Frobenius automorphism}.
For $\alpha=\alpha_0+2\alpha_1$ ($\alpha_0,\alpha_1\in\mathcal{T}\cup\{0\}$), we have
\[
\alpha^f=\alpha_0^2+2\alpha_1^2,
\]
and the trace of $\alpha\in\mathcal{R}$ is defined by
\[
\mathcal{S}(\alpha)=\alpha+\alpha^f+\cdots+\alpha^{f^{e-1}}.
\]
Note that $\mathcal{S}$ is an $\mathcal{A}$-linear mapping from $\mathcal{R}$ to $\mathcal{A}$.
The additive characters of $\mathcal{R}$ are given by $\alpha\mapsto\chi(\alpha\beta)$
($\beta\in\mathcal{R}$), where
\[
\chi(\alpha)=i^{e\mathcal{S}(\alpha)},
\]
and $i^2=-1$. Define
\[
\lambda_{\alpha}(A)=\sum_{\beta\in A}\chi(\alpha\beta),
\]
for $\alpha\in\mathcal{R}$ and $A\subset\mathcal{R}$.
Set $\lambda=\lambda_1$.

Observe
\begin{align}
\mathcal{E}&=H\cup(-H), \nonumber \\
\alpha\mathcal{T}&=\mathcal{P}\setminus\{0\} \quad \text{for any $\alpha\in\mathcal{P}\setminus\{0\}$}, \label{eq:aT} \\
\mathrm{Tr}(\zeta^j)&\equiv \mathcal{S}(\xi^j)\bmod{2}. \nonumber
\end{align}
Since $e$ is odd, we obtain
\begin{align}
\chi(\alpha)&=
\begin{cases}
1 &  \ \text{if} \ \alpha\in\mathcal{P}_0, \\
-1 &  \ \text{if} \ \alpha\in\mathcal{P}\setminus\mathcal{P}_0,
\end{cases}\label{chival} \\
\chi(\alpha)&=i \quad \text{for} \ \alpha\in H. \label{chiH}
\end{align}

\begin{lem}  \label{lem:aP0}
For any $\alpha\in\mathcal{R}^{\ast}\setminus\mathcal{E}$ we have the following.
\begin{align}
\lambda(\alpha\mathcal{P}_0)&=0, \label{eq:aP0} \\
\lambda(\alpha H)&=0. \label{eq:aH}
\end{align}
\end{lem}
\begin{proof}
Note that the correspondence $2\xi^j\mapsto \zeta^j$ extends to an
isomorphism from $\mathcal{P}$ to $K$ as additive groups. Under this isomorphism,
$\mathcal{P}_0$ is mapped to $\mathrm{Tr}^{-1}(0)$. Since
$\alpha\notin\mathcal{E}$, the image $a=\alpha+\mathcal{P}$ of $\alpha$ in $K=\mathcal{R}/\mathcal{P}$
is not $1$. By the non-degeneracy of the trace function on $K$, we have
$a\mathrm{Tr}^{-1}(0)\neq\mathrm{Tr}^{-1}(0)$. This implies
$\alpha\mathcal{P}_0\neq\mathcal{P}_0$.
Then
\[A=\alpha\mathcal{P}_0\cap\mathcal{P}_0\]
is a subgroup of index $4$ in the additive group $\mathcal{P}$.
By \eqref{chival} we have
\[
\chi(\beta)=
\begin{cases}
    1 & \text{if} \quad \beta\in A, \\
    -1 & \text{if} \quad \beta\in\alpha\mathcal{P}_0\setminus A.
\end{cases}
\]
Therefore $\lambda(\alpha\mathcal{P}_0)=0$.
Since $\alpha H=\alpha+\alpha\mathcal{P}_0$, we have
$\lambda(\alpha H)=\lambda(\alpha)\lambda(\alpha\mathcal{P}_0)
=0$.
\end{proof}

Recall $b=|\mathcal{P}_0|=|H|$.

\begin{lem} \label{lem:lamH}
We have the following.
\begin{itemize}
\item[\rm{(i)}] $\lambda(H)=bi$,
\item[\rm{(ii)}] $\lambda(\mathcal{T} H)=bi$.
\end{itemize}
\end{lem}
\begin{proof}
(i) Obvious from \eqref{chiH}.

(ii) Since $\mathcal{T}\setminus\{1\}\subset\mathcal{R}^{\ast}\setminus\mathcal{E}$,
we have
$\lambda((\mathcal{T}\setminus\{1\})H)=0$
by \eqref{eq:aH}. Then the result follows from  (i).
\end{proof}

\begin{lem} \label{lem:lamaH}
We have
\[
\lambda_{\alpha}(H)=
\begin{cases}
b & \text{if} \quad \alpha\in\mathcal{P}_0, \\
-b & \text{if} \quad \alpha\in\mathcal{P}\setminus\mathcal{P}_0.
\end{cases}
\]
\end{lem}
\begin{proof}
Let $\alpha\in\mathcal{P}$. Then $\alpha H=\{\alpha\}$.
Hence $\lambda_{\alpha}(H)=|H|\chi(\alpha)$, and
the result follows from \eqref{chival}.
\end{proof}

\begin{lem}  \label{lem:lamaTH}
For any $\alpha\in\mathcal{P}\setminus\{0\}$ we have
\[ \lambda_{\alpha}(\mathcal{T} H)=-b. \]
\end{lem}
\begin{proof}
By \eqref{eq:aT} we have
\begin{align*}
\lambda_{\alpha}(\mathcal{T} H)&=\sum_{\beta\in\mathcal{P}\setminus\{0\}}\lambda_{\beta}(H) \\
&=\sum_{\beta\in\mathcal{P}}\lambda_{\beta}(H)-\lambda_0(H) \\
&=-\lambda_0(H) && (\text{by Lemma~\ref{lem:lamaH}}) \\
&=-b.
\end{align*}
\end{proof}

Then our second main theorem is the following.
\begin{thm} \label{thm:1}
Let
\begin{align*}
S_0&=\{0\}, \\
S_1&=(\mathcal{T}\setminus\{1\})H, \\
S_2&=(\mathcal{T}\setminus\{1\})(-H), \\
S_3&=H, \\
S_4&=-H, \\
S_5&=\mathcal{P}_0\setminus\{0\}, \\
S_6&=\mathcal{P}\setminus\mathcal{P}_0.
\end{align*}
Then $\lambda_{\alpha}(S_j)$ is constant for any $\alpha\in S_i$ {\rm($i=0,1,\ldots,6$)}.
Put $p_{i,j}=\lambda_{\alpha}(S_j)$ for $\alpha\in S_i$ {\rm($i,j=0,1,\ldots,6$)}.
Then $P=(p_{i,j})_{\substack{0\leq i\leq 6 \\ 0\leq j\leq 6}}$
 is given by {\rm(\ref{P6})}.
In particular, if we set
\[
R_j=\{(\alpha,\beta)\in\mathcal{R}\times\mathcal{R} \mid \alpha-\beta\in S_j \} \quad (j=0,1,\ldots,6),
\]
then $\mathfrak{X}=(\mathcal{R},\{R_i\}_{i=0}^6)$ is a nonsymmetric association scheme of class $6$
with the first eigenmatrix \eqref{P6}.
\end{thm}
\begin{proof}
It is easy to check that $|S_1|=|S_2|=2b(b-1)$, $|S_3|=|S_4|=b$, $|S_5|=b-1$, and $|S_6|=b$.
So we have $p_{0,j}$ ($j=0,1,\ldots,6$) as in \eqref{P6}.
It is trivial that $\lambda_{\alpha}(S_0)=1$ for any $\alpha\in S_0\cup S_1\cup \cdots\cup S_6$.
Hence $p_{j,0}=1$ for $j=0,1,\ldots,6$ as in \eqref{P6}.

We claim that
\begin{align}
\lambda_{\alpha}(S_1)&=\begin{cases}
bi & \text{if} \quad \alpha\in S_1,\\
0 & \text{if} \quad \alpha\in S_3,\\
-2b & \text{if} \quad \alpha\in S_5,\\
0 & \text{if} \quad \alpha\in S_6,
\end{cases} \label{case:1}
\displaybreak[0]\\
\lambda_{\alpha}(S_3)&=\begin{cases}
0 & \text{if} \quad \alpha\in S_1,\\
bi & \text{if} \quad \alpha\in S_3,\\
b & \text{if} \quad \alpha\in S_5,\\
-b & \text{if} \quad \alpha\in S_6,
\end{cases} \label{case:2}
\displaybreak[0]\\
\lambda_{\alpha}(S_5)&=\begin{cases}
-1 & \text{if} \quad \alpha\in S_1,\\
b-1 & \text{if} \quad \alpha\in S_3,\\
b-1 & \text{if} \quad \alpha\in S_5,\\
b-1 & \text{if} \quad \alpha\in S_6.
\end{cases} \label{case:3}
\end{align}

First we prove \eqref{case:1}.
Assume that $\alpha\in S_1$.
Then $\lambda_{\alpha}(\mathcal{T} H)=\lambda(\mathcal{T} H)=bi$ by Lemma~\ref{lem:lamH} (ii),
and $\lambda_{\alpha}(H)=\lambda(\alpha H)$.
Since $S_1\subset\mathcal{R}^{\ast}\setminus\mathcal{E}$, we have $\lambda(\alpha H)=0$ by
\eqref{eq:aH}.
Hence $\lambda_{\alpha}(S_1)=bi$.
Next assume that $\alpha\in S_3$.
Similarly, we have $\lambda_{\alpha}(\mathcal{T} H)=\lambda(\mathcal{T} H)$ and
$\lambda_{\alpha}(H)=\lambda(H)$.
By Lemma~\ref{lem:lamH} (i), (ii) we have $\lambda_{\alpha}(S_1)=0$.
Next assume that $\alpha\in S_5\cup S_6$.
Then by Lemma~\ref{lem:lamaH} and Lemma~\ref{lem:lamaTH},
we have $\lambda_{\alpha}(S_1)=-2b$ or $0$ depending on
$\alpha\in S_5$ or $\alpha\in S_6$.
This completes the proof of \eqref{case:1}.

Next we prove \eqref{case:2}.
Assume that $\alpha\in S_1$.
Then $\lambda_{\alpha}(S_3)=\lambda(\alpha S_3)$.
Since $S_1\subset\mathcal{R}^{\ast}\setminus\mathcal{E}$, we have $\lambda(\alpha S_3)=0$ by
\eqref{eq:aH}.
Next assume that $\alpha\in S_3$.
Then $\lambda_{\alpha}(S_3)=\lambda(S_3)=bi$ by Lemma~\ref{lem:lamH} (i).
Next assume that $\alpha\in S_5\cup S_6$.
Then by Lemma~\ref{lem:lamaH}, $\lambda_{\alpha}(S_3)=b$ or $-b$ depending on
$\alpha\in S_5$ or $\alpha\in S_6$.
This completes the proof of \eqref{case:2}.

Finally we prove \eqref{case:3}.
Assume that $\alpha\in S_1$.
Since $S_1\subset\mathcal{R}^{\ast}\setminus\mathcal{E}$,
we have $\lambda_{\alpha}(S_5)=\lambda(\alpha(\mathcal{P}_0\setminus\{0\}))=-1$ by \eqref{eq:aP0}.
Next assume that $\alpha\in S_3$.
Since $\alpha\mathcal{P}_0=\mathcal{P}_0$, we have $\lambda_{\alpha}(S_5)
=\lambda(\mathcal{P}_0\setminus\{0\})=|S_5|=b-1$.
Next assume that $\alpha\in S_5\cup S_6$.
Then $\alpha(\mathcal{P}_0\setminus\{0\})=\{0\}$.
Hence $\lambda_{\alpha}(S_5)=|S_5|=b-1$.
This completes the proof of \eqref{case:3}.

From \eqref{case:1}--\eqref{case:3}
we have $p_{1,1}$, $p_{3,1}$, $p_{5,1}$, $p_{6,1}$,
$p_{1,3}$, $p_{3,3}$, $p_{5,3}$, $p_{6,3}$,
$p_{1,5}$, $p_{3,5}$, $p_{5,5}$, and $p_{6,5}$ as in \eqref{P6}.
Since $\{S_j\}_{j=0}^6$ is a partition of $\mathcal{R}$,
we have $\sum_{j=0}^6\lambda_{\alpha}(S_j)=0$ for $\alpha\not=0$.
Thus, it is enough to check that
$\lambda_{\alpha}(S_j)$ is a constant independent of $\alpha\in S_i$ for $i=1,\ldots,6$ and $j=1,\ldots,5$.
Since $S_2=-S_1$, we have $p_{2,j}=\overline{p_{1,j}}$ and $p_{j,2}=\overline{p_{j,1}}$ for $j=1,\ldots,6$.
Since $S_4=-S_3$, we have $p_{4,j}=\overline{p_{3,j}}$ and $p_{j,4}=\overline{p_{j,3}}$ for $j=1,\ldots,6$.
Therefore we find all $p_{i,j}$ as in \eqref{P6}
from \eqref{case:1}--\eqref{case:3}.
\end{proof}

Let $\mathfrak{X}=(X,\{R_i\}_{i=0}^d)$ be a commutative association scheme.
A partition $\Lambda_0, \Lambda_1, \ldots, \Lambda_e$ of the index set $\{0,1,\ldots,d\}$
of the association scheme $\mathfrak{X}$
is said to be {\it admissible} if $\Lambda_0=\{0\}$, $\Lambda_i\not=\emptyset$
$(1\leq i\leq e)$ and
$\Lambda_i'=\Lambda_j$ for some $j$ ($1\leq j\leq e$),
where $\Lambda_i'=\{\alpha' \mid \alpha\in\Lambda_i\}$, $R_{\alpha'}=\{(x,y)\mid (y,x)\in R_{\alpha}\}$.
Let $R_{\Lambda_i}=\bigcup_{\alpha\in\Lambda_i}R_{\alpha}$.
If $\mathfrak{Y}=(X,\{R_{\Lambda_i}\}_{i=0}^e)$ becomes an association scheme,
then it is called a {\it fusion} scheme of $\mathfrak{X}$.

\medskip\par\noindent
{\bf Bannai--Muzychuk criterion} (\cite{B, Muzychuk}).
Let $\mathfrak{X}$ be a commutative association scheme
with the first eigenmatrix $P$.
Let $\{\Lambda_j\}_{j=0}^e$ be an admissible partition of the index set $\{0,1,\ldots,d\}$.
Then $\{\Lambda_j\}_{j=0}^e$ gives rise to a fusion scheme $\mathfrak{Y}$
if and only if there exists a partition $\{\Delta_i\}_{i=0}^e$ of $\{0,1,\ldots,d\}$
with $\Delta_0=\{0\}$ such that
each $(\Delta_i,\Lambda_j)$ block of the first eigenmatrix $P$ has a constant row sum.
Moreover, the constant row sum of the $(\Delta_i,\Lambda_j)$ block
is the $(i,j)$ entry of the first eigenmatrix of $\mathfrak{Y}$.

\medskip

Let $\mathfrak{X}$ be an association scheme given in Theorem~\ref{thm:1}.
Fusion schemes of $\mathfrak{X}$ with at least three classes
are listed in Table~\ref{table}.
\begin{table}[htb]
 \begin{center}
  \begin{tabular}{|c|c|c|c|}
  \hline
  &fused relations & class & nonsymmetirc or symmetric \\
  \hline\hline
  $\mathfrak{X}_1$&$\{1,2\}$ & $5$ & nonsymmetric \\
  \hline
  $\mathfrak{X}_2$&$\{3,4\}$ & $5$ & nonsymmetric \\
  \hline
  $\mathfrak{X}_3$&$\{1,2\},\{3,4\}$ & $4$ & symmetric \\
  \hline
  $\mathfrak{X}_4$&$\{3,4,6\}$ & $4$ & nonsymmetric \\
  \hline
  $\mathfrak{X}_5$&$\{1,2\},\{3,4\},\{5,6\}$ & $3$ & symmetric \\
  \hline
  $\mathfrak{X}_6$&$\{1,2,3,4\}$ & $3$ & symmetric \\
  \hline
  $\mathfrak{X}_7$&$\{1,3\},\{2,4\},\{5,6\}$ & $3$ & nonsymmetric \\
  \hline
  $\mathfrak{X}_8$&$\{1,4\},\{2,3\},\{5,6\}$ & $3$ & nonsymmetric \\
  \hline
  \end{tabular}
  \caption{Fusion schemes of $\mathfrak{X}$}
  \label{table}
 \end{center}
\end{table}
The first eigenmatrix of $\mathfrak{X}_7$ and $\mathfrak{X}_8$ in Table~\ref{table}
are the matrix (\ref{d3P1}) by putting $a=b$. Put $b=4^p$.
We can verify that $\mathfrak{X}_7$ and $\mathfrak{X}_8$ for $p=3$ are isomorphic,
and $\mathfrak{X}_7$ and $\mathfrak{X}_8$ for $p=5$ are not isomorphic.
The computation needed to verify these facts was done with the
help of Magma \cite{magma}.

\section{Complex Hadamard matrices and Galois rings}\label{sec:5}

We continue to use the notation introduced in Section~\ref{sec.Galoisring}.
Let $\{A_j\}_{j=0}^6$ be the set of the adjacency matrices of $\mathfrak{X}$.
Then $A_1^T=A_2$, $A_3^T=A_4$, and $A_5,A_6$ symmetric.
We call the algebra $\mathfrak{A}=\langle A_0,A_1,\ldots,A_6\rangle$
the Bose--Mesner algebra of $\mathfrak{X}$.
Our next theorem gives a classification of
hermitian complex Hadamard matrices belonging to $\mathfrak{A}$.

\begin{thm}\label{thm:2}
Let $w_0=1$ and $w_j$ {\rm(}$1\leq j\leq6${\rm)} be complex numbers of absolute value $1$.
Set
\[ W=\sum_{j=0}^6w_jA_j\in\mathfrak{A}, \]
and assume that $W$ is hermitian.
Then, $W$ is a complex Hadamard matrix if and only if
\begin{align}
W=&A_0+\epsilon_1i(A_1-A_2)+\epsilon_2i(A_3-A_4)+A_5+A_6, \quad \text{or} \label{eq:W1}\\
W=&A_0+\epsilon_1i(A_1-A_2)+\epsilon_2(A_3+A_4)+A_5-A_6,\label{eq:W2}
\end{align}
for some $\epsilon_1,\epsilon_2\in\{\pm1\}$.
\end{thm}

\begin{proof}
Recall that $\{A_j\}_{j=0}^6$ is the set of the adjacency matrices of $\mathfrak{X}$ given in Theorem~\ref{thm:1}.
Notice that $A_1^T=A_2$, $A_3^T=A_4$, while $A_5$ and $A_6$ are symmetric.
Since $W$ is hermitian, we have $w_1w_2=1$ and $w_3w_4=1$.

Suppose that the matrix $W$ is a complex Hadamard matrix.
Then $(w_i)_{i=1}^6$ is a common zero of the polynomials $e_k$ ($k=1,\ldots,6$) defined in \eqref{eq:ek}.
Since
\begin{equation}\label{01-7}
e_2-e_1=4ib(X_1-X_2)(X_5-1),
\end{equation}
we have $w_1=w_2$ or $w_5=1$.

Suppose first that $w_1=w_2$. After specializing $X_1=X_2$ and $X_1^2=1$,
we have
$$e_1=(X_5-(2b+1))(X_5+2b-1).$$
Then $|w_5|=|2b\pm1|\geq 2b-1>1$, which is a contradiction.
Therefore, we must have $w_5=1$.
After specializing $X_5=1$, we have
$$e_1=-b^2(X_1-X_2-2i)(X_1-X_2+2i).$$
Hence $w_1=-w_2\in\{\pm i\}$.
Moreover, after specializing $X_5=1$, $X_1=-X_2$ and $X_1^2=-1$, we have
\begin{align}
e_4-e_3&=4ib^2(X_6-1)(X_3-X_4),\label{01-8}\\
e_5-e_6&=4b^2(X_6+1)(X_3+X_4).\label{01-9}
\end{align}
If $w_6=1$, then $w_4=-w_3\in\{\pm i\}$ by \eqref{01-9}.
Therefore we have \eqref{eq:W1}.
If $w_3=w_4$, then $w_3\in\{\pm1\}$ and $w_6=-1$ by \eqref{01-9}.
Therefore we have \eqref{eq:W2}.

Conversely, assume that $W$ is one of the matrices \eqref{eq:W1}, \eqref{eq:W2}.
We show that $W$ is a complex Hadamard matrix.
By Lemma~\ref{lem:equiv}, it suffices to show that
$(w_i)_{i=1}^6$ is common zero of the polynomials \eqref{eq:ek}, and
it is easy to do this.
\end{proof}

We note that the matrix \eqref{eq:W1}
belongs to the Bose--Mesner algebra of  $\mathfrak{X}_7$ or $\mathfrak{X}_8$  in Table~\ref{table},
depending on $\epsilon_1=\epsilon_2$ or $\epsilon_1=-\epsilon_2$.
Also, the matrix \eqref{eq:W2}
belongs to the Bose--Mesner algebra of $\mathfrak{X}_2$ or $\mathfrak{X}_4$ in Table~\ref{table},
depending on $\epsilon_2=1$ or $\epsilon_2=-1$.
Therefore, no proper fusion scheme of $\mathfrak{X}$ contains
all the matrices \eqref{eq:W1} and \eqref{eq:W2} in its Bose--Mesner algebra $\mathfrak{A}$.
We also note that Ma \cite{Ma} considered association schemes which are
invariant under the multiplication by $\mathcal{T}$.
The only association schemes
in Table~\ref{table} which are invariant under the multiplication by $\mathcal{T}$
is $\mathfrak{X}_7$, and it has the first eigenmatrix described by
\cite[Theorem 7]{Ma}.

redpajama_set_name: RedPajamaArXiv